GridSPT: Grid-based calculation for perturbation theory of large-scale structure

# GridSPT: Grid-based calculation for perturbation theory of large-scale structure

Atsushi Taruya Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8583, Japan    Takahiro Nishimichi Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8583, Japan    Donghui Jeong Department of Astronomy and Astrophysics and Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA
July 14, 2019
###### Abstract

Perturbation theory (PT) calculation of large-scale structure has been used to interpret the observed non-linear statistics of large-scale structure at the quasi-linear regime. In particular, the so-called standard perturbation theory (SPT) provides a basis for the analytical computation of the higher-order quantities of large-scale structure. Here, we present a novel, grid-based algorithm for the SPT calculation, hence named GridSPT, to generate the higher-order density and velocity fields from a given linear power spectrum. Taking advantage of the Fast Fourier Transform, the GridSPT quickly generates the nonlinear density fields at each order, from which we calculate the statistical quantities such as non-linear power spectrum and bispectrum. Comparing the density fields (to fifth order) from GridSPT with those from the full N-body simulations in the configuration space, we find that GridSPT accurately reproduces the N-body result on large scales. The agreement worsens with smaller smoothing radius, particularly for the under-dense regions where we find that 2LPT (second-order Lagrangian perturbation theory) algorithm performs well.

cosmology, large-scale structure
###### pacs:
98.80.-k, 98.62.Py, 98.65.-r
preprint: YITP-18-73

## I Introduction

Large-scale matter inhomogeneities in the Universe, as partly traced by the spatial distribution of galaxies and clusters, are thought to have evolved from tiny fluctuations under the influence of gravity and cosmic expansion. The statistical properties of large-scale structure therefore contains rich cosmological information. This is why the large-scale structure observations has been playing a major role to improve our understanding of the cosmology, and there are ongoing/upcoming observations aiming at precisely measuring the power spectrum and the two-point correlation function of large-scale structure over a gigantic cosmological volume. After an idealistic observation of cosmic microwave background anisotropies by Planck Planck Collaboration et al. (2016), large-scale structure would be the best cosmological probe especially to pin down the late-time evolution of the Universe.

In confronting with future precision observations, an accurate modeling of the large-scale structure is crucial for the unbiased estimation of the cosmological parameters. While the cosmological -body simulation is an essential and standard tool to deal with dark matter and halo clustering even at small scales, the analytic treatment with perturbation theory (PT) is powerful and indispensable to the statistical prediction at quasi-linear scales Bernardeau et al. (2002). It is at the quasi-linear scales where the measurements of baryon acoustic oscillations and redshift-space distortions, which probe the nature of cosmic acceleration as well as gravity on cosmological scales, are undertaken with large-volume galaxy redshift surveys. This is why numerous works on PT treatment have been done (e.g., Jeong and Komatsu (2006); Crocce and Scoccimarro (2008); Bernardeau et al. (2008); Taruya and Hiramatsu (2008); Taruya et al. (2010); Matsubara (2008); Pietroni (2008)). One advantage of the PT calculation is that it provides a way to directly compute the statistical quantities. This is indeed the standard PT (SPT) treatment, in which the higher-order PT kernels of the density and velocity fields (the building blocks of PT calculations defined in the Fourier space) are analytically constructed with recursion relations (e.g., Bernardeau et al. (2002); Goroff et al. (1986)), and are used to give analytical expressions for the higher-order corrections to the power spectrum or bispectrum (but see for Refs. Taruya (2016); Bose and Koyama (2016) for the cases when the analytical construction of PT kernels is intractable). Although the SPT expansion is known to have a poor convergence in predicting statistical quantities (e.g., Crocce and Scoccimarro (2006a); Taruya et al. (2009)), the framework of the SPT treatment is still useful and in fact used in the re-summed PT scheme which improves the convergence of PT expansion Crocce et al. (2012); Taruya et al. (2012).

So far, most of the works with the SPT takes advantage of its analytical basis, and has focused on the direct calculation of the statistical quantities (but see Ref. Roth and Porciani (2011); Tassev (2014)) such as nonlinear matter power spectrum and bispectrum. In this paper, we present a novel GridSPT method to generate the order-by-order nonlinear density fields from a realization of the Gaussian linear density field. Because the Gaussian random field is used to generate the initial conditions for -body simulations, this method makes possible a face-to-face comparison between SPT and the -body simulation. We study in detail the properties of the generated density fields from both statistical and morphological point-of-view.

This paper is organized as follows. We begin by briefly reviewing the basis of SPT treatment in Sec. II. We then present a Fast-Fourier-Transform based method to compute the higher-order density fields from the SPT on grids starting with random density field (Sec. II.3). The results of the grid-based PT calculations up to the fifth order are shown in Sec. III, and the structure and statistics of the density fields are compared with -body simulations in detail. In particular, we study the cross-correlation properties of the grid-based PT and -body density fields in both real and Fourier space, and discuss the pros and cons of the SPT prediction. Finally, Sec. IV is devoted to conclusion and discussion on the possible application of grid-based PT calculations.

## Ii Standard perturbation theory

### ii.1 Basic equations

The main objective of this study is to find an accurate model for the large-scale matter inhomogeneities in the cold dark matter (CDM) dominated Universe. Rigorously speaking, the gravitational evolution of such a system is described by the Vlasov-Poisson equation in a cosmological background. Starting with a cold initial condition, the motion of dark matter distribution basically follows a single-stream flow at an early phase of structure formation, and, in such a regime, the system is described by the pressureless fluid equations coupled with the Poisson equation. Then the basic equations for mass density field become (e.g., Bernardeau et al. (2002))

 ∂δ∂t+1a∇⋅[(1+δ)\boldmathv]=0, (1) ∂\boldmathv∂t+H\boldmath% v+1a(\boldmathv⋅∇)⋅\boldmathv=−1a∇ψ, (2) 1a2∇2ψ=4πGρmδ. (3)

