<figure><img src="image_1.png" /></figure>

HAL
open science

# Lagrangian theory of structure formation in relativistic cosmology. V. Irrotational fluids

Yong-Zhuang Li, Pierre Mourier, Thomas Buchert, David L. Wiltshire

<figure><img src="image_2.png" /></figure>

► To cite this version:

Yong-Zhuang Li, Pierre Mourier, Thomas Buchert, David L. Wiltshire. Lagrangian theory of structure formation in relativistic cosmology. V. Irrotational fluids. Physical Review D, 2018, 98 (4), pp.043507. 10.1103/PhysRevD.98.043507. hal-01827979

HAL Id: hal-01827979

https://hal.science/hal-01827979v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Lagrangian theory of structure formation in relativistic cosmology.
V. Irrotational fluids

Yong-Zhuang Li,¹ Pierre Mourier,² Thomas Buchert,²,\* and David L. Wiltshire¹

¹School of Physical & Chemical Sciences, University of Canterbury,
Private Bag 4800, Christchurch 8140, New Zealand

²Univ Lyon, Ens de Lyon, Univ Lyon1, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574,
F-69007 Lyon, France

<figure><img src="image_4.png" /></figure>

(Received 15 June 2018; published 7 August 2018)

We extend the general relativistic Lagrangian perturbation theory, recently developed for the formation of cosmic structures in a dust continuum, to the case of model universes containing a single fluid with a single-valued analytic equation of state. Using a coframe-based perturbation approach, we investigate evolution equations for structure formation in pressure-supported irrotational fluids that generate their rest-frame spacetime foliation. We provide master equations to first order for the evolution of the trace and traceless parts of barotropic perturbations that evolve in the perturbed space, where the latter describes the propagation of gravitational waves in the fluid. We illustrate the trace evolution for a linear equation of state and for a model equation of state describing isotropic velocity dispersion, and we discuss differences to the dust matter model, to the Newtonian case, and to standard perturbation approaches.

DOI: 10.1103/PhysRevD.98.043507

I. INTRODUCTION

Relativistic cosmological perturbation theory is based on evolving the Einstein equations with a global foliation of the spacetime metric, via the 3 + 1 formalism [1,2]. In the standard approach a spatially homogeneous and isotropic Friedmann-Lemaître-Robertson-Walker (FLRW) geometry is assumed as the unperturbed global background spacetime, and Einstein's equations are then solved to some order on this predefined background [3]. The standard approach is Eulerian in the sense that perturbations are represented and propagate on this background that corresponds, in the Newtonian limit, to Eulerian perturbation theory. In this latter case, a perturbation method for the density and velocity fields is used to solve the Euler-Poisson system of equations that governs the fluid evolution [4]. Cosmological structure formation in the nonperturbative regime is also generally modeled within the Newtonian framework.

An alternative approach to structure formation has also
been developed, principally in the Newtonian regime, which
is directly tied to fluid elements. It is consequently known as
Lagrangian perturbation theory [5–22] to distinguish it from
the Eulerian approach based on coordinates on an assumed
global background. The Lagrangian approach uses a single
perturbation variable, the fluid’s deformation field. This gives
it the advantage of also applying in the nonlinear regime,
where Eulerian density perturbations are large. In recent
years, Lagrangian perturbation theory has been generalized

to general relativistic cosmologies with a dust continuum
[23–26]; see also [27–33].

In the Newtonian regime, an extension of Lagrangian perturbation theory to fluids with dynamic pressure was considered first in terms of isotropic pressure [34]. The resulting Lagrangian perturbation equations have been solved up to second order for a polytropic fluid [35,36]. For third-order perturbative solutions in Newtonian Lagrangian perturbation theory with pressure, see Ref. [37]. Models with isotropic pressure can also be considered as phenomenological models for the generally anisotropic pressure originating from the velocity dispersion of dust particles [38–40], by taking velocity moments of the collisionless Vlasov equation [41,42]. For a sequence of modeling assumptions used in nonperturbative extensions of Lagrangian perturbation theory, see the summary [43].

In this paper we will extend relativistic Lagrangian
perturbation theory for a dust matter model [23–26] to
the case of irrotational perfect fluids, and also to cases that
are relevant for the modeling of multistream regimes where
the dust approximation breaks down. This will provide a
framework not only to deal with a relativistic generalization
of Newtonian Lagrangian perturbation theory with pressure
at late epochs, but also to the fully relativistic situation of
the early Universe.

A primary motivation for such an investigation is to estab-
lish a framework which is better suited to studies of the
backreaction of inhomogeneities in cosmology as compared
to standard perturbation theory. In particular, standard cos-
mological perturbation theory conventionally assumes that
average cosmic evolution is *exactly* described by a solution to

<footnote>*Corresponding author.
buchert@ens-lyon.fr</footnote>
Einstein's equations with a prescribed energy-momentum tensor on a global hypersurface irrespective of the scale of coarse-graining of the matter fields. No fundamental physical principle demands such an outcome [44].

The scalar averaging scheme introduced in [45–48] is an example of an approach to backreaction of inhomogeneities in cosmology, in which the Einstein equations are assumed to hold on small scales, where they are well tested, but not for the average cosmic evolution on arbitrarily large spatial scales. This is a consequence of the fact that a generic averaging operation includes nonlocal fluctuation terms, and it should not be confused with modified gravity approaches which change the Einstein-Hilbert action. A variety of phenomenological interpretations of the Buchert scheme are possible [49–58], since additional ingredients are required to relate statistical quantities to physical observables determined from our own cosmological observations.

To date, no phenomenological approach to backreaction has fully utilized the general scalar averaging framework for perfect fluids [46]. In the timescape scenario [51–53], solutions to the corresponding system of averaged scalar equations have been given with matter and radiation [59] extending smoothly into the early radiation-dominated epoch in the early Universe. However, in deriving these solutions it was assumed that backreaction is insignificant before photon-electron decoupling, so that backreaction involving pressure terms was neglected.

Neglecting backreaction in the primordial plasma may seem to be a reasonable approximation for the evolution of the background universe to leading order, given that it is extremely close to being spatially homogeneous and isotropic at early times. However, backreaction can nonetheless make a significant difference when considering the growth of perturbations. In particular, even if the difference from the Friedmann equation is of order $10^{-5}$ as a fraction of energy density at decoupling, this is nonetheless of the same order as the density perturbations. A recent study of cosmic microwave background (CMB) anisotropies in the timescape model found that neglecting such small differences in initial conditions at last scattering leads to systematic uncertainties of 8–13% for particular cosmological parameters at the present epoch [60]. This remark applies to the conservative assumption that the background universe does not already contain backreaction arising from earlier epochs that could be compatible with large-scale homogeneity and isotropy [57].

For these reasons, we desire a new approach to cosmological perturbation theory which is intrinsic to the fluid and not anchored to an embedding space. Relativistic Lagrangian perturbation theory represents a promising avenue, as it is intimately tied to physical particles. To proceed to a fully realistic theory will require important steps beyond those which we investigate in this paper. Such steps will include

(i) An extension from one fluid to the many fluids pertinent to the early Universe, which requires

considering a tilt between various fluid flow vectors and the normal to the spatial hypersurfaces<fnref n="1" />;

(ii) Identifying relevant physical scales and volume partitioning the model universe into regions whose average evolves by averaged dynamical equations, rather than by global Friedmann equations;

(iii) Aiming at a background-free description. While perturbations are still formulated in the present paper as deviations from a fixed background cosmology, a general volume partitioning can be implemented without referring to a background [54,62].

As a first step towards these goals, in the present paper we will firstly consider relativistic Lagrangian perturbation theory for the same system that was considered in Ref. [46], namely a single component perfect fluid with barotropic equation of state. We will also include an explicit cosmological constant term.

This paper is organized as follows. In Sec. II we employ a 3 + 1 formalism [1,2] with Lagrangian spatial coordinates, presenting the general framework and foliation structure for a general irrotational matter model. We then restrict our attention to a barotropic fluid and discuss in detail the fluid variables and their equation of state. In this context, in Sec. III we introduce Cartan's coframe formalism, proceeding with the relativistic Lagrangian perturbation approach. We develop the first-order Lagrangian scheme and derive master equations for the trace and trace-free parts of the perturbation field. In Sec. IV we apply the first-order Lagrangian scheme to particular matter models, allowing us to explicitly derive solutions for the trace part, and we illustrate and discuss the results. Particular solutions for the gravitoelectric traceless part are studied in the Appendix. We summarize our main results in Sec. V.

## II. SPACETIME FOLIATION STRUCTURE AND 3+1 EINSTEIN EQUATIONS

In this paper we will consider a model universe containing a single irrotational fluid, so that a foliation of spacetime into flow-orthogonal hypersurfaces can be introduced.

### A. Decomposition of Einstein's equations for flow-orthogonal hypersurfaces

The irrotationality assumption on the fluid amounts to the existence of two scalar functions, $N$ and $t$, such that the 1-form dual to the normalized 4-velocity vector $u^\mu$ of the fluid can be written as<fnref n="2" />

<footnote n="1">¹Note that in the standard approach, the same FLRW frame is used for different matter components. (Even in this idealized case there are important differences to be respected for the different matter components [61].)</footnote>

<footnote n="2">²In the convention we use here, greek letters $\mu, \nu, ...$ are spacetime indices running from 0 to 3, while lowercase latin letters $i, j, ...$ are spatial indices running from 1 to 3. We use units in which $c = 1$, if not otherwise stated.</footnote>
$$u_{\mu} = -N \partial_{\mu} t; \quad N := (-\partial^{\mu} t \partial_{\mu} t)^{-1/2}. \quad (1)$$

The level sets of $t$ then define flow-orthogonal hypersurfaces, labeled $\Sigma_t$, which foliate spacetime, with unit normal $u^\mu$, $u^\mu u_\mu = -1$. We will now follow the 3+1 formalism [1,2] and define our time coordinate as coinciding with this function $t$. In this case, $N(t, x^i)$ is the lapse function.

In addition, we choose the spatial coordinates to be spatial Lagrangian (or comoving) coordinates, denoted $X^i$, that are assumed to be constant along each flow line. In the set of coordinates $(X^\mu) = (t, X^i)$, the components of the fluid 4-velocity vector and its dual are then, respectively,

$$u^{\mu} = \frac{1}{N}(1, 0, 0, 0); \quad u_{\mu} = (-N, 0, 0, 0), \quad (2)$$

while the line element can be written as

$$ds^2 = g_{\mu\nu} dX^\mu dX^\nu = -N^2 dt^2 + g_{ij} dX^i dX^j. \quad (3)$$

Here, $g_{ij}$ corresponds both to the spatial coefficients of the 4-metric $g_{\mu\nu}$ and to the components of the 3-metric that it induces on the hypersurfaces $\Sigma_t$. Introducing the projector onto $\Sigma_t$, $h_{\mu\nu} = g_{\mu\nu} + u_\mu u_\nu$, this 3-metric is indeed

$$h_{ij} := g_{\mu\nu} h^\mu_i h^\nu_j = g_{ij}. \quad (4)$$

The spatial metric and the lapse function $N$ together encode the inhomogeneities. (We will later use the more elementary coframe coefficients instead of the 3-metric coefficients). We use $\mathcal{R}^i_j$ to denote the Ricci tensor coefficients of this spatial metric, with $\mathcal{R}$ the corresponding Ricci scalar.

Without loss of generality, the energy-momentum tensor of the fluid is given by

$$T_{\mu\nu} = (\epsilon + p)u_{\mu}u_{\nu} + pg_{\mu\nu} + \pi_{\mu\nu} + q_{\mu}u_{\nu} + q_{\nu}u_{\mu}, \quad (5)$$

where $\pi_{\mu\nu}$ is an anisotropic pressure, with $\pi_{[\mu\nu]} = 0$, $u^\mu \pi_{\mu\nu} = 0$ and $\pi^\mu_\mu = 0$, and $q_\mu$ the heat flux, with $q_\mu u^\mu = 0$.

Introducing the expansion tensor (as minus the extrinsic curvature) of the hypersurfaces,

$$\Theta_{ij} := \nabla_{\nu} n_{\mu} h^{\mu}_{i} h^{\nu}_{j} = \frac{1}{2N} \partial_{i} g_{ij}, \quad (6)$$

Einstein's equations with a cosmological constant may be cast into a set of constraint and evolution equations. The constraint equations are the energy and momentum constraints<fnref n="3" />:

$$\begin{gathered} \mathcal{R} + \Theta^2 - \Theta^i_j \Theta^i_j = 16\pi G\epsilon + 2\Lambda; \\ \Theta^i_{j\parallel i} - \Theta_{j\parallel} = -8\pi G q_j. \end{gathered} \quad (7)$$

The propagation equations are given by

$$\begin{gathered} \Theta^i_j = \frac{1}{2N} g^{ik} \partial_t g_{kj}; \\ N^{-1} \partial_t \Theta^i_j = -\Theta \Theta^i_j - \mathcal{R}^i_j + \mathcal{A}^i_j \\ + 4\pi G [(\epsilon - p)\delta^i_j + 2\pi^i_j] + \Lambda \delta^i_j, \end{gathered} \quad (8)$$

where $a_\mu := u^\nu \nabla_\nu u_\mu = N^{-1} N_{\parallel\mu}$ is the covariant acceleration of the fluid (with $\nabla$ denoting the 4-covariant derivative), and $\mathcal{A}^i_j := a^i_{\parallel j} + a^i_j a_j = N^{-1} N^i_{\parallel j}$. Combining the trace of the second equation with the energy constraint yields the Raychaudhuri equation:

$$N^{-1} \partial_t \Theta = -\frac{1}{3} \Theta^2 - 2\sigma^2 - 4\pi G(\epsilon + 3p) + \mathcal{A} + \Lambda, \quad (9)$$

$$\text{where } \mathcal{A} := \mathcal{A}^i_i = \nabla_\mu a^\mu = N^{-1} N^i_{\parallel i}.$$

With the spacetime described by the given metric, the energy-momentum conservation laws are expressed as follows:

$$\begin{align} \partial_t \epsilon + N \Theta(\epsilon + p) &= -N(q^\mu_{\parallel\mu} + 2q^\mu a_\mu + \sigma_{\mu\nu} \pi^{\mu\nu}); \tag{10} \\ (\epsilon + p)a_\mu + p_{\parallel\mu} &= -(\pi_{\mu\nu}^{\parallel\nu} + a^\nu \pi_{\mu\nu}) \nonumber \\ &\quad - \left(\frac{4}{3}\Theta q_\mu + q^\nu \sigma_{\mu\nu} + u^\nu \nabla_\nu q_\mu - q^\nu a_\nu u_\mu\right). \tag{11} \end{align}$$

In what follows, we will specialize to the case of isotropic pressure, $\pi_{\mu\nu} = 0$, and vanishing heat flux, $q_\mu = 0$. Note that with these assumptions we do still allow for some nonperfect fluids, since $p$ is not necessarily the local thermodynamic equilibrium pressure [63]. Such a restriction is required here since both extra terms in general create vorticity, which cannot be covered by the class of flow-orthogonal foliations considered in this work.

Let us illustrate this by considering more closely the irrotationality condition for a fluid with negligible heat flux, $q_\mu = 0$, to see how this condition constrains the equation of state and the anisotropic pressure. The vanishing of the vorticity 2-form implies vanishing of the antisymmetrized projected gradient of the acceleration, $a_{[\nu]\mu} = 0$, since $a_\mu = (\ln N)_{\parallel\mu}$ from (1), being a consequence of the existence of the fluid-orthogonal foliation. From this, one obtains through (11) the following constraint on the energy-momentum components:

$$\begin{align} \epsilon_{\parallel\mu} p_{\parallel\nu} + (\epsilon + p)_{\parallel\mu} h^{\rho}_{\nu} \nabla_{\sigma} \pi^{\sigma}_{\rho} \nonumber \\ - (\epsilon + p) h^{\rho}_{[\mu} h^{\sigma}_{\nu]} \nabla_{\rho} \nabla_{\tau} \pi^{\tau}_{\sigma} = 0. \tag{12} \end{align}$$

<footnote n="3">³The symbol $\parallel$ denotes the covariant derivative with respect to the 3-metric $h_{ij}$. When applied to scalars it reduces to a partial derivative, denoted $|$, with respect to the Lagrangian coordinates, $X^i$.</footnote>
Since $\nabla_{\mu}\pi^{\mu}_{\nu} = 0$ would imply the vanishing of the right-hand sides of (10)–(11), an anisotropic pressure that does contribute to the dynamics will satisfy $\nabla_{\mu}\pi^{\mu}_{\nu} \neq 0$ and thus will not fulfill the above condition in general, producing a vortical flow. Conversely, a barotropic fluid flow with $\pi_{\mu\nu} = 0$ and an effective equation of state of the form $p = \beta(\epsilon)$ automatically satisfies the above constraint. Moreover, for such a fluid, (11) allows one to write the acceleration as a flow-orthogonal projected gradient, and it will indeed obey the relativistic equivalent of the Kelvin-Helmholtz theorem, so that irrotationality will be preserved along the flow lines [63,64].

## B. Barotropic perfect fluid spacetimes

For the remainder of this paper we will only consider fluids with $q_{\mu} = 0$ and $\pi_{\mu\nu} = 0$. The energy-momentum tensor (5) then reduces to perfect fluid form:

$$T_{\mu\nu} = (\epsilon + p)u_{\mu}u_{\nu} + pg_{\mu\nu}, \quad (13)$$

while its conservation equations (10)–(11) become, respectively

$$\partial_t \epsilon + N\Theta(\epsilon + p) = 0; \quad (14)$$

$$a_{\mu} = - \frac{p_{\parallel \mu}}{\epsilon + p}. \quad (15)$$

As a further restriction we will assume that the fluid flow is barotropic; i.e., we assume a local relation of the form $p = \beta(\epsilon)$ to effectively hold throughout the entire fluid,<fnref n="4" /> that we will henceforth call the “equation of state” or EoS. As noted earlier, such a relation will ensure that the flow remains irrotational. For such a fluid, setting some reference constant energy and rest mass density values $\epsilon_1$, $q_1$, we may use the EoS to define a formal rest mass density $\rho(\epsilon)$ and a related specific enthalpy $h(\epsilon)$—as an injection energy per fluid element and unit formal rest mass [65]— respectively, by

$$\rho := F(\epsilon) := q_1 \exp \int_{\epsilon_1}^{\epsilon} \frac{dx}{x + \beta(x)}; \quad (16)$$

$$h(\epsilon) := \frac{\epsilon + \beta(\epsilon)}{F(\epsilon)} = \frac{\epsilon + p}{\rho}. \quad (17)$$

The energy-momentum conservation equations (10) and (11) then, respectively, provide a conservation law for $\rho$,

$$\partial_t \rho + N\Theta \rho = 0, \quad (18)$$

and a relation between the specific enthalpy (17) and the lapse,

$$\frac{N_{\parallel \mu}}{N} = a_{\mu} = - \frac{h_{\parallel \mu}}{h} : (Nh)_{|i} = 0. \quad (19)$$

By an appropriate choice of the hypersurface-labeling function $t$, the lapse can thus be rescaled so that [46,63]

$$N = \frac{1}{h} = \frac{F(\epsilon)}{\epsilon + \beta(\epsilon)}. \quad (20)$$

If we assume that the fluid remains in thermodynamic equilibrium locally, and if it has a nonvanishing rest mass density, then this density will follow the same evolution law (18) as $\rho = F(\epsilon)$, by rest mass conservation. This formal $\rho$ and the actual rest mass density will then coincide up to a possible different spatial dependence (cf., footnote 4). These two quantities may be made equal by a suitable choice of initial conditions for the rest mass density or local thermodynamic equilibrium assumptions.<fnref n="5" /> This would then ensure the validity of the interpretation of $\rho$ and $h$ as the physical rest mass density (or particle number density) and specific enthalpy of the fluid, respectively. We will not, however, make such assumptions in the following Sec. III, to keep its level of generality. This will allow us to consider the case of a zero rest mass fluid [for which $F(\epsilon) \neq 0$ and $h(\epsilon)$ are still well defined], as well as that of a nonzero rest mass density with less constrained initial conditions. It will also allow us to consider the variable $p$ as an effective pressure term—e.g., modeling velocity dispersion— instead of the local thermodynamic equilibrium pressure.

