# The Zel’dovich approximation

###### Abstract

This year marks the anniversary of the birth of Yakov Zel’dovich. Amongst his many legacies is the Zel’dovich approximation for the growth of large-scale structure, which remains one of the most successful and insightful analytic models of structure formation. We use the Zel’dovich approximation to compute the two-point function of the matter and biased tracers, and compare to the results of N-body simulations and other Lagrangian perturbation theories. We show that Lagrangian perturbation theories converge well and that the Zel’dovich approximation provides a good fit to the N-body results except for the quadrupole moment of the halo correlation function. We extend the calculation of halo bias to order and also consider non-local biasing schemes, none of which remove the discrepancy. We argue that a part of the discrepancy owes to an incorrect prediction of inter-halo velocity correlations. We use the Zel’dovich approximation to compute the ingredients of the Gaussian streaming model and show that this hybrid method provides a good fit to clustering of halos in redshift space down to scales of tens of Mpc.

###### keywords:

gravitation; galaxies: haloes; galaxies: statistics; cosmological parameters; large-scale structure of Universe^{†}

^{†}pagerange: The Zel’dovich approximation–A

## 1 Introduction

This year marks the anniversary of the birth of Yakov Zel’dovich, who was a pioneer in the study of large-scale structure and introduced the approximate dynamics that bears his name (Zel’dovich, 1970). The Zel’dovich approximation provides an intuitive way to understand the emergence of the beaded filamentary structure which has become known as the cosmic web and a fully realized (though approximate) model of non-linear structure formation (Peebles, 1980; Coles & Lucchin, 1995; Peacock, 1999). The Zel’dovich approximation predicts the rich structure of voids, clusters, sheets and filaments observed in the Universe (Doroshkevich et al., 1980; Pauls & Melott, 1995), and indeed it provides a reasonably good match to N-body simulations on large scales (Coles, Melott & Shandarin, 1993; Tassev & Zaldarriaga, 2012a, c). For a discussion of why the Zel’dovich approximation works so well, see Buchert (1989); Pauls & Melott (1995); Yoshisato et al. (2006); Tassev (2014a). For reviews of the Zel’dovich approximation, see the textbooks referenced above and Shandarin & Zeldovich (1989); Sahni & Coles (1995); Coles & Sahni (1996); Gurbatov, Saichev & Shandarin (2012); Hidding, Shandarin & van de Weygaert (2014).

The last few years have seen a resurgence of interest in the Zel’dovich approximation. It has been applied to understanding the effects of non-linear structure formation on the baryon acoustic oscillation feature in the correlation function (Padmanabhan & White, 2009; McCullagh & Szalay, 2012; Tassev & Zaldarriaga, 2012a) and to understanding how “reconstruction” (Eisenstein, et al., 2007) removes those non-linearities (Padmanabhan, White & Cohn, 2009; Noh, White & Padmanabhan, 2009; Tassev & Zaldarriaga, 2012b). It has been used as the basis for an effective field theory of large-scale structure (Porto, Senatore & Zaldarriaga, 2014). It has been compared to “standard” perturbation theory (Tassev, 2014a), extended to higher orders in Lagrangian perturbation theory (Matsubara, 2008a, b; Okamura, Taruya, & Matsubara, 2011; Carlson, Reid & White, 2013) and to higher order statistics (Tassev, 2014b). Despite the more than 40 years since it was introduced, the Zel’dovich approximation still provides one of our most accurate models for the distribution of cosmological objects.

In this paper we investigate to what extent the Zel’dovich approximation can be used as a quantitatively accurate model of the low-order clustering of objects in cosmology. The outline is as follows. After some background and review to establish notation in Section 2 we present a derivation of the 2-point function within the Zel’dovich approximation (see also Carlson, Reid & White, 2013; Tassev, 2014a, b) both for matter (Section 3) and for biased tracers (Section 4). In these sections we show that the principle ingredient to the calculation, the Lagrangian correlator, can be inverted analytically and thus the correlation function expressed as a simple quadrature. All of the ingredients to the Zel’dovich approximation involve only one dimensional integrals of the linear theory power spectrum, and these can be efficiently precomputed and tabulated, making numerical evaluation fast and efficient. We compare the Zel’dovich calculation to some other Lagrangian perturbation theory schemes, and to the results of N-body simulations, in Section 5. We extend previous calculations of the effects of bias to higher order and include non-local, Lagrangian bias in Section 6. Finally we introduce the Zel’dovich Streaming Model (ZSM) as a hybrid method for accurately computing the redshift-space correlation function of biased tracers in Section 7 and end with a discussion of our results in Section 8.

