Fourth-order perturbative equations in Lagrangian perturbation theory for a cosmological dust fluid

We have derived fourth-order perturbative equations in Lagrangian perturbation theory for a cosmological dust fluid. These equations are derived under the supposition of Newtonian cosmology in the Friedmann-Lema\^{i}tre-Robertson-Walker Universe model. Even if we consider the longitudinal mode in the first-order perturbation, the transverse mode appears in the third-order perturbation. Furthermore, in this case, six longitudinal-mode equations and four transverse-mode equations appear in the fourth-order perturbation. The application of the fourth-order perturbation leads to a precise prediction of the large-scale structure.

Recently, resummation theories have been proposed; 26), 27), 29), 28), 30) these theories reorganize an infinite series of diagrams into the form of standard perturbation theory and partially resum them. Resummation theories can precisely describe physical quantities such as power spectra and correlation functions. In particular, theories based on Lagrangian perturbation have several merits: they make is easy to describe physical quantities in not only real space but also in redshift space 28), 32), 33) and nonlocal bias is included. 31), 32) A formula for arbitrary-order (m-loop) correction has been proposed. 33) If we have greater knowledge of higher-order Lagrangian perturbations, we can further improve resummation theories that are based on Lagrangian perturbation. For example, recently, the matter bispectrum with higher-order Lagrangian perturbations was computed. 34), 35) Although perturbation theory cannot describe strongly nonlinear regions, it still plays an important role. Cosmological N -body simulations have been applied to describe the evolution of strongly nonlinear regions. The initial condition for the simulations is given by Lagrangian perturbation. It is known that the accuracy of the initial condition affects time evolution in strongly nonlinear regions. 36), 37), 38), 39) In the future, this knowledge would be important for precision cosmology.
The application of Lagrangian perturbation theory (LPT) to observational results will also become important. For example, the Sloan Digital Sky Survey (SDSS) Data Release 7 (DR7) contains spectroscopic data on approximately one million galaxies. 40) In other words, the SDSS project provides a three-dimensional map of the large-scale structure. Because the magnitude limit of the project is approximately 22, the survey is carried out within the low-z region. In future surveys, high-z (z > 1) galaxies are expected to be observed. It is expected that the results of these surveys will provide us with knowledge on the evolution of the large-scale structure. Furthermore, because the high-z region is still quasi-nonlinear, higher-order LPT will be valid. We expect that very precise predictions regarding the large-scale structure can be made by incorporating fourth-order perturbation into resummation theory. This paper is organized as follows. In Sec. 2, we briefly introduce the basic equations for Lagrangian perturbation; in this study, we consider only a dust fluid. Next, we present the perturbative equations in Sec. 3. The fourth-order perturbative equations consist of six longitudinal-mode (scalar) equations and four transversemode (vector) equations. We also derive the perturbative solutions in the Einsteinde Sitter (E-dS) Universe model. Every mode grows at the same rate; thus, it is not possible to ignore any part of these modes. Next, the application of these equations to our results is considered. In Sec. 4, we discuss the applications of the higher-order perturbations. In Sec. 4.1, we describe one possible outcome of applying our result to resummation theories. In Sec. 4.2, we consider initial condition problems in cosmological N -body simulations. Finally, we summarize our findings in Sec. 5. In Appendix A, we describe the relation between the determinant of the displacement matrix and the combination of triplet terms. Appendix B presents a detailed derivation of the fourth-order perturbative equations. Appendix C presents higher-order loop correction for the power spectrum. §2. Basic equations In this section, we briefly introduce the Lagrangian description of a Newtonian cosmological fluid. When the scale of objects is smaller than that of the horizon, this description is valid. The effect of cosmic expansion is reflected by the scale factor a. The solution a is given by the Friedmann equations (or alternative equations). In the comoving coordinates, the basic equations for a cosmological dust fluid are In Eulerian perturbation theory, the density fluctuation δ is regarded as a perturbative quantity. However, in LPT, displacement from a homogeneous distribution is considered: where x and q are the comoving Eulerian coordinates and the Lagrangian coordinates, respectively. s is the displacement vector, which is regarded as a perturbative quantity. From Eq. (2 . 6), we can solve the continuous equation (2 . 1) exactly. Then, the density fluctuation is given in the formally exact form J refers to the Jacobian of the coordinate transformation from Eulerian x to Lagrangian q. Therefore, when we derive the solutions for s, we can determine the evolution of the density fluctuation. The peculiar velocity is given as follows: Then, we introduce the Lagrangian time derivative Taking the divergence and rotation of Eq. (2 . 2), we obtain the evolution equations for the Lagrangian displacement: where (˙) refers to the Lagrangian time derivative (Eq. (2 . 9)). To solve the perturbative equations, we decompose the Lagrangian perturbation into its longitudinal and transverse modes: (2 . 13) where ∇ refers to the Lagrangian spatial derivative.

