text:
# Improving the Accuracy of Cosmic Magnification Statistics

## 1 Introduction

Gravitational lensing by large-scale structures magnifies sources and distorts their images. The systematic distortion of faint background galaxies near matter overdensities, the *cosmic shear*, has been measured by several groups in the past few years (Bacon et al. 2000, 2002; Hämmerle et al. 2002; Hoekstra et al. 2002; Kaiser et al. 2000; Maoli et al. 2001; Réfrégier et al. 2002; Rhodes et al. 2001; Van Waerbeke et al. 2000, 2001, 2002; Wittman et al. 2000). It was found to be in remarkable agreement with theoretical predictions based on the Cold Dark Matter model, and has already provided new constraints on cosmological parameters (Van Waerbeke et al. 2001).

In a similar way, systematic magnifications of background sources near foreground matter overdensities, the *cosmic magnification*, can be measured and can provide largely independent constraints on cosmological parameters (Ménard & Bartelmann 2002, Ménard et al. 2002). Gravitational magnification has two effects: first, the flux received from distant sources is increased, and the solid angle in which they appear is stretched, thus their density is diluted. The net result of these competing effects depends on how the loss of sources due to dilution is balanced by the gain of sources due to flux magnification. Sources with flat luminosity functions, like faint galaxies, are depleted by cosmic magnification, while the number density of sources with steep luminosity functions, like quasars, is increased. Thus, cosmic magnification gives rise to apparent angular cross-correlations between background sources and foreground matter overdensities which are physically completely uncorrelated with the sources. These overdensities can be traced by using the distribution of foreground galaxies.

Numerous studies have confirmed the existence of quasar-galaxy correlations on angular scales ranging from one arc minute to about one degree, as expected from cosmic lensing (for a review, see Bartelmann & Schneider 2001; also Guimarães et al. 2001). However, the measured amplitude of these correlations has been systematically higher than usually theoretically predicted.

While cosmic shear can directly be related to observable quantities like image ellipticities, the theoretical interpretation of cosmic magnification involves several approximations:

* the luminosity function of the sources is described by a power-law over the range probed by the flux limit of the observation; and
* the magnification is assumed to fall into the weak lensing regime, i.e. to deviate weakly from unity. Thus, the magnification can with sufficient accuracy be approximated by its first-order Taylor expansion and its deviation from unity becomes proportional to the lensing convergence alone.

While the first assumption is comfortably satisfied, in particular for quasars, the validity of the second needs to be verified. This is the goal of the present paper.

Our paper is structured as follows: first, we introduce the formalism of the effective magnification and its Taylor expansion in Sect. 2. We then describe a number of statistics related to the lensing convergence, and evaluate the amplitude of the second-order terms which appear in the Taylor expansion. In Sect. 3, we describe the numerical simulations we use to test our analytical results and estimate the accuracy of several approximations for the magnification. As an application, we investigate second-order effects on quasar-galaxy correlations in Sect. 4, and we summarise our results in Sect. 5.

## 2 Formalism

### 2.1 Expanding the magnification

Cosmic magnification can be measured statistically through characteristic changes in the number density of the background sources. Along a given line-of-sight, this effect depends on two quantities:

* the magnification factor $`\mu `$, which describes whether sources are magnified or demagnified, depending on whether the matter along their lines-of-sight is preferentially over- or underdense compared to the mean,
* and the logarithmic slope $`\alpha `$ of the source counts as a function of flux, which quantifies the amplitude of source number-count modifications due to flux magnification. As mentioned in the introduction, magnification by gravitational lensing not only increases the observed flux, but also stretches the sky, thus the number density of sources on a magnified patch of the sky is reduced. The net magnification effect, called *magnification bias*, depends on the balance between the number of sources lost by dilution and gained by flux magnification. The steeper the number-count function of the sources is, the more pronounced is the magnification bias.

If the number-count function of the background sources can be described as a power law in a sufficiently wide range around the flux limit of the observation, the magnification bias is quantified by the *effective magnification* $`\mu ^{\alpha 1}`$. It directly expresses the changes of the background source density caused by lensing through the relation

$$n(>S,\theta )=\mu ^{\alpha 1}(\theta )n_0(>S)$$
(1)
where $`n_0(>S)`$ is the intrinsic number-count function of sources whose observed flux exceeds $`S`$ in the absence of lensing, and $`n(>S)`$ is the corresponding number-count function in presence of lensing.

The local properties of the gravitational lens mapping are characterised by the convergence $`\kappa `$, which is proportional to the surface mass density projected along the line-of-sight, and the shear $`\gamma `$, which is a two-component quantity and describes the gravitational tidal field of the lensing mass distribution. The effective magnification is related to $`\kappa `$ and $`\gamma `$ through

$$\mu ^{\alpha 1}=\left[(1\kappa )^2|\gamma |^2\right]^{1\alpha },$$
(2)
where $`|\gamma |=(\gamma _1^2+\gamma _2^2)^{1/2}`$ is taken as the absolute value of the shear. In the weak-lensing regime, both $`\kappa `$ and $`|\gamma |`$ are small compared to unity, and the previous expression can be expanded in a Taylor series:

$$\mu ^{\alpha 1}=1+(\alpha 1)\left[2\kappa +(2\alpha 1)\kappa ^2+|\gamma |^2\right]+𝒪(\kappa ^3,|\gamma |^3).$$
(3)
Previous studies using analytical formulae for magnification statistics focused only on the first-order term of this expansion, i.e. they used the approximation $`\mu ^{\alpha 1}1+2(\alpha 1)\kappa `$, which potentially causes the amplitude of the effect to be underestimated. In this section, we investigate the second-order terms in the expansion and estimate their contribution.

In doing so, we first note that $`\kappa ^2(\theta )`$ and $`|\gamma |^2(\theta )`$ share the same statistical properties (e.g. Blandford et al. 1991), because both $`\kappa `$ and $`\gamma `$ are linear combinations of second-order derivatives of the lensing potential. The identity of their statistics is most easily seen in Fourier space. Since we will only deal with ensemble averages of the magnification later on, $`\kappa ^2`$ and $`|\gamma |^2`$ can be combined into a single variable, which we denote by $`\kappa `$ for simplicity. Thus, we can write for our purposes,

$$\mu ^{\alpha 1}=1+2(\alpha 1)\left[\kappa +\alpha \kappa ^2\right]+𝒪(\kappa ^3).$$
(4)
Observable effects are due to departures from the mean value of the magnification. Therefore, the relevant quantity to correlate is $`\delta \mu ^{\alpha 1}=\mu ^{\alpha 1}\mu ^{\alpha 1}`$. Then, up to second order in $`\kappa ^2`$, the autocorrelation function of the effective magnification is

$`\delta \mu ^{\alpha 1}(\varphi )\delta \mu ^{\alpha 1}(\varphi +\theta )`$ $`=`$ $`4(\alpha 1)^2[\kappa (\varphi )\kappa (\varphi +\theta )`$ (5)
$`+`$ $`2\alpha \kappa (\varphi )\kappa ^2(\varphi +\theta )],`$
and the corresponding power spectrum can be expanded in a similar way,

$$P_{\mu ^{\alpha 1}}(s)=4(\alpha 1)^2\left[P_\kappa (s)+2\alpha P_{\mu ,2}(s)\right];$$
(6)
the power spectrum $`P_{\mu ,2}(s)`$ will be defined in Eq. (16) below. The last two equations show that the importance of the second-order terms in the expansion (3) increases as the number-count function of the background sources steepens, i.e. as $`\alpha `$ increases. In the following, we will use $`\alpha =2`$ unless stated otherwise. This value applies, for instance, to the number counts of bright quasars with $`m_\mathrm{B}<19.5`$ (Pei 1995). For simplicity, we abbreviate $`\delta \mu \delta \mu `$ by $`\mu \mu `$.

### 2.2 Second and Third-Order Correlations

We will now estimate several $`\kappa `$-related statistical quantities needed in the Taylor expansion of the magnification. For this purpose, we first introduce the $`\kappa `$ projector such that

$$\kappa (\theta )=_0^{w_\mathrm{H}}dwp_\kappa (w)\delta [\theta f_K(w),w]$$
(7)
can be written as a weighted line-of-sight projection of the density contrast $`\delta `$ from the observer to the Hubble distance $`w_\mathrm{H}`$. The projector is

$`p_\kappa (w)`$ $`=`$ $`{\displaystyle \frac{3}{2}}\mathrm{\Omega }_0\left({\displaystyle \frac{H_0}{c}}\right)^2`$ (8)
$`\times `$ $`{\displaystyle _w^{w_\mathrm{H}}}{\displaystyle \frac{\mathrm{d}w^{}}{a(w)}}n_\mathrm{S}(w^{}){\displaystyle \frac{f_K(w)f_K(w^{}w)}{f_K(w^{})}},`$
where $`w`$ is the radial coordinate distance, $`f_K(w)`$ is the comoving angular-diameter distance, $`n_\mathrm{S}(w)`$ is the normalised distance distribution of the sources, and $`a(w)`$ is the cosmological scale factor. Using Limber’s equation, we can then relate the autocorrelation function of $`\kappa `$ to the dark-matter power spectrum $`P_\delta `$,

$`\kappa (\varphi )\kappa (\varphi +\theta )`$ $`=`$ $`{\displaystyle dw\frac{p_\kappa ^2(w)}{f_K^2(w)}}`$ (9)
$`\times `$ $`{\displaystyle \frac{s\mathrm{d}s}{2\pi }P_\delta (\frac{s}{f_k(w)},w)\mathrm{J}_0(s\theta )},`$
where J0 is the zeroth-order Bessel function, and the power spectrum $`P_\kappa `$ corresponding to this correlation function is

$$P_\kappa (s)=dw\frac{p_\kappa ^2(w)}{f_K^2(w)}P_\delta (\frac{s}{f_k(w)},w).$$
(10)
As indicated by Eq. (5), the estimation of second-order terms requires the computation of the cross-correlation between $`\kappa `$ and $`\kappa ^2`$. We do this by first introducing a three-point correlation function for $`\kappa `$ and then identifying two of its three points. As usual, we define the three-point function by

$$z_\kappa (\theta _1,\theta _2)=\kappa (\varphi )\kappa (\varphi +\theta _1)\kappa (\varphi +\theta _2).$$
(11)
Using the $`\kappa `$ projector defined in (8), we can then write

$`z_\kappa (\theta _1,\theta _2)`$ $`=`$ $`{\displaystyle dw_1p_\kappa (w_1)dw_2p_\kappa (w_2)dw_3p_\kappa (w_3)}`$ (12)
$`\times `$ $`\delta [f_K(w_1)\varphi ,w_1]\delta [f_K(w_2)(\varphi +\theta _1),w_2]`$
$`\times `$ $`\delta [f_K(w_3)(\varphi +\theta _2),w_3].`$
Next, we employ the approximation underlying Limber’s equation, which asserts that the coherence length of the density fluctuation field is much smaller than the scales on which the projector $`p_\kappa `$ varies appreciably. Finally, we insert the expression for the bispectrum of the dark-matter fluctuations detailed in Appendix A, and find

$`z_\kappa (\theta _1,\theta _2)`$ $`=`$ $`{\displaystyle dwp_\kappa ^3(w)\frac{\mathrm{d}^2k_1}{(2\pi )^2}\mathrm{e}^{ik_1\theta _1f_K(w)}\frac{\mathrm{d}^2k_2}{(2\pi )^2}}`$ (13)
$`\times `$ $`B_\delta (k_1,k_2,k_1k_2,w)\mathrm{e}^{ik_2\theta _2f_K(w)},`$
where $`B_\delta (k_1,k_2,k_3)`$ is defined by