Note that the single-stream flow approximation is eventually violated in the nonlinear regime where the multi-stream flow must be included. Although this can basically happen at small scales, and hence the single-stream approximation is expected to give an accurate description of the non-linear mode coupling in the weakly non-linear regime, recent studies Bernardeau et al. (2014); Blas et al. (2014); Nishimichi et al. (2016) have advocated that the small-scale multi-stream flows can give an impact on the prediction at large scales, and their effects need to be accounted for (see Taruya and Colombi (2017); McDonald and Vlah (2018) for an attempt in 1D). We shall later discuss the quantitative impact of the breakdown of the single-stream flow on the prediction of the large-scale matter distribution.

Adopting Eqs. (1)-(3) as our basic equations, we further impose the irrotationality of the mass flow consistent with the standard picture of structure formation. Then the velocity field follows the potential flow, and is characterized by the scalar quantity. We define with being the linear growth rate defined by with being the linear growth factor. Further, we introduce a new time variable and rewrite Eqs. (1)-(3) with evolution equations for and . We obtain

 (4)

where the quantity is the reduced velocity field defined by . Under the irrotational flow assumption, it is related to the velocity divergence through

 \boldmathu(\boldmathx)=∇[∇−2θ(\boldmathx)]=∫d3\boldmathk(2π)3(−i\boldmathkk2)θ(\boldmath% k)ei\boldmathk⋅\boldmathx. (5)

In Eq. (4), the quantity is a time-dependent matrix:

 Ωab(η)=⎛⎜ ⎜ ⎜ ⎜⎝0−1−4πGρmf2H21f(2+˙HH2+dfdη)⎞⎟ ⎟ ⎟ ⎟⎠. (6)

Note that in the Einstein-de Sitter (flat, matter-dominated) Universe, this matrix is reduced to

 Ωab(η)⟶ΩEdSab=⎛⎜ ⎜ ⎜⎝0−1−3212⎞⎟ ⎟ ⎟⎠. (7)

### ii.2 Real-space formalism

Eq. (4) with (6) is the basis to develop a systematic perturbative expansion. Since we are interested in the late-time evolution of matter fluctuations, we may consider the perturbations dominated by the linear growing-mode solution. Further, we shall adopt the so-called Einstein-de Sitter approximation, by which the matrix in Eq. (4) is replaced with the one in the Einstein-de Sitter Universe, Eq. (7). The PT calculation with Einstein-de Sitter approximation is shown to give a sufficiently accurate prediction in the cosmological model close to the CDM model (e.g., Pietroni (2008); Hiramatsu and Taruya (2009)), and it ensures that the perturbative quantities and can be expanded by a power series of growth factor . We thus have

 δ(\boldmathx)=∑nenηδn(\boldmathx),θ(\boldmathx)=∑nenηθn(\boldmathx),\boldmathu(\boldmathx)=∑nenη\boldmathun(\boldmathx). (8)

Substituting these into Eq. (4), the order-by-order calculation leads to (for )

 (9)

This gives the recursion relation for perturbative quantities and :

 ⎛⎜⎝δn(\boldmathx)θn(\boldmathx)⎞⎟⎠=2(2n+3)(n−1)⎛⎜ ⎜ ⎜ ⎜ ⎜⎝n+12132n⎞⎟ ⎟ ⎟ ⎟ ⎟⎠n−1∑m=1⎛⎜⎝(∇δm)⋅\boldmathun−m+δmθn−m[∂j(\boldmathum)k][∂k(% \boldmathun−m)j]+\boldmathum⋅(∇θn−m)⎞⎟⎠, (10)

for . For the linear-order quantities (), the growing-mode initial condition implies

 ⎛⎜⎝δ1(\boldmathx)θ1(\boldmathx)⎞⎟⎠=⎛⎜⎝11⎞⎟⎠δ0(\boldmathx), (11)

where is the initial density field.

The real-space recursion relation given above contains the gradient and vector fields in the nonlinear source terms. To analytically evaluate these nonlinear terms, a standard way is to go to Fourier space, and obtain the local expression for the Fourier kernels of perturbations. This is what has been done in the statistical predictions of large-scale structure. That is, we express the Fourier transform of the density and velocity-divergence fields, and , as

 δn(\boldmathk)=∫d3\boldmathk% 1⋯d3\boldmathkn(2π)3(n−1)δD(%\boldmath$k$−\boldmathk1⋯n)Fn(\boldmathk1,⋯,\boldmathkn)δ0(\boldmathk1)⋯δ0(\boldmathkn), (12) θn(\boldmathk)=∫d3\boldmathk% 1⋯d3\boldmathkn(2π)3(n−1)δD(%\boldmath$k$−\boldmathk1⋯n)Gn(\boldmathk1,⋯,\boldmathkn)δ0(\boldmathk1)⋯δ0(\boldmathkn), (13)

with . Here, the field is the initial linear density field, and the functions and are the Fourier-space PT kernels. Substituting these expressions into Eq. (10), one obtains the Fourier-space recursion relation for PT kernels. Writing , we have

 F(n)a(\boldmathk1,⋯,% \boldmathkn)=n−1∑m=1σab(n)γbcd(% \boldmathk1⋯m,\boldmathk(m+1)⋯n)F(m)c(\boldmathk1,⋯,\boldmathkm)F(n−m)d(\boldmathkm+1,⋯,\boldmathkn) (14)