For plots and numerical comparisons we assume a CDM cosmology with , , , , and . Our simulation data are derived from a suite of 20 N-body simulations run with the TreePM code described in White (2002). Each simulation employed equal mass () particles in a periodic cube of side length Gpc as described in (Reid & White, 2011; White et al., 2011). Halos are found using the friends-of-friends method, with a linking length of times the mean inter-particle spacing.

## 2 Background and review

In this section we provide a brief review of cosmological perturbation theory,
focusing on the Lagrangian formulation^{1}^{1}1See Bernardeau et al. (2002) for a
comprehensive (though somewhat dated) review of Eulerian perturbation theory.
(Buchert, 1989; Moutarde et al., 1991; Hivon et al., 1995; Taylor & Hamilton, 1996).
This material should be sufficient to remind the reader of some essential
terminology, and to establish our notational conventions.
Our notation and formalism follows closely that in
Matsubara (2008a, b); Carlson, Reid & White (2013); Wang, Reid & White (2013)
to which we refer the reader for further details.

In terms of the fractional density perturbation, , we can write the 2-point correlation function,

(1) |

and its Fourier transform, the power spectrum , as

(2) |

Here denotes the 3-dimensional Dirac delta function, and we use the Fourier transform convention

(3) |

Angle brackets around a cosmological field, e.g. , signify an ensemble average of that quantity over all possible realizations of our universe; in most cases of interest, ergodicity allows us to replace these ensemble averages with spatial averages over a sufficiently large cosmic volume.

In the Lagrangian approach to cosmological fluid dynamics, one traces the trajectory of an individual fluid element through space and time. For a fluid element located at position at some initial time , its position at subsequent times can be written in terms of the Lagrangian displacement field ,

(4) |

where . Every element of the fluid is uniquely labeled by its Lagrangian coordinate and the displacement field fully specifies the motion of the cosmological fluid. Lagrangian Perturbation Theory (LPT) finds a perturbative solution for the displacement field,

(5) |

The first order solution is the Zel’dovich approximation (Zel’dovich, 1970), which shall be the focus of this paper. Henceforth we shall denote simply as :

(6) |

This formalism makes it particularly easy to include redshift space distortions. In this work we adopt the standard “plane-parallel” or “distant-observer” approximation, in which the line-of-sight direction to each object is taken to be the fixed direction . This has been shown to be a good approximation at the level of current observational error bars (e.g., Figure 10 of Samushia, Percival, & Raccanelli 2012 or Figure 8 of Yoo & Seljak 2014). Under this assumption, the position of an object located at true comoving position , will be mis-identified due to its peculiar velocity along the line-of-sight, as

(7) |

Thus including redshift-space distortions requires only a simple additive offset of the displacement field. The peculiar velocity of a fluid element, labeled by its Lagrangian coordinate , is so in redshift space the apparent displacement of the fluid element is

(8) |

To a good approximation the time dependence of the th order term in Eq. (5) is given by . Therefore , where is the growth rate, often approximated as . Thus, within the Zel’dovich approximation, the mapping to redshift space is achieved via the matrix

(9) |

which simply multiplies the -component of the vector by .

Finally we must consider biased tracers of the density field. We begin by considering a local Lagrangian bias model, which posits that the locations of discrete tracers at some late time are determined by the overdensities in the initial matter density field, specifically:

(10) |

Here is the mean comoving number density of our tracer and the function is called the Lagrangian bias function. Matsubara (2011) provides an extensive discussion of non-local Lagrangian bias.

The correlation function within the Zel’dovich approximation then follows by elementary manipulations (Bond & Couchman, 1988; Fisher & Nusser, 1996; Matsubara, 2008a, b; Carlson, Reid & White, 2013; Wang, Reid & White, 2013; Tassev, 2014b). We begin by writing

(11) |