$$\widehat{\delta }(k_1)\widehat{\delta }(k_2)\widehat{\delta }(k_3)=\delta _D(k_1+k_2+k_3)B_\delta (k_1,k_2,k_3).$$
(14)
Then, using Eq. (12) and identifying two points of the three-point correlation function $`\theta _1\theta _2\theta `$, we find

$`\kappa (\varphi )\kappa ^2(\varphi +\theta )`$ $`=`$ $`{\displaystyle dwp_\kappa ^3(w)\frac{\mathrm{d}^2k_1}{(2\pi )^2}\mathrm{e}^{ik_1\theta f_K(w)}}`$ (15)
$`\times `$ $`{\displaystyle \frac{\mathrm{d}^2k_2}{(2\pi )^2}B_\delta (k_1,k_2,k_1k_2,w)}.`$
The term $`\kappa (\varphi )\kappa ^2(\varphi +\theta )`$ is a function of $`\theta `$ only. Its contribution $`P_{\mu ,2}(s)`$ to the power spectrum of the magnification is given by the inverse Fourier transform of Eq. (15), which reads

$`P_{\mu ,2}(s)`$ $`=`$ $`{\displaystyle 𝑑w\frac{p_\kappa ^3(w)}{f_K^4(w)}\frac{d^2s^{}}{(2\pi )^2}}`$ (16)
$`\times `$ $`B_\delta ({\displaystyle \frac{s^{}}{f_k(w)}},{\displaystyle \frac{s}{f_k(w)}},{\displaystyle \frac{s^{}s}{f_k(w)}},w).`$
### 2.3 Results and predictions

We can now numerically evaluate the first two contributions to the Taylor expansion of the magnification autocorrelation function defined in Eq. (5). As mentioned before, we use $`\alpha =2`$ here.

For evaluating the correlation functions, we use a CDM power spectrum in a spatially flat Universe parameterised with $`\mathrm{\Omega }_0=0.3`$, $`\sigma _8=0.9`$, $`h=0.7`$ and $`\mathrm{\Gamma }=0.21`$. The non-linear evolution of the power spectrum and the bispectrum are computed according to the formalisms developed by Peacock & Dodds (1996) and Scoccimarro et al. (2000), see Appendix A. The upper panel of Fig. 1 shows the first- and second-order contributions (dashed and dotted lines, respectively) to the Taylor expansion of the magnification for a fixed source redshift of $`z_\mathrm{s}=1`$. The sum of the two contributions is shown by the solid line. The figure shows that the contribution of the second-order term reaches an amplitude of more than $`30\%`$ of the first-order term on angular scales smaller than one arc minute. According to Eq. (5) which describes the Taylor expansion of the magnification autocorrelation, we define the relative contribution of the second-order compared to the first-order term as

$$R_{\mu \mu }(\theta )=\frac{2\alpha \kappa (\varphi )\kappa ^2(\varphi +\theta )}{\kappa (\varphi )\kappa (\varphi +\theta )}.$$
(17)
The lower panel shows this ratio in per cent for different source redshifts as a function of angular scale. From the lower to the upper curves, the source redshifts are $`1`$, $`1.5`$, $`2`$ and $`3`$. For each source redshift, the contribution of the second term exhibits a similar dependence on angular scale:

* on scales larger than a few degrees, the contribution drops to negligible values,
* effects become relevant on smaller scales, with a fairly constant amplitude from a few degrees down to around 10 arc minutes,
* on yet smaller scales, the second-order contribution increases steeply, due to the non-linear evolution of the density field. For sources at redshift $`2`$, the amplitude of the second term reaches *half* of the amplitude of the first term below one arc minute.

Thus, given the amplitude of $`R_{\mu \mu }`$, the correcting term introduced in Eq. (5) is relevant and must be taken into account for describing the magnification autocorrelation with an accuracy better than $`30\%50\%`$ on scales smaller than a few degrees.

So far, we have only investigated the amplitude contributed by the second-order term. In order to estimate the remaining contributions of all missing terms of the magnification expansion, we will now use numerical simulations allowing a direct computation of $`\mu `$ as a function of the convergence $`\kappa `$ and the shear $`\gamma `$.

## 3 Magnification statistics from numerical simulations

### 3.1 The ray-tracing simulation

For testing the theoretical predictions we performed ray-tracing experiments in a Very Large $`N`$-body Simulation (VLS) recently carried out by the Virgo Consortium (Jenkins et al. 2001, and see also Yoshida et al. 2001 for simulation details)111The ray-tracing data are available from T. Hamana on request, hamanatk@cc.nao.ac.jp The simulation was performed using a parallel P3M code (MacFarland et al. 1998) with a force softening length of $`l_{\mathrm{soft}}30h^1\mathrm{kpc}`$. The simulation employed $`512^3`$ CDM particles in a cubic box of $`479h^1\mathrm{Mpc}`$ on a side. It uses a flat cosmological model with a matter density $`\mathrm{\Omega }_0=0.3`$, a cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and a Hubble constant $`h=0.7`$. The initial matter power spectrum was computed using CMBFAST (Seljak & Zaldarriaga 1996) assuming a baryonic matter density of $`\mathrm{\Omega }_\mathrm{b}=0.04`$. The particle mass ($`m_{\mathrm{part}}=6.86\times 10^{10}h^1M_{}`$) of the simulation is sufficiently small to guarantee practically no discreteness effects on dark-matter clustering on scales down to the softening length in the redshift range of interest for our purposes (Hamana, Yoshida & Suto 2002).

The multiple-lens plane ray-tracing algorithm we used is detailed in Hamana & Mellier (2001; see also Bartelmann & Schneider 1992 and Jain, Seljak & White 2000 for the theoretical basics); we thus describe only aspects specific to the VLS $`N`$-body data in the following. In order to generate the density field between $`z=0`$ and $`z3`$, we use a stack of ten snapshot outputs from two runs of the $`N`$-body simulation, which differ only in the realisation of the initial fluctuation field. Each cubic box is divided into 4 sub-boxes of $`479^2\times 119.75h^3\mathrm{Mpc}^3`$ with the shorter box side being aligned with the line-of-sight direction. The $`N`$-body particles in each sub-box are projected onto the plane perpendicular to the shorter box side and thus to the line-of-sight direction. In this way, the particle distribution between the observer and $`z3`$ is projected onto $`38`$ lens planes separated by $`119.75h^1\mathrm{Mpc}`$. Note that in order to minimise the difference in redshift between a lens plane and an output of $`N`$-body data, only one half of the outputs (i.e. two sub-boxes) at $`z=0`$ are used.

The particle distribution on each plane is converted into the surface density field on either a $`1024^2`$ or $`2048^2`$ regular grid using the triangular shaped cloud (TSC) assignment scheme (Hockney & Eastwood 1988). The two grid sizes are adopted for the following reasons:

* the $`1024^2`$ grid is chosen to maintain the resolution provided by the $`N`$-body simulation and removing at the same time the shot noise due to discreteness in the $`N`$-body simulation. Its computation follows the procedure described in Hamana & Mellier (2001) and Jain et al. (2000). The corresponding outputs will be labelled with *large-scale smoothing* in the following.
* the $`2048^2`$ grid is also chosen to examine effects of small-scale nonlinear structures which are smoothed in the *large-scale smoothing* simulation. We should, however, note that in this case the shot noise is not sufficiently removed. Actually, the shot-noise power spectrum amplitude exceeds the convergence power spectrum on scales below $`1`$ arcmin. In the following, therefore, we will only consider measured correlation functions on scales larger than 1 arcmin. The corresponding outputs will be labelled with *small-scale smoothing* below.

Having produced surface density fields on all lens planes, $`1024^2`$ rays are traced backwards from the observer’s point using the multiple-lens plane algorithm (e.g. Schneider, Ehlers & Falco 1992). The initial ray directions are set on $`1024^2`$ grids with a grid size of $`0.25\mathrm{arcmin}`$, thus the total area covered by rays is $`4.27^2`$square degrees. We produced $`36`$ realizations of the underlying density field by randomly shifting the simulation boxes in the direction perpendicular to the line-of-sight using the periodic boundary conditions of the $`N`$-body boxes. Note that the lens planes coming from the same box are shifted in the same way in order to maintain the clustering of matter in the box.

We point out that second and higher-order statistics of point-source magnifications are generally ill-defined in presence of caustic curves because the differential magnification probability distribution asymptotically decreases as $`\mu ^3`$ for large $`\mu `$ (see Fig. 2). This is a generic feature of magnification near caustics and is thus independent of the lens model. Strong lensing effects on point sources near caustic curves give rise to rare, but arbitrarily high magnification values in the simulations, and therefore the variance of the measured statistics of $`\mu `$ cannot be defined. However, the smoothing procedure introduced above allows this problem to be removed because it smoothes out high density regions in the dark matter distribution and thus the fractional area of high magnification decreases. In reality, infinite magnifications do not occur, for two reasons. First, each astrophysical source is extended and its magnification (given the surface brightness-weighted point-source magnification across its solid angle) remains finite. Second, even point sources would be magnified by a finite value since for them, the geometrical-optics approximation fails near critical curves and a wave-optics description leads to a finite magnification (Schneider et al. 1992, Chap. 7).

### 3.2 Filtering

The computation of correlation functions from numerical simulations is mainly affected by two effects; on large scales by the finite box size of the dark matter simulation, and on small scales by the grid size used for computing the surface density field from the particle distribution. These boundaries set the limits for the validity of correlation functions measured in numerical simulations. In other words, this means that measuring a correlation function on a given scale is relevant only if this scale falls within the range of scales defined by the simulation. As shown in the previous section, our method for computing the cross-correlation between $`\kappa `$ and $`\kappa ^2`$ consists of first computing a three-point correlation function $`\kappa (\varphi )\kappa (\varphi +\theta _1)\kappa (\varphi +\theta _2)`$, and then identifying two of its three points. In such a case, one of the correlation lengths of the triple correlator becomes zero, thus necessarily smaller than the smallest relevant scales in any simulation. This prevents us from using any numerical simulation for directly comparing the results.

In order to avoid this problem, and for comparing our analytical with numerical results, we will introduce an effective smoothing into the theoretical calculations, such that each value of $`\kappa `$ at a given position $`\theta `$ is evaluated by averaging the $`\kappa `$-values in a disk of radius $`\theta _\mathrm{S}`$ centred on $`\theta `$. Indeed, the limit imposed by the grid size of the simulation gives rise to an unavoidable smoothing-like effect which cancels all information coming from scales smaller than a corresponding smoothing scale $`\theta _\mathrm{S}`$. For this purpose, we introduce a smoothed three-point correlator,

$$z_{\mathrm{smooth}}(\theta _2\theta _1,\theta _3\theta _1)=\kappa (\theta _1)\kappa (\theta _2)\kappa (\theta _3)_{\theta _\mathrm{S}}$$
(18)
$`=`$ $`{\displaystyle d\theta _1^{}d\theta _2^{}d\theta _3^{}\kappa (\theta _1^{})\kappa (\theta _2^{})\kappa (\theta _3^{})}`$
$`\times `$ $`W_{\theta _\mathrm{S}}(\theta _1^{}\theta _1)W_{\theta _\mathrm{S}}(\theta _2^{}\theta _2)W_{\theta _\mathrm{S}}(\theta _3^{}\theta _3),`$
where the function $`W_{\theta _\mathrm{S}}(\theta ^{})`$ is a normalised top-hat window of radius $`\theta _\mathrm{S}`$. Introducing this smoothing scheme into the expression for $`\kappa (\varphi )\kappa ^2(\varphi +\theta )`$ yields

