text:
# Evolution of level density step structures from 56,57Fe to 96,97Mo

## Abstract

Level densities have been extracted from primary $`\gamma `$ spectra for 56,57Fe and 96,97Mo nuclei using (3He,$`\alpha \gamma `$) and (3He,3He$`{}_{}{}^{}\gamma `$) reactions on 57Fe and 97Mo targets. The level density curves reveal step structures above the pairing gap due to the breaking of nucleon Cooper pairs. The location of the step structures in energy and their shapes arise from the interplay between single-particle energies and seniority-conserving and seniority-non-conserving interactions.

PACS number(s): 21.10.Ma, 25.55.Hp, 27.40.+z, 27.60.+j

The group at the Oslo Cyclotron Laboratory (OCL) has recently developed a new method (the so-called Oslo method) to extract the level density and radiative strength function from primary $`\gamma `$ spectra . The method can be characterized as a further development of the sequential extraction method . The Oslo method has been extensively tested in the rare-earth mass region  and has been successfully extended to the very light 27,28Si nuclei . The robustness of the Oslo method in the mid sd shell where the level density is rather low and the single-particle levels are widely spaced prompted the present investigations in mass regions around the $`f_{7/2}`$ shell (the iron nuclei) and around the fpg shell (the molybdenum nuclei). The 56Fe nucleus was especially chosen due to astrophysical interest in this isotope. While Ref.  discusses the $`e`$-process and the direct production of 56Fe, a modern view in terms of a hierarchy of restricted nuclear statistical equilibria (NSE) on different time scales predicts the direct production of 56Fe only on the highest level of NSE, namely in the presence of weak equilibrium . These rather extreme conditions are envisioned only in the core of massive stars (where ejection into space is unlikely) and in the core of Type Ia supernovae (where ejection into space seems possible). In the equations of NSE, the nuclear partition function, i.e., the Laplace transform of the nuclear level density, plays a prominent role. The nucleus 96Mo is of special interest in the investigation of the $`|NZ|`$ dependence of level densities, since the mass 96 isobars are the only ones in the nuclear chart where one can find three different stable nuclei with $`|NZ|`$ varying by eight units from 96Zr to 96Ru. Also the $`|NZ|`$ dependence of level densities can have significant astrophysical importance , since many reactions take place on unstable isotopes and therefore most of their reaction cross sections must be estimated by Hauser-Feshbach type calculations .

The purpose of this Rapid Communication is to report on experimental level densities of 56,57Fe and 96,97Mo nuclei and to provide a schematic explanation of the observed step structures in these curves.

The experiments were performed at the MC-35 cyclotron at the OCL using a $`2`$ nA beam of 45 MeV 3He particles. The self-supporting targets were isotopically enriched to 94.7% and 94.2% and had thicknesses of 3.4 mg/cm2 and 2.1 mg/cm2 for 57Fe and 97Mo, respectively. Each experiment ran for approximately five days each, and about 200,000 relevant particle-$`\gamma `$ coincidences were recorded in each analyzed reaction channel. The charged ejectiles were identified and their energies measured in a ring of eight collimated Si $`\mathrm{\Delta }E`$-$`E`$ telescopes placed at 45 with respect to the beam direction. The thicknesses of the front and end detectors were 140 and 3000 $`\mu `$m, respectively and shielding against $`\delta `$ electrons was achieved with a 19-$`\mu `$m-thick Al foil. The distance from the target was 5 cm, giving a total solid angle coverage of 0.3% of $`4\pi `$, and the energy resolution was $``$0.3 MeV over the entire spectrum. The $`\gamma `$-rays were detected in 28 collimated 5”x5” NaI(Tl) detectors called CACTUS  surrounding the target and particle detectors. The total efficiency of CACTUS was $``$15% of $`4\pi `$ and the resolution was $``$6% of the deposited energy at 1.3 MeV. In addition, one 60% Ge(HP) detector was used in the setup to monitor the selectivity and populated spin distribution of the reactions. Raw $`\alpha `$-particle data are shown in Fig. 1 for the case of the 97Mo(3He,$`\alpha `$)96Mo reaction, where discrete transfer peaks are observed up to $``$6 MeV in excitation energy indicating the opening of the neutron $`g_{9/2}`$ shell.

From the known $`Q`$-value and the reaction kinematics, the ejectile energy can be transformed into initial excitation energy of the residual nuclei. Using the particle-$`\gamma `$ coincidence technique, each $`\gamma `$ ray can be assigned to a cascade depopulating a certain initial excitation energy in the residual nucleus. The data are therefore sorted into total $`\gamma `$-ray spectra originating from different initial excitation energy bins. Every spectrum is then unfolded using a Compton-subtraction method which preserves the fluctuations in the original spectra and does not introduce further, spurious fluctuations . From the unfolded spectra, a primary $`\gamma `$ matrix is constructed using the subtraction method of Ref. . The basic assumption behind this method is that the $`\gamma `$-ray decay pattern from any excitation energy bin is independent of whether states in this bin are populated directly via the (3He,$`\alpha `$) or (3He,3He) reactions or indirectly via $`\gamma `$ decay from higher excited levels following the initial nuclear reaction. This assumption is trivially fulfilled if one populates the same levels with the same weights within any excitation energy bin, since the decay branchings are properties of the levels and do not depend on the population mechanisms. On the other hand, if one populates, e.g., vastly different spin distributions within one excitation energy bin via the two different population mechanisms, the $`\gamma `$-decay patterns should be different as well. As an example of the data analysis discussed in this paragraph, the raw, unfolded and primary $`\gamma `$ spectra are shown for the 57Fe(3He,$`\alpha \gamma `$)56Fe and 57Fe(3He,3He$`{}_{}{}^{}\gamma `$)57Fe reactions in Fig. 2.

Finally, the primary $`\gamma `$ matrix is factorized using the generalized Brink-Axel hypothesis . The original hypothesis states that the giant dipole resonance (GDR) can be built on every excited state, and that the properties of the GDR do not depend on the temperature of the nuclear state on which it is built. This hypothesis can be generalized to include not only the GDR but any type of nuclear excitation and results in the assumption that primary $`\gamma `$ spectra originating from the excitation energy $`E`$ can be factorized into a $`\gamma `$-ray transmission coefficient $`𝒯(E_\gamma )`$ which depends only on the $`\gamma `$-transition energy $`E_\gamma `$ and into the level density $`\rho (EE_\gamma )`$ at the final energy. This factorization is determined by a least $`\chi ^2`$ fit to the primary $`\gamma `$ matrix, using no a priori assumptions about the functional form of either the level density or the $`\gamma `$-ray transmission coefficient . An example to illustrate the quality of the fit is shown in Fig. 3, where for the 97Mo(3He,3He$`{}_{}{}^{}\gamma `$)97Mo reaction we compare the experimental primary $`\gamma `$ spectra from two different initial excitation energies to the least $`\chi ^2`$ fit. Unfortunately, the mathematical structure of the relevant equations in the least $`\chi ^2`$ fit does not allow us to find a unique solution for the level density and $`\gamma `$-ray transmission coefficient. However, it has been shown that all solutions with the same $`\chi ^2`$ can be obtained by the transformation of one randomly chosen solution according to

$`\stackrel{~}{\rho }(EE_\gamma )`$ $`=`$ $`A\mathrm{exp}[\alpha (EE_\gamma )]\rho (EE_\gamma )`$ (1)
$`\stackrel{~}{𝒯}(E_\gamma )`$ $`=`$ $`B\mathrm{exp}(\alpha E_\gamma )𝒯(E_\gamma ).`$ (2)
The three free parameters $`A`$, $`B`$, and $`\alpha `$ have to be determined to give the physically most relevant solution to the least $`\chi ^2`$ fit using independent experimental information. The most common way is to count the number of discrete levels at low excitation energies and to use the neutron resonance spacing at $`B_n`$ to find values for $`A`$ and $`\alpha `$. The remaining parameter $`B`$ is then determined using the average total radiative width of neutron resonances . Unfortunately, in the case of 56Fe, there are no data on neutron resonances and thus, the information about the level density around $`B_n`$ in 56Fe has to be obtained by different means. In order to do so, we calculate the level density at $`B_n`$ in 57Fe using a backshifted Fermi-gas expression with the parameterization of von Egidy et al. , where we apply an additional overall re-normalization factor to match the level density determined from neutron resonance spacings. Then, we use the same level-density formula (including the same re-normalization factor but with the other parameters for 56Fe) to calculate the level density at $`B_n`$ in 56Fe. Using this data point instead of the unknown neutron resonance spacing, we proceed in the same way as for the other three nuclei .

In Fig. 4 we show the physically most relevant solutions for the level densities in 56,57Fe and 96,97Mo from the least $`\chi ^2`$ fit. The most striking feature in those curves are the steps starting at 2.9 MeV and 1.8 MeV in 56Fe and 57Fe. There are also less pronounced but still statistically significant step structures at 2.0 MeV and 1.2 MeV in 96Mo and 97Mo. It has been established in the rare-earth region that such low-lying step structures are connected to the breaking of the first nucleon Cooper pair . It is therefore natural to assume that the same is true for the lighter mass regions. The additional step at 6 MeV in the level density curve of 96Mo might be due to the opening of the $`g_{9/2}`$ shell which is indicated to happen at this excitation energy (see Fig. 1). A future systematic investigation of Mo nuclei might shed more light on this issue.

The difference in binding energy between the even and odd systems is a measure for pairing correlation and can be calculated from the three-mass indicator of Dobaczewski et al. . This indicator yields 1.2 MeV and 0.9 MeV for the Fe and Mo nuclei, respectively, which agrees very well with the differences in excitation energy for the first steps in the level density curves. Thus we know that the steps in neighboring nuclei appear at the same ’effective’ excitation energy, i.e., the excitation energy corrected for the contribution to the binding energy due to neutron-pairing correlations. Further, the average proton pairing energies are 0.8 MeV and 1.1 MeV for the 57Fe and 97Mo nuclei, respectively. These energies should now correspond to the excitation energies of the steps in the two odd nuclei. However, the steps are delayed in excitation energy by 1.0 MeV and 0.1 MeV for these two nuclei. This might be explained by the fact that one not only has to invest the energy to break a proton Cooper pair but also at least one of the unpaired protons has to be promoted to the next unoccupied single-particle level. The spacing of those levels can also be calculated using the Dobaczewski three-mass indicator, giving 2.0 MeV and 0.9 MeV in the two cases. Clearly, the higher single particle spacing for the 56Fe nucleus leads to a larger delay in excitation energy for the step structure to appear compared to the 97Mo nucleus, but obviously the exact excitation energy for the steps in the level density curves cannot be estimated from binding energies, since it will depend on the exact location of the Fermi energy within the single-particle level scheme.

Another complication might arise from the effect of seniority non-conserving interactions, which mix configurations of different seniority and smooth out the step structures in the level densities. To investigate this, we have performed a model calculation. In the model we assume a system of eight particles scattered into an equidistant single-particle level scheme with eight doubly-degenerate levels. As residual interactions we consider a pairing interaction and a seniority non-conserving interaction. Thus, the model Hamiltonian is written as

$`\widehat{H}=ϵ{\displaystyle \underset{i=1}{\overset{8}{}}}ia_i^{}a_i{\displaystyle \frac{1}{2}}G{\displaystyle \underset{i,j=1}{\overset{8}{}}}a_i^{}a_{\overline{ı}}^{}a_{\overline{ȷ}}a_j`$ (4)
$`{\displaystyle \frac{1}{2}}\kappa {\displaystyle \underset{i,j,k,l=1}{\overset{8}{}}}W_{ijkl}a_i^{}a_j^{}a_ka_l,`$
where $`a^{}`$ and $`a`$ are Fermion creation and annihilation operators and the labels with bars stand for time reversed orbits. The single particle level spacing $`ϵ`$, the strength $`G`$ of the pairing interaction and the strength $`\kappa `$ of the seniority non-conserving interaction $`W`$ are the only macroscopic parameters of the model. This model with good seniority, i.e., the case of $`\kappa =0`$ has already been diagonalized in Ref. , and the dotted line in Fig. 5 gives the distribution of eigenvalues with excitation energy, i.e., the exact level density of the model. The individual bumps contain mainly levels with the same seniority, thus the step structures in the level density can be explained by the consecutive breaking of nucleon Cooper pairs. However, the experimental data show that in general the steps are much smoother than in such a simple model calculation indicating the need for seniority non-conserving terms in the Hamiltonian.