with the initial condition, . The explicit functional form of and can be found in e,g., Refs. Goroff et al. (1986); Bernardeau et al. (2002); Crocce and Scoccimarro (2006a); Nishimichi et al. (2007). One advantage of the Fourier-space formulation in Eqs. (12)-(14) is that cosmology dependence are entirely separated out, and the structure of the PT kernels is determined irrespective of the initial conditions and background cosmology (but see Refs. Taruya (2016); Bose and Koyama (2016); Fasiello and Vlah (2016) for generalized cosmologies). However, if one wants to compare the PT prediction with a particular realization of the -body simulation and/or observed large-scale structure at field level, the Fourier-space formulation Roth and Porciani (2011) becomes impractical to evaluate the higher-order density and velocity fields because of the multi-dimensional convolution integrals. We therefore implement the right hand side of Eq. (10) directly in the real space.

### ii.3 GridSPT: Generating higher-order density fields on grids

In this subsection, based on the real-space recursion relations, Eqs. (10) and (11), we present a grid-based PT calculation, called GridSPT, which enables us to systematically evaluate the higher-order PT solutions at the field level. Making use of the Fast-Fourier Transform (FFT), the basic procedure to construct the density and velocity fields on grids are as follows. Note that one can find in the literature some studies closely related to this work, in which the FFT technique has been applied to directly solve Eqs. (1)-(3) Gouda (1995, 1994). Here, we rather stick to a perturbative calculation, and give a recipe to compute the PT solutions on grids:

1. Generate the initial density field on Fourier-space grids drawn from the Gaussian random distribution specified by the linear power spectrum .

2. Perform the inverse FFT to obtain the following quantities on real-space grids after multiplying the relevant factors to in Fourier space:

 δ0(\boldmathk)i\boldmathkδ0(\boldmathk)(−i\boldmathkk2)δ0(\boldmathk)(kikjk2)δ0(% \boldmathk)⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭\lx@stackrelinverseFFT⟹⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩δ1(\boldmathx)=θ1(\boldmathx)∇δ1(\boldmathx)=∇θ1(\boldmathx)\boldmathu1(\boldmathx)∂i(\boldmathu1)j (15)
3. Substitute the quantities obtained in the step (2) to the recursion relation (10) to construct the second-order solutions and on real-space grids.

4. Perform the forward FFT to have and .

5. Perform the inverse FFT to obtain the following quantities on real-space grids:

 i\boldmathkδ2(% \boldmathk)i\boldmathkθ2(\boldmathk)(−i\boldmathkk2)θ2(\boldmathk)(kikjk2)θ2(% \boldmathk)⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭\lx@stackrelinverseFFT⟹⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∇δ2(\boldmathx)∇θ2(\boldmathx)\boldmathu2(\boldmathx)∂i(\boldmathu2)j (16)
6. Use the recursion relation and quantities obtained so far (i.e., , , , , , and for and ) to construct the third-order solutions and on real-space grids.

Repeating the last three steps , this procedure can be generalized to an arbitrary higher-order order in PT solutions. That is, provided the solutions up to the -th order, the -th order solutions, and , are constructed using FFT. In this paper, we will explicitly demonstrate the grid-based PT calculations up to fifth order. Note that in computing -th order, we do not necessarily store all the field data set up to -th order. What is needed is the density and velocity-divergence data up to the -th order, from which we can generate the gradient and vector/tensor fields in Fourier space.

Although the above procedure is simple and thus the code implementation is rather straightforward, one needs to take care of the aliasing effect in practice. The aliasing effect arises when we Fourier-transform the nonlinear terms in the right hand side of Eq. (10). For simplicity, consider the two fields, and , in one-dimensional grid space of with a grid number . The discrete Fourier transform is described by

 a(xj)=N/2−1∑n=−N/2a(kn)eiknxj,a(kn)=1NN/2−1∑j=−N/2a(xj)e−iknxj (17)

with and for . Then the Fourier transform of the product, , leads to

 c(kn)=1NN/2−1∑j=−N/2c(xj)e−iknxj=N/2−1∑l,m=−N/2δKl+m,na(kl)b(km)+N/2−1∑l,m=−N/2δKl+m,n±Na(kl)b(km) (18)

with being the Kronecker delta. In the right hand side of Eq. (18), the first term is the signal that we want to calculate, and second term is the aliasing contribution coming from the discrete sampling. A simple way to eliminate the aliasing effect is to discard the high-frequency modes that produces the aliasing contribution. For example, if we force the fields and to zero at , the non-vanishing modes that appear in Eq. (18) are restricted to , and hence the aliasing term never appears. This zero-padding method, referred to as the rule, can be also applied to the three-dimensional case (e.g., Orszag (1971)). The method suits well for our purpose because the validity of the PT calculation is basically limited at large-scale modes (i.e., low- modes).

In what follows, we will demonstrate the grid-based PT calculations with side length of the boxsize Mpc and number of grids . We then eliminating the aliasing effect by applying the isotropic low-pass filter (called sharp- filter) in the Fourier space to all variables. Note that strict rule would suggest to set the modes with  Mpc to zero. After carefully investigating the impact of cutoff scales on the statistical results in Appendix A, we find that employing the two different cutoff wavenumbers performs better. As the default setup, we apply the sharp- filter of  Mpc to the linear density fields, and then apply the same sharp- filter with  Mpc to the higher-order density fields. Though this exceeds the critical wavenumber , the resultant PT calculations give the best performance. A possible reason for this is discussed in Appendix A.