$`z_{\mathrm{smooth}}(\theta _\mathrm{S})`$ $`=`$ $`{\displaystyle dw\frac{p_\kappa ^3(w)}{f_K^4(w)}\frac{\mathrm{d}^2s_1}{(2\pi )^2}\frac{\mathrm{d}^2s_2}{(2\pi )^2}}`$ (19)
$`\times `$ $`\mathrm{I}(s_1\theta _\mathrm{S})\mathrm{I}(s_2\theta _\mathrm{S})\mathrm{I}(|s_1+s_2|\theta _\mathrm{S})`$
$`\times `$ $`B_\delta ({\displaystyle \frac{s_1}{f_k(w)}},{\displaystyle \frac{s_2}{f_k(w)}},{\displaystyle \frac{s_1s_2}{f_k(w)}},w)e^{is_1\theta },`$
where $`\mathrm{I}(x)=2\frac{\mathrm{J}_1(x)}{x}`$. Similarly, introducing the smoothing scheme into the two-point correlation function gives

$`w_{\mathrm{smooth}}(|\theta _2\theta _1|)=\kappa (\theta _1)\kappa (\theta _2)_{\theta _\mathrm{S}}`$ (20)
$`=`$ $`{\displaystyle dw\frac{p_\kappa ^2(w)}{f_K^2(w)}\frac{\mathrm{d}^2s}{(2\pi )^2}P(\frac{s}{f_K(w)},w)}`$
$`\times `$ $`|\mathrm{I}(s\theta _\mathrm{S})|^2e^{is\theta }.`$
The effective smoothing scale depends on two parameters:

* the evolution of the apparent grid size of the simulation as a function of redshift, and
* the radial selection function of the dark-matter field whose correlation function has to be measured.

These quantities are plotted in Fig. 3. In order to use a unique smoothing scale valid on the final convergence map, we define the effective angular smoothing scale by

$$\theta _\mathrm{S}=dzW(z)\theta _{\mathrm{grid}}(z),$$
(21)
where $`W(z)`$ is the relevant normalised selection function along the line-of-sight. Measuring $`w_\mathrm{s}`$ means probing the power spectrum along the line-of-sight, weighted by $`p_\kappa ^2(z)`$. Therefore, we will use $`W(z)=D_+^2(z)p_\kappa ^2(z)`$, where $`D_+(z)`$ is the growth factor. In a similar way, we will use $`W(z)=D_+^4(z)p_\kappa ^3(z)`$ for measuring $`z_{\mathrm{smooth}}`$. The numerical values of the corresponding effective angles are presented in Table 1.

The second important difference between analytical calculations and measurements in numerical simulations is the finite box size effect. Indeed, the analytical correlation functions presented above were computed taking into account all modes in the power spectrum. However, the finite size of the box used in the simulation introduces an artificial cutoff in the power spectrum since wavelengths larger than the box size are not sampled by the simulation. This effect can also be taken into account in the analytical calculations by simply cancelling all the power on wavelengths with wave number $`k<k_{\mathrm{min}}`$. The boxes we use have a comoving size of $`480h^1\mathrm{Mpc}`$ which corresponds to $`k_{\mathrm{min}}=0.013h\mathrm{Mpc}^1`$.

### 3.3 Comparing $`\kappa \kappa `$ and $`\kappa \kappa ^2`$

With the help of the filtering schemes introduced in the previous section, we can now compare our theoretical predictions with correlation functions measured from the numerical simulations. We first compare the amplitude and angular variation of the two first terms of the Taylor expansion of the magnification separately. In the next section, we will then compare their sum to the total magnification fully computed from the simulation.

In Fig. 4, we overplot analytical and numerical results. The upper curve shows the autocorrelation function of $`\kappa `$ as a function of angular scale. We plot in circles the average measurement from $`36`$ realisations of the simulation, and the corresponding 1-$`\sigma `$ error bars to show the accuracy of the numerical results as a function of angular scale. The solid line shows the analytical prediction, including effective smoothing and an artificial cut of the power at scales below $`k_{\mathrm{min}}`$. The agreement is good on all scales. For comparison, the dotted line shows the result if we do not impose the large-wavelength cut, and the dashed line is the result if no cut and no smoothing are applied. In both cases, the deviations from the fully filtered calculation remain small since we are probing angular scales within the range allowed by the simulation.

The lower curves in Fig. 4 show a quantity proportional to the second-order correction of the Taylor expansion, namely the correlation function $`\kappa \kappa ^2`$. In the same way as before, the circles show average measurements from $`36`$ realisations, and the error bars denote the corresponding 1-$`\sigma `$ deviation. The prediction including smoothing and small-wavelength cut (solid line) shows a relatively good agreement given the expected accuracy of the bispectrum fitting formula, which is approximately 15% (Scoccimarro & Couchman 2000). This time, including smoothing changes the amplitude dramatically, and this effect affects all scales (see the dashed line). As discussed before, this is expected since we are measuring a three-point correlator on triangles which have one side length smaller than the angular grid size of the simulation. Finally, as shown by the difference between the dotted and solid lines, cancelling the power on scales where $`k<k_{\mathrm{min}}`$ again improves the agreement on large scales.

The agreement between our analytical and numerical computations of $`\kappa \kappa `$ and $`\kappa \kappa ^2`$ demonstrates the validity of the formalism introduced in Sect. 2 as well as the choice of the effective smoothing scale (Eq. ) for describing the second-order term in the Taylor expansion of the magnification.

### 3.4 Deviations from simulated $`\mu `$-statistics

We now want to investigate how well the second-order expansion describes the full magnification expression (2) which can be computed using maps of $`\kappa `$, $`\gamma `$ and $`\omega `$ (a net rotation term which arises from lens-lens coupling and the lensing deflection of the light ray path; see Van Waerbeke et al. 2001b) obtained from the simulations (see Hamana et al. 2000 for more detail).

Before doing so, we recall that the amplitude of the magnification autocorrelation measured from the simulation depends on the smoothing scale, as seen in Sect. 3.2, since $`\mu `$ is nonlinear in the density field. Therefore, all the following comparisons are valid for a given effective smoothing length only.

We further emphasise that two problems will complicate this comparison. First, our analytical treatment is valid in the weak-lensing regime only, i.e. as long as convergence and shear are small compared to unity, $`\kappa 1`$, $`|\gamma |1`$. While most light rays traced through the numerical simulations are indeed weakly lensed, a non-negligible fraction of them will experience magnifications well above two, say. Such events are restricted to small areas with high overdensities and thus affect the magnification statistics only at small angular scales. Second, a separate problem sets in if and where caustics are formed. The magnification of light rays going through caustics is infinite, and the magnification probability distribution near caustics drops like $`\mu ^3`$ for $`\mu \mathrm{}`$. As noted above, second- or higher-order statistics of $`\mu `$ then become meaningless because they diverge.

Departures of the numerical from the analytical results will thus have two distinct reasons, viz. the occurrence of non-weak magnifications which causes the analytical to underestimate the numerical results on small angular scales; and the formation of caustics, which causes second-order magnification statistics to break down entirely. Both effects will be demonstrated below. They can be controlled or suppressed in numerical simulations by smoothing, which makes lensing weaker, or by masking highly magnified light rays or regions containing caustics.

In Fig. 5, we plot with circles the autocorrelation function $`\mu (\varphi )\mu (\varphi +\theta )`$ measured from the large- and small-scale smoothing simulations in the left and right panels, respectively. The presence of caustics is more pronounced in the case of small-scale smoothing than in the large-scale smoothing simulations. The dotted line shows the theoretical prediction given by the first-order term of the Taylor expansion, namely $`4\kappa (\varphi )\kappa (\varphi +\theta )`$. This yields a low estimate of the correlation, with a discrepancy of order 10% on large scales, and more than 20% below a few arc minutes.

As expected from the preceding discussion, this level of discrepancy also depends on the effective smoothing scale and can increase if simulations with a smaller grid size are used. Estimating the contribution of the two lowest-order terms of $`\mu ^{\alpha 1}`$, we computed in Sect. 2.3 a lower bound to this discrepancy for a real case without smoothing, and found it to reach a level of 25% at large scales, and above 30% below a few arc minutes. The smoothed results taking the additional contribution of the second-order term into account are plotted as solid lines, and give a much better agreement, as expected. To quantify this in more detail, the lower panels of the figure show several contributions compared to the first-order term, i.e. to $`4\kappa \kappa `$.

* The symbols show the additional amplitude of the magnification statistics measured from the simulation, compared to the first-order term also obtained from the simulation,

$$R_{\mu \mu }^{\mathrm{num}}=\frac{\mu (\varphi )\mu (\varphi +\theta )_{\mathrm{num}}}{4\kappa (\varphi )\kappa (\varphi +\theta )_{\mathrm{num}}}1.$$
(22)
The error bars indicate the 1-$`\sigma `$ deviation across $`36`$ realisations.
* The solid line shows the contribution of the second-order relative to the first-order term computed from the analytical expression including the effective smoothing,

$$R_{\mu \mu }^{\mathrm{smooth}}=\frac{2\alpha \kappa (\varphi )\kappa ^2(\varphi +\theta )_{\theta _\mathrm{S}}}{\kappa (\varphi )\kappa (\varphi +\theta )_{\theta _\mathrm{S}}}$$
(23)
with $`\alpha =2`$.

In each case, we use the appropriate reference for $`\kappa (\varphi )\kappa (\varphi +\theta )`$, i.e. the numerical measurement in the first and the analytical estimation in the second case. Indeed, the measurement of $`\kappa \kappa `$ from the simulation agrees with the analytical estimation within some uncertainty, which is due to numerical effects like the finite number of dark-matter boxes used for simulating the light cone. It introduces a bias into our comparisons which is impossible to separate from the real offset due to all higher-order terms of the Taylor expansion that were not taken into account. The two contributions plotted in Fig. 5 are thus of different nature, but are suitable for a relative comparison.

As the lower panel of the large-scale smoothing simulation shows, the simple $`4\kappa \kappa `$ estimate of the magnification misses 20% of the real amplitude near one arc minute. This discrepancy almost vanishes after adding the contribution of the second-order term, which gives at all scales a final agreement on the per cent level: the additional amplitude reaches 19% at the smallest scales of the figure, compared to a value of 20% given by the simulation, and agrees within better than one per cent on larger scales. Therefore, taking into account the $`2\alpha \kappa \kappa ^2`$ correction allows the accuracy to be increased by a factor of $`20`$ compared to the approximation 4$`\kappa \kappa `$, in the case of our *large-scale smoothing* simulation. On the largest scales, between $`6`$ and $`30`$ arc minutes, the agreement even improves. Above these scales, the numerical results do not allow any relevant comparison because the number of available independent samplings corresponding to a given separation decreases. On scales below a few arc minutes, the offset between the measured points and the analytical estimate gives the amplitude of all higher-order terms neglected in the Taylor expansion of the magnification. As we can see, their contribution is on the one per cent level for the large-scale smoothing simulation.

