Modelling the emergence of cosmic anisotropy from non-linear structures

Astronomical observations suggest that the Universe may be anisotropic on the largest scales. In order to model this situation, we develop a new approach to cosmology that allows for large-scale anisotropy to emerge from the growth of non-linear structure. This is achieved by decomposing all relevant fields with respect to a preferred space-like direction, and then averaging the resulting scalar quantities over spatial domains. Our approach allows us to derive a set of large-scale effective field equations that govern the dynamics of any emergent large-scale anisotropy, and which (up to back-reaction terms) take the form of the field equations of the locally rotationally symmetric Bianchi cosmologies. We apply our approach to the dust-filled Farnsworth solutions, which are an interesting set of exact cosmological models that allow for both anisotropic expansion and large-scale bulk flow.


Introduction
The standard paradigm for interpreting cosmological observations treats the Universe as being homogeneous and isotropic on large scales. This is justified in part by the near perfect isotropy of the CMB [1], but has led to a cosmological model that appears to be in tension due to the nonvanishing of some direction-dependent observables [2] as well as due to self-inconsistencies [3][4][5][6][7]. In particular, the values of cosmological parameters inferred from the CMB appear to show some directional dependence [8,9], an apparent hemi-spherical power asymmetry [10,11] and the alignment of low-l multipoles [12,13]. In the more nearby Universe there are asymmetries in galaxy cluster scaling relationships [14], observations of radio galaxies [15], quasars [16] and perhaps even supernovae [17]. Beyond this, there is a long history of larger-than-expected bulk flows on large scales, with controversial results on a "dark flow" [18] through to more recent results [19].
These results are certainly not without controversy, sometimes with regard to the specifics of the approach used when inferring them, and often with regard to the statistical significance of finding what appear to be anomalies in large and varied collections of data sets. We do not aim to argue for or against the veracity or significance of any of these results, but merely to note that if there exist data that suggest the Universe may not be isotropic on large scales, then it should also be of interest to know how large-scale anisotropy may or may not be able to emerge within an inhomogeneous universe that contains non-linear structure, but which itself started off close to homogeneous and isotropic at early times. This is our present purpose.

arXiv:2302.05715v2 [gr-qc] 23 Jun 2023
From a theoretical perspective, the possible emergence of anisotropy at late times is itself highly intriguing. It was suggested long ago by Ellis & Baldwin [20] that any misalignment between the frame in which the galaxy/quasar distribution is at rest, and the frame in which the dipole of the CMB vanishes, could be used as a test of large-scale isotropy. Indeed, the idea that observers may be "tilted" with respect to a preferred foliation of homogeneous spaces has itself recently been the subject of much study [21][22][23][24][25], with suggestions that the consequences for some observables might even be large enough to nullify the evidence for dark energy entirely [26][27][28]. These proposals have so far gone without a concrete mechanism by which such large-scale misalignments could be generated. We take this as further motivation for the development of a framework for studying emergent anisotropy, which could then be used to evaluate the possibility of these scenarios being realised.
An important aspect in our approach is the concept of "emergence", which is intended to describe the possibility of the properties in which we are interested being realised on average, from a space-time which is highly inhomogeneous on small scales. This concept is vital for the goal we have outlined so far, as anisotropy in exactly homogeneous cosmologies typically decays rapidly in time [29], meaning that any anisotropy in the late universe would mean enormous anisotropy at early times. This is clearly not physically realistic, and so the only other possibility would be to have large-scale anisotropy realised as an emergent property of a more complicated inhomogeneous space-time. We expect that anisotropy generated in this way should weaken the tight bounds that could otherwise be imposed from the CMB [30][31][32], and might even allow for the intruiging fits of Bianchi VII h templates to the observational data to be reinterpreted [33][34][35][36].
The approach we will use in our construction will follow the general philosophy of the scalar averaging formalism developed by Buchert [37], which works by covariantly averaging locallydefined scalars over three-dimensional spatial volumes. This formalism is designed to include the consequences of the non-commutativity of averaging and evolution under Einstein's equations, and results in a set of effective field equations that govern the behaviour of large-scale averages. These equations take the same form as the locally defined versions that are more frequently used in the standard approach to cosmological modelling, but with "back-reaction" terms that account for the effects of non-linear structure on the cosmological dynamics. This is exactly what we require for our current purpose, although we will need to extend and generalize Buchert's approach in a number of regards in order to make it suitable for the task at hand.
First, if we are to model situations in which matter is to have a bulk motion then we will need to perform our averaging procedure on a set of hypersurfaces that are not necessarily orthogonal to the flow of matter [38]. This will require generalizing the standard approach, in which the averaging surfaces are usually chosen such that the matter is comoving. Second, we will need to introduce a preferred spatial direction, and covariantly decompose all geometric, kinematic and matter field quantities with respect to it. This will allow us to average all resultant scalar quantities on our generalized foliation, and identify back-reaction terms by comparing to the field equations of the locally rotationally symmetric (LRS) Bianchi or Kantowski-Sachs cosmologies. Our many back-reaction terms will generalize the single term that occurs when considering only the isotropic part of the expansion, and will allow us to study if and how anisotropic cosmological expansion and bulk flows could emerge from inhomogeneous space-times. This paper is arranged as follows: in Section 2 we introduce our mathematical formalism, then in Section 3 we discuss how this formalism can be applied. The resulting equations are presented and discussed in Section 4, before being applied to an illustrative example in Section 5. Finally, we offer some concluding remarks in Section 6. We choose units such that c = 8πG = 1 throughout, and use metric signature (−, +, +, +). Bold typeface is used to denote vectors.

Mathematical formalism
Let us now introduce the mathematical formalism that will be necessary to identify all of the required kinematic, curvature and matter quantities, to relate these quantities to those that would be measured by a set of observers, and to construct suitable, well-defined spatial averages of those local quantities, from which our concept of an emergent cosmological model can be realised.

