On average properties of inhomogeneous fluids in general relativity III: general fluid cosmologies

Abstract

We investigate effective equations governing the volume expansion of spatially averaged portions of inhomogeneous cosmologies in spacetimes filled with an arbitrary fluid. This work is a follow-up to previous studies focused on irrotational dust models (Paper I) and irrotational perfect fluids (Paper II) in flow-orthogonal foliations of spacetime. It complements them by considering arbitrary foliations, arbitrary lapse and shift, and by allowing for a tilted fluid flow with vorticity. As for the first studies, the propagation of the spatial averaging domain is chosen to follow the congruence of the fluid, which avoids unphysical dependencies in the averaged system that is obtained. We present two different averaging schemes and corresponding systems of averaged evolution equations providing generalizations of Papers I and II. The first one retains the averaging operator used in several other generalizations found in the literature. We extensively discuss relations to these formalisms and pinpoint limitations, in particular regarding rest mass conservation on the averaging domain. The alternative averaging scheme that we subsequently introduce follows the spirit of Papers I and II and focuses on the fluid flow and the associated $1+3$ threading congruence, used jointly with the $3+1$ foliation that builds the surfaces of averaging. This results in compact averaged equations with a minimal number of cosmological backreaction terms. We highlight that this system becomes especially transparent when applied to a natural class of foliations which have constant fluid proper time slices.

Appendices

A $3+1$ evolution equations along the congruence of the fluid

The evolution equations for $h_{ij}$ and ${\mathcal {K}}^i_{ j}$ along the congruence of the fluid are obtained from expressions (3.1) and (3.2), by relating the derivative $\partial _t |_{x^i}$ to ${\mathrm{d}}/{\mathrm{d}}t$ with the help of (2.25). They read:

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}}t} h_{ij} =&- 2 N {\mathcal {K}}_{ij} + N_{i || j} + N_{j || i} + V^k \partial _k h_{ij} \ ; \end{aligned}$$

(A.1)

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}}t} {\mathcal {K}}^i_{ j} =&\; N \, \Big ( {\mathcal {R}}^i_{ j} + {\mathcal {K}}{\mathcal {K}}^i_{ j} + 4 \pi G \, \big [ \left( S - E \right) \delta ^i_{ j} - 2\, S^i_{ j} \big ] - {\Lambda }\,\delta ^i_{ j} \Big ) \nonumber \\&- N^{|| i}_{\ \; || i} + N^k {\mathcal {K}}^i_{ j || k} + {\mathcal {K}}^i_{ k} N^k_{\ || j} - {\mathcal {K}}^k_{ j} N^i_{\ || k} + V^k \partial _k {\mathcal {K}}^i_{ j} \ . \end{aligned}$$

(A.2)

Comoving coordinate system. In the comoving picture, as described in Sect. 2.4, we have ${\varvec{N}} = N {\varvec{v}}$, or equivalently ${\varvec{V}} = {\varvec{0}}$. Equations (A.1) and (A.2) hence read:

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}}t} h_{ij} =&- 2 N {\mathcal {K}}_{ij} + N ( v_{i || j} + v_{j || i}) + v_{i} N_{ || j} + v_{j} N_{ || i} \ ; \end{aligned}$$

(A.3)

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}}t} {\mathcal {K}}^i_{ j} =&\; N \, \Big ( {\mathcal {R}}^i_{ j} + {\mathcal {K}}{\mathcal {K}}^i_{ j} + 4 \pi G \, \big [ \left( S - E \right) \delta ^i_{ j} - 2\, S^i_{ j} \big ] - {\Lambda }\,\delta ^i_{ j} \Big ) - N^{|| i}_{\ \; || i} \nonumber \\&+ N v^k {\mathcal {K}}^i_{ j || k} + N {\mathcal {K}}^i_{ k} v^k_{\ || j} + {\mathcal {K}}^i_{ k} v^k N_{|| j} - N {\mathcal {K}}^k_{ j} v^i_{\ || k} - {\mathcal {K}}^k_{ j} v^i N_{|| k} \ . \end{aligned}$$

(A.4)

For these equations, focusing on the local foliation-related variables ${\mathbf {h}}$, ${\mathcal {K}}^i_{\ j}$ rather than fluid variables, the additional assumption of a fluid proper time foliation and of the choice $t=\tau$ (implying $N=\gamma$), corresponding to the Lagrangian picture, does not affect the comoving picture equations (A.3)–(A.4) above, except for the possibility of replacing each occurrence of N by $\gamma$. We shall accordingly not rewrite the equations for this case.

The local equations of evolution along the fluid flow for arbitrary slices and coordinates (A.1)–(A.2) allow for an alternative derivation of the coordinate-time derivative of the extrinsic volume, Eq. (3.17):

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}}t} {\mathcal {V}}_{\mathcal {D}}= \int _{{\mathcal {D}}_{{\varvec{x}}}} \left( -N {\mathcal {K}}+ \left( N v^i \right) _{|| i} \right) \sqrt{h} \, {\mathrm{d}}^3 x \ , \end{aligned}$$