The curves shown in the right panel demonstrate how the use of a smaller smoothing scale increases the discrepancy between the analytical and the numerical results. The fraction of non-weakly magnified light rays increases, and caustics appear which give rise to a power-law tail in the magnification probability distribution. We investigate the impact of the rare highly magnified light rays by masking pixels where the simulated magnification exceeds 4 or 8, and show that caustics have no noticeable effect on the amplitude of the magnification autocorrelation function determined from these simulated data. Note, however, that the impact of the caustics depends on the source redshift. The higher the redshift, the more caustics appear, and the larger is their impact on the correlation amplitude.

Imposing lower masking thresholds removes a significant fraction of the area covered by the simulation, changing the spatial magnification pattern and thus the magnification autocorrelation function. The corresponding measurements are represented by the dashed error bars in the lower right panel of Fig. 5. We note that the error bars of $`R_{\mu \mu }^{\mathrm{num}}`$ computed with the small-scale smoothing simulation become larger at small scales compared to the lower left panel. This reflects the fact that second-order magnification statistics are ill-defined once caustics appear. In the next section, we will investigate similar smoothing effects on cross-correlations between magnification and dark matter fluctuations. These quantities are not affected by problems of poor definition when the smoothing scale becomes small, and therefore do not show larger error bars at small scales when the smoothing scale decreases.

These comparisons show that the approximation $`\mu 1+2\kappa `$ misses a non-negligible part of the total amplitude of weak-lensing magnification statistics. The formalism introduced in Sect. 2 allows second-order corrections to be described with or without smoothing of the density field. This provides a better description of the correlation functions, but still gives a lower amplitude than the simulation results. As we noticed, the analytic computation based on the Taylor expansion is sufficiently accurate only in the weak lensing regime. In reality, however, the strong lensing, which can not be taken into account in the analytic formalism, has a significant impact on the magnification correlation especially at small scales as shown in the small-scale smoothing simulation. Therefore, one should carefully take the strong lensing effect into consideration when one interprets the magnification related correlation functions. However, we will see in the next section that counts-in-cells estimators are less affected by the strong lensing than correlation functions and thus enable better comparisons of observations with results from simulations.

## 4 Applications to quasar-galaxy correlations

As a direct application of the formalism introduced previously, we now investigate the effects of second-order terms on a well-known magnification-induced correlation, namely the quasar-galaxy cross-correlation (the results can also be applied to galaxy-galaxy correlations induced by magnification; Moessner & Jain 1998). In order to estimate cosmological parameters from this kind of correlations, we then suggest the use of a more suitable estimator using counts-in-cells rather than two-point correlation functions. It has the advantage of making the observational results more easily reconciled with the ones from numerical simulations.

### 4.1 Formalism and correcting terms

The magnification bias of large-scale structures, combined with galaxy biasing, leads to a cross-correlation of distant quasars with foreground galaxies. The existence of this cross-correlation has firmly been established (e.g. Benítez & Martínez-González 1995; Williams & Irwin 1998; Norman & Impey 1999; Norman & Williams 2000; Benítez et al. 2001; Norman & Impey 2001). Ménard & Bartelmann (2002) showed that the Sloan Digital Sky Survey (York et al. 2000) will allow this correlation function to be measured with a high accuracy. Its amplitude and angular shape contain information on cosmological parameters and the galaxy bias factor. Thus, it is important to accurately describe these magnification-related statistics in order to avoid a biased estimation of cosmological parameters as well as the amplitude of the galaxy bias.

As shown in Bartelmann (1995), the lensing-induced cross-correlation function between quasars and galaxies can be written as

$`w_{\mathrm{QG}}(\theta )`$ $``$ $`\delta _{\mathrm{QSO}}(\varphi )\delta _{\mathrm{gal}}(\varphi +\theta )`$ (24)
$`=`$ $`\delta \mu ^{\alpha 1}(\varphi )\delta _{\mathrm{gal}}(\varphi +\theta ).`$
Using the above formalism, we can expand the effective magnification fluctuation $`\delta \mu ^{\alpha 1}`$ up to second order and find the correcting term:

$$w_{\mathrm{QG}}(\theta )=2(\alpha 1)\left[\kappa \delta _{\mathrm{gal}}+\alpha \kappa ^2\delta _{\mathrm{gal}}\right].$$
(25)
The second term is proportional to $`\alpha `$ (contrary to the factor $`2\alpha `$ in Eq. ), since there is only one contribution of the magnification. Therefore, the expected effects will be roughly half of those on the autocorrelation of the effective magnification seen in the previous section. Assuming a linear bias $`b`$ between galaxies and dark matter, the cross-correlation between $`\delta _{\mathrm{gal}}`$ and $`\kappa ^2`$ can be written as

$`\delta _{\mathrm{gal}}(\varphi )\kappa ^2(\varphi +\theta )=b\delta _{DM}(\varphi )\kappa ^2(\varphi +\theta )`$ (26)
$`=`$ $`{\displaystyle dw\frac{p_\kappa ^2(w)p_\delta (w)}{f_K^4(w)}\frac{\mathrm{d}^2s_1}{(2\pi )^2}\frac{\mathrm{d}^2s_2}{(2\pi )^2}}`$
$`\times `$ $`B_\delta ({\displaystyle \frac{s_1}{f_k(w)}},{\displaystyle \frac{s_2}{f_k(w)}},{\displaystyle \frac{s_1s_2}{f_k(w)}},w)\mathrm{e}^{is_1\theta },`$
where $`p_\delta (w)`$ is the normalised distance distribution of the galaxies. For this example, we will use

$$p_\delta (z)\mathrm{d}z=\frac{\beta z^2}{z_0^3\mathrm{\Gamma }(3/\beta )}\mathrm{exp}\left[\left(\frac{z}{z_0}\right)^\beta \right]\mathrm{d}z,$$
(27)
with $`\beta =1.5`$ and $`z_0=0.3`$.

The results are shown in Fig. 6. As we can see, previous estimates using the approximation $`\mu 2\kappa `$ missed approximately $`15\%`$ of the amplitude on small scales for quasars at redshift unity. Using quasars at redshift $`2`$, these effects reach up to 25%. These offsets, which are only lower limits, would lead to biased estimates of $`\mathrm{\Omega }_0`$ or $`b`$, for example.

As for the magnification autocorrelation, we can compare our theoretical estimates against numerical estimations. We can first introduce a coefficient $`R_{\delta \mu }`$ describing the accuracy of our second-order correction:

$$R_{\delta \mu }(\theta )=\frac{\alpha \delta (\varphi )\kappa ^2(\varphi +\theta )}{\delta (\varphi )\kappa (\varphi +\theta )}.$$
(28)
We plot the results in Fig. 7. Note that contrary to the magnification autocorrelation, this quantity does not suffer from poor definition, even without smoothing. The difference can be seen by the same size of the error bars between the two simulation results at small scales, whereas they were larger in the case of $`\mu \mu `$ for the small-scale smoothing simulation (Fig. 5). The results for $`R_{\delta \mu }`$ are very similar those obtained for $`R_{\mu \mu }`$: for the large-scale smoothing ray-tracing we find very good agreement which reaches the one percent level on small scales. However, when the smoothing length decreases, we see from the small-scale smoothing outputs that we are missing a part of the total amplitude on small scales, which shows that higher-order terms play a non negligible role on those scales.

### 4.2 Using count-in-cells estimators

For precisely estimating cosmological parameters as well as the amplitude of the galaxy bias, it is necessary to employ theoretical magnification statistics that closely describe the observables. However, we have seen in Sect. 3 that analytical estimates as well as numerical simulations have intrinsic limitations and prevent us from accurately describing usual $`n`$-point correlation functions related to magnification statistics.

Besides, it is possible to focus on another estimator closely related to correlation functions, namely a count-in-cells estimator, which naturally smoothes effects originating from the density field and can thus more easily be reconciled with numerical simulations. So far, quasar-galaxy or galaxy-galaxy correlations have been quantified measuring the excess of background-foreground pairs at a given angular separation. Instead, we can correlate the amplitude of the background and foreground fluctuations, both measured inside a given aperture. We will therefore introduce a count-in-cells estimator,

$`\overline{w}_{\mathrm{QG}}(\theta )`$ $`=`$ $`\delta _{\mathrm{QSO}}(\varphi )\delta _{\mathrm{gal}}(\varphi )_\theta `$ (29)
$`\delta \mu ^{\alpha 1}(\varphi )\delta _{\mathrm{gal}}(\varphi )_\theta ,`$
where the subscript $`\theta `$ indicates averaging of $`\delta _{\mathrm{QSO}}(\varphi )`$ and $`\delta _{\mathrm{gal}}(\varphi )`$ inside a cell of radius $`\theta `$. In practice, this estimator is intended to be applied to galaxy-galaxy rather than to quasar-galaxy correlations, since the average angular separation between bright distant quasars is of order one degree for current surveys, thus averaging the source counts inside cells with radii of several arc minutes will not be relevant. Using galaxies as background sources, this limitation occurs at much smaller scales.

Using a first-order Taylor expansion for the magnification, the new estimator $`\overline{w}_{\mathrm{QG}}(\theta )`$ can be written

$`{\displaystyle \frac{\overline{w}_{\mathrm{QG}}(\theta )}{2(\alpha 1)b}}`$ $`=`$ $`{\displaystyle dw\frac{p_\kappa (w)p_\delta (w)}{f_K^2(w)}\frac{s\mathrm{d}s}{2\pi }}`$ (30)
$`\times `$ $`P({\displaystyle \frac{s}{f_k(w)}},w)|\mathrm{I}(s\theta )|^2,`$
where $`\mathrm{I}(x)=2\frac{\mathrm{J}_1(x)}{x}`$. This expression differs from the 2-point correlation function (9) by its Fourier-space filtering of the power spectrum. The additional smoothing wipes out the power on scales smaller than the physical scale corresponding to the angular smoothing scale $`\theta `$. For any observational result to be compared to a numerical simulation, $`\theta `$ and the smoothing scale used in the simulation will have to be carefully adapted to each other and to the redshift distribution of the foreground galaxy distribution.

Note that gravitational lensing by the foreground galaxies themselves is entirely irrelevant here. The angular scale on which galaxies act as efficient lenses is on the order of one arc second and below, much smaller than the angular scales we are concerned with. Moreover, the probability for a quasar to be strongly lensed by a galaxy is well below one per cent. Bartelmann & Schneider (1991) demonstrated this point explicitly by including galaxies into their numerical simulations and showing they had no noticeable effect.

## 5 Conclusion

As surveys mapping the large-scale structure of the Universe become wider and deeper, measuring cosmological parameters as well as the galaxy bias with cosmic magnification will become increasingly efficient and reliable. Therefore, an accurate theoretical quantification of magnification statistics becomes increasingly important.

Previous estimates of cosmic magnification relied on the assumption that the magnification deviates sufficiently little from unity that it can be accurately approximated by its first-order Taylor expansion about unity, i.e. $`\mu 1+2\kappa `$. In this paper, we have tested the validity of this assumption in the framework of magnification statistics, by investigating the second-order terms in the Taylor expansion of $`\mu `$. We have shown that:

* Second-order terms can be related to the cross-correlation between $`\kappa `$ and $`\kappa ^2`$,
* their importance increases as the number-count function of the background sources steepens, i.e. as $`\alpha `$ increases,
* their amplitude is *not* negligible: for the magnification autocorrelation, their contribution is typically on the order of 30%-50% at scales below one degree. Therefore, previous estimates of cosmic magnification were systematically low.