First-order perturbative equations
The first-order perturbative equation in the longitudinal mode was derived by Zel'dovich. 9) The evolution equation is given as follows: The solution can be divided into its spatial and temporal parts: In the E-dS Universe model (a(t) ∝ t 2/3 ), the solutions to the temporal part are given as Hereafter, we consider only the growing solution (g 1 (t) = t 2/3 ). The spatial part ψ (1) (q) is given by the initial condition. The effect of the transverse mode was examined by Buchert 17) and Barrow and Saich: 18) We choose a reasonable boundary condition for this equation. The spatial differential operator can be removed, leaving In the E-dS Universe model, the solution is given as Because there are no growing solutions for first-order perturbation in the transverse mode, we consider only the longitudinal mode in the first-order perturbation.

Second-order perturbative equations
The second-order perturbative equation in the longitudinal mode was derived by Bouchet et al. 12) and Buchert and Ehlers. 13) The evolution equation for the second-order longitudinal mode is derived from Eq. (2 . 10). This equation is given as follows:S (2) ,ii + 2ȧ aṠ We decompose the perturbation into its spatial and temporal parts: The evolution equation is then given as follows: In the E-dS Universe model, the special solution for the temporal part is given as follows: The equation in the transverse mode is given as follows: The evolution equation for the second-order perturbation coincides with that of the first-order perturbation. In other words, there are no perturbative solutions for the second-order transverse mode.
The transform between the determinant and the triplet terms is described in Appendix A. We decompose the perturbation into two spatial and two temporal parts: The evolution equations are given as follows: In the E-dS Universe model, the special solutions to the temporal parts are given as follows: (3 . 21) Sasaki and Kasai 19) derived third-order perturbative solutions to both the longitudinal and transverse modes. In this study, we only consider the effect of the firstorder longitudinal growing mode. The evolution equation for the third-order transverse mode is derived from Eq. (2 . 11). By using first-and second-order longitudinalmode equations, (3 . 1)-(3 . 3) and (3 . 9)-(3 . 11), respectively, the equation for the thirdorder transverse mode is given as follows: ,jl ψ ,kl − g 1 g 2 ε ijk ψ ,jl ψ ,kl = 0. (3 . 22) We decompose the perturbation into its spatial and temporal parts: The evolution equations are given as follows: Even if we consider only the longitudinal mode in the first-order perturbation, a third-order perturbation appears in the transverse mode. The meaning of the transverse mode will be examined in Sec. 5.
In the E-dS Universe model, the solution for the temporal part is given as follows: (3 . 26)

Fourth-order perturbative equations
Although fourth-order perturbative solutions have been derived in the past (for example, 41), 42) ), these solutions were reported before the publication of. 19) Therefore, these derivations ignored the effect of transverse modes. In this paper, we consider the transverse modes.
The fourth-order longitudinal mode is also affected by the third-order transverse mode, which is generated by the first-and second-order longitudinal modes.
The derivation of the fourth-order perturbative equations is presented in Appendix B. The fourth-order perturbative equations for the longitudinal mode can be derived. In this section, we decompose the perturbation into twelve different parts, six spatial and six temporal: The evolution equations are given as follows: In the E-dS Universe model, the special solutions to the temporal parts are given as follows: Next, we consider the evolution equation in the transverse mode. The evolution equation for the fourth-order transverse mode is derived from Eq. (2 . 11). We can decompose the perturbation into eight parts, four spatial and four temporal: (3 . 46) The fourth-order perturbative equations for the transverse mode are decomposed as follows:g In the E-dS Universe model, the special solutions to the temporal parts are given as follows: The fourth-order perturbative equations in this subsection are the main result of this paper. In E-dS Universe model, g 4 is proportional to a 4 .
For the generic Friedmann Universe model, we must solve the equations for the temporal parts. According to the behavior of g 2 and g 3 , g 4 would be proportional to a 4 mostly. In the matter-dominant era, the solutions almost coincide with those in E-dS Universe model. §4. Implication of fourth-order solutions for practical applications

