{'fraction_non_alphanumeric': 0.03383671010017323, 'fraction_numerical': 0.008623936130149883, 'mean_word_length': 4.697779208239459, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Perturbation theory for cosmologies with non-linear structure

The next generation of cosmological surveys will operate over unprecedented scales, and will therefore provide exciting new opportunities for testing general relativity. The standard method for modelling the structures that these surveys will observe is to use cosmological perturbation theory for linear structures on horizon-sized scales, and Newtonian gravity for non-linear structures on much smaller scales. We propose a two-parameter formalism that generalizes this approach, thereby allowing interactions between large and small scales to be studied in a self-consistent and well-defined way. This uses both post-Newtonian gravity and cosmological perturbation theory, and can be used to model realistic cosmological scenarios including matter, radiation and a cosmological constant. We find that the resulting field equations can be written as a hierarchical set of perturbation equations. At leading-order, these equations allow us to recover a standard set of Friedmann equations, as well as a Newton-Poisson equation for the inhomogeneous part of the Newtonian energy density in an expanding background. For the perturbations in the large-scale cosmology, however, we find that the field equations are sourced by both non-linear and mode-mixing terms, due to the existence of small-scale structures. These extra terms should be expected to give rise to new gravitational effects, through the mixing of gravitational modes on small and large scales - effects that are beyond the scope of standard linear cosmological perturbation theory. We expect our formalism to be useful for accurately modelling gravitational physics in universes that contain non-linear structures, and for investigating the effects of non-linear gravity in the era of ultra-large-scale surveys.


I. INTRODUCTION
There is hope that the next generation of astronomical surveys [1-3], which will collect data on scales comparable to the cosmological horizon, will have sufficient precision to provide a new testing ground for non-linear relativistic gravity [4][5][6][7]. This is a particularly exciting prospect as, to date, non-linear gravitational effects have only been observed in the solar system [8][9][10][11], binary pulsar systems [12,13], and the newly discovered binary black hole mergers observed using LIGO [14][15][16]. The observation of corresponding effects in cosmology would allow general relativity to be investigated on entirely new length and time scales, as well as in a totally different physical environment. There is no telling what this might reveal about Einstein's theory.
To develop a mathematical formalism for investigating the non-linear properties of gravity in cosmology is, however, a highly non-trivial task. There is now a substantial literature dedicated to developing different approaches to modelling non-linear gravitational physics in cosmology. The most common approach is a direct implementation of second-order cosmological perturbation theory [17][18][19], which allows relativistic gravitational perturbations around a homogeneous and isotropic background to be modelled in the presence of linear density contrasts. Other approaches, however, have started to import techniques from post-Newtonian gravity [20][21][22], where gravitational fields are assumed to be slowly varying and where non-linear density contrasts can be consistently modelled.
These two different approaches are the obvious non-linear extensions of the standard method that has been so successful when considering linear gravity in cosmology; using cosmological perturbation theory to model linear structures on scales above the homogeneity scale ( 100Mpc), and Newtonian gravity to model non-linear structures on smaller scales. This looks very natural at linear order in the gravitational fields, partly because the equations of Newtonian gravity can be recovered as the quasi-static limit of cosmological perturbation theory (when the gravitational fields slowly vary in time). Indeed, this pleasing feature extends to non-linear gravitational fields in cosmology [20].
However, if one wants to consider non-linear gravity in a universe that simultaneously contains linear structures on large scales and non-linear structures on small scales, then one must adopt a more sophisticated approach. This is because, when considering the quasi-static limit, the terms that have been relegated to higher-order can no longer be entirely forgotten; they can and should be expected to appear in the next-to-leading-order gravitational field equations. This could be at second-order on small scales, but could in principle be at what is usually thought of as first-order on large scales. The question of how to construct a perturbative expansion that can systematically perform the required re-ordering, and produce a self-consistent and well-motivated set of field equations on all scales, is the purpose of this paper.
This work is a development of the two-parameter perturbative expansion proposed and investigated in a previous paper [23]. We extend and apply this work to the case of realistic cosmological models that contain relativistic fluids with barotropic equations of state, as well as a cosmological constant, Λ, and non-relativistic dust-like matter that can be used to model dark matter and baryons. The result is a set of equations that can be used to calculate the effect of small-scale structure on the leading-order perturbations on large scales. These equations contain terms that are quadratic in short-scale potentials and can be written as an effective fluid [24], as well as mode-mixing terms 1 that couple scalar, vector and tensor perturbations in the large-scale cosmology. Both of these two types of terms offer exciting possibilities for testing non-linear gravity with upcoming surveys.
In Section II we give the final field equations for the two-parameter perturbative expansion constructed during this paper, written in their effective fluid form. In Section III we discuss the construction of the two-parameter perturbation expansion. Section IV contains a study of infinitesimal gauge transformations in the context of twoparameter expansions, as well as the construction of gauge-invariant variables. The application of these variables to the field equations is then discussed in Section V. We use Greek indices to represent spacetime coordinates, and Latin indices for spatial coordinates. Dots refer to differentiation with respect to cosmic timeφ ≡ ∂ϕ ∂t , and dashes refer to differentiation with respect to conformal time ϕ ′ ≡ ∂ϕ ∂τ . Single spatial derivatives are given by ϕ ,i ≡ ∂ i ϕ ≡ ∇ϕ, and ∇ 2 refers to the Laplacian associated with spatial partial derivatives with respect to the comoving coordinates. Throughout, we use units in which G = c = 1.

II. PERTURBED FIELD EQUATIONS
In this section we will present the perturbed field equations that result from simultaneously considering linear structures on large scales, and non-linear single-stream structures on small scales. The quantities that appear in these equations will then be explained further in the sections that follow, with technical details reserved for the appendices. We hope this will allow the busy reader to see the most physically interesting aspects of this work without first having to undergo the lengthy mathematical considerations that went into their derivation. It will, however, require the reader to temporarily suspend their satisfaction of knowing exactly what all the quantities involved represent, as we then move on to more carefully lay out their specific technical definitions and properties. The multi-stream generalization should follow straightforwardly.
The field equations we present will be expressed in terms of the following set of gauge-invariant gravitational fields: as well as a corresponding set of gauge-invariant matter perturbations: where ρ, p and v i correspond to the energy-density, pressure and peculiar velocity, respectively. These quantities all represent fluctuations about a spatially-flat FLRW geometry, which in a particular choice of coordinates can be written as where a is the scale factor and τ is a conformal time coordinate. In what follows we will also require the Hubble rate in conformal time, which we write as H ≡ a ′ /a. There are a couple of things that require some explanation before we can proceed further. First, U has been removed from the scalar potentials φ and ψ because we intend it to correspond to the Newtonian gravitational potential. A precise explanation of what we mean by this will come in the following sections. For now it sufficient to understand U as the leading-order part of the gravitational field produced by non-linear density contrasts. The potentials φ and ψ, as well as h ij and S i , contain information about both the large-scale cosmological potentials and the small-scale post-Newtonian potentials. Again, the precise meaning of this will become clearer in subsequent sections. Likewise, the Newtonian density contrast is δρ N , and the cosmological and post-Newtonian density contrast is given by δρ. The former of these is allowed to be arbitrarily large, while the latter is required to be small (again, to be made precise later on). Similar comments apply to v Ni and v i .
After simultaneously expanding the field equations in post-Newtonian and cosmological perturbation theories we find the leading-order parts are given by the effective Friedmann equations whereρ =ρ M +ρ R andp =p R are the leading order parts of the spatial averages of the energy density and pressure, respectively. They have both radiation (ρ R andp R ), and dark and baryonic matter (ρ M ) contributions. The Newtonian gravitational field equation occurs at the same order in our expansion, and is given by where η is the expansion parameter for the post-Newtonian expansion. In this case it is used to characterise the size of structure on scales of order the homogeneity scale; the largest-scale at which the post-Newtonian expansion can sensibly be performed. The reader may note that only dark matter and baryonic matter contribute to δρ N , and not radiation. Subsequent orders of the perturbation expansion in the field equations yield the following two equations for the scalar part of the gravitational field: (6) and where D ij ϕ ≡ ϕ ,(ij) − 1 3 δ ij ∇ 2 ϕ is the trace-free second derivative operator on any field ϕ, and where perturbations in radiation, and dark and baryonic matter contribute to both δρ and δp. The reader will note that these equations contain extra terms when compared to standard cosmological perturbation theory. First, there are effective energy density and pressure terms, δρ eff and δp eff . These are solely due to the presence of non-linear structures on small scales, and are given explicitly in Eqs. (11) and (12), below. In other words, by writing the field equations in an effective fluid description, one can clearly identify that small-scale non-linearities lead to, amongst other things, an effective pressure on large-scales. Second, in the above equations, the potential U couples to h ij and there are extra source terms on the right-hand-side of these equations that are linear in φ and ψ. These interaction terms do not exist in standard cosmological perturbation theory and vanish in the limit in which non-linear small-scale structures vanish. In general, the interaction terms should be expected to produce mode-mixing between scalar, vector and tensor parts of the gravitational field on cosmological scales and coupling between different Fourier modes in Fourier-space.
The remaining parts of the gravitational field are the vector and tensor modes. For the vectors we find that we can write the following single equation to describe S i , accurate up to order O(η 5 ): We can take the leading-order part of this equation, at O(η 3 ), and write it as the following simple Poisson equation The leading-order part of the vector gravitational field, given by the solution to Eq. (9), is only sourced by small-scale quantities, and is a hundred times greater than might naively be expected from cosmological perturbation theory. This is the equation that was identified in the post-Friedmann approach of [20], and solved for numerically in [25]. For the full vector equation (8), accurate up to O(η 5 ), it can be seen that there exists sources on both small and large scales and mode-mixing (which is missing from [27,28]). This equation has an effective energy flux, Q eff i , which is due to small scale potentials. It also has extra source terms on the right-hand-side that are linear in S i . Both of these vanish when small-scale structures are absent. The explicit expression for Q eff i is given in Eq. (13), below, along with the other effective fluid quantities.
The final field equations we require, in order to complete our set to the desired accuracy, is given as follows: where angle brackets around indices indicate a symmetric and trace-free operation has been used, so that T ij ≡ T (ij) − 1 3 δ ij T kk for any field T ij . This equation can be used to determine the tensor part of the gravitational field, h ij . It also has an effective fluid source, Π eff ij , which this time acts as an effective anisotropic stress and is formed from the quadratic contractions of the lower-order small-scale potentials, see Eq. (14). Again, the non-linear structure on small scales couples the large-scale scalar and tensor parts of the cosmological gravitational fields, and again we have additional terms on the right-hand-side that are linear in h ij , resulting in mode-mixing.
As promised, the effective fluid quantities in the perturbation equations above are given as follows: It can be seen that each of these quantities was constructed only from variables that correspond to small-scale gravitational fields, or background quantities (which will be shown later to be calculated from the average of smallscale quantities). We therefore have a hierarchy of equations that can be solved order-by-order: first, the Friedmann and Newtonian equations (4) and (5), and then the equations that contain large-scale perturbations (6)- (10). The former of these sets are already calculated routinely in modern cosmological N-body simulations. The latter are modified versions of the usual cosmological perturbation equations on large scales, and can be used to find post-Newtonian equations on small scales (as recently solved for numerically in [22,[26][27][28]). Finally, note that the above effective quantities, in Eqs. (11)- (14), contain terms that would normally only be included in second or third order in cosmological perturbation theory. In particular, the term δρ N v N i v Nj in Eq. (14) would appear at third order in standard perturbation theory, but here should be expected to source a gravitational "slip" in the leading-order part of the large-scale physics. Our approach can be compared to the effective fluid approach studied previously in [24,29], as well as the large and small wavelength split used in [30,31]. As a final comment, before moving on to explain the origin of these equations and give detailed explanations of the gauge invariant quantities involved, we note that the usual trick of separating equations like (8) and (10) into scalar, vector and tensor parts is much more difficult to apply here. This is due to the fact terms like (D ij U)(φ + ψ) do not have scalar, vector and tensor parts that are easy to identify. This term, for example, is a scalar multiplied by a tensor, and in general should be expected to contain scalar, vector and tensor parts. This does not mean that such a separation is impossible -indeed we very much expect it to be possible. It just means that the resulting equations are very messy to write down, which is the reason why we have chosen to present these equations without such a decomposition. The reader should also be warned that manipulation of these equations is considerably more difficult than in either cosmological perturbation theory or standard post-Newtonian theory. This is due to different derivative operators changing the order to the terms they operate on in different ways. This will be made clearer in the sections that follow, and expanded upon in more detail in a subsequent publication [32].

III. PERTURBATIVE FRAMEWORKS
In order to explain the origin of the equations presented in Section II we first need to outline the perturbative expansions used in their derivation. These are the post-Newtonian expansion and cosmological perturbation theory, which will be used to describe perturbations on small and large scales, respectively. For the former of these expansions we use the expansion parameter η ≪ 1, while for the latter we use ǫ ≪ 1. We assume that any field Q can be expanded in both ǫ and η as follows: where Q (n,m) is a quantity of order O(ǫ n η m ). The characteristic length scales on which post-Newtonian and cosmological perturbations exist and vary on will be labelled L N and L C , respectively. Let us now briefly outline the essential features of these two different expansion schemes, before considering how they can be performed simultaneously.

A. Post-Newtonian theory
Post-Newtonian gravity is a slow-motion and weak-field expansion of the field equations, usually applied to systems that are much larger than their Schwarzschild radius but also much smaller than the Hubble radius [33]. Characteristic dimensionless velocities on such scales are usually small, from which it follows that time derivatives of fields are also small in comparison to spatial derivatives:φ where L N is the spatial length scale over which these perturbations are expected to vary, and ϕ represents metric potentials. The perturbative order-of-smallness of metric potentials is given by noting that the magnitude of peculiar velocities are small, |v (0,1)i | ∼ η, and assigning terms with time derivatives an extra order-of-smallness in η, compared to spatial derivatives. The magnitude of the gravitational potentials, in terms of η, can then be determined through the geodesic equations for freely falling particles and the field equations. This process can also be used to assign orders of smallness to the sources of stress-energy [33,34]. At this stage orders-of-smallness in η can be assigned to all matter fields and gravitational potentials. The perturbed energy density and pressure for a perfect fluid are given by the expansions ρ = ρ (0,2) + 1 2 ρ (0,4) + . . . and p = p (0,4) + . . ., respectively, where ellipses denote higher-order terms. The quantity ρ (0,2) is the Newtonian energy density, while ρ (0,4) and p (0,4) are the post-Newtonian contributions to energy density and pressure. The reader may note that there is no pressure term at O(η 2 ). This is no accident -if such a term were to be included then the stress-energy conservation equations would require it to be spatially homogeneous. This means that barotropic fluids with p = wρ and w = 0 do not fit into post-Newtonian gravity in a natural way, unless they are diffuse enough to be considered post-Newtonian (i.e. to occur only at O(η 4 ) or above).
Similarly, order-of-smallness can be assigned to the metric components. The minimal set of metric perturbations required to consistently expand a time-dependent background metric, g (0,0) µν (t), using the post-Newtonian expansion, is therefore given by It is important to note that the post-Newtonian formalism outlined above was formulated to include matter, but not radiation or a cosmological constant. This will be discussed further in Section III D below, in the context of incorporating post-Newtonian expansions into cosmology. For further details about post-Newtonian expansions, the reader is referred to the textbooks by Will [33], and Poisson and Will [34].

B. Cosmological perturbation theory
Cosmological perturbation theory is a weak-field expansion about the spatially homogeneous and isotropic class of Friedmann-Lemaître-Robertson-Walker (FLRW) geometries: where a(t) is the scale factor and k is the spatial curvature (we take k = 0 for the rest of this paper). All perturbations to all matter and gravitational fields in this perturbative expansion are taken to occur at the same order: where L C is the characteristic length scale gravitational fields vary on, and is necessary above to compare the dimensionless expansion parameter, peculiar velocity and gravitational potentials to dimensionful quantities like the perturbed energy density and pressure. Cosmological perturbation theory also assumes that time and space derivatives do not change the order-of-smallness of a quantity, so we havė This result occurs because the perturbations can be as large as, or larger than, the horizon, and is in stark contrast with the post-Newtonian expansion outlined previously. It means that orders at which terms appear in the two expansions cannot in general be expected to be comparable, a difference that changes the character of the resulting field equations. For our present purpose it is important to note that standard cosmological perturbation theory is designed in such a way that it can be applied to many different epochs in the Universe. Specifically, given a fluid with equation of state p = wρ then both the background and perturbed equations can be written down straightforwardly. This allows the theory to be applied to both the radiation-dominated and matter-dominated stages of the Universe's evolution, as well as to the current cosmological-constant-dominated epoch. This is a versatility that is absent from the standard approach in post-Newtonian gravity, as radiation and Λ are completely negligible for the study of gravity in the Solar System, binary pulsars, and other such very-small-scale astrophysical environments. If one wants to apply such expansions to super-clusters in a cosmological context, however, then more care may be required. For further explanation of cosmological perturbation theory, the reader is referred to the review by Malik and Wands [35].

C. Two-parameter perturbation theory
In reality, both post-Newtonian and cosmological perturbations should be expected to exist in any realistic model of the Universe [36]. In this section we review and extend the two-parameter framework developed in [23] that simultaneously performs a perturbative expansion in both sectors. In the next section we will extend this formalism to include radiation and a cosmological constant. Our two-parameter expansion in both ǫ and η will be constructed around an FLRW geometry, corresponding to the line-element in Eq. (17). This is the standard background for cosmological perturbation theory, but is so far little used for post-Newtonian gravity (see however [20,21]). Nevertheless, it can be shown that both expansions can be performed in such a background in an entirely self-consistent and well-posed way [23,37,38].
The first step in doing this is to expand the total energy density and pressure in both ǫ and η: The terms ρ (0,0) and p (0,0) can be considered as the background energy density and pressure, as they are not perturbed in either ǫ or η. All other terms correspond to perturbations at the order indicated by the superscript, but we have neglected to include δ symbols before them to keep the notation as compact as possible. To be even more precise, the orders-of-magnitude of these perturbed quantities are given by where {m, n} ∈ N + , and again L C and L N are the characteristic length scales of the cosmological and post-Newtonian sytems, respectively. A similar expression holds for the expansion of p. The length scales are necessary in the denominators of these expressions, as ρ is a quantity with dimension L −2 , and because it only makes sense to compare the magnitude of quantities with the same dimensions. The first thing to notice about Eq. (20) is that the mixed-order terms do not always appear at the same order as the product of post-Newtonian and cosmological terms (i.e. we have included ρ (1,1) , even though there is no O(η) term in the post-Newtonian expansion). The reason for this is that such terms are necessarily generated by arbitrary gauge transformations, and so cannot be assumed to vanish in general, even if they are assumed to do so in one particular coordinate system. Secondly, we have included a background energy density, which was absent in Ref. [23]. This will be useful for including radiation.
We also need to expand the metric in both ǫ and η, which we do as follows: (1,2) 00 (1,2) 00 where in the second line of each of these equations we have chosen our background metric g (0,0) µν to be the flat FLRW metric from Eq. (17), and simultaneously defined the perturbations h µν . The orders of magnitude of each of the perturbations to each of the components of this metric are the minimal set required to self-consistently account for the gravitational fields of the two-parameter perturbed perfect fluid discussed above, in any arbitrary coordinate system. We find that the inclusion of radiation and a cosmological constant does not require the introduction of any new metric potentials at any new order, so the form of Eqs. (23)-(25) is the same as in Ref. [23].
The final ingredient of the field equations that must be perturbed is the peculiar velocity, v i . This is split into post-Newtonian and cosmological parts such that which leads to the following components of the reference four-velocity u µ : (1,0) 00 which are derived using the normalization condition u µ u µ = −1, and Eqs. (23)- (25). The components of the two-parameter perturbed energy-momentum tensor that arise from these equations are given in Appendix A. The components of the Ricci tensor are unchanged from Ref. [23], and can be found in the appendix of that paper.
The reader should note that within the context of the two-parameter formalism, time derivatives are taken to add an extra order-of-smallness, η, compared to spatial derivatives whenever they act on an object that contains any non-zero perturbation in its post-Newtonian sector. So, for example, we takė As in Eq. (22), the purpose of this is to reflect the expectation that quantities perturbed in the post-Newtonian sector should be slowly varying in time and change over spatial length scales L N , while quantities that are perturbed only in the cosmological sector should vary equally over both time and length scales L C . This is explained further in [23].

D. Including radiation and Λ
Let us now consider how to add radiation and Λ to our two-parameter expansion. For radiation this can be achieved by writing where ρ M and p M are the energy density and pressure of non-relativistic matter, ρ R and p R are the energy density and pressure of radiation, and v i M and v i R are the peculiar velocities of the matter and radiation fluids. We then want to expand each of these new quantities in ǫ and η, which we do according to and These equations can, of course, be compared to Eqs. (20), (21) and (26) to read off values for the perturbations to the total energy density, pressure and v i . They can also be seen to generalise those from Ref. [23] by the inclusion of ρ R and p R , as well as by the inclusion of an extra mixed-order term p The reader will note that the expansions of the matter and radiation fluids have not been performed in an identical way: We have omitted (i) a time-dependent background-level contribution to the matter energy density and pressure, and (ii) a Newtonian-level contribution to the radiation energy density and pressure, so that The former of these is neglected because it corresponds to a term that could otherwise be as large as the Newtonian rest mass energy density ρ (0,2) M , which we consider to be highly unphysical. In the real universe there is no time-dependent background matter component to the energy density, ρ (0,0) M (t). This is because the leading-order contribution to ρ M is in fact dominated by the (inhomogeneous) rest mass of galaxies, dust etc., which is exactly what ρ Let us now consider the expansion of ρ R and p R given in Eqs. (32) and (33). For this purpose it is useful to consider the stress-energy conservation equation for the total stress-energy tensor T µν : where T Mµν and T Rµν are the matter and radiation contributions to the total stress-energy tensor, respectively. This implies ∇ µ T Mµν = Q ν and ∇ µ T Rµν = −Q ν , where Q ν = 0 for interacting fluids and Q ν = 0 for non-interacting fluids.
In either case, the lowest-order part of Eq. (35) is given by which implies p is a function of time only. If we now take p R = 1 3 ρ R , then this result implies that the leading-order part of the energy density in radiation must also be spatially homogeneous, such that ρ This is, in fact, exactly what is required for a background-level contribution to the energy density in an FLRW model.
A similar argument can now be used to understand why it would be inappropriate to include a term ρ (t) would be Newtonian in size, and such a term would be highly unusual in normal post-Newtonian gravity. We therefore find that the lowest order at which inhomogeneous perturbations in radiation fit into our two-parameter expansion is at order O(p (1,0) R ) ∼ O(ǫL −2 C ), which corresponds to a cosmological-scale perturbation. The reader may also note that there is no term ρ (1,1) R = 0 in any coordinate system, and the same result will hold in any other coordinate system related by an infinitesimal gauge transformation.
Finally, let us consider the cosmological constant Λ. We assign an order of magnitude and dimensions to the cosmological constant in the following way: This choice is motivated by the fact that the cosmological constant in the standard model of cosmology must be of background order, in order for it to be influential in the Friedmann equations at late times. There is also no point in perturbing it in either ǫ or η, as it is a constant, and the Taylor expansion is trivial. The cosmological constant therefore fits naturally into our two-parameter expansion at lowest-order, as a cosmological background quantity with corresponding scale L −2 C . The full field equations, in arbitrary coordinates, are given in terms of the perturbed quantities introduced in this section in Appendix B. In the sections that follow we will perform gauge transformations in order to determine how these quantities transform between different coordinate systems. We will then construct a set of gauge-invariant quantities that obey the same equations in any coordinate system, before introducing combinations of potentials that can be taken in order to write the gauge-invariant field equations in terms of an effective fluid.

IV. CONSTRUCTING GAUGE-INVARIANT VARIABLES
A general infinitesimal gauge transformation between coordinate systems can be written as where ξ µ is the gauge generator. All tensors, T , are taken to transform under the gauge transformation in Eq. (38) asT where L ξ denotes the Lie derivative operator with respect to ξ µ . This exponential map results in an invertible transformation, and can be applied to both the metric and the stress-energy tensor. We must now expand the components of the gauge generator in terms of ǫ and η, which we do as follows: These non-vanishing components of the gauge generator have been chosen so that no new components of the metric or the stress-energy tensor are generated by this transformation, which is an important condition to ensure the problem is being treated in a self-consistent manner.

A. Infinitesimal coordinate transformations
In order to perform infinitesimal coordinate transformations it is useful to decompose the perturbed gauge generator into a scalar and a divergence-free vector. Omitting superscripts, these can be written as where δx i ,i = 0. The transformation of the metric perturbations due to a gauge transformation of this type are unchanged from the dust-only case, and are given in Section V of Ref. [23]. In the remainder of this section we will outline how the presence of radiation affects the transformation properties of the matter fields {ρ, p, v i , Λ}. This is done using the result from Eq. (39), and by solving for the decomposed matter variables.
In order to present these results in a form that can be used for cosmology we choose to take L N /L C ∼ η. This means that we are restricting the post-Newtonian sector of our expansion to apply on scales below about 100Mpc, which is realistic also about the size of the homogeneity scale. This is ideal for considering the influence of galaxies, clusters and super-clusters on large-scale linear cosmological perturbations. We also choose, without loss of generality, to express our results in terms of L N . Given this, the total energy density transforms as follows: while the total pressure transforms as The transformations in Eqs. (44), (46) and (49) remain exactly the same as the dust-only case studied in Ref. [23], while all other transformations are affected by the presence of the radiation. The term ρ (0,0) can be seen to transform in the same was as the Newtonian energy density, ρ (0,2) . This is not unexpected, as both quantities have magnitude Similarly, ρ (0,0) appears alongside ρ (0,2) in the transformation given in Eq. (45). With the inclusion of radiation, we find that p (0,0) is automatically gauge invariant. Furthermore, as can be seen in Eq. (48), the inclusion of radiation means that the transformation of the cosmological and mixed-order perturbations to the pressure are no longer gauge invariant (as they were in the dust-only case). The reader may note that these results differ from the quasi-static limit of cosmological perturbation theory, as space and time derivatives are treated on a different footing, and because velocities come in at different orders [41].
Meanwhile, the peculiar velocities transform in the following way: These transformations are the same as in the dust-only case studied in Ref. [23]. Note particularly that in Eq. (50) the small-scale Newtonian velocity contributes to the transformation of the large-scale velocity -this is a by-product of our two-parameter expansion, and is an effect that would otherwise only appear at second order in standard cosmological perturbation theory. Finally, we find that the cosmological constant Λ (0,0) does not transform under the gauge transformation in Eq. (38), as it is a constant in space and time: The transformations above will now be used to construct gauge-invariant quantities.

B. Gauge invariant quantities
Let us now create gauge-invariant quantities for the matter degrees of freedom in the presence of radiation and Λ. Such variables isolate and remove superfluous degrees of freedom, as well as allowing the field equations to be written in a greatly simplified way. To do this it is useful to perform an irreducible decomposition on the metric. Omitting superscripts for simplicity, and without loss of generality, we can do this as follows: where B i , F i ,ĥ ij are divergenceless andĥ ij is trace-free. Applying the gauge transformation (39) to the metric components (23)-(25) this gives the transformation rules for the irreducibly decomposed components, and allows gaugeinvariant gravitational perturbations to be constructed (see Sections V and VI of Ref. [23]). The presence of radiation does not affect the construction of gauge-invariant gravitational perturbations, but does affect the construction of gauge-invariant quantities for the matter variables, which is what we will elaborate upon here. The method we use to calculate gauge-invariant quantities is as follows: we choose gauge generators δx, δx i and δt such that the gauge transformed metric potentialsẼ =B =F i = 0. We then substitute these gauge generators, now written in terms of E, B and F i , back into the expressions for all of the transformed perturbations presented in Section IV. Because the original gauge transformations were written down in a completely arbitrary coordinate system, these new results are automatically gauge invariant [35]. All such quantities also reduce to metric perturbations in longitudinal gauge when E = B = F i = 0, and have been explicitly checked to be truly gauge invariant.
To construct gauge-invariant matter perturbations we require the transformation laws for E (1,0) , B (1,0) , F (1,0)i , E (0,2) and F (0,2)i under Eq. (38). These are given in [23] and arẽ For the total energy density perturbations it can then be seen that the following quantities are gauge invariant: Correspondingly, for the pressure perturbations we find the following gauge-invariant quantities: (1,0) (64) and for the peculiar velocity we construct These last two quantities can be separated into scalar and divergenceless vector parts in a straightforward way. Finally, the gauge-invariant cosmological constant is trivial to construct: There are no further quantities to consider in the stress-energy tensor, so, when combined with the set of gaugeinvariant metric potentials {Φ, Ψ, B i , h ij }, constructed in Section VI of Ref. [23], this gives us a full set of gaugeinvariant quantities in our two-parameter perturbative expansion. The field equation in terms of these gauge-invariant variables are given in Appendix C.

V. CONSTRUCTING THE FIELD EQUATIONS
The two-parameter expansion described in the previous sections could in principle be applied to numerous different physical systems. While the perturbed metric and stress-energy tensor can be written down without specifying any specific relationship between either ǫ and η or L C and L N , we must choose how to express these quantities in terms of one another if we want to be able to solve a hierarchical set of field equations. In order to model a realistic universe that has non-linear structure on scales up to ∼ 100Mpc, as well as linear structure on large scales, we choose L N /L C ∼ η. On the other hand, to model a realistic universe, gravitational potentials must have similar magnitude on both small and large scales, so we choose ǫ ∼ η 2 , see Section III of Ref. [23] for justification of this. Both of these requirements are therefore satisfied by the choice where 10 −5 is the typical depth of a potential on both cosmological and post-Newtonian scales. With these relations we can translate our two-parameter expansion into effectively a single-parameter expansion in η, and write the field equations order-by-order in η. We further choose to express the field equations in units of L −2 N . This last choice has no particular physical significance, and is purely for expediency. We use Eq. (69) in the work that follows, as well in Section II and Appendices B and C.

A. Background and Newtonian cosmological equations
Within the formalism outlined above, the Friedmann-like equations that govern the evolution of the scale factor a(t), and hence the large-scale expansion of the Universe, are not independent of the perturbations. This can be seen explicitly in Eqs. (C3) and (C4) from Appendix C, where the Newtonian mass density and gravitational potential act as sources for the cosmological expansion. This is in some sense a very pleasing result; the large-scale expansion of space is driven by the same Newtonian mass that governs the leading-order part of the gravitational field on small scales. On the other hand, it means that our "background" is not by itself an exact solution of Einstein's equations. This stretches the meaning of what is usually implied by the phrase "perturbation theory" in Einstein's theory 3 . Nevertheless, both the fundamental objects being perturbed and the field equations themselves are being consistently expanded in the perturbative parameters ǫ and η, and we see no reason to expect this expansion should not converge. Indeed the present expansion seems to have much better convergence properties than the standard approach to cosmological perturbation theory, in the presence of non-linear structures [39]. Furthermore, a change of coordinates on a sub-horizon-sized region of space can be shown to be isometric to perturbed Minkowski space, with the cosmological expansion arising from boundary conditions at the edge of the region [37]. In this sense, the cosmological expansion can be considered an emergent property, and the background on small-scales could equally well be considered to be either a Friedmann model or Minkowski space (which definitely is a solution when ǫ = η = 0).
In any case, we can now proceed in a similar manner to Section VIIA of Ref. [23] to find the simplest way in which to express the equations that govern the large-scale expansion of space. In order to do this, it us useful to calculate the average mass density and radiation density on distances above the homogeneity scale, L hom ∼ 100Mpc [42]. These are given by where V hom indicates the spatial volume associated with the homogeneity scale. Of course, we know from Eq. (36) that there can be no small-scale inhomogeneities in the radiation fluid. For the matter fluid, on the other hand, small-scale fluctuations most definitely do exist and are of order unity. To accommodate these fluctuations we define This equation implies that the leading-order inhomogeneous part of the matter energy density, δρ (0,2) , is formally of the same order as the background component of the matter fields, ρ (0,2) , both being O(η 2 L −2 N ). These quantities can now be used to write Eqs. (C3) and (C4) into a more useful form.
To derive a set of effective Friedmann equations we first integrate Eq. (C4) over the volume corresponding to the homogeneity scale: where H ≡ȧ/a. Using Gauss' theorem this can be written as If we now assume that on the homogeneity scale there is no net flux of ∇Φ (0,2) into or out of the surface S hom , then the second term in Eq. (73) vanishes. This leaves us with which is exactly the same form as the standard Friedmann equation in the presence of matter, radiation and a cosmological constant. What is more, the lowest-order parts of the stress-energy conservation equations yields the results [40] ρ M ∝ a −3 and which are again exactly as expected from Friedmann cosmology. Substituting these results back into Eq. (C4) gives which is identical to the standard equation used in Newtonian N-body simulations for cosmology. In summary, we find that the leading-order parts of the field equations, in the context of our two-parameter expansion, reproduce exactly the same results as standard Friedmann cosmology with dust, radiation and a cosmological constant (although the meaning of the equations is slightly different). In the following section we will find that this is not the case when the non-linear aspects of Einstein's equations become important, on large scales. Finally, we note that the Friedmann equation is recovered from our expansion when δρ (0,2) = 0 (i.e. when the leading-order contribution to the energy density is homogeneous). This is not the same condition as setting ǫ = η = 0, which would correspond to an empty space within our framework.

B. Leading-order cosmological perturbation equations
The equations presented in Appendix C constitute a hierarchy of field equations, where the equations from Section V A are the leading-order parts. Once the Friedmann equation (74) and the Newtonian equation (76) have been solved, then their solutions can be substituted into the remaining higher-order equations to gain a set of solutions for the leading-order cosmological perturbations. This latter set of solutions, at O(η 4 L −2 N ), contain linear-order cosmological large-scale potentials and post-Newtonian potentials from small scales. With this in mind, we therefore seek to recast the O(η 4 L −2 N ) equations in the form of the equations of standard first-order cosmological perturbation theory, modified by the addition of terms related to the existence of inhomogeneity on the length scale L N . These terms will be then be treated as forming the components of an effective fluid on large scales, whose characteristics and behaviour is determined by the small-scale gravitational physics. Such an approach has similarities to the effective fluid approaches in for example [24,29], but in our case is also required to reduce the number of gravitational degrees of freedom to be no more than the available number of field equations.
In the end, we want to reduce to a set of six perturbed field equations for six degrees of freedom (i.e. the 10 degrees of freedom in the metric minus the four coordinate freedoms). At present, Eqs. (C5)-(C14) from Appendix C contain a total of sixteen degrees of freedom: six scalars (Φ (1,0) , Φ (1,2) , Φ (0,4) , Ψ (1,0) , Ψ (1,2) and Ψ (0,4) ), six in the tensors (h ) and four in the vectors (B (1,0) and B (1,2) ). Taking into account the four degrees of freedom removed by gauge fixing implies that we need to remove six degrees of freedom. This is achieved by defining new sets of variables as follows: and which are exactly the sets of variables used in Section II. A number of these new variables could be considered to be "composite quantities", as they contain a number of different perturbative orders in the same variable. For example, the variable ψ is dominated by O(ǫ) terms on cosmological length scales L C , but contains smaller terms at O(η 4 ) on small-scales L N . This is quite atypical in perturbation theory. However, the way in which these quantities arise together in the field equations suggest that they should be solved for together. The equations displayed in Section II can now be seen to be given by those from Appendix C, but written in terms of the variables in Eqs. (77)-(86) and in terms of conformal time dτ ≡ a −1 dt. To solve this system, one must solve the Newton-Poisson equation and Friedmann equations for U and a(t), respectively. This should not be too challenging, as they are identical to the standard equations used in Newtonian N-body simulations for Friedmann cosmology. From these results one can then solve for the effective fluid parameters, and then solve the cosmological perturbation equations. This will inevitably be complicated by the additional "mode-mixing" terms in the cosmological perturbation equations, which will require more sophisticated techniques than at leading-order in standard cosmological perturbation theory. These terms suggest that it may in fact be possible to generate vector and tensor modes from scalar fluctuations, which is already well known in second-order cosmological perturbation theory [17,19], but is not usually seen at first order. One should also note that these terms, for example 8πa 2 3 δρ N (ψ − φ), also mean that Fourier modes no longer decouple in a trivial way as they do in standard first-order perturbation theory, even if no mode-mixing occurs. This is because the Fourier transforms of such terms are expressible only in terms of a convolution integral over all Fourier modes.

VI. CONCLUSIONS
We construct a two-parameter perturbation expansion around an FLRW background that simultaneously describes non-linear structures on small-scales and linear structures on large scales. Moreover, it includes matter, radiation and a cosmological constant. In doing so we use both cosmological perturbation theory and the post-Newtonian expansion. As this expansion is able to model large density contrasts and different matter components it therefore both contains the essential features of the real Universe and has a number of potential advantages over standard cosmological perturbation theory. We derived the two-parameter perturbed field equations valid for structure on the order of a fraction on the horizon size, the two-parameter gauge transformations of the matter sector, and construct gauge-invariant quantities in this sector. The consistency of the gauge transformations requires not only gravitational potentials and matter perturbations at the orders expected from post-Newtonian gravity and cosmological perturbation theory alone, but also a number of others at orders in perturbation which may not naively have been expected. We have therefore identified a minimal set of perturbations that are required for mathematical consistency of the problem, and written down gauge-invariant versions of the field equations that contain them all.
We find that the small-scale Newton-Poisson equation for the scalar gravitational potential occurs at the same order in perturbations as the Friedmann equation, but that they can be separated after the introduction of a suitable homogeneity scale. At leading order, this results in a small-scale Newton-Poisson equation sourced by the inhomogeneous part of the Newtonian energy density, and large-scale Friedmann equations sourced by the spatial average of the leading-order parts of the energy density, pressure, and the cosmological constant. Our results give no indications that the effects of small-scale non-linearities should be expected to cause acceleration of the large-scale Universe, but we do find that they should be expected to affect large-scale perturbations. This is because the higher-order field equations include quadratic Newtonian potentials within the effective fluid terms. They therefore contain valuable information about non-linear gravity, and could potentially be used to identify relativistic effects in observations of large-scale structure.
By presenting the higher-order field equations in terms of an effective fluid we are able to highlight the similarities and differences between our formalism and regular cosmological perturbation theory. We expect this to aid further application of our equations by allowing some standard techniques from cosmological perturbation theory to be imported. Our effective fluid description also enables an easier physical interpretation of the effects of non-linearities in the field equations, which clearly lead to (for example) a large-scale effective pressure and anisotropic stress. Since the effective fluid terms are all constructed from the solution to the short-scale Newtonian gravitational potential, their properties should be able to be determined from N-body simulations. Once the form of these effective fluids has been identified, one can proceed to solve the cosmological equations for the long-wavelength perturbations. This method of solution is available to us because of the hierarchical nature of the perturbation equations -short-scale fluctuations appear at lower-order compared to cosmological perturbations, and so can be solved for before cosmological perturbations. Within this prescription we observe a mixing of scales, as well as mode-mixing at what would normally be considered to be linear order in cosmological potentials. Understanding the consequences of these relativistic effects on the formation of non-linear structure in the Universe is of importance not only for removing sources of observational bias, but also because it has the potential to offer new ways of probing Einstein's theory with cosmology.
The terms on the right-hand side of Eq. (A1) are given by This appendix contains the field equations in terms of the variables introduced in Section III, with the choice of relations between ǫ, η, L C and L N given in Eq. (69).
Note that a ∼ 1 andä ∼ 1/L 2 C , as the time variation of a(t) is over cosmological scales. This equation is a combination of the Raychaudhuri equation from Friedmann cosmology, and the Newton-Poisson equation from post-Newtonian gravity. The leading-order contribution to the trace of the ij-field equation at O(η 2 L −2 N ) is given by ȧ a This equation is the same as given in Ref. [23] because neither the cosmological constant nor radiation contribute trace-free components. The other equations contain contributions from both radiation and the cosmological constant.

Vector potentials
The 0i-field equations give the governing equations for the vector gravitational potentials. The leading-order contribution is O(η 3 L −2 N ) and is given by 0j,ij − aḣ This is the equation for the small-scale post-Newtonian vector potential, responsible for phenomena such as the Lense-Thirring effect, and is the one studied in Ref. [20]. Interestingly, this field equation implies the gravitomagnetic potential is ∼ 100 times larger than second-order perturbation theory predicts [23,46]. The next-to-leading-order 0i-field equation occurs at O(η 4 L −2 N ), and is given by The 00-field equation at O(η 2 L −2 N ) can be written as a a + 1 6a 2 ∇ 2 Φ (0,2) = − 4π 3 ρ (0,0) + ρ (0,2) + 3p (0,0) + 1 3 Λ , where we have substituted in the results from Eq. (C2). These equations govern the leading-order part of the gravitational field, at O(η 2 L −2 N ).

Vector potentials
The 0i-field equations at order O(η 3 L −2 N ) give (1,0) i , which can also be split into scalar and divergenceless vector part using a derivative. The reader may note that the quadratic term, which includes the lower-order potential Φ (0,2) , does not source the vector part of Eq. (C6).