{'fraction_non_alphanumeric': 0.041890585733616856, 'fraction_numerical': 0.015580900831239126, 'mean_word_length': 4.679108573937864, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Closure Theory for Non-linear Evolution of Cosmological Power Spectra

We apply a non-linear statistical method in turbulence to the cosmological perturbation theory and derive a closed set of evolution equations for matter power spectra. The resultant closure equations consistently recover the one-loop results of standard perturbation theory and beyond that, it is still capable of treating the non-linear evolution of matter power spectra. We find the exact integral expressions for the solutions of closure equations. These analytic expressions coincide with the renormalized one-loop results presented by Crocce&Scoccimarro (2006,2007). By constructing the non-linear propagator, we analytically evaluate the non-linear matter power spectra based on the first-order Born approximation of the integral expressions and compare it with those of the renormalized perturbation theory.


INTRODUCTION
Cosmology now enters the era of precision cosmology. With large data set from the precision measurements of the cosmic microwave background anisotropies as well as the matter density fluctuations in the large-scale structure, the standard cosmological model has been fully established (e.g., Spergel et al. 2006;Tegmark et al. 2006). The associated cosmological parameters are well-determined with an accuracy of less than 10% level. With improved sensitivity and higher precision of future observations, the modern picture of the universe will be further reinforced and one can even explore a tiny signature of new physics beyond the standard cosmological model.
In fact, several ambitious missions for galaxy redshift survey are planed in order to reveal the nature of dark energy (e.g., Albrecht et al. 2006;Peacock et al. 2006, and references therein). Among these, Wide-field Fiber-fed Multi-Object Spectrograph (WFMOS) may be one of the best facility capable of achieving the percent-level measurement of baryon acoustic oscillations (BAOs) imprinted in the matter power spectrum (Meiksin, White & Peacock 1999). The recent observations from Sloan Digital Sky Survey (SDSS) and Two-degree-Field Galaxy Redshift Survey (2dF GRS) showed that the characteristic scale of BAOs can be used as the cosmic standard ruler to determine the distanceredshift relation of high-redshift galaxies (Cole et al. 2005;Eisenstein et al. 2005;Hütsi 2006;Percival et al. 2007). Since the distance-redshift relation is sensitive to the cosmic expansion history, details of the accelerated expansion can be clarified from the accurate measurement of BAOs (Seo & Eisenstein 2003). With percent-level measurement of characteristic scale of BAOs, the determination of dark energy equation-of-state, parametrized by w ≡ P de /ρ de , will achieve a few percent accuracy, where P de and ρ de are the pressure and the energy density of dark energy, respectively.
On the other hand, pursuit of the nature of dark energy highlights various fundamental problems which are potentially crucial for the accurate determination of dark energy equation-of-state. For example, the observation of BAOs requires high-precision theoretical template for the matter power spectrum in the relevant wavenumber, k ∼ 0.1 − 0.3hMpc −1 . To achieve the required accuracy in the determination of w, several systematic effects must be incorporated into the theoretical predictions. Among known systematic effects, the non-linear gravitational clustering is one of the most fundamental building block in the theory of structure formation. The recent N -body simulations showed that the non-linear growth of matter distribution significantly alters the shape of the power spectrum and the acoustic signature of BAOs tends to be erased, where linear theory prediction of matter power spectrum is no longer valid (e.g., Meiksin, White & Peacock 1999;Seo & Eisenstein 2005). To tackle the issue, the perturbation theory for gravitational clustering has been revived and has been applied to the study of BAOs (e.g., Suto & Sasaki 1991;Makino, Sasaki & Suto 1992;Jain & Bertschinger 1994;Scoccimarro & Frieman 1996;Jeong & Komatsu 2006;Nishimichi et al. 2007). The inclusion of leading-order correction to the non-linear clustering effect somehow improves a performance and reproduces the N -body results very well (Jeong & Komatsu 2006). At lower redshifts z < 2, however, next-to-leading order effect becomes important and the theoretical prediction with leading-order correction is insufficient to reproduce the N -body simulations.
Going beyond the perturbation theory, existing theoretical tools dealing with the non-linear gravitational clustering are the N -body simulation and the fitting formula for matter power spectrum (e.g., Peacock & Dodds 1996;Smith et al. 2003), as well as the phenomenological approach based on the halo model (Cooray & Sheth 2002, for a review). Currently, however, none of the reliable methods to ensure the percent-level precision exist. While the N -body simulation has a potential to provide a high-precision prediction, at present, one cannot blindly trust the N -body results unless reliable and comparable counterpart is established and is fully reconciled with N -body results. In this respect, development of new analytical method beyond the perturbation theory is necessary and essential for the progress on precision cosmology.
In this paper, we present a non-linear statistical method to predict the time evolution of matter power spectrum. Very recently, there appear several works on the statistical treatment going beyond the perturbation theory (Valageas 2004;Crocce & Scoccimarro 2006a,b;McDonald 2007;Valageas 2007a;Matarrese & Pietroni 2007;Izumi & Soda 2007). Based on the field-theoretical approach, the perturbation theory has been reformulated by improving the summation of the naive perturbative expansion. The so-called renormalized perturbation theory (RPT) developed by Crocce & Scoccimarro (2006a) seems viable theoretical tool alternative to the N -body simulation, suited to a high-precision prediction. Using RPT, some attempts to predict the non-linear evolution of BAOs has been reported (Crocce & Scoccimarro 2007). Here, we consider an alternative statistical method accepted widely in the subject of statistical theory of turbulence (e.g., Kraichnan 1959;Leslie 1973;Kaneda 1981;Kida & Goto 1997). In contrast to the sophisticated treatment based on the field-theoretical approach, our method is rather primitive in the sense that the effect of non-linearity of matter power spectrum is simply described by a systematic expansion and a truncation of the naive perturbation. After applying the so-called reversed expansion, the perturbative expansion is effectively re-organized and a class of higher-order corrections is systematically re-summed (Wyld 1961). With this treatment, some non-perturbative effects are also incorporated. We derive a closed set of moment equations characterizing the non-linear evolution of power spectra. The resultant evolution equations consistently recover the leading-order results of standard perturbation theory (one-loop PT). The solutions for these closure equations have exact integral expressions, which coincides with the one-loop results of RPT. Constructing the non-linear propagator, we attempt to evaluate the non-linear matter power spectra analytically. Based on the first-order Born approximation of the exact integral expressions, we find that the power spectra from the closure theory reasonably agree with those from the RPT.
The outline of the paper is as follows. In Section 2, governing equations for the dynamics of cosmological gravitational clustering is presented on the basis of the fluid description of the Vlasov equation. The closure problem, i.e., theoretical issue on the self-consistent treatment of the hierarchy of moment equations, is briefly mentioned. In Section 3, nonlinear statistical method in turbulence called direct interaction approximation is introduced and is applied to the present cosmological situation. Then, a closed set of evolution equations for matter power spectrum is obtained by the systematic perturbative expansion and the so-called reversion procedure. In Section 4, important properties of the resultant closure equations, i.e., recovery of standard PT and exact integral solutions, are discussed, based on detailed mathematical calculations presented in Appendices A and B. In Section 5, analytical treatment for our closure system is presented for illustrative purpose. Constructing approximate solutions for non-linear propagator, power spectrum of density fluctuations is calculated based on the Born approximation of the integral solutions. The results are compared with those obtained from RPT, particularly focusing on the non-linear evolution of BAOs. Finally, Section 6 is devoted to discussion and conclusions.

Basic equations
Throughout the paper, we consider the evolution of mass distribution in the flat universe, neglecting the tiny contribution from the massive neutrinos. To evaluate the non-linear growth of the density perturbations, hydrodynamic description of mass distribution is useful. Strictly speaking, this treatment is not exact and is often called the single-stream approximation of the Vlasov equation. However, at least in a statistical sense, it would be the best approximation if the scale of our interest is sufficiently large so that one can safely ignore the effect of shell-crossing. Denoting the fluctuation of mass distribution consisting of the cold dark matter and the baryon fluid by δ, we have ∂δ ∂t Here, a is the scale factor of the universe, H is the Hubble parameter and ρ m is the homogeneous mass density field. Assuming the irrotationality of the fluid flow, the above equations can be recast as (e.g., Taruya 2000): where the quantity θ is the velocity divergence defined by The quantity Ω w is the density parameter of dark energy satisfying the equation of sate, P de = wρ de , defined by Ω w = 8π G ρ de /(3H 2 ). Note that the relation Ω m + Ω w = 1 holds in the flat cosmology.
To treat the non-linear evolution of the matter power spectrum, the Fourier representation of equations (4) and (5) is useful. To do this, we introduce the Fourier transform of the perturbed quantities: The assumption of irrotational flow implies Then, the fluid equations (4) and (5) can be written as (e.g., Taruya 2000): where the kernels in the Fourier integrals, α and β, are respectively given by For later analysis, it is convenient to introduce the vector-field notation and rewrite the equations (9) and (10) in more compact form. For this purpose, we define where the function f (t) is given by f (t) ≡ d ln D(t)/d ln a and the quantity D(t) being the linear growth rate. Then, in terms of the new time variable η ≡ ln D(t), the evolution equation for the vector quantity Φ a (k; t) becomes (e.g., Crocce & Scoccimarro 2006a;Valageas 2007a): where δ D is the Dirac delta function.
Here and in what follows, we use the summation convention that the repetition of the same subscripts indicates the sum over the whole vector components. The time-dependent matrix Ω ab (η) is given by Each component of the vertex function γ abc becomes Note that the vertex function γ abc has the following symmetric properties: γ abc (k 1 , k 2 ) = γ acb (k 2 , k 1 ), γ abc (−k 1 , −k 2 ) = γ abc (k 1 , k 2 ), γ abc (k 1 , −k 2 ) = γ abc (−k 1 , k 2 ) and γ abc (k, −k) = 0. Equation (13) with (14) (15) is the basic equation for our subsequent analysis.

Moment equations
Before addressing the non-linear statistical method, it is worthwhile to mention the closure problem for the dynamics of statistical quantities.
First of all, we define the two kinds of the power spectra for fluid field u i (k; t): where the bracket · stands for the ensemble average. The quantity P ab is the ordinary power spectra which we are interested in and we have the symmetry, P ab = P ba . On the other hand, the quantity R ab represents the cross power spectrum between different times, which will be later important when we derive a closed set of equations. Note that R ab = R ba , in general.
Since we are specifically concerned with the time evolution of statistical quantities P ab (k; η) and R ab (k; η, η ′ ), rather than focusing on equation (13), it seems convenient to treat the moment equations for these quantities. To derive the moment equation for P ab , we first note that With a help of equation (13), we eliminate the time derivative ∂Φ a /∂η. Then, we obtain Similarly, moment equations for R ab becomes Here, we have introduced the two kinds of operators, Σ abcd and Λ ac : Equations (18) and (19) are not yet closed because they contain the higher-order correlation functions (or bi-spectra). In order to obtain the closed set of evolution equations, it is necessary to derive the evolution equations for higher-order correlation functions. However, the repetition of this treatment produces an infinite number of evolution equations and one cannot obtain a closed set of equations. This is the so-called closure problem for dynamics of statistical quantities. Note that the closure problem considered here is very close to the concept of BBGKY hierarchy, but slightly different in some sense. The BBGKY hierarchy arises from the many-body system characterized by the Loiuville equations and it also appears in the linear system. On the other hand, the origin of the closure problem essentially comes from the non-linearity of equation (13). Hence, to derive a closed set of moment equations, one must devise to introduce some truncation procedures by approximately treating the non-linear interaction in a self-consistent manner. The self-consistent truncation procedure is referred to as the closure theory (or closure approximation) in the statistical theory of turbulence and various closure theories have been so far exploited. In what follows, we especially consider the so-called direct-interaction approximation as one of the reliable closure theories.
Before closing this subsection, we define the propagator G ab (k, η|k ′ , η ′ ), which will later play a key role in deriving a closed set of equations: where δ stands for a functional derivative. It represents the influence on Φ a (k; η) at time η due to an infinitesimal disturbance for Φ b (k ′ ; η ′ ) (η ≥ η ′ ). Taking a functional derivative of equation (13), we obtain the governing equation for the propagator as with the boundary condition:

CLOSURE THEORY
3.1. Direct-interaction approximation In the subject of fluid mechanics, statistical characterization of turbulence for incompressible fluid flows is one of the major goals to understand the non-linear dynamics of Navier-Stokes equations. Among various attempts o construct a statistical theory of turbulence, there are approaches on systematic renormalized expansion (e.g., Kraichnan 1959;Wyld 1961;Leslie 1973). Direct-interaction approximation (DIA) is one of the best-known approximations and provides a simple truncation procedure (Kraichnan 1977;Kaneda 1981;Kida & Goto 1997). DIA has several desirable properties such as the local energy conservation and the realizability of energy spectrum. Also, the agreement with numerical simulations of isotropic turbulence is excellent even at high Reynolds number. Although we are especially concerned with the dynamics of compressible and irrotational fluid flow, the non-linearity in the cosmic fluid system essentially come from the same advection terms as in the Navier-Stokes equations. In this respect, DIA is a promising method to give a quantitative prediction of matter power spectrum.
To derive a closed set of evolution equations in the DIA, one transparent and intuitive way is to decompose the true field Φ a into the direct-interaction (DI) and the non-direct-interaction (NDI) parts. Let us consider particular Fourier modes (k, p, q) and especially focus on the time evolution of Φ a for a specific Fourier mode k. Through the non-linear term in the right-hand side of equation (13), the time evolution of Φ a (k; η) is determined by the infinite sum of three Fourier modes. We can formally decompose the quantity Φ a (k; η) into In the above expression, the quantity Φ (DI) a denotes the DI field, whose time evolution is determined by the direct interaction with the particular Fourier modes, p and q. On the other hand, the quantity Φ (NDI) a is defined as a fictitious field without the direct interaction between three modes k, p and q.
In the system governed by dimensionless equation (13), the strength of the non-linearity is characterized by the number of interacting Fourier modes. In the non-linear regime, one naturally expects that NDI field Φ (NDI) a plays the most dominant part in the time evolution of the quantity Φ a (k; η). Hence, in the DIA, we treat the DI field as a small perturbed quantity, relying on the assumption Φ . In addition, we put the following assumptions: (i) Gaussianity of the NDI field; (ii) statistical independence among the modes without direct interaction; (iii) statistical independence between the NDI field and the propagator. These assumptions basically come from the physical intuition for fully developed turbulence. In the presence of the infinite sum of the quadratic interaction, the fluid fields are expected to be nearly Gaussian, as naively indicated from the central-limit theorem. Also, it seems plausible that the initial correlation between different modes or quantities tends to be lost and the fluid quantities become statistically independent along the course of the non-linear interaction. Then, relying on these assumptions, we systematically expand the moment equations and the governing equation for propagator. Evaluating the ensemble average by using the formal solution of DI field in terms of the NDI fields, the closed set of equations for NDI fields can be finally derived (Kraichnan 1959;Kida & Goto 1997;Goto & Kida 1998).
In the following, we shall derive a closed set of equations in an alternative route. As it has been shown by Goto & Kida (1998), the resultant closure equations by DIA is identical with those obtained from the so-called reversed expansion procedure by introducing the fictitious parameter 1 (see also Leslie 1973;Kraichnan 1977;Kaneda 1981). The reversed expansion procedure seems rather straightforward and the assumptions made in this procedure are relevant for the present cosmological situations.

Derivation
Let us first introduce the fictitious parameter λ, which represents the strength of the non-linearity, and consider the weakly non-linear regime. We rewrite the field Φ a and the propagator G ab with In terms of these, the basic equations (13) and (22) respectively become Our task is to derive the consistent closure equations from the moment equations (18) and (19) with a help of the factitious parameter λ. To do this, we regard λ as a small expansion parameter (i.e., λ ≪ 1) and put the following assumptions: Assumptions (i) The evolution of field Φ a (k; η) is started from a tiny fluctuation and the development of non-linearity is mild.
(ii) The field Φ a (k; η) is a homogeneous random field and the statistical property of Φ a (k; η) is approximately described by the Gaussian statistics.
(iii) At the leading order, the field Φ a and the propagator G bc are statistically independent from each other.
Then, based on the perturbative calculation, we first evaluate the higher-order correlation terms in the moment equations (18) and (19). Inverting the perturbative expansion by a formal replacement of the perturbed quantities and setting the fictitious parameter to λ = 1 at a final step, we obtain a closed set of evolution equations. The set of equations derived here may be regarded as the result of renormalization and/or resummation of the perturbative expansion, which will be applicable to the nonlinear regime of gravitational clustering beyond the standard perturbation theory.

Naive perturbation
Based on the assumption (i), let us evaluate the nonlinear terms by a perturbative expansion of the small parameter λ. To do this, we expand the quantities Φ a and G ab as: The boundary condition for propagator in each order is Also, the boundary condition for the first-order quantity Φ (1) Since the equations for zeroth-order quantities are source-free, no mode-mode coupling occurs. Thus, the zeroth-order propagator G ab satisfying the boundary condition (29) can be expressed in the following form: Note that the function G ab (k|η, η ′ ) is a dimensionless quantity. Using this functional form, we formally write down the first-order solutions to the quantities Φ (1) i and G (1) . From (26) and (27), we obtain for η > η ′ > η 0 . According to the assumption (ii), one may treat Φ (0) a as a Gaussian random variable. Then, applying the perturbative expansion (28), we calculate the lowest-order non-vanishing contribution to the three-point correlation in equations (18) and (19). Let us first deal with the right-hand side of equation (18). At the lowest-order contribution, the three-point correlation becomes which are of the order of O(λ 1 ). Substituting the formal solution (32) into the above, we obtain The above expression is further reduced if we use the perturbative expression for the power spectra (16): The leading terms P (0) ab (k; η, η ′ ). Then, after some algebra, equation (35) is finally reduced to the following form: The function F (2) apq is explicitly written in terms of the two-time correlation R (0) ab : In deriving equation (38), we have used the symmetric properties of the vertex function, i.e., γ abc (k 1 , k 2 ) = γ acb (k 2 , k 1 ) and γ abc (k, −k) = 0.
In similar manner, we perturbatively evaluate the three-point correlation in equation (19). The resultant expression becomes Here, the function Θ(t) is the Heaviside step function. Now, summing up the perturbative expressions of three-point correlations, the moment equations (18) and (19) respectively become

Non-linear propagator
Next, we evaluate the higher-order terms in the governing equation for propagator (27). To do this, we first notice that the propagator G ab is no longer deterministic. Because of the interaction with the random field Φ a , the time evolution of G ab also exhibits stochastic nature. We thus treat the equation (27) statistically: The higher-order term in the right-hand side of equation (43) is perturbatively evaluated as In the above equation, the first two terms at the right-hand side become vanishing because of the assumption (iii).
The only non-vanishing contribution comes from the last term, which can be recast as In the last equality, we have used equation (31) and the definition of R ab (see Eq. [36] below). Hence, the perturbative evaluation of equation (27) becomes up to the contribution of O(λ 2 ).

Reversed expansion
We are in position to employ the procedure of the so-called reversion to rewrite the perturbative expressions (41), (42) and (46).
Recall from the perturbative expansion (36) that the O(λ 4 ) and the higher-order terms of the power spectra can be expressed in terms of the P ab . Also, the ensemble average of the higher-order propagator is formally expressed in terms of the zeroth-order quantities, G ab and R (0) ab . We then regard the expansion (36) as equations for P ab and R ab , the solutions of which are written in powers of λ as Similarly, we may write This procedure is called the reversion and it corresponds to the resummation of the perturbation series. Thus, at the leading order, equations (41), (42) and (46) are written in terms of the true field variables P ab , R ab and G ab . Now, we do not necessarily treat λ as the small parameter. This is just a book-keeping parameter and we finally set it to unity. Dropping the tilde over the quantities P ab , R ab and G ab , we at last reach the closure equations: The explicit expressions for the kernels F apq and K cpq are summarized as F apq (k, k 1 , k 2 ; η) = 2 η η0 dη ′′ 2 G ql (k 2 |η, η ′′ ) γ lrs (k, k 1 ) R ar (k; η, η ′′ )R ps (k 1 ; η, η ′′ ) where we have used the fact that the functions F apq and K cpq always appear as the product of γ bpq F apq and γ apq K cpq .

PROPERTIES OF CLOSURE EQUATIONS
In this section, we discuss some important properties of the closure equations derived in previous section.
4.1. Recovery of one-loop perturbations The closure equations in previous section have been derived in somewhat non-trivial manner by the reversed expansion procedure. While we employ the perturbative approach when evaluating the moment equations, the final governing equations for matter power spectrum become the non-linear coupled system and it seems unclear whether the power spectra calculated from the closure equations consistently recover the results of the standard perturbation theory at some level. In this respect, the recovery of the perturbation calculation may be a fast important check for the usefulness of the closure approximation.
In the standard treatment of the perturbation theory, the field Φ a is assumed to be a small perturbed quantity and is expanded as Substituting the above expansion into the evolution equation (13), we systematically derive the perturbation equations and through the order-by-order treatment, the solutions for higher-order quantities are expressed in terms of the linearorder quantity Φ (1) a . Further assuming the Gaussianity of the linear-order quantity Φ (1) a , the power spectra can be summarized as P ab (k) = P where the first term in the right-hand-side of equation is the linear power spectra and the quantities in the curly bracket is the so-called one-loop corrections to the power spectra, given by ab (k; η), ab (k; η).
In Appendix A, we show that the linear plus one-loop power spectra satisfy the following evolution equations: where the function F apq exactly coincides with the definition (52), with the replacement of the non-linear propagator and power spectra, G ab and R ab , with those of the linear counterparts, g ab and R ab : and Thus, in the weakly non-linear regime, the closure approximation faithfully recover the one-loop results of the standard perturbation theory. This will be manifestly apparent in next section by calculating the power spectra. A great emphasis is that among several non-perturbative approaches, the closure equations as non-linear coupled system also have the potential to go beyond the perturbation theory. In fact, our closure theory is basically equivalent to the the 2PI effective action approach by Valageas (2007a) and the one-loop level of the RPT by Crocce & Scoccimarro (2006a). This will be explicitly shown in next subsection when we obtain the exact integral expressions.

Exact integral solutions
The closure equations derived before seem rather complicated and analytically intractable because of its non-linearity and non-locality. In practice, a sophisticated numerical treatment is required to get the exact solutions for closure system. Note, however, that the closure equations possess the exact integral expressions for the power spectra P ab and R ab , which are formal solutions of the closure equations: The above expressions contain the function Φ(k; η 1 , η 2 ), which represents the non-linear mode-coupling between different Fourier modes, given by Note that the mode-coupling function Φ possesses the following symmetry: Φ ab (k; η 1 , η 2 ) = Φ ba (k; η 2 , η 1 ). In Appendix B, the integral expressions (60) and (61) are indeed compatible with the closure equations if the mode-coupling function Φ is given by equation (62).

The integral expressions given above have been also derived based on the RPT by Crocce & Scoccimarro (2006a) and/or through the path-integral formulation by Matarrese & Pietroni (2007) (see also Valageas 2007a;
Crocce & Scoccimarro 2007), although their derivations are quite formal. In contrast to their formal expressions, our integral solutions have the explicit functional dependence of the mode-coupling function Φ ab on the power spectra R ab and the propagators G ab . In the language of the RPT, this corresponds to the renormalized expressions for the mode-coupling power up to the one-loop order 2 . In this respect, the closure equations is a non-perturbative description of the power spectra going beyond the perturbation theory, and have an ability to predict the matter power spectra accurately, the result of which will be comparable to the one-loop results from RPT or the path-integral approach. Note, however, that our closure system has time evolution of the non-linear propagator, whose governing equation has been also derived by the self-consistent truncation of the higher-order corrections. On the other hand, in the RPT, no such truncation is considered in deriving the integral expressions. This may cause a major difference in the prediction of matter power spectra, which we will address in detail in next section.

ANALYTICAL TREATMENT OF NON-LINEAR POWER SPECTRUM
In this section, we compute the power spectrum of mass density fluctuation from the closure equations. Based on the exact integral expressions, we employ the Born approximation to obtain the analytic expressions for power spectrum. Further, the approximate expressions for the non-linear propagator is obtained by matching the two asymptotic behaviors. Combining these results, we evaluate the power spectrum of mass fluctuations and compare it with the one obtained from the RPT (Crocce & Scoccimarro 2007). In what follows, the following cosmological parameters are adopted to compute the power spectra : Ω m,0 = 0.27, Ω b,0 = 0.043, Ω w,0 = 0.73, w = −1, h = 0.7, σ 8 = 0.8, and n s = 1.

Born approximation
In practice, numerical treatment to directly solve the closure equations would be essential for an accurate evaluation of non-linear power spectrum. Nevertheless, analytical evaluation of the power spectrum is instructive and very helpful to understand the behavior of the non-linear corrections incorporated into the closure system. Crocce & Scoccimarro (2007) recently applied the first-order Born approximation to the exact integral expressions for the power spectra and derived the analytic expressions correctly up to the two-loop order. In the following, adopting the same approximation, we will analytically evaluate the power spectrum.
Here, the term, Born approximation, means the iterative approximation scheme to evaluate the integral equations like (60) and (61). Applying the Born approximation to the integral solutions, the power spectra P ab and R ab are first evaluated by substituting the linear-order quantities into the right-hand side of the expressions (60) and (61). This is the first-order Born approximation and it can be further improved by repeating the iterative substitution of the leading-order solutions into the right-hand side of the integral solutions. This treatment can be accurate in principle and the higher-order corrections become negligible as long as the contribution from the mode-coupling function is small, compared to the first term in right-hand side of integral equations.
In next subsection, we will present the approximate expression for the non-linear propagator, G approx ab . Here, assuming the analytic form of G ab , we derive analytical expressions for the quantity P ab (k; η) by substituting the iterative solutions of the different-time power spectra R ab (k; η, η ′ ). Let us denote the linear power spectra given at initial time by P lin ab (k; η 0 ). For a sufficiently small value of η 0 , the late-time evolution is dominated by the growing-mode solution. We thus put P lin ab (k; η 0 ) = e 2η0 P 0 (k) with the quantity P 0 (k) being the linearly extrapolated spectrum at the present time. Then, from equation (61), different-time power spectra are iteratively evaluated as follows: with G a ≡ G a1 + G a2 . Substituting the iterative solutions of R ab into the integral expression (60), one obtains Here, the source functions I a and J a are respectively given by Compared to the results given by Crocce & Scoccimarro (2007), we find that the first and second terms in equation (64) (2007), but the factor 2 is different. This discrepancy may be caused by the fact that the term P (III) given above is associated with the higher-order contributions to the Born approximation of the mode-coupling function (62), while the expression for P 2loop MC has been properly derived as the leading-order contribution to the mode-coupling term at two-loop order. For the influence of the higher-order corrections, we leave the discussions in section 5.3.

Approximate solution for non-linear propagator
Having obtained the analytic expressions for power spectra, our remaining task is to get the approximate solution for the non-linear propagator, G approx ab , from the closure equation (51). Here, we just follow the procedure suggested by Crocce & Scoccimarro (2006b) and construct the approximate solutions restricting their validity to the low-k or the high-k regions. Matching these asymptotic solutions appropriately at an intermediate regime, we obtain the global solutions for the non-linear propagator.

Solutions at one-loop order
Let us first consider the low-k limit of the non-linear propagator, where the perturbative treatment is safely applied. As it has been shown in section 4.1, our closure system consistently reproduces the one-loop results of the perturbation theory for power spectra P ab . Indeed, this is also true for the propagator G ab . Here, we explicitly write down the perturbative solutions at the one-loop level, δG 1-loop ab (k|η, η ′ ). Using the linear propagator g ab , perturbative solution of equation (51) is generally expressed as with the quantity R lin ab being the linear-order solution for different-time power spectra. In the approximation that the time-dependent matrix Ω ab (η) given by equation (14) is replaced with the constant matrix with Ω w = 0 and f = 1, the analytical solution for linear propagator g ab satisfying the evolution equation (58) is obtained and is given by (e.g., Crocce & Scoccimarro 2006a;Valageas 2007a): with Θ(x) being the Heaviside step function. For the different-time spectrum R lin ab , we have R lin ab (k; η 1 , η 2 ) = e η1+η2 1 1 where we have neglected the contribution from the decaying mode. Substituting the quantities (66) and (67) into the next-to-leading order solution (65), we first perform the time integrals over η 1 and η 2 . As for the three-dimensional Fourier integral, we can write d 3 q as 2π dq q 2 dx, where we have performed the integral over azimuthal angle. The variable x is the cosine of the angle between k and q, i.e., x = k · q/(kq), and the integral over x is performed analytically. A straightforward but lengthy calculation leads to (Crocce & Scoccimarro 2006b): where the matrices X ab and Y ab are given by The above expressions contain time-dependent functions α, β g,d , γ g,d and δ and scale-dependent functions f , g, h and i, whose explicit expressions are summarized in Appendix C. Note that the large-k limit of the above functions satisfies with the quantity σ v being the velocity dispersion for linear fluctuation defined by

Solutions in the high-k limit
Turn next to consider the high-k limit of the non-linear propagator. In the evolution equation (51), we take the limit k → +∞, while keeping q finite. In this limit, the vertex functions behave like To estimate the leading-order behavior analytically, the different-time spectrum R pr (q; η 1 , η 2 ) in equation (51) is also treated approximately by replacing it with the linear-order quantity R lin pr , given by equation (67). Then, the governing equation for non-linear propagator (51) is greatly simplified and we obtain where the quantity σ v is the rms fluctuation of the linear velocity given by equation (72). To solve the above equation, we adopt the following ansatz: with the initial condition, f (k|0) = 1. Using the basic property of the linear propagator, g ab (η, η ′ )g bc (η ′ , η ′′ ) = g ab (η, η ′′ ), equation (73) is rewritten with Here, for convenience, we introduced the new time variable τ = e η − e η ′ . The above equation has the analytical solution (Valageas 2007a). Writing the Laplace transform of the function f as f (k|s) = +∞ 0 ds e −sτ f (k|τ ), we have The solution of this equation satisfying the limit, f → 0 for s → +∞, becomes The inverse Laplace transform of the above expression is well-known and can be read off from the mathematical table: with the function J 1 (x) being a Bessel function of the first kind. Hence, the non-linear propagator in the high-k limit finally becomes

Matching the two solutions
The asymptotic behaviors obtained in section 5.2.1 and 5.2.2 have an overlapping region in which both of the approximations are applied. Therefore, matching these two solutions, one can obtain a global solution which would be a good approximation for the full propagator G. Let us recall from equation (71) that the propagator including the one-loop correction has the following asymptotic form: where we have only considered the dominant terms at η → ∞. On the other hand, the propagator in the high-k limit, (76), is perturbatively expanded as Comparing these two expressions, the approximate solution smoothly matching these asymptotic behaviors at η → ∞ may be Here, the matrices P ab and Q ab are defined as with X ab ≡ |2X ab | 1/2 and Y ab ≡ |2Y ab | 1/2 (see Eqs. [69] and [70] for definitions of X ab and Y ab ). In the weakly nonlinear regime, the propagator G approx ab correctly reproduces the one-loop results. In the large-k limit, the function (77) asymptotically approaches the solution (76).
Note that the approximate propagator (77) is derived in the same way as done in RPT of Crocce & Scoccimarro (2006b), although the functional dependence is somewhat different because of the different high-k behavior. In RPT, the matrices P ab and Q ab defined above should be replaced with which lead to the asymptotic behavior, G (approx) ab → g ab exp(−x 2 /2), in the high-k limit. Figure 1 shows the propagators G 1 = G 11 +G 12 (left) and G 2 = G 21 +G 22 (right) multiplying the factor D(z int )/D(z), for specific redshifts z = 0, 2 and 5, and with initial redshift z init = 35. The solid lines represent the approximate solutions obtained by matching the two asymptotic solutions. While the dotted lines show the results from the oneloop perturbation, the dashed lines indicate the non-linear propagators obtained from the RPT (Crocce & Scoccimarro 2006b). As increasing the wave number k, all the results exhibit the decaying behavior and the characteristic scale of the decay is shifted to low-k as decreasing the redshift. A closer look at small scale (high-k) reveals that the oneloop propagators, g 1, 2 + δ G 1-loop 1, 2 , show unphysical behavior, which eventually become negative and tend to diverge. On the other hand, the approximate solutions G approx 1, 2 obtained from the closure theory and RPT asymptotically approach zero as k → ∞. These damping behaviors are regarded as the non-perturbative effect as a result of the renormalization and/or self-consistent closure, which effectively takes account of the infinite series of higher-order corrections. Nevertheless, there exist some differences in the damping behaviors of the propagators. While the propagators in the RPT exhibit the exponential damping, the approximate solutions in closure theory show a damping oscillation. These differences may affect the final result of power spectrum. This point will be carefully discussed in next subsection.

Results and discussion
Having provided the basic ingredients for calculating the power spectrum, we now present the analytic results for power spectrum of density fluctuations, i.e., P 11 (k; z), and compare those with the results obtained from the RPT, particularly focusing on the characteristic scale of the BAOs. Figure 2 illustrates the overall behaviors of non-linear power spectrum of density fluctuations P (k; z) ≡ P 11 (k; z) given at z = 1, based on the Born approximation (64). Here, the contributions to the total power spectrum up to the first-order Born approximation, i.e., P (I) (k) and P (II) (k), are separately plotted. Thin and thick lines represent the results from RPT and the closure theory, respectively. The result from RPT is basically the same one as presented by Crocce & Scoccimarro (2007), although they further considered the higher-order contribution coming from the two-loop correction. Due to the damping behavior in the non-linear propagators, each contribution to the total power spectrum rapidly falls off in both predictions and their amplitudes become significantly lower than the linear power spectrum on small scales (labeled by P linear (k)), where the differences between the two predictions become manifest. Turning to focus on the scales larger than the damping scale, the contribution coming from P (II) becomes maximum around k ∼ 0.2h Mpc −1 , where the non-linear enhancement of the power spectrum P (k) can be seen and the differences between closure theory and RPT become fairly small. In particular, for the scale of our interest on BAOs around k < ∼ 0.01 − 0.3h Mpc −1 , one cannot clearly distinguish between both predictions from Figure 2.
To enlarge the differences between the two predictions and to clarify the non-linear behaviors of the BAOs, in Figure  3, snapshots of the power spectra divided by the smooth linear spectrum, P (k)/P no-wiggle (k) are plotted for specific redshifts z = 3, 2, 1 and 0.5 together with the N -body results kindly provided by Jeong & Komatsu (2006) (except for z=0.5), while in Figure 4, we present the logarithmic derivative of the power spectra, d ln P (k)/d ln k. All the results are plotted in linear scales. The smooth power spectra, P no-wiggle (k), was calculated from the linear transfer function without BAOs based on the fitting formula of Eisenstein & Hu (1998) (see Eq.[29] in their paper). Note that the power spectra calculated from the closure theory and RPT are the sum of the leading-order contributions, P (I) (k) and P (II) (k), not including the higher-order term P (III) (k). For comparison, the one-loop predictions from the standard perturbation theory are also depicted as dashed lines.
On large scales (low-k), the predictions both from the closure theory and RPT reasonably match the one-loop results of standard perturbation theory, as anticipated. This is just the quantitative check for the non-perturbative methods discussed in section 4.1. On the other hand, on smaller scales (high-k), the deviations from the one-loop perturbation become manifest and the amplitude of the predictions both from closure theory and RPT is suppressed compared to the one-loop predictions. The reduction of the power spectrum amplitude is the natural outcome of the damping behaviors appearing in the non-linear propagator and it qualitatively explains the behaviors seen in the N-body simulations (e.g., Jeong & Komatsu 2006;Matarrese & Pietroni 2007;Crocce & Scoccimarro 2007). However, at lower redshifts z = 1 and z = 0.5, the suppression of the amplitude is so significant that the predictions eventually become lower than the [64]). The contributions to the total power spectrum are separately plotted as indicated by P (I) (k) and P (II) (k) in the panel and the total power spectrum, P (I) (k) + P (II) (k), is depicted as the dashed lines. Note that in evaluating the power spectrum, the approximate solutions for the non-linear propagators G approx ab were used. Thick and thin lines indicate the results using the approximate solutions G approx ab from the closure theory and RPT, respectively. linear theory prediction. Compared to the prediction from RPT, the suppression of the amplitude is even larger for the prediction from the closure theory.
The strong suppression seen in the low redshifts seems somewhat unphysical, indicating the failure of our present analysis. In the panels with redshift z = 1 and z = 0.5 of Figures 3 and 4, the maximum wave number for the limitation of one-loop perturbation theory is estimated according to the condition, ∆ 2 (k) ≡ k 3 P (k)/(2π 2 ) < ∼ 0.4, given by Jeong & Komatsu (2006) and is depicted as vertical arrows, above which the perturbation prediction fails to reproduce the N-body results within the 1% accuracy. Comparaed to those critical scales, the validity range of the non-perturbative predictions seems rather narrow. While this statement seems rather controversial, our present analysis actually relies on several approximate treatments. Hence, the inclusion of the higher-order terms such as P (III) (k) might be essential for the correct non-linear behaviors. Crocce & Scoccimarro (2007) reported that the contributions higher than the two-loop correction become important on small scales and the prediction including the two-loop correction does improve the prediction, which reasonably agrees with N-body results within an accuracy of the N-body calculation. As we noted previously, however, the two-loop correction evaluated by the first-order Born approximation is comparable to the higher-order Born approximation coming from the P (III) (k) term, and at a level of the present analysis, it is not clear whether the higher-loop corrections rather than the higher-order Born approximation of the one-loop correction are essential or not. Rather, one obvious thing is that a further improvement of the approximate treatment is necessary to faithfully predict the non-linear evolution of BAOs. In the line of this, self-consistent treatment to solve the evolution equations (49), (50) and (51), would be one plausible approach, which we will address this issue in next task.

