text:
# References

CHARMED HADRONS PRODUCTION IN HIGH-ENERGY $`\mathrm{\Sigma }^{}`$ BEAM

A.K. Likhoded1 and S.R. Slabospitsky

State Research Center

Institute for High Energy Physics,

Protvino, Moscow Region 142284,

RUSSIA

Abstract

We present the calculation of the inclusive $`x_F`$-distributions of charmed hadrons, produced in high-energy $`\mathrm{\Sigma }^{}`$-beam. Our calculation is based on the modified mechanism of charmed quarks fragmentation as well as on the mechanism of $`c`$-quark recombination with the valence quarks from initial hadrons.

1E–mail: LIKHODED$`\mathrm{@}`$mx.ihep.su,

Perturbative QCD provides a reasonable description of the experimental data on the inclusive cross sections of open charm and beauty production on fixed targets . However, it is well known that the experiments indicate that there is a substantial difference in the production of the charmed and anticharmed hadrons in the fragmentation region of the initial hadrons (a leading particle effect). This leading particle asymmetry $`A`$ is defined as follows:

$`A{\displaystyle \frac{\sigma (\mathrm{leading})\sigma (\mathrm{nonleading})}{\sigma (\mathrm{leading})+\sigma (\mathrm{nonleading})}}.`$ (1)
It worth to note, that perturbative QCD calculation  is unable to reproduce this effect. Indeed, in the quark parton model framework the production of hadrons containing a heavy quark proceeds via two subsequent stages

* heavy quark pair $`Q\overline{Q}`$ is produced as a result of the hard collision of the partons from the initial hadrons (e.g. the subprocesses $`ggc\overline{c}`$ and $`q\overline{q}c\overline{c}`$ in Born approximation);
* transition of the heavy quarks $`c`$ into charmed hadrons (“hadronization”).

The standard way of describing heavy quark hadronization is to use the fragmentation function $`D(z)`$ of the heavy $`c`$-quark into the charmed hadron ($`D`$-meson or baryon) (here $`z=|\stackrel{}{p}_D|/|\stackrel{}{p}_c|`$ is the fraction of the heavy quark momentum carried away by the charmed hadron $`D`$).

It should be noted that the usage of the fragmentation function assumes the absence of the interaction of the produced heavy quark $`Q`$ with the remnants of the initial hadrons. Therefore, it should be no difference between the spectra of charmed and anticharmed hadrons. Moreover, any modification of the fragmentation mechanism can not reproduce the production asymmetry (the leading particle effect).

Note, that the fragmentation mechanism can be apply for the production of the $`c\overline{c}`$ pair in the color–singlet state or for high $`p_T`$ production of the open charm. On the other hand, for the case of the hadronic production of color $`c\overline{c}`$ pair with small $`p_T`$ one should takes into account the possibility of charmed $`c`$ and $`\overline{c}`$ quarks interaction with the initial hadron remnants. Therefore, due to the different valence quarks in the initial hadrons one may expect the different inclusive spectra of the final charmed hadrons.

In the parton model framework, a heavy $`c`$–quark should interact with a high probability with its nearest neighbor in the rapidity space able to form a color-singlet state with it. In some cases, the heavy antiquark may find itself close (in rapidity space) to a valence light quark from the initial hadron. This would result in the formation of a fast heavy meson in the fragmentation region of the initial hadron. Alternatively, the proximity of a heavy quark to a valence diquark results in the fast charmed $`B(cq_1q_2)`$-baryon production.

Therefore, the ”hard” part of charmed hadron spectra is very sensitive to the form of valence quark distribution in the initial hadrons.

In this note we consider the charmed hadron production in high-energy beam of $`\mathrm{\Sigma }^{}`$-hyperon. We may expect the different behavior of the distributions of the valence $`d`$\- and $`s`$-quarks. As a result we should observe the different $`x_F`$-dependence of spectra of charmed hadrons with $`d`$\- or $`s`$-quarks, namely, $`D^{}(\overline{c}d)`$ and $`D_s^{}(\overline{c}s)`$, $`\mathrm{\Xi }_c^0(cds)`$ and $`\mathrm{\Sigma }_c^{}(cdd)`$, etc.

Indeed, very roughly, the distribution of valence quark in the baryon $`B(q_1q_2q_3)`$ can be presented as follows :

$$V_{q_1}^B(x)x^{\alpha _1}(1x)^{\gamma _b\alpha _2\alpha _3},$$
(2)
where $`\alpha _i`$ is the intercept of the leading Regge-trajectory for $`q_i`$-quark, while $`\gamma _B4`$. Note, that due to violation flavor $`SU(N)`$-symmetry, we have different intercepts for $`d(u)`$\- and $`s`$-quarks :

$`\alpha _u=\alpha _d={\displaystyle \frac{1}{2}},\alpha _s0,\alpha _c2.2.`$ (3)
As a result, the $`x`$-dependence of the valence $`d`$ and $`s`$-quark in the $`\mathrm{\Sigma }^{}(sdd)`$-hyperon has the form as follows :

$$V_d^\mathrm{\Sigma }\frac{1}{\sqrt{x}}(1x)^{3.5},V_s^\mathrm{\Sigma }(1x)^3$$
(4)
It is seen from the (4) that the valence $`s`$-quark in the $`\mathrm{\Sigma }^{}`$-hyperon has slightly harder $`x`$-distribution than that for $`d`$-quark.

We use the model  to describe the production asymmetry for charmed hadrons. In this model the interaction of the charmed quarks with valence quarks from the initial hadrons describes with the help of the recombination function . The detail description of this mechanism is given elsewhere . The recombination of $`q_V`$ and $`\overline{c}`$ quarks into $`D`$-meson is described by the function of $`R_M(x_V,z;x)`$:

$`R_M(x_q,z;x)=Z_M\xi _q^{(1\alpha _q)}\xi _c^{(1\alpha _c)}\delta (1\xi _q\xi _c),`$ (5)
where $`\xi _q=x_q/x`$ and $`\xi _c=z/x`$, while $`x_q`$, $`z`$, and $`x`$ are the fractions of the initial-hadron c.m. momentum that are carried away by the valence $`q`$-quark, charmed $`c`$-quark, and the meson $`M_{\overline{c}q}`$, respectively. The corresponding recombination of three quarks into baryon can be described by means of the similar recombination function:

$`R_B(x_1,x_2,z;x)=Z_M\xi _1^{(1\alpha _1)}\xi _2^{(1\alpha _2)}\xi _c^{(1\alpha _c)}\delta (1\xi _1\xi _2\xi _c).`$ (6)
These functions take into account the momentum conservation and the proximity of partons in the rapidity space. Actually, the recombination function is the modulus squared of the heavy meson wave function in momentum space, being considered in the infinite momentum frame in the valence quark approximation.

As a result, the total differential cross section for the production of $`H_c`$–hadron looks as follows:

$`{\displaystyle \frac{d\sigma (H_c)}{dx}}={\displaystyle \frac{d\sigma ^R(H_c)}{dx}}+{\displaystyle \frac{d\sigma ^F(H_c)}{dx}},`$ (7)
where the first term in the r.h.s. is the cross section for $`H_c`$–hadron production due to the recombination of the charmed $`c`$–quark with valence quarks from initial hadron, while the second term is the cross section for $`H_c`$ production due to charmed quark fragmentation.

Note, that this model provides more or less successful description of the charmed $`D`$-meson production in $`\pi ^{}N`$-interactions (see Fig. 1, 2 and  for details).

We use LO formulas for the cross sections of quark-antiquark and gluon-gluon annihilation into charmed quark pair. We set $`m_c=1.25`$ GeV, $`\alpha _s=0.3`$ and find the following cross section value of $`\mathrm{\Sigma }^{}p`$ interaction at $`P_{LAB}=600`$ GeV:

$$\sigma (\mathrm{\Sigma }^{}pc\overline{c}X)8\mu b$$
(8)
In our calculations we do not pretend to reproduce the absolute value of this cross section (see , for detail consideration of this problem). We concentrate on the description of $`x_F`$-distribution of charmed mesons and baryons. The corresponding distributions (integrated over $`p_{}`$) are presented in Fig.3–5. We may see from these figures, that the considered charmed quark interaction in the final state (recombination) leads, indeed, to noticeable differences in $`x_F`$-spectra. These differences can be explicitly seen in Fig.6, where we present the corresponding asymmetry $`A`$ (see (1 for definition). The most non-trivial prediction of the proposed model is presented in two lower plots in Fig.6, where we present the ratio of the inclusive spectra of $`D_s^{}(\overline{c}s)`$ and $`D^{}(\overline{c}d)`$ mesons. It is evident from this figure the difference of $`x_F`$-spectra of the these two mesons, which is a result of of different $`x`$-distributions of the valence $`d`$-and $`s`$-quarks in the initial $`\mathrm{\Sigma }^{}`$-beam (see (4)).

Conclusion
In the present note we wish to stress once more, that the source of the observed asymmetry in charmed hadron production is the interaction of produced charmed quarks with valence quarks from initial hadrons. Note, the model under consideration provides also the additional method to measure the valence quark distribution functions of $`K`$-meson and $`\mathrm{\Sigma }`$-baryons.


format: text

text:
# Viscosity critical behaviour at the gel point in a 3𝑑 lattice model

## Abstract

Within a recently introduced model based on the bond-fluctuation dynamics we study the viscoelastic behaviour of a polymer solution at the gelation threshold. We here present the results of the numerical simulation of the model on a cubic lattice: the percolation transition, the diffusion properties and the time autocorrelation functions have been studied. From both the diffusion coefficients and the relaxation times critical behaviour a critical exponent $`k`$ for the viscosity coefficient has been extracted: the two results are comparable within the errors giving $`k1.3`$, in close agreement with the Rouse model prediction and with some experimental results. In the critical region below the transition threshold the time autocorrelation functions show a long time tail which is well fitted by a stretched exponential decay.
PACS: 05.20.-y, 82.70.Gg, 83.10.Nn

Polymeric materials are characterized by a rich and complex phenomenology, which goes from non newtonian dynamic behaviour in polymer liquids to viscoelastic properties of polymer gels and the glass transition of polymer melts, intensively investigated in the last decades. The interest in such systems has been recently further increased by several possibilities of technological applications in many different fields but still the non trivial viscoelastic behaviour of polymeric systems does not have a completely satisfying description.
The gelation transition which transforms the polymeric solution, i.e. the sol, in a polymeric gel is characterized by a dramatic change in the viscoelastic properties: this is usually described in terms of the divergence of the viscosity coefficient and the appearance of an elastic modulus which characterizes the gel phase . In the experiments both the viscosity coefficient and the elastic modulus dependence on the polymerization reaction extent are well fitted by a power law but the experimental determination of these critical exponents is quite controversial: the results are in fact scattered, probably because of the practical difficulties in obtaining the gelation transition in a reproducible manner, and could only be interpreted on the basis of a better comprehension of the relevant mechanisms in the transition. Recent experimental measurements of the viscosity critical exponent $`k`$ give values ranging from $`0.7`$ in diisocyanate/triol to $`1.5`$ in epoxy resins whereas for the elastic modulus critical exponent $`t`$ the values are even more scattered, ranging from $`1.9`$ in diisocianate/triol to $`3.0`$ in polyesters (see references ).
With simple statistical mechanics models it is possible to analyse the essential aspects of the transition and its critical properties. From the Flory model this has led to the description of the gelation process in terms of a percolation transition. The percolation model has turned out to be the satisfactory model for the sol-gel transition , it is able to describe the role of connectivity and gives the critical exponents for all the geometrical properties, which are in perfect agreement with the experimental results. On the other hand the viscoelastic dynamic behaviour is not simply obtained in terms of the connectivity transition.
The difficulties in studying the viscosity critical behaviour at the gel point come from the determination of the viscosity of a very complex medium, which is the sol at the transition threshold, a highly polydisperse polymeric solution at high concentration. The complex polymer dynamics is characterized by the relaxation processes over many different time scales and compete with the increasing connectivity to produce the observed viscoelastic behaviour. The simplest approach consists in considering the sol as a polydisperse suspension of solid spheres neglecting the cluster-cluster interactions and generalizing the Einstein formula for the viscosity of a monodisperse suspension of solid spheres , which corresponds to a highly diluted regime. Within the Flory classical theory of gelation the viscosity remains finite  or diverges at most logarithmically. Using instead the Rouse model for the polymer dynamics , which neglects the entanglement effects and the hydrodynamics interactions, the viscosity in a solution of polymeric clusters, expressed in terms of the macroscopic relaxation time, grows like $`<R^2>`$ as the cluster radius $`R`$ grows in the gelation process. The contribution of the $`n_s`$ molecules of size $`s`$ and gyration radius $`R_s`$ to the average $`<>`$ is of the order of $`sn_sR_s^2`$ leading to the critical exponent $`k=2\nu \beta `$ , where $`\nu `$ is the critical exponent of the correlation length diverging at the gel point and $`\beta `$ is the critical exponent describing the growth of the gel phase. With the random percolation exponents the value $`k1.35`$ is found, that agrees quite well with the experimental measurements for silica gels of ref.. Actually this Rouse exponent could be considered as an upper limit due to the complete screening of the hydrodynamic interactions and the entanglement effects. The Zimm approach , where the monomer correlation due to the hydrodynamic interactions are not completely screened, would give a smaller exponent.
Another approach has been proposed by de Gennes  using an analogy between the viscosity at the gelation threshold and the diverging conductivity in the random superconducting network model. Following this analogy an exponent $`k0.7`$ is obtained in $`3d`$, according to the determination of the conductivity critical exponent in the random superconducting model . This result is in good agreement with the values experimentally obtained in gelling solutions of polystyrene/divinylbenzene and diisocyanate/triol .
Our approach consists in directly investigating the viscoelastic properties at the sol-gel transition, introducing within the random percolation model the bond fluctuation ($`BF`$) dynamics which is able to take into account the polymer conformational changes . We study a solution of tetrafunctional monomers at concentration $`p`$ and with a probability $`p_b`$ of bonds formation. In terms of these two parameters one has different cluster size distributions and eventually a percolation transition. Monomers interact via excluded volume interactions and can diffuse with local random movements. The monomer diffusion process produces a variation of the bond vectors and is constrained by the excluded volume interaction and the SAW condition for polymer clusters: these two requirements can be satisfied if the bond lengths vary within an allowed range. This dynamics results to take into account the main dynamic features of polymer molecules.
We have performed numerical simulations of the model on the cubic lattice of different lattice sizes ($`L=24,32,40`$) with periodic boundary conditions. The eight sites which are the vertices of a lattice elementary cell are simultaneously occupied by a monomer, with the constraint that two nearest neighbour ($`nn`$) monomers are always separated by an empty elementary cell, i.e. two occupied cells cannot have common sites. The lattice of cells, with double lattice spacing, has been occupied with probability $`p`$, which coincides with the monomer concentration on the main lattice in the thermodynamic limit . Monomers are randomly distributed on the main lattice via a diffusion process, then between two $`nn`$ or next nearest neighbour ($`nnn`$) monomers bonds are instantaneously created with probability $`p_b`$ along lattice directions. Since most of the experimental data on the gelation transition refer to polymers with monomer functionality $`f=4`$ we have considered this case allowing the formation of at most four bonds per monomer.
First, a qualitative phase diagram has been determined, studying the onset of the gel phase varying $`p`$ and $`p_b`$: on this basis we have then fixed $`p_b=1`$ and let $`p`$ vary in the interesting range from the sol to the gel phase. The percolation transition has thus been studied via the percolation probability $`P`$ and the mean cluster size $`\chi `$ on lattices of different size. From their finite size scaling behaviours we have obtained the percolation threshold $`p_c0.718\pm 0.005`$, the critical exponents $`\nu 0.89\pm 0.01`$ for the percolation correlation length $`\xi `$ and $`\gamma 1.8\pm 0.05`$ for the mean cluster size $`\chi `$ (fig.(1)). These results do agree with the random percolation critical exponents . ***The same agreement with the random percolation critical exponents has already been obtained in the $`d=2`$ version of the model . Here it is worth to mention that this case on a cubic lattice with monomer functionality $`f=4`$ is rather the problem of a restricted valence percolation. This is a percolation on a lattice where the number of bonds emanating from the same site is restricted or no site may have more than a fixed number of nn. It does reproduce the occurrence of valence saturation for monomers in the gelation process. If the number of allowed bonds per site is greater than $`2`$ this problem is expected to belong to the same universality class as the random percolation . As the restriction on the number of bonds per monomer clearly introduces a correlation effect in the process of bond formation we have optimized our algorithm to minimize this effect and actually obtained a good agreement with the random percolation exponents.
The system evolves according to the bond fluctuation dynamics, i.e. the monomers diffuse with random local movements along lattice directions within the excluded volume constraint and produce the bond length fluctuation among the allowed values. In fact the $`BF`$ dynamics can be easily expressed in terms of a lattice algorithm and on a cubic lattice it can be shown that the bond lengths which guarantee the SAW condition are $`l=2,\sqrt{5},\sqrt{6},3,\sqrt{10}`$ in lattice spacing units . In fig.(2) a simple example of time evolution of a cluster is shown.
We remind that the percolation properties are not modified during the dynamic evolution of the system, as it is not possible to break or form bonds, which would change the cluster size distribution.
In order to determine the viscosity critical behaviour we use two independent ways, based respectively on the diffusion behaviour of the clusters and on the relaxation times. Within the study of diffusion properties in the system an interesting picture is obtained with a simple scaling argument on the diffusion coefficients as presented in : the sol at the sol-gel transition is a heteregenous medium formed by the solvent and all the other clusters of different sizes, with the mean cluster size rapidly growing near the percolation threshold. A cluster with gyration radius $`R`$ can be seen as a probe diffusing in such a medium: as long as its radius is much greater than the value of the percolation correlation length in the sol the Stokes-Einstein relation is expected to be valid, and the diffusion coefficient $`D(R)`$ of the probe will decrease proportionally to the inverse of the viscosity coefficient of the medium, $`D(R)1/R^{d2}\eta `$. Then the generic probe of size $`R`$ diffuses in a medium with a viscosity coefficient depending on $`R`$, $`\eta (R)`$, and a Stokes-Einstein generalized relation $`D(R)1/R^{d2}\eta (R)`$ should be expected to hold.
As the percolation threshold is approached the probe diffuses in a medium where a spanning cluster appears with a self-similar structure of holes at any length scale. At the percolation threshold the viscosity coefficient of the sol (the bulk viscosity coefficient) diverges as $`\eta (p_cp)^k`$ and for the viscosity coefficient depending on $`R`$ the scaling behaviour $`\eta (R)R^{k/\nu }`$ should be expected. When $`R`$ is of the order of the correlation length then $`\eta (R)\eta `$. Following this scaling argument $`D(R)1/R^{d2+k/\nu }`$ at $`p=p_c`$.
Within this description the use of the Rouse model would consist in taking $`D(s)1/s`$ for a cluster of size $`s`$. Then for large enough cluster sizes $`s`$ in percolation $`sR^{d_f}`$ where $`d_f`$ is the fractal dimensionality of the percolating cluster. Taking $`D(R)1/R^{d_f}`$ leads to $`k(d_f+2d)\nu `$, which again gives $`k=2\nu \beta `$, the Rouse exponent given in .
We then study the diffusion of monomers and clusters via the mean square displacement of the center of mass. For a cluster of size $`s`$ (an $`s`$-cluster) this quantity is calculated from the coordinates of its center of mass $`\stackrel{}{R}_s(t)`$ as

$$\mathrm{\Delta }R_s^2(t)=\frac{1}{N_s}\underset{\alpha =1}{\overset{N_s}{}}(\stackrel{}{R}_s^\alpha (t)\stackrel{}{R}_s^\alpha (0))^2$$
(1)
where the index alpha refers to the $`\alpha `$th $`s`$-cluster and $`N_s`$ is the number of $`s`$-clusters so that this quantity is averaged over all the $`s`$-clusters. All the data refer to a lattice size $`L=32`$ and the calculated quantities have been averaged over $`30`$ different site and bond configurations of the same $`(p,p_b)`$ values .
On the basis of the theory of Brownian motion and a simple Rouse model approach  the center of mass of a polymeric molecule is expected to behave as a Brownian particle after a sufficiently long time. Indeed we find a linear dependence on time in the long time behaviour of $`\mathrm{\Delta }R_s^2(t)`$ and determine the diffusion coefficient of an $`s`$-cluster in the environment formed by the solvent and the other clusters. In fig.(3) it is shown how the asymptotic diffusive behaviour is reached after a time which increases with the cluster size $`s`$, due to the more complicated relaxation mechanism linked to the inner degrees of freedom in the molecule. As $`p`$ increases towards $`p_c`$ the $`s`$-clusters move in a medium which is more viscous and whose structure is more and more complex. As a consequence, we observe an immediate increase in the time necessary to reach the asymptotic diffusive behaviour of the center of mass motion for all the cluster sizes.
At $`p=p_c`$ we have calculated the diffusion coefficients of clusters of size $`s`$ and radius of gyration $`R`$ in order to obtain the dependence $`D(R)`$. The diffusion coefficients decrease with the increasing size of the cluster, as it is expected, but after a gradual decrease their values dramatically go to zero. This behaviour which is due to the block of the diffusion for finite cluster size in a finite system has allowed to consider only cluster size up to $`s=30`$. In fig.(4) we have plotted $`D(R)`$ for different cluster sizes ($`s=5`$ to $`30`$): the data results to be well fitted by a power law behaviour. Using this scaling argument which gives the prediction $`D(R)1/R^{1+k/\nu }`$ we obtain for the viscosity coefficient a critical exponent $`k1.3\pm 0.1`$.
We here briefly mention that on the other side the diffusion coefficients of very small clusters are not expected to be linked to the macroscopic viscosity, and in fact the diffusion coefficient of monomers does not go to zero at $`p_c`$, but has a definitively non-zero value for $`p>p_c`$ (fig.(5)). It is also interesting to notice that the data seems to agree with a dependence $`e^{1/1p}`$, suggesting some cooperative mechanism in the diffusion process.
In order to study the viscosity in the system independently from the scaling hypothesis given before, we have studied the relaxation times via the density time autocorrelation functions. We have calculated the time autocorrelation function $`g(t)`$ of the number of pairs of $`nn`$ monomers $`\epsilon (t)`$, defined as

$$g(t)=\frac{\overline{\epsilon (t^{})\epsilon (t^{}+t)}\overline{\epsilon (t^{})}^2}{\overline{\epsilon (t^{})^2}\overline{\epsilon (t^{})}^2}$$
(2)
where the bar indicates the average over $`t^{}`$ (of the order of $`10^3`$ time intervals) and the brackets indicate the average over about $`30`$ different initial site and bond configurations.
At different $`p`$ values in the critical region after a fast transient $`g(t)`$ decays to zero but can not be fitted by a simple time exponential behaviour. This is a sign of the existence of a distribution of relaxation times which cannot be related to a single time and this behaviour had already been observed in the $`d=2`$ study of the model. It is a typical feature of polymeric systems where the relaxation process always envolves the rearrangement of the system over many different length scales . This idea of a complex relaxation behaviour is further confirmed by the good fit of the long-time decay of the $`g(t)`$ with a stretched exponential law (fig.(6)). This behaviour of the relaxation functions is considered typical of complex materials and usually interpreted in terms of a very broad distribution of relaxation times or eventually an infinite number of them and it is in fact experimentally observed in a sol in the gelation critical regime . The picture we obtain via the density time autocorrelation functions is coherent and very close to the experimental characterization of the sol at the gelation threshold.
For $`pp_c`$ we have then fitted the $`g(t)`$ with a stretched exponential behaviour $`e^{(t/\tau _0)^\beta }`$ (fig.(7)), where $`\beta 0.3`$. This value is quite lower than the ones experimentally obtained in ref. for a gelling solution or analitically predicted in ref. for randomly branched polymers: this discrepancy could be due to the fact that our data refer to a quite narrow region near the gelation threshold where $`\beta `$ is expected to assume the lowest value. On the other hand this value of $`\beta `$ agrees with the experimental results in ref. and with the asymptotic value of the stretched exponential exponent $`\beta `$ in both Ising spin glasses and polymer melts close to the freezing point in ref..
The characteristic time $`\tau _0`$ varies with $`p`$ and increases as $`p_c`$ is approached. Plotting $`\tau _0`$ as a function of $`(p_cp)`$ the data can be fitted by a power law dependence, with a critical exponent $`1.27\pm 0.05`$ (fig.(8)), very close to the viscosity critical exponent obtained from the diffusion properties. This characteristic time extracted from the long time relaxation behaviour diverges at the percolation threshold because of the passage between two different viscoelastic regimes.
The most immediate way to characterize the distribution of the relaxation times in the system is the average characteristic time defined as

$$\tau (p)=\frac{_0^tt^{}g(t^{})𝑑t^{}}{_0^tg(t^{})𝑑t^{}}$$
(3)
which is a typical macroscopic relaxation time and can be directly linked to the viscosity coefficient.
Numerically in the eq.(3) $`t`$ has been chosen by the condition $`g(t^{})0.001`$ for $`t^{}t`$. This characteristic time grows with $`p`$ and diverges at the percolation threshold according to the critical behaviour

$$\tau (p_cp)^k$$
(4)
with an exponent $`k1.3\pm 0.03`$ (fig.9), which gives the critical exponent for the viscosity at the sol-gel transition. From an immediate comparison between $`k`$ and the critical exponent for $`\tau _o`$ there seems to be a unique power law characterizing the divergence of the relaxation times in the system approaching from below the transition threshold.
For a more detailed description of the distribution of relaxation times at the sol-gel transition together with a study of the frequency dependence of the viscoelastic properties see references .
The critical exponent obtained for the viscosity critical behaviour is then $`k1.3`$ and it well agrees with the value of $`k`$ previously given, although independently obtained. This value is quite close to the one experimentally measured in silica gels  and interestingly these systems are characterized by a polyfunctional condensation mechanism with tetrafunctional monomers, which is actually very similar to the case simulated here. This result is also close to the value obtained by recent accurate measurements in PDMS .
It does not agree with the random superconducting network exponent $`k=0.7`$ according to the de Gennes’ analogy, whereas it is quite close to the Rouse exponent discussed above. We have already mentioned that a Rouse-like description of a polymer solution corresponds to a complete screening of the entanglement effects and the hydrodynamic interactions and is usually considered not realistic enough. Actually the entanglement effects could be not so important in the relaxation mechanism in the sol on the macroscopic relaxation time scale: due to the fractal structure of the gel phase the system is in fact quite fluid, probably there is no blocking entanglement yet and such temporary entanglements relaxe on a smaller time scale, not really affecting the macroscopic relaxation time . Furthermore the screening effect of the hydrodynamic interactions in a polymeric solution in the semidiluted regime can be quite strong, drastically reducing the range of the interactions so that the Rouse model results to be in fact very satisfactory . This could reasonably be the case of the sol at the gelation threshold too, and the deviation of the real critical exponent from the Rouse value would turn out to be actually very small.
The numerical simulations have been performed on the parallel Cray-t3e system of CINECA (Bologna,Italy), taking about 15000 hours of CPU time. The authors aknowledge partial support from the European TMR Network-Fractals c.n.FMRXCT980183 and from the MURST grant (PRIN 97). This work was also supported by the INFM Parallel Computing and by the European Social Fund.


text:
# Limits on a Composite Higgs Boson

(March, 2000)
Precision electroweak data are generally believed to constrain the Higgs boson mass to lie below approximately 190 GeV at 95% confidence level. The standard Higgs model is, however, trivial and can only be an effective field theory valid below some high energy scale characteristic of the underlying non-trivial physics. Corrections to the custodial isospin violating parameter $`T`$ arising from interactions at this higher energy scale dramatically enlarge the allowed range of Higgs mass. We perform a fit to precision electroweak data and determine the region in the $`(m_H,\mathrm{\Delta }T)`$ plane that is consistent with experimental results. Overlaying the estimated size of corrections to $`T`$ arising from the underlying dynamics, we find that a Higgs mass up to 500 GeV is allowed. We review two composite Higgs models which can realize the possibility of a phenomenologically acceptable heavy Higgs boson. We comment on the potential of improvements in the measurements of $`m_t`$ and $`M_W`$ to improve constraints on composite Higgs models.

Precision electroweak data are generally believed to constrain the Higgs boson mass to lie below approximately 190 GeV at 95% confidence level . The standard Higgs model is, however, trivial  and can only be an effective field theory valid below some high energy scale $`\mathrm{\Lambda }`$ characteristic of the underlying non-trivial physics. Additional interactions coming from the underlying theory, and suppressed by the scale $`\mathrm{\Lambda }`$, give rise to model-dependent corrections to measured electroweak quantities. When potential corrections from physics at higher energy scales are included, the limit on the Higgs boson mass becomes weaker111See also, Langacker and Erler in . .

In the context of the triviality of the standard model and given the relatively weak (logarithmic) dependence of electroweak observables on the Higgs boson mass , the typical size of corrections to $`T`$ arising from custodial symmetry violating  non-trivial underlying physics can dramatically enlarge222In contrast, in theories lacking a custodial symmetry the contributions to $`S`$ are relatively small  and do not have a significant effect on Higgs mass bounds. the allowed Higgs mass range . In this note we perform a fit to precision electroweak data and determine the region in the $`(m_H,\mathrm{\Delta }T)`$ plane that is consistent with experimental results. Overlaying the predicted size of corrections to $`T`$ arising from the underlying dynamics, we find that a Higgs mass up to 500 GeV is allowed. We review two composite Higgs models which can realize the possibility of a phenomenologically acceptable heavy Higgs boson.

For a given Higgs boson mass, an upper bound on the scale $`\mathrm{\Lambda }`$ is given by the position of the Landau pole  of the Higgs boson self-coupling $`\lambda `$. As the Higgs boson mass is proportional to $`\lambda (m_H)`$, the larger the Higgs boson mass the smaller the upper bound on the scale $`\mathrm{\Lambda }`$. We may estimate333While this estimate is based on perturbation theory, non-perturbative calculations yield essentially the same result . this upper bound by integrating the one-loop beta function for the self-coupling $`\lambda `$, which yields

$$\mathrm{\Lambda }\stackrel{<}{}m_H\mathrm{exp}\left(\frac{4\pi ^2v^2}{3m_H^2}\right),$$
(1.1)
where $`m_H`$ is the Higgs boson mass and $`v246`$ GeV is the vacuum expectation value of the Higgs boson.

The leading corrections to electroweak observables from the underlying theory are encoded in dimension six operators  which contribute to the Peskin-Takeuchi $`S`$ and $`T`$ parameters . Given the scale of the underlying non-trivial physics, dimensional analysis  may be used to estimate the size of effects from these dynamics in the low-energy Higgs theory . If the underlying theory does not respect custodial symmetry , the contribution to $`T`$ is dominant and is estimated to be

$$|\mathrm{\Delta }T|\frac{b\kappa ^2v^2}{\alpha _{em}(M_Z^2)\mathrm{\Lambda }^2},$$
(1.2)
or larger . Here $`\alpha _{em}`$ is the electromagnetic coupling renormalized at $`M_Z^2`$, $`b`$ is a model-dependent coefficient of order 1, and $`\kappa `$ is a measure of the size of dimensionless couplings in the effective Higgs theory and is expected to lie between 1 and $`4\pi `$ . Combining eqn. 1.2 with the bound on $`\mathrm{\Lambda }`$ shown in eqn. 1.1, we find

$$|\mathrm{\Delta }T|\stackrel{>}{}\frac{b\kappa ^2v^2}{\alpha _{em}(M_Z^2)m_H^2}\mathrm{exp}\left(\frac{8\pi ^2v^2}{3m_H^2}\right).$$
(1.3)
Since the Higgs model is trivial, the potential effects of the underlying non-trivial dynamics must be included when establishing constraints on the Higgs mass . As the contributions to $`T`$ are expected to dominate, we have performed a fit to electroweak measurements  and have determined the region in the $`(m_H,\mathrm{\Delta }T)`$ plane that is consistent with these results. In addition to measurements at the $`Z`$-pole from LEP and SLD, we include measurements of $`M_W`$ from LEP and the Tevatron, and measurements of $`m_t`$ from the Tevatron. In performing these fits, we have used ZFITTER 6.21  to generate the standard model predictions for a given value of the $`Z`$ mass, Higgs mass, top-quark mass, and strong ($`\alpha _s`$) and electromagnetic ($`\alpha _{em}`$) couplings, and have introduced the effect of non-zero $`\mathrm{\Delta }T`$ linearly . We have included the determinations of $`\alpha _{em}`$

$$\alpha _{em}^1(M_Z^2)=128.905\pm 0.036$$
(1.4)
and the (non-electroweak) determinations of $`\alpha _s`$

$$\alpha _s(M_Z^2)=0.119\pm 0.002$$
(1.5)
as observations, i.e. we have included deviations from the listed central values in our computation of $`\chi ^2`$. The correlation matrices listed in ref.  are incorporated in our calculation of $`\chi ^2`$.

The result of our fit is summarized in Figure 1. The best-fit value444This best fit value is, of course, below the direct experimental lower bound  of order 108 GeV. is shown and occurs at a Higgs boson mass of 90 GeV; it corresponds to a minimum value of $`\chi ^2=21.7`$ for 21 observables while varying 5 fit parameters ($`M_Z`$, $`m_H`$, $`m_t`$, $`\alpha _s`$, and $`\alpha _{em}`$). For two degrees of freedom, the 68% and 95% CL bounds correspond to $`\mathrm{\Delta }\chi ^2`$ of 2.30 and 6.17 respectively. The two degree555Note that the one degree of freedom 95% CL upper bound on $`m_H`$, $`\mathrm{\Delta }\chi ^2=4`$ and $`\mathrm{\Delta }T=0`$, is approximately 190 GeV in agreement with . of freedom 95% CL upper bound on $`m_H`$ is 243 GeV for $`\mathrm{\Delta }T=0`$.

Extending this bound to non-zero $`\mathrm{\Delta }T`$, we see that the region in the $`(m_H,\mathrm{\Delta }T)`$ plane which fits the observed data as well as the “standard model” at 95% CL extends to large Higgs masses for a positive value of $`\mathrm{\Delta }T`$.

It is not possible, however, for this entire region to be realized consistent with the constraints of triviality. For example, motivated by the models we consider below, the area excluded by eqn. 1.3 with $`b\kappa ^2=4\pi `$ is shown as the light region on the right in Figure 1. Overlaying the constraints, we see that Higgs masses above 500 GeV would likely imply the existence of new physics at such low scales ($`\mathrm{\Lambda }\stackrel{<}{}12`$ TeV from eqn. 1.1) as to give rise to a contribution to $`T`$ which is too large .

We emphasize that these estimates are based on dimensional arguments, and we are not arguing that it is impossible to construct a composite Higgs model consistent with precision electroweak tests with $`m_H`$ greater than 500 GeV. Rather, barring accidental cancellations in a theory without a custodial symmetry, contributions to $`\mathrm{\Delta }T`$ consistent with eqn. 1.3 are generally to be expected. This expectation is illustrated in the two models which we now discuss.

The top-quark seesaw theory of electroweak symmetry breaking  provides a simple example of a model with a potentially heavy composite Higgs boson consistent with electroweak data. In this case, electroweak symmetry breaking is due to the condensation, driven by a strong topcolor  gauge interaction, of the left-handed top-quark with a new right-handed singlet fermion $`\chi `$. Such an interaction gives rise to a composite Higgs field at low energies, and the mass of the top-color gauge boson sets the scale of the Landau pole $`\mathrm{\Lambda }`$ . The weak singlet $`\chi _L`$ and $`t_R`$ fields are introduced so that the $`2\times 2`$ mass matrix,

$$\left(\begin{array}{cc}0& m_{t\chi }\\ m_{\chi t}& m_{\chi \chi }\end{array}\right)$$
(1.6)
is of seesaw form and has a light eigenvalue corresponding to the observed top quark. The value of $`m_{t\chi }`$ is related to the weak scale, and its value is estimated to be 600 GeV .

The coupling of the top-quark to $`\chi `$ violates custodial symmetry in the same way that the top-quark mass does in the standard model. The leading contribution to $`T`$ from the underlying top seesaw physics arises from contributions to $`W`$ and $`Z`$ vacuum polarization diagrams involving the $`\chi `$. This contribution is positive and is calculated to be

$$\mathrm{\Delta }T=\frac{N_c}{16\pi ^2\alpha _{em}(M_Z^2)}\frac{m_{t\chi }^4}{m_{\chi \chi }^2v^2}\frac{0.7}{\alpha _{em}}\left(\frac{\mathrm{\Lambda }^2}{m_{\chi \chi }^2}\right)\left(\frac{v^2}{\mathrm{\Lambda }^2}\right),$$
(1.7)
which is of the form of eqn. 1.2 with $`b\kappa ^2(\mathrm{\Lambda }/m_{\chi \chi })^2`$. Note that $`\mathrm{\Lambda }/m_{\chi \chi }`$ cannot be small since top-color gauge interactions must drive $`t\chi `$ chiral symmetry breaking. Taking $`\mathrm{\Lambda }/m_{\chi \chi }4`$, we reproduce the positive branch of the boundary of the light region excluded by triviality shown in Figure 1. By varying $`\mathrm{\Lambda }`$ and $`m_{\chi \chi }`$, the entire allowed $`(m_H,\mathrm{\Delta }T)`$ region with positive $`\mathrm{\Delta }T`$ and to the left of the triviality constraint can be obtained. In particular, we note that it is possible to obtain a light Higgs boson in this context as well.

The fact that contributions to $`T`$ greatly expand the region of allowed Higgs mass in the top seesaw model is discussed in detail in ref. . Here we see that the running of the Higgs self-coupling encoded in the constraints of eqn. 1.3 prevent Higgs masses higher than about 500 GeV from being realized .

“Composite Higgs Models”  also provide examples of theories with a potentially heavy composite Higgs boson. In the simplest of these models, one introduces three new fermions which couple to a vectorial “ultracolor” $`SU(N)`$ gauge interaction. Two of these fermions ($`\psi `$) transform as a vectorial doublet under $`SU(2)_W`$, while the third ($`\sigma `$) is assumed to be a singlet. Dirac mass terms can be introduced for all of these fermions and, as so far described, chiral symmetry breaking driven by the ultracolor interactions leaves the vectorial $`SU(2)_W\times U(1)_Y`$ unbroken. Extra chiral interactions are then introduced to misalign the vacuum by a small amount, causing a nonzero $`\overline{\psi }\sigma `$ condensate and breaking the weak interactions.

The octet of pions which result from ultracolor chiral symmetry breaking include a set, the analogs of the kaons, which form a composite Higgs boson. Models can be constructed  in which the Higgs boson can formally be as heavy as a TeV (i.e. at tree-level), while the other four pions have masses controlled by the ultracolor scale and can be much heavier.

This simplest model does not have a custodial symmetry. A direct calculation of the $`W`$ and $`Z`$ masses yields the positive contribution

$$\mathrm{\Delta }T=\frac{v^2}{4\alpha _{em}(M_Z^2)f^2}.$$
(1.8)
Here $`f`$ is the pion decay constant for ultracolor chiral symmetry breaking, the analog of $`f_\pi `$ in QCD. The ultracolor chiral symmetry breaking scale, estimated  to be $`𝒪(4\pi f)`$, sets the compositeness scale $`\mathrm{\Lambda }`$ of the Higgs boson. Comparing eqns. 1.2 and 1.8, we see that the contribution to $`T`$ is of the same form with $`b\kappa ^24\pi ^2`$, excluding the light and dark shaded regions to the right in Figure 1. From this we see that phenomenologically acceptable composite Higgs models can be constructed with Higgs masses up to approximately 450 GeV. Again, in this case by varying the Dirac masses of the fermions and adjusting the size of the chiral interaction, it is possible to construct models that realize any $`(m_H,\mathrm{\Delta }T)`$ to the left of the triviality constraint for positive $`\mathrm{\Delta }T`$.

Finally, we briefly consider the prospects for improving these indirect limits over the next few years. The measurements of $`M_W`$ and $`m_t`$ are likely to be greatly improved during Run II of the Fermilab Tevatron. With an integrated luminosity of 10 fb-1, it may be possible to reduce the uncertainty in the top mass to 2 GeV and in the $`W`$ mass to 30 MeV . To illustrate the potential of these measurements, in Figure 2 we plot the 68% and 95% CL bounds in the $`(m_H,\mathrm{\Delta }T)`$ plane which would be allowed if $`M_W`$ and $`m_t`$ assumed their current “best-fit” values while the uncertainties dropped as projected. Note that although the 95% CL region is somewhat smaller than in Fig. 1 (e.g. the two degree of freedom upper bound on the “standard model” Higgs boson mass – $`\mathrm{\Delta }T=0`$ – drops to $`𝒪`$(180 GeV)), there would still be composite Higgs models consistent with electroweak data with a Higgs boson mass up to 500 GeV for positive $`\mathrm{\Delta }T`$.

In a forthcoming publication , we will detail the calculation of corrections to precisely measured electroweak quantities in the two composite Higgs models we reviewed above and consider the complementary constraints arising from bounds on $`Zb\overline{b}`$, flavor-changing neutral currents , and CP-violation.

Acknowledgments

We thank Gustavo Burdman, Aaron Grant, Marko Popovic, and Elizabeth Simmons for useful discussions. NE is grateful to PPARC for the sponsorship of an Advanced Fellowship. This work was supported in part by the Department of Energy under grant DE-FG02-91ER40676.


text:
# Lithium abundances in main-sequence F stars and sub-giants

## References

Burkhart, C., Coupry, M.F. 2000, preprint

Charbonnel, C., Talon, S. 1999, A&A 351, 635

Charbonnel, C., Talon, S. 2000, this conference

Dias do Nascimento Jr., J.D., Charbonnel, C., Lèbre, A., de Laverny, P., de Medeiros, J.R., 1999, preprint

Maeder, A., & Zahn, J.-P. 1998, A&A, 334, 1000

Mestel L., 1953, M.N.R.A.S. 113, 716

Théado and Vauclair 2000, this conference

Vauclair, S. 1991, in IAU Symp. 145, Evolution of Stars : the Photospheric Abundance Connection ( Michaud, G., Tutukov, A. eds.) p. 327

Vauclair, S. 1999, A&A, 351, 973

Vauclair, S. 2000, this conference

Zahn, J.-P. 1992, A&A, 265, 115


text:
# The diagonalization of quantum field Hamiltonians

## 1 Introduction

Most computational work in non-perturbative quantum field theory and many body phenomena rely on one of two general techniques, Monte Carlo or diagonalization. These methods are nearly opposite in their strengths and weaknesses. Monte Carlo requires relatively little storage, can be performed using parallel processors, and in some cases the computational effort scales reasonably with system size. But it has great difficulty for systems with sign or phase oscillations and provides only indirect information on wavefunctions and excited states. In contrast diagonalization methods do not suffer from fermion sign problems, can handle complex-valued actions, and can extract details of the spectrum and eigenstate wavefunctions. However the main problem with diagonalization is that the required memory and CPU time scales exponentially with the size of the system.

In view of the complementary nature of the two methods, we consider the combination of both diagonalization and Monte Carlo within a computational scheme. We propose a new approach which takes advantage of the strengths of the two computational methods in their respective domains. The first half of the method involves finding and diagonalizing the Hamiltonian restricted to an optimal subspace. This subspace is designed to include the most important basis vectors of the lowest energy eigenstates. Once the most important basis vectors are found and their interactions treated exactly, Monte Carlo is used to sample the contribution of the remaining basis vectors. By this two-step procedure much of the sign problem is negated by treating the interactions of the most important basis states exactly, while storage and CPU problems are resolved by stochastically sampling the collective effect of the remaining states.

In our approach diagonalization is used as the starting point of the Monte Carlo calculation. Therefore the two methods should not only be efficient but work well together. On the diagonalization side there are several existing methods using Tamm-Dancoff truncation , similarity transformations , density matrix renormalization group , or variational algorithms such as stochastic diagonalization . However we find that each of these methods is either not sufficiently general, not able to search an infinite or large dimensional Hilbert space, not efficient at finding important basis vectors, or not compatible with the subsequent Monte Carlo part of the calculation. The Monte Carlo part of our diagonalization/Monte Carlo scheme is discussed separately in a companion paper . In this paper we consider the diagonalization part of the scheme. We introduce a new diagonalization method called quasi-sparse eigenvector (QSE) diagonalization. It is a general algorithm which can operate using any basis, either orthogonal or non-orthogonal, and any sparse Hamiltonian, either real, complex, Hermitian, non-Hermitian, finite-dimensional, or infinite-dimensional. It is able to find the most important basis states of several low energy eigenvectors simultaneously, including those with identical quantum numbers, from a random start with no prior knowledge about the form of the eigenvectors.

Our discussion is organized as follows. We first define the notion of quasi-sparsity in eigenvectors and introduce the quasi-sparse eigenvector method. We discuss when the low energy eigenvectors are likely to be quasi-sparse and make an analogy with Anderson localization. We then consider three examples which test the performance of the algorithm. In the first example we find the lowest energy eigenstates for a random sparse real symmetric matrix. In the second example we find the lowest eigenstates sorted according to the real part of the eigenvalue for a random sparse complex non-Hermitian matrix. In the last example we consider the case of an infinite-dimensional Hamiltonian defined by $`1+1`$ dimensional $`\varphi ^4`$ theory in a periodic box. We conclude with a summary and some comments on the role of quasi-sparse eigenvector diagonalization within the context of the new diagonalization/Monte Carlo approach.

## 2 Quasi-sparse eigenvector method

Let $`|e_i`$ denote a complete set of basis vectors. For a given energy eigenstate

$$|v=\underset{i}{}c_i|e_i,$$
(1)
we define the important basis states of $`|v`$ to be those $`|e_i`$ such that for fixed normalizations of $`|v`$ and the basis states, $`\left|c_i\right|`$ exceeds a prescribed threshold value. If $`|v`$ can be well-approximated by the contribution from only its important basis states we refer to the eigenvector $`|v`$ as quasi-sparse with respect to $`|e_i`$.

Standard sparse matrix algorithms such as the Lanczos or Arnoldi methods allow one to find the extreme eigenvalues and eigenvectors of a sparse matrix efficiently, without having to store or manipulate large non-sparse matrices. However in quantum field theory or many body theory one considers very large or infinite dimensional spaces where even storing the components of a general vector is impossible. For these more difficult problems the strategy is to approximate the low energy eigenvectors of the large space by diagonalizing smaller subspaces. If one has sufficient intuition about the low energy eigenstates it may be possible to find a useful truncation of the full vector space to an appropriate smaller subspace. In most cases, however, not enough is known a priori about the low energy eigenvectors. The dilemma is that to find the low energy eigenstates one must truncate the vector space, but in order to truncate the space something must be known about the low energy states.

Our solution to this puzzle is to find the low energy eigenstates and the appropriate subspace truncation at the same time by a recursive process. We call the method quasi-sparse eigenvector (QSE) diagonalization, and we describe the steps of the algorithm as follows. The starting point is any complete basis for which the Hamiltonian matrix $`H_{ij}`$ is sparse. The basis vectors may be non-orthogonal and/or the Hamiltonian matrix may be non-Hermitian. The following steps are now iterated:

1. Select a subset of basis vectors $`\{e_{i_1},\mathrm{},e_{i_n}\}`$ and call the corresponding subspace $`S`$.
2. Diagonalize $`H`$ restricted to $`S`$ and find one eigenvector $`v`$.
3. Sort the basis components of $`v`$ according to their magnitude and remove the least important basis vectors.
4. Replace the discarded basis vectors by new basis vectors. These are selected at random according to some weighting function from a pool of candidate basis vectors which are connected to the old basis vectors through non-vanishing matrix elements of $`H`$.
5. Redefine $`S`$ as the subspace spanned by the updated set of basis vectors and repeat steps 2 through 5.

If the subset of basis vectors is sufficiently large, the exact low energy eigenvectors will be stable fixed points of the QSE update process. We can show this as follows. Let $`|i`$ be the eigenvectors of the submatrix of $`H`$ restricted to the subspace $`S`$, where $`S`$ is the span of the subset of basis vectors after step 3 of the QSE algorithm. Let $`|A_j`$ be the remaining basis vectors in the full space not contained in $`S`$. We can represent $`H`$ as

$$\left[\begin{array}{cccccc}\lambda _1& 0& \mathrm{}& 1\left|H\right|A_1& 1\left|H\right|A_2& \mathrm{}\\ 0& \lambda _2& \mathrm{}& 2\left|H\right|A_1& 2\left|H\right|A_2& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ A_1\left|H\right|1& A_1\left|H\right|2& \mathrm{}& E\lambda _{A_1}& A_1\left|H\right|A_2& \mathrm{}\\ A_2\left|H\right|1& A_2\left|H\right|2& \mathrm{}& A_2\left|H\right|A_1& E\lambda _{A_2}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right].$$
(2)
We have used Dirac’s bra-ket notation to represent the terms of the Hamiltonian matrix. In cases where the basis is non-orthogonal and/or the Hamiltonian is non-Hermitian, the meaning of this notation may not be clear. When writing $`A_1\left|H\right|1`$, for example, we mean the result of the dual vector to $`|A_1`$ acting upon the vector $`H|1`$. In (2) we have written the diagonal terms for the basis vectors $`|A_j`$ with an explicit factor $`E`$. We let $`|1`$ be the approximate eigenvector of interest and have shifted the diagonal entries so that $`\lambda _1=0.`$ Our starting hypothesis is that $`|1`$ is close to some exact eigenvector of $`H`$ which we denote as $`|1_{\text{full}}`$. More precisely we assume that the components of $`|1_{\text{full}}`$ outside $`S`$ are small enough so that we can expand in inverse powers of the introduced parameter $`E.`$

We now expand the eigenvector as

$$|1_{\text{full}}=\left[\begin{array}{c}1\\ c_2^{}E^1+\mathrm{}\\ \mathrm{}\\ c_{A_1}^{}E^1+\mathrm{}\\ c_{A_2}^{}E^1+\mathrm{}\\ \mathrm{}\end{array}\right]$$
(3)
and the corresponding eigenvalue as

$$\lambda _{\text{full}}=\lambda _1^{}E^1+\mathrm{}.$$
(4)
In (3) we have chosen the normalization of $`|1_{\text{full}}`$ such that $`1|1_{\text{full}}=1`$. From the eigenvalue equation

$$H|1_{\text{full}}=\lambda _{\text{full}}|1_{\text{full}}$$
(5)
we find at lowest order

$$c_{A_j}^{}=\frac{A_j\left|H\right|1}{\lambda _{A_j}}.$$
(6)
We see that at lowest order the component of $`|1_{\text{full}}`$ in the $`|A_j`$ direction is independent of the other vectors $`|A_j^{}`$. If $`|1`$ is sufficiently close to $`|1_{\text{full}}`$ then the limitation that only a fixed number of new basis vectors is added in step 4 of the QSE algorithm is not relevant. At lowest order in $`E^1`$ the comparison of basis components in step 3 (in the next iteration) is the same as if we had included all remaining vectors $`|A_j`$ at once. Therefore at each update only the truly largest components are kept and the algorithm converges to some optimal approximation of $`|1_{\text{full}}`$. This is consistent with the actual performance of the algorithm as we will see in some examples later. In those examples we also demonstrate that the QSE algorithm is able to find several low energy eigenvectors simultaneously. The only change is that when diagonalizing the subspace $`S`$ we find more than one eigenvector and apply steps 3 and 4 of the algorithm to each of the eigenvectors.

## 3 Quasi-sparsity and Anderson localization

As the name indicates the accuracy of the quasi-sparse eigenvector method depends on the quasi-sparsity of the low energy eigenstates in the chosen basis. If the eigenvectors are quasi-sparse then the QSE method provides an efficient way to find the important basis vectors. In the context of our diagonalization/Monte Carlo approach, this means that diagonalization does most of the work and only a small amount of correction is needed. This correction is found by Monte Carlo sampling the remaining basis vectors, a technique called stochastic error correction . If however the eigenvectors are not quasi-sparse then one must rely more heavily on the Monte Carlo portion of the calculation.

The fastest and most reliable way we know to determine whether the low energy eigenstates of a Hamiltonian are quasi-sparse with respect to a chosen basis is to use the QSE algorithm and look at the results of the successive iterations. But it is also useful to consider the question more intuitively, and so we consider the following example.

Let $`H`$ be a sparse Hermitian $`2000\times 2000`$ matrix defined by

$$H_{jk}=\mathrm{log}(j)\delta _{jk}+x_{jk}M_{jk},$$
(7)
where $`j`$ and $`k`$ run from $`1`$ to $`2000`$, $`x_{jk}`$ is a Gaussian random real variable centered at zero with standard deviation $`x_{\text{rms}}=0.25`$, and $`M_{jk}`$ is a sparse symmetric matrix consisting of random $`0`$’s and $`1`$’s such that the density of $`1`$’s is $`5\%`$. The reason for introducing the $`\mathrm{log}(j)`$ term in the diagonal is to produce a large variation in the density of states. With this choice the density of states increases exponentially with energy. Our test matrix is small enough that all eigenvectors can be found without difficulty. We will consider the distribution of basis components for the eigenvectors of $`H`$. In Figure 1 we show the square of the basis components for a given low energy eigenvector $`|v.`$ The basis components are sorted in order of descending importance. The ratio of $`\mathrm{\Delta }E`$, the average spacing between neighboring energy levels, to $`x_{\text{rms}}`$ is $`0.13`$. We see that the eigenvector is dominated by a few of its most important basis components. In Figure 2 we show the same plot for another eigenstate but one where the spacing between levels is three times smaller, $`\mathrm{\Delta }E/x_{\text{rms}}=0.041.`$ This eigenvector is not nearly as quasi-sparse. The effect is even stronger in Figure 3, where we show an eigenvector such that the spacing between levels is $`\mathrm{\Delta }E/x_{\text{rms}}=0.024`$.

Our observations show a strong effect of the density of states on the quasi-sparsity of the eigenvectors. States with a smaller spacing between neighboring levels tend to have basis components that extend throughout the entire space, while states with a larger spacing tend to be quasi-sparse. The relationship between extended versus localized eigenstates and the density of states has been studied in the context of Anderson localization and metal-insulator transitions . The simplest example is the tight-binding model for a single electron on a one-dimensional lattice with $`Z`$ sites,

$$H=\underset{j}{}d_j|jj\left|+\underset{jj^{}}{}t_{jj^{}}\right|jj^{}|.$$
(8)
$`|j`$ denotes the atomic orbital state at site $`j,`$ $`d_j`$ is the on-site potential, and $`t_{jj^{}}`$ is the hopping term between nearest neighbor sites $`j`$ and $`j^{}`$. If both terms are uniform ($`d_j=d,`$ $`t_{jj^{}}=t`$) then the eigenvalues and eigenvectors of $`H`$ are

$`Hv_n`$ $`=(d+2t\mathrm{cos}\frac{2\pi n}{Z})v_n,`$ (9)
$`v_n`$ $`=\frac{1}{\sqrt{Z}}{\displaystyle \underset{j}{}}e^{i{\scriptscriptstyle \frac{2\pi nj}{Z}}}|j,`$ (10)
where $`n=1,\mathrm{},Z`$ labels the eigenvectors. In the absence of diagonal and off-diagonal disorder, the eigenstates of $`H`$ extend throughout the entire lattice. The eigenvalues are also approximately degenerate, all lying within an interval of size 4$`t`$. However, if diagonal and/or off-diagonal disorder is introduced, the eigenvalue spectrum becomes less degenerate. If the disorder is sufficiently large, the eigenstates become localized to only a few neighboring lattice sites giving rise to a transition of the material from metal to insulator.

We can regard a sparse quantum Hamiltonian as a similar type of system, one with both diagonal and general off-diagonal disorder. If the disorder is sufficient such that the eigenvalues become non-degenerate, then the eigenvectors will be quasi-sparse. We reiterate that the most reliable way to determine if the low energy states are quasi-sparse is to use the QSE algorithm. Intuitively, though, we expect the eigenstates to be quasi-sparse with respect to a chosen basis if the spacing between energy levels is not too small compared with the size of the off-diagonal entries of the Hamiltonian matrix.

## 4 Finite matrix examples

As a first test of the QSE method, we will find the lowest four energy states of the random symmetric matrix $`H`$ defined in (7). So that there is no misunderstanding, we should repeat that diagonalizing a $`2000\times 2000`$ matrix is not difficult. The purpose of this test is to analyze the performance of the method in a controlled environment. One interesting twist is that the algorithm uses only small pieces of the matrix and operates under the assumption that the space may be infinite dimensional. A sample MATLAB program similar to the one used here has been printed out as a tutorial example in .

The program starts from a random configuration, 70 basis states for each of the four eigenvectors. With each iteration we select $`10`$ replacement basis states for each of the eigenvectors. In Figure 4 we show the exact energies and the results of the QSE method as functions of iteration number. In Figure 5 we show the inner products of the normalized QSE eigenvectors with the normalized exact eigenvectors. We note that all of the eigenvectors were found after about 15 iterations and remained stable throughout successive iterations. Errors are at the $`5`$ to $`10\%`$ level, which is about the theoretical limit one can achieve using this number of basis states. The QSE method has little difficulty finding several low lying eigenvectors simultaneously because it uses the distribution of basis components for each of the eigenvectors to determine the update process. This provides a performance advantage over variational-based techniques such as stochastic diagonalization in finding eigenstates other than the ground state.

As a second test we consider a sparse non-Hermitian matrix with complex eigenvalues. This type of matrix is not amenable to variational-based methods. We will find the four eigenstates corresponding with eigenvalues with the lowest real part for the random complex non-Hermitian matrix

$$H_{jk}^{}=(1+ic_{jk})H_{jk}.$$
(11)
$`H_{jk}`$ is the same matrix used previously and$`c_{jk}`$ is a uniform random variable distributed between $`1`$ and 1. As before the program is started from a random configuration, 70 basis states for each of the four eigenvectors. For each iteration $`10`$ replacement basis vectors are selected for each of the eigenvectors. In Figure 6 the exact eigenvalues and the results of the QSE run are shown in the complex plane as functions of iteration number. In Figure 7 we show the inner products of the QSE eigenvectors with the exact eigenvectors. All of the eigenvectors were found after about 20 iterations and remained stable throughout successive iterations. Errors were again at about the $`5`$ to $`10\%`$ level.

## 5 $`\varphi ^4`$ theory in $`1+1`$ dimensions

We now apply the QSE method to an infinite dimensional quantum Hamiltonian. We consider $`\varphi ^4`$ theory in $`1+1`$ dimensions, a system that is familiar to us from previous studies using Monte Carlo  and explicit diagonalization . The Hamiltonian density for $`\varphi ^4`$ theory in $`1+1`$ dimensions has the form

$$=\frac{1}{2}\left(\frac{\varphi }{t}\right)^2+\frac{1}{2}\left(\frac{\varphi }{x}\right)^2+\frac{\mu ^2}{2}\varphi ^2+\frac{\lambda }{4!}\text{:}\varphi ^4\text{:},$$

where the normal ordering is with respect to the mass $`\mu `$. We consider the system in a periodic box of length $`2L`$. We then expand in momentum modes and reinterpret the problem as an equivalent Schrödinger equation . The resulting Hamiltonian is

$`H`$ $`=\frac{1}{2}{\displaystyle \underset{n}{}}\frac{}{q_n}\frac{}{q_n}+\frac{1}{2}{\displaystyle \underset{n}{}}\left(\omega _n^2(\mu )\frac{\lambda b(\mu )}{8L}\right)q_nq_n`$ (12)
$`+\frac{\lambda }{4!2L}{\displaystyle \underset{n_1+n_2+n_3+n_4=0}{}}q_{n_1}q_{n_2}q_{n_3}q_{n_4}`$
where

$$\omega _n(\mu )=\sqrt{\frac{n^2\pi ^2}{L^2}+\mu ^2}$$
(13)
and $`b(\mu )`$ is the coefficient for the mass counterterm

$$b(\mu )=\underset{n}{}\frac{1}{2\omega _n(\mu )}.$$
(14)
It is convenient to split the Hamiltonian into free and interacting parts with respect to an arbitrary mass $`\mu ^{}`$:

$$H_{free}=\frac{1}{2}\underset{n}{}\frac{}{q_n}\frac{}{q_n}+\frac{1}{2}\underset{n}{}\omega _n^2(\mu ^{})q_nq_n,$$
(15)
$`H`$ $`=H_{free}+\frac{1}{2}{\displaystyle \underset{n}{}}\left(\mu ^2\mu ^2\frac{\lambda b(\mu )}{8L}\right)q_nq_n`$ (16)
$`+\frac{\lambda }{4!2L}{\displaystyle \underset{n_1+n_2+n_3+n_4=0}{}}q_{n_1}q_{n_2}q_{n_3}q_{n_4}.`$
$`\mu ^{}`$ is used to define the basis states of our Fock space. Since $`H`$ is independent of $`\mu ^{}`$, we perform calculations for different $`\mu ^{}`$ to obtain a reasonable estimate of the error. It is also useful to find the range of values for $`\mu ^{}`$ which maximizes the quasi-sparsity of the eigenvectors and therefore improves the accuracy of the calculation. For the calculations presented here, we set the length of the box to size $`L=5\pi \mu ^1`$. We restrict our attention to momentum modes $`q_n`$ such that $`\left|n\right|N_{\mathrm{max}}`$, where $`N_{\mathrm{max}}=20`$. This corresponds with a momentum cutoff scale of $`\mathrm{\Lambda }=4\mu .`$

To implement the QSE algorithm on this infinite dimensional Hilbert space, we first define ladder operators with respect to $`\mu ^{}`$,

$`a_n(\mu ^{})`$ $`=\frac{1}{\sqrt{2\omega _n(\mu ^{})}}\left[q_n\omega _n(\mu ^{})+\frac{}{q_n}\right]`$ (17)
$`a_n^{}(\mu ^{})`$ $`=\frac{1}{\sqrt{2\omega _n(\mu ^{})}}\left[q_n\omega _n(\mu ^{})\frac{}{q_n}\right].`$ (18)
The Hamiltonian can now be rewritten as

$`H`$ $`={\displaystyle \underset{n}{}}\omega _n(\mu ^{})a_n^{}a_n+\frac{1}{4}(\mu ^2\mu ^2\frac{\lambda b}{8L}){\displaystyle \underset{n}{}}\frac{\left(a_n+a_n^{}\right)\left(a_n+a_n^{}\right)}{\omega _n(\mu ^{})}`$ (19)
$`+\frac{\lambda }{192L}{\displaystyle \underset{n_1+n_2+n_3+n_4=0}{}}\left[\frac{\left(a_{n_1}+a_{n_1}^{}\right)}{\sqrt{\omega _{n_1}(\mu ^{})}}\frac{\left(a_{n_2}+a_{n_2}^{}\right)}{\sqrt{\omega _{n_2}(\mu ^{})}}\frac{\left(a_{n_3}+a_{n_3}^{}\right)}{\sqrt{\omega _{n_3}(\mu ^{})}}\frac{\left(a_{n_4}+a_{n_4}^{}\right)}{\sqrt{\omega _{n_4}(\mu ^{})}}\right].`$
In (19) we have omitted constants contributing only to the vacuum energy. We represent any momentum-space Fock state as a string of occupation numbers, $`|o_{N_{\mathrm{max}}},\mathrm{},o_{N_{\mathrm{max}}}`$, where

$$a_n^{}a_n|o_{N_{\mathrm{max}}},\mathrm{},o_{N_{\mathrm{max}}}=o_n|o_{N_{\mathrm{max}}},\mathrm{},o_{N_{\mathrm{max}}}.$$
(20)
From the usual ladder operator relations, it is straightforward to calculate the matrix element of $`H`$ between two arbitrary Fock states.

Aside from calculating matrix elements, the only other fundamental operation needed for the QSE algorithm is the generation of new basis vectors. The new states should be connected to some old basis vector through non-vanishing matrix elements of $`H`$. Let us refer to the old basis vector as $`|e`$. For this example there are two types of terms in our interaction Hamiltonian, a quartic interaction

$$\underset{n_1,n_2,n_3}{}\left(a_{n_1}+a_{n_1}^{}\right)\left(a_{n_2}+a_{n_2}^{}\right)\left(a_{n_3}+a_{n_3}^{}\right)\left(a_{n_1n_2n_3}+a_{n_1+n_2+n_3}^{}\right),$$
(21)
and a quadratic interaction

$$\underset{n}{}\left(a_n+a_n^{}\right)\left(a_n+a_n^{}\right).$$
(22)
To produce a new vector from $`|e`$ we simply choose one of the possible operator monomials

$`a_{n_1}a_{n_2}a_{n_3}a_{n_1n_2n_3},a_{n_1}^{}a_{n_2}a_{n_3}a_{n_1n_2n_3},\mathrm{},`$ (23)
$`a_na_n,a_n^{}a_n^{},\mathrm{}`$
and act on $`|e`$. Our experience is that the interactions involving the small momentum modes are generally more important than those for the large momentum modes, a signal that the ultraviolet divergences have been properly renormalized. For this reason it is best to arrange the selection probabilities such that the smaller values of $`\left|n_1\right|`$, $`\left|n_2\right|`$, $`\left|n_3\right|`$ and $`\left|n\right|`$ are chosen more often.

For each QSE iteration, $`50`$ new basis vectors were selected for each eigenstate and $`250`$ basis vectors were retained. The results for the lowest energy eigenvalues are shown in Figure 8. The error bars were estimated by repeating the calculation for different values of the auxiliary mass parameter $`\mu ^{}`$.

From prior Monte Carlo calculations we know that the theory has a phase transition at $`\frac{\lambda }{4!}2.5\mu ^2`$ corresponding with spontaneous breaking of the $`\varphi \varphi `$ reflection symmetry. In the broken phase there are two degenerate ground states and we refer to these as the even and odd vacuum states. In Figure 8 we see signs of a second order phase transition near $`\frac{\lambda }{4!}2.5\mu ^2`$. Since we are working in a finite volume the spectrum is discrete, and we can track the energy eigenvalues as functions of the coupling. Crossing the phase boundary, we see that the vacuum in the symmetric phase becomes the even vacuum in the broken phase while the one-particle state in the symmetric phase becomes the odd vacuum. The energy difference between the states is also in agreement with a Monte Carlo calculation of the same quantities. The state marking the two-particle threshold in the symmetric phase becomes the one-particle state above the odd vacuum, while the state at the three-particle threshold becomes the one-particle state above the even vacuum. These one-particle states should be degenerate in the infinite volume limit. One rather unusual feature is the behavior of the first two-particle state above threshold in the symmetric phase. In the symmetric phase this state lies close to the two-particle threshold. But as we cross the phase boundary the state which was the two-particle threshold is changed into a one-particle state. Thus our two-particle state is pushed up even further to become a two-particle state above the even vacuum and we see a pronounced level crossing.

We note that while the one-particle mass vanishes near the critical point, the energies of the two-particle and three-particle thresholds reach a minimum but do not come as close to zero energy. It is known that this model is repulsive in the two-particle scattering channel. In a large but finite volume the ground state and one-particle states do not feel significant finite volume effects. The two-particle state at threshold, however, requires that the two asymptotic particles be widely separated. In our periodic box of length 2$`L`$ the maximal separation distance is $`L`$ and we expect an increase in energy with respect to twice the one-particle mass of size $`V(L)`$, where $`V`$ is the potential energy between particles. Likewise a three-particle state will increase in energy an amount $`3V(2L/3)`$. Our results indicate that finite volume effects for the excited states are significant for this value of $`L`$.

## 6 Summary

We have proposed a new approach which combines both diagonalization and Monte Carlo within a computational scheme. The motivation for our approach is to take advantage of the strengths of the two computational methods in their respective domains. We remedy sign and phase oscillation problems by handling the interactions of the most important basis states exactly using diagonalization, and we deal with storage and CPU problems by stochastically sampling the contribution of the remaining states. We discussed the diagonalization part of the method in this paper. The goal of diagonalization within our scheme is to find the most important basis vectors of the low energy eigenstates and treat the interactions among them exactly. We have introduced a new diagonalization method called quasi-sparse eigenvector diagonalization which achieves this goal efficiently and can operate using any basis, either orthogonal or non-orthogonal, and any sparse Hamiltonian, either real, complex, Hermitian, non-Hermitian, finite-dimensional, or infinite-dimensional. Quasi-sparse eigenvector diagonalization is the only method we know which can address all of these problems.

We considered three examples which tested the performance of the algorithm. We found the lowest energy eigenstates for a random sparse real symmetric matrix, the lowest eigenstates (sorted according to the real part of the eigenvalue) for a random sparse complex non-Hermitian matrix, and the lowest energy eigenstates for an infinite-dimensional Hamiltonian defined by $`1+1`$ dimensional $`\varphi ^4`$ theory in a periodic box.

We regard QSE diagonalization as only a starting point for the Monte Carlo part of the calculation. Once the most important basis vectors are found and their interactions treated exactly, a technique called stochastic error correction is used to sample the contribution of the remaining basis vectors. This method is introduced in .

### Acknowledgments

We thank P. van Baal, H. S. Seung, H. Sompolinsky, and M. Windoloski for useful discussions. Support provided by the National Science Foundation.


text:
# Interaction-assisted propagation of Coulomb-correlated electron-hole pairs in disordered semiconductors

## Abstract

A two-band model of a disordered semiconductor is used to analyze dynamical interaction induced weakening of localization in a system that is accessible to experimental verification. The results show a dependence on the sign of the two-particle interaction and on the optical excitation energy of the Coulomb-correlated electron-hole pair.

The problem of two interacting particles (TIP) in a random potential is an excellent paradigm for the general question of the interplay of disorder and interactions in many-body systems. First addressed in a 1990 paper by Dorokhov , the subject has been especially well studied since Shepelyansky’s 1994 publication . Considering the TIP localization length $`l_2`$, most authors  obtain an interaction-induced increase $`l_2>l_1`$ over the single-particle localization length $`l_1`$ independent of the sign of the interaction, with $`l_2/l_1l_1^a`$ and $`a=1`$ or $`a=0.65`$ . Here, $`l_1`$ and $`l_2`$ are measured in units of the lattice constant of a one-dimensional Anderson chain. Similar results have been obtained for TIP in a quasiperiodic chain . The independence of the predicted effect on the sign of the interaction is an especially intriguing feature.

Early work approached this problem using a wide variety of theoretical techniques and focused on establishing the existence of the TIP-effect, while more recent work has dealt with quantitative details like scaling behaviour and the influence of interaction range, strength and sign .

Existing works comprise purely theoretical case studies since the model of just two particles in a single band does not correspond to any real physical situation. Obviously, experimental study is needed to promote further understanding of the TIP problem, and would put the presently rather academic discussion on a firm physical basis. Exploiting the fact that the coherent spatio-temporal dynamics of Coulomb-correlated electron-hole pairs is strongly influenced by the two-particle interaction, we show in the present paper that the TIP localization properties should be accessible to modern ultrafast optical techniques. The corresponding spatial and temporal scales are on the order of sub-$`\mu `$m and 1 ps. Our numerical studies are based on integration of the Semiconductor Bloch Equations and include no a priori assumptions about energy hierarchies or interaction matrix elements. The calculations were performed for a disordered 1D semiconductor quantum wire, where localization effects are most important. Despite its simplicity, this model system contains already all essential ingredients to describe the dynamics following optical excitation even in disordered systems of higher dimension . The model parameters have been given values that resemble those of realistic disordered semiconductor quantum wires. Additional physical parameters (excitation energy, spectral pulse width, screening length, different masses of the two particles) allow the study of a wide variety of observable phenomena.

In our numerical calculations we investigate the spreading of an electron-hole wave packet after local excitation by an optical pulse. Here the interaction is given by the long-range Coulomb potential which, besides producing bound states (excitons) near the edges of the excitation spectrum, also correlates the electrons and holes in the pair continuum. Previous theoretical studies of the spatial-temporal dynamics of wave packets formed from excitons show that their motion is rather limited in the presence of scattering . This result is recovered by our present calculations. Here we focus our interest on the dynamics of optically generated wave packets in the pair continuum. We find that the excitation conditions in the presence of particle-particle interaction influence the carrier dynamics dramatically. In addition, and in contrast to some previous claims in the literature, we find that the sign of the interaction has a pronounced effect on the spatio-temporal dynamics.

We consider a 1D array of sites $`i`$ with diagonal disorder in both the valence band ($`vb`$) and the conduction band ($`cb`$). The site energies $`\epsilon _{vi}`$ and $`\epsilon _{ci}`$ corresponding to the $`vb`$ and the $`cb`$, respectively, are randomly distributed over the interval $`[W/2,W/2]`$ and are uncorrelated. The nearest neighbor $`cb`$ levels are coupled by the tunneling term $`J^c`$, the $`vb`$ levels by $`J^v`$. We use the Coulomb interaction in its monopole-monopole form  with matrix elements

$$V_{ij}=\frac{U}{4\pi ϵϵ_0}\frac{e^2}{r_{ij}+\alpha }$$
(1)
which has been regularized in order to cope with the pathological singularity in one dimension. The constant $`\alpha `$ has been chosen to be 5 times the lattice constant . The total Hamiltonian is written as $`\widehat{H}=\widehat{H_0}+\widehat{H_I}+\widehat{H_C}`$, where $`\widehat{H_0}`$ describes the $`vb`$ and $`cb`$ band structure, $`\widehat{H_I}`$ represents the semiclassical interaction with the electromagnetic field $`E_i(t)`$ in dipole approximation, and $`\widehat{H_C}`$ defines the electron-electron interaction term .

In the following, we assume a local initial excitation at the central site $`i=0`$, which is modeled by setting $`\mu _i=\delta _{i,0}`$ for the local dipole matrix element in $`\widehat{H_I}`$. The optical polarization $`p_{ij}(t)`$ is obtained from the equation of motion for the polarization operator $`\widehat{p}_{ij}=\widehat{d}_i\widehat{c}_j`$, which is coupled to the equation of motion of the electron and hole intraband quantities $`\widehat{n}_{ij}^e=\widehat{c}_i^+\widehat{c}_j`$ and $`\widehat{n}_{ij}^h=\widehat{d}_i^+\widehat{d}_j`$, respectively, where the operators $`\widehat{c}_i^+,\widehat{c}_i`$ ($`\widehat{d}_i^+,\widehat{d}_i`$) describe the electron (hole) creation and annihilation operators at site $`i`$. The equation of motion for the expectation values $`p_{ij}(t)`$ and $`n_{ij}^{e,h}(t)`$ is treated using the well-known Semiconductor Bloch Equations for $`p_{ij}(t)`$ and $`n_{ij}^{e,h}(t)`$ in the real-space representation . Detailed derivations of the Semiconductor Bloch Equations are given in  and the textbook . As we are interested in small excitation densities, we write only the equation for $`p_{ij}`$ in the lowest (linear-response) order in the exciting field

$`_tp_{ij}=`$ $``$ $`i\left(\epsilon _i^e\epsilon _j^hV_{ij}\right)p_{ij}+i{\displaystyle \underset{l=1}{\overset{N}{}}}\left(J^ep_{il}+J^hp_{lj}\right)`$ (2)
$`+`$ $`i\mu _jE_j(t)\delta _{ij}.`$ (3)
Using the conservation laws $`n_{ij}^e=_lp_{lj}p_{li}^{}`$ and $`n_{ij}^h=_lp_{jl}p_{il}^{}`$ valid in this lowest order , we obtain the intraband quantities.

Instead of studying a rather academic localization length which describes only the asymptotic behavior of wave functions, we calculate the experimentally more relevant participation number $`\mathrm{\Lambda }(t)=(_in_{ii}^2)^1`$. Here $`n_{ii}`$ stands for either $`n_{ii}^e`$ or $`n_{ii}^h`$. With the packet localized at site $`0`$, $`n_{ii}=\delta _{i0}`$ and $`\mathrm{\Lambda }=1`$, while for an excitation uniformly extended over the sample of $`N`$ sites, $`n_{ii}=1/N`$ and $`\mathrm{\Lambda }=N`$. Our calculations were performed for chains containing $`N=240`$ sites. Boundary effects can easily be identified in the temporal evolution of $`\mathrm{\Lambda }`$ and do not play any role as long as $`\mathrm{\Lambda }<N/2`$. All the data presented are free of finite-size effects.

The transform-limited optical pulse is defined by its mean energy $`\mathrm{}\omega `$ and the temporal width $`\tau `$ of the gaussian envelope $`exp\{(t/\tau )^2\}`$. We define an excitation energy $`E_{exc}`$ refered to the bottom of the (ordered) absorption band, i.e. $`E_{exc}=\mathrm{}\omega E_{gap}`$. All results are given for $`\tau =100`$ fs, which corresponds to an energetic width (FWHM) of 22 meV.

To make contact with the previous work where two particles in a single band were placed initially at a single site, we first consider the situation of a symmetric band structure with $`J^e=J^h=20`$ meV. The absorption spectra with and without Coulomb interaction are shown in Fig. 1 for the ordered case. The peak structure near the absorption is due to the excitonic resonances. Upon changing the sign of the Coulomb interaction, the bound state resonances are shifted from the bottom to the top of the absorption spectrum.

As the dynamics of electrons and holes are the same for the assumed symmetric band structure, we restrict our discussion to the electrons. We first discuss the situation in the absence of Coulomb interaction. Fig. 2 shows the corresponding $`\mathrm{\Lambda }_e(t)`$ for two different disorder parameters $`W`$ after excitation by a pulse at $`E_{exc}=80`$ meV. The excitation is centered in the absorption spectrum as indicated in Fig. 1. $`\mathrm{\Lambda }_e(t)`$ evolves exponentially with rise time less than 1 ps. Here and below, we take the saturation value as a measure of localization. As expected, it decreases rapidly with increasing disorder. We find $`\mathrm{\Lambda }W^{1.3}`$ as $`W`$ is varied over the range 40 meV to 240 meV for $`J=20`$ meV. A discussion of related exponents can be found in .

Fig. 2 contrasts the interacting and noninteracting behavior for two values of disorder and reveals three remarkable features. i) The interaction clearly leads to a reduction of the localization of the particles. We have carefully checked that the saturation value of $`\mathrm{\Lambda }_e(t)`$ at long times is not due to a finite size effect; values $`N/2`$ are fully converged with respect to the sample size. ii) While the participation number in the noninteracting situation evolves exponentially and saturates quickly ($`<1`$ ps), the interacting wave packets evolve diffusively and reach their saturation values at much longer times. iii) The sign of the Coulomb interaction ($`U=\pm 1`$) has virtually no influence on the propagation of the particles in this case. The same is true if we apply a very short excitation pulse which spectrally covers the whole band. The spectral position of the central pulse frequency within the band is then completely irrelevant. In this situation the excited particle pair-wave packet is initially situated exclusively at site $`i=0`$. These observations are not new . However, iii) has been questioned  and ii) remained unexplained.