(A.5)

by starting over from (3.14) and expanding its integrand as

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}}t} {\mathcal {V}}_{\mathcal {D}}= \int _{{\mathcal {D}}_{{\varvec{x}}}} \left( \frac{1}{2} h^{ij} \frac{\mathrm{d}}{{\mathrm{d}}t} h_{ij} + J^{-1} \frac{\mathrm{d}}{{\mathrm{d}}t} J \right) \sqrt{h} \, {\mathrm{d}}^3 x \ . \end{aligned}$$

(A.6)

The trace of (A.1) can then be used together with (3.8) to obtain:

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}}t} {\mathcal {V}}_{\mathcal {D}}= \int _{{\mathcal {D}}_{{\varvec{x}}}} \left( - N {\mathcal {K}}+ N^i_{\ || i} + \frac{1}{2} h^{ij} V^k \partial _k h_{ij} + \partial _k V^k \right) \sqrt{h} \, {\mathrm{d}}^3 x \ . \end{aligned}$$

(A.7)

This expression then allows to catch up with the end of the derivation given in Sect. 3.2.2, so that a similar use of relations (3.16) and (2.10) again gives the evolution of the volume (A.5).

B Extrinsic averaging operator in fluid-intrinsic variables

The system of equations for extrinsic averages on ${\mathcal {D}}$ derived in Sect. 3 (Theorem 1) is mostly expressed in terms of geometric variables of the ${\varvec{n}}$-orthogonal hypersurfaces, such as their extrinsic curvature. We here present an alternative formulation of the same equations focusing instead on the intrinsic, rest-frame kinematic and dynamical quantities of the fluid (see Sect. 2.2.2) which do not depend on the foliation choice.

We can first rewrite the volume expansion rate (3.19) and the commutation rule (3.22) in terms of the intrinsic local expansion rate of the fluid $\varTheta$ by re-expressing the three-divergence of $N {\varvec{v}}$ as

$$\begin{aligned} \big ( N v^i \big )_{|| i} = N v^i_{\ || i} + v^i N_{|| i} = N \nabla _\mu v^\mu \ , \end{aligned}$$

(B.1)

where we have employed (2.6) for the last equality. Noticing that ${\mathcal {K}}= - \nabla _\mu n^\mu$ and making use of expression (2.7), we get:

$$\begin{aligned} - N {\mathcal {K}}+ \big ( N v^i \big )_{|| i} = \frac{N}{\gamma } \varTheta - \frac{1}{\gamma } \frac{{\mathrm{d}} \gamma }{{\mathrm{d}} t}= {\tilde{\varTheta }}- \frac{1}{\gamma } \frac{{\mathrm{d}} \gamma }{{\mathrm{d}} t}:= {\tilde{\varTheta }}^T \ , \end{aligned}$$

(B.2)

where we have defined a tilted and scaled expansion rate${\tilde{\varTheta }}^T$ out of the scaled rate ${\tilde{\varTheta }}= (N/\gamma ) \, \varTheta$. The factor $N / \gamma$ in ${\tilde{\varTheta }}^T$ above adjusts the local clock rates between the proper time of the fluid and the coordinate time. This can also be seen upon writing:

$$\begin{aligned} {\tilde{\varTheta }}^T = \frac{N}{\gamma } \varTheta - \frac{1}{\gamma } \frac{\mathrm{d} \gamma }{{\mathrm{d}}t} = \frac{N}{\gamma } \left( \varTheta - \frac{1}{\gamma } \frac{\mathrm{d} \gamma }{{\mathrm{d}}\tau } \right) = \frac{{\mathrm{d}}\tau }{{\mathrm{d}}t} \left( \varTheta - \frac{1}{\gamma } \frac{\mathrm{d} \gamma }{{\mathrm{d}}\tau } \right) \ , \end{aligned}$$

(B.3)

where we have used the relation (2.26) between $\mathrm {d}/\mathrm {d}t$ and $\mathrm {d}/\mathrm {d}\tau$. The additional tilt term $- \gamma ^{-1} \, \mathrm {d}\gamma / \mathrm {d}t$ can be understood as the effect of the evolving mutual tilt (Lorentz boost) between the hypersurfaces in which ${\mathcal {D}}$ is embedded, and the fluid flow. This affects the local volume measure so that the evolution of the volume is not only due to the fluid’s intrinsic expansion.

The above rewriting (B.2) allows us to recast the volume and scale factor evolution rates, and the commutation rule, respectively, into the following expressions:

$$\begin{aligned} \frac{3}{a_{\mathcal {D}}} \frac{{\mathrm{d}} a_{\mathcal {D}}}{{\mathrm{d}} t}&{} = \frac{1}{{\mathcal {V}}_{\mathcal {D}}} \frac{{\mathrm{d}} {\mathcal {V}}_{\mathcal {D}}}{{\mathrm{d}} t} = \left\langle {{\tilde{\varTheta }}^T} \right\rangle _{{\mathcal {D}}} \; ; \end{aligned}$$

(B.4)