Applications to resummation theory
Although the precision of the perturbation improves with increasing order, the effect of higher-order perturbation is quite small. When the density fluctuation evolves in the nonlinear region, it is difficult to describe the evolution using perturbation theory. Recently, a set of remarkable methods that substantially improve the precision has been studied. These methods constitute the resummation theory. 26), 27), 29), 28) There are many methods of the resummation theory, and still new methods have been also developed. Hiramatsu and Taruya 43) have developed numerical algorithm to solve closure evolutions in both perturbative and nonperturbative regimes. For time evolution of matter power spectra, the calculation result in nonperturbative regime almost coincides to that from N -body simulations. Crocce, Scoccimarro and Bernardeau 44) have developed fast renormalized perturbative scheme. The feature of the scheme is to apply multi-point propagators. They can derive nonlinear power spectra quickly. Pietroni 45) proposed a different method. The method is similar to the familiar BBGKY hierarchy. 1) To solve group of evolution equations, we must introduce some kind of truncation. In this paper. as one assumption, a part of fourpoint function is truncated. In other words, four-point correlator is consist from two-point correlators. Of course, these works consider generic Friedmann Universe models. In these works, although the growth rate of n-th order perturbation is given by g n 1 , the power spectrum is calculated with high precision. As we mentioned in Sec. 2, the relation between the density fluctuation and the Lagrangian perturbation is nonlinear. In other words, Lagrangian linear perturbation includes a nonlinear effect in the density fluctuation (Eulerian perturbation). Therefore, if we choose Lagrangian perturbation as the model, a subset of the infinite series of diagrams in Eulerian perturbation would be included. 29) Furthermore, the resummation based on LPT can compute the power spectrum in redshift space. The mapping from real space and redshift space is given analytically.
Recently, Okamura et al. 33) proposed a generic formula for m-loop resummation via LPT. The formula is shown in Appendix C. Their formula can apply to generic Friedmann Universe models. They described the evolution of baryon acoustic oscillation (BAO). 46), 47) Additionally, they compared improved LPT to N -body simulations in real space. Their two-loop correction can predict the power spectrum with a higher accuracy than when using a one-loop correction. However, because one-loop correction is sufficiently accurate for describing the nonlinearity, the effect of the two-loop correction is shown to be less significant for two-point correlation functions.
Their generic formula is applied to not only the power spectrum in real space but also in redshift space. Furthermore, their formula is adopted in the generic Friedmann Universe model. In this paper, we show the apparent form of the perturbative solutions in a flat Universe model without a cosmological constant.
When we consider the transverse mode in third-order LPT (3LPT), at least the power spectrum in redshift space would change. The transverse mode in 3LPT does not affect the growth of the density fluctuation (2 . 7). However, the displacement whose direction is different from that of the longitudinal mode is given. Details of the difference between the directions of the longitudinal and transverse modes will be discussed in next subsection.
For application of the resummation theory, higher-order perturbative solutions are the point. For two-loop correction for the power spectrum, we require the perturbative solutions up to 3LPT. For more improvement of the power spectrum, we will consider three-loop correction. As we show in Appendix C, 5LPT is required for three-loop correction.
Making a comparison between one-or two-loop corrections and N -body simulations, one finds that two-loop LPT is valid for k ≤ 0.3h −1 Mpc −1 at z = 2 with 1% accuracy.
On the other hand, when we consider bispectrum, 4LPT solutions are enough to compute at one-loop order in RPT . 34), 35) The contribution of the one-loop correction is closed by 4LPT.
Given the effects of higher-order LPT solutions, RPT offers a more effective prediction method for future deep-space surveys.