In all cases where $`J^e`$ and $`J^h`$ are of comparable magnitude we find that the participation number is enhanced by the interaction. In a mean field picture, it is the temporal fluctuations of the field originating from the partner particle which destroy the coherence necessary to produce localization. This explanation in terms of a dynamic-correlation-induced weakening of the influence of disorder can be nicely corroborated by a number of case studies.

We note that contrary to previous statements, the independence of the sign of the Coulomb interaction is not a general feature, but is a consequence of the imposed electron-hole symmetry. In particular, displacing the central frequency of excitation pulses from the center of the absorption band, the situation changes completely. Note that this choice of the excitation frequency corresponds to the realistic situation where electron-hole pairs are excited close to the absorption edge in semiconductors. Fig. 3 shows the participation number $`\mathrm{\Lambda }`$ for light electrons and heavy holes, i.e., $`J^e=2J^h=20`$ meV. The central excitation energy of the pulse is placed in the lower part of the pair continuum at $`E_{exc}=40`$ meV. Results averaged over 60 realizations are shown for $`W=80`$ meV and $`U=0,\pm 1`$. The results are invariant under reflection of the excitation frequency through band center and simultaneous switching of the sign of the interaction. This reflects the approximate symmetry (within fluctuations in the site energy distribution) of the Hamiltonian.

