text:
# Magnetoresistance of a 2D electron gas caused by electron interactions in the transition from the diffusive to the ballistic regime.

## Abstract

On a high-mobility 2D electron gas we have observed, in strong magnetic fields $`(\omega _c\tau >1)`$, a parabolic negative magnetoresistance caused by electron-electron interactions in the regime of $`k_BT\tau /\mathrm{}1`$, which is the transition from the diffusive to the ballistic regime. From the temperature dependence of this magnetoresistance the interaction correction to the conductivity $`\delta \sigma _{xx}^{ee}(T)`$ is obtained in the situation of a long-range fluctuation potential and strong magnetic field. The results are compared with predictions of the new theory of interaction-induced magnetoresistance.

Electron-electron interaction (EEI) corrections to the Drude conductivity $`\sigma _0`$ of 2D systems have been intensively studied over two decades. These studies were based on the theory of interactions in the diffusive regime, $`k_BT\tau /\mathrm{}<1`$ Altshuler . Physically this condition implies that the effective interaction time, $`\mathrm{}/k_BT`$, is larger than the momentum relaxation time $`\tau `$ and therefore the two interacting electrons experience scattering by many impurities. In the ballistic regime, $`k_BT\tau /\mathrm{}>1`$, electrons interact when scattered by a single impurity. A theory of the interaction correction for such a case was only recently developed Aleiner1 , and there have already been several experimental attempts to apply it to the conductance of high-mobility (large $`\tau `$) semiconductor structures Yuri2 ; Shashkin ; Pudalov ; Vitkalov ; Coleridge ; Noh . An essential feature of this theory is that the impurities are treated as point-like scatterers – the condition which is not satisfied in structures where the impurities are separated from the 2D channel by an undoped spacer (unless the spacer is thick enough for the background impurities to dominate the scattering). There is then a question of how the interaction correction in the ballistic regime manifests itself in a smooth fluctuation potential.

Introducing a long-range scattering potential is expected to suppress the interaction correction in the ballistic regime considered in Aleiner1 . This correction is caused by electron back-scattering, but in the case of a smooth potential the backscattering is significantly reduced. However, as shown in Rudin ; Gornyi , applying a strong magnetic field increases the probability of an electron to return back and restores the interaction correction.

Experimentally, the interaction correction $`\delta \sigma _{xx}^{ee}(T)`$ is usually obtained from the temperature dependence of the conductance, where it has to be separated from the interference correction $`\delta \sigma _{xx}^{WL}(T)`$ caused by the weak localisation (WL) effect Lee , as well as a possible classical contribution from phonon scattering. It has been shown however that in the diffusive regime there is an elegant method of detecting $`\delta \sigma _{xx}^{ee}(T)`$ from the perpendicular magnetoresistance Paalanen . The interaction correction gives rise to the parabolic negative magnetoresistance (NMR), expressed by the following relation at $`\omega _c\tau >1`$:

$$\rho _{xx}=\frac{1}{\sigma _0}+\frac{1}{\sigma _0^2}(\mu ^2B^2)\delta \sigma _{xx}^{ee}(T),$$
(1)
where $`\mu `$ is the electron mobility. This relation is derived by converting the conductivity tensor into the resistivity tensor and using the fact that in the diffusive regime the Hall conductivity is not affected by interactions: $`\delta \sigma _{xy}^{ee}(T)=0`$. Another essential specific of the diffusive regime used in its derivation is that $`\delta \sigma _{xx}^{ee}(T)`$ and $`\delta \sigma _{xy}^{ee}(T)`$ are not changed with strong magnetic field applied Houghton .

There have been several experiments where this method of extracting the interaction correction was used Paalanen ; Minkov ; Choi . However, apart from the experiment Minkov  on low-mobility structures, these experiments were in fact performed not in the diffusive but ballistic regime, where Eq. 1 is not justified. The extracted quantity $`\delta \sigma _{xx}^{ee}(T)`$ should be neither a logarithmic correction in the diffusive regime Altshuler , nor a linear correction in the ballistic regime Aleiner1  derived for classically small magnetic fields.

The new theory Gornyi  considers the interactions in the ballistic and intermediate regimes in strong fields $`\omega _c\tau >1`$. It shows that interactions in systems with a long-range potential will also produce a parabolic magnetoresistance described by Eq. 1. According to this theory, strong magnetic field not only restores the interaction correction, but also provides the required condition for Eq. 1: $`\delta \sigma _{xy}^{ee}(T)/\sigma _{xy}<\delta \sigma _{xx}^{ee}(T)/\sigma _{xx}`$. The theory gives a distinct prediction for the magnitude of the magnetoresistance $`\delta \sigma _{xx}^{ee}(T)`$ in Eq. 1. The aim of this work is to study the magnetoresistance caused by electron-electron interactions in the intermediate regime, in a structure with long-range fluctuation potential, and to compare the results with the prediction of Gornyi .

The sample is a standard modulation doped n-GaAs heterostructure. The doped layer is separated from the conducting channel by an undoped spacer $`d=20`$ nm. The wafer has been mesa etched to a Hall bar pattern. Measurements have been performed by a standard 4-terminal method with a current of 0.4 - 2 nA. The electron mobility changes in the range $`0.42\times 10^55.5\times 10^5`$cm2/Vs with increasing carrier density, which is lower than that of the samples in Paalanen . The range of electron densities is from $`0.46\times 10^{11}`$ cm-2 to $`2\times 10^{11}`$ cm-2. The parameter $`k_BT\tau /\mathrm{}`$ in our experiment is varied from 0.04 to 3.8. (This value is 3.3– 33 in Paalanen , 1.7–5.6 in Choi  and 0.004–0.18 in Minkov ).

The magnetoresistance of the 2DEG with $`n=6.8\times 10^{10}`$cm-2 is shown in Fig. 1(a). To analyse the data in terms of theory Gornyi , we have to prove first that the experimental conditions satisfy the theoretical approximations. Firstly, in the measured electron density range the $`k_Fd`$ value varies from 1.2 to 2.2, which proves that the fluctuation potential, with the correlation length $`d`$, is indeed long-range. This is further supported by an estimation of the ratio of the momentum relaxation time to the quantum time found from the magnitude of the Shubnikov–de Haas oscillations: $`\tau /\tau _q`$ varies from 24 at the highest to 4 at the lowest electron density. Secondly, the magnetoresistance is analysed in the magnetic field range $`\omega _c\tau >1`$. Thirdly, the gated structure is suitable for the study of interaction effects as the gate is separated from the 2DEG by 70 nm, while the separation between electrons is in the range from 13 nm to 26 nm and therefore interactions are not screened by the metallic gate.

One can see in Fig. 1 that the negative magnetoresistance exhibits a sharp change in small fields, followed by a parabolic dependence in higher fields. The sharp change is caused by the WL effect which is suppressed by magnetic field. We will analyse the parabolic magnetoresistance in the range of fields well above the ‘transport’ magnetic field $`B_{tr}=\mathrm{}/4De\tau 0.013`$T, to make sure that the WL contribution to the negative magnetoresistance is negligible. On the other hand, the magnetic field should not be too large, in order to prevent the development of the positive magnetoresistance caused by the Zeeman effect on the interaction correction Altshuler . This condition is also satisfied in our experiments where $`g\mu _BB/k_BT1`$. The contribution of this effect to the measured mangetoresistance can be estimated, from the theory of interactions in the diffusive regime Altshuler  and in the ballistic regime with point scatterers Aleiner1 , as being less than $`1\%`$ in the entire range of studied fields and thus can be neglected.

After the initial rapid change, the magnetoresistance develops a parabolic dependence, Fig. 1 (b). It is temperature dependent and shows immediately the qualitative features of the model Gornyi , Fig. 1(c). A flat region can be seen at small fields, which is a clear indication that the long-range potential at small field suppresses the interaction correction which is then restored by larger fields. In accordance with the theory, the flat region is better seen at higher temperatures. Another feature is seen in Fig. 1(d) - with increasing temperature the magnetoresistance becomes positive (in the model this is the result of the Hartree term becoming dominant over the exchange term at higher temperatures).

Before proceeding with the analysis of the strength of the magnetoresistance at different temperatures and extracting the temperature dependent EEI correction, we would like to note that there is also a classical contribution to the parabolic negative magnetoresistance in high-mobility structures which can also be seen in strong fields Dmitriev ; Polyakov . This NMR originates from the fact that the magnetic field makes it possible for electrons to ‘cycle’ around impurities and become localised, and hence do not contribute to the conductance. The magnitude of this effect is very different for short-range and long-range fluctuation potentials. In the case of short range scatterers the classical magnetoresistance has the form Baskin ; Bobylev ; Dmitriev :

