Splitting the spacetime: a systematic analysis of foliation dependence in cosmic averaging

It is a fundamental unsolved question in general relativity how to unambiguously characterize the effective collective dynamics of an ensemble of fluid elements sourcing the local geometry, in the absence of exact symmetries. In a cosmological context this is sometimes referred to as the averaging problem. At the heart of this problem in relativity is the non-uniqueness of the choice of foliation within which the statistical properties of the local spacetime are quantified, which can lead to ambiguity in the formulated average theory. This has led to debate in the literature on how to best construct and view such a coarse-grained hydrodynamic theory. Here, we address this ambiguity by performing the first quantitative investigation of foliation dependence in cosmological spatial averaging. Starting from the aim of constructing slicing-independent integral functionals (volume, mass, entropy, etc.) as well as average functionals (mean density, average curvature, etc.) defined on spatial volume sections, we investigate infinitesimal foliation variations and derive results on the foliation dependence of functionals and on extremal leaves. Our results show that one may only identify fully foliation-independent integral functionals in special scenarios, requiring the existence of associated conserved currents. We then derive bounds on the foliation dependence of integral functionals for general scalar quantities under finite variations within physically motivated classes of foliations. Our findings provide tools that are useful for quantifying, eliminating or constraining the foliation dependence in cosmological averaging.


Introduction
The formulation of the theory of relativity came with the remarkable insight that proper time varies between observers, and, consequently, the Newtonian notion of a unique time parameterization of physical phenomena was abandoned.In concrete applications of relativity to the modelling of physical systems, it is nevertheless practical to introduce a time parametrization.Foliations of spacetime into a set of spatial leaves that are labelled by a time coordinate, also known as '3 + 1' decompositions, appear very commonly in applications of general relativity (and alternative theories of gravity), and especially in the cosmological context.Such foliations allow for an initial value formulation of general relativity [1,2].
Due to the non-uniqueness of time parameterization in general relativity, there is a broad freedom in choosing the foliation of a given spacetime.The non-uniqueness of the foliation is not a problem per se: it may be seen as an advantage that there is the freedom to consider a foliation where the physical phenomena of interest are more easily described.Once one applies operations that are tied to a particular slicing, such as the integration or averaging of a field over the leaves of the foliation, the choice of foliation is important.In this paper, it will be our goal to systematically examine and quantify the impact of the choice of foliation in spatial averaging with a focus on applications for cosmology, and we thus give a brief review of the cosmological averaging problem here.
Cosmology typically aims at describing overall, statistical, or average properties of the Universe as a whole as a function of cosmic time.Thus, the concept of 3+1 foliations is inherent in the questions posed in cosmology.The formulation of a large-scale effective cosmological theory can be approached either (i) by making an ansatz for the large-scale metric (that must be assumed to be a meaningful object) and for the large-scale matter content and seeking solutions within this ansatz -possibly also allowing for smaller-scale perturbations thereof; or (ii) by explicitly averaging over small-scale dynamics to derive the large-scale evolution laws of the locally defined spacetime.
The route in (i) is most often applied.It works well when the small and intermediate scales of gravitational phenomena can effectively be ignored in the large-scale evolution of the Universe, and it becomes simple when there is a notion of large-scale symmetries.When assuming homogeneity and isotropy over spatial sections of the Universe at its largest scales, one arrives at a Friedman-Lemaître-Robertson-Walker (FLRW) metric description by this route.The Lambda Cold Dark Matter (ΛCDM) paradigm of cosmology is founded on this approach, where structures are described as perturbations around a large-scale (and background) FLRW model.
The alternative route in (ii) is more involved, since it starts with the local spacetime description -which must necessarily be complicated by accounting for the hierarchy of scales and nonlinearity of gravitational phenomena in the Universe -as the basis for deriving the dynamics of the largest scales.This step also involves the construction of an appropriately defined coarse-graining or averaging operation.If the assumptions of decoupling of the physics at scales comparable to the size of the visible Universe from physics at smaller scales hold, and if large-scale homogeneity and isotropy apply, then one should, by such a procedure, arrive at the same FLRW metric description as for the approach in (i).However, if these conditions are not satisfied, then the results of the route in (ii) is expected to differ from the usual applications of the procedure in (i) -such a difference is called a backreaction effect of the dynamics of inhomogeneous structures.Despite of its empirical success, the ΛCDM paradigm is subject to model anomalies on a wide range of scales [3][4][5][6][7], and continues to face the fundamental challenges of the nature of the dark energy and dark matter sources.There is a debate regarding whether exploring the route in (ii) could help resolve the anomalies and interpretational challenges of the ΛCDM paradigm; see for instance [8,9].
The problem of how to formally approach the route in (ii) may be denoted as the averaging problem or the fitting problem, and was first discussed in detail in early works by [10][11][12].The most widely studied procedure for cosmological averaging is the Buchert averaging scheme for scalar quantities [13][14][15][16], but there have been numerous contributions on complementary/alternative procedures, e.g., [17][18][19][20][21][22].A variety of applications, mainly following the Buchert averaging scheme, and generalisations thereof, have directly addressed the impact of the choice of foliation in these schemes, either quantitatively or in qualitative discussions [23][24][25][26][27][28][29].Very briefly summarized, these papers point towards a variation of the cosmological averages performed in a given spacetime with the choice of foliation, which is indeed natural: the choice of foliation is defining the averaging domains and therefore also the final values of the computed averages.Such a variation may be worsened by additional dependences of the averaging scheme in the choice of slices, e.g. if the quantities to be averaged or the spatial boundaries of the domain depend on the foliation too.Without such extra dependences, it may be possible to select from physical constraints a broad class of foliations within which averaged observables remain approximately invariant, as argued qualitatively in [28].The foliation dependence in cosmological averaging is not in itself a problem, but if left uncontrolled, it can introduce ambiguity in the average cosmological theory.
Foliation dependence in cosmological averaging has not yet been quantified systematically in the literature.In this paper, we present covariant and broadly applicable methods for quantifying, eliminating and/or constraining the foliation dependence in cosmological averaging.Although the methods are investigated with the cosmological averaging problem in mind, they have broad applications to foliation studies in general relativity and differential geometry.
In section 2, we present the general scalar averaging formalism that we use to quantify foliation dependence throughout this paper.We then consider foliation (in)dependence in an exact way, in section 3, through calculus of variation, for integral and average functionals.We also investigate how foliations might be singled out uniquely from their extremal properties.In section 4, we consider bounds on foliation dependence for finite variations, which are relevant when considering physically motivated restricted classes of cosmic foliations.Finally, in section 5, we summarize and discuss our results and their potential application to the foliation dependence of large-scale backreaction terms and effects in cosmological averaging.

Covariant averaging over spatial foliations
We now define an averaging scheme relevant for averaging scalar functions over 3-dimensional (spatial or null) slices as embedded within the four-dimensional spacetime manifold M. We use the 3 + 1 averaging formalism presented in [22], building upon [20,21,30] -see also the related development in [31].This scheme is very general for 3 + 1 averaging of scalars1 , as it allows for the selection of an arbitrary averaging domain within an arbitrary slicing, and -as an extension to [20,21,30] -it further allows for any volume measure and weighting function for the integrand.
Following [22], we consider spacetime domains that are selected by a window function of the form: where H is the Heaviside step-function and δ D the Dirac delta function.Above, the foliation scalar A, with an everywhere nonzero gradient, defines the foliation, each level set {A = cst.}defining a hypersurface.The gradient of A is assumed to be either time-like or null, making the hypersurfaces space-like or light-like, accordingly.The constant A 0 parametrizes the hypersurface within the foliation.The boundary scalar B is set to have a space-like, outwards-pointing gradient, hence one may think of B as a radial coordinate in a suitable coordinate system.Together with the (fixed) constant B 0 , it defines the spatial boundaries of the finite domain of averaging within each hypersurface.We shall restrict ourselves to the case where B is independent of A and its derivatives, i.e., by selecting the region where B ≤ B 0 we consider a fixed spatially bounded (3 + 1)-dimensional tube in spacetime.
We illustrate this setup and these notations on Fig. 1.
The volume measure vector V , on the other hand, determines the integration measure on the selected 3-dimensional surfaces.For instance, for space-like hypersurfaces, taking V to be their unit normal vector n, i.e., reduce to the Riemannian volume element within the hypersurfaces, inside the integration domain, where g is the modulus of the determinant of the spacetime metric g (the components of which are noted g µν ), g ≡ | det(g µν )| = − det(g µν ).We shall, however, allow V to be any time-like vector field, and the integration measure will thus generally not coincide with the Riemannian volume measure on the leaves.Allowing for a non-Riemannian volume measure is convenient for certain types of fluidintrinsic averaging [15,28].Those formalisms are recovered by setting V to be the 4-velocity u of a Representation of the foliation and integration domain setup defined by the scalar functions A and B in the window function (2.1).The figure shows several hypersurfaces from the foliation defined by the level sets of A, the spacetime tube specified by {B ≤ B 0 } bounding the spatial extent of the region of interest, and the compact integration domain defined within each {A = cst.}slice from the intersection of the slice with the above tube.Note that we here illustrate the case of space-like slices, but a function A will light-like level sets could be chosen as well.
certain fluid flow, independently of the foliation set by A, and the volume measure then differs from the Riemannian one on the slices (if ∇A is time-like) by the Lorenz factor γ ≡ −u • n.The generality of V shall also be convenient in the present paper for making explicit which properties of the averaged expressions are related to the domain selection and which are related to the volume measure.We finally note that V is not necessarily requested to be unitary.This allows for weighted averageswhere the norm |V | of V can be considered as the local weight and the normalized V as defining the volume measure -, such as mass-weighted averages when |V | is a mass density [22,32].
For the examples of applications considered in this paper we mainly have in mind spatial foliations (∇A • ∇A < 0), which are relevant for formulating averaged evolution equations in time, i.e., viewing the averaging problem as an initial value problem.We shall however also be interested in considering null foliations (∇A • ∇A = 0) in some cases.In such situations, we will sometimes be interested in allowing B to rather have a time-like gradient.Then, {B = B 0 } would become a space-like hypersurface rather than a time-like tube boundary, and the regions selected by W A,A 0 ,B,B 0 ,V for a range of values of A 0 would then rather be interpreted as a set of light-cones truncated at that fixed hypersurface.
For the unit Heaviside step-function H, we use the right-continuous convention H(0) = 1 throughout.In the following we shall omit the subscripts on W A,A 0 ,B,B 0 ,V except for those relating to the foliation and refer to the window function (2.1) as W A A 0 , leaving the choice of boundaries and volume measure implicit.
We define the integral over a scalar S over the spacetime domain singled out by W A A 0 in the following way We then define the volume V A A 0 of the domain and the volume average ⟨S⟩ A A 0 of an arbitrary scalar S, respectively, as 3) The derivative of the integral (2.2) with respect to A 0 is given by (see [22]): and the analogous derivative of the average (2.3) is given by: When A is a time-function, meaning that ∇A • ∇A < 0, these derivatives may be thought of as describing the time-evolution of the integral/average.In the above expressions, we re-introduced the volume weighting factor V • ∇A to explicitly separate the window function and the scalar to be integrated (generally foliation dependent due to its V σ ∇ σ A denominator).For instance, given our H(0) = 1 convention, Eq. (2.4) may be rewritten as,