<footnote n="4">⁴Considering the local dynamical solution for these variables, there is always a freedom of integration constant that depends on the Lagrangian coordinates, i.e., on the particular fluid element. We assume here that the same relation holds for all fluid elements. Only this assumption makes the dynamical relation an apparent equation of state that is valid throughout the fluid flow. All related variables then also depend on this assumption, which is a restriction imposed on initial data.</footnote>

<footnote n="5">⁵Let us take the local state of the fluid to belong to a thermodynamic Gibbs space admitting the equation of state $u(s, v)$, where $s$ is the specific entropy, $v$ is the specific volume and $u = \epsilon v$ is the specific internal energy. If we now assume that $p$ is the local thermodynamic equilibrium pressure of the fluid, it can then be expressed as $p(s, v) = -\partial u/\partial v$. Provided that a specific equation of state does not render the above relations degenerate, then these relations may be inverted to provide $v(\epsilon, p)$. Within a barotropic flow satisfying $p = \beta(\epsilon)$, the actual rest mass density $v^{-1}$ thus only depends on the energy density $\epsilon$, which fully determines its initial conditions. From the conservation equations of both quantities, $\partial_t \epsilon/(\epsilon + \beta(\epsilon)) = -N\Theta = \partial_t (v^{-1})/v^{-1}$, this dependency must be $v^{-1} = F(\epsilon)$, for $\Theta$ not identically vanishing, up to a constant prefactor which can be absorbed in the choice of $q_1$. Hence, in this case, $F(\epsilon)$ is indeed the rest mass density of the fluid with no further loss of generality. Also note that under the same assumptions, $s$ is also a function of $\epsilon$, preserved along the flow lines as the flow is adiabatic: $\partial_t s = 0 = (ds/dc)\partial_t \epsilon$, while $\partial_t \epsilon$ is not identically vanishing. The flow is thus isentropic, $s$ being a constant $s_1$ that depends neither on time nor on the fluid element. The barotropic relation then corresponds to the equation $p(\epsilon, s)$ deduced from the thermodynamic equation of state, and taken at constant $s$, $\beta(\epsilon) = p(\epsilon, s = s_1)$ (see [63–68]).</footnote>
For the general treatment including these cases it will suffice to formally define $\varrho$ and $h$ from Eqs. (16) and (17) using the single barotropic assumption $p = \beta(\epsilon)$. We follow the notation of Ref. [46] here.

### III. LAGRANGIAN PERTURBATION SCHEME

In this section we will introduce the coframe formalism to describe spacetime, which is a set of four deformation 1-form fields dual to a generally noncoordinate basis of vectors at every point of the manifold [69–71]. A general relativistic version of a coframe-based perturbative approach for an irrotational dust continuum has been proposed in Ref. [28], developed further in Ref. [29] and in final form, featuring only the coframes as the single perturbation variable in Ref. [23].

#### A. Coframe formulation

Following [24–26], we construct a set of three spatial coframes $\eta^a$ such that the spatial metric can be rewritten in the form

$$^{(3)}\mathbf{g} = G_{ab} \eta^a \otimes \eta^b: g_{ij} = G_{ab} \eta^a{}_i \eta^b{}_j. \quad (21)$$

Here $G_{ab}(X)$ is the Gram matrix that encodes all the initial spatial metric perturbations, $G_{ab}(X) := \delta_a^i \delta_b^j G_{ij}(X)$, with the initial metric coefficients, $G_{ij}(X) := g_{ij}(t_i, X)$. On the other hand we can also include the temporal component into the matrix and rewrite it as

$$\tilde{G}_{\alpha\beta} = \begin{pmatrix} -1 & 0 \\ 0 & G_{ab} \end{pmatrix}. \quad (22)$$

With this we introduce a full set of four spacetime coframes $\eta^a$ to describe the 4-metric $^{(4)}\mathbf{g}$:

$$^{(4)}\mathbf{g} = \tilde{G}_{\alpha\beta} \eta^a \otimes \eta^b, \quad (23)$$

by defining the coframe components as

$$\eta^0_\mu = (-N, 0, 0, 0); \quad \eta^a_\mu = (0, \eta^a_i). \quad (24)$$

We now define the transformation between coordinate and noncoordinate bases as $\mathcal{J} = \sqrt{-g}/\sqrt{\tilde{G}} = \sqrt{-\bar{g}}/\sqrt{G}$ [the signature adopted here being $(-1, 1, 1, 1)$, and using the notation $g := \det(^{(4)}\mathbf{g})$, $\tilde{G} := \det(\tilde{G}_{\alpha\beta})$ and $G := \det(G_{ab})$]. This corresponds to $\mathcal{J} = -\det(\eta^a_\mu)$, or

$$\begin{aligned} & \frac{1}{4!} \epsilon_{\alpha\beta\gamma\delta} \eta^\alpha \wedge \eta^\beta \wedge \eta^\gamma \wedge \eta^\delta \\ & = -\frac{1}{4!} \mathcal{J} \epsilon_{\mu\nu\rho\sigma} \mathrm{d}X^\mu \wedge \mathrm{d}X^\nu \wedge \mathrm{d}X^\rho \wedge \mathrm{d}X^\sigma. \end{aligned} \quad (25)$$

From Eq. (24), in terms of the spatial components of the coframes, $\mathcal{J}$ becomes

$$\mathcal{J} = \frac{1}{3!} N \epsilon_{abc} \epsilon^{ijk} \eta^a{}_i \eta^b{}_j \eta^c{}_k = N \det(\eta^a{}_i), \quad (26)$$

while correspondingly, the dual vector basis can be described by the four frames $e_\alpha = e_\alpha^\mu \partial/\partial X^\mu$:

$$\begin{aligned} e_\alpha^\mu \eta^a_\nu &= \delta^\mu_\nu; & e_\alpha^\mu \eta^b_\mu &= \delta^b_\alpha; \\ e_\alpha^\mu &= -\frac{1}{6\mathcal{J}} \epsilon_{\alpha\beta\gamma\delta} \epsilon^{\mu\nu\rho\sigma} \eta^b_\nu \eta^\gamma_\rho \eta^\delta_\sigma; \\ e_a^i &= \frac{1}{2\mathcal{J}} N \epsilon_{abc} \epsilon^{ijk} \eta^b{}_j \eta^c{}_k; \\ e_0^\mu &= \frac{1}{N} (-1, 0, 0, 0); & e_a^\mu &= (0, e_a^i). \end{aligned} \quad (27)$$

With this choice, the evolution equations for $\mathcal{J}$ and the expansion tensor coefficients $\Theta^i{}_j$ read

$$\begin{aligned} \partial_t \mathcal{J} &= \frac{\partial_t N}{N} \mathcal{J} + \mathcal{J} N \Theta; \\ \Theta^i{}_j &= \frac{1}{2\mathcal{J}} \epsilon_{abc} \epsilon^{ikl} (\partial_t \eta^a{}_j) \eta^b{}_k \eta^c{}_l; \\ \frac{\partial_t \Theta^i{}_j}{N} &= -\Theta \Theta^i{}_j + \frac{1}{2\mathcal{J}} \epsilon_{abc} \epsilon^{ikl} \partial_t \left( \frac{1}{N} \partial_t \eta^a{}_j \right) \eta^b{}_k \eta^c{}_l \\ &\quad + \frac{1}{N \mathcal{J}} \epsilon_{abc} \epsilon^{ikl} (\partial_t \eta^a{}_j) (\partial_t \eta^b{}_k) \eta^c{}_l. \end{aligned} \quad (28)$$

From the constraint and evolution equations (7)–(9), together with the definition of $\mathcal{J}$ and Eq. (28), the Lagrange-Einstein system of an irrotational barotropic fluid model is cast into the following form:

$$\begin{aligned} G_{ab} \partial_t \eta^a{}_{[i} \eta^b{}_{j]} &= 0; \quad (29) \\ \frac{1}{2\mathcal{J}} \epsilon_{abc} \epsilon^{ikl} \partial_t \left( \frac{1}{N} (\partial_t \eta^a{}_j) \eta^b{}_k \eta^c{}_l \right) &= A^i{}_j - \mathcal{R}^i{}_j + [4\pi G(\epsilon - p) + \Lambda]\delta^i{}_j; \quad (30) \\ \epsilon_{abc} \epsilon^{ijk} (\partial_t \eta^a{}_i) (\partial_t \eta^b{}_j) \eta^c{}_k &= (16\pi G\epsilon + 2\Lambda - \mathcal{R})N\mathcal{J}; \quad (31) \end{aligned}$$

$$\left( \frac{1}{\mathcal{J}} \epsilon_{abc} \epsilon^{ikl} (\partial_t \eta^a{}_j) \eta^b{}_k \eta^c{}_l \right) \Big|_{\parallel i} = \left( \frac{1}{\mathcal{J}} \epsilon_{abc} \epsilon^{ikl} (\partial_t \eta^a{}_i) \eta^b{}_k \eta^c{}_l \right) \Big|_{\parallel j}; \quad (32)$$

$$p = \beta(\epsilon). \quad (33)$$

Equations (29)–(32) are not closed unless an EoS, here (33), is specified. Recall that the lapse appearing above can be replaced by its expression in terms of $\epsilon$, $N = (\epsilon + \beta(\epsilon))^{-1}F(\epsilon)$. The evolution equation (30) may be split into a trace part, which we then combine with the energy constraint (31) to obtain the Raychaudhuri equation, and a traceless part, yielding respectively

$$\frac{1}{2\mathcal{J}} \epsilon_{abc} \epsilon^{ikl} \partial_t \left( \frac{1}{N} \partial_t \eta^a{}_i \right) \eta^b{}_k \eta^c{}_l = \mathcal{A} - 4\pi G(\epsilon + 3p) + \Lambda; \quad (34)$$
$$ \frac{1}{2J} \left[ \epsilon_{abc} \epsilon^{ikl} \partial_i \left( \frac{1}{N} (\partial_t \eta^a_j) \eta^b_k \eta^c_l \right) - \frac{1}{3} \epsilon_{abc} \epsilon^{mkl} \partial_i \left( \frac{1}{N} (\partial_t \eta^a_m) \eta^b_k \eta^c_l \right) \delta^i_j \right] = \xi^i_j - \tau^i_j, \quad (35) $$

where $\tau^i_j := R^i_j - \frac{1}{3}R\delta^i_j$ are the coefficients of the traceless part of the spatial Ricci tensor, and $\xi^i_j := A^i_j - \frac{1}{3}A\delta^i_j$.

The Lagrange-Einstein system, Eqs. (29)–(33), is closed and provides the components $\eta^a_i$ of coframes, from which one can calculate the evolution of the perturbations. The system comprises 14 equations, where 9 equations describe the evolution for the coefficient functions of 3 spatial Cartan coframe fields, and the remaining 5 equations originate from the 4 constraints and the EoS defining the properties of the fluid.

## B. Perturbation ansatz

### 1. Background

We will choose a spatially flat, homogeneous and isotropic model universe as the background spacetime, with the same barotropic EoS, and including a possible constant curvature term into the first-order perturbations (cf., e.g., [25]). Accordingly, the spatial metric coefficients of the background will be $a^2(t)\delta_{ij}$, $a(t)$ being the background scale factor. We prescribe a homogeneous lapse $N_H(t)$ for this homogeneous and isotropic background, by setting its relation to the background energy density $\epsilon_H$, formal rest mass density $\rho_H := F(\epsilon_H)$ and pressure $p_H = \beta(\epsilon_H)$ as being the same relations as those for the inhomogeneous quantities,

$$ N_H = \frac{\rho_H}{\epsilon_H + p_H} = \frac{F(\epsilon_H)}{\epsilon_H + \beta(\epsilon_H)}. \quad (36) $$

We may then write the background line element as

$$ ds_H^2 = -N_H^2(t)dt^2 + a^2(t)\delta_{ij}dX^i dX^j. \quad (37) $$

Note that the evolution of the background lapse function $N_H(t)$ will be given by its definition (36) and the EoS, making it time-dependent for $p_H \neq 0$. One should keep in mind that our choice of time coordinate $t$ will consequently not coincide in general with the usual “cosmic time” coordinate for the background, and it will evolve at a different rate. The usual cosmic time $\tilde{t}$ would rather be defined by $d\tilde{t} = N_H(t)dt$, so that the background line element (37) would take the usual Friedmannian form for homogeneous and isotropic model universes$^6$:

$$ ds_H^2 = -d\tilde{t}^2 + a^2[\tilde{t}]\delta_{ij}dX^i dX^j. \quad (38) $$

With this time variable, the standard Friedmann equations would indeed be recovered:

$$ \begin{aligned} 3 \frac{\partial_t^2 a}{a} &= -4\pi G(\epsilon_H + 3p_H) + \Lambda; \\ 3 \left( \frac{\partial_t a}{a} \right)^2 &= 8\pi G \epsilon_H + \Lambda; \\ \partial_t \epsilon_H + 3 \frac{\partial_t a}{a} (\epsilon_H + p_H) &= 0. \end{aligned} \quad (39) $$

However, for consistency with the lapsed foliation used for the full inhomogeneous spacetime, in what follows we include the homogeneous lapse $N_H$ into the background and use the coordinate $t$. In terms of this variable, the acceleration and Friedmann equations are respectively

$$ \begin{aligned} \frac{3}{N_H^2} \frac{\partial_t^2 a}{a} &= -4\pi G(\epsilon_H + 3p_H) + \Lambda + 3 \frac{\partial_t a}{a} \frac{\partial_t N_H}{N_H^3}; \\ \frac{3}{N_H^2} \left( \frac{\partial_t a}{a} \right)^2 &= 8\pi G \epsilon_H + \Lambda, \end{aligned} \quad (40) $$

while the energy-momentum conservation equation is formally unchanged:

$$ \partial_t \epsilon_H + 3 \frac{\partial_t a}{a} (\epsilon_H + p_H) = 0. \quad (41) $$

### 2. Coframes decomposition

It is important to express the full set of equations in terms of a single perturbation variable, the coframes, so that the Lagrangian perturbation approach is well defined. Although this is not made fully explicit in the Lagrange-Einstein system (29)–(33), it is implicitly the case as the Ricci tensor and covariant derivatives are functionals of the metric, and hence of the coframes, and $\epsilon, p, N$ and $A^i_j$ can be expressed in terms of the coframes and initial energy density data. The latter relations are obtained via the conservation equation (18) for $\rho = F(\epsilon)$ and the evolution equation for $J := J/N = \det(\eta^a_i)$ from the first equation in (28):

$$ N\Theta = -\frac{\partial_t F(\epsilon)}{F(\epsilon)} = \frac{\partial_t J}{J}; \quad \epsilon = F^{-1}\left(\frac{F(\epsilon_i)}{J}\right), \quad (42) $$

where for any quantity $A$, $A_i$ denotes the quantity at initial time $t_i$. Here $J_i = 1$ as a result of the choice of initial conditions for the coframes. The barotropic EoS and choice of $N$ then allow us to determine $p, N$ and $A^i_j = N^{-1}N^{\parallel i}_{\parallel j}$, and to express these fields as functions of $J = \det(\eta^a_i)$.

We then follow the previous papers [23–26] and decompose the coframes into a FLRW coframe set and deviations thereof,

$$ \eta^a = \eta^a_i dX^i = a(t)(\delta^a_i + P^a_i)dX^i. \quad (43) $$

<footnote>$^6$The notation $a[\tilde{t}]$ signifies that the scale factor still takes the same values, $a[\tilde{t}] := a(t)$, but has a different functional dependence on the alternative time coordinate.</footnote>
At this nonperturbative level, the metric coefficients are then related to the deformation field by

$$g_{ij} = a^2(t)(G_{ij} + 2P_{(ij)} + G_{ab}P^a_i P^b_j), \quad (44)$$

where we have defined

$$P^i_j := \delta_a^i P^a_j; \quad P := P^k_k = \delta_a^k P^a_k; \quad P_{ij} := G_{ai} P^a_i. \quad (45)$$

Recall that the Gram matrix coefficients $G_{ab}$ have been defined to encode the initial metric inhomogeneities, so that the coefficients $P^a_i$ can be set to zero initially. Also recall that this coframe split is made with respect to a FLRW background with a nontrivial lapse included, and that the functional dependence of $a$, or of the deformation field, on the time coordinate $t$ will be affected accordingly.

We then expand the deformation fields $P^a_i$ into a perturbative sum, so that the coframes are given by

$$\eta^a = a(t) \left( \delta^a_i + \sum_{m=1}^{\infty} P^a_i^{(m)} \right) dX^i, \quad (46)$$

where the m $th$-order deformation field coefficients $P^a_i^{(m)}$ are of order $\epsilon^m$ for some bookkeeping parameter $\epsilon \ll 1$. In this paper we will only focus on first-order deformations.

### 3. Initial conditions

We will follow the steps of Refs. [25,26] to prescribe the initial data. The deformation field and its time derivatives are given at some initial time $t_i$ by

$$\begin{aligned} P^a_i(t_i) &= 0; \\ (\partial_t P^a_i)(t_i) &=: U^a_i; \\ (\partial_t^2 P^a_i)(t_i) &=: W^a_i - 2H_i U^a_i, \end{aligned} \quad (47)$$

where $H := \partial_t a/a$ is the Hubble function. Hereafter, we will normalize the scale factor as $a_i = 1$. The six 1-form fields $U^a = U^a_i dX^i$ and $W^a = W^a_i dX^i$ are 1-form generalizations of the initial Newtonian peculiar-velocity and peculiar-acceleration gradient fields, respectively.

The Lagrange-Einstein system with its split into trace and traceless parts according to (29)–(35) then translates into constraints on the initial data:

$$U_{[ij]} = 0; \quad W_{[ij]} = 0; \quad (48)$$

$$\begin{aligned} W - U \left( \frac{\partial_t N}{N} \right)_i = & 3H_i \left[ \left( \frac{\partial_t N}{N} \right)_i - \left( \frac{\partial_t N_H}{N_H} \right)_i \right] \\ & + \Lambda(N_i^2 - N_{Hi}^2) + N_i^2 A_i \\ & - 4\pi G \left[ (\epsilon_i + 3p_i) N_i^2 - (\epsilon_{Hi} + 3p_{Hi}) N_{Hi}^2 \right]; \end{aligned} \quad (49)$$

$$\begin{aligned} & {}^{\mathrm{t}}W^a_j \delta_a^i + \left( H_i - \left( \frac{\partial_i N}{N} \right)_i \right) {}^{\mathrm{t}}U^a_j \delta_a^i \\ & \qquad + U {}^{\mathrm{t}}U^a_j \delta_a^i - \left( U^a_k \delta_a^i U^b_j \delta_b^k - \frac{1}{3} U^a_l \delta_a^k U^b_k \delta_b^l \delta^i_j \right) \\ & = N_i^2 (\xi^i_j(t_i) - \tau^i_j(t_i)); \end{aligned} \quad (50)$$

$$\begin{aligned} & U^2 - U^a_i \delta_a^j U^b_j \delta_b^i + 4H_i U \\ & = 16\pi G (\epsilon_i N_i^2 - \epsilon_{Hi} N_{Hi}^2) + 2\Lambda(N_i^2 - N_{Hi}^2) - \mathcal{R}_i N_i^2; \end{aligned} \quad (51)$$

$$(N_i^{-1} U^a_j \delta_a^i)_{\parallel i} = (N_i^{-1} U)_{\parallel j} + 2H_i (N_i^{-1})_{\parallel j}; \quad (52)$$

$$p_i = \beta(\epsilon_i); \quad p_{Hi} = \beta(\epsilon_{Hi}). \quad (53)$$

The abbreviations $U := U^a_k \delta_a^k$, $W := W^a_k \delta_a^k$, and ${}^{\mathrm{t}}W^a_i := W^a_i - (1/3)W \delta^a_i$, ${}^{\mathrm{t}}U^a_i := U^a_i - (1/3)U \delta^a_i$ are used for the trace and traceless parts, respectively.