$$\rho _{xx}=\rho _0\frac{1p}{1+p^2/(\omega _c\tau )^2},$$
(2)
where $`p=exp(2\pi /\omega _c\tau )`$ is the fraction of cycling electrons, and $`\omega _c`$ is the cyclotron frequency. Fig. 2 shows the comparison of the prediction by Dmitriev  with experimental results at the intermediate electron densities $`n=5.7\times 10^{10}`$ cm-2 and $`n=9.0\times 10^{10}`$ cm-2. At high densities the magnitude of the classical effect expected for short-range scatterers is much stronger than the measured magnetoresistance. Furthermore, the direction of its change with varying density (different $`\tau `$) would be opposite to the experimental situation. This proves again that in the studied samples we are dealing with long-range rather than short-range scattering. It is the presence of a smooth potential which significantly decreases the classical NMR Polyakov , since the ‘cycling’ electron trajectories will now become the trajectories ‘wandering’ in the potential landscape. For the moment we assume that the contribution of the classical magnetoresistance is negligible and attribute all magnetorestance to the interaction correction $`\delta \sigma _{xx}^{ee}(T)`$ (later we will come back to this question).

We plot the resistivity as a function of $`B^2`$ and from the slope of the straight line obtain $`\delta \sigma _{xx}^{ee}(T)`$. Fig. 3 shows the temperature dependence of $`\delta \sigma _{xx}^{ee}`$ at different electron densities, where experimental points concentrate around one curve. This curve becomes rather close to the interaction correction in the exchange channel Gornyi  if one makes a vertical shift of the theoretical curve by $`\mathrm{\Delta }\sigma =0.07e^2/h`$. We want to emphasise that apart from this small shift there are no adjustable parameters involved in the analysis. It is interesting to note that the interaction correction $`\delta \sigma _{xx}^{ee}(T)`$ found from the temperature dependence of the conductance at $`B=0`$ is defined with an accuracy of a constant contributing to a renormalised value of the momentum relaxation time. In the method of quadratic magnetoresistance, this constant does not contribute to the magnetoresistance and no shift is allowed in comparing the results with the theory. We will discuss below a possible physical origin of this additional, temperature independent contribution to the quadratic magnetoresistance.

The theoretical curve plotted in Fig. 3 is given by the following expression that describes electron interactions in the exchange channel:

$$\frac{\delta \rho _{xx}(B)}{\rho _0}=\frac{(\omega _c\tau )^2}{\pi k_Fl}G_F(T\tau ),$$
(3)
where $`k_F`$ is the Fermi wave number and $`l`$ is the mean free path. The function $`G_F(x)`$ has the asymptotes $`G_F(x1)lnx+const`$ and $`G_F(x1)(c_0/2)x^{1/2}`$, with $`c_00.276`$.

Measurements at different electron densities have enabled us to cover a broad range of the parameter $`T\tau `$ and detect the specific features of this dependence. It shows a logarithmic behaviour at small $`T\tau `$, followed by a rapid disappearance of the interaction correction at higher temperatures. The results show good agreement with the expected ‘saturation’ at $`T\tau 0.5`$, which means that the turning effect of magnetic field is reduced with increasing temperature and electrons have less chance to return back in the long-range potential.

Note that if one naively compared the obtained dependence $`\sigma _{xx}^{ee}(T)`$ with the one calculated in Aleiner1  (for point-like scatterers in zero field), a striking difference would be seen, Fig. 3 (dotted line). The latter has a linear asymptote at $`T\tau 1`$ and does not show any saturation.

It is important to emphasise that the comparison was made with the contribution from the exchange channel only, however it is known that there is another (Hartree) term in interactions controlled by the Fermi liquid interaction parameter $`F_0^\sigma `$. Comparing the total interaction correction (exchange plus Hartree Gornyi ) with the experiment shows that the Hartree contribution is much smaller than the exchange contribution. It can be seen in Fig. 3 that within the experimental error the magnitude of the parameter $`F_0^\sigma `$ in our case cannot be larger than 0.1–0.2. The value of $`F_0^\sigma `$ depends on the parameter $`r_s=1/(\pi n)^{1/2}a_B`$, which is the ratio of the Coulomb and kinetic energy of electrons. The value of the Fermi liquid parameter is only known for $`r_s<1`$: $`F_0^\sigma =\frac{1}{2\pi }\frac{r_s}{\sqrt{2r_s^2}}\mathrm{ln}(\frac{\sqrt{2}+\sqrt{2r_s^2}}{\sqrt{2}\sqrt{2r_s^2}})`$ Aleiner1 . We plot this value in Fig. 4, together with results of two recent experiments where it was determined, in different 2D structures, at large $`r_s`$. (The latter results were obtained on systems with little effect of the smooth potential, using analysis based on theory Aleiner1 .) The overall trend on Fig. 4 indicates that in our range of $`r_s=1.22.6`$ it is reasonable to expect the value of $`F_0^\sigma `$ to be $`0.15`$.

Let us now return to the observed shift in Fig. 3. We believe that it is caused by a contribution from the classical NMR Polyakov , known to be $`T`$-independent, which we completely ignored earlier because the long-range potential suppresses the mechanism described by Eq. 2. Quantitatively, the parabolic classical negative magnetoresistance depends on the ratio of the short- and long-range scattering times, $`\tau _s/\tau _L`$. For the case of $`\tau _L>\tau _s`$, the magnetoresistance is given by the following relation Polyakov :

$$\delta \rho _{xx}/\rho _0=\omega _c^2/\omega _0^2,$$
(4)
where $`\omega _0=(2\pi n_s)^{1/2}v_F(2l_s/l_L)^{1/4}`$, with $`n_s`$ the concentration of the strong scatterers, $`l_s`$ and $`l_L`$ the mean free paths for strong and smooth potential scatterers, respectively. Increasing the contribution of long-range scattering (increasing this ratio) significantly suppresses the classical magnetoresistance. We can estimate the value of $`\tau _L`$ using the expression for $`\tau `$ in the case of remote donor scattering Davies : $`\frac{1}{\tau }=\frac{(\pi \mathrm{})n_s^{2D}}{8m(k_Fd)^3}`$, where $`n_s^{2D}`$ is the concentration of scattering donors which is approximately equal to the electron density. For short-range scatterers the estimation of $`\tau `$ can be done using the relation $`\frac{1}{\tau }=\frac{ha_Bn_s^{3D}}{m(k_Fa_B)^3}`$ for background impurity scattering, assuming that $`n_s^{3D}3\times 10^{14}`$ cm-3. For $`n=9.0\times 10^{10}`$cm-2 these estimations give $`\tau _L6`$ ps and $`\tau _s10`$ ps, which are close to the momentum relaxation time $`\tau =6.8`$ ps. An estimate using Eq. 4 shows that this effect can indeed account for the experimentally observed shift. A more accurate comparison is, however, complicated due to uncertainties in the values of $`\tau _s`$ and $`\tau _L`$ and the fact that Eq. 4 is, strictly speaking, only valid for $`\tau _L\tau _s`$.

In conclusion, we have studied the magnetoresistance of a high-mobility 2D electron gas in a GaAs/AlGaAs heterostructure where the electron scattering is determined by a long-range fluctuation potential. In classically strong magnetic fields we have observed negative magnetoresistance, which is parabolic and temperature dependent. We have shown that this magnetoresistance originates from the electron-electron correction to the conductance. We have extracted this correction and demonstrated that it is well described by the recent theory of interactions in the regime which is intermediate between the diffusive and ballistic regimes.

We are grateful to I. V. Gornyi and A. D. Mirlin for stimulating discussions, and I. L. Aleiner for reading the manuscript and valuable comments. L.L. and Y.Y.P. are grateful to the ORS scheme for the financial support.

Figure captions:
Fig. 1. (a) Longitudinal resistivity versus magnetic field for electron density $`n=6.8\times 10^{10}`$cm-2 at different temperatures: $`T=0.2`$, 0.8, 1.2 K. (b) The same data presented as a function of $`B^2`$. (c) Zoomed-in region of $`\rho _{xx}(B)`$ from (a) at $`T=1.2`$ K, showing a flat region at small fields. (d) $`\rho _{xx}`$ versus $`B^2`$ for another density, $`n=9\times 10^{10}`$cm-2, at $`T=0.4`$, 1.2, 4.2 K, showing a transition from negative to positive magnetoresistance.
Fig. 2. Solid lines: Negative magnetoresistance for the momentum relaxation time $`\tau =6.8`$ ps and $`\tau =2.3`$ ps. Dashed lines: NMR introduced by the classical effect Dmitriev , with $`\tau =6.8`$ and 2.3 ps.
Fig. 3. Conductivity correction due to interactions, obtained experimentally at different electron densities $`n=0.46,0.57,0.68,0.90,1.2,1.4,2.0\times 10^{11}`$ cm-2 (different symbols for different densities, circles correspond to the lowest density). Solid line - theoretical prediction for the correction due to exchange interaction, shifted by $`0.07e^2/h`$. Dashed line - theory for the total correction which includes exchange (Fock) and Hartree interactions with Fermi-liquid parameter $`F_0^\sigma =0.15`$. Dotted line shows the result of the interaction theory for point-like scatterers in the transition from the diffusive to the ballistic regime (exchange term at zero magnetic field). Inset: the same results presented in the logarithmic scales.
Fig. 4. The dependence of the Fermi liquid parameter on $`r_s`$. Dashed box indicates an approximate range of $`r_s`$ and $`F_0^\sigma `$ for the structures in this work. Open squares: results from Ref.Yuri2  for hole gas in p-GaAs heterostructures. Solid squares: results from Ref.Shashkin  for electron gas in Si MOSFETs. Solid line is the theoretical curve for small $`r_s`$ Aleiner1 .