An important contribution to seniority breaking comes from quadrupole collectivity, which can change rapidly with mass number for the nuclei under study. Strong, residual quadrupole-quadrupole interactions can lead to deformation of the nucleus and to a modification of the single-particle level spacing. Unfortunately, our simple model cannot incorporate a realistic quadrupole-quadrupole interaction; the smoothing of the level density, however, should not depend on the details of the seniority non-conserving residual interaction. We choose therefore to model the quadrupole-quadrupole interaction by a random two-body interaction  of roughly equivalent strength where all the pairing-like terms have been set to zero i.e., the $`W_{ijkl}`$ in Eq. (4) are Gaussian random numbers of mean zero and width equal to one except for the cases of $`j=\overline{ı}`$ and $`k=\overline{l}`$ where $`W_{i\overline{ı}\overline{l}l}=0`$.

Since in this more general case Eq. (4) does not have good seniority, exact diagonalization of the model Hamiltonian has to be performed within the full model space which is computationally very demanding, as one requires all levels, not just the lowest. For eight particles in eight doubly-degenerate states there are 12870 states (counting all magnetic substates). The number of levels (not counting magnetic substates) equals the number of states with spin projection $`J_z=0`$, thus yielding 4900 levels. For the present calculations, we used $`ϵ=0.25`$ MeV, $`G=0.5`$ MeV, and $`\kappa 0`$ MeV (pure pairing) and $`\kappa =0.14`$ MeV (pairing+random interaction). In Fig. 5, we show the resulting level densities of the two calculations. The resolved bumps with definite seniority in the pure pairing case are smeared out by the random interaction. The gaps between the bumps are rapidly being filled, while the bumps themselves are degraded into a step structure quite similar to the step structure seen in the experimental level density curves of the iron nuclei. The occurrence of the step structure in the calculation is actually quite sensitive to the exact choice of the strength of the random interaction. A weak random interaction ($`\kappa 0.1`$ MeV) does not fill the gaps between the bumps, a strong random interaction ($`\kappa 0.2`$ MeV) produces a more smeared out step structure in the level density curve similar to the Mo data, indicating the presence of relatively stronger seniority non-conserving interactions in the Mo nuclei. The range of values which produces good qualitative agreement with the Fe data is about 0.13 MeV$`\kappa 0.16`$ MeV. However, due to the simplicity of the model, no attempt has been made to achieve quantitative agreement between the model and the experimental data.

In conclusion, we have presented new experimental data on level densities below $`B_n`$ in 56,57Fe and 96,97Mo. Step structures in the level density curves have been related to the breaking of nucleon Cooper pairs. The occurrence of the step structures in energy depends on the exact location of the Fermi energy in the single particle level scheme. The smoothness of the step structures depends on the strength of the seniority non-conserving interactions, e.g., the quadrupole-quadrupole interaction compared to the pairing interaction and the average single-particle level spacing. The more pronounced step structures in the level density data of Fe nuclei compared to Mo indicate relatively stronger seniority non-conserving interactions in the latter case. Our schematic calculations where random two-body interactions were used to model residual seniority non-conserving interactions show that level density curves are very sensitive to the relative strengths of the different terms in the model Hamiltonian. The present experimental data therefore break new ground for the possible theoretical investigation of two of the most important residual interactions in atomic nuclei, namely the pairing and the quadrupole-quadrupole interaction

Part of this work was performed under the auspices of the U.S. Department of Energy by the University of California, Lawrence Livermore National Laboratory under Contract No. W-7405-ENG-48. Financial support from the Norwegian Research Council (NFR) is gratefully acknowledged. G.M. acknowledges support by a U.S. Department of Energy grant with No. DE-FG02-97-ER41042.


format: text

text:
# Disorder–Induced Effects in III–V Semiconductors with Mn

## I Introduction

The properties of the diluted magnetic semiconductors (DMS), such as II-VI and III-V materials doped with Mn, are mostly determined by the presence of open d-shells of Mn. The d-electrons form the local magnetic moments, and the hybridization of the d-states with the band states leads to strong magnetooptical and magnetoelectric phenomena Furdyna88 ; Ohno99 . In addition, the magnetic polarization of the band states of the host semiconductor results in kinetic exchange interaction between the moments.

The change of the crystal potential due to the Mn substitution has also an important direct effect on the band states. In the simplest description, the substitution is represented by replacing the atomic s- and p-levels $`E_s`$ and $`E_p`$ of the host atom by the atomic levels of the impurity. In traditional II-VI DMS with Mn, the differences in $`E_s`$ are of order of 1 - 2 eV tables . This, with respect to the typical band structure of the II-VI semiconductors, represents a perturbation of an intermediate strength Larson88 , which can not be treated within the virtual crystal approximation (VCA).

In the III-V semiconductors, however, the atomic levels of Mn differ from the atomic levels of the group-III cations even more. The change of $`E_s`$ in GaAs and GaP doped with Mn is approximately 4 eV tables  and also the difference of the p-levels is remarkable. Strong impurity potential, comparable to the bandwidth, is expected to cause a strong reconstruction of the electronic structure of the mixed crystal. This prediction of the tight-binding theory has been recently confirmed also by our recent calculations Masek01 ; Maca02 .

We expect a particularly strong effect of the Mn substitution in the lowest conduction bands, because it is mostly composed of the cationic s–states. This means that the impurity scattering should have a pronounced implication on the near-edge absorption due to the violation of the momentum conservation and the appearance of the non-direct transitions Spicer67  from the valence to the conduction band.

The situation is, however, more complex. It is well established that Mn doped III-V DMS reveal an almost complete compensation Beschoten99 . This means that the concentration of either intentional donors or native compensating defects increases proportionally to the concentration of Mn Maca02 . These donors represent an additional source of the disorder scattering.

We present here the first ab-initio study of the interband absorption in the III-V DMS. It is combined with tight-binding CPA calculations suitable for the interpretation of the disorder-induced features in the electronic structure of the dilute mixed crystals.

We concentrate on the interband absorption. That is why we consider only the limiting case of a complete compensation, where the IR transition inside the valence band are absent.

## II Absorption coefficient

Although the main interest is paid to (Ga,Mn)As, we applied our density functional study to (Ga,Mn)P with a wider bandgap that is less sensitive to narrowing of the band gap due to the local density approximation. We take SeP as the compensating donor. The mixed crystal is represented by a periodic system with a large unit cell (LUC) containing 16 atoms. One molecular unit (Ga–P) of the LUC is replaced by a close Mn–Se pair. The electronic structure and the imaginary part of the dielectric function, $`\epsilon _2(\omega )`$, were calculated by the full-potential linearized augmented-plane-waves (FPLAPW) method WIEN97 . The spin–orbit coupling was not included in our spin–polarized calculations. The generalized gradient approximation to the exchange–correlation potential was used. The results are summarized in Figs. 1 and 2.

Fig. 1 shows the spin-polarized total density of states (DOS) of the Ga7MnP7Se crystal in the low–temperature, ferromagnetic state altogether with the DOS of pure GaP. The bands of Ga7MnP7Se are strongly magnetically polarized due to high concentration (12.5 at. %) of magnetic impurities. The spin splitting of the valence band (1.2 eV) corresponds to the exchange constant $`J_{pd}2`$ eV. The bandgaps for both majority and minority spin electrons, 0.57 eV and 1.29 eV, respectively, are smaller than the bandgap obtained for the pure GaP (1.67 eV). In contrast to (Ga,Mn)As Maca02 , the present material with a small overlap (0.16 eV) of both bandgaps is a semiconductor. The width of the common bandgap is expected to increase for lower Mn doping.

Although the uppermost part of the majority spin DOS at the top of the valence band can be viewed as an impurity band formed by the acceptor states, no gap states are created due to Se impurities. They act as donors with zero activation energy. The presence of Se has no important effect, neither on the distribution of the Mn d-states (see Fig. 1) nor on the total DOS around the bandgap. It means that most of the changes in the absorption should be attributed to Mn. The Mn s-states contribute mostly to the broad maximum of the conduction band DOS around 4.5 eV and they are almost completely absent at the edge of the conduction band, so that its bottom is locally eroded around Mn. Also the empty Mn d-states substantially modify the nature of the bottom of the conduction band for the minority-spin electrons.

In Fig. 2 we compare the calculated $`\epsilon _2(\omega )`$ for Ga7MnP7Se with the results for the GaP crystal from near-edge to the UV range. The absorption edge of pure GaP is followed by a steep increase of the absorption at 2.85 eV (E1 transition) and by a main sharp peak at 4.4 eV (E2 transition).

The change due to the Mn (and Se) doping is remarkable. The onset of the absorption is shifted to lower energies by $`1.1`$ eV. Instead of well developed absorption edge, a smooth increase of the absorption reminds an Urbach tail Urbach53 . The interband absorption does not show any structure typical for the zinc-blende semiconductors, only a broadened E2 peak survives. The E1 transition, more sensitive to the arrangement of the bands along the $`\mathrm{\Lambda }`$ line in the Brillouin zone, disappears completely and the absorption gradually increases between 1.5 eV and 4 eV eV.

## III Spectral density

To describe the DMS with lower impurity concentrations we use the coherent potential description (CPA) Velicky68  of the electronic structure of the DMS. The coherent potential method, resulting in configurationally averaged quantities, can also easily handle the room-temperature, magnetically disordered phase of the DMS. The avaraging over all random distributions of the impurities restores the full translational symmetry of the lattice, so that the wavevector $`k`$ can be used as a quantum number. The spectral density $`A(k,E)`$, i.e. the DOS decomposed into contributions from the various points of the Brillouin zone, represents a detailed information about the dispersion of the electron states.

We apply a tight-binding version of the CPA Masek87  to GaAs co-doped with 4 at. % of Mn and 4 at. % of Se. The parametrization of the tight-binding Hamiltonian provides a correct band gap for the pure GaAs crystal Talwar82 , as well as appropriate exchange splitting of the Mn d–states. As shown in Fig. 3 for minority spin electrons, the influence of the dopants on the spectral density is strong. Similar results were obtained for the majority spin electrons and also for the paramagnetic phase of Ga0.92Mn0.04Se0.04.

The linewidth of the Bloch states in the highest valence band and in the lowest conduction band is, in all cases, approximately 0.5 eV. The linewidth of the corresponding excitations, in particular of the E1 transitition is $`1.0`$ eV, large enough to smear out the E1 peak in $`\epsilon _2(\omega )`$. However, it is important to note that the effect of the impurity broadening depends on the orbital composition of the bands and on the position in the Brillouin zone, as seen in Fig. 3. The narrow peaks of $`A(k,E)`$ indicate those parts of the bands where the effect of the impurities is rather weak.

## IV Comparison with other donors

The almost complete compensation Beschoten99  in the III-V DMS with Mn means that the total concentration of donors, either intentional or native defects, is almost equal to the concentration of Mn, i.e. a few atomic percent. In such case, we have mixed rather than doped semiconductors and the modification of their band structure due to the donors must be considered. We compare the electronic structure of GaAs with 4 at. % of Mn and with a corresponding concentration of the compensating donors. We consider donors in both anion and cation sublattices (Se, Sn), and As antisite defects.

The most favorite donors are the As antisite defects. They are double donors and create an impurity band deep in the bandgap, centered at 0.8 eV above the valence band. Also the co-doping with tin results in an impurity band around 1.7 eV, which is, however, resonant with the bottom of the conduction band. Both impurities give an inhomogeneous broadening ($`0.7`$ eV) of the conduction–band states.

Finally, the intentional compensation with Se, considered in Sections 2 a 3, can be viewed as an alloying of the III-V materials with MnSe, without formation of the impurity states.

## V Summary

The substitution of Mn in the GaAs crystal is connected with a scattering of electrons on the dopants. It results in broadening of the spectral lines and in smearing out the electron spectra. This is particularly seen in the calculated absorption coefficient. We found that the absorption edge of the compensated (Ga,Mn)P co–doped with Se is shifted by $``$ 1.0 eV with respect to pure GaP crystal. The absorption is smooth, without any structure typical to the E1 transitions in the III-V materials.

Three aspects of the electron scattering on impurities can be distiguished: (i) Hybridization of the host states with the Mn d–states, acting locally on the p–states of As atoms bonded to Mn. This effect is particularly important in the upper part of the valence band. (ii) Impurity potential acting on the cationic s- and- p–states that modifies mostly the conduction band states. (iii) Scattering on compensating donors that are present in real III-V materials in a concentration comparable to the concentration of Mn.

While (i) and (ii) are inherent to the Mn impurity, (iii) may be optimized by the choice of the compensating donor. We found that the effect of Se is weak in comparison to the influence of Mn. The modification of the band states by both As antisite defects and substitutional Sn is much stronger. These two impurities also form impurity bands in (Ga,Mn)As while Se does not.

