Multi-Scale Perturbation Theory II: Solutions and Leading-Order Bispectrum in the $\Lambda$CDM Universe

Two-parameter perturbation theory (2PPT) is a framework designed to include the relativistic gravitational effects of small-scale nonlinear structures on the large-scale properties of the Universe. In this paper we use the 2PPT framework to calculate and study the bispectrum of matter in a spatially-flat $\Lambda$CDM cosmology. This is achieved by deploying Newtonian perturbation theory to model the gravitational fields of quasi-nonlinear structures, and then subsequently using them as source terms for the large-scale cosmological perturbations. We find that our approach reproduces some of the expected relativistic effects from second-order cosmological perturbation theory, but not all. This work therefore provides a first step in deploying a formalism that can simultaneously model the weak gravitational fields of both linear and nonlinear structures in a realistic model of the Universe.


Introduction
The Universe we inhabit contains nonlinear structures on small scales ( 100 Mpc) [1], and linear structures on large scales ( 100 Mpc) [2]. The gravitational physics in these two different regimes can both be understood by using weak-field expansions of Einstein's field equations, assuming one does not wish to resolve compact objects like neutron stars and black holes [3]. However, the precise nature of the approach that is required to accurately model each regime is not the same.
On small scales cosmologists routinely appeal to the post-Newtonian limit of Einstein's theory, as Newton's law of gravity is capable of modelling the gravitational fields produced by arbitrarily large density constrasts (assuming gravity remains weak, and the speed of matter fields remain non-relativistic) [4]. However, on large-scales the Newtonian approach cannot be applied, and cosmological perturbation theory must be used to model the gravitational fields of matter [5]. Cosmological perturbation theory is extremely versatile, but is only applicable in the limit where density constrasts and the peculiar velocity of matter fields are both small (of the order of magnitude as the gravitational potentials themselves, i.e. ∼ 10 −5 ).
The usual approach to dealing with the different mathematical expansions required in each of these two limits is to assume that each can be performed independently of the other, and that the results of both should overlap in some intermediary regime. While this is plausible at linear order, it is unlikely to be a viable approach at second and higher order in either of the two expansions. This is because second-(and higher-) order field equations are known to mix scales. That is, functions that have support only on small scales only can source large-scale fluctuations when combined quadratically, or in higher powers, with other perturbations.
The "2-Parameter Perturbation Theory" (2PPT) was recently introduced in order to provide a mathematical framework in which to understand and study this problem [6,7]. This approach simultaneously performs a (post-)Newtonian expansion and an expansion in cosmological perutrbation theory around the same Friedmann background. It makes explicit which features of each expansion can be understood as being independent of the other, as well as the ways in which the two expansions are linked and interact. In particular, this approach makes explicit how the Friedmann equations of the background cosmology are sourced by the average of the Newtonian masses that exist within it, and how the linear perturbation equations for large-scale perturbations can be sourced by quadratic and cubic combinations of the Newtonian quantities from the small-scale sector of the expansion [8].
In a previous paper we examined the 2PPT approach in the test case of a universe with vanishing cosmological constant (i.e. with Λ = 0), and during the epochs in which the gravitational effects of relativistic matter are negligible [9]. The leading-order equations for describing the small-scale physics in this case are identical to the usual equations used to describe Newtonian gravity on an expanding background [10], and were solved using Newtonian perturbation theory [11]. These solutions were then used to solve the equations for the leading-order part of the gravitational fields on large scales, order-by-order in the Newtonian perturbation theory expansion.
It was found that when Λ = 0 the leading-order contribution to the 2PPT equations precisely reproduces the expected equations from cosmological perturbation theory and Newtonian gravity, within their respective domains of applicability. The next-to-leading order contributions to the 2PPT equations, which describe the evolution of long wavelength linear perturbations on top of an FLRW background that has been allowed to develop short-scale Newtonian nonlinear inhomogeneities, were then investigated. These contributions provide terms that appear in the calculation of the bispectrum of matter fluctuations, and were shown to produce sizeable corrections to the Newtonian bispectrum on scales k 10 −2 Mpc [9]. They include most (though not all) of the terms that occur in second-order cosmological perurbation theory, and give the leading-order relativistic effect of small-scale nonlinearities on the matter bispectrum.
In this paper we generalize the study performed in Ref. [9] to Friedmann cosmologies with non-zero cosmological constant, Λ = 0. Such situations are of obvious significance for observational cosmology, and require separate study due to the extra complexity that a non-zero Λ introduces into the evolution of perturbations. The mathematical foundations of including Λ = 0 in the 2PPT approach have already been investigated in Ref. [7], where the gauge depenencies and field equations that result were thoroughly investigated. Here we will use these results to calculate the leading-order relativistic corrections to the matter bispectrum within the 2PPT formalism, and with a non-vanishing cosmological constant.
We find that the time-dependence that the presence of a non-zero cosmological constant induces into the first-order gravitational potentials significantly complicates the application of the 2PPT equations, and increases the disparity between this approach and the more traditional use of second-order cosmological perturbation theory. While the 2PPT equations faithfully include the effects of small-scale nonlinearities on the relativisitic bispectrum, they do not include the second-order effects of large-scale linear fluctuations, which are nonnegligible on very large scales. This demonstrates that the 2PPT approach may need to be supplemented by additional terms in order to reliably include all relativistic effects in the calculation of at least some cosmological observables on very large scales.
In Section 2 we explain the 2PPT formalism, as required for use in a ΛCDM universe.
In Sections 3 and 4 we then solve these equations using Newtonian perturbation theory. This is followed by a calculatation of the bispectra of gravitational potentials and dark matter in Section 5, and a conclusion in Section 6. We use Latin letters for spatial indices and dashes to refer to differentiation with respect to conformal time. Spatial derivatives are written ∂ i , and ∇ 2 refers to the Laplacian operator associated with these derivatives. We also choose to work in geometrized units G = c = 1, and in longitudinal gauge [12], throughout.
2 Two-parameter perturbation theory in ΛCDM Two-parameter perturbation theory is constructed by performing post-Newtonian and cosmological perturbation theory expansions around a single Friedmann background, which results in a line-element that can be written in longitudinal gauge as [7] In what follows, we will briefly outline the equations that must be obeyed in the presence of a non-zero Λ by the Newtonian potential U , the scale factor a, and the large-scale cosmological perturbations φ and ψ. For further details, and more thorough mathematical treatment, justification and introduction, the reader is referred to Ref. [7].