For testing our theoretical calculations, we have compared our results to magnification statistics found in numerical simulations by performing ray-tracing experiments in a very large $`N`$-body simulation. We have first checked the validity of our formalism describing the correlation $`\kappa \kappa ^2`$, and demonstrated the importance of including an effective smoothing into the analytical calculations. Indeed, $`\mu `$ is nonlinear in the density field and the amplitude of magnification statistics measured from numerical simulations depends therefore on the available resolution.

Using a simulation with an effective smoothing scale of 0.8 arc minutes, we found that our second-order formalism is accurate to the percent level for describing magnification autocorrelations. Compared to previous estimates, this improves the accuracy by a factor of $`20`$. For smaller effective smoothing scales, the contribution of third- and higher-order terms becomes important on scales below a few arc minutes.

Finally we have applied our formalism to observed correlations, like quasar-galaxy and galaxy-galaxy correlations due to lensing. We have shown that second-order corrections increase their amplitude by 15% to 25% on scales below one degree. These correlations are valuable tools to probe cosmological parameters as well as the galaxy bias. However, even including our correcting terms, analytical or numerical estimates of magnification statistics can only provide lower bounds to the real amplitude of the correlation functions in the weak-lensing regime. Thus, we propose using count-in-cells estimators rather than correlation functions since the intrinsic smoothing in determining counts-in-cells allows the observational results to be more directly related to those obtained in numerical simulations.

## Acknowledgements

We thank Francis Bernardeau and Stéphane Colombi for helpful discussions. This work was supported in part by the TMR Network “Gravitational Lensing: New Constraints on Cosmology and the Distribution of Dark Matter” of the EC under contract No. ERBFMRX-CT97-0172.

## Appendix A Bispectrum and Non-Linear evolution

The bispectrum can be estimated using second-order perturbation theory. Indeed, an expansion of the density field to second nonlinear order as

$$\delta (x)\delta ^{(1)}(x)+\delta ^{(2)}(x),$$
(31)
where $`\delta ^{(2)}`$ is of order $`(\delta ^{(1)})^2`$ and represents departures from Gaussian behaviour, yields the bispectrum

$`\delta _1\delta _2\delta _3`$ $``$ $`\delta _1^{(1)}\delta _2^{(1)}\delta _3^{(1)}+\delta _1^{(1)}\delta _2^{(1)}\delta _3^{(2)}`$ (32)
$`+`$ $`\text{cyclic terms (231, 312)}.`$
The first term in Eq. (32) vanishes because the density fluctuation field is Gaussian to first order, hence the third moment of $`\delta ^{(1)}`$ is zero. Thus, the leading term in Eq. (32) is of the order of $`\delta _1^{(1)}\delta _2^{(1)}\delta _3^{(2)}`$ and can be quantified using second-order perturbation theory.

The bispectrum $`B_\delta (𝐤_1,𝐤_2,𝐤_3)`$ is defined only for closed triangles formed by the wave vectors $`𝐤_1,𝐤_2,𝐤_3`$. It can be expressed as a function of the second-order kernel $`F(𝐤_1,𝐤_2)`$ and the power spectrum

$`B_\delta (𝐤_1,𝐤_2,𝐤_3)`$ $`=`$ $`2F(𝐤_1,𝐤_2)P(k_1)P(k_2)`$ (33)
$`+`$ $`2F(𝐤_2,𝐤_3)P(k_2)P(k_3)`$
$`+`$ $`2F(𝐤_1,𝐤_3)P(k_1)P(k_3).`$
For describing the bispectrum on all scales, we use the fitting formula derived by Scoccimarro & Couchman (2000) for the non-linear evolution of the bispectrum in numerical simulations of CDM models, extending previous work for scale-free initial conditions. In that case, we have