Assuming that the concentration of Mn determines the total number of the donors we deduce that the presence of the intentional donors can reduce the number of native antisite defects. If it is so, the co-doping of the III-V DMS with Se–like impurities will be quite important.

Acknowledgment The financial support was provided by the Academy of Sciences of the Czech Republic (Grant No. A1010214) and by RTN project ”Computational Magnetoelectronics” of the European Commission (HPRN-CT-2000-00143).


text:
# Universal crossing probabilities and incipient spanning clusters in directed percolation

## 1 Introduction

Most of the properties of the critical state at a second order phase transition are related to the existence of self-similar correlated clusters. In the percolation problem, the critical point is characterized by the spontaneous occurence of infinite clusters of connected sites, crossing the system in the thermodynamic limit.

The study of crossing probabilities for standard percolation on finite systems at the percolation threshold has been an active field of research during the last decade \[1–44\] (for a recent review see cardy01 ).

In two dimensions, the crossing probability $`\pi `$ can be defined as the probability to find at least one cluster joining two opposite edges of a rectangular-shaped system with length $`L_{}`$ and width $`L_{}`$. At the percolation threshold, in the finite-size scaling (FSS) limit, i.e., when $`L_{}\mathrm{}`$, $`L_{}\mathrm{}`$, while the aspect ratio $`r=L_{}/L_{}`$ remains fixed, the crossing probability is a scale-invariant universal function $`\pi (r)`$ of the aspect ratio langlands92 ; cardy92 .

Following the extensive numerical study of Langlands et al. langlands92 , an analytical expression for $`\pi (r)`$ at criticality was derived by Cardy cardy92 . Using the relation between percolation and the cluster formulation of the $`q`$-state Potts model in the limit $`q1`$ kasteleyn69 ; kasteleyn72  as well as the techniques of boundary conformal field theory, he was able to calculate the probability for two non-overlapping segments on the edge of the half-plane to be connected. The crossing probability in the rectangular geometry could then be deduced from the appropriate conformal mapping. One may notice that the scale invariance of $`\pi (r)`$ is a non-trivial property, resulting from the occurence of a vanishing scaling dimension $`x(q)`$ associated with a boundary condition changing operator of the Potts model in the limit $`q1`$ for percolation.

Other exact crossing formulas were later obtained for different geometries pinson94 ; watts96 . Since then, mathematically rigorous proofs have been provided for some of these results smirnov01 ; schramm01 .

Crossing probabilities on same-spin Ising clusters langlands00 ; lapalme01 ; arguin02a  and Potts clusters arguin02b  in two dimensions have been also considered.

Recently, some attention has been paid to the properties of incipient spanning clusters \[31–44\]. The probability to find $`n`$ incipient spanning clusters, connecting two disjoint segments at the boundary of a finite critical system, has been studied through Monte Carlo simulations \[31–39\]. Rigourous bounds on the spanning probability aizenman97 ; aizenman98  have been obtained as well as analytical expressions through conformal and Coulomb-gas methods cardy98 ; aizenman99 ; cardy02 .

In two recent works turban02 ; turban03 , we started a study of the critical crossing probability (in space-time) for off-equilibrium systems in 1+1 dimensions. In these strongly anisotropic systems, when the temperature-like scaling field $`\mathrm{\Delta }`$ is non-vanishing, the correlation functions generally display an exponential decay, in the space and time directions, characterized by a correlation length $`\xi `$ and a relaxation time $`\tau `$ diverging respectively as binder89

$$\xi =\widehat{\xi }\mathrm{\Delta }^\nu ,\tau =\widehat{\tau }\mathrm{\Delta }^{z\nu },$$
(1.1)
when $`\mathrm{\Delta }0`$. The prefactors $`\widehat{\xi }`$ and $`\widehat{\tau }`$ are non-universal amplitudes and $`z`$ is the dynamical (anisotropy) exponent.

The directed percolation problem broadbent57 ; kinzel83 ; hinrichsen00  was studied numerically through Monte Carlo simulations in turban02  and analytical results were obtained for an exactly solvable diffusion-limited coagulation process in turban03 . In these sytems, the crossing probability in the time direction $`\pi _t`$ is defined as the probability that the system of size $`L`$ remains active at time $`t`$. Anisotropic scaling binder89 ; hucht02  then implies that the appropriate aspect ratio is $`r=t/L^z`$, i.e., the rescaled time. We found that in the FSS limit, as for isotropic critical systems, the critical crossing probability is a scale-invariant universal function of an effective aspect ratio which is the product of $`r`$ by a non-universal constant $`c`$ r .

In the present work, we continue the examination of the universal properties of the crossing probability $`\pi _t`$ in directed percolation which, in a dynamical interpretation, essentially corresponds to the survival probability of an epidemic spreading without immunization harris74 . Section 2 is devoted to a detailed study of the influence of the initial state on the universal behaviour of $`\pi _t`$. In Section 3 we present some results concerning the scaling behaviour of the probability to find $`n`$ incipient spanning clusters, i.e., $`n`$ disjoint centres of infection surviving in the FSS limit.

## 2 Critical crossing probability

We study the critical crossing probability $`P_t`$, in the time direction, for bond and site directed percolation in $`1+1`$ dimensions through Monte Carlo simulations on the square lattice. The time axis is oriented in the diagonal direction. The space and time coordinates take alternatively integer and half-integer values on successive spatial rows. On a site located at $`(x,t)`$, with $`1xL`$, two directed bonds are leaving, which terminate on the nearest-neighbour sites at $`(x\pm 1/2,t+1/2)`$ as shown in Figure 1. We use either free boundary conditions (FBC) or periodic boundary conditions (PBC) in the spatial direction. In the later case, $`x`$ is defined modulo $`L`$.

In the epidemic language, occupied (empty) sites correspond to infected (healthy) individuals. In the bond problem, the bonds are independently open with probability $`p`$. A site is occupied at $`t+1/2`$ when it is connected via an open bond to a site which was occupied at time $`t`$. In the site problem all the bonds are open and a site is occupied with probability $`p`$ at time $`t+1/2`$ when at least one of its first neighbours at time $`t`$ is occupied.

A directed percolation cluster is a collection of connected occupied sites starting from some source at $`t=0`$. We consider successively the three types of initial states shown in Figure 1. A sample contributes to the crossing probability when at least one cluster survives at time $`t`$. The simulations are performed at the percolation threshold, $`p_\mathrm{c}^{\mathrm{bond}}=0.644700185(5)`$ for the bond problem and $`p_\mathrm{c}^{\mathrm{site}}=0.70548522(4)`$ for the site problem. These values, as well as the dynamical exponent $`z=1.580745(10)`$, are taken from jensen99 . The crossing probability $`P_t(L,t,l)`$ is determined as a function of the aspect ratio, $`r=t/L^z`$, and the fraction $`f=l/L`$ of occupied sites in the initial state, using $`10^6`$ samples with sizes up to $`L=2^{10}`$. When $`f`$ is non-vanishing, $`P_t`$ converges to a scale-invariant function, $`\pi _t(r,f)`$, in the FSS limit where $`L`$, $`t`$ and $`l\mathrm{}`$, while keeping $`r`$ and $`f`$ fixed.

### 2.1 Initial state with a fully occupied lattice

When the $`L`$ sites are randomly occupied with probability $`p_i>0`$, the crossing probability converges to a scale-invariant function in the FSS limit. This is illustrated in Figure 2 in the case of site percolation with FBC. When $`p_i=1`$ the convergence is from above and the results obtained with $`L=8`$ are already quite close to the FSS limit. Thus $`p_i`$ is an irrelevant variable when $`f>0`$. In the following we always take $`p_i=1`$ in order to obtain a better convergence towards the scale-invariant behaviour.

We now consider site and bond directed percolation with either FBC or PBC in order to study the universality of the crossing probability. The lattice is fully occupied in the initial state. The behaviour of the crossing probability as a fonction of the aspect ratio is shown in Figure 3. The raw data converge quickly to scale-invariant functions as the system size is increased. They are shown in the inset for the largest size, $`L=256`$. Different functions are obtained for the two types of boundary conditions and the curves for the bond and site problems are shifted with respect to one another, along the horizontal axis, by an amount which is the same for FBC and PBC. The main figure shows the data collapse on universal curves, one for each type of boundary condition, after the results for site percolation have been shifted.

In order to explain this behaviour let us consider the scaling behaviour of the crossing probability for an off-critical system with a deviation $`\mathrm{\Delta }=|pp_\mathrm{c}|`$ from the percolation threshold. Under an anisotropic change of the length scale by a factor $`b`$ the length transforms as $`L^{}=L/b`$, the time as $`t^{}=t/b^z`$ and $`\mathrm{\Delta }`$ as $`\mathrm{\Delta }^{}=b^{1/\nu }\mathrm{\Delta }`$, where $`\nu `$ is the correlation length exponent defined in (1.1). Assuming that, as indicated by the Monte Carlo data, the crossing probability is scale invariant, then it satisfies the scaling relation

$$P_t(L,t,\mathrm{\Delta })=P_t(\frac{L}{b},\frac{t}{b^z},b^{1/\nu }\mathrm{\Delta }).$$
(2.1)
With $`b=\mathrm{\Delta }^\nu `$ one obtains

$$P_t(L,t,\mathrm{\Delta })=P_t(\frac{L}{\mathrm{\Delta }^\nu },\frac{t}{\mathrm{\Delta }^{z\nu }},1)=g(\frac{L}{\xi },\frac{t}{\tau }),$$
(2.2)
where $`g(u,v)`$ is a universal fonction of its dimensionless arguments. The FSS limit corresponds to $`\mathrm{\Delta }=0`$ and $`b=L`$ in (2.1), which gives

$$P_t(L,t,0)=P_t(1,\frac{t}{L^z},0)=\pi _t(cr).$$
(2.3)
For each type of boundary conditions, the crossing probability is a scale-invariant universal function of the effective aspect ratio $`r_{\mathrm{eff}}=cr`$. The non-universal amplitude $`c`$ depends on the choice of the length and time units. It can be expressed as a function of $`\widehat{\xi }`$ and $`\widehat{\tau }`$ hucht02 . Comparing (2.2) and (2.3), $`r`$ appears through the dimensionless ratio $`v/u^z`$ so that

$$r_{\mathrm{eff}}=\frac{t/\tau }{(L/\xi )^z}=\frac{\widehat{\xi }^z}{\widehat{\tau }}r,c=\frac{\widehat{\xi }^z}{\widehat{\tau }}.$$
(2.4)
The non-universal constant $`c`$ is different for the two percolation problems. According to (2.3), for a given type of boundary condition, the crossing probabilities $`\pi _t`$ are identical for the site and bond problems when the values of the aspect ratio, $`r_\mathrm{s}`$ and $`r_\mathrm{b}`$, satisfy the relation $`c_\mathrm{s}r_\mathrm{s}=c_\mathrm{b}r_\mathrm{b}`$. Thus the data collapse in Figure 3 is obtained through a shift of the site percolation data by

$$\delta _{\mathrm{sb}}^r=\mathrm{ln}r_\mathrm{b}\mathrm{ln}r_\mathrm{s}=\mathrm{ln}\left(\frac{c_\mathrm{s}}{c_\mathrm{b}}\right)=.2476(8).$$
(2.5)
The numerical value of $`\delta _{\mathrm{sb}}^r`$ was estimated in turban02  through a least-square fit.

The crossing probability can be calculated using a transfer matrix in the time direction $`𝖳`$, working with a restricted basis of states corresponding to configurations with at least one occupied site derrida80 . The matrix element $`T_{\alpha \beta }`$ gives the probability that the configuration in state $`|\alpha `$ at time $`t`$ leads to the configuration in state $`|\beta `$ at time $`t+1`$. Given an initial state $`|i`$, the crossing probability then reads

$$\pi _t(cr)=\underset{\alpha }{}i|𝖳^t|\alpha .$$
(2.6)
When $`r=t/L^z1`$, the leading contribution comes from the largest eigenvalue $`\mathrm{\Lambda }_{\mathrm{max}}<1`$ of $`𝖳`$ and

$$\pi _t(cr)\mathrm{\Lambda }_{\mathrm{max}}^t=\mathrm{exp}\left(a\frac{ct}{L^z}\right)=\mathrm{exp}(acr)$$
(2.7)
since the relaxation time $`|\mathrm{ln}(\mathrm{\Lambda }_{\mathrm{max}})|^1`$ scales as $`L^z`$. In the argument of the exponential, $`a`$ is some universal constant. Thus $`\mathrm{ln}\pi _t`$ is a universal linear function of $`cr`$ when $`r1`$. This universal linear variation is shown in Figure 4 where $`r`$ for the site problem has been multiplied by $`c_\mathrm{s}/c_\mathrm{b}=1.281(1)`$ in order to eliminate the influence of the non-universal constant $`c`$ on the data collapse. Since the argument leading to (2.7) is quite general, the same behaviour is expected for the other types of initial states considered below.