format: text

text:
# Damage spreading in the Bak-Sneppen model: Sensitivity to the initial conditions and equilibration dynamics

## I Introduction

The concept of self-organized criticality was introduced to describe the tendency of large dynamical systems to organize themselves in an out of equilibrium state exhibiting spatio-temporal complexity. Such complexity is reflected by the presence of correlations between events separated in space and time over a wide range of length and time scales. The emergence of spatio-temporal complexity in dynamical systems driven by rules based on extremal principles is in the origin of the widespread occurrence in nature of fractal structures, noise with $`1/f`$ power spectrum and punctuated equilibrium.

The simplest mathematical model of a driven dynamical system, which exhibits self-organized criticality, was introduced by Bak and Sneppen as a toy model of biological evolution of an ecology of interacting species. In this model random numbers $`f_i`$ (fitness), uniformly distributed on the interval , are assigned to the sites of a $`d`$-dimensional lattice. At each time step the site with the smaller random number is located and new random numbers from a uniform distribution are assigned to that site and also to its first neighbors.

The above model system achieves a statistically stationary state in which the density of random numbers in the system vanishes for $`f<f_c`$ and is uniform above $`f_c`$, with $`f_c=0.667`$  in a chain and $`f_c=0.38`$ in a square lattice geometry . Once the stationary state is achieved the system exhibits punctuated equilibrium as real biology, characterized by intermittent co-evolutionary avalanches of all sizes. The complexity of this stationary regime can be revealed by the show up of spatio-temporal power-law distribution of events. For example, the distribution $`C(x)`$ of the distance $`x`$ between two subsequent extremal sites scales as $`C(x)x^\pi `$, with $`\pi =3.23`$ in one dimension . Further the temporal long-range correlations can be also observed in the distribution of first return times $`P(t)`$ which scales as $`t^{\tau _{first}}`$, with $`\tau _{first}=1.58`$ in $`1d`$ and $`\tau _{first}=1.28`$ in $`2d`$ .

The propagation of local perturbations in the self-organized state of Bak-Sneppen model was recently investigated using the damage spreading algorithm. By measuring how the difference between two initially close configurations evolves in time under the same noise, it was shown that an initial power-law divergence sets up before the distance saturates in a finite size-dependent plateau. This slow short-time power-law dynamics is characteristics of systems poised at critical states such as at the critical point of second-order phase transitions and the onset of chaos in non-linear dynamical systems.

However, usual critical phenomena exhibit two dynamical regimes characterized by distinct critical exponents. The first one is related to a short-time slow dynamics, which governs the power-law increase of the order parameter when its value in an initial non-equilibrium configuration is small and non-zero. In finite systems this power-law increase persists up to a characteristic time, after which the order parameter slowly decays towards its equilibrium value.

The same dynamical regimes have been identified in the onset of chaos of low-dimensional non-linear dynamical maps. The distance between two nearby orbits was shown to diverge following a power-law whose exponent is directly related to geometric exponents characterizing the extremal sets of the critical dynamical attractor. On the other hand, the long-time relaxation towards the dynamical attractor is governed by a new exponent, which seems to be related to the fractal dimension of the support of the dynamical attractor in phase-space.

Damage spreading studies of critical phenomena have also reported two dynamical regimes with distinct and independent exponents. Therefore, it is natural to expect this to be the general trend of systems poised at criticality as well as systems exhibiting self-organized criticality such as the Bak-Sneppen model. The fact that previous damage spreading works have not identified the long-time relaxation towards equilibrium in the BS model just reflects the insensitivity of the proposed measure to capture the long-time dynamics.

In this work we introduce a new Hamming distance measure that is sensitive to the long-time dynamics of the Bak-Sneppen model. Following its temporal evolution we are going to show that, besides displaying a power-law short-time dynamics with the same exponent previously reported, it also exhibits a long-time relaxation dynamics. We will report the scaling properties of this regime and discuss its possible relation with the geometric properties of the dynamical attractor in phase space.

## II Model system and numerical simulation

In our simulations we implement the Bak-Sneppen algorithm in linear chains with periodic boundary conditions. Random numbers $`f_i`$ (fitness) are initially assigned to all sites from a uniform distribution in the interval $`[0,1]`$. The system evolution is based on the standard re-assignment of the extremal site, i.e., the site with the minimum fitness and of its first neighbors. In our simulations we worked with chains of up to $`N=1000`$ sites and left $`t_{trans}=20N`$ collective time steps for the system to achieve the statistically stationary state. Our collective time unit corresponds to $`N`$ elementary time steps so that each site is going to be updated only once in average during a collective time step.

After achieving the statistically stationary state we produce a copy of the system’s configuration. A small damage is introduced in the copy by interchanging the position of the site with the minimum fitness with a randomly chosen one. After this, we follow the temporal evolution of both configurations using always the same set of random numbers to update both replicas. After a characteristic time $`\tau `$, both configurations will be composed by the same sequence of random numbers just shifted by a random distance. Once periodic boundary conditions are used, these two configurations are indeed indistinguishable and should be considered as identical. In figure 1 we report the finite size scaling of the average time needed for two initially random configurations to become identical when updated following the Bak-Sneppen dynamics and under the same set of random numbers. This characteristic time, measured in units of collective time steps, scales as $`\tau N^\varphi `$, with $`\varphi =1.46(4)`$. We would like to stress that, within the error bars, this is the same scaling exponent obtained by the size dependence of the average self-organization time to reach the critical attractor \[$`D1=1.43(1)`$\] and is independent of the initial state of the system.

In order to accomplish our task, we introduce a new Hamming distance $`D(t)`$ defined at each time $`t`$ as the smallest among the $`N`$ possible values of $`D_j(t)`$ ($`j=1,2,\mathrm{}N`$) defined as

$$D_j(t)=\frac{1}{N}\underset{i=1}{\overset{N}{}}|f_i^1f_{i+j}^2|.$$
(1)
where $`\mathrm{}`$ represents configurational average over distinct runs. Those runs at which the Hamming distance vanishes are not considered in the averaging process and therefore our simulations are limited in size by the characteristic time $`\tau `$. Notice that $`D_0(t)`$ corresponds to the measure used in previous damage spreading studies of Bak-Sneppen model. After the initial power-law growth $`D_0(t)t^\alpha `$ with $`\alpha =0.32`$, it saturates at a constant value which can be shown to be $`1/3`$ of the width of the fitness distribution in equilibrium. This saturation is due to the fact that by just measuring $`D_0(t)`$ one cannot identify that both replicas have converged to the same random sequence configuration.

The procedure proposed here is to measure the Hamming distance as the minimum distance between the replicas taking in account all possible shifts between them. In respect to the short-time dynamics it is equivalent to measuring $`D_0(t)`$ itself and shall give the same power-law exponent for the initial Hamming distance growth. However, after achieving a maximum, it must decrease to vanishing small values as both replicas converge to the same sequence, capturing then the long-time dynamics.

In figure 2 we show our results for $`D(t)`$ for chains up to $`N=1000`$ sites. As larger chain sizes are considered the initial power-law regime extends for longer periods. The dashed line corresponds to the power law $`D(t)t^\alpha `$ with $`\alpha =0.32`$ in full agreement with the value previously reported in the literature. After a size dependent characteristic time, it reaches a maximum and starts to monotonically decrease. However, much longer runs than these presently reported would be needed to precisely estimate the scaling behavior of the Hamming distance decay.

In order to investigate in detail the long-time dynamical regime, it is computationally more efficient to use a slightly different approach. We start with two uncorrelated equilibrium configurations and follow the time evolution of the distance between them when the same set of random numbers is used. In this way, the initial growth is absent and the relaxation regime can be observed using shorter runs. Usually the critical power-law relaxation dynamics takes place after some transient time. This typical behavior means that the Hamming distance at criticality satisfies the non-linear differential equation