∂I(S)
We do require that B be parametrized in such a way that space-like or null sections of {B ≤ cst.} tubes are compact.In the case where the spatial sections of M are infinite, integrals or averages over the entire {A = cst.}slices may be obtained from a B 0 → ∞ limit, the existence of which would however set some convergence conditions that S must obey; we will not directly consider such a limit in this work.In the case of a spatially closed topology, on the other hand, the above formalism can directly encompass the special case of integrating over the entire compact slices.This is obtained by choosing a large enough B 0 ; in that case, the corresponding equations may as well be simplified by dropping the Heaviside factors H(B 0 − B).Boundary conditions arising from derivatives of this factor, i.e. δ D (B − B 0 ) terms, then disappear as B < B 0 everywhere within M.

Infinitesimal variation of the foliation
In this section we consider local extrema of the integrals and averages defined above, when they are viewed as functionals of the foliation.We thus compute stationarity conditions of the integrals/averages under variation of the foliation scalar A. Such stationarity conditions have at least two useful applications.Firstly, we are interested in identifying foliation independent quantities.Foliation-independence is equivalent to stationarity of the given functional for all foliations.Secondly, we may be interested in finding a foliation (or leaf of a foliation) that extremise a specific integral or average functional.Extremal leaves of such functionals can be thought of as generalisations of paths of shortest distance, and they may be thought of defining natural hypersurfaces in specific contexts.Thus, such extremals may provide an interesting way of defining leaves/foliations uniquely.Moreover, the corresponding extremum value of the functional would be a natural foliation independent measure of this functional, generalising the shortest distance (minimised over all possible paths) as a preferred measure of distance between two points.
We leave the foliation scalar unconstrained in the variation.Thus the leaves associated with the solutions to the resulting stationarity conditions may in principle be space-like, time-like, null, or of a varying nature (depending on the point) amid these three possibilities; and the nature of the solution will have to be checked in each case.Analogous constrained stationarity conditions can be derived by imposing local or global constraints of physical interest, such as constraining the foliation scalar to have an everywhere null gradient (e.g. when discussing light cones) or to be a proper-time function for a given 4-velocity field.We avoid such considerations in this paper for the sake of simplicity.
We shall consider stationarity conditions either for a single leaf or for an entire foliation, the latter being more restrictive.We consider cases where the integrated scalar S of interest is possibly dependent on ∇A, but not A or higher-order derivatives, and where the volume measure vector V is possibly dependent on A and ∇A2 .

Freedom of parameterisation of the foliation
Before considering the variation of the average/integral functionals in generality, we shall consider the trivial subset of variations that define a map of the foliation onto itself.We shall require that the functional is invariant under such mappings.Consider a given foliation F = {Σ A=A 0 }, where A is a spacetime scalar representing the foliation, A 0 is a parameter defined over some range A 0,1 ≤ A 0 ≤ A 0,2 selecting a leaf of the foliation, and where Σ A=A 0 represents a leaf of the foliation.The transformation where f is a strictly monotonic function, defines a map from the foliation onto itself F → F. We thus denote the choice of f the freedom of parameterisation of the foliation.We can define an averaging operation independent of this parameterisation by requiring that a change of parametrisation leaves the integral (2.2) invariant when S is itself independent of the parametrisation of the foliation.Thus, we require that transformations A → f (A), A 0 → f (A 0 ) lead to mappings of the integral onto itself I(S) A A 0 → I(S) , where S is an arbitrary scalar independent of the parametrisation of the foliation F. Demanding that I(S) A A 0 remains invariant under the gauge transformations A → f (A), A 0 → f (A 0 ) is equivalent to requiring stationarity of I(S) A A 0 for all representations A of the foliation F under variations A → Ã = A+δf (A) , A 0 → Ã0 = A 0 +δf (A 0 ), where δf (A) (and its derivatives) is an infinitesimal scalar function of A.
For a generic integral operation (2.2) the integral evaluated at the representation Ã of the foliation can be expressed in terms of the integral evaluated in the representation A as follows where the second line follows from We have made the functional dependence of V explicit by writing V µ [ Ã, ∇ Ã].Performing the first order functional expansion of V µ ( Ã) around the representation A and plugging it into (3.2) we have where δI(S) A A 0 F denotes the first order variation of δI(S) A A 0 under variations of A, A 0 that maps the foliation F onto itself.We require (3.3) to vanish for any choice of S. Since δf (A 0 ) and its derivative can be chosen arbitrarily and independently on a given leaf A = A 0 , the following two conditions, must be independently satisfied.Loosely speaking, (3.4) states that the volume measure (as determined by V ) can depend only on the direction defined by ∇A, but not on its norm or the values of A. Besides any foliation-independent V , these conditions allow in particular for taking V as the unit normal to the hypersurfaces as in [20,21], in the case where those surfaces are space-like.

Variation of integral quantities with respect to the foliation
We shall now derive the stationarity conditions for the integral functional (2.2) under variation of the hypersurface scalar3 A. Physical integral functionals of interest may for instance be volume or mass functionals.
We write the first order variation of the integral I(S) A A 0 (2.2) as a function of the variation A → A + δA as with We note that local stationarity requirements cannot yet be directly extracted from the variation as written in (3.6), because δ(∇ µ A) = ∇ µ (δA) may not be considered fully independent of δA.It is only the component of ∇ µ (δA) along a direction pointing away from the surface, that may be considered independent from δA as evaluated on the surface.Hence, we make the following rewriting, first using partial integration of (3.6) to obtain and then introducing an arbitrary, fixed vector field Z chosen such that Z •∇A ̸ = 0 (this is for instance guaranteed for a time-like Z, if A has an everywhere time-like or null gradient) to rewrite the first term in (3.7).This results in the following useful form of the variation: where the requirement of gauge-invariance of the averaging operation (3.4) has been used, and with In Eq. (3.8) above, the variation δA on a given surface can now be considered independent from its derivative as the latter is taken along the direction Z, away from the surface.Accordingly, the first line of (3.8) gives rise to a first local constraint equation, and the last two lines give rise to a second one.Altogether, we arrive at the two independent constraints for stationarity of I(S) A A 0 around the surface A = A 0 : where we have made use of the assumption V ν ∇ ν A ̸ = 0.The first condition above is nothing but the requirement for S to be invariant under a change of the foliation F = {Σ A=A 0 }.Note that the conditions for stationarity in Eq. (3.9) are independent on the arbitrary choice of Z despite the apparent dependence on Z in the second condition.This can be seen by expanding Z in a component proportional to ∇A and a component orthogonal to ∇A, and noting that all contributions orthogonal to ∇A vanish due to the first condition of (3.9) and to its spatial derivatives along the {A = A 0 } hypersurface.
Requiring that I(S) A A 0 is stationary for all surfaces A 0 results in the following conditions: Note that the stationarity requirements (3.9) and (3.10) for the integral I(S) A A 0 are local conservation equations.Due to the first condition, when S is dependent on the representation A of the foliation such that ∇ µ A ∂S/∂(∇ µ A) ̸ = 0, no stationary foliations exist for I(S) A A 0 , when allowing for all possible variations of the foliation5 .Due to the derivative of the spatial boundary selection function H(B 0 −B), the second condition in either Eq.(3.9) or Eq.(3.10) may itself be split into a condition inside the B ≤ B 0 domain and a condition on the B = B 0 boundary.For instance for Eq.(3.10), the second condition is expanded as and is thus equivalent to: It can be shown from (2.4) and the conditions (3.9) that stationarity of I(S) A A 0 implies the vanishing of its derivative with respect to A 0 : These results are expected since (infinitesimal) constant translations of A at fixed A 0 , which are part of the class of variations δA, can also be seen as a translation in time within the original foliation.In fact, the results (3.13) and (3.14) are general and apply to integral functionals with arbitrary functional dependence on ∇A, ∇∇A, .., ∇ (n) A.
We now consider the special case where S and V are independent of the foliation.In particular, we let S and V be independent on the direction vector ∇A of the foliation, and it follows that in the case of extremisation for a single leaf; or, for extremisation for the entire foliation.In these cases, the first condition of (3.9) and of (3.10) respectively are automatically satisfied.Moreover, the condition for stationarity for a single hypersurface simplifies to δI(S) and the stationarity condition for all surfaces reduces to The requirements (3.15) or (3.16) are indeed natural in many cases 6 .In the case where V is independent of the foliation, the stationarity requirement (3.17) for a single surface is only dependent on the foliation through the surface of evaluation, and the requirement (3.18) for stationarity for an entire foliation does not depend on the foliation at all.In this case, if (3.18) is satisfied for a particular foliation, then it is satisfied for any possible foliation.
As in the general case of Eqs.(3.9)-(3.10), the above stationarity conditions may be split into a bulk and a boundary condition, to be simultaneously satisfied either on the single {A = A 0 } slice or on all slices.
We remark that the stationarity conditions for integral functionals are restrictive.In particular, the stationarity requirement (3.18) for an entire foliation is extremely restrictive.The foliation is stationary only when SV is a conserved current.The stationarity requirement (3.17) for a single leaf is more flexible: obtaining stationarity in a given foliation amounts to being able to collect points for which ∇(SV ) = 0 is satisfied to construct a space-like (or null-like) leaf.We emphasise that the stationarity conditions (3.17) and (3.18) are derived under assumptions.For instance, they are derived under the assumption of no functional dependence of I(S) A A 0 on second or higher order derivatives of A (thus, for integrals over foliation-adapted curvature degrees of freedom the results derived in the present section do in general not apply).We also emphasise, that extrema which are not stationary points can exist in the form of infimums or in the form of local extrema introduced by a "boundary" in the solution space.