The 1+1+2 covariant decomposition
In order to describe a cosmology with a preferred space-like direction, such as the symmetry axis of an LRS model, it is useful to perform a semi-tetrad decomposition of the field equations with respect to both a preferred time-like vector n and a preferred space-like vector m (satisfying m a n a = 0, n a n a = −1, and m a m a = 1). One can then decompose all vector and tensor quantities with respect to n and m, while computing the Ricci and Bianchi identities for n and m gives rise to a set of equations which is physically equivalent to Einstein's equations. This 1 + 1 + 2 formulation was first introduced by Greenberg [39], developed in Refs. [40][41][42][43][44][45][46], and worked out in full by Clarkson [47] and Keresztes et al. [48]. We will introduce the salient features of this approach here, and explain why they are of use for our approach.
Let us suppose that we have a time-like vector n. One can then define the projection tensor f ab = g ab + n a n b , as well as derivative operators along and orthogonal to n: with obvious generalizations to tensors of other ranks. These allow us to decompose the covariant derivative of n into kinematic quantities as [49] ∇ a n b = −n aṅb + 1 3 whereṅ a = n b ∇ b n a , Θ = D a n a , σ ab = D a n b and ω a = 1 2 η abc D b n c . These are, respectively, the acceleration of observers with 4-velocity n, the isotropic expansion of the space-like 3-surfaces orthogonal to n, and the volume-preserving shear and vorticity. In writing these expressions down, we have introduced angled brackets around indices to denote a tensor that is projected orthogonally to n and made symmetric and trace-free, such that e.g.
We have also used η abc ≡ n d η dabc as the alternating tensor on the 3-surfaces orthogonal to n.
Let us now further suppose that we have a space-like vector m, which is orthogonal to n. The covariant derivative of m can be decomposed as and a projection tensor onto the two-dimensional "screen spaces" orthogonal to both n and m can be defined as It can be seen that the projected time and space derivatives of m can be written aṡ where ab = η abc m c is the alternating tensor on the screen space. The new quantities introduced in Eq. (6) are α a = M ac n b ∇ b m c and A = −n aṁ a = m aṅ a , which are the projections ofṁ a orthogonal and parallel to n, and which are the non-geodesy of m, and the expansion, shear and twist of the screen space along m. We use the names "expansion" and "shear" at this point in analogy with the quantites defined for n, despite m being a space-like vector. This is done for linguistic ease.
To go further, we note that all 3-vectors and 3-tensors defined with respect to n can be further projected parallel and orthogonal to m.
where V = v b m b is the projection of v a parallel to m a , and V a = M b a v b is the projection of v a into the screen spaces. Similarly, for projected, symmetric and tracefree 3-tensors t ab = t ab , we can write where are the projections of t ab twice parallel, once parallel and once orthogonal, and twice orthogonal to m, respectively.
Using this approach, we can decompose the acceleration and vorticity 3-vectors associated with n in the following way: n a = Am a + A a and ω a = Ωm a + Ω a .
The only other kinematic quantity is the shear of n, which can be decomposed as The quantities A, Ω and Σ are covariantly defined scalars, while A a , Ω a , Σ a and Σ ab are covariantly defined 2-vectors and tensors defined in the screen spaces orthogonal to n and m.
We have so far considered the decompositions of the derivatives of n and m. For a complete system we will also require quantities associated with the Ricci and Weyl curvature of the spacetime. The former of these can be related to the matter content of the space-time using Einstein's equations, and therefore written in terms of a stress-energy tensor that can be decomposed with respect to n as T ab = µn a n b + pf ab + 2q (a n b) + π ab , where µ = T ab n a n b , p = 1 3 T ab f ab , q a = −T bc n b f c a , and π ab = T ab are the energy density, isotropic pressure, momentum density and anisotropic stress, respectively. The last two of these can be further decomposed with respect to m to get where all quantities are defined as in Eqs. (7) and (8), above. Finally, the Weyl tensor can be decomposed with respect to n into an "electric" and a "magnetic" part, as E ab = C acbd n c n d and which can be further decomposed with respect to the space-like vector m as Again, all quantities in these expressions are defined as in Eqs. (7) and (8). This completes the decomposition of all required quantities with respect to both n and m. The full and irreducible set of quantities describing the kinematics of n and m, the matter variables and the Weyl curvature is therefore given as follows: The reader is referred to Refs. [47] and [48] for the further details of these decompositions, and to Appendix A for a full list of all quantities defined above. For the remainder of this paper we will take Ω = 0 = Ω a , so that n is hypersurface-orthogonal. This is required for having a foliation of space-like 3-surfaces over which averaging can be performed. It is not a restriction on the properties of matter.