We now replace the delta function with its Fourier representation, and also
introduce the Fourier^{2}^{2}2An alternative route to the same final
expressions is to expand in powers of and use the
properties of Gaussian integrals. We will use the Fourier methodology since
it was also used in Matsubara (2008b); Carlson, Reid & White (2013); Wang, Reid & White (2013). transform
of , so the expression for becomes

(12) | ||||

(13) |

The 2-point correlation function for the biased tracer is then given by

(14) | |||||

where , , and . By statistical homogeneity, the expectation value above depends only on the difference in Lagrangian coordinates, . The change of variables then leads to

(15) |

where we have defined

(16) |

and . This expression, derived in Carlson, Reid & White (2013), is the exact configuration space analog of Eq. (9) in Matsubara (2008b) and is analogous to the power spectrum derived in Fisher & Nusser (1996).

We can expand the expectation value in Eq. (16) in terms of cumulants. Since and are Gaussian only the second cumulant is non-zero

(17) | |||||

where we have defined

(18) | ||||||

(19) |

Eq. (16) then evaluates to

(20) |

This expression is exact, within the Zel’dovich approximation. The quantity is simply the variance of the smoothed linear density field, while is the corresponding smoothed linear correlation function. The matrix may be decomposed as

(21) | |||||

(22) |

where is the 1-D dispersion of the displacement field, and and are the transverse and longitudinal components of the Lagrangian 2-point function, . The vector is the cross-correlation between the linear density field and the Lagrangian displacement field. In the Zel’dovich approximation these quantities are given by

(23) | |||

(24) | |||

(25) | |||

(26) |

and shown in Fig. 1. Up to factors of 2 and , these expressions are identical to the Eulerian velocity correlators in linear theory (e.g. Gorski, 1988; Fisher, 1995; Reid & White, 2011), which is not surprising since in the Zel’dovich approximation. It is also useful to define , the pairwise velocity dispersion at an angle to the line-of-sight with the line-of-sight and perpendicular components and .

## 3 Correlation function – matter

For the unbiased case (i.e. the matter field) we can write our expression for in closed form by carrying out the Gaussian integral

(27) | ||||

(28) |

Further discussion of this expression, and the approach to linear theory, can be found in Carlson, Reid & White (2013).

Evaluation of involves a numerical integral of a Gaussian function. The inversion of which is required can be done analytically by use of the Sherman-Morrison formula which states that for matrices and vectors and ,

(29) |

Writing (see Eq. 22) we have

(30) | |||||

(31) |

where , and the combination . To compute the determinant we make use of det for scalar and matrix and that det. Then

(32) |

as expected.
The integrand can thus be expressed analytically in terms of the 2-point
functions defined previously and evaluated by simple quadratures^{3}^{3}3We
use the midpoint rule in and Gauss-Legendre integration in
. in .
The integral is dominated by .
By tabulating the functions and in advance these integrals
can be performed very quickly.
In redshift space we replace

(33) | ||||

(34) |

This corresponds simply to dividing the -components of the inverse of by and multiplying by .

## 4 Perturbative expansion for biased tracers

Returning to the case of biased tracers, consider again Eq. (20). In the unbiased case the integration in Eq. (15) took the form of a Gaussian integral, which we carried out analytically. In the biased case, we can achieve the same thing if we first partially expand Eq. (20) as

(35) | |||||

We may justify this choice of expansion by noting that both and vanish in the large-scale limit , while and approach non-zero values. At this stage we note one difference between doing this expansion within Lagrangian perturbation theory (e.g. Matsubara, 2008b; Okamura, Taruya, & Matsubara, 2011; Carlson, Reid & White, 2013) and performing it within the context of the Zel’dovich expansion. In the Zel’dovich approximation all the terms are just multiples of 2-point functions and we can go to arbitrarily high order without the need to evaluate any high dimensional mode coupling integrals or numerically difficult terms. Here we have gone to cubic order in the 2-point function and indicated how the quartic terms appear. We shall show later that the expansion seems to be converging quickly.

The and integrations give , the expectation value of the th derivative of the Lagrangian bias function (Matsubara, 2008b). In order to make the expressions more readable we shall write . Then

(37) |

The integration reduces to a series of multi-variate Gaussian integrals which can be done with the formulae in Carlson, Reid & White (2013). In the end we obtain

(39) |

where

(40) |

so that

(41) |

Our final expression for the correlation function is

(42) |

The remaining integration over must be performed numerically as before.