CONCLUSIONS
In this paper, motivated by a forthcoming experiment on the precision measurement of the BAOs imprinted on the matter power spectrum, a new theoretical tool for predicting the non-linear gravitational evolution of power spectrum was presented. In particular, we have applied the non-linear statistical method in turbulence to the cosmological perturbation theory and derived a closed set of matter power spectrum. The resultant evolution equations (49)-(51) with Fourier kernels (52) and (53) are the non-linear coupled system of integro-differential equations and these consistently recover the one-loop results of standard perturbation theory. Further, the exact integral expressions for the solutions of closure equations were obtained (see Eqs. [60] and [61] with mode-coupling function [62]), whose analytical expressions coincide with the renormalized one-loop results of the theory developed by Crocce & Scoccimarro (2006a), apart from the corrections coming from the vertex renormalization.
Based on the Born approximation to the exact integral expressions, we next tried to evaluate the non-linear power spectrum analytically. Constructing an approximate solution for non-linear propagator, which smoothly matches the two asymptotic solutions valid at low-k and high-k, power spectrum of density fluctuations was computed and the results are compared with those obtained from the RPT. Due to the non-linear damping behavior of the non- Fig. 3.-Ratio of non-linear power spectrum to smoothed linear spectrum, P (k)/P no-wiggle (k), given at specific redshifts, z = 3, 2, 1 and 0.5. The error bar represents the N-body results taken from Jeong & Komatsu (2006), in which different color indicates the results with different box size (see their paper in detail). Here, smoothed linear spectra P no-wiggle (k) were calculated from the linear transfer function without baryon acoustic oscillation according to the fitting formula of Eisenstein & Hu (1998) (Eq.[29] of their paper). The non-linear power spectra are obtained from the first-order Born approximation to the integral solution (Eq.[64]), with approximate solutions of the non-linear propagator given by closure theory (thick) and RPT (thin). For comparison, one-loop predictions from the standard perturbation theory are plotted in dashed lines. Also, in panels with z = 1 and 0.5, maximum wave number for limitation of one-loop perturbation is indicated by vertical arrows, according to the criterion, ∆ 2 (k) ≡ k 3 P (k)/(2π 2 ) < ∼ 0.4 (Jeong & Komatsu 2006).
linear propagator, the resultant amplitude of the power spectra up to the first-order Born approximation is strongly suppressed on small scales and this effect becomes significant as decreasing redshift. As a result, predicted spectra at high-k both from the closure theory and RPT fail to reproduce the N-body trends and the inclusion of the higherorder corrections is required. Nevertheless, at the intermediate scales, where the damping behavior of the non-linear propagator is rather mild, the closure theory and RPT both predict the deviation from the one-loop results of standard perturbation theory, which qualitatively agree with the N-body results. Although the analytical treatment presented here is still primitive and the range of the applicability is severely restricted by the validity of the approximations, the results indicate that our closure theory is a promising nonperturbative approach comparable to the RPT and a more elaborated study will provide an accurate prediction for non-linear power spectra going beyond standard perturbation theory. In this respect, a direct numerical treatment of the closure equations (49)-(51) is our urgent task. As it has been reported by Valageas (2007a), thanks to the full numerical treatment, the evolved result of the power spectrum shows several desirable properties and it qualitatively matches the N-body trends even at high-k. In future publication, these points will be further investigated from a quantitative point-of-view, particularly focusing on the non-linear evolution of BAOs.
Finally, one criticism on the present approach including the currently existing non-perturbative methods is that the present methodology heavily relies on the single-stream approximation of the Vlasov equation. In a strongly non-linear regime at the high-k region, the shell-crossing eventually occurs and the fluid description would be broken. A preliminary investigation suggests that the break down of the single-stream approximation may arise at k ∼ 1 − 2 hMpc −1 (Scoccimarro 2000; Afshordi 2007). Therefore, for the high-k region beyond the critical scale, the present approach cannot be applied and a more delicate treatment based on the Vlasov equation must be developed. This issue would be particularly crucial for accurate theoretical predictions to the cosmic shear statistics.
We are grateful to Keisuke Izumi, Shun Saito, Misao Sasaki, Jiro Soda and Yasushi Suto for useful discussions. We also thank Eiichiro Komatsu for providing the N-body data of Jeong & Komatsu (2006). This work is supported in Fig. 4.-Logarithmic derivative of non-linear power spectrum, d ln P (k)/d ln k, at redshifts z = 3, 2, 1 and 0.5. Line types and labels are the same as in Fig.3. part by a Grant-in-Aid for Scientific Research from the Japan Society for Promotion of Science (No. 18740132).
Now, let us consider the time evolution of power spectrum, P (22) (k; η). Acting the operator Σ abcd (η) defined by

INTEGRAL SOLUTIONS
In this appendix, we show that the integral expressions given in Section 4.2 are compatible with the closure equations (49) and (50). Here, we shall particularly focus on the integral expression (61) and explicitly derive the closure equation (50) from equation (61). As for the integral expression (60), it is straightforward to show the compatibility between equations (49) and (60), just repeating the same procedure as presented below.