Changing foliation
We have so far not specified anything about the time-like vector n, other than it being hypersurfaceforming. In the sections that follow we will consider specific choices of this vector, which are orthogonal to foliations with certain special properties. Here, we will simply note that it will be of interest to be able to transform between foliations, and present some of the mathematics that is required to do so. In order to understand this, let us consider a boost from n to some second time-like vector u. The components of these vectors are related by and where Here the boost vectors v a and w a exist respectively in the surfaces orthogonal to n and u, such that n a v a = u a w a = 0. The projection tensor into the hypersurfaces orthogonal to u (if they exist) can then be defined as h ab = g ab +u a u b , which allows us to write w a = γ −1 h b a v b . The same projection tensor also allows us to define a new vectorm with components where k = 1 + (γm c v c ) 2 −1/2 andk = 1 + (γm c w c ) 2 −1/2 . These are the components of the preferred space-like vector m projected orthogonal to u, as illustrated in Figure 1. The whole system of 1 + 1 + 2-equations from Section 2.1 can now be set up with respect to either n or u, with the preferred space-like vector projected according to Eq. (18). After we boost from n to u, we wish to express the scalars in the new foliation in terms of those in the old one, plus any terms arising from the boost 3-velocity v that connects them. If |v| 1, which is expected if, for example, n corresponds to the canonical foliation of FLRW spacetime, and u is the matter 4-velocity in that spacetime, then one can expand these expressions order by order in v, leading to vast simplifications. However, in the general case, one need not make this approximation. Performing the full calculations, one finds the following results: The isotropic expansionΘ = h ab ∇ a u b is related to Θ = f ab ∇ a n b bỹ where κ = D a n a and β ab = D a v b . This equation shows how the expansion of the rest spaces defined with respect either n or u are related to each other, and explicitly provides the additional terms that contribute to the expansion after performing a boost.
and the scalar acceleration,Ã =m a u b ∇ a u b , and the scalar vorticityΩ = 1 2m a η dabc u d ∇ b u c , arẽ where we have used that Ω = 1 2 m a η dabc n d ∇ b n c = 0. One sees immediately that if the vorticity W a associated with the relative velocity v vanishes, thenΩ = 0 as well.
Meanwhile, the twist and expansion of the screen spaces alongm,ξ = 1 2˜ abM acM bdD cmd and φ =M abD amb , are related to ξ and φ bỹ For the matter variables, the energy densityμ = T ab u a u b , isotropic pressurep = 1 3 T ab h ab , and scalar heat fluxQ = −T bc u b h c am a are given bỹ and the scalar anisotropic Finally, one can compute the scalar parts of the Weyl curvature for the new congruence. For the scalar parts of the electric and magnetic Weyl tensors this gives [50] This completes the transformation rules for all of our covariantly defined scalars. The 2-vector and tensor quantities that result from the 1 + 1 + 2 decomposition can also be transformed between foliations, but only appear in our averaged description via "backreaction" terms. We therefore omit presenting their transformation rules here.
In what follows, we will use often use u to refer to the flow lines of an irrotational dust fluid. In the presence of such matter, this vector is uniquely defined and therefore also provides a unique foliation on which to perform averaging. This is the foliation that is most often used in the literature on this subject, and is the one chosen for the standard approach to formulating Buchert's equations [37]. We will then use n to refer to any other foliation, chosen (for example) for observational or geometric reasons. Such freedom is required if we are to allow for the possibility of bulk flows, as the 3-velocity of matter vanishes identically if we choose the foliation to be induced by the flow of dust u. It also allows us the freedom to re-foliate in situations in which the description of space-time might break down if it were specified by the fluid flow, as would happen, for example, in perturbed FLRW cosmologies that contain non-linear structures, or if one wished to consider fluids with non-zero vorticity or caustics in their flow lines.

Scalar averaging
Averaging in general relativity is a notoriously difficult problem to define, primarily as the theory is defined in terms of tensors, which cannot be straightforwardly compared at different points in a curved space-time [51]. This makes defining averages difficult, as quantities that exist at different points in an averaging domain first have to be transported to the same point for comparison, which is a non-unique and path-dependent process. The problem is further exacerbated by the non-linearity of Einstein's equations, which means that the acts of averaging and evolution do not commute. The evolution of an averaged quantity cannot therefore be guaranteed to be the same as an average taken at some future time.
The first of these problems has been the focus of a large number of different studies, using a variety of different approaches (see e.g. Refs. [52][53][54][55][56][57]). The benefits and drawbacks of each of these approaches are a subject of much debate, but in order to make progress one can consider a more limited case: the averaging of scalars. Such quantities can be readily compared at different points in averaging domains, and it is therefore much more straightforward to compute their averages. Indeed, attempts have been made to construct averaging schemes that completely characterize space-times only in terms of the averages of scalars [58]. In what follows we will proceed by constructing an averaging formalism for the scalars that result from our 1+1+2 decomposition.
We will be interested in spatial averages of scalar quantities, over the 3-dimensional spatial volumes orthogonal to the irrotational flow lines n. Following Buchert [37], we will define these as follows: where S is the scalar being averaged over the spatial domain D, which has volume V D = D dV , and where dV = d 3 x √ f is the spatial volume element. In these expressions, we have taken x i to be coordinates on the spatial hypersurface orthogonal to n, and written f = det(f ab ) where f ab = g ab + n a n b is the induced metric. This definition of averaging depends on our choice of foliation and on the selected averaging domain, but is otherwise unambiguous.
The rate of change of this average can be determined by considering that the factor √ f represents the volume of an infinitessimal region of space. Assuming the spatial coordinates are fixed by the flow of n, we then have ( √ f ) · = Θ √ f , where the dot represents the time derivative operator n a ∇ a and where Θ = ∇ a n a is the expansion scalar. This allows one to derive Buchert's commutation rule for time derivative and averaging operations [37]: where we have defined the covariance as Cov(S 1 , S 2 ) := S 1 S 2 − S 1 S 2 . This result relies on the averaging domain D being propagated by the flow of n. In what follows we will make repeated use of the covariance of scalar quantities, as well as the variance, Var(S) := Cov(S, S). In Eq. (32) it is assumed that the lapse function N , between the leaves of the chosen foliation, is independent of spatial position. This is the case if (for example) the time-like 4-vector n is identified with the 4-velocity of a field of irrotational dust u, as in this caseṅ a =u a = 0. Each leaf in the foliation can then be identified with a single value of the proper time along the world-lines of this dust, and will be the same at all points in each spatial hypersurface. More generally, we will haveṅ a = D a ln N = 0, meaning that the proper time of a family of observers following n is not constant across D. For the purposes of calculating cosmological averages it is desirable to have quantities defined over an averaging domain to have a single value for the time coordinate t. This can be achieved by writing n a = −N ∂ a t where N =ṫ −1 , such that a derivative with respect to t can be written as ∂ t S = NṠ. The commutation rule (32) then becomes which allows us to calculate the evolution of scalar averages as a function of the coordinate time t, once the lapse function N has been determined from the acceleration of n.
One can apply the averaging operation from Eq. (31) to each of the terms in any scalar equation that is defined locally at each point in space-time. After suitable application of the commutation rule from Eq. (32) or (33), one can then obtain from that averaged scalar equation a new one that is written entirely in terms averages and their derivatives with respect to time. Equations derived in this way can be thought of as governing the behaviour of averages, rather than the behaviour of locally-defined quantities (as usually occurs in theories of gravity). This is ideal for describing the emergent large-scale properties of a space-time, as is relevant for the study of cosmology. Indeed, the averaged equations that result from this procedure can often be written in a form identical to that of the local gravitational equations of an homogeneous cosmology, with any extra terms that occur being considered to describe the "back-reaction" effect that results from averaging small-scale structures.
In the standard approach to this problem, three scalar equations are usually considered: the Raychaudhuri equation for the evolution of the expansion, the Hamiltonian constraint equation that relates expansion and energy density, and a contracted Bianchi identity that corresponds to energy conservation. The averages of these are all evaluated in the foliation orthogonal to the flow of dust, and a set of Friedmann-like equations emerge that govern the average expansion [37]. If one wishes to consider the emergence of any anisotropic quantity, however, then this necessarily requires further scalar equations, such as those governing the scalar shear Σ and momentum density Q. This is what we will pursue in the following sections. We note that an interesting step in this direction has already been taken by extending the Buchert formalism to include an evolution equation for the magnitude of the shear, σ 2 = 1 2 σ ab σ ab [59]. We extend this result further by deriving the full set of averaged scalar equations in the 1 + 1 + 2 formalism. We find a set of equations that can be used to interpret that large-scale cosmological average in terms of an emergent LRS Bianchi model with tilt, and with back-reaction terms that give the effect of small-scale structure on all emergent quantities.