$$\frac{dD(t)}{dt}=\lambda D^q$$
(2)
whose solution is

$$D(t)=D(0)/[1+\lambda (q1)t]^{1/(q1)}.$$
(3)
and reproduces the above trend. In numerical simulations, the critical relaxation dynamics is always followed by an exponential relaxation due to finite size effects. A pure exponential relaxation implies that $`D(t)`$ must asymptotically obey a linear differential equation. Therefore in order to consider simultaneously the critical and exponential relaxation we consider $`D(t)`$ to obey a more general non-linear differential equation such as

$$\frac{dD(t)}{dt}=(\lambda \mu )D^q\mu D.$$
(4)
In the above equation, $`\mu `$ represents the finite size correction and vanishes in the thermodynamic limit. The general solution of the above differential equation is in the form

$$D(t)=D(0)/[1(\gamma /\mu )+(\gamma /\mu )e^{(q1)\mu t}]^{1/(q1)}.$$
(5)
where $`\gamma =\mu (\mu \lambda )/D_0^{(1q)}`$. It is worth to mention that similar equations have also been used to describe experimental data on the re-association in folded proteins, quantitative linguistics and fluxes of cosmic rays  within the framework of nonextensive thermostatistics . Indeed, the above form is expected to generally describe the sensitivity to initial conditions of systems at the vicinity of critical states.

In Figure 3 we report our results for $`D(t)`$ starting with two uncorrelated equilibrium configurations for several chain sizes. The average initial distance equals to $`D(0)=0.111`$ which corresponds to $`1/3`$ of the width of the interval $`[f_c,1]`$ where $`f_c=0.667`$ for the $`1d`$ Bak-Sneppen model. The coefficient $`\mu `$ determines characteristic exponential decay which sets up in the long-time regime as can be seen in the inset. The solid lines correspond to best fits to the form of Eq. (5). We recall that the fitting parameters govern the typical behavior at distinct time regimes. In the table we summarize the parameters obtained. Notice that $`\mu `$ decreases monotonically when chain size is increased, once the exponential relaxation is postponed for large chains. Also $`\gamma `$ decreases with $`N`$ so that the transient period diverges in the thermodynamic limit. However, $`\mu \gamma `$ and therefore the pre-asymptotic decay can be observed just at intermediate times satisfying $`1/\gamma t1/\mu `$.

In figure 4 we show the size dependence of $`\mu `$ and $`\gamma `$. Once $`1/\mu `$ determines the time scale needed to finite size effects become relevant, it shall present the same scaling behavior as the average time $`\tau `$ for both replicas to become identical. We found $`1/\mu N^{1.49(4)}`$, which is consistent with the above conjecture. On the other hand, $`1/\gamma `$ defines the time scale when the crossover between the transient and the pre-asymptotic regime takes place. It is governed by a distinct exponent $`1/\gamma N^\kappa `$, with $`\kappa =1.11(3)`$ and is directly related to the average time needed to update all active sites in the initial avalanche.

For the system sizes simulated, the intermediate regime satisfying the condition $`1/\gamma t1/\mu `$ is achieved only in a narrow time interval which grows slowly as the chain size is increased. This feature introduces a large uncertainty in the estimates of the non-linear exponent $`q`$, particularly the ones from small chain data. The best fitting values for $`q`$ were found to be much larger than one, specially the more confident estimate from chain size $`N=1000`$. This may indicate a very slow logarithmic pre-asymptotic decay of the Hamming distance.

At this point it is worth to explore the similarities between the damage spreading analysis of the dynamical properties of extended critical systems and recent results concerning the sensitivity to initial conditions of low-dimensional dissipative maps poised at criticality. More precisely, the numerical analysis of such systems (i.e., logistic map, logistic-like maps, circular maps, asymmetric logistic maps, single-site map, Henon map) has shown that, at critical points such as the chaos threshold, tangent bifurcations etc, the standard type of the sensitivity to the initial conditions given by the sensitivity function

$$\xi (t)\underset{\mathrm{\Delta }x(0)0}{lim}\frac{\mathrm{\Delta }x(t)}{\mathrm{\Delta }x(0)}=\mathrm{exp}\left(\lambda _1t\right)$$
(6)
(where $`\lambda _1`$ is the standard Lyapunov exponent) can be replaced by a power-law type of sensitivity

$$\xi (t)=\left[1+(1q)\lambda _qt\right]^{1/(1q)},$$
(7)
which indicates a weak sensitivity to the initial conditions. Here, $`\lambda _q`$ is the generalized Lyapunov exponent, and satisfies the generalized Pesin equality $`K_q=\lambda _q`$ if $`\lambda _q>0`$ and $`K_q=0`$ otherwise, where $`K_q`$ is the generalized Kolmogorov-Sinai entropy .

Now we are in a position to discuss the connection between this scenario and the BS model. Indeed, the BS model belongs to another class of dynamical systems, namely high-dimensional dissipative systems, for which this scenario is also needed to be tested. As a matter of fact, the BS model evolves spontaneously (without a tuning parameter) towards a critical state, whereas a one-dimensional map approaches to a critical point (e.g. chaos threshold point) with a fine tuning of a map parameter. Moreover, in both cases, the standard Lyapunov exponent would be zero. Therefore, the properties of the sensitivity to the initial conditions of the BS model are expected to be similar to those of the low-dimensional dissipative maps at the edge of chaos. Sensitivity to the initial conditions, which is reflected by the sensitivity function for the low-dimensional dissipative maps, can be conveyed by the damage spreading (and Hamming distance) for the BS model. It is clear from ref. and from our results using a different definition of the Hamming distance that a power-law spread of damage emerges for the BS model in the short-time dynamics regime. Then, one can easily determine the proper $`q`$ value of the BS model as $`q^{}2.1`$.

In what follows we discuss the picture seen from the analysis of long-time dynamics of the BS model. Firstly, we should recall the results of low-dimensional dissipative maps. Recently, it has been shown  that, for the long-time dynamics of such dynamical systems, a new class of $`q`$ values emerges. This new class appeared to be related to the equilibration (or long-time dynamics). It is seen from the results of  that, whenever the critical attractor of the dynamical system has a fractal support (or in other words, the critical attractor is not dense), the system (e.g., logistic-like maps) exhibits a power-law long-time relaxation with a $`q`$ value greater than unity. On the other hand, if the critical attractor of the dynamical system (e.g., circular maps) is dense, then a slow logarithmic decay has been observed. When we look at the BS model, the situation is similar to that of the circular maps in the sense that the attractor in phase space is fully occupied in the interval $`[f_c,1]`$, which makes the critical attractor dense. The observed slow (logarithmic) pre-asymptotic dynamics seems to indicate that it shall take place in general system models whose critical dynamical attractor is dense.

## III summary and conclusions

In summary, we have applied the damage spreading algorithm to study the dynamics of the Bak-Sneppen model for biological evolution which is one of the prototypes models exhibiting self-organized criticality. By defining a proper Hamming distance between two system’s configurations we investigated both the short-time and long-time critical dynamics. A power-law short-time dynamics was found in agreement with previous reports. In the relaxation regime, the Hamming distance was shown to follow the solution of a generalized non-linear equation, first proposed in the context of nonextensive thermostatistics, which includes in a single expression the transient, pre-asymptotic and asymptotic exponential relaxation regimes. We found that the pre-asymptotic regime is very close to a slow logarithmic decay in a close relationship with the relaxation dynamics of non-linear dissipative maps at the onset of chaos with a dense dynamical critical attractor. It would be interesting to investigate the relaxation dynamics of self-organized critical systems with a fractal dynamical attractor in phase-space. In this case, based on previous studies of the critical dynamics of low-dimensional deterministic non-linear maps at the onset of chaos, one would expect a power-law relaxation. This point will be the subject of a future contribution.

## IV acknowledgements

This work has been supported by the Turkish Academy of Sciences, in the framework of the Young Scientist Award Program (UT/TUBA-GEBIP/2001-2-20). M.L.L. would very much appreciate to acknowledge the hospitality of the Physics Department of Ege University where this work was done and thanks TUBITAK for financially supporting his visit via NATO-D Program. Partial financial support from the Brazilian research agencies CNPq, CAPES and FAPEAL is also acknowledged.

## V Figure and Table Captions

Figure 1 - The finite size scaling of the average time for two initially uncorrelated random sequences to become identical when updated following the Bak-Sneppen dynamical rules using the same set of random numbers. The solid line corresponds to a power-law fit $`\tau N^\varphi `$ with $`\varphi =1.46(4)`$. The number of experiments used in calculations for different chain sizes is 100.