$$\begin{aligned} \frac{\mathrm{d}}{{\mathrm{d}}t} \left\langle {\psi } \right\rangle _{{\mathcal {D}}}&{} = \left\langle { \frac{\mathrm{d}}{{\mathrm{d}}t} \psi } \right\rangle _{{\mathcal {D}}} - \left\langle { {\tilde{\varTheta }}^T } \right\rangle _{{\mathcal {D}}} \left\langle {\psi } \right\rangle _{{\mathcal {D}}} + \left\langle { {\tilde{\varTheta }}^T \psi } \right\rangle _{{\mathcal {D}}} \; . \end{aligned}$$

(B.5)

We notice that, even for the general configuration we are investigating (see Fig. 1), the commutation rule, as well as the domain volume expansion rate, can be cast into a simple form with respect to the fluid quantities, although this extrinsic averaging framework requires the explicit additional contribution from the evolving tilt.

The use of the Raychaudhuri equation (4.14) and the energy constraint (4.16) (instead of the scalar parts of the extrinsic $3+1$ Einstein equations (3.2), (3.4)), together with the above alternative form of the commutation rule, allows for a rewriting of the evolution equations for the extrinsic effective scale factor $a_{\mathcal {D}}$. This yields the following equivalent formulation of Theorem 1, in terms of rescaled fluid-intrinsic kinematic and dynamical variables, ${\tilde{\sigma }}^2 = (N^2/\gamma ^2) \, \sigma ^2$, ${\tilde{\omega }}^2 = (N^2/\gamma ^2) \, \omega ^2$, $\tilde{{\mathscr {R}}}= (N^2/\gamma ^2) \, {\mathscr {R}}$, ${\tilde{\epsilon }}= (N^2/\gamma ^2) \, \epsilon$, ${\tilde{p}} = (N^2/\gamma ^2) \, p$, and $\tilde{{\mathcal {A}}} = (N^2/\gamma ^2) \, {\mathcal {A}}$ (with ${\mathcal {A}}= \nabla _\mu a^\mu$), as well as ${\tilde{{\Lambda }}} := (N^2/\gamma ^2) \,{\Lambda }$:

Corollary 1a

(extrinsically averaged evolution equations in fluid variables) The averaged evolution equations for the extrinsic effective scale factor $a_{\mathcal {D}}$ can also be written under the following form:

$$\begin{aligned} 3 \left( \frac{1}{a_{\mathcal {D}}} \frac{{\mathrm{d}} a_{\mathcal {D}}}{{\mathrm{d}} t} \right) ^2= & {} 8 \pi G \left\langle {{\tilde{\epsilon }}} \right\rangle _{{\mathcal {D}}} + \left\langle {{\tilde{{\Lambda }}}} \right\rangle _{{\mathcal {D}}} - \frac{1}{2} \left\langle {{\tilde{{\mathscr {R}}}}} \right\rangle _{{\mathcal {D}}} - \frac{1}{2} {\tilde{{\mathcal {Q}}}}^{\mathrm{T}}_{\mathcal {D}}\ ; \end{aligned}$$

(B.6)

$$\begin{aligned} \frac{3}{a_{\mathcal {D}}} \frac{{\mathrm{d}}^2 a_{\mathcal {D}}}{{\mathrm{d}}t^2}= & {} -4 \pi G \left\langle {({\tilde{\epsilon }} + 3 {\tilde{p}})} \right\rangle _{{\mathcal {D}}} + \left\langle {{\tilde{{\Lambda }}}} \right\rangle _{{\mathcal {D}}} + {\tilde{{\mathcal {Q}}}}^{\mathrm{T}}_{\mathcal {D}}+ {\tilde{{\mathcal {P}}}}^{\mathrm{T}}_{\mathcal {D}}\ , \end{aligned}$$

(B.7)

with alternative, ‘tilted’ kinematical and dynamical backreactions, reading respectively:

$$\begin{aligned} {\tilde{{\mathcal {Q}}}}^{\mathrm{T}}_{\mathcal {D}}&:= \frac{2}{3} \left[ \left\langle {\left( {\tilde{\varTheta }}^T \right) ^2} \right\rangle _{{\mathcal {D}}} - \left\langle {{\tilde{\varTheta }}^T} \right\rangle _{{\mathcal {D}}}^2 \right] - 2 \left\langle {{\tilde{\sigma }}^2} \right\rangle _{{\mathcal {D}}} + 2 \left\langle {{\tilde{\omega }}^2} \right\rangle _{{\mathcal {D}}} + \frac{2}{3} \left\langle { 2 \, {\tilde{\varTheta }}^T \frac{1}{\gamma } \frac{{\mathrm{d}} \gamma }{{\mathrm{d}} t}+ \left( \frac{1}{\gamma } \frac{{\mathrm{d}} \gamma }{{\mathrm{d}} t}\right) ^2} \right\rangle _{{\mathcal {D}}} \; ; \end{aligned}$$

(B.8)