Application of the formalism
Let us now consider how to apply the formalism from Section 2 to a cosmological space-time. This will require a choice of the time-like and space-like vectors n and m, and some consideration of how to understand and interpret the quantities that result.

Choice of foliation
An application of the averaging procedure outlined above requires a choice of foliation, or equivalently a choice of the irrotational time-like 4-vector n. The averages obtained using Eq. (31) will then correspond to the large-scale properties of the 3-dimensional spaces that constitute the leaves of this foliation. Different foliations will mean that one is considering different 3-dimensional spaces, and hence different averages will be obtained. It is therefore necessary to make sure an appropriate choice of foliation is made, for the situation being considered. This is given further importance by the fact that observers in different frames will infer different cosmological parameters from the Hubble flow around them [21][22][23][24][25][26][27][28].
In general, one might be interested in choosing a foliation that is expected to give results that can be associated with a particular observable [60], that have a particular mathematical or physical meaning associated with them [61], or that are perhaps convenient in some other way [38]. For example, one may wish to construct Hubble diagrams in a frame comoving with the flow of matter [37], or in a frame of 'most uniform Hubble flow' [62][63][64]. Of course, in order to relate observables corresponding to quantities calculated on different foliations one will need to be able to transform between frames ‡, as considered in Section 2.2.
Let us now consider some specific choices of foliation that one might use: The comoving foliation exists when the Universe is filled with irrotational dust, and n is chosen to be coincident with the 4-velocity of that dust, u. This choice shares some properties with the "comoving synchronous" gauge that may be familiar from studies of perturbed Friedmann-Lemaître-Robertson-Walker (FLRW) models. It has the distinct benefit of being tied to a physical quantity: the flow of matter. It therefore takes a particularly privileged position in the pantheon of possible choices, but does require the matter flow to be vorticity-free, which is not expected to hold for realistic astrophysical structures. It also cannot account for the existence of any bulk flow, as it corresponds to the choice of frame in which matter is at rest. Nevertheless, this is the standard choice of foliation in much of the literature on mathematical cosmology [49], as well as in many studies of perturbed FLRW.
The gravitational rest frame was defined in Ref. [38] to be that in which the magnetic part of the Weyl tensor vanishes: so that the emergent cosmological model is as Newtonian in character as possible [66]. Cosmologies that obey this condition have been described as "silent" [67]. Up to second order in perturbations around a Robertson-Walker geometry this choice corresponds to the Poisson gauge [38], but when applied non-perturbatively ends up being a strong restriction on the allowed solutions [68,69]. It therefore appears to be particularly useful for studies of perturbed FLRW models, as the Poisson gauge is one of the few that is expected to be valid into the regime of non-linear structure formation, where post-Newtonian expansions are required to describe weak gravitational fields [70].
The gravitational wave frame is a more conservative choice defined by the demand that H ab is divergence-free, such that This is a covariant way of stating that the only gravito-magnetic contributions to the curvature measured by an observer comoving with n would be those coming from gravitational waves [71,72]. It is therefore less restrictive than the condition H ab = 0, and typically corresponds to a congruence which is distinct from the flow of dust [73]. It is well-suited to numerical-relativistic cosmological simulations [74], but in a general space-time it is not guaranteed that the frames it picks out are hypersurface forming (though one may be able to take the irrotational part, which would be).
The constant mean curvature foliation has a long history [61], and is defined by the condition that the spatial gradient of the expansion scalar (also known as the "mean curvature") vanishes: This foliation has the property of being unique (under certain circumstances), and having a monotonic variation in the expansion scalar between leaves [75]. It therefore provides a plausible candidate for the provision of a universal, cosmological arrow of time. The literature on this particular choice of foliation has been largely resticted to the literature on mathematical cosmology, but it also exists in perturbed FLRW models under the name of "uniform expansion" gauge.
The list of choices presented above is by no means exhaustive. Other interesting choices of foliation are specified by the "zero-shear" condition, σ ab = 0. This is closely related to the gravitational rest frame, in that the Ricci identities provide a constraint that sets H ab equal to the covariant curl of the shear tensor. Vanishing of the shear therefore corresponds to the vanishing of H ab , which explains why the Poisson gauge (sometimes referred to as the "zero shear" gauge) corresponds to the gravitational rest frame. Another is the "constant density" foliation, in which the spatial gradient of the energy density vanishes, D a µ = 0. As with the comoving foliation, this has the benefit of being directly tied to the matter content of the space-time. However, it also shares the weakness that it will be highly problematic in the presence of non-linear structures.