Figure 2 - The time evolution of the Hamming distance between an equilibrium configuration and its slightly modified replica. After the initial power-law growth $`D(t)t^\alpha `$ with $`\alpha =0.32`$, it monotonically decreases as both configurations converge to the same sequence of random numbers. The number of experiments used in calculations varies from 100 up to 500 for different chain sizes.

Figure 3 - The time evolution of the Hamming distance between two uncorrelated equilibrium configurations. After the initial transient where the Hamming distance is close to its initial value $`D(0)=0.111`$, $`D(t)`$ exhibits a non-trivial decay followed by an exponential relaxation due to finite size effects. The solid lines correspond to best fits to the form of Eq. (5) that contains the above three regimes. The fitting parameters are reported in the table. The number of experiments used in calculations varies from 100 up to 500 for different chain sizes.

Figure 4 - The finite size dependence of the fitting parameters $`\mu `$ and $`\gamma `$. The solid lines are best fits to power-laws. We found that $`1/\mu N^{1.49(4)}`$ and therefore it is, within the error bar, proportional to the average time needed to both configurations coincide. On the other hand $`1/\gamma N^\kappa `$ with $`\kappa =1.11(3)`$ has, within the error bar, the same scaling behavior of the average time needed to update all sites in the initial avalanche.

Table Caption - The size dependence of the best fitting parameters. Notice that $`\mu `$ decreases as the chain size $`N`$ increases, as expected. The large values of $`q`$ indicate a possible slow logarithmic relaxation in the pre-asymptotic regime as discussed in the text.

Table

| $`N`$ | $`\gamma `$ | $`\mu `$ | $`q`$ |
| --- | --- | --- | --- |
| $`\mathrm{\hspace{0.33em}100}`$ | $`\mathrm{\hspace{0.33em}0.088}`$ | $`\mathrm{\hspace{0.33em}0.012}`$ | $`\mathrm{\hspace{0.33em}8.7}`$ |
| $`\mathrm{\hspace{0.33em}200}`$ | $`\mathrm{\hspace{0.33em}0.039}`$ | $`\mathrm{\hspace{0.33em}0.0043}`$ | $`\mathrm{\hspace{0.33em}8.7}`$ |
| $`\mathrm{\hspace{0.33em}500}`$ | $`\mathrm{\hspace{0.33em}0.013}`$ | $`\mathrm{\hspace{0.33em}0.00094}`$ | $`\mathrm{\hspace{0.33em}10.0}`$ |
| $`\mathrm{\hspace{0.33em}1000}`$ | $`\mathrm{\hspace{0.33em}0.0062}`$ | $`\mathrm{\hspace{0.33em}0.00034}`$ | $`\mathrm{\hspace{0.33em}14.1}`$ |


text:
# 1 Introduction

## 1 Introduction

Cosmic ray observations have determined that the spectrum of the highest energy cosmic rays extends beyond $`10^{20}`$ eV . Observations have also indicated that the highest energy spectrum is dominated by protons rather than photons . Above $`5\times 10^{19}`$ eV, protons can interact with cosmic background photons at the $`\mathrm{\Delta }`$-resonance generating pions. Above this energy, called the GZK (Greisen-Zatsepin-Kuzmin) cutoff , the proton energy loss length is near 50 Mpc, thus requiring semi-local sources to produce the observed flux. The lack of any such known sources has spawned a great deal of speculation as to the origin of these particles. A common class of models, called top-down scenarios, involve supermassive particles which decay or annihilate generating the highest energy cosmic rays .

The decay of superheavy particles has been studied in some detail . In particular, it has been demonstrated that a significant amount of the initial energy of such a particle can be emitted in the form of ultra-high energy supersymmetric particles . In most models, the lightest supersymmetric particle is a neutralino. This neutralino, weakly interacting and stable by virtue of $`R`$parity, can travel cosmological distances without absorption or scattering. In this paper, we discuss the prospects for observing ultra-high energy cosmic neutralinos in future very large area, satellite-borne air shower experiments.

## 2 Ultra-High Energy Fragmentation To Neutralinos

In the general framework of top-down scenarios, one has to consider the decay of super-heavy $`X`$ particles with a mass of the order of $`10^{21}`$ to $`10^{25}`$ eV, and a lifetime comparable to or longer than the age of the universe. Such a long lifetime can be ensured by “storing” the $`X`$particles in cosmological defects, which can survive into the present epoch . Alternatively, free $`X`$particles might be long–lived since their decay is suppressed, e.g. by (approximate) symmetries . For a review of different candidates, see . Such particles could be produced in the very early times of the universe, e.g. at the end of inflation . The typical decay modes of the $`X`$ particles are generally unknown and/or quite model dependent. However, if there is no additional “new physics” scale between $`M_{\mathrm{SUSY}}1`$ TeV and $`M_X`$This hypothesis, known as the “desert hypothesis”, is well motivated by the fact that the existence of new physics between the GUT scale and the SUSY breaking scale would destroy the very impressive feature of “natural” unification of the gauge couplings at $`M_{GUT}10^{16}`$ GeV occuring in the MSSM ., the $`X`$ particles should decay into $`N`$ “known” particles of the minimal supersymmetric standard model (MSSM), and usual particle physics allows us to study in detail the shower generated by the primary products of the initial $`X`$ decay.

A detailed computation of the spectra of stable particles (protons, photons, neutralino LSPs, electrons and neutrinos of the three species) obtained in such decay showers has been described in . We recall here that at the energies we are considering, it is necessary to take into account all the gauge couplings of the MSSM; indeed, at the scale of unification, they are all of the same strength, so that electroweak (and some Yukawa) interactions can be as relevant as the QCD ones. The perturbative part of the shower was computed by solving numerically the complete set of evolution equations  for the relevant fragmentation functions of the MSSM. We carefully modeled the decays of unstable particles with mass near $`M_{\mathrm{SUSY}}1`$ TeV, as well as the hadronization process at the GeV scale for light quarks and gluons. We found that the LSP flux depends only mildly on the spectrum of superparticles, as long as the LSP is a bino–like neutralino . Some sample spectra are shown in Fig. 1.The primary 10–body decay $`X5q5\stackrel{~}{q}`$ has been modeled using phase space only, i.e. ignoring any possible dependence of the matrix element on external momenta.

Here we have conservatively assumed that $`X`$ particles have an overdensity of $`10^5`$ in the vicinity of our galaxy, as expected  for $`X`$ particles that move freely under the influence of gravity.§§§The exact profile of the halo of $`X`$ particles does not affect our results as long as most UHECR events originate at distances well below one GZK interaction length. This minimizes the expected neutralino flux, since all scenarios are normalized by matching  the predicted proton spectrum at $`E10^{20}`$ eV to the highest energy cosmic ray observations. Data  indicate that most UHE events have protons (or heavier ions, which however cannot be produced in top–down models), not photons, as primary particles. Specifically, Yakutsk, Haverah Park and, most recently, AGASA observe more muons in the events than expected from photon primaries. The longitudinal development of the most energetic Fly’s Eye event indicates that it is not due to a photon. Finally, the highest energy AGASA events should show a North–South asymmetry if they were due to photons (which, at that energy, initiate an electromagnetic cascade already in the Earth’s magnetosphere), but no such asymmetry is observed. On the other hand, top–down models predict the photon flux at source to exceed the proton flux .

Following ref. we are thus compelled to assume that most UHE photons are absorbed somewhere between source and Earth. According to ref. this is in fact expected if most sources are at cosmological distances. At $`E10^{20}`$ eV the attenuation length for photons is at least one order of magnitude shorter than that for protons. This ratio is closer to 50, if the extragalactic $`\stackrel{}{B}`$ fields average at least $`10^9`$ G, as is widely expected. On the other hand, according to standard estimates  the density of radio photons in our galaxy is only about one order of magnitude higher than the average value in intergalactic space, leading to an attenuation length of $`10^{20}`$ eV photons of at least 100 kpc. In order to obtain near–complete absorption of such photons in our galaxy the actual density of radio photons would have to be at least ten times higher than the accepted value. This may not be impossible. The galactic (and extragalactic) radio background is reasonably well known only at frequencies above 100 MHz. However, we are most interested in photons with frequency of only a few MHz, which have the highest cross section for $`e^+e^{}`$ pair production with $`10^{20}`$ eV photons. This frequency band is difficult to observe on Earth, due to strong absorption in the ionosphere.

## 3 Signatures of Ultra-High Energy Neutralinos

Ultra–relativistic neutralinos interact with quarks by $`t`$channel $`Z`$ and $`W^\pm `$ exchange, as well as by the exchange of squarks in the $`s`$ or $`u`$channel. These interactions either directly yield an LSP, or produce a heavier neutralino or chargino which quickly decays to the lightest neutralino (except, perhaps, in the case of near–degenerate masses). Either interaction generates a shower which can be observed by air shower experiments.

