text:
# ON LOCAL APPROXIMATIONS TO THE NONLINEAR EVOLUTION OF LARGE-SCALE STRUCTURE

## 1 INTRODUCTION

Large-scale structures are believed to have formed from the gravitational amplification of primordial perturbations. At its first stages, the process of gravitational clustering can be investigated using linear perturbation theory. However, as the universe evolves, nonlinear concentrations of mass arise. Many structures we see today correspond to fluctuations several orders of magnitude higher than the mean density of the universe; for example, clusters of galaxies have typically $`\rho _{\mathrm{cluster}}/\rho _{\mathrm{universe}}10^210^3`$. For larger scales this ratio decreases, approaching unity in the largest structures.

As there is no analytical treatment for the nonlinear regime, $`N`$-body simulations are often resorted to. The numerical simulations had an enormous development in the last decade (see Bertschinger 1998 and references therein), being able to reproduce many features of the large scale structure. However they do not always provide a clear insight of the physics of nonlinear gravitational collapse. Moreover, they are usually very time-consuming, making it difficult to scan a large part of the parameter space of the cosmological models.

For this reason, semi-analytical methods have been devised to tackle such a complex problem. The first approximation developed to study the nonlinear regime was introduced by Zel’dovich (1970). There are now various approximation schemes to analyze different aspects of non-linear clustering, including extensions of the Zel’dovich approximation (for a review see Sahni & Coles 1995). Among them, the so-called local Lagrangian approximations have been introduced rather recently. The basic feature of these local approximations is that the kinematical parameters in each fluid element evolve independently of those of other elements. Thus the time evolution of a self gravitating fluid is replaced by a set of ordinary differential equations. This comes at the expense of losing information about the positions of each fluid element. Only local quantities, such as the density contrast, shear and expansion rate can be determined.

Due to their handy applicability compared to the numerical simulations, as seen in the case of the widely used Zel’dovich approximation, they deserve a closer investigation. For some of these methods, certain aspects of their performance and applicability have already been discussed. However, to the authors’ knowledge, no systematic comparison among them has ever been done. We consider it worthwhile to analyze them in a unified way in order to exploit general properties of these approximations, clarifying their similarities and differences. It is also important to compare their performance in some practical applications. In this paper, we discuss the following four approximations, in addition to the Zel’dovich approximation: the Local Tidal Approximation (Hui & Bertschinger 1996), the Deformation Tensor Approximation (Audit & Alimi 1996), the Complete Zel’dovich Approximation (Betancort-Rijo & López-Corredoira 2000) and the Modified Zel’dovich Approximation (Reisenegger & Miralda-Escudé 1995). All of them intend to be applicable to the highly non-linear regime. To the best of our knowledge, these comprise all existent local approximations in the literature, that are exact for planar, spherical, and cylindrical symmetries (except for the Zel’dovich approximation).

The paper is organized as follows: In section 2, we briefly review various local approximations in a unified way. In section 3, these methods are applied to several cases. First, we discuss the homogeneously collapsing ellipsoid. We then study their behavior under general initial conditions. Finally, we apply some of these approximations to the calculation of the cosmological mass function. We sum up our results and present conclusions in section 4. In the appendices we present useful formulae for the calculation of the mass function together with fitting formulae for the collapse time in the approximations considered here.

## 2 LOCAL APPROXIMATIONS

Throughout this paper we will only consider the case of cold dark matter (CDM), which is assumed to be collisionless, at least on large scales. This is well justified since 80 to 90% of the matter that clusters is composed by CDM (Turner 2000, Durrer & Novosyadlyj 2000). Furthermore, as long as the trajectories do not intersect, we can treat the CDM as a pressureless fluid.

We will be working in a matter-dominated flat universe (the Einstein-de Sitter universe, hereafter EdS). Recent observational evidences are consistent with a zero curvature universe (De Bernardis et al. 2000; Hanany et al 2000). Even if we had a non-flat universe we would only require that the curvature be negligible in the scales of interest. The assumption of matter dominance may seem unrealistic since the observations indicate that the universe is now dominated by a repulsive homogeneous cosmological term (Perlmuter et al. 1998; Riess et al. 1998; Zehavi & Dekel 1999). However, the energy density of this term decays more slowly than the matter density. In the case of a cosmological constant we would have $`\rho _\mathrm{\Lambda }=const`$. whereas for matter we have $`\rho _Ma^3`$, where $`a`$ is the scale factor of the universe. Since most structures form at a time when $`\rho _\mathrm{\Lambda }\rho _M`$, we can safely ignore the effect of the cosmological term on the collapse process.

The peculiar motions in the universe are much smaller than the speed of light. For perturbations on scales smaller than the Hubble radius, we can use the Newtonian approximation to describe the gravitational clustering. The basic equations for nonrelativistic pressureless matter in a perturbed EdS universe are the Euler, the continuity and the Poisson equations (Bertschinger 1996):

$$\frac{1}{a}\frac{dv_i}{d\tau }+\frac{\dot{a}}{a^2}v_i=\frac{\varphi }{x^i},$$
(1)
$$\frac{d\delta }{d\tau }+a\left(1+\delta \right)\theta =0,$$
(2)
$$\frac{^2\varphi }{x^ix_i}=4\pi Ga^2\overline{\rho }\delta ,$$
(3)
where $`\delta =\left(\rho \overline{\rho }\right)/\overline{\rho }`$ is the density contrast $`v_i=\left(dx_i/d\tau \right)/a`$ is the peculiar velocity, $`\theta =v^i/x^i`$ is the expansion, $`\varphi `$ is the peculiar gravitational potential, and the time variable $`\tau `$ is related to the cosmic time $`t`$ (also known as proper time) by $`d\tau =dt/a^2`$. The comoving coordinate $`x_i`$ is given in terms of the position $`r_i`$ by $`x_i=r_i/a`$. The left hand side of equation (1) is simply $`\left(d^2x_i/d\tau ^2\right)/a^2`$, so that it looks like the usual Euler equation (apart form the factor $`a^2`$). In an Einstein-de Sitter background the scale factor is proportional to $`\tau ^2`$. We set $`a=\tau ^2`$ such that $`4\pi Ga^2\overline{\rho }=6\tau ^2=6/a`$. The present value of the scale factor $`a_0`$ is fixed to be unity.

The Lagrangian coordinates $`q_i`$ are often used instead of the position $`x_i`$ in nonlinear analyses. In terms of $`q_i`$ the convective derivative $`d/d\tau =/\tau |_x+v_i/x_i`$ is simply given by the time derivative at fixed $`q`$: $`d/d\tau =/\tau |_q`$. The Lagrangian coordinates are chosen to be the initial comoving positions: $`q_i=lim_{a0}r_i/a`$.

The Jacobian matrix of the transformation $`x_iq_i`$

$$J_{ij}=\frac{x_i}{q^j},$$
(4)
is known as the deformation tensor. The velocity gradient $`v_i/x^j`$ can be expressed in terms of $`J_{ij}`$ as

$$\frac{v_i}{x^j}=\frac{1}{a}J_{kj}^1\frac{dJ_i^k}{d\tau }.$$
(5)
The density is given by $`\rho (x,t)=\overline{\rho }J`$, where $`J`$ is the determinant of $`J_{ij}`$. It is easy to see that the continuity equation (2) is solved exactly with $`\delta =J^11`$.

Differentiating equation (1) with respect to $`x_j`$ we find

$$J_{jk}^1\frac{d^2J_i^k}{d\tau ^2}\frac{1}{a^2}=\frac{\varphi }{x^ix^j},$$
(6)
whose trace furnishes Raychaudhuri’s equation

$$J_{ij}^1\frac{d^2J^{ji}}{d\tau ^2}=4\pi Ga^4\overline{\rho }\left(J^11\right).$$
(7)
This is a local equation for $`J_{ij}`$ in the sense that it has no spatial derivatives, although it is not sufficient for determining the nine components of the deformation tensor. Usually this equation is written in terms of the kinematical parameters, $`\theta ,`$ $`\sigma _{ij}`$ (shear) and $`\omega _i`$ (vorticity), defined by:

$$\frac{v_j}{x_i}=\frac{1}{3}\theta \delta _{ij}+\sigma _{ij}+\omega _{ij}\sigma _{ij}=\sigma _{ji},\omega _{ij}=ϵ_{ijk}\omega ^k=\omega _{ji}.$$
(8)
If the initial conditions have no vorticity, then we have $`\omega ^i=0`$ during all the evolution, as long as the trajectories do not intersect. Here we will consider only the case of vanishing vorticity.

Equations (1) to (3) form a set of nonlinear partial differential equations. However, for certain specific configurations the time evolution of the deformation tensor $`J_{ij}`$ behaves as if each space point evolves independently from the others. One might then expect that for more general situations the locality may hold, at least approximately, for these variables. Accordingly, several methods have been introduced which are known as local approximations. In their framework, the influence of the neighbors may enter only through the initial conditions.

In addition to the solution of the continuity and Euler equations, the local approximations discussed here will replace the essentially nonlocal exact equation (3) either by some Ansatz inspired on equation (7), or by local evolution equations for the second derivative of the peculiar gravitational potential $`\varphi `$ (see also Kofman & Pogosian 1995 for a discussion).

One of the basic features of local approximations is that the eigenvectors of the deformation tensor do not change with time. Thus, once diagonalized, $`J_{ij}`$ remains diagonal in the same frame, along all the evolution. This condition is either assumed from the beginning or appears as a consequence of the approximation introduced in the evolution equations. Actually, this assumption is not strictly consistent with the evolution of the mapping $`q_ix_i`$, so that the reconstruction of space coordinates in these local approximations is not possible (cf. subsection 2.6).

In the basis where $`J_{ij}`$ is diagonal,

$$J_{ij}=\left(1+w_i\right)\delta _{ij},$$
(9)
Raychaudhuri’s equation (7) is written as

$$\underset{i=1}{\overset{3}{}}\frac{\ddot{w}_i}{\left(1+w_i\right)}=4\pi Ga^4\overline{\rho }\left(\frac{1}{\left(1+w_1\right)\left(1+w_2\right)\left(1+w_3\right)}1\right).$$
(10)
The local approximations discussed here are required to be exact for planar, spherical and cylindrical symmetries. In the spherical case we have $`w_1=w_2=w_3`$; for a cylindrical perturbation $`w_1=w_2`$, and $`w_3=0;`$ whereas for planar symmetry $`w_2=w_3=0`$. In these three cases, as we have only one independent eigenvalue $`w_i`$ of the deformation tensor, this equation can be solved for $`w_i`$.

### 2.1 Zel’dovich Approximation

The Zel’dovich approximation (Zel’dovich 1970), hereafter ZA, can be viewed as a solution of the linearized form of equation (10):

$$\underset{i=1}{\overset{3}{}}\ddot{w}_i=4\pi Ga^4\overline{\rho }\underset{i=1}{\overset{3}{}}w_i.$$
(11)
Zel’dovich used the solution of the linearized equations (1 to 3) $`x_i=q_iD\left(t\right)\mathrm{\Psi }_i\left(q\right)`$ and extrapolated it into the nonlinear regime. The eigenvalues of the deformation tensor are thus given by

$$w_i=D\left(t\right)\lambda _i^0\left(q\right),$$
(12)
where $`\lambda _i^0`$ are the eigenvalues of$`\mathrm{\Psi }_i/x^j`$. Substituting this expression into the eq. (11) we find two solutions for $`D`$, known as the growing and decaying modes. For an EdS universe we have

$$w_i^+=^+\lambda _i^0a\text{ and }w_i^{}=^{}\lambda _i^0a^{3/2}.$$
(13)
Since the decaying mode becomes negligible very quickly, only the growing mode is relevant for our discussion. The initial conditions are specified in terms of the $`\lambda _i^0`$, which are functions of the initial positions $`q`$. The principal axes of $`\mathrm{\Psi }_i/x^j`$ are generally different for each point. We will denote the linear growing mode solution by $`\lambda _i`$:

$$\lambda _i=\lambda _i^0a.$$
(14)
In the linear regime the density contrast $`\delta `$ will be given by $`\delta _L=\left(\lambda _1^0+\lambda _2^0+\lambda _3^0\right)a=\delta _0a`$.

The gist of the Zel’dovich approximation is that the linearized trajectories can lead to nonlinear density perturbations. Analogous ideas have been applied in many approximations. An example is the higher-order Lagrangian expansions, where the perturbed quantity is the displacement field. In an EdS universe the solution may be written in the form $`x_i=q_i+_{n=1}a^n\mathrm{\Psi }_i^{(n)}\left(q\right).`$ The first order solution $`\mathrm{\Psi }_i^{(1)}`$ is the Zel’dovich approximation. The determination of the higher order $`\mathrm{\Psi }_i^{(n)}`$ follows from the lower order ones through the solution of Poisson equations. The second order solution is known as Post-Zel’dovich approximation (Moutarde at al. 1991; Buchert 1992; Lachièze-Rey 1993), and the third order is called Post-post-Zel’dovich (Juszkiewicz, Bouchet, & Colombi 1993; Buchert 1994).

The Zel’dovich approximation is widely used for the weakly non-linear regime, and for generating initial conditions for numerical simulations. It gives the exact solution for the case of planar symmetry.

### 2.2 Modified Zel’dovich Approximation

In the ZA the time factor in eq. (12) is independent of the initial conditions, and it is valid only for the linearized limit in $`w_i`$ (eq. 11). Reisenegger & Miralda-Escudé (1995) have proposed a generalization of the Zel’dovich approximation where $`D`$ may depend on the position through the initial conditions $`\lambda _i^0`$. The Ansatz $`w_i=D(\tau ,\lambda _i^0)\lambda _i^0\left(q\right)`$ is substituted in equation (10) to give

$$\frac{d^2D}{d\tau ^2}=4\pi G\overline{\rho }a^4\frac{\eta _1D\eta _2D^2+\eta _3D^3}{\eta _12\eta _2D+3\eta _3D^2},$$
(15)
where $`\eta _1=\lambda _1^0+\lambda _2^0+\lambda _3^0`$, $`\eta _2=\lambda _1^0\lambda _2^0+\lambda _1^0\lambda _3^0+\lambda _2^0\lambda _3^0`$ e $`\eta _3=\lambda _1^0\lambda _2^0\lambda _3^0`$. This equation, which must be solved numerically, determines completely the function $`D(\tau ,\lambda _i^0)`$ and defines the Modified Zel’dovich Approximation, hereafter MZA. It is exact for spherical, planar and cylindrical symmetries. However, for underdense regions ($`\delta _0<0`$), the MZA may not work, as pointed out by Reisenegger & Miralda-Escudé (1995). This is due to the fact that, when not all the three eigenvalues $`\lambda _i`$ have the same sign, the denominator in the right hand side of eq. (15) will eventually vanish. Thus MZA cannot be used with this kind of initial conditions.

### 2.3 Deformation Tensor Approximation

In the two local approximations discussed above, the time dependence of the three $`w_i`$ is the same and it can be completely determined from equation (10). The next two approximations will provide an equation for each of the three $`w_i`$ and an analytical solution of $`w_i`$ in terms of the linear solution $`\lambda _i`$ (eq. 14). Due to the symmetry among the axes, both the equation for $`w_i`$ and the explicit solution in terms of $`\lambda _i`$ should be invariant under any exchange of indices, $`(i,j,k)`$.

Equation (10) may be written in the form

$$\underset{i=1}{\overset{3}{}}\left[\left(1+w_j+w_k+w_jw_k\right)\ddot{w}_i4\pi Ga^4\overline{\rho }\left(1+\frac{w_j+w_k}{2}+\frac{w_jw_k}{3}\right)w_i\right]=0,$$
(16)
where $`(i,j,k)`$ is a permutation of $`(1,2,3)`$. Audit & Alimi (1996) have, as an Ansatz, split this equation into three equations for each $`w_i`$:

$$\left(1+w_j+w_k+w_jw_k\right)\ddot{w}_i=4\pi Ga^4\overline{\rho }\left(1+\frac{w_j+w_k}{2}+\frac{w_jw_k}{3}\right)w_i\text{.}$$
(17)
This equation defines the deformation tensor approximation, hereafter DTA. Another motivation for the above equation is that it is exact for planar, spherical and cylindrical perturbations. We have thus a set of local equations that allows to determine each $`w_i`$ completely. Of course this spliting of equation (16) is not unique and we could add more local terms in equation (17) which would obey the symmetry requirement.

### 2.4 Complete Zel’dovich Approximation

The Complete Zel’dovich Approximation, CZA (Betancort-Rijo & López-Corredoira 2000) assumes that the $`w_i`$ can be expanded in terms of the linear solution $`\lambda _i`$ (eq. 14). To satisfy the symmetries required, the power series must have the following expression:

$$r_i(\lambda _i,\lambda _j,\lambda _k)=1+\underset{l,m,n=0}{\overset{\mathrm{}}{}}C_{l,m,n}^p(\lambda _j+\lambda _k)^l(\lambda _j\lambda _k)^{2n}\lambda _i^m,$$
(18)
where

$$w_i=\lambda _ir_i,$$
(19)
and $`C_{l,m,n}^p`$ are the coefficients of the $`p`$-th order terms, with $`pl+2n+m`$. The Zel’dovich approximation corresponds to $`r_i=1`$. The second order term

$$w_i^{(2)}=\lambda _i\frac{3}{14}\left(\lambda _j+\lambda _k\right)$$
(20)
coincides with that of the DTA.

For planar configurations one should have $`r_i=1,`$ thus $`C_{0,m,0}^p=0`$. The other coefficients of the expansion are determined from equations (1) and (3) through a recursive scheme. Betancort-Rijo & López-Corredoira (2000) calculated explicitly the coefficients $`C_{l,m,n}^p`$ up to the terms of fourth order in $`\lambda `$ in an EdS universe.

When the higher-order terms become important, all of them contribute roughly the same. Thus Betancort-Rijo & López-Corredoira have chosen to truncate the series at the fourth order and approximate the rest by a function $`R(\lambda _i,\lambda _j,\lambda _k)`$. This function is parametrized in such a way that the result is in agreement with the exact planar, spherical and cylindrical dynamics. Their expression for $`R`$ is:

$`R(\lambda _i,\lambda _j,\lambda _k)`$ $`=`$ $`\left[19\left(\lambda _i{\displaystyle \frac{\lambda _j+\lambda _k}{2}}\right)\left(1{\displaystyle \frac{\lambda _i+\lambda _j+\lambda _k}{1.3}}\right)\right]`$ (21)
$`\left(R_{\mathrm{sp}}(\lambda _i+\lambda _j+\lambda _k)R_{\mathrm{sp}}(\lambda _i)+R_{\mathrm{sp}}\left({\displaystyle \frac{\lambda _j+\lambda _k}{2}}\right)\right),`$
where $`R_{\mathrm{sp}}`$ is the correction term $`R`$, corresponding to the spherical symmetry. By comparing the numerical results for overdense perturbations with the truncated series solution, they fitted $`R_{\mathrm{sp}}`$ as

$$R_{\mathrm{sp}}\left(x\right)=2.58\times 10^3x^5\left(1\frac{x}{2.06}\right)^1.$$
(22)
The expansion (18) up to $`p=4`$, together with expression (21), gives nearly exact results for spherical and cylindrical overdense perturbations, and the exact result in the planar case. Indeed, the CZA predicts that a spherical perturbation with $`\delta _0=1`$ will collapse at $`a_c=1.72`$, whereas the exact solution gives $`a_c=1.69`$.

The CZA does not apply for perturbations with negative values of $`\lambda _i^0`$. For example, when all the three $`\lambda _i^0`$ are negative the volume element should expand indefinitely, hence the $`\lambda _i`$ will approach infinity and therefore the series expansion breaks down. It can be easily seen that, as all $`C_{l,m,n}^p`$ are positive, if we truncate the series in an odd power of $`\lambda ,`$ $`r`$ will change sign, and the fluid element will eventually collapse.

### 2.5 Local Tidal Approximation

The general relativistic equations for the kinematical parameters in the projection formalism are very akin to the Newtonian ones (see Ellis 1973). The analog of the Poisson equation is obtained from the equation for the Weyl tensor. Barnes & Rowlinson (1989) pointed out that by neglecting the magnetic part of Weyl tensor $`H_{\mu \nu }`$, the evolution equation for the electric part $`E_{\mu \nu }`$ becomes local. The dynamics of kinematical parameters is then reduced to a closed set of local equations. This result was first applied to structure formation by Matarrese, Pantano, & Saez (1993). Since the magnetic part of the Weyl tensor has no Newtonian analogue, Bertschinger & Jain (1994) introduced the non-magnetic approximation, by simply discarding the magnetic part $`H_{\mu \nu }`$ in the equation for $`E_{\mu \nu }`$ in the application to the Newtonian cosmology. This approximation is exact for spherical and planar configurations, but fails for cylindrical symmetry. Also it was not able to reproduce the dynamics of the collapse even for a homogeneous ellipsoid. Thus we will not consider this approximation further in this work.

Bertschinger & Hamilton (1994) pointed out that, in a “Newtonian” limit, the role of the magnetic part $`H_{ij}`$ is not altogether negligible (see also Ellis & Dunsby 1997). Within this framework, Hui & Bertschinger (1996) have proposed the local tidal approximation (LTA), which consists in discarding some terms in the evolution equation for $`E_{ij}`$, to get

$$\frac{dE_{ij}}{d\tau }+\frac{1}{a}\frac{da}{d\tau }E_{ij}=4\pi Ga^3\overline{\rho }\sigma _{ij},$$
(23)
where $`E_{ij}`$ is the Newtonian limit of $`E_{\mu \nu }`$, which gives the tidal field:

$$E_{ij}=\frac{^2\varphi }{x^ix^j}\frac{1}{3}\frac{^2\varphi }{x^kx_k}\delta _{ij}=\frac{^2\varphi }{x^ix^j}\frac{4\pi Ga^2\overline{\rho }\delta }{3}\delta _{ij}.$$
(24)
Equations (23) and (6) written in terms of $`E_{ij}`$ form a closed set of local equations. It was shown that the LTA is exact for spherical, planar and cylindrical symmetries. In general, it is exact whenever the orientation and axis ratios of the gravitational and velocity equipotentials are equal and constant for the mass element under consideration (Hui & Bertschinger 1996).

It is possible to show that in the LTA, once the velocity gradient is diagonalized, it will remain diagonal (Hui & Bertschinger 1996) and so will the deformation tensor. We may write equation (23) in terms of $`w_i`$ by using the equation (24) together with equation (6). We will then have a set of three third-order equations for $`w_i`$ which completely determines their evolution, once appropriate initial conditions are provided. Alternatively, equation (23) can be solved in terms of the kinematical parameters. In this case, the $`w_i`$ are calculated by (eqs. (5) and (8))

$$\frac{dw_i}{d\tau }=a\left(\sigma _i+\frac{1}{3}\theta \right)\left(1+w_i\right),$$
(25)
where $`\sigma _i`$ are the eigenvalues of the shear $`\sigma _{ij}`$.

### 2.6 General Features

The local approximations discussed above are either a system of ordinary differential equations or explicit expressions in terms of the linear solution. In these approximations each point evolves independently of the others. The influence of the other fluid elements enters only through the initial conditions. They give the time evolution of the deformation tensor, and thus the kinematical parameters for each volume element.