Consider a conserved local rest mass current
where ϱ is a rest mass density of a fluid with 4-velocity u.We might seek to define a total rest mass of spatial domains as volume integrals over ϱ.
Let us consider the case where we take V to coincide with the fluid 4-velocity: V = u.This is a natural choice when averaging local quantities intrinsic to the fluid [14,28], and it allows for the integrated ϱ, I(ϱ) A A 0 , to properly define a total rest mass [22].We will assume that we can define boundaries B intrinsic to the fluid flow through parallel transport, u • ∇B = 0, such that the spatial boundary is comoving with the fluid flow.We may, alternatively, consider the case of a boundaryfree domain, i.e., for a spatially closed manifold and an integration domain coinciding with the whole hypersurfaces.Either way, this ensures the preservation of the fluid content of the averaging domain over time, and will accordingly allow for the preservation of its total rest mass.
Accordingly, we define a total rest-mass associated with the (fluid-comoving) averaging domain in a given foliation A as with The mass (3.21) can be shown to be conserved over time (∂M A A 0 /∂A 0 = 0) for any choice of foliation A [22].Furthermore, we have from (3.8) that and we recover the condition (3.18), ∇ µ (ϱu µ H(B 0 − B)) = 0, for the stationarity of M A A 0 with respect to the foliation.Using the conservation (3.20) of the rest mass current and the comoving boundary (or no-boundary) assumption, always vanishes, implying that the mass is stationary for all foliations and thus foliation independent.
The rest mass definition (3.21) and its stationarity still hold if ∇A is light-like7 .Thus, this mass is not only invariant with respect to a choice of spatial hypersurface; it is also invariant under a change from spatial to light-like hypersurfaces, providing a potentially interesting correspondence between light-cone and spatial hypersurface averaging.

Suppose we want to consider extremal foliations for the volume V A
A 0 ≡ I(1) A A 0 of a domain lying within the hypersurfaces.We consider variations A → A + δA that vanish (δA = 0) on the boundary B = B 0 (or, as above, that the domain is boundary-free, i.e.B < B 0 everywhere on M).This implies slightly weaker stationarity requirements, removing the conditions imposed at the domain boundary since only the interior region has to be constrained.
The problem of finding the extremal volume enclosed by a fixed spatial boundary, can be thought of as a higher-dimensional generalisation of finding the shortest path between two fixed spacetime events.This is a well-studied problem in the literature, at least for the Riemannian volume measure (see e.g.[33] and references therein for the case of Riemannian geometry and [34] and references therein for Lorentzian manifolds).It is nevertheless worth recalling here as an illustrative example for the present discussion, as the volume functional is of great importance in cosmology.
We first consider the case where the hypersurfaces are space-like and V is their unit normal.In this case, we recover the Riemannian volume measure on the surfaces.We write where n is the future-pointing unit normal to the hypersurfaces, and N −1 is the associated lapse function.In this case the condition (3.16) is satisfied, and the stationarity requirement reduces to (3.18).This simplifies further due to the fixed value of A at the boundary (or the lack of a boundary), removing the condition at It follows that a foliation extremises I(1) A A 0 if and only if the extrinsic curvature scalar of each hypersurface vanishes inside the domain, ∇ • V = ∇ • n = 0. We thus recover the well-known condition for the stationarity of this volume functional.
Suppose that we can find a foliation for which this condition is satisfied.We want to determine whether such an extremal foliation is a maximum, minimum, or a saddle-point.For this purpose we use the identities to compute the second variation of the volume where the second last equality follows from partial integration of the term involving N δA 2 (∂ 2 δ D (A 0 − A)/∂A 2 ), the condition for stationarity ∇ • V = 0, and (3.25).The last equality follows from (3.25).Since h µν is positive semi-definite, we have δ 2 V A A 0 ≤ 0, with equality when the perturbation δA depends solely on A (with ∇δA ∝ ∇A).Such perturbations simply map the foliation onto itself -and in fact, they have to vanish when the domain does have a boundary, due to the boundary condition δA| B=B 0 = 0 that we set in this case for the current example.Thus δ 2 V A A 0 < 0 for any infinitesimal actual change of foliation, and hence the extremal foliation maximises the volume.This is an expected consequence of the Lorentzian signature of the metric, while a Riemannian signature would induce a minimisation of the volume.
We next consider the case where V is independent of the foliation.For variations of A that are zero on the boundary B = B 0 and in the boundary-free case, all foliations are extrema of the volume if and only if ∇ • V = 0, and it follows that the volume is foliation independent if this condition is satisfied.Furthermore, we trivially have δ 2 V A A 0 = 0 for all foliations in this class.The volume has no local extremum, on the other hand, if ∇ • V ̸ = 0.This is the case for a fluid proper volume measure, i.e.V = u where u is the 4-velocity of a fluid source [22,28], if any expansion or contraction of the fluid occurs within the integration domain.
Choosing V as a conserved current defines a foliation-independent "volume".We can consider again the example of a conserved rest mass current M = ϱu (3.20) from the above subsection, and set V = ϱu.In this example, the conserved "volume" is simply the total rest mass as defined in (3.21), and averages with such a window function are mass-weighted [22,32].The foliation-independence is still, in principle, restricted to the case δA| B=B 0 = 0 or that of a boundary-free domain, but is recovered for any deformation and more generic domains with the additional requirement of a domain boundary comoving with the rest mass current, u • ∇B = 0.
A few remarks are in order on extremal leaves.
In this section we have commented on volume extremising foliations.However, we remark that a single leaf that extremises volume -or another integral functional of physical interest -could be used as a preferred surface for the initial value problem in cosmology.A preferred initial surface may in turn be extrapolated to form a foliation, by for instance propagating the initial surface along a physically motivated 4-velocity field.Since stationarity requirements for a single leaf are easier to satisfy than for a full foliation, single leaf extrema may be explored in cases where it is not possible to identify a full extremal foliation.