Choice of preferred spatial direction
In order to extract scalars from anisotropic quantities, we also need to choose a preferred spacelike direction m. All 3-vector and tensor quantities can then be decomposed with respect to this direction as prescribed by Eqs. (7)- (8), and the scalar parts of the decomposed quantities can be averaged according to Eq. (31). Clearly, in order to make sense, this vector must (in some sense) point in the same direction at every point in an averaging domain. Just as in the case of the time-like vector n this presents a choice; there is in general not going to be any single uniquely preferred space-like direction m, and we will need to select a suitable way to define this direction based on the situation at hand. This choice will be important for the outcome of our averaging process.
As in the case of choosing n, we may wish to apply geometric or observational considerations when selecting a preferred direction m. In space-times with symmetry, it may be that a preferred direction is selected by the existence of some Killing vectors. For example, in LRS geometries m can be unambiguously taken to correspond to the spatial axis around which the rotational symmetry exists. This can be identified by searching for the non-degenerate eigenvector of any non-vanishing 3-tensor, such as σ ab or E ab [47,76]. Alternatively, in algebraically special geometries it may be possible to pick out a unique direction from the projections of the canonical null tetrad into our chosen foliation. These will be specified by the properties of the Weyl tensor.
An alternative method for selecting a preferred spatial direction would be to use observational methods at each point within an averaging domain. This method would be better adapted to situations in which anisotropy in a particular observable is being considered, or where the space-time has no explicit symmetries or special algebraic properties (as is the case for the real Universe). For example, one might choose to take the axis of CMB parity asymmetry [77,78], the direction of greatest asymmetry in the galaxy distribution [14,15], the quasar or Type Ia supernovae dipole [16,17], or the direction of greatest variation of a cosmological parameter or coupling constant such as H 0 or the fine structure constant α [79][80][81]. As long as these directions line up at different points in space, as one would expect in an anisotropic universe, then they should provide a well-motivated choice of preferred direction.

Interpretation as an LRS Bianchi cosmology
As well as considering the choice of foliation and space-like direction, it seems worthwhile to also consider what the model that emerges from averaging the scalars in our 1+1+2-decomposition will represent. These will be precisely the set of scalars represented in Eq. (14). It is an extended set compared to the scalars that result from a 1+3-decomposition, and include (among other things) the scalar parts of 3-vectors such as the momentum density Q, and the scalar part of 3-tensors such as the shear Σ.
The averaged scalars that result from the 1+3-decomposition, as provided in Buchert's groundbreaking approach [37], are relatively straightforward to interpret. They result in a set of equations that govern the averages of the expansion scalar, the energy density of dust, and the curvature scalar of the 3-spaces: { θ , µ , (3) R } , with everything else being collected together into "back-reaction" terms. Defining an effective scalar factor a D and an effective curvature parameter k D of the averaging domain D by the equations that result bear a striking similarity to the dust-dominated Friedmann equations that govern the behaviour of homogeneous and isotropic cosmologies [37]. This is no accident, as the set of averaged scalars listed above are exactly those required to fully characterise an homogeneous and isotropic universe. All other quantities will be 3-vectors or 3-tensors, but these must vanish in geometries that are isotropic around every point in space. So, after discarding such quantities in the scalar averaging approach, we are left with exactly the set that is required to describe an FLRW model. One could say that the geometry has been averaged to a Friedmann cosmology.
If we now consider the scalars that result from the 1+1+2-decomposition, as given in Eq. (14), then we have a more complicated situation. It can, however, be understood in terms of a well studied class of cosmological models: the Locally Rotationally Symmetric (LRS) models pioneered by Ellis [82,83]. A space-time is locally rotationally symmetric if it admits everywhere a onedimensional continuous isotropy group, which essentially corresponds to rotations about a local axis of symmetry m. The existence of such a symmetry means that all 3-vectors are proportional to m, and all projected symmetric and trace-free 3-tensors are proportional to m a m b − 1 2 M ab , where M ab is defined in Eq. (5). In the language of the 1+1+2-decompostions described in Section 2.1, this means that all 2-vectors and 2-tensors must vanish. The entire dynamics of an LRS space-time is therefore described purely by the set of scalars given in Eq. (14). For this reason, we expect to recover a set of equations for our averaged scalars that can be written in a form similar to the equations that govern LRS cosmologies, with additional back-reaction terms due to the averaging of small-scale structure. We will therefore say that we are averaging to an LRS Bianchi cosmology.
The Bianchi classification contains all spatially homogeneous cosmological models (except Kantowski-Sachs) and divides them into classes labelled I to IX. The classes that admit LRS symmetry are I, II, III, V, VII and IX. Of these, the tilted LRS cosmologies can be found in classes V and VII h , and the fully isotropic cosmologies can be found in classes I, V, VII and IX. Due to their homogeneity, the only non-zero derivative of any of the scalars S within these models will be a time derivatice. All the governing equations will then necessarily have the form where S i labels the allowed scalars from Eq. (14), b ij and c ijk are constants, and a i = 0 or 1 depending on whether the equation is a constraint or an evolution equation. The exponent p i of the lapse function N depends on the scalar being considered: for the kinematic scalars {Θ, Σ, A, φ, ξ} we have p i = 1, and for the matter and Weyl curvature scalars, {µ, p, Q, Π, E, H} we have p i = 2. These equations are a significant simplification of the full set of equations of motion, and will be the form we expect for the equations that govern the average of the scalars from our 1+1+2-decompostion. The only difference will be that we will find back-reaction terms on the right-hand side of each equation, instead of zero.

Averaged equations for anisotropic cosmologies
In this section we will present the explicit form of the equations governing the averaged scalars from Eq. (14). These equations are valid for any hypersurface-orthogonal n, and any choice of m.

