{'fraction_non_alphanumeric': 0.031328695738898446, 'fraction_numerical': 0.010733787119455458, 'mean_word_length': 4.823230490018148, 'pattern_counts': {'":': 0, '<': 2, '<?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}}
Does the growth of structure affect our dynamical models of the universe? The averaging, backreaction and fitting problems in cosmology

Structure occurs over a vast range of scales in the universe. Our large-scale cosmological models are coarse-grained representations of what exists, which have much less structure than there really is. An important problem for cosmology is determining the influence the small-scale structure in the universe has on its large-scale dynamics and observations. Is there a significant, general relativistic, backreaction effect from averaging over structure? One issue is whether the process of smoothing over structure can contribute to an acceleration term and so alter the apparent value of the cosmological constant. If this is not the case, are there other aspects of concordance cosmology that are affected by backreaction effects? Despite much progress, this 'averaging problem' is still unanswered, but it cannot be ignored in an era of precision cosmology.


I. INTRODUCTION
A crucial feature of the real universe is that structure occurs on many scales. Large-scale models of the universe, such as the standard models of cosmology, are coarse-grained representations of what is actually there, which has much more structure than is represented by those models. An important issue for cosmology is the question of whether the smaller scale structures influence the dynamics of the universe on larger scales: is there a significant back-reaction effect from the small scales to the large scales?
The standard Friedmann-Lemaître-Robertson-Walker (FLRW) models of cosmology are a resounding success. They can account for all observations to date with just a handful of parameters. These are the simplest reasonably realistic universe models possible within General Relativity: homogeneous, isotropic, and flat to a first approximation, with a scale-invariant spectrum of Gaussian perturbations from inflation added on top to account for structures down to the scale of clusters of galaxies. Only the physical motivation for the value of the cosmological constant and the associated coincidence problem is required for a complete understanding at late times. Cosmology in the future will be refining this picture.
But is it as simple as this? The universe may well be statistically homogeneous and isotropic above a certain scale, but on smaller scales it is highly inhomogeneous, quite unlike a FLRW universe. General Relativity is a theory in which spacetime itself is the dynamical field with no external reference space. Yet it is ubiquitous in cosmology to talk of a 'background' which is exactly homogeneous and isotropic, on which galaxies and structure exist as perturbations. Is this the same as starting with a more detailed truly inhomogeneous metric of spacetime, and progressively smoothing it -probably by a non-covariant process -until we get to this background?
It is commonly assumed that whichever way one goes about it the results of these two processes must be the same; but the non-linearity of the field equations ensures that they are not. The problem for cosmologists is whether this difference is important. This has been the subject of controversy, with opinions ranging from the suggestion it could completely explain the recently observed accelerating expansion of the Universe without a need for any dark energy [1][2][3][4][5][6][7] to the claim that it is completely negligible [8][9][10][11][12][13][14][15][16][17][18][19], with a middle-ground of papers claiming it is enough to disturb the cosmic concordance of the standard model [20][21][22][23][24][25][26][27][28].
Averaging has kinematic, dynamical, and observational aspects. One estimate of the importance of this effect would be to note that within the standard model, variations in the expansion rate of around a few percent occur on scales of order 100 Mpc [20,22,23,25,[28][29][30][31]. At the very least, then, the averaging problem may be important for observations in precision cosmology; indeed it is obviously built into standard observational practice.
But does it have dynamic consequences? Green and Wald [19] claim to derive in a rigorous and systematic way the effects of small scale inhomogeneities on large scale dynamics, and thereby prove that matter inhomogeneities produce no new effects on the dynamics of the background metric: they can only mimic radiation. However others reach an opposite conclusion: Buchert's approach [1,2] shows backreaction from inhomogeneity can potentially mimic dark energy; Kolb and colleagues concur [32], and even have suggested inhomogeneity effects can fully explain the acceleration of the universe [5,33], as has Wiltshire [34,35]. Recently Baumann et al. [18], using second-order perturbation theory, have claimed that virialised structures freeze out of the cosmic expansion and only affect the background dynamics through a renormalised mass. Clearly there is no agreement, but even from a conservative viewpoint there is evidence that backreaction might have an impact on precision cosmology: Baumann et al [18] conclude that inter alia it significantly affects the Baryon Acoustic Oscillations (BAO).

This review
Our aim in this paper is to present the different approaches to the dynamical aspects of backreaction effects in recent eras in cosmology, highlighting open questions and so future research agendas. While we will briefly touch on them, we will not deal comprehensively with observational aspects of the fitting problem: this raises a number of further issues we do not have space to deal with here [137]. We will also not deal with possible backreaction effects in the very early universe, which raise a rather different set of issues (see e.g. [36,37]).
There have been quite a number of attempts to estimate the backreaction effect. One set is based on starting with an inhomogeneous model and using some averaging process to produce an approximately FLRW model. This has been developed both in terms of general formalisms, and in terms of detailed model building. Alternatively, by averaging over structure in the standard model, the backreaction effect can be calculated perturbatively; this should then provide an estimate of its size, and, in principle, allows us to evaluate the fluctuations which become of the order of the mean, and the background one starts with no longer correctly describes the averages.
We first discuss non-perturbative attempts to estimate backreaction. Then we provide an overview of backreaction in the standard model, where it appears the effect is small -at least to second-order in perturbation theorybut nevertheless is at least sufficient to affect precision cosmology results.

II. AVERAGING, BACKREACTION, AND FITTING
A. The Lumpy universe and the Background Any description of physical reality embodies a scale of averaging that is usually not made explicit. The world looks completely different depending on the scale of description chosen: for example a fluid looks wholly different on an everyday scale (say 1 m) as opposed to atomic scale (say 10 −11 m). The same is true for astronomy and cosmology. The background models in standard cosmology are the FLRW models, which ignore all local structural details. More realistic models include linearised perturbations about this background, in principle representing the largest scale growing structures down to the scale of clusters of galaxies (10's of Mpc) from early times to today; but they do not represent the non-linear smaller scale structures, such as our galaxy or the solar system, at later times. The universe looks far from homogeneous when viewed on any scale from 1 AU to 1 Mpc, and only approaches statistical homogeneity past 100 Mpc (though even this is in dispute [38,39]) -see Fig. 1.
To describe the Universe on its largest scales, one has to make approximations that postulate or derive a high degree of symmetries for the metric of space-time. Prac-tically, this means that one wants to calculate the large scale observables using a background geometry, i.e. a geometry that ignores, on average, the details present on small scales and that are not probed by the observables. Such a background is usually found in a FLRW solution to the Einstein Field Equations (EFE). This issue is referred to as the fitting problem [41]: what is the best-fit FLRW model to the lumpy Universe? In the standard concordance cosmology, the existence of this background is postulated, and no smoothing mechanism is provided to obtain it from the real lumpy Universe. Despite the success of this approach in fitting the observations there still remains the problem of properly defining the background geometry in relation to the real lumpy universe.
For example, which set of observers (worldlines) are associated to the background, i.e. are supposed to measure an homogeneous and isotropic Universe? It is clear that such a fitting procedure needs an explicit method of averaging or smoothing. This averaging can be performed with different techniques and on different quantities. For example, one can argue that homogeneity is a spatial property and average quantities that define particular spatial hypersurfaces, such as densities, pressures, or expansion rate of geodesics bundles. In that case, the background is defined by surfaces of constant density, pressure or expansion rate and the average quantities are used to fit the background. Another possibility is to fit the model via averaged observable relations such as the magnitude-redshift or number count-redshift relations; then the average has to be performed in some sense on the past-null cone of the observer. Of course, these procedures will, in general, give different results and the FLRW fitting model reconstructed will depend on which method has been used.
The construction of the background is a crucial issue: if the wrong background is compared to the data, it will imply the existence of a backreaction that may disappear if a better background is chosen. This emphasizes the fact that a gauge choice is always part of an averaging procedure. As such, all the approaches listed and commented on in what follows somehow intend to clarify the (usually unstated) way this is handled in the standard approach to cosmological modeling in most papers on cosmology.
We have three closely related but distinct problems to consider: Averaging Coarse-graining of structure, such that small-scale effects are hidden to reveal large scale geometry and dynamics.
Backreaction Gravity gravitates, so local gravitational inhomogeneities may affect the cosmological dynamics. How this is calculated depends on the degree of coarse graining.
Fitting How do we appropriately fit an idealized model to observations made from one location in a lumpy universe, given that this 'background' does not in fact exist?
FIG. 1: Structure in the Millennium simulation [40] (from [26]). Can we describe the universe as smooth on scales of order 150Mpc, shown here in the black and white boxes (top panel)? The averaging problem is shown in the bottom row: how do we go from left to right? Does this process give us corrections to the 'background', or is it the 'background' itself? How does it relate to the 'background' left at the end of inflation?

B. Averaging and Backreaction
The basic issue is the non commutativity of averaging and the field equations. Start with a realistic description of the universe on a small scale, with metric g (local) ab (e.g. this might describe individual stars and planets in the universe, and the vacuum between them). Average it by a smoothing procedure to a metric g (gal) ab with an averaging scale where galaxies are well represented but individual stars are invisible. Average this in turn to a metric g (lss) ab with an averaging scale where large scale structures are well represented but individual galaxies are invisible (there are many possible scales that are omitted in this a description). The largest scale (completely smoothed) model will have a FLRW metric g (cos) ab , where all traces of inhomogeneity have been removed. There will similarly be averaged stress energy tensors T representing the matter present at each of these scales. Now the Einstein equations may be assumed to hold at the 'local' scale: after all, this is the scale where they have been exquisitely checked, so But the averaging process: ab , (2) doe not commute with evaluating the inverse metric, connection coefficients, Ricci tensor, and Ricci scalar, for example: Hence if the EFE hold at scale 'local' (i.e. equation (1) is true), they will not hold at scales 'gal' or 'lss'; for example one will find where the extra term E (gal) ab = 0 is due to this noncommutativity. It is the effective matter source term representing the effect of averaging out smaller scale structures, which is then an effective source term for averaged EFE at scale 'gal'. Similarly there will be such an effective such source term E (lss) ab at that scale -the scale usually represented by perturbed FLRW models -and E (cos) ab at the cosmological scale. This is the backreaction from the small scales to the larger scales.
In essence, it is an assumption that Einstein's equations also hold for an averaged geometry, as well as a local one. In fact, it is not clear that the whole machinery of GR holds after averaging -e.g., concepts such as spacetime and objects such as tensors need to be assumed to make sense after coarse graining.
The classic example of this effect is Isaacson's calculation of the effective backreaction of small scale gravitational radiation on an averaged large scale metric [42,43]. Now this affect applies to cosmology: if there is significant gravitational radiation at early times, this will affect the dynamics at later times, according to Isaacson's calculation (which can also be derived from a variational principle, see [44]). However that effect will be exceedingly small at late times (it may be significant at early times, see [36,37]).
We are concerned however with the dynamic effects of non-linear structure on cosmology at late times (i.e. after decoupling of matter and radiation). In principle, the same effect will occur in this context [45]; the question is whether this is a significant effect or not. We will look at general formalisms and specific models in the next section.

C. Fitting and observations
One can do a spacetime fitting, asking which FLRW model is best if we average invariant quantities in a spatial or spacetime volume: e.g., the energy density of particles and their velocities (when the matter averaging may be represented by kinetic theory); one may choose to smooth the metric or scalar invariants on the geometry side. Alternatively, one can do a null fitting, where one in effect averages astronomical observations.

Fitting and the past null cone
As we determine our best fit model by null cone observations, this is in the end what we'd like to do (in effect the standard observational cosmology approach is this type, but is not usually phrased this way). So a key issue is how this all relates to cosmological observations. Here we have to take into account not only the effect of averaging on the geometry and dynamics (as represented by (5)) but also the effect of lumps on null geodesics and on observational relations.
The basic question is: What are observables in these approaches? What do observations of averaged quantities really tell us about our Universe? Cosmological observations probe quantities such as the redshift, the angular diameter distance, the luminosity distance, and the image distortion, i.e. quantities related to light emitted by a distance source and propagating on the past null cone. The general way this works was described in a perturbative framework in the pioneering work of Kristian and Sachs [46] (extended to more general cases in [47][48][49]). This involves expanding the geodesic devi-ation vector in orders of the affine parameter along the past null cone. Here one relates the geometric quantities to the observational quantities using the relation between the redshift and the four-velocity of the observer u a : 1 + z = uak a |e u b k b |o , and the relation between the intrinsic cross-sectional area, dA of a distant source and the measured solid angle subtended by the observer dΩ: dA = r 2 A dΩ, where r 2 A is the angular-diameter distance. This expansion as given by [47,48] is exactly covariant and can be applied to any space time. In the present context one would like to compare the general result with that of the FLRW background. This is done by matching the background FLRW results with that from the spacetime indicated by observations [9,50].
However there are complications. In the real universe, as pointed out in [46,48,51,52], observations take place via null geodesics lying in the empty spacetime between galaxies, which are focused only by the curvature actually inside the beam, not the matter that would be there in a completely uniform model. The effect on the observational relations of introducing inhomogeneities into a given background spacetime is twofold: it alters the redshift, and it changes area distances. This should also be taken into account in any fitting procedure.
Others have emphasized the importance of smoothing the past lightcone. A sketch of how averaged Raychauhduri equation on the past null cone will look like was given in [53], while Räsänen in [54,55] gave the past null cone averaged equations for the scalars in a statistically homogeneous and isotropic space time. Work in this important area is still at its infant stage with no concrete quantifiable physical result that can be directly fitted to observations.

Averaging on the past light cone
Averaging is in some respects a fitting process, but does not necessarily correspond to any actual observational procedure. Can one propose an average model of the universe based on the past null cone? A very recent attempt at a comprehensive approach is Gasperini et al. [56]. However this should be approached with caution: observations to cosmological distances mix spatial and time variation, and it does not make sense to simply average today's state of the universe with what it was like in the past. One would expect any averaging operation to leave the background invariant, and it's not obvious that this can happen for FLRW somehow averaged on its light cones. So averaging based on observations would need to involve comparing the universe today with earlier times by use of dynamical equations relating variables at these different times: a very model dependent procedure, and not 'averaging' in a normal sense.
What is clear is that backreaction won't be important along the light cone, because the key causal effects in cosmology propagate in a timelike fashion [57]. Nevertheless it is obviously important to relate the results of averaging and backreaction effects to observations. The exact and perturbation approaches that follow attempt to do this.

III. NON-PERTURBATIVE BACKREACTION
One approach is to build a model of the universe 'bottom up'. Advocated by Buchert [1,2,58], Zalaletdinov [59,60], Wiltshire [34,35,61], Räsänen [6,54,55], among others [62][63][64][65][66], these approaches dispute the idea that one needs a background to work from: rather the background model and its dynamics should emerge as a large-scale approximation to a more detailed inhomogeneous model, which can be compared with a FLRW model for the universe assumed ab initio as in the standard approach. As explained in previous sections, the field equations for such a model may be expected to be different from those in a standard FLRW model: we want to understand that difference.
Non-perturbative approaches fall into two main categories. One are generic averaging formalisms, which aim to understand the nature of the backreaction terms in general, a bit like deriving and understanding the macroscopic Maxwell Equations. The other approach is to create fully relativistic inhomogeneous or even N-body models by use of simplifying assumptions; then, by comparing observables in these models with their averaged FLRW counterparts one can hope to quantify non-perturbatively the backreaction effect and the magnitude of the fitting problem. Let us consider each approach in turn, with some of the main attempts in the literature.

Buchert's approach
Alongside early attempts by [20,[67][68][69][70], Buchert [1,2] builds on the Newtonian averaging by Buchert and Ehlers [71,72] to provide a bare-bones approach to the problem, concentrating on averaging scalars on spatial hypersurfaces. The kinematic scalar equations for vorticity-free perfect fluid are averaged, to give evolution equations for the averaged expansion and shear scalars. Several authors [2,14,25,54,73,74] have generalized this approach to any arbitrary spacetime, but we will illustrate with the original Buchert proposal.
Let us assume a dust spacetime, and observers and coordinates at rest with respect to the dust. The average of a scalar quantity S may be (non-covariantly) defined as simply its integral over a region of a spatial hypersurface D of constant proper time divided by the Riemannian volume: Taking the time derivative of Eq. (6) yields the commutation relation where Θ is the expansion of the dust, and we assume the domain is comoving with the dust. The dimensionless volume scale factor is defined as a D ∝ V D 1/3 , which ensures Θ D = 3∂ t ln a D . Then, the second derivative of the scale factor is given by the averaged Raychaudhuri equation: where is the kinematic backreaction term and σ 2 = 1 2 σ ab σ ab is the magnitude of the shear tensor. The non-local variance of the local expansion rate can act in the same way as the cosmological constant, causing the average expansion rate to speed up, even if the local expansion rate is slowing down. Even more tantalisingly, if this were the cause of the observed acceleration, the coincidence problem would be solved in the most natural way: as structure forms the variance in the expansion rate grows, as matter coalesces and virialises [3]. This is a truly remarkable possibility in moving from local to non-local quantities on a non-trivial geometry, and is the reason for the recent excitement in the averaging problem.
One can see how this counter-intuitive idea works as follows [75]: If the average scale factor of a universal domain, a D , can be written as a union of locally homogeneous and isotropic regions, each with its scale factor a i , then the acceleration of the universal domain D is given by [6,76,77]: where a i represents the locally defined scale factor in the i-th sub-region, and, a D ≡ (a 3 1 + a 3 2 + · · · ) 1/3 . Acceleration of the universal domain,ä D > 0, can easily be achieved, for example, for a two disjointed dust filled FLRW sub-region in which one might be expanding while the other is contracting at a time t, i.eȧ 1 = −ȧ 2 (assuming same sized sub-regions a ≡ a 1 = a 2 ), one obtains Here one easily obtain an acceleration for the universal domain,ä D > 0, even when the two sub-regions are deceleratingä 1 < 0, a 2 < 0; i.e all observers see only deceleration. However, it has been argued [8] that acceleration found in this toy model does not necessary imply that the physical universe is accelerating, since this model has not been shown to satisfy other rigorous observational tests. It also ignores problems due to matching/junction conditions for any two regions.
In Buchert's scheme, all tensor contributions appear as scalars in these equations, and are collected into unknown source terms. Of course, the system of scalar equations is not closed, so one has to make an ansatz about the effect of averaging the shear terms; so it is very difficult to say how big the backreaction effect is. One can derive an evolution equation for the averaged shear scalar; but that would be sourced by products of Weyl curvature tensors, amongst other things, and one quickly sees that the system of equations can never close. So the method of averaging only scalars reaches this limitation quickly. This feature can further be understood by considering the integrability condition: where R D is the average local curvature. This coupling between the curvature and the volume scale factor implies that if R D ∼ a −2 D as in FLRW cosmology, the kinematic backreaction term will scale as Q D ∼ a −6 D , which mimics the behavior of some kind of dark fluid.
Having said that, [78] has suggested that, if the scalar curvature invariants can uniquely characterise any spacetime, a scalar averaging scheme can work in general by averaging these invariants. He thus arrives at a complete, closed way of averaging spacetime using only scalars. Whether it is practical remains to be seen.
The averaged quantities in Buchert's formalism do not have a clear observational meaning. Nevertheless, it is worth noticing that [54,55] argues that in a statistically homogeneous and isotropic Universe, these average quantities are exactly the ones that describe observations along the past lightcone. It will be interesting to see if such an argument can be made rigorous.

Zalaletdinov's Macroscopic Gravity
A comprehensive approach to covariantly averaging tensors was initiated and explored by Zalaletdinov [59,60,[79][80][81] and later by others [65,66,82,87]. It is a foundational attempt to average the complete set of Cartan structure equations, in order to define Einstein field equations for the averaged quantities. This approach is directly inspired by the way in which a macroscopic theory of electromagnetism can be obtained from the microscopic Lorentz-Maxwell theory [60]. Once averaged, the 'macroscopic field equations' resemble Einstein's, but with a source term -analogously to the polarization term in macroscopic Maxwell equations [83].
The core issue of averaging tensors covariantly is managed by using bi-local extensions of tensors, so that they transform as tensors at some point of interest, x, but as scalars in a neighbourhood of the point. Because of that, they can be averaged over that region Σ. The covariant spacetime average is defined as [59] T is the bi-local transport operator. The backreaction generated by the smoothing procedure takes the form of a correlation tensor for the gravitational degrees of freedom which appears as an effective source in the Einstein field equations. Alternatives to this definition, and the resulting averaging scheme were initiated in [66]. Implementing this operation results in a completely covariant smoothing procedure provided the transport operators satisfy certain conditions. The result is formally independent of the averaging scale and as such can be seen as generating a universal description of the collective behaviour of local gravitational degrees of freedom when only their very large scale properties matter (much the same way as thermodynamics encompasses the collective behaviour of particles on very large scales compared to the particles themselves).
This approach makes an extensive use of transport of tensorial quantities along geodesics, but by using the natural parallel transport bitensor the metric is invariant so that no smoothed metric is obtained. Hence, it relies on a specific choice of bi-tensor that satisfies a set of differential equations and conditions, in order to imply the correct properties of the average. This bi-tensor is used to evaluate integrals on finite domains, and it is not clear how the formalism is affected by the choice of this bitensor [65]. Additionally his averaged Einstein equations rely on some "splitting rules" (see eqns (45) and (48) of [79]) which can be questioned (although they are consistent with an analysis of high frequency gravitational waves, see eqn (68) of the same paper).
In [84][85][86][87] it is shown that in a flat FLRW macroscopic background, the correlation tensor is of the form of a spatial curvature, while [88] showed that Zalaletdinov's MG reduces to Buchert's equations with corrections in an appropriate limit. It has further been employed to evaluate the backreaction effect in a perturbed FLRW model [89,90], but the requirement that a FLRW background exists makes this attempt fall under the category described in Sec. IV A. While the amplitude of backreaction is similar to that obtained by simpler means, the details of using such a covariant approach will be important for an effective fluid description of perturbations at second-order, such as [18].

B. Other approaches
Another rigorous approach is based on the deformation of the spatial metric of initial data sets along its Ricci flow [91,92]. In principle, this method is a nice, natural way of smoothing a space-time that can be linked to the standard renormalization group approach of effective field theories [67], but the non-linearity of the Ricci flow equations is a serious complication that can lead to the development of singularities along the flow; that makes its use in a cosmological context particularly difficult. A renormalisation group approach to coarse graining in the very early universe is given in [37].
C. Model building Approaches

Timescape Cosmology
Another very original viewpoint called the Timescape cosmology, has been proposed and investigated by Wiltshire [34,35,93]. It is a brave but contentious attempt to seriously look at the status of bound regions and their interaction with an expanding cosmological model.
The idea is to separate the Universe into expanding, underdense regions whose boundaries are overdense regions enclosing virialized regions such as the one we, as observers, live in. An average is then performed spatially (using Buchert's formalism [58]) to define a reference cosmic background. Interestingly, the amount of backreaction in these models is at most of order a few percent (when normalized as a fraction of the energy density) and is not solely responsible for explaining the apparent cosmic acceleration: the non standard effects principally come from the desynchronization of local clocks (in the virialized regions) with respect to cosmic clocks defined via the average background. Indeed, the gravitational redshift effects imply different ticking rates for clocks inside the voids and in the virialized regions. Wiltshire argues that the effect is cumulative when an average is performed to define the background clocks. A possible interpretation of the model is that the extra redshift effects change the observable relation in the effective FLRW background. As such, the model is not actually accelerating, but the extra redshift accounts fully for the dimming of supernovae because they appear to be at a higher redshift than expected. A detailed discussion of the observational consequences of the Timescape cosmology can be found in [94,95].
Wiltshire's proposal is more original than simply viewing the problem as one of backreaction of structures on the overall global dynamics. It also recognises that the position of the observer (in virialized structures as opposed to voids) may be important to the fitting problem when the variance in local geometries becomes large. Nevertheless, it suffers from its own problems, among which one is the use of a pure two-zones model to describe the Universe, without proper junction conditions between the zones (this problem is emphasized in [96,97]).

Swiss-Cheese Models
The Swiss-Cheese model consists of one or more spherically symmetric vacuum regions, each described by a Schwarzchild metric, joined across spherical boundaries to a FLRW model [98][99][100][101][102]. This set-up represents a very natural way to model the part of the universe that we see, for example the boundary of a galaxy and intergalactic space, the lack of effect of the expansion of the universe on the motion of planets, etc.
The Swiss-cheese type construction can be adopted to study the effect of inhomogeneities on cosmological observations in a fully non-linear and relativistic manner.
The algorithm commonly in use is to start with FLRW metric and cut out comoving holes and fill them with a Schwarzschild metric, while making sure through the matching conditions that the mass in the holes equals the mass that was removed. Hence, Swiss Cheese models do not affect the global dynamics of a lumpy universe, i.e there is no dynamical backreaction effect.
Another more interesting construction is to modify the homogeneous FLRW metric by the introduction of spherical regions of the spherically symmetric Lemaître-Tolman-Bondi (LTB) dust spacetime [103][104][105][106][107][108][109][110][111]; a quasispherical Szekeres model can also be introduced [112]. One has to ensure through the matching conditions that the spherical region of inhomogeneities is comoving and mass compensating, to ensure that the LTB and FLRW regions evolve independently. Thus again there are no dynamical backreaction effects.
There are however significant observational differences from standard cosmology. These spacetimes model precisely the difference between Weyl and Ricci focussing of null geodesics. The null geodesics in the void regions are focussed only by shear induced by the Weyl tensor. One of the critical problems in this area is the relation between Ricci curvature and Weyl curvature, and Swisscheese models are a very interesting way to study this. We refer to other papers for a discussion of these observational effects: see [103][104][105][106][107][108][109]111] for differing viewpoints.

Lindquist-Wheeler type models
All the preceding models relied on the hypothesis that a cosmological fluid can be employed to model the distribution of matter in the Universe. On the contrary, Lindquist-Wheeler type models are a genuine approach at constructing an expanding universe model out of locally static domains. It consists in modelling the Universe by approximately paving its compact spatial sections [138] (topologically homomorphic to S 3 ) with Schwarzchild domains that stand for the static regions constituting the 'particles' of cosmology (such as galaxies) [113]. Of course, the matching cannot be exact, and shells have to be introduced at the boundaries between the cells. These boundaries then obey equations of motion that produce an overall expanding and recollapsing model that closely mimics a k = 1 FLRW Universe.
This approach fundamentally differs from a Swiss-Cheese model in that no reference to a FLRW metric is needed in addition to the static regions: the dynamical properties really emerge from the interaction between the static cells encoded in the motion of their boundaries. As such, the fluid approximation is not required and the fluid-like behaviour only appears when the dynamics is coarse-grained over the detailed, local, structure. An interesting point is that it leads to a solution that is similar to the equivalent FLRW k = +1 solution, but different in the details (for example, the relation between the total mass and the maximum radius is modified).
It has recently been extended interestingly by Clifton and Ferreira [114,115] in their Archipelagean cosmology. They showed that the optical properties of such a model are different from those of the "equivalent" FLRW model and can lead to a correction in the fitted value of Ω Λ of order 10%. However, the solution is not, as such, self consistent: it is an approximation that neglects the interaction between neighbouring cells, interaction that would result in the deformation of the geometry around each vertex and the appearance of anisotropies. But the approximation is plausible and is worth exploring because it is a neat way to explore the emergence of a collective expanding Universe formed by locally static regions.

IV. PERTURBATIVE BACKREACTION
All these approaches are enlightening, but do not yet result in detailed models that can be directly compared with precision cosmology observations. For that we need to turn to the standard perturbed FLRW models, which enable us to comprehensively investigate backreaction effects in the linear and weak non-linear regime.
Whether this adequately reflects the dynamics of full non-linear models can be debated. Paranjape [116] emphasizes that the background scale factor, needed to compute the backreaction is affected by the backreaction itself, hence making it impossible to calculate the backreaction precisely as circularity occurs. We assume that structure formation can be described by perturbing a well behaved background (since it seems to give the right power spectrum and indeed everything that has to do with observations of structure formation). The background that has to be fitted to observations, however, is not the one we first thought of with scale factor a(t), but the one with coarse-grained scale factor a D (t). Thus fluctuations change the background that has to be compared with large scale observables, such as the luminosity distances. We return to the issue of the adequacy of these approaches below.

A. Backreaction in the standard model
The standard model of cosmology ignores all the complexity of smoothing the spacetime and assumes that on 'large' scales (say larger than a few Mpc) we can model the universe as homogeneous and isotropic, with linear fluctuations describing structure propagating as smooth fields on this background. On smaller scales we can jump to Newtonian gravity, and model the universe as discrete particles in simulations. Because it is the only model we have where we can calculate anything realistically at all, it is the perfect arena to study backreaction in detail. We shall give a rough overview of the issues involved, which illustrate more about the backreaction problem in cosmology.
The fields that are propagating on the background alter its dynamics through the non-linearity of the field equations. Could the gravity of the gravitational poten-tial be important? According to observers at rest with the gravitational field, the potential itself is small everywhere outside objects less dense than neutron stars, and so if we write g = g 0 + h where h 1, how might the backreaction of this perturbation h add up to something large? The affine connection or Ricci rotation coefficients determine the dynamics of the spacetime, and are generically O(∂h); the field equations are O(∂ 2 h). A perturbation of wavelength λ = a/k, where k is the comoving wavenumber and a is the scale factor, much less than the Hubble scale can give rise to large fluctuations in the field equations O(k 2 h), even though the change from the background metric is small. Such terms describe density fluctuations, and these can be large even though the metric potentials are small. Can this change the background significantly?
Let us consider this in some detail, in the simplest cosmology which agrees with observations: a flat LCDM model with Gaussian scalar perturbations. Averaging FLRW perturbations has been discussed in different guises in the literature (mostly for an Einstein-de Sitter model): some authors investigate specifically the modification to the Hubble expansion rate or other variables [3, 20-25, 28, 64, 69, 106, 117, 118]; others reformulate the average of the backreaction into an effective fluid [14,18,27,82,88,[119][120][121][122][123], while one of the first attempts considered the important problem of how to calculate the averaged metric [70]. Rather than summarize these approaches, let us discuss what happens in a general way.
In the Poisson gauge to second-order in scalar perturbations the metric reads [18,124] (12) where ∂ i V i = 0, h i i = 0, ∂ i h ij = 0 (because scalar, vector, and tensor modes interact at second order, it is inconsistent to include scalar modes alone). The background evolution of the scale factor a(t) at late times is determined by the Friedmann equation where the Hubble constant H 0 is the present day expansion rate, and Ω m the normalised matter content today. The first-order scalar perturbations are given by Φ, Ψ (and are all that is required for observations at the moment), and the second-order by Φ (2) , Ψ (2) (which are needed for a consistent analysis of backreaction). In this gauge we have the metric in its Newtonian-like form, which we may think of as the local rest-frame of the gravitational field because it is the frame in which the magnetic part of the Weyl tensor vanishes for vanishing vector and tensor perturbations [25]. For a single fluid with zero pressure and no anisotropic stress Ψ = Φ, and Φ obeys the 'master' equation For a LCDM universe the solution in time to this is equation has Φ constant until Λ becomes important, and then starts to decay as Λ suppresses the growth of structure on all scales by about a factor of 2. There is no scale dependence in the equation, which all comes from the initial conditions -usually a nearly scale-invariant Gaussian spectrum from frozen quantum fluctuations during inflation -and subsequent evolution during the radiation era. In Fourier space, assuming scale invariant initial conditions from inflation, the power spectrum of Φ, P Φ , is independent of scale for modes larger than the equality scale, k eq ≈ 0.07Ω m h 2 Mpc −1 , and ∼ k −4 (ln k) 2 for modes much smaller than it, up to some non-linear scale k N L k eq . The change in behaviour at the equality scale, arising from modes which enter the Hubble radius before matter-radiation equality, is important for backreaction because it is the modes larger than the equality scale which are primarily responsible for any backreaction at all. In essence, the equality scale determines the size of the backreaction effect.
All first-order quantities can be derived from Φ; for example, is the first-order velocity perturbation, which governs the peculiar velocity between the matter flow and the restframe of the gravitational field. Meanwhile, the gaugeinvariant density perturbation is The second-order solutions for Ψ (2) and Φ (2) are given by [124]. These are complicated expressions involving time integrals over products of Φ and its derivatives. For backreaction, however, the important thing is that these contain terms of the form ∂ i Φ ∂ i Φ, as well as non-local terms such as What is the backreaction of perturbations onto the the expansion rate? At second-order, after substituting for v (1) i from Eq. (15), the expansion rate looks something like [3, 14, 21-23, 25, 28, 82, 88, 117, 118, 121, 122]: and time derivatives thereof.
Provided the domain is small, this quantity can also correspond to the sky-averaged Hubble rate, measured from observations on the past nullcone [4,118]. In principle, then, the perturbative terms give the backreaction to the local expansion rate. To evaluate H we can use a realisation of Φ given an inflationary model. Alternatively, we can assume a spectrum for Φ and evaluate the statistics of H. This allows us to calculate the expectation value of H, H, as well as its variance, in terms of integrals over the power spectrum of Φ multiplied by powers of k. The reason we must go to second-order now becomes clear when we calculate the expectation value: for Gaussian perturbations from inflation, the ensemble average of Φ is zero, which implies -assuming ergodicity -that when averaged on the background over a large (strictly, infinite) domain they are zero too. Thus, the second-order terms provide the main backreaction effect; the first-order terms give the variance.

Averages used
Before going further, it is worth summarizing the different types of averaging that are used in this kind of analysis: Riemannian averaging this is the 'correct' way to average scalars in a non-homogeneous and nonisotropic spacetime. In perturbation theory it can be expanded in terms of the Euclidean average, introducing products of averaged quantities. It doesn't commute with the time derivative, implying more non-connected terms. Acting on a secondorder quantity, this is the same as Euclidean averaging at that order.
Euclidean averaging or Smoothing This is spatial averaging on the background. When thinking of perturbations in the metric as fields on the background this is the natural averaging to use. When thinking of the perturbed metric as an approximate spacetime in its own right, quantities should be smoothed intrinsically using the Riemannian averaging. Smoothing is the same, but often applied directly to Φ to remove structure below a given scale. It was used in [18] to integrate out small scales, where it was applied to tensor components as well as scalars.

Ensemble averaging
The key tool for statistical evaluation, assuming Gaussian perturbations from inflation. A Euclidean average over an infinite domain is the same as ensemble averaging, assuming ergodicity.
It is often assumed that these are interchangeable, but they are not. In particular, the scale dependence of averaging is very different depending on the procedure used, and this could have an impact observationally [22,25,28]. Furthermore, Euclidean averaging -or often even Riemannian -is often replaced with an ensemble average, loosing the non-connected terms which are responsible for scale-dependence. This gives only the super-Hubble contribution to backreaction. The exact relation between ensemble averaging and Riemanning averaging is another open important question.

Averaging
Consider now the average of H over a spatial domain D, where the spatial surfaces are defined in the coordi-nate frame. The Riemannian average of a quantity Υ, can be expanded in terms of the Euclidean average defined on the background space slices, Υ = D d 3 x Υ D d 3 x, as: where Υ (0) , Υ (1) and Υ (2) denote respectively the background, first order and second order parts of the scalar function Υ = Υ (0) + Υ (1) + Υ (2) . Note the important term in square brackets, which encapsulates the relativistic part of the averaging procedure. Distinct types of terms now appear in the averaged expansion rate: This smoothed Hubble parameter may be the closest variable to what is measured in practise [22][23][24], at least for small domains where the linear Hubble law applies. Again, the expectation value (ensemble average) and variance can be calculated using the statistics of Φ. Once the ensemble average is taken, the first-order terms drop out, the 2nd-order connected terms become independent of scale (as is the case for all terms in H), while only the non-connected terms retain scale dependence. Ensemble averages of pure divergence terms, such as ∂ i v i 2 (which actually contains (∂ 2 Φ) 2 terms -see below) also drop out [25]; these are usually dropped in backreaction studies because they are boundary terms, so by statistical homogeneity must be small [117].
A similar analysis has been carried out for average of the deceleration parameter q = −(1 +Ḣ/H 2 ), which quantifies the rate of change of the local Hubble rate; it can be extended by replacing H → H D to measure the deceleration of the average Hubble rate. These are different things, and have either ∂ 2 Φ∂ 2 Φ or ∂ 2 Φ ∂ 2 Φ terms. Such terms also appear in the deceleration parameter defined via a series expansion of the distance-redshift relation, which is a physical observable [4,118].
The relations for determining the scaling behaviour for the backreaction terms are where W is an appropriate window function specifying the domain. Note that the connected terms have no dependence on the domain size at all. Given the approximate behaviour of P Φ ∼ ∆ 2 R ∼ 10 −9 for k k eq , and P Φ ∼ ∆ 2 R (k eq /k) 4 ln(k/k eq ) 2 for k k eq , the scaling of the backreaction terms may be estimated. For this we need to replace ∞ 0 → k U V k IR . Then, the key scalings for backreaction become: Terms like these appear in the average of the Hubble rate, as well as its ensemble average. The second type of term, Φ∂ 2 Φ ∼ k 2 eq , arising from peculiar velocity terms, is the term primarily responsible for setting the fundamental amplitude of the backreaction in the Hubble rate. It is quite small, for the concordance model. (The overall effect is somewhat larger than this due to the contribution of several such terms.) So, we're looking at sub-percent changes to the Hubble rate from backreaction, though non-connected terms make it larger on small scales. Yet, we can observe that the backreaction is small because the equality scale is large in our universe, which is because the temperature of matter-radiation equality is very low. Modes which enter the Hubble radius during the radiation era are significantly damped compared to those which remain outside until after equality; so, the longer the radiation era, the less power there is on small scales to cause a significant backreaction effect. For a scaleinvariant spectrum, then, it may be considered that the long-lived radiation era is the reason that the dynamical backreaction is small. The temperature will have to drop by several orders-of-magnitude before backreaction in the Hubble rate becomes significant.
The first type of term, Φ 2 , is nominally much smaller, O(10 −10 ). However, it also tells us that the IR divergence in Φ 2 must be cut-off by hand. This is interesting because it implies that a scale invariant spectrum can not go on forever, and so in this sense the universe cannot be infinite. This implies that backreaction could in principle be used to place limits on the start on inflation, which governs the largest mode k IR , which is the first mode to leave the Hubble radius during inflation [4,21,33,125]. There has been speculation that this might lead to very important effects and even mimic dark energy [4,5,21], though this has been criticised [9][10][11]. However, for this to be significant we require k IR /k eq ∼ 10 ∼10 4 (where Φ 2 would compete with the other terms), which is quite large. Alternatively, given that we only measure the primordial power spectrum to be nearly scale-invariant over a comparatively narrow range of scales, the appearance of Φ 2 implies that it can't be too tilted to the red on super-Hubble scales. A red spectrum would convert the logarithmic divergence into a power-law one, and constraints on the largest mode would be much stronger (e.g., k IR /k H 10 ∼80 for n s = 0.95).
In both the variance of the Hubble rate, and the backreaction on q, much larger terms appear: ∂ 2 Φ∂ 2 Φ, which are of order the density fluctuation squared. (While they appear in the Hubble rate through the second-order velocity perturbation, the ensemble average conspires to cancel them out.) This now behaves as, ∂ 2 Φ∂ 2 Φ ∼ ∆ 2 R k 4 eq × [divergent integral], and so: which is pretty significant in size, and diverges as k U V → ∞. The function F is roughly F (x) ∼ 0.5x 2.1 for 1 x 10, ∼ 70x −0.1 (log 10 x) 4.75 for x 1, and approaches ∼ 53 ln 3 x as x → ∞. Importantly, F overcomes the pre-factor around k U V ∼ 10k eq , so these terms are big, and difficult to know what to do with: here, the UV cutoff is really a measure of our ignorance. Within linear perturbation theory it should be set by the end of inflation and the reheating temperature, as well as the small-scale physics of dark matter, both of which are sub-pc scales today. Replacing the UV cutoff with a smoothing function in Φ implies that we might do better to smooth order-by-order, and calculate second-order terms from smoothed first-order ones, rather than average directly at a given order. Even for domains much larger than the non-linear scale, where linear perturbation theory breaks down (somewhere around a few Mpc), F is quite sizeable, and so we have backreaction terms O (1). From this, we also recover that the variance in the Hubble rate is O(1) on scales of Mpc.
The divergencies arise not only because of the high derivative terms which appear at second-order but also because of the scale invariant initial conditions and power-law suppression of modes below the equality scale. Without this suppression, even the Φ 2 terms would lead to significant backreaction, as modes all the way up to the inflationary cutoff would contribute. In such a scenario it is clear that the canonical approach to perturbation theory would fail completely. As we discuss below, a reformulation of perturbation theory might avoid these type of problems [18].
Would higher-order perturbation theory affect these conclusions? On the one hand it seems clear that higherorder perturbations should be suppressed. Provided Φ is Gaussian, only even orders will be important, once ensemble averages are taken. We might expect the largest terms at any order n to behave like Φ (n) ∼ (∂Φ) n (e.g., from relativistic corrections to the peculiar velocity), the ensemble average of which goes like ∆ n R k n eq . Terms which appear in the Hubble rate at order n of the form ∂ 2 Φ (n) do not have enough derivatives to overcome the suppression from (∂Φ) (n−2) terms. By this argument, secondorder should be as large as it gets, and backreaction from structure is irrelevant. On the other hand, others [6,117,126] have argued that at higher order terms such as (∂Φ) 2 (∂ 2 Φ) n−2 are the norm; in this case, from 4th order on perturbation theory formally diverges -at . The heavy black lines shows the deviation from the background, and the shaded region the ensemble variance as a function of averaging scale (that is, it is the variance one gets randomly throwing down Gaussian averaging balls of radius R). For the deceleration parameter there is some subtlety in the definition [118], which give different variances and scales: from outside to in we have the deceleration of the averaged scale factor (Buchert's definition, dashed), the spatially averaged deceleration parameter, and the ensemble average of the local deceleration parameter, now as a function of the UV cutoff (here evaluated using a smoothing function on Φ). The UV cutoff completely dominates the amplitude.
least as far as calculating averages is concerned. Even if not divergent, if (∂ 2 Φ) n−2 ∼ 1 then higher-order terms are at least as large as at second-order and must be included to evaluate backreaction properly, and so correctly identify the background [118]. But do these terms cancel out? This is definitely an open important problem. It is likely that an extension of the successful methods of [127], which have worked so well for Newtonian perturbation theory, will be important, but they must be extended to include vector and tensor degrees of freedoma significant challenge.
It is intriguing that the amplitude of backreaction can be dominated by the UV cutoff, and so difficult to quantify. In some ways this seems rather unphysical: we know that spherical systems cannot alter the expansion rate, which follows from Birkhoff's theorem. Recent work by Kolb [32] argues that this in fact sets the cutoff and so backreaction is dominated by that scale. Alternate work by Baumann et al [18], which we discuss further in a following subsection, has argued that when a system is smoothed over a scale larger than its size the virial theorem holds, in a pseudo-Newtonian limit. This holds in the Newtonian case where sources can be thought of as 'localised', in equilibrium, and so freeze out of the cosmic expansion. Of course, such a general theorem cannot exist in general relativity because energy is radiated to infinity, and only the stationary part of a system virialises [128]. Nevertheless, in the same way as for spherical structures this suggests the Newtonian virial scale is reasonable for the cutoff, as scales below this are shielded from the expansion. However, even a conservative cutoff at a few Mpc gives backreaction effects of order unity (see Fig. 2), the precise amplitude of which is completely dominated by how the cutoff is performed (a UV cutoff, vs a formal smoothing procedure, say). It is an open problem as to how this should be achieved in practise in a precise way: the cutoff is there to provide an ad hoc adjustment to the linear spectrum to simulate the non-linear, which has not been calculated.
One can argue that such arguments are in fact irrelevant. As a referee has commented to us: "As voids form they occupy a larger and larger fraction of the volume, since they expand faster than the mean expansion. The increase with time of the fraction of volume in which the expansion is faster than average means a decreasing average deceleration, and with smaller smoothing volumes fully developed voids can occupy larger fractions of the smoothing volume, increasing the variance. However, voids cannot expand faster than the average expansion indefinitely, since they collide with neighboring voids. Calculations [. . . ] which fail to take void collisions into account, cannot hope to give a correct account of the large-scale average deceleration once voids on any scale are fully developed. Since most of the matter ends up in the walls, filaments, and clusters between voids, which contain almost all the galaxies, what is relevant for the observed average deceleration on large scales is the deceleration of these walls and filaments with respect to each other, which has nothing to do with a volume-averaged expansion rate dominated by void interiors. The large scale deceleration is not a question which can be addressed by perturbation theory to the extent that backreaction coupling small-scale nonlinearities to the large scale dynamics has any importance." Though it is not clear how one envisages the wall/void interaction, the void dynamics will clearly be important. This alternative view highlights the kinds of interesting issues that remain unresolved.
We now discuss two recent approaches which may provide motivation around some of the problems we have discussed.

B. Short Wavelength Approximation
The problem associated with the ill-definition of an average for tensors was examined by Green and Wald in [19], generalizing earlier work by Burnett [129]. They replaced averages by the notion of weak limits in the weak field approximation, and thereby obtained strong restrictions on back-reaction effects (in this context) through a mathematically precise point limit process. In this approximation, they show that if the small-scale motions of matter inhomogeneities are non-relativistic, the effect of small scale inhomogeneities on large scale dynamics can be written as an effective trace-free (radiation) stress energy tensor, and hence cannot lead to acceleration via a negative active gravitational mass.
However the degree to which this analysis captures the physically relevant degrees of freedom is debatable because of the nature of the ultra-local limiting process, which ignores all details of the actual clustering of matter; but we shall see below that whether backreaction effects are significant or not depends crucially on the nature of that clustering. Furthermore, this analysis is based on the work in [129], which centers on handling a singular limit at the origin (a vanishing energy momentum tensor at all points away from the origin has a non-zero limit at the origin) arising through behaviour of the 'x sin(1/x)' variety. But does that realistically represent any real physical matter distribution in cosmology? Indeed it seems likely that the crucial quantity representing this discontinuity (µ mnabcs in [129], µ abcdef in [19]) will vanish for any realistic matter distribution: does the kind of ultra-local backreaction mechanism envisaged in these papers occur in physical reality? It is not generically the same as the backreaction due to averaging over finite volumes that is the concern of this paper, although it might be a limit of such a mechanism in specific singular geometric circumstances.

C. Effective Fluid Approach
In order to avoid the divergence problems caused when δ exceeds unity, [18] have considered a reorganization of the perturbative expansion, using a coordinate based Euclidean smoothing to separate long and short wavelength modes. Further, as we have discussed, v 2 ∼ (∂Φ) 2 ∼ Φ in magnitude, and on small scales when δ ∼ 1 a natural expansion variable is v 2 , provided each spatial derivative reduces the order by v. The field equations are linear in the matter variables, so there is no need to expand δ. Using this, averaging over suitable scales, one can derive effective pressure and densities which, when averaged, obey Newtonian-like equations [123] for the kinetic and potential energies on small scales and proving a viral theorem for local systems imbedded in the expanding universe. From this they argue that backreaction is always small, even in a model with no radiation era, and (in agreement with Green and Wald) cannot generate a negative active gravitational mass (ρ + 3p 0 always).
In effect, the potential for a negative gravitational mass to be generated by Buchert's averaging process is vitiated because of the existence of local equilibrium states characterised by their virial theorem. The main idea behind the effective fluid approach proposed in [18] is to re-write the Einstein equations into the background, forms linear in X, and those non-linear in X:Ḡ and to assume that the background equations,Ḡ ab =T ab , and the linearized Einstein equations, (G ab ) L = (T ab ) L , are defined in the standard way. Then the Einstein equations may be written in a form that is very similar to the linear equations, (G ab ) L = (τ ab −T ab ), where the effective stress-energy pseudo-tensor τ ab may then be defined as, τ ab ≡ T ab − (G ab ) NL . The second part in this process requires that the perturbation on the right hand side of the field equation be performed in orders of the peculiar velocity instead of density, since the density contrast is ill-defined at non-linear scales. Then at some scale Λ −1 each linear term is split into short wavelength modes and the long wavelength modes as, X = X + X s , and the non-linear splits as, After smoothing, the effective energy momentum pseudo tensor becomes, The superscripts s, and ∂ 2 denote the short wavelength, the long wavelength and suppressed higher derivative parts, and Λ is the cut-off for the effective theory. The tensor τ a b is conserved by virtue of the linearized Bianchi identity, and can be re-written into the form of a fluid with density, pressure and anisotropic stress, whereũ a is the renormalized matter 4-velocity and overline denotes ensemble average. Non-linear terms in τ ab may be re-written in the form of the kinetic energy κ and the potential energy ω, which evolve as d(κ + ω)/dt + H(2κ + ω) = 0. The effective density and pressure are given by ρ eff =ρ m (1+κ+ω) and 3p eff =ρ m (2κ+ω) , and its equation of state becomesw eff ≡p eff ρ eff = 1 3 (2κ + ω) . They argue that structure below the virial scale decouples from the effective long wavelength expansion of the universe and effective pressure vanishes: 2κ+ω = 0. The more detailed version of this proof in [18] relies on the fact that within the sub-horizon region, one can safely ignore the expansion of the universe (so one can set a = 1) and that the smoothing domain is much larger than the size of the system. Setting a = 1 and stating that a system is localised is equivalent to imposing a stationary orbit condition as was done in [128] for a general relativistic version of the virial theorem.
The effective theory they develop holds on large scales k Λ and contains a set of effective fluid parameters which must be determined from N-body simulations or directly from observations, or from higher-order perturbation theory.
This paper is a significant and interesting study of the back reaction issue, calculating some of the effects of averaging and taking the relevant scales into account, to see when it may be important. It brings together ideas of both Buchert and Zalaletdinov of deriving effective field equations and hence an effective fluid which hold on macroscopic scales. They also take care to motivate how virialised regions are cut off from the global expansion, as Wiltshire has emphasised. They approach the problem from the point of view of the standard model, allowing more quantitative predictions, but consequently suffer from relying on the assumptions of the existence of a global background geometry.
One of the major differences between [18] and other works lies in the fact that they are unconcerned with the problems of relativistic averaging, which is emphasised as critical by many authors discussed in Sec. III. All their smoothing is performed on the background, as all fields are considered as fields on that background. Averages of tensors are performed on tensor components. We can see from the general definition of the Riemannian average of a scalar, Eq. (18), that if Riemannian averaging were used instead extra terms such as Φδ ∼ v 2 δ would appear in their effective fluid which could change the nature of the result (though probably not the overall amplitude). Would this require just a re-definition of the effective fluid parameters, or something more significant? Related issues are that their results rely on working in a "good" gauge, and the assumption that coordinates exist which cover the whole spacetime on both small and large scales adequately; this is of course in common with most perturbative approaches.

D. Discussion
What can we conclude from this? One could conclude that dynamical backreaction is small by good fortune: the universe is so hot and had such a long radiation era that small-scale power is significantly reduced over its scale-invariant initial conditions. The backreaction terms in the Hubble expansion, Φ∂ 2 Φ, are the largest ones which appear in the lhs of the Einstein Field Equations, because the Einstein tensor has at most two derivatives of the metric in it. When backreaction is re-formulated as an effective fluid, these are also the largest terms which appear [18]. In some respects then this settles it: backreaction is small by virtue of there being a very small hierarchy of scales between the Hubble scale at equality and the Hubble scale today (they are only a factor of 50 apart in comoving terms). In this evaluation of backreaction, then, what happens on scales smaller than the equality scale is actually of little relevance. This is perhaps surprising given how we normally think of backreaction arising from small-scale structure -it is really power on very large scales which are responsible for the backreaction effect.
On the other hand, however, we should be able describe the universe using an orthonormal tetrad version of the Field Equations. This is entirely equivalent to the EFE, but reformulates gravity as a system of first-order PDE's in the Ricci Rotation coefficients and Weyl curvature tensor. Because of this, metric variables occur at up to two derivative levels higher than in the field equations. As we have seen with H, taking its derivative in perturbation theory can result in divergent terms appearing. Consequently, terms such as (∂ 2 Φ) 2 will occur frequently, and as we have seen, in q at least, can appear in both their connected and non-connected forms giving rise to large backreaction terms. That they can appear in purely kinematical quantities, and not just in a perhaps erroneous expansion of δ, was not addressed in [18]. The non-connected terms appear from the commutation relation for the time evolution of the spatial domain, so may be an artefact of the non-covariant averaging procedure. The connected form, however, is much more subtle to interpret as it appears even if we just calculate the expectation value of q -there's nothing really to do with averaging here. Where they appear in their connected form, they depend on how we cut-off the non-linear scales (assuming ergodicity). Rather intriguingly, a cutoff at the virial scale yields significant changes to the background. In this reading of backreaction, then, backreaction could be very significant indeed, and is large precisely because of the large hierarchy of scales between the non-linear scale and the Hubble scale, and is dominated by the UV cutoff. The large equality scale dampens the effect, but not enough for it to be insignificant.
Observationally, can there be any signature of backreaction? When measuring the Hubble rate, perturbations are significant in the variance of the Hubble rate on subequality scales [22,24]. Perturbations affect the whole distance-redshift relation which has been calculated to first-order by [4,12,13]. Corrections to the luminosity distance include corrections ∂ 2 Φ, just like for the Hubble rate. This allows us to see that the variance includes divergent terms like (∂ 2 Φ) 2 -the lensing term -which appear in the variance of the all-sky average of the luminosity distance -see [12], although this problem was not discussed. It is an important open question to find out what happens if their results are extended to secondorder, where an overall ensemble averaged offset to the luminosity distance will be present. Terms of the form (∂ 2 Φ) 2 appear in the series expansion of the distanceredshift relation, which can be used to define an observational deceleration parameter [4,118], but is this just an artefact of using a power series?
Finally, Baumann et al [18] claim that there are significant effects on the BAO. This means that backreaction effects are of significance to precision cosmology at these scales, even if they are not significant on the largest scales.

V. THE ALTERNATIVES
We can observe why the averaging problem is so difficult in GR, straight from the definition of the average of a scalar Eq. (6): to average or smooth a scalar quantity S it is not sufficient to just know S(t, x) as it would be in other areas of physics which have a spacetime prescribed; we also need to know h ij (t, x) too which requires the solution to the field equations in the nearby domain of interest. To prescribe a distribution of matter with mean energy density, say, requires us to know the full solution for the spacetime -before we can state its mean. Many of the problems we have discussed stem from this simple fact: describing an averaged spacetime accurately is just as difficult as modelling the full lumpy spacetime. Combined with the fact that averaging tensors covariantly is ill defined, the non-linearity of the field equations and so on, this is all quite a problem. That it is difficult, everyone agrees; what its effects are, on the other hand, depends on who is asking and what method they're using. Let us summarise some of the different viewpoints.

A. The Skeptic
The view that backreaction is negligible has been strongly argued by Peebles and Ratra [123,130], Ishibashi and Wald [8], among quite a few others [9-15, 18, 19]. The backbone of their argument is based on the fact that the gravitational potential Φ is small everywhere in the universe except in the immediate vicinity of black holes, and so the average of Φ, and all relevant physical quantities derived from it must consequently be small. While the density contrast can fluctuate by many tens of orders-of-magnitude, it is not this that causes backreaction because the field equations are linear in the density. Relativistic backreaction is not caused by density fluctuations per se; rather, it arises due to peculiar velocity contributions, which at the order relevant for backreaction are v 2 ∼ (∂Φ) 2 , and are ∼ 10 −5 in dimensionless form. While density fluctuations behave as ∼ ∂ 2 Φ on small scales these can only cause a very large variance in cosmological parameters, and this is only on scales which are small compared to the Hubble scale, about the non-linear scale of a few Mpc.
Although there are non-linear interactions of the gravitational potential with itself, this general relativistic effect of 'gravity gravitating' is tiny, ∼ Φ 2 ∼ O(10 −10 ). So backreaction and all perturbative effects on the background, while interesting, are essentially irrelevant for cosmology.

B. The Enthusiast
Some argue that the view of the Skeptic misses the point in its entirety: backreaction is a general relativistic effect arising from the fully non-linear field equations [4,5,35,131]. Arguments which claim that backreaction is small rely on perturbed FLRW which are inherently quasi-Newtonian. Such models miss the main possible backreaction effect both because they are Newtonian, and because they don't take inhomogeneity seriously. Even though they can be corrected to give some relativistic effects [132], this still only brings them into line with linear relativistic perturbation theory, which does not adequately capture the reality of a vast network of walls and filaments forming structures around expanding voids. Furthermore, Newtonian N-body simulations remove backreaction from the start because they employ periodic boundary conditions enforced on the background [71]. While such a condition is enforced, there can't be a net flow of particles into or out of the boxexactly what backreaction looks for -and the expansion is simply put in by hand. The divergent terms from perturbation theory (∂ 2 Φ) 2 cancel neatly in the backreaction terms because of their Newtonian nature.
Because the averaged field equations do not easily close, the backreaction terms have to be estimated from a realistic, fully non-linear, solution of the field equations, which are hard to make realistic -but that must be the aim.

C. The Fence-Sitter
While backreaction is clearly an important mechanism to understand, it seems difficult to believe it can really cause significant effects on the largest cosmological scales, both because of the large scale differences involved and the small value of the gravitational potential. However, that does not in itself imply that backreaction is insignificant in cosmology: perturbations about FLRW do not necessarily give small effects on smaller scales. Depending on the definition of the deceleration parameter, for example, (∂ 2 Φ) 2 terms either appear or cancel. They cancel in the deceleration parameter defined as the deceleration of the averaged scale-factor; but in the average value of the local deceleration parameter, they appear and give 10% effects even when power is artificially cutoff at 10 Mpc, and diverge if power on smaller scales is allowed to contribute [25,118]. This is hardly an unreasonable thing to calculate, yet perturbation theory doesn't give a sensible answer. Similar terms may appear at higher order in perturbation theory. Consequently, it seems too early to say that backreaction and the associated fitting problem is negligible purely from perturbative arguments. Even a conservative calculation indicates it is significant on BAO scales [18]. That means it cannot be ignored in precision cosmology.

D. Open issues and future directions
Perhaps we should ask ourselves the following questions, which may be fruitful avenues for further research: The role of (∂ 2 Φ) 2 at second-order How should such terms be dealt with? Do they signify that something unphysical has been calculated, or do they signify a breakdown in perturbation theory? Do they appear in the luminosity distance at second-order and so give a large effect there? Or is it that the ensemble average of such terms is actually irrelevant? Should perturbations be smoothed on a given scale at each order before calculating the next, which can then go to smaller scales?
Higher-order perturbation theory Given the uncertainty in guessing what higher-order perturbation theory might look like, it is timely to calculate directly, and to see if the perturbative expansion remains small, and effectively truncates at secondorder, or if it does indeed diverge and require renormalisation methods.
The validity of the coordinates How large a domain can be covered by the perturbative coordinates in a spacetime with major voids -which are genuine vacuum and so static -combined with a crowd of near-delta-functions for the density wherever there is any structure? It is crucial to the usual arguments.
Description of matter as dust Cosmology doesn't deal with many well defined particles. Smoothing over 10 ∼70 CDM particles is well defined in kinetic theory to give a fluid description, but once a structure freezes out of the expansion, it becomes our 'point particle', in an ill defined sense. Then we have 10 ∼12 virialised structures -significantly less than Avagadro's number -which are large compared to the Hubble volume, and large compared to the total mass, compared to normal fluids. Are these correctly modelled as dust? Within perturbation theory it appears consistent [18]; can this be made more general [133]? How does this relate to the radiation-reaction problem [134]?
Relativistic simulations Is a fully general relativistic simulation out of the question? This would have to be attempted very differently from current N-body simulations, but could be attempted on the expectation of improving computing power. Rather than work with point particles, one could start with the ADM equations with different fluids, plus perturbations laid down early. A simulation could consist of perturbative modes up to a maximum wavenumber, rather than N particles of a given mass. The resolution would be determined by the maximum wavenumber, so simulations would start big and get smaller as computing power improved.
Averaging observables The role of observables and the lightcone is still very unclear in the averaging scenario. Observable quantities such as distances are usually implicitly sky-averaged quantities. While it is not clear what an averaged lightcone might mean, sky-averages of observables such as the redshift-distance relation have a reasonably clear meaning, even in a very non-homogeneous universe. How does this fit into other approaches?
Covariant formulation The present detailed analyses are highly gauge dependent. Initial work on gauge invariant formalism is interesting [36,74] but ideally one would use a 1+3 covariant and gauge invariant geometric formalism for the purpose. Alternatively, a tetrad approach may be useful as all the Ricci rotation coefficients, Weyl tensors and so on are then scalar invariants. How does one average the tetrad?

VI. CONCLUSION
Averaging, even though it implicitly defines the standard concordance cosmology, is not understood in general relativity. Because extra terms in the field equations seem inevitable, then it is surely timely to understand how to deal with it -or if we can't, then decide what limitations that imposes on cosmology.
We have seen in this article some of the large variety of ways that this problem is being investigated. Some argue it's utterly irrelevant, suppressed by orders of magnitude below the background model which consequently holds genuine physical significance -analogous to being able to say "The Earth is approximately a sphere". Others take the opposite view and contend that it holds the key to the whole dark energy issue. More moderate views prevail between these two extremes, pointing out there may be significant effects in intermediate scales of cosmological interest. It is also fair to point to studies where simplified but fully non-linear spacetimes have observational properties significantly different from their spatially averaged counterparts [135,136]. And it is not clear that perturbation theory is genuinely convergent, and that it is a well behaved process. Is it properly converged at second order, where claims for tiny backreaction are made? Why are we so lucky that the equality scale is so large, suppressing what could be a huge backreaction effect? Why, in short, is cosmology so easy?
At the very least, then, these considerations surely tell us that it is important to understand the averaging, fitting and backreaction problems to see what effects there may be on cosmology. There are some scales where backreaction may be important -probably not the largest scales relevant to the cosmic acceleration, but others where precision cosmology is significant. In investigating this, we must get a clearer distinction between dynamical and observational effects -the latter not covered here, but certainly relevant to null fitting, which is the core of observational cosmology.