Initial conditions for cosmological N -body simulations
Although Lagrangian nonlinear perturbations describe quasi-nonlinear evolution well, they cannot describe strongly nonlinear evolution. After the formation of caustics, the hydrodynamical description no longer has physical meaning. For strongly nonlinear evolution, cosmological N -body simulations have been used. With improvements to algorithms for the computation of gravitational interactions and ad-vancements in computer technology, large, detailed simulations can be carried out.
Although the growth of the density fluctuations starts in the recombination era (z ≃ 1000), N -body simulations are begun from a later time, because of the difficulty in performing the simulation with small fluctuations. In the standard scheme, the small fluctuations evolve following the perturbative approach until they reach the quasi-nonlinear regime (z ≃ 30). 48) In the perturbative approach, ZA has long been used. In ZA, although the perturbation remains in the linear regime, quasi-nonlinear density fluctuations can be described. However, ZA does not reproduce the higher-order statistics, because the acceleration is always parallel to the peculiar velocity. In the quasilinear regime, effects that cannot be expressed using only linear perturbations appear. If we set up the initial condition for use with a cosmological N -body simulation at the quasinonlinear stage, ZA seems insufficient for describing the initial condition.
To solve this problem, Crocce, Pueblas, and Scoccimarro 37) proposed an improvement in which different initial conditions are adopted. Essentially, the initial conditions are based on approximations that are valid up to second-order LPT (2LPT) and reproduce the exact value of the skewness in the weakly nonlinear region. They set up the initial conditions for ZA and 2LPT at z = 49. The simulations were run using the Gadget2 code 49) and contained N = 512 3 particles. The cosmological parameters used were Ω m = 0.27, Ω Λ = 0.73, Ω b = 0.046, h = 0.72, and σ 8 (z = 0) = 0.9. The box size was L box = 512, 1024h −1 [Mpc].
Crocce, Pueblas, and Scoccimarro showed the evolution of higher-order cumulants such as skewness and kurtosis. First, they computed these quantities for the reference runs (2LPT, z i = 49). Next, they compared the quantities for reference runs with those for several simulations performed using other initial conditions. The difference between the skewness and the kurtosis values for the reference runs and the simulations using ZA is several percent for 2 < R < 60h −1 [Mpc] at z = 0, 1, 2, 3. Simulations using ZA underestimated the power spectrum in the nonlinear regime by approximately 2%, 4%, and 8% at z = 0, 1, and 3, respectively.
Thereafter, the effect of 3LPT was verified by Tatekawa and Mizuno. The conditions they considered differed from those in the former study. The simulations were run using the P 3 M code 50) and contained N = 128 3 particles. The cosmological parameters were taken from the WMAP three-year results. 51) They considered only the longitudinal mode in 3LPT. Their initial condition was given by z ≃ 23, 33, 84. In their results, the difference of both the skewness and the kurtosis appeared between the simulations using ZA and 2LPT. Even if the initial condition is given at z ≃ 84, the difference appeared. The result is consistent with that by Crocce, Pueblas, and Scoccimarro. On the other hand, if the initial condition was given by z ≃ 33, the difference between both the skewness and the kurtosis for the simulations using 2LPT and 3LPT is less than 2% for R = 2h −1 [Mpc] at late times. According to their results, 2LPT is sufficiently accurate when several percentage points of precision are needed for the higher-order cumulants. In other words, the initial conditions given by 2LPT have the enough accuracy in simulations with precision of one percent.
In the future, when higher precision is needed, we should select the initial conditions for N -body simulations using 3LPT or 4LPT. The transverse mode of 3LPT was ignored in previous analyses. Before estimating the effect of 4LPT, it is necessary to estimate the effect of the transverse mode in 3LPT. Here we suppose simple perturbation in first-order LPT (1LPT): This perturbation gives the concentration of the matter. The displacement vector is shown in Fig. 1. The longitudinal mode in 2LPT and 3LPT promotes the concentration. The direction of the displacement in both 2LPT and the longitudinal mode in 3LPT is parallel to that in 1LPT. However, the displacement of the transverse mode in 3LPT differs from that of the others. This difference in the direction of the displacement would especially affect the statistical quantities in redshift space.
Here, we roughly estimate the effect of 4LPT. The difference between the higherorder cumulants in the cases of ZA and 2LPT was approximately 10%. In contrast, the difference between the higher-order cumulants in the cases of 2LPT and 3LPT was approximately 1%. When the order of the perturbation is raised, the difference is 1/10. This tendency suggests that when the desired precision is less than 1%, we should select the initial condition using 4LPT. §5. Summary We have derived Lagrangian fourth-order perturbative equations for a cosmological fluid. The equations consist of six scalar equations and four vector equations. The spatial components of the perturbation are described by Laplace's equation.
In this paper, although we assumed that the linear perturbation only had a longitudinal mode, a transverse mode appears in the third-order perturbation; we consider the effect of this transverse mode. The density fluctuation is described by Eq. (2 . 7). In 3LPT, even if the transverse mode is present in the third-order perturbation, the mode does not affect the density fluctuation, because the mode disappears in Eq. (2 . 14); this mode affects only the displacement. The simple case was shown in Sec. 4.2.
In 4LPT, the doublet term consists of ZA, and the third-order transverse mode appears. When the order of the approximation is raised to fourth order, it is necessary to consider the effect of the transverse mode on not only the displacement but also the density fluctuation.
The appearance of the transverse mode does not imply the occurrence of vorticity. In LPT, the vorticity is described exactly: 17), 18), 19) ω(q, t 0 ) refers to the initial vorticity. Therefore, if vorticity does not exist in the initial condition, it never appears. This result is consistent with Kelvin's circulation theorem.
Here, we note a problem in higher-order perturbation. If the primordial density fluctuation is negative, the fluctuation δ approaches −1. In Eulerian perturbation, the density fluctuation exceeds −1, i.e., negative mass density appears. However, in LPT, δ decreases gently and converges to −1. However, the first-order perturbation remains the best approximation for describing the late-time evolution of voids. In particular, if we stop expansion when an even order (second, fourth, sixth, . . . ) is reached, the perturbative solution describes the contraction of a void. 25), 52), 53) Therefore, when significantly advancing the time evolution, we should apply the fourth-order perturbation.
Resummation theory significantly improves the precision of the perturbation. Even if we consider fourth-order Lagrangian perturbation, we could still avoid the above problem. In particular, the Lagrangian perturbation easily computes physical quantities not only in real space but also in redshift space. In the future, the resummation of the Lagrangian perturbation up to the higher-order will contribute to the prediction methods used in precision cosmology.
Rampf and Buchert 34) have been derived 4LPT equations independently. Although their peculiar-system approach corresponds to our approach, their 4LPT equations are different from our equations. The origin of the difference are treatment of temporal components and spatial derivative. For temporal components, they introduced conformal time for the equations. Then they considered both growing and decaying modes. For the spatial derivative, they considered inverse matrix for the Jacobian matrix between Eulerian coordinate and Lagrangian coordinate. At present, especially in the spatial components, both equations are inconsistent. The meaning of the difference in the spatial terms and its effects should be confirmed in near future.
For more higher-order correction, 4LPT would not be enough. Enough higherorder correction to resummation, fifth-order LPT (5LPT) has been achieved. Based on the calculation of 4LPT, 5LPT equations will likely consist of at least 19 scalar and 18 vector equations with our approach. Although the derivation of 5LPT equations appears to be quite a challenging problem, useful predictions for future deep-space observation might be possible; therefore, it is worth wrestling with this problem. Determinant of displacement We consider rewriting the determinant of the Lagrangian displacement matrix with triplet terms.
First, we define D as follows:  where A, B, and C are scalar functions of q. In this paper, we consider fourth-order perturbation. Therefore, we must consider a determinant that consists of first-and second-order perturbations: The determinant is rewritten in triplet terms as ,jk S ,jk + ,jk S ,ki , (A . 5)