Variation of averaged quantities with respect to the foliation
We will now derive stationaity conditions for the average functional (2.3) under variations of the hypersurface scalar A, analogous to the above results for the integral functional.Examples of physical average functionals of interest in cosmology are average density, expansion rate, and spatial curvature degrees.
We write the first order variation of the average (2.3) under the variation A → A + δA as The variation of the average can be expressed through the variation of integral quantities in the following way and we can plug in (3.8) to obtain stationarity conditions for averaged quantities.The conditions for demanding stationarity for the entire foliation are The condition for stationarity of an average functional for a single leaf can be similarly derived from (3.8), providing an analogous criterion to (3.9) for integral functionals.We do not include this condition for simplicity given that such a situation will not be of interest for the investigations discussed below.We also note that, as for stationary integrals above, the second condition in Eq. (3.29) -or its singleslice equivalent -may be further split into a bulk (B ≤ B 0 ) and a boundary (B = B 0 ) conditions stemming from the domain selection factor H(B − B 0 ) and its derivative.The global constraint equation must be satisfied in order to obtain stationarity for a single slice of a foliation selected by A = A 0 and for stationarity for the entire foliation respectively.The necessary condition (3.30) is analogous to the integral condition (3.14).As for (3.14), the condition (3.30) can be shown to hold for averages arising from a general window function with arbitrary functional dependence on ∇A, ∇∇A, . . ., ∇ (n) A.
We now consider the special case where S and V are independent of the foliation.In this case, the constraints (3.16) are satisfied.Stationarity conditions for the entire foliation, Eq. (3.29), accordingly reduce to 30), we see that the above condition is equivalent to ⟨S⟩ A A 0 being constant (independent of A 0 ) together with S − ⟨S⟩ A A 0 V σ being a conserved current comoving with the boundaries of the domain, such that The latter boundary condition can be neglected if we consider variations that are fixed on the boundary, or a global averaging domain on a spatially closed manifold.
We remark that the conditions for stationarity of average functionals are very restrictive.The existence of extremal foliations are conditioned on the existence of a locally-conserved current.Considering the case V = u, where u is a fluid four velocity field with an associated rest mass density ϱ, we note that a natural conserved current is ϱ u: ∇ µ (ϱu µ ) = 0.However, since ⟨ϱ⟩ A A 0 ≥ 0 with equality only when ϱ = 0 everywhere, it follows that setting S − ⟨S⟩ A A 0 ∝ ϱ cannot generate stationary solutions to (3.31), except for the trivial case of taking S as a constant of spacetime.We note that ϱ might itself be such a constant, but in addition to a homogeneous rest-mass distribution, this would imply a non-expanding fluid flow, ∇ µ u µ = 0, unless ϱ = 0, due to the conservation equation.
We note that solutions found in this section are valid only for the specified functional dependence of the average functional on the foliation.As discussed for integral functionals in section 3.2, there might exist extremals that are not stationary points, occuring as infimums or as local extrema on the boundary of a set of allowed foliations.

Example: Entropy
The study of entropy is a rich topic in gravitational physics, and in cosmology it has found various applications, for instance in the characterisation of initial conditions and inflationary scenarios [35], cyclic universe models [36], and structure formation [37].Here, we focus on the following entropy measure, inspired by the Kullback-Leibler relative information entropy [38,39]: where S must be a field in which gravity induces an increased clustering or inhomogeneity.Physically relevant substitutions for S include rest mass densities ϱ as in [39], expansion scalars θ (in case of positive expansion everywhere), and proper time measures τ .The variation of S with the foliation reads We restrict our investigation to the case where V is independent of the foliation.In this case, the condition (3.16) applies, and using (3.8), the stationarity condition becomes We have here explicitly split the condition into its boundary and bulk components (respectively from δ D (B − B 0 ) and H(B 0 − B) terms).Accordingly, the second requirement above is to be satisfied everywhere inside the B ≤ B 0 tube; the first requirement, as a boundary condition, is dropped in the case of a boundary-free averaging domain.Assuming that S is hypersurface-forming, the foliation into constant-S hypersurfaces, A = S, satisfies (3.34) since S is then by construction homogeneous over each slice, and S = ⟨S⟩ A A 0 .It follows that S is stationary with value S = 0 for this foliation.This point of stationarity is guaranteed to be a unique global minimum when S is a strictly positive function (as is the case for a rest mass density) since where the inequality follows from Jensen's inequality for the function x → x ln(x), with equality only if S is constant over the A = A 0 hypersurface.For a hypersurface-forming S, such an extremal-foliation exists irrespective of how V and B are defined.However, to satisfy non-singularity of the integration measure (V • ∇A ̸ = 0) for this foliation, we would have to demand V • ∇S ̸ = 0. Note that in case ∇S • ∇S < 0 is not fulfilled everywhere, the constant-S solution does not define a spatial foliation.This situation may occur for the choices of S suggested above, S ∈ {ϱ, θ, τ }, for which the hypersurface-forming property of S must always be checked.Moreover, in general, the constant-S foliation needs not be the unique solution to (3.34), i.e., there may be other local extrema.
We now consider the special case where SV is a conserved current.In this case, we have ∇ µ (SV µ ) = 0.This includes the physical example where V is a 4-velocity field and S is the associated rest mass density.In this case, the stationarity conditions (3.34) reduce to For ∇ • V = 0 the second condition of (3.36) is automatically satisfied for any foliation.Hence, in this case, the entropy is foliation independent and has a constant value (which is zero if and only if S is a constant of spacetime), up to the boundary condition, i.e. the first condition of (3.36).The latter can be accounted for by either imposing V • ∇B = 0, by considering the particular case of a global boundary-free integration domain, or by simply keeping the foliation fixed at the domain's boundaries, δA| B=B 0 = 0. Note that this case of ∇ • V = 0 also implies V • ∇S = 0, and thus a foliation defined by A = S (where S would realise its global minimum at S = 0) would result in a singular volume element.
For ∇ • V ̸ = 0 everywhere, (3.36) is equivalent to S = ⟨S⟩ A A 0 , i.e., the global minimum of S corresponding to the constant-S foliation is the only local extremum.It follows for instance that, for an everywhere expanding fluid with 4-velocity V = u and rest mass density S = ϱ, the constant-ϱ foliation is the unique minimiser for the entropy.

Example: Minimally differing frames
Suppose that we have a physical time-like vector field u in our cosmological theory in the frame of which averaged quantities would be desirable.This would for instance mean averaging in the rest frame of a fluid source if u represents its 4-velocity.This 4-velocity field can have vorticity, which will prevent defining hypersurfaces that are orthogonal to its flow lines.In this case we may ask whether there is a unique space-like foliation (defined by a scalar A with time-like gradient), or a family of foliations, such that their normal vector field n, given by is maximally close to u by some measure.In cases where such a foliation could be defined, this would provide a natural frame for definining averages as close to the frame of u as possible.The tilt between the two normalised time-like vector fields is a natural local scalar measure of their closeness, where γ = 1 if and only if n = u.We define the measure of "statistical closeness" of the vector fields n and u over the domain defined by {A = A 0 , B ≤ B 0 }, as with being the window functions with volume measure defined with respect to the normal vector n (3.37) and the vector field u, respectively.Plugging in V = n and S = γ in (3.29), we have that the first condition of (3.29) is automatically satisfied while the second condition of (3.29), split into its bulk and boundary components, becomes The second condition in (3.41), i.e., the boundary condition, needs to be satisfied only if the domain has a boundary and we include variations of the foliation at that boundary as well as on the domain's interior.A necessary condition for (3.41) is ∂ A 0 ⟨γ⟩ A A 0 = 0.In solving (3.41), we can thus consider ⟨γ⟩ A A 0 as a constant parameter.
Investigations of the general solution to Eq. (3.41) is beyond the scope of this paper.However, we consider here the special case where u is divergence-free, ∇ • u = 0.The first condition of (3.41) then reduces to ∇ • n = 0. Finding a solution to the first condition of Eq. (3.41) in this case amounts to examining the existence of zero extrinsic scalar curvature foliations.The problem of extremising the averaged tilt (3.39) thus becomes equivalent to extremising the volume as in section 3.2.2.This is because M d 4 x √ g W A A 0 ,u becomes foliation independent, so that we are extremising the inverse of the Riemannian volume (as defined by M d 4 x √ g W A A 0 ,n ).From the results in section 3.2.2we know that stationary points for the Riemannian volume are maximal, meaning that the stationary points for the averaged tilt (3.39) are minimal in this case.

Finite foliation changes and quantitative bounds on foliation dependence
The results of section 3 show that, while foliation independent statements can be made for special cases, most integral or average functionals are foliation dependent.Nevertheless, it may be possible in many cases of interest to quantify the level of foliation dependence.The aim of this section is to determine quantitative bounds on the foliation dependence in scenarios relevant for cosmological models.

Correspondence between hypersurfaces of different foliations
In the following, we will consider two different foliations F and F ′ corresponding to the respective level sets of two scalars with past-pointing time-like or null gradients8 , A and A ′ , where the transformation A → A ′ needs not be infinitesimal.
In order to make comparisons of leaves Σ 0 ≡ Σ A=A 0 and Σ ′ 0 ≡ Σ A ′ =A ′ 0 of the two foliations in a meaningful way, we must ensure that the two slices correspond to the "same time" in some sense.We shall specifically require that Σ 0 and Σ ′ 0 intersect at at least one event within the bounded region determined by the tube T B 0 ≡ {x ∈ M / B(x) ≤ B 0 }.This ensures a notion of synchronisation at at least one point within the domain of interest.It prevents in particular artificial differences due to the comparison of two different parametrisations of the same foliation.
We can always choose a parametrization such that This can be achieved by using the freedom of reparametrization of F ′ as per transformations of A ′ of the class (3.1).The requirement of intersection of the pairs of corresponding slices from the foliations will suffice in what follows, even though it does not always uniquely specify the parametrization9 A ′ for F ′ from the parametrization choice A for F. It already ensures, in particular, that A ′ must be chosen as equal to A (i.e., the transformation A → A ′ = f (A) reduces to the identity) in the case F = F ′ mentioned above.

