The Cosmological Trajectories Method: Modelling cosmic structure formation in the non-linear regime
MMNRAS , 1–15 (2021) Preprint 4 February 2021 Compiled using MNRAS L A TEX style file v3.0
The Cosmological Trajectories Method: Modelling cosmic structureformation in the non-linear regime
F. C. Lane, ★ A. N. Taylor, † and D. Sorini ‡ Scottish Universities Physics Alliance, Institute for Astronomy, School of Physics and Astronomy, University of Edinburgh, Royal Observatory,Blackford Hill, Edinburgh, EH9 3HJ, U.K.
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We introduce a novel approach, the Cosmological Trajectories Method (CTM), to model nonlinear structure formation inthe Universe by expanding gravitationally-induced particle trajectories around the Zel’dovich approximation. A new BeyondZel’dovich approximation is presented, which expands the CTM to leading second-order in the gravitational interaction andallows for post-Born gravitational scattering. In the Beyond Zel’dovich approximation we derive the exact expression for thematter clustering power spectrum. This is calculated to leading order and is available in the CTM Module. We comparethe Beyond Zel’dovich approximation power spectrum and correlation function to other methods including 1-loop StandardPerturbation Theory (SPT), 1-loop Lagrangian Perturbation Theory (LPT), and Convolution Lagrangian Perturbation Theory(CLPT). We find that the Beyond Zel’dovich approximation power spectrum performs well, matching simulations to within ± 𝑧 = 𝑧 = ± Key words:
Cosmology – methods: data analysis,statistical – cosmological parameters – large-scale structure of Universe
Deciphering how the cosmic web and large-scale structure is formedin our Universe is an essential part of understanding cosmology.Better knowledge of large-scale structure formation will allow usto extract more information from current (e.g. Planck Collaborationet al. 2020, Hildebrandt et al. 2017, Abbott et al. 2018) and futureobservations of our Universe. Gathering more statistical informa-tion from current and upcoming surveys such as the Dark EnergySpectroscopic Instrument (Levi et al. 2019), the Vera Rubin Obser-vatory (LSST Science Collaboration et al. 2009) and
Euclid (Raccaet al. 2016) will lead to tighter constraints on viable cosmological,gravity and structure formation models.Modelling the cosmic web involves knowing how structures formunder the influence of gravity. In the first approximation, the equa-tions governing the evolution of density perturbations can be lin-earised. While this approach is accurate enough to describe thelarge-scale modes, it inevitably breaks down on small scales, wherethe local density field can become much larger than the averagebackground density of the Universe. The breakdown of linear the-ory was found to occur around Fourier modes with wavenumber 𝑘 > . − (Sugiyama 2014, McQuinn & White 2016), hencethis regime is generally referred to as “non-linear regime”. A furtherdegree of complexity comes into play when considering the impactof baryonic effects on galactic scales, such as winds ejected due to ★ E-mail: fl[email protected] † E-mail: [email protected] ‡ E-mail: [email protected] supernovae explosions or jets from active galactic nuclei (see e.g.Somerville & Davé 2015, for a review).Because of the non-linear and interconnected nature of the physicalprocesses driving structure formation, the current preferred methodfor investigating structure formation in the non-linear regime is torun large cosmological simulations. N-body (dark matter only) andhydrodynamic (dark matter and baryons) simulations can be usedto simulate the gravitational evolution of structure in the Universe.One of the first large N-body simulations, the Millennium simu-lation (Springel et al. 2005), modelled the evolution of around amillion dark matter particles from 𝑧 =
127 to 𝑧 =
0. Recent hy-drodynamic simulations such as EAGLE (Crain et al. 2015, Schayeet al. 2015), Illustrius-TNG (Weinberger et al. 2017, Pillepich et al.2018), the New Horizon runs (Kim et al. 2011, Dubois et al. 2020)and Simba (Davé et al. 2019) have furthered our understanding ofstructure formation and baryonic effects.The large volume and high precision of data from forthcomingsurveys demands at least comparable accuracy in theoretical modelsof structure formation. For this reason, simulations would need toprobe a wide range of scales, while retaining high enough resolutionto properly capture small-scale physics. However, the consequentcomputational cost in terms of memory and computer time hindersthe exploration of a wide parameter space. This represents an issuewhen testing multiple theories of gravity and cosmological models,which typically requires obtaining predictions for several choicesof the underlying parameters. Thus, there is clearly an interest forsearching alternative and less costly methods.Cosmological emulators provide a way of predicting the non-linear growth of structure for a range of cosmological parameters © a r X i v : . [ a s t r o - ph . C O ] F e b F. C. Lane et al. and some modified gravity theories. Emulators are generally trainedon large sets of high-resolution simulation runs but once they havebeen trained on the simulation output they can be made publiclyavailable for the community to utilise. In this paper, we will utilisethe Euclid Emulator (Knabenhans et al. 2019) which was devel-oped for the
Euclid space telescope and was trained on a sample of100 input runs of PKDGRAV3 (Stadel et al. 2002, Potter et al. 2017).Other examples of emulators include CosmicEmu (Heitmann et al.2010, 2009, Lawrence et al. 2010, Heitmann et al. 2013, Lawrenceet al. 2017) trained on the Coyote Universe simulations, mgemu (Ra-machandra et al. 2020) an emulator that can model the ratio betweenthe Λ -CDM power spectrum and the Hu-Sawiki 𝑓 ( 𝑅 ) gravity (Hu& Sawicki 2007) power spectrum, the AEMULUS (DeRose et al.2019) project and the Dark Quest (Nishimichi et al. 2019) project.Although both simulations and emulators allow us to probe the non-linear regime accurately, analytic techniques can allow us to see howdensity field correlations arise more easily.One technique that does not require the running of simulationsand can give a more in-depth insight into structure formation is per-turbation theory. Standard Perturbation Theory (SPT) or EulerianPerturbation Theory (EPT) can reproduce cosmological observa-tions up to scales of around 𝑘 ∼ . − , where non-lineargravitational and baryonic effects become important (Peebles 1980,Bertschinger 1995, Bouchet 1996, Bernardeau et al. 2002, Carlsonet al. 2009, Bernardeau 2013). SPT can be extended to model smallerscales using IR-resummation (accounting for the physical effects ofbulk flows) and loop corrections (next to leading order corrections),as demonstrated in Crocce & Scoccimarro (2006a,b), Taruya & Hi-ramatsu (2008), Bernardeau et al. (2008, 2012, 2014) and Blas et al.(2014). Lagrangian Perturbation Theory (LPT) is another techniquethat can be used to model the formation of structure (Moutarde et al.1991, Catelan 1995, Buchert 1992, Buchert & Ehlers 1993, Bouchet1996, Tatekawa 2004, Rampf & Buchert 2012). LPT, unlike SPT,follows the motion of particles through a system and has been foundto be more accurate at equal orders than SPT (Matsubara 2008a,b,Bouchet et al. 1995, Carlson et al. 2013, White 2014, Catelan 1995).The Zel’dovich approximation (Zel’dovich 1970, Taylor 1993,Schneider & Bartelmann 1995, Taylor & Hamilton 1996, White2014), a first-order LPT, is unique in that in 1D it is exact up untilshell-crossing (the point at which streams of matter from differentdirections intersect) occurs. In 3D it behaves competitively with EPTand higher-order LPT. It is an intuitive method for describing howparticles form the structures we see in the cosmic web (McQuinn &White 2016).As discussed in McQuinn & White (2016), although techniquesthat aim to address the breakdown of perturbation theory on smallscales have made an improvement (matching simulations up to 𝑘 = . − as discussed in Sugiyama 2014), fundamental failingson these scales remain. For example, it is well known that when theoverdensity field becomes large ( 𝛿 (cid:29)
1) these schemes are no longervalid. As mentioned above, both EPT and LPT also breakdown aftershell-crossing occurs. However, we know that virialised structuresin our Universe, such as dark matter haloes, are formed after shell-crossing occurs.The Effective Field Theory of Large Scale Structure (EFTofLSS)is a method that aims to fix these issues by integrating out small wave-lengths (Carrasco et al. 2012, 2014a,b, Carroll et al. 2014, Porto et al.2014, Senatore & Zaldarriaga 2015, Vlah et al. 2015a,b, 2016a,b),to reduce the uncontrolled small-scale perturbative effects affectinglarge scales. EFTofLSS generally requires either data from simu-lations or observations to fix free parameters in the theory. Othermethods for extending perturbation theory into the non-linear regime including semi-classical propagators and methods based on field the-ory have been suggested in Seljak & Vlah (2015), Taruya & Colombi(2017), McDonald & Vlah (2018), Friedrich & Prokopec (2017),Friedrich & Prokopec (2018), Uhlemann et al. (2019), Friedrich &Prokopec (2019) and Halle et al. (2020).A statistical mechanics approach to modelling gravitational in-teraction into the non-linear regime was introduced in Bartelmannet al. (2014a) and further developed in subsequent works (Bartel-mann et al. 2014b, Fabis et al. 2014, Kozlikin et al. 2014, Viermannet al. 2015, Bartelmann et al. 2017, Sorini 2017, Lilow et al. 2019,Bartelmann et al. 2019). This theory is called Kinetic Field Theory(KFT). The theory was re-derived using particle trajectories in Ali-Haïmoud (2015). We will focus on the trajectories implementationof the technique. The initial results for the matter power spectrum(Bartelmann et al. 2014a) hinted that the method could match cur-rent simulations. Another advantage of this method is that it hasthe potential to be easily adapted to multiple cosmological models,therefore allowing predictions to be made without running numeroussimulations.In this paper, we will introduce the Cosmological TrajectoriesMethod (CTM), which expands the trajectory around the Zel’dovichapproximation. We present exact results for the matter power spec-trum to leading second-order in the displacement field, in the BeyondZel’dovich approximation, and show how an expanded version of thepower spectrum can be calculated numerically in Section 2. Finally,in Section 3.3 we will compare the Beyond Zel’dovich approxima-tion to other approximations including SPT 1-loop and LPT 1-loop.We find that the Beyond Zel’dovich approximation power spectrummatches the Euclid Emulator (Knabenhans et al. 2019) more ac-curately than the Zel’dovich approximation above 𝑧 =
1. We also findthat the Beyond Zel’dovich approximation captures the BAO peak inthe two-point correlation function more accurately than SPT 1-loop,LPT 1-loop and CLPT at 𝑧 = The fundamental idea behind KFT (Bartelmann et al. 2014a,b) isthat an ensemble of dark matter particles moves dictated by someinitial conditions until some redshift, 𝑧 ∗ , when a gravitational in-teraction term is “switched on” as in an N-body simulation. Thisgravitational interaction term is then expanded perturbatively. Thistranslates to an initial particle trajectory set by a Zel’dovich prop-agator, thus capturing the decaying velocity, with the addition of agravitational correction term, the size of which is controlled by an ex-pansion parameter, 𝜖 . The formalism in KFT is based on field theoryand therefore involves functional integrals, which is one motivationbehind the work presented in Ali-Haïmoud (2015).In Ali-Haïmoud (2015), the KFT results were re-derived in termsof particle trajectories and the Zel’dovich approximation. The trajec-tory is found from the solution to the particle equations of motion, x ( q , 𝑧 ) = x ( q , 𝑧 ∗ ) + ∫ 𝑧 ∗ 𝑧 𝑑𝑧 (cid:48) 𝑎 (cid:48) 𝐻 ( 𝑧 (cid:48) ) u ( q , 𝑧 ∗ )− 𝜖 ∫ 𝑧 ∗ 𝑧 𝑑𝑧 (cid:48) 𝑎 (cid:48) 𝐻 ( 𝑧 (cid:48) ) ∫ 𝑧 ∗ 𝑧 (cid:48) 𝑑𝑧 (cid:48)(cid:48) 𝐻 ( 𝑧 (cid:48)(cid:48) ) ∇ 𝜙 (cid:0) x (cid:0) q , 𝑧 (cid:48)(cid:48) (cid:1) , 𝑧 (cid:48)(cid:48) (cid:1) , (1)where 𝑧 ∗ is the redshift when the gravitational terms are “switchedon”, 𝑧 is redshift, 𝑎 is the scale factor, 𝐻 is the Hubble parameter, 𝜙 isthe gravitational potential, x and q are the Eulerian and Lagrangianpositions and u = 𝑎 v where v is the proper velocity. In Ali-Haïmoud MNRAS , 1–15 (2021) TM (2015) the initial position x ( q , 𝑧 ∗ ) and velocity u ( q , 𝑧 ∗ ) at time 𝑧 ∗ are set by the Zel’dovich approximation. The gravitational interactionterm, the final term in Equation (1), is found via the Poisson equationand the overdensity field, calculated in the Zel’dovich approximation.The final result derived in Ali-Haïmoud (2015) is the power spectrumto first-order in this interaction term but is exact for the Zel’dovichdensity field.The result and method presented in Ali-Haïmoud (2015) is verypromising. However in both that approach, and in KFT, the correctlinear growth is not recovered on large scales from the expansionof Equation (1). This is because to the lowest-order gravitationalinteraction takes too long to overcome the damping effect of theexpansion, leading to an underestimate of the displacement of parti-cles and growth of structure. In both cases a re-normalisation of thepredicted linear matter-density power spectrum is required.This problem motivates us to introduce the Cosmological Trajec-tories Method (CTM), which expands the gravitational interactionaround the free-field Zel’dovich approximation. This guarantees thatto first order in the displacement we match linear theory, and tosecond-order we include the effects of gravitational scattering. Thisexpansion avoids the non-local gravitational terms that appear atsecond-order in LPT (Matsubara 2008a,b, Buchert & Ehlers 1993).One issue is that the free-field Zel’dovich approximation is already anapproximation to gravitational collapse, so to avoid double-countingterms we remove a linear term from the gravitational interaction. Afull derivation of the CTM is given in Appendix A and the CTMtrajectory is presented in Equation (A2).In this paper, we will focus on the implementation where thegravitational interaction term in the CTM trajectory is expanded tosecond order in the displacement field. We shall refer to this as theBeyond Zel’dovich approximation, where the trajectory is given by(see Appendix A for a full derivation) 𝑥 𝑖 ( q , 𝑧 ) = 𝑞 𝑖 + 𝐴 ( 𝑧 ) Ψ 𝑖 ( q , 𝑧 𝑖 )+ 𝐵 𝜖 ( 𝑧 ) Ψ 𝑗 ( q , 𝑧 𝑖 ) (cid:18) 𝐸 𝑖 𝑗 ( q , 𝑧 𝑖 ) + 𝛿 ( ) ( q , 𝑧 𝑖 ) 𝛿 𝑖 𝑗 (cid:19) . (2)Here 𝑞 𝑖 is the initial Lagrangian position of the particles, Ψ 𝑖 is thelinear displacement field, 𝛿 ( ) is the linear overdensity field, 𝐸 𝑖 𝑗 = (cid:18) ∇ 𝑖 ∇ 𝑗 ∇ − − 𝛿 𝑖 𝑗 (cid:19) 𝛿 ( ) , (3)is a dimensionless, trace-free linear tidal tensor, and 𝑧 𝑖 is some ini-tial redshift. Equation (2) is the leading lowest-order gravitationalcorrection to the Zel’dovich approximation, describing post-Borngravitational deflections from the unperturbed trajectory.As we are expanding around the Zel’dovich approximation, thelinear time-dependence function in Equation (2) is 𝐴 ( 𝑧 ) = 𝐷 ( 𝑧 ) 𝐷 ( 𝑧 𝑖 ) .We note there is freedom to choose other time-dependencies, butthis ensures the lowest-order theory matches linear growth on largescales. The second time-dependent function 𝐵 𝜖 ( 𝑧 ) in Equation (2) isderived in Appendix A to match the gravitational field, and is givenby 𝐵 𝜖 ( 𝑧 ) = − 𝜖 CTM 𝜔 ∫ 𝑧𝑧 𝑖 𝑑𝑧 (cid:48) 𝑎 (cid:48) 𝐻 ( 𝑧 (cid:48) ) ∫ 𝑧 (cid:48) 𝑧 𝑖 𝑑𝑧 (cid:48)(cid:48) 𝐻 ( 𝑧 (cid:48)(cid:48) ) (cid:18) 𝐷 ( 𝑧 (cid:48)(cid:48) ) 𝐷 ( 𝑧 𝑖 ) (cid:19) . (4)where 𝜖 CTM controls the size of the higher-order gravitational term(the tidal tensor) and 𝜔 = 𝐻 Ω 𝑚 . Note that in principle one coulduse the linear growth factor for a scale-independent modified gravitytheory (Clifton et al. 2012, Nojiri et al. 2017) instead. Comparing the Beyond Zel’dovich approximation time dependence to that obtainedin Bartelmann et al. (2014a,b) and Ali-Haïmoud (2015), we see thatas we do not include the initial, decaying velocity term. Instead theparticle follows a Zel’dovich trajectory and the displacement fieldis therefore proportional to the linear growth factor. In Appendix Dmore detail on the application of second-order CTM to KFT is given.There are two free parameters in the second-order CTM trajectory;the initial redshift 𝑧 𝑖 and the expansion parameter 𝜖 CTM . The expan-sion parameter, 𝜖 CTM , controls the size of the gravitational terms. Ifone considers 𝜖 CTM as a perturbative parameter then by definition itshould be small ( 𝜖 CTM (cid:28) 𝜖 CTM val-ues will increase the impact the tidal field has on non-linear structureformation.
In this section, we give details on the calculation of the matter-density power spectrum for the Beyond Zel’dovich trajectory. Wefind that, assuming Gaussian initial conditions, we can calculate anexact expression for the matter power spectrum in this approximation.In order to explore the numerical implementation of this result, weexpand around the exact solution. Our numerical results are availableusing the CTM Module. We begin with the statistical properties ofthe linear fields.
The linear displacement field, Ψ 𝑖 , the tidal field, 𝐸 𝑖 𝑗 , and the linearoverdensity field, 𝛿 ( ) , are correlated Gaussian fields at the initialredshift, 𝑧 𝑖 . As we shall show, we can calculate the matter powerspectrum in the Beyond Zel’dovich approximation using the statis-tics of Gaussian fields (Bardeen et al. 1986, van de Weygaert &Bertschinger 1996, Taylor & Watts 2000). As the fields are Gaus-sian, they are fully specified by their covariance matrix, C , whichcontains the correlation of the fields with each other at two differentLagrangian points, q and q , at the initial redshift.We define a vector of the fields at each position, X = (cid:16) Ψ 𝑖 ( q ) , Ψ 𝑖 ( q ) , 𝐸 𝑖 ( q ) , 𝐸 𝑖 ( q ) , 𝛿 ( ) ( q ) , 𝛿 ( ) ( q ) (cid:17) , (5)where 𝐸 𝑖 = vec ( E ) is the 6-dimensional vectorisation of the distinctterms in the symmetric tidal tensor 𝐸 𝑖 𝑗 . The covariance matrix ofthe vector, C = (cid:104) XX 𝑇 (cid:105) , is given by C = 𝐶 Ψ 𝑖 Ψ 𝑗 𝐶 Ψ 𝑖 Ψ 𝑗 𝐶 𝐸 𝑖 Ψ 𝑗 𝐶 𝐸 𝑖 Ψ 𝑗 𝐶 𝛿 ( ) Ψ 𝑖 𝐶 𝛿 ( ) Ψ 𝑖 𝐶 Ψ 𝑖 Ψ 𝑗 𝐶 Ψ 𝑖 Ψ 𝑗 𝐶 𝐸 𝑖 Ψ 𝑗 𝐶 𝐸 𝑖 Ψ 𝑗 𝐶 𝛿 ( ) Ψ 𝑖 𝐶 𝛿 ( ) Ψ 𝑖 𝐶 Ψ 𝑖 𝐸 𝑗 𝐶 Ψ 𝑖 𝐸 𝑗 𝐶 𝐸 𝑖 𝐸 𝑗 𝐶 𝐸 𝑖 𝐸 𝑗 𝐶 𝛿 ( ) 𝐸 𝑖 𝐶 𝛿 ( ) 𝐸 𝑖 𝐶 Ψ 𝑖 𝐸 𝑗 𝐶 Ψ 𝑖 𝐸 𝑗 𝐶 𝐸 𝑖 𝐸 𝑗 𝐶 𝐸 𝑖 𝐸 𝑗 𝐶 𝛿 ( ) 𝐸 𝑖 𝐶 𝛿 ( ) 𝐸 𝑖 𝐶 Ψ 𝑖 𝛿 ( ) 𝐶 Ψ 𝑖 𝛿 ( ) 𝐶 𝐸 𝑖 𝛿 ( ) 𝐶 𝐸 𝑖 𝛿 ( ) 𝐶 𝛿 ( ) 𝛿 ( ) 𝐶 𝛿 ( ) 𝛿 ( ) 𝐶 Ψ 𝑖 𝛿 ( ) 𝐶 Ψ 𝑖 𝛿 ( ) 𝐶 𝐸 𝑖 𝛿 ( ) 𝐶 𝐸 𝑖 𝛿 ( ) 𝐶 𝛿 ( ) 𝛿 ( ) 𝐶 𝛿 ( ) 𝛿 ( ) , (6)where the numerical ‘1’ and ‘2’ indicate the position. The correlatorsof the linear density and displacement fields are given by MNRAS000
The linear displacement field, Ψ 𝑖 , the tidal field, 𝐸 𝑖 𝑗 , and the linearoverdensity field, 𝛿 ( ) , are correlated Gaussian fields at the initialredshift, 𝑧 𝑖 . As we shall show, we can calculate the matter powerspectrum in the Beyond Zel’dovich approximation using the statis-tics of Gaussian fields (Bardeen et al. 1986, van de Weygaert &Bertschinger 1996, Taylor & Watts 2000). As the fields are Gaus-sian, they are fully specified by their covariance matrix, C , whichcontains the correlation of the fields with each other at two differentLagrangian points, q and q , at the initial redshift.We define a vector of the fields at each position, X = (cid:16) Ψ 𝑖 ( q ) , Ψ 𝑖 ( q ) , 𝐸 𝑖 ( q ) , 𝐸 𝑖 ( q ) , 𝛿 ( ) ( q ) , 𝛿 ( ) ( q ) (cid:17) , (5)where 𝐸 𝑖 = vec ( E ) is the 6-dimensional vectorisation of the distinctterms in the symmetric tidal tensor 𝐸 𝑖 𝑗 . The covariance matrix ofthe vector, C = (cid:104) XX 𝑇 (cid:105) , is given by C = 𝐶 Ψ 𝑖 Ψ 𝑗 𝐶 Ψ 𝑖 Ψ 𝑗 𝐶 𝐸 𝑖 Ψ 𝑗 𝐶 𝐸 𝑖 Ψ 𝑗 𝐶 𝛿 ( ) Ψ 𝑖 𝐶 𝛿 ( ) Ψ 𝑖 𝐶 Ψ 𝑖 Ψ 𝑗 𝐶 Ψ 𝑖 Ψ 𝑗 𝐶 𝐸 𝑖 Ψ 𝑗 𝐶 𝐸 𝑖 Ψ 𝑗 𝐶 𝛿 ( ) Ψ 𝑖 𝐶 𝛿 ( ) Ψ 𝑖 𝐶 Ψ 𝑖 𝐸 𝑗 𝐶 Ψ 𝑖 𝐸 𝑗 𝐶 𝐸 𝑖 𝐸 𝑗 𝐶 𝐸 𝑖 𝐸 𝑗 𝐶 𝛿 ( ) 𝐸 𝑖 𝐶 𝛿 ( ) 𝐸 𝑖 𝐶 Ψ 𝑖 𝐸 𝑗 𝐶 Ψ 𝑖 𝐸 𝑗 𝐶 𝐸 𝑖 𝐸 𝑗 𝐶 𝐸 𝑖 𝐸 𝑗 𝐶 𝛿 ( ) 𝐸 𝑖 𝐶 𝛿 ( ) 𝐸 𝑖 𝐶 Ψ 𝑖 𝛿 ( ) 𝐶 Ψ 𝑖 𝛿 ( ) 𝐶 𝐸 𝑖 𝛿 ( ) 𝐶 𝐸 𝑖 𝛿 ( ) 𝐶 𝛿 ( ) 𝛿 ( ) 𝐶 𝛿 ( ) 𝛿 ( ) 𝐶 Ψ 𝑖 𝛿 ( ) 𝐶 Ψ 𝑖 𝛿 ( ) 𝐶 𝐸 𝑖 𝛿 ( ) 𝐶 𝐸 𝑖 𝛿 ( ) 𝐶 𝛿 ( ) 𝛿 ( ) 𝐶 𝛿 ( ) 𝛿 ( ) , (6)where the numerical ‘1’ and ‘2’ indicate the position. The correlatorsof the linear density and displacement fields are given by MNRAS000 , 1–15 (2021)
F. C. Lane et al. 𝐶 Ψ 𝑖 Ψ 𝑗 = (cid:104) Ψ 𝑖 ( q ) Ψ 𝑗 ( q )(cid:105) = 𝜎 𝑖 𝑗 ( 𝑞 ) , (7a) 𝐶 𝛿 ( ) Ψ 𝑖 = (cid:104) 𝛿 ( ) ( q ) Ψ 𝑖 ( q )(cid:105) = Π 𝑖 ( 𝑞 ) , (7b) 𝐶 𝛿 ( ) 𝛿 ( ) = (cid:104) 𝛿 ( ) ( q ) 𝛿 ( ) ( q )(cid:105) = 𝜉 ( 𝑞 ) , (7c)where 𝑞 = | q − q | is the distance between points, while the corre-lations of the vectorised tidal field are 𝐶 𝐸 𝑖 Ψ 𝑗 = (cid:104) 𝐸 𝑖 ( q ) Ψ 𝑗 ( q )(cid:105) , (8a) 𝐶 𝐸 𝑖 𝐸 𝑗 = (cid:104) 𝐸 𝑖 ( q ) 𝐸 𝑗 ( q )(cid:105) , (8b) 𝐶 𝐸 𝑖 𝛿 ( ) = (cid:104) 𝐸 𝑖 ( q ) 𝛿 ( ) ( q )(cid:105) . (8c)These can be written in terms of the correlations of the tensor tidalfield, (cid:104) 𝐸 𝑖 𝑗 ( q ) Ψ 𝑘 ( q )(cid:105) = Φ 𝑖 𝑗𝑘 ( 𝑞 ) , (9a) (cid:104) 𝐸 𝑖 𝑗 ( q ) 𝐸 𝑘𝑙 ( q )(cid:105) = 𝜂 𝑖 𝑗𝑘𝑙 ( 𝑞 ) , (9b) (cid:104) 𝐸 𝑖 𝑗 ( q ) 𝛿 ( ) ( q )(cid:105) = Σ 𝑖 𝑗 ( 𝑞 ) . (9c)The correlation functions given in Equations (7) and Equations (9)are defined in Appendix B. The matter-density power spectrum, P ( 𝑘 ) , is defined by the correlatorof the Fourier modes of the density field, (cid:104) 𝛿 ( k ) 𝛿 ( k )(cid:105) = ( 𝜋 ) P ( 𝑘 ) 𝛿 D ( k + k ) , (10)where the expectation value is calculated by an ensemble average.The Fourier transform of the overdensity field is given by ( 𝜋 ) 𝛿 D ( k ) + 𝛿 ( k ) = ∫ 𝑑 𝑞 e 𝑖 k · x ( q ,𝑧 ) , (11)where x is the trajectory defined in Equation (2). Therefore, the powerspectrum for the second-order CTM trajectory is given byP ( 𝑘, 𝑧 ) = ∫ 𝑑 𝑞 e 𝑖 k · q × (cid:20)(cid:28) e 𝑖𝑘 𝑖 Ψ 𝑗 ( q ,𝑧 𝑖 ) (cid:16) 𝐴 ( 𝑧 ) 𝛿 𝑖𝑗 + 𝐵 𝜖 ( 𝑧 ) 𝐸 𝑖𝑗 ( q ,𝑧 𝑖 )+ 𝐵 𝜖 ( 𝑧 ) 𝛿 ( ) ( q ,𝑧 𝑖 ) 𝛿 𝑖𝑗 (cid:17) e − 𝑖𝑘 𝑖 Ψ 𝑗 ( q ,𝑧 𝑖 ) (cid:16) 𝐴 ( 𝑧 ) 𝛿 𝑖𝑗 + 𝐵 𝜖 ( 𝑧 ) 𝐸 𝑖𝑗 ( q ,𝑧 𝑖 )+ 𝐵 𝜖 ( 𝑧 ) 𝛿 ( ) ( q ,𝑧 𝑖 ) 𝛿 𝑖𝑗 (cid:17) (cid:29) − (cid:21) . (12)We can simplify this by introducing a new vector, K = 𝐴 ( 𝑧 ) ( 𝑘 𝑖 , − 𝑘 𝑖 , , , , ) , (13)with the same dimensionality as X . If we define a new matrix M ,with the same dimensionality as C , M = 𝑖 𝐵 𝜖 𝑘 𝑖 − 𝛿 𝑗𝑘 − 𝛿 𝑗𝑘 − 𝛿 𝑗𝑘 𝛿 𝑗𝑘 − , (14) the ensemble average in Equation (12) can be rewritten in the multi-variate Gaussian form (cid:104) e 𝑖𝑘 𝑖 Ψ 𝑗 ( q ) (cid:16) 𝐴𝛿 𝑖𝑗 + 𝐵 𝜖 𝐸 𝑖𝑗 ( q )+ 𝐵 𝜖 𝛿 ( ) ( q ) 𝛿 𝑖𝑗 (cid:17) × e − 𝑖𝑘 𝑖 Ψ 𝑗 ( q ) (cid:16) 𝐴𝛿 𝑖𝑗 + 𝐵 𝜖 𝐸 𝑖𝑗 ( q )+ 𝐵 𝜖 𝛿 ( ) ( q ) 𝛿 𝑖𝑗 (cid:17) (cid:105) = ( 𝜋 ) ∫ 𝑑 𝑋 | det C | − / e − X 𝑇 C − X e 𝑖 K · X e − X 𝑇 MX . (15)Equation (15) can be integrated, resulting in an exact expression forthe matter power spectrum for the second-order CTM trajectory,P ( 𝑘, 𝑧 ) = ∫ 𝑑 𝑞 e 𝑖 k · q (cid:104) | det ( + MC )| − / e − K 𝑇 C [ + MC ] − K − (cid:105) . (16)This expression is the main result of the paper. While Equation (16) is exact, and the matrix manipulation can inprinciple be carried out numerically, the integration is highly oscil-latory and can be numerically unstable. To explore the features ofthe Beyond Zel’dovich approximation we shall expand the solutionin such a way as the take advantage of existing algorithms to treat theintegration, and to compare to other methods.The argument of the exponential in Equation (16) can be expanded; K 𝑇 C [ + MC ] − K ≈ K 𝑇 CK − K 𝑇 CMCK . (17)The first term here is K 𝑇 CK = 𝐴 ( 𝑧 ) (cid:104) 𝑘 𝜎 𝜓 ( 𝑧 𝑖 ) − 𝑘 𝑖 𝑘 𝑗 𝜎 𝑖 𝑗 ( 𝑞, 𝑧 𝑖 ) (cid:105) , (18)where 𝜎 𝑖 𝑗 ( , 𝑧 𝑖 ) = 𝜎 𝜓 𝛿 𝑖 𝑗 , while the second term vanishes. We canexpand the determinant in Equation (16) asdet ( + MC ) = exp ( tr ln ( + MC )) ≈ exp (cid:18) −
12 tr
MCMC (cid:19) (19)where tr ( MC ) = ( MCMC ) = 𝐵 𝜖 𝑘 𝑖 𝑘 𝑗 (cid:2) Π 𝑖 ( 𝑞 ) Π 𝑗 ( 𝑞 ) + 𝜎 𝑖𝑛 ( 𝑞 ) Σ 𝑛 𝑗 ( 𝑞 )+ 𝜉 ( 𝑞 ) 𝜎 𝑖 𝑗 ( 𝑞 ) − 𝜉 ( ) 𝜎 𝜓 𝛿 𝑖 𝑗 (cid:105) . (20)In this approximation the power spectrum isP ( 𝑘, 𝑧 ) = ∫ 𝑑 𝑞 e 𝑖 k · q (cid:104) e − [ K 𝑇 CK − tr ( MCMC ) ] − (cid:105) . (21)To lowest order this reduces to the Zel’dovich power spectrum (Taylor1993, Schneider & Bartelmann 1995, Taylor & Hamilton 1996).Both of the terms in the exponential in equation (21) have a factor 𝑘 𝑖 𝑘 𝑗 , so the function is Gaussian. The covariance of this Gaussianis the differential displacement covariance. Hence, we can interpretthe extra term as the lowest-order correction to the displacementcovariance matrix due to gravitational scattering. MNRAS , 1–15 (2021) TM It is useful to define the correlation function ¯ Σ 𝑖 𝑗 ( 𝑞 ) = (cid:104) ¯ 𝐸 𝑖 𝑗 ( q ) 𝛿 ( ) ( q )(cid:105) , which can be related to the un-barred corre-lation function,¯ Σ 𝑖 𝑗 ( 𝑞 ) = Σ 𝑖 𝑗 ( 𝑞 ) + 𝜉 ( 𝑞 ) 𝛿 𝑖 𝑗 . (22)The correlations 𝜎 𝑖 𝑗 , Π 𝑖 and ¯ Σ 𝑖 𝑗 can be split into irreducible compo-nents (Vlah et al. 2015a, Catelan et al. 2000, Crittenden et al. 2001).The method used to split 𝜎 𝑖 𝑗 , Π 𝑖 and ¯ Σ 𝑖 𝑗 is shown in Appendix B1.The correlations Π 𝑖 and ¯ Σ 𝑖 𝑗 can be expanded as¯ Σ 𝑖 𝑗 ( 𝑞 ) = 𝐷 ( 𝑞 ) 𝛿 𝑖 𝑗 + 𝐹 ( 𝑞 ) ˆ 𝑞 𝑖 ˆ 𝑞 𝑗 , (23) Π 𝑖 ( 𝑞 ) = 𝐺 ( 𝑞 ) ˆ 𝑞 𝑖 . (24)and 𝐷, 𝐹 and 𝐺 are defined as 𝐷 ( 𝑞 ) = 𝜋 ∫ ∞ 𝑑𝑘 [ 𝑗 ( 𝑘𝑞 ) + 𝑗 ( 𝑘𝑞 )] 𝑘 P L ( 𝑘 ) , (25a) 𝐹 ( 𝑞 ) = − 𝜋 ∫ ∞ 𝑑𝑘 𝑗 ( 𝑘𝑞 ) 𝑘 P L ( 𝑘 ) , (25b) 𝐺 ( 𝑞 ) = − 𝜋 ∫ ∞ 𝑑𝑘 𝑗 ( 𝑘𝑞 ) 𝑘 P L ( 𝑘 ) . (25c)Finally, the correlation of the displacement field can be decomposedas, 𝜎 𝑖 𝑗 = (cid:18) 𝜎 𝜓 − 𝑋 (cid:48) ( 𝑞 ) (cid:19) 𝛿 𝑖 𝑗 − 𝑌 (cid:48) ( 𝑞 ) ˆ 𝑞 𝑖 ˆ 𝑞 𝑗 , (26)with 𝑋 (cid:48) ( 𝑞 ) = 𝜋 ∫ ∞ 𝑑𝑘 (cid:20) − 𝑗 ( 𝑘𝑞 ) 𝑘𝑞 (cid:21) P L ( 𝑘, 𝑧 𝑖 ) 𝑌 (cid:48) ( 𝑞 ) = 𝜋 ∫ ∞ 𝑑𝑘 (cid:20) 𝑗 ( 𝑘𝑞 ) 𝑘𝑞 − 𝑗 ( 𝑘𝑞 ) (cid:21) P L ( 𝑘, 𝑧 𝑖 ) . (27)Substituting the decomposed correlations into Equation (20) thensplitting the integral into 𝑘 and 𝑘 𝜇 parts using the method fornumerically calculating the Zel’dovich power spectrum described inSchneider & Bartelmann (1995), Carlson et al. (2013), Sugiyama(2014) and Vlah et al. (2015a). We can use the kth moment of theintegral to calculate the angular integral 𝐼 𝑘 = ∫ − 𝑑𝜇 𝜇 𝑘 e 𝑖𝑎𝜇 e 𝑏 𝜇 (28)which can be solved using the general prescription (Schneider &Bartelmann 1995, Vlah et al. 2015a), 𝐼 𝑘 = (− 𝑖 ) 𝑘 e 𝑏 ∞ ∑︁ 𝑛 = (− 𝑏 ) 𝑛 (cid:18) 𝑑𝑑𝑘 (cid:19) 𝑘 𝑎 − 𝑛 𝑗 𝑛 ( 𝑎 ) , (29)where 𝑗 𝑛 is a spherical Bessel function. The angular parts of Equa-tion (16) are calculated using an identity resulting in the power spec-trum becoming,P ( 𝑘, 𝑧 ) ≈ 𝜋 ∫ ∞ 𝑑𝑞 𝑞 ∫ − 𝑑𝜇 e 𝑖𝑘𝑞𝜇 (cid:20) e − 𝑘 𝐴 ( 𝑋 (cid:48) + 𝜇 𝑌 (cid:48) ) e 𝐵 𝜖 𝑘 ( 𝑊 (cid:48) + 𝜇 𝑍 (cid:48) ) − e − 𝑘 𝜎 𝜓 (cid:16) 𝐴 + 𝐵 𝜖 𝜂 𝐸 (cid:17) (cid:21) , (30) where the second exponential term is a Dirac delta function at the ori-gin and has been added to cancel oscillations as described in Schnei-der & Bartelmann (1995). The Beyond Zel’dovich power spectrumto second-order is finally given byP ( 𝑘, 𝑧 ) = 𝜋 ∫ ∞ 𝑑𝑞 𝑞 e − 𝑘 [ 𝐴 ( 𝑋 (cid:48) + 𝑌 (cid:48) )− 𝐵 𝜖 ( 𝑊 (cid:48) + 𝑍 (cid:48) ) ]× ∞ ∑︁ 𝑛 = 𝑘 (cid:16) 𝐴 𝑌 (cid:48) − 𝐵 𝜖 𝑍 (cid:48) (cid:17) 𝑞 𝑛 𝑗 𝑛 ( 𝑘𝑞 ) (31)where 𝑊 (cid:48) = − 𝜎 𝜓 𝜂 𝐸 + (cid:18) 𝜎 𝜓 − 𝑋 (cid:48) (cid:19) (cid:18) 𝐷 − 𝜉 (cid:19) , (32a) 𝑍 (cid:48) = 𝐺 + (cid:18) 𝜎 𝜓 − 𝑋 (cid:48) (cid:19) 𝐹 − 𝑌 (cid:48) (cid:18) 𝐷 + 𝐹 − 𝜉 (cid:19) . (32b)In the above expressions for 𝑊 (cid:48) and 𝑍 (cid:48) all functions apart from 𝜎 𝜓 and 𝜂 𝐸 are evaluated at 𝑞 . The power spectra presented in the remainder of this paper have beencalculated using the CTM Module . The initial power spectra andcosmological parameters are calculated using classylss and thespherical Bessel integrals are calculated using mcfit . The powerspectra are calculated using Planck18 (Planck Collaboration et al.2020) cosmology ( Ω 𝑚 = . , ℎ = . , 𝑛 𝑠 = . 𝜎 = . 𝑛 =
32 spherical-Bessel functionswhen calculating the power spectra. See Appendix C for more detailson the numerical integration tools used in the CTM Module.There are two free parameters in the second-order CTM trajectory:the initial redshift, 𝑧 𝑖 , and the expansion parameter, 𝜖 CTM . The choiceof the initial redshift does not make a noticeable difference to the finalpower spectrum unless a very low value such as 𝑧 𝑖 =
10 is chosen.Since we assume that the fields are initially Gaussian a sufficientlyhigh value of 𝑧 𝑖 must be chosen to not invalidate the method. In thispaper, we will set 𝑧 𝑖 = 𝜖 CTM values, 𝜖 CTM = . , . , 𝑧 = 𝑧 = 𝑧 = 𝑧 = (Kn-abenhans et al. 2019) are shown in dashed grey lines. In Knabenhanset al. (2019), the emulator is found to be ±
1% accurate compared tosimulations at 𝑧 = ±
1% up to 𝑘 = − at 𝑧 =
1. Above 𝑧 = ≈
3% accurate. The emulator was built on a sampleof 100 input runs of PKDGRAV3 (Stadel et al. 2002, Potter et al.2017).The power spectra shown in Figure 1 have been truncated at 𝑘 = . − . The second-order CTM trajectory is only appli-cable until this 𝑘 -value as the method suffers from numerical issues https://github.com/franlane94/CTM https://classylss.readthedocs.io/en/stable/ https://github.com/eelregit/mcfit https://github.com/miknab/EuclidEmulator/wiki/III)-UsageMNRAS000
3% accurate. The emulator was built on a sampleof 100 input runs of PKDGRAV3 (Stadel et al. 2002, Potter et al.2017).The power spectra shown in Figure 1 have been truncated at 𝑘 = . − . The second-order CTM trajectory is only appli-cable until this 𝑘 -value as the method suffers from numerical issues https://github.com/franlane94/CTM https://classylss.readthedocs.io/en/stable/ https://github.com/eelregit/mcfit https://github.com/miknab/EuclidEmulator/wiki/III)-UsageMNRAS000 , 1–15 (2021) F. C. Lane et al. − − − ∆ ( k , z ) z = 0LinearBeyond Zel dovich , (cid:15) CTM = 1Beyond Zel dovich , (cid:15) CTM = 0 . dovich , (cid:15) CTM = 0 . z = 110 − − k [h Mpc − ]10 − − − − ∆ ( k , z ) z = 2 10 − − k [h Mpc − ] z = 3 Figure 1.
The dimensionless power spectrum for linear theory (black solid line), Beyond Zel’dovich with 𝜖 CTM = 𝜖 CTM = . 𝜖 CTM = .
01 (dotted blue line) at 𝑧 = 𝑧 = 𝑧 = 𝑧 = beyond this point and it is difficult to disentangle these from physicaleffects. This is addressed in Appendix C. At all redshifts, the BeyondZel’dovich approximation calculated with 𝜖 CTM = . 𝜖 CTM = . 𝜖 CTM parameter controls the size of the gravitational correctionto the Zel’dovich trajectory. The Beyond Zel’dovich approximationwith 𝜖 CTM = 𝑧 = 𝜖 CTM parameter on the Beyond Zel’dovich powerspectrum is shown in more detail in Figure 2. The maximum 𝑘 -valuereached before the difference, Δ diff = P calc ( 𝑘 ) − P emu ( 𝑘 ) P emu ( 𝑘 ) , (33) where P emu is the power spectrum obtained using the Euclid Emu-lator exceeds Δ diff = ± .
05 is shown.One can see more clearly that the Beyond Zel’dovich approxi-mation calculated with smaller values of 𝜖 CTM (shown in blue plussigns and orange crosses) converges to the Zel’dovich approxima-tion (shown in black circles). This validates our approximationas if 𝜖 CTM = 𝜖 CTM = To reduce the impact of the small-scale breakdown on larger scales,we will introduce a Gaussian damped initial power spectrum defined
MNRAS , 1–15 (2021) TM − − k m a x [ h M p c − ] Beyond Zel dovich , (cid:15) CTM = 1Beyond Zel dovich , (cid:15) CTM = 0 . dovich , (cid:15) CTM = 0 . dovich Figure 2.
The maximum 𝑘 -value reached, 𝑘 max , before the difference betweenthe Beyond Zel’dovich approximation and the Euclid Emulator exceeds ±
5% versus redshift. The purple diamonds represent the Beyond Zel’dovichapproximation calculated with 𝜖 CTM =
1, the orange crosses are calculatedwith 𝜖 CTM = . 𝜖 CTM = .
01. The black circlesrepresent the Zel’dovich approximation. asP damped ( 𝑘, 𝑧 ) = e − (cid:16) 𝑘𝑘𝑐 (cid:17) P L ( 𝑘, 𝑧 ) (34)where 𝑘 𝑐 is the cut-off scale.In Figure 3, the maximum 𝑘 − value reached before the differencebetween the calculated Beyond Zel’dovich power spectrum and theemulator power spectrum becomes larger than Δ diff = ±
5% is shownversus redshift. The Beyond Zel’dovich power spectra were calcu-lated using 𝜖 CTM = 𝑘 𝑐 =
50 h Mpc − (blue plus signs), 𝑘 𝑐 = − (purple di-amonds) and 𝑘 𝑐 = . − (orange crosses). The highest cut-offvalue of 𝑘 𝑐 =
50 h Mpc − has no noticeable effect and the lowest cut-off value of 𝑘 𝑐 = . − does not counteract the over dampingon small scales until 𝑧 =
0. The largest and smallest cut-off valuesare too stringent and either restrict structure formation too much ortoo little in the desired regime. The value of 𝑘 𝑐 = − , how-ever, appears to effectively remove the influence of the breakdownon larger scales at a wide range of redshifts.Therefore, in Figure 4 a range of cut-off values centered around 𝑘 𝑐 = − are tested. As was the case previously small cut-off values ( 𝑘 𝑐 < − ) have a detrimental effect on structureformation at high redshifts. We will choose to set 𝑘 𝑐 = − toremove the effect of the breakdown on the Beyond Zel’dovich powerspectrum for as wide a redshift range as possible. Hence, in all futureplots the Beyond Zel’dovich power spectra are calculated with aninitial Gaussian damped power spectrum with 𝑘 𝑐 = − .Although not shown in this paper, we tested the dependence ofthe cut-off parameter 𝑘 𝑐 on the cosmology chosen. We comparedthe value of 𝑘 max reached when the damped power spectrum wascalculated using a value of Ω 𝑚 = . Ω 𝑚 = . ∼
5% on average,implying that the cut-off scale is likely only weakly dependent onboth cosmology and redshift. As we were comparing our results tothe Euclid Emulator we were limited in the range of Ω 𝑚 values we − k m a x [ h M p c − ] Undamped k c = 0 . − k c = 5 hMpc − k c = 50 hMpc − Figure 3.
The maximum 𝑘 − value reached before the difference between theBeyond Zel’dovich approximation power spectrum, calculated with 𝜖 CTM = ±
5% is shown vs. redshift. k c [h Mpc − ]10 − k m a x [ h M p c − ] z = 0 z = 1 z = 2 z = 3 z = 4 z = 5 Figure 4.
The maximum 𝑘 − value reached before the difference between theBeyond Zel’dovich approximation power spectrum, calculated with 𝜖 CTM = ±
5% is shown vs. the cut-off value 𝑘 𝑐 . could choose. We leave it to future work to obtain simulation datafor a wider range of cosmological parameters to more stringently testthe dependence of the cut-off on cosmology. In this section, we will investigate the performance of the Be-yond Zel’dovich approximation for modelling the two-point cor-relation function. Specifically, we are interested in modelling themildly non-linear regime, as this is where we find the BAO signal
MNRAS000
MNRAS000 , 1–15 (2021)
F. C. Lane et al. ( 𝑟 ≈
100 Mpc h − or 𝑘 ≈ .
01 h Mpc − ) first detected in Eisen-stein et al. (2005) and Cole et al. (2005). In Figure 5, the solid linesshow the scaled correlation functions for the Beyond Zel’dovich ap-proximation (upper-left panel), SPT 1-loop (upper-right panel) asdiscussed in Crocce & Scoccimarro (2006a,b), Taruya & Hiramatsu(2008), Bernardeau et al. (2008, 2012, 2014), Blas et al. (2014), LPT1-loop (lower-left panel) as described in Matsubara (2008a,b), Carl-son et al. (2009), Vlah et al. (2015a), Sugiyama (2014), Carlson et al.(2013), McQuinn & White (2016), Vlah et al. (2015a) and CLPT(lower-right panels) as described in Carlson et al. (2013), Wang et al.(2013), Vlah et al. (2015a,b, 2016b). The different colours representredshifts 𝑧 = 𝑧 = 𝑧 = 𝑧 =
0. The BeyondZel’dovich approximation matches the Euclid Emulator moreclosely in the mildly non-linear regime than CLPT. The accuracyof the Beyond Zel’dovich approximation in the BAO peak regime isdue to the inclusion of the tidal field term in Equation (2), as thespatial deformation responsible for the non-linear evolution of theBAO peak, is encoded within the tidal tensor.
In this section, we will compare the Beyond Zel’dovich approxima-tion (with 𝜖 CTM = 𝑘 𝑐 = − ) power spectrum calculatedusing the CTM Module to other methods. The first method we com-pared the Beyond Zel’dovich approximation to is the Zel’dovich ap-proximation which has also been computed with an initial Gaussiandamped power spectrum. In Figure 7, the top-left panel shows the dif-ference between the Beyond Zel’dovich (solid lines) and Zel’dovich(dashed) power spectra at different redshifts. Above 𝑧 = 𝑧 = 𝑧 = 𝑧 =
5, the Beyond Zel’dovich approxi-mation performs as well LPT 1-loop until around 𝑘 = . − .Finally, in the bottom-right panel the Beyond Zel’dovich approxi-mation is compared to 3-point Convolution Lagrangian PerturbationTheory (CLPT) and computed using CLEFT (the CLPT correlationfunction presented in Section 3.2 was also calculated this way). AgainCLPT matches the emulator results more stringently in the non-linearregime at all redshifts compared to the Beyond Zel’dovich approxi-mation. In summary, the Beyond Zel’dovich approximation is moreaccurate than the Zel’dovich approximation above redshift 𝑧 = 𝑘 = . − . In this paper, we have introduced the Cosmological TrajectoriesMethod (CTM). The leading second-order CTM trajectory, the Be-yond Zel’dovich approximation is comprised of the Zel’dovich ap-proximation with a gravitational correction term given by the productof the linear displacement field and a tidal tensor. This post-Bornapproximation to the Zel’dovich approximation should capture non-linear effects such as gravitational deflection. We then introduceda special case of second-order CTM called the Beyond Zel’dovichapproximation in which the linear order terms were proportional tothe linear growth factor, 𝐷 . We have calculated the exact expressionfor the Beyond Zel’dovich matter power spectrum, assuming Gaus-sian initial conditions. A numerical implementation of this expandsaround this solution, for stability. The Beyond Zel’dovich approxi-mation computed with an initial Gaussian damped power spectrumoutperformed the Zel’dovich approximation (also computed with aninitial Gaussian damped power spectrum) when compared to thepower spectrum obtained using the Euclid Emulator at redshiftsabove 𝑧 =
1. The Beyond Zel’dovich approximation also matchedthe emulator correlation function between 𝑟 = − and 𝑟 =
10 Mpc h − as well as SPT 1-loop, LPT 1-loop and CLPTat 𝑧 =
0. As demonstrated in Figure 6 the Beyond Zel’dovich ap-proximation further models the BAO peak in the correlation functionmore accurately than the Zel’dovich approximation, SPT 1-loop andCLPT at 𝑧 = 𝑘 = . − when compared to the emulator power spectrum above 𝑧 =
1. Al-though LPT 1-loop, SPT 1-loop and CLPT outperformed the Be-yond Zel’dovich approximation on small scales ( 𝑘 ≥ . − ),on mildly non-linear scales the Beyond Zel’dovich approximationmatched the emulator power spectrum. This suggests that the CTMcould be implemented to produce mock observables for future BAOobservations taken by instruments such as DESI in its current state.We investigated and fixed the two free parameters of the theory, 𝜖 CTM which controls the size of the correction term, and 𝑧 𝑖 the initialredshift at which the correlations of the fields are calculated at to be 𝜖 CTM = 𝑧 𝑖 = 𝑘 ≥ − ) numericalintegration issues were encountered, resulting in the distrust of resultsbeyond this point. The numerical issues and possible solutions werediscussed in Appendix C. We introduced an initial Gaussian damped https://github.com/modichirag/CLEFTMNRAS , 1–15 (2021) TM r ξ ( r ) [ M p c h − ] Beyond Zel dovich z = 0 z = 1 z = 2 z = 3 SPT 1 − loop20 40 60 80 100 120 r [Mpc h − ]010203040 r ξ ( r ) [ M p c h − ] LPT 1 − loop 20 40 60 80 100 120 r [Mpc h − ] CLPT Figure 5.
The scaled two-point correlation function calculated using the Beyond Zel’dovich approximation (upper-left panel), SPT 1-loop (upper-right panel),LPT 1-loop (lower-left panel) and CLPT (lower-right panel) is shown for four redshifts ( 𝑧 = 𝑧 = 𝑧 = 𝑧 = power spectrum with a cut-off scale of 𝑘 𝑐 = − to reduce theexcessive damping on mildly non-linear scales due to shell-crossingshowed in Figure 1. We leave it to future work to investigate thedependence of this Gaussian cut-off on both cosmology and redshiftas it may have an impact on the application of the CTM to modifiedgravity theories or large deviations from Λ -CDM cosmology.In conclusion, both the CTM and the Beyond Zel’dovich approxi-mation appear to be valuable tools for studying nonlinear clusteringof matter and galaxies in the Universe. The Beyond Zel’dovich ap-proximation can be used to interpret future BAO observations byinstruments such as DESI and
LSST . Furthermore, both approxi-mations could be used in conjunction with Lyman- 𝛼 observations,re-ionisation studies and other high redshift surveys, as we verifiedthat our approach performs particularly well at higher redshift. Thus,the CTM promises to enable placing even tighter statistical con-straints on viable models of dark matter and dark energy, as well ason modified gravity theories. ACKNOWLEDGEMENTS
The authors would like to thank Matthias Bartelmann, Yacine Ali-Haimoud, Zvonimir Vlah and Yanchuan Cai for useful discussions.F.C. Lane acknowledges the support of the UK Science and Technol-ogy Facilities Council and the Scottish Universities Physics Alliance.A.N. Taylor thanks the Royal Society for the support of a WolfsonResearch Merit Award and a STFC Consolidated Grant. D. Soriniis supported by the European Research Council, under grant no.670193.
REFERENCES
Abbott T. M. C., et al., 2018, MNRAS, 480, 3879Ali-Haïmoud Y., 2015, Phys. Rev. D, 91, 103507Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, ApJ, 304, 15MNRAS000
Abbott T. M. C., et al., 2018, MNRAS, 480, 3879Ali-Haïmoud Y., 2015, Phys. Rev. D, 91, 103507Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, ApJ, 304, 15MNRAS000 , 1–15 (2021) F. C. Lane et al.
80 90 100 110 r [Mpc h − ] − . − . . . . . . ξ c a l c ( r ) − ξ e m u ( r ) ξ e m u ( r ) Zel dovichBeyond Zel dovichLPT 1 − loopCLPT Figure 6.
The difference, Δ diff , between the Beyond Zel’dovich (purple dia-monds), Zel’dovich approximation (black circles), LPT 1-loop (orange trian-gles) and CLPT (blue stars) correlation functions and the Euclid Emulatorat 𝑧 =
0. The grey shaded region shows Δ diff ± . arXiv:astro-ph/9603013 )Bouchet F. R., Colombi S., Hivon E., Juszkiewicz R., 1995, A&A, 296, 575Buchert T., 1992, MNRAS, 254, 729Buchert T., Ehlers J., 1993, MNRAS, 264, 375Carlson J., White M., Padmanabhan N., 2009, Phys. Rev. D, 80, 043531Carlson J., Reid B., White M., 2013, MNRAS, 429, 1674Carrasco J. J. M., Hertzberg M. P., Senatore L., 2012, Journal of High EnergyPhysics, 2012, 82Carrasco J. J. M., Foreman S., Green D., Senatore L., 2014a, Journal ofCosmology and Astro-Particle Physics, 2014, 056Carrasco J. J. M., Foreman S., Green D., Senatore L., 2014b, Journal ofCosmology and Astro-Particle Physics, 2014, 057Carroll S. M., Leichenauer S., Pollack J., 2014, Phys. Rev. D, 90, 023518Catelan P., 1995, MNRAS, 276, 115Catelan P., Porciani C., Kamionkowski M., 2000, MNRAS, 318, L39Clifton T., Ferreira P. G., Padilla A., Skordis C., 2012, Physics Reports, 513,1–189Cole S., et al., 2005, Monthly Notices of the Royal Astronomical Society,362, 505–534Crain R. A., et al., 2015, MNRAS, 450, 1937Crittenden R. G., Natarajan P., Pen U.-L., Theuns T., 2001, ApJ, 559, 552Crocce M., Scoccimarro R., 2006a, Phys. Rev. D, 73, 063519Crocce M., Scoccimarro R., 2006b, Phys. Rev. D, 73, 063520Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, Monthly Notices of the Royal Astronomical Society,486, 2827–2849DeRose J., et al., 2019, ApJ, 875, 69Dubois Y., et al., 2020, arXiv e-prints, p. arXiv:2009.10578Eisenstein D. J., et al., 2005, The Astrophysical Journal, 633, 560–574Fabis F., Berg D., Kozlikin E., Bartelmann M., 2014, arXiv e-prints, p.arXiv:1412.2572Fang X., Blazek J. A., McEwen J. E., Hirata C. M., 2017, Journal of Cosmol-ogy and Astroparticle Physics, 2017, 030–030Friedrich P., Prokopec T., 2017, Phys. Rev. D, 96, 083504Friedrich P., Prokopec T., 2018, Phys. Rev. D, 98, 025010Friedrich P., Prokopec T., 2019, Phys. Rev. D, 100, 103527Halle A., Nishimichi T., Taruya A., Colombi S., Bernardeau F., 2020, arXive-prints, p. arXiv:2001.10417Hamilton A. J. S., 2015, FFTLog: Fast Fourier or Hankel transform(ascl:1512.017)Heitmann K., Higdon D., White M., Habib S., Williams B. J., Lawrence E.,Wagner C., 2009, The Astrophysical Journal, 705, 156–174Heitmann K., White M., Wagner C., Habib S., Higdon D., 2010, The Astro-physical Journal, 715, 104–121Heitmann K., Lawrence E., Kwan J., Habib S., Higdon D., 2013, The Astro-physical Journal, 780, 111Hildebrandt H., et al., 2017, MNRAS, 465, 1454Hu W., Sawicki I., 2007, Physical Review D, 76Kim J., Park C., Rossi G., Lee S. M., Gott J. Richard I., 2011, Journal ofKorean Astronomical Society, 44, 217Knabenhans M., et al., 2019, MNRAS, 484, 5509–5529Kozlikin E., Fabis F., Lilow R., Viermann C., Bartelmann M., 2014, arXive-prints, p. arXiv:1412.2715LSST Science Collaboration et al., 2009, LSST Science Book, Version 2.0( arXiv:0912.0201 )Lawrence E., Heitmann K., White M., Higdon D., Wagner C., Habib S.,Williams B., 2010, The Astrophysical Journal, 713, 1322–1331Lawrence E., et al., 2017, The Astrophysical Journal, 847, 50Levi M. E., et al., 2019, The Dark Energy Spectroscopic Instrument (DESI)( arXiv:1907.10688 )Lilow R., Fabis F., Kozlikin E., Viermann C., Bartelmann M., 2019, Journalof Cosmology and Astro-Particle Physics, 2019, 001Matsubara T., 2008a, Phys. Rev. D, 77, 063530Matsubara T., 2008b, Phys. Rev. D, 78, 083519McDonald P., Vlah Z., 2018, Phys. Rev. D, 97, 023508McEwen J. E., Fang X., Hirata C. M., Blazek J. A., 2016, Journal of Cosmol-ogy and Astroparticle Physics, 2016, 015–015McQuinn M., White M., 2016, Journal of Cosmology and Astro-ParticlePhysics, 2016, 043Mehrem R., 2011, The Plane Wave Expansion, Infinite Integrals and Identitiesinvolving Spherical Bessel Functions ( arXiv:0909.0494 )Moutarde F., Alimi J. M., Bouchet F. R., Pellat R., Ramani A., 1991, ApJ,382, 377Nishimichi T., et al., 2019, ApJ, 884, 29Nojiri S., Odintsov S., Oikonomou V., 2017, Physics Reports, 692, 1–104Peebles P. J. E., 1980, The large-scale structure of the universePillepich A., et al., 2018, MNRAS, 473, 4077Planck Collaboration et al., 2020, A&A, 641, A6Porto R. A., Senatore L., Zaldarriaga M., 2014, Journal of Cosmology andAstro-Particle Physics, 2014, 022Potter D., Stadel J., Teyssier R., 2017, Computational Astrophysics and Cos-mology, 4, 2Racca G. D., et al., 2016, Space Telescopes and Instrumentation 2016: Optical,Infrared, and Millimeter WaveRamachandra N., Valogiannis G., Ishak M., Heitmann K., 2020, Mat-ter Power Spectrum Emulator for f(R) Modified Gravity Cosmologies( arXiv:2010.00596 )Rampf C., Buchert T., 2012, Journal of Cosmology and Astroparticle Physics,2012, 021–021Schaye J., et al., 2015, MNRAS, 446, 521Schneider P., Bartelmann M., 1995, MNRAS, 273, 475Seljak U., Vlah Z., 2015, Phys. Rev. D, 91, 123516MNRAS , 1–15 (2021) TM − . − . − . − . − . . . P c a l c ( k ) − P e m u ( k ) P e m u ( k ) Zel dovich z = 0 z = 1 z = 2 z = 3 z = 4 z = 5 SPT 1 − loop10 − − k [h Mpc − ] − . − . − . − . − . . . P c a l c ( k ) − P e m u ( k ) P e m u ( k ) LPT 1 − loop 10 − − k [h Mpc − ]CLPT Figure 7.
The difference between a given theory and the emulator results is shown for redshifts 𝑧 = , , , , Δ diff ± . arXiv:astro-ph/0412025 )Taylor A. N., 1993, in Bouchet F., Lachieze-Rey M., eds, Cosmic VelocityFields. p. 585Taylor A. N., Hamilton A. J. S., 1996, MNRAS, 282, 767Taylor A. N., Watts P. I. R., 2000, MNRAS, 314, 92 Uhlemann C., Rampf C., Gosenca M., Hahn O., 2019, Phys. Rev. D, 99,083524Viermann C., Fabis F., Kozlikin E., Lilow R., Bartelmann M., 2015, Phys.Rev. E, 91, 062120Vlah Z., Seljak U., Baldauf T., 2015a, Phys. Rev. D, 91, 023508Vlah Z., White M., Aviles A., 2015b, Journal of Cosmology and Astro-ParticlePhysics, 2015, 014Vlah Z., Seljak U., Yat Chu M., Feng Y., 2016a, Journal of Cosmology andAstro-Particle Physics, 2016, 057Vlah Z., Castorina E., White M., 2016b, Journal of Cosmology and Astro-Particle Physics, 2016, 007Wang L., Reid B., White M., 2013, MNRAS, 437, 588–599Weinberger R., et al., 2017, MNRAS, 465, 3291White M., 2014, MNRAS, 439, 3630Zel’dovich Y. B., 1970, A&A, 5, 84van de Weygaert R., Bertschinger E., 1996, MNRAS, 281, 84MNRAS000
The difference between a given theory and the emulator results is shown for redshifts 𝑧 = , , , , Δ diff ± . arXiv:astro-ph/0412025 )Taylor A. N., 1993, in Bouchet F., Lachieze-Rey M., eds, Cosmic VelocityFields. p. 585Taylor A. N., Hamilton A. J. S., 1996, MNRAS, 282, 767Taylor A. N., Watts P. I. R., 2000, MNRAS, 314, 92 Uhlemann C., Rampf C., Gosenca M., Hahn O., 2019, Phys. Rev. D, 99,083524Viermann C., Fabis F., Kozlikin E., Lilow R., Bartelmann M., 2015, Phys.Rev. E, 91, 062120Vlah Z., Seljak U., Baldauf T., 2015a, Phys. Rev. D, 91, 023508Vlah Z., White M., Aviles A., 2015b, Journal of Cosmology and Astro-ParticlePhysics, 2015, 014Vlah Z., Seljak U., Yat Chu M., Feng Y., 2016a, Journal of Cosmology andAstro-Particle Physics, 2016, 057Vlah Z., Castorina E., White M., 2016b, Journal of Cosmology and Astro-Particle Physics, 2016, 007Wang L., Reid B., White M., 2013, MNRAS, 437, 588–599Weinberger R., et al., 2017, MNRAS, 465, 3291White M., 2014, MNRAS, 439, 3630Zel’dovich Y. B., 1970, A&A, 5, 84van de Weygaert R., Bertschinger E., 1996, MNRAS, 281, 84MNRAS000 , 1–15 (2021) F. C. Lane et al.
APPENDIX A: THE CTM TRAJECTORY CALCULATION
The CTM trajectory is an expansion of the gravitationally inducedtrajectory around the Zel’dovich approximation, given by x ( q , 𝑡 ) = q + 𝐴 ( 𝑡 ) 𝚿 ( q , 𝑡 𝑖 ) + 𝜖 CTM ∫ 𝑡𝑡 𝑖 𝑑𝑡 (cid:48) 𝑎 (cid:48) ∫ 𝑡 (cid:48) 𝑡 𝑖 𝑑𝑡 (cid:48)(cid:48) Δ g (cid:0) x (cid:0) q , 𝑡 (cid:48)(cid:48) (cid:1) , 𝑡 (cid:48)(cid:48) (cid:1) (A1)where 𝑎 is the scale factor, q is the initial position, 𝚿 is the linearorder displacement field, 𝑡 𝑖 is some initial time. The gravitationalfield is given by g = −∇ x 𝜙 where 𝜙 is the gravitational potentialand 𝜖 CTM is an expansion parameter used to control the size of thehigher-order gravitational term. As the Zel’dovich approximationalready extrapolates the effects of the linear gravitational field, weadd the differential gravitational field Δ g = g − g L , where we haveremoved the linear field to avoid double-counting forces. The CTMtrajectory in Equation (A1) then describes a particle moving underfree motion given by the Zel’dovich approximation with the additionof a gravitational correction term.The CTM approach is a hybrid of KFT (Bartelmann et al. 2014a,b,Ali-Haïmoud 2015) and LPT (Moutarde et al. 1991, Catelan 1995,Buchert 1992, Buchert & Ehlers 1993, Bouchet 1996, Tatekawa 2004,Rampf & Buchert 2012). In the KFT approach, the trajectory is basedon the formal solution to the particle equations of motion, and so havethe free particle motion is damped by the expansion. Here we havechosen to have the trajectory in Equation (A1) to be defined by theZel’dovich approximation plus a higher-order gravitational term weare free to pick the time-dependent function, 𝐴 ( 𝑡 ) , to be the lineargrowth factor, 𝐷 , which will allow us to avoid re-normalisation onlarge scales. In the remainder of this section, we will demonstrate howthe second-order CTM trajectory, given in Equation (2) is obtainedby solving the gravitational field using the Zel’dovich approximationas a basis. This allows us to avoid the inclusion of non-local terms,which arise when considering second-order LPT.We can write the trajectory (A1) as x ( q , 𝑡 ) = q + x ( q , 𝑡 ) + x ( q , 𝑡 ) (A2)with x ( q , 𝑡 ) = 𝐴 ( 𝑡 ) 𝚿 ( q , 𝑡 𝑖 ) , (A3) x ( q , 𝑡 ) = 𝜖 CTM ∫ 𝑡𝑡 𝑖 𝑑𝑡 (cid:48) 𝑎 (cid:48) ∫ 𝑡 (cid:48) 𝑡 𝑖 𝑑𝑡 (cid:48)(cid:48) Δ g (cid:0) x (cid:0) q , 𝑡 (cid:48)(cid:48) (cid:1) , 𝑡 (cid:48)(cid:48) (cid:1) (A4)The overdensity field, 𝛿 ( x , 𝑡 ) , is given by 𝛿 ( x , 𝑡 ) = ∫ 𝑑 𝑞 𝛿 D ( x − q − x − x ) − is ( 𝜋 ) 𝛿 D ( k ) + 𝛿 ( k , 𝑡 ) = ∫ 𝑑 𝑞 e 𝑖 k · q e 𝑖 k ·( x + x ) . (A6)Using the Poisson equation, ∇ 𝜙 = 𝐻 Ω 𝑚 𝑎 − 𝛿 = 𝜔 𝑎 − 𝛿 , with 𝜔 = 𝐻 Ω 𝑚 the gravitational field g ( x , 𝑡 ) can be written as g ( x , 𝑡 ) = 𝑖𝜔 𝑎 − ∫ 𝑑 𝑘 ( 𝜋 ) e − 𝑖 k · x k 𝑘 𝛿 ( k , 𝑡 ) (A7) Our Fourier transform convention is 𝑓 ( k ) = ∫ 𝑑 𝑥 e 𝑖 k · x 𝑓 ( x ) and 𝑓 ( x ) = ( 𝜋 ) ∫ 𝑑 𝑘 e − 𝑖 k · x 𝑓 ( k ) . where 𝛿 ( k , 𝑡 ) is the non-linear overdensity field. The overdensityfield is given by the Zel’dovich approximation, 𝛿 ( x , 𝑡 ) = ∫ 𝑑 𝑞 𝛿 D ( x − q − x ) − , (A8)and the linear displacement is given by 𝚿 = ∫ 𝑑 𝑘 ( 𝜋 ) e − 𝑖 k · x k 𝑘 𝛿 ( ) ( k , 𝑡 ) . (A9)These equations define our Cosmological Trajectories Method.The Beyond Zel’dovich approximation calculates the gravitationalcorrection to second order in the displacement field. To this order thegravitational field is g = 𝜔 𝑎 𝐴 ( 𝑡 ) (cid:16) 𝚿 − 𝐴 ( 𝑡 ) (cid:104) ( 𝚿 · ∇ ) 𝚿 + ∇ ∇ − (cid:16) ¯ 𝐸 − | 𝛿 ( ) | (cid:17)(cid:105)(cid:17) , (A10)where ¯ 𝐸 = ¯ 𝐸 𝑖 𝑗 ¯ 𝐸 𝑗𝑖 , the tidal tensor is ¯ 𝐸 𝑖 𝑗 = ∇ 𝑖 ∇ 𝑗 ∇ − 𝛿 ( ) and 𝛿 ( ) is the linear overdensity field. The first term in equation (A10) isthe linear gravitation field, while the second term represents a localchange from Lagrangian to Eulerian coordinates. The last term isthe force generated by nonlinear, second-order growth of structure.This term is non-local, depending on all points in the density andtidal field through the inverse Laplacian. However, we can keep ouranalysis local as we can expect the second term, proportional tothe displacement field which we will extrapolate, to be larger thanthe third term. In addition, we can expect the cancellation between¯ 𝐸 and | 𝛿 ( ) | in the third terms to reduce its effects. Hence thedifferential gravitational field, to leading second-order, is Δ g ≈ − 𝜔 𝑎 𝐴 ( 𝑡 ) [ 𝚿 ( q , 𝑡 𝑖 ) · ∇ ] 𝚿 ( q , 𝑡 𝑖 ) . (A11)Therefore, we define the Beyond Zel’dovich approximation as 𝑥 𝑖 ( q , 𝑡 ) = 𝑞 𝑖 + Ψ 𝑖 ( q , 𝑡 𝑖 ) (cid:2) 𝐴 ( 𝑡 ) 𝛿 𝑖 𝑗 + 𝐵 𝜖 ( 𝑡 ) ¯ 𝐸 𝑖 𝑗 ( q , 𝑡 𝑖 ) (cid:3) , (A12)where the tidal term describes the effects of gravitational scattering.The time-dependent function 𝐵 𝜖 ( 𝑡 ) is, 𝐵 𝜖 ( 𝑡 ) = − 𝜖 CTM 𝜔 ∫ 𝑡𝑡 𝑖 𝑑𝑡 (cid:48) 𝑎 (cid:48) ∫ 𝑡 (cid:48) 𝑡 𝑖 𝑑𝑡 (cid:48)(cid:48) 𝑎 (cid:48)(cid:48) 𝐴 (cid:0) 𝑡 (cid:48)(cid:48) (cid:1) , (A13)which after using 𝑑𝑡 = − 𝑎𝐻 𝑑𝑧 can be written as, 𝐵 𝜖 ( 𝑡 ) = − 𝜖 CTM 𝜔 ∫ 𝑧𝑧 𝑖 𝑑𝑧 (cid:48) 𝑎 (cid:48) 𝐻 ( 𝑧 (cid:48) ) ∫ 𝑧 (cid:48) 𝑧 𝑖 𝑑𝑧 (cid:48)(cid:48) 𝐻 ( 𝑧 (cid:48)(cid:48) ) 𝐴 (cid:0) 𝑧 (cid:48)(cid:48) (cid:1) . (A14)The Beyond Zel’dovich approximation used in the main body of thepaper is the second-order CTM trajectory given in Equation (A12)with the following time-dependent functions 𝐴 ( 𝑧 ) = 𝐷 ( 𝑧 ) 𝐷 ( 𝑧 𝑖 ) , (A15) 𝐵 𝜖 ( 𝑧 ) = − 𝜖 CTM 𝜔 ∫ 𝑧𝑧 𝑖 𝑑𝑧 (cid:48) 𝑎 (cid:48) 𝐻 ( 𝑧 (cid:48) ) ∫ 𝑧 (cid:48) 𝑧 𝑖 𝑑𝑧 (cid:48)(cid:48) 𝐻 ( 𝑧 (cid:48)(cid:48) ) (cid:18) 𝐷 ( 𝑧 (cid:48)(cid:48) ) 𝐷 ( 𝑧 𝑖 ) (cid:19) . (A16)This time dependence reproduces linear growth on large scales with-out the need for re-normalisation. MNRAS , 1–15 (2021) TM APPENDIX B: CORRELATION FUNCTIONS
In Section 2.1.1 we defined the covariance matrix, C = (cid:104) XX 𝑇 (cid:105) ,where 𝑋 𝛼 = (cid:16) Ψ 𝑖 , Ψ 𝑖 , 𝐸 𝑖 , 𝐸 𝑖 , 𝛿 ( ) , 𝛿 ( ) (cid:17) . In order to evaluate theCTM power spectrum it is useful to have the correlation functionsdefined in Equations (7a, 8a, 8b, 8c, 7b, 7c) re-expressed in of ¯ 𝐸 𝑖 𝑗 .For example, Φ 𝑖 𝑗𝑘 ( 𝑞 ) = (cid:104) 𝐸 𝑖 𝑗 ( q ) Ψ 𝑘 ( q )(cid:105) = (cid:28)(cid:18) ¯ 𝐸 𝑖 𝑗 ( q ) − 𝛿 ( ) ( q ) (cid:19) Ψ 𝑘 ( q ) (cid:29) = ¯ Φ 𝑖 𝑗𝑘 ( q ) − Π 𝑘 ( q ) 𝛿 𝑖 𝑗 (B1)where¯ Φ 𝑖 𝑗𝑘 ( 𝑞 ) = (cid:104) ¯ 𝐸 𝑖 𝑗 ( q ) Ψ 𝑘 ( q )(cid:105) = 𝑖 ∫ 𝑑 𝑘 ( 𝜋 ) e 𝑖 k · q 𝑘 𝑖 𝑘 𝑗 𝑘 𝑘 𝑘 P L ( 𝑘, 𝑧 𝑖 ) , (B2a) Π 𝑖 ( 𝑞 ) = (cid:104) Ψ 𝑖 ( q ) 𝛿 ( ) ( q )(cid:105) = 𝑖 ∫ 𝑑 𝑘 ( 𝜋 ) e 𝑖 k · q 𝑘 𝑖 𝑘 P L ( 𝑘, 𝑧 𝑖 ) , (B2b)with P L ( 𝑘, 𝑧 𝑖 ) being the linear power spectrum evaluated at theinitial redshift. Similarly, 𝜂 𝑖 𝑗𝑘𝑙 ( 𝑞 ) = ¯ 𝜂 𝑖 𝑗𝑘𝑙 ( 𝑞 ) + 𝜉 ( 𝑞 ) 𝛿 𝑖 𝑗 𝛿 𝑘𝑙 − (cid:0) ¯ Σ 𝑖 𝑗 ( 𝑞 ) 𝛿 𝑘𝑙 + ¯ Σ 𝑘𝑙 ( 𝑞 ) 𝛿 𝑖 𝑗 (cid:1) , (B3) Σ 𝑖 𝑗 ( 𝑞 ) = ¯ Σ 𝑖 𝑗 ( 𝑞 ) − 𝜉 ( 𝑞 ) 𝛿 𝑖 𝑗 , (B4)with¯ 𝜂 𝑖 𝑗𝑘𝑙 ( 𝑞 ) = (cid:104) ¯ 𝐸 𝑖 𝑗 ( q ) ¯ 𝐸 𝑘𝑙 ( q )(cid:105) = ∫ 𝑑 𝑘 ( 𝜋 ) e 𝑖 k · q 𝑘 𝑖 𝑘 𝑗 𝑘 𝑘 𝑘 𝑙 𝑘 P L ( 𝑘, 𝑧 𝑖 ) , (B5a)¯ Σ 𝑖 𝑗 ( 𝑞 ) = (cid:104) ¯ 𝐸 𝑖 𝑗 ( q ) 𝛿 ( ) ( q )(cid:105) = ∫ 𝑑 𝑘 ( 𝜋 ) e 𝑖 k · q 𝑘 𝑖 𝑘 𝑗 𝑘 P L ( 𝑘, 𝑧 𝑖 ) , (B5b) 𝜉 ( 𝑞 ) = (cid:104) 𝛿 ( ) ( q ) 𝛿 ( ) ( q )(cid:105) = ∫ 𝑑 𝑘 ( 𝜋 ) e 𝑖 k · q P L ( 𝑘, 𝑧 𝑖 ) . (B5c) B1 An Example of the splitting of a correlation function
We will implement the identities presented in Catelan et al. (2000)and Crittenden et al. (2001) to calculate ¯ Σ 𝑖 𝑗 and Π 𝑖 . For example,¯ Σ 𝑖 𝑗𝑘 given in Equation (B5b) can be written as¯ Σ 𝑖 𝑗 ( 𝑞, 𝑧 𝑖 ) = ( 𝜋 ) ∫ − 𝑑𝜇 e 𝑖𝑘𝑞𝜇 ∫ ∞ 𝑑𝑘 𝑘 𝑖 𝑘 𝑗 P L ( 𝑘, 𝑧 𝑖 ) (B6)which after performing the angle integral and substituting in 𝑖𝑘 𝑖 = ∇ 𝑖 results in¯ Σ 𝑖 𝑗 ( 𝑞, 𝑧 𝑖 ) = − 𝜋 ∫ ∞ 𝑑𝑘 P L ( 𝑘, 𝑧 𝑖 ) ∇ 𝑖 ∇ 𝑗 𝑗 ( 𝑘𝑞 ) . (B7)Defining the following ∇ 𝑖 = 𝑞 𝑖 𝑞 𝑑𝑑𝑞 = 𝑞 𝑖 𝐷 𝑞 and 𝐷 𝑛𝑟 𝑗 ( 𝑟 ) = (− ) 𝑛 𝑟 − 𝑛 𝑗 𝑛 ( 𝑟 ) (B8) where 𝑟 = 𝑘𝑞 in our case allows us to decompose ¯ Σ 𝑖 𝑗 as¯ Σ 𝑖 𝑗 ( 𝑞, 𝑧 𝑖 ) = 𝐷 ( 𝑞, 𝑧 𝑖 ) 𝛿 𝑖 𝑗 + 𝐹 ( 𝑞, 𝑧 𝑖 ) ˆ 𝑞 𝑖 ˆ 𝑞 𝑗 (B9)with 𝐷 ( 𝑞, 𝑧 𝑖 ) and 𝐹 ( 𝑞, 𝑧 𝑖 ) as defined in Equations (25a) and (25b). APPENDIX C: COMMENTS ON NUMERICALINTEGRATION
In Section 3 it is mentioned that we only trust the CTM power spec-trum up until 𝑘 = . − . After this point, there are numericaluncertainties due to the highly oscillatory spherical Bessel integralsinvolved in the calculation. We investigated multiple techniques toremedy these numerical issues for large- 𝑘 values.We first implemented an alternative numerical integration tech-nique to the one introduced in Section 2.1.4. This alternate techniquewas introduced in Vlah et al. (2015a) and involves a generalisationof the plane wave expansion (Mehrem 2011). We will briefly sum-marise this alternative integration technique here but refer the readerto Vlah et al. (2015a) for a full derivation. The plane wave expansionise 𝑖 k · q = ∞ ∑︁ 𝑙 = 𝑖 𝑙 ( 𝑙 + ) P 𝑙 ( cos 𝜃 ) 𝑗 𝑙 ( 𝑘𝑟 ) (C1)using this we can write that ( 𝑖𝑥 ) 𝑛 = ∞ ∑︁ 𝑙 = 𝑖 𝑙 ( 𝑙 + ) P 𝑙 ( 𝑥 ) (cid:18) 𝑑 𝑛 𝑗 𝑙 ( 𝛼 ) 𝑑𝛼 𝑛 (cid:19) 𝛼 = . (C2)This is simply the Taylor expansion of the spherical Bessel functionaround zero. Comparing this Taylor expansion with another wellknown representation of spherical Bessel function , 𝑗 𝑙 ( 𝛼 ) = 𝛼 𝑙 ∞ ∑︁ 𝑘 = (− ) 𝑘 𝑘 𝑘 ! 𝛼 𝑘 ( 𝑙 + 𝑘 + ) !! , (C3)we can write Equation (C1) as ( 𝑖𝑥 ) 𝑛 = ∞ ∑︁ 𝑙 = 𝑖 𝑙 ( 𝑙 + ) P 𝑙 ( 𝑥 ) 𝑏 𝑙𝑛 (C4)where 𝑏 𝑙𝑛 = 𝑖 𝑛 − 𝑙 𝑛 ! √ 𝑛 − 𝑙 (cid:16) ( 𝑛 − 𝑙 ) (cid:17) ! , ( 𝑛 + 𝑙 + ) !! , if 𝑛 ≥ 𝑙 and 𝑛 and 𝑙 are both even or odd0 , otherwise. (C5)Therefore, ∫ − 𝑑𝜇 e 𝑖𝐴𝜇 e 𝐵𝜇 = ∞ ∑︁ 𝑛 = ( 𝑛 ) !2 𝑛 𝑛 ! 𝐵 𝑛 × 𝑛 ∑︁ 𝑝 = (− ) 𝑝 𝑝 + ( 𝑛 − 𝑝 ) ! ( 𝑛 + 𝑝 + ) !! 𝑗 𝑝 ( 𝐴 ) . (C6) https://dlmf.nist.gov MNRAS000
In Section 3 it is mentioned that we only trust the CTM power spec-trum up until 𝑘 = . − . After this point, there are numericaluncertainties due to the highly oscillatory spherical Bessel integralsinvolved in the calculation. We investigated multiple techniques toremedy these numerical issues for large- 𝑘 values.We first implemented an alternative numerical integration tech-nique to the one introduced in Section 2.1.4. This alternate techniquewas introduced in Vlah et al. (2015a) and involves a generalisationof the plane wave expansion (Mehrem 2011). We will briefly sum-marise this alternative integration technique here but refer the readerto Vlah et al. (2015a) for a full derivation. The plane wave expansionise 𝑖 k · q = ∞ ∑︁ 𝑙 = 𝑖 𝑙 ( 𝑙 + ) P 𝑙 ( cos 𝜃 ) 𝑗 𝑙 ( 𝑘𝑟 ) (C1)using this we can write that ( 𝑖𝑥 ) 𝑛 = ∞ ∑︁ 𝑙 = 𝑖 𝑙 ( 𝑙 + ) P 𝑙 ( 𝑥 ) (cid:18) 𝑑 𝑛 𝑗 𝑙 ( 𝛼 ) 𝑑𝛼 𝑛 (cid:19) 𝛼 = . (C2)This is simply the Taylor expansion of the spherical Bessel functionaround zero. Comparing this Taylor expansion with another wellknown representation of spherical Bessel function , 𝑗 𝑙 ( 𝛼 ) = 𝛼 𝑙 ∞ ∑︁ 𝑘 = (− ) 𝑘 𝑘 𝑘 ! 𝛼 𝑘 ( 𝑙 + 𝑘 + ) !! , (C3)we can write Equation (C1) as ( 𝑖𝑥 ) 𝑛 = ∞ ∑︁ 𝑙 = 𝑖 𝑙 ( 𝑙 + ) P 𝑙 ( 𝑥 ) 𝑏 𝑙𝑛 (C4)where 𝑏 𝑙𝑛 = 𝑖 𝑛 − 𝑙 𝑛 ! √ 𝑛 − 𝑙 (cid:16) ( 𝑛 − 𝑙 ) (cid:17) ! , ( 𝑛 + 𝑙 + ) !! , if 𝑛 ≥ 𝑙 and 𝑛 and 𝑙 are both even or odd0 , otherwise. (C5)Therefore, ∫ − 𝑑𝜇 e 𝑖𝐴𝜇 e 𝐵𝜇 = ∞ ∑︁ 𝑛 = ( 𝑛 ) !2 𝑛 𝑛 ! 𝐵 𝑛 × 𝑛 ∑︁ 𝑝 = (− ) 𝑝 𝑝 + ( 𝑛 − 𝑝 ) ! ( 𝑛 + 𝑝 + ) !! 𝑗 𝑝 ( 𝐴 ) . (C6) https://dlmf.nist.gov MNRAS000 , 1–15 (2021) F. C. Lane et al.
In Vlah et al. (2015a) this integration technique was used to calculateboth the LPT 1-loop and CLPT power spectra. It was found that therewas only a difference between this method and the method used inthis paper for high 𝑘 -values. We also reached the same conclusionin regards to the CTM power spectrum. The method in Section 2.1.4has numerical advantages as it contains only one infinite sum, henceit was this method that we implemented in the CTM Module.In order to calculate the infinite sum numerically in Equation (31)we truncate the sums at 𝑛 =
32. To reduce the impact of the higher-order spherical Bessel functions on the summation we investigatedthe impact of truncating the sums at 𝑛 =
10 instead. We found thatthis removed some of the numerical noise, however, did not impactthe maximum 𝑘 -value reached before we dropped below 5% of theEuclid Emulator.The spherical Bessel functions in this paper have been calculatedusing the publicly available mcfit. This software is based on the FFT-Log algorithm (Hamilton 2015) and the FFTLog code . Althoughthese codes can be fully optimised to calculate the Zel’dovich powerspectrum, we encountered issues when the correction term in theCTM trajectory becomes large for either large 𝑘 -values or low red-shifts. We leave it to future work to implement an original integrationroutine, fully optimised for the CTM power spectrum. APPENDIX D: APPLICATION OF THE CTM TO KFT
In Ali-Haïmoud (2015) a more detailed computation of the powerspectrum to first-order in the gravitational interaction is given. Thecomputation is also described in Bartelmann et al. (2014a,b) fromthe statistical mechanics perspective. We will simply summarise theresults here so that we may compare our power spectrum to thatpresented in Bartelmann et al. (2014a) and Bartelmann et al. (2014b).Using the definition of the overdensity field (A5) the Dirac delta canbe expanded such that1 + 𝛿 ( ) ( x ) = ∫ 𝛿 D ( x − x ) 𝑑 𝑞, (D1a) 𝛿 ( ) ( x ) = − ∫ x ( q ) · ∇ 𝛿 D ( x − x ) 𝑑 𝑞. (D1b)To first-order in the gravitational interaction (with P ( ) ∝ (cid:104) 𝛿 ( ) 𝛿 ( ) (cid:105) and P ( ) ∝ (cid:104) 𝛿 ( ) 𝛿 ( ) (cid:105) ) the power spectrum is given by,P ( 𝑘 ) ≈ P ( ) ( 𝑘 ) + 𝜖 P ( ) ( 𝑘 ) . (D2)It is then noted that x is equal to the Zel’dovich approxima-tion with the exception of the time dependent function 𝛼 ( 𝑡 ) = + 𝑎 𝑖 (cid:164) 𝐷 ( 𝑡 𝑖 ) 𝐷 ( 𝑡 𝑖 ) ∫ 𝑡𝑡 𝑖 𝑑𝑡 (cid:48) 𝑎 (cid:48) . Therefore, the power spectrum to zeroth orderin the interaction is,P ( ) ( 𝑘 ) = ∫ 𝑑 𝑞 e 𝑖 k · q (cid:104) e − 𝑘 𝛼 𝜎 𝜓 + 𝑘 𝑖 𝑘 𝑗 𝛼 𝜎 𝑖𝑗 ( q ) − (cid:105) . (D3)Calculating P ( ) is more involved as it requires taking the correlationof 𝛿 ( ) and 𝛿 ( ) . It is calculated in full in Ali-Haïmoud (2015),Bartelmann et al. (2014a) and Bartelmann et al. (2014b). However, https://jila.colorado.edu/ ajsh/FFTLog/ let us focus on the result obtained if one expands in terms of thelinear power spectrumP ( ) L ( 𝑘, 𝑧 ) = 𝛼 ( 𝑧 ) P L ( 𝑘, 𝑧 ∗ ) , (D4a)P ( ) L ( 𝑘, 𝑧 ) = 𝜔 𝛼 ( 𝑧 ) ∫ 𝑡𝑡 ∗ 𝑑𝑡 (cid:48) 𝑎 (cid:48) ∫ 𝑡 (cid:48) 𝑡 ∗ 𝑑𝑡 (cid:48)(cid:48) 𝑎 (cid:48)(cid:48) 𝛼 (cid:0) 𝑧 (cid:48)(cid:48) (cid:1) P L ( 𝑘, 𝑧 ∗ ) , (D4b)where 𝜔 = 𝐻 Ω 𝑚 . There are a number of issues raised in Ali-Haïmoud (2015) concerning the results presented in Bartelmann et al.(2014a) and Bartelmann et al. (2014b). One such issue is with theexpansion carried out to calculate the power spectrum. Expandingthe Zel’dovich power spectrum in the usual way gives (Crocce &Scoccimarro 2006a)P zel ( 𝑘 ) ≈ P L ( 𝑘 ) − 𝑘 𝜎 𝜓 P L ( 𝑘 ) + P − loop ( 𝑘 ) (D5)withP − loop ( 𝑘 ) = ∫ 𝑑 𝑘 (cid:48) ( 𝜋 ) P L ( 𝑘 (cid:48) ) P L ( 𝑘 (cid:48)(cid:48) ) 𝑘 (cid:48) 𝑘 (cid:48)(cid:48) (cid:0) k · k (cid:48) (cid:1) (cid:0) k · k (cid:48)(cid:48) (cid:1) . (D6)where k (cid:48)(cid:48) = k − k (cid:48) . The following expansion however, is chosen inBartelmann et al. (2014a) and Bartelmann et al. (2014b)P zel , B14 ( 𝑘 ) ≈ P L ( 𝑘 ) + P − loop + 𝑘 𝜎 𝜓 . (D7)This may lead to enhancement of power on small scales, which iswhat one would expect from the non-linear power spectrum. Notethat there are two free parameters in this approach. There is the timeat which the new trajectory is “switched on”, 𝑧 ∗ and there is thebook-keeping parameter, 𝜖 . In Ali-Haïmoud (2015) 𝑧 ∗ is chosen tobe 𝑧 ∗ =
99 and 𝜖 = 𝐴 ( 𝑧 ) and 𝐵 𝜖 ( 𝑧 ) defined in Equation (A14) for KFT are 𝐴 ( 𝑧 ) = 𝛼 ( 𝑧 ) + 𝜖 CTM 𝜔 𝛽 ( 𝑧 ) , (D8a) 𝐵 𝜖 ( 𝑧 ) = − 𝜖 CTM 𝜔 𝛾 ( 𝑧 ) . (D8b)with 𝛼 ( 𝑧 ) = + 𝑎 𝑖 𝐻 ( 𝑧 𝑖 ) 𝐷 (cid:48) ( 𝑧 𝑖 ) 𝐷 ( 𝑧 𝑖 ) ∫ 𝑧𝑧 𝑖 𝑑𝑧 (cid:48) 𝑎 (cid:48) 𝐻 ( 𝑧 (cid:48) ) , (D9a) 𝛽 ( 𝑧 ) = ∫ 𝑧𝑧 𝑖 𝑑𝑧 (cid:48) 𝑎 (cid:48) 𝐻 ( 𝑧 (cid:48) ) ∫ 𝑧 (cid:48) 𝑧 𝑖 𝑑𝑧 (cid:48)(cid:48) 𝐻 ( 𝑧 (cid:48)(cid:48) ) 𝐷 ( 𝑧 (cid:48)(cid:48) ) 𝐷 ( 𝑧 𝑖 ) , (D9b) 𝛾 ( 𝑧 ) = ∫ 𝑧𝑧 𝑖 𝑑𝑧 (cid:48) 𝑎 (cid:48) 𝐻 ( 𝑧 (cid:48) ) ∫ 𝑧 (cid:48) 𝑧 𝑖 𝑑𝑧 (cid:48)(cid:48) 𝐻 ( 𝑧 (cid:48)(cid:48) ) 𝐷 ( 𝑧 (cid:48)(cid:48) ) 𝐷 ( 𝑧 𝑖 ) 𝛼 (cid:0) 𝑧 (cid:48)(cid:48) (cid:1) . (D9c)The difference between the Beyond Zel’dovich approximation (solidlines) and KFT calculated using the CTM (dashed lines) and theemulator is shown in Figure D1 for 𝑧 = , , , , ,
5. Both theBeyond Zel’dovich approximation and KFT were calculated using 𝜖 CTM = 𝑧 𝑖 =
100 and 𝑘 𝑐 = − . At redshifts, 𝑧 = MNRAS , 1–15 (2021) TM − − k [h Mpc − ] − . − . − . − . − . . . P c a l c ( k ) − P e m u ( k ) P e m u ( k ) z = 0 z = 1 z = 2 z = 3 z = 4 z = 5 Figure D1.
The difference between a given theory and the emulator results isshown for redshifts 𝑧 = , , , , Δ diff ± . using the CTM) when compared to the Euclid Emulator. This islikely due to the time-dependent functions 𝐴 ( 𝑧 ) and 𝐵 𝜖 ( 𝑧 ) beingmarginally larger in the Beyond Zel’dovich approximation.Part of the motivation for the introduction of the CTM and theBeyond Zel’dovich approximation was that KFT does not regainlinear growth on large scales as expected. Thus, in Figure D1 the KFTresults have been re-normalised by a factor of 𝐴 − ( 𝑧 ) (cid:16) 𝐷 ( 𝑧 𝑖 ) 𝐷 ( 𝑧 ) (cid:17) . This paper has been typeset from a TEX/L A TEX file prepared by the author. MNRAS000