The Ricci identities for n
The Ricci identities taken with respect the time-like vector n can be written as Taking the trace of S abc over its first and third indices, and contracting with n, gives the Raychaudhuri equation (n b S a ba = 0): where ∂ t S = NṠ andṠ = n a ∇ a S. This is an evolution equation for the expansion, written in terms of coordinate time, and is identical to the corresponding equation from LRS Bianchi cosmologies up to the back-reaction term B 1 , and the presence of the lapse function N , which cannot in general be set to unity. We find that the back-reaction term can be written as which encodes all information about the influence of small-scale inhomogeneities on the acceleration of the expansion of space. This is the only equation that can be derived from S abc that has a counterpart in the standard approach of averaging to isotropic cosmology. As this is an equation for the acceleration of expansion, a sufficiently large (and positive) B 1 would lead to an accelerating universe, at least within the spatial domain being considered. By again taking the trace of S abc over the first and last indices, but this time contracting the middle index with respect to the space-like vector m we can obtain a momentum constraint equation (f b a S c bc m a = 0): where the back-reaction term is given by This term describes the direct contribution from inhomogeneity to the large-scale momentum density in the direction of m.
There are two remaining equations that can be derived from S abc . The first of these is an evolution equation for the scalar part of the shear (f d a f e b n c S dec m a m b = 0): where the back-reaction term is given by This drives the generation of anisotropy in the expansion of space. The last equation to be derived from S abc is a constraint that is usually used to set either the acceleration of n (i.e. A) or the twist of m (i.e. ξ) to zero. In this case, however, we find (f d a f e b n c S dec η abf m f = 0): where the back-reaction term is given by (N A, N ξ), In this case we see that any non-zero effect from the small-scale inhomogeneity will lead to a more complicated result, with N A = 0 = N ξ .

The Ricci identities for m
We can also consider the Ricci identities defined with respect to the preferred space-like direction m: This time, the only equation that has a counterpart in isotropic cosmology is obtained by contracting over the last two indices with the screen space projection tensor, and contracting the first index with m. This gives the Hamiltonian constraint equation (m a M bc R abc = 0): The back-reaction term in this case is given by the following expression: This term acts in the same way as an additional effective energy density, in what would reduce to the Friedmann equation in isotropic cosmology. Comparing (39) to the trace of the Gauss embedding equation, (3) R = 2µ − 2 3 Θ 2 + 2σ 2 + 2Λ, and using σ 2 = 3 4 Σ 2 + Σ a Σ a + 1 2 Σ ab Σ ab , one finds that the average (lapse-weighted) Ricci curvature of the space-like hypersurfaces is given by The second equation to be gained from contracting R abc with the screen space projector is to be found using the time-like vector n, and is an evolution for the expansion of m in the preferred space-like direction (n a M bc R abc = 0): where we find This equation has no counterpart in isotropic cosmological modelling. The first of the three further equations that can be derived from R abc are a constraint on the magnetic part of the Weyl tensor ( ab n c R abc = 0): where we find For our averaged cosmology to be "silent" we therefore require B 7 to vanish along with either N ξ or N Σ , or to conspire to cancel 3 N ξ N Σ . The second equation is an evolution equation for the twist of m (n a bc R abc = 0): with back-reaction term The final equation involving R abc is a contraint on the product of the expansion of m in the preferred space-like direction and its twist (m a bc R abc = 0): The back-reaction term in this equation is As in Eq. (38), the back-reaction term in this last equation prevents one obtaining the usual result that either N φ or N ξ vanishes, with N φ = 0 = N ξ being required if B 9 = 0.

The Bianchi identities
Our final set of equations are derived from the Bianchi identities: The only equation that could be derived from these identities in the case of isotropic cosmologies would be the energy conservation equation, which in our case takes the form (n a W bc abc = 0): where the back-reaction term is given by The back-reaction term here provides an additional source that drives the rate of change of the average energy density, N 2 µ . All other equations correspond to anisotropic cosmologies only, and vanish in the isotropic limit. It can be seen that if one chooses the foliation such that n is the 4-velocity of pressureless dust, then B 10 = 0. For a generic foliation, however, this need not be the case, as an observer that is not comoving with the dust will typically measure non-zero pressure, momentum density and anisotropic stress. The first of the remaining equations gives conservation of momentum density (m a W bc abc = 0): This equation contains the back-reaction term This term drives the evolution of the momentum density, and if non-zero can therefore source a bulk flow in the averaged cosmology.
There are four remaining equations to be derived from W abcde . Two of these are evolution equations for the scalar parts of the electric and magnetic Weyl tensor, and two are constraints for these quantities. The evolution equation for the electric part of the Weyl tensor is given by (f de n c W cdabe m a m b = 0): with back-reaction The evolution equation for the magnetic part of the Weyl tensor is (η cd a n e n f W ecdf b m a m b = 0): which has back-reaction term The remaining two equations are a constraint equation for the electric Weyl tensor (f b a n c n d W e db ec m a = 0): which we find to have a back-reaction term given by and a constraint for the magnetic part of the Weyl tensor (η bc a n d W e db ec m a = 0): The final back-reaction term for this equation is given by Eqs. (35)-(50) provide a complete set that can be used to describe the large-scale properties of an anisotropic cosmological model after averaging. These equations are a considerable complication on those that govern exact LRS cosmologies, due to the presense of the back-reaction terms, B i . If all of these terms are sufficiently small when calculated for a domain D, then the expansion of D should be expected to be well approximated by an exact LRS Bianchi model. If this is not the case, then we will be in a situation where back-reaction of small-scale structures on the large-scale properties of the cosmology is no longer negligible. The reader may note that the equations above are all written in terms of the coordinate time t, which (by definition) takes the same value at every spatial averaging domain D, as required. In the next section we will consider a family of exact cosmological models, to help develop our understanding of these terms by explicitly calculating them for an example geometry.