Appendix B Derivation of fourth-order perturbative equations
In this section, we provide a detailed description of the derivation of fourth-order perturbative equations. First, we present a derivation of the longitudinal mode. The left-hand-side terms of the evolution equation (2 . 10) for the fourth-order longitudinal mode are given as follows: ,ii − s i,j + 2ȧ aṡ ,ki + 2ȧ aṠ The right-hand-side terms of the evolution equation (2 . 10) for the fourth-order longitudinal mode are given as follows: ,ii + 2S (1) ,ii S ,jj + S ,ii S ,jj − ,ii S ,jj − S ,ij S ,jj − S ,ij s ,jj S (2) ,kk − S (1) ,jk S (2) ,jk + S ,jj S ,kk S ,ll − S ,kl S ,jj S ,kk S ,ll .

Appendix C Higher-order loop correction for power spectrum
The generic formula about the power spectrum with N -loop corrections in both real space and redshift space has been derived. 33) Here we suppose Gaussian initial fluctuation. In real space, the power spectrum is described as δ D is Dirac's delta function. C i 1 ···i 2n is a 2n−order cumulant of the displacement field in Fourier space Ψ : For derivation of the power spectrum, N -loop contribution to the power spectrum predicted by SPT is required. In N -loop correction of Lagrangian resumption, the series in the exponent is expanded up to O (P L ) N . The remaining factor is expanded up to O (P L ) N +1 . Following the formula (C . 1), three-loop correction for the power spectrum is derived. Here we extend correction terms.