### 2.2 Initial state with a sequence of occupied sites

Let us consider now the behaviour of the crossing probability when $`l`$ consecutive sites are occupied in the initial state, as shown in Figure 1b. The results for PBC and FBC are presented in Figures 5 and 6, respectively, for three values of the occupation ratio $`f=l/L`$, from top to bottom $`f=1/2`$, $`1/4`$, $`1/16`$. The raw data for the largest size, $`L=256`$, are shown in the insets. As above, different curves are obtained for the site and bond problems. Shifting the site results by $`\delta _{\mathrm{sb}}^r`$, the data collapse on three universal curves, one for each value of $`f`$, as shown on the main figures.

The problem involves the new length scale, $`l`$, and the crossing probability scales as

$$P_t(L,t,l,\mathrm{\Delta })=P_t(\frac{L}{b},\frac{t}{b^z},\frac{l}{b},b^{1/\nu }\mathrm{\Delta }).$$
(2.8)
With $`b=\mathrm{\Delta }^\nu `$ one obtains

$$P_t(L,t,l,\mathrm{\Delta })=P_t(\frac{L}{\mathrm{\Delta }^\nu },\frac{t}{\mathrm{\Delta }^{z\nu }},\frac{l}{\mathrm{\Delta }^\nu },1)=h(\frac{L}{\xi },\frac{t}{\tau },\frac{l}{\xi }),$$
(2.9)
where $`h(u,v,w)`$ is a universal function of its arguments. In the FSS limit, at the critical point, equation (2.8) gives

$$P_t(L,t,l,0)=P_t(1,\frac{t}{L^z},\frac{l}{L},0)=\pi _t(cr,f).$$
(2.10)
The dimensionless occupation ratio $`f`$ does not introduce any new non-universal constant in the scaling limit where $`L`$ and $`l\mathrm{}`$. Thus, as shown in Figures 5 and 6, the crossing probability is a scale-invariant universal function of the effective aspect ratio $`r_{\mathrm{eff}}`$, which depends on the value of $`f`$ as well as on the boundary conditions. One may notice a slightly slower convergence to the limiting behaviour when $`f`$ decreases.

### 2.3 Initial state with a single occupied site

When a single site (located in the middle of the system for FBC) is occupied in the initial state, the crossing probability $`\pi _t^{(1)}`$ is no longer scale invariant. It acquires a dimension which is the scaling dimension of the order parameter,

$$x_m=\frac{\beta }{\nu }=.252072(8),$$
(2.11)
corresponding to the values

$$\beta =.276486(8),\nu =1.096854(4),$$
(2.12)
taken from jensen99 . This behaviour follows from the scaling of $`P_t^{(1)}`$, the probability that a cluster grown from a single central occupied site survives after $`t`$ time steps,

$$P_t^{(1)}(L,t,\mathrm{\Delta })=b^{x_m}P_t^{(1)}(\frac{L}{b},\frac{t}{b^z},b^{1/\nu }\mathrm{\Delta }).$$
(2.13)
On the infinite system with $`b=t^{1/z}`$, one obtains

$$P_t^{(1)}(\mathrm{},t,\mathrm{\Delta })=t^{x_m/z}P_t^{(1)}(\mathrm{},1,t^{1/z\nu }\mathrm{\Delta }),$$
(2.14)
in agreement with equation (94) in reference hinrichsen00 , where the exponent $`\beta ^{}/\nu _{}=\beta /\nu _{}=x_m/z`$ for directed percolation with our notations. Taking $`b=L`$ at criticality ($`\mathrm{\Delta }=0`$) leads in the FSS limit to

$$P_t^{(1)}(L,t,0)=L^{x_m}P_t^{(1)}(1,\frac{t}{L^z},0)=\pi _t^{(1)}(cr,L),$$
(2.15)
As shown in Figure 7 a good data collapse is obtained for $`L^{x_m}\pi _t^{(1)}`$ with the different system sizes, $`L=2^6`$ to $`2^{10}`$. In this figure, the data for the site percolation problem have been shifted by $`\delta _{\mathrm{sb}}`$ given in (2.5) in order to take into account the non-universal constants $`c_\mathrm{s}`$ and $`c_\mathrm{b}`$.

For small values of $`r`$, $`L^zt`$, the boundary conditions are irrelevant since, with a single occupied site in the initial state, the crossing cluster issued from the middle of the sample does not explore the boundary regions. For the same reason $`\pi _t^{(1)}`$ is independent of the size $`L`$ in this regime. Different curves are obtained for the site and bond problems for a given type of boundary conditions, meaning that a new non-universal constant is involved. In a log-log plot, it is clear that a single curve can be obtained through a vertical shift $`\delta _{\mathrm{sb}}^\pi `$ of the data for the site problem. The resulting data collapse is shown in Figure 8. Thus the new non-universal constant $`A`$ is a multiplicative one, either $`A_\mathrm{s}`$ or $`A_\mathrm{b}`$, and the shift is given by:

$$\delta _{\mathrm{sb}}^\pi =\mathrm{ln}\left(\frac{A_\mathrm{b}}{A_\mathrm{s}}\right)=.178(4).$$
(2.16)
It was calculated by first shifting the data of the site problem by $`\delta _{\mathrm{sb}}^r`$. Due to this horizontal shift of the site data, a direct calculation of the vertical shift is not possible. Thus, considering $`y=\mathrm{ln}(L^{x_m}\pi _t^{(1)})`$ as a function $`\mathrm{ln}r`$, the data of the bond problem were fitted to a quadratic polynomial in the linear region, for small values of $`\mathrm{ln}r`$. Then $`\delta _{\mathrm{sb}}^\pi `$ was obtained as the average of $`y_\mathrm{b}y_\mathrm{s}`$, evaluated at small enough values of $`\mathrm{ln}r`$ for the site problem.

Alternatively, the scaling behaviour of $`\pi _t^{(1)}`$ in (2.15) can be recovered by assuming that, asymptotically when $`f1`$, the scale-invariant crossing probability in (2.10) behaves as

$$\underset{f1}{lim}\pi _t(cr,f)=f^{x_m}\phi (cr)=\left(\frac{l}{L}\right)^{x_m}\phi (cr),$$
(2.17)
so that $`L^{x_m}\pi _t^{(1)}=l^{x_m}\phi (cr)`$.

Furthemore the non-universal amplitude, $`A_\mathrm{s}`$ or $`A_\mathrm{b}`$, which appears only in $`\pi _t^{(1)}`$, i.e., when the number $`l`$ of occupied sites in the initial state is O(1), can be explained if one associates to these $`l`$ sites an effective length $`l_{\mathrm{eff}}`$ given by

$$l_{\mathrm{eff}}=l+\delta l,\delta l=\delta l_\mathrm{s},\delta l_\mathrm{b},$$
(2.18)
where the increment $`\delta l`$ is O(1) and non-universal. Thus, replacing $`f`$ by $`f_{\mathrm{eff}}=(l+\delta l)/L`$, there is no change in the FSS limit when $`f>0`$, but when $`l=1`$ we obtain

$$L^{x_m}\pi _t^{(1)}=(1+\delta l)^{x_m}\phi (cr)$$
(2.19)
with a non-universal amplitude given by

$$A=(1+\delta l)^{x_m}.$$
(2.20)
According to (2.16) we have $`\delta l_\mathrm{s}>\delta l_\mathrm{b}`$.

Since, as noticed above, $`\pi _t^{(1)}`$ becomes independent of $`L`$ for small values of $`r`$, the scaling function $`\phi (w)`$ with $`w=cr`$ has to behave as $`w^{x_m/z}`$ when $`w1`$. This leads to a slope $`x_m/z=\beta /z\nu =.159464(5)`$ for the linear part of $`\mathrm{ln}(L^{x_m}\pi _t^{(1)})`$ as a function of $`\mathrm{ln}r`$, in good agreement with the numerical results in Figure 8.

## 3 Critical spanning clusters

In this Section we present a preliminary study of the scaling behaviour of the probability $`\pi _{tn}(r)`$ to find exactly $`n`$ independent clusters surviving after $`t=rL^z`$ time steps on a system with size $`L`$ in the FSS limit. In the epidemic interpretation $`\pi _{tn}`$ is the probability to find $`n`$ disjoint centres of infection surviving in the FSS limit. This problem presents some similarity with the statistics of family names first studied by Galton and Watson watson74  (see also broeker99 ; janssen01 ).

The spanning probability was estimated by generating $`10^6`$ samples with $`L=64,128`$ and $`256`$ for the site problem with FBC or PBC and $`L`$ occupied sites in the initial state. As the size increases, it converges to a scale-invariant function $`\pi _{tn}`$ of the effective aspect ratio $`r_{\mathrm{eff}}`$. A better convergence is obtained with FBC, thus we shall only analyse the results obtained with this type of boundary conditions. An example of the raw data for $`\pi _{tn}`$ at the largest sizes studied is shown in Figure 9.

The scaling behaviour can be obtained by extending to the case of a strongly anisotropic system a simple argument proposed earlier for isotropic percolation aizenman97 ; cardy98 .

Let $`P_t(n,r)`$ denote the probability to have exactly $`n`$ spanning clusters at time $`t`$, on a system in $`d+1`$ dimensions, with volume $`L^d`$. The scale invariance of $`P_t`$ implies that the size dependence is wholly contained in the aspect ratio $`r`$. Dividing the system in $`b^d`$ subsystems with the same volume $`(L/b)^d`$, as shown schematically in Figure 10 for $`d=1`$ and $`b=n`$, the dominant event will correspond to $`n^{}=n/b^d`$ spanning clusters in space-time for each subsytem. Since the aspect ratio of the large system is $`r=t/L^z`$, it becomes $`r^{}=t/(L/b)^z=rb^z`$ in the smaller subsystems. Thus, up to prefactors, one expects the following behaviour:

$$P_t(n,r)\left[P_t(b^dn,b^zr)\right]^{b^d}.$$
(3.1)
According to the transfer-matrix argument of Section 2.1, a leading exponential dependence on $`t`$ (through $`r`$) is expected. Thus, with $`b=n^{1/d}`$, we obtain

$$P_t(n,r)\left[P_t(1,n^{z/d}r)\right]^n\mathrm{exp}(\alpha rn^\omega ).$$
(3.2)
Applying the first relation to the exponential form allows us to identify the exponent $`\omega =1+z/d`$ so that the scale-invariant crossing probability behaves as

$$\pi _{tn}(r)\mathrm{exp}(\alpha rn^{1+z/d}),$$
(3.3)
where $`\alpha `$ is a non-universal prefactor. Actually, the crossing probability $`\pi _{tn}`$, like $`\pi _t`$, is a universal function of the effective aspect ratio $`cr`$. This scaling behaviour generalizes the form conjectured by Aizenman for isotropic percolation in $`D=d+1`$ dimensions aizenman97 , which was recently tested numerically with success shchur02 . The isotropic limit corresponds to $`z=1`$ and $`d=D1`$ in (3.3). The anisotropic result should remain valid below the upper critical dimension $`d_\mathrm{c}=4`$ for directed percolation kinzel83 . At and above $`d_\mathrm{c}`$, as in the isotropic case, one expects a proliferation of spanning clusters, linked with the breakdown of hyperscaling aizenman97 ; aharony84 .

The Monte Carlo results in Figures 11 and 12 for the site and bond percolation problems with FBC at the largest size, $`L=256`$, are in good agreement with this scaling behaviour. The insets show the linear variation of $`\mathrm{ln}\pi _{tn}`$ when $`r`$ is sufficiently large. Different slopes are obtained for different values of $`n`$. As shown in the main figures, after multiplication by $`n^{(z+1)}`$, a good data collapse is obtained.

The universality of $`\pi _{tn}(r_{\mathrm{eff}})`$ is illustrated in Figure 12. The dotted line on the main figure corresponds to the site problem with $`n=2`$ after rescaling $`r`$ by the non-universal factor $`c_\mathrm{s}/c_\mathrm{b}`$. The remaining shift between the site and bond data is likely to be due to finite-size corrections.

## 4 Conclusion

For the directed percolation problem in the FSS limit, the crossing probability in the time direction, $`\pi _t`$, is a scale-invariant universal function of the effective aspect ratio $`r_{\mathrm{eff}}=ct/L^z`$, appropriate for a strongly anisotropic system. It also depends on the fraction $`f=l/L`$ of occupied sites in the initial state, where $`l`$ is the length of the sequence of occupied sites. This function vanishes with $`f`$ as $`f^{x_m}`$, where $`x_m=\beta /\nu `$ is the scaling dimension of the order parameter. When the $`L`$ sites are occupied with probability $`p_i>0`$ at $`t=0`$, the crossing probability is asymptotically the same as for a fully occupied system in the initial state.