Implementing this filtering technique, the code of grid-based PT calculations, which we hereafter call GridSPT, is written in c++ using FFTW library111version 2.1.5; http://www.fftw.org. With the thread-parallerized FFT calculation, the code can quickly generates higher-order density fields222To be precise, with the CPU of Xeon E5-2695 2.1GHz (36 cores) and using the Intel compiler and 36 threads for FFT, it takes roughly two minutes to generate density fields up to fifth order..

## Iii Properties of SPT density fields

Here, we present the results of grid-based PT calculations up to fifth order, and in comparison with -body simulations, we investigate in detail the properties of generated density fields. The cosmological parameters used to generate linear density field as well as to run the -body simulation are determined by WMAP5 results Komatsu et al. (2009) assuming the flat-CDM model: for matter density, for dark energy with equation-of-state parameter , for baryon fraction, for Hubble parameter, for scalar spectral index, and finally, for the normalization of the fluctuation amplitude at Mpc. The cosmological -body simulation is carried out by publicly available code, GADGET-2 Springel (2005), with particles in comoving periodic cubes of Mpc, starting with the initial density field calculated from the 2LPT code Crocce et al. (2006) at redshift . Note that we use the same initial density field ( large-scale modes) for the GridSPT calculation.

### iii.1 The morphology of non-linear density fields

Let us first present the morphology of density fields generated by GridSPT, and investigate their properties, comparing with -body results. Figs. 2 and 2 show the 2D slices of the density fields at from GridSPT (with the order label) and from -body results (bottom right). In Fig. 2, the density fields over the entire box are shown with the amplitude plotted in linear scale. On the other hand, Fig. 2 shows the local density fields particularly taken from Fig. 2 over the Mpc-sized region enclosed by dashed line in bottom right panel. In both figures, we apply the Gaussian filter of the radius Mpc to the results, and on each grid, the density fields are averaged over the Mpc depth. While the top left panel in Figs. 2 and 2 are the linear density field (labeled as ), the four successive panels (i.e., top middle, top right, bottom left, and bottom middle) are the GridSPT results summing up higher-order density fields, i.e., with the number indicated in each panel. For comparison, the last panel (bottom right) is the -body result obtained from the same initial density field. Adding higher-order PT corrections, the resultant density fields tend to show a clearer contrast, and resemble the -body density field. At the fifth order, the PT density field seems to almost match the -body result from a visual inspection. It is indeed hard from Fig. 2 to discriminate between the N-body and GridSPT, but a closer look at high-density region seen in Fig. 2 reveals that there is a slight mismatch between the PT prediction (at the fifth order) and the simulation result.

To scrutinize the difference between the two, we select representative regions from the local patch in Fig. 2 (indicated in horizontal and vertical dashed lines), and plot the 1D density fields in Figs. 4 and 4. The left panels plot the GridSPT density fields, with (indicated in the panel), while the right panels compare the result of with those obtained from the -body simulation and second-order Lagrangian PT with 2LPT code. As increasing the order of perturbative expansion, the density peak becomes steepened, and the amplitude gets increased at high density regions, and the low-density regions are flattened conversely. A notable point may be that the convergence of the PT predictions tends to be slow at low density regions. This implies that the PT prediction is generally poor to describe the bulk matter flows, consistent with what has been found in Ref. Roth and Porciani (2011); Tassev (2014). On the other hand, the convergence at high-density regions looks faster, and at the fifth order, the PT prediction reproduces the -body density fields remarkably well. Although there is a general trend that the predicted peak amplitude slightly overestimates the -body results, the overall structure of density peaks is better described by the fifth-order SPT than second-order Lagrangian PT. Nevertheless, this does not imply that the statistical correlation of SPT with -body density field is better than that of Lagrangian PT. This point will be discussed in detail in Sec. III.2.2 and III.2.3.

Finally, we note here that the resultant density fields obtained from the GridSPT calculation is rather sensitive to the choice of the smoothing scale. Fig. 5 plots the 1D density fields at the same regions as shown in Figs. 4 and 4, but the fifth-order PT results at different smoothing scales are depicted as different colors. A slight decrease of the smoothing scale results in a spurious wiggle structure at low density region, and this can also affect small density peaks, leading to an un-physical behavior of . We have checked that this behavior appears persistently in the GridSPT density fields smoothed with the same scales, regardless of the box size and the aliasing correction. The fact that the density field is less than is thus not simply due to the numerical artifact, if any, but rather a drawback of Eulerian perturbation theory; for the region with , continuity equation yields the violation of mass conservation, namely, density contrast increases when the velocity field is divergent. The explicit violation of mass conservation exhibits a rather strong mode coupling between long and short modes. We will discuss this UV-sensitive behavior from the statistical point-of-view in the next section.

### iii.2 Statistical properties

In this subsection, based on the density fields generated with GridSPT in Sec. III.1, we measure the statistical quantities, and study their properties in comparison with the result of the -body simulation. In particular, we measure the cross-correlation of the SPT density field with -body results, and discuss on how well the PT predictions reproduce the N-body results.

#### iii.2.1 Power spectrum and bispectrum

Before discussing the cross-correlation with -body results, let us first check if the GridSPT calculations properly reproduce the analytic prediction computed with Fourier-space recursion relation in Eq. (14). Fig. 6 shows the measured results of the power spectrum from the GridSPT calculations, which are compared with both analytic PT predictions and -body results. Here, the results at are particularly shown. All the results are multiplied by () for the power spectrum (bispectrum).