These local approximations are exact under some geometrical symmetries. In particular, they are exact whenever $`w_1=w_2=w_3`$; or $`w_1=w_2`$, with $`w_3=0;`$ or $`w_2=w_3=0`$ are satisfied locally. They are nonperturbative, i.e., valid, in principle, for any $`\delta `$ or $`\lambda _i`$.

The second order solution of the CZA and DTA (eq. 20) is in agreement with the result of second order Lagrangian perturbation theory for the density contrast (Sahni & Coles 1995, Betancort-Rijo & López-Corredoira 2000). The ZA fails at second order.

The local approximations are not appropriate for recovering the positions. To see this, let us consider an initial configuration such that the deformation tensor is diagonal at every point

$$J_{ij}=\frac{x_i}{q^j}=\left(1+w_i\right)\delta _{ij}.$$
(26)
If this holds initially, it will be valid throughout the evolution, according to the local approximations. In this case $`x_i`$ would be given by

$$x_i=\left(1+w_i\right)𝑑q_i.$$
(27)
If $`w_i`$ had an explicit dependence on $`q_j`$ or $`q_k,`$ nondiagonal terms would arise in equation (26); hence $`w_i`$ must be a function of $`q_i`$ only. Consequently, for this particular choice of initial conditions, each $`\lambda _i^0`$ must depend only on the coordinate $`q_i`$. However, as the $`w_i`$ evolve, they will in general depend on the three $`\lambda _i^0`$ and, ultimately, on the three coordinates: $`w_i=w_i(\lambda _1^0\left(q_1\right),\lambda _2^0\left(q_2\right),\lambda _3^0\left(q_3\right),\tau )`$. Thus equation (26) can no longer be satisfied. This shows that the local approximations in general violate the integrability of the deformation tensor. In other words, we cannot recover the actual positions in the local approximations (except when $`w_i`$ is independent of the initial position).

Another way of seeing that the integrability is violated is as follows: If it were possible to reconstruct $`x_1`$ from $`J_{1j}`$, for example, then $`J_{1j}`$ ought to be a gradient field in lagrangian space. Therefore its curl should vanish. In the particular case of equation (26), this implies that $`w_1/q^2=0`$ and $`w_1/q^3=0`$. These conditions are satisfied, in general, only by the ZA for which $`w_i=a\lambda _i^0`$ (provided that $`\lambda _i^0=\lambda _i^0\left(q_i\right)`$). Thus, the only approximation which always permits the direct computation of the positions is the Zel’dovich approximation. In spite of the non-integrability, these methods offer an approximate solution for the deformation tensor, allowing to calculate local quantities, such as the kinematical parameters.

If any eigenvalues of the deformation tensor approach $`1`$, the density contrast $`\delta `$ will diverge. Since they are functions only of $`a`$ and the initial conditions, $`\lambda _i^0`$, we can expand them near the collapse time $`a_c`$, as

$$w_i=1\frac{dw_i}{da}|_{a=a_c}\left(a_ca\right)+\mathrm{}.$$
(28)
Therefore, the density contrast $`\delta `$ behaves, at the collapse time, as

$$\delta \left(a_ca\right)^\gamma ,$$
(29)
where $`\gamma `$ is the dimensionality of the collapse ($`\gamma =1`$ for the collapse in only one axis, $`\gamma =2`$ for the collapse in two axes simultaneously, and $`\gamma =3`$ for the collapse in three axes). On the other hand, for the expansion $`\theta `$, we have from equation (2)

$$\theta \frac{2\sqrt{a_c}\gamma }{\left(a_ca\right)}$$
(30)
for $`aa_c`$. The asymptotic behavior of the other kinematical parameters can also be determined in a similar fashion.

## 3 APPLICATIONS

In order to compare the performance of these approximation schemes, we apply them to some specific situations in the following subsections.

### 3.1 The homogeneous ellipsoid

An initially homogeneous ellipsoid in an expanding universe develops in such a way that the homogeneity is almost preserved during all the evolution. Therefore, the homogeneously collapsing ellipsoid model (HCE) is considered to be very accurate (Eisenstein & Loeb 1995; Hui & Bertschinger 1996). The result of local approximations have been compared to this model. Such a comparison is useful since it offers the possibility of checking these approximations in a less symmetrical situation (Hui & Bertschinger 1996; Audit & Alimi 1996; Betancort-Rijo & López-Corredoira 2000). It is worthwhile to compare theses analyses including the MZA.

The equation of motion for the HCE is given by (Icke 1973; White & Silk 1979)

$$\frac{d^2Y_i}{d\tau ^2}=\frac{2}{9}aY_i\left(X_1X_2X_3Y_1Y_2Y_3\right)C_D(Y_k^2,Y_j^2,Y_i^2),$$
(31)
where, as before, $`(i,j,k)`$ are permutations of $`(1,2,3)`$, $`Y_i`$ represent the axes of the ellipsoid in comoving coordinates, and $`X_i`$ are their asymptotic values for $`a0`$. The function $`C_D`$ is the degenerate case of Carlson’s integral of the third kind (Carlson 1977; Press et al. 1992):

$$C_D(x,y,z)=\frac{3}{2}_0^{\mathrm{}}\frac{ds}{\left(z+s\right)^{3/2}\sqrt{\left(x+s\right)\left(y+s\right)}}.$$

The linear growing mode is

$$Y_i\underset{a0}{}X_i\left(1\frac{1}{3}X_1X_2X_3C_D(X_k^2,X_j^2,X_i^2)\delta _0a\right).$$
(32)
As there is no rotation, the orientation of the principal axes does not change. Thus, the position of each element will be proportional to the expansion in each direction. If an element inside the ellipsoid has initial position $`q_i`$, then its coordinates at a later time will be given by

$$x_i=\frac{Y_i}{X_i}q_i.$$
(33)
With this expression we may compute the kinematical parameters which will not depend on the position. The same holds for the tidal field. Hence we can compute the evolution of a fluid element according to the local approximations, and compare with the evolution of $`E_{ij}`$, $`\sigma _{ij}`$, $`\delta `$ and $`\theta `$ as derived from the ellipsoid solution with the same initial conditions.

From expression (33) we can see that the deformation tensor does not depend on $`q_i`$ inside the ellipsoid:

$$J_{ij}=\frac{Y_i}{X_i}\delta _{ij}.$$
(34)
As discussed in the previous section, the positions of the fluid elements may not necessarily be recoverable in the local approximations. However, the choice of the same $`w_i`$ for any fluid element ($`w_i`$ independent of $`q_i`$) is consistent with the HCE. In this case, we may recover the positions from the $`w_i`$ as

$$x_i=\left(1+w_i\right)q_i.$$
(35)
In Figure 1 we compare the time evolution of the axes $`R_i=aY_i`$ of an ellipsoid in the five approximations discussed in this paper. Here, the initial values $`X_i`$ of the axes were arbitrarily chosen to be $`1:1.25:1.5`$ with $`\delta _0=1`$. The general conclusion does not depend substantially on the choice of these values, as will be seen in the next section. We see that the results of these approximations, except for the ZA, are very close to the one given by the homogeneous ellipsoid model. The ZA overestimates the collapse time, showing that a simple extrapolation of the linear trajectories underestimates the nonlinear effects. The common feature we observe in the local approximations is that the collapse occurs a little bit earlier in the directions of the two initially larger axes than the HCE case, whereas the collapse in the direction of the shortest axis is slightly delayed compared to the HCE. In other words, in the local approximations, the tidal forces are reduced compared with the HCE model.

Concerning the collapse time $`a_c`$, all these approximations give very similar results as shown in Table 1. The differences are less than $`5\%`$. The MZA gives the closest value to that of the HCE. Considering, however, that the HCE itself neglects the effect of the interaction of the background with the ellipsoid, this will not necessarily indicate that the MZA has the better performance among the other local approximations. In fact, for larger shear, the MZA deviates from the others as will be seen in the next section.

### 3.2 Generic Initial Conditions

Following Bertschinger & Jain (1994) we will parametrize the initial conditions in the following way

$$\lambda _i^0=\frac{2}{3}\epsilon _0Q_i\left(\alpha _0\right)+\frac{1}{3}\delta _0,$$
(36)
where $`Q_i\left(\alpha \right)`$ are the diagonal terms of the traceless quadrupole matrix

$$Q_{ij}\left(\alpha \right)=\mathrm{diag}[\mathrm{cos}\left(\frac{\alpha +2\pi }{3}\right),\mathrm{cos}\left(\frac{\alpha 2\pi }{3}\right),\mathrm{cos}\left(\frac{\alpha }{3}\right)].$$
(37)
It is easy to show that $`\epsilon _0`$ is related to the magnitude of the shear and tide, $`Q_i\left(\alpha _0\right)`$ gives ratios of the eigenvalues of $`E_{ij}`$ and $`\sigma _{ij}`$ (note that in the linear regime $`E_{ij}\sigma _{ij}`$), and $`\delta _0`$ the density contrast. The parameter $`\epsilon _0`$ varies from $`0`$ to $`\mathrm{}`$, $`\alpha _0`$ from $`0`$ to $`\pi `$, and $`\delta _0`$ can go from $`\mathrm{}`$ to $`+\mathrm{}`$. However, it is sufficient to study the dynamics for $`\delta _0=+1`$ and $`\delta _0=1`$, as we shall see below.

The initial perturbation $`\delta _0`$ is defined as the ratio between $`\delta `$ and the growth factor $`D`$ in the linear regime:

$$\delta _0=\underset{a0}{lim}\frac{\delta }{D}.$$
(38)
In an Einstein-de Sitter universe we have $`D=a`$. Thus choosing different values of $`\delta `$ is equivalent to rescaling $`a`$. This is so for all the kinematical parameters. Therefore, the equations of motion in the local approximations are invariant under the following scaling

$$\delta _0c\delta _0,\epsilon _0c\epsilon _0,\mathrm{and}ac^1a.$$
(39)
Due to this invariance we can express the collapse time $`a_c`$ as (Audit, Teyssier, & Alimi 1997):

$$a_c(\delta _0,\epsilon _0,\alpha _0)=\{\begin{array}{cc}\left|\delta _0^1\right|a_c^+(\epsilon _0/\delta _0,\alpha _0),\hfill & \text{if }\delta _0>0,\hfill \\ \left|\delta _0^1\right|a_c^{}(\epsilon _0/\left|\delta _0\right|,\alpha _0),\hfill & \text{if }\delta _0<0,\hfill \end{array}$$
(40)
where $`a_c^\pm (\epsilon _0/\delta _0,\alpha _0)=a_c(\pm 1,\epsilon _0/\left|\delta _0\right|,\alpha _0)`$. Hence we just need to compute the two functions $`a_c^+`$ and $`a_c^{}`$ which depend on $`\epsilon _0/\delta _0`$ and $`\alpha _0`$ only.

In Figures 2 and 4, we plot the collapse time $`a_c`$ as a function of $`\epsilon _0`$ and $`\alpha _0`$ for overdense and underdense perturbations, respectively. Since the MZA and CZA do not apply for some underdense regions, we have not displayed the results of these approximations in figure 4. We also show the signs of $`\lambda _i^0`$ corresponding to the initial conditions in these two figures.

The parameter space of initial conditions that can be spanned by an ellipsoid with any axes ratios is equivalent to having the three $`\lambda _i^0`$ positive. The region corresponding to the homogeneous ellipsoid is limited to relatively small shear, and all the local approximations, except for the ZA, agree significantly well in this region. The ZA overestimates the collapse time for spherical configurations. We see that they still quite similar for overdense perturbations in general. The MZA substantially deviates from the others for high shear. In all the cases the shear accelerates the collapse, which is a well-known nonlinear effect. Thus the first regions to collapse are not necessarily those with higher density. We can also see that oblate initial configurations (for which $`\mathrm{cos}\alpha _0>0`$) collapse first. Thus planar collapse is favored by these approximations.

For negative perturbations the difference among the approximations is enhanced. The LTA systematically gives slightly larger collapse times than the DTA. The collapse time given by the ZA is the shortest among the three. It is important to notice that in the local approximations underdense regions may also collapse, due to the effects of the shear.

The relevance of the shear in the nonlinear phase of gravitational clustering is in agreement with $`N`$-body simulations (Katz, Quinn, & Gelb 1993), yet it is sometimes ignored in structure formation studies. Any model based on the spherical collapse would miss this effect.

An interesting aspect of the local approximations is that the collapse time has the same asymptotic behavior

$$a_c(\delta _0,\epsilon _0,\alpha _0)\frac{C}{\epsilon _0}$$
(41)
for high initial shear ($`\epsilon _0\delta _0`$) in all the approximations, where $`C`$ is a (slowly varying) function of $`\alpha _0`$ only (see Appendix A). That is, the collapse time for large shear does not depend on $`\delta _0`$ and is inversely proportional to the initial shear $`\epsilon _0`$.

### 3.3 The Cosmological Mass Function

The mass function $`n\left(M\right)`$ is defined such that $`n\left(M\right)dM`$ gives the number density of collapsed dark matter clumps with masses between $`M`$ and $`M+dM`$. These clumps are associated with proto-galactic haloes, and with galaxy groups and clusters. Comparing theoretical mass functions with observations provides important constraints on the cosmological parameters (Bahcall & Cen 1993; Girardi et al. 1998; Rahman & Shandarin 2000) and the spectrum of primordial perturbations (Lucchin & Matarrese 1988; Ribeiro, Wuensche, & Letelier 2000). The approach of Press & Schechter (1974) to calculate the mass function, hereafter PS, has been extended to nonspherical collapse and applied to some local approximations (Monaco 1995 for the ZA, and Audit et al. 1997 for the DTA). Here, we extend such analysis to the LTA and compare them.

Let $`F(M;a_0)`$ be the fraction of collapsed objects at $`a_0`$ with mass higher than $`M`$; then the mass function is given by

$$n\left(M\right)=\frac{\overline{\rho }}{M}\frac{dF}{dM}.$$
(42)
The fraction $`F`$ may be calculated as an integral over all the possible initial conditions weighted by their probabilities:

$$F=\frac{1}{F_0}_0^\pi _0^{\mathrm{}}_{\mathrm{}}^+\mathrm{}s(a_0;\delta _0,\epsilon _0,\alpha _0)P_M(\delta _0,\epsilon _0,\alpha _0)𝑑\delta _0𝑑\epsilon _0𝑑\alpha _0.$$
(43)
The function $`s`$ is equal to one if an element with parameters $`\delta _0,\epsilon _0,\alpha _0`$ has already collapsed at $`a_0`$, and is zero otherwise; $`F_0`$ is a normalization factor. The collapse time $`a_c`$ of a fluid element with initial perturbations parametrized by $`\delta _0,\epsilon _0,\alpha _0`$ can be computed in the local approximations. As mentioned in subsection 2.6, the collapse is characterized by the divergence of the density, which is equivalent to the first axis collapse. Beyond this point the Lagrangian formalism breaks down. Some authors (Audit et al. 1997; Lee & Shandarin 1998; Sheth, Mo & Tormen 1999) have suggested other alternatives for the definition of collapse in the calculation of the mass function. Here we will prefer to keep the simplest assumption of first axis collapse, since it does not introduce any free parameter.

What we need now is the probability distribution function $`P_M(\delta _0,\epsilon _0,\alpha _0)`$ for the initial conditions. Assuming Gaussian initial fluctuations, Doroshkevich (1970) derived the joint probability for the three eigenvalues of the deformation tensor $`\lambda _1^0,\lambda _2^0`$ and $`\lambda _3^0`$. Using this result, $`P_M(\delta _0,\epsilon _0,\alpha _0)`$ is given by the product of three independent probabilities for each parameter $`\delta _0`$, $`\epsilon _0`$ and $`\alpha _0`$:

$`P_\nu \left(\delta _0\right)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi \mathrm{\Delta }^2}}}\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\delta _0}{\mathrm{\Delta }}}\right)^2\right],`$ (44)
$`P_\chi \left(\epsilon _0\right)`$ $`=`$ $`{\displaystyle \frac{50}{3}}\sqrt{{\displaystyle \frac{5}{2\pi \mathrm{\Delta }^2}}}\left({\displaystyle \frac{\epsilon _0}{\mathrm{\Delta }}}\right)^4\mathrm{exp}\left[{\displaystyle \frac{5}{2}}\left({\displaystyle \frac{\epsilon _0}{\mathrm{\Delta }}}\right)^2\right],`$ (45)
$`P_\alpha \left(\alpha _0\right)`$ $`=`$ $`\mathrm{sin}\left({\displaystyle \frac{\alpha _0}{3}}\right)\left[{\displaystyle \frac{3}{2}}2\mathrm{sin}^2\left({\displaystyle \frac{\alpha _0}{3}}\right)\right].`$ (46)
The variance $`\mathrm{\Delta }`$ is related to the mass $`M`$ and the power spectrum of the primordial density field $`\sigma _k`$ through

$$\mathrm{\Delta }^2\left(R\right)=_0^{\mathrm{}}\frac{dk}{2\pi ^2}W_k^2\left(R\right)k^2\sigma _k^2,$$
(47)
where $`M=(4\pi /3)f_WR^3\overline{\rho }`$ and $`W_k\left(R\right)`$ is the Fourier transform of a filter with width $`R`$ in physical space. The factor $`f_W`$ depends on the shape of the filter function, for a top-hat filter we have $`f_W=1`$, whereas for a sharp$`k`$ filter $`f_W=9\pi /2`$. The mass function can now be written in the form

$$n\left(M\right)=\frac{\overline{\rho }}{M}\frac{d\mathrm{\Delta }}{dM}\mathrm{\Phi }\left(\mathrm{\Delta }\right),$$
(48)
where $`\mathrm{\Phi }\left(\mathrm{\Delta }\right)=dF(\mathrm{\Delta })/d\mathrm{\Delta }`$. The function $`\mathrm{\Phi }\left(\mathrm{\Delta }\right)`$ contains all the influence of the dynamics and depends neither on the particular form of the power spectrum nor of the filter $`W`$; it is referred to as the universal mass function (Audit et al. 1997).

We calculate the universal mass function for the ZA, LTA and DTA but not for the MZA and CZA since they do not apply for negative density perturbations. In appendix B we show the detailed calculation.

In Figure 5, we show the mass functions for those approximations. For comparison, we also display in this figure the fit to $`N`$-body simulations obtained by Jenkins et al. (2000), together with the standard PS mass function.

We see that the results of the DTA and LTA are very similar. Furthermore, in the high-mass tail ($`\mathrm{\Delta }0.5`$), they reproduce well the results of the $`N`$-body simulations. However, we can see that these approximations overestimate the concentration of masses near $`\mathrm{\Delta }=1`$. The right-end tail of the distribution decays more rapidly compared to the $`N`$-body simulations. This tendency is still enhanced in the ZA. However, in these approximations, the position of the maximum of the distribution is close to that of the $`N`$-body simulations, giving a better estimate than that of the PS; in particular, the LTA and DTA give nearly the same value as the $`N`$-body results.

As for the normalization factor $`F_0`$, there exists an extensive discussion on its origin (see, for example, Peacock & Heavens 1990; Bond et al. 1991; Jedamzik 1995; Yano, Nagashima, & Gouda 1996). The normalization factors for the local approximations are close to one ($`1/0.92`$ for the DTA, and $`1/0.89`$ for the LTA) whereas in the original PS derivation the normalization factor needed is $`F_0=1/0.5`$. This is due to the fact that, in the spherical collapse model, only overdense regions collapse.

The fact that around $`\mathrm{\Delta }=1`$ the number of objects is overestimated in the local approximations implies, due to the normalization of the mass function, that they should provide a lower estimate than the $`N`$-body simulations for large enough $`\mathrm{\Delta }`$. There, the contribution from the low-mass objects is dominant; in any realistic process, they may also arise from the fragmentation of larger clusters. The criterion for the formation of a clump from the direct collapse of an initially perturbed region does not account for these complex processes of fragmentation. Therefore, the discrepancy of the mass function for high $`\mathrm{\Delta }`$ might be attributed to the use of expression (43) rather than to the definition (42).

## 4 DISCUSSION

We have investigated local Lagrangian approximations to the nonlinear dynamics of pressureless dark matter. We have selected the modified Zel’dovich approximation (MZA), the deformation tensor approximation (DTA), the complete Zel’dovich approximation (CZA), and the local tidal approximation (LTA), in addition to the original Zel’dovich approximation (ZA). These four approximations were designed to improve the ZA, and are in fact exact for planar, spherical and cylindrical symmetries, whereas the ZA is only exact for the planar case. They are semi-analytic and easy to be implemented in any application where local quantities are involved, such as the calculation of the mass function in the PS approach.

All the local approximations discussed here, except for the ZA, provide quite a similar evolution for an ellipsoid, reproducing the results of the homogeneous ellipsoid model. Thus, for this kind of positive density perturbations, these methods work fairly well. However, the MZA turns out to deviate substantially for large values of the shear as was shown in section (3.2), reflecting the fact that it does not give the correct second order solution. Furthermore, the MZA cannot deal with initially underdense regions that will eventually collapse. Therefore, its applicability is rather limited when compared to the other approximations.

We note that the second order expansions of the CZA, LTA and DTA coincide with the second order Lagrangian perturbation theory, whereas those of the MZA and ZA do not. It is interesting to recall that the LTA and DTA have very different origins from the CZA, but still they give the correct second order result.

The CZA, LTA and DTA give quite analogous results for generic initial conditions. However, at least in its original form, the CZA cannot be used for negative values of $`\lambda _i`$. One possible solution to this problem might be achieved through an expansion such as

$$r_i=\frac{1+_{l,m,n=0}^NE_{l,m,n}^p(\lambda _j+\lambda _k)^l(\lambda _j\lambda _k)^{2n}\lambda _i^m}{1+_{l,m,n=0}^MD_{l,m,n}^p(\lambda _j+\lambda _k)^l(\lambda _j\lambda _k)^{2n}\lambda _i^m}.$$
(49)
The coefficients $`E`$ and $`D`$ should be appropriately chosen to adjust the asymptotic behavior; in particular, we could use the numerical solution for underdense cylindrical and spherical perturbations to fit some of these coefficients, as done for the overdense case in the CZA. Besides, to agree with the perturbative solution (18) we should have

$$E_{l,m,n}^1D_{l,m,n}^1=C_{l,m,n}^1.$$
(50)
For higher orders, the determination of these coefficients is rather complicated. Further investigations on this possibility should be pursued.

Concerning the mass function, it is found that the LTA and DTA give an accurate result for large masses as compared to the $`N`$-body simulations. The position of the peak is also in good agreement, whereas its amplitude is overestimated by a factor 2 compared to the $`N`$-body results. Since the mass function is normalized to unity, this means that the local approximations, together with the PS formalism, underestimates the density of low-mass clusters. However, this might be a consequence of the criterion for the formation of a collapsed object based only on the collapse time.

It is interesting to notice that the collapse time, as a function of $`\alpha _0`$ and $`\epsilon _0,`$ has an approximate scaling property (see eq. A1), which is very precise for $`\epsilon _0/\delta _01`$. We conclude that this may be a general feature of gravitational collapse in local approximations, whose validity is worth checking in a more general setting.

While there is still not a clear theoretical understanding or support to the local approximations, they proved to be very accurate in the situations investigated here. They reproduce some well known features of nonspherical collapse, such as the possibility of collapse of some initially underdense regions, and the fact that the shear accelerates the collapse (see Sahni & Coles 1995). The main limitation of these approximations is that they only provide information about the internal state of a given mass element, but do not determine its position. Even so, their simplicity is highly expedient for practical applications, such as the calculation of nonlinear corrections to the microwave background anisotropies, and the Gunn-Peterson effect. In particular, they are suitable for obtaining statistical properties of the present fields as a function of the primordial ones, as in the case of the mass function. Further studies on the validity of these approximations, based on a comparison to $`N`$-body simulations, are required. Such a comparison would allow to fully test the approximations described in this paper and more generically, the locality hypothesis. If they still provide accurate results in this case, the local approximations could represent good alternatives to the computer simulations, taking much less computational time, allowing thereby a larger scanning of initial conditions. They could give complementary information to the $`N`$-body simulations and would provide a better physical understanding of the nonlinear dynamics of self gravitating systems.

Most of the results of this paper may be extended to more general backgrounds. The influence of any smooth component only alters the behavior of $`a(t)`$, and of the growing mode growth factor $`D\left(t\right)`$, that will not be equal any more. It would be interesting to study the relativistic analogue of the LTA. Another interesting extension would be to include vorticity in the local approximations as was done for the ZA in (Buchert 1992; Barrow & Saich 1993). We could use them to test the effects of a possible primeval vorticity on large scales (Li Xin-Li 1998).

MM acknowledges the participants of the “Pequeno Seminário” at CBPF, and was partially supported by a CNPq fellowship (contract no 142338/97-4). TK would like to acknowledge the partial support of the Brazilian sponsoring agencies CNPq (contract no 300962/86-0) and FAPERJ (contract no E-26/150.942/99), and a PRONEX grant (contract no 41.96.0886.00). M. O. C. acknowledges financial support from FUJB-UFRJ.

## Appendix A FITTING FORMULAE FOR THE COLLAPSE TIME

In order to avoid repeated numerical integrations of differential equations in the LTA and DTA, we have parametrized the collapse time as a function of initial conditions. For both of these local approximations, the following scaling property is approximately satisfied:

$$a_c^\pm (x,\alpha _0)H^\pm \left(f\right),\mathrm{with}f=xg^\pm \left(\alpha _0\right),$$
(A1)
where $`H^\pm `$ and $`g^\pm `$ are functions to be fitted for each approximation, and $`x=\epsilon _0/\left|\delta _0\right|`$. This relation becomes more accurate for increasing $`x`$. In any case the error of the fit is less than a few percent.

For $`g\left(\alpha \right)`$ we found that a kind of truncated Fourier series can be used to a very good approximation:

$`g\left(\alpha \right)`$ $`=`$ $`c_1\mathrm{cos}\left({\displaystyle \frac{\alpha }{2}}\right)+c_2\mathrm{cos}\alpha +c_3\mathrm{cos}\left({\displaystyle \frac{3\alpha }{2}}\right)+c_4\mathrm{cos}\left(2\alpha \right)+`$ (A2)
$`c_5\mathrm{cos}\left({\displaystyle \frac{5\alpha }{2}}\right)+\left(1c_1c_2c_3c_4c_5\right)\mathrm{cos}\left(3\alpha \right).`$
The values of the parameters for each case are shown in Table 2.

For the dependence of $`a_c`$ on $`x`$ for fixed $`\alpha _0`$, we fitted the function $`f=H^1`$ rather than $`H`$ itself, because this is the quantity we will need to compute the mass function. Taking into account the boundary value and asymptotic behavior, we parametrize $`f`$ by

$$f^+=\frac{d_1z\left(1+d_2z+d_4z^2\right)}{\left(1+d_3z+d_5z^2\right)a_c},\mathrm{with}z=d_0a_c,$$
(A3)
for overdense regions ($`\delta _0>0`$), and

$$f^{}=d_0+\frac{d_1a_c^{1/2}\left(1+d_2a_c^1+d_4a_c^{5/2}\right)}{1+d_3a_c^1+d_5a_c^2},$$
(A4)
for underdense regions ($`\delta _0<0`$). The parameters for the LTA and DTA are given in Table 3.

Notice that for high shear ($`x1`$) we have $`a_c1`$ such that:

$$a_c^+\underset{x\mathrm{}}{}\frac{d_1d_0\left(1+d_2d_0+d_4d_0^2\right)}{\left(1+d_3d_0+d_5d_0^2\right)}\frac{1}{xg^+\left(\alpha _0\right)}.$$
(A5)
and

$$a_c^{}\underset{x\mathrm{}}{}\frac{d_1d_4}{d_5}\frac{1}{xg^{}\left(\alpha _0\right)},$$
(A6)
In the Zel’dovich approximation we have an analytical expression for $`a_c`$. As the collapse occurs when the greatest $`\lambda _i`$ reaches the value $`1`$, using the parametrization (36) we get:

$$a_c^{ZA}=\frac{3}{\delta _0+2\epsilon _0\mathrm{cos}\left(\alpha _0/3\right)}.$$
(A7)
In this case we can clearly see the features of $`a_c`$:

$`a_c^{ZA}(\delta _0,\epsilon _0,\alpha _0)`$ $`=`$ $`{\displaystyle \frac{1}{\left|\delta _0\right|}}{\displaystyle \frac{3}{\pm 1+2\left(\epsilon _0/\left|\delta _0\right|\right)\mathrm{cos}\left(\alpha _0/3\right)}}`$ (A8)
$`=`$ $`{\displaystyle \frac{1}{\left|\delta _0\right|}}a_c^{ZA}(\pm 1,\epsilon _0/\left|\delta _0\right|,\alpha _0).`$
Note that the property (A1) is satisfied exactly for the ZA.

## Appendix B CALCULATION OF THE MASS FUNCTION

With the integral (43) we may write the universal mass function $`\mathrm{\Phi }\left(\mathrm{\Delta }\right)=dF/d\mathrm{\Delta }`$ in the form

$$\mathrm{\Phi }\left(\mathrm{\Delta }\right)=\frac{1}{F_0}\frac{d}{d\mathrm{\Delta }}_0^\pi _0^{\mathrm{}}_{\mathrm{}}^+\mathrm{}s(\nu \mathrm{\Delta },\chi \mathrm{\Delta },\alpha _0)P(\nu ,\chi ,\alpha _0)𝑑\nu 𝑑\chi 𝑑\alpha _0,$$
(B1)
where $`\nu =\delta _0/\mathrm{\Delta }`$ and $`\chi =\epsilon _0/\mathrm{\Delta }`$. With these new variables the dependence on $`\mathrm{\Delta }`$ will be present only in the function $`s`$, which may be written as

$$s=\mathrm{\Theta }\left(1a_c(\delta _0,\epsilon _0,\alpha _0)\right).$$
(B2)
Note that in the case of the Press & Schechter original approach this function is given by: $`s=\mathrm{\Theta }\left(\delta _0\delta _c\right),`$ where $`\delta _c=1.686`$ is the value at which a spherical perturbation collapses at $`a=1`$.

To calculate $`ds/d\mathrm{\Delta }`$ one uses the relation (40), obtaining:

$$\frac{ds}{d\mathrm{\Delta }}=\delta _D\left(1a_c(\nu \mathrm{\Delta },\chi \mathrm{\Delta },\alpha _0)\right)\frac{a_c^\pm (\chi /\nu ,\alpha _0)}{\left|\nu \right|\mathrm{\Delta }^2},$$
(B3)
where $`\delta _D`$ is the Dirac delta function. Therefore we may eliminate one of the integrals in expression (B1), with the mass function being calculated over the surface $`a_c(\nu \mathrm{\Delta },\chi \mathrm{\Delta },\alpha _0)=\left|\nu \mathrm{\Delta }\right|^1a_c^\pm (\chi /\nu ,\alpha _0)`$ $`=1`$. For this sake we need to write one of the three variables in terms of the others on this surface.

Let us assume the we have $`\chi `$ as a functions of $`\nu \mathrm{\Delta }`$ and $`\alpha `$: $`\chi =\chi _a(\nu \mathrm{\Delta },\alpha )`$, where the subscript $`a`$ indicates that $`\chi `$ is calculated over the surface $`a_c=1`$. The integral in $`\chi `$ in equation (B1) is thus eliminated using the relation

$$\delta _D\left(1a_c(\nu \mathrm{\Delta },\chi \mathrm{\Delta },\alpha _0)\right)=\delta _D\left(\chi \chi _a\right)\left|\frac{a_c}{\chi }\right|^1.$$
(B4)
The mass function (B1) is now given by:

$$\mathrm{\Phi }\left(\mathrm{\Delta }\right)=\frac{1}{F_0}_0^\pi _{\mathrm{}}^+\mathrm{}\frac{1}{\mathrm{\Delta }}\left|\frac{a_c}{\chi }\right|^1P_\nu \left(\nu \right)P_\chi \left(\chi _a\right)P_\alpha \left(\alpha _0\right)𝑑\nu 𝑑\alpha _0,$$
(B5)
where $`a_c/\chi `$ is evaluated in $`\chi _a(\nu \mathrm{\Delta },\alpha ).`$ We can simplify this expression further if the collapse time $`a_c(x,\alpha )`$ is only a function of the product $`xg\left(\alpha \right),`$ which we have seen is an excellent approximation for the local approximations studied here (see appendix A). Using the property (A1) we have

$$\chi _a=\frac{\nu }{g\left(\alpha _0\right)}H^1\left(\left|\nu \mathrm{\Delta }\right|\right),$$
(B6)
and

$$\left|\frac{a_c}{\chi }\right|=\left|\frac{1}{\nu \mathrm{\Delta }}\frac{H}{f}\frac{g\left(\alpha _0\right)}{\nu }\right|,$$
(B7)
where the superscript $`+`$ is implied for positive $`\nu `$, and the $``$ for negative $`\nu `$. Replacing these results in expression (B5) we get finally:

$$\mathrm{\Phi }\left(\mathrm{\Delta }\right)=\frac{1}{F_0}_{\mathrm{}}^+\mathrm{}_0^\pi \nu ^2\left|\left(\frac{H}{f}\right)^1\frac{1}{g\left(\alpha _0\right)}\right|P_\nu \left(\nu \right)P_\chi \left(\chi _a\right)P_\alpha \left(\alpha _0\right)𝑑\alpha _0𝑑\nu ,$$
(B8)
where $`\chi _a=\left|\nu \right|f_a/g\left(\alpha _0\right)`$, with $`f_a\left(\nu \mathrm{\Delta }\right)=H^1\left(\left|\nu \mathrm{\Delta }\right|\right)`$ and $`H/f`$, is calculated in $`f_a`$. This is why we have chosen to fit the function $`H^1`$, instead of its inverse.

As $`P\left(\nu \right)`$ and $`H/f`$ are independent of $`\alpha _0,`$ we integrate first in this variable:

$$I_1(\nu ,\mathrm{\Delta }):=_0^\pi \frac{1}{g^5\left(\alpha _0\right)}\mathrm{exp}\left[\frac{5}{2}\left(\frac{\nu f}{g\left(\alpha _0\right)}\right)^2\right]P_\alpha \left(\alpha _0\right)𝑑\alpha _0.$$
(B9)
The universal mass function will now be given by

$$\mathrm{\Phi }\left(\mathrm{\Delta }\right)=\frac{1}{F_0}N_{\mathrm{}}^+\mathrm{}I_1(\nu ,\mathrm{\Delta })\nu ^6f_a^4\left|\left(\frac{H}{f}|_{f_a}\right)^1\right|\mathrm{exp}(\frac{\nu ^2}{2})d\nu .$$
(B10)
where $`N=50\sqrt{5}/6\pi `$ is the product of the normalizations for $`P_\nu `$ and $`P_\chi `$ (equations (44) and (45)). Note that the above integral is limited for positive values of $`\nu `$, as $`H^1\left(\left|\nu \mathrm{\Delta }\right|\right)=0`$ for $`\nu \mathrm{\Delta }>f_0`$ ($`f_0=1.686`$ for the DTA and LTA, and $`f_0=3`$, for the ZA). Although our fitting formulas (sec. A) become less accurate for $`\nu \mathrm{}`$ and $`\nu f_0/\mathrm{\Delta },`$ the mass function is not affected, since the integrand in (B10) goes to zero in these regions. Here it is clear that the underdense regions do contribute to the mass function, as pointed out by Audit et al. (1997).

As an example, let us consider the Zel’dovich approximation. In this case we have

$$H_\pm \left(f\right)=\frac{3}{\pm 1+2f},\mathrm{and}g\left(\alpha \right)=\mathrm{cos}\left(\frac{\alpha }{3}\right).$$
(B11)
Using the following transformation of variables

$$x=\mathrm{sin}\left(\frac{\alpha }{3}\right)P_\alpha d\alpha =x\left[\frac{3}{2}2x^2\right]3\frac{dx}{\sqrt{1x^2}},$$
(B12)
we find an analytical expression for the integral (B9):

$$I_1:=\frac{1}{\left(\nu f\right)^2}\mathrm{exp}\left[\frac{5}{2}\left(\nu f\right)^2\right]\left[\frac{3}{25}\frac{1}{\left(\nu f\right)^2}\left(\mathrm{exp}\left[\frac{15}{2}\left(\nu f\right)^2\right]1\right)+\frac{9}{10}\right],$$
(B13)
where $`f^\pm =\left(3/\left(\left|\nu \right|\mathrm{\Delta }\right)1\right)/2`$. The mass function will be given by:

$`\mathrm{\Phi }\left(\mathrm{\Delta }\right)`$ $`=`$ $`{\displaystyle \frac{1}{F_0}}{\displaystyle \frac{15\sqrt{5}}{8\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{1}{\mathrm{\Delta }^4}}\mathrm{exp}({\displaystyle \frac{\nu ^2}{2}})\mathrm{exp}[{\displaystyle \frac{5}{8}}{\displaystyle \frac{\left(3\nu \mathrm{\Delta }\right)^2}{\mathrm{\Delta }^2}}]\times `$ (B14)
$`\left[{\displaystyle \frac{12}{5}}\mathrm{\Delta }^2\left(\mathrm{exp}\left[{\displaystyle \frac{15}{8}}{\displaystyle \frac{\left(3\nu \mathrm{\Delta }\right)^2}{\mathrm{\Delta }^2}}\right]1\right)+{\displaystyle \frac{9}{2}}\left(3\nu \mathrm{\Delta }\right)^2\right]d\nu .`$


format: text

text:
# Considerations on the quantum double-exchange Hamiltonian

## I Introduction

Introduced by Zener  in the early 1950s the notion of double exchange together with mixed-valency manganites $`R_{1x}A_x`$MnO3 (where $`R`$ = La, Pr, Nd and $`A`$ = Sr, Ca, Ba, Pb) attracted renewed attention when a colossal magnetoresistive effect was discovered in these compounds some years ago . The magnetic and electronic properties of manganese oxides, to some extend, are believed to arise from the large Coulomb and Hund’s rule interaction of the manganese $`d`$ shell electrons. Due to the almost octahedral coordination within the perovskite structure the $`d`$ levels split into two subbands labeled according to their octahedral symmetry, $`e_g`$ and $`t_{2g}`$. In the case of zero doping ($`x=0`$) there are four electrons per Mn site which fill up the three $`t_{2g}`$ levels and one $`e_g`$ level, and by Hund’s rule, form a $`S=2`$ spin state. Doping will remove the electron from the $`e_g`$ level, and by hopping via bridging oxygen sites these holes acquire mobility. However, this hopping acts in a background of local spins $`S=3/2`$ formed by the $`t_{2g}`$ electrons and its amplitude depends on the overlap of the spin states at neighbouring sites (or, in a classical language, on their relative angle), it is largest if the total bond spin is maximal and vice versa .

Another ingredient, that is assumed to significantly influence the physical properties of manganites, is electron-lattice interaction. Namely the two $`e_g`$ orbitals, which are degenerate in a perfect cubic environment, will couple to lattice vibrations of the same symmetry, giving rise to a Jahn-Teller effect and polaronic behaviour in some regions of the phase diagram. It is this close interplay of three different subsystems (electrons in degenerate orbitals, background of localized spins, and lattice vibrations) that makes the physics of manganites both, rich and complicated.

In the present work we concentrate on the double exchange (DE) part of the interactions and consider different possibilities for an approximate treatment of the exact DE Hamiltonian on a lattice in terms of effective electronic one- or two-band models. These can be used in a more elaborate modelling of the real materials (see our forthcoming work ). It turns out that quantum double exchange on a lattice is most suitably derived and described with the help of Schwinger bosons. We therefore include a detailed and pedagogic derivation of the quantum DE Hamiltonian using Schwinger bosons. Although this approach has been used before , we feel a comprehensive presentation of the subject is still missing. In two appendices we reexamine the derivation for two sites, and consider the semiclassical limit ($`S\mathrm{}`$). In addition, by means of numerical experiments, we illustrate how this limit evolves from the quantum case.

## II Schwinger boson representation of double exchange

To derive the quantum DE Hamiltonian on a lattice, as a starting point we take the Kondo lattice model including on-site Coulomb repulsion,

$`H`$ $`=`$ $`t{\displaystyle \underset{ij\sigma }{}}\left[c_{i\sigma }^{}c_{j\sigma }+\text{H.c.}\right]J_H{\displaystyle \underset{i\sigma \sigma ^{}}{}}(𝐒_i𝝈_{\sigma \sigma ^{}})c_{i\sigma }^{}c_{i\sigma ^{}}`$ (2)
$`+U{\displaystyle \underset{i}{}}n_in_i,`$
where summation is over nearest neighbour bonds $`ij`$ or sites $`i`$, respectively. For clarity and since it can be included easily in the final result, here we have neglected the orbital degeneracy of the $`e_g`$ electrons. That is, $`c_{i\sigma }^{()}`$ denote electrons in a single band, which interact with some localized spin $`𝐒_i`$ via the Hund’s coupling $`J_H`$. In the real materials this localized spin corresponds to the remaining $`t_{2g}`$ electrons, that tend to form a high spin state with an electron in the $`e_g`$ shell.

In the manganites the situation is such, that $`UJ_H>t`$ (cf. Refs.  ). Hence, we first take the limit $`U\mathrm{}`$, resulting in

$$H=t\underset{ij\sigma }{}\left[\stackrel{~}{c}_{i\sigma }^{}\stackrel{~}{c}_{j\sigma }+\text{H.c.}\right]J_H\underset{i\sigma \sigma ^{}}{}(𝐒_i𝝈_{\sigma \sigma ^{}})\stackrel{~}{c}_{i\sigma }^{}\stackrel{~}{c}_{i\sigma ^{}}$$
(3)
with restricted fermions $`\stackrel{~}{c}_{i\sigma }=c_{i\sigma }(1n_{i\sigma })`$. Next, following Kubo and Ohata , the exchange term in Eq. (3) is solved while the hopping term is considered as a small perturbation. For positive $`J_H`$ the ground-state of the exchange term is a free spin $`S`$, if there is no electron at site $`i`$, or a coupled spin $`\overline{S}=S+1/2`$ otherwise (note, that we use $`S`$ for the length of the localized spin $`𝐒_i`$ formed by $`t_{2g}`$ electrons). To describe the effective hopping we therefore need a projection operator, which restores these conditions,

$$(P_i^+)_{\sigma \sigma ^{}}=\frac{(𝐒_i𝝈_{\sigma \sigma ^{}})+(S+1)\delta _{\sigma \sigma ^{}}}{2S+1}.$$
(4)
Then the DE-Hamiltonian (Eq. (2.3) in Ref. ) is given in terms of spin and restricted fermion operators by

$$H_{\mathrm{el}}^{\mathrm{DE}}=t\underset{ij\sigma \sigma ^{}}{}\left[\stackrel{~}{c}_{i\sigma }^{}(P_i^+P_j^+)_{\sigma \sigma ^{}}\stackrel{~}{c}_{j\sigma ^{}}+\text{H.c.}\right].$$
(5)
However, this expression turns out to be unwieldy for analytic as well as numeric calculations. Although in principle the electronic spin is absorbed into the total spin at each site, the spin index $`\sigma `$ is still present in Eq. (5). Here the advantages of Schwinger bosons come into play, namely the possibility to describe spins of arbitrary amplitude with the same set of boson operators $`a_i`$ and $`b_i`$,

$`S_i^+`$ $`=`$ $`a_i^{}b_i,S_i^{}=b_i^{}a_i,`$ (6)
$`S_i^z`$ $`=`$ $`(a_i^{}a_ib_i^{}b_i)/2,`$ (7)
$`|S_i|`$ $`=`$ $`(a_i^{}a_i+b_i^{}b_i)/2.`$ (8)
Using these operators, we can rewrite the projection operators $`P_i^+`$,

$`(P_i^+)_{\sigma \sigma ^{}}`$ $`=`$ $`{\displaystyle \frac{1}{2S+1}}\left[\begin{array}{cc}(S+1)+S_i^z& S_i^{}\\ S_i^+& (S+1)S_i^z\end{array}\right]`$ (11)
$`=`$ $`{\displaystyle \frac{1}{2S+1}}\left[\begin{array}{cc}a_ia_i^{}& a_ib_i^{}\\ b_ia_i^{}& b_ib_i^{}\end{array}\right],`$ (14)
where we can keep $`S`$ in the denominator, because it is conserved. The last matrix can be decomposed easily,

$$(P_i^+)_{\sigma \sigma ^{}}=\frac{1}{2S+1}\left[\begin{array}{c}a_i\\ b_i\end{array}\right]\left[\begin{array}{cc}a_i^{}& b_i^{}\end{array}\right],$$
(15)
which leads to

$$H_{\mathrm{el}}^{\mathrm{DE}}=\frac{t}{2S+1}\underset{ij}{}\left[(R_i^+)^{}(a_i^{}a_j+b_i^{}b_j)R_j^++\text{H.c.}\right]$$
(16)
with the projectors

$$R_i^+=\frac{\stackrel{~}{c}_ia_i^{}+\stackrel{~}{c}_ib_i^{}}{\sqrt{2S+1}}.$$
(17)
If we restrict the Hilbert space of the problem by fixing the spin length at each site to the value $`S+n_i/2`$ (where $`n_i=\stackrel{~}{n}_i+\stackrel{~}{n}_i`$), we can simply replace the projectors $`R_i^+`$ by spinless fermion operators $`c_i`$. This can be seen, by taking a closer look on the operators $`R_i^+`$. Let us assume, site $`i`$ is occupied by an electron, and therefore the total spin at the site is $`S+1/2`$. The corresponding state $`|n_e;S,m`$ can be written as

$`|1;S+\frac{1}{2},m`$ $`=`$ $`\sqrt{{\displaystyle \frac{S+\frac{1}{2}+m}{2S+1}}}||S,m\frac{1}{2}`$ (18)
$`+`$ $`\sqrt{{\displaystyle \frac{S+\frac{1}{2}m}{2S+1}}}||S,m+\frac{1}{2}.`$ (19)
Using the representation of an $`S^z`$ eigenstate in terms of Schwinger bosons,

$$|S,m=\frac{(a^{})^{S+m}(b^{})^{Sm}}{\sqrt{(S+m)!(Sm)!}}|0,$$
(20)
we find, that applying the operator $`R^+`$ the state $`|1;S+1/2,m`$ in Eq. (18) is transformed into the corresponding Schwinger boson representation of the coupled spin $`S+1/2`$, while the electron is annihilated

$$R^+|1;S+\frac{1}{2},m=|0;S+\frac{1}{2},m.$$
(21)
Backwards, the operator $`(R^+)^{}`$ creates the decomposition of the coupled spin into electronic and localized spin, i.e., it produces appropriate Clebsch-Gordan coefficients.

Now it is straightforward to omit the electronic spin index, using only spinless fermions, $`c_i^{()}`$, and Schwinger bosons to describe double-exchange. The corresponding Hamiltonian is given by

$$H_{\mathrm{el}}^{\mathrm{DE}}=\frac{t}{2S+1}\underset{ij}{}\left[(a_i^{}a_j+b_i^{}b_j)c_i^{}c_j+\text{H.c.}\right],$$
(22)
where, for every site $`i`$, the Hilbert space is constrained to

$$a_i^{}a_i+b_i^{}b_i=2S+c_i^{}c_i.$$
(23)
In the case of low doping usually it is more appropriate and natural to consider holes instead of electrons. Here “hole” denotes a fermion and a spin $`S`$ moving together in a background of spins $`\overline{S}=S+1/2`$. For the transformation of the electronic Hamiltonian, Eq. (16), the operator $`R^+`$ needs to be replaced by a suitable counterpart involving restricted hole operators $`\stackrel{~}{h}_{i\sigma }`$. In analogy to Eqs. (18) to (21), the state

$`|n_h=1;S,m`$ $`=`$ $`\sqrt{{\displaystyle \frac{\overline{S}m+\frac{1}{2}}{2\overline{S}+1}}}||\overline{S},m\frac{1}{2}`$ (24)
$``$ $`\sqrt{{\displaystyle \frac{\overline{S}+m+\frac{1}{2}}{2\overline{S}+1}}}||\overline{S},m+\frac{1}{2}`$ (25)
is transformed into its corresponding Schwinger boson representation by the operator

$$R_i^{}=\frac{\stackrel{~}{h}_ib_i\stackrel{~}{h}_ia_i}{\sqrt{2\overline{S}+1}}.$$
(26)
Hence, in hole representation, the DE Hamiltonian reads

$$H_{\mathrm{hole}}^{\mathrm{DE}}=\frac{t}{2\overline{S}}\underset{ij}{}\left[(R_i^{})^{}(a_ia_j^{}+b_ib_j^{})R_j^{}+\text{H.c.}\right].$$
(27)
The Hamiltonian for spinless fermions, Eq. (22), changes only little, becoming

$$H_{\mathrm{hole}}^{\mathrm{DE}}=\frac{t}{2\overline{S}}\underset{ij}{}\left[(a_ia_j^{}+b_ib_j^{})h_i^{}h_j+\text{H.c.}\right]$$
(28)
with the constraint

$$a_i^{}a_i+b_i^{}b_i=2\overline{S}h_i^{}h_i.$$
(29)
## III Effective transport Hamiltonian

To obtain an effective Hamiltonian for the spin-dependent hole-hopping, the spin part of the DE interaction is considered within mean field approximation. However, there are two representations of the DE Hamiltonian to start from: Eqs. (27) and (28). The resulting effective Hamiltonians describe carriers with or without spin, respectively. Given Eq. (28), a mean hopping of the spinless carriers, $`\stackrel{~}{t}^{(\mathrm{b})}`$, is obtained by considering each bond $`ij`$ separately, and taking the expectation value of the DE term for all values of the total bond spin $`S_T`$ and $`S_T^z`$ in an effective ordering field $`\lambda =\beta g\mu _BH_{\mathrm{eff}}^z`$. This reproduces the result of Kubo and Ohata . Having spin $`\overline{S}`$ at site $`i`$ and spin $`S`$ at site $`j`$, the coupled state $`S_T`$ with maximal $`S_T^z`$ is given by

$`|S_T,S_T_{(\overline{S}S)}=`$ (30)
$`=`$ $`C(b_i^{}a_j^{}b_j^{}a_i^{})^{2S+\frac{1}{2}S_T}(a_i^{})^{S_T+\frac{1}{2}}(a_j^{})^{S_T\frac{1}{2}}|0,`$ (31)
where $`C`$ is a normalization factor,

$$C=\sqrt{\frac{(2S_T+1)!}{(S_T+\frac{1}{2})!(S_T\frac{1}{2})!(2S+\frac{1}{2}S_T)!(2S+\frac{3}{2}+S_T)!}}.$$
(32)
Applying $`(a_j^{}a_i+b_j^{}b_i)`$ we arrive at

$`(a_j^{}a_i+b_j^{}b_i)|S_T,S_T_{(\overline{S}S)}=`$ (33)
$`=`$ $`C(S_T+\frac{1}{2})(b_i^{}a_j^{}b_j^{}a_i^{})^{2S+\frac{1}{2}S_T}(a_i^{})^{S_T\frac{1}{2}}(a_j^{})^{S_T+\frac{1}{2}}|0`$ (34)
$`=`$ $`(S_T+\frac{1}{2})|S_T,S_T_{(S\overline{S})}.`$ (35)
Hence, we rederived the effective matrix element for a single bond (cf. Appendix A or Ref. ),

$$t^{(\mathrm{b})}=\frac{S_T+1/2}{2\overline{S}}t,$$
(36)
which is averaged over all values and directions of $`S_T`$,

$$\stackrel{~}{t}^{(\mathrm{b})}=\gamma _{\overline{S}}[\overline{S}\lambda ]t,$$
(37)
where

$`\gamma _{\overline{S}}[\overline{S}\lambda ]={\displaystyle \frac{_{S_T=1/2}^{2\overline{S}1/2}_{M=S_T}^{S_T}\frac{S_T+1/2}{2\overline{S}}\text{ e}^{M\lambda }}{_{S_T=1/2}^{2\overline{S}1/2}_{M=S_T}^{S_T}\text{ e}^{M\lambda }}}`$ (38)
$`=`$ $`\frac{1}{2}+\frac{\overline{S}}{2\overline{S}+1}\mathrm{coth}\left(\frac{2\overline{S}+1}{2}\lambda \right)\left[\mathrm{coth}(\overline{S}\lambda )\frac{1}{2\overline{S}}\mathrm{coth}\left(\frac{\lambda }{2}\right)\right].`$ (39)
The effective Hamiltonian describing spinless fermionic holes in an averaged background of ordered spins reads

$$H_{\mathrm{hole}}^{\mathrm{eff},\mathrm{I}}=\stackrel{~}{t}^{(\mathrm{b})}\underset{ij}{}\left[h_i^{}h_j+\text{H.c.}\right].$$
(40)
In Appendix B we compare the classical limit of both, this Hamiltonian and the exact expression. At least the bandwidth turns out to be represented very well.

Another possible way to obtain an effective Hamiltonian is based on the picture of itinerant carriers of spin $`\frac{1}{2}`$ moving in the background of localized spins, the correlations of which change on a large time scale compared with the hopping frequency. Then an effective hopping Hamiltonian is obtained averaging $`H_{\mathrm{hole}}^{\mathrm{DE}}`$, Eq. (27), over free spins $`\overline{S}`$ in a homogeneous field $`\lambda `$. According to Eqs. (6)-(7) and the fact that $`S^\pm =0`$ and $`S^z=\overline{S}B_{\overline{S}}[\overline{S}\lambda ]`$ (where $`B_{\overline{S}}[z]`$ denotes the Brillouin function), only two terms contribute. The resulting Hamiltonian involves two effective hopping matrix elements, one for each spin channel,

$$H_{\mathrm{hole}}^{\mathrm{eff},\mathrm{II}}=\underset{ij}{}\left[\stackrel{~}{t}_{}\stackrel{~}{h}_i^{}\stackrel{~}{h}_j+\stackrel{~}{t}_{}\stackrel{~}{h}_i^{}\stackrel{~}{h}_j+\text{H.c.}\right]$$
(41)
with

$`\stackrel{~}{t}_{}`$ $`=`$ $`{\displaystyle \frac{t}{(2\overline{S})(2\overline{S}+1)}}\left[\overline{S}(1B_{\overline{S}}[\overline{S}\lambda ])\right]^2,`$ (42)
$`\stackrel{~}{t}_{}`$ $`=`$ $`{\displaystyle \frac{t}{(2\overline{S})(2\overline{S}+1)}}\left[\overline{S}(1+B_{\overline{S}}[\overline{S}\lambda ])\right]^2.`$ (43)
In a fully polarized background ($`\lambda \mathrm{}`$) only holes with anti-parallel spin can hop ($`\stackrel{~}{t}_{}\frac{2\overline{S}}{2\overline{S}+1}t`$), while holes with parallel spin are blocked ($`\stackrel{~}{t}_{}0`$). In general, the situation is complicated by the fact that Eq. (41) involves restricted fermion operators (Hubbard operators) forbidding an exact solution of the model. However, to a good approximation, in the polarized phase the up-band can be neglected, while the down-band is taken into account using unrestricted operators $`h_i^{()}`$. On the other hand, in a disordered phase ($`\lambda 0`$) both bands are equivalent ($`\stackrel{~}{t}_{}=\stackrel{~}{t}_{}\frac{\overline{S}/2}{2\overline{S}+1}t`$), making an approximate treatment less evident.

## IV Conclusions

In the present work we review the subject of double exchange and derive an effective Hamiltonian for the spin-dependent hopping of holes in an averaged background of local spins. This is used in a subsequent publication within a two-phase scenario for the description of colossal magnetoresistant manganites, alternatively to the effective Hamiltonian for spinless carriers derived in Ref. . Besides, we illustrate that on the level of quantum spins all features of double exchange - namely its derivation, limiting cases, as well as approximations - are probably most clearly represented in terms of Schwinger bosons.

Acknowledgements: This work was supported by the Deutsche Forschungsgemeinschaft and the Czech Academy of Sciences under Grant No. 436 TSE 113/33.

## A Reexamination of the 2-site problem

As noted above, the DE matrix element, Eq. (36), was first derived by Anderson and Hasegawa  by considering a system of two spins $`S`$ at neighbouring sites and a mobile electron whose spin $`\sigma `$ is coupled to the local spin at the same site. Recently the problem was reexamined by Müller-Hartmann and Dagotto , who found a “nontrivial” phase factor, which the authors interpreted in terms of a Berry phase within the limit $`S\mathrm{}`$. We argue that the phase factor of their quantum DE Hamiltonian merely compensates for the phase introduced by the permutation of neighbouring spin states, whereas independently the Berry phase evolves in the classical ($`S\mathrm{}`$) limit of the correct DE Hamiltonian, Eq. (22). Moreover, on a lattice double exchange can not be described in terms of spin and permutation operators, which do not take into account fermionic commutation relations.

To be specific we briefly summarize the two site problem. Assuming the electron initially to be at site 1 we start with coupling the local spin $`𝐒_1`$ and the electron spin $`𝝈`$, to give an on-site spin $`\overline{𝐒}_1=𝐒_1+𝝈`$. We note, that the construction of the corresponding wave function is not unique, as is known from textbooks on spin algebra . In more detail

$`|\overline{S}_1\overline{m}_1_{(\sigma S_1)}`$ $`=`$ $`{\displaystyle \underset{\mu m_1}{}}\left(\begin{array}{ccc}\sigma & S_1& \overline{S}_1\\ \mu & m_1& \overline{m}_1\end{array}\right)|\sigma \mu |S_1m_1,`$ (A3)
$`|\overline{S}_1\overline{m}_1_{(S_1\sigma )}`$ $`=`$ $`(1)^{\sigma +S_1\overline{S}_1}|\overline{S}_1\overline{m}_1_{(\sigma S_1)},`$ (A4)
where we give the order of the constituting spin states as a subscript and denote Clebsch-Gordan coefficients by $`()`$. Adding this spin to the unchanged core spin at site 2, we arrive at an initial state of total spin $`S_T`$

$`|\phi _{\mathrm{in}}`$ $`=`$ $`|S_Tm_T_{((\sigma S_1)S_2)}`$ (A5)
$`=`$ $`{\displaystyle \underset{\overline{m}_1m_2}{}}\left(\begin{array}{ccc}\overline{S}_1& S_2& S_T\\ \overline{m}_1& m_2& m_T\end{array}\right)|\overline{S}_1\overline{m}_1_{(\sigma S_1)}|S_2m_2.`$ (A8)
Now, if the electron moves to site 2, we can couple $`𝝈`$ and $`𝐒_2`$ to give $`\overline{𝐒}_2`$, which together with $`𝐒_1`$ yields another state of total spin $`S_T`$. However, as can be seen above, there are different ways to connect the three spins. One possible final state is

$`|\phi _{\mathrm{fi}}^A`$ $`=`$ $`|S_Tm_T_{(S_1(\sigma S_2))}`$ (A9)
$`=`$ $`{\displaystyle \underset{m_1\overline{m}_2}{}}\left(\begin{array}{ccc}S_1& \overline{S}_2& S_T\\ m_1& \overline{m}_2& m_T\end{array}\right)|S_1m_1|\overline{S}_2\overline{m}_2_{(\sigma S_2)}.`$ (A12)
The corresponding effective hopping matrix element $`t_A^{(\mathrm{b})}`$ is proportional to the overlap $`\phi _{\mathrm{fi}}^A|\phi _{\mathrm{in}}`$, which can be evaluated by rewriting the multiple sum over the product of four Clebsch-Gordan coefficients with a 6-j-symbol,

$`\phi _{\mathrm{fi}}^A|\phi _{\mathrm{in}}=`$ (A16)
$`(1)^{\overline{S}_1+S_2+S_T}\sqrt{(2\overline{S}_1+1)(2\overline{S}_2+1)}\left\{\begin{array}{ccc}\sigma & S_1& \overline{S}_1\\ S_T& S_2& \overline{S}_2\end{array}\right\}.`$
We obtain the value of the matrix relevant for strong Hund’s coupling by setting $`|\overline{𝐒}_i|=S_i+1/2\overline{S}`$ and $`|𝐒_i|=S=\overline{S}1/2`$:

$$t_A^{(\mathrm{b})}=\frac{S_T+1/2}{2\overline{S}}t,$$
(A17)
i.e., there is no $`S_T`$-dependent phase factor. However, Müller-Hartmann and Dagotto derive an effective double-exchange Hamiltonian, where the hopping is expressed as the permutation of the spin $`\overline{S}`$ “particle” at site 1 and the spin $`(\overline{S}1/2)`$ “hole” at site 2. Therefore the final state they have to consider for the matrix element is

$`|\phi _{\mathrm{fi}}^B`$ $`=`$ $`|S_Tm_T_{((\sigma S_2)S_1)}`$ (A18)
$`=`$ $`{\displaystyle \underset{\overline{m}_2m_1}{}}\left(\begin{array}{ccc}\overline{S}_2& S_1& S_T\\ \overline{m}_2& m_1& m_T\end{array}\right)|\overline{S}_2\overline{m}_2_{(\sigma S_2)}|S_1m_1.`$ (A21)
Obviously the permutation of the indices in the Clebsch-Gordan coefficients yields a different phase factor in the overlap

$`\phi _{\mathrm{fi}}^B|\phi _{\mathrm{in}}=`$ (A25)
$`(1)^{S_1+S_2+\overline{S}_1+\overline{S}_2}\sqrt{(2\overline{S}_1+1)(2\overline{S}_2+1)}\left\{\begin{array}{ccc}\sigma & S_1& \overline{S}_1\\ S_T& S_2& \overline{S}_2\end{array}\right\},`$
and consequently the authors obtain a $`S_T`$-dependent phase factor for $`t_B^{(\mathrm{b})}`$,

$$t_B^{(\mathrm{b})}=(1)^{2\overline{S}S_T1/2}\frac{S_T+1/2}{2\overline{S}}t.$$
(A26)
To recover the effective Hamiltonian Eq.(6) of Ref. ,

$$H^{\mathrm{eff}}=tP_{12}Q_{\overline{S}}(y),$$
(A27)
we simply have to express $`S_T`$ in $`t_B^{(\mathrm{b})}`$ with the help of $`y=𝐒_1𝐒_2/(\overline{S}(\overline{S}1/2))`$. In a direct way we can use interpolation polynomials in the form of Lagrange, i.e.

$`Q_{\overline{S}}(y)`$ $`=`$ $`{\displaystyle \underset{S_T=1/2}{\overset{2\overline{S}1/2}{}}}t_B^{(\mathrm{b})}{\displaystyle \underset{\genfrac{}{}{0.0pt}{}{j=1/2}{jS_T}}{\overset{2\overline{S}1/2}{}}}{\displaystyle \frac{yy_j}{y_{S_T}y_j}}`$ (A28)
$`y_j`$ $`=`$ $`{\displaystyle \frac{j(j+1)\overline{S}(\overline{S}+1)(\overline{S}1/2)(\overline{S}+1/2)}{2\overline{S}(\overline{S}1/2)}},`$ (A29)
or the recursive formula given in Eq.(7) of Ref. . Of course, taking the matrix element $`t_A^{(\mathrm{b})}`$ for the construction of $`Q_{\overline{S}}(y)`$ is wrong, as this does not account for the phase factor due to the permutation $`P_{12}`$. Concerning the limit of classical spins, we can see no connection between the above $`S_T`$-dependent phase factor and the Berry phase.

The peculiarities concerning the phase factor indicate that the above derivation is not suitable for the generalization to a lattice. Moreover, the spin-$`\overline{S}`$ “particles” still obey fermion commutation relations, which can not be expressed by conventional permutation or spin operators. That means, the expression given in Eq.(6) of Ref. ,

$$H^{\mathrm{eff}}=t\underset{ij}{}P_{ij}Q_{\overline{S}}(y),$$
(A30)
is not the correct quantum DE-Hamiltonian on a lattice.

## B Classical limit of the DE model

The limit $`S\mathrm{}`$ of $`H_{\mathrm{el}}^{\mathrm{DE}}`$, Eq. (22), is easily derived by taking its expectation value with spin coherent states ,

$$|\mathrm{\Omega }(S,\theta ,\varphi )=\frac{(ua^{}+vb^{})^{2S}}{\sqrt{(2S)!}}|0,$$
(B1)
where $`u=\mathrm{cos}(\theta /2)\text{ e}^{i\varphi /2}`$ and $`v=\mathrm{sin}(\theta /2)\text{ e}^{i\varphi /2}`$. Using the properties of coherent states,

$`a|\mathrm{\Omega }(S,\theta ,\varphi )`$ $`=`$ $`\sqrt{2S}u|\mathrm{\Omega }(S\frac{1}{2},\theta ,\varphi )`$ (B2)
$`b|\mathrm{\Omega }(S,\theta ,\varphi )`$ $`=`$ $`\sqrt{2S}v|\mathrm{\Omega }(S\frac{1}{2},\theta ,\varphi ),`$ (B3)
for a given spin configuration $`\{\theta _k,\varphi _k\}`$ and two electronic states $`|\psi _1`$ and $`|\psi _2`$,

$$|\psi _j=\underset{k}{}|n_{j,k}|\mathrm{\Omega }(S+\frac{n_{j,k}}{2},\theta _k,\varphi _k)$$
(B4)
(where $`|n_{j,k}=(c_k^{})^{n_{j,k}}|0`$ with numbers $`n_{j,k}\{0,1\}`$), we find the average

$$\psi _1|H_{\mathrm{el}}^{\mathrm{DE}}|\psi _2=\underset{k}{}n_{1,k}|\left(\underset{ij}{}\left[t_{ij}c_i^{}c_j+\text{H.c.}\right]\right)\underset{k}{}|n_{2,k}$$
(B5)
with the matrix element

$`t_{ij}`$ $`=`$ $`\mathrm{cos}\left({\displaystyle \frac{\theta _i}{2}}\right)\mathrm{cos}\left({\displaystyle \frac{\theta _j}{2}}\right)\text{ e}^{i(\varphi _i\varphi _j)/2}`$ (B6)
$`+`$ $`\mathrm{sin}\left({\displaystyle \frac{\theta _i}{2}}\right)\mathrm{sin}\left({\displaystyle \frac{\theta _j}{2}}\right)\text{ e}^{i(\varphi _i\varphi _j)/2}.`$ (B7)
Hence, the classical Hamiltonian should read

$$H_{\mathrm{class}}^{\mathrm{DE}}=\underset{ij}{}\left[t_{ij}c_i^{}c_j+\text{H.c.}\right],$$
(B8)
which is equivalent to the results obtained in Refs. . In our case, however, the classical limit followed from the quantum Hamiltonian, Eq. (22), in an obvious and more straightforward way. Note also, that averaging the operators $`(R_i^+)^{}`$ and $`R_i^+`$ over coherent states yields the unitary transformation onto rotated electrons $`d_j^{()}`$ (compare Refs. ),

$$d_j=\mathrm{cos}\left(\frac{\theta _j}{2}\right)\text{ e}^{i\varphi _j/2}\stackrel{~}{c}_j+\mathrm{sin}\left(\frac{\theta _j}{2}\right)\text{ e}^{i\varphi _j/2}\stackrel{~}{c}_j,$$
(B9)
i.e., naturally Eq. (16) has the same classical limit.

To check, whether the description of double-exchange in terms of classical spins is appropriate, we compared the (canonical) density of states (DOS) for a fixed number of carriers on a small cluster, which interact with quantum (Eq. (22)) or classical (Eq. (B8)) spins, respectively. Using Chebychev expansion and maximum entropy methods , for the quantum case the spectra can be obtained numerically for rather large $`S`$. The classical DOS is found by averaging the eigenvalues of $`H_{\mathrm{class}}^{\mathrm{DE}}`$, Eq. (B8), over a large number of spin configurations.

In Fig. 1 we consider two electrons on a ring of four sites. Comparing the running average (bold dot-dashed) over the discrete spectrum (thin dashed) and the classical limit (bold solid), we find good convergence already for a moderate spin length $`S=10`$. But even for the case $`S=3/2`$, which is realized in the manganites, the classical description appears to be acceptable. If we consider four electrons on a ring of eight sites, the spectrum is much more dense, making similarities to be recognized easily, see Fig. 2. Of course, from the density of states we can learn nothing about correlations or other more involved features. Note, that in both figures we subtracted the peak at $`E=0`$ consuming a large fraction of spectral weight.

Another interesting check concerns the effective hopping matrix element $`\stackrel{~}{t}^{(\mathrm{b})}`$, Eq. (37). With classical spins and the Chebychev expansion methods mentioned above, one can calculate the grand-canonical DOS of the tight-binding model, Eq. (B8), for rather large clusters (here $`64^3`$ sites on a simple cubic lattice) without much effort. By considering a thermalized ensemble of classical spins in a homogeneous field we compare the resulting bandwidth with the limit $`S\mathrm{}`$ of $`\stackrel{~}{t}^{(\mathrm{b})}`$, where the limit of $`\gamma _{\overline{S}}[\overline{S}\lambda ]`$ is obtained omitting $`\overline{S}`$ in the argument and setting $`\overline{S}\mathrm{}`$ in the index,

$$\gamma _{\overline{S}\mathrm{}}[\lambda ]=\frac{1}{2}\left(1+\mathrm{coth}(\lambda )\left[\mathrm{coth}(\lambda )\frac{1}{\lambda }\right]\right).$$
(B10)
The agreement is rather satisfactory, as can be seen in Fig. 3. Naturally, the precise shape of the DOS (see inset) will not be reproduced by the effective electronic model, it remains simple cubic tight-binding for all fields and temperatures.