Simplifying assumptions and notations
In this section, we shall only consider cases where S and V are invariant under deformations A → A ′ of the foliation.As has been argued in the above sections, defining V as a physical vector field independent of the foliation is natural for many purposes -a natural choice of V could for instance be the 4-velocity field of a physical matter source.In most cases we will also be interested in averaging physical scalars that are independent of the foliation10 , as for instance the rest mass density or expansion rate of a physical matter source.
The boundaries of the domain and their propagation between slices, as determined by the scalar B, are already considered to be set independently of the foliation as part of our averaging scheme.Here, we shall make the simplifying assumption that the domain propagation follows the flow of the volume-measure vector: V • ∇B = 0 -unless a boundary-free domain is being considered so that boundary terms can already be discarded, with H(B 0 − B) = 1 and δ D (B − B 0 ) = 0 everywhere, and T B 0 = M.We set V to be unitary, V • V = −1; a non-normalized vector field (corresponding to weighted volume averages) could formally be absorbed into the scalar S to be averaged, The above assumption and normalization convention are again compatible with the choice of V as a source fluid's 4-velocity; in this case with a fluid-comoving domain propagation.This is indeed one of the main applications we have in mind for this section (see, e.g., [14,28]).

The difference of integral functionals between leaves of two foliations
We now consider the difference between integrals of a given scalar S over leaves of two foliations F and F ′ : ∆I(S) ≡ I(S) A ′ A 0 − I(S) A A 0 , where S is any foliation-independent scalar.Note that we are using the short hand notation for the difference, ∆I(S), where the dependence on the foliations and the leaves selected by A 0 is implicit.We shall rewrite this difference in a way that will be convenient for defining upper bounds for its norm.

∆I(S) in terms of covariant 4-integration
We consider a given A 0 and the corresponding pair of intersecting slices Σ 0 = Σ A=A 0 , Σ ′ 0 = Σ A ′ =A 0 from the two foliations F, F ′ , obeying the above assumptions.While the roles played by F and F ′ are formally symmetric, we will consider F as the reference foliation in which the integrals I(ξ) A A 0 or averages ⟨ξ⟩ A A 0 of various scalars ξ are known.F ′ , on the other hand, will be considered as an arbitrary other space-like or null foliation that may be subjected to certain conditions, such as having a small tilt everywhere with respect to F in the case where both foliations are space-like.We will keep F ′ fully general in the present subsection.We assume that M is a path-connected manifold, hence so are Σ 0 and Σ ′ 0 according to the global hyperbolicity assumption.In the following, we will make use of the flowlines of V as a diffeomorphism between the domains of Σ 0 and Σ ′ 0 that are within the tube T B 0 .While this mapping is covariantly defined, it will be convenient to introduce an associated set of spatial coordinates.We do so by arbitrarily choosing a coordinate basis (X i , i = 1, 2, 3) on one of the slices Σ A=A 1 of F (or some open subset of Σ A=A 1 containing Σ A=A 1 ∩T B 0 ).The three coordinates X i , assumed to span R 3 without loss of generality, can then be extended into an incomplete (spatial) set of coordinates in spacetime by requiring them to be constant along the flow lines of V : V • ∇(X i ) = 0, ∀i ∈ {1, 2, 3}.The spacetime tube T B 0 = {x ∈ M / B(x) ≤ B 0 } with V • ∇B = 0, then corresponds to a given compact domain in the space of the spatial coordinates (X i ); in other words H(B 0 − B) is a function of (X i ).(This also holds in the case of a spatially closed manifold with a global, boundary-free averaging domain.)One may then introduce any time coordinate to complete (X i ) into a spacetime coordinate set -again for convenience in the below calculations.We shall complete it into a synchronous-comoving coordinate set adapted to the 4−vector field V by introducing a proper time τ of V , i.e., a function satisfying V • ∇τ = 1 (see section 2.4.2 of [14]).We can use the residual freedom in the definition of τ to demand that τ = 0 on Σ 0 -even if Σ 0 is where we again made use of two relations holding in these coordinates: The first and last sides of Eq. (4.5), and hence their equality, are nevertheless independent of the choice of the time coordinate.The evolution rate of the volume measure √ b, which is additionally invariant where the operator (d/dτ )| X i corresponds to a derivative along V with respect to τ .Using the coordinates (τ, X i ), we can now explicitly integrate Eq. (4.6) into for any value τ 1 of τ .The main integrand in expression (4.4) for ∆I(S) is then rewritten as: with With this rewriting, the expression for ∆I(S) finally reduces to where ψ(X i ) is extended into a spacetime scalar ψ by defining the latter as equalling ψ(X i ) on some Cauchy hypersurface (say Σ 0 ) and satisfying V • ∇ ψ = 0; i.e., ψ(τ,X i ) = ψ(X i ) ∀τ .With a slight abuse of notation, one may simply write ∆I(S) = I ψ(X i ) A A 0 .We give in Appendix A.2 an alternative form of ψ(X i ) re-expressed in terms of the local current ∇ µ (SV µ ).This form allows for a more direct connection with our results on exact foliation-(in)dependence of Sec.3.2 and for alternative bounds on finite variations of spatial integrals and averages to the ones presented below.

Bounds for foliations with small relative tilts
Cosmology is typically studied under the assumption of small (nonrelativistic) relative velocities between relevant observers in spacetime.In this section, we shall therefore consider bounds on integrals relevant for a class of space-like foliations which are close to the reference field V (and to each other) in terms of their relative Lorentz factors.

The small tilts assumption
In this subsection, we restrict our attention to space-like foliations.We denote as n the future-pointing unit time-like normal vector field to the foliation F, i.e. satisfying where N is the lapse function associated with A. We similarly define n ′ as the unit time-like normal associated with F ′ and N ′ as the lapse associated with A ′ , with n ′ = −N ′ ∇A ′ .We can then define the local Lorentz factors between V , n and n ′ , as follows: and introduce the local decomposition of n ′ with respect to V : where v ′ automatically satisfies The key assumption that we shall use in this subsection is that both n and n ′ are close to V , that is, that the tilt velocities 1 − γ −2 V ,n and 1 − γ −2 V ,n ′ ( = v ′µ v ′ µ ) are small, i.e., globally bounded by a small parameter v 1 ≪ 1.This implies in particular that the relative tilt between slices of the two foliations also remain small and globally bounded: ) −1/2 ≃ 1, and we can introduce the global small parameter as the main characteristic spatial velocity to be used in the below bounds.