In left panel of Fig. 6, the contribution to the power spectrum at each perturbative order, (red), (green), and (blue), are shown. These are measured from the GridSPT density fields analogously to the standard procedure on the snapshots of -body simulations:

 Plin(k) =D2+1Nk∑|\boldmathk|=k|δ1(\boldmathk)|2, (19) P1-loop(k) (20) P2-loop(k) (21)

with being the number of Fourier modes in a -bin. The corresponding analytic predictions, taking account of the finite box and high- cutoff333To be precise, in analytic PT calculations, we introduced the cutoff scales in the linear power spectrum, given by  Mpc and  Mpc. , are depicted as solid lines. We expect that the GridSPT results would converge to the analytic curves after we take an ensemble average over many different random realizations of the GridSPT density field.

Overall, the measured results reasonably agree with the analytic predictions. However, a closer look at small scales reveals a small discrepancy between GridSPT and analytic calculations in both and . In particular, the discrepancy is manifest in the two-loop correction, and GridSPT calculations are prone to overestimate the analytic prediction. This could happen probably due to the accumulation of small numerical flaw at each order in the GridSPT calculations. Indeed, in SPT, higher-order correction of the power spectrum is known to have a heavy cancellation among multiple diagrams at the same order with positive and negative contributions, and the cancellation becomes more significant as we go to higher-order. Thus, even a small error on each diagram at lower order may result in a big systematics at higher-order corrections through an imperfect cancellation.

In right panel of Fig. 6, summing up all the PT corrections, the GridSPT results are shown in green and blue filled circles, which are compared with analytic PT predictions depicted as solid lines. As anticipated, discrepancy is manifest at two-loop order, while the one-loop results show a tiny amount of error at high-, which apparently looks insignificant. Although this point has to be kept in mind in our subsequent analyses, the discrepancy is large only at the scales where the deviation from the -body result, depicted as filled red circles, is significant. Further, the discrepancy remains mild, and is smaller than a large underestimation found in the prediction based on 2LPT (filled gray circles).

Indeed, such a discrepancy is not clearly seen in the case of the bispectrum. Fig. 8 presents the measured results of the PT contribution at each order, and , (left) and their total amplitudes (right), which are compared with -body and 2LPT results. Here, the measurement of the bispectrum are done in the equilateral configuration, , and the results are plotted as function of . The PT corrections and are defined as:

 Btree(k1,k2,k3) =D4+1N123∑\boldmathk1,\boldmathk2,\boldmathk3δK% \boldmathk1+\boldmathk2+\boldmathk3,% \boldmath0[δ1(\boldmathk1)δ1(% \boldmathk2)δ2(\boldmathk3)+(2 perm.)], (22) B1-loop(k1,k2,k3) =D6+1N123∑\boldmathk1,\boldmathk2,\boldmathk3δK% \boldmathk1+\boldmathk2+\boldmathk3,% \boldmath0[{δ1(\boldmathk1)δ2(\boldmathk2)δ3(\boldmathk3)+(5 perm.)} +δ2(\boldmathk1)δ2(% \boldmathk2)δ2(\boldmathk3)+{δ1(\boldmathk1)δ1(\boldmathk2)δ4(% \boldmathk3)+(2 perm.)}]. (23)

Note that in actual calculation of these expressions, we use a fast estimator based on FFT (e.g., Baldauf et al. (2015); Sefusatti et al. (2016)). Only with the single realization, the resultant bispectrum is rather noisy, however, increasing the number of realizations up to , shown in Fig. 8, the tree- and one-loop corrections are found to reproduce the analytic PT results (solid lines) quite well. Also, in the right panel, the overall behavior of the one-loop prediction better agrees with -body simulation than the 2LPT prediction (gray filled circles).

#### iii.2.2 Cross correlation

In order to systematically compared the GridSPT density fields with -body simulations, we calculate the cross-correlation between them. First, we consider the density field at each PT order, and compute the cross correlation with the density field constructed from -body simulations in Fourier space. Fig. 10 shows the cross-correlation coefficients, , defined by

 r(n)corr(k)≡∑|\boldmathk|=kRe[δn(\boldmathk)δN-body(\boldmathk)]√∑|\boldmathk|=k|δn(\boldmathk)|2∑|\boldmathk|=k|δN-% body(\boldmathk)|2. (24)

Note that . The measured results at (left) and (right) are shown up to the fifth order of GridSPT density fields.

As shown in Fig. 10, all the correlations tend to get suppressed at high-, and, at , the suppression of correlation appears prominent even on large scales. This is indeed what is expected. An interesting point may be that the correlation of the higher-order PT fields with -body is not monotonic, and exhibits anti-correlation () for a certain range of .

Indeed, these non-monotonic behaviors, together with a strong damping at high-, are predicted by a resummed PT treatment. In Fig. 10, we also show (solid lines) the predictions of the same quantity based on the multi-point propagator expansion Bernardeau et al. (2008) with regularized propagators, called RegPT Taruya et al. (2012). The RegPT predictions are made by assuming that the evolved density field in -body simulation is described with the multi-point propagator expansion at two-loop order. In Appendix B, we present a recipe to compute , and derive the analytic expressions for the cross-power spectrum between SPT and nonlinear density fields, valid at the two-loop order.