$$\begin{aligned} {\tilde{{\mathcal {P}}}}^{\mathrm{T}}_{\mathcal {D}}&:= \left\langle {\tilde{{\mathcal {A}}}} \right\rangle _{{\mathcal {D}}} + \left\langle {\frac{\gamma }{N} \frac{{\mathrm{d}} }{{\mathrm{d}} t} \left( \frac{N}{\gamma } \right) {\tilde{\varTheta }}^T } \right\rangle _{{\mathcal {D}}} - \left\langle { 2 \, {\tilde{\varTheta }}^T \frac{1}{\gamma } \frac{{\mathrm{d}} \gamma }{{\mathrm{d}} t}+ \left( \frac{1}{\gamma }\frac{{\mathrm{d}} \gamma }{{\mathrm{d}} t}\right) ^2 + \frac{N}{\gamma } \frac{{\mathrm{d}} }{{\mathrm{d}} t} \left( \frac{\gamma }{N} \frac{1}{\gamma } \frac{{\mathrm{d}} \gamma }{{\mathrm{d}} t} \right) } \right\rangle _{{\mathcal {D}}} .\nonumber \\ \end{aligned}$$

(B.9)

We recall that, as in Theorem 1a, the left-hand sides in the above equations should be seen as derivatives with respect to the chosen parameter t, and be interpreted according to the physical meaning of the latter. In particular, the term $3 \, a_{\mathcal {D}}^{-1} \, \mathrm {d}^2 a_{\mathcal {D}}/ \mathrm {d}t^2$ in equation (B.7) corresponds in a Lagrangian picture to the proper time scale factor acceleration, but not in general.

Under this form, only two backreaction terms appear, ${\tilde{{\mathcal {Q}}}}^{\mathrm{T}}_{\mathcal {D}}$ and ${\tilde{{\mathcal {P}}}}^{\mathrm{T}}_{\mathcal {D}}$, as the tilt only contributes under these combinations. These backreaction terms will not in general be directly related to the terms ${\mathcal {Q}}_{\mathcal {D}}$ and ${\mathcal {P}}_{\mathcal {D}}$ appearing in Theorem 1a, as they do not collect the same local terms in their expressions. They do coincide, however, for a fluid-orthogonal foliation assuming an irrotational fluid (with in this case ${\tilde{{\mathcal {Q}}}}^{\mathrm{T}}_{\mathcal {D}}= {\mathcal {Q}}_{\mathcal {D}}$ and ${\tilde{{\mathcal {P}}}}^{\mathrm{T}}_{\mathcal {D}}= {\mathcal {P}}_{\mathcal {D}}$, while ${\mathcal {T}}_{\mathcal {D}}= 0$).

Note that there is no explicit non-perfect fluid contribution to these effective evolution equations, although the non-perfect fluid components of the energy-momentum tensor do have an influence on the dynamics via the local (and average, see below) evolution of $\epsilon$ and p.

As before, this set of equations goes along with an integrability condition, and must be complemented by the evolution equation for the averaged energy density and pressure.

Corollary 1b

(integrability and energy balance conditions to Corollary 1a) The integrability condition corresponding to Corollary 1a reads:

$$\begin{aligned}&\frac{{\mathrm{d}} }{{\mathrm{d}} t} {\tilde{{\mathcal {Q}}}}^{\mathrm{T}}_{\mathcal {D}}+ \frac{6}{a_{\mathcal {D}}} \frac{{\mathrm{d}} a_{\mathcal {D}}}{{\mathrm{d}} t} \, {\tilde{{\mathcal {Q}}}}^{\mathrm{T}}_{\mathcal {D}}+ \frac{{\mathrm{d}} }{{\mathrm{d}} t} \left\langle {{\tilde{{\mathscr {R}}}}} \right\rangle _{{\mathcal {D}}} + \frac{2}{a_{\mathcal {D}}} \frac{{\mathrm{d}} a_{\mathcal {D}}}{{\mathrm{d}} t} \left\langle {{\tilde{{\mathscr {R}}}}} \right\rangle _{{\mathcal {D}}} + \frac{4}{a_{\mathcal {D}}} \frac{{\mathrm{d}} a_{\mathcal {D}}}{{\mathrm{d}} t} \, {\tilde{{\mathcal {P}}}}^{\mathrm{T}}_{\mathcal {D}}\nonumber \\&\quad = 16 \pi G \left( \frac{{\mathrm{d}} }{{\mathrm{d}} t} \left\langle {{\tilde{\epsilon }}} \right\rangle _{{\mathcal {D}}} + \frac{3}{a_{\mathcal {D}}} \frac{{\mathrm{d}} a_{\mathcal {D}}}{{\mathrm{d}} t} \left\langle {{\tilde{\epsilon }} + {\tilde{p}}} \right\rangle _{{\mathcal {D}}} \right) + 2 \, \frac{{\mathrm{d}} }{{\mathrm{d}} t} \left\langle {{\tilde{{\Lambda }}}} \right\rangle _{{\mathcal {D}}} \ , \end{aligned}$$

(B.10)

while the associated averaged conservation equation for the scaled energy density ${\tilde{\epsilon }}$ and pressure ${\tilde{p}}$ becomes:

$$\begin{aligned}&\frac{{\mathrm{d}} }{{\mathrm{d}} t} \left\langle {{\tilde{\epsilon }}} \right\rangle _{{\mathcal {D}}} + \frac{3}{a_{\mathcal {D}}} \frac{{\mathrm{d}} a_{\mathcal {D}}}{{\mathrm{d}} t} \left\langle {{\tilde{\epsilon }} + {\tilde{p}}} \right\rangle _{{\mathcal {D}}} = \left( \left\langle {{\tilde{\varTheta }}^T} \right\rangle _{{\mathcal {D}}} \left\langle {{\tilde{p}}} \right\rangle _{{\mathcal {D}}} - \left\langle {{\tilde{\varTheta }}^T {\tilde{p}}} \right\rangle _{{\mathcal {D}}}\right) - \left\langle {\frac{1}{\gamma } \frac{{\mathrm{d}} \gamma }{{\mathrm{d}} t} \, \left( {\tilde{\epsilon }}+ {\tilde{p}} \right) } \right\rangle _{{\mathcal {D}}} \nonumber \\&\quad - \left\langle {\frac{N^3}{\gamma ^3} \Big (a_\mu q^\mu + \nabla _\mu q^\mu + \pi ^{\mu \nu } \sigma _{\mu \nu } \Big )} \right\rangle _{{\mathcal {D}}} + 2 \left\langle {{\tilde{\epsilon }} \, \frac{\gamma }{N} \frac{{\mathrm{d}} }{{\mathrm{d}} t} \left( \frac{N}{\gamma } \right) } \right\rangle _{{\mathcal {D}}} . \quad \;\; \end{aligned}$$

(B.11)

(From the commutation rule, the expression $\langle {{\tilde{\varTheta }}^T}\rangle _{\mathcal {D}}\left\langle {{\tilde{p}}} \right\rangle _{{\mathcal {D}}} - \langle {{\tilde{\varTheta }}^T {\tilde{p}}}\rangle _{\mathcal {D}}$ can also be written as $\left\langle {{\mathrm{d}} {\tilde{p}} / \mathrm {d}t} \right\rangle _{{\mathcal {D}}} - \mathrm {d}\left\langle {{\tilde{p}}} \right\rangle _{{\mathcal {D}}}/\mathrm {d}t$.)

Table 1 Summary of the main differences between the various generalization proposals discussed in Sects. 3.5.4–3.5.6 for the scalar averaging of the $3+1$ Einstein equations in arbitrary foliations. This table is split into three parts, considering respectively the setup, the equations presented (and the corresponding effective Hubble parameter), and the main variables used

C Summarized literature comparison

We present in Table 1 a comparative overview of the various formalisms used in the existing generalization proposals of the system of averaged scalar equations of Papers I and II to general foliations, discussed in Sects. 3.5.4–3.5.6, including the extrinsic averaging formalism of Sect. 3 of this work.

In this table we express all notations in terms of those used in this work to make comparisons easier. In particular, when considering the 4-scalar expressions of [53, 74, 82], we define the lapse N as $(- \partial _\mu T \, \partial ^\mu T)^{-1/2}$, where T is the scalar function which labels the hypersurfaces. This quantity (noted ${\Gamma }$ in [74]) satisfies $n_\mu = - N \partial _\mu T$, thus playing an analogous role to the $3+1$ lapse, and it indeed coincides with the usual lapse if the 4-scalar formalism is split into a $3+1$ description with the natural choice of T as the time coordinate. Since the domain propagation varies between these papers, in terms of the present notations the averaging domain should be noted ${\mathcal {E}}_{{\varvec{\partial }}_t}$ or ${\mathcal {E}}_{{\varvec{n}}}$, and ${\mathcal {D}}$ for the present paper. Rather than indicating the corresponding domain for each average and backreaction term, we remove any domain-labelling subscript for these objects in the table for notational ease.

Additional specificities of some of the papers compared here:

• Tanaka and Futamase [87]: The commutation rule obtained in the corresponding framework is not explicitly provided or derived. It thus remains implicit as a necessary intermediate step for the derivation of the evolution equations for the effective scale factor from averages of the local dynamical equations. We note, moreover, that the main aim of the paper, as for [60], is to set up a link to a background-dependent perturbation approach.

• Räsänen [74]: The velocity field ${\varvec{u}}$ that is introduced in addition to ${\varvec{n}}$ is fully general and is not related either to ${\varvec{n}}$ nor to the content of the model universe (it could be chosen to be the 4-velocity of some fluid as in the present approach, but this would be a restriction of generality). It is supposed to represent the 4-velocity field of the observers. In the application paper [75], this field is restricted to be everywhere very close to ${\varvec{n}}$ (and so has a small vorticity), whereas ${\varvec{n}}$ is assumed to be chosen such that it builds hypersurfaces of statistical homogeneity and isotropy. These restrictions are already both suggested in the original paper [74] but the equations are kept general.

