Estimating cosmic velocity fields from density fields and tidal tensors
Francisco-Shu Kitaura, Raul E. Angulo, Yehuda Hoffman, Stefan Gottl"ober
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 17 September 2018 (MN L A TEX style file v2.2)
Estimating cosmic velocity fields from density fields andtidal tensors
Francisco-Shu Kitaura , (cid:63) , Raul E. Angulo , Yehuda Hoffman and Stefan Gottl¨ober Leibniz-Institut f¨ur Astrophysik (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany Max-Planck Institut f¨ur Astrophysik (MPA), Karl-Schwarzschildstr. 1, D-85748 Garching, Germany Racah Institute of Physics, Hebrew University, Jerusalem, Israel
17 September 2018
ABSTRACT
In this work we investigate the nonlinear and nonlocal relation between cosmologicaldensity and peculiar velocity fields. Our goal is to provide an algorithm for the recon-struction of the nonlinear velocity field from the fully nonlinear density. We find thatincluding the gravitational tidal field tensor using second order Lagrangian perturba-tion theory (2LPT) based upon an estimate of the linear component of the nonlineardensity field significantly improves the estimate of the cosmic flow in comparison tolinear theory not only in the low density, but also and more dramatically in the highdensity regions. In particular we test two estimates of the linear component: the log-normal model and the iterative Lagrangian linearisation. The present approach relieson a rigorous higher order Lagrangian perturbation theory analysis which incorpo-rates a nonlocal relation. It does not require additional fitting from simulations beingin this sense parameter free, it is independent of statistical-geometrical optimisationand it is straightforward and efficient to compute. The method is demonstrated toyield an unbiased estimator of the velocity field on scales > ∼ h − Mpc with closelyGaussian distributed errors. Moreover, the statistics of the divergence of the peculiarvelocity field is extremely well recovered showing a good agreement with the true onefrom N -body simulations. The typical errors of about 10 km s − (1 sigma confidenceintervals) are reduced by more than 80% with respect to linear theory in the scalerange between 5 and 10 h − Mpc in high density regions ( δ >
Key words: (cosmology:) large-scale structure of Universe – galaxies: clusters: gen-eral – catalogues – galaxies: statistics
Gravitational instability is one of the key ingredients to ex-plain the rich hierarchy of structures we observe today in theUniverse. It has amplified small mass fluctuations producedafter inflation to give rise from large galaxy clusters to hugevoids. Simultaneously, the same local gravitational field im-printed “peculiar” velocities in galaxies; deviations from theoverall expansion of the Universe which is responsible forthe Hubble flow.The peculiar velocity of galaxies is a valuable quantityin cosmology since it contains similar but complementaryinformation to that enclosed in the galaxies position. Forinstance, by requiring isotropy in the measured galaxy clus-tering, the cosmological mass density parameter and even (cid:63)
E-mail: [email protected], Karl-Schwarzschild fellow the nature of gravity can be constrained (see e.g. Daviset al. 1996; Willick & Strauss 1998; Branchini et al. 2001;Guzzo et al. 2008). In addition, these motions can be usedto reconstruct the properties of the universe at an earliertime, in principle, even at recombination where perturba-tions were completely linear (Nusser & Dekel 1992; Gra-mann 1993a; Croft & Gaztanaga 1997; Narayanan & Wein-berg 1998; Monaco & Efstathiou 1999; Frisch et al. 2002).A method able to accurately determine the peculiarvelocity field can be used in many different applications;ranging from bias studies combining galaxy redshift surveyswith measured peculiar velocities (see e.g. Fisher et al. 1995;Zaroubi et al. 1999; Courtois et al. 2011), Baryon AcousticOscillations reconstructions (see e.g. Eisenstein et al. 2007),determination of the growth factor, to estimates of the initialconditions of the Universe which in turn can be used for con-strained simulations (see e.g. Gottloeber et al. 2010; Klypin c (cid:13) a r X i v : . [ a s t r o - ph . C O ] J un Kitaura et al.
Works θ - δ relation parameters methodology Yahil et al. (1991) θ = δ { p } linear theory: LPTNusser et al. (1991) θ = δ/ (1 + αδ ) α , { p } empirical approximationBernardeau (1992) θ = α (cid:104) (1 + δ ) β − (cid:105) α, β , { p } PTGramann (1993b) θ = δ − αD µ (2) ( φ g ) α, { p } approx. 2LPTWillick et al. (1997) θ = (cid:2) (1 + α σ δ ) δ + ασ δ (cid:3) / (1 + αδ ) α, σ δ , { p } empirical approximationChodorowski et al. (1998) θ = (cid:2) δ − α (cid:0) δ − σ δ (cid:1) + βδ (cid:3) α, β, σ δ , { p } PT+ N -bodyBernardeau et al. (1999) θ = (cid:2) γδ + α (cid:0) δ − σ δ (cid:1) + βδ (cid:3) α, β, γ, σ δ , { p } PT+ N -bodyKudlicki et al. (2000) θ = α (cid:104) (1 + δ ) /α − (cid:105) + ( α − / (2 α ) σ δ α, σ δ , { p } PT+Eulerian grid-based codeBilicki & Chodorowski (2008) θ = γ (cid:104) (1 + δ ) /α − (1 + δ ) β (cid:105) α, β, γ , { p } spherical collapse modelthis work θ = D δ (1) − D f /f µ (2) ( φ (1) ) { p } Table 1.
The parameters α, β, γ are either derived from first principles (Bernardeau 1992; Gramann 1993b; Bilicki & Chodorowski2008, this work) or from fitting to simulations being different for each case. These parameters also depend on the scale of interest. Theparameters which are in parenthesis are fully determined by the chosen cosmology and the theoretical model. The variance of the densityfield is given by σ δ = (cid:104) δ (cid:105) . PT stands for perturbation theory and is applied only in the univariate case (local relation). 2LPT standsfor second order Lagrangian perturbation theory and yields the only nonlocal (and nonlinear) relation from the list. Hivon et al. (1995);Monaco & Efstathiou (1999) proposed 2LPT to correct for redshift distortions. Let us mention here the least-action principle methodsto determine the peculiar velocities from mass tracer objects (galaxies) introduced by Peebles (1989) and further extended by Nusser &Branchini (2000); Branchini et al. (2002); Lavaux et al. (2008). et al. 2003). A particularly well suited application regardsthe topological methods to detect the kinematic Sunyaev-Zeldovich effect.There is in addition recent interest in the measurementof large-scale flows. Some authors claim to have detecteda so-called “dark flow” caused by super-horizont perturba-tions (see e.g. Kashlinsky et al. 2009, 2011). Others discusssuch flows as a challenge to the standard cosmological modelas a whole (see e.g. Watkins et al. 2009).Unfortunately, the apparent shift in spectral features ofa galaxy is also affected by the expansion of the universe,therefore it is not possible to directly measure the peculiarvelocity field in spectroscopic surveys. For this reason, onehas to resort to indirectly infer it from the mass density fluc-tuations (but see Nusser et al. 2011, for a recent alternativemethod). However, this is not a trivial procedure due thehighly nonlinear state of the density fluctuations today anddue to its nonlocal relationship with the velocity field.The simplest approach is given by the linear continu-ity equation, which is routinely used in clustering stud-ies. However, it has a range of applicability only limitedto very large scales (e.g. Angulo et al. 2008). Alternativemethods devised to improve upon linear theory can be sep-arated into three areas. The first one consists on developingnonlinear relationships with higher-order perturbation the-ory (Bernardeau 1992; Chodorowski et al. 1998; Bernardeauet al. 1999; Kudlicki et al. 2000), with spherical collapsemodels (Bilicki & Chodorowski 2008) or based on empiricalrelations found in cosmological N-body simulations (Nusseret al. 1991; Willick et al. 1997).Another idea is to solve the boundary problem of find-ing the initial positions of a set of particles governed by theEulerian equation of motion and gravity with the least ac-tion principle (see Peebles 1989; Nusser & Branchini 2000;Branchini et al. 2002). A similar approach consists on relat-ing the observed positions of matter tracers (e.g. galaxies)in a geometrical way to a homogeneous distribution by min-imizing a cost function, which combined with the Zeldovichapproximation (Zel’dovich 1970) provides an estimation of the velocity field (see Mohayaee & Tully 2005; Lavaux et al.2008).The diversity of strategies and approximations for ob-taining the velocity from the density field hint at the diffi-culty of the problem. Some approaches are simply not accu-rate enough and some are computationally very expensive.This sets the agenda for an improvement in the field. Anynew method should be accurate, unbiased, computationallyefficient and applicable to observational data.A further shortcoming of the existing approaches is thatthey are mostly particle-based, which is not applicable formore general matter tracers like the Lyman-alpha forest orthe 21-cm line, nor they can be combined with optimal den-sity field estimators (see Kitaura et al. 2010; Jasche & Ki-taura 2009; Kitaura et al. 2010).In this paper we investigate a different approach basedon higher order Lagrangian perturbation theory, and it ismotivated by the pioneering work of Gramann (1993b) andfurther extended by Hivon et al. (1995); Monaco & Ef-stathiou (1999). The theoretical basis for LPT was care-fully worked out by Moutarde et al. (1991); Buchert &Ehlers (1993); Buchert (1994); Bouchet et al. (1995); Cate-lan (1995) (for further references see Bernardeau et al. 2002).Contrary to classic applications of LPT, in which theproperties of an evolved distribution are predicted from alinear density field in Lagrangian coordinates (e.g. in thegeneration of initial conditions for N -body simulations orof galaxy mock catalogues: Jenkins (2010); Scoccimarro &Sheth (2002)), our starting point is an evolved field in Eule-rian coordinates (e.g. the present-day galaxy distribution).The key realisation of our approach is that it is possibleto decompose a nonlinear density field into a Gaussian andNon-Gaussian term, which are related to each other throughLPT (see similar approaches Gramann 1993a,b; Monaco &Efstathiou 1999). In other words, it is possible to find aclosely Gaussian field which would evolve, under the as-sumption of LPT, into the measured density field today.This Gaussian field can then naturally be used to predictthe corresponding velocity field today in LPT.Our method combines i) the equations of motion and c (cid:13) , 000–000 osmic velocity fields continuity for a fluid under self gravity in 2-nd (3-rd) orderLagrangian Perturbation theory (LPT) with ii) the idea thatthe present-day galaxy distribution can be expanded into aclosely linear-Gaussian field and a highly skewed nonlinearcomponent consistent with 2LPT (or 3LPT) (see Kitaura& Angulo 2011). The former aspect makes our approachphysically motivated and also captures the nonlocal natureof the density-velocity relation via the gravitational tidalfield tensor. The latter aspect self-consistently minimises theimpact of the approximations of a 2nd order formulation ofgravitational evolution, but more importantly, it enables theapplication of LPT to nonlinear fields.We note that the use of the lognormal transformation(including the subtraction of a mean field) to obtain an esti-mate of the linear field was proposed by Kitaura et al. (2010)and has been recently confirmed to give already a good esti-mate of the divergence of the displacement field (Falck et al.2011). To estimate the velocity field one needs, however, togo to higher order perturbation theory as we show in thiswork.In the next section we recap Lagrangian perturbationtheory and derive the velocity-density relation to second andthird order. In section 2.2 we will present the method to com-pute the peculiar velocity field from the nonlinear densityfield. In section 3 we will present our numerical tests basedon the Millennium Run simulation. Here we will analyze thegoodness of the recovered velocity divergence as comparedto the true one and the same for each velocity component. Astudy of the errors in our method is also presented. Finallywe present our conclusions. The first part of this section presents the relation betweendensity and velocity fields in 2LPT, and how it can be ap-plied to an evolved field. In the second part, we outline apractical implementation of this method.
The basic idea in Lagrangian perturbation theory is thatan initial homogeneous field expressed in Lagrangian coordi-nates q can be related to a final field in Eulerian coordinates x trough a unique mapping: x = q + Ψ ( q ) determined by adisplacement field Ψ (see e.g. Bernardeau et al. 2002).The linear solution for this expression is the well-knownZeldovich approximation, in which the displacement field isgiven by the Laplacian of the gravitational potential at q .This result has been successfully applied to many aspects ofcosmology, but it fails to describe the dynamics of a non-linear field where shell crossings and curved trajectories arecommon.An improvement is found by expanding the displace-ment field and considering higher order terms (see e.g.Buchert et al. 1994; Melott et al. 1995; Bouchet et al. 1995).For instance, the displacement field to second order is givenby Ψ ( q ) = − D ∇ q φ (1) ( q ) + D ∇ q φ (2) ( q ) , (1) where D is the linear growth factor normalised to unity to-day, D the second order growth factor, which is given by D (cid:39) − / − / D (see Buchert & Ehlers 1993; Bouchetet al. 1995). The potentials φ (1) ( q ) and φ (2) ( q ) are obtainedby solving a pair of Poisson equations: ∇ φ (1) ( q ) = δ (1) ( q ),where δ (1) is the linear over-density, and ∇ φ (2) ( q ) = δ (2) ( q ).It is important to realise that these terms are not in-dependent of each other. The second order nonlinear term δ (2) is fully determined by the linear over-density field δ (1) through the following quadratic expression δ (2) ( q ) ≡ µ (2) ( φ (1) ( q )) = (cid:88) i>j (cid:16) φ (1) ( q ) ,ii φ (1) ( q ) ,jj − [ φ (1) ( q ) ,ij ] (cid:17) , (2)where we use the following notation φ ,ij ( q ) ≡ ∂ φ ( q ) /∂q i ∂q j , and the indices i, j run over the threeCartesian coordinates.Similarly, the particle co-moving velocities v are givento second order by: v ( q ) = − D f H ∇ q φ (1) ( q ) + D f H ∇ q φ (2) ( q ) , (3)where f i = d ln D i / d ln a , H is the Hubble constant and a isthe scale factor. For flat models with a non-zero cosmolog-ical constant, the following relations apply f ≈ Ω / , f ≈ / (see Bouchet et al. 1995), where Ω( z ) is the matterdensity at a redshift z .We note that going to third order in Lagrangian pertur-bation theory provides modest improvements (see Buchert& Ehlers 1993; Melott et al. 1995; Catelan 1995; Bouchetet al. 1995; Bernardeau et al. 2002).To apply the Lagrangian framework to a density fieldin the Eulerian frame δ M ( x ) one must be careful. Under theassumption that there is no shell-crossing, one can write theinverse equation relating Eulerian to Lagrangian coordinates q = x − Ψ ( x ). Mass conservation leads to the followingequation (see Nusser et al. 1991):1 + δ M ( x ) = ˜ J ( x ) , (4)with ˜ J ( x ) ≡ (cid:12)(cid:12)(cid:12)(cid:12) ∂ q ∂ x (cid:12)(cid:12)(cid:12)(cid:12) . (5)Expanding the Jacobian ˜ J ( x ) one finds (Kitaura & An-gulo 2011); δ M ( x ) = | − ∇ x · Ψ ( x ) | − (cid:39) −∇ x · Ψ ( x ) + µ (2) (Θ( x )) + µ (3) (Θ( x )) . (6)The displacement or velocity field derived from this ex-pression will automatically be expressed in Eulerian coor-dinates as we will show below. We note that the same ex-pression is found as a function of the Lagrangian coordinate q when expanding the inverse of the corresponding Jacobian J ( q ) ≡ | ∂ x /∂ q | (see Kitaura & Angulo 2011).Using the displacement field (Eq. 1) in Eulerian coordi-nates, the final density field can be written in terms of thelinear and nonlinear fields: c (cid:13)000
The basic idea in Lagrangian perturbation theory is thatan initial homogeneous field expressed in Lagrangian coordi-nates q can be related to a final field in Eulerian coordinates x trough a unique mapping: x = q + Ψ ( q ) determined by adisplacement field Ψ (see e.g. Bernardeau et al. 2002).The linear solution for this expression is the well-knownZeldovich approximation, in which the displacement field isgiven by the Laplacian of the gravitational potential at q .This result has been successfully applied to many aspects ofcosmology, but it fails to describe the dynamics of a non-linear field where shell crossings and curved trajectories arecommon.An improvement is found by expanding the displace-ment field and considering higher order terms (see e.g.Buchert et al. 1994; Melott et al. 1995; Bouchet et al. 1995).For instance, the displacement field to second order is givenby Ψ ( q ) = − D ∇ q φ (1) ( q ) + D ∇ q φ (2) ( q ) , (1) where D is the linear growth factor normalised to unity to-day, D the second order growth factor, which is given by D (cid:39) − / − / D (see Buchert & Ehlers 1993; Bouchetet al. 1995). The potentials φ (1) ( q ) and φ (2) ( q ) are obtainedby solving a pair of Poisson equations: ∇ φ (1) ( q ) = δ (1) ( q ),where δ (1) is the linear over-density, and ∇ φ (2) ( q ) = δ (2) ( q ).It is important to realise that these terms are not in-dependent of each other. The second order nonlinear term δ (2) is fully determined by the linear over-density field δ (1) through the following quadratic expression δ (2) ( q ) ≡ µ (2) ( φ (1) ( q )) = (cid:88) i>j (cid:16) φ (1) ( q ) ,ii φ (1) ( q ) ,jj − [ φ (1) ( q ) ,ij ] (cid:17) , (2)where we use the following notation φ ,ij ( q ) ≡ ∂ φ ( q ) /∂q i ∂q j , and the indices i, j run over the threeCartesian coordinates.Similarly, the particle co-moving velocities v are givento second order by: v ( q ) = − D f H ∇ q φ (1) ( q ) + D f H ∇ q φ (2) ( q ) , (3)where f i = d ln D i / d ln a , H is the Hubble constant and a isthe scale factor. For flat models with a non-zero cosmolog-ical constant, the following relations apply f ≈ Ω / , f ≈ / (see Bouchet et al. 1995), where Ω( z ) is the matterdensity at a redshift z .We note that going to third order in Lagrangian pertur-bation theory provides modest improvements (see Buchert& Ehlers 1993; Melott et al. 1995; Catelan 1995; Bouchetet al. 1995; Bernardeau et al. 2002).To apply the Lagrangian framework to a density fieldin the Eulerian frame δ M ( x ) one must be careful. Under theassumption that there is no shell-crossing, one can write theinverse equation relating Eulerian to Lagrangian coordinates q = x − Ψ ( x ). Mass conservation leads to the followingequation (see Nusser et al. 1991):1 + δ M ( x ) = ˜ J ( x ) , (4)with ˜ J ( x ) ≡ (cid:12)(cid:12)(cid:12)(cid:12) ∂ q ∂ x (cid:12)(cid:12)(cid:12)(cid:12) . (5)Expanding the Jacobian ˜ J ( x ) one finds (Kitaura & An-gulo 2011); δ M ( x ) = | − ∇ x · Ψ ( x ) | − (cid:39) −∇ x · Ψ ( x ) + µ (2) (Θ( x )) + µ (3) (Θ( x )) . (6)The displacement or velocity field derived from this ex-pression will automatically be expressed in Eulerian coor-dinates as we will show below. We note that the same ex-pression is found as a function of the Lagrangian coordinate q when expanding the inverse of the corresponding Jacobian J ( q ) ≡ | ∂ x /∂ q | (see Kitaura & Angulo 2011).Using the displacement field (Eq. 1) in Eulerian coordi-nates, the final density field can be written in terms of thelinear and nonlinear fields: c (cid:13)000 , 000–000 Kitaura et al. -5 0 5 10 1510 -6 -5 -4 -3 -2 -1 r S =5 h − Mpc P D F θ k NbodyLINGRAMLOGLOG-2LPT2LPT -2 0 2 4 610 -6 -5 -4 -3 -2 -1 r S =10 h − Mpc θ k Figure 1.
Statistics of the scaled velocity divergence ( θ k : true from the simulation θ Nbody and different reconstructed ones θ rec asexplained below) with different smoothing scales: Left: r S = 5 h − Mpc, Right: 10 h − Mpc. Curves from left to right in the order ofappearance: yellow: approximate 2LPT solution from the nonlinear density field (GRAM) (see G93 in table 1), cyan: logarithm of thedensity field (LOG), blue dashed: true scaled velocity divergence at z = 0 the Millennium Run (Nbody), black: 2LPT estimate fromthe logarithm of the density field (LOG-2LPT) (Eq. 11), red: 2LPT estimate from the iterative solution (2LPT) (see § δ M ( x )= D δ (1) ( x ) − D µ (2) ( φ (1) ( x )) + µ (2) (Θ( x )) + µ (3) (Θ( x ))= δ L ( x ) + δ NL ( x ) , (7)with Θ( x ) being the potential associated to the divergenceof the displacement field: Θ( x ) = −∇ − ∇ · Ψ ( x ), δ L ( x ) = D δ (1) ( x ) and δ NL ( x ) = − D µ (2) ( φ (1) ( x )) + µ (2) (Θ( x )) + µ (3) (Θ( x )). The third order term in the Jacobian expansionis given by: µ (3) (Θ( x )) = det (Θ( x ) ,ij ) . (8)From now on the Eulerian coordinate dependence ( x ) is leftout. Assuming that any primordial vorticity has no growingmodes associated (the first vorticity terms appear in 3rd or-der Lagrangian perturbation theory, see appendix) impliesthat the velocity field today is fully characterised by its di-vergence ( ∇ · v ), or, for convenience, by the scaled velocitydivergence, defined as: θ ≡ − f (Ω , Λ , z ) − ∇ · v , (9)with f (Ω , Λ , z ) ≡ ˙ D /D = f (Ω , Λ , z ) H ( z ).Combining Eqs. 7 and 3 truncated to quadratic termsin Φ (1) : δ NL = (cid:0) D − D (cid:1) δ (2) , one gets θ = δ M − (cid:20) D + (cid:18) f f − (cid:19) D (cid:21) δ (2) . (10)This expression is very similar to the one found by Gra-mann (1993b) (see Tab. 1). Note, however, that using di-rectly the evolved field as the source for the second orderterm δ (2) is a good approximation for the true velocity field only on very large scales (where δ M is close to unity), as wewill see in § h − Mpc for both estimations of the nonlinear field and ofthe velocity divergence.In this paper we follow a different approach, and solveiteratively the following equation; θ = D δ (1) − D f f δ (2) . (11)which results from taking the divergence of Eq. 3. For this,we rely on an estimation of the linear component of thepresent-day density field, which in turn can be calculatedby solving iteratively Eq. (7). Note that the Eulerian de-scription forces one to expand the inverse of the Jacobian(see Eqs. 4-7 and Kitaura & Angulo 2011). This is the maindifference with respect to the work of Monaco & Efstathiou(1999) in which first the Jacobian is expanded and then theinverse of it is taken and could explain why we avoid prob-lems caused by artificial Lagrangian caustics in low densityregions as reported in their work.In practice, a good approximation for the linear term, δ (1) , is simply given by the logarithm of the density field; δ L = D δ (1) = ln(1 + δ M ) − µ , with µ = (cid:104) ln(1 + δ M ) (cid:105) , asshown by (Neyrinck et al. 2009; Kitaura & Angulo 2011).Note that this expression is essentially the lognormal ap-proximation for the matter field (Coles & Jones 1991). Thislocal transformation has the advantage of being computa-tionally inexpensive.In summary, approach finds a linear field which whenplugged into 2LPT expressions produces the observed mat-ter field (or third order, see appendix). If gravity workedonly at a second order level, then this linear field would beidentical to the actual linear field that originated δ M , but c (cid:13) , 000–000 osmic velocity fields naturally in reality it is just an approximation. Thus, it isimportant to characterise the performance and accuracy ofthe method, which we do in §
3. But first, we discuss a prac-tical implementation of our method in the next subsection.
The method to estimate the peculiar velocity field from thenonlinear density field is straightforward and fast to com-pute. We now outline the steps to be followed for its im-plementation. For this, we have assumed that there is anunbiased and complete estimation of the matter field δ M .The extra layer of complication introduced by shot noise,a survey mask, biasing and redshift space distortions is outof the scope of this paper and will be addressed in a futurepublication.(i) Linear density field
One starts by computing an estimate of the linear com-ponent of the density field. We propose two alternatives forthis:(a) Lognormal model δ LLOG = D δ (1) = ln(1 + δ M ) − µ (12)(b) Gaussianisation with LPT (Kitaura & Angulo2011) δ LLPT = D δ (1) = δ M − δ NL (13)(ii) Linear potential
Then the Poisson equation is solved to obtain the linearpotential: φ (1) = ∇ − δ (1) (14)(iii) Nonlinear second order density field
The tidal field contribution to second order is calculatedfrom the linear potential: δ (2) (cid:104) φ (1) (cid:105) = (cid:88) i>j (cid:16) φ (1) ,ii φ (1) ,jj − ( φ (1) ,ij ) (cid:17) , (15)(iv) Scaled velocity divergence
One can now construct the 2nd order divergence of thevelocity field taking both the linear and the second ordercontribution: θ = D δ (1) − D f f δ (2) , (16)(v) Peculiar velocity field
Finally, one obtains the 3D velocity field: v = − f (Ω , Λ , z ) ∇∇ − θ . (17)Please note that the equations in steps (ib), (ii), (iii) and(v) can be solved with FFTs. Details of the Gaussianisationstep with LPT can be found in Kitaura & Angulo (2011). σ ( ∇ × v | x ) / σ ( ∇ · v ) [ % ] r S [ h − Mpc]
Figure 5.
Standard deviation of the x -component of the curl σ ( ∇ × v | x ) divided by the standard deviation of the velocity di-vergence σ ( ∇ · v ) in % as a function of the scale r S . The meanof both the curl and the divergence of the peculiar field are veryclose to zero. N -BODYSIMULATIONS In this section we test the performance of the method out-lined above by comparing the velocity field directly ex-tracted from a N -body simulation with our estimation basedon the respective nonlinear density field.With this purpose, we employ the Millennium simu-lation which tracks the nonlinear evolution of more than10 billion particles in a box of comoving side-length 500 h − Mpc (Springel et al. 2005). In particular, we use theoutput corresponding to redshift 0, which is where the mostnonlinear structures are present.At such output we first compute the velocity and den-sity field by mapping the information of dark matter par-ticles onto a grid of 256 cells using the nearest-grid-point(NGP) assignment scheme, which gives a spatial resolutionof about 2 h − Mpc. We then apply the algorithm presentedin § h − Mpc on which linear theory is usually assumed toperform well, and 5 h − Mpc which enters into the mildlynonlinear regime.
The first ingredient in our algorithm is the linear compo-nent of the density field. We stress that this field are not the“initial conditions” of the universe, since structures havemoved from its Lagrangian position to the Eulerian ones at z = 0. Contrarily to the linear term, the probability distri-bution function (PDF) of the nonlinear component is highlyskewed.Fig. 1 compares the PDF of the velocity divergence θ as given by different estimators, with the one directly mea- c (cid:13)000
The first ingredient in our algorithm is the linear compo-nent of the density field. We stress that this field are not the“initial conditions” of the universe, since structures havemoved from its Lagrangian position to the Eulerian ones at z = 0. Contrarily to the linear term, the probability distri-bution function (PDF) of the nonlinear component is highlyskewed.Fig. 1 compares the PDF of the velocity divergence θ as given by different estimators, with the one directly mea- c (cid:13)000 , 000–000 Kitaura et al. r S =5 h − MpcLIN θ r ec L I N r S =10 h − MpcLIN θ r ec L O G − L P T r S =5 h − MpcLOG-2LPT r S =10 h − MpcLOG-2LPT θ r ec L P T r S =5 h − Mpc2LPT θ Nbody r S =10 h − Mpc2LPT θ Nbody
Figure 2.
Cell-to-cell comparison between the true and the reconstructed scaled velocity divergence, θ Nbody and θ rec respectively. Leftpanels show results on scales of 5 h − Mpc and right panels on scales of 10 h − Mpc. Upper panels: LIN, middle panels: LOG-2LPT,lower panels: 2LPT. The dark colour-code indicates a high number and the light colour-code a low number of cells.c (cid:13) , 000–000 osmic velocity fields -0.5 0.0 0.5 1.0 1.5 2.0 Y [ h − M p c ] -0.1 -0.0 0.1 0.2 0.3 0.4 0.5 Y [ h − M p c ] Y [ h − M p c ] X [ h − Mpc] X [ h − Mpc]
Figure 3.
Difference between the true and the estimated scaled velocity divergence. On scales of r S =5 h − Mpc (left panels) and 10 h − Mpc (right panels) with different colour ranges. Upper panels: LIN, middle panels: LOG-2LPT, lower panels: 2LPT.c (cid:13)000
Difference between the true and the estimated scaled velocity divergence. On scales of r S =5 h − Mpc (left panels) and 10 h − Mpc (right panels) with different colour ranges. Upper panels: LIN, middle panels: LOG-2LPT, lower panels: 2LPT.c (cid:13)000 , 000–000
Kitaura et al. -0.5 0.0 0.5 1.0 1.5 2.0 Y [ h − M p c ] -0.5 0.0 0.5 1.0 Y [ h − M p c ] Y [ h − M p c ] X [ h − Mpc] X [ h − Mpc]
Figure 4.
Difference between the true and the estimated speeds (velocity magnitude). On scales of r S =5 h − Mpc (left panels) and 10 h − Mpc (right panels) with different colour ranges. Upper panels: LIN, middle panels: LOG-2LPT, lower panels: 2LPT. The colour-codeis in units of 10 − s − km. c (cid:13) , 000–000 osmic velocity fields sured from the simulation (blue dashed line). The left panelshows the fields smoothed with a 5 h − Mpc whereas theright panels do so with a larger smoothing of 10 h − Mpc.In both cases, the predictions of linear perturbation theory(i.e. θ = δ ) displays the worst performance of all – it overes-timates systematically the number of volume elements withlarge value of θ and underestimates the ones with low val-ues of θ (green curve). This is a consequence of linear theorybreaking down even in mildly under- or overdense regions.Using linear theory together with the linear componentof the nonlinear density field given by the logarithm of thefield (ln(1 + δ M ) − µ , see cyan curve in Fig. 1) is still apoor description of the PDF, suggesting the need of higherorder corrections. We note that this approximation corre-sponds to linear Lagrangian perturbation theory (Zeldovichapproximation).While standard linear (Eulerian) theory overestimatesthe values of the divergence of the peculiar velocity field, lin-ear Lagrangian perturbation theory underestimates them.In the former case we are assuming that the scaled diver-gence of the peculiar velocity field is given by the densityfield at (final) Eulerian coordinates ( θ = δ ( x )), i. e. by thefull nonlinear density field; whereas in the latter case weassume that it is given by the density field at (initial) La-grangian coordinates transformed to Eulerian coordinates( θ = D δ (1) ( x )), i. e. by the first linear (Gaussian) termin Lagrangian perturbation theory. We note that in generalthe underestimation of the Zeldovich approximation is lesssevere (cyan curve) than the overestimation of linear theory(green curve) in high density regions (see Fig. 1). Meaningthat the Zeldovich approximation performs better than lin-ear theory being more conservative, as it has been repeatedlyshown in the literature (see e.g. Nusser and Branchini 2000).Velocity reconstruction approaches like PIZA or MAK arebased on the Zeldovich approximation (see Croft & Gaz-tanaga 1997; Lavaux et al. 2008, respectively). We thus ex-pect that the approach presented here is more accurate thanthese methods.The first second order estimation we consider is thatgiven by Eq. (10) which is closely related to the one pro-posed by Gramann and uses the nonlinear field as a proxy forthe linear density field (yellow curve). In the regime wherethis approximation is valid ( δ ∼
0) this approach performsremarkably well. However, there is a clear and rapid degra-dation for volume elements with larger deviations of homo-geneity. For instance, this solution yields to values even of θ = −
140 at scales of 5 h − Mpc! This behaviour is dueto a complete misestimation of the nonlinear term δ (2) andtherefore of the nonlinear corrections to the velocity field.Finally, black and red lines indicate the results of themethod proposed in this paper: Lagrangian perturbationtheory based upon an estimate of the linear component ofthe density field using the lognormal model (LOG-2LPT:black) or based on the iterative Lagragian linearisation(2LPT: red). In both panels, the predicted PDF very closelyfollows that measured in the Millennium simulation, evenon the extreme tails (especially with LOG-2LPT). The onlyappreciable difference with the lognormal model is a slightoverestimation for low values of θ (static regions), we notehowever, that this could be potentially improved by higherorder expressions. Indeed, 3rd order PT appears to performbetter than 2LPT for underdense regions (see appendix). In spite of this, on both 5 and 10 h − Mpc our method isclearly superior to any of the other methods we investigatedhere, as far as the PDF is concerned and for any value of θ . The iterative 2LPT solution yields a moderate overesti-mation for high values of θ due to the laminar flow approxi-mation used in Lagrangian perturbation theory which doesnot fully capture nonlinear structure formation. Neverthe-less, the PDF of θ using this solution is clearly superior tothe linear approximation.We now continue with a more detailed testing of ourmethod. In Fig. 2 we plot the predicted velocity divergence, θ recLOG − and θ rec2LPT (based on δ LLOG and δ LLPT , respec-tively), for each of the 256 cells in our mesh, as a functionof the value measured in the simulation. As in the previousplot, we display results on two different scales; 5 h − Mpc onthe left panels and 10 h − Mpc on the right panels. For com-parison we also provide the results using linear perturbationtheory θ recLIN = δ .The values of θ are remarkably well predicted by La-grangian perturbation theory. In fact, measurements liearound the 1:1 line in the LOG-2LPT case, implying thatthere are no appreciable biases in our estimation over all therange probed by the Millennium simulation (with the excep-tion of the low values of θ , for an improvement on this seeappendix). In contrast, the linear theory prediction presentsoverestimations of up to a factor of 3 for the 5 h − Mpcsmoothing and of 2 for 10 h − Mpc.The iterative solution (2LPT) produces smaller disper-sions but also a slight overestimation of θ for high values, aswe already mentioned before.The distribution of differences in our method is well ap-proximated by a Gaussian function, whereas in linear theorythere are significant extended tails, we will return to thisin more detail in § θ in agiven volume element.Although not displayed by the figure, our method alsoperforms better than the other methods shown in Fig. 1. Inparticular, the classic application of 2LPT recovers θ quitewell for the range − . < θ <
2, as shown in Gramann(1993b). However, outside this range it displays an erraticbehaviour as it could have been anticipated from Fig. 1.A complementary visual assessment of our method isprovided in Fig. 4. In these images we display the relativedifference θ rec − θ Nbody ( θ recLIN − θ Nbody in the upper panels, θ rec2LPT − θ Nbody in the central panels and θ recLOG − − θ Nbody in the lower ones) projected on a slice 2 h − Mpc thick. Aspreviously, we explore 5 and 10 h − Mpc and also displaythe linear theory predictions for comparison. We see thatthis difference field is not uniformly distributed across thesimulation but there are well defined regions in which theprediction is very accurate and others where the predictionis somewhat worse. Not surprisingly the latter coincide withhigh density regions. Nevertheless, and consistently withprevious plots, we see that areas where linear theory failsdramatically are much better handled in our approach.As a final crucial check, we compute the power-spectraof the scaled velocity divergence according to the N -body We have also checked that this is true on 3,4,6,7 and 8 h − Mpcc (cid:13) , 000–000 Kitaura et al. r S =10 h − Mpc (dashed) r S = 5 h − Mpc (solid)
LINLOG-2LPT2LPT k [ h Mpc − ] P r ec θ ( k ) / P N b o d y θ ( k ) Figure 6.
Ratio of the true P Nbody θ ( k ) and reconstructed P rec θ ( k ) power-spectra of the scaled velocity divergence θ forgreen: linear theory (LIN), black: 2LPT estimate from the log-arithm of the density field (LOG-2LPT) (Eq. 11), red: 2LPT es-timate from the iterative solution (2LPT) (see § h − Mpc smoothing we have multiplied the linesby a factor of 2 for clarity. simulation and our reconstruction with both the LPT andthe lognormal linearisation approaches. This is shown inFig. 6. Linear theory, as expected, increasingly overestimatesthe power towards small scales, whereas the 2LPT solutionspeform remarkably well in a wide range going down to scalesof k (cid:39) . h /Mpc and k (cid:39) . h /Mpc for r S =10 h − Mpcand r S =5 h − Mpc, respectively. We can also see that thereis a systematic deviation originated by the lognormal trans-formation. The LPT estimate of the linear component cor-rects these and the results are extremely close to the actualpower-spectrum over most of the k -range shown. We have calculated each component of the peculiar velocityfield ( v x , v y , v z ) from the inferred velocity divergence as-suming ∇ × θ = 0. In this case, the Fourier transform of thevelocity along a direction is given by ˆ v = k · ˆ θ with k beingthe k-vector. The approximation that the velocity field is ir-rotational is actually a good one for the scales and redshiftswe consider here. In fact, the curl of the velocity fields is onaverage more than 25 times smaller than its divergence inthe z=0 output of the Millennium simulation on scales of 3 h − Mpc and even smaller on larger scales (see Fig. 5).In Fig. 7 we compare the velocity field along the x-axisin the simulation and in our method (the two other axisshow essentially the same features). The contours show thevariation of number of cells (the darker colours representhigher numbers). By comparing both columns, it is clearthat linear theory presents biases for large speeds, which areremoved by our 2LPT approach, on both 5 and 10 h − Mpc.We can see that the distribution gets sharper with 2LPTshowing a clear decrease in the number of outliers. We quantified the uncertainty in the recovered velocityin Fig. 8. The x-axis corresponds to the errors in the re-construction, defined as (cid:15) v,x ≡ v rec x − v x . As in the previousplot, we show only results for the x-axis since the other threeCartesian coordinates provide consistent results.We find that the errors in our method are closely Gaus-sian distributed – the skewness and kurtosis are dramati-cally smaller than for linear theory. This property is veryimportant when applying the method to real data, since theobservational uncertainties can be added to the model un-certainties within a Gaussian likelihood without the needof introducing complex error models. The typical errors arealso significantly reduced with smaller standard deviations.The errors in the reconstruction using linear theory havevery long tails. Such outliers are not present in our method(see Fig. 8).As we have discussed along this paper, the larger im-provement of our method concerns velocities in high densityregions. For regions with δ > δ > h − Mpcand cells with δ > − and are reduced with 2LPT to ∼ − , i.e. 81% smaller. The corresponding 2 sigma confi-dence intervals are about 160 and 28 km s − , respectively,i.e. an error reduction of about 83%. One can see that the 2sigma confidence intervals are about double as large as the1 sigma confidence intervals for the 2LPT estimation. How-ever, this is not the case for the linear estimates as these arenot Gaussian distributed. In this paper we presented an improved method to recon-struct the peculiar velocity field from the density field. Itbuilds from 2nd order Lagrangian perturbation theory andthe realisation that the density field can be split into a linearplus a nonlinear term. The latter is the key concept, whichenables the application of Lagrangian perturbation theoryto an evolved field in Eulerian coordinates. This in turn, cre-ates an approach that is nonlinear and nonlocal by includingthe information of the gravitational tidal field tensor.We have shown that this approach is efficient and accu-rate over the dynamical range probed by the Millenniumsimulation. When the reconstruction is carried out on 5 h − Mpc, each component of the velocity field can be recov-ered to an accuracy of about 10 km s − . On 10 h − Mpc thisfigure is reduced to about 7 km s − . If we consider high den-sity regions, the typical uncertainty is of 13 km s − , whichimproves dramatic over linear perturbation theory; typicalerrors are 81% smaller. An accurate description of the veloc-ity divergence, both in terms of its PDF and on a point-by-point basis, is also achieved. In addition, we have shown thatthe 1- and 2-point statistics of the scaled velocity divergenceare extremely well recovered, being almost indistinguishablefrom the true ones. Contrarily, linear theory dramaticallyover-estimates the velocity divergence. This especially onthe mildly nonlinear scale of 5 h − Mpc where our methodshows more clearly its advantages. Finally, we highlight thatour method does not require calibration nor free parametersto predict the velocity divergence field.There are a number of simplifications and assumptions c (cid:13) , 000–000 osmic velocity fields -1500-1000-500050010001500-1500 -1000 -500 0 500 1000 1500 r S =5 h − MpcLIN v r ec x [ s − k m ] -1000 -500 0 500 1000 -1000-50005001000 r S =10 h − MpcLIN -1500-1000-500050010001500 v r ec x [ s − k m ] r S =5 h − MpcLOG-2LPT -1000-50005001000 r S =10 h − MpcLOG-2LPT -1500 -1000 -500 0 500 1000 1500-1500-1000-500050010001500 v r ec x [ s − k m ] r S =5 h − Mpc2LPT v Nbody x [s − km] -1000 -500 0 500 1000 -1000-50005001000 r S =10 h − Mpc2LPT v Nbody x [s − km] Figure 7.
Cell-to-cell comparison between the true velocity v Nbody x and the reconstructed one v rec x for the x -component. Left panels onscales of 5 h − Mpc and right panels on scales of 10 h − Mpc. Upper panels: LIN, middle panels: LOG-2LPT, lower panels: 2LPT (see § (cid:13)000
Cell-to-cell comparison between the true velocity v Nbody x and the reconstructed one v rec x for the x -component. Left panels onscales of 5 h − Mpc and right panels on scales of 10 h − Mpc. Upper panels: LIN, middle panels: LOG-2LPT, lower panels: 2LPT (see § (cid:13)000 , 000–000 Kitaura et al. -8 -7 -6 -5 -4 -3 -2 -1 -1500 -1000 -500 0 500 1000 1500 r S =5 h − Mpc P D F LINLOG-2LPT2LPT -8 -7 -6 -5 -4 -3 -2 -1 -400 -200 0 200 400 r S =10 h − Mpc -6 -5 -4 -3 -2 -1 P D F δ > r S =5 h − Mpc -5 -4 -3 -2 -1 δ > r S =10 h − Mpc -1500 -1000 -500 0 500 1000 150010 -8 -6 -4 -2 P D F (cid:15) v,x [s − km] δ < . r S =5 h − Mpc -400 -200 0 200 40010 -8 -6 -4 -2 δ < . (cid:15) v,x [s − km] r S =10 h − Mpc
Figure 8.
Statistics of the errors in the reconstructed velocity for the x -component (green: LIN; black: LOG-LPT, red:2LPT) withdifferent smoothing scales r S = 5 (left panels) and 10 (right panels) h − Mpc. Upper panels in the entire δ range, middle panels: for therange δ >
2, lower panels: for the range δ < .
5. c (cid:13) , 000–000 osmic velocity fields that we have adopted throughout our paper. First, our anal-yses focused on the peculiar velocity averaged over a volume.Another aspect is that we have neglected the rotational com-ponent of the velocity field. This however, does not seem tobe relevant at the scales we have investigated (larger than5 h − Mpc). Another simplification is performing our com-parison assuming that there is an unbiased estimation of theunderlying real-space density field. But, of course, densitiesmeasured in a survey are in redshift-space. The transfor-mation can be done, but not trivially. One alternative is toapply the Gibbs-sampling method suggested by Kitaura &Enßlin (2008); Kitaura et al. (2010) to correct for redshift-space distortions. In this, the Gaussian distribution of er-rors in our method is highly convenient, since it permitsto model the uncertainties in the peculiar velocity field in-cluding observational errors in a straightforward way. Oneshould consider also classical iterative approaches pioneeredby Yahil et al. (1991) based on linear theory. In a similarway they have been applied with different approximationsby different groups (see Croft & Gaztanaga 1997; Nusser &Branchini 2000; Monaco et al. 2000; Branchini et al. 2002;Mohayaee & Tully 2005; Lavaux et al. 2008; Wang et al.2012). Here it is crucial to have an accurate descriptionrelating the peculiar velocity field to the density field (orgalaxy distribution) which subject of study is central in thework presented here. We expect that the improved relationfound in this work with respect to linear Eulerian or linearLagrangian perturbation theory yields better estimates ofthe peculiar velocity field. Further studies need to be doneto test the performance including redshift-space distortions.We would like to note that our comparison and the un-certainties quoted here, were based on the present-day out-put of the Millennium Run. Such simulation was carriedout with a value for σ about 10% higher than the currentbest estimations (see Angulo & White 2010, for a methodto correct for this). Therefore, our uncertainties should beregarded as an upper limit of the reachable uncertainties fora hypothetical spectroscopic survey, which should contain aless nonlinear underlying dark matter distribution.We finalise this paper by emphasising that the methodpresented here can potentially be used in many differentapplications, and should be further developed and testedto perform bias studies combining galaxy redshift surveyswith measured peculiar velocities, Baryon Acoustic Oscilla-tion reconstructions, determination of the growth factor, toestimates of the initial conditions of the Universe. ACKNOWLEDGEMENTS
FSK and REA thank Simon White for encouraging conver-sations. The authors appreciate Francisco Prada’s commentson the manuscript. We are indebted to the German Astro-physical Virtual Observatory (GAVO) and the MPA facili-ties for providing us the Millennium Run simulation data.The work of REA was supported by Advanced Grant 246797”GALFORMOD” from the European Research Council. YHhas been partially supported by the ISF (13/08).
REFERENCES
Angulo R. E., Baugh C. M., Frenk C. S., Lacey C. G., 2008,MNRAS, 383, 755Angulo R. E., White S. D. M., 2010, MNRAS, 405, 143Bernardeau F., 1992, ApJ, 390, L61Bernardeau F., Chodorowski M. J., (cid:32)Lokas E. L., Stompor R.,Kudlicki A., 1999, MNRAS, 309, 543Bernardeau F., Colombi S., Gazta˜naga E., Scoccimarro R., 2002,Phys.Rep., 367, 1Bilicki M., Chodorowski M. J., 2008, MNRAS, 391, 1796Bouchet F. R., Colombi S., Hivon E., Juszkiewicz R., 1995,Astr.Astrophy., 296, 575Branchini E., Eldar A., Nusser A., 2002, MNRAS, 335, 53Branchini E., Freudling W., Da Costa L. N., Frenk C. S., Gio-vanelli R., Haynes M. P., Salzer J. J., Wegner G., Zehavi I.,2001, MNRAS, 326, 1191Buchert T., 1994, MNRAS, 267, 811Buchert T., Ehlers J., 1993, MNRAS, 264, 375Buchert T., Melott A. L., Weiss A. G., 1994, Astr.Astrophy.,288, 349Catelan P., 1995, MNRAS, 276, 115Chodorowski M. J., (cid:32)Lokas E. L., Pollo A., Nusser A., 1998, MN-RAS, 300, 1027Coles P., Jones B., 1991, MNRAS, 248, 1Courtois H. M., Hoffman Y., Tully R. B., Gottlober S., 2011,ArXiv e-printsCroft R. A. C., Gaztanaga E., 1997, MNRAS, 285, 793Davis M., Nusser A., Willick J. A., 1996, ApJ, 473, 22Eisenstein D. J., Seo H.-J., Sirko E., Spergel D. N., 2007, ApJ,664, 675Falck B. L., Neyrinck M. C., Aragon-Calvo M. A., Lavaux G.,Szalay A. S., 2011, ArXiv e-printsFisher K. B., Lahav O., Hoffman Y., Lynden-Bell D., ZaroubiS., 1995, MNRAS, 272, 885Frisch U., Matarrese S., Mohayaee R., Sobolevski A., 2002, Na-ture, 417, 260Gottloeber S., Hoffman Y., Yepes G., 2010, ArXiv e-printsGramann M., 1993a, ApJ, 405, 449Gramann M., 1993b, ApJ, 405, L47Guzzo L., Pierleoni M., Meneux B., Branchini E., Le F`evre O.,Marinoni C., Garilli B., et al 2008, Nature, 451, 541Hivon E., Bouchet F. R., Colombi S., Juszkiewicz R., 1995,Astr.Astrophy., 298, 643Jasche J., Kitaura F. S., 2009, ArXiv e-printsJenkins A., 2010, MNRAS, 403, 1859Kashlinsky A., Atrio-Barandela F., Ebeling H., 2011, ApJ, 732,1Kashlinsky A., Atrio-Barandela F., Kocevski D., Ebeling H.,2009, ApJ, 691, 1479Kitaura F., Jasche J., Metcalf R. B., 2010, MNRAS, 403, 589Kitaura F. S., Angulo R. E., 2011, ArXiv e-printsKitaura F. S., Enßlin T. A., 2008, MNRAS, 389, 497Kitaura F.-S., Gallerani S., Ferrara A., 2010, ArXiv e-printsKlypin A., Hoffman Y., Kravtsov A. V., Gottl¨ober S., 2003, ApJ,596, 19Kudlicki A., Chodorowski M., Plewa T., R´o˙zyczka M., 2000,MNRAS, 316, 464Lavaux G., Mohayaee R., Colombi S., Tully R. B., BernardeauF., Silk J., 2008, MNRAS, 383, 1292Melott A. L., Buchert T., Weib A. G., 1995, Astr.Astrophy., 294,345Mohayaee R., Tully R. B., 2005, ApJ, 635, L113Monaco P., Efstathiou G., 1999, MNRAS, 308, 763Monaco P., Efstathiou G., Maddox S. J., Branchini E., FrenkC. S., McMahon R. G., Oliver S. J., Rowan-Robinson M., et al2000, MNRAS, 318, 681c (cid:13)000
Angulo R. E., Baugh C. M., Frenk C. S., Lacey C. G., 2008,MNRAS, 383, 755Angulo R. E., White S. D. M., 2010, MNRAS, 405, 143Bernardeau F., 1992, ApJ, 390, L61Bernardeau F., Chodorowski M. J., (cid:32)Lokas E. L., Stompor R.,Kudlicki A., 1999, MNRAS, 309, 543Bernardeau F., Colombi S., Gazta˜naga E., Scoccimarro R., 2002,Phys.Rep., 367, 1Bilicki M., Chodorowski M. J., 2008, MNRAS, 391, 1796Bouchet F. R., Colombi S., Hivon E., Juszkiewicz R., 1995,Astr.Astrophy., 296, 575Branchini E., Eldar A., Nusser A., 2002, MNRAS, 335, 53Branchini E., Freudling W., Da Costa L. N., Frenk C. S., Gio-vanelli R., Haynes M. P., Salzer J. J., Wegner G., Zehavi I.,2001, MNRAS, 326, 1191Buchert T., 1994, MNRAS, 267, 811Buchert T., Ehlers J., 1993, MNRAS, 264, 375Buchert T., Melott A. L., Weiss A. G., 1994, Astr.Astrophy.,288, 349Catelan P., 1995, MNRAS, 276, 115Chodorowski M. J., (cid:32)Lokas E. L., Pollo A., Nusser A., 1998, MN-RAS, 300, 1027Coles P., Jones B., 1991, MNRAS, 248, 1Courtois H. M., Hoffman Y., Tully R. B., Gottlober S., 2011,ArXiv e-printsCroft R. A. C., Gaztanaga E., 1997, MNRAS, 285, 793Davis M., Nusser A., Willick J. A., 1996, ApJ, 473, 22Eisenstein D. J., Seo H.-J., Sirko E., Spergel D. N., 2007, ApJ,664, 675Falck B. L., Neyrinck M. C., Aragon-Calvo M. A., Lavaux G.,Szalay A. S., 2011, ArXiv e-printsFisher K. B., Lahav O., Hoffman Y., Lynden-Bell D., ZaroubiS., 1995, MNRAS, 272, 885Frisch U., Matarrese S., Mohayaee R., Sobolevski A., 2002, Na-ture, 417, 260Gottloeber S., Hoffman Y., Yepes G., 2010, ArXiv e-printsGramann M., 1993a, ApJ, 405, 449Gramann M., 1993b, ApJ, 405, L47Guzzo L., Pierleoni M., Meneux B., Branchini E., Le F`evre O.,Marinoni C., Garilli B., et al 2008, Nature, 451, 541Hivon E., Bouchet F. R., Colombi S., Juszkiewicz R., 1995,Astr.Astrophy., 298, 643Jasche J., Kitaura F. S., 2009, ArXiv e-printsJenkins A., 2010, MNRAS, 403, 1859Kashlinsky A., Atrio-Barandela F., Ebeling H., 2011, ApJ, 732,1Kashlinsky A., Atrio-Barandela F., Kocevski D., Ebeling H.,2009, ApJ, 691, 1479Kitaura F., Jasche J., Metcalf R. B., 2010, MNRAS, 403, 589Kitaura F. S., Angulo R. E., 2011, ArXiv e-printsKitaura F. S., Enßlin T. A., 2008, MNRAS, 389, 497Kitaura F.-S., Gallerani S., Ferrara A., 2010, ArXiv e-printsKlypin A., Hoffman Y., Kravtsov A. V., Gottl¨ober S., 2003, ApJ,596, 19Kudlicki A., Chodorowski M., Plewa T., R´o˙zyczka M., 2000,MNRAS, 316, 464Lavaux G., Mohayaee R., Colombi S., Tully R. B., BernardeauF., Silk J., 2008, MNRAS, 383, 1292Melott A. L., Buchert T., Weib A. G., 1995, Astr.Astrophy., 294,345Mohayaee R., Tully R. B., 2005, ApJ, 635, L113Monaco P., Efstathiou G., 1999, MNRAS, 308, 763Monaco P., Efstathiou G., Maddox S. J., Branchini E., FrenkC. S., McMahon R. G., Oliver S. J., Rowan-Robinson M., et al2000, MNRAS, 318, 681c (cid:13)000 , 000–000 Kitaura et al.
Moutarde F., Alimi J.-M., Bouchet F. R., Pellat R., Ramani A.,1991, ApJ, 382, 377Narayanan V. K., Weinberg D. H., 1998, ApJ, 508, 440Neyrinck M. C., Szapudi I., Szalay A. S., 2009, ApJ, 698, L90Nusser A., Branchini E., 2000, MNRAS, 313, 587Nusser A., Branchini E., Davis M., 2011, ArXiv e-printsNusser A., Dekel A., 1992, ApJ, 391, 443Nusser A., Dekel A., Bertschinger E., Blumenthal G. R., 1991,ApJ, 379, 6Peebles P. J. E., 1989, ApJ, 344, L53Scoccimarro R., Sheth R. K., 2002, MNRAS, 329, 629Springel V., White S. D. M., Jenkins A., Frenk C. S., Yoshida N.,Gao L., Navarro J., Thacker R., Croton D., Helly J., PeacockJ. A., Cole S., Thomas P., Couchman H., Evrard A., ColbergJ., Pearce F., 2005, Nature, 435, 629Wang H., Mo H. J., Yang X., van den Bosch F. C., 2012, MN-RAS, 420, 1809Watkins R., Feldman H. A., Hudson M. J., 2009, MNRAS, 392,743Willick J. A., Strauss M. A., 1998, ApJ, 507, 64Willick J. A., Strauss M. A., Dekel A., Kolatt T., 1997, ApJ,486, 629Yahil A., Strauss M. A., Davis M., Huchra J. P., 1991, ApJ, 372,380Zaroubi S., Hoffman Y., Dekel A., 1999, ApJ, 520, 413Zel’dovich Y. B., 1970, Astr.Astrophy., 5, 84
APPENDIX A: THIRD ORDER LAGRANGIANPERTURBATION THEORY (3LPT)
For the sake of completeness we investigate here third orderLPT. Following Buchert (1994) and Catelan (1995) one canwrite the displacement field to third order as Ψ = − D ∇ φ (1) + D ∇ φ (2) + D ∇ φ (3)a + D ∇ φ (3)b + D ∇ × A (3) (A1)where { D , D , D } are the 3rd order growth factorscorresponding to the gradient of two scalar potentials( φ (3)a , φ (3)b ) and the curl of a vector potential ( A (3) ). Particu-lar expressions for the irrotational 3rd order growth factors( D , D ) can be found in Bouchet et al. (1995), the growthfactor of the rotational term ( D ) is given in Catelan (1995).We assume that the fields are curl-free on scales > ∼ h − Mpc (see Bouchet et al. 1995, and the comparison be-tween the velocity divergence and the curl in the MillenniumRun in Fig. 5). We therefore consider only the scalar poten-tial terms φ (3)a and φ (3)b .The first term is cubic in the linear potential δ (3)a ≡ µ (3) ( φ (1) ) = det (cid:16) φ (1) ,ij (cid:17) , (A2)and the second term is the interaction term between thefirst- and the second-order potentials: δ (3)b ≡ µ (2) ( φ (1) , φ (2) ) = 12 (cid:88) i (cid:54) = j (cid:16) φ (2) ,ii φ (1) ,jj − φ (2) ,ij φ (1) ,ji (cid:17) , (A3)(see Buchert 1994; Bouchet et al. 1995; Catelan 1995). Inorder to ensure that the term δ (3)b is curl-free one has tointroduce some proper weights in the expression A3 (seeBuchert 1994; Catelan 1995). From the displacement field one can derive the expres-sion for the velocity field v = − D f H ∇ φ (1) + D f H ∇ φ (2) (A4)+ D f H ∇ φ (3)a + D f H ∇ φ (3)b , with f i = d ln D i / d ln a (particular expressions for f and f can be found in Bouchet et al. 1995). As we can see fromEq. A4 one can construct all components from the linearpotential φ (1) .We consider here two models for the linear potential.The first model relies on the local lognormal estimate (see § θ = D δ (1) − f f D µ (2) ( φ (1) ) − f f D µ (3) ( φ (1) ) − f f D µ (2) ( φ (1) , φ (2) ) . (A5)The second model yields a non-local estimate of φ (1) to3rd order given by δ L = δ − δ NL (A6)with δ NL = − D µ (2) ( φ (1) ) − D µ (3) ( φ (1) ) − D µ (2) ( φ (1) , φ (2) ) + µ (2) (Θ) + µ (3) (Θ). Note that thepotential Θ is also different and is determined by Eq. A1recalling the relation to the displacement field Ψ = −∇ Θ.Using the latter expression one can write the θ - δ rela-tion to 3rd order in LPT as θ = δ − (cid:18) f f − (cid:19) D µ (2) ( φ (1) ) − (cid:18) f f − (cid:19) D µ (3) ( φ (1) ) − (cid:18) f f − (cid:19) D µ (2) ( φ (1) , φ (2) ) − µ (2) (Θ) − µ (3) (Θ) . (A7)Note that Eq. A6 should be solved iteratively as we doin the 2LPT case (see § φ (1) into Eq. A5 yielding a non-local estimate ofthe linear field in one iteration. To minimize the deviationbetween LPT and the full nonlinear evolution we have addi-tionally smoothed the density δ in Eq. A7 with a Gaussiankernel of 2 h − Mpc radius.We do not find an improvement in the determination ofthe velocity divergence with respect to 2LPT including bothcurl-free terms using any of the estimates of the linear com-ponent. This could be due to an inaccurate determinationof the interaction term φ (3)b , since uncertainties in the linearcomponent φ (1) propagate more dramatically than in termswhich depend only on the initial conditions ( φ (2) , φ (3)a ).One may neglect the interaction term δ (3)b and consideronly the cubic term δ (3)a to 3rd order. Such a truncated 3LPTmodel includes the main body of the perturbation sequencewith the rest of the sequence being made up of interactionterms (Buchert 1994).Our calculations show in this case better results than2LPT for low values of the velocity divergence (see Fig. A1). c (cid:13) , 000–000 osmic velocity fields r S =5 h − MpcLOG-3LPT θ r ec L O G − L P T r S =10 h − MpcLOG-3LPT θ r ec L P T r S =5 h − Mpc3LPT θ Nbody r S =10 h − Mpc3LPT θ Nbody
Figure A1.
Cell-to-cell correlation between the true and the reconstructed scaled velocity divergence using 3LPT. Left panels on scalesof 5 h − Mpc and right panels on scales of 10 h − Mpc. Upper panels: LOG-3LPT (based on the truncated 3LPT expansion computingEq. A5 with the local model for the linear component), lower panels: 3LPT (based on the truncated 3LPT expansion computing Eq. A7with the non-local model for the linear component).
However, for large values of ∇· v we obtain larger dispersionswhich could be also due to numerical errors in the estimateof the linear density component δ L . The errors in the velocityestimation are only moderately reduced with respect to the2LPT case (see § c (cid:13)000