One can treat the as fitting parameters, or attempt to compute them from a bias model. One such model is the peak-background split, which begins with the unconditional multiplicity function

(43) |

which can be fit with

(44) |

where , gives the Press-Schecter mass function (Press & Schechter, 1974), while , yields the Sheth-Tormen mass function (Sheth & Tormen, 1999). Within the assumption of the peak-background split, the conditional multiplicity function is given by the substitution,

(45) |

where is the background density and is the critical overdensity for collapse. The Lagrangian bias parameters then follow from Taylor expanding the (appropriately normalized) conditional multiplicity function as a function of , yielding or

(46) |

(47) |

and

(48) | |||||

For the halo sample shown in Fig. 2 we have , , and from the Sheth-Tormen mass function and we have used these values in the relevant figures. There is some evidence (Baldauf et al., 2012) that simplistic bias models such as the above are not quantitatively accurate when compared to N-body simulations. On large scales the level of agreement is quite insensitive to the value assumed for as long as is not much larger than because the terms involving are numerically small compared to the terms. Thus assuming the peak-background split value for is perfectly adequate and in matching the Zel’dovich theory to observations there is only one free parameter ( or ). This is very well constrained by the overall amplitude of .

## 5 Results

Figs. 2 and 3 compare the predictions of the Zel’dovich approximation to a number of other Lagrangian perturbation theory schemes and to the clustering of halos in an N-body simulation. The agreement between the predictions of the Zel’dovich approximation and N-body simulations is very good on large scales. For the real-space correlation function and the redshift-space monopole the agreement extends down to Mpc for both the matter and halos. However, the theory does much less well for the halo quadrupole (a problem shared with all of the local Lagrangian bias schemes shown in Fig. 2). One major difference between the halo calculation and the matter calculation is the Taylor series expansion of the bias terms.

While the perturbative schemes plotted in Fig. 2 stop at , we have extended the Zel’dovich calculation one additional order to see if the higher order terms could contribute to the quadrupole on intermediate and small scales. We find that the terms cubic in and contribute negligibly to the real-space correlation function and the monopole and quadrupole moments of the redshift-space correlation function for Mpc. Since the terms do not contribute significantly to any of the statistics for this halo sample it appears that the Taylor series expansion is not the source of the discrepancy. This also suggests that truncating the expansion at is a good approximation and that the expansion is converging. Henceforth we shall drop the terms.

The contribution to the correlation functions is not shared equally among the terms. Working above Mpc we find that the real-space correlation function is dominated by the , “” and terms in the square brackets of Eq. (39). The redshift-space monopole is dominated by the same three terms. Fig. 4 shows how the different terms contribute to the final peak in near Mpc and to the quadrupole. Above Mpc the terms “” and in the square brackets in Eq. (39) contribute the vast majority of the total quadrupole signal, with approximately equal contributions. By Mpc the other terms contribute about 10 per cent of the total, with the , and terms contributing the remainder in decreasing order of importance (the term provides a negligible contribution).

Fig. 5 shows the degree to which the and terms depend on scale differently than the matter terms. In all cases, above Mpc the level of scale-dependence is quite small. The actual shape of the real-space correlation function and the redshift-space monopole correlation function differ at the ten per cent level near the peak, but this difference is due to the impact of redshift space distortions on the correlation function and not due to scale-dependent bias in the sense that we mean it here. Note that this relative scale-independence does not need to hold in Fourier space. As a trivial example, has a scale-independent configuration-space bias if the transform of is arbitrarily small on the scales of interest. Taking , or the convolution of two halo profiles can satisfy this criterion. Depending on the size of this could lead to a large (but smoothly varying) scale-dependent bias in Fourier space (see e.g. the discussion in Schulz & White, 2006). Conversely, one finds that the Fourier transform of the Zel’dovich correlation function has almost no power beyond . To use the Zel’dovich approximation in Fourier space requires the addition of other terms which provide the missing power but affect the correlation function only at small scales.

Interestingly, the Zel’dovich approximation predicts that halos which are locally biased in Lagrangian space will have approximately the same small-scale quadrupole to large-scale quadrupole ratio as the matter, while the halos in N-body simulations display a significantly larger small-scale quadrupole when scaled to the same large-scale quadrupole. The term involving gives the desired increase in small-scale quadrupole, though at too small an amplitude. Making large and positive helps slightly, but the quadrupole still has the wrong overall shape.