Apart from a small discrepancy in , RegPT predictions agree well with measured results of cross-correlation coefficient. The discrepancy in the case is presumably due to the lack of higher-loop corrections, although the breakdown of the single-stream PT treatment might possibly play a role. In any case, an important consequence of this agreement is that the suppression of the measured correlation at high- comes from the randomness of the linear displacement field, which can be described with RegPT by a (partial) resummation of an infinite series of SPT expansion in the high- limit. This explains why the GridSPT density field with the naive SPT expansion looses quickly the correlations with the fully nonlinear density field at small scales. On the other hand, the non-monotonic behaviors at lower- is originated from the SPT kernels, , in Eq. (12), which can can be both positive and negative. In other words, the non-monotonic behavior of directly manifests the poor convergence of the SPT expansion, and this explains why summing up each order of perturbation improves very slowly, as shown in Fig. 10, where we plot the cross-correlation coefficient defined by

 R(n)corr(k)≡∑|\boldmathk|=kRe[(∑ni=1Dn+δi(\boldmath% k))δN-body(\boldmathk)]√∑|\boldmathk|=k|∑ni=1Dn+δi(% \boldmathk)|2∑|\boldmathk|=k|δN-body(\boldmathk)|2. (25)

As we can see from Fig. 10, adding higher-order density fields recover the correlation at large scales to some extent, and the RegPT predictions, depicted as solid line, explain the measured results quantitatively up to . At , however, higher-order corrections rather worsen the correlation with -body simulation at intermediate scales around  Mpc. While this behavior might be possibly originated from small numerical flows as observed in the auto-power spectrum (see Fig. 6), we see in Sec. III.1 that no severe cancellation of the higher-order terms has occurred at the density fields. Though the impact of numerical flaws inferred from Fig. 6 may result in a % systematic error, we rather suspect that substantial de-correlation found in cases is related to the UV-sensitive features of higher-order density fields seen in Fig. 5.

To clarify this point, we measure the cross-moments of the real-space density fields smoothed with a Gaussian filter. Varying the smoothing scales, we calculate the cross-correlation coefficients, which is defined similar to Eq. (25) in real-space, as a function of the smoothing scales in Fig. 11. The scale-dependent behaviors seen in Fig. 11 are what is anticipated from the Fourier-space cross correlation coefficients. That is, adding higher-order corrections does not improve the correlation at , and rather leads to a de-correlation at small scales. The notable point is its characteristic scale, i.e.,  Mpc at and  Mpc at . The latter indeed matches the scales where we found a spurious wiggle feature in Fig. 5. We thus think that the behaviors seen in Figs. 10 and 11 are real, and capture the limitation of SPT. This point will be further discussed in detail in Sec. III.2.3.

Finally, as a reference, we compare the results with those between second-order Lagrangian PT (2LPT) and -body simulation, which are plotted in magenta filled circles in Figs. 10 and 11. We find that the correlation obtained from 2LPT is much better than those from SPT, indeed, at every redshift and on every scale. This is again due to the fact that the density fields generated with Lagrangian PT is constructed based on the displacement of the particles which respects the mass conservation (that is, Lagrangian PT ensures that ). It can therefore describe the large-scale matter flow, at a certain precision, in the single-stream regime. In contrast, the SPT requires a (partial) re-summation of an infinite series of PT expansion to describe such an effect. Nevertheless, as shown in Figs. 6-8, a better performance in the cross-correlation coefficients does not necessarily imply that the Lagrangian PT always gives a better statistical prediction for the mass distribution. Indeed, the density fields generated with 2LPT generally under-predicts the amplitude of peaks (see Fig. 4), and this can result in the underestimation of the power spectrum and bispectrum amplitudes. In this respect, SPT suits better for the statistical predictions sensitive to the high-density regions. In other words, at high-density regions, the statistical correlation of SPT with -body simulation may be better than that of 2LPT. We shall show it below by the direct measurement of the correlation in the point-by-point manner.

#### iii.2.3 Joint probability distribution of GridSPT and N-body simulation

Our last measurement is the joint probability distribution function (PDF) for the local density contrast obtained both from the GridSPT calculation, , and -body simulation, . This quantity, denoted by , characterizes the probability that at a given position in the space, the overdensity from -body simulation is , and the GridSPT calculation predicts . It thus gives the point-by-point comparison of the two density fields.

Fig. 12 shows the measured result of the joint PDF for smoothed density fields with Gaussian filter at . Here, we show the results with the smoothing radius Mpc. In each panel, the results with GridSPT are presented up to the fifth order (i.e., ). For reference, the joint PDF between 2LPT and -body density fields is also shown (lower right). Compared to the results with 2LPT, the lower-order results of the GridSPT calculations show not only a large scatter but also a declined correlation feature in the joint PDF. Here, the magenta line in each panel represents the conditional mean for a given -body density field, , defined by:

 ¯¯¯δSPT(δN-body)=∫dδSPTδSPTProb(δSPT,δN-body)Prob(δN-body). (26)

As increasing (i.e., the order of perturbation), the correlation gets tighter, and the conditional mean tends to follow a linear relation depicted as black solid line, meaning that the GridSPT density field gets closer to the -body density field. Overall, the result with 2LPT still looks much better as the correlation is much tighter and closer to the linear relation. However, a closer look at the high-density region reveals a small tilt i.e., , where the GridSPT results at are apparently better correlated with the -body results. Beyond this regime, on the other hand, the conditional mean of the GridSPT becomes steeper, and thus the GridSPT tends to over-predicts the amplitude of density fields. These behaviors are indeed what we saw in the power spectrum and bispectrum (Figs. 6-8), and partly explains why SPT gives a better prediction.