The background for this signal consists of showers generated by ultra-high energy cosmic neutrinos. The neutrino interaction length becomes comparable to the radius of the earth around $`10^5`$ GeV. By $`10^9`$ GeV, only about one out of 1000 neutrinos passes through the Earth without interaction (see figure 2). A neutralino, however, depending on the choice of SUSY parameters, will have a different interaction cross section and, therefore, different absorption properties. The size of this cross section depends sensitively on the neutralino eigenstate, which in general is a composition of bino, wino and neutral higgsinos. A wino– or higgsino–like neutralino has couplings to $`W`$ and/or $`Z`$ bosons that resemble or even exceed those of neutrinos. In contrast, a bino–like neutralino has very small couplings to gauge boson, because its superpartner, the $`U(1)_Y`$ gauge boson, does not couple to other gauge bosons. The couplings of bino–like neutralinos to squarks are of full $`U(1)_Y`$ gauge strength, but squark searches at the Tevatron  tell us that first and second generation squarks must be at least three times heavier than $`W`$ bosons. Note also that models with radiative breaking of the electroweak gauge symmetry prefer the lightest neutralino to be bino–like in most of parameter space . Typical parameter choices therefore predict neutralino-nucleon cross sections one or two orders of magnitude smaller than neutrino-nucleon cross sections . With a significantly smaller cross section, very high energy cosmic neutralinos may travel through the Earth producing upgoing events at much higher energies than neutrinos. Upgoing showers with energy above 100 PeV or so would be a smoking gun for cosmic neutralinos.

Furthermore, by virtue of $`R`$parity, neutralinos will generate less energetic neutralinos in each interaction, thus not depleting their number. Tau neutrinos also display this property , but not as dramatically. The difference comes from the fact that high energy tau leptons lose energy in propagation whereas charginos decay quickly enough to lose very little energy in propagation. Also, phase space arguments indicate that a larger fraction of a decaying chargino’s energy goes into the resulting (massive) neutralinos than a decaying tau’s energy goes into the new (essentially massless) tau neutrino. Together, these effects indicate that tau regeneration is largely ineffective above about $`10^8`$ GeV. On the other hand, for even moderately smaller neutralino cross sections, the Earth can remain effectively transparent to cosmic neutralinos at much higher energies.

Our calculations of tau neutrino and neutralino regeneration in the Earth were done with a Monte Carlo simulation which, at each interaction, calculated the energy lost in the interaction and following propagation . Our treatment of $`\tau `$ propagation includes $`e^+e^{}`$ pair production, photonuclear interactions, bremsstrahlung and ionization energy losses. As stated earlier, any unstable superparticle produced in LSP interactions is too short–lived to lose energy prior to its decay. We estimate that each interaction, if necessary followed by superparticle decay, will reduce the energy of the LSP by slightly more than a factor of two; this effect is included in our treatment of LSP regeneration. Our code demonstrated the appearance of a ‘pile-up’ of outgoing particles at an energy corresponding to an interaction length equal to the size of the Earth. For tau neutrinos, this occurs at PeV energies, but can be considerably higher for neutralinos, due to their smaller cross section.

## 4 Prospects For Detection In Air Shower Experiments

The flux of very high energy neutralinos from top-down scenarios can be calculated assuming that this is the mechanism which generates the highest energy cosmic rays . Given a sufficient cosmic flux, these neutralinos may be detected in future air shower experiments. The challenge, however, is not merely observing the showers generated in neutralino interactions but in differentiating these cosmic neutralinos from neutrinos.

We have calculated the number of neutralino events predicted for a variety of top-down models associated with the highest energy cosmic rays in a future experiment such as EUSO  or OWL . EUSO and OWL are proposed satellite experiments which observe fluorescence in the Earth’s atmosphere generated in very high energy showers. Such experiments are expected to observe on the order of 150,000 square kilometers of surface area on the Earth. Particles which pass through the Earth can interact in the shallow Earth or atmosphere generating upgoing showers observable by fluorescence or Cerenkov radiation. Ultra-high energy showers reach a maximum near a slant depth of 850 $`\mathrm{g}/\mathrm{cm}^2`$, corresponding to a depth of 8.5 meters in water. Including the effective slant depth of the lower atmosphere extends this to $`0.015`$ km, thus providing a water equivalent effective volume of $`150,000\times 0.0152250`$ cubic kilometers, a truly enormous volume. Such an experiment will be capable of measuring both the energy and the direction of an observed particle.

Estimating the rate of neutrino–induced “background” events is difficult at present since the neutrino flux at $`E>10^9`$ GeV is not known. The flux of atmospheric neutrinos is completely negligible at these energies. However, most proposed explanations of the UHECR events also predict a significant UHE neutrino flux. We therefore use the neutrino flux predicted by top–down models  to estimate the neutrino background. Fig. 3 compares signal and background at $`E1`$ EeV for one such model, where we assume a galactic distribution of $`X`$ particles, with primary $`Xq\overline{q}`$ decay and $`M_X=210^{12}`$ GeV. We see that signal and background clearly have very different angular distributions even for the larger LSP–nucleon cross section of $`\sigma _\nu /10`$. Regeneration effects are included, but they cannot produce neutrino events at large energy and large angle. Requiring the events to emerge more than $`5^{}`$ below the horizon removes almost all the background, with little loss of signal; in the case at hand, we expect about 2 signal events per year, compared to 0.1 background event. If the LSP–nucleon cross section is smaller, a somewhat stronger angular cut may be advantageous; on the other hand, at even higher energies it might be better to use a slightly weaker cut. However, this variation of the angular cut has negligible effect on the predicted signal rate, compared to the uncertainty inherent in our estimates. In the following we therefore apply a fixed angular cut of $`5^{}`$ on the signal in all cases. This cut will have to be optimized once the angular resolution of the experiment is known. Moreover, measurements at neutrino telescopes as well as AUGER should soon greatly improve our knowledge of the neutrino flux at very high energies. Finally, this figure also shows that a measurement of the angular distribution of the signal will allow to determine the LSP scattering cross section: for the larger cross section shown, there will be very few vertically upgoing events. The dependence of the angular distribution of the signal on the cross section becomes even more pronounced at higher energies.

Table 1 shows signal event rates for two choices of energy threshold, $`E_{\chi ^0}`$ 1 EeV and 100 EeV. We also show results for the stronger cut on energy in order to illustrate that at least in some cases the LSP spectrum should be measurable over a significant range of energies. The first case shown in the table corresponds to the situation depicted in Fig. 3. Of course, the choice of a 100 EeV threshold is even more effective in reducing the background, to the level of $`10^3`$ events per year. From the physics point of view an energy threshold of 100 EeV should only be necessary in the unlikely case that the total background of ultra–high energy neutrinos is dominated by some mechanism not related to the observed UHECR events. Regarding the energy threshold which can be achieved experimentally, it has been argued that for upgoing events, the threshold could be as small as a PeV .

The rates shown in table 1 are for a variety of primary $`X`$ decay modes, and for “galactic” and homogeneous distributions of $`X`$ particles. It seems highly unlikely that $`X`$ particles will indeed be distributed homogeneously, but it is conceivable that the majority of sources contributing to the LSP flux is at cosmological distances (e.g. if the $`X`$ particles are embedded in topological defects); the homogeneous distribution is meant to be representative for such models. Our results show that the $`X`$ distribution throughout the universe has significant impact on the expected size of our signal. For a full description of these models, see our previous paper . We note that the neutralino signal is more sensitive to the primary $`X`$ decay mode than the neutrino signal analyzed in  is. Not surprisingly, scenarios with (at least) one superparticle in the primary decay produce a higher neutralino flux than models where $`X`$ only decays into quarks. Moreover, leptonic $`X`$ decays increase the predicted neutralino flux by another order of magnitude, since in this case relatively few protons are produced, leading to a higher source density required to explain the observed UHECR events. On the other hand, choosing $`M_X=210^{25}`$ eV rather than $`210^{21}`$ eV significantly reduces the predicted flux. Note, however, that in this case $`X`$ decays can only describe the UHECR flux above $`10^{20}`$ eV ; events at a few times $`10^{19}`$ eV then have to be produced by an as yet unknown source.