• Beltrán Jiménez et al. [8]: The main objective of using a general $T_{\mu \nu }$ in this paper is to account for theories beyond General Relativity whose differences are transferred into $T_{\mu \nu }$ as effective terms. Note also that this paper features an additional average equation giving the evolution of the average squared rate of shear $\partial _t \left\langle {\sigma ^2} \right\rangle$, as well as the corresponding local equation; these equations are absent from the other papers quoted in this Appendix, including the present work (the reason being that the resulting system is still not closed by adding this equation; work about looking deeper into the hierarchy of equations will be published in [15], where the next level of the hierarchy is explored through a topological approach).

• Smirnov [82]: Not only the choice of hypersurfaces (or of ${\varvec{n}}$) and the choice of the time that parametrizes them obey specific criteria in this paper, but this is also the case for the averaging domain, although this is not explicit in the averaged equations and it could as well be any domain propagated along the chosen ${\varvec{n}}$. Indeed, the domain is there chosen as a ‘sphere’ in some ${\varvec{n}}$-comoving coordinates $Z^i$ on the hypersurfaces, as defined by $H_{i j} Z^i Z^j \le r_0$ for some $r_0 > 0$ and with $H_{i j}$ the components of the spatial metric in these coordinates. This choice was a response to the series of papers of Gasperini et al. and Marozzi [52, 53, 69] to show how one may determine the boundary of the domain via a scalar function (here in the sense that the $Z^i$ are fixed a priori without any link to the actual spacetime coordinates choices).

• In the present work, we also introduce, in Sect. 4, a different (intrinsic) averaging formalism that measures scalar quantities and volume in the local rest frames of the fluid, even if they are integrated over a domain lying in the not necessarily fluid-orthogonal hypersurfaces. We then obtain the corresponding commutation rule and averaged dynamics under rather simple forms in terms of the intrinsic kinematic and dynamical quantities of the fluid —for instance, the effective Hubble parameter, still defined as 1/3 of the volume expansion rate, can be simply expressed as the average of $(N/\gamma ) \, (\varTheta /3)$— and only two backreactions, kinematical and dynamical, distinct in general from the terms ${\mathcal {Q}}_{\mathcal {D}}$, ${\mathcal {P}}_{\mathcal {D}}$ of Sect. 3. This formalism and this system of equations clearly differ from the literature compared in this Appendix (including our Sect. 3, although it otherwise follows the same setup), due to the different volume and averaging operator definition.

  We also give in the present work, in Appendix B, a re-expression of the averaged equations arising from the extrinsic averaging operator of Sect. 3, in terms of averages of fluid-intrinsic kinematic and dynamical variables. This is obtained at the expense of the appearance of additional contributions from the evolution of the tilt factor.

D Remarks on volume 3-forms and manifestly covariant rewritings

Let us consider the two different volume 3-forms obtained as the respective Hodge duals $\star {{{{\underline{\varvec{n}}}}}}$ , $\star {{{{\underline{\varvec{u}}}}}}$ to the 1-forms ${{{{\underline{\varvec{n}}}}}} = n_\mu \, {\mathbf {d}}x^\mu$ (metric-dual to the normal vector to the slices ${\varvec{n}}$) and ${{{{\underline{\varvec{u}}}}}} = u_\mu \, {\mathbf {d}}x^\mu$ (metric-dual to the fluid 4-velocity vector ${\varvec{u}}$). We can decompose these 3-forms as follows in the exact basis $({\mathbf {d}}x^\mu ) = ({\mathbf {d}}t,{\mathbf {d}}x^i)$, associated with the reference coordinates $(t,x^i)$:

$$\begin{aligned} \star {{{{\underline{\varvec{n}}}}}}&= \frac{1}{6} n^\mu \varepsilon _{\mu \nu \rho \sigma } \, {\mathbf {d}}x^\nu \wedge {\mathbf {d}}x^\rho \wedge {\mathbf {d}}x^\sigma \nonumber \\&= \frac{1}{6} \sqrt{h} \,\epsilon _{ijk} ({\mathbf {d}}x^i + N^i {\mathbf {d}}t) \wedge ({\mathbf {d}}x^j + N^j {\mathbf {d}}t) \wedge ({\mathbf {d}}x^k + N^k {\mathbf {d}}t) \; ; \end{aligned}$$

(D.1)

$$\begin{aligned} \star {{{{\underline{\varvec{u}}}}}}&= \frac{1}{6} u^\mu \varepsilon _{\mu \nu \rho \sigma } \, {\mathbf {d}}x^\nu \wedge {\mathbf {d}}x^\rho \wedge {\mathbf {d}}x^\sigma \nonumber \\&= \frac{1}{6} \sqrt{b} \,\epsilon _{ijk} ({\mathbf {d}}x^i - V^i {\mathbf {d}}t) \wedge ({\mathbf {d}}x^j - V^j {\mathbf {d}}t) \wedge ({\mathbf {d}}x^k - V^k {\mathbf {d}}t) \; , \end{aligned}$$

(D.2)