Finally, as we discussed in Sec. III.1, the measured result of joint PDF is also sensitive to the choice of filter scales. Fig. 13 shows the same plot as in Fig. 12, but the results of the Gaussian filter of the radii (top) and Mpc (bottom) are shown. Decreasing the filter scale, a sizable amount of scatter appears in the joint PDF of the GridSPT calculation at higher-order. In particular, the scatter becomes developed at the low-density region, and it extends to the un-physical region, . This is another manifestation of what we see in Fig. 5, This is presumably originated from the UV-sensitive mode-coupling behavior in the SPT, which explains why the correlation between GridSPT and -body results becomes substantially reduced at higher-order, as seen in Figs. 10 and 11.

## Iv Conclusion

In this paper, we have presented a novel, grid-based calculation for standard perturbation theory (SPT) of large-scale structure, with which a face-to-face comparison of the PT prediction with -body simulations is made possible. With the FFT, the code, GridSPT, written in c++, can quickly compute the higher-order density fields based on the real-space recursion formula in Eq. (10). Then, we have demonstrated the grid-based PT calculation up to the fifth order, and studied the morphological and statistical properties of the SPT density fields in comparison with -body simulations and Lagrangian PT predictions.

Our findings are summarized as follows:

• Increasing the PT order up to more than third order, the SPT reproduces the structure of density fields in -body simulation reasonably well on large scales ( Mpc). In particular, the convergence of the PT expansion is found to be faster at higher-density region, while it conversely gets slower at lower-density regions.

• On the other hand, on small scales, the SPT is prone to produce spurious and un-physical structure in the density field. The statistical correlation with -body simulations is generally poor on small scales ( Mpc or Mpc at ), and including the higher-order corrections does not improve the cross correlation at all.

• In contrast, the second-order Lagrangian PT (2LPT) prediction gives a better correlation with -body density field even at the second order. In particular, it reproduces the structure at low-density regions remarkably well. However, the Lagrangian PT systematically underestimates the amplitude at higher-density regions, and this can lead to a poor statistical description of the power spectrum and bispectrum. Rather, the SPT gives a better statistical prediction, and albeit the lack of precision, it can qualitatively capture the trend of nonlinear growth.

The deficiencies of the SPT predictions certainly come from the facts that the present expansion scheme does not respect the mass conservation at finite order (that is, regions with exist), and the single-stream approximation in the PT treatment is invalid on small scales. Our findings are thus regarded as a direct manifestation of these deficiencies both at field and statistical levels. To remedy these drawbacks, resummed PT scheme helps the convergence of PT expansion (e.g., Crocce and Scoccimarro (2006a); Crocce and Scoccimarro (2008); Bernardeau et al. (2008); Taruya and Hiramatsu (2008); Pietroni (2008); Taruya et al. (2009); Crocce et al. (2012); Taruya et al. (2012)), and can mitigate the impact of violating the mass conservation. Further, the effective-field-theory approach (e.g., Baumann et al. (2012); Carrasco et al. (2012); Hertzberg (2014); Baldauf et al. (2015)) would be crucial to remedy the UV sensitive behavior in SPT arising from the violation of single-stream flows. Although various works have been so far done, most of the works has been made at the level of statistical quantities, not at the field level. Development of resummation scheme and effective-field-theory treatment at the field level is thus rather interesting subject.

Finally, although the present GridSPT treatment is still primitive and needs to be improved for the practical purposes along this line, possible application of GridSPT calculation may be worthy of consideration. One interesting application is obviously a direct comparison with observed galaxy distributions, after implementing appropriate prescriptions for galaxy bias and redshift-space distortions. There has been recently a general description of the galaxy bias exploited on the basis of the SPT Desjacques et al. (2018a); Mirbabayi et al. (2015); Desjacques et al. (2018b), and thus the implementation of galaxy bias is straightforward at the field level. Also, the redshift-space distortions are described by a simple mapping formula, and thus easy to handle at the field level. Taking one step further, the method may be applied to the reconstruction problem to infer the initial conditions of the universe (e.g., Jasche and Wandelt (2013); Wang et al. (2013); Kitaura (2013); Wang et al. (2014); Schmittfull et al. (2017); Seljak et al. (2017); Shi et al. (2018); Modi et al. (2018) for recent works). Since the generation of higher-order density fields is rather fast, with the help of a sophisticated algorithm for optimization, the grid-based PT approach may have a potential to efficiently extract the cosmological information at weakly nonlinear regime, rather than the statistical PT calculation.

###### Acknowledgements.
This work was supported in part by MEXT/JSPS KAKENHI Grant Number JP15H05899 and JP16H03977 (AT), and JP17K14273 (TN). TN also acknowledges financial support from Japan Science and Technology Agency (JST) CREST Grant Number JPMJCR1414. DJ acknowledges support from NSF grant (AST-1517363) and NASA 80NSSC18K1103. Numerical computation was partly carried out at the Yukawa Institute Computer Facility.

## Appendix A Sensitivity of GridSPT calculations to high-k cutoff

In this Appendix, we investigate the impact of the high- cutoff for the sharp- filter on the GridSPT calculations, particularly paying an attention to the systematic error associated with the aliasing effect.

As we discussed in Sec. II.3, high- modes of  Mpc can produce the non-vanishing spurious contribution (aliasing error) to the grid-based PT calculations through the non-linear interaction. These modes are to be subtracted, and in this paper, we adopt the isotropic sharp-k filter, with which modes with are set to zero at each order of PT calculation. In general, a smaller value of the cutoff wavenumber is preferred for a secure removal of the aliasing error, but this would change the mode-coupling behavior, also affecting the power spectrum and bispectrum predictions. It is thus important to find an optimal cutoff scale.