Example: Farnsworth dust solutions
In order to understand our formalism, it is illustrative to apply it to a class of exact cosmological models. For this we choose the anisotropic cosmologies found by Farnsworth [84]. These are exact solutions to Einstein's equations. They are of Bianchi type V , admitting a four-parameter group of isometries, including local rotational symmetry. The three-dimensional space-like surfaces of transitivity of this isometry group are in general not coincident with the hypersurfaces orthogonal to the dust 4-velocity: these are tilted cosmologies. The k = −1 FLRW metric is contained within this wider class of metrics as a special case, and although they do not contain anything that could be considered as non-linear structure, they do provide us with a precise example geometry to illustrate the application of our formalism.
The metric for the Farnworth solutions can be written as [84] where the functional dependence of the metric functions X and Y has been fixed by the presence of the Killing vectors and where C is a constant. The 1-dimensional isotropy group can be seen to be generated by X 4 . The dust 4-velocity can be written in these coordinates as u = ∂ t , for which the corresponding axis of rotational symmetry is given bym = X −1 ∂ r . Substituting the metric from Eq. (51) into Einstein's equations, one obtains a relationship between X and Y , where a prime denotes a derivative with respect to t + Cr and k is a positive constant, as well as the following Friedmann-like equation: where D is a positive constant related to the energy density and where we have set Λ = 0. The solution to Eq. (54) can be written in parametric form as where we have defined W := D/3k 3 . The constant C controls the size of the bulk flow, and in the limit C → 0 we recover FLRW and η becomes the usual conformal time coordinate. In what follows, when required, we will choose C = 2 and W = 125 in order to display numerical results.
Let us now consider re-foliating this space-time by boosting along the symmetry axis from the matter frame u into some new frame n, such that and γ = (1 − v 2 ) −1/2 takes its usual form in terms of the 3-velocity v. The functional dependence of v is motivated by the symmetries of the space-time, and it can be seen that this choice of boost means that all 2-vectors and tensors vanish in the new frame. This will greatly simplify the calculation of our backreaction scalars. After the boost, the kinematic scalars can be written as while the matter scalars become where the rest mass density is given by ρ = G ab u a u b . As the space-time is of Petrov type D, the electric Weyl scalar E is invariant under boosts, and takes the following form for all possible v: With these results in place, we can now consider re-foliations of the Farnsworth cosmology, and calculate the 2-scalars and their associated backreaction terms for any of them.

Homogeneous foliation
Let us first consider the foliation composed of leaves that coincide with the surfaces of transitivity of the Killing vectors from Eq. (52). These are spatially homogeneous 3-dimensional spaces with a time-like normal n that must satisfy n a X a 3 = 0. This immediately implies v = −C/X, and results in a set of spaces of constant η. The set of non-vanishing scalars in this foliation is larger than in the foliation orthogonal to the matter flow, as it contains a bulk flow that gives rise to non-zero pressure, heat flux and anisotropic stress. The full set of scalars is {Θ, Σ, E, φ, µ, p, Q, Π}, which can be immediately obtained from Eqs. (58)- (62). For each of these scalars (S) we also know that their projected covariant derivative (D a S) must vanish, due to homogeneity. This immediately implies that this foliation is both constant-density and constant-mean-curvature.
The line-element in Eq. (51) can be recast into coordinates adapted to this foliation by making the transformation [85] t = Cx + dT β(T ) and where α = C/(X 2 − C 2 ) and β 2 = 1 − C 2 /(X 2 − C 2 ). This results in where X 2 = C 2 + A and Y 2 = Be 4x exp 2 dT α/β . All quantities that were functions of t + Cr are now purely functions of T , and moreover we have that n = ∂ T and X 3 = −∂ x + y∂ y + z∂ z . This makes clear the orthogonality of n with the homogeneous spaces spanned by the Killing vectors from Eq. (52). The lapse function is equal to unity in these coordinates. The behaviour of the expansion rates of space, and the expansion-normalised variables Σ (n) , φ (n) , E, µ (n) , are shown in Figures 2 and 3. Figure 2 shows the expansion rates parallel and orthogonal to the direction picked out by the bulk flow, and displays significant anisotropy at early times (small η). As the bulk flow decays at late times (large η), the difference in expansion rates in different directions also decays, demonstrating that space isotropises in this foliation. Figure 3 shows the other kinematic scalars as a function of η, normalised by the appropriate power of the expansion rate Θ to make the resultant quantity dimensionless. One sees shear Σ and and electric Weyl curvature E both increase at early times, consistent with growing anisotropy in this limit. In both plots η is made dimensionless by normalising it relative to η eq , its value  when Y = −D/3k 2 (such that the two terms on the right-hand side of Eq. (54) are equal). In Figure 2 the expansion rates are made dimensionless by normalising them with respect to the mean expansion rateH (n) = 1 3 Θ (n) at the equality time η eq . In this foliation it can be seen that all of back-reaction scalars B i from Section 3.3 must be equal to zero. This happens because all variances and covariances must vanish (as D d 3 x √ f S(T ) = S(T )V D ), and because all terms containing projected covariant derivatives (such as m a D a S) must also vanish. The cosmology one observes in this foliation can therefore be entirely described by the locally defined scalar quanitites, without any need for averaging. This is as exactly as expected for a space-time that admits a homogeneous foliation. Of course, in most geometries such a simplifying set of symmetries will not exist, and in more realistic cosmologies one would need to perform explicit averages in order to gain a set of quantities that could be used to describe the large-scale properties of space. Let us now consider an inhomogeneous foliation of the Farnsworth solutions, to show how this would work in this simple example space-time.