### C. First-order Lagrange-Einstein system

We now expand the above Lagrange-Einstein system and its initial conditions to first order<fnref n="7" /> in the only dynamical variable in this Lagrangian perturbation approach, namely the deformation field $P^a_i$. In what follows we omit the index (1) for the first-order deformation field and the associated initial conditions $U_{ij}$, $W_{ij}$, but keep the index for the other variables, as functionals of $P^a_i$. We first need to express these functionals explicitly at first order.

#### 1. Dependent variables at first order

In order to express the first-order Ricci tensor and scalar curvature in terms of the coframes, we expand the initial metric coefficients to first order as $G_{ij}(X) = \delta_{ij} + G_{ij}^{(1)}(X)$ since they reduce to $\delta_{ij}$ at the unperturbed zero-order level. Introducing the first-order quantities $G^{(1)ij} := \delta^{ik}\delta^{jl}G_{kl}^{(1)}$, $P^{ij} := \delta^{ik}\delta^{jl}P_{kl}$ for the inverse metric, we can then substitute the metric and its inverse, truncated to first order,

$$g_{ij} = a^2(\delta_{ij} + G_{ij}^{(1)} + 2P_{(ij)}); \quad (54)$$

$$g^{ij} = a^{-2}(\delta^{ij} - G^{(1)ij} - 2P^{(ij)}), \quad (55)$$

into the definitions of the spatial Christoffel symbols and of the spatial Ricci tensor. We thereby obtain

$$\begin{aligned} \Gamma^k_{ij}{}^{(1)} = & \frac{1}{2} \delta^{kl} (G^{(1)}_{il|j} + G^{(1)}_{lj|i} - G^{(1)}_{ij|l}) \\ & + \delta^{kl} (P_{(il)|j} + P_{(lj)|i} - P_{(ij)|l}); \end{aligned} \quad (56)$$