Anisotropic scaling, together with the generalization of a simple geometrical argument due to Aizenman for isotropic percolation, leads to the scaling behaviour for the probability $`\pi _{tn}`$ to find $`n`$ critical directed percolation clusters surviving at time $`t`$ on a system with size $`L^d`$. The probability $`\pi _{tn}`$ is a scale-invariant universal function which decays exponentially with $`n^{1+z/d}r_{\mathrm{eff}}`$. The numerical data in $`1+1`$ dimensions support the conjectured expression although further work is needed with larger system sizes, in order to be closer to the true asymptotic behaviour. The scaling behaviour should be also examined in higher dimensions. In both cases one should use a more efficient simulation method, like the “go with the winner” strategy introduced by Grassberger grassberger00 ; grassberger02 ; shchur02 . It would be particularly interesting to study the expected crossover of $`\pi _{tn}`$, due to the proliferation of clusters, at and above the upper critical dimension shchur02 . It should be easier to do this for directed percolation, where $`d_\mathrm{c}=4`$, than for isotropic percolation, where $`D_\mathrm{c}=6`$.

The present work can be also extended by considering the crossing probability on same-spin clusters in strongly anisotropic spin systems at equilibrium, like the three-dimensional uniaxial ANNNI model at its Lifshitz point, for which the critical parameters are now accurately known pleimling01 .


text:
# 𝑆𝑈(2) Landau gluon propagator on a 140³ lattice

## I Introduction

The study of the infra-red (IR) limit of QCD is of central importance for understanding the mechanism of quark confinement and the dynamics of partons at low energy. Despite being non-gauge-invariant, the gluon propagator is a powerful tool in this (non-perturbative) investigation . In particular, it would be interesting to express it in a closed form for recovering the phenomenology of Pomeron exchange from first principles .

Studies of the coupled set of Dyson-Schwinger equations for gluon and ghost propagators in Landau gauge predict for the gluon propagator an IR behavior of the form $`D(p)p^{4\kappa 2}`$ \[implying $`D(0)=0`$ if $`\kappa >0.5`$\]. The available predictions for the IR exponent are $`\kappa [0.52,1.00]`$ in the four-dimensional case  and $`\kappa 0.648`$ or $`\kappa =0.75`$ in three dimensions .

Furthermore, in the minimal Landau gauge, the gauge-fixed configurations belong to the region $`\mathrm{\Omega }`$ of transverse configurations, for which the Faddeev-Popov operator is non-negative. This implies a rigorous inequality  for the Fourier components of the gluon field $`A_\mu (x)`$ and a strong suppression of the gluon propagator in the IR limit. In particular, for dimension $`d`$ and infinite volume, it is proven that the (unrenormalized) gluon propagator is less singular than $`p^{2d}`$ and that, very likely, it vanishes in the IR limit . A vanishing gluon propagator at $`p=0`$, given by the form $`p^2/(p^4+\lambda ^4)`$, was also obtained by Gribov . Here the mass scale $`\lambda `$ arises when the configuration space is restricted to the region $`\mathrm{\Omega }`$. A generalization of this expression has been introduced in Ref.  as an Ansatz for a non-perturbative solution of the gluon Dyson-Schwinger equation.

Numerical studies  have now established that the gluon propagator in lattice Landau gauge shows a turnover in the IR region and attains a finite value for $`p=0`$. Evidences of a decreasing gluon propagator for small $`p`$ have been obtained in the $`4d`$ $`SU(2)`$ and $`SU(3)`$ cases (but only in the strong-coupling regime)  in the $`3d`$ $`SU(2)`$ case (also in the scaling region) , in the $`3d`$ $`SU(2)`$ adjoint Higgs model , in the $`4d`$ $`SU(2)`$ case at finite temperature  and for the equal-time three-dimensional transverse gluon propagator in $`4d`$ $`SU(2)`$ Coulomb gauge . In this last case, one obtains an excellent fit of the transverse propagator by a Gribov-like formula.

This work aims to verify the possibility of using Gribov-like formulae to fit data of the gluon propagator also in Landau gauge. At the same time, we will try to obtain a value for the IR critical exponent $`\kappa `$ to be compared to the analytic determinations mentioned above. In order to probe the infinite-volume limit and the IR region we consider the three-dimensional case and the $`SU(2)`$ group, using lattice sizes up to $`140^3`$. Note that the study of the gluon propagator in three dimensions is also of interest in finite-temperature QCD .

## II Numerical simulations

We consider the standard Wilson action for $`SU(2)`$ lattice gauge theory in three dimensions with periodic boundary conditions. The numerical code is entirely parallelized using MPI. (Technical details and a study of the code performance are left for a subsequent work .) For the construction of staples we follow Ref. . For the random number generator we use a double-precision implementation of RANLUX (version 2.1) with luxury level set to 2. Computations were performed on the PC cluster at the IFSC-USP. The system has 16 nodes and a server with 866 MHz Pentium III CPU and 256/512 MB RAM memory. The machines are connected with a 100 Mbps full-duplex network. The total computer time used for the runs was about 80 days on the full cluster.

In Table I we report, for each coupling $`\beta `$ and lattice volume $`V`$, the parameters used for the simulations. All our runs start with a random gauge configuration. For thermalization we use a hybrid overrelaxed (HOR) algorithm . Each HOR iteration consists of one heat-bath sweep over the lattice followed by $`m`$ micro-canonical sweeps. We did not try to find the best tuning for $`m`$; we use $`m=4`$ for $`V=40^3,\mathrm{\hspace{0.17em}80}^3`$ and $`m=5`$ for $`V=140^3`$. In order to optimize the heat-bath code, we implement two different $`SU(2)`$ generators, namely methods 1 and 2 described in \[20, Appendix A\] with $`h_{cutoff}=2`$.

For the numerical gauge fixing we use the stochastic overrelaxation (SOR) algorithm  with even/odd update. In Table I we report the value of the tuning parameter $`p_{sor}`$ used for each pair $`(\beta \text{,}V)`$. We stop the gauge fixing when the average value of $`[(A)^b(x)]^2`$ is smaller than $`10^{12}`$. (For a definition of the lattice gauge field $`A_\mu ^b(x)`$ and of the lattice divergence $``$ we refer to .) We note that half of the configurations for $`V=80^3`$ were done using the so-called Cornell method , with tuning parameters $`\alpha _{corn}=0.325,0.32,0.316`$ respectively at $`\beta =4.2,5.0,6.0`$. In fact, the Cornell method is somewhat faster than the SOR algorithm and leads to a gauge fixing of comparable quality if one uses an even/odd update . A good estimator of the quality of the gauge fixing is the quantity $`\mathrm{\Sigma }_Q`$ (see eq. 6.8 in Ref. ), which should be zero when the configuration is gauge-fixed. By averaging over the gauge-fixed configurations, we find (at $`\beta =4.2`$) that the ratio between the final and the initial values of $`\mathrm{\Sigma }_Q`$ is (in $`95\%`$ of the cases) about $`5.3\times 10^{10}`$ with the Cornell method and $`2.4\times 10^{11}`$ with the SOR algorithm. At the same time, the average CPU-time needed for updating each site variable is about $`11\%`$ smaller for the Cornell method. In any case, the CPU-time for gauge fixing was quite significant for the large lattices. In order to go to even larger lattices one should probably implement a global algorithm such as the Fourier acceleration method  with the multigrid or conjugate gradient implementations introduced in Ref. , which are highly parallelizable.

We ran the $`40^3`$ lattices on a single node, the $`80^3`$ lattices on two nodes and the $`140^3`$ lattices on four nodes. The parallelization of the code worked well. In fact, for runs on 2 and 4 nodes, we obtain speed-up factors 1.82 and 3.41 for the heat-bath link update. For the micro-canonical link update the factors are respectively 1.87 and 3.77 and for the SOR site update we get 1.90 and 3.72.

For each $`\beta `$ we evaluate the average plaquette $`W_{1,1}`$ (see Table II). Results are in agreement with the data reported in Ref. , but we now have smaller statistical errors. (The data have been analyzed using various methods, described in footnote 4 of Ref. . Here we always report the largest error found.) We also evaluate the tadpole-improved coupling $`\beta _I\beta W_{1,1}`$. In this way, by using the fit given in eq. 2 and Table IV of Ref.  we calculate the string tension $`\sqrt{\sigma }`$ in lattice units (see Table II) and the inverse lattice spacing $`a^1`$ using the input value $`\sqrt{\sigma }=0.44`$ GeV. The fit is valid for $`\beta \mathrm{}>3.0`$, i.e. the couplings $`\beta `$ considered here are well above the strong-coupling region. Let us notice that, if we compare the data for the string tension (in lattice units) with data obtained for the $`SU(2)`$ group in four dimensions \[25, Table 3\], then our values of $`\beta `$ correspond to $`\beta 2.28,2.345,2.41`$ in the four dimensional case. Finally, in the same table we report the lattice spacing $`a`$ in fm and the smallest non-zero momentum (in MeV) that can be considered for each $`\beta `$. Thus, with the lattice volumes and the $`\beta `$ values used here we are able to consider momenta as small as 59 MeV (in the deep IR region) and physical lattice sides almost as large as 25 fm.

In this work we did not do a systematic study of Gribov-copy effects for the gluon propagator. However, we compared data obtained using the SOR and the Cornell gauge-fixing methods. In principle, different gauge-fixing algorithms — or even the same algorithm with different values of the tuning parameter  — can generate different Gribov copies starting from the same thermalized configuration. Thus, this comparison provides an estimate of the bias (Gribov noise) introduced by the gauge-fixing procedure. We found that, in most cases, the difference between the two sets of gluon-propagator data is within 1 standard deviation and that, in all cases, it is smaller than 2 standard deviations. Moreover, this difference did not show any systematic effect, suggesting that the influence of Gribov copies on the gluon propagator (if present) is of the order of magnitude of the numerical accuracy. This is in agreement with previous studies in Landau gauge for the $`SU(2)`$ and $`SU(3)`$ groups in four dimensions . A similar result has also been obtained for the Coulomb gauge . Note that, in the $`U(1)`$ case , Gribov copies can affect the behavior of the photon propagator, making it difficult to reproduce the known perturbative behavior in the Coulomb phase. The situation is very different for the $`SU(2)`$ and $`SU(3)`$ cases, at least when considering the lattice Landau gauge. In fact, as said in the Introduction, in the non-Abelian case the minimal-Landau-gauge condition implies  the positiveness of the Faddeev-Popov matrix and a strict bound for the Fourier components of the gluon field $`A_\mu (k)`$. The bound applies to all Gribov copies obtained with the numerical gauge fixing. Thus, if the behavior of the gluon propagator is determined by this bound , then this behavior should be the same for all lattice Gribov copies. This would explain why we do not see Gribov-copy effects here. Clearly, the same result does not apply to the $`U(1)`$ theory, since in this case the Faddeev-Popov matrix is independent of the gauge field.

## III Results and Conclusions

We evaluate the lattice gluon propagator $`D(k)`$ and study it as a function of the lattice momentum $`p^2(k)`$ (see Ref.  for definitions). In our simulations we consider, for each gauge-fixed configuration, all vectors $`k(k_x,k_y,k_t)`$ with only one component different from zero and average over the three directions. For the gluon propagator we analyze the data by estimating the statistical error with three different methods: standard deviation, jack-knife with single-data elimination and bootstrap (with 10000 samples). We found that the results obtained are in agreement in all cases. Here we always use the standard-deviation error.

Let us recall that in the three-dimensional case the coupling $`g^2`$ has dimensions of mass. Thus, in order to obtain a dimensionless lattice coupling we have to set $`\beta =4/(ag^2)`$. Then, with our notation , the quantity $`aD(k)`$ approaches $`g^2D^{(cont)}(k)/4`$ in the continuum limit, where $`D^{(cont)}(k)`$ is the unrenormalized continuum gluon propagator.