It is at first sight counterintuitive that the enhancement of the participation number is larger for attractive ($`U=1`$) than for repulsive ($`U=+1`$) interaction. This behavior can be attributed to the fact that for attractive interaction and positive masses (i.e. for excitation into the lower half of the excitation continuum) the electron-hole pair tends to stay closer together. The fluctuating field due to the accompanying particle is then more pronounced as compared to the case of repulsive interaction, where the mutually repulsive particle pair tends to be separated. Hence the dynamic-correlation-induced weakening of the influence of disorder is less effective for repulsive than for attractive interaction.

Completely different behavior is found for a static field. We consider an infinitely heavy hole, $`J^h=0`$, which now produces a static field, and excitation at the (interaction-free) band center. For both attractive and repulsive interactions, the participation number is decreased with respect to the noninteracting case. This result is easily understood without invoking fluctuating fields since at band center electron states have maximal extent. In the presence of interaction, off-center states are admixed leading to greater confinement. The effect of the static interaction is thus opposite to that of a fluctuating field .

The strong retardation of the saturation in the interacting case can also be understood in our picture. Whether with or without interaction, the electron and hole wave packets spread over a range given by the single-particle levels involved in the optical transition just after the short excitation pulse. The fluctuating Coulomb field due to the partner particle then leads to an increase of the spread of the wave packets. As a consequence, the average fluctuating field acting on a given particle is reduced, which in turn tends to slow further spreading, eventually leading to the observed saturation at long times. The neglect of phonon interactions in our model is justified a posteriori. Fig. 2 makes it clear that the time scales between 100 fs and $``$ 3 ps are fully sufficient for experimental observations while near-band-edge acoustic phonon scattering occurs on longer time scales .