where $\varepsilon _{\mu \nu \rho \sigma }$ denote the components of the four-dimensional Levi-Civita tensor, while $\epsilon _{ijk}$ is simply the Levi-Civita symbol with 3 (spatial) indices. We have used above the components expressions (2.1) and (2.14) of ${\varvec{n}}$ and ${\varvec{u}}$, respectively, and we have related the modulus g of the 4-metric determinant to the spatial determinants h and bvia$h = g/N^2$ (see, e.g., [55, sec. 5.2.3], or [53, sec. 3]) and $b = \gamma ^2 h = (\gamma /N)^2 g$ (from Eq. (3.23)). Note that in the expression for $\star {{{{\underline{\varvec{n}}}}}}$, ${\mathbf {d}}x^i + N^i {\mathbf {d}}t$ is the hypersurface-projected element $h^i_{\ \mu }\, {\mathbf {d}}x^\mu$, related to the rewriting of the four-dimensional line element (2.5) in all generality as $\mathrm {d}s^2 = - N^2 \mathrm {d}t^2 + h_{ij} (\mathrm {d}x^i + N^i \mathrm {d}t) (\mathrm {d}x^j + N^j \mathrm {d}t)$.

The restrictions of these 3-forms to the tangent spaces to the spatial hypersurfaces thus read $( \star {{{{\underline{\varvec{n}}}}}} )_{|{\Sigma }} = \sqrt{h} \, \mathrm {d}^3 x$ and $( \star {{{{\underline{\varvec{u}}}}}} )_{|{\Sigma }} = \sqrt{b} \, \mathrm {d}^3 x$, respectively, since the integration element on the slices $\mathrm {d}^3 x$ previously introduced corresponds to $(1 / 6) \, \epsilon _{ijk} \, ({\mathbf {d}}x^i \wedge {\mathbf {d}}x^j \wedge {\mathbf {d}}x^k )_{|{\Sigma }}$. We thus get, for any scalar $\psi$, using here again the label ${}^h$ for the extrinsic volume and averages for a clearer distinction:

$$\begin{aligned} {\mathcal {V}}_{\mathcal {D}}^h \left\langle \psi \right\rangle _{\mathcal {D}}^h= & {} \int _{\mathcal {D}}\psi \, \sqrt{h} \,\mathrm {d}^3 x = \int _{\mathcal {D}}\psi \, \left( \star {{{{\underline{\varvec{n}}}}}} \right) _{|{\Sigma }} \; , \;\; \text {and,} \nonumber \\ {\mathcal {V}}_{\mathcal {D}}^b\left\langle {\psi } \right\rangle ^b_{{\mathcal {D}}}= & {} \int _{\mathcal {D}}\psi \,\sqrt{b} \,\mathrm {d}^3 x = \int _{\mathcal {D}}\psi \, \left( \star {{{{\underline{\varvec{u}}}}}} \right) _{|{\Sigma }} \; , \end{aligned}$$

(D.3)

the $\psi = 1$ case giving the volumes ${\mathcal {V}}_{\mathcal {D}}^h$ and ${\mathcal {V}}_{\mathcal {D}}^b$. We thus recover under this form the extrinsic and intrinsic operators of Sects. 3 and 4 as natural averages (and volumes) definitions based, respectively, on the hypersurface normal, or on the fluid 4-velocity.

The similar roles played by ${\varvec{n}}$ in the extrinsic averaging formalism and by ${\varvec{u}}$ in the intrinsic formalism can alternatively be highlighted by rewriting the corresponding volumes and averages under a manifestly 4-covariant form using a 4-scalar window function [53, 57]. The framework of Gasperini et al. [53] indeed gives a manifestly covariant form of the extrinsic volume and of the associated averaging operator as (for any scalar $\psi$):

$$\begin{aligned} {\mathcal {V}}_{\mathcal {D}}^h = \int _{{\mathcal {M}}^4} {\mathfrak {W}}_{\mathcal {D}}^h \, \sqrt{g} \, \mathrm {d}^4 x \quad ; \quad \left\langle \psi \right\rangle _{\mathcal {D}}^h = \frac{1}{{\mathcal {V}}_{\mathcal {D}}^h} \int _{{\mathcal {M}}^4} {\mathfrak {W}}_{\mathcal {D}}^h \, \psi \sqrt{g} \, \mathrm {d}^4 x \; , \end{aligned}$$

(D.4)

respectively. Here, ${\mathfrak {W}}_{\mathcal {D}}^h$ is the window function selecting the averaging domain:

$$\begin{aligned} {\mathfrak {W}}_{\mathcal {D}}^h(x^\alpha ) := n^\mu \nabla _\mu \big ({{\mathcal {H}}}(A(x^\alpha ) - A_0) \big ) {{\mathcal {H}}}(r_0 - B(x^\alpha )) \; , \end{aligned}$$

(D.5)