In order to compare lattice data at different $`\beta `$’s, we apply the matching technique described in \[30, Sec. V.B.2\]. (Note that we have already determined the lattice spacing $`a`$, as described above.) We start by checking for finite-size effects, comparing data at different lattice sizes and same $`\beta `$ value. In this way, we find (for each $`\beta `$) a range of ultra-violet (UV) momenta for which the data are free from finite-volume corrections. We then perform the matching using data for these momenta and $`V=40^3`$, since for this lattice volume the errors are smallest (about $`1\%`$). In particular, when matching data obtained at two different values of $`\beta `$, we first interpolate the data for the larger $`\beta `$ (the fine lattice) using a spline. Finally, we find the multiplicative factor $`R_Z=Z(a_f)/Z(a_c)`$ corresponding to the best fit of the (multiplied) coarse-lattice data to the same spline. The error of $`R_Z`$ is estimated using a procedure similar to the one described in . The method works very well (see Figure 1). Notice that we did not fix the remaining global factor $`Z`$ imposing a renormalization condition, as done for example in Ref. . Our case is equivalent to setting $`Z(a)=1`$ at $`\beta =6.0`$.

The data obtained after the matching are shown for $`V=80^3`$ in Figure 1. Clearly, we find that the gluon propagator decreases in the IR limit for momenta smaller than $`p_{dec}`$, which corresponds to the mass scale $`\lambda `$ in a Gribov-like propagator. From the plot we can estimate $`p_{dec}=350_{50}^{+100}\text{MeV}`$, in agreement with Ref. .

In Figure 2 we plot the rescaled gluon propagator at zero momentum, namely $`aD(0)/Z(a)`$, as a function of the inverse lattice side $`L^1=1/(aN)`$ in physical units (fm-1). We see that $`aD(0)/Z(a)`$ decreases monotonically as $`L`$ increases, in agreement with Ref. . It is interesting to notice that these data can be well fitted using the simple Ansatz $`d+b/L^c`$ both with $`d=0`$ and $`d0`$ (see Figure 2). In order to decide for one or the other result one should go to significantly larger lattice sizes. We plan  to extend these simulations to $`\beta =3.4`$ and lattice sizes up to $`260^3`$, allowing us to consider a value $`L^10.017`$ fm-1. (This requires running in parallel on all nodes of our PC cluster.)

Following Ref.  we fit the data using a Gribov-like (or Stingl-like) formula

$$D(p)=\frac{s+zp^{2\alpha }}{y^2+(p^2+x)^2}\text{,}$$
(1)
where $`z,s,\alpha ,x`$ and $`y`$ are fitting parameters. For non-negative $`\alpha `$ this implies a finite gluon propagator in the IR limit, with a behavior given by $`D(p)(s+zp^{2\alpha })`$. If $`y^2>0`$ this form corresponds to a propagator with poles at $`m_\pm ^2=x\pm iy`$, while if $`y^20`$ the poles are real. Finally, note that the IR exponent $`\kappa `$ considered in studies of Dyson-Schwinger equations is given in terms of $`\alpha `$ by $`\kappa =(1+\alpha )/2`$, assuming $`D(0)=0`$. Results of our fits using the un-rescaled lattice data $`D(k)`$ are reported in Table III. We see that $`\alpha `$ decreases when the physical lattice volume increases. Also, we always get $`y^2>0`$ and $`x<\left|y\right|`$, which seems to support the scenario of purely imaginary poles found in . Let us recall that the poles of the gluon propagator are gauge-invariant (at all orders in perturbation theory) .

Notice that this fitting form leads in general to the wrong UV behavior, namely $`D(p)p^{2\alpha 4}`$. This is not a serious problem, since the largest momentum we can consider is about $`3.5`$ GeV, i.e. we are not really exploring the UV limit. On the other hand, the exponent $`\alpha `$ in (1) plays a role both in the IR and in the UV regimes. Thus, the values obtained for $`\alpha `$ probably correspond to averaging the behavior of the gluon propagator in these two regions. In particular, since we expect $`D(p)p^2`$ in the UV limit, it is likely that the IR behavior of the propagator be given by a smaller exponent $`\alpha `$ than the ones reported in Table III. To check this we set $`x=0`$ and introduce an anomalous dimension $`\gamma `$:

$$D(p)=\frac{s+zp^{2\alpha }}{y^2+p^{2(1+\gamma )}}\text{.}$$
(2)
Results for this fitting form are also reported in Table III. Indeed, we get values of $`\alpha `$ smaller than in the previous case (and still decreasing with increasing physical lattice volume). For the anomalous dimension we obtain $`\gamma 0.65`$, with small volume- and $`\beta `$-dependence. The problem with eq. (2) is that the introduction of $`\gamma `$ compromises the pole interpretation.

Finally, one can also try the form

$$D(p)=\frac{s+zp^{2\alpha }}{(y^2+p^4)^\gamma }\text{.}$$
(3)
This corresponds to a propagator with purely imaginary poles $`m_\pm ^2=\pm iy`$ and at the same time allows the data to select the IR and the UV behaviors separately. The problem in this case is that the fit is unstable for small lattice volumes. On the contrary, for $`V=140^3`$ the fit works well and we obtain $`\alpha =0.27(6)`$, $`0.29(7)`$, $`0.38(8)`$ respectively for $`\beta =4.2,5.0,6.0`$ and $`\gamma 0.72`$. Hence, our $`\alpha `$ values are of the order of $`0.3`$, corresponding to $`\kappa 0.65`$. (Again, there is a decrease of $`\alpha `$ when the physical lattice volume increases.)

In order to check for possible effects from the breaking of rotational invariance  we redid our fits substituting $`p^2(k)`$ by $`\stackrel{~}{p}^2(k)p^2(k)+p^{[4]}(k)/12`$ (see Ref. ). One expects this modification to play an important role in the UV limit and to have a small effect in the IR case. For all $`\beta `$ values and lattice volumes we obtain good fits and results similar to those reported above. In particular, using eq. (1) one still finds $`y^2>0`$ and $`x<\left|y\right|`$, supporting the scenario of purely imaginary poles. At the same time, for the three fitting functions, $`\alpha `$ decreases when the physical lattice volume increases, but its value is about $`2030\%`$ larger than the one reported in Table (III) and above. In particular, using eq. (3) and data for $`V=140^3`$ we obtain $`\alpha =0.32(7)`$, $`0.39(8)`$, $`0.59(12)`$ respectively for $`\beta =4.2,5.0,6.0`$ and $`\gamma 0.7`$. This implies slightly larger values for $`\kappa =(1+\alpha )/2`$. The analysis of these discretization effects and their influence on the extrapolation of $`\alpha `$ and $`\kappa `$ to the continuum limit will be done elsewhere .

We have confirmed, by numerical simulations in the scaling region, that the transverse gluon propagator in $`3d`$ $`SU(2)`$ Landau gauge is a decreasing function for momenta $`p\mathrm{}<350\text{MeV}`$ (attaining a finite value at $`p=0`$). Also, the data are well fitted by the Gribov-like formulae (1)–(3) and we obtain an IR critical exponent $`\kappa `$ in agreement with recent analytic results. In order to eliminate remaining discretization and finite-volume effects \[and in particular to check if $`D(0)=0`$\], we need to simulate at larger lattice volumes and at other values of $`\beta `$.

## Acknowledgments

We thank Reinhard Alkofer, Terry Goldman, Kurt Langfeld, Axel Maas and Dan Zwanziger for helpful comments and suggestions and Martin Lüscher for sending us the latest double-precision version of the RANLUX random number generator. The research of A.C. and T.M. is supported by FAPESP (Project No. 00/05047-5). A.T. thanks CNPq for financial support and the IFSC-USP for the hospitality.


text:
# Multiparticle correlations from momentum conservation

## I Multiparticle cumulants

Consider a collision in which a total of $`N`$ particles are emitted (throughout this paper, $`N`$ is assumed to be large). Let $`𝐩_1`$, …, $`𝐩_N`$, be their $`D`$-dimensional momenta, where the space dimension $`D`$ is kept arbitrary for sake of generality. For a given $`M<N`$, we denote by $`f(𝐩_{j_1},\mathrm{},𝐩_{j_M})`$ the $`M`$-particle probability distribution of $`𝐩_{j_1}`$, …, $`𝐩_{j_M}`$.

To derive multiparticle correlations of arbitrary order in a systematic way, we introduce a generating function of the distributions, which we define as

$`G(x_1,\mathrm{},x_N)`$ $``$ $`1+x_1f(𝐩_1)+x_2f(𝐩_2)+\mathrm{}`$ (2)
$`+x_1x_2f(𝐩_1,𝐩_2)+\mathrm{},`$
and so on for every order $`M`$. Thus, the coefficient of the term $`x_{j_1}x_{j_2}\mathrm{}x_{j_M}`$ is the $`M`$-particle probability distribution $`f(𝐩_{j_1},𝐩_{j_2},\mathrm{},𝐩_{j_M})`$.

We shall mostly be interested in the connected part of the $`M`$-particle distribution, that is, the term which cannot be expressed as a product of correlations between less than $`M`$ particles. For instance, the two-particle distribution can be decomposed into two parts:

$$f(𝐩_1,𝐩_2)=f(𝐩_1)f(𝐩_2)+f_c(𝐩_1,𝐩_2),$$

where the first term in the right-hand side (rhs) is the mere product of the one-particle probability distributions, while the second one, the “connected part,” reflects correlations between $`𝐩_1`$ and $`𝐩_2`$. In particular, $`f_c(𝐩_1,𝐩_2)`$ vanishes if $`𝐩_1`$ and $`𝐩_2`$ are uncorrelated. This genuine correlation, which is also called the cumulant of the distribution in probability theory, is obtained by taking the logarithm of the generating function of distributions vanKampen . The cumulant $`f_c(𝐩_{j_1},\mathrm{},𝐩_{j_M})`$ is the coefficient of $`x_{j_1}\mathrm{}x_{j_M}`$ in the expansion of $`\mathrm{ln}G(x_1,\mathrm{},x_N)`$:

$`\mathrm{ln}G(x_1,\mathrm{},x_N)`$ $``$ $`x_1f_c(𝐩_1)+x_2f_c(𝐩_2)+\mathrm{}`$ (4)
$`+x_1x_2f_c(𝐩_1,𝐩_2)+\mathrm{}.`$
If the system splits into independent subsystems or, more generally, if correlations in the system are short-range, the probability distributions of each subsystem add up: $`f(\{𝐩_j\})=_A(N_A/N)f_A(\{𝐩_j\})`$, where the sum runs over subsystems with sizes $`N_A`$. It follows that $`G`$ can be factorized into the product of functions for each subsystem: $`G(\{x_j\})=_Ag_A(\{N_Ax_j/N\})`$. In the limit where each individual particle is a subsystem, $`G(\{x_j\})=[g(\{x_j/N\})]^N`$. Therefore, $`\mathrm{ln}G`$, and the cumulants, are the sums of the corresponding quantities for the subsystems. As a result, inspecting the term in $`x_{j_1}\mathrm{}x_{j_M}`$ shows that, while the $`M`$-particle probability distribution $`f(𝐩_{j_1},\mathrm{},𝐩_{j_M})`$ is independent of $`N`$, the cumulant $`f_c(𝐩_{j_1},\mathrm{},𝐩_{j_M})`$ scales like $`1/N^{M1}`$. In Sec. II, we shall show that the cumulants from the correlations due to global momentum conservation follow the same scaling to leading order in $`1/N`$. More precisely, we shall demonstrate that the generating function of distributions reads

$$G(x_1,\mathrm{},x_N)e^{Ng(\{x_j/N\})}\left(1+\frac{\{x_j\}^l}{N^q}\right),$$

where the sum runs over terms with $`ql`$.

## II Cumulants in the large-$`N`$ limit

In this Section, we derive a compact expression of the cumulants of multiparticle correlations due to global momentum conservation to leading order in $`1/N`$. In particular, we show that the $`M`$-particle cumulants scale like $`1/N^{M1}`$.

Let us assume for simplicity that the only source of correlation between the particles is global momentum conservation, $`𝐩_1+\mathrm{}+𝐩_N=\mathrm{𝟎}`$. Under this assumption, the $`M`$-particle distribution of $`𝐩_1`$, …, $`𝐩_M`$, with $`M<N`$, is defined as

$$f(𝐩_1,\mathrm{},𝐩_M)\frac{\left({\displaystyle \underset{j=1}{\overset{M}{}}}F(𝐩_j)\right){\displaystyle \delta ^D(𝐩_1+\mathrm{}+𝐩_N)\underset{j=M+1}{\overset{N}{}}\left[F(𝐩_j)d^D𝐩_j\right]/𝒩_D^{NM}}}{{\displaystyle \delta ^D(𝐩_1+\mathrm{}+𝐩_N)\underset{j=1}{\overset{N}{}}\left[F(𝐩_j)d^D𝐩_j\right]/𝒩_D^N}},$$
(5)
where $`F(𝐩)`$ is the one-particle momentum distribution unrenormalized for the momentum conservation constraint, and $`𝒩_DF(𝐩)d^D𝐩`$ is a normalization constant. In the following, we denote by $`\mathrm{}`$ an $`F(𝐩)`$-weighted average, that is, $`g(𝐩)g(𝐩)F(𝐩)d^D𝐩/𝒩_D`$, for any function of momentum $`g(𝐩)`$.