As stated earlier, we normalize the LSP flux by assuming that (almost) all UHE photons are absorbed between source and Earth, as indicated by experiment . Since absorption of UHE photons in our galaxy is speculative, we comment on how the expected signal is changed if this evidence is ignored, i.e. if the observe UHECR spectrum is normalized to the sum of photon and proton fluxes. The predicted LSP event rate for galactic models with $`M_X=210^{25}`$ eV would go down by about a factor of 4. If $`M_X=210^{21}`$ eV, the predicted event rate would go down by a factor of 2 to 3 for hadronic primary $`X`$ decays, and by about an order of magnitude for purely leptonic primary $`X`$ decay.The UHE neutrino background from $`X`$ decay would be reduced by the same factor. Note that this ”uncertainty” in the predicted event rate from taking refs. seriously or not is comparable to the variation between different primary $`X`$ decay modes. Finally, we remind the reader that the UHECR spectra measured by AGASA and HiRes differ significantly in the post–GZK region, where we normalize our fluxes, leading to a corresponding uncertainty in our predicted signal.

## 5 Conclusions

The cosmic neutralino flux predicted in top–down scenarios could possibly provide an interesting test of both supersymmetry and GUT scale particle physics. To identify any showers generated in future experiments as being generated by cosmic neutralinos, they will need to occur at energies and from directions at which neutrinos would be absorbed by the Earth. We have calculated the event rates for a variety of such models for a large area air shower experiment such as OWL or EUSO. We find that for many scenarios, the event rate is large enough to be observable in principle. We should mention here that our estimates of annual event rates assume 100% duty cycle. This is clearly not realistic for any experiment based on optical observations. However, planning for the kind of space–based experiment we envision is still in its early stage; a smaller duty cycle might be compensated by a larger area and/or a longer period of observation.

We believe that searching for UHE LSPs is very important, since it is the only measurement that can qualitatively distinguish between “top–down” and the more conventional “bottom–up” explanations for the observed UHE events: in bottom–up models superparticles can only be produced in the collision of accelerated protons, so the UHE LSP flux will be a tiny fraction \[typically $`𝒪(10^6)`$ or less\] of the UHE neutrino flux, much too small to be observed in any currently conceivable experiment. In contrast, a sizable UHE LSP flux is a generic prediction of top–down models. Moreover, the neutralino event rate turns out to be a far more sensitive probe of details of the model than the flux of neutrinos with energy exceeding $`1`$ PeV . We therefore find it encouraging that the observation of UHE LSPs along the lines suggested in this paper, while certainly not easy, should at least be possible.

Acknowledgments: We thank Andrew Strong for a useful discussion. This work was supported in part by a DOE grant No. DE-FG02-95ER40896 and in part by the Wisconsin Alumni Research Foundation. The work of M.D. was partially supported by the SFB375 of the Deutsche Forschungsgemeinschaft.


text:
# A meron cluster solution for the sign problem of the two-dimensional O(3) model

## 1 Introduction

The numerical simulation of many interesting physical systems, including in particular quantum chromodynamics (QCD) at non-zero chemical potential or at non-zero theta vacuum angle, suffers from a complex action problem or a sign problem (see  and references therein). The numerical simulation of these systems is impossible with standard Monte Carlo methods. In this article we investigate the two-dimensional O(3) model at a vacuum angle $`\theta =\pi `$. This model has a severe sign problem. Unlike gauge theories, for which no reliable cluster algorithms are available , the O(3) model can be simulated with the help of cluster algorithms and the meron cluster concept provides a solution of the sign problem. The results obtained in our studies are of interest for the understanding of certain systems in statistical mechanics. Haldane e.g. conjectured  that integer and half-integer one-dimensional antiferromagnetic quantum spin chains behave qualitatively differently. While integer spin chains have a mass gap, half-integer spin chains should be gapless. In the large $`s`$ limit ($`2s+1`$ being the number of possible spin orientations) there is a mapping from the one-dimensional antiferromagnetic quantum spin chain on the two-dimensional classical O(3) model with a $`\theta `$ vacuum term . Integer spin chains correspond to $`\theta =0`$ and half-integer spin chains correspond to $`\theta =\pi `$. The O(3) model at $`\theta =0`$ has a mass gap in agreement with Haldane’s conjecture. On the other hand, Haldane’s conjecture together with the (non-rigorous) mapping of spin chains on the O(3) model implies that the mass gap disappears at $`\theta =\pi `$. This corresponds to a phase transition governed by the vacuum angle $`\theta `$.

###

The paper is organized as follows. In section 2 the Wolff cluster algorithm and the sign problem of the two-dimensional O(3) model at $`\theta =\pi `$ are described. The meron clusters are defined. They are used for the construction of an improved estimator for the correlation function of two spins at $`\theta =\pi `$ (section 3). In order to calculate this correlation function a meron cluster algorithm is developed (section 4). The improved estimator together with the meron cluster algorithm will be shown to solve the sign problem. The correlation function at $`\theta =\pi `$ is calculated numerically. The results of these calculations are presented in section 5.

## 2 The Wolff cluster algorithm and the sign problem of the two-dimensional O(3) model

The lattice action of the two-dimensional O(3) model is given by

$$S=\beta \underset{xy}{}e_xe_y,\beta >0.$$
(1)
The summation is over all pairs $`xy`$ of nearest neighbor lattice sites. The three-component vectors of unit length are denoted by $`e_x`$ and called *spins*. The action has a global O(3) rotation symmetry.

In numerical simulations of this model one has to generate spin configurations $`\{e\}`$ that are distributed with a probability given by the statistical weight factor $`\mathrm{exp}(S(\{e\}))`$. In the Wolff cluster algorithm  some spins are combined and then updated together in one step. For every pair of nearest neighbor lattice sites $`x`$ and $`y`$ there is a *bond variable* $`b_{xy}`$ which only can take the values 0 and 1. If $`b_{xy}=1`$, one says that $`x`$ and $`y`$ are connected by a bond. The bond variables define the *clusters*. A cluster is a set of lattice sites. Lattice sites connected by bonds with each other belong to the same cluster.

The Wolff cluster algorithm consists of the following steps:

1. A vector $`r`$ of unit length is chosen whose direction (Wolff direction) is completely random.
2. To generate the clusters, nearest neighbor sites $`x`$ and $`y`$ are connected with probability

$$p(e_x,e_y)=1\mathrm{exp}(\mathrm{min}[0,2\beta (e_xr)(e_yr)]).$$
(2)
(bond activation probability) by a bond. The probability is largest if both spins $`e_x`$ and $`e_y`$ are aligned parallel or anti-parallel to the Wolff direction $`r`$. It is zero if the projections of the spins on the Wolff direction have different signs.
3. A cluster is flipped with probability $`1/2`$ independent of all the other clusters. Cluster flip means that all the spins belonging to the cluster are reflected on the plane perpendicular to the vector $`r`$ (Wolff plane)

$$e_xe_x^{}=e_x2(e_xr)r.$$
(3)
$`e_x^{}`$ denotes the reflected spin.

It has been shown that this algorithm fulfills detailed balance and is ergodic .

###

The spins of the two-dimensional O(3) model define a mapping from the lattice with periodic boundary conditions to the sphere $`S^2`$. These mappings fall into different classes characterized by an integer valued *topological charge* $`Q`$. For this characterization the spins must be interpolated between the lattice points. The charge $`Q`$ counts, how many times the sphere is covered by the mapping. The reason for the existence of the integer valued topological charge $`Q`$ is that the homotopy group of continuous mappings from the torus $`T^2`$ (the lattice) to the sphere $`S^2`$ is isomorphic to the group of integers.

We consider the two-dimensional O(3) model with an additional topological term in its action

$$SS\mathrm{i}\theta Q.$$
(4)
with $`\theta `$ real. A complex probability distribution $`\mathrm{exp}(S)`$ does not make sense. For the calculation of expectation values the additional factor $`\mathrm{exp}(\mathrm{i}\theta \mathrm{Q})`$ must be included in the observable

$$F\mathrm{exp}(\mathrm{i}\theta Q)F\stackrel{\theta =\pi }{=}(1)^QF.$$
(5)
For the special case $`\theta =\pi `$ the additional factor becomes $`(1)^Q`$. The factor $`(1)^Q`$ causes a *sign problem*. In a Monte Carlo calculation the alternating sign $`(1)^Q`$ leads to cancellations between the contributions from different configurations to the expectation value. The numerical calculation of the expectation value is difficult because it has a large variance if $`F`$ is one of the usual physical observables.

###

Using a Wolff cluster algorithm the flip of a cluster changes the topological charge from $`Q`$ to $`Q^{}`$. We can assign to every cluster a *cluster charge*

$$q=\frac{QQ^{}}{2}$$
(6)
that is either integer or half-integer. Clusters with half-integer charge are called *meron clusters*. Note that the flip of a meron cluster changes the sign $`(1)^Q^{}=(1)^{2q}(1)^Q=(1)(1)^Q`$. This property of the meron clusters is exploited for the construction of improved estimators at $`\theta =\pi `$ (see next section).

## 3 An improved estimator for the correlation function