Small-scale Newtonian equations
The U in this expression corresponds to the Newtonian gravitational potential, which is the leading-order contribution to the line-element from the post-Newtonian expansion, and which satisfies the usual Newton-Poisson equation on an expanding background: whereρ = ρ N is the cosmological average of the Newtonian mass density (taken over a suitable domain), and δ N = ρ N /ρ N is the density contrast in this quantity. The evolution equations for δ N is given by the energy conservation equation, where v Ni is the 3-velocity of the Newtonian masses (assumed to be of order v 2 ∼ U ), which has its own evolution equation as follows: 4) where θ N = ∂ i v Ni is the corresponding velocity divergence. The relative size of the terms δ N , θ N and H are determined in this expansion by noting that time derivatives add an order of smallness ∼ v compared to spatial derivatives (see Ref. [6]). These equations are all identical to the usual expressions for Newtonian gravity and dynamics on an expanding background [10], but here are formally derived as the leadingorder parts of a relativistic post-Newtonian expansion on an expanding background [13], which is used to describe gravitational fields on small scales. The reader will note that there is no restriction implied on the magnitude of the quantity δ N in these equations, though the average of the density is found to be of orderρ ∼ U . This means that δ N can be used to describe highly nonlinear structures. It may also be noted that the cosmological constant, Λ, does not appear in these equations directly, and only influences U through the effect it has on the scale factor, a(τ ).

Dynamical equations for the background
The background dynamics in the 2PPT approach are given by equations that are obtained by averaging the Newtonian-level contributions to the field equations, which results in the effective Friedmann equations [7] where H = a /a, and where primes denote differentiation with respect to conformal time τ .
The quantityρ that appears in these equations is the average of the Newtonian mass density, which obeys the conservation equation [6] ρ + 3Hρ = 0 , where time derivatives have again been used to assign an extra order of smallness of order ∼ v, and where we have taken Λ ∼ v 2 .
The reader may note that the densityρ that appears in these equations is the direct result of integrating the nonlinear Newtonian mass density ρ N over large cosmological domains, and is not a separate component of the energy density (as it is in cosmological perturbation theory). Otherwise, these equations are identical to those of Friedmann cosmology containing dust and a cosmological constant, and correspondingly have the same solutions.