Previous work  on the TIP problem suggests a scaling of the two-particle localization length $`l_2/l_1l_1^a(U/J)^2`$, with $`a=1`$ or $`a=0.65`$ . Our results for the participation number do not obey such a scaling law as far as the dependence on $`U`$ is concerned. We obtain for electrons and holes, both for attractive and repulsive interaction, $`\mathrm{\Lambda }(U=\pm 1)/\mathrm{\Lambda }(U=0)\mathrm{\Lambda }(U=0)^b`$ with $`b=0.65\pm 0.3`$. For electrons and attractive interaction the present model predicts $`\mathrm{\Lambda }_eW^c`$ with a larger exponent $`c2.2`$ compared to $`c=1.3`$ for the noninteracting case.

In conclusion, we have studied the localization of a pair of interacting particles in a situation which, in principle, is accessible to experiments. Optical excitation in the pair continuum of a disordered one-dimensional semiconductor with long-range Coulomb interaction has been considered. Starting from a tight-binding description, the temporal evolution of the participation numbers of the electron and the hole wave packets has been calculated by a direct solution of the equation of motion of the correlated material excitation within linear response with regard to the exciting laser field. The participation number increases with interaction for both attractive and repulsive interaction. We find that in general the degree of delocalization depends strongly on the sign of the interaction, in contrast to previously published predictions. The sign of the interaction becomes irrelevant (even if the masses of electrons and holes are different) only for two special situations: excitation in the center of the pair continuum, or excitation of the whole band. We have checked that this result is independent of the assumed form of the interaction and that it remains true also for the short-range interactions studied in the literature. Compared to the single-band models treated in the past, the present semiconductor model admits a richer variety of phenomena, which can be qualitatively explained within a mean-field picture. We emphasize that the enhancement of the participation number is clearly not due to a finite size effect, and that it should be experimentally observable. Ultra-short time-of-flight experiments on arrays of semiconductor quantum wires in the coherent limit using pump-probe techniques are a promising option. The enhancement should also be observable in disordered semiconductor quantum wells. In this case we expect the enhancement to be even more pronounced, since, in contrast with one-dimensional systems, only states close to the band edge are essentially affected by the disorder in two dimensions, so that the interaction will lead to coupling with rather extended states.

This work is supported by DFG, SFB 383 and 341, the Leibniz Prize, OTKA (T021228, T024136, F024135), SNSF (2000-52183.97), and the A. v. Humboldt Foundation. Discussions with A. Knorr and F. Gebhard are gratefully acknowledged.