where ${\mathcal {H}}$ is the Heaviside step function. The scalar functions A, B respectively define the foliation and the four-dimensional tube spanned by the domain, with the parameters $A_0$, $r_0$ selecting respectively a specific slice as $\{ x^\alpha \vert \; A(x^\alpha ) = A_0 \}$ and the domain boundaries from the condition $B(x^\alpha ) \le r_0$. (In [53], a hypersurface-orthogonal propagation is then set for the averaging domain; this is done by requiring $n^\mu \partial _\mu B = 0$.) As pointed out and analyzed in more detail in [57], the fluid volume (4.5) and the intrinsic averager (4.6) can be written under the same form by simply replacing ${\varvec{n}}$ by ${\varvec{u}}$ in the window function, i.e., by replacing the window function ${\mathfrak {W}}_{\mathcal {D}}^h$ by a window function ${\mathfrak {W}}_{\mathcal {D}}^b$ in (D.4), with

$$\begin{aligned} {\mathfrak {W}}_{\mathcal {D}}^b(x^\alpha ) := u^\mu \nabla _\mu \big ({{\mathcal {H}}}(A(x^\alpha ) - A_0) \big ) {{\mathcal {H}}}(r_0 - B(x^\alpha )) \; . \end{aligned}$$

(D.6)

With the addition of the constraint $u^\mu \partial _\mu B = 0$, which defines a comoving domain propagation, this yields the same intrinsic averaging framework as that considered in our Sect. 4 above.

E Remark on the threading lapse

We here add a technical remark. In the fluid-intrinsic approach we can borrow one element from the $1+3$ formalism to foliate spacetime, the so-called spacetime threading, although spatial volume averaging only makes sense on hypersurfaces. We recall that in the $1+3$ decomposition, and in comoving spatial coordinates, the four-dimensional line element (2.33) reads (see, e.g., [59, sec. 10]):

$$\begin{aligned} {\mathrm{d}}s^2 = - {\mathcal {M}}^2 \; {\mathrm{d}}t^2 + 2 \, {\mathcal {M}}^2 {\mathcal {M}}_i \; {\mathrm{d}}t \, {\mathrm{d}}X^i + ( b_{ij} - {\mathcal {M}}^2 {\mathcal {M}}_i {\mathcal {M}}_j) \; {\mathrm{d}}X^i {\mathrm{d}}X^j , \end{aligned}$$

(E.1)

with ${\mathcal {M}}$ the threading lapse and $\underline{\varvec{{\mathcal {M}}}}$ the threading shift, which relate to the lapse and to the Eulerian velocity dual 1-form ${{{{\underline{\varvec{v}}}}}}$ as follows:

$$\begin{aligned} {\mathcal {M}}:= \frac{N}{\gamma } \;\; ; \qquad \underline{\varvec{{\mathcal {M}}}} := \frac{\gamma ^2}{N} \left( {{{{\underline{\varvec{v}}}}}} + v^k v_k \, {{{{\underline{\varvec{n}}}}}} \right) \; : \; {\mathcal {M}}{\mathcal {M}}_\mu = \gamma (0, v_i) = (0, u_i) . \end{aligned}$$

(E.2)

In the Lagrangian description we have:

$$\begin{aligned} {\mathcal {M}}= 1 \;\; ; \qquad \underline{\varvec{{\mathcal {M}}}} = \gamma \left( {{{{\underline{\varvec{v}}}}}} + v^k v_k \, {{{{\underline{\varvec{n}}}}}} \right) \; : \; {\mathcal {M}}_\mu = \gamma (0, v_i) = (0, u_i) . \end{aligned}$$

(E.3)

Note that in the intrinsic form of the general averaged equations (see Sect. 4.2), we only deal with occurrences of $N / \gamma ={\mathcal {M}}$.

F Erratum

We wish to point out a small mistake in Paper II [19] (repeated in the appendix of [22] after equations (A23) and (A28) therein). For this we recall the spatial components of the 4-acceleration, $a_i = N_{|| i} / N$, its 4-divergence ${\mathcal {A}} := \nabla _{\mu }a^{\mu } = a^i_{|| i} + a^i a_i$, and the expression of the latter in terms of the lapse N or the injection energy per fluid element h (related to the relativistic enthalpy),

$$\begin{aligned} {\mathcal {A}} = \left( \frac{N^{|| i}}{N}\right) _{|| i} + \frac{N^{||i} N_{||i}}{N^2} = \frac{N^{|| i}_{\ \ || i}}{N} = h \left( \frac{1}{h} \right) ^{||i}_{\ \ ||i} = - \frac{h^{||i}_{\ \ ||i}}{h} + 2 \frac{h^{||i} h_{||i}}{h^2}\;\;, \end{aligned}$$

which are correctly written (for such a fluid-orthogonal framework). However, the first equality in equation (10a) of Paper II (and of its review in [22]) is incorrect, ${{{\mathcal {A}}}} \ne \left( N^{|| i} /N \right) _{|| i}$, due to an omission of the $a^i a_i$ contribution to ${\mathcal {A}}$.

There is also an imprecise statement: in Paper II, in footnote 3, it is stated that for scalars the operator || amounts to a partial derivative. This statement is only true for spatial components (for a scalar, ${ }_{||i} = \partial _i$, but ${ }_{||0} {} = N^i \partial _i \ne \partial _t$; ${ }_{||0}$ was identically zero for scalars in Paper II due to the vanishing shift).