Large-scale cosmological perturbations
The gravitational potentials φ and ψ, which appear in (2.1), correspond to the gravitational fields generated by the low amplitude density contrasts δ 1 that appear on scales k −1 10 2 Mpc. Given the existence of the small-scale Newtonian perturbations on a Friedmann background, and assuming U ∼ φ ∼ ψ, gives the following evolution equation for these fields [7]: and the constraints and where v i is the 3-velocity associated with the large-scale perturbations δ (not to be confused with the nonlinear small-scale density contrast δ N ). These equations have been derived under the assumption that vector and tensor modes can be neglected, and exist at a higher order in the 2PPT perturbative expansion than the previous sets presented in Sections 2.1-2.2. They combine both large and small scale gravitational potentials, though the small-scale potentials U are already specified by the Newton-Poisson constraint equation (2.2). They can also be seen to reduce to the linearized scalar equations from cosmological perturbation theory [5], in the case where U vanishes and quadratic products of density constrasts and 3-velocities are neglected.
In the absense of their neglect, there can be seen to be terms that are quadratic and cubic in quantities that would normally be considered to be small in cosmological perturbation theory. These terms appear in these equations because they are derived from the post-Newtonian sector of the theory, and because they are required to not be small in order to recover the appropriate Newtonian limit. Consistent application of the exact same logic forces these terms to appear in Equations (2.8)-(2.11), where they can act as sources for the large-scale gravitational potentials φ and ψ.
While equations (2.8)-(2.11) are intended to give the leading-order part of the field equations for the large-scale gravitational potentials in 2PPT, it is instructive to be able to compare them to the equations that one obtains from performing an expansion of the field equations up to second order in cosmological perturbation theory. In this case, it can readily be seen that there are both extra terms in the equations above, when compared to the second-order cosmological perturbation theory equations, as well as missing terms. In particular, equation (2.8) does not contain the terms 1 2 (U ) 2 , 4HU U and 2a 2 ΛU 2 , and equation (2.10) does not contain the terms −2H 2 U 2 and − 1 2 (U ) 2 . This is because, in every case, these terms appear at order ∼ v 6 in 2PPT expansion, which is smaller than every other term in these equations (which are all ∼ v 4 ). This difference originates from the fact that time-derivatives add an order of smallness in this approach, as is required for self-consistency of the post-Newtonian limit of this theory (see Ref. [6] for details).
On the other hand, equations (2.8)-(2.11) contain terms that would not only normally appear at second-order in cosmological perturbation theory, but also a handfull of terms that would not appear until third order. These include the term 4πa 2ρ . The appearance of these terms at leading-order in the 2PPT approach demonstrates the fact that δ N is allowed to be non-perturbatively large in this approach, as is required for consistent inclusion of non-linearities in the density field. We refer the reader to Refs. [6,7] for a more detailed explanation of the derivation of these equations, as well as for mathematical and physical justification of the assumptions that go into them.
Here we will proceed to find solutions to equations (2.8)-(2.11) by using Newtonian perturbation theory to appropriately expand them, as we will now describe. Once we have found solutions, we will compare these to the corresponding equations from second-order cosmological perturbation theory, which we have verified are obtained precisely by re-introducing the terms discussed above into equations (2.8)-(2.11), and which can be found in Ref. [16].

Solutions of the two-parameter equations
The order in which the 2PPT equations from Section 2 need to be solved goes as follows: (i) Solve the effective Friedmann equations for the background evolution to find a(τ ); i.e. solve equations (2.5)-(2.7) from Section 2.2.
(ii) Use the scale factor a(τ ) from step (i) to solve the equations for Newtonian gravity on an expanding background; equations (2.2)-(2.4) from Section 2.1.
(iii) Use the solutions from steps (i) and (ii) to solve for the large-scale gravitational potentials φ and ψ; using equations (2.8)-(2.11) from Section 2.3.
The first step here is straightforward, with known exact solutions available for spatially flat FRW cosmologies with dust and Λ [14]. The second step is more complicated, and could be approached in a number of ways [10,11]. In this paper we use Newtonian perturbation theory to model the solutions to these equations, in order to get approximate solutions for step (ii). This will largely follow the treatments given in Refs. [9,11,15]. With solutions for the U , δ N and v Ni in hand, we will then complete step (iii) by solving equations (2.8)-(2.11) order-by-order in the Newtonian perturbation theory expansion. This will involve borrowing techniques from second-order cosmological perturbation theory, as explained in Ref. [9].
Each step in the method described above will require the specification of appropriate initial conditions, which we will discuss as we proceed. Once we have solutions on both large and small scales, we will proceed to calculate the matter bispectrum of perturbations in Section 5. This is the simplest statistic that can be constructed that will have non-trivial contributions from the extra term in equations (2.8)-(2.11). For the rest of this section, we will proceed to find solutions at first and second order in the Newtonian perturbation theory expansion, which we will refer to as the 'first' and 'second' approximations.

