GridSPT: Gridbased calculation for perturbation theory of largescale structure
Abstract
Perturbation theory (PT) calculation of largescale structure has been used to interpret the observed nonlinear statistics of largescale structure at the quasilinear regime. In particular, the socalled standard perturbation theory (SPT) provides a basis for the analytical computation of the higherorder quantities of largescale structure. Here, we present a novel, gridbased algorithm for the SPT calculation, hence named GridSPT, to generate the higherorder 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 nonlinear power spectrum and bispectrum. Comparing the density fields (to fifth order) from GridSPT with those from the full Nbody simulations in the configuration space, we find that GridSPT accurately reproduces the Nbody result on large scales. The agreement worsens with smaller smoothing radius, particularly for the underdense regions where we find that 2LPT (secondorder Lagrangian perturbation theory) algorithm performs well.
pacs:
98.80.k, 98.62.Py, 98.65.rI Introduction
Largescale 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 largescale structure therefore contains rich cosmological information. This is why the largescale 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 twopoint correlation function of largescale structure over a gigantic cosmological volume. After an idealistic observation of cosmic microwave background anisotropies by Planck Planck Collaboration et al. (2016), largescale structure would be the best cosmological probe especially to pin down the latetime evolution of the Universe.
In confronting with future precision observations, an accurate modeling of the largescale 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 quasilinear scales Bernardeau et al. (2002). It is at the quasilinear scales where the measurements of baryon acoustic oscillations and redshiftspace distortions, which probe the nature of cosmic acceleration as well as gravity on cosmological scales, are undertaken with largevolume 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 higherorder 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 higherorder 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 resummed 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 orderbyorder 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 facetoface comparison between SPT and the body simulation. We study in detail the properties of the generated density fields from both statistical and morphological pointofview.
This paper is organized as follows. We begin by briefly reviewing the basis of SPT treatment in Sec. II. We then present a FastFourierTransform based method to compute the higherorder density fields from the SPT on grids starting with random density field (Sec. II.3). The results of the gridbased 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 crosscorrelation properties of the gridbased 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 gridbased PT calculations.
Ii Standard perturbation theory
ii.1 Basic equations
The main objective of this study is to find an accurate model for the largescale matter inhomogeneities in the cold dark matter (CDM) dominated Universe. Rigorously speaking, the gravitational evolution of such a system is described by the VlasovPoisson equation in a cosmological background. Starting with a cold initial condition, the motion of dark matter distribution basically follows a singlestream 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))
(1)  
(2)  
(3) 
Note that the singlestream flow approximation is eventually violated in the nonlinear regime where the multistream flow must be included. Although this can basically happen at small scales, and hence the singlestream approximation is expected to give an accurate description of the nonlinear mode coupling in the weakly nonlinear regime, recent studies Bernardeau et al. (2014); Blas et al. (2014); Nishimichi et al. (2016) have advocated that the smallscale multistream 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 singlestream flow on the prediction of the largescale 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
(5) 
In Eq. (4), the quantity is a timedependent matrix:
(6) 
Note that in the Einsteinde Sitter (flat, matterdominated) Universe, this matrix is reduced to
(7) 
ii.2 Realspace formalism
Eq. (4) with (6) is the basis to develop a systematic perturbative expansion. Since we are interested in the latetime evolution of matter fluctuations, we may consider the perturbations dominated by the linear growingmode solution. Further, we shall adopt the socalled Einsteinde Sitter approximation, by which the matrix in Eq. (4) is replaced with the one in the Einsteinde Sitter Universe, Eq. (7). The PT calculation with Einsteinde 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
(8) 
Substituting these into Eq. (4), the orderbyorder calculation leads to (for )
(9) 
This gives the recursion relation for perturbative quantities and :
(10) 
for . For the linearorder quantities (), the growingmode initial condition implies
(11) 
where is the initial density field.
The realspace 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 largescale structure. That is, we express the Fourier transform of the density and velocitydivergence fields, and , as
(12)  
(13) 
with . Here, the field is the initial linear density field, and the functions and are the Fourierspace PT kernels. Substituting these expressions into Eq. (10), one obtains the Fourierspace recursion relation for PT kernels. Writing , we have
(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 Fourierspace 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 largescale structure at field level, the Fourierspace formulation Roth and Porciani (2011) becomes impractical to evaluate the higherorder density and velocity fields because of the multidimensional convolution integrals. We therefore implement the right hand side of Eq. (10) directly in the real space.
ii.3 GridSPT: Generating higherorder density fields on grids
In this subsection, based on the realspace recursion relations, Eqs. (10) and (11), we present a gridbased PT calculation, called GridSPT, which enables us to systematically evaluate the higherorder PT solutions at the field level. Making use of the FastFourier 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:

Generate the initial density field on Fourierspace grids drawn from the Gaussian random distribution specified by the linear power spectrum .

Perform the inverse FFT to obtain the following quantities on realspace grids after multiplying the relevant factors to in Fourier space:
(15) 
Substitute the quantities obtained in the step (2) to the recursion relation (10) to construct the secondorder solutions and on realspace grids.

Perform the forward FFT to have and .

Perform the inverse FFT to obtain the following quantities on realspace grids:
(16) 
Use the recursion relation and quantities obtained so far (i.e., , , , , , and for and ) to construct the thirdorder solutions and on realspace grids.
Repeating the last three steps , this procedure can be generalized to an arbitrary higherorder 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 gridbased 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 velocitydivergence 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 Fouriertransform the nonlinear terms in the right hand side of Eq. (10). For simplicity, consider the two fields, and , in onedimensional grid space of with a grid number . The discrete Fourier transform is described by
(17) 
with and for . Then the Fourier transform of the product, , leads to
(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 highfrequency modes that produces the aliasing contribution. For example, if we force the fields and to zero at , the nonvanishing modes that appear in Eq. (18) are restricted to , and hence the aliasing term never appears. This zeropadding method, referred to as the rule, can be also applied to the threedimensional case (e.g., Orszag (1971)). The method suits well for our purpose because the validity of the PT calculation is basically limited at largescale modes (i.e., low modes).
In what follows, we will demonstrate the gridbased PT calculations with side length of the boxsize Mpc and number of grids . We then eliminating the aliasing effect by applying the isotropic lowpass 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 higherorder 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 gridbased PT calculations, which we hereafter call GridSPT, is written in c++ using FFTW library^{1}^{1}1version 2.1.5; http://www.fftw.org. With the threadparallerized FFT calculation, the code can quickly generates higherorder density fields^{2}^{2}2To be precise, with the CPU of Xeon E52695 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 gridbased 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 flatCDM model: for matter density, for dark energy with equationofstate 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, GADGET2 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 ( largescale modes) for the GridSPT calculation.
iii.1 The morphology of nonlinear 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 Mpcsized 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 higherorder 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 higherorder 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 Nbody and GridSPT, but a closer look at highdensity 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 secondorder 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 lowdensity 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 highdensity 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 fifthorder SPT than secondorder 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 fifthorder 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 unphysical 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 UVsensitive behavior from the statistical pointofview 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 crosscorrelation of the SPT density field with body results, and discuss on how well the PT predictions reproduce the Nbody results.
iii.2.1 Power spectrum and bispectrum
Before discussing the crosscorrelation with body results, let us first check if the GridSPT calculations properly reproduce the analytic prediction computed with Fourierspace 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:
(19)  
(20)  
(21) 
with being the number of Fourier modes in a bin. The corresponding analytic predictions, taking account of the finite box and high cutoff^{3}^{3}3To 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 twoloop 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, higherorder 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 higherorder. Thus, even a small error on each diagram at lower order may result in a big systematics at higherorder 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 twoloop order, while the oneloop 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:
(22)  
(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 oneloop corrections are found to reproduce the analytic PT results (solid lines) quite well. Also, in the right panel, the overall behavior of the oneloop 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 crosscorrelation 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 crosscorrelation coefficients, , defined by
(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 higherorder PT fields with body is not monotonic, and exhibits anticorrelation () for a certain range of .
Indeed, these nonmonotonic 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 multipoint 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 multipoint propagator expansion at twoloop order. In Appendix B, we present a recipe to compute , and derive the analytic expressions for the crosspower spectrum between SPT and nonlinear density fields, valid at the twoloop order.
Apart from a small discrepancy in , RegPT predictions agree well with measured results of crosscorrelation coefficient. The discrepancy in the case is presumably due to the lack of higherloop corrections, although the breakdown of the singlestream 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 nonmonotonic behaviors at lower is originated from the SPT kernels, , in Eq. (12), which can can be both positive and negative. In other words, the nonmonotonic 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 crosscorrelation coefficient defined by
(25) 
As we can see from Fig. 10, adding higherorder 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, higherorder 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 autopower spectrum (see Fig. 6), we see in Sec. III.1 that no severe cancellation of the higherorder 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 decorrelation found in cases is related to the UVsensitive features of higherorder density fields seen in Fig. 5.
To clarify this point, we measure the crossmoments of the realspace density fields smoothed with a Gaussian filter. Varying the smoothing scales, we calculate the crosscorrelation coefficients, which is defined similar to Eq. (25) in realspace, as a function of the smoothing scales in Fig. 11. The scaledependent behaviors seen in Fig. 11 are what is anticipated from the Fourierspace cross correlation coefficients. That is, adding higherorder corrections does not improve the correlation at , and rather leads to a decorrelation 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 secondorder 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 largescale matter flow, at a certain precision, in the singlestream regime. In contrast, the SPT requires a (partial) resummation of an infinite series of PT expansion to describe such an effect. Nevertheless, as shown in Figs. 68, a better performance in the crosscorrelation 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 underpredicts 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 highdensity regions. In other words, at highdensity 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 pointbypoint manner.
iii.2.3 Joint probability distribution of GridSPT and Nbody 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 pointbypoint 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 lowerorder 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:
(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 highdensity 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 overpredicts the amplitude of density fields. These behaviors are indeed what we saw in the power spectrum and bispectrum (Figs. 68), 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 higherorder. In particular, the scatter becomes developed at the lowdensity region, and it extends to the unphysical region, . This is another manifestation of what we see in Fig. 5, This is presumably originated from the UVsensitive modecoupling behavior in the SPT, which explains why the correlation between GridSPT and body results becomes substantially reduced at higherorder, as seen in Figs. 10 and 11.
Iv Conclusion
In this paper, we have presented a novel, gridbased calculation for standard perturbation theory (SPT) of largescale structure, with which a facetoface comparison of the PT prediction with body simulations is made possible. With the FFT, the code, GridSPT, written in c++, can quickly compute the higherorder density fields based on the realspace recursion formula in Eq. (10). Then, we have demonstrated the gridbased 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 higherdensity region, while it conversely gets slower at lowerdensity regions.

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

In contrast, the secondorder Lagrangian PT (2LPT) prediction gives a better correlation with body density field even at the second order. In particular, it reproduces the structure at lowdensity regions remarkably well. However, the Lagrangian PT systematically underestimates the amplitude at higherdensity 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 singlestream 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 effectivefieldtheory 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 singlestream 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 effectivefieldtheory 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 redshiftspace 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 redshiftspace 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 higherorder density fields is rather fast, with the help of a sophisticated algorithm for optimization, the gridbased 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 (AST1517363) and NASA 80NSSC18K1103. Numerical computation was partly carried out at the Yukawa Institute Computer Facility.Appendix A Sensitivity of GridSPT calculations to high 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 nonvanishing spurious contribution (aliasing error) to the gridbased PT calculations through the nonlinear interaction. These modes are to be subtracted, and in this paper, we adopt the isotropic sharpk 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 modecoupling 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 twoloop order and bispectrum prediction at oneloop 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 interparticle distance in our body simulations with Mpc and (that is, these predictions slightly differ from those shown in Figs. 68, but the differences are subtle).
As we see, a small enhances the smallscale 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 Crosscorrelation 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:
(27)  
(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:
(29)  
(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 nonvanishing expressions relevant up to (e.g., Refs. Fry (1994); Carlson et al. (2009); Taruya et al. (2009)):
(31)  
(32)  
(33)  
(34) 
On the other hand, to derive the analytic expression of , we need several steps. First, we use Eq. (12) to express
(35) 
To rewrite the ensemble average at righthand side, we identify with the nonlinear density field which can be dealt with singlestream treatment based on Eqs. (1)(3). While this is the assumption valid in the singlestream regime, the density field is not necessarily treated as small quantity. We then introduce the multipoint propagator Bernardeau et al. (2008). The point propagator, denoted by , is defined as