The numerical calculation of the correlation function $`(1)^Qe_xe_y`$ of two spins at $`\theta =\pi `$ suffers from large statistical errors because of the sign problem. However the meron clusters can be used to construct an improved estimator (see eq. (3)) for the correlation function. This construction is based on a proposal by U.-J. Wiese.

Because of the O(3) rotation symmetry it is sufficient to consider the product of the parallel components of the spins (parallel to the Wolff direction). There are two possibilities: Either the lattice sites $`x`$ and $`y`$ belong to the same cluster $`𝒞`$ or they belong to different clusters

$`(1)^Qe_x^{}e_y^{}=`$ (7)
$`(1)^Qe_x^{}e_y^{}\mathrm{\Theta }(𝒞_x=𝒞_y)+(1)^Qe_x^{}e_y^{}\mathrm{\Theta }(𝒞_x𝒞_y).`$
The symbol $`\mathrm{\Theta }()`$ is 1 or 0, if the statement in brackets is true or false respectively. $`𝒞_x`$ denotes the cluster that contains the site $`x`$. The flip of a meron cluster changes the sign $`(1)^Q`$. Therefore when averaging over cluster flips the first term only gets a non-vanishing contribution if the cluster configuration contains no merons. If $`x`$ is contained in a meron cluster $`_x`$ the flip of the meron cluster changes $`(1)^Q`$ and $`e_x^{}`$ simultaneously. Therefore when averaging over cluster flips the second term only gets non-canceling contributions if $`x`$ and $`y`$ are both in meron clusters and there are no other meron clusters

$`(1)^Qe_x^{}e_y^{}=`$
$`(1)^Qe_x^{}e_y^{}\mathrm{\Theta }(𝒞_x=𝒞_y)\mathrm{\Theta }(N_{}=0)+(1)^Qe_x^{}e_y^{}\mathrm{\Theta }(_x_y)\mathrm{\Theta }(N_{}=2).`$
The factor $`(1)^Q`$ can be removed in both terms on the right hand side. Every cluster configuration can be thought of as being generated by cluster flips from a configuration in which all spins have been on the same side of the Wolff plane. For this configuration $`(1)^Q=1`$, because $`Q=0`$, and $`e_x^{}e_y^{}>0`$. The flip of a non-meron cluster does not change $`(1)^Q`$ and $`e_x^{}e_y^{}`$. The flip of a meron cluster simultaneously changes $`(1)^Q`$ and the sign of $`e_x^{}e_y^{}`$ in the case $`_x_y`$. Therefore

$`(1)^Qe_x^{}e_y^{}=`$
$`\underset{=g_{0\mathrm{M}}}{\underset{}{|e_x^{}e_y^{}|\mathrm{\Theta }(𝒞_x=𝒞_y)\mathrm{\Theta }(N_{}=0)}}+\underset{=g_{2\mathrm{M}}}{\underset{}{|e_x^{}e_y^{}|\mathrm{\Theta }(_x_y)\mathrm{\Theta }(N_{}=2)}}.`$
With equation (3) the sign problem has turned into the problem of finding cluster configurations with only 0 or 2 merons. Usually cluster configurations contain many merons and configurations with only 0 or 2 merons are quite rare. Therefore one aims at constructing an algorithm that generates only configurations with 0 and 2 merons. Such an algorithm is described in the next section. The improved estimator on the right hand side of equation (3) has no sign problem and therefore enables reliable numerical calculations.

The first term $`g_{0\mathrm{M}}`$ only gets non-vanishing contributions from cluster configurations in the zero-meron sector while the second term $`g_{2\mathrm{M}}`$ only gets non-vanishing contributions from the two-meron sector. In the term $`g_{0\mathrm{M}}`$ the sites $`x`$ and $`y`$ must be in the same cluster. Therefore, if the distance $`r=r(x,y)`$ between the lattice sites $`x`$ and $`y`$ is much larger than the cluster size, there is no contribution to $`g_{0\mathrm{M}}(r)`$. The decay of $`g_{0\mathrm{M}}(r)`$ is known to be mainly governed by the cluster size distribution which is exponential . On the other hand in the term $`g_{2\mathrm{M}}`$ the sites $`x`$ and $`y`$ must be in two different meron clusters. The two meron clusters can be far apart. Then $`g_{2\mathrm{M}}(r)`$ gets non-vanishing contributions at large distances $`r`$ (large compared to the cluster size). In this way long range correlations can arise.

In figure 1 the correlation functions at $`\theta =0`$ and at $`\theta =\pi `$ are compared. This figure shows the slower decay of correlation function at $`\theta =\pi `$.

## 4 A meron cluster algorithm

Assume that we have a cluster configuration with at most 2 merons. In the Wolff cluster algorithm the Wolff direction, used in the definition of the clusters, is not fixed. After a change of the Wolff direction, however, all bonds become invalid and have to be calculated anew. The new cluster configuration likely contains more than 2 merons.

Thus, a gradual update (that is an update which changes only a few bonds in each step) of the allowed cluster configuration is only possible if the Wolff direction is fixed. In this case the spins are only reflected on a single plane (perpendicular to the Wolff direction). These reflections do not change the spin components perpendicular to the Wolff direction. Therefore the resulting algorithm is not ergodic. To make it ergodic we must combine it with another update step effecting the perpendicular spin components. In this additional step only a few spins and bonds should be updated because we want to make only small changes on the allowed cluster configuration.

###

For the numerical calculations we use an algorithm with two Wolff directions. Wolff direction 1 is fixed during the whole simulation. The clusters and their charges are defined with respect to this Wolff direction 1. The cluster flips on the corresponding Wolff plane are not ergodic. To establish ergodicity we use a single cluster algorithm with a Wolff direction 2 that is chosen randomly anew for every update step. In one update step a single cluster is calculated and flipped with respect to the Wolff direction 2. After the single cluster flip the bonds in configuration 1 are calculated anew and the meron number constraint for configuration 1 is checked. If the meron number constraint is violated the single cluster flip is undone and the old bond configuration 1 is restored (accept reject step).

After the single cluster flip only those bonds in configuration 1 are calculated anew whose corresponding spins got changed. Thereby the change of the configuration of spins and bonds induced by the update step is restricted to a relatively small part of the lattice and the acceptance ratio for the update step becomes independent of the lattice size.

With the introduction of a second Wolff direction, the obvious problems with ergodicity are removed. There remains however still the posibility that independent of the details of the algorithm non-ergodic behaviour could arise if certain components of the zero- and two-meron sector are not accessible to our algorithm. It could be that configurations in the zero- and two-meron sector exist which could be reached only from sectors with more than two merons. This would invalidate any approach which deals with selected meron sectors only. We have not been able to rule out this possibility analytically. However we have not found any dependencies on the initial configuration in our results for the correlation function.

## 5 Numerical results

The meron cluster algorithm is used to investigate the large distance behavior of the correlation function of two spins at $`\theta =\pi `$. We calculate the projected correlation function at $`\theta =\pi `$

$$g_\mathrm{p}(r)=(1)^Q\left(\frac{1}{\sqrt{L_2}}\underset{x_2=0}{\overset{L_21}{}}e_{x_1,x_2}^{}\right)\left(\frac{1}{\sqrt{L_2}}\underset{y_2=0}{\overset{L_21}{}}e_{y_1,y_2}^{}\right)$$
(10)
(see figure 2) for various values of the lattice size $`L=L_1=L_2`$ at fixed $`\beta =0.1`$. The spins are summed over one space direction (the 2-direction) in order to suppress the contributions from excited states. This corresponds to a projection on the zero mode. The data points are fitted with a hyperbolic cosine

$$h(r)=A\mathrm{cosh}\left(m_\pi \left(\frac{L}{2}r\right)\right).$$
(11)
$`A`$ and $`m_\pi `$ are the fit parameters which must be determined. The $`L`$ dependence of the masses $`m_\pi `$ obtained in this way is shown in figure 3. The $`L`$ dependence of $`m_\pi `$ seems to be approximately of the form

$$m_\pi =\frac{C}{L},C=3.138.$$
(12)
The value of the constant $`C`$ is obtained from a fit. These data suggest that in the infinite volume limit the mass $`m_\pi `$ vanishes and thus long range correlations are present. The numerical data confirm the conjecture that the correlation function at $`\theta =\pi `$ is long ranged (cf. section 3). Also Haldane’s conjecture about the behavior of quantum spin chains is supported by the results (cf. section 1).

## 6 Summary

An improved estimator has been defined for the correlation function of two spins in the two-dimensional O(3) model at $`\theta =\pi `$. A meron cluster algorithm for the numerical calculation of the improved estimator has been developed. The numerical results indicate the occurrence of long range correlations at $`\theta =\pi `$.

## Acknowledgement

The author thanks J. Cox, O. Jahn, F. Lenz and U.-J. Wiese for help and discussions.