Bounding the distance between two tilted slices
From the expression of ψ above, Eq. (4.9), or its rewriting in Eq.(A.4), it is clear that in addition to assuming upper limits on the local S-and V -based variables |S| (on Σ 0 ), |V µ ∂ µ S|, |∇ µ (SV µ )| and/or |∇ µ V µ |, one also needs to be able to bound the proper-time parametrization function |τ 0 (X i )| of Σ ′ 0 for all points.This quantity provides a measure of the distance between the two slices Σ 0 and Σ ′ 0 , and in this section we shall provide a bound of |τ 0 (X i )| in terms of the small velocity parameter v 0 .
Figure 2. The spatial slices Σ 0 and Σ ′ 0 , the spatial curve C ⊂ Σ 0 , and the main points and vector fields used for expressing and bounding the time-like distance τ 0 between the two slices.P is an arbitrary point of Σ 0 , while P 0 , at coordinates (τ = 0, X i 0 ), belongs to the intersection Σ 0 ∩ Σ ′ 0 of both slices under consideration: τ 0 (X i 0 ) = 0. Note that for this schematic representation which is not specifically concerned with causality, we use a Riemannian picture of orthogonality for easier visualisation.
Consider any given point P -of coordinates (τ = 0, X i ) for a certain (X i ) -within the integration domain on Σ 0 .One can draw a geodesic spatial curve C within Σ 0 joining P to a reference point P 0 -say of coordinates (τ = 0, X i 0 ) -within the integration domain on Σ 0 where τ 0 (X i 0 ) = 0.That is, P 0 is taken as an intersection point of the two slices: P 0 ∈ Σ 0 ∩ Σ ′ 0 ∩ T B 0 which we assumed to be non-empty.(See Fig. 2 for an illustration of this geometric setup.)C can be parametrized by its unit space-like n-orthogonal tangent vector K, and the associated affine parameter λ.Setting λ = 0 at P 0 , λ then runs from 0 at P 0 to L at P , where L is the total proper length of the curve within Σ 0 .The coordinates of the point at parameter λ along C can then be parametrized as x µ (λ) = (τ = 0, X i (λ)).These definitions allow us to perform spatial integrations along C , using P 0 as a reference point where τ 0 = 0, writing for instance In the following we will use the short-hand notation τ 0 (λ) ≡ τ 0 (X i (λ)).
We may extend the tangent vector K of C into a vector field L on the congruence generated by V from C -i.e., along all worldlines of V that intersect Σ 0 at a point on C -by Lie dragging K along V : L V L ≡ ∇ V L − ∇ LV = 0, i.e., (d/dτ ) Lν | X i = Lµ ∂ µ V ν , with L = K at τ = 0, along any given C -intersecting worldline of V .Along these same worldlines, we can then introduce a unit space-like vector L built from L: The V -orthogonal projection b µ ν Lν is indeed nonvanishing, given that L is not parallel to V at τ = 0 (as K • K = 1 while V is time-like) and that this property is preserved by the Lie dragging.Note that, on points of C , K coincides with L (by construction), but a priori not with the projected vector L, since K needs not be orthogonal to the local V .
The derivative of τ 0 along C involved in Eq. (4.14) is then obtained as It can also be re-expressed as follows: where ) is V • K evaluated at the current point on C and where for any τ 1 .The derivation of these results, Eqs.(4.16)-(4.19), is given in Appendix A.3.In the above expressions for F and G, we have used the decomposition of the covariant derivative of V with respect to V [40]: This defines the expansion tensor of V , with components and its acceleration vector a, with components a µ ≡ V ν ∇ ν V µ .We can right away use the above global smallness assumption on the relative tilts between n, n ′ and V to derive bounds on the terms |V • K| and , where h µν ≡ g µν + n µ n ν are the components of the spatial projector (which is also the induced metric) on the leaves of F. Using this, the second above term obeys the following inequality: With the small tilt velocities v 1 , v 0 satisfying everywhere both of the above terms from Eq. (4.17) are everywhere smaller than v 0 .
We shall now assume the existence of a global bound on the norm of the expansion tensor.That is, we assume that there exists a bound on (Θ µ ν Θ ν µ ) 1/2 that applies throughout the part of M under consideration (i.e., a certain portion of T B 0 ).This is equivalent to assuming a global bound H on all eigenvalues of the (diagonalizable) matrix Θ µ ν over the local tangent spaces orthogonal to V , so that for any point P in the spacetime region of interest, for each eigenvalue θ k (P ) (k = 1, 2, 3) of Θ µ ν at P over the tangent space orthogonal to V at P , |θ k (P )| ≤ H.Note that these eigenvalues are real, and covariantly defined.This implies, in particular, that the norm of the volume expansion rate µ is also globally bounded: |∇ µ V µ | ≤ 3 H.We also set the acceleration of V to zero (see the discussion below for the case of a nonzero acceleration), so that G(τ 0 (λ), λ) = 0.
(for any τ ), and with G(τ 0 (λ), λ) = 0, Eq. (4.17) gives remembering that F > 0 by definition.One moreover has where we have used again the orthogonality of V and its expansion tensor, implying Injecting the above inequality into Eq.(4.18) implies F (τ 0 (λ), λ) ≤ exp H τ 0 (λ) for any λ ∈ [0, L], given that τ 0 (λ) ≥ 0. Eq. (4.22) then becomes, This may then be integrated between λ 0 , where τ 0 (λ 0 ) = 0, and any λ ≥ λ 0 , to give If the above right-hand side is nonnegative, which is ensured for all λ ∈ [λ 0 , L] provided v 0 is small enough such that v 0 H(L − λ 0 ) < ln(2), we then obtain the following bound on τ 0 (X i ) = τ 0 (λ = L) as a function of the tilt velocity v 0 : As the integration domain Σ 0 ∩T B 0 is compact, v 0 H(L−λ 0 ) admits a maximum; we can thus introduce a small parameter η (with still η < ln(2) to satisfy the above assumption) such that v 0 H(L − λ 0 ) ≤ η 13 This requires the integration domain Σ0 ∩ TB 0 to be path-connected, and geodesically convex or at least star-shaped with respect to P0.If this is not satisfied but Σ0 as a whole does satisfy these conditions, spatial paths may be drawn within an extended flow-tube TB 0 of V , encompassing TB 0 , such that Σ0 ∩ TB 0 satisfies these properties while keeping its diameter as small as possible (e.g., Σ0 ∩ TB 0 could be taken as the smallest sphere containing Σ0 ∩ TB 0 ).The assumed global bounds on tilts, Θµν , or the integration domain's diameter, then simply are understood to hold on a portion of TB 0 rather than only the corresponding one of TB 0 , while the domain of spatial integration is unchanged (Σ0 ∩ TB 0 or Σ ′ 0 ∩ TB 0 ).If Σ0 is itself not geodesically convex e.g.due to punctures, when referring to spatial curves within Σ0 the meaning of 'geodesic' has to be extended to a non-necessarily unique path taken as short as possible within Σ0 (for instance, we do not need K to satisfy the spatial geodesic equation within Σ0).
for any P .We can then use the convexity of the function x → − ln(2e . This then gives a bound on τ 0 (X i ) that is linear in v 0 and independent of H (apart from setting the above condition The (positive) numerical factor α = −η −1 ln(2e −η − 1) yields for instance α ≃ 2.11 for η = 1/10, and converges to α = 2 for η → 0.
The factor L − λ 0 (≤ L) in the above bound and condition, is simply the length of a reduced curve C , the (still spatially geodesic) part of C parametrized by λ ∈ [λ 0 , L], i.e., the part of C joining the point of coordinate x µ (λ 0 ) to P .It is in fact clear that this bound may be written in terms of the length of the shortest among all curves on Σ 0 that join P to a point of Σ 0 ∩ Σ ′ 0 .Since this simply amounts to a redefinition (if necessary) of C and of its length L, we may simply rewrite the above bound as where the previous assumed bound on v 0 (≪ 1) and L becomes: where L is now interpreted as the length of the shortest curve on Σ 0 joining P to Σ 0 ∩ Σ ′ 0 as above.This length may be significantly smaller than that of a geodesic on Σ 0 joining P to a given arbitrary point P 0 ∈ Σ 0 ∩ Σ ′ 0 in cases where the Σ 0 and Σ ′ 0 slices intersect multiple times (e.g., periodically).The symmetric case τ 0 (X i ) < 0, corresponding to Σ ′ 0 lying to the past of the V -comoving observer at P , can be handled in a similar way.One may first assume that τ 0 (λ) < 0 ∀λ ∈]0, L] with τ 0 (λ = 0) = 0 up to a suitable redefinition of C , L (and the parameter λ) as above.The bound , along with Eq. (4.23) (writing this time Θ µν L µ L ν ≥ − H), can then again be used to bound the derivative of τ 0 (λ) as given by Eq. (4.17), this time from below (cf.Eqs.(4.22) and (4.24)): Proceeding similarly to the τ 0 (X i ) > 0 case above, this results in still assuming v 0 HL ≤ η < ln(2), with the same numerical factor α = −η −1 ln(2e −η − 1) as above.In this bound, L is to be interpreted in the same way as discussed for Eqs.(4.29)-(4.30)above.Combining Eqs.(4.29), for the τ 0 (X i ) > 0 case, and (4.32), for the τ 0 (X i ) < 0 case, implies that in all cases 14 under the condition (4.30), with α = α(η) as above and still the same interpretation for L as discussed above.τ 0 (X i ) corresponds to the distance between the slices Σ 0 and Σ ′ 0 as measured along the flow line of V passing through P .The spatial curve length L depends on the point P ∈ Σ 0 ∩T B 0 considered, but it may itself be globally bounded for all points P in this domain by a finite length L. The latter may for instance be defined as the diameter of the domain within Σ 0 , L = L1 ≡ max or as the maximal distance along Σ 0 of a point in the domain to the intersection of the two slices within the domain, L = L2 ≡ max Above, d(P, P ′ ) is the spatial distance along Σ 0 between the points P and P ′ , that is, the proper length of the shortest curve on Σ 0 joining P and P ′ .The existence of the above maxima is guaranteed by the compactness of the spatial domain Σ 0 ∩ T B 0 (and consequently of Σ 0 ∩ T B 0 ∩ Σ ′ 0 , which is non-empty).L2 is always a smaller (or equal) length than L1 , and may be much smaller, but it could be harder to determine in practice, and it generally depends on Σ ′ 0 .We assume that the bounds on the tilt, v 0 ≪ 1, on the expansion rate, H, and on the distances on the reference slice as defined by L, can be taken as small enough to ensure a relatively small v 0 H L. We denote as η the best upper limit that may be set a priori on this value, and as above we assume that this limit is strictly smaller than ln(2).We thus have and consequently the condition (4.30) is satisfied with this same η for all points P ∈ Σ 0 ∩T B 0 .Eq. (4.33) then implies that for all such points, of coordinates (τ = 0, X i ), This global bound means in particular that, within the spacetime tube T B 0 delineating the domain of interest, the distance between the slices Σ 0 and Σ ′ 0 (along V ) is everywhere much smaller than the spatial size of the domain as measured in the reference slice Σ 0 , if v 0 ≪ 1.
Case of a nonvanishing 4-acceleration / velocity dispersion.In the more general case, the 4acceleration a of V cannot be everywhere neglected.When for instance V represents the 4-velocity of a source fluid, a can correspond to non-gravitational accelerations, which we would consider negligible apart from extremely small fractions of the volume of the domain considered, but it can also arise from the effective modelling of velocity dispersion within the source fluid by effective pressure terms.We can then assume some global bound ā (with dimension time −1 or length −1 ) on √ a • a = b µν a µ a ν and follow a similar derivation as above from Eq. (4.17), noting that |G(τ 0 (λ), λ)| ≤ ā τ 0 (λ) 0 exp ( Hτ ) dτ = ā/ H exp ( Hτ 0 (λ)) − 1 .For a small enough ā so that ā L is at most of order unity and under a possibly more stringent, ā-dependent constraint on v 0 HL than previously, bounds on |τ 0 | now involving ā L can be similarly obtained.
For instance, let us assume that ζ(ā L − v 0 H L) ≤ K for some constant K > 0, where ζ is the nonnegative, nondecreasing function defined by ζ(x) ≡ (e x − 1)/x.The inverse function ζ −1 can be defined and is nondecreasing as well, and the above condition is equivalently rewritten as Note that the above requirements imply again that η < ln 2, and thus K > 1/(2 ln 2).The prefactor in the above expression is always larger than 2 (corresponding to its Kη ≪ 1 limit), but remains of order unity if 2Kη is not too close to 1. Consequently, for the above bound on |τ 0 | to be significant, K should not be very large, and ā L should accordingly be at most of order unity (with ζ −1 (K) ∼ ln K for large K).For instance, ā L ≃ 1 would still allow for taking K ≲ 2. Note that one can take e.g.
In practice, the above requirements may be too restrictive, or the bound too large, to be used directly. 15This can be considered as being due to accumulating proper-time differences along Σ 0 when 4-accelerations have a specific consistent spatial orientation (corresponding to making the Cauchy- • a in G into an equality, with no changes of sign).However, especially when V models a physical fluid 4-velocity, and regardless of a corresponding to non-gravitational forces or to effective velocity dispersion, a residual 4-acceleration is only expected around localized overdensities, hence small parts of the domain.The acceleration vector is moreover expected to have radial orientations around the centers of those overdensities, leading to compensating signs in L • a over paths crossing such regions.In such a case, we may thus assume that the fluid's 4-acceleration only contributes small corrections to the above (a = 0) bounds.It can moreover be argued that the small spatial regions where a may be non-negligible could be avoided by the spatial path C up to a small increase in its total length, while the cases where the endpoint P of the path, spanning the whole integration domain, falls within such a dense region may be neglected for their small contribution to the volume-weighted spatial integrals I(S) in either slice.When V arises from a different, more geometric construction, we shall simply assume it to be a geodesic vector field for simplicity.Accordingly, in the following, we will adopt either of these simplifying assumptions and use the bounds obtained in the vanishing-acceleration case above -with small corrections being e.g.encompassed into taking a slightly pessimistic value for η or v 0 , if necessary.
The range τ ∈ [0, τ 0 (X i )] for the maximum simply corresponds geometrically to taking a maximum over the segment of a flow line of V that joins the two slices Σ 0 and Σ ′ 0 .The interval [0, τ 0 (X i )] 15 Considering for instance, in the late Universe, the natural case of V representing the 4-velocity field of a nonrelativistic (nearly dust) matter fluid, the main contribution to its 4-acceleration over most of the volume is expected to arise from velocity dispersion within the fluid acting as effective pressure forces which contribute to oppose gravitational collapse in bound structures.The magnitude of such a 4-acceleration from effective pressure can then be estimated as being at most of the same order as the (Newtonian) gravitational acceleration within those virialized domains.A typical present-day value for these accelerations in the outskirts of galaxies or within galaxy cluster haloes, for instance, would then be of the order of 10 −10 m.s −2 (e.g., [41,42]), or a little smaller.This value of ā accumulated over a large cosmological domain with a diameter L of a Hubble length, would correspond to ā L of about 10 −1 .A more pessimistic estimate of the maximum ā of the 4-acceleration local amplitude, e.g.accounting for the central regions of dense clusters or for the haloes of very massive elliptic galaxies, and still with L ≃ 1/H0, may then lead to an ā L of order unity or beyond.The regions where this may occur, however, would occupy a very small volume fraction at such scales and may accordingly be avoided by the spatial path and/or neglected in the volume-weighted spatial integrals considered.
(which should be read as [τ 0 (X i ), 0] in case τ 0 (X i ) < 0) is part of the larger interval [−α v 0 L, α v 0 L] from Eq. (4.33), and we may then use to get rid of the remaining dependence of Eq. (4.41) on the specific foliation F ′ .Setting again as η the best upper threshold that may be set on the global v 0 H L, and assuming that it can obey the constraint (4.36), η < ln(2), we can moreover replace the two X i -dependent factors L in Eq. (4.41) by the global length bound L. With these two remarks, inserting the above Eq.(4.41) into Eq.(4.10) provides the following bound on the variation of the spatial integral of the scalar S between the two slices under the condition (4.36): Above, we have kept the X i -dependent distance L in the second spatial integral, but it may as well be replaced by the global size L (since ) to give a simpler, though a priori weaker bound.We also note that, from Eq. (4.9), we could instead have written and thus the second spatial integral in Eq. (4.43) could alternatively be replaced by: , or by: where we have omitted the prefactor α(1 + η α) v 0 L in front of this integral.Another variant of the above bound, Eq. (4.43), is provided in Appendix A.2 as Eq.(A.7).It is based on an alternative writing of ψ(X i ) and requires some knowledge about the local variable ∇ µ (SV µ ), rather than S and V µ ∂ µ S as in the above.
The bound obtained here, Eq. (4.43), ensures that |∆I(S)|/I(|S|) A A 0 ≪ 1 provided the 'tiltweighted' domain size v 0 L (≪ L) can be considered much smaller than the characteristic lengths (or times) associated with H−1 and with S (0,X i ) / |(V µ ∇ µ S) (τ,X i ) |.At least in the cases where S does not change sign on Σ 0 (e.g., for an energy or mass density, or for S = 1, or as a known property of Σ 0 ), this expresses a direct constraint on the relative variation of its integral, |∆I(S)/I(S) A A 0 |.In general, this bound, as well as its alternative form in Eq. (A.7), require that the global upper limit H on the local directional expansion rates exists over the appropriate spacetime region, and that it and/or the domain size and tilt velocity be sufficiently small for the condition (4.36) to be satisfied.As mentioned, the bounds become more stringent if one can even ensure that v 0 H L ≪ 1 (i.e., η ≪ 1).In a cosmological setup, and in the case where V can be associated with the 4-velocity of a matter fluid source component, one can expect the typical expansion rate within the domain to be of the order of the Hubble parameter at a characteristic time of Σ 0 .If local fluctuations of the expansion rate do not substantially exceed this value, then v 0 H L ≪ 1 is ensured up to a domain size of the order of the associated Hubble length (and possibly larger for v 0 small enough).This assumption might be violated in rapidly collapsing, overdense regions.However, as discussed earlier about the 4-acceleration, which is also expected to only have potentially non-negligible values in strong overdensities, such regions typically occupy very small spatial volumes and may accordingly be avoided and neglected in the spatial integrals, allowing for a tighter H bound in the remaining domain.This may as well be applied to the tilt velocities (interpreted as peculiar velocities of the source fluid in each slice), considering that a tighter threshold may be imposed on the tilts everywhere outside small overdense regions.This would allow for an even smaller v 0 to be picked while neglecting the contributions of the remaining overdensities as small corrections to the result.
As for the variation of integrals ∆I(S) in Sec.4.4.3above, these bounds are expressed as a function of v 0 H L and require this factor to be small enough to obey the constraint (4.36); hence, similar remarks on setting H and/or v 0 apply.The bounds on the variation ∆⟨S⟩ still depend on the local evolution rate V µ ∇ µ S or non-conserved current ∇ µ (SV µ ); but we note in this case that their additional dependence is on the local fluctuations of S on the reference slice, rather than just ⟨S⟩ A A 0 .