First approximation
The Newtonian perturbation theory expansion we wish to use decomposes δ N , θ N and U into linear, second-order, and higher-order parts as where the label in brackets is denoting the order of a quantity in this expansion. These expressions are then subsituted into Equations (2.3) and (2.4), which can be manipulated to obtain and The first of these two equations is solved by δ 0 (x), where D(a) is the standard growth factor and δ (1) 0 (x) is the spatial initial condition. The second equation is solved by U (1) = ϕ = g(τ )ϕ 0 (x), where g(τ ) is a time-dependent factor and ϕ 0 (x) is the initial condition for the gravitational potential itself. These equations are identical to the ones used in standard Newtonian perturbation theory at linear in a ΛCDM universe, and we will not dwell on them further here.
We can now find the first approximation to equations (2.8)-(2.11), using the linearized solutions to the Newtonian equations given above. This begins by decomposing the cosmological variables into the different orders that will be solved for one-by-one using the solutions from the Newtonian perturbation theory applied to the short-scale domain, and results in the following: where numbers in brackets again refer to the order of approximation in the Newtonian perturbation theory. These expression can be subsituted into Equations (2.8)-(2.11), which immediately result in The evolution equation for the one degree of freedom that remains in the scalar sector of the gravitational theory can then be found to be This equation is identical in form to Equation (3.5), and is therefore solved by function with the same time dependency, g(τ ). This fact leads us to consider ψ (1) + U (1) as one single variable, with support on all length scales, which we will henceforth label as ϕ. Similarly absorbing the parts of θ (1) and δ (1) that behave like the first-approximations to the corresponding Newtonian quantities into θ N and δ N gives us that only remaining non-Newtonian part of these quantities is given by This process of re-defining Newtonian quantities to absorb the long-wavelength parts of the corresponding large-scale perturbations is described in more detail in Equations (5.12)-(5.16) of Ref. [9], and proceeds in exactly the same way here. These results show that the first approximation to the 2PPT equations gives identical results to standard first-order cosmological perturbation theory in a ΛCDM universe. Let us now consider the second approximation, which gives more interesting results.

Second approximation
After applying Newtonian perturbation theory to equations (2.2)-(2.4), we can extract the following evolution equation at second approximation: is the growth rate of structure, and where g and ϕ 0 are as in the equations presented above. This equation is solved by the usual expression for second-order Newtonian solutions in a ΛCDM universe [16]: where we have used the Newtonian kernel To find the large-scale potentials φ and ψ at second approximation it is useful to isolate the combination ψ (2) − φ (2) , using the trace-free ij field equation. The cosmological constant does not affect this equation, so the result is identical to the corresponding equation in an Einstein-de Sitter universe [9]: where (3.18) and where P = P i i . Using this allows us to write the second-approximation to the large-scale equations (2.8)-(2.11) as where we have in addition made use of the results If one were to compare Equation (3.19) to the corresponding equation in standard secondorder cosmological perturbation theory [17], it would be immediately apparent that there are a number of differences. In particular the terms 8Hϕϕ , (ϕ ) 2 , and 4a 2 Λϕ 2 are all absent here. This difference is due to the two-parameter counting scheme, which at leading-order neglects the effects of terms that are quadratic in large-scale perturbations. Such a difference does not occur in dust-only Einstein-de Sitter cosmologies [9], as in that case the potentials are independent of time.
The solutions to Equation (3.19) can be written as where the b n satisfy The solution for b 2 is the same as in second-order cosmological perturbation theory [16], and is given by The solution for b 1 , on the other hand, is given by where I(a) and W(a) are the source and Wronskian functions This result can be seen to reduce to the Einstein-de Sitter limit as Λ → 0, by noting that I → 0 in this case. While this result is cumbersome to deal with analytically, it straightforward to evaluate numerically.