<footnote n="7">⁷Note that initial data can be assumed, without loss of generality, to be first order.</footnote>
$$ \mathcal{R}_{ij}^{(1)} = \mathcal{R}_{ij} + P_{[j}^{k}[k|i} + P_{[j}^{k}|k|i} + P_{(ik)}|j|^{k} - P_{(ij)}|^{k}|_{k}; \quad (57) $$

$$ \mathcal{R}^{(1)} = a^{-2}\mathcal{R} + 2a^{-2}(P^{ki}|_{i|k} - P^{k}|_{k}), \quad (58) $$

where $\mathcal{R}_{ij} := G_{i[k|j]}^{(1)|k} + G_{[j|k]|i}^{k(1)}$ and $\mathcal{R} := \delta^{ij}\mathcal{R}_{ij} = 2G_{[k|k]}^{l}|^{k(1)}$ are the initial conditions for the curvature tensor coefficients and their trace, respectively.

An important difference from the dust case is that here the spatial Ricci scalar will in general not be constrained to evolve as $\mathcal{R}(X)a(t)^{-2}$ at first order, due to the contributions from the lapse in the momentum constraints. As will be shown below, these contributions give rise to a nonzero evolution for the (initially vanishing) second term $(P^{ki}|_{i|k} - P^{k}|_{k})$, or equivalently a nonconserved scalar curvature, $\partial_t \mathcal{R}^{(1)} + 2H\mathcal{R}^{(1)} = a^{-2}\partial_t(a^2\mathcal{R}^{(1)}) \neq 0$, in contrast to the dust case.

Using the barotropic EoS and the corresponding solution (42) to the energy conservation equation (14), we can also expand $\epsilon$, $p$, $N$ and $\mathcal{A}_j^i$ in terms of the first-order deformation field. We write $\epsilon_i := \epsilon_{Hi}(1 + \delta\epsilon_i)$ at first order, and expand $J^{-1} = a^{-3}\det(\delta a_i^j + P^a_i)^{-1}$ at the same order as $a^{-3}(1 - P)$. The solution (42) for $\epsilon$ as a function of $J$ can then be expanded to first order in the perturbation as

$$ \begin{aligned} \epsilon &= F^{-1} \left( \frac{F(\epsilon_{Hi}) + F'(\epsilon_{Hi})\epsilon_{Hi}\delta\epsilon_i - F(\epsilon_{Hi})P}{a^3} \right) \\ &= F^{-1} \left( \frac{F(\epsilon_{Hi})}{a^3} \right) \\ &\quad + \left[ \frac{1}{a^3}\epsilon_{Hi}F'(\epsilon_{Hi})\delta\epsilon_i - P\frac{F(\epsilon_{Hi})}{a^3} \right] (F^{-1})' \left( \frac{F(\epsilon_{Hi})}{a^3} \right). \end{aligned} \quad (59) $$

The energy-momentum conservation equation (42) still holds for background quantities, giving

$$ F(\epsilon_H) = \frac{F(\epsilon_{Hi})}{a^3}. \quad (60) $$

This can be substituted into (59) to give

$$ \epsilon = \epsilon_H \left[ 1 + \frac{F(\epsilon_H)}{\epsilon_H F'(\epsilon_H)} \left( \frac{\epsilon_{Hi} F'(\epsilon_{Hi})}{F(\epsilon_{Hi})} \delta \epsilon_i - P \right) \right]. \quad (61) $$

The further use of the definition of $F$, Eq. (16), allows us to simplify the above to

$$ \epsilon = \epsilon_H \left[ 1 - \left( 1 + \frac{p_H}{\epsilon_H} \right) \bar{P} \right], \quad (62) $$

which we have written for convenience in terms of a shifted deformation trace,

$$ \bar{P} := P - \alpha_{Hi} \delta \epsilon_i, \quad (63) $$

where $\alpha_{Hi} := (\epsilon_{Hi} + \beta(\epsilon_{Hi}))^{-1}\epsilon_{Hi}$ is a constant, and $\delta\epsilon_i$ is the initial energy perturbation.

The pressure can in turn be expanded to first order as $p = \beta(\epsilon)$, yielding

$$ p = p_H - \beta'(\epsilon_H)(\epsilon_H + p_H)\bar{P}. \quad (64) $$

Note that the factor $\beta'(\epsilon_H)$ corresponds to the (generally time-dependent) dimensionless ratio of the background speed of sound to speed of light squared, $\beta'(\epsilon_H) := c_S^2(t)/c^2$, if $p_H$ is the thermodynamic equilibrium pressure for the background fluid.

We then expand the lapse $N = (\epsilon + p)^{-1}F(\epsilon)$ as

$$ N = N_H [1 + \beta'(\epsilon_H) \bar{P}] \quad (65) $$

at first order in the deformation field. At this order, one will then have (with $\partial_t P = \partial_t \bar{P}$)

$$ \begin{aligned} \frac{\partial_t N}{N} = \frac{\partial_t N_H}{N_H} &+ \beta'(\epsilon_H) \partial_t \bar{P} \\ &- 3H(\epsilon_H + \beta(\epsilon_H)) \beta''(\epsilon_H) \bar{P}, \end{aligned} \quad (66) $$

with

$$ \frac{\partial_t N_H}{N_H} = 3H\beta'(\epsilon_H). \quad (67) $$

This also allows one to obtain the first-order expression for $\mathcal{A}_j^i = N^{-1}N^{\parallel j}$:

$$ \mathcal{A}_j^i(1) = a^{-2}\beta'(\epsilon_H)\delta^{ik}\bar{P}_{|j,k}. \quad (68) $$

## 2. First-order system

Using the above expansions, the Lagrange-Einstein system (29)–(32) can be rewritten at first order in the deformation field as follows:

$$ \partial_t P_{[ij]} = 0; \quad (69) $$

$$ \begin{aligned} & \partial_t^2 P^i_j + 3H[1 - \beta'(\epsilon_H)]\partial_t P^i_j \\ & \qquad + H[1 - \beta'(\epsilon_H) - \mathcal{V}(t)]\partial_t \bar{P} \delta^i_j \\ & = N_H^2 \mathcal{A}_j^{(1)} - N_H^2 (\mathcal{R}_j^{(1)} - \frac{\mathcal{V}(t)}{4} \mathcal{R}^{(1)} \delta^i_j); \end{aligned} \quad (70) $$

$$ \partial_t(P^i_{j|i} - \bar{P}_{|j}) = -2H\beta'(\epsilon_H)\bar{P}_{|j}, \quad (71) $$

$$ \begin{aligned} & H \partial_t \bar{P} + 4\pi G[\epsilon_H + p_H - (2\epsilon_H + \tilde{\Lambda})\beta'(\epsilon_H)] N_H^2 \bar{P} \\ & = -\frac{1}{4} N_H^2 \mathcal{R}^{(1)}, \end{aligned} \quad (72) $$

with $\partial_t P = \partial_t \bar{P}$, and where $\mathcal{A}_j^{(1)}, \mathcal{R}_j^{(1)} = a^{-2}\delta^{ik}\mathcal{R}_{kj}^{(1)}$ and $\mathcal{R}^{(1)}$ are expressed as functions of $P^a_i$ according to the
formulas given above, $\tilde{\Lambda} := \Lambda/(4\pi G)$, and we introduce the abbreviation

$$ \mathcal{V}(t) := [\epsilon_H + p_H - (2\epsilon_H + \tilde{\Lambda})\beta'(\epsilon_H)]^{-1} \times \{\epsilon_H + p_H - (3\epsilon_H - p_H + 2\tilde{\Lambda})\beta'(\epsilon_H) + (2\epsilon_H + \tilde{\Lambda})(\epsilon_H + p_H)\beta''(\epsilon_H)\}. \quad (73) $$

Equation (70) has been obtained from the first-order expansion of the extrinsic curvature evolution equation (30) by combining it with the first-order energy constraint (72). The EoS (33) has already been used to expand $\epsilon, p$ and $N$ in terms of the first-order deformation field.

## D. First-order master equations

Following the approach of Ref. [26] the above system can be reexpressed by decomposing the deformation fields into trace, trace-free symmetric and antisymmetric parts:

$$ P^i_j = \frac{1}{3} P \delta^i_j + \Pi^i_j + \mathfrak{P}^i_j, \quad (74) $$

where $\Pi_{ij} = P_{(ij)} - \frac{1}{3}P\delta_{ij}$ and $\mathfrak{P}_{ij} = P_{[ij]}$.

We will now derive the governing equations for these parts, named “master equations.” For the trace part we use the new variable $\bar{P}$ from Eq. (63). Accordingly, (69)–(71) become

$$ \partial_t \mathfrak{P}_{ij} = 0: \quad \mathfrak{P}_{ij} = \mathfrak{P}_{ij}(t_i) = 0; \quad (75) $$

$$ \begin{aligned} \partial_t^2 \bar{P} + 3H[2 - 2\beta'(\epsilon_H) - \mathcal{V}(t)]\partial_t \bar{P} \\ = N_H^2 \mathcal{A}^{(1)} - N_H^2 \left(1 - \frac{3}{4}\mathcal{V}(t)\right) \mathcal{R}^{(1)}; \end{aligned} \quad (76) $$

$$ \partial_t^2 \Pi^i_j + 3H[1 - \beta'(\epsilon_H)]\partial_t \Pi^i_j = N_H^2 (\xi^i_j^{(1)} - \tau^i_j^{(1)}); \quad (77) $$

$$ \partial_t \left( \Pi^i_{j|i} - \frac{2}{3} \bar{P}_{|j} \right) = -2H\beta'(\epsilon_H)\bar{P}_{|j}. \quad (78) $$

Once again the first-order quantities $\mathcal{A}^{(1)}, \xi^i_j^{(1)}, \mathcal{R}^{(1)}$ and $\tau^i_j^{(1)}$ are used as shorthand notations but are meant to be expressed in terms of the deformation field. These expressions are obtained from the results above, Eqs. (56), (58), (68), as follows$^8$:

$$ a^2 \mathcal{A}^{(1)} = \beta'(\epsilon_H) \Delta_0 \bar{P}; \quad (79) $$

$$ a^2 \xi_j^{(1)} = \beta'(\epsilon_H) \left( \bar{P}_{|j}^i - \frac{\delta^i_j}{3} \Delta_0 \bar{P} \right); \quad (80) $$

$$ a^2 \mathcal{R}^{(1)} = \mathcal{R} + 2 \left( \Pi^{ki}_{k|i} - \frac{2}{3} P^i \right) \quad (81) $$

$$ \begin{aligned} a^2 \tau^i_j^{(1)} = \mathcal{F}^i_j + 2\Pi^i_{k|j}{}^k - \Pi^i_{j|k}{}^k \\ - \frac{1}{3} \left( 2\Pi^k_{|k}{}^l \delta^i_j + P^i_{|j} - \frac{1}{3} \Delta_0 P \delta^i_j \right), \end{aligned} \quad (82) $$

with $\mathcal{F}^i_j := \mathcal{R}^i_j - \frac{1}{3}\mathcal{R}\delta^i_j = \tau^i_j^{(1)}(t_i)$, and with $\Delta_0$ the coordinate Laplacian operator in the Lagrangian coordinates $\{X^i\}$, $\Delta_0 := \delta^{ij}\partial_i\partial_j$.

## 1. Master equation for the trace

Contracting the momentum constraints (78) with a spatial derivative $|_j$ yields the first-order evolution equation for the nontrivial part of the scalar curvature:

$$ \begin{aligned} \partial_t (P^{ki}_{|k|i} - P^k{}_{|k}) &= \partial_t \left( \Pi^{ki}_{k|i} - \frac{2}{3} \bar{P}^k{}_{|k} \right) \\ &= -2H\beta'(\epsilon_H)\Delta_0 \bar{P}. \end{aligned} \quad (83) $$

From the respective expressions (58) and (79) for $\mathcal{R}^{(1)}$ and $\mathcal{A}^{(1)}$, this amounts to the following evolution for $\mathcal{R}^{(1)}$:

$$ \begin{aligned} \partial_t \mathcal{R}^{(1)} + 2H\mathcal{R}^{(1)} &= -4H a^{-2} \beta'(\epsilon_H) \Delta_0 \bar{P} \\ &= -4H \mathcal{A}^{(1)}, \end{aligned} \quad (84) $$

which unlike the case of dust does remain coupled to the dynamics of the inhomogeneous perturbation.

Combining this evolution equation with the linearized energy constraint (72) and its time derivative one then obtains the master equation for the evolution of the trace (63) of the first-order deformation field$^9$:

$$ \begin{aligned} \partial_t^2 \bar{P} + 2H(1 - 3\beta'(\epsilon_H))\partial_t \bar{P} - \mathcal{W}(t) N_H^2 \bar{P} \\ = a^{-2} N_H^2 \beta'(\epsilon_H) \Delta_0 \bar{P}, \end{aligned} \quad (85) $$

where $p_H = \beta(\epsilon_H)$ and $N_H = F(\epsilon_H)/(\epsilon_H + p_H)$ still, and

$$ \begin{aligned} \mathcal{W}(t) &:= 4\pi G[\epsilon_H + p_H - (2\epsilon_H + \tilde{\Lambda})\beta'(\epsilon_H)][4 - 3\mathcal{V}(t)] \\ &= 4\pi G[\epsilon_H + p_H + (\epsilon_H - 3p_H + 2\tilde{\Lambda})\beta'(\epsilon_H)] \\ &\quad - 12\pi G(2\epsilon_H + \tilde{\Lambda})(\epsilon_H + p_H)\beta''(\epsilon_H). \end{aligned} \quad (86) $$

<footnote>$^8$The expression given for $\tau^i_j^{(1)}$ makes use of the momentum constraints (78), which imply, through their spatial derivative, $\partial_t \Pi_{k|i|j}^i = 0$, and thus $\Pi_{k|i|j}^i = \Pi_{k|i|j}^i(t_i) = 0$. Also note that since $P$ and $\bar{P}$ differ by an initial spatial function, we can express (79)–(82) in terms of either variable. Here we have adopted the most compact possibility, noting that the initial value of $\bar{P}$ is nonzero, whereas (81) and (82) involve the initial curvature which is independent of the initial perturbation field.</footnote>

<footnote>$^9$This equation can also be derived by combining the energy constraint (72) with the trace (76) of the evolution equation to eliminate $\mathcal{R}^{(1)}$, or equivalently by directly expanding the Raychaudhuri equation (34) to first order. In both cases, the master equation for the trace would then be recovered after replacing the first-order acceleration divergence $\mathcal{A}^{(1)}$ by its explicit expansion (79).</footnote>
To avoid potential confusion, since the time coordinate $t$ used in this paper has a different rate as compared to the conventional cosmic time, it will sometimes be convenient for further applications to use the (time-coordinate-independent) background scale factor $a$ as the time variable instead. With this change of parametrization, the energy constraint (72) and the master equation for the trace (85) may be rewritten as follows:

$$a \frac{\partial \bar{P}}{\partial a} + \alpha_0 \bar{P} = - \frac{N_H^2}{4H^2} \mathcal{R}^{(1)}; \quad (87)$$

$$\frac{\partial^2 \bar{P}}{\partial a^2} + \frac{\alpha_1}{a} \frac{\partial \bar{P}}{\partial a} - \frac{\alpha_2}{a^2} \bar{P} = \frac{\alpha_3}{a^4} \Delta_0 \bar{P}, \quad (88)$$

respectively, with time-dependent coefficients,

$$\begin{align*}
\alpha_0 &:= 4\pi G \frac{N_H^2}{H^2} [\epsilon_H + p_H - (2\epsilon_H + \bar{\Lambda})\beta'(\epsilon_H)]; \\
\alpha_1 &:= \alpha_0 + 4\pi G \frac{N_H^2}{H^2} [\bar{\Lambda} - 2p_H]; \\
\alpha_2 &:= N_H^2 \mathcal{W}(t)/H^2; \quad \alpha_3 := N_H^2 \beta'(\epsilon_H)/H^2, \quad (89)
\end{align*}$$

where we recall that from the background Friedmann equation we have $H^2/N_H^2 = 4\pi G(2\epsilon_H + \bar{\Lambda})/3$.

From Eq. (88) we can introduce a time-dependent background Jeans wave number $k_J(\epsilon_H)$ by$^{10}$

$$k_J(\epsilon_H) := \frac{1}{c} \sqrt{\frac{\alpha_2}{\alpha_3}} = \frac{1}{c} \sqrt{\frac{\mathcal{W}(t)}{\beta'(\epsilon_H)}}, \quad (90)$$

provided that the term in the square root is positive. Pressure should be positive for sound waves to resist gravitational collapse, and the existence of the Jeans length is intimately related to the energy conditions satisfied by the matter field.

A remark is in order here. In general, one would expect the evolution of the inhomogeneous deformation to be affected by the local, inhomogeneous speed of sound and density, so that a nonperturbative Lagrangian realization would rather feature a local Jeans wave number $k_J(\epsilon)$ [43]. The dynamics in the presence of a significant density contrast will thus only be partially captured by the above first-order equation, where $\epsilon$ has been expanded in $P^a_i$ and, accordingly, only zero-order background factors such as $k_J(\epsilon_H)$ survive in front of the first-order $\bar{P}$.

As in the dust case, the advantages of the Lagrangian approach are only fully realized via nonlinear extrapolation, e.g., by computing the energy density as a full nonlinear functional from the first-order deformation. This is part of the Relativistic Zel'dovich Approximation scheme, as defined for dust fluids in [23]. As in the dust

case and in contrast to standard Eulerian linear perturbation schemes, applying this procedure to compute the energy density out of the solution to first-order equations such as (85) will already capture part of the nonlinear features. This is due to the nonlinear extrapolation and to the use of Lagrangian spatial coordinates which follow the fluid propagation in an exact manner. Further nonlinear effects of inhomogeneous pressure will, however, still be missed due to the absence of local Jeans length contributions in the equation used for $\bar{P}$, compared to what should appear in the nonperturbative evolution equation.

We will not go beyond this procedure in the present paper. Let us nonetheless suggest here a possible direction for improvement. It would require properly defining the local Jeans length in the relativistic context as a functional of the deformation. This quantity would then replace the background Jeans length in the trace master equation. The corresponding nonlinear master equation could then be solved in an iterative manner, by computing at each step the local Jeans length via functional extrapolation out of the previous estimate for the deformation field. Note that each step would also involve a search for the traceless part of the deformation, as all of its components would be required for the extrapolation.

The evolution equation (85) may be rewritten in an alternative form via a time-dependent rescaling of the variable $\bar{P} \mapsto \bar{P}/N_H(t)$. Using the variation rate (67) of the background lapse one finds the more transparent form:

$$\begin{align}
\partial_t^2 \left( \frac{\bar{P}}{N_H} \right) + 2H \partial_t \left( \frac{\bar{P}}{N_H} \right) - 4\pi G (\epsilon_H + p_H) N_H^2 \left( \frac{\bar{P}}{N_H} \right) & \nonumber \\
&= a^{-2} N_H^2 \beta'(\epsilon_H) \Delta_0 \left( \frac{\bar{P}}{N_H} \right). \tag{91}
\end{align}$$

Dust limit.—Setting $p_H = \beta(\epsilon_H) = 0$, we find $\mathcal{W}(t) = 4\pi G\epsilon_H = 4\pi GQ_H = 4\pi GQ_{Hi}a^{-3}$ and $N_H(t) = (\epsilon_H + p_H)^{-1}Q_H = 1$, and consequently both $t$-variable forms of the trace master equation, Eqs. (85) and (91), reduce to the dust deformation trace evolution equation of [23–26]. The trace master equation becomes

$$\partial_t^2 P + 2H \partial_t P - 4\pi G Q_{Hi} a^{-3} P = -4\pi G Q_{Hi} a^{-3} \delta \epsilon_i. \quad (92)$$

With $N_H = 1$ the time variable used is the standard FLRW time coordinate $\tilde{t} = t$, so that the above time derivatives coincide with those used in [23–26] (denoted there by overdots). Finally, as evaluating Eq. (92) at the initial time gives $W = -4\pi G Q_{Hi} \delta \epsilon_i$, its right-hand side can always be rewritten as $Wa^{-3}$, and the dust-case master equation for the trace [e.g., Eq. (41) of [26]] is thus recovered.

Newtonian limit.—The Newtonian limit is obtained by the joint application of the Minkowski Restriction (MR) for the deformation field, as introduced for dust in [23,24], and of

<footnote>$^{10}$We include the factor $c$ explicitly so that the dimensional content of this relation is clear. The right-hand side of (86) must be divided by $c^2$ if units $c \neq 1$ are restored.</footnote>
the $c \to \infty$ limit together with the assumption of a non-relativistic pressure.

The latter two assumptions imply that the pressure is no longer a source of the gravitational field, as the energy density is then $\epsilon \simeq \rho c^2 \gg p$ (where the constant $c$ has been temporarily restored), so that all source terms reduce to the contribution of $\rho$. Note that $\rho$ can be considered as equal to the actual rest mass density in this limit. A further consequence of this is that the lapse becomes trivial, $N = \rho c^2/(\epsilon + p) \simeq 1$, consistent with its spatial variation rate, $N^{-1}N_{|i} = -(\epsilon + p)^{-1}p_{|i} \simeq -(\rho c^2)^{-1}p_{|i} \to 0$ when $c \to \infty$, for any pressure spatial gradient. It is also the case for the (already homogeneous, but generally time-dependent) background lapse that $N_H \simeq 1$. Consequently, the fluid-orthogonal hypersurface time label $t$ now coincides with the fluid's proper time $\tau$ (since $1 \simeq N = \partial_t \tau$) as well as with the standard background cosmic (proper) time $\tilde{t}$. All these notions thus consistently define a time reference that can be used as the Newtonian absolute time. We will denote the corresponding Lagrangian time-derivative operator by an overdot.

With $N=1$ the line element (3) reduces to the one used in [23,24] for irrotational dust, and one can thus directly use the corresponding definition of the MR in this context.<fnref n="11" /> This restriction amounts to assuming that the initial metric is Euclidean and that the spatial coframes are exact in the three-dimensional hypersurfaces, i.e., that there exist spatial coordinates $x^a = f^a(X^i, t)$ such that $G_{ab} = \delta_{ab}$ and

$$ \eta^a_i = a(t)(\delta^a_i + P^a_i) = f^a_{|i}. \quad (93) $$

In any $t$ = const hypersurface, the spatial line element then reads $ds^2 = \delta_{ab}dx^a dx^b$. The coordinates $x^a$ thus define Cartesian-type Eulerian coordinates in which the metric coefficients are manifestly Euclidean at each time, and they can be used to define a Newtonian spatial reference frame. Through its second equality, Eq. (93) also implies that the deformation 1-forms $P^a$ are also exact and accordingly define a deformation vector $\mathbf{P}$, with components $P^a$,

$$ \mathbf{x} = a(t)[\mathbf{X} + \mathbf{P}(\mathbf{X}, t)], \quad P^a_i := P^a_{|i}. \quad (94) $$

With these two assumptions the master equation (91) on the trace $P = \delta_a^i P^a_i$ becomes an equation on the Lagrangian divergence $\nabla_0 \cdot \mathbf{P} := \delta_a^i P^a_{|i}$ of $\mathbf{P}$:

$$ \begin{aligned} \nabla_0 \cdot \ddot{\mathbf{P}} + 2H \nabla_0 \cdot \dot{\mathbf{P}} - 4\pi G_{QH} (\nabla_0 \cdot \mathbf{P} - \delta_{Q_i}) \\ = a^{-2} \frac{d p_H}{d Q_H} \Delta_0 (\nabla_0 \cdot \mathbf{P} - \delta_{Q_i}), \end{aligned} \quad (95) $$

with $Q_H = a^{-3}Q_{Hi}$ still, and $q_i := Q_{Hi}(1 + \delta_{Q_i})$. Note that, although the pressure itself no longer contributes as a source of gravitation, its spatial gradient still produces an acceleration (as obviously expected in a Newtonian framework), which is why it still affects the dynamics of $\nabla_0 \cdot \mathbf{P}$ above through the sound speed squared factor $dp_H/dQ_H$ in front of its Laplacian.

The above Eq. (95) matches<fnref n="12" /> the corresponding equation for the deformation vector obtained in the Newtonian Lagrangian framework, Eq. (24b) in [34] or Eq. (45) in [42] written for the longitudinal part of the deformation vector. By definition, this part obeys the same evolution equation as the Lagrangian divergence of the vector, as can be seen in the unnumbered equations involving that divergence before Eq. (24a) in [34]. Note that in this reference, the Laplacian term features a local sound speed squared (related to the local Jeans length) $dp/dQ$, but it is already noted there that it should actually be replaced by the background value for consistency with the first-order expansion, and it is indeed replaced by the corresponding background expression in [42].

## 2. Master equation for the traceless part

The first-order evolution of the traceless symmetric part $\Pi^i_j$ is given by Eq. (77), with $\xi^i_j$ and $\tau^i_j$ replaced by their expressions (80) and (82), respectively. Eliminating the initial traceless curvature $\mathcal{T}^i_j$ by evaluation of the evolution equation at the time corresponding to the initial condition (114), then first yields the following evolution equation for the traceless symmetric part:

$$ \begin{aligned} \partial_t^2 \Pi^i_j + 3H[1 - \beta'(\epsilon_H)] \partial_t \Pi^i_j \\ + \frac{N_H^2}{a^2} \left( 2\Pi^i_{k|j}{}^k - \Pi^i_{j|k}{}^k - \frac{2}{3} \Pi^k_{l|k}{}^l \delta^i_j \right) \\ = \frac{N_H^2}{3a^2} \left( [1 + 3\beta'(\epsilon_H)] \mathcal{D}^i_j \bar{P} - [1 + 3\beta'(\epsilon_{Hi})] \mathcal{D}^i_j \bar{P}_i \right) \\ + \frac{N_H^2}{a^2 N_{Hi}^2} ({}^tW^i_j + H_i[1 - 3\beta'(\epsilon_{Hi})]^tU^i_j). \end{aligned} \quad (96) $$

<footnote n="11">¹¹Note that the Minkowski Restriction introduced for the dust case is in principle independent of a possible $c \to \infty$ limit and can still otherwise be applied in a Minkowskian regime, as the name suggests. In the present case, when $c$ is still finite, this procedure would need to be extended to the presence of pressure and consequently of an inhomogeneous lapse. We believe, however, that such an extension to this case would require a modification of the perturbation framework used so far in this paper, through the use of a spacetime foliation better adapted to this purpose, and we will consequently not attempt to provide such a generalization here.</footnote>

<footnote n="12">¹²Equation (95) features additional contributions from the initial density perturbations $\delta_{Q_i}$ as compared to the original Newtonian result obtained in [34]. These perturbations were actually neglected there, by assuming $q_i = q_{Hi}$, as is also assumed in Zel'dovich's original work for the dust case [22]. However, as is demonstrated in Appendix A of [34], such an assumption can be made without loss of generality in the Newtonian context within the first-order perturbation scheme in the deformation vector, through a suitable change of Lagrangian coordinates, making both approaches equivalent.</footnote>
Here $\bar{P}_i = -\alpha_{Hi}\delta\epsilon_i$ due to the vanishing of the initial spatial perturbation field, and we have introduced the coordinate traceless spatial Hessian operator $\mathcal{D}_j^i := \delta^{ik}\partial_k\partial_j - (1/3)\delta_j^i\Delta_0$.

This equation still explicitly features the trace, but it can be fully expressed in terms of $\Pi^i_j$ by making use of the momentum constraints (78). This can be achieved by rewriting (78) as

$$ \frac{1}{N_H} \partial_t \Pi^i_{j|i} = \frac{2}{3} \partial_t \left( \frac{\bar{P}_{|j}}{N_H} \right). \qquad (97) $$

A time-integration and spatial differentiation of this equation allows one to express $\mathcal{D}_j^i \bar{P}$ as

$$ \frac{\mathcal{D}_j^i \bar{P}}{N_H} = \frac{\mathcal{D}_j^i \bar{P}_i}{N_{Hi}} + \frac{1}{2} \int_{t_i}^t \partial_t (3\Pi^k_{j|k}{}^{|i} - \Pi^k_{|l|k}{}^{|l} \delta^i_j) dt'. \quad (98) $$

The pair of equations { (96), (98) } together comprise the master equation for the traceless part. When $p_H = 0$, one simply has $N_H(t) = 1$ and $\beta'(\epsilon_H) = 0$ so that this master equation reduces to the corresponding one in the dust case, Eq. (43) in [26].

### 3. Master equations for free and scattered gravitational waves

Following the approach developed in [25,26], we can gain more insight into the evolution of $\Pi^i_j$ by splitting the full master equation for the traceless variable into gravito-electric and gravitomagnetic parts.

To this end, we first define a corresponding split of the initial conditions for the traceless variables:

$$ {}^{\mathrm{t}}U^i_j = {}^{\mathrm{t}}E U^i_j + {}^{\mathrm{t}}H U^i_j; \quad {}^{\mathrm{t}}W^i_j = {}^{\mathrm{t}}E W^i_j + {}^{\mathrm{t}}H W^i_j; \quad (99) $$

$$ {}^{\mathrm{t}}H U^i_{j|i} = 0; \quad {}^{\mathrm{t}}H W^i_{j|i} = 0; \quad (100) $$

$$ 2\Delta_0 {}^{\mathrm{t}}E U^i_j + {}^{\mathrm{t}}E U^k_{|l}{}^{|l} \delta^i_j - 3 {}^{\mathrm{t}}E U^i_{k|j}{}^{|k} = 0; \quad (101) $$

$$ 2\Delta_0 {}^{\mathrm{t}}E W^i_j + {}^{\mathrm{t}}E W^k_{|l}{}^{|l} \delta^i_j - 3 {}^{\mathrm{t}}E W^i_{k|j}{}^{|k} = 0. \quad (102) $$

These conditions can be jointly required because of the following geometric identity (taking its first two time derivatives and evaluating them at the initial time):

$$ (2\Delta_0 \Pi^i_j + \Pi^k_{|l}{}^{|l} \delta^i_j - 3\Pi^i_{k|j}{}^{|k})|_i = 0. \quad (103) $$

This in turn is due to $\Pi^k_{[i|j]k} = 0$, which is a consequence of the momentum constraints (see footnote 8).

We can then define the gravitoelectric and gravitomagnetic traceless parts, respectively, $E^i \Pi^i_j$ and $H^i \Pi^i_j$, from their vanishing initial values and their respective initial first time derivatives ${}^{\mathrm{t}}E U^i_j$ and ${}^{\mathrm{t}}H U^i_j$, requiring them to obey the following evolution equations:

$$ \begin{aligned} \partial_t^2 \Pi^i_j + 3H[1 - \beta'(\epsilon_H)]\partial_t^H \Pi^i_j - \frac{N_H^2}{a^2} \Delta_0^E \Pi^i_j \\ &= \frac{N_H^2}{a^2 N_{Hi}^2} ({}^{\mathrm{t}}H W^i_j + H_i [1 - 3\beta'(\epsilon_{Hi})]{}^{\mathrm{t}}H U^i_j); \end{aligned} \quad (104) $$

$$ \begin{aligned} \partial_t^2 \Pi^i_j + 3H[1 - \beta'(\epsilon_H)]\partial_t^E \Pi^i_j + \frac{N_H^2}{3a^2} \Delta_0^E \Pi^i_j \\ &= \frac{N_H^2}{3a^2} ([1 + 3\beta'(\epsilon_H)]\mathcal{D}_j^i \bar{P} - [1 + 3\beta'(\epsilon_{Hi})]\mathcal{D}_j^i \bar{P}_i) \\ &\quad + \frac{N_H^2}{a^2 N_{Hi}^2} ({}^{\mathrm{t}}E W^i_j + H_i [1 - 3\beta'(\epsilon_{Hi})]{}^{\mathrm{t}}E U^i_j). \end{aligned} \quad (105) $$

Equation (104) is the master equation for free gravitational waves, while Eq. (105), after elimination of the coupling to the trace, is the master equation for the gravitational wave part that is scattered at the fluid source. We will discuss the coupling to the trace of this latter equation in more detail below.

The above evolution equations ensure that we indeed get a decomposition of the traceless deformation field obeying (96) at all times:

$$ \Pi^i_j = E^i \Pi^i_j + H^i \Pi^i_j. \quad (106) $$

They will also propagate the initial constraints (99)–(102) that define the split of ${}^{\mathrm{t}}U^i_j$ and ${}^{\mathrm{t}}W^i_j$. This will ensure the preservation at all times of the divergence-free nature of free gravitational waves as well as the geometric identity on their scattered part, similar to the dust case (cf., [25,26]):

$$ {}^H \Pi^i_{j|i} = 0; \quad (107) $$

$$ 2\Delta_0^E \Pi^i_j + E^i \Pi^k_{|l}{}^{|l} \delta^i_j - 3 E^i \Pi^i_{k|j}{}^{|k} = 0. \quad (108) $$

The (also propagating) momentum constraints (97) split as follows:

$$ {}^H \Pi^i_{j|i} = 0; \quad \frac{1}{N_H} \partial_t E^i \Pi^i_{j|i} = \frac{2}{3} \partial_t \left( \frac{\bar{P}_{|j}}{N_H} \right). \quad (109) $$

Observe that ${}^H \Pi^i_j$ decouples from the trace in both the momentum constraints and the evolution equation, while $E^i \Pi^i_j$ remains coupled to the trace in both cases.

Alternatively, using a time integral of the momentum constraints,

$$ E^i \Pi^i_{j|i} = \frac{2}{3} \int_{t_i}^t N_H \partial_t \left( \frac{\bar{P}_{|j}}{N_H} \right) dt', \quad (110) $$

the geometric constraint (108) on $E^i \Pi^i_j$ can be expressed as follows:

$$ \Delta_0^E \Pi^i_j = \mathcal{D}_j^i \left( \int_{t_i}^t N_H \partial_t \left( \frac{\bar{P}_{|j}}{N_H} \right) dt' \right). \quad (111) $$
This is to be compared to the dust-case relation, Eq. (51) in [26], to which it reduces when $p_H = 0$ and accordingly $N_H(t) = 1$: $\Delta_0 \bar{E} \Pi_i^j = \mathcal{D}_i^j (\bar{P} - \bar{P}_i) = \mathcal{D}_i^j P$. Hence, in the presence of pressure, in contrast to the dust case, the gravitoelectric traceless part and the trace, although still tightly coupled, will in general have different time behaviors.

With the antisymmetric part vanishing at all times, the evolution equations for the trace and for the gravitoelectromagnetic split of the traceless symmetric part, coupled through the momentum constraints, characterize the behavior of the first-order Lagrangian deformation field for this general barotropic single fluid. These evolution equations have yet to be complemented by the set of initial constraints (48)–(53), to which we turn now.

### E. First-order initial conditions

The constraints on the initial conditions for the deformation field, the density and the spatial curvature are expressed at the first-order level as follows:

$$U_{[ij]} = 0; \quad W_{[ij]} = 0; \qquad (112)$$

$$W - 6H_i \beta'(\epsilon_{Hi}) U = -N_{Hi}^2 \alpha_{Hi} [W(t_i) \delta\epsilon_i + \beta'(\epsilon_{Hi}) \Delta_0(\delta\epsilon_i)]; \quad (113)$$

$$\begin{aligned} & {}^{\mathrm{d}}W^i_j + H_i[1 - 3\beta'(\epsilon_{Hi})]{}^{\mathrm{d}}U^i_j \\ & = -N_{Hi}^2 \mathcal{F}^i_j - N_{Hi}^2 \alpha_{Hi} \beta'(\epsilon_{Hi}) \left[ (\delta\epsilon_i)^{i|_j} - \frac{1}{3} \Delta_0(\delta\epsilon_i) \delta^i_j \right]; \end{aligned} \qquad (114)$$

$$H_i U = -\frac{1}{4} \mathcal{R} N_{Hi}^2 + 4\pi G N_{Hi}^2 \alpha_{Hi} \delta\epsilon_i \times [\epsilon_{Hi} + p_{Hi} - (2\epsilon_{Hi} + \bar{\Lambda})\beta'(\epsilon_{Hi})]; \quad (115)$$

$$U^i_{j|i} - U_{lj} = 2H_i \alpha_{Hi} \beta'(\epsilon_{Hi}) (\delta\epsilon_i)_{lj}; \qquad (116)$$

$$p_i = p_{Hi} + \epsilon_{Hi} \beta'(\epsilon_{Hi}) \delta\epsilon_i; \quad p_{Hi} = \beta(\epsilon_{Hi}). \qquad (117)$$

This set of initial conditions can also be obtained by evaluating the linearized Lagrange-Einstein system at the initial time. It can be complemented by the requirements (99)–(102) which define the initial split into gravitoelectric and gravitomagnetic parts of the traceless deformation field.

Note that the above set keeps more variables coupled than the corresponding ones in [26]. This is to be expected, since in the dust case a vanishing pressure and a constant lapse allowed for the elimination of $\epsilon$ and $\Lambda$ between the first two constraints, leaving only a relation among $U$, $W$ and $\mathcal{R}$. Here, we also have contributions from $p$, $\Lambda$ (due to the lapse factor in the $\Lambda$ term) and the nonvanishing $A_i^{(1)}$. Accordingly, the dependence on the initial energy density $\epsilon_i$ and its spatial derivatives can no longer be explicitly removed in general. However, as in the dust case, the scalar

constraints (113) and (115), together with the initial EoS (117), show that only two independent first-order initial conditions need to be given for the scalar variables $U$, $W$, $\mathcal{R}$, $\epsilon_i$, and $p_i$. One could for instance only specify $U$ and $W$ as can be done in the dust case, fully determining the other scalar initial conditions. In contrast to the dust case, however, determining $\epsilon_i$ in this situation would involve solving for the Laplacian differential equation (113).

## IV. APPLICATION TO SPECIFIC EQUATIONS OF STATE

Concrete results can be obtained by looking at special cases of the barotropic EoS. In this section, we will first consider the family of linear relations between the pressure and the energy density. We then proceed to a special nonlinear polytropic EoS that allows one to model the isotropic part of a velocity dispersion field up to late epochs of nonlinear structure formation.

### A. Case of a linear equation of state: $p=w\epsilon$

In the previous section we have derived the evolution equations for the first-order deformation field, sourced by a general barotropic fluid. In this section we will consider as an example the simplest barotropic EoS, $p = \beta(\epsilon) = w\epsilon$ with $w$ a constant parameter obeying the dominant energy condition, $-1 \le w \le 1$. In addition to the radiation fluid, with $w = 1/3$, other interesting cases include a “stiff fluid” corresponding to a free scalar field source, with $w = 1$, and a “curvature” or “string gas” equation of state, with $w = -1/3$. For this class of linear EoS we can readily apply the procedure suggested in [25,26] to find the relativistic Lagrangian first-order solutions.

The formal rest mass density $F(\epsilon)$ and the lapse are found to be as follows:

$$F(\epsilon) = Q_1 \left(\frac{\epsilon}{\epsilon_1}\right)^{1/(1+w)}; \quad N = \frac{Q_1}{\epsilon_1(1+w)} \left(\frac{\epsilon}{\epsilon_1}\right)^{-w/(1+w)}, \qquad (118)$$

if $w \neq -1$. (The case $w = -1$ for a “vacuum energy equation of state” can be treated separately by the explicit cosmological term.)

The solution (42) of the energy conservation law then yields the energy density, and the lapse as deduced from (118), as the following respective functionals of the coframes, with $J = \det(\eta^{a_i})$:

$$\epsilon = \epsilon_i J^{-(1+w)}; \quad N = N_i J^w. \qquad (119)$$

Similar equations hold for the background spacetime,

$$\epsilon_H = \epsilon_{Hi} a^{-3(1+w)}; \quad N_H = N_{Hi} a^{3w}; \quad \frac{\partial_t N_H}{N_H} = 3wH. \qquad (120)$$
Given the linear barotropic EoS, the pressure and background pressure are immediately deduced from the expression of the corresponding energy densities, and will share their functional dependencies.

## 1. First-order equations

With the linear EoS $\beta(\epsilon) = w\epsilon$, $\beta'(\epsilon_H)$ reduces to the constant value $w$, $\beta''(\epsilon_H)$ vanishes at all times, and $\alpha_{Hi} = (1+w)^{-1}$. Consistent with a first-order evaluation of the exact formulas above, the first-order expressions (62)–(66) for $\bar{P}$, $\epsilon$, $p$, $F(\epsilon)$, $N$ (and its rate of evolution) thus simplify to

$$ \bar{P} = P - (1+w)^{-1} \delta \epsilon_i; $$

$$ \epsilon = \epsilon_H [1 - (1+w)\bar{P}]; \quad p = p_H - w(1+w)\epsilon_H \bar{P}; $$

$$ F(\epsilon) = F(\epsilon_H)[1 - \bar{P}]; $$

$$ N = N_H[1 + w\bar{P}]; \quad \frac{\partial_t N}{N} = 3wH + w\partial_t \bar{P}. \quad (121) $$

Equation (73) reduces to

$$ \mathcal{V}(t) = \frac{\epsilon_H (1-w)^2 - 2w \tilde{\Lambda}}{\epsilon_H (1-w) - w \tilde{\Lambda}}, \quad (122) $$

and the first-order Lagrange-Einstein system (72), (75)–(78) becomes<fnref n="13" />

$$ \begin{aligned} \partial_t \mathfrak{P}_{ij} &= 0: \quad \mathfrak{P}_{ij} = \mathfrak{P}_{ij}(t_i) = 0; \\ \partial_t^2 \bar{P} + 3H \frac{\epsilon_H (1-w)^2 + 2w^2 \tilde{\Lambda}}{\epsilon_H (1-w) - w \tilde{\Lambda}} \partial_t \bar{P} &= N_{Hi}^2 a^{6w} \left[ \mathcal{A}^{(1)} - \frac{\epsilon_H (1-w)(1+3w) + 2w \tilde{\Lambda}}{4\epsilon_H (1-w) - 4w \tilde{\Lambda}} \mathcal{R}^{(1)} \right]; \end{aligned} \quad (123) $$

$$ \partial_t^2 \Pi^i_j + 3H(1-w)\partial_t \Pi^i_j = N_{Hi}^2 a^{6w} (\xi_j^{(1)} - \tau^i_j^{(1)}); \quad (124) $$

$$ \begin{aligned} & H \partial_t \bar{P} + 4\pi G [\epsilon_H (1-w) - w \tilde{\Lambda}] N_{Hi}^2 a^{6w} \bar{P} \\ & = -\frac{1}{4} N_{Hi}^2 a^{6w} \mathcal{R}^{(1)}; \end{aligned} \quad (125) $$

$$ \partial_t \left( \Pi^i_{ji} - \frac{2}{3} \bar{P}_{lj} \right) = -2wH\bar{P}. \quad (126) $$

The acceleration gradient and its trace and traceless parts are expressed in terms of the deformation field at first order according to Eqs. (68), (79), and (80), yielding

$$ \partial_t^2 \bar{P} + 3H(1-w)\partial_t \bar{P} = N_{Hi}^2 a^{6w} \left[ \mathcal{A}^{(1)} - \frac{1+3w}{4} \mathcal{R}^{(1)} \right]. $$

$$ \mathcal{A}_j^{(1)} = a^{-2w} \bar{P}^i{}_j; \quad (127) $$

$$ \mathcal{A}^{(1)} = a^{-2w} \Delta_0 \bar{P}; \quad (128) $$

$$ \xi_j^{(1)} = a^{-2w} \mathcal{D}_j \bar{P}, \quad (129) $$

while the first-order expressions (58), (81), and (82) of the Ricci tensor and its trace/traceless split are formally unchanged. Since for the chosen EoS, $\mathcal{W}(t)$ yields

$$ \begin{aligned} \mathcal{W}(t) &= 4\pi G [\epsilon_H (1-w)(1+3w) + 2w\tilde{\Lambda}] \\ &= 4\pi G [\epsilon_{Hi} a^{-3(1+w)} (1-w)(1+3w) + 2w\tilde{\Lambda}], \end{aligned} \quad (130) $$

the master equation (85) for the trace of the perturbation now reads

$$ \begin{aligned} \partial_t^2 \bar{P} + 2H(1-3w)\partial_t \bar{P} & - 4\pi G N_{Hi}^2 [\epsilon_{Hi}(1-w)(1+3w)a^{3(w-1)} + 2w\tilde{\Lambda} a^{6w}] \bar{P} \\ &= w N_{Hi}^2 a^{6w-2} \Delta_0 \bar{P}. \end{aligned} \quad (131) $$

In turn, the master equation (96) for the traceless symmetric part of the deformation field becomes

$$ \begin{aligned} & \partial_t^2 \Pi^i_j + 3H(1-w)\partial_t \Pi^i_j \\ & \qquad + N_{Hi}^2 a^{6w-2} \left\{ 2\Pi^i_{k|j}{}^k - \Pi^i_{j|k}{}^k - \frac{2}{3}\Pi^k_{|l|k}{}^l \delta^i_j \right. \\ & \qquad \left. - \frac{1}{3}(1+3w)\mathcal{D}^i_j (\bar{P} - \bar{P}_i) \right\} \\ & = a^{6w-2} [\mathcal{d}W^i_j + (1-3w)H_i^{\mathcal{d}} U^i_j], \end{aligned} \quad (132) $$

with, from the momentum constraints (126),

$$ a^{-3w} \mathcal{D}^i_j \bar{P} = \mathcal{D}^i_j \bar{P}_i + \int_{t_i}^{t} \frac{\partial_t (3\Pi^k_{j|k}{}^i - \Pi^k_{|l|k}{}^l \delta^i_j)}{2a^{3w}} dt'. \quad (133) $$

We can finally rewrite the set of initial conditions (112)–(117) for the linear EoS:

$$ U_{[ij]} = 0; \quad W_{[ij]} = 0; \quad (134) $$

$$ W - 6wH_i U = - \frac{N_{Hi}^2}{1+w} \left( w\Delta_0(\delta\epsilon_i) + 4\pi G [\epsilon_{Hi}(1-w)(1+3w) + 2w\tilde{\Lambda}]\delta\epsilon_i \right); \quad (135) $$

$$ \mathcal{d}W^i_j + (1-3w)H_i^{\mathcal{d}} U^i_j = -N_{Hi}^2 \left[ \mathcal{F}^i_j + \frac{w}{1+w} \mathcal{D}^i_j (\delta\epsilon_i) \right]; \quad (136) $$

<footnote n="13">¹³It is worth noting in the case when $\Lambda = 0$, $\mathcal{V}(t)$ simplifies further and reduces to the constant $1-w$, so that (123) becomes</footnote>
$$H_i U = -\frac{1}{4} \mathcal{R} N_{Hi}^2 + \frac{4\pi G N_{Hi}^2}{1+w} [\epsilon_{Hi}(1-w) - w\tilde{\Lambda}]\delta\epsilon_i; \qquad (137)$$

$$U_{j|i}^i - U_{|j} = 2 \frac{w}{1+w} H_i (\delta\epsilon_i)_{|j}; \qquad (138)$$

$$p_i = p_{Hi} + w\epsilon_{Hi}\delta\epsilon_i; \qquad p_{Hi} = w\epsilon_{Hi}. \qquad (139)$$

## 2. Solutions for the trace of the deformation field

Similarly to [23,24], we will now further investigate the behavior of the trace P of the first-order deformation. For simplicity, we will restrict attention to the case of a vanishing cosmological constant, $\Lambda = 0$, as may be reasonably assumed during the radiation-dominated era. In this case Eqs. (87)–(89) reduce to

$$a \frac{\partial \bar{P}}{\partial a} + \frac{3}{2} (1-w) \bar{P} = \frac{-3}{32\pi G \epsilon_{Hi}} a^{3(1+w)} \mathcal{R}^{(1)}; \qquad (140)$$

$$\frac{\partial^2 \bar{P}}{\partial a^2} + \frac{\alpha_1}{a} \frac{\partial \bar{P}}{\partial a} - \frac{\alpha_2}{a^2} \bar{P} = \alpha_{3i} a^{3w-1} \Delta_0 \bar{P}, \qquad (141)$$

with the constant parameters

$$\alpha_1 = \frac{3(1-3w)}{2}; \qquad \alpha_2 = \frac{3(1-w)(3w+1)}{2}; \\ \alpha_{3i} = \frac{3w}{8\pi G \epsilon_{Hi}}. \qquad (142)$$

If $w > 0$ (implying $\alpha_{3i} > 0$), as we will assume in the following, then Eq. (141) is a second-order hyperbolic partial differential equation (PDE).<fnref n="14" /> This equation is formally analogous to the standard Eulerian propagation equations for a linearized density contrast [63,72,73] once those are reexpressed in terms of the variable $a$.<fnref n="15" /> In the Eulerian case, assuming global flat-space spatial coordinates, one can find the analytical general solution using a Fourier transformation. A discussion of the differences between the Eulerian and Lagrangian approaches has been given in [26]. (See also the related discussion in [21].) Reference [26] also elucidated a procedure for finding general-relativistic Lagrangian first-order solutions for the deformation field in the dust case. We show here that this procedure can be readily extended to the presence of

pressure and apply it to the determination of a Lagrangian solution for the trace part.<fnref n="16" />

First, we can use the formal identity of Eq. (141), written in Lagrangian coordinates on the nontrivial spacetime manifold, with an equation written in Euclidean space. We can thus work within this flat space with its effective “Eulerian” Cartesian spatial coordinates $x^i$ and solve Eq. (141) with $\Delta_0 \mapsto \delta^{ij}\partial_{x^i}\partial_{x^j}$ for the unknown $\bar{P}(a, x)$. On this space we can then apply an inverse Fourier transformation

$$\bar{P}(a, x) = \iiint \bar{P}_k(a, k) e^{-i\mathbf{k}\cdot\mathbf{x}} d^3\mathbf{k}, \qquad (143)$$

and thus get a second-order linear ordinary differential equation:

$$\frac{d^2 \bar{P}_k}{da^2} + \frac{\alpha_1}{a} \frac{d \bar{P}_k}{da} - (\alpha_2 a^{-2} - \alpha_{3i} k^2 a^{3w-1}) \bar{P}_k = 0, \qquad (144)$$

where we have used $\mathbf{k} \cdot \mathbf{x} := \delta_{ij}k^i x^j$ and $k := (\delta_{ij}k^i k^j)^{1/2}$.

In this case the background Jeans wave number (90) satisfies

$$k_J(\epsilon_H)^2 = \frac{\alpha_2}{\alpha_{3i}} a^{-3(1+w)} \\ = 4\pi G \epsilon_{Hi} \frac{(1-w)(3w+1)}{wa^{3(1+w)}}, \qquad (145)$$

where we recall that $0 < w \le 1$ is assumed. The behavior of the solution to Eq. (141) will then depend on the relative values of $k$ and $ak_J(\epsilon_H)$.

One can first proceed by investigating the extreme cases, as is commonly done in the Eulerian analyses. When $k \ll ak_J(\epsilon_H)$, Eq. (144) may be solved as

$$\bar{P}_k = a^{1+3w} C_{k,1} + a^{3(w-1)} C_{k,2}, \qquad (146)$$

where $C_{k,1(2)}$ are two functions of $\mathbf{k}$ encoding the initial conditions. This corresponds, as expected, to the unstable regime since the term with coefficient $C_{k,1}$ is a growing mode.

In the opposite situation when $k \gg ak_J(\epsilon_H)$, the solution reads

$$\bar{P}_k = a^{\frac{9w-1}{4}} [J_\hat{\nu}(Ba^{\frac{1+3w}{2}}k)C_{k,1} + Y_\hat{\nu}(Ba^{\frac{1+3w}{2}}k)C_{k,2}]; \\ B := \frac{2\sqrt{\alpha_{3i}}}{1+3w}; \qquad \hat{\nu} := \frac{9w-1}{2+6w}, \qquad (147)$$

with different $\mathbf{k}$-dependent coefficients $C_{k,1(2)}$, and where $J_\nu(x)$ and $Y_\nu(x)$ denote the Bessel functions of the first and

<footnote n="14">¹⁴It would be an elliptic PDE for $w < 0$ (i.e., $\alpha_{3i} < 0$), while for the parabolic case $w = 0$ (and consequently $\alpha_{3i} = 0$) it reduces, as expected, to the evolution equation for the dust case, with decoupled time and space variables.</footnote>

<footnote n="15">¹⁵Note that in terms of the conventional cosmic time $\tilde{t}$ introduced in (38), Eq. (141) reduces to $\partial_t^2 \bar{P} + (2-3w)a^{-1} \times \partial_\tau a \partial_\tau \bar{P} - 4\pi G[(1-w)(1+3w)\epsilon_H + 2w\tilde{\Lambda}]\bar{P} = wa^{-2}\Delta_0 \bar{P}$. This is formally equivalent to the linearized Eulerian equation (3.2.17) of Ref. [73] in that the coefficients agree, but both the dependent and (spatial) independent variables differ.</footnote>

<footnote n="16">¹⁶A complementary picture of an equivalent procedure is shown in Appendix A 2 and applied to the search for a particular solution for the traceless part.</footnote>
second kind, respectively. This corresponds to a “stable” regime of acoustic oscillations, although their amplitude will grow over time [as $a^{(3w-1)/2}$ for large $a$] for an unusual EoS with $w > 1/3$. The latter remark includes the “stiff fluid” EoS $w = 1$, for which the above solution is exact at all times, since it corresponds to $k_J(\epsilon_H) = 0$.

From the expression (145) of $k_J(\epsilon_H)$, the noncomoving Jeans wave number $ak_J(\epsilon_H)$ decreases over time, so that even an initially unstable solution will eventually enter the stable regime. Such a solution will cross the threshold $k \approx ak_J(\epsilon_H)$ and it may be useful to be able to describe this transition period as well.

As in the Newtonian case in the Eulerian approach, with different coefficients (see, e.g., [74]), the Bessel functions actually allow for an explicit solution of Eq. (144) for any mode at all times. The general solution is the same as (147) up to a change of the order of the Bessel functions:

$$ \begin{aligned} \bar{P}_k &= a^{\frac{9w-1}{4}} [J_\nu(Ba^{\frac{1+3w}{2}}k)C_{k,1} + Y_\nu(Ba^{\frac{1+3w}{2}}k)C_{k,2}]; \\ B &= \frac{2\sqrt{\alpha_{3i}}}{1+3w}; \quad \nu = \frac{5+3w}{2+6w}. \end{aligned} \quad (148) $$

The integration constants $C_{k,1(2)}$ are derived from the initial conditions on $\bar{P}$ and its time derivative, $\bar{P}_i(\mathbf{x})$ and $U(\mathbf{x})$. To this end, one formally replaces these quantities by functions of the “Eulerian” coordinates $x^i$ on the Euclidean space, with the same functional dependence, $\bar{P}_i(\mathbf{x})$ and $U(\mathbf{x})$. One is then working on flat space, and the respective Fourier transforms $\bar{P}_k(a = a_i = 1, \mathbf{k})$ and $(\partial_t \bar{P}_k)(a = 1, \mathbf{k}) = H_i(\partial_a \bar{P}_k)(a = 1, \mathbf{k})$ can be computed, from which $C_{k,1(2)}(\mathbf{k})$ are deduced. Knowing these, $\bar{P}(a, \mathbf{k})$ is expressed as the full solution given by Eq. (148) and its inverse Fourier transform (143) gives $\bar{P}(a, \mathbf{x})$ in Euclidean space.

Finally, one can formally replace the Eulerian spatial coordinates by the Lagrangian ones in $\bar{P}(a, \mathbf{x})$ while preserving the functional form. The resulting Lagrangian function $\bar{P}(a, \mathbf{X})$ then gives a solution to the evolution equation (141) in the nonconstant curvature spatial sections, thanks to the algebraic identity of this equation with its Euclidean space counterpart. It is now a Lagrangian solution, however, and must be interpreted as such: the coordinates $X^i$ are comoving with the inhomogeneous fluid flow. They are local coordinates on the perturbed manifold; thus the solution $P(a, \mathbf{X})$ describes perturbations as they evolve in the perturbed space. This perturbed space is in general not isometric to Euclidean space. Note that the Fourier modes $\bar{P}(a, \mathbf{k})$ are only an intermediate resolution step as they only correspond to modes in the ancillary Euclidean space. As the inversion of the solution (148) does not allow for an explicit general analytic expression, it requires the specification of the initial conditions and will usually involve numerical integration with the given $C_{k,1(2)}(\mathbf{k})$ to realize this solution procedure.

## B. Case of a polytropic equation of state: $p = \kappa Q^\gamma$

As a second class of models we will now turn to the nonlinear case of polytropic equations of state.

### 1. Equation of state and resolution procedure

We consider the polytropic EoS, $p = \kappa Q^\gamma$, $Q = F(\epsilon)$, where $\kappa$ is the polytropic constant, and $\gamma > 1$ the polytropic exponent. For such flows the pressure and the energy density obey the relation [63,75]

$$ \epsilon = \beta^{-1}(p) = \frac{1}{\gamma - 1} p + A p^{1/\gamma} = \frac{1}{\gamma - 1} \kappa Q^{\gamma} + A \kappa^{1/\gamma} Q, \quad (149) $$

where $A$ is a constant parameter. We will assume in this section that the formal $Q = F(\epsilon)$ actually coincides with the rest mass density of the fluid, e.g., via suitable initial conditions. For $A = 0$, we again obtain the (nondust) linear case $p = w\epsilon$ with $w := \gamma - 1 > 0$. In the following, we will instead consider the case $A\kappa^{1/\gamma} = 1$ (in particular $A > 0$), corresponding to an EoS of the type of a nonrelativistic adiabatic ideal gas, the energy density being the sum of the rest mass density and an internal energy density equal to $p/(\gamma - 1)$.

As a relevant example, we will focus on the case $\gamma = 5/3$, which is an exact solution for a locally isotropic distribution with velocity dispersion, derived from the relativistic kinetic theory of collisionless matter [76]. (See also [77] and references therein.) This EoS also coincides with the corresponding exact solution in Newtonian cosmology derived from kinetic theory [41,42]. In these latter papers it is also shown that this particular EoS arises in the inhomogeneous case by closing the hierarchy of kinetic equations through truncation of the third and higher reduced moments. In the inhomogeneous case this law is, however, phenomenological, since there is a nonvanishing anisotropic part. Neglecting this part strictly results in shear-free motion confirming the exactness of the law in the homogeneous case.

The conservation law (18), combined with $p = \kappa Q^\gamma$, gives for the evolution of $p$

$$ \partial_t p + \gamma N \Theta p = 0; \quad \gamma = \frac{5}{3}. \quad (150) $$

The same relation holds within the background spacetime, so that $p_H a^5 = p_{H_i} a_i^5$. The assumption of the background sources following the same EoS also gives, for $\gamma = 5/3$,

$$ \begin{aligned} \epsilon_H &= \beta^{-1}(p_H) = \frac{3}{2} p_H + A p_H^{3/5}; \\ \beta'(\epsilon_H) &= \frac{2}{3} \frac{5}{5 + 2A p_H^{-2/5}}; \quad \beta''(\epsilon_H) = \frac{80 A p_H^{-7/5}}{9(5 + 2A p_H^{-2/5})^3}. \end{aligned} \quad (151) $$
The procedure outlined in the last subsection for solving the trace master equation, Eq. (88), in terms of Fourier transformation within a set of coordinates formally equivalent to Eulerian spatial coordinates, is still applicable in this case. We can thus substitute (151) and (89) in the Eulerian coordinate analogue of (88), and solve the corresponding ordinary differential evolution equation for each Fourier mode. This has to be performed by numerical integration as the more complicated time evolution of the coefficients prevents an explicit analytic solution. Once initial conditions are specified we can then numerically compute the inverse Fourier transform and formally replace the (Eulerian) spatial coordinates by the Lagrangian coordinates $X^i$ (see Sec. IVA 2) to obtain the solution for $\bar{P}(t, X^i)$.

## 2. Behavior of the first-order trace for a model overdense region

As an instructive toy model, we will now consider the evolution of an initial spherical Gaussian deformation:

$$ -\bar{P}_i = \alpha_{Hi} \delta \epsilon_i = c_i \exp \left( -\frac{R^2}{2\sigma^2} \right), \quad (152) $$

where $\sigma$ and $c_i$ respectively define the characteristic scale and maximum amplitude of the initial perturbation, and $R := (\delta_{ij}X^i X^j)^{1/2}$ is a Lagrangian coordinate “radius.”$^{17}$ We will take $c_i > 0$ and $c_i \ll 1$. The perturbation can then be seen to describe a small initial local overdensity, since the initial rest mass density contrast,

$$ \delta_i := \frac{\rho_i}{\rho_{Hi}} - 1 = \frac{F(\epsilon_i)}{F(\epsilon_{Hi})} - 1 = \frac{F(\epsilon_{Hi}[1 + \delta\epsilon_i]) - F(\epsilon_{Hi})}{F(\epsilon_{Hi})}, \quad (153) $$

is well approximated by $\alpha_{Hi}\delta\epsilon_i = -\bar{P}_i$ for $c_i \ll 1$.

The actual value of the amplitude $c_i$ is irrelevant for the evolution of $\bar{P}$ itself, since it obeys a linear equation. However, it will matter for the nonlinear evaluation of any physical quantity such as $\rho$ determined by the first-order solution for $\bar{P}$ through the extrapolation procedure mentioned above from the Relativistic Zel'dovich Approximation. To best illustrate the effect of this procedure, we choose a rather large overdensity with the arbitrary amplitude $c_i = 10^{-3}$ at an initial time that corresponds to the epoch of last scattering. As we will see, this will let the unstable perturbations enter the mildly

nonlinear regime (where $|\bar{P}| < 1$ but is of order 1) around the present epoch, i.e., around $a = a_0 \approx 1090$ since we set $a_i = 1$.

The other independent initial condition amounts to specifying the first time derivative $(\partial_t \bar{P})(t_i)$. For this we simply consider an initially stationary deformation and set $(\partial_t \bar{P})(t_i) = U = 0$.

The present formalism focuses on the description of a single fluid source, as it allows for a description in terms of a single velocity field and a single EoS. We will consequently make the simplifying assumption of a model universe filled with a single-component matter fluid and a cosmological constant. The description of model universes with multicomponent fluids is beyond the scope of the present paper and is left to future work. The background density parameters $\Omega_m, \Omega_\Lambda$ for the matter component and the cosmological constant, respectively, satisfy $\Omega_m + \Omega_\Lambda = 1$. We will take the present epoch value $\Omega_\Lambda^0 = 0.692$ in agreement with the best-fit $\Lambda$CDM parameters from the Planck Collaboration [78].

The background is also affected by the polytropic EoS (149) of the source fluid. As noted above, our polytrope is exact for the background and is parametrized by the arbitrary constant $\kappa$, or equivalently $A$ as we set $A\kappa^{1/\gamma} = 1$. Specifying its value amounts to choosing the initial instability scale as determined by $k_J(\epsilon_{Hi})$. It also controls the ratio between pressure and rest mass density at a given time, and hence the deviation of the background from a dust-fluid $\Lambda$CDM model. The value we adopt for our examples below, $A\rho_{Hi}^{-2/5} = 3/2$, requires the background fluid pressure to be relativistic (and radiationlike) at the initial time, $\rho_{Hi} = \epsilon_{Hi}/3$, with $\rho_{Hi}/\rho_{Hi} = 2/3$. However, it subsequently quickly becomes negligible as $\rho_H/\rho_H \propto a^{-2}$, keeping the late-time dynamics of the background very close to that of the $\Lambda$CDM model.

We choose to make the lengths $R$, $\sigma$ dimensionless by setting the initial instability scale $k_J(\epsilon_{Hi})^{-1}$ [as derived from substituting (151) into (90) at the initial time] to be our length unit. Thus $\sigma < 1$ means that the scale of the initial perturbation is below the Jeans scale $k_J(\epsilon_{Hi})^{-1}$ and above it for $\sigma > 1$. For the value of $A$ adopted in the present example and estimating $\rho_{Hi}$ from $\Lambda$CDM background parameters [78], this length unit is approximately 98 kpc. This would correspond to a large background comoving initial overdensity size of $a_0 k_J(\epsilon_{Hi})^{-1} \approx 107$ Mpc.$^{18}$

<footnote>$^{17}$We have chosen the set of Lagrangian coordinates $X^i$ such that the components of the spatial metric at initial time, $G_{ij}$, are approximately $\delta_{ij}$ (at leading order) in these coordinates. They can thus be considered as Cartesian-like coordinates, and $R$ is thus a fluid-comoving radial coordinate. It does not, however, coincide with the spatial metric distance between the fluid elements of the respective Lagrangian coordinates ($X^i$) and (0,0,0). [This is true irrespective of a possible normalization by $a(t)$ to make it a background comoving distance.]</footnote>

<footnote>$^{18}$Note that $k_J(\epsilon_{Hi})^{-1}$ defines an initial instability “scale” only in terms of Lagrangian coordinates, e.g., in terms of $R$. This means that the corresponding “background comoving” distance, $a(t)k_J(\epsilon_{Hi})^{-1}$ evaluated at present time, does not coincide with the present-day physical size of an object that would initially have been of this scale, as such a size must be evaluated using the actual, deformed, spatial metric. (See previous footnote.) $a_0 k_J(\epsilon_{Hi})^{-1}$ may be seen as a rough estimate of this physical size, as obtained by fully neglecting the deformations $G_{ab}^{(1)}$, $P^a_i$, in the evaluation of the integrated spatial line element.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 1. Numerical solution for the first-order trace $-\bar{P}$ in Lagrangian space, for an initial spherical Gaussian overdensity with a peak amplitude of $10^{-3}$ at $R=0$ and a standard deviation $\sigma$ such that $k_J(\epsilon_{H_i})\sigma = 10$. (a) Evolution of $-\bar{P}$ as a function of $a$ for fixed values of the Lagrangian radius $R$. From top to bottom, $R=0, 10, 20$ and $30$. (b) Spatial variation of $-\bar{P}$ with $R$, for several values of the background scale factor. From bottom to top, $a=1, 10, 200, 500$ and $1000$. The perturbation strongly grows over time, corresponding to a collapsing structure.</figcaption></figure>

Figures 1–3 show the numerical results for $\bar{P}$ with the procedure, initial conditions and parameters given above, for three different values of $\sigma$.

The first case, $\sigma = 10$ (Fig. 1), corresponds to a super-Jeans length, and hence unstable, initial perturbation. Figures 1(a) and 1(b) show the numerical results for the evolution of the perturbation $-\bar{P}$ as a function of the scale factor at several values of $R$, and over the whole range of radii $R$ for increasing values of $a$, respectively. As expected, this perturbation is unstable and remains so by growing at all times, the pressure gradient being insufficient

to prevent the collapse of the structure. The evolution is similar to the dust case with the fast onset of a linear growth of the perturbation with $a$ before a late-time slowdown due to the presence of $\Lambda$.

The second case, $\sigma = 0.2$ (Fig. 2), illustrates the opposite situation of an initially sub-Jeans length perturbation. Figures 2(a) and 2(b) show the numerical solution for $-\bar{P}$ in this situation along the same reasoning as for Figs. 1(a) and 1(b). At the early stage, the pressure gradient dominates and opposes the gravitational collapse. The perturbation behaves as an acoustic wave and is damped

<figure><img src="image_6.png" /><figcaption>FIG. 2. Numerical solution for the first-order trace $-\bar{P}$ in Lagrangian space, for an initial spherical Gaussian overdensity with a peak amplitude of $10^{-3}$ at $R=0$ and a standard deviation $\sigma$ such that $k_J(\epsilon_{H_i})\sigma = 0.2$. (a) Evolution of $-\bar{P}$ as a function of $a$ at fixed distance $R$. From top to bottom at $a = 1000, R=3, R=4, R=1$ and $R=0$. The inset panel shows a detail of the early evolution (small values of $a$), where only the $R=0$ (solid line) and $R=1$ (dashed line) are visibly nonzero. (b) Spatial variation of $-\bar{P}$ with the Lagrangian radius, for several values of the background scale factor. The structure is first damped and spread out by the Lagrangian pressure gradient, before starting to grow back after the critical wave number $ak_J(\epsilon_H)$ has increased, as the perturbation enters the unstable regime.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 3. Numerical solution for the first-order trace $-\bar{P}$ in Lagrangian space, for an initial spherical Gaussian overdensity with a peak amplitude of $10^{-3}$ at $R=0$ and a standard deviation $\sigma$ such that $k_J(\epsilon_{H1})\sigma = 1$. (a), (b) Evolution of $-\bar{P}$ as a function of $a$ at a given distance $R$, for late and early times, respectively. From top to bottom at $a=1000$ for (a), $R=3$, $R=1$, $R=0$, $R=4$, $R=5$; same order for (b) at $a=20$. (c) Spatial variation of $-\bar{P}$ with $R$, for fixed values of the background scale factor. From top to bottom at $R=0$, $a=1000$, $a=500$, $a=200$, $a=1$, $a=20$, $a=4$. The behavior is rather similar to the previous case of $k_J(\epsilon_{H1})\sigma = 0.2$; as expected, the unstable regime is, however, reached sooner, and the perturbation then grows similarly to the case of $k_J(\epsilon_{H1})\sigma = 10$, up to much above its initial amplitude.</figcaption></figure>

as it propagates away from the initial peak at $R=0$. However, the instability wave number $ak_J(\epsilon_H)$ quickly starts increasing over time (cf., Fig. 4). That is why around $a=50$ to $100$ the perturbation starts to grow as its typical wave number (estimated by $\sigma^{-1}=5$) ends up below the critical value, with $ak_J(\epsilon_H) = 5$ for $a \approx 94$, and it enters the unstable regime. The peak of this growing structure remains at a mostly stationary Lagrangian position, at $R \approx 3.7$, while its increasing amplitude still remains small and below the initial value $-\bar{P}(a=1, R=0) = 10^{-3}$ up to present time ($a \approx 1090$).

For comparison we also consider the special case where the initial scale lies at the stability threshold, $\sigma = 1$. The evolution of the corresponding solution for $-\bar{P}$ with $a$ at several radii is shown in Figs. 3(a) and 3(b), with the latter highlighting the early evolution ($1 \le a \le 20$). Figure 3(c) shows the spatial dependence of $-\bar{P}$ with $R$ at some values

of the scale factor. The behavior of the perturbation in this case is as expected intermediate, with an initial acoustic damping and propagation away from $R=0$ similarly to the $\sigma = 0.2$ case, but more rapidly entering an unstable regime, after $a \approx 5$. The amplitude of the perturbation then starts growing with a dustlike behavior up to beyond 20 times its initial value at present time, with a shifted peak as in the $\sigma = 0.2$ case, that stays around $R \approx 2.5$.

### 3. Evaluating the nonlinear density contrast

As we recalled above, even the first-order Lagrangian perturbation scheme allows one to probe part of the nonlinear regime in the evaluation of observable quantities. This involves extrapolating these observables as exact, nonlinear functionals of the deformation field, the latter being evaluated as a solution to its first-order evolution equations and constraints.
<figure><img src="image_8.png" /><figcaption>FIG. 4. Evolution of the instability wave number $ak_J(\epsilon_H)$ with the scale factor $a$ for the polytropic EoS considered here, with the unit of length convention $k_J(\epsilon_H) = 1$. As this wave number only depends on the background by construction, this result applies to all examples considered in Sec. IV B. After a small initial dip, $ak_J(\epsilon_H)$ exceeds its initial value around $a \approx 4$ and enters the increasing power law regime $ak_J(\epsilon_H) \propto \sqrt{a}$ [valid as long as $(\Omega_N/\Omega_m)(a/a_i)^{-2} \ll 1$, which is satisfied up to the present epoch] as expected from the large $a$ expansion of its expression for the present polytropic EoS.</figcaption></figure>

Adopting this procedure for the rest mass density we evaluate it as the exact integral to the rest mass conservation equation (18):

$$Q = \frac{Q_i}{J}; \quad J = \det(\eta^a_i) = a^3 \det(\delta^a_i + P^a_i), \quad (154)$$

where $P^a_i$ are the components of the deformation field. The density contrast $\delta$ is then deduced from the above:

<figure><img src="image_9.png" /></figure>

$$\delta := \frac{Q - Q_H}{Q_H} = \frac{Q_i}{Q_{Hi} a^{-3} J} - 1; \quad a^{-3} J = \det(\delta^a_i + P^a_i), \quad (155)$$

and it is evaluated by replacing $P^a_i$ by the first-order solution.

Using the polytropic EoS and the parameters adopted here, the lapse may then be computed from

$$N = \frac{Q}{\epsilon + p} = \frac{Q}{Q + \frac{\gamma}{\gamma-1} \kappa Q^{\gamma}} = \frac{1}{1 + \frac{5}{3}(1 + \delta)^{2/3} a^{-2}}, \quad (156)$$

with $\delta$ expressed from the deformation field as above. This formula shows that the lapse is 1 in pressure-free (here empty) regions ($\delta = -1$) and decreases with increasing density contrast at a given time. The deviation $(1 - N)$ rapidly decreases over time as $\propto a^{-2}$, with late time values of order $10^{-6}$ (when $a \approx 1000$), as long as $\delta$ remains at most of order unity.

We will now illustrate this process for the density contrast with two examples using the same polytropic EoS as above. Note that this evaluation requires the knowledge of all components of the deformation field, including the traceless part. We specify procedures in Appendix to obtain a particular (gravitoelectric) solution for the first-order traceless part from the initial conditions for the trace in specific cases. These procedures have been used to determine a consistent solution for the full deformation field in the examples below. We have also made use of the fact that the initial density $Q_i = F(\epsilon_{Hi}[1 + \delta\epsilon_i])$ is well approximated by $F(\epsilon_{Hi})(1 + \alpha_{Hi}\delta\epsilon_i) = Q_{Hi}(1 - \bar{P}_i)$ for a small, still linear, initial density perturbation (with $\alpha_{Hi} = 3/4$ for the chosen EoS parameters) for the evaluation of $\delta$.

<figure><img src="image_10.png" /><figcaption>FIG. 5. Numerical evaluation of the nonlinear density contrast $\delta$ as extrapolated from the first-order Lagrangian perturbation, where the initial $-\bar{P}$ is the same spherical Gaussian field as for Fig. 1, with peak value of $10^{-3}$ and $k_J(\epsilon_{Hi})\sigma = 10$. (a) Evolution of $\delta$ with the background scale factor at fixed distances $R$. From top to bottom, $R = 0, 10, 20$ and $30$. (b) Spatial variation of $\delta$ with the Lagrangian radius, for given values of $a$. From bottom to top, $a = 1, 10, 200, 500$ and $1000$. The overall behavior of $\delta$ is similar to the results of Fig. 1 for the first-order $-\bar{P}$ in the same situation, but the extrapolated density contrast grows faster at late times near the $R = 0$ maximal overdensity. Additional nonlinear effects concerning the comparison with a standard perturbation approach, not studied here, could also be revealed by using instead as the $x$-axis for (b) the actual spatial metric distance to the $R = 0$ fluid element (as an “Eulerian radius”), altering the spatial dependence. (See the discussion in Sec. IV B 4.)</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 6. Comparison of the extrapolated nonlinear density contrast $\delta$ (dashed line) with the first-order solution for the sign-inverted deformation trace $-\bar{P}$ (solid line) within the same setting as Figs. 1 and 5. (a) Comparison of the evolution of both quantities as a function of $a$ at the center of the overdensity ($R=0$). (b) Comparison of the spatial variation of both quantities with $R$ at a late time ($a=1000$). In this situation, the perturbation grows large enough to enter the nonlinear regime and to render the time evolution and spatial behavior of the extrapolated $\delta$ clearly deviating from those of $-\bar{P}$.</figcaption></figure>

**Localized overdensity.**—Let us first retain the “spherical” initial overdensity example studied thus far in this section, with the initial conditions for the trace given by (152), with $c_i = 10^{-3}$, and $U=0$. The first-order solution for the trace in this situation has been determined above, and it is complemented by a gravitoelectric solution for the first-order traceless part through the use of the procedure given in Appendix A 2 that directly applies to this case. The determinant $J$ is then computed from this solution as in Appendix A 4, giving $\delta$ from Eq. (155).

Note that when all components of the deformation field are very small, i.e., when it lies fully in the linear regime, then the extrapolated $\delta$ remains quantitatively close to $-\bar{P}$, which corresponds to its expansion at first order in the deformation field. This is the case in the initially stable or marginally stable cases $\sigma = 0.2$ and $\sigma = 1$, where the initial acoustic damping of the perturbation keeps its amplitude small up to the present time despite the late-time growth. In both of these cases, the resulting density contrast indeed remains indistinguishable from the value of $-\bar{P}$ already depicted above (Figs. 2–3).

We will consequently focus, from now on, on the case $\sigma = 10$, where the unstable deformation reaches into the mildly nonlinear regime before the present time, as can be seen for the trace (whose amplitude reaches about 0.5 at the present epoch).

Figures 5(a) and 5(b) show the result of the nonlinear evaluation of the density contrast in this situation, as a function of $a$ at given radii $R$, and as a function of the radius at several moments in its evolution, respectively. Although the general behavior is roughly similar to that of $-\bar{P}$ (cf., Fig. 1), nonlinear effects are visible in the amplified growth of $\delta$ at late times near $R=0$, with a maximal overdensity reaching about 0.7 at present.

This nonlinear deviation of the density contrast functional with respect to its first-order estimate $-\bar{P}$ is made explicit by the direct comparison of the peak ($R=0$) amplitude evolution of $\delta$ and $-\bar{P}$ as a function of the background scale factor in Fig. 6(a). The spatial dependence on $R$ of both quantities at late times, compared in Fig. 6(b) at $a=1000$, is also visibly affected by the amplified growth of the density contrast where $\bar{P}$ is no longer small, i.e., around $R=0$.

**Lagrangian monochromatic wave.**—The second toy model we consider is that of a single Lagrangian monochromatic wave deformation. The choices of background parameters and the length unit [$k_J(\epsilon_{Hi})=1$] are unchanged. The initial perturbation is now chosen to be

$$-\bar{P}_i = c_i \cos(KX); \quad U=0, \quad (157)$$

where we will again take $c_i = 10^{-3}$ as an initial amplitude. This situation corresponds to an initially stationary monochromatic wave in the given Lagrangian coordinate set,<fnref n="19" /> $-\bar{P}_i = c_i \cos(\delta_{ij}K^i X^j + \phi_0)$ with $\phi_0 = 0$, and a

<footnote n="19">¹⁹Similarly to the interpretation of $R$ for the previous example, it is important to keep in mind that the perturbation we are considering here only has a sinusoidal dependence in the chosen Lagrangian coordinates $X^i$. It would have a different functional dependence in terms of actual physical (metric) spatial distance between two points on a given hypersurface $t=\text{const}$. One expects for instance, at a given late time $t$ and along a given spatial geodesic line, the distance between the successive perturbation nodes at $KX = -\pi/2$ and $KX = \pi/2$ (surrounding a collapsing overdensity) to be shorter than the distance between the nodes at $KX = \pi/2$ and $KX = 3\pi/2$ (surrounding an expanding underdensity), despite all nodes being equally separated in terms of the Lagrangian coordinate $X$.</footnote>
<figure><img src="image_12.png" /><figcaption>FIG. 7. Numerical evaluation of the nonlinear density contrast $\delta$ as extrapolated from the first-order Lagrangian perturbation. The first-order deformation trace is taken as a plane wave in Lagrangian coordinates of wave vector $\mathbf{K}$ (of norm $K$) along the $X$ coordinate, $-\bar{P} \propto \cos(KX)$, of initial amplitude $10^{-3}$. The result is shown at a given time as a function of $KX$ for three possible values of $K$, which is expressed in units $k_J(\epsilon_{Hi}) = 1$. (a) At $a = 10$, for $K = 0.1$ ($K^{-1} = 10$), $K = 5$ ($K^{-1} = 0.2$) and $K = 1$ by order of decreasing amplitude. (b) At $a = 1000$, for $K = 0.1$ ($K^{-1} = 10$), $K = 1$ and $K = 5$ ($K^{-1} = 0.2$) by order of decreasing amplitude. The side panel displays the (otherwise barely visible) latter two curves on a different vertical scale. The most unstable perturbation, for $K^{-1} = 10$, displays a nonsinusoidal asymmetric shape at late times as it reaches the mildly nonlinear regime. This shape would be further nonlinearly modified, via a different $x$-axis dependence, if this axis were expressed alternatively in terms of an Eulerian-type, regularly spaced (in terms of spatial metric distances) $x$ coordinate.</figcaption></figure>

Lagrangian wave vector **K** along the first coordinate $X$, with components $K^i = (K, 0, 0)$.

The first-order trace solution then remains in this monochromatic mode form in the Lagrangian coordinates at all times, $\bar{P} = \hat{P}_{\mathbf{K}}(t) \cos(KX)$. The amplitude $\hat{P}_{\mathbf{K}}(t)$ evolves according to the ordinary differential equation (A2) which is solved by numerical integration for a given wave number $K$. A gravitoelectric solution for the traceless part is then determined along the lines of Appendix A 1, where the relevant amplitude $\hat{Q}_{\mathbf{K}}(t)$ is again numerically evaluated, knowing $\hat{P}_{\mathbf{K}}(t)$, through its defining time integral formula (A3). From these, one can calculate the density contrast in the same way as in the previous example, with the determinant $J$ evaluated as detailed in Appendix A 4.

Here we again study three cases distinguished by their wave number in direct analogy to the previous example, with $K^{-1}$ playing the role of the characteristic length $\sigma$. We accordingly choose $K^{-1} = 0.2$, $K^{-1} = 1$ and $K^{-1} = 10$, which at the initial time are stable, marginally stable and unstable, respectively. The corresponding spatial dependence of $\delta$ as a function of $KX$ for the three wave number choices is shown at an early time ($a = 10$) in Fig. 7(a) and at a late time ($a = 1000$) in Fig. 7(b).

In this situation, in the first two cases the components of the deformation field again remain small at all times, due to initial acoustic oscillations, and the density contrast thus follows the sinusoidal shape of $-\bar{P}$ at all times. This is also the case for the unstable mode $K^{-1} = 10$ at $a = 10$ when it is still in the linear regime. At $a = 1000$, however, this mode clearly deviates from this behavior as its amplitude is

no longer linear. In particular, an asymmetry develops between the under- and overdensity magnitudes as the latter is sharply amplified by the nonlinear evolution of $\delta$.

## 4. Discussion

In both examples above, the Lagrangian scheme and the proposed extrapolation procedure exhibit nonlinear effects on the overdensity for unstable perturbations when they become large enough. The amplitude of large overdensities in these examples is clearly underestimated when they are approximated by the first-order expression $-\bar{P}$ instead of using the nonlinear extrapolation for $\delta$.

An even higher initial overdensity amplitude could actually lead to a vanishing determinant $a^{-3}J$ at the maximum overdensity at a late enough time, implying $q \to \infty$ with deformation coefficients still of order 1. This situation corresponds to a shell crossing, beyond which the first-order Lagrangian scheme in no longer valid.

The presence of pressure can delay its occurrence by damping the perturbation. An improvement of the perturbative scheme to account for further local nonlinear effects in the dynamical evolution, e.g., allowing for a nonlinear coefficient to define the Jeans length, is needed, however, to fully circumvent this problem. Velocity dispersion effects may in principle allow us to model the multistream regime, and the stabilization of structure formation in the form of virialization, which may help to avoid shell crossings [42,79].

We emphasize that the current Lagrangian perturbation scheme already contains another effect of nonlinear structure evolution, which lies in the exact propagation of
the spatial coordinates used along the fluid flow lines. This is analogous to the inclusion of quadratic convection terms within linear Lagrangian time derivatives in the Newtonian framework.<fnref n="20" />

Let us suggest a procedure that would be required to make these effects explicit also in the relativistic context; its concrete application is beyond the scope of this paper.

Eulerian-like coordinates could first be recovered, at least along a given spatial geodesic direction, by labeling points at equal intervals of spatial metric distances. This would involve solving for the initial metric components $G_{ab}$ such that their Ricci tensor is consistent with the initial conditions (114)–(115) for given initial deformation field data, and then functionally evaluating and integrating the line element as given by (A18) from the first-order solution for $P^a_i$. The resulting length, as a function of a Lagrangian coordinate, could then be used as an estimate of the Eulerian coordinate distance. Finally, this relation would have to be numerically inverted so that a given Lagrangian function obtained through the Relativistic Zel'dovich Approximation, such as $\varphi(X^i)$, could be expressed as a function of the Eulerian coordinate $x$ along the chosen line.

A different functional dependence on this spatial distance [which may be normalized by $a(t)$ to become a background comoving distance], as compared to the fluid-comoving coordinates $X^i$, would thus include nonlinear effects of the fluid-propagation-dependent coordinate transformation.

Recall, however, that a three-dimensional family of Eulerian observers generally does not exist in a relativistic (intrinsic) description. Strictly, a coordinate transformation to Eulerian space can only be conducted after the Minkowski Restriction of the relativistic solution has been executed.

## V. CONCLUSION

In this paper we have generalized the Lagrangian perturbation approach to the nonlinear evolution of inhomogeneous general relativistic model universes containing a single irrotational fluid obeying a general barotropic relation.

By choosing a suitable set of coframes, we obtained the master partial differential equations for the evolution of the trace and traceless parts of the first-order deformation field that reduce to the corresponding equations in the dust case. The trace part also matches the Newtonian limit of the corresponding Lagrangian perturbation problem.

We discussed the procedure proposed in previous papers of how to find the solution for perturbations that propagate in the perturbed space, and we applied this procedure to specific toy models, illustrating the mildly nonlinear evolution of the density contrast. We also discussed the limits of a first-order Lagrangian scheme, and we proposed ideas for a nonperturbative generalization, which is needed especially in application to cases where the pressure term is taken to model multistreaming beyond the mildly nonlinear regime.

## ACKNOWLEDGMENTS

This work is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (ERC Advanced Grant No. 740021-ARTHUS, PI: T. B.). It is also supported by Catalyst Grant No. CSG-UOC1603 administered by the Royal Society of New Zealand. Y. L. is supported by the China Scholarship Council. P. M. is supported by a *spécifique Normalien* Ph.D. grant, and acknowledges support and hospitality for a visit to the University of Canterbury. P. M. also acknowledges support by the National Science Centre, Poland, under Grant No. 2014/13/B/ST9/00845. Y. L. and D. L. W. acknowledge support and hospitality for visits to CRAL-ENS, Lyon. We wish to thank Asta Heinesen for helpful discussions.

## APPENDIX: EXAMPLES OF SOLUTIONS FOR THE GRAVITOELECTRIC TRACELESS PART

In this paper we will not attempt to find the general solution of Eqs. (96) and (97) for the traceless part. We will, however, discuss a procedure for finding one possible solution for suitably chosen traceless-part initial conditions. For any barotropic EoS, this yields one example of a full gravitoelectric solution for all components of the deformation field $P^a_i$. It can then be substituted into exact nonlinear formulas to extrapolate functionals of the coframes such as metric distances or the rest mass density.

To find such an example solution, we will focus on the gravitoelectric part which is directly coupled to the trace, and accordingly we set the gravitomagnetic part to zero.

### 1. Case of a Lagrangian monochromatic wave

Let us first assume that the first-order trace solution can be written as a single monochromatic wave mode in the given set of Lagrangian spatial coordinates $X^i$:

$$ \bar{P}(t, X^i) = \varphi(\mathbf{K} \cdot \mathbf{X}) \hat{P}_{\mathbf{K}}(t), \quad (A1) $$

with constant Lagrangian wave vector $\mathbf{K}$, where $\mathbf{K} \cdot \mathbf{X} := \delta_{ij}K^i X^j$, and $\varphi(\mathbf{K} \cdot \mathbf{X}) = \cos(\mathbf{K} \cdot \mathbf{X} + \phi_0)$, for some constant phase $\phi_0$. This form is a solution of the first-order trace master equation, if and only if $\hat{P}_{\mathbf{K}}(t)$ is a solution of the ordinary differential equation

<footnote n="20">²⁰In addition to the time derivatives being taken at different fixed spatial coordinates, a difference also comes from the spatial derivative operators, such as the Laplacian $\Delta_0$ appearing in the trace master equation (85), being expressed in terms of Lagrangian coordinates and thus differing from the corresponding Eulerian operators. (See [34] for the explicit transformation in the Newtonian case.)</footnote>
$$ \frac{d^2}{dt^2} \hat{P}_{\mathbf{K}} + 2H(1 - 3\beta'(\epsilon_H)) \frac{d}{dt} \hat{P}_{\mathbf{K}} - W(t)N_H^2 \hat{P}_{\mathbf{K}} \\ = -a^{-2}N_H^2 \beta'(\epsilon_H) K^2 \hat{P}_{\mathbf{K}}, \qquad (A2) $$

with $K := (\delta_{ij}K^i K^j)^{1/2}$. Then $P = \bar{P} - \bar{P}_i = \varphi(\mathbf{K} \cdot \mathbf{X}) \times (\hat{P}_{\mathbf{K}}(t) - \hat{P}_{\mathbf{K}}(t_i))$. Setting

$$ \begin{aligned} \hat{Q}_{\mathbf{K}}(t) &:= \int_{t_i}^{t} N_H(t') \partial_t \left( \frac{\hat{P}_{\mathbf{K}}}{N_H} \right) (t') dt' \\ &= \hat{P}_{\mathbf{K}}(t) - \hat{P}_{\mathbf{K}}(t_i) - 3 \int_{t_i}^{t} H(t') \beta'(\epsilon_H)(t') \hat{P}_{\mathbf{K}}(t') dt', \end{aligned} \qquad (A3) $$

the time integral of the momentum constraints (97) is

$$ \Pi^{i}_{jl} = \frac{2}{3} \hat{Q}_{\mathbf{K}}(t) K_j \varphi'(\mathbf{K} \cdot \mathbf{X}). \qquad (A4) $$

We now take $\Pi^i_j$ to be a purely longitudinal mode and get the following solution to the momentum constraints (with $K_j := \delta_{jl}K^l$):

$$ \Pi^i_j = \left( \frac{K^i K_j}{K^2} - \frac{1}{3} \delta^i_j \right) \hat{Q}_{\mathbf{K}}(t) \varphi(\mathbf{K} \cdot \mathbf{X}) \qquad (A5) $$

$$ = \left( \frac{K^i K_j}{K^2} - \frac{1}{3} \delta^i_j \right) \left( \frac{\hat{Q}_{\mathbf{K}}(t)}{\hat{P}_{\mathbf{K}}(t) - \hat{P}_{\mathbf{K}}(t_i)} \right) P(t, X^i). \qquad (A6) $$

Substituting this form into the master equation (96) shows that it is consistently a solution of both equations for the traceless part. It is straightforward to show from the above formula that $2\Delta_0\Pi^i_j + \Pi^k_{\quad j|k} |^l \delta^i_j - 3\Pi^i_{\quad k|j} |^k = 0$; i.e., this $\Pi^i_j$ obeys the defining relation (108) for the gravitoelectric part and evolves according to (105). This solution is thus a pure gravitoelectric one, amounting to setting the grav- itomagnetic part to zero by the choice of vanishing gravitomagnetic traceless part of the initial deformation: $\Pi^i_j = E \Pi^i_j$.

Choosing this solution amounts to specifying the follow- ing (gravitoelectric) initial conditions:

$$ {}^t U^i_j = \left( \frac{K^i K_j}{K^2} - \frac{1}{3} \delta^i_j \right) (U + 3H_i \beta'(\epsilon_{Hi}) \alpha_{Hi} \delta \epsilon_i); \qquad (A7) $$

$$ {}^t W^i_j = \left( \frac{K^i K_j}{K^2} - \frac{1}{3} \delta^i_j \right) \left( W + 3H_i \beta'(\epsilon_{Hi}) U + 3[\partial_t(H\beta'(\epsilon_H))(t_i) + 2H_i^2 \beta'(\epsilon_{Hi})]\alpha_{Hi} \delta \epsilon_i \right). \qquad (A8) $$

This is compatible with the set of constraints on the initial conditions given in Sec. III E, in particular the initial momentum constraints (116) and Eq. (114), provided that the latter is used to specify the traceless part of the initial first-order Ricci tensor $\mathcal{T}^i_j$.

The corresponding full perturbation field $P^i_j = \Pi^i_j + \frac{1}{3} \delta^i_j P$ then reads

$$ P^i_j = \frac{K^i K_j}{K^2} \left( \frac{\hat{Q}_{\mathbf{K}}(t)}{\hat{P}_{\mathbf{K}}(t) - \hat{P}_{\mathbf{K}}(t_i)} \right) P + \frac{1}{3} \delta^i_j \left( 1 - \frac{\hat{Q}_{\mathbf{K}}(t)}{\hat{P}_{\mathbf{K}}(t) - \hat{P}_{\mathbf{K}}(t_i)} \right) P. \qquad (A9) $$

Note that the corresponding deformation 1-forms $\mathbf{P}^a = \delta^a_k P^k_i dX^i$ are not exact due to the different time evolution of the trace and gravitoelectric traceless parts. This con- trasts with the dust case where a purely gravitoelectric perturbation would lead to integrable coframes [26], so that only the nonflat initial metric would prevent one obtaining a Euclidean spatial metric at all times in that situation.

By linearity of the equations, a solution for $\Pi^i_j$ can also be obtained when the trace is a finite sum of such mono- chromatic waves, or the sum of the two time-evolution modes' solutions of the evolution equation (A2) for a given wave vector $\mathbf{K}$, simply by summing the corresponding solutions as given by (A5).

## 2. Case of a spatially localized solution

We assume here either that the spatial slices are globally diffeomorphic to the Euclidean space $\mathbb{R}^3$, i.e., that they can be covered by a single chart, or that the deformation field can be assumed to vanish outside a given chart. In either case it suffices to work within the Euclidean space spanned by the spatial coordinates in a given chart.

Let us now consider a spatially localized solution for the trace, e.g., a local overdensity evolving from an initial Gaussian perturbation in terms of the given set of spatial Lagrangian coordinates, as studied in the numerical exam- ples of Sec. IV. More specifically, we require the solution for the trace to always be a square-integrable function of the spatial coordinates in the chart, so that its Fourier transform in these coordinates can be performed and inverted. We can thus write

$$ \bar{P}(t, X^i) = \iiint e^{-i\mathbf{K} \cdot \mathbf{X}} \hat{P}(t, \mathbf{K}) d^3\mathbf{K}, \qquad (A10) $$

where $\hat{P}(t, \mathbf{K})$ is a solution of the evolution equation (A2) at fixed $\mathbf{K}$, with the initial conditions set by the forward Fourier transform in the chart coordinates:

$$ \hat{P}(t_i, \mathbf{K}) = - \frac{1}{(2\pi)^3} \alpha_{Hi} \iiint e^{i\mathbf{K} \cdot \mathbf{X}} \delta\epsilon_i(\mathbf{X}) d^3\mathbf{X}; \qquad (A11) $$

$$ (\partial_t \hat{P})(t_i, \mathbf{K}) = \frac{1}{(2\pi)^3} \iiint e^{i\mathbf{K} \cdot \mathbf{X}} U(\mathbf{X}) d^3\mathbf{X}. \qquad (A12) $$

Note that the above approach represents an alternative and complementary formulation of the method of solution
presented in [26] which formally replaces the Lagrangian coordinates by Eulerian ones. In the present paper it is applied in Secs. IVA 2 and IV B. The reformulation suggested here allows us to be more explicit about the required assumptions, as well as expressing the coordinate components of tensors such as $\Pi^i_j$ in a more convenient form. In both formulations, the use of plane-wave modes and flat-space Fourier transformations is sufficient since the Lagrangian first-order master equations to be solved only involve the metric-independent coordinate spatial derivatives $|_i$ and Laplacian $\Delta_0 = |_i|_j \delta^{ij}$ as spatial derivative operators.

By linearity of the equations, a solution for the (gravitoelectric) traceless part is obtained by summation of the plane-wave solutions for all Fourier modes:

$$ \Pi^i_j = E^i \Pi^i_j = \iiint e^{-i\mathbf{K} \cdot \mathbf{X}} \frac{K^i K_j}{K^2} \hat{Q}(t, \mathbf{K}) d^3\mathbf{K} \\ -\frac{1}{3} \delta^i_j \iiint e^{-i\mathbf{K} \cdot \mathbf{X}} \hat{Q}(t, \mathbf{K}) d^3\mathbf{K}, \quad (A13) $$

with

$$ \hat{Q}(t, \mathbf{K}) := \int_{t_i}^{t} N_H(t') \partial_t \left( \frac{\hat{P}(t, \mathbf{K})}{N_H(t)} \right) (t') dt'. \quad (A14) $$

Using this solution again implies a specific choice of initial conditions for the traceless deformation field (in particular taking it to be gravitoelectric) and for the traceless part of the spatial Ricci tensor.

In the case of spherically symmetric initial conditions in the chart coordinates, i.e., when $\delta_{ij}(X^i)$ and $U(X^i)$ only depend on $R := (\delta_{ij}X^i X^j)^{1/2}$, their Fourier transform will also depend only on $K$. From the evolution equation (A2), this feature is preserved over time, so that one can write $\hat{P}(t, \mathbf{K})$ as $\hat{P}(t, K)$ and consequently $\hat{Q}(t, \mathbf{K})$ as $\hat{Q}(t, K)$ and $\hat{P}(t, X^i)$ as $\hat{P}(t, R)$. The above solution for $\Pi^i_j$ can then be computed as

$$ \Pi^i_j = \left( \frac{X^i X_j}{R^2} - \frac{1}{3} \delta^i_j \right) q(t, R), \quad (A15) $$

with $X_j := \delta_{jk}X^k$ and

$$ q(t, R) := \frac{4\pi}{R} \int_{0}^{\infty} K \sin(RK) \hat{Q}(t, K) dK \\ - \frac{4\pi}{R^3} \int_{0}^{\infty} \left( \frac{\sin(RK)}{K} - R \cos(RK) \right) \hat{Q}(t, K) dK. \quad (A16) $$

### 3. Time integral of the gravitoelectric evolution equation

The above procedure gives a way of obtaining a traceless part consistent with the momentum constraints and

evolution equations in particular situations, and when only initial conditions on the trace part (or on the energy density) are explicitly specified. Alternatively, and still focusing on a purely gravitoelectric traceless part, a solution can be derived from the gravitoelectric traceless evolution equation (105), if the trace part and the (gravitoelectric) traceless initial conditions are known. It can be achieved by rewriting this evolution equation as follows:

$$ \begin{aligned} & \partial_t \left( \frac{a^3}{N_H} \partial_t E^i \Pi^i_j \right) \\ &= -\frac{aN_H}{3} \mathcal{D}^i_j \left( \int_{t_i}^t N_H \partial_t \left( \frac{\bar{P}}{N_H} \right) dt' \right) \\ &\quad + \frac{aN_H}{3} \left( [1 + 3\beta'(\epsilon_H)] \mathcal{D}^i_j \bar{P} - [1 + 3\beta'(\epsilon_{Hi})] \mathcal{D}^i_j \bar{P}_i \right) \\ &\quad + \frac{aN_H}{N_H^2} \left( \mathrm{t}^{\mathrm{l},\mathrm{E}} W^i_j + H_i [1 - 3\beta'(\epsilon_{Hi})]^{\mathrm{t},\mathrm{E}} U^i_j \right), \quad (A17) \end{aligned} $$

after replacing $\Delta_0 E^i \Pi^i_j$ by its integral expression (111) in terms of $\bar{P}$. It can be readily time integrated twice to give $E^i \Pi^i_j$. This yields the full $\Pi^i_j$ if the initial conditions are chosen such that the gravitomagnetic part vanishes.

In contrast to the previous subsections, this procedure can be applied in general, allowing the gravitoelectric initial conditions for the traceless part to be freely set. However, this requires the initial conditions $\mathrm{t}^i U^i_j = \mathrm{t}^{\mathrm{l},\mathrm{E}} U^i_j$ and $\mathrm{t}^i W^i_j = \mathrm{t}^{\mathrm{l},\mathrm{E}} W^i_j$ to be explicitly specified. While the trace parts relate to the energy density and spatial scalar curvature, the trace-free parts are related to properties of the gravitational wave components at the initial time. The latter have to be set in such a way as to fulfill the momentum constraints and their time derivative at the initial time, as well as the geometric constraints (101)–(102) for the gravitoelectric parts.

### 4. On the evaluation of physical quantities

From given solutions for the trace and traceless parts, the full deformation field is straightforwardly obtained as $P^i_j = \Pi^i_j + (1/3)P \delta^i_j$, with $P = \bar{P} - \bar{P}_i$. This expression can then be inserted into the Lagrangian functional expressions for various physical quantities in terms of the deformation field. They can then be directly evaluated without any further linearization. This extrapolation is a crucial part of the Relativistic Zel'dovich Approximation as defined in [23], and it generally requires the knowledge of all components of the deformation field.

One would for instance directly compute a spatial distance from the line element

$$ ds^2 = a(t)^2 G_{ab} (\delta^a_i + P^a_i) (\delta^b_j + P^b_j) dX^i dX^j, \quad (A18) $$

where knowledge of $G_{ab}(X^k)$ is also required. In turn, the rest mass density [with initial conditions set in such a way that it does coincide with $q = F(\epsilon)$] would be computed as
$$Q = \frac{Q_i}{J} = \frac{Q_{Hi}(1 + \alpha_{Hi}\delta_{ci})}{a^3 \det(\delta_i^a + P_i^a)}. \quad (A19)$$

For the evaluation of the latter, note that in the case of a monochromatic wave (with one or both time-evolution modes), the deformation field components can be written as follows:

$$P_j^i = \lambda_1 \frac{K^i K_j}{K^2} + \lambda_2 \delta_j^i, \quad (A20)$$

and similarly in the case of a localized spherically symmetric perturbation,

$$P_j^i = \lambda_1 \frac{X^i X_j}{K^2} + \lambda_2 \delta_j^i. \quad (A21)$$

The coefficients $\lambda_1(t, X^k)$, $\lambda_2(t, X^k)$ for the monochromatic case are directly deduced from (A9) or from a sum of two such solutions, while in the localized spherically symmetric case, $\lambda_1(t, X^k) = q(t, R)$ and $\lambda_2(t, X^k) = (P(t, R) - q(t, R))/3$. The determinant of the spatial coframe coefficients,

from which $Q$ is evaluated, is then expressed in both cases by

$$J = a^3(1 + \lambda_2)^2(1 + \lambda_1 + \lambda_2), \quad (A22)$$

leading to an infinite rest mass density (from shell crossing) whenever $\lambda_2 \to -1$ or $\lambda_1 + \lambda_2 \to -1$.

Such an extrapolation procedure provides the exact metrical distances, density and other physical properties as produced by the deformation field at a given order. In particular, this gives powerful approximations for the Ricci and Weyl curvatures that are not available in standard perturbation theory. It is, however, clear that the resulting expressions are approximations that must be controlled.

We can further combine the exact functionals for a given deformation with exact averages of Einstein's equations. An example was given in [24] that also showed that the resulting prescription can even lead to exact results. For example, the combination of the first-order Lagrangian dust model with exact averages led to an exact formula for the kinematical backreaction within a class of averaged Lemaître-Tolman-Bondi solutions [24].

[1] R. L. Arnowitt, S. Deser, and C. W. Misner, The dynamics of general relativity, in *Gravitation: An Introduction to Current Research*, edited by L. Witten (Wiley, New York, 1962), pp. 227–265, republished in *Gen. Relativ. Gravit.* 40, 1997 (2008).

[2] É. Gourgoulhon, 3+1 Formalism in General Relativity: Bases of Numerical Relativity, Lecture Notes in Physics (Springer-Verlag, Berlin Heidelberg, 2012).

[3] H. Kodama and M. Sasaki, Cosmological perturbation theory, Prog. Theor. Phys. Suppl. 78, 1 (1984).

[4] F. Bernardeau, S. Colombi, E. Gaztanaga, and R. Scoccimarro, Large scale structure of the universe and cosmological perturbation theory, Phys. Rep. 367, 1 (2002).

[5] H. Asada, Lagrangian description of fluid flow with pressure in relativistic cosmology, Phys. Rev. D 62, 127301 (2000).

[6] H. Asada and M. Kasai, Lagrangian description of the fluid flow with vorticity in the relativistic cosmology, Phys. Rev. D 59, 123515 (1999).

[7] S. Bildhauer, T. Buchert, and M. Kasai, Solutions in Newtonian cosmology—The pancake theory with cosmological constant, Astron. Astrophys. 263, 23 (1992).

[8] F. R. Bouchet, S. Colombi, E. Hivon, and R. Juszkiewicz, Perturbative Lagrangian approach to gravitational instability, Astron. Astrophys. 296, 575 (1995).

[9] T. Buchert, A class of solutions in Newtonian cosmology and the pancake theory, Astron. Astrophys. 223, 9 (1989).

[10] T. Buchert, Lagrangian theory of gravitational instability of Friedmann-Lemaître cosmologies and the “Zel’dovich approximation”, Mon. Not. R. Astron. Soc. 254, 729 (1992).

[11] T. Buchert, Lagrangian perturbation theory—A key-model for large-scale structure, Astron. Astrophys. 267, L51 (1993).

[12] T. Buchert, Lagrangian theory of gravitational instability of Friedmann-Lemaître cosmologies: Generic third-order model for nonlinear clustering, Mon. Not. R. Astron. Soc. 267, 811 (1994).

[13] T. Buchert, Lagrangian perturbation approach to the formation of large-scale structure, in *Proceedings of the International School of Physics Enrico Fermi. Course CXXXII: Dark matter in the Universe*, edited by S. Bonometto, J. Primack, and A. Provenzale (Società Italiana di Fisica, Varenna sul Lago di Como 1996).

[14] T. Buchert and J. Ehlers, Lagrangian theory of gravitational instability of Friedmann-Lemaître cosmologies—Second order approach: An improved model for nonlinear clustering, Mon. Not. R. Astron. Soc. 264, 375 (1993).

[15] T. Buchert and G. Götz, A class of solutions for self-gravitating dust in Newtonian gravity, J. Math. Phys. 28, 2714 (1987).

[16] J. Ehlers and T. Buchert, Newtonian cosmology in Lagrangian formulation: Foundations and perturbation theory, Gen. Relativ. Gravit. 29, 733 (1997).

[17] S. Nadkarni-Ghosh and D. F. Chernoff, Extending the domain of validity of the Lagrangian approximation, Mon. Not. R. Astron. Soc. 410, 1454 (2011).

[18] S. Nadkarni-Ghosh and D. F. Chernoff, Modelling nonlinear evolution using Lagrangian perturbation theory expansions, Mon. Not. R. Astron. Soc. 431, 799 (2013).

[19] C. Rampf, The recursion relation in Lagrangian perturbation theory, J. Cosmol. Astropart. Phys. 12 (2012) 004.
[20] C. Rampf and T. Buchert, Lagrangian perturbations and the matter bispectrum I: Fourth-order model for nonlinear clustering, J. Cosmol. Astropart. Phys. 06 (2012) 021.

[21] E. Villa, S. Matarrese, and D. Maino, Cosmological dynamics: From the Eulerian to the Lagrangian frame. Part I. Newtonian approximation, J. Cosmol. Astropart. Phys. 06 (2014) 041.

[22] Ya. B. Zel'dovich, Gravitational instability: An approximate theory for large density perturbations, Astron. Astrophys. 5, 84 (1970).

[23] T. Buchert and M. Ostermann, Lagrangian theory of structure formation in relativistic cosmology. I. Lagrangian framework and definition of a nonperturbative approximation, Phys. Rev. D 86, 023520 (2012).

[24] T. Buchert, C. Nayet, and A. Wiegand, Lagrangian theory of structure formation in relativistic cosmology. II. Average properties of a generic evolution model, Phys. Rev. D 87, 123503 (2013).

[25] A. Alles, T. Buchert, F. Al Roumi, and A. Wiegand, Lagrangian theory of structure formation in relativistic cosmology. III. Gravitoelectric perturbation and solution schemes at any order, Phys. Rev. D 92, 023512 (2015).

[26] F. Al Roumi, T. Buchert, and A. Wiegand, Lagrangian theory of structure formation in relativistic cosmology. IV. Lagrangian approach to gravitational waves, Phys. Rev. D 96, 123538 (2017).

[27] G. F. R. Ellis and C. G. Tsagas, Relativistic approach to nonlinear peculiar velocities and the Zel'dovich approximation, Phys. Rev. D 66, 124015 (2002).

[28] M. Kasai, Tetrad-based perturbative approach to inhomogeneous universes: A general relativistic version of the Zel'dovich approximation, Phys. Rev. D 52, 5605 (1995).

[29] S. Matarrese and D. Terranova, Post-Newtonian cosmological dynamics in Lagrangian coordinates, Mon. Not. R. Astron. Soc. 283, 400 (1996).

[30] C. Rampf and G. Rigopoulos, Zel'dovich approximation and general relativity, Mon. Not. R. Astron. Soc. 430, L54 (2013).

[31] C. Rampf and A. Wiegand, Relativistic Lagrangian displacement field and tensor perturbations, Phys. Rev. D 90, 123503 (2014).

[32] H. Russ, M. Morita, M. Kasai, and G. Börner, Zel'dovich-type approximation for an inhomogeneous universe in general relativity: Second-order solutions, Phys. Rev. D 53, 6881 (1996).

[33] D. S. Salopek, J. M. Stewart, and K. M. Croudace, The Zel'dovich approximation and the relativistic Hamilton-Jacobi equation, Mon. Not. R. Astron. Soc. 271, 1005 (1994).

[34] S. Adler and T. Buchert, Lagrangian theory of structure formation in pressure supported cosmological fluids, Astron. Astrophys. 343, 317 (1999).

[35] M. Morita and T. Tatekawa, Extending Lagrangian perturbation theory to a fluid with velocity dispersion, Mon. Not. R. Astron. Soc. 328, 815 (2001).

[36] T. Tatekawa, M. Suda, K. Maeda, M. Morita, and H. Anzai, Perturbation theory in Lagrangian hydrodynamics for a cosmological fluid with velocity dispersion, Phys. Rev. D 66, 064014 (2002).

[37] T. Tatekawa, Third-order perturbative solutions in the Lagrangian perturbation theory with pressure, Phys. Rev. D 71, 044024 (2005).

[38] R. Maartens, Causal thermodynamics in relativity, arXiv: astro-ph/9609119.

[39] R. Maartens and J. Triginer, Density perturbations with relativistic thermodynamics, Phys. Rev. D 56, 4640 (1997).

[40] R. Maartens, J. Triginer, and D. R. Matravers, Stress effects in structure formation, Phys. Rev. D 60, 103503 (1999).

[41] T. Buchert and A. Domínguez, Modeling multistream flow in collisionless matter: Approximations for large scale structure beyond shell crossing, Astron. Astrophys. 335, 395 (1998).

[42] T. Buchert and A. Domínguez, Adhesive gravitational clustering, Astron. Astrophys. 438, 443 (2005).

[43] T. Buchert, The non-perturbative regime of cosmic structure formation, Astron. Astrophys. 454, 415 (2006).

[44] D. L. Wiltshire, Cosmic structure, averaging and dark energy, in *Cosmology and Gravitation: Proc. XVth Brazilian School*, edited by M. Novello and S. E. Perez Bergliaffa (Cambridge Scientific Publishers, Cambridge, 2014), pp. 203–244.

[45] T. Buchert, On average properties of inhomogeneous fluids in general relativity. I. Dust cosmologies, Gen. Relativ. Gravit. 32, 105 (2000).

[46] T. Buchert, On average properties of inhomogeneous fluids in general relativity: Perfect fluid cosmologies, Gen. Relativ. Gravit. 33, 1381 (2001).

[47] T. Buchert, Dark energy from structure: A status report, Gen. Relativ. Gravit. 40, 467 (2008).

[48] T. Buchert, P. Mourier, and X. Roy, Cosmological backreaction and its dependence on spacetime foliation, arXiv:1805.10455.

[49] T. Buchert, J. Larena, and J.-M. Alimi, Correspondence between kinematical backreaction and scalar field cosmologies: The “morphon field”, Classical Quantum Gravity 23, 6379 (2006).

[50] J. Larena, J.-M. Alimi, T. Buchert, M. Kunz, and P.-S. Corasaniti, Testing backreaction effects with observations, Phys. Rev. D 79, 083011 (2009).

[51] D. L. Wiltshire, Cosmic clocks, cosmic variance and cosmic averages, New J. Phys. 9, 377 (2007).

[52] D. L. Wiltshire, Exact Solution to the Averaging Problem in Cosmology, Phys. Rev. Lett. 99, 251101 (2007).

[53] D. L. Wiltshire, Average observational quantities in the timescape cosmology, Phys. Rev. D 80, 123512 (2009).

[54] A. Wiegand and T. Buchert, Multiscale cosmology and structure-emerging dark energy: A plausibility analysis, Phys. Rev. D 82, 023523 (2010).

[55] X. Roy and T. Buchert, Chaplygin gas and effective description of inhomogeneous universe models in general relativity, Classical Quantum Gravity 27, 175013 (2010).

[56] X. Roy, T. Buchert, S. Carloni, and N. Obadia, Global gravitational instability of FLRW backgrounds—Interpreting the dark sectors, Classical Quantum Gravity 28, 165004 (2011).

[57] T. Buchert and N. Obadia, Effective inhomogeneous inflation: Curvature inhomogeneities of the Einstein vacuum, Classical Quantum Gravity 28, 162002 (2011).

[58] W. Zimdahl, Cosmic bulk viscosity through backreaction, Gen. Relativ. Gravit. 48, 51 (2016).
[59] J. A. G. Duley, M. A. Nazer, and D. L. Wiltshire, Timescape cosmology with radiation fluid, *Classical Quantum Gravity* **30**, 175006 (2013).

[60] M. A. Nazer and D. L. Wiltshire, Cosmic microwave background anisotropies in the timescape cosmology, *Phys. Rev. D* **91**, 063519 (2015).

[61] W. Valkenburg and F. Villaescusa-Navarro, Accurate initial conditions in mixed dark matter–baryon simulations, *Mon. Not. R. Astron. Soc.* **467**, 4401 (2017).

[62] T. Buchert and M. Carfora, On the curvature of the present-day Universe, *Classical Quantum Gravity* **25**, 195001 (2008).

[63] G. F. R. Ellis, R. Maartens, and M. A. H. MacCallum, *Relativistic Cosmology* (Cambridge University Press, Cambridge, 2012).

[64] J. Ehlers, Akad. Wiss. Lit. (Mainz); Abh. Math.-Nat. Kl. **11**, 793 (1961) [English Translation, Contributions to the relativistic mechanics of continuous media, *Gen. Relativ. Gravit.* **25**, 1225 (1993)].

[65] W. Israel, Nonstationary irreversible thermodynamics: A causal relativistic theory, *Ann. Phys. (N.Y.)* **100**, 310 (1976).

[66] H. Friedrich, Evolution equations for gravitating ideal fluid bodies in general relativity, *Phys. Rev. D* **57**, 2317 (1998).

[67] W. Israel and J. M. Stewart, Transient relativistic thermodynamics and kinetic theory, *Ann. Phys. (N.Y.)* **118**, 341 (1979).

[68] H. Stephani, D. Kramer, M. MacCallum, C. Hoenselaers, and E. Herlt, *Exact Solutions of Einstein's Field Equations*, Cambridge Monographs on Mathematical Physics (Cambridge University Press, Cambridge, 2003).

[69] G. F. R. Ellis, Dynamics of pressure free matter in general relativity, *J. Math. Phys.* **8**, 1171 (1967).

[70] G. F. R. Ellis and H. van Elst, Cosmological models: Cargèse lectures 1998, *Theoretical and Observational Cosmology* NATO Science Series Vol. 541 (Springer, Dordrecht, 1999).

[71] A. F. Zakharov, V. A. Zinchuk, and V. N. Pervushin, Tetrad formalism and reference frames in general relativity, *Phys. Part. Nuclei* **37**, 104 (2006).

[72] P. Peter and J.-P. Uzan, *Primordial Cosmology* (Oxford University Press, New York, 2009).

[73] C. G. Tsagas, A. Challinor, and R. Maartens, Relativistic cosmology and large-scale structure, *Phys. Rep.* **465**, 61 (2008).

[74] R. M. Gailis and N. E. Frankel, Two-component cosmological fluids with gravitational instabilities, *J. Math. Phys.* **47**, 062505 (2006).

[75] L. Rezzolla and O. Zanotti, *Relativistic Hydrodynamics* (Oxford University Press, New York, 2013).

[76] J. Ehlers and W. Rienstra, The locally isotropic solution of the Liouville and Poisson equations, *Astrophys. J.* **155**, 105 (1969).

[77] R. Treciokas and G. F. R. Ellis, Isotropic solutions of the Einstein-Boltzmann equations, *Commun. Math. Phys.* **23**, 1 (1971).

[78] P. A. R. Ade et al. (Planck Collaboration), Planck 2015 results XIII. Cosmological parameters, *Astron. Astrophys.* **594**, A13 (2016).

[79] T. Buchert, Stabilization of large-scale structure by adhesive gravitational clustering, in *From Stars to the Universe*, 1993 Shanghai (PR China) 1998, Ann. Shanghai Obs., Acad. 1994 Sin., No. 21 (2000), pp. 85–92, arXiv:astro-ph/9901002.