Bounds for constant-proper time foliations
In cosmology, the age of the Universe is typically measured with respect to the proper time of fundamental observers.However, the choice of proper time function is associated with a calibration freedom, and a family of proper time functions in general exist for a given 4-velocity field.In the present section, we shall investigate bounds of integrals formulated within classes of proper time foliations.In particular, we give as an example proper time foliations that are calibrated within the epoch of last scattering.

The family of constant-proper time foliations
We consider proper time foliations of a single 4-velocity field, as an example of a set of foliations that are natural to compare.Let us consider the class of proper time functions of a given 4-velocity field.We shall use this (unit) 4-velocity field as our volume measure vector V , and as above use this V to define the propagation of the boundaries of the tube T B 0 .We say that τ is a proper time function of V if it satisfies the equation where ξ is a given function satisfying the transport rule V µ ∇ µ ξ = 0.If the function ξ, describing the distance of a given solution τ to the reference time function τ ref , can be bounded on a single hypersurface, then it can be bounded throughout the tube T B 0 .In the following, we shall consider bounds that are relevant for when we can assume that we can bound ξ on the initial surface Σ init such that |ξ| ≤ δT everywhere on Σ init .It then follows immediately from the above transport rule that |ξ| ≤ δT globally within T B 0 .This scenario is thus substantially simpler than the setup considered in subsection 4.4 above where the main difficulty was bounding the proper time distance between the two slices considered.With such a bound already ensured, we do not need to assume small tilts between the slices of the constant-proper time foliations.