Initital conditions and kernels
Let us now turn to the question of how to set initial conditions at τ in , for the expressions derived in the previous section. For this we choose to proceed by using standard cosmological perturbation theory to model the growth of structure in the early universe (before nonlinear structures developed), and then mapping those perturbations onto our 2PPT variables as some crossover time, τ cross . We therefore take g in = g(τ cross ) . Now, the function g(z) is known to approach a constant value in matter-domination, which means we are free to select any moment of time as the crossover time, provided that it is (a), calculated through numerical integration. The two curves can be seen to deviate when Λ becomes significant.
well within the matter-dominated era. We choose a redshift of z = 19 for z cross , which is equivalent to a in = 0.05.
In the case of the quantity b 1 , we can note that we are free to add a homogeneous solution to the differential equation it must obey. We use this fact to ensure continuity of the metric at τ in by setting for comparison. The two quantities can be seen to differ at late times due to the Λ-induced time-dependence of the first-order solutions, and the presence of extra terms at second-order in cosmological perturbation theory that contain their time derivatives.

Gravitational potential kernel
In order to calculate the bispectrum of the gravitational potential, we write the solutions from the previous section as where These are the same as the corresponding quantities in second-order perturbation theory [16], apart from C 1 which is modified due to the structure of the 2PPT equations. Taking the Fourier transform of this expression we obtain which immediately allows us to write down an expression for the kernel Let us know consider the kernel for matter perturbations.

Dark matter kernel
Given all the constituent parts of our solution in Equation (3.21), we can proceed by using the following expression to determine δ (2) : (4.10) After some manipulation, we find that this equation implies where ) , (4.12) Dg in H 2 0 Ω m0 (1 + 10a NL ) + 10 3 Dg (4.14) (4.16) The functions J 1 , J 3 and J 5 differ from the corresponding quantities determined using secondorder cosmological perturbation theory, and are plotted in Figures 2, 3 and 4 for comparison.
It may be noted that the expected Einstein-de Sitter behaviour is recovered at high redshift, up to a mathematical normalization that is physically unimportant. The solution from Equation (4.11) can be expressed as N (q 2 )K 2 (q 1 , q 2 , τ ) , (4.17)  where the kernel K 2 (k 1 , k 2 ) is given by and where These expressions are identical to those from second-order cosmological perturbation theory [18], except that in this case the functions J n are given by Equations (4.12)-(4.16). These expressions enable the calculation of the correlation functions of the dark matter overdensity.

Calculation of the matter bispectrum
In this section, we will present the results of the leading order corrections to the bispectra of the metric potential ψ (2) + U (2) and the dark-matter overdensity δ (2) + δ N .

Gravitational potentials
The leading order contribution to the dimensionless bispectrum of the scalar metric potential ψ + U is calculated as follows: where P ψ (k) is the dimensional linear power spectrum of the gravitational potential, as predicted by CLASS [19][20][21][22]. The implication of this result is that we can write the dimensionless bispectrum B ψ+U in terms of the kernel K ψ+U where ∆ ψ = 2π 2 k 3 P ψ (k). We plot the dimensionless bispectrum of gravitational potentials at redshift z = 0 in various different types of triangular configurations in Figures 5, 6 and 7. These results are plotted together with the linear power spectrum ∆ ψ in order to give the reader a reference point for the magnitudes of the quantities involved.
These plots demonstrate how the 2PPT solution for the gravitational potential deviates from the Newtonian prediction at small values of k. This behaviour can be anticipated from the form of the integral kernel; Newtonian terms come with an extra k 2 and hence dominate at large k, whereas relativistic corrections are dimensionless, and hence are only visible at small values of k. The difference with cosmological perturbation theory solution is maximised in the equilateral configuration, whilst it is significantly reduced in the flattened case.
In the case of the squeezed limit shown in Figure 6, it is clear that the differences between the 2PPT solution and the second-order CPT solution is negligible even at scales of k ∼ 10 −4 . This is because the kernel is dominated by the terms proportional toq 1 ·q 2 in this limit, and these terms are the same in both the cosmological perturbation theory and 2PPT case. We have verified that differences between 2PPT and CPT become negligible at higher redshifts, even in the equilateral configuration, as in this case we move into the matterdominated era where the differences between the gravitational potentials in each approach become negligible.