Some authors (e.g. Melott, Pellman & Shandarin, 1994) have suggested that the Zel’dovich approximation performs better when small-scale power is filtered out of the linear theory power spectrum. We have tested this by Gaussian filtering the input on scales Mpc. We find that none of these scales improves the agreement of the quadrupole of the redshift-space correlation function on smaller scales.

At this point it is unclear whether the discrepancy above is due to our assumption of local Lagrangian bias or the Zel’dovich dynamics predicting the wrong velocity field for halos. To explore this issue further, we have run another set of 8 simulations with a simplified set-up. Each simulation started with initial conditions generated with the Zel’dovich approximation at (where the rms displacement was about 10 per cent of the inter-particle spacing). To ensure numerical convergence we smoothed the linear power spectrum with a Gaussian of Mpc. Again particles in a Gpc box were employed and for each particle, the value of the initial density, evaluated on the grid, was stored. The particles were integrated to either using a particle-mesh (PM) code (with a mesh) or the Zel’dovich approximation. Particles were then selected if their density field in the initial conditions (extrapolated to using linear theory) exceeded some threshold. In this manner the simulations mimic the analytic calculation closely.

In the analytic calculation we also used a linear theory power spectrum smoothed with a Mpc Gaussian, and we set assuming , i.e.

(49) | |||||

(50) |

(Szalay, 1988; Matsubara, 2011) where the limits shown are for and can be compared to the leading order behavior in Eqs. (46, 47) with . We have chosen such that the large-scale bias is approximately , as for the halo sample in Fig. 2.

Fig. 6 shows the real-space correlation function and the monopole and quadrupole moments of the redshift-space correlation function for all three methods, focusing on intermediate scales. The lower curves/points show the results for the matter field while the upper curves/points show the results for all particles with initial above a threshold. Note that there are no free parameters in this comparison! The agreement between the N-body results, the Zel’dovich simulations and the theory is excellent for the real-space correlation function and the monopole of the redshift-space correlation function. The agreement remains good (though not perfect) for the quadrupole of the matter field, but is less good for particles selected by initial density. In particular the PM results show the same qualitative difference from the Zel’dovich results as was found in the TreePM runs with halos. This suggests that the mismatch in the halo quadrupole that we are seeing in Figs. 2 and 3 is at least partly due to inadequacies in the Zel’dovich prediction for the inter-halo relative velocities (see also Seto & Yokoyama, 1998; Tassev & Zaldarriaga, 2012b, c), though some may be due to our assumption of local Lagrangian bias. We discuss non-local bias next.

## 6 Beyond local bias

The failure of our model to match the quadrupole moment on small and intermediate scales may be due in part to our assumption of local Lagrangian bias. While this approximation has received support from N-body simulations (Roth & Porciani, 2011; Baldauf et al., 2012; Chan, Scoccimarro, & Sheth, 2012; Wang & Szalay, 2012) it must break down at some level.

Perhaps the simplest modification to our formalism would be to allow a -dependence to the bias coefficients, . For example we could consider . For suitably chosen , such a modification can improve the agreement of the quadrupole at intermediate scales (Mpc) but it changes the monopole in a manner qualitatively similar to the overshoot seen in iPT in Fig. 2. In general a large enough modification to create agreement for the quadrupole removes the good agreement with the monopole. However it is possible to adjust so that the disagreement for the monopole is on small scales (Mpc) while the agreement for the quadrupole is improved non-negligbly over the range Mpc. Indeed, by “softening” the bias using a form like both the monopole and quadrupole can be made to agree with the N-body results to better than 5 per cent for Mpc, though the theory rapidly departs from the N-body results for smaller scales.

We also expect that terms involving e.g. the tidal tensor, can become important for high mass halos (Sheth, Chan & Scoccimarro, 2012). Such terms are naturally quadrupolar in nature and may affect the predictions. While shear terms are naturally induced by gravitational evolution, here we are interested in any dependence in the initial conditions.

Suppose we extend to also include a tidal-shear dependence, , as discussed for example by McDonald & Roy (2009)? We can Fourier transform on as well, generating a term . Since , with

(51) |