Matter-rest-space foliation
Let us now foliate the exact same space-time into hypersurfaces orthogonal to the matter flow u, with induced metric h ab = g ab + u a u b . This corresponds to v = 0 in Eqs. (58)- (62), and to surfaces of constant coordinate time t. The lapse is again equal to unity, so covariant and coordinate time derivatives are equivalent. The only non-zero quantities in this case are the scalars {Θ, Σ, E, µ, φ}, which in terms of the parameters from the solution in Eq. (55) are and the dust density The electric Weyl curvature is given by Eq. (63). The normalised expansion rates as a function of η are shown in Figure 4, and the other expansion-normalised scalars are shown in Figure 5. These plots, displaying quantities calculated in the foliation generated by the matter flow u, may be directly compared with Figures 2 and 3.
In order to average these quantities we need an averaging domain D on each of the surfaces orthogonal to u. If we take this domain to have coordinate extents {∆y, ∆z, ∆r}, then the spatial volume is given by where J(η) = (cosh η − 1) 2 W (cosh η − 1) 2 − 2C sinh η exp {−W (sinh η − η) /C}. We can then where ∆r = r max − r min . In choosing the integration domain (i.e. r min/max ), it is important to avoid the big bang singularity that is located at r BB (t) such that t + Cr BB (t) = 0. For any C > 0 this is satisfied by r min > 0. The value of r max needs only to be sufficiently large that η min and η max are relatively far apart, but once this is ensured to be true the backreaction effects are only weakly dependent on its value. If we choose r min = 1, r max = 20 and C = 2 as an example, then one obtains the integration domain shown in Figure 6. At early times the distance from r min to r max corresponds to a large separation in η, while at late times the integration domain shrinks (as the space isotropises). The evolution with t of our averaged scalars, normalised by the averaged expansion Θ , is shown in Figure 7. The inhomogeneity in these spaces has been averaged out to give a set of purely time-dependent scalar functions. These averages provide a cosmological description of the large-scale properties of the spaces orthogonal to the matter flow u, including anisotropy. The magnitude of the back-reaction scalars that produce them are displayed in Figure 8, where we focus on {B 1 , B 3 , B 6 , B 12 } as they appear in evolution rather than constraint equations. These quantities  are shown as ratios of the time derivative of the relevant average, as per the equations in Section 4, in order to show the fractional contribution they make to their evolution. In each case, the averaging domain used is the one displayed in Figure 6, and we have again used {C, W } = {2, 125}.
It can be seen from Figure 8 that the back-reaction effect becomes large at early times, and dissipates to zero at late times. This is as expected in the Farnsworth solution, as anisotropy (and hence inhomogeneity in the matter-rest-space foliation) is entirely due to the bulk flow v = −C/X, which becomes small as the scale factor X grows and the space isotropises. This example demonstrates that the large-scale averaged properties of a space depend on the chosen foliation, and that one requires in addition a choice of averaging domain (and in general a spatial direction m). The kinematic and matter quantities that result can display different behaviours when calculated on different foliations, and so can the magnitude of back-reaction scalars. This clearly demonstrates the need for these choices to be made carefully, and for the interpretation of scalar averaging results in cosmology to be understood in terms of these choices.

Discussion
We have introduced an approach for modelling the emergence of cosmic anisotropy from inhomogeneous space-times. This has been achieved by performing a 1 + 1 + 2-decomposition of all relevant fields [47], and averaging the covariantly defined scalars that result [37]. The equations that result describe a locally rotationally symmetric Bianchi cosmological model, with additional source terms due to the back-reaction from inhomogeneity. Our approach is very general, allowing for averaging surfaces and the direction of anisotropy to be freely chosen, and does not require any special properties from the matter content of the space-time.One can freely convert between different choices of foliation, according to the transformation rules set out in Section 2.2, in order to accommodate a variety of different types of cosmological model. A set of models of particular interest will be near-FLRW geometries with small tilt (i.e. with v c), in which case one may expect our general expressions to simplify considerably.
We expect our formalism to be of use for modelling situations in which astronomical observations suggest large-scale cosmological anisotropy [2], as well as for assessing the significance of potential inconsistencies between different data sets [3][4][5]. In the former case, the existence of large-scale anisotropy has been claimed by some to be enough to nullify at least some of the evidence for dark energy [26][27][28], while in the latter case the Hubble tension points to a possible failure of a single FLRW model to consistently account for all cosmological observables [94,95]. Our formalism provides a way to model deviations from FLRW that may be of interest for understanding these observational anomalies, within a fully relativistic and covariant framework, and without introducing any exotic new physics. It remains to be seen whether any particular model will be able to produce an emergent anisotropy that would be compatible with any of the anomalies claimed in the literature, or whether any of the additional degrees of freedom in this approach would be sufficient to alleviate any tensions. Our framework provides a mechanism by which such questions can be investigated.
The set of equations presented in Section 4 for our emergent scalars is, in general, not closed, although there do exist integrability conditions that relate the backreaction scalars B i to one another. This non-closure is to be expected on mathematical grounds, as they do not contain evolution or constraint equations for the 2-vectors or tensors that result from the 1 + 1 + 2decomposition. On more physical grounds, one would not expect equations that govern the largescale average of a wide class of inhomogeneous equations to be closed, as this would suggest that all inhomogeneous space-times would have the same average properties, which cannot possibly be true. In our formalism, the set of large-scale cosmological properties is dependent on the presence of non-linear structure through the B i . Once these are specified, the averaged cosmological equations become closed, and can be integrated.
The magnitude of back-reaction terms will depend on the geometry of the underlying spacetime. In the case of nearly FLRW cosmologies, which have been extensively studied, most authors find the relative size of back-reaction terms to be small in the Buchert equations of isotropic cosmology, such that B/H 2 ∼ 10 −5 (see e.g. [93]). While such contributions are indeed small, they are not as small as one may have naively assumed, given that B itself has leading-order contributions at second-order in cosmological perturbation theory. The reason why B is not even smaller is that it contains terms ∼ Φ D 2 Φ ∼ Φ δ, where δ is a density constrast that can become of order unity (or larger) in the presence of non-linear structures. In the present case we have back-reaction terms appearing not only in the averaged versions of the Friedmann equations, but also (for example) in the momentum evolution equation (46). Using the same logic, the presence of non-linear structure might be expected to result in back-reaction terms of order B 11 /H 2 ∼ 10 −5 . In this case, however, such a contribution would be of the same size as the usual terms on the lefthand side of Eq. (46), meaning that the back-reaction terms could potentially be more significant overall. Whether these expectations are, in fact, realised will require a detailed study.
In Section 5 we applied our formalism to a simple exact cosmological solution: the Farnsworth solution [84]. This explicitly demonstrated the dependence that arises from choice of foliation and averaging domain, as well as how we envisage the formalism to be applied. It shows the necessity of choosing space-like surfaces carefully when considering concepts such as the expansion of space and bulk flows in non-FLRW cosmological models, and the effect that back-reaction from inhomogeneity can have on these quantities. In future work, we will consider more general inhomogeneous cosmological models, as well as how our formalism should be applied in the real Universe. This will necessarily require developing an understanding of how to interpret observations in terms of averages, and how to most appropriately choose both the foliation and the preferred spatial direction.