Dark matter overdensity
The second approximation to the dark-matter overdensity leads to the following expression for the leading order bispectrum in our approach: This quantity is plotted at redshift z = 0 in Figures 8, 9 and 10, in the equilateral, squeezed and flattened configurations respectively. Figures 8 9 and 10 display more significant deviation from the results of second-order CPT than in the Einstein-de Sitter case examined in Ref. [9]. This is to be expected, as the differences in the field equations are more significant. In the Einstein-de Sitter case, time derivatives of the first-order solution vanished, meaning that the omission of time-dependent terms like ϕ 2 and a 2 H 2 ϕ 2 were not a factor. However, these terms are non-zero in ΛCDM. The fact that they are present in the second-order CPT equations, but not the 2PPT field equations, therefore leads to the deviations shown in the plots.

Discussion
We have studied the application of the 2-parameter perturbation theory (2PPT) approach to a ΛCDM universe. This has involved solving for the short-scale Newtonian gravitational potentials using the leading-order part of the 2PPT-expanded field equations, and then using these solutions to determine the effects of nonlinear structure on long-wavelength cosmological fluctuations. We have, in particular, calculated the dark matter bispectrum, and have compared the results to standard second-order cosmological perturbation theory.
The method we used to solve the small-scale equations is Newtonian perturbation theory. This method allows analytically tractable solutions to be found in the quasi-nonlinear regime, where the equations of Newtonian physics on an expanding background are expected to accurately describe the physics, but where perturbation theory is still expected to give a realiable approximation to the true behaviour of the system. This approach has the very considerable benefit of allowing solutions to be written down in a relatively straightforward fashion, but suffers from the drawback that it does not accurately model the gravitational fields of matter on scales where it is highly nonlinear (i.e. on scales 10 Mpc).
With solutions for the Newtonian part of the gravitational field in hand, we then solved for the long-wavelength cosmological perturbations. We found that the leading-order contributions to the field equation that governs the evolution of these perturbations, (2.8), is      similar in form to the corresponding equation in second-order cosmological perturbation theory, though not identical. This is in contrast to the case where Λ = 0, where no such difference exists in the evolution equation [9]. The 2PPT equations are missing terms that would have occurred if calculating the corresponding quantities in second-order cosmological perturbation theory, and the bispectra that result are consequently modified (as shown in Figs. 5 -10).
Our findings illustrate both the advantages and disadvantages of using the 2PPT approach, as compared with other approaches to modelling structure in the Universe. First of all, our results show that the effects of the quasi-nonlinear Newtonian regime on cosmological fluctuations do indeed faithfully reproduce some of the known effects that are included in second-order cosmological perturbation theory. However, our study also shows that there exist large-scale relativistic effects that are not accounted for in the leading-order 2PPT equations that we have studied here, and that the absense of the terms that correspond to these effects can have considerable consequences for the calculation of observables such as the matter bispectrum on large scales.
The occurence of deviations from second-order cosmological perturbation theory is expected. This is because second-order perturbation theory is a theory of second-order corrections to the Einstein field equations in a regime where all fluctuations are small with respect to the background. The 2PPT field equations instead describe the evolution of first-order perturbations on top of a background which itself contains nonlinear Newtonian structure on small scales. As such, 2PPT will necessarily not contain the same gravitational selfinteractions as a second-order theory. Such terms could be included by going to higher order within the two-parameter expansion, or by selectively including terms which can be shown to have significant consequences for a given observable. Our calculation of the bispectrum shows that this may be necessary for at least some observables.
Instead, the strength of the 2PPT approach is that it can include the consequences of nonlinear much more readily. One way of doing this, in the context of the present study, would be to go to higer-order approximations in the Newtonian perturbation theory. This would be significantly easier than going to higher than second-order in cosmological perturbation theory, and could well provide a way of probing deeper into the regime of nonlinear structure formation. A more direct route, though less analytically tractible, would be to use the numerical output of Newtonian cosmological n-body simulations as the solutions to the Newtonian sector of the 2PPT system. This would allow for the full nonlinear regime to be included in ways that would be simply impossible in cosmological perturbation theory, and could lead to significant departures from the predictions of second-order perturbation theory on all scales, not just in the non-perturbative regime. We leave the study of these ideas for future work.