$`F(𝐤_1,𝐤_2)`$ $`=`$ $`{\displaystyle \frac{5}{7}}a(n,k_1)a(n,k_2)`$ (34)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{𝐤_1𝐤_2}{k_1k_2}}\left({\displaystyle \frac{k_1}{k_2}}+{\displaystyle \frac{k_2}{k_1}}\right)b(n,k_1)b(n,k_2)`$
$`+`$ $`{\displaystyle \frac{2}{7}}\left({\displaystyle \frac{𝐤_1𝐤_2}{k_1k_2}}\right)^2c(n,k_1)c(n,k_2),`$
with

$`a(n,k)`$ $`=`$ $`{\displaystyle \frac{1+\sigma _8^{0.2}(z)\left[0.7Q_3(n)\right]^{1/2}(q/4)^{n+3.5}}{1+(q/4)^{n+3.5}}}`$ (35)
$`b(n,k)`$ $`=`$ $`{\displaystyle \frac{1+0.4(n+3)q^{n+3}}{1+q^{n+3.5}}}`$ (36)
$`c(n,k)`$ $`=`$ $`{\displaystyle \frac{1+4.5/\left[1.5+(n+3)^4\right](2q)^{n+3}}{1+(2q)^{n+3.5}}},`$ (37)
and $`qk/k_{\mathrm{NL}}(z)`$, where $`4\pi k_{\mathrm{NL}}^3P_L(k_{\mathrm{NL}})=1`$, and $`P_\mathrm{L}(k)`$ is the linear power spectrum at the desired redshift. The effective spectral index is taken from the linear power spectrum as well. The function $`Q_3(n)`$ is given by

$$Q_3(n)=\frac{(42^n)}{(1+2^{n+1})},$$
(38)
For more detail, see Scoccimarro & Couchman (2000).


format: text

text:
# 1 Introduction

## 1 Introduction

Glassy dynamics is often described as motion in a complex energy landscape , with many valleys separated by barriers. Valleys have been given various definitions , which are not equivalent from a dynamical point of view . One common feature of all these approaches is the exponential growth with the system size of the number of valleys at fixed energy density $`E`$:

$$𝒩(N;E)\mathrm{exp}(NS_{\mathrm{ap}}(E)),$$
(1.1)
where $`S_{\mathrm{ap}}(E)`$ is the ‘a priori’ configurational entropy, or complexity.

As valleys are so numerous, an important issue concerns the possibility of giving an a priori statistical description of them. The key question to address is then the following. Assuming the initial configuration be random, does the system sample all the possible valleys with equal statistical weights, i.e., with a uniform or flat measure, or, to the contrary, is the size of the basin of attraction of each valley relevant? The assumption that an a priori statistical description with a flat measure could hold was first formulated by Edwards  in the context of granular materials, and is therefore referred to as the Edwards hypothesis.

One-dimensional spin systems at zero temperature provide an adequate setting for the investigation of the dynamical role of valleys, and especially for testing the Edwards hypothesis. Indeed, the blocked configurations reached by zero-temperature dynamics are truly metastable states with infinite lifetimes. In an earlier work  we have addressed this question for Ising chains evolving under irreversible zero-temperature dynamics, where each spin may flip at most once before a blocked state is reached. These dynamical models can be exactly mapped onto random sequential adsorption problems , for which analytical tools are available.

In the present work we pursue the investigation of metastable states in Ising chains, by considering the richer situation of zero-temperature Kawasaki (spin-exchange) dynamics . We allow the rate $`W_0`$ of diffusion processes (constant-energy moves) to vary . If the kinetic constraint $`W_0=0`$ is imposed , only irreversible processes are present . For the unconstrained dynamics $`(W_0>0)`$, where diffusion processes are allowed, the dynamics is partly reversible. Each spin only flips a finite number of times, before the system globally reaches a blocked state after a finite time. The statistics of this blocking time is studied in Section 2. Our main goal is then again to test the Edwards hypothesis, i.e., to evaluate various observables in the blocked configurations reached by the dynamics, starting from a random non-magnetized initial configuration. We shall systematically compare the data of extensive numerical simulations, presented in Section 4, with the predictions of the a priori approach, derived in Section 3. Finally, some aspects of persistence, in particular the distribution of the number of flips of a given spin, are discussed in Section 5.

To be more specific, we consider a ferromagnetic chain of Ising spins $`\sigma _n=\pm 1`$, whose Hamiltonian reads

$$=\underset{n}{}\sigma _n\sigma _{n+1}.$$
(1.2)
In Kawasaki dynamics , only pairs of opposite spins may be flipped $`(++)`$, so that the magnetization is locally conserved.

For simplicity, we limit ourselves to Monte-Carlo dynamics with random sequential updating, and we assume that the flipping rate only depends on the energy difference $`\delta `$ involved in the proposed move. Any given pair of opposite spins thus has a probability $`W_\delta `$ of flipping per unit time, with

$$\delta =2(\sigma _{n1}\sigma _n+\sigma _{n+1}\sigma _{n+2})\{4,0,4\}.$$

The requirement of detailed balance with respect to the Hamiltonian (1.2) at temperature $`T=1/\beta `$ yields a single condition,

$$\frac{W_4}{W_4}=\mathrm{e}^{4\beta },$$

upon the rates $`W_4`$, $`W_0`$, $`W_4`$.

We furthermore restrict ourselves to zero-temperature dynamics, hence $`W_4=0`$. We choose time units such that $`W_4=1`$, keeping $`W_0`$ as a free parameter. The allowed moves and the corresponding rates are listed in Table 1. The zero-temperature limits of the Metropolis and heat-bath rules correspond respectively to $`W_0=1`$ and $`W_0=1/2`$. Hereafter we focus our attention on the range $`0W_01`$.

## 2 Statistics of blocking time

Let us first recall that for the kinetically constrained model $`(W_0=0)`$ a finite system consisting of $`N`$ spins reaches a blocked configuration after a finite blocking time $`T_N\mathrm{ln}N`$, which is the jamming time of the equivalent problem of random sequential adsorption of hollow trimers . The blocked configurations are characterized by the property that the spin patterns $`++`$ and $`++`$ are absent. Equivalently, there are at most two consecutive unsatisfied bonds.

In the present case $`(W_0>0)`$ the system still gets trapped in a blocked configuration. However the diffusive motion of free $`+`$ spins in domains of $``$ spins, and vice-versa, is allowed. Each free spin will eventually be annihilated, by meeting either another free spin or one of the boundaries of the domain. Blocked configurations of Kawasaki dynamics are therefore characterized by the property that the patterns $`++`$ and $`+`$ are absent. Equivalently, isolated spins are absent, or unsatisfied bonds are isolated.

In order to understand the statistics of the blocking time $`T_N`$, we consider first the regime $`W_01`$, where the time scales of condensation and diffusion are well separated . The fast part of the dynamics, which occurs with unit rate, is identical to the constrained dynamics considered in Refs. . For intermediate times of order $`1t1/W_0`$, the system is therefore approximately left in one of the final configurations of the constrained dynamics. The slow, diffusive part of the dynamics then takes place at rate $`W_01`$. The late stages of the dynamics are governed by large domains, on which a single free spin diffuses. We assume that such large domains, of size $`L1`$, occur with an exponentially small probability

$$f_{\mathrm{dif}}(L)\mathrm{exp}(L/\xi _{\mathrm{dif}}),$$

where $`\xi _{\mathrm{dif}}`$ is the relevant characteristic length. The spin diffusion constant reads $`D=W_0`$ in our units. For a random initial point, the survival probability is known  to decay as $`S(t;L)(8L/\pi ^2)\mathrm{exp}(\pi ^2W_0t/L^2)`$. The mean density of free spins at time $`t`$ can therefore be estimated as

$$S(t)\underset{L}{}f_{\mathrm{dif}}(L)S(t;L)_0^{\mathrm{}}\mathrm{exp}\left(\frac{L}{\xi _{\mathrm{dif}}}\frac{\pi ^2W_0t}{L^2}\right)LdL.$$
(2.1)
Evaluating the integral by the method of steepest descent, we thus obtain stretched exponential decay for the density of free spins:

$$S(t)t^{1/2}\mathrm{exp}\left(\frac{(W_0t)^{1/3}}{A_{\mathrm{dif}}}\right),$$
(2.2)
with

$$A_{\mathrm{dif}}=\left(\frac{4\xi _{\mathrm{dif}}^2}{27\pi ^2}\right)^{1/3}.$$
(2.3)
This behavior, already emphasized in Ref. , is a general characteristic feature of diffusion processes in the presence of random traps . The exponent $`1/3`$ is related to the one-dimensional geometry; it would read $`d/(d+2)`$ for trapping problems in higher spatial dimension $`d`$. The above result should hold for the late stages $`(W_0t1)`$ of Kawasaki dynamics, for any finite value of $`W_0`$. The characteristic length $`\xi _{\mathrm{dif}}`$ is expected to depend smoothly on $`W_0`$.

For a finite system of $`N`$ spins, the last free spin will be annihilated at a time $`T_N`$ such that $`NS(T_N)1`$, hence $`(W_0T_N)^{1/3}A_{\mathrm{dif}}\mathrm{ln}N`$. More precisely, as the histories of spins diffusing on different domains are statistically independent, it can be argued along the lines of  that $`T_N`$ is given according to extreme-value statistics  as

$$(W_0T_N)^{1/3}A_{\mathrm{dif}}\left(\mathrm{ln}\left(N(\mathrm{ln}N)^{3/2}\right)+b+X_N\right).$$
(2.4)
In this expression, the factor $`(\mathrm{ln}N)^{3/2}`$ takes account of the $`t^{1/2}`$ prefactor in the survival probability (2.2), while the effective constant $`b`$ encompasses all subleading effects that have been neglected, and the random variable $`X_N`$ is distributed according to the Gumbel law

$$f(X)=\mathrm{exp}(X\mathrm{e}^X).$$
(2.5)
These predictions have been checked against extensive numerical simulations, performed according to the zero-temperature Kawasaki dynamics summarized in Table 1, with random sequential updates and periodic boundary conditions, starting from a random initial configuration with zero magnetization. The dynamics is run until the system gets trapped in a blocked configuration. The blocking time $`T_N`$ is recorded for each sample. This measurement has been performed for many samples ($`10^8`$ spins in total for each value of $`W_0`$ and of $`N`$). The mean $`(W_0T_N)^{1/3}`$ is found to follow an almost perfect linear law when plotted against $`\mathrm{ln}(N(\mathrm{ln}N)^{3/2})`$, at least for $`N`$ ranging from 50 to 3200, whereas the same data plotted against $`\mathrm{ln}N`$ are bent in a significant way. Using (2.3), the slope of the latter plot yields the value of the characteristic length $`\xi _{\mathrm{dif}}`$, which is plotted in Figure 1, against $`W_0`$. The error on this estimate, containing a systematic and a statistical part, is roughly comparable to the symbol size. The length $`\xi _{\mathrm{dif}}`$ exhibits a rather weak dependence on $`W_0`$, decreasing from $`\xi _{\mathrm{dif}}=1.77`$ in the $`W_00`$ limit to $`\xi _{\mathrm{dif}}=1.61`$ for $`W_0=1`$.

Another confirmation of the above picture is provided by Figure 2, showing a histogram plot of the variable $`X_N`$ defined by (2.4), for 125,000 samples of size $`N=800`$, with $`W_0=1`$. The parameter $`A_{\mathrm{dif}}=0.339`$, i.e., $`\xi _{\mathrm{dif}}=1.61`$, is taken from the data of Figure 1, while the constant $`b`$ is chosen by fitting the average $`X_N`$. A convincing agreement is found with the limit law (2.5).

## 3 Blocked states: a priori ensemble

As already stated, the blocked states of zero-temperature Kawasaki dynamics are the configurations where isolated spins are absent. Equivalently, unsatisfied bonds are isolated.

This section is devoted to the statistical description of the blocked configurations taken with equal a priori weights (Edwards approach). We shall distinguish the full ensemble of all the blocked configurations, irrespective of their energy, the restricted ensemble of blocked configurations with prescribed energy density $`E`$, and the canonical ensemble of blocked configurations, obtained by fixing the parameter $`\beta `$ conjugate to the energy density $`E`$ (see below).

### 3.1 Statistics of energy and configurational entropy

For a finite chain of $`N`$ spins, we first investigate the number $`𝒩(N;E)`$ of blocked configurations with prescribed energy density $`E`$. This number can be evaluated by an elementary combinatorial reasoning . For the sake of generality, we prefer to resort to the transfer-matrix formalism. We introduce the partition function $`Z_N(\beta )`$, defined as a sum over all the blocked configurations $`𝒞=\{\sigma _n\}`$ of the Boltzmann weight associated with the Hamiltonian (1.2). We have

$$Z_N(\beta )=\underset{𝒞}{}\mathrm{e}^{\beta (𝒞)}𝒩(N;E)\mathrm{e}^{\beta NE}dE.$$
(3.1)
In this framework, the parameter $`\beta `$ can be positive or negative, and it is not related to physical temperature. The transfer matrix is then a very useful tool. Indeed the partition functions $`Z_N^\alpha `$, labeled by the prescribed value $`\alpha =(\sigma _{N1},\sigma _N)`$ of the last two spins, obey the recursion relation

$$\left(\begin{array}{c}Z_{N+1}^{++}\\ Z_{N+1}^+\\ Z_{N+1}^+\\ Z_{N+1}^{}\end{array}\right)=𝒯\left(\begin{array}{c}Z_N^{++}\\ Z_N^+\\ Z_N^+\\ Z_N^{}\end{array}\right),$$

where the $`4\times 4`$ transfer matrix $`𝒯`$ reads

$$𝒯=\left(\begin{array}{cccc}\mathrm{e}^\beta & 0& \mathrm{e}^\beta & 0\\ \mathrm{e}^\beta & 0& 0& 0\\ 0& 0& 0& \mathrm{e}^\beta \\ 0& \mathrm{e}^\beta & 0& \mathrm{e}^\beta \end{array}\right).$$

The characteristic polynomial of $`𝒯`$ factors as $`det(\lambda \mathrm{𝟏}𝒯)=(\lambda ^2\mathrm{e}^\beta \lambda 1)(\lambda ^2\mathrm{e}^\beta \lambda +1)`$. The first (resp. the second) factor yields eigenvalues $`\lambda _1`$, $`\lambda _2`$ (resp. $`\lambda _3`$, $`\lambda _4`$), with

$$\lambda _{1,2}=\frac{1}{2}\left(\mathrm{e}^\beta \pm \sqrt{\mathrm{e}^{2\beta }+4}\right),\lambda _{3,4}=\frac{1}{2}\left(\mathrm{e}^\beta \pm \sqrt{\mathrm{e}^{2\beta }4}\right).$$

The left and right eigenvectors $`L_a|`$ and $`|R_a`$ are even (resp. odd) for $`a=1,2`$ (resp. $`a=3,4`$) under the spin symmetry $`+`$. We assume that they are normalized so that $`L_a|R_b=\delta _{ab}`$. Their explicit expressions will not be needed hereafter.

For large $`N`$, we have $`Z_N(\beta )\lambda _1^N`$, as $`\lambda _1`$ is the largest eigenvalue. Using (3.1), we obtain an exponential law of the form (1.1) for $`𝒩(N;E)`$, where the a priori configurational entropy $`S_{\mathrm{ap}}(E)`$ is related to $`\mathrm{ln}\lambda _1(\beta )`$ by a Legendre transform:

$$S_{\mathrm{ap}}(E)\mathrm{ln}\lambda _1(\beta )=\beta E,E=\frac{\mathrm{d}\mathrm{ln}\lambda _1}{\mathrm{d}\beta },\beta =\frac{\mathrm{d}S_{\mathrm{ap}}}{\mathrm{d}E}.$$

Explicitly, we have

$$E=\frac{\mathrm{e}^\beta }{\sqrt{\mathrm{e}^{2\beta }+4}},\mathrm{e}^\beta =\frac{2E}{\sqrt{1E^2}},$$
(3.2)
and the a priori entropy reads

$$S_{\mathrm{ap}}(E)=E\mathrm{ln}(2E)+\frac{1E}{2}\mathrm{ln}(1E)\frac{1+E}{2}\mathrm{ln}(1+E).$$

This entropy is non-zero for $`1<E<0`$. It takes its maximal value

$$S{}_{}{}^{}=\mathrm{ln}\mathrm{\Phi }=0.481212,$$

where $`\mathrm{\Phi }=(1+\sqrt{5})/2`$ is the golden mean, for

$$E{}_{}{}^{}=\frac{1}{\sqrt{5}}=0.447214,$$
(3.3)
corresponding to $`\beta =0`$. Equation (3.3) therefore yields the typical a priori energy density of a blocked configuration.

The result (1.1) can be recast as follows. Consider the full ensemble of all the blocked configurations, irrespective of their energy. The probability of observing, in that ensemble, a blocked configuration with energy density $`E`$ reads

$$P_{\mathrm{ap}}(E)\mathrm{exp}(N\mathrm{\Sigma }_{\mathrm{ap}}(E)),\mathrm{\Sigma }_{\mathrm{ap}}(E)=S{}_{}{}^{}S_{\mathrm{ap}}(E).$$
(3.4)
The function $`\mathrm{\Sigma }_{\mathrm{ap}}(E)`$ vanishes quadratically as

$$\mathrm{\Sigma }_{\mathrm{ap}}(E)c(EE{}_{}{}^{})^2,c=\frac{5\sqrt{5}}{8}.$$

The bulk of the a priori distribution of $`E`$ is therefore a narrow Gaussian around $`E^{}`$, whose rescaled variance asymptotically reads

$$NVarE=N(E^2E{}_{}{}^{}{}_{}{}^{2})\underset{N\mathrm{}}{}\frac{1}{2c}=\frac{4\sqrt{5}}{25}=0.357771.$$

The above results also allow to determine the higher cumulants of the energy in the canonical ensemble (at fixed parameter $`\beta `$). One has indeed

$$\frac{1}{N}\mathrm{ln}\mathrm{e}^{sNE}_\beta =\frac{1}{N}\mathrm{ln}\frac{Z_N(\beta s)}{Z_N(\beta )}\underset{N\mathrm{}}{}\mathrm{ln}\frac{\lambda _1(\beta s)}{\lambda _1(\beta )}.$$

By expanding this result as a power series in $`s`$, we obtain explicit expressions for the cumulants $`E^k`$ as a function of $`\beta `$.

As far as mean quantities are concerned, the microcanonical or restricted ensemble (fixed energy density $`E`$) and the canonical one (fixed conjugate parameter $`\beta `$) are equivalent. It is therefore justified to recast canonical results in terms of the prescribed value $`E`$ of the mean energy, using (3.2). Generalizing this equivalence prescription to the cumulants of the energy, we obtain

$`NVarE=NE^2`$ $`\underset{N\mathrm{}}{}`$ $`E(1E^2),`$
$`N^2E^3`$ $`\underset{N\mathrm{}}{}`$ $`E(1E^2)(3E^21),`$ (3.5)
$`N^3E^4`$ $`\underset{N\mathrm{}}{}`$ $`E(1E^2)(15E^412E^2+1).`$
The results (3.5) can be given the following interpretation. Fixing $`\beta `$ amounts to fixing the extensive part $`NE`$ of the energy, while the equivalence prescription is a natural Ansatz to describe the fluctuations of its non-extensive part.

### 3.2 Spin correlation function

The spin correlation function $`C_n=\sigma _0\sigma _n`$ in the canonical a priori ensemble can also be evaluated by the transfer-matrix method. In the bulk of an infinitely long chain, and for $`n0`$, we have

$$C_n=\frac{L_1|𝒮𝒯^n𝒮|R_1}{\lambda _1^n}=\underset{a}{}L_1|𝒮|R_aL_a|𝒮|R_1\left(\frac{\lambda _a}{\lambda _1}\right)^n,$$
(3.6)
where

$$𝒮=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)$$

is the spin operator.

Because of symmetry, only the eigenvectors with $`a=3,4`$ contribute to the above sum. The values $`C_0=1`$ and $`C_1=E`$ allow to determine the products of matrix elements entering (3.6), without knowing the eigenvectors explicitly. The following alternative reasoning can also be used. Being a linear combination of $`(\lambda _3/\lambda _1)^n`$ and $`(\lambda _4/\lambda _1)^n`$, $`C_n`$ can be shown to obey the three-term recursion relation

$$(1E)C_{n+2}+2EC_{n+1}+(1+E)C_n=0.$$

The initial values $`C_0=1`$ and $`C_1=E`$ are therefore again sufficient to determine the correlation function for all values of the distance $`n`$.

For $`1/\sqrt{2}E0`$, which contains the range of final energies reached by the dynamics, $`\lambda _{3,4}=\mathrm{exp}(\pm \mathrm{i}Q)`$ are complex numbers with unit modulus, with

$$\mathrm{tan}Q=\frac{\sqrt{12E^2}}{E}(0Q\pi /2).$$
(3.7)
We are thus led to the expression

$$C_n=\left(\frac{1+E}{1E}\right)^{n/2}\left(\mathrm{cos}nQ+\frac{E^2}{\sqrt{12E^2}}\mathrm{sin}nQ\right),$$
(3.8)
for $`n0`$. In the a priori ensemble, the spin correlation function therefore exhibits an exponential decay, of the form $`\mathrm{exp}(n/\xi _{\mathrm{spin}})`$, modulated by oscillations at wavevector $`Q`$. Both the correlation length

$$\xi _{\mathrm{spin}}=\frac{2}{\mathrm{ln}{\displaystyle \frac{1E}{1+E}}}$$
(3.9)
and the wavevector $`Q`$, given by (3.7), depend continuously on energy.

We finally quote for further reference the value of the reduced susceptibility

$$\widehat{\chi }=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}C_n.$$
(3.10)
After some algebra, we obtain the simple result

$$\widehat{\chi }=\frac{E(1E)}{1+E}.$$
(3.11)
### 3.3 Distribution of domain sizes

Another characteristic feature of blocked states is the distribution of domain sizes $`f(\mathrm{})`$, defined as the probability that a given domain consists exactly of $`\mathrm{}`$ consecutive parallel spins. Since isolated spins are absent, domains have at least size two ($`\mathrm{}2`$).

In the a priori ensemble, the distribution $`f(\mathrm{})`$ can again be evaluated by the transfer-matrix method. Indeed $`\rho (\mathrm{})`$, the density (per unit length) of domains consisting of exactly $`\mathrm{}`$ spins, admits an expression similar to the middle side of (3.6), namely

$$\rho (\mathrm{})=\frac{L_1|𝒫\mathrm{e}^{(\mathrm{}2)\beta }|R_1}{\lambda _1^{\mathrm{}}}=L_1|𝒫|R_1\frac{1+E}{1E}\left(\frac{2E}{1E}\right)^\mathrm{}2,$$

where

$$𝒫=(++)+(++)=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\end{array}\right)$$

is the appropriate domain boundary operator. The explicit expression of the matrix element is not needed, as it can be fixed by normalization. The probability distribution $`f(\mathrm{})`$ of domain sizes indeed reads $`f(\mathrm{})=\rho (\mathrm{})/\rho `$, where $`\rho `$ is the total density of domains (or equivalently, of domain walls).

We thus obtain the geometric (i.e., discrete exponential) probability distribution

$$f(\mathrm{})=\frac{1+E}{1E}\left(\frac{2E}{1E}\right)^\mathrm{}2$$
(3.12)
for $`\mathrm{}2`$, and consistently

$$\rho =\underset{\mathrm{}=2}{\overset{\mathrm{}}{}}\rho (\mathrm{})=\frac{1}{\mathrm{}}=\frac{1+E}{2}=L_1|𝒫|R_1.$$
(3.13)
The characteristic length of (3.12),

$$\xi _{\mathrm{dom}}=\frac{1}{\mathrm{ln}{\displaystyle \frac{1E}{2E}}},$$
(3.14)
is in general different from the spin correlation length $`\xi _{\mathrm{spin}}`$ of (3.9), except for the typical value of energy (3.3), where $`\xi _{\mathrm{spin}}=\xi _{\mathrm{dom}}=1/S^{}`$.

The mean of the domain size distribution (3.12) agrees with (3.13), while its variance reads

$$Var\mathrm{}=\underset{\mathrm{}=2}{\overset{\mathrm{}}{}}\mathrm{}^2f(\mathrm{})\mathrm{}^2=\frac{2E(1E)}{(1+E)^2}.$$

Consider now a large sample of $`N`$ spins, in the canonical a priori ensemble (fixed parameter $`\beta `$). The number $`M`$ of domains in the sample is such that $`N=\mathrm{}_1+\mathrm{}+\mathrm{}_M`$ is the sum of $`M`$ independent variables distributed according to $`f(\mathrm{})`$, neglecting boundary effects. The expected number of domains and its variance are given by

$$\frac{M}{N}\underset{N\mathrm{}}{}\frac{1}{\mathrm{}}=\frac{1+E}{2},\frac{VarM}{N}\underset{N\mathrm{}}{}\frac{Var\mathrm{}}{\mathrm{}^3}=\frac{E(1E^2)}{4}.$$
(3.15)
One has therefore

$$NVarE=\frac{4VarM}{N}\underset{N\mathrm{}}{}\frac{4Var\mathrm{}}{\mathrm{}^3},$$
(3.16)
so that the second result of (3.15) agrees with expression (3.5) for $`NVarE`$.

## 4 Blocked states: dynamics

In this section we compare the predictions of the a priori approach, derived in Section 3, to the results of numerical simulations concerning the blocked configurations reached by the dynamics. We have used the rules of the zero-temperature Kawasaki dynamics summarized in Table 1, with random sequential updates, starting from a random non-magnetized initial configuration.

### 4.1 Mean energy

The first and the simplest quantity to be measured is the mean energy $`E`$ of the blocked configurations reached by the dynamics.

Let us start with a reminder of the kinetically constrained model $`(W_0=0)`$. For a random initial configuration, the two kinds of allowed defects, namely isolated unsatisfied bonds (domain walls) and domains of two unsatisfied bond (isolated spins) occur respectively with the following densities in blocked configurations :

$$q_1=\frac{1}{2}\left(1\mathrm{e}^{5/4}\mathrm{e}^{9/4}_1^{3/2}\mathrm{e}^{y^2}dy\right)=0.219704,q_2=\frac{\mathrm{e}^{5/4}}{4}=0.071626,$$

so that the mean energy of blocked configurations reads

$$E_0=1+2q_1+4q_2=\mathrm{e}^{9/4}_1^{3/2}\mathrm{e}^{y^2}dy=0.274087.$$
(4.1)
In the present case $`(W_00)`$, the mean energy $`E`$ of blocked states is expected to be below this number. Indeed the diffusive moves can only help relaxing more efficiently the energy excess of the disordered initial state.

Figure 3 shows a plot of the mean final energy $`E`$ against the diffusive rate $`W_0`$. Each data point corresponds to $`10^8`$ spins in total. We have checked that no appreciable size dependence is to be observed. The final energy is found to be well below (4.1), and below the typical a priori value (3.3). It only exhibits a very weak dependence on the rate $`W_0`$, increasing from the extrapolated minimum value

$$E(0)=0.5279$$
(4.2)
in the $`W_00`$ limit to the maximum value

$$E(1)=0.51633$$
(4.3)
for $`W_0=1`$.

The regime $`W_01`$ again deserves some more attention. The mean energy has a discontinuity at $`W_0=0`$. Indeed the $`W_00`$ limit of the energy, $`E(0)`$ given by (4.2), is different from that of the constrained dynamics, $`E_0`$ given by (4.1). This discontinuity can be analyzed as follows. In the diffusive part of the dynamics, free spins annihilate by meeting either each other or domain walls. The collision of a single free spin with a domain wall relaxes 4 units of energy, while the coalescence between $`n2`$ spins within a domain relaxes $`4(n1)`$ units of energy, i.e., $`4(n1)/n`$ per spin. If all the free spins were annihilated in meeting domain walls, the final energy would assume the value $`E=E_04p_2`$. The efficiency of the diffusive relaxation mechanism can therefore be characterized by the ratio

$$\eta =\frac{E_0E(0)}{4p_2}=1\underset{n2}{}\frac{\mathrm{\Pi }_n}{n},$$

where $`\mathrm{\Pi }_n`$ is the probability that a given free spin gets annihilated in a coalescence of $`n2`$ spins. The extrapolated value (4.2) of the final energy yields a rather high efficiency: $`\eta =0.886`$.

In view of the weak dependence of the final energy on the rate $`W_0`$ (see Figure 3), hereafter we restrict the numerical analysis to the case $`W_0=1`$ of zero-temperature Metropolis dynamics.

### 4.2 Higher cumulants of energy

We now turn to the statistics of the energy of the blocked configurations, besides its mean value studied above. In analogy with the a priori estimate (3.4), the final energy is expected to obey a large-deviation formula of the type

$$P(E)\mathrm{exp}(N\mathrm{\Sigma }(E)),$$

with $`\mathrm{\Sigma }(E)`$ being the dynamical entropy. This formula implies that the cumulants of the total energy scale as $`(NE)^kN`$, just as in usual equilibrium situations.

Instead of measuring the whole function $`\mathrm{\Sigma }(E)`$, which would require a quite extensive numerical effort, we have measured the first four cumulants of the energy of the blocked configurations. Table 2 gives the measured values of the scaled energy cumulants for $`W_0=1`$ and a random initial configuration. For comparison we also list the predictions of the full a priori ensemble (3.5) for $`E=E^{}`$ of (3.3), and of the restricted a priori ensemble, obtained by inserting into (3.5) the observed value (4.3) of the mean energy. The simulations have been performed on samples of various sizes ranging from $`N=50`$ to 200, having $`10^{10}`$ spins in total. No systematic size dependence is observed. Statistical errors can be estimated to be of the order of one unit of the least significant digit. Both a priori schemes perform very unequally in predicting the energy cumulants. The variance is rather accurately predicted by both schemes, which perform equally poorly for the third cumulant, while the restricted scheme is definitely better for the fourth cumulant.

### 4.3 Spin correlation function

We now turn to the spin correlation function $`C_n=\sigma _0\sigma _n`$ in the blocked configurations. This quantity has two remarkable properties. First, since isolated spins are not allowed in blocked configurations, we have $`C_1=12\rho =E`$ and $`C_2=14\rho `$, where $`\rho `$ is the density of domain walls (3.13), hence

$$C_2=2C_11.$$
(4.4)
Second, the total magnetization $`M=_n\sigma _n`$ is exactly conserved in Kawasaki dynamics. Furthermore, for a homogeneous non-magnetized state, we have $`M^2=N\widehat{\chi }`$, with the definition (3.10). The quantity $`\widehat{\chi }`$ therefore has the same value in the initial and final configurations, i.e., for a random initial state,

$$\widehat{\chi }=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}C_n=1.$$
(4.5)
Figure 4 shows a comparison between numerical data (circles and full line), corresponding to $`310^9`$ spins in total, and the prediction (3.8) of the restricted a priori ensemble at the observed mean energy (4.3) (triangles and dashed line). The a priori prediction and numerical data coincide both for $`n=1`$ (by construction, as $`E`$ has been imposed) and for $`n=2`$ (as a consequence of (4.4)). For larger distances, the spin correlations oscillate in sign and fall off exponentially, in qualitative agreement with the a priori prediction. Both the period of oscillations and the decay length are observed to be slightly larger than those of the a priori ensemble. These observations remain however at a qualitative level.

We also notice that the a priori ensemble fails to reproduce the identity (4.5). Indeed $`\widehat{\chi }`$ of (3.11) equals unity for $`E=1\sqrt{2}=0.414214`$, a value which neither agrees with the most probable a priori energy $`E^{}`$ of (3.3), nor with the observed mean energy (4.3). Conversely, the prediction for the correlation function of the restricted a priori ensemble at energy (4.3), plotted in Figure 4, has $`\widehat{\chi }=1.6192`$, which is significantly different from unity.

### 4.4 Distribution of domain sizes

We have determined the distribution of domain sizes $`f(\mathrm{})`$ by means of extensive numerical simulations, with a random initial state. Figure 5 shows our data corresponding to $`310^9`$ spins in total. The observed distribution is not exponential, at variance with the prediction (3.12) of the a priori ensemble. The full line suggests an exponential asymptotic fall-off of the distribution, with a characteristic length $`\xi _{\mathrm{dom}}=1.75`$. This length is close to the characteristic length of domains with a single diffusive spin, $`\xi _{\mathrm{dif}}=1.61`$, plotted in Figure 1. As a matter of fact, the equality $`\xi _{\mathrm{dom}}=\xi _{\mathrm{dif}}`$ is expected to hold in the $`W_00`$ limit. The prediction (3.12) of the restricted a priori ensemble at the observed mean energy (4.3) is shown as a straight dashed line, whose inverse slope is the prediction $`\xi _{\mathrm{dom}}=2.603`$ (3.14) of the a priori approach.

Another noticeable difference with the prediction of the a priori approach is that the sizes of successive domains in the final states of the dynamics are not independent. Indeed, if they were so, the relation (3.16) would hold, while the data yield $`4Var\mathrm{}/\mathrm{}^3=0.2790`$, a number significantly below the observed variance $`NVarE=0.3446`$, listed in Table 2. Hence domain sizes exhibit a weak but definite positive correlation. Domain sizes which are neither statistically independent nor exponentially distributed have also been observed recently  in spin chains undergoing tapping dynamics.

## 5 Features of persistence

We now turn to features related to the whole history of a single given spin, say $`\sigma _0`$, the spin situated at the origin. This kind of problems belongs to the realm of persistence. In the present context, the central quantity is $`\nu `$, the total number of times $`\sigma _0`$ flips during the history of the sample, before a blocked configuration is reached. The number $`\nu `$ of spin flips is finite with unit probability in the limit of a large system. Zero-temperature Kawasaki dynamics is therefore of type $``$ in the classification of . Furthermore $`\nu `$ is random, as it depends both on the initial configuration and on the history of the system. We are therefore interested in the distribution $`p_\nu `$ of the number of flips, which is expected to have a well-behaved limit when $`\sigma _0`$ is deep inside a large enough sample.

Before presenting numerical data, we first predict the main salient features of the distribution $`p_\nu `$, following the lines of Section 2. Consider again a large domain of $`L1`$ spins, containing the origin, on which a single free spin diffuses. The spin $`\sigma _0`$ flips twice each time the free spin traverses the origin. We are thus led to the following effective problem.

Consider a random walker in the interval $`L_1<n<L_2`$, with absorbing boundaries at $`n=L_1`$ and $`n=L_2`$. The walker starts from the origin ($`n=0`$). The probability that the walker returns to the origin before being absorbed by either boundary reads

$$P_{\mathrm{ret}}(L_1,L_2)=1\frac{1}{2}\left(\frac{1}{L_1}+\frac{1}{L_2}\right).$$
(5.1)
This formula relies on a well-known result in the gambler’s ruin problem . Indeed, suppose that the walker’s first jump is to the right (resp. to the left), and consider $`0nL_2`$ (resp. $`0nL_1`$) as the gambler’s wealth. Then the ruin probability reads $`P=11/L_2`$ (resp. $`P=11/L_1`$). The expression (5.1) is the arithmetical mean of both ruin probabilities.

The probability that $`\sigma _0`$ flips an even number $`\nu =2k1`$ of times therefore approximately reads $`p_{2k}(L_1,L_2)(P_{\mathrm{ret}}(L_1,L_2))^k(1P_{\mathrm{ret}}(L_1,L_2))`$. Hence the distribution of the number of flips can be estimated, in analogy with (2.1), as

$`p_{2k}`$ $``$ $`{\displaystyle \underset{L_1,L_2}{}}(L_1+L_21)f_{\mathrm{dif}}(L_1+L_21)P_{\mathrm{ret}}(L_1,L_2)^k(1P_{\mathrm{ret}}(L_1,L_2))`$
$``$ $`\left({\displaystyle _0^{\mathrm{}}}\mathrm{exp}\left({\displaystyle \frac{L}{\xi _{\mathrm{dif}}}}{\displaystyle \frac{k}{2L}}\right)dL\right)^2.`$
The integral entering this expression closely resembles that entering (2.1). For $`k1`$ it is legitimate to use the steepest-descent method. The saddle point lies at $`L_c=\sqrt{k\xi _{\mathrm{dif}}/2}`$. We thus obtain the stretched exponential law

$$p_{2k}\sqrt{k}\mathrm{exp}\left(2\sqrt{\frac{2k}{\xi _{\mathrm{dif}}}}\right)$$
(5.2)
for the distribution of the number of spin flips, provided $`\nu =2k`$ is a large even number.

The occurrence of odd numbers $`\nu =2k+1`$ of spin flips can also be explained in the above framework, if the spin $`\sigma _0`$ is either situated at an endpoint of a domain, or involved in a coalescence event between two free spins. Both effects are expected to scale as the inverse of the domain size $`L`$. This leads us to predict that odd values of $`\nu `$ are suppressed by a factor of order $`1/L_c`$, i.e.,

$$\frac{p_{2k+1}}{p_{2k}}\frac{a}{\sqrt{k\xi _{\mathrm{dif}}}}.$$
(5.3)
Figure 6 shows a logarithmic plot of the distribution $`p_\nu `$ of the number of flips. The simulations again concern samples of various sizes having $`210^{10}`$ spins in total, with $`W_0=1`$. Even numbers of spin flips $`\nu =2k`$ (full symbols) are clearly more frequent than odd numbers $`\nu =2k+1`$ (open symbols), especially for large values of $`k`$. From a quantitative viewpoint, the full lines on Figure 6 show a common fit of the numerical data for $`\nu >10`$ according to the asymptotic predictions (5.2), (5.3). We thus obtain $`2/\sqrt{\xi _{\mathrm{dif}}}1.55`$, in agreement with the data of Figure 1, $`\xi _{\mathrm{dif}}1.61`$, i.e., $`2/\sqrt{\xi _{\mathrm{dif}}}1.57`$. We also obtain $`a1.5`$, albeit with a large uncertainty.

Besides the above analysis of the regime of large numbers of spin flips, due to long surviving free spins, our data yield yet other interesting informations. First, the persistence probability, i.e., the probability for a spin to never flip, has the following value:

$$p_0=0.44739.$$

Then, the probability that a spin flips an even or an odd number of times reads

$$P_{\mathrm{even}}=\underset{k0}{}p_{2k}=0.72302,P_{\mathrm{odd}}=\underset{k0}{}p_{2k+1}=1P_{\mathrm{even}}=0.27698.$$

These figures can be related to the overlap between a random initial configuration and the corresponding final one. We have indeed

$$Q=\sigma _0(0)\sigma _0(\mathrm{})=\underset{\nu 0}{}(1)^\nu p_\nu =P_{\mathrm{even}}P_{\mathrm{odd}}=0.44604.$$

The mean number of spin flips reads

$$\nu =\underset{\nu 0}{}\nu p_\nu =2.06916.$$

This quantity can be used to determine the fractions $`f_{\mathrm{cond}}`$ of condensation moves and $`f_{\mathrm{dif}}`$ of diffusive moves in a typical history, starting from a random initial configuration. Indeed, on the one hand, a condensation move lowers the total energy by four units, while a diffusive moves leaves it unchanged (see Table 1). On the other hand, any move involves exactly two spin flips. We have therefore, using (4.3),

$$f_{\mathrm{cond}}=1f_{\mathrm{dif}}=\frac{E}{2\nu }=0.12477.$$

To close up, let us compare the above results to the case of constrained Kawasaki dynamics ($`W_0=0`$, where only condensation moves are allowed. In this situation, the above quantities can be simply related to the final energy $`E_0`$ of (4.1). Indeed every spin flips at most once, so that only $`p_0`$ and $`p_1=1p_0`$ are non-zero, and $`f_{\mathrm{cond}}=1`$. Starting again from a random configuration, we are left with

$$2(1p_0)=1Q=2\nu =E_0=0.274087.$$

## 6 Discussion

First, we wish to emphasize the richness of the zero-temperature dynamics of the ferromagnetic Ising chain. There are indeed four different natural kinds of dynamics, summarized in Table 3. Only Glauber dynamics gives rise to a bona fide coarsening dynamics, obeying dynamical scaling with a typical domain size growing as $`L(t)\sqrt{t}`$. With the three other dynamics, the system is left in a metastable configuration after a relatively short blocking time. The constrained Glauber and Kawasaki dynamics have been analyzed in our previous work . Each spin may flip at most once, before a global blocked state is reached. The blocking time scales with the number of spins as $`\mathrm{ln}N`$. These models can be exactly mapped onto the random sequential addition problem of dimers and hollow trimers, respectively, hence allowing an analytical treatment.

The last case not considered so far from the viewpoint of metastability, Kawasaki dynamics, has been the subject of this work. Because of the diffusive moves, the dynamics is only partly irreversible. This novel feature makes the Kawasaki problem both richer and hopefully closer to more realistic situations. The number of flips of a given spin, although finite with probability one, may be arbitrarily large. The blocking time grows as $`(\mathrm{ln}N)^3`$. On the other hand, analytical tools being no longer available, we needed to have recourse to numerical simulations.

The present work demonstrates that there is no good a priori statistical description of the metastable states reached by Kawasaki dynamics. In other words, the Edwards hypothesis is invalidated, at least as an exact prescription. Systematic differences are indeed observed between numerical results and a priori predictions, especially in the pattern of spin correlations (Figure 4) and in the distribution of domain sizes (Figure 5). The latter are neither statistically independent nor exponentially distributed.

We also want to underline that the effective trapping description of the late stages of the zero-temperature Kawasaki dynamics, already emphasized in Ref. , appears to yield quantitative predictions for several novel physical quantities, including the statistics of the blocking time ((2.4), (2.5), Figure 2), and the distribution of the number of spin flips ((5.2), (5.3), Figure 6).

#### Acknowledgements

Interesting discussions with Silvio Franz are gratefully acknowledged.