The denominator in Eq. (5) is a constant, which we denote by $`1/𝒞_D`$, independent of $`M`$. The actual value will not influence the forthcoming discussion.

Consider next the numerator of Eq. (5). Introducing a Fourier representation of the Dirac distribution, it reads

$$\frac{d^D𝐤}{(2\pi )^D}\left(\underset{j=1}{\overset{M}{}}F(𝐩_j)e^{i𝐤𝐩_j}\right)e^{i𝐤𝐩}^{NM}=\frac{d^D𝐤}{(2\pi )^D}e^{i𝐤𝐩}^N\left(\underset{j=1}{\overset{M}{}}F(𝐩_j)\frac{e^{i𝐤𝐩_j}}{e^{i𝐤𝐩}}\right).$$
(6)
Inserting in Eq. (2) the expression of the $`M`$-particle distribution, Eq. (5), with the numerator replaced by the rhs of Eq. (6), we obtain

$`G(x_1,\mathrm{},x_N)`$ $`=`$ $`𝒞_D{\displaystyle \frac{d^D𝐤}{(2\pi )^D}e^{i𝐤𝐩}^N\underset{j=1}{\overset{N}{}}\left(1+x_jF(𝐩_j)\frac{e^{i𝐤𝐩_j}}{e^{i𝐤𝐩}}\right)}`$ (7)
$``$ $`𝒞_D{\displaystyle \frac{d^D𝐤}{(2\pi )^D}e^{i𝐤𝐩}^N\mathrm{exp}\left(\underset{j=1}{\overset{N}{}}x_jF(𝐩_j)\frac{e^{i𝐤𝐩_j}}{e^{i𝐤𝐩}}\right)}.`$ (8)
In passing from the first line to the second one, we have used the fact that we shall only consider the coefficient of $`x_{j_1}\mathrm{}x_{j_M}`$ where the $`M`$ indices $`j_1`$, …, $`j_M`$ are all different \[see Eq. (4)\]. One easily checks that the only difference between the two forms of the generating function in Eq. (7) comes from terms in which at least one $`x_j`$ is raised to some power $`m2`$ (corresponding to the “autocorrelation” of particle $`j`$ with itself). Therefore, as far as we are concerned, Eq. (7) really is an identity, not an approximation.

Please note that in Eq. (7), the variable $`x_j`$ is always multiplied by a factor $`F(𝐩_j)`$. We may thus rescale $`x_j`$ by this factor, and drop it in the following to simplify expressions. This is quite satisfactory, since it means that the measurable multiparticle correlations $`f(𝐩_{j_1},\mathrm{},𝐩_{j_M})`$ or $`f_c(𝐩_{j_1},\mathrm{},𝐩_{j_M})`$ will not depend on the non-measurable distribution $`F(𝐩)`$—this justifies our previously calling it “unrenormalized”.

To evaluate the integral in Eq. (7), we rely on the fact that $`N`$ is large, and use a saddle-point approximation. To simplify the discussion, we introduce the notation

$$(𝐤)\mathrm{ln}e^{i𝐤𝐩}+\underset{j=1}{\overset{N}{}}\frac{x_j}{N}\frac{e^{i𝐤𝐩_j}}{e^{i𝐤𝐩}}.$$
(9)
Thus, the integrand in Eq. (7) is $`e^{N(𝐤)}`$. We call $`𝐤_0`$ the position of the saddle point. Please note that $``$ and its successive derivatives, which we shall denote by $`^{}`$, $`^{\prime \prime }`$, $`^{(3)}`$, …, depend on $`x`$ only through $`x/N`$, where $`x`$ stands for any of the $`x_j`$. Therefore, $`𝐤_0`$, which is of course the solution of $`^{}(𝐤)=0`$, also depends on $`x/N`$ only. One should pay attention to the fact that the saddle point $`𝐤_0`$ is not merely the origin $`𝐤=\mathrm{𝟎}`$.

We shall now demonstrate that to leading order in $`1/N`$, all multiparticle cumulants are determined by the saddle-point value $`N(𝐤_0)`$. More precisely, we show that the $`M`$-particle cumulant is of order $`1/N^{M1}`$, and that corrections to the saddle-point calculation only yield subleading terms, suppressed by (positive) powers of $`1/N`$.

Using a Taylor expansion of $`(𝐤)`$ around the saddle point, the generating function (7) reads

$`G(x_1,\mathrm{},x_N)`$ $`=`$ $`𝒞_De^{N(𝐤_0)}{\displaystyle \frac{d^D𝐤}{(2\pi )^D}e^{N^{\prime \prime }(𝐤_0)(𝐤𝐤_0)^2/2}\mathrm{exp}\left[N\underset{m=3}{\overset{+\mathrm{}}{}}\frac{^{(m)}(𝐤_0)}{m!}(𝐤𝐤_0)^m\right]}`$ (10)
$`=`$ $`𝒞_De^{N(𝐤_0)}{\displaystyle \frac{d^D𝐤}{(2\pi )^D}e^{N^{\prime \prime }(𝐤_0)(𝐤𝐤_0)^2/2}\left(1+\underset{n=1}{\overset{+\mathrm{}}{}}\frac{1}{n!}\left[N\underset{m=3}{\overset{+\mathrm{}}{}}\frac{^{(m)}(𝐤_0)}{m!}(𝐤𝐤_0)^m\right]^n\right)}.`$ (11)
As recalled in Sec. I, the cumulants are given by the logarithm of the generating function. Now, the logarithm of Eq. (10) will split into three parts. First, there is $`\mathrm{ln}𝒞_D`$, which does not depend on $`x`$, and thus does not influence the values of the multiparticle cumulants.

Next, there is $`N(𝐤_0)`$. As noted above, both $``$ and $`𝐤_0`$ only involve powers of $`x/N`$. Thus, the coefficient of $`x_{j_1}\mathrm{}x_{j_M}`$ in $`(𝐤_0)`$ contains a factor $`1/N^M`$. Multiplying by the overall factor of $`N`$, we find that the contribution of $`N(𝐤_0)`$ to the cumulant $`f_c(x_{j_1}\mathrm{}x_{j_M})`$ scales as $`1/N^{M1}`$.

Finally, the cumulants involve the contribution of the logarithm of the integral in Eq. (10). We shall use the fact that the integrand is the sum of a Gaussian function and of combinations of its moments. After integration, this yields a common factor, $`1/[2\pi N^{\prime \prime }(𝐤_0)]^{D/2}`$, which multiplies a sum $`\mathrm{\Sigma }`$. Let us show that $`\mathrm{\Sigma }`$ only involves terms in $`x^l/N^q`$ with $`ql`$. Noting that the values which contribute to the Gaussian integral are of order $`(𝐤𝐤_0)^21/N`$, one sees that each term $`N^{(m)}(𝐤_0)(𝐤𝐤_0)^m`$ is of order $`N^{1m/2}`$, multiplied by a function of $`x/N`$. Since the sum runs over $`m3`$, $`N^{1m/2}N^{1/2}`$: therefore, a power $`x^l`$ goes with at most a factor $`1/N^{l+1/2}`$. Even after raising the sum to the power $`n`$, and taking the logarithm, this will always remain of the form $`x^l/N^q`$ with $`q>l`$. Finally, the logarithm of the overall factor $`1/[2\pi N^{\prime \prime }(𝐤_0)]^{D/2}`$ involves $`x`$ only through the form of powers of $`x/N`$ in $`\mathrm{ln}^{\prime \prime }(𝐤_0)`$. All in all, the contribution to the $`M`$-particle cumulant of the integral is therefore at most of order $`1/N^M`$, subleading with respect to the contribution of $`N(𝐤_0)`$.

Therefore, we have shown that the leading contribution to the cumulants of multiparticle correlations due to momentum conservation is the saddle-point value: $`\mathrm{ln}G(x_1,\mathrm{},x_N)=N(𝐤_0)`$. This means first that the mere knowledge of $`𝐤_0`$ gives access to all cumulants at once, at least to leading order. One easily checks that $`𝐤_0`$ is pure imaginary, i.e., $`i𝐤_0`$ has real-valued components. From Eq. (9), it follows that $`(𝐤_0)`$ is real-valued as well, and so are the cumulants, as should be. We shall illustrate these points by computing explicitly the two- and three-particle cumulants in the following Section. Our result also means that genuine $`M`$-particle correlations arising from momentum conservation, which is a long-range effect, scale in the same way as correlations from short-range sources, namely as $`𝒪(1/N^{M1})`$. This was certainly not obvious a priori.

## III Two- and three-particle cumulants

As an example of the order of magnitude derived in the previous Section, let us compute the two- and three-particle cumulants. According to the discussion in Sec. I, they are given by the terms in $`x^2`$ and $`x^3`$ in the generating function of cumulants $`\mathrm{ln}G(x_1,\mathrm{},x_N)`$, that is, following Sec. II, the corresponding terms in $`N(𝐤_0)`$. We shall for simplicity perform calculations assuming that $`F(𝐩)`$ only depends on the modulus $`|𝐩|`$, not on the azimuthal angle of $`𝐩`$. As a consequence, the average momentum $`𝐩`$ vanishes. Departures from this assumption are discussed at the end of this Section.

Now a straightforward calculation shows that the saddle point $`𝐤_0`$ is given by

$$\left(\underset{j=1}{\overset{N}{}}\frac{x_j}{N}\frac{e^{i𝐤_0𝐩_j}}{e^{i𝐤_0𝐩}}1\right)𝐩e^{i𝐤_0𝐩}=\underset{j=1}{\overset{N}{}}\frac{x_j}{N}𝐩_je^{i𝐤_0𝐩_j}.$$
(12)
As mentioned above, $`𝐤_0`$ is a function of $`x/N`$. Equation (12) can be solved order by order in $`x/N`$, using the fact that one term (in the left-hand side) is of order $`0`$ in $`x/N`$ while the other two are linear in $`x/N`$. Inspecting Eq. (9), and using the fact that the first term in the rhs is even in $`𝐤`$, we see that the calculation of $`(𝐤_0)`$ to order $`x^3`$ requires our knowing $`𝐤_0`$ to order $`x^2`$ (while the calculation to order $`x^2`$ only requires $`𝐤_0`$ to order $`x`$).

We can then compute $`𝐤_0`$, expanding Eq. (12) in powers of $`𝐤_0`$, which gives

$$i𝐤_0=\left[1_D\left(X_01_D\frac{D}{𝐩^2}X_2\right)\right]^1\frac{D}{𝐩^2}𝐗_\mathrm{𝟏},$$
(13)
where $`1_D`$ denotes the unit $`D\times D`$ matrix and we have introduced

$$X_0\underset{j=1}{\overset{N}{}}\frac{x_j}{N},𝐗_\mathrm{𝟏}\underset{j=1}{\overset{N}{}}\frac{x_j}{N}𝐩_j,X_2\underset{j=1}{\overset{N}{}}\frac{x_j}{N}𝐩_j𝐩_j.$$

Note that $`i𝐤_0`$ is real-valued, as expected. Reporting its value in Eq. (9), we obtain

$`(𝐤_0)`$ $`=`$ $`X_0{\displaystyle \frac{D}{2𝐩^2}}(𝐗_\mathrm{𝟏})^2`$
$`{\displaystyle \frac{D}{2𝐩^2}}𝐗_\mathrm{𝟏}\left(X_01_D{\displaystyle \frac{D}{𝐩^2}}X_2\right)𝐗_\mathrm{𝟏}.`$
Multiplying this result by $`N`$, we finally obtain, to leading order in $`1/N`$:

$`\mathrm{ln}G(x_1,\mathrm{},x_N)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}x_j{\displaystyle \frac{D}{2N𝐩^2}}{\displaystyle \underset{j,k}{}}x_jx_k(𝐩_j𝐩_k)`$ (15)
$`{\displaystyle \frac{D}{2N^2𝐩^2}}{\displaystyle \underset{j,k,l}{}}x_jx_kx_l\left[𝐩_j𝐩_l{\displaystyle \frac{D}{𝐩^2}}(𝐩_j𝐩_k)(𝐩_k𝐩_l)\right]+𝒪(x^4).`$
Hence the coefficients of $`x_1x_2`$ and $`x_1x_2x_3`$:

$`f_c(𝐩_1,𝐩_2)`$ $`=`$ $`{\displaystyle \frac{D𝐩_1𝐩_2}{N𝐩^2}},`$ (16)
$`f_c(𝐩_1,𝐩_2,𝐩_3)`$ $`=`$ $`{\displaystyle \frac{D}{N^2𝐩^2}}(𝐩_1𝐩_2+𝐩_1𝐩_3+𝐩_2𝐩_3)`$ (18)
$`+{\displaystyle \frac{D^2}{N^2𝐩^2^2}}\left[(𝐩_1𝐩_2)(𝐩_1𝐩_3)+(𝐩_1𝐩_2)(𝐩_2𝐩_3)+(𝐩_1𝐩_3)(𝐩_2𝐩_3)\right].`$
We recover, in the case $`D=2`$, the expression of the two-particle correlation due to transverse momentum conservation already derived in Refs Borghini:2000cm ; Danielewicz:1988in : the correlation is back-to-back, and stronger between high-momenta particles.

Finally, let us comment on the assumption that $`F(𝐩)`$ only depends on the modulus of $`𝐩`$, not on its azimuthal angle.

Our expansion of Eq. (12) relies on both $`𝐩=\mathrm{𝟎}`$ and on the identity $`(𝐤𝐩)^2=𝐤^2𝐩^2/D`$, which are no longer valid if $`F`$ is non-isotropic. In the general case, $`𝐤^2𝐩^2/D`$ will be replaced by another quadratic form in $`𝐤`$, namely $`𝐤𝐩^{}𝐩^{}𝐤`$, where $`𝐩^{}𝐩𝐩`$. To perform explicit calculations, it becomes necessary to introduce the principal-axes frame, in which $`𝐩^{}𝐩^{}`$ is diagonal. For instance, the two-particle cumulant, Eq. (16), will read

$$f_c(𝐩_1,𝐩_2)=\underset{i=1}{\overset{D}{}}\frac{(𝐩_1)_i(𝐩_2)_i}{N(𝐩)_i^2},$$
(19)
where the sum runs over the coordinates along the principal axes. The distinction between the various directions is for example relevant in a high-energy collision: while the transverse momentum distribution can be isotropic (unless there is some anisotropy, as, e.g., a correlation to the impact parameter direction), the isotropy will not extend to the beam direction. In that case, one may want to use Eq. (16), with $`D=2`$, when studying two-particle correlations in the transverse plane, but turn to Eq. (19) if interested in $`3`$-dimensional correlations.

## IV Discussion

We have shown how it is possible to calculate the multiparticle correlations arising from momentum conservation between any number of particles using a generating function formalism:

$$\mathrm{ln}G(x_1,\mathrm{},x_N)=N(𝐤_0),$$
(20a)
where $`(𝐤)`$ is given by Eq. (9) and $`𝐤_0`$ by
$$^{}(𝐤_0)=0.$$
(20b)
In particular, we have seen that the $`M`$-particle cumulant scales as $`1/N^{M1}`$, where $`N`$ is the total number of emitted particles. That means that the correlation scales in the same way as correlations arising from short-range final interactions or from resonance decays, although the underlying reason is less obvious. In the latter cases, the scaling is a simple consequence of combinatorics: the $`M`$ particles are required either to be altogether in a small phase-space region, or to originate from a single decay, hence the factor $`1/N^{M1}`$. However, in the case of global momentum conservation, the same arguments do not apply a priori.

The scaling of the genuine $`M`$-particle correlation due to momentum conservation has several consequences. First, a good feature: in the context of collective flow in heavy-ion collisions JY:QM97 , a new method of analysis has been proposed Borghini:2001vi , which relies on the idea that a cumulant expansion allows one to separate genuine collective phenomena from trivial short-range correlations, while they interfere in the data. Since correlations from momentum conservation behave as short-range correlations, they can be removed as efficiently as them, to leave a clean collective flow signal.

Then, some bad news. Since the $`M`$-particle correlation from momentum conservation is of the same order as that from short-range interactions, as for instance quantum correlations, it may bias measurements of these other correlations. It is thus worth checking that momentum conservation does not contribute significantly to the correlations measured by particle interferometry and attributed to quantum (anti)symmetrization of the wave-function. This is especially true when studying the dependence of HBT parameters on the average momentum of the particles since the correlation from momentum conservation increases with momentum \[see Eqs. (16) and (18)\].

###### Acknowledgements.

Fruitful discussions with François Gelis and Jean-Yves Ollitrault are warmly acknowledged.


text:
# The Clairvoyant Demon Has a Hard Task

## 1. Statement of the result

### 1.1. Introduction

Let $`x=(x(0),x(1),\mathrm{})`$ be an infinite sequence and $`u=(u(0),u(1),\mathrm{})`$ be a binary sequence with elements in $`\{0,1\}`$. Let $`s_n=_{i=0}^{n1}u(i)`$. We define the delayed version $`x^{(u)}`$ of $`x`$, by

$$x^{(u)}(n)=x(s_n).$$

Thus, if $`i=s_n`$ then $`x^{(u)}(n)=x(i)`$, and $`x^{(u)}(n+1)=x(i)`$ or $`x(i+1)`$ depending on whether $`u(n)=0`$ or 1. If $`u(n)=0`$ then we can say that $`x^{(u)}`$ is delayed at time $`n+1`$. For two infinite sequences $`x`$, $`y`$ we say that they do not collide if there is a delay sequence $`u`$ such that for each $`n`$ we have

$$x^{(u)}(n)y^{(1u)}(n).$$

Here, $`1u`$ is the delay sequence complementary to $`u`$: thus, $`y^{(1u)}`$ is delayed at time $`n`$ if and only if $`x^{(u)}`$ is not.

For a given $`m>1`$, suppose that $`X=(X(0),X(1),\mathrm{})`$ is an infinite sequence of independent random variables, and $`Y=(Y(0),Y(1),\mathrm{})`$ is another such sequence, also independent of $`X`$, where all variables are uniformly distributed over $`\{0,\mathrm{},m1\}`$. The following theorem has been conjectured in .

###### Proposition 1.1 (See ).

If $`m`$ is sufficiently large then with positive probability, $`X`$ does not collide with $`Y`$.

The sequences $`X,Y`$ can be viewed as two independent random walks on the complete graph $`K_m`$, and then the problem is whether a “clairvoyant demon”, i.e. a being who knows in advance both infinite sequences $`X`$ and $`Y`$, can introduce delays into these walks in such a way that they never collide.

### 1.2. Graph reformulation

We define a graph $`G=(V,E)`$ as follows. $`V=_{>0}^2`$ is the set of points $`(i,j)`$ where $`i,j`$ are positive integers. Let us define the distance of two points $`(i,j)`$, $`(k,l)`$ as $`|ki|+|lj|`$ ($`L_1`$ distance). When representing the set $`V`$ of points $`(i,j)`$ graphically, the right direction is the one of growing $`i`$, and the upward direction is the one of growing $`j`$. The set $`E`$ of edges consists of all pairs of the form $`((i,j),(i+1,j))`$ and $`((i,j),(i,j+1))`$.

Given $`X,Y`$ as in the theorem, let us say that a point $`(i,j)`$ has color $`k`$ if $`X(i)=Y(j)=k`$. Otherwise, it has color $`1`$, which we will call white. It is easy to see that $`X`$ and $`Y`$ do not collide if and only if there is an infinite directed path in $`G`$ starting from $`(0,0)`$ and proceeding on white points. Indeed, each path corresponds to a delay sequence $`u`$ such that $`u(n)=1`$ if and only if the edge is horizontal. Thus, the two sequences do not collide if and only if the graph of white points “percolates”. We will say that there is percolation if the probability that there is an infinite path is positive.

It has been shown independently in  and  that if the graph is undirected then there is percolation even for $`m=4`$. In traditional percolation theory, when there is percolation then typically (unless the probability of blocking is at a “critical point”) the probability that there is percolation to a distance $`n`$ but no percolation to infinity is exponentially small in $`n`$. This is the case also in the papers cited above, but it is not true for the directed percolation problem we are facing.

###### Theorem 1.

If there is percolation from the origin to infinity with positive probability then the probability of percolating from the origin to distance $`n`$ but not to infinity is at least $`Cn^\alpha `$ for some constants $`C,\alpha >0`$ depending only on $`m`$.

## 2. The proof

Let $`b_m=(0,1,2,\mathrm{},m1)`$ be called the basic color sequence of length $`m`$: it is simply the list of all different colors. Let $`b_m^{}`$ be the reverse of $`b_m`$, i.e. $`b_m^{}(i)=b_m(mi1)`$. Let $`_{n,k}`$ be the event that for all $`i[0,k1]`$, $`j[0,m1]`$ we have

$$Y(n+im+j1)=b_m^{}(j),$$

i.e. starting with the index $`n1`$, the sequence $`Y`$ has $`k`$ consecutive repetitions of $`b_m^{}`$. We say that $`i`$ is an index of the occurrence of $`b_m`$ in the sequence $`X`$ if $`X(i+j)=b_m(j)`$ for $`j[0,m1]`$. For $`i>0`$, let $`\tau _i`$ be the $`i`$-th index of occurrence of $`b_m`$ in $`X`$. Let $`_{n,k}`$ be the event that for all $`i[1,n]`$ we have

$$\tau _{i+1}m\tau _ik1,$$
(1)
and also $`\tau _1k1`$.

###### Lemma 2.1.

If for some integer $`n>0`$ both $`_{n,k}`$ and $`_{n,k}`$ hold then there is no white infinite directed path.

###### Proof.

Let us assume that there is a white infinite path. Since there are $`n`$ consecutive copies of $`b_m`$ in $`X`$, the $`i`$-th one starting at index $`\tau _i`$, for each $`p[1,n]`$ there must be a vertical step in the path with an $`x`$ projection in $`[\tau _p,\tau _p+m1]`$. Therefore the path ascends to the segment $`[0,\tau _n+m1]\times \{n1\}`$, before its $`x`$ projection reaches $`t_n`$.

For each $`1pn`$, $`0q<k`$ there is a diagonally descending sequence of $`m`$ colored points

$$\{(\tau _p+j,n+(q+1)mj2):j[0,m1]\}.$$

For a fixed $`p`$ there are $`k`$ such diagonal barriers stacked above each other, forming an impenetrable column of height $`km`$. The path would have to ascend between two of these columns, say between column $`i`$ and $`i+1`$. The distance of two consecutive columns from each other is

$$\tau _{i+1}m\tau _ik1.$$

Since there are $`k`$ consecutive copies of $`b_m^{}`$ in $`Y`$ starting at index $`n1`$, for each $`q[0,k1]`$ there must be a horizontal step in the path with a height in $`n+qm1+[0,m1]`$. Since the distance of the two columns is at most $`k1`$, it is not possible for the path to pass between the two columns. (The same holds for the space before the first column.) ∎

###### Lemma 2.2.

There is a constant $`\alpha `$ such that for all $`s>0`$, there is a $`k`$ with

$$Prob(_{n,k})m^mn^{\alpha (s+1)},Prob(_{n,k})1n^s.$$

###### Proof.

With $`p_1=m^m`$, we have $`Prob(_{n,k})=p_1^k`$. Let us estimate the probability of $`_{n,k}`$. The probability that $`\tau _1>k1`$ is upperbounded by the probability that a copy of $`b_m`$ does not begin at $`i`$ in $`X`$ for any $`i`$ in $`\{0,m,\mathrm{},(k1)/mm\}`$, which is

$$(1p_1)^{(k1)/m+1}(1p_1)^{k/m}<e^{p_1k/m}.$$

The same estimate holds for the probability of (1) assuming that the similar conditions for smaller $`i`$ have already been satisfied. Hence

$$Prob(_{n,k})>1ne^{p_1k/m}.$$

Let us choose

$$k=(s+1)(m\mathrm{log}n)/p_1$$
(2)
for some $`s>0`$, then $`Prob(_{n,k})>1n^s`$, while

$$Prob(_{n,k})p_1n^{\frac{m\mathrm{log}p_1}{p_1}(s+1)}.$$

###### Proof of Theorem 1.

Assume that there is an infinite path with some positive probability $`P_0`$. Choose $`s`$ such that $`n^s<0.5P_0`$ and choose $`k`$ as a function of $`s`$ as in (2). Then the probability that $`_{n,k}`$ holds and there is an infinite path is at least $`0.5P_0`$. Let $`𝒢_n`$ be the event that there is a path leaving the rectangle $`[0,\tau _n]\times [0,n1]`$. Then $`Prob(_{n,k}𝒢_n)0.5P_0`$. Since $`_{n,k}𝒢_n`$ is independent of $`_{n,k}`$, we have

$$Prob(_{n,k}_{n,k}𝒢_n)0.5P_0m^mn^{\alpha (s+1)}.$$

The author is thankful to János Komlós and Endre Szemerédi for telling him about the problem.