To see how the choice of cutoff scale affects the GridSPT results, we plot in Figs. 15 and 15 respectively the power spectrum prediction at two-loop order and bispectrum prediction at one-loop order at , varying the cutoff scales. In each figure, left panels examine the cases with , while right panels show the results with , whose amplitudes are plotted in logarithmic scales. For references, gray dashed lines are the analytic PT predictions with the high- cutoff  Mpc, corresponding to the Nyquist frequency determined by the inter-particle distance in our -body simulations with Mpc and (that is, these predictions slightly differ from those shown in Figs. 6-8, but the differences are subtle).

As we see, a small enhances the small-scale amplitudes of the power spectrum and bispectrum. This is expected from the fact that the mode transfer generally occurs from large to small scales in the CDM spectrum, and introducing the high- cutoff makes the mode transfer at smaller scales ineffective compared to that at larger scales (see, e.q., Ref. Nishimichi et al. (2016)). On the other hand, the results with larger exhibit a strong enhancement in amplitudes at large scales, and a slight change of the cutoff scale largely alters the low- behaviors in power spectrum and bispectrum. While this behavior is expected and is basically explained by the spurious aliasing contributions, a notable point is that the effect appears prominent at  Mpc, which is larger than the critical wavenumber  Mpc.

A possible reason for this may be that we use the isotropic filter, and even if the cutoff wavenumber slightly exceeds , majority of the high- modes that can produce the aliasing error is set to zero except the modes close to the axes. Although the isotropic sharp- filter ceases to be effective at , a part of the aliasing error is still suppressed if the cutoff scale is chosen below . Based on this consideration, one may adopt the cutoff scale of the isotropic filter larger than . Indeed, we find empirically that applying the sharp- filter with  Mpc for the linear density field and  Mpc for the PT density fields higher than second order, the GridSPT calculations give the best performance to the power spectrum and bispectrum predictions amongst those we examined, depicted as filled black symbols in Figs. 15 and 15. We have also compared other quantities such as cross correlation coefficients, and , and joint PDF with those obtained with the cutoff scale , and checked that all the results are hardly changed. Hence, we decided to adopt this choice as default setup and presented the GridSPT results.

## Appendix B Cross-correlation coefficient from RegPT

In this Appendix, based on the resummed PT calculation with RegPT, we describe a prescription to compute and shown in solid lines of Figs. 10 and 10 is presented.

Let us first rewrite Eqs. (24) and (25) with those in the continuum limit:

 r(n)corr(k)=P(n)cross(k)P(nn)SPT(k)PN-body(k), (27) R(n)corr(k)=∑ni=1P(i)cross(k)∑ni≤jP(ij)SPT(k)PN-body(k). (28)

In the above, is the power spectrum of the evolved density field in -body simulations, for which we use the measured data. On the other hand, the quantities and , which involve the PT density fields, are defined by:

 ⟨δn(\boldmathk)δN-body(\boldmathk′)⟩=(2π)3δD(% \boldmathk+\boldmathk′)P(n)cross(k), (29) ⟨δi(\boldmathk)δj(% \boldmathk′)⟩=(2π)3δD(\boldmathk+\boldmathk′)P(ij)SPT(k), (30)

These are what we want to deal with based on the analytic PT treatment. Below, we present their analytic expressions up to .

The quantity, , appears frequently in the statistical calculation of power spectrum in SPT. With the expression given at Eq. (12), the Gaussianity of the initial density field leads to the following non-vanishing expressions relevant up to (e.g., Refs. Fry (1994); Carlson et al. (2009); Taruya et al. (2009)):

 P(11)SPT(k) =P0(k) (31) P(22)SPT(k) =2∫d3\boldmathp(2π)3{F2(% \boldmathk,\boldmathk−\boldmathp)}2P0(p)P0(|\boldmathk−\boldmathp|), (32) P(13)SPT(k) =6P0(k)∫d3\boldmathp(2π)3F3(\boldmathk,\boldmathp,−\boldmathp)2P0(p), (33) P(33)SPT(k) =6∫d3\boldmathpd3\boldmathq(2π)6{F3(\boldmathp,\boldmathq,\boldmathk−\boldmathp−\boldmathq)}2P0(p)P0(q)P0(|% \boldmathk−\boldmathp−\boldmathq|)+9P0(k)[∫d3\boldmathp(2π)3F3(\boldmathk,\boldmathp,−\boldmathp)P0(p)]2. (34)

On the other hand, to derive the analytic expression of , we need several steps. First, we use Eq. (12) to express

 ⟨δn(\boldmathk)δN-body(\boldmathk′)⟩ =∫d3\boldmathp1⋯d3% \boldmathpn(2π)3(n−1)Fn(\boldmathp1,⋯,\boldmathpn)⟨δ0(\boldmathp1)⋯δ0(\boldmathpn)δN-body(% \boldmathk′)⟩. (35)

To rewrite the ensemble average at right-hand side, we identify with the nonlinear density field which can be dealt with single-stream treatment based on Eqs. (1)-(3). While this is the assumption valid in the single-stream regime, the density field is not necessarily treated as small quantity. We then introduce the multi-point propagator Bernardeau et al. (2008). The -point propagator, denoted by , is defined as

 1p!⟨δpδN-body(\boldmathk,η)δδ0(\boldmathk1)⋯δδ0(\boldmathkp)⟩=δD(% \boldmathk−\boldmathk1⋯p)1(2π)3(p−1)Γ