Resulting bounds on the norm of ∆I(S)
We now consider an arbitrary scalar function S, and two intersecting leaves Σ 0 and Σ ′ 0 of the foliations F and F ′ corresponding to the level sets of τ ref and τ ′ , where τ ′ is an arbitrary member of the class of solutions (4.52).We may choose τ ref (Σ 0 ) = 0 = τ ′ (Σ ′ 0 ) without loss of generality, using the gauge freedom of shifting τ ′ by an additive constant if necessary.Let (X i ) again be a set of three V -comoving spatial coordinates.The values τ 0 (X i ) taken by τ ref on the Σ ′ 0 hypersurface are equivalently given by the values of −ξ on that same surface since τ ref = τ ′ − ξ; and it follows that |τ 0 (X i )| ≤ δT within T B 0 .From this global upper bound on |τ 0 (X i )| we can formulate an upper bound on |ψ(X i )| from Eq. (4.9), which reads, in the (τ ref , X i ) coordinate system: Similarly to subsection 4.4, the exponential terms above may as well be replaced by affine expressions in H δT as e 3 H δT ≤ 1 + η −1 (e 3η − 1) H δT by considering an upper threshold value η that can be set on H δT .The above local constraints on |ψ(X i )| then bound the variation of the integral of S, ∆I(S) = I(ψ(X i )), as |∆I(S)| ≤ I(|ψ(X i )|).This translates as well into a bound on the variation of averages, ∆⟨S⟩, along the same lines as in Sec.4.4.4above.These bounds are useful when it is possible to constrain (V µ ∂ µ S) (τ,X i ) over the bounding time span ±δT between the leaves.When δT is much smaller than typical values of the time scales set by within the subdomain of T B 0 defined by −δT ≤ τ ref ≤ δT , then |∆I(S)|/I(|S|) A A 0 is much smaller than 1.This directly constrains the relative variation |∆I(S)/I(S) A A 0 | at least when S has a constant sign along the reference slice.The full width at half maximum of the visibility function in ΛCDM cosmology is ∼ 10 5 years, and can be used to define the duration of the last scattering epoch.Consider the sub-class of proper time functions (4.52) of u which can all be initialized within the epoch of last scattering.That is, we set τ ref as defined from an initialization hypersurface Σ init that is fully contained within the spacetime region of last scattering, and all remaining proper time functions considered are also required to have an ("initial") level set contained within the region of last scattering.These proper time functions accordingly satisfy the global constraint |ξ| ≤ δT with δT ∼ 10 5 years.We might think of this class of synchronisation hypersurfaces as defining a set of equally preferred calibrations of a cosmic age function.
Let us investigate how integral quantities computed at the present epoch, defined as the τ = t 0 slice in each foliation for the fixed present-day age t 0 , differ between foliations of this set.Let typical values of the matter fluid's expansion rate |∇ µ u µ | around the present-epoch Universe be of the order of 3H 0 with the Hubble constant H 0 ∼ 10 −10 years −1 , and let its local fluctuations reach at most a factor of a few times this value (say, at most ∼ 10H 0 ), such that

Summary and discussion
The 3 + 1 foliation of spacetime in general relativity is a powerful tool for casting Einstein's equations into an initial value problem and for constructing coarse-grained variables by integration operations defined on the leaves.This is in particular used in approaches for cosmological averaging, as in [13][14][15][16][17][18][19][20][21].It has been noted that average quantities that appear in such formalisms have dependence on the foliation within which they are formulated [23][24][25][26][27][28][29].
In the present paper, we go beyond these previous investigations and consider the 3 + 1 foliation problem in a broad context.We are hence not addressing a particular metric model nor a perturbative setting, but are treating the problem of foliation dependence in relativistic integrals and averages of scalar variables over arbitrary bounded regions of 3-dimensional hypersurfaces, in presence of a generic nonsingular spacetime metric.This setup also encompasses the case of integrals and averages over the entire hypersurfaces when those have a closed topology, such as that of a 3-torus or a 3-sphere.We view the integral and averaged quantities as functionals of the general scalar function defining the foliation.Our systematic analysis of these functionals considering infinitesimal variations of foliations revealed that globally foliation-independent functionals do exist but must be associated with a locally-conserved current.Thus, the only physically relevant exactly preserved quantities are total rest masses within a bounded domain or other quantities that are per construction preserved within the individual volume elements.This is not in contradiction with the gauge-independent scheme recently proposed in [29] in the context of perturbation theory, since there, the coordinate transformations considered are not affecting the spacetime foliation, and thus actual foliation changes are by construction absent in the authors' scheme.
We additionally examined choosing foliations specifically to leave a certain functional invariant under infinitesimal deformations, i.e., foliations or individual slices which extremize a given functional, as in the well-known case of extremal-volume slices.Such extremals provide examples of ways to uniquely specify a foliation and thus to eliminate the ambiguity of foliation choice.We briefly discussed, as specific applications, the extremization of entropy functionals and the selection of slices with a minimal average tilt with respect to a given vortical (hence not hypersurface-orthogonal) fluid flow.
Since strictly foliation-independent integral/average quantities are rare, as might have been expected, in the second half of our paper we investigated the bounding of foliation dependence of integrals and averages under finite changes of foliations.There we showed that we can in some cases set upper limits on this dependence if we consider foliations with space-like leaves that can be bounded in terms of their relative distance.In particular we have considered classes of space-like foliations with associated normal vectors that have a small relative tilt.We have also considered families of constant-proper time foliations (for a given family of observers defining a 4-velocity field) that are bounded in terms of distance between their respective initial synchronization slices.In both cases, we derived bounds on the relative variation of integrals of scalar quantities between different foliations, under the assumption that the local volume expansion rate also remains bounded within a certain spacetime region.These bounds depend on the ability to set constraints on the evolution rate of the scalar quantity considered or on the associated non-conserved current.This implies, as a special case, bounds on the volume within the integration domain.The results on integral functionals of scalars also directly imply bounds on the variation of the associated averages, which we explicited in the small-tilt case, and which depend on the local fluctuations of the corresponding scalar.These bounds are consistent with the qualitative discussion for similarly restricted classes of space-like foliations in [28].
The foliation dependence of cosmological backreaction terms as defined in [13,14] and extended to arbitrary spatial foliations in [15,28] turns out to be complicated to bound rigorously in general.This is due to the appearance of the threading lapse M ≡ (V • ∇A) −1 = (∇A • ∇A) −1/2 γ −1 V ,n , or of M 2 , as a factor in the corresponding integrands, inducing a dependence on the local foliation scalar's gradient 16 .We have accordingly not derived general analytical expressions for the bounding of variations of backreaction terms.We note, however, that our results in the case of constant-proper time foliations (V • ∇A = 1) are directly applicable to the averages and backreaction terms appearing in [15,28] when the same class of foliations is considered there.Our results in the case of spacelike foliations with small relative tilts also remain applicable to such backreaction terms under the additional assumption of geodesic slicings, for which the slicing lapse N ≡ (∇A • ∇A) −1/2 can be set to 1 -up to the small extra corrections induced by the presence of the Lorentz factor γ V ,n ∼ 1.Moreover, since we have succeeded in bounding the volume of spatial sections, and since the important implication of backreaction is precisely its impact on the growth of cosmic volume over time, our results still provide implicit bounds on the foliation-dependent contributions to backreaction.In a setting where the foliation dependence of the volume is tightly bounded at all times, while the foliation dependence of a particular backreaction term might be larger, the backreaction effect on the volume (as the combined impact of all backreaction terms accumulated over some time evolution) also has to remain tightly bounded in this scenario.
Our results indicate that while the averaged properties of a given region of spacetime are generally going to depend on the reference time-slicing, there are nevertheless tight bounds that can be constructed within physically motivated classes of foliations that are close to each other by some suitable measure.In particular, for cosmological purposes where most of the matter in the Universe is thought to have non-relativistic relative speeds, physically meaningful space-like foliation frames that approximately trace the matter of the Universe will have a relative tilt velocity much smaller than unity.There are also epochs in the early Universe that constitute natural choices for setting a synchronisation of relevant cosmological foliations.Since, for instance, the time duration of the epoch of last scattering is very short relative to the present-day Hubble time scale, a synchronisation within this epoch provides a class of natural foliations which have small separations, compared to the characteristic time scale of expansion.The results that we have obtained are thus directly applicable to the averaging problem in cosmology.
.51) Let Σ init be an initial reference hypersurface chosen at convenience, and let τ = τ ref be the proper time function satisfying τ ref (Σ init ) = constant.All solutions to Eq. (4.51) can be expressed as .53) Assume that we can define a global upper bound 3 H on the volume expansion rate, time-averaged along V for τ ref spanning [−δT, 0] or [0, δT ]:

4. 5 . 3
Example: proper time foliations synchronised near the last scattering epoch We consider a class of proper time foliation scalars (4.52) of a physical matter congruence V = u which are synchronised at the epoch of last scattering.The epoch of last scattering defines a natural initialisation epoch in cosmology, and this epoch may thus naturally be used for synchronizing the proper time of physical observers.This epoch does however cover a spacetime region with a finite width in cosmic time.In the physics of recombination, the visibility function quantifies the probability of the streaming of photons as protons and electrons combine into hydrogen.The visibility function might then be used to quantify an epoch of last scattering, after which most photons are freely streaming.

±δTτ
=0 |∇ µ u µ | (τ,X i ) dτ ≤ 3 H δT ≲ 10 −4 , with 3 H ≲ 10H 0 .Let us see for instance how the present-epoch volume of the Universe is bounded within members of this set of proper time foliations.We accordingly take S(τ, X i ) = 1, and the bound (4.53) becomesψ(X i ) ≤ exp max ±δT ±δT τ =0 |∇ µ V µ | (τ,X i ) dτ − 1 ≲ exp 3 H δT − 1 ≈ 3 H δT ≲ 10 −4 .(4.56)Thus, the difference in the present-epoch volume V 0,τ ≡ I(1) τ t 0 , as measured in two members of the class of proper time foliations synchronised at the last scattering epoch, is bounded as|∆V 0 | = |∆I(1)| ≤ I(|ψ(X i )|) τ ref t 0 ≲ 3 H δT I(1) τ ref t 0 = 3 H δT V 0,τ ref ≲ 10 −4 × V 0,τ ref .(4.57)Hence, the volume of the present-epoch Universe is well-defined at the level of ∼ 10 −4 in the class of proper time foliations calibrated within the region of last scattering.This is due to the time scale set by the duration of last scattering through the width of the visibility function being much smaller than the typical time scale associated with volume expansion in the present-epoch Universe.
.54) This may arise still as a consequence of 3 H holding as a bound on the local expansion rate everywhere over the spacetime range considered, as in subsection 4.4 above; but in the present case only the less restrictive time-averaged bound is required.With this assumption, Eq. (4.53) becomes