is already quadratic in , in the exponential it appears^{4}^{4}4When
using the cumulant theorem, bear in mind that is not Gaussian. only
multiplied by other expectation values.
Throughout we shall subtract the mean,
, from .
Working to lowest order in the additional terms to be added inside
the in Eq. (35) go as

(52) |

since the other terms involve the expectation value of an odd number of Gaussian fields. The relevant formulae can be found in the appendix. The terms are shear-density correlations: ; shear-shear terms: and ; shear-displacement terms: and cross-terms . These last terms include contributions with and into Eq. (37) and appear to be the best bet for influencing the quadrupole.

The terms involving are very small, as expected since they are higher order in . Even allowing for arbitrary prefactors in front of the terms, the halo quadrupole cannot be matched without spoiling the agreement with the redshift-space monopole and real-space correlation function. This is because the terms which enter contribute approximately the same amount to the monopole as to the quadrupole (as was the case with the non-shear terms). It thus appears that the lack of shear terms in the bias function is not the reason for the discrepancy seen in Fig. 2.

One may take a more general approach. Within the context of the Zel’dovich approximation, the terms appearing in the square brackets in Eq. (37) will be functions of and will be contracted with various factors of . There will be scalars, like , vectors, like , and tensors of various ranks. The vectors must be proportional to . The rank-2 tensors must go as a sum of terms like and , and similarly for higher rank objects. The most general biasing scheme would therefore consist of all such terms, with general dependence on . We have not undertaken an exploration of this large parameter space, but our experience above suggests that any such terms will contribute approximately equally to the monopole and the quadrupole moment, making it difficult to substantially adjust the quadrupole on small scales without spoiling the agreement seen for the monopole. As we saw above when we modified , it is possible to improve the level of agreement in some cases, though we do not posses a theory which predicts the required functional form at present.

On the basis of these investigations it appears that the disagreement between the Zel’dovich prediction for the quadrupole moment of the redshift-space, halo correlation function and that measured in N-body simulations may be due to simplifications inherent in the Zel’dovich approximation itself. Perhaps the relative velocities predicted from the lowest order displacement field are not as accurate for larger . Consistent with this view, we note that CLPT (which goes to next order in Lagrangian perturbation theory) does perform (very slightly) better than the Zel’dovich approximation in the range Mpc. Higher order LPT calculations are also known to give more accurate values for the moments (e.g. Munshi, Sahni & Starobinsky, 1994). However, the convergence appears to be very slow at best.

## 7 Zel’dovich streaming model (ZSM)

If the pure Zel’dovich calculation cannot match the small-scale quadrupole moment of the halo correlation function, can it be part of an extended model which can? The tests above suggest that the failure of the Zel’dovich approximation lies in the pairwise velocity distribution predicted by the model. Reid & White (2011) showed that the halo pairwise velocity distribution was quite well approximated by a Gaussian. What if we enforce this functional form (the “Gaussian streaming model”), using the Zel’dovich approximation to compute the ingredients: the real-space clustering of biased tracers, the mean infall velocity and the velocity dispersion? Specifically we assume

(53) |

with , and from the analytic theory.
This expression simply enforces pair counting, assuming that the functional
form of the velocity distribution is Gaussian, centered at ;
the mean LOS velocity between a pair of tracers as a function of their real
space separation.
We have just shown that the Zel’dovich approximation works well for the
real-space correlation function of biased tracers, such as halos.
The scale-dependence of the velocity dispersions^{5}^{5}5If our goal is to
model the clustering of galaxies, then we must include a phenomenological
model for the finger-of-god effect. Reid et al. (2012) showed that a single extra
parameter – an isotropic velocity dispersion – sufficed to model
fingers-of-god on large scales.
is well predicted by linear theory (Reid & White, 2011; Wang, Reid & White, 2013).
Thus the only missing ingredient is the mean infall velocity.

We can use the method described in Wang, Reid & White (2013) to compute within the Zel’dovich approximation. One simply adds a term to the exponent in Eq. (16) and computes derivatives with respect to , This is a subset of the calculation presented in Wang, Reid & White (2013):

(54) |

with^{6}^{6}6There is a typographical error in the subscripts in Eqs. (31, 32)
of Wang, Reid & White (2013) that is corrected here. The combination
should have been and
should have been . None of
the results in Wang, Reid & White (2013) were affected.

(55) | |||||