# Lagrangian theory of structure formation in relativistic cosmology. IV. Lagrangian approach to gravitational waves

Fosca Al Roumi, Thomas Buchert, Alexander Wiegand

► To cite this version:

Fosca Al Roumi, Thomas Buchert, Alexander Wiegand. Lagrangian theory of structure formation in relativistic cosmology. IV. Lagrangian approach to gravitational waves. Physical Review D, 2017, 96 (12), pp.123538. 10.1103/PhysRevD.96.123538. hal-01704711

HAL Id: hal-01704711

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

Submitted on 4 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Lagrangian theory of structure formation in relativistic cosmology.
## IV. Lagrangian approach to gravitational waves

Fosca Al Roumi,<fnref n="1" /> Thomas Buchert,<fnref n="1" />,* and Alexander Wiegand<fnref n="2" />,<fnref n="3" />,<fnref n="4" />

The relativistic generalization of the Newtonian Lagrangian perturbation theory is investigated. In previous works, the perturbation and solution schemes that are generated by the spatially projected gravitoelectric part of the Weyl tensor were given to any order of the perturbations, together with extensions and applications for accessing the nonperturbative regime. We here discuss more in detail the general first-order scheme within the Cartan formalism including and concentrating on the gravitational wave propagation in matter. We provide master equations for all parts of Lagrangian-linearized perturbations propagating in the perturbed spacetime, and we outline the solution procedure that allows one to find general solutions. Particular emphasis is given to global properties of the Lagrangian perturbation fields by employing results of Hodge-de Rham theory. We here discuss how the Hodge decomposition relates to the standard scalar-vector-tensor decomposition. Finally, we demonstrate that we obtain the known linear perturbation solutions of the standard relativistic perturbation scheme by performing two steps: first, by restricting our solutions to perturbations that propagate on a flat unperturbed background spacetime and, second, by transforming to Eulerian background coordinates with truncation of nonlinear terms.

DOI: 10.1103/PhysRevD.96.123538

# I. INTRODUCTION

In this series of papers we generalize the Lagrangian perturbation theory to relativistic cosmology by employing spatial Cartan coframes as the only perturbation variable. In [1] we investigated 3 + 1 first-order solutions for the trace and antisymmetric parts of perturbations of a dust fluid, and we proposed an extrapolation in the spirit of Zel'dovich's approximation [2-4].<fnref n="1" /> In [10] the resulting nonperturbative scheme of structure formation was then applied to quantify average properties of inhomogeneous cosmologies in relation to the dark energy and dark matter problems. In [11] we gave the general perturbation and solution schemes at any order of the perturbations that cover the full hierarchy of the Newtonian Lagrangian perturbation and solution schemes.

In this paper we proceed with a detailed discussion of the general first-order scheme concentrating on the trace-free symmetric parts of the perturbations and, thus, focusing on gravitational waves propagating in continuous matter.

The recent detections of gravitational waves by the LIGO/VIRGO teams [12] has heated up this research subject, and we here provide further steps toward a nonlinear comprehension of the theory of gravitational waves; for a historical account and references see [13]. As in previous work we restrict our attention to irrotational dust continua for simplicity. The generalization to more general matter models and general foliations of spacetime is scheduled.

We recall the conceptual differences of our framework compared with the standard approach in relativistic cosmology (for a list of key references on the standard approach as well as work that relates to our approach see [11]). We choose the (spatially diffeomorphism invariant) Cartan coframes as a single perturbation variable, not the metric; we seek to describe the evolution of the perturbation fields intrinsically, thus operating on the physical Riemannian manifold, not on a global background manifold. As a consequence, perturbations cannot be expressed in terms of inertial (Eulerian) coordinates of the background, but in terms of local coordinates in the tangent spaces at each point of the physical manifold. These local coordinates take the role of Lagrangian coordinates in a flow-orthogonal (comoving) foliation of spacetime, adopted in this paper. With this approach we expect to capture nonlinear properties of gravitational waves, since a perturbative solution of these latter propagate in the perturbed spacetime, not in the background spacetime. This holds true even for linearized solutions for the Cartan coframes that we will consider

<footnote n="1">¹Univ Lyon, Ens de Lyon, Univ Lyon 1, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574,
F-69007 Lyon, France</footnote>

<footnote n="2">²Max-Planck-Institut für Gravitationsphysik, Albert-Einstein-Institut, Am Mühlenberg 1,
D-14476 Potsdam, Germany</footnote>

<footnote n="3">³Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, Massachusetts 02138, USA</footnote>

<footnote n="4">⁴Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Strasse 1, D-85741 Garching, Germany
(Received 6 November 2017; published 29 December 2017)</footnote>

<footnote>*Corresponding author.
buchert@ens-lyon.fr

¹In Newtonian cosmology in its Lagrangian form [5,6] this extrapolation is obtained by employing general functional expressions of the perturbed variable (see also the original proposal [7] and the discussions in [8,9]; recall that a relativistic generalization of Zel'dovich's approximation can only be identified with a Lagrangian perturbation solution, if the full system of equations can be expressed in terms of a single variable [1]).</footnote>
throughout this paper. As for the gravitoelectric part of the relativistic solutions, studied in [11], we here will also demonstrate how the solutions of the gravitomagnetic part of the standard perturbation approach are recovered. The reader may consult the introduction to [11] for further explanation of this approach.

This paper is structured as follows. Section II recalls the basic system of equations and provides the general perturbation scheme up to the first order; we rewrite this scheme in a couple of steps by employing decompositions in terms of symmetry and a split into gravitoelectric and gravitomagnetic parts; we then summarize the equations for each part in terms of master equations. Section III provides an alternative global view on the perturbation scheme by employing Hodge-de Rham theory. Section IV performs a detailed comparison with the standard first-order approach to gravitational waves. Section V discusses the results with an outlook. In Appendix A we provide the Lagrangian propagation equations for the gravitoelectric and gravitomagnetic parts of the Weyl tensor.

## II. FOUNDATIONS AND FIRST-ORDER PERTURBATION SCHEME

First, we will provide a brief recapitulation of the basic system of equations and the perturbation ansatz.

### A. Foundations

As in previous work, we employ a flow-orthogonal spacetime foliation admitting the following bilinear metric forms:

$$^{(4)}\mathbf{g} = -\mathbf{d}t \otimes \mathbf{d}t + ^{(3)}\mathbf{g}, \quad \text{with} \quad ^{(3)}\mathbf{g} = g_{ij} \mathbf{d}X^i \otimes \mathbf{d}X^j, \qquad (1)$$

where we denote by $X^i$ Gaussian normal (or Lagrangian) coordinates. The spatial metric form is decomposed into Cartan coframes,

$$^{(3)}\mathbf{g} = G_{ab} \eta^a \otimes \eta^b, \qquad (2)$$

where $G_{ab}(\mathbf{X})$ is Gram's matrix that encodes the initial metric perturbations [see Eq. (9) below]. The dynamical variable is composed of a $t$ parametrization of three spatial one-form fields.<fnref n="2" />

The system of 3 + 1 coefficient equations for the Einstein equations of an irrotational dust model, projected on the exact cotangent basis $\mathbf{d}X^i$ reads<fnref n="3" />

$$G_{ab} \dot{\eta}_{[i}^{a} \eta_{j]}^{b} = 0; \qquad (3)$$

$$\frac{1}{2J} \epsilon_{abc} \epsilon^{ikl} \partial_l (\dot{\eta}^a_j \eta^b_k \eta^c_l) = -\mathcal{R}^i_j + (4\pi G\mathcal{Q} + \Lambda) \delta^i_j; \qquad (4)$$

$$\frac{1}{2J} \epsilon_{abc} \epsilon^{mjk} \dot{\eta}_m^a \eta_j^b \eta^c_k = -\frac{\mathcal{R}}{2} + (8\pi G\mathcal{Q} + \Lambda); \qquad (5)$$

$$\left(\frac{1}{J} \epsilon_{abc} \epsilon^{ikl} \dot{\eta}_j^a \eta^b_k \eta^c_l\right)_{\parallel i} = \left(\frac{1}{J} \epsilon_{abc} \epsilon^{ikl} \dot{\eta}_i^a \eta^b_k \eta^c_l\right)_{\parallel j}; \qquad (6)$$

an overdot denotes a partial derivative with respect to the coordinate time (being equivalent to the covariant time derivative in the present setting); a double vertical slash denotes the covariant spatial derivative with respect to the 3-metric, the connection is assumed to be Levi Cività; the dust density is given by $\mathcal{Q} = Q_i J^{-1}$, $J \ge 0$, where the index $i$ marks initial data; $J$ defines the coefficient function of the 3-volume form, $J \equiv \sqrt{g}/\sqrt{G}$, with $\sqrt{g} \mathbf{d}^3 X$ the 3-volume form on the exact basis, $g := \det(g_{ij}(\mathbf{X}, t))$ and $G := \det(G_{ij}) = \det(g_{ij}(\mathbf{X}, t_i))$; $\mathcal{R}_{ij}$ denote the coefficients of the spatial Ricci tensor, with $\mathcal{R}$ its trace; the constant $\Lambda$ is the cosmological constant.

The set of equations { (3)–(6) } is composed of 13 equations, where the first 9 equations are the needed evolution equations for the 9 coefficient functions of the 3 Cartan coframe fields, while the remaining 4 equations originate from the energy and momentum constraints.

A key insight underlying previous work is the possibility of paraphrasing formalisms and results of the Newtonian description in Lagrangian form [5,6,14] that allowed construction of a part of the relativistic perturbation solutions in [11].<fnref n="4" /> This insight is based on the formal analogy between the Lagrange-Newton system of equations and the following set of 4 equations, in which the Ricci curvature is eliminated. It is composed of the first 3 Eqs. (3), and Raychaudhuri's equation [derived from the trace of the equation of motion (4) combined with the energy constraint (5)],<fnref n="5" />

$$\begin{gathered} G_{ab} \dot{\eta}_{[i}^{a} \eta_{j]}^{b} = 0; \\ \frac{1}{2J} \epsilon_{abc} \epsilon^{ikl} \dot{\eta}_{i}^{a} \eta_{j}^{b} \eta^{c}_{\ell} = \Lambda - 4\pi G\mathcal{Q}. \end{gathered} \qquad (7)$$

<footnote n="2">²Note that in place of orthonormal coframes we employ more general *adapted coframes*; see [10,11] for more details. Restricting the orthonormal spatial coframes $\eta^a_{on}$ to exact forms $\eta^a_{on} = \mathbf{d}f^a$ implies for the 3-metric a spatial diffeomorphism equivalence with a Euclidean space. This statement is obvious for orthonormal coframes, but it also remains true for our adapted coframes and the requirement of being exact, as long as we assume exactness for the orthonormal coframes: the reader may consult Appendix A in [11] for the proof of equivalence of this property for both types of coframes. Requiring exactness of the coframes also implies that the counterindices (denoted by $a, b, c, \dots$) can then be used as coordinate indices (denoted by $i, j, k, \dots$); the functions $f^a$ then define coordinates $x^i = f^{a \to i}$.</footnote>

<footnote n="3">³We here only give the equations in coordinate components; for their representation in terms of differential forms the reader is directed to [1,11].</footnote>

<footnote n="4">⁴This part, however, completely covers all Newtonian solutions at any order, whose explicit expressions in this representation were given up to the fourth order [15–19].</footnote>

<footnote n="5">⁵Equations (3) imply $G_{ab} \dot{\eta}_{[i}^{a} \eta_{j]}^{b} = 0$.</footnote>
For integrable coframe coefficients (or exact forms) the above system is closed and equivalent to the Lagrange-Newton system. A short discussion of the interrelation between the Lagrange-Newton system and the above gravitoelectric part of the Lagrange-Einstein system can be found in [20].

## B. Perturbation ansatz

We follow the prescription of [11] and decompose the coframes into a Friedmann-Lemaître-Robertson-Walker (FLRW) coframe and deviations thereof, the perturbations being locally developed to a given order $n$,

$$ \eta^a = \eta^a_i \mathbf{d}X^i = a(t) \left( \delta^a_i + \sum_n P^a_i^{(n)} \right) \mathbf{d}X^i, \quad (8) $$

in the local basis $\mathbf{d}X^i$. We use the normalization $a(t_i) = a_i = 1$. The local metric coefficients are given by$^6$

$$ g_{ij} = G_{ab}\eta^a_i\eta^b_j; \quad G_{ij} := g_{ij}(t_i), \quad (9) $$

with $G_{ij} = G_{ab}\eta^a_i(t_i)\eta^b_j(t_i) = G_{ab}\delta^a_i\delta^b_j$.

As the coframes are the only dynamical variables, we are entitled to employ functional definitions of all other variables such as the density, the metric, or the curvature in terms of the deformation of the manifold at a given order. The strategy followed in previous work of this series was to inject the $n$-order coframes into these functional definitions without further truncation, furnishing nonperturbative approximations that may improve iteratively by going to higher order deformations. In this paper, however, we do not exploit this strategy and concentrate on strictly linearized equations.

The deformation coefficients $P^a_i$ only appear summed over the noncoordinate index in the equations, so we introduce the following tensor coefficients and their trace:

$$ P^i_j \equiv \delta^a_i P^a_j \quad \text{and} \quad P \equiv P^k_k = \delta^k_a P^a_k, \quad (10) $$

and we use this notation throughout this paper. Confining ourselves to first-order deviations, we have for the coefficients of the perturbation one-forms,

$$ \eta^a_i = a(t)(\delta^a_i + P^a_i), \quad (11) $$

and we have for the metric coefficients (9),$^7$

$$ g_{ij} = a^2(t)(G_{ij} + 2P_{(ij)} + P_{ai}P^a_j); \quad (12) $$

we have defined

$$ P_{ij} := G_{ai}P^a_j. \quad (13) $$

For further details, especially concerning the initial data, the reader is directed to [11]. However, we here briefly recall the notations for the general initial data set.

We have chosen to prescribe initial data (being considered, without loss of generality, first order) in terms of the 6 one-form fields $\mathbf{U}^a = U^a_i \mathbf{d}X^i$ and $\mathbf{W}^a = W^a_i \mathbf{d}X^i$, being nonintegrable generalizations of the Newtonian peculiar-velocity and -acceleration gradients, and the initial values of the perturbations $[P^a_i(t_i) := P^a_i]$ together with their time derivatives,$^8$

$$ \left\{ \begin{array}{l} \mathcal{P}^a_i = 0; \\ \dot{\mathcal{P}}^a_i = U^a_i, U_{[ij]} = 0; \\ \dot{\mathcal{P}}^a_i = W^a_i - 2H_i U^a_i, W_{[ij]} = 0, \end{array} \right. \quad (14) $$

with the Hubble function $H = \dot{a}/a$ at initial time, $H_i = \dot{a}_i$ (we recall the normalization of the scale factor to be $a_i = 1$). We use the abbreviations $\delta^k_a U^a_k := U$, $\delta^k_a W^a_k := W$ for the trace expressions. The energy and momentum constraints have to be obeyed for admissible initial data,

$$ \left\{ \begin{array}{l} H_i U = -\frac{\mathcal{R}}{4} - W; \\ (U^a_j \delta^a_i)_{\parallel i} = (U^a_i \delta^a_j)_{\parallel j}, \end{array} \right. \quad (15) $$

with the initial Ricci curvature as found from the equation of motion (4),

$$ \mathcal{R}_{ij}^{(1)}(t_i) := \mathcal{R}_{ij} = -(W_{ij} + H_i U_{ij}) - \delta_{ij}(W + H_i U). \quad (16) $$

These initial data determine the problem completely: in our 3+1 split, the system of equations of motion for the coframes $\{(3)-(6)\}$ is composed of the 6 second-order differential Eqs. (4). As $U_{ij}$ and $W_{ij}$ are symmetric 3-matrices [originating from Eqs. (3)], they provide the 12 initial conditions needed for the first and second derivatives of the coframes. For both, the initial data and the coframes there are 4 constraints.

## C. First-order equations

The first-order Lagrange-Einstein system [the index (1) is omitted for the perturbation fields but kept for the Ricci curvature] read

$$ \dot{P}_{[ij]} = U_{[ij]}a^{-2} = 0; \quad (17) $$

<footnote>$^6$Recall that with this ansatz we choose to perturb a zero-curvature FLRW model, since it is possible to encode an initial first-order constant curvature in the coefficient functions $G_{ab}$ [11].</footnote>

<footnote>$^7$Note that only the zero-order metric tensor will appear in the linearized equations.</footnote>

<footnote>$^8$$\dot{P}^a_i$ or $\ddot{P}^a_i$ is a shorthand of $\dot{P}^a_i(t_i)$ or $\ddot{P}^a_i(t_i)$, respectively.</footnote>
$$ \ddot{P}_{ij} + 3H\dot{P}_{ij} = -a^{-2}\left(\mathcal{R}_{ij}^{(1)} - \frac{\mathcal{R}^{(1)}}{4}a^2\delta_{ij}\right); \quad (18) $$

$$ H\dot{P} + 4\pi G_{QH_i} a^{-3} P = -\frac{\mathcal{R}^{(1)}}{4} - a^{-3} W; \quad (19) $$

$$ \dot{P}^i_{[i|j]} = 0. \quad (20) $$

A more transparent representation of these equations decomposes the deformation fields into trace, trace-free symmetric, and antisymmetric parts,

$$ P_{ij} = P_{(ij)} + P_{[ij]} = \frac{1}{3} P \delta_{ij} + \Pi_{ij} + \mathfrak{P}_{ij}, \quad (21) $$

where we defined $\Pi_{ij} := P_{(ij)} - 1/3P\delta_{ij}$, $\mathfrak{P}_{ij} := P_{[ij]}$, and we introduce the trace-free symmetric part of the Ricci tensor, $\tau_{ij}^{(1)} := \mathcal{R}_{ij}^{(1)} - 1/3\mathcal{R}^{(1)}\delta_{ij}$. The first-order system for the deformation coefficients now reads

$$ \dot{\mathfrak{P}}_{ij} = U_{[ij]} a^{-2} = 0; \quad (22) $$

$$ \ddot{P} + 3H\dot{P} = -\frac{\mathcal{R}^{(1)}}{4}; \quad (23) $$

$$ \ddot{\Pi}_{ij} + 3H\dot{\Pi}_{ij} = -a^{-2}\tau_{ij}^{(1)}; \quad (24) $$

$$ H\dot{P} + 4\pi G_{QH_i} a^{-3} P = -\frac{\mathcal{R}^{(1)}}{4} - a^{-3} W; \quad (25) $$

$$ \frac{1}{3}\dot{P}_{|j} - \frac{1}{2}(\dot{\Pi}_{j}^{i})_{|i} = 0. \quad (26) $$

In order to solve the first-order trace equation and the traceless symmetric equation, it is necessary to express the first-order scalar curvature and the traceless Ricci tensor $^{(1)}\tau_{ij}$. To do so, we inject the metric and its inverse, truncated to first-order,

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

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

into the definitions of the spatial Christoffel symbol and spatial Ricci tensor, also making use of the momentum constraints, to obtain

$$ \Gamma^k_{ij}{}^{(1)} = \frac{1}{2} \delta^{kl} (G^{(1)}_{li|j} + G^{(1)}_{lj|i} - G^{(1)}_{ij|l}) \quad (29) $$

$$ + \delta^{kl} (P_{(li)|j} + P_{(lj)|i} - P_{(ij)|l}); \quad (30) $$

$$ \mathcal{R}_{ij}^{(1)} = G_{i[k|j]}^{(1)} + G_{[j|k]i}^{(1)} + P_{i[k|j]}^{k} + P_{j[k|i]}^{k}; \quad (31) $$

$$ \mathcal{R}^{(1)} = 2a^{-2}G_{[k|l]}^{(1)|k} := \frac{\mathcal{R}}{a^2}. \quad (32) $$

Using the split into parts with different symmetries, we also express the curvature through the parts P and $\Pi_{ij}$. A further manipulation, using again the momentum constraints, allows one to rewrite some terms in the curvature to obtain

$$ \mathcal{R}_{ij}^{(1)} = \mathcal{R}_{ij} + P_{|ij} - \frac{1}{3} P^{|k}_{|k} \delta_{ij} - \Pi_{ij}^{|k}_{|k}; \quad (33) $$

$$ \mathcal{R}^{(1)} = a^{-2}\mathcal{R}. \quad (34) $$

We can now write the first-order traceless part of the Ricci curvature tensor as

$$ \tau_{ij}^{(1)} = \mathcal{T}_{ij} + P_{|ij} - \frac{1}{3} P^{|k}_{|k} \delta_{ij} - \Pi_{ij}^{|k}_{|k}, \quad (35) $$

where we defined $\mathcal{T}_{ij} := \tau_{ij}(t_i) := \mathcal{R}_{ij} - 1/3\mathcal{R}\delta_{ij}$. These expressions for the intrinsic curvature show that first-order perturbations only generate traceless curvature in the course of time, while the trace is a conformal rescaling of the initial scalar curvature.

Knowing the first-order scalar curvature and the traceless Ricci tensor, we can inject them into the first-order system for the deformation coefficients (22)–(26). In this step we also perform the time integration of the (trivial) antisymmetric part and of the momentum constraints, and apply the constraints on initial data (14) and (15),<fnref n="9" />

$$ \mathfrak{P}_{ij} = \mathfrak{P}_{ij}(t_i), \quad \mathfrak{P}_{ij}(t_i) = 0; \quad (36) $$

$$ \ddot{P} + 3H\dot{P} = -\frac{a^{-2}\mathcal{R}}{4}; \quad (37) $$

$$ \ddot{\Pi}_{ij} + 3H\dot{\Pi}_{ij} - a^{-2}\Pi_{ij}^{|k}_{|k} = -a^{-2}\left(\mathcal{T}_{ij} + P_{|ij} - \frac{1}{3}P^{|k}_{|k}\delta_{ij}\right); \quad (38) $$

$$ H\dot{P} + 4\pi G_{QH_i} a^{-3} P = -\frac{a^{-2}\mathcal{R}}{4} - a^{-3} W; \quad (39) $$

$$ \frac{2}{3} P_{|j} = \Pi_{j|k}^{k}. \quad (40) $$

In the following subsection, we will separate the symmetric traceless deformation field $\Pi_{ij}$ into a part that is (through the momentum constraints) determined to have the same time dependence as the trace P, and a part that encodes free gravitational waves.

<footnote n="9">⁹Note that there is no constant of integration appearing in Eq. (40), because we chose our initial data in (14) such that $P^a = 0$.</footnote>
## D. First-order master equations

In this subsection, evolution and constraint equations will be combined into a set of equations for the symmetric parts of the first-order perturbation fields (the antisymmetric part is trivially solved). The result is a set of equations that we will call master equations.

### 1. Master equations for the trace and traceless symmetric parts

Inserting the energy constraint (39) into the evolution Eq. (37) provides the master equation for the trace part (Raychaudhuri's equation),

$$ \ddot{P} + 2H\dot{P} - 4\pi G_{Hij}a^{-3}P = a^{-3}W. \quad (41) $$

For the example of an Einstein-de Sitter universe model (EdS), the solution reads

$$ P = \frac{3}{5} \left[ \left( U t_i + \frac{3}{2} W t_i^2 \right) a - \left( U t_i - W t_i^2 \right) a^{-\frac{3}{2}} - \frac{5}{2} W t_i^2 \right]. \quad (42) $$

In the same spirit, we combine the momentum constraints (40) with the traceless part of the evolution Eq. (38). We eliminate the initial traceless curvature via the initial data set (16) and obtain the following master equation for the traceless symmetric part:

$$ \ddot{\Pi}_{ij} + 3H\dot{\Pi}_{ij} - a^{-2}(W_{ij}^{tl} + H_i U_{ij}^{tl}) \\ = a^{-2} \left( \Pi_{ij}^{|k} + \frac{1}{2}\delta_{ij}\Pi_{l|k}^k - \frac{3}{2}\Pi_{j|ki}^k \right), \quad (43) $$

with $^{tl}$ denoting the traceless part of the initial tensor fields. We note that the trace and the traceless part of the perturbation fields are coupled via the momentum constraints (40). Thus, their dynamics is not independent. Understanding how they are coupled is in the focus of the next subsection.

### 2. The gravitoelectric and gravitomagnetic parts

We now split the traceless symmetric part into a component that is coupled to the trace, $^E\Pi_{ij}$, and one that is decoupled from it, $^H\Pi_{ij}$. As was explicitly shown in [11], $^E\Pi_{ij}$ can be obtained from a generalization of the Newtonian trace solution. We will now show that it is separable in space and time and has the same time dependence as the trace P. Furthermore, we will show that $^E\Pi_{ij}$ and $^H\Pi_{ij}$ are closely related to the gravitoelectric and gravitomagnetic parts of the spatially projected Weyl tensor (hence our index notations E and H).

From the momentum constraints (40) we infer that a decoupling of $^H\Pi_{ij}$ from P implies that it must be divergence-free. We write (40) as

$$ \frac{2}{3} P_{|j} = \Pi_{j|k}^k := ({}^E\Pi_{j|k}^k + {}^H\Pi_{j|k}^k), \quad (44) $$

which can then be divided into two separate constraints,

$$ \frac{2}{3} P_{|j} = {}^E\Pi_{j|k}^k; \quad {}^H\Pi_{j|k}^k := 0. \quad (45) $$

We split the initial data accordingly,

$$ U_{ij} := {}^E U_{ij} + {}^H U_{ij}; \quad W_{ij} := {}^E W_{ij} + {}^H W_{ij}. \quad (46) $$

The initial traceless part of the Ricci curvature reads

$$ W_{ij}^{tl} + H_i U_{ij}^{tl} := ({}^E W_{ij}^{tl} + H_i {}^E U_{ij}^{tl}) + ({}^H W_{ij}^{tl} + H_i {}^H U_{ij}^{tl}). \quad (47) $$

The propagation of the constraint Eqs. (44) guarantees the preservation of the decomposition (45) in time.

This split can be carried through to the master Eq. (43). Inserting the superposition $\Pi_{ij} = {}^E\Pi_{ij} + {}^H\Pi_{ij}$, we first extract the divergence-free part, which obeys the master equation for the gravitomagnetic part,

$$ {}^H \ddot{\Pi}_{ij} + 3H^H \dot{\Pi}_{ij} - a^{-2} {}^H \Pi_{ij}^{|k}{}_k = a^{-2} ( {}^H W_{ij}^{tl} + H_i {}^H U_{ij}^{tl} ). \quad (48) $$

This equation describes the propagation of gravitational waves and has the form of d'Alembert's equation with a damping term due to expansion—note the important difference to the standard perturbation approach, discussed in Sec. IV: this equation assumes this form in the local coordinates of the tangent spaces at each point of the manifold.

The gravitoelectric part $^E\Pi_{ij}$ is the solution of

$$ \begin{aligned} & {}^E\ddot{\Pi}_{ij} + 3H^E \dot{\Pi}_{ij} - a^{-2}( {}^E W_{ij}^{tl} + H_i {}^E U_{ij}^{tl}) \\ & = a^{-2} \left( {}^E\Pi_{ij}^{|k}{}_k + \frac{1}{2}\delta_{ij} {}^E\Pi_{l|k}^k{}_l - \frac{3}{2} {}^E\Pi_{j|ki}^k \right). \end{aligned} \quad (49) $$

As was explicitly shown in [11], and as can be obtained from Eq. (45), $^E\Pi_{ij}$ has the same time dependence as the trace part P. It is therefore the solution of the master equation for the gravitoelectric part, which is formally analogous to (37) [using (15)],

$$ {}^E\ddot{\Pi}_{ij} + 3H^E \dot{\Pi}_{ij} = a^{-2} ({}^E W_{ij}^{tl} + H_i {}^E U_{ij}^{tl}). \quad (50) $$

The equivalence of (49) and (50) implies with (45)

$$ \Delta_0 {}^E\Pi_{ij} = \mathcal{D}_{ij} P, \quad (51) $$

where $\Delta_0$ denotes the ordinary Laplacian with respect to the local Lagrangian coordinates. The traceless spatial double derivative operator is defined as follows:

$$ \mathcal{D}_{ij} = \partial_i \partial_j - \frac{1}{3} G_{ij} \Delta_0. $$
(Here, $G_{ij} \approx \delta_{ij}$.) In the next subsection we are going to discuss some important issues of interpretation of the solutions to the above equations. In particular, we recall the relation between the Newtonian Lagrangian perturbation theory and the gravitoelectric part $E\Pi_{ij}$, we discuss the role of the gravitomagnetic part $H\Pi_{ij}$ for the dynamical evolution of the traceless Ricci tensor, and we motivate the need for global considerations.

### E. Discussion

Interesting insights into the gravitoelectric part of the perturbation field, denoted by $E^P_{ij}$, can be obtained from a formal analogy between the Lagrange-Newton system and the gravitoelectric part of the Einstein equations—see [11]. This set of equations, which comprises (7) and which does not explicitly involve the spatial Ricci curvature tensor, can be expressed in terms of identities on the trace and the antisymmetric part of the gravitoelectric part of the Weyl tensor [cf. Eqs. (33) of [11]]. In the comoving synchronous foliation of spacetime considered here, this nonclosed set of equations becomes identical to the closed Lagrange-Newton system after the execution of the “Minkowski restriction” (MR). This mathematical operation sends Cartan coframe fields $\eta^a$ to exact forms $d\mathbf{f}^a$ (the reader may recall its definition in [11] and note that the perturbation field $E^P_{ij} = 1/3P\delta_{ij} + E\Pi_{ij}$ coincides with the gravitoelectric solution of [11]). Reverting the MR allows us to build the relativistic counterpart of Newtonian solutions to the Lagrangian perturbation theory,

$$d\mathbf{f}^i = a(t)(d\mathbf{X}^i + d\mathbf{P}^i) \rightarrow \eta^a = a(t)(d\mathbf{X}^a + E^P \mathbf{P}^a). \quad (52)$$

Although $E^P_{ij}$ features a Newtonian time behavior, the spatial coefficient functions are a priori nonintegrable in the relativistic case, so that its trace and its trace-free part separately generate curvature. However, we emphasize a subtle property of the gravitoelectric perturbation field: by comparing Eqs. (49) and (50), leading to the constraint (51), with the defining equation for the Ricci tensor (33), we observe that the gravitoelectric part $E\Pi_{ij}$ exactly compensates the nonconstant part of the traceless Ricci curvature generated by $P$. Thus, only the gravitomagnetic perturbation $H\Pi_{ij}$ encodes the dynamics of the Ricci curvature or its traceless part (35), respectively.<fnref n="10" />

$$\Delta_0 H_{ij} (E^P_{ij}) = 0; \quad (53)$$

i.e., there is in general a harmonic global part of the solution, an issue that we will address in the next section.

Furthermore, we point out that the master equations for both perturbation fields require global conditions to be solved, which leads us to consider the topology of the spatial hypersurfaces. In the next section, we will therefore adopt a global description of the spatial sections. Results of the Hodge-de Rham theory will be used to investigate the global properties of the Cartan coframes and to relate their Hodge decomposition to the local decomposition $P_{ij} = E^P_{ij} + H^P_{ij}$. Local considerations do not fully determine the physical content of $E^P_{ij}$ and $H^P_{ij}$, as the elliptic equations we encounter need topological boundary conditions to be solved.

### III. GLOBAL CONSIDERATIONS: TOPOLOGY AND HODGE-DE RHAM THEORY

In the last section we considered the projection of the first-order equations on a local coordinate basis {$d\mathbf{X}^i$}, hence constraining the validity of the results to the regions of the manifold that can be covered by a single coordinate chart. Topological considerations enable us to specify the boundary conditions needed to globally solve the governing equations. Note that the need for topological considerations is not specific to the intrinsic description of perturbation fields defined on perturbed spatial sections; boundary conditions have to be specified too in standard perturbation theory where the perturbation fields are described to evolve on flat space sections. In this latter case the topology usually adopted is that of a 3-torus (periodic boundary conditions).

In the following subsection we will consider spatial hypersurfaces that have a closed topology. Under this assumption it is possible to apply the Hodge decomposition of the Cartan coframes. As we will see, this decomposition will provide a new understanding of the global properties of perturbations.

#### A. Hodge theorem and Thurston's geometrization program

In order to set the stage for the Hodge decomposition of forms, some additional formalism is necessary and summarized in this subsection. For introductions into differential geometry and Hodge's theorem [21], the reader may refer to [22–25].

##### 1. Laplace-de Rham operator and the Weizenböck formula for a 1-form

Let $\mathcal{M}$ be a closed Riemannian $n$-dimensional manifold, and let us introduce two $p$-forms $\alpha$ and $\beta$. We define a global positive definite inner product,

$$(\alpha, \beta) = \int_{\mathcal{M}} \alpha \wedge *\beta, \quad (54)$$

where the dual of the form $\beta$ in an $n$-dimensional space is $*\beta$,

<footnote n="10">¹⁰This holds up to a global contribution of the gravitomagnetic part that is generated by $E\Pi_{ij}$. To understand this latter remark we note that the gravitomagnetic part generated by $E^P_{ij}$ solves the following equation (cf. [11], Appendix C):</footnote>
$$ *\beta = \frac{1}{(n-p)!} \epsilon_{\nu_1 \cdots \nu_p} \mu_1 \cdots \mu_{n-p} \beta^{[\nu_1 \cdots \nu_p]} \eta^{\mu_1} \wedge \cdots \wedge \eta^{\mu_{n-p}}; \qquad (55) $$

$\{\eta^i\}$, $i = 1 \cdots n$ are the Cartan coframes defined on an $n$-dimensional manifold $M$.

The Laplace-de Rham operator $\Delta^{\text{dR}}$ is a mapping from $p$-forms to $p$-forms; it generalizes the simple Laplacian from flat to curved manifolds. It is defined by

$$ \Delta^{\text{dR}} := \mathbf{d}\mathbf{d}^* + \mathbf{d}^*\mathbf{d} = (\mathbf{d} + \mathbf{d}^*)^2, \qquad (56) $$

where the codifferential operator $\delta = \mathbf{d}^*$ acting on the $p$-form $\beta$ is defined, on closed Riemannian manifolds, as follows:

$$ \mathbf{d}^*\beta = (-1)^{n(p+1)+1} * \mathbf{d} * \beta. \qquad (57) $$

According to the Weizenböck formula, $\Delta^{\text{dR}}$ applied to three 1-form fields $\eta^a$, which is the case of interest for the Cartan coframes, yields in components, projected onto the exact basis $\{\mathbf{d}X^i\}$,

$$ (\Delta^{\text{dR}}\eta^a)_i = (\Delta\eta^a)_i + \eta^a_k \mathcal{R}_i^k, \qquad (58) $$

where $\mathcal{R}_i^k$ are the 3-Ricci curvature tensor components. The first term features the “rough Laplacian”,<fnref n="11" /> which implicitly depends on the geometry through the covariant derivatives, and the second term explicitly takes into account the local geometry through the Ricci curvature.

## 2. The Hodge theorem

The Hodge theorem [21] asserts that on a closed, oriented, and smooth Riemannian $n$-dimensional manifold, equipped with a smooth metric, the vector space of harmonic $p$-forms<fnref n="12" /> $\mathcal{H}^p$ is of finite dimension. The generalized Poisson equation for $p$-forms $\alpha^p$ and $\rho^p$ on a curved manifold,

$$ \Delta^{\text{dR}}\alpha^p = \rho^p, \qquad (59) $$

has a solution $\alpha^p$ if and only if $\rho^p$ is orthogonal to $\mathcal{H}^p$.

For any $p$-form $\gamma^p$ there exists a $p$-form $\theta^p$ such that $\gamma^p$ can be Hodge-decomposed as follows:

$$ \gamma^p = \mathbf{d}\alpha^{p-1} + \mathbf{d}^*\beta^{p+1} + \mathbf{h}^p, \qquad (60) $$

where $\alpha^{p-1} = \mathbf{d}^*\theta^p$, $\beta^{p+1} = \mathbf{d}\theta^p$, and $\mathbf{h}^p$ is harmonic. In other words, the vector space of $p$-forms can be decomposed as

$$ \Lambda^p = \mathbf{d}\Lambda^{p-1} \oplus \mathbf{d}^*\Lambda^{p+1} \perp \mathcal{H}^p, \qquad (61) $$

where $\mathbf{d}\Lambda^{p-1}$, $\mathbf{d}^*\Lambda^{p+1}$, and $\mathcal{H}^p$ are the $p$-dimensional exact, coexact, and harmonic vector spaces.

In view of what follows, we remark that the manifold need not be simply connected for employing the Hodge decomposition.

## 3. Harmonic forms, and the geometry and topology of the Universe

In order to apply the Hodge decomposition (60), we first have to address the question of the harmonic part that, if nonvanishing, would make a decomposition into exact and coexact forms nonunique. In general relativity, for regular Einstein flows, it is generally assumed that the topology of spacetime is given by the topology of the initial (Cauchy) hypersurface $\Sigma$, being conserved in time (hyperbolicity of Einstein's equations), $\Sigma \times \mathbb{R}$. Hence, we have to specify the topology of the Cauchy surface.

If this initial hypersurface is flat, then there is a simple reasoning to remove the harmonic part [18]: in the case of a flat space, the harmonic part is a vector field that obeys the simple Laplace equation, $\Delta h^a = 0$. Solutions of the Poisson equation on a flat space are unique, if we impose periodic boundary conditions, i.e., we impose a 3-torus topology to close the spatial section. Harmonic vector fields are then spatially constant, and this constant can be set to zero without loss of generality by exploiting the translational invariance of the equations on a flat space (e.g., Newton's equations have the same physical content if we add to the velocity model a spatially constant, eventually time-dependent vector field at each point).

For curved manifolds the vanishing of harmonic forms cannot be based on simple arguments. There are theorems (see, e.g., [22]) stating that, if the oriented closed Riemannian manifold $M$ has everywhere a positive Ricci tensor, then a harmonic 1-form will vanish identically. This situation is, however, nongeneric in cosmology, since the Ricci curvature will alternate between positive and negative values. If generic inhomogeneities are present, one may (in line with Thurston's geometrization program [26], see below, and also [27]) conceive the idea of using regionally different backgrounds in our perturbative analysis. Here, “regional background” refers to length scales set by an average over a given region containing inhomogeneities governing the background; see [28]. Such a formulation would allow us to construct a link to Thurston's geometrization program that asserts that every closed three-dimensional manifold can be decomposed into a connected sum of 3-manifolds modeled after the 8 model geometries listed by Thurston. Thurston's program was proven by Perelman in 2003 using Ricci flow with surgery [29]. For its realization one has to devise a way to match the different “homogeneity patches,” a possibly difficult task. In such a case we would have different situations for the harmonic

<footnote n="11">¹¹We here denote by $\Delta$ the rough Laplacian (as the negative Bochner Laplacian) whose action on the coframe components can be written $-\eta_i^a \partial_i \partial^a$, where $\partial_i$ denotes the Lagrangian covariant derivative with respect to the spatial metric.</footnote>

<footnote n="12">¹²Since $(\Delta\alpha^p, \alpha^p) = \|\mathbf{d}^*\alpha\|^2 + \|\mathbf{d}\alpha\|^2 \ge 0$, $\Delta\alpha^p = 0$ if and only if $\mathbf{d}^*\alpha = 0$ and $\mathbf{d}\alpha = 0$. Harmonic forms are both closed and coclosed.</footnote>
forms (depending on the regional background). Perelman's work also implies a proof of Poincaré's conjecture: Thurston's geometrization program also asserts that closed and simply connected 3-manifolds are homeomorphic to the hypersphere $\mathbb{S}^3$. Since $\mathbb{S}^3$ admits the round metric (which has positive Ricci curvature), it can be shown that the only harmonic 1-form is the null-form.<fnref n="13" />

The assumption of simple connectivity has obvious limitations. For instance, it is valid for a constant-curvature space with positive curvature, or for infinite models with constant curvature (and trivial fundamental groups), but already for a toroidal model with zero curvature or negative constant curvature spaces that are finite, it fails. An alternating and inhomogeneous curvature is the generic situation. Notwithstanding, we henceforth assume an $\mathbb{S}^3$ topology for the spatial sections, such that the harmonic part vanishes. This assumption is to be considered as a representative example of space forms for which the vanishing of the harmonic part can be rigorously proved. (Recall that the vanishing of the harmonic part on a curved manifold generalizes the assumption made in Newtonian cosmology or standard perturbation theory.)

## B. Hodge decomposition of the Cartan coframes

As we have seen in the previous subsection, the Hodge theorem allows us to decompose the Cartan coframes into exact, coexact, and harmonic forms. In the case of a closed and simply connected manifold, e.g., $\mathbb{S}^3$, the harmonic 1-form in this decomposition is set to zero. Furthermore, since any closed three-dimensional manifold is parallelizable, the Cartan coframes are continuously defined on the manifold while the coordinate charts undergo singularities. This case allows us to extend the definition of the coframes globally and apply to them the Hodge decomposition.<fnref n="14" />

With these assumptions, the three Cartan coframes, which are three 1-forms, can be decomposed as follows:

$$ \eta^a = d\alpha^a + d^*\beta^a = \Delta^{dR}\gamma^a, \qquad (62) $$

where $\alpha^a$ are three scalars, $\beta^a$ are three 2-forms, and $\gamma^a$ are three 1-forms.

For later reference we project Eq. (62) onto our local Lagrangian basis $\{dX^i\}$, $g := \det(g_{ij})$,

$$ \eta^a = \alpha^a|_j dX^j + (\beta^a_{km} e^{pkm} g_{rp} \sqrt{g})|_l g_{jn} e^{nrl} \sqrt{g} dX^j, \qquad (63) $$

where the 2-form coefficient matrix $\beta^a_{km}$ occurs in the projection of the 2-form $\beta^a$ as $\beta^a = \beta^a_{km}dX^k \wedge dX^m$.

In three dimensions, each antisymmetric coefficient matrix $\beta^a_{km}$ has only three independent components. Writing these components as $B^a_v$, we can choose to express $\beta^a_{km}$ by $\beta^a_{km} = 1/2g^{vv} \sqrt{g} \epsilon_{wkm} B^a_v$. With this substitution the coefficient form of $\eta^a$ reads

$$ \eta^a = \alpha^a|_j dX^j + B^a_{u|l} \sqrt{\det(g)} g_{kj} e^{ulk} dX^j. \qquad (64) $$

## 1. Consequences for the intrinsic spatial curvature

The Hodge decomposition bears interesting consequences in the context of our decomposition into gravito-electric and gravitomagnetic parts. From (53) we conclude that, for topologies with a null-dimensional harmonic space of 1-forms, the gravitomagnetic part generated by $E^P a$ is zero. For this we note that the structure coefficients associated with the Cartan coframes,

$$ d\eta^a = -\frac{1}{2} C^a_{bc} \eta^b \wedge \eta^c, \qquad (65) $$

are related, at first order, to the (spatially projected) gravitomagnetic part of the Weyl tensor as follows:

$$ H^i_j = \frac{1}{2} \epsilon^i_{jl} \dot{C}^j_{lk}. \qquad (66) $$

Since the gravitoelectric perturbation $E^P a$ generates only a harmonic gravitomagnetic part (that vanishes according to our global assumption), and since from (66) the structure coefficients calculated from $E^P a$ are constant in time and are initially zero $[\eta^a(t_i) = dX^a]$, we then conclude from (65) that $E^P a$ represents the integrable part $d\alpha^a$ of the perturbation fields. This is in accord with our finding that the curvatures generated separately by the trace and the trace-free gravitoelectric parts exactly compensate each other; see Sec. II E. Thus, this global consideration confirms that the only curvature-generating part of the perturbation fields is $H^i \Pi^a$ (at first order).

Note that the scalar part of the curvature (32) is at first order not part of the perturbation fields, but included in the

<footnote n="13">¹³This is a direct consequence of the theorem by de Rham [30,31]: If $v$ is a harmonic vector field on a compact Riemannian manifold and if the Ricci curvature $Ric(v, v) \ge 0$, then $v$ is parallel and $Ric(v, v) = 0$. This in turn implies that there exists no nonzero harmonic field on a compact Riemannian manifold of positive Ricci curvature. Hence, since $\mathbb{S}^3$ admits the round metric, the first de Rham cohomology group $\mathcal{H}_{\mathrm{DR}}^1(\mathcal{M}) = 0$ [Hodge theory indeed states that the space of the harmonic 1-forms $\mathcal{H}_{\mathrm{DR}}^1(\mathcal{M})$ is isomorphic to the first de Rham cohomology group $\mathcal{H}_{\mathrm{DR}}^1(\mathcal{M})$] [32,33].</footnote>

<footnote n="14">¹⁴The $\mathbb{S}^3$ topology is a known example in cosmology for finite universe models with constant positive curvature such as the classical Einstein or Eddington universe models. Speaking in terms of constant-curvature template geometries for the Universe, this assumption may be linked to observations [a marginally spherical model [34] has been favored by Cosmic Microwave Background measurements of a reported first Doppler peak that is shifted by a few percent toward larger angular scales with respect to the peak predicted by the standard cold dark matter model [35]]. For the opposite case of a hyperbolic universe model with small negative curvature, see [36], and for recent constraints on the constant-curvature universe model from the Planck mission, see [37].</footnote>
normalization of the coframes $G_{ab}$. This defines the curved reference space for our first-order perturbation fields. In [11] we explicitly demonstrated that, with respect to this space, the time evolution of the perturbations correspond to the Newtonian solutions. With the decomposition (62) this can now be made more rigorous: the exact part can, as in the Newtonian case, be absorbed into a redefinition of global coordinates, leaving us this time not with a reparametrization of flat space, but rather of the initially perturbed metric $G_{ij} = g_{ij}(t_i)$. The additional component in the nonexact part that had not been considered in [11] encodes the nonintegrable geometrical deviations from the zero-order solution, featuring gravitational waves.

## 2. Obtaining the Hodge fields from the perturbed coframes

The Hodge decomposition of the Cartan coframes (62) can equivalently be defined for the perturbation fields. In our case these are taken to be of first order. Thus, we face the decomposition (not changing the notation for $\alpha^a$ and $\beta^a$ for simplicity),

$$ \mathbf{P}^a = \mathbf{d}\alpha^a + \mathbf{d}^*\beta^a. \quad (67) $$

We project these fields onto the local Lagrangian basis $\{\mathbf{d}X^i\}$, and using again the fact that we work at first order, we find

$$ P^a_j = \alpha^a_{\ |j} + 2\beta^a_{\ |kj} \delta^{jk}, \quad (68) $$

where $\alpha^a$ are three scalar fields and $\beta^a_{kj}$ the three 2-form coefficient matrices $\beta^a = \beta^a_{km} \mathbf{d}X^k \wedge \mathbf{d}X^m$. Note that we used the fact that we work at first order and so $P^a_j, \alpha^a$, and $\beta^a_{kj}$ are first-order quantities and $g_{ij} \approx a^2 \delta_{ij}$. In three dimensions, each antisymmetric coefficient matrix $\beta^a_{km}$ has only three independent components. Writing these components as $B^a_v$, we can express the matrix by $\beta^a_{km} = 1/2 g^{vw} \sqrt{g} \epsilon_{wkm} B^a_v$. With this substitution the coefficient form of $P^a_j$ reads

$$ P^a_j = \alpha^a_{\ |j} + \delta_{rj} \epsilon^{ulr} \mathbf{B}^a_{\ u|l}, \quad (69) $$

again at first order. Formally the coefficients $B^a_v$ are attributed to the 1-form fields $\mathbf{B}^a$ and related to $\beta^a$ through the duality of 2- and 1-forms in three dimensions.

Since $E^\mathbf{P}^a$ encodes the integrable part of the perturbations, $E^\mathbf{P}^a = \mathbf{d}\alpha^a$, it is straightforwardly related to the scalar modes of the standard perturbation approach (cf. Sec. IV A). The complementary part, $H^\Pi^a$, is equal to $\delta\beta^a$.

From (69) we conclude that the Hodge fields $\alpha_i$ and $B_{ij}$ can be obtained from the perturbation fields by solving the following Poisson equations in local (Lagrangian) coordinates:

$$ \Delta_0 \alpha_j = P_{ij}^{|j|}, \quad (70) $$

$$ \Delta_0 B_{ij} = \epsilon_j^{mr} P_{ir|m}, \quad (71) $$

where we have chosen the gauge condition on $B_{ij}$ to be $B^m_{i|m} = 0$.<fnref n="15" />

The spatial gravitomagnetic part of the Weyl tensor can be expressed in terms of the expansion tensor [cf. Eq. (70) of [1]],

$$ H^i_j = -\frac{1}{j} \epsilon^{ikl} \Theta_{jk|l}. \quad (72) $$

At first order, this yields

$$ H^i_j = -\frac{1}{a} \epsilon^{ikl} \dot{P}_{jk|l} \quad \text{and} \quad H_{ij} = -a \epsilon_i^{kl} \dot{P}_{jk|l}. \quad (73) $$

We conclude that the first-order gravitomagnetic part is given by

$$ H_{ij} = -a \Delta_0 \dot{B}_{ij}. \quad (74) $$

# IV. COMPARISON WITH THE STANDARD PERTURBATION SCHEME

In this section we will explain the link between our approach and the standard one that uses the scalar-vector-tensor (S-V-T) decomposition of the perturbed metric. In order to do so, we compare our results to the ones obtained in the comoving synchronous gauge being the choice of coordinates that initially coincides with the Lagrangian coordinates in the 3 + 1 foliation we consider. (At places where we refer to explicit solutions, we will adopt an EdS background solution as an example, without further indication in the text.)

## A. S-V-T decomposition in the comoving synchronous gauge

In most of the literature, the perturbative approach is applied to the metric instead of the coframes (as the more elementary variable). The synchronous line element is decomposed as

$$ ds^2 = -dt^2 + a^2(t) \gamma_{ij}(\mathbf{x}, t) dx^i dx^j. \quad (75) $$

The first-order corrections to a flat FLRW background are then written as

<footnote n="15">¹⁵This gauge degree of freedom comes from the fact that the physical content of $\mathbf{B}^a$ is determined by its relation to $\mathbf{P}^a$, given by Eq. (69). In other words, due to the Levi-Civita symbol, the term $B^m_{i|m}$ does not intervene in Eq. (69); we are free to choose its value. Whatever the value of $B^m_{i|m}$, it will equally contribute to $\mathbf{P}^a$.</footnote>
$$ \gamma_{ij} = \delta_{ij} + \gamma_S \delta_{ij}. \qquad (76) $$

In the standard perturbation framework, it is common to split the perturbations into scalar, vector, and tensor contributions, rather than to consider a split into gravito-electric and gravitomagnetic parts. The perturbation to the background can be decomposed as

$$ \gamma_{S \ ij} = -2\phi_S\delta_{ij} + \mathcal{D}_{ij}\chi_S^{\parallel} + \partial_i\chi_S^{\perp j} + \partial_j\chi_S^{\perp i} + \chi_{ij}^{\top}, \quad (77) $$

where we follow the notation of [38] and denote the traceless derivative by $\mathcal{D}_{ij} = \partial_i\partial_j - \frac{1}{3}\Delta\delta_{ij}$, where $\Delta$ denotes the ordinary flat-space Laplacian in Eulerian background coordinates, and

$$ \partial^i \chi_{S i}^{\perp} = 0; \quad \chi^{\top i}_i = 0; \quad \partial^i \chi_{ij}^{\top} = 0. \qquad (78) $$

As we consider only first-order corrections to the metric in this subsection, $\phi_S$, $\chi_S^{\parallel}$, $\chi_S^{\perp}$, and $\chi_{ij}^{\top}$ are first-order quantities as well.

The decomposition (77) simplifies for irrotational fluids. In this case the vector perturbations represent gauge modes that can be set to zero [38]. This result is in accord with the simplification of our coframe decomposition (21). Note that Eq. (21) does not contain any vector fields anyway, as we have a local coordinate system (in the tangent space at a point in the Riemannian manifold) in which vector fields simply are not present in a flow-orthogonal foliation. However, the part that corresponds to the vector perturbations, the antisymmetric tensor part, also vanishes for the irrotational case.

With the vectors removed, we only need the synchronous gauge expressions for $\phi_S$, $\mathcal{D}_{ij}\chi_S^{\parallel}$, and $\chi_{ij}^{\top}$ as given in [38]. Let us start with the tensor part. It satisfies the equation

$$ \ddot{\chi}_{ij}^{\top} + 3H\dot{\chi}_{ij}^{\top} - \frac{1}{a^2}\Delta\chi_{ij}^{\top} = 0. \qquad (79) $$

This equation corresponds to our Eq. (48) in the sense explained further below. The solution of (79) reads [38]

$$ \chi_{ij}^{\top}(\mathbf{x}, t) = \frac{1}{(2\pi)^3} \int d^3\mathbf{k} \exp(i\mathbf{k} \cdot \mathbf{x}) \chi_{\sigma}(\mathbf{k}, t) e_{ij}^{\sigma}(\hat{\mathbf{k}}), \qquad (80) $$

with $\hat{\mathbf{k}}$ a normalized $\mathbf{k}$-vector, and where $e_{ij}^{\sigma}(\hat{\mathbf{k}})$ introduces the polarization tensor; $\sigma$ is ranging over the polarization components + and $\times$. The $\chi_{\sigma}(\mathbf{k}, t)$ are the corresponding amplitudes. The time evolution is obtained from (79) and provides the solution in terms of the first spherical Bessel function, $j_1(x) = \sin x/x^2 - \cos x/x$,

$$ \chi_{\sigma}(\mathbf{k}, t) = A(k) a_{\sigma}(\mathbf{k}) \left( \frac{3j_1(3kt_i^{2/3}t^{1/3})}{3kt_i^{2/3}t^{1/3}} \right). \qquad (81) $$

$a_{\sigma}(\mathbf{k})$ denotes a zero mean random variable, and $A(k)$ encodes the form of the spectrum of primordial gravitational waves.

Second, we turn to the scalar sector. The solutions given in [38] are restricted to the growing mode only. In addition, the authors use the residual gauge freedom of synchronous gauge to fix $\chi_S^{\parallel}$ such that $\Delta\chi_{Si}^{\parallel} = -2\delta_i$. Then, their solutions for the scalar sector are given by

$$ \mathcal{D}_{ij}\chi_S^{\parallel} = -3t_i^2 \left(\frac{t}{t_i}\right)^{\frac{2}{3}} \left(\partial_{ij}\varphi - \frac{1}{3}\delta_{ij}\Delta\varphi\right), \qquad (82) $$

where $\varphi$ is time independent and defined by its relation to $\delta_i$ in the cosmological flat-space Poisson equation $\Delta\varphi(\mathbf{x}) = 2/(3t_i^2)\delta_i(\mathbf{x})$, and the solution for the second scalar mode

$$ \phi_S(\mathbf{x}, t) = \frac{5}{3}\varphi(\mathbf{x}) + \frac{1}{2}t_i^2 \left(\frac{t}{t_i}\right)^{\frac{2}{3}} \Delta\varphi(\mathbf{x}). \qquad (83) $$

## B. Comparing the S-V-T decomposition to the intrinsic perturbation fields

In our case, $P_{ij}$ as well as the initial metric $G_{ij}$ appear in the decomposition of the metric perturbation that we can formally compare,

$$ \gamma_{Sij} \leftrightarrow G_{ij} + 2P_{(ij)}. $$

The relation between the two sets of perturbation fields can be formally obtained by a replacement of the dependent variables,

$$ \phi_S \leftrightarrow -\frac{1}{3} \left( P + \frac{1}{2} G \right), \\ \frac{1}{2} \left( \mathcal{D}_{ij} \chi_S^{\parallel} + 2 \partial_i \chi_{Sj}^{\perp} + \chi_{ij}^{\top} \right) \leftrightarrow \left( \Pi_{ij} + \frac{1}{2} G_{ij}^{\top} \right). $$

Note, however, that this correspondence is formal in the sense that the right-hand-side expressions are given in terms of local (Lagrangian) coordinates on the manifold that are constant along the solution, while the left-hand-side expressions are given in terms of global (Eulerian) background coordinates.

In order to better compare the standard results to ours, we modify the restrictions imposed in [38]. First, we use a different gauge choice, which is $\Delta\chi_{Si}^{\parallel} = 0$ instead of $\Delta\chi_{Si}^{\parallel} = -2\delta_i$. Second, we include the decaying mode. Then, the potential is no longer independent of time and $\Delta\varphi(\mathbf{x}, t) = 2/(3t_i^2a)\delta(\mathbf{x}, t)$. Therefore, $\varphi$ now has two components $\varphi_1$ and $\varphi_2$ that are determined by $\Delta\varphi_1(\mathbf{x}) + a^{-5/2}\Delta\varphi_2(\mathbf{x}) = 2/(3t_i^2a)\delta(\mathbf{x}, t)$. Finally, we include the prefactor on the right-hand-side into the definition of a renormalized potential $\psi(\mathbf{x}, t) = \psi_1(\mathbf{x}) + a^{-5/2}\psi_2(\mathbf{x})$, which then satisfies $\Delta\psi_1(\mathbf{x}) + a^{-5/2}\Delta\psi_2(\mathbf{x}) = -a^{-1}\delta(\mathbf{x}, t)$.
With these changes, the resulting metric perturbations read

$$ \mathcal{D}_{ij} \chi_S^\parallel = 2((t/t_i)^{\frac{2}{3}} - 1) \mathcal{D}_{ij} \psi_1 + 2((t/t_i)^{-1} - 1) \mathcal{D}_{ij} \psi_2 $$

and

$$ \phi_S = -10/(9t_i^2)\psi_1 - \frac{1}{3}((t/t_i)^{\frac{2}{3}} - 1)\Delta\psi_1 - \frac{1}{3}((t/t_i)^{-1} - 1)\Delta\psi_2. $$

The scalar and tensor fields are formally linked, in the comoving synchronous gauge, to the first-order metric and intrinsic perturbation field components of $^{E}\mathbf{P}^a$ and $^{H}\mathbf{\Pi}^a$, providing a dictionary between our formalism and the standard notation,

$$ \begin{aligned} {}^E P_{ij} &\leftrightarrow [(t/t_i)^{\frac{2}{3}} - 1]\psi_{1|ij} + [(t/t_i)^{-1} - 1]\psi_{2|ij}; \\ {}^H \Pi_{ij} &\leftrightarrow \frac{1}{2}(\chi_{ij}^{\top}(t) - \chi_{ij}^{\top}(t_i)), \end{aligned} \quad (84) $$

together with the correspondence with our initial data fields,

$$ \begin{aligned} G_{ij} &\leftrightarrow \frac{20}{9t_i^2} \psi_1 \delta_{ij} + \chi_{ij}^{\top}(t_i), \\ U_{ij} &\leftrightarrow \frac{2}{3} \frac{1}{t_i} \psi_{1|ij} - \frac{1}{t_i} \psi_{2|ij} + \frac{1}{2} \dot{\chi}_{ij}^{\top}(t_i), \\ W_{ij} &\leftrightarrow \frac{2}{3} \frac{1}{t_i^2} [\psi_{1|ij} + \psi_{2|ij}] + \frac{1}{2} \ddot{\chi}_{ij}^{\top}(t_i) + H_i \dot{\chi}_{ij}^{\top}(t_i). \end{aligned} \quad (85) $$

These three initial fields are not independent. Combining Eq. (16) with Eq. (31) at $t = t_i$, gives a differential equation that links $G_{ij}$ to $U_{ij}$ and $W_{ij}$. We have checked that this equation is consistent with the three relations (85). Taking the trace of the third of Eqs. (85) and using $\Delta\psi_1(\mathbf{x}) + \Delta\psi_2(\mathbf{x}) = -\delta(\mathbf{x}, t_i)$ provide another consistency check. We find $W^k_k \leftrightarrow -(3/2)H_i^2\delta(\mathbf{x}, t_i)$ [which is Eq. (49) in [11]].

### C. Discussion of the difference between standard (Eulerian) and the intrinsic (Lagrangian) approaches

Let us first illustrate the difference between a Lagrangian and an Eulerian approach in a flat-space situation. We take the example of inertial motion that, for a continuum, is governed by the following evolution equations for the three velocity field components $v^a$, indexed by the counter $a = 1, 2, 3$ (a comma denotes a partial derivative with respect to time $t$ or with respect to Eulerian coordinates $x^i$, respectively):

$$ \frac{d}{dt} v^a = v^a,_{t} + \frac{dx^i}{dt} v^a,_{,x^i} = 0; \quad \frac{dx^i}{dt} = v^i. \quad (86) $$

$d/dt$ denotes the total or Lagrangian time-derivative operator. The linearized version of this equation on the Euclidean background $v^a,_{t} = 0$ is solved by a constant velocity field in Eulerian coordinates, $v^a(x^i, t) = v^a(x^i, t_i)$. Introducing Lagrangian coordinates that follow the solution curves, $dX^i/dt = 0$, the total (Lagrangian) time derivative reduces to a partial time derivative at fixed Lagrangian coordinates, so that Eq. (86) reduces to a set of decoupled ordinary and linear differential equations, $v^a,_{t}(X^i, t) = 0$. This allows us to find the general solution of the coupled set of partial and nonlinear differential Eqs. (86) in Eulerian space, as long as the mapping between the Eulerian and the Lagrangian coordinates, provided by the Galilean solution curves $f^a := X^a + v^a(X^i, t_i)(t - t_i)$, is invertible,

$$ v^a(X^i, t) = v^a(X^i, t_i); \quad (87) $$

$$ \Rightarrow v^a(x^i, t) = v^a(X^i = (f^a)^{-1}[x^i, t], t_i). \quad (88) $$

We see in this example that the linearized Eulerian and the linear Lagrangian equations have the same algebraic form, but the solution of the latter includes Eulerian nonlinearities, here the convective term. Thus, in the close vicinity of the initial time, both solutions coincide, but the Lagrangian solution moves away from the Eulerian solution, and this motion is dictated by the solution itself. For the discussion in the next subsection we note that these deviations have an *recursive* structure [5],

$$ v^a(t, x^i) = v^a(t_i, X^i + (t - t_i)v^a(t_i, x^i = X^i + \dots)). \quad (89) $$

This example gives an intuition for the flat-space case. Working on nonflat spaces with Cartan coframes, the local coordinate representation of the coframes $\eta^a(X^i, t)$ exists as a projection on the local cotangent spaces with exact basis $\{dX^i\}$. Nonetheless, the diffeomorphism $f^a$ to the flat-space coordinates $x^i = f^a=_{i}$, as used above, does not exist, since the coframes are not exact forms, $\eta^a \neq d\mathbf{f}^a$. This is the reason why we cannot map the standard solution to an intrinsic solution by a coordinate transformation. Notwithstanding, we can find a solution of the intrinsic equation as we are going to outline in the following subsection.

### D. How to find relativistic Lagrangian first-order solutions for gravitational waves

In order to find explicit first-order Lagrangian solutions of, e.g., the master equation for gravitational waves, Eq. (48), we can exploit the algebraic similarity of the standard linear equation and the Lagrangian master equation. To determine the solution of (48) on a curved manifold we need to prescribe eigenfunctions other than flat-space Fourier modes. This is possible for constant-curvature spaces, but in general situations of an inhomogeneous curvature the solution procedure is not obvious. In the cases where an equation on flat space has algebraically the same
form as an equation on a curved manifold, it is possible to find a solution. In what follows, we outline the procedure of how to find the general solution of the master Eq. (48),

$$H^H \ddot{\Pi}_{ij} + 3H^H \dot{\Pi}_{ij} - a^{-2H} \Pi_{ij}^{|k} = a^{-2} (H^W W_{ij}^{tl} + H_i^H U_{ij}^{tl}).$$

(i) We employ the dictionary, i.e., the last equation of (84), to redefine the dependent variable,

$$H^H \Pi_{ij} \leftrightarrow \frac{1}{2} (\chi_{ij}^{\top}(\mathbf{x}, t) - \chi_{ij}^{\top}(\mathbf{x}, t_i)).$$

Equation (48) is then cast into the standard linear propagation Eq. (79), except that the former is an equation for functions of local (Lagrangian) coordinates, and the latter is an equation for functions of (Eulerian) background coordinates.

(ii) The linear propagation equation on flat space in the standard perturbation theory, Eq. (79), can be solved. The solution in the Fourier space is provided by Eq. (81). From this we have to prescribe initial data and perform the inverse Fourier transformation to obtain the solution as a function of Eulerian coordinates.

(iii) The fact that our master Eq. (48) is algebraically identical (up to the redefinition of the dependent variable) to the Eulerian linear Eq. (79) entitles us to take this Eulerian solution and replace the background coordinates by Lagrangian ones. The resulting solution is transported along geodesics, which are themselves solutions of this equation. They therefore exhibit a nonlinear behavior when seen from the background space reference frame. In practice, we would prescribe the same initial data in the Lagrangian solution, since, initially, the Eulerian and Lagrangian coordinates coincide by definition.

The above outlined solution procedure is not explicitly carried out in examples. We shall return to similar equations in the follow-up paper that includes pressure.

## V. DISCUSSION AND OUTLOOK

As perturbations of space propagate in the perturbed physical spacetime, the equations that describe their dynamics are intrinsically recursive. The standard linear solution for gravitational waves does not capture the nonlinear nature of these perturbations; see the criticisms advanced in [39] that provide arguments why the standard linear solution is at odds with the physical phenomenon, and see also [40]. The common way out is to go to second- and higher-order solutions within the standard approach [38,41,42]. However, as we discussed in the previous section, the intrinsic Lagrangian approach enjoys more generality, even if the equations were linearized. This fact

can be exploited to combine the simplicity of a linear solution with the nonlinear character of the physical problem at hand. Thanks to the present Lagrangian intrinsic approach we are in the situation to uncover nonlinear features of the gravitational waves. In other words, the linear intrinsic solution contains components that would be higher-order contributions in the standard linear approach.

In order to find out the wealth of information revealed by the Lagrangian approach, we now explore some avenues toward a more general solution than the one outlined in the last section.

### A. Extrapolation of the solution into the nonperturbative regime

We can certainly go beyond the linear solution by considering higher-order Lagrangian perturbations, as done in the standard Eulerian approach. However, we can already employ the linear solution supplemented by extrapolation procedures in the spirit of the Newtonian Zel'dovich approximation (see [1] for the discussion of the relativistic generalization of this extrapolation idea).

We now discuss different possibilities to understand nonperturbative properties of gravitational waves in terms of different extrapolation strategies.

#### 1. Functional extrapolation

The first strategy follows the line of thoughts in [1,10]. The idea is to insert the linear Lagrangian solution into exact functionals of the fields of interest [1], and we may then study their average properties [10]. With this strategy we obtain nonperturbative expressions from a linear perturbation. Consider, as an example, the spatial metric. This is a bilinear form of the Cartan coframes, and a linearization of the metric would truncate a quadratic term, thus, leading to a functional expression that is not the metric corresponding to the Cartan coframe that describes the deformation to a given order. The proposal is to keep this quadratic term, since the basic system of equations is written in terms of the coframes only and a consistent perturbation theory already determines the approximation for the coframes to a given order. Another example is the curvature that can be functionally defined in terms of coframes. Looking at the first-order deformation of the coframes, the corresponding curvature is the full functional expression, not the truncated one (see again [10] for an explicit demonstration). (A related paper on averaging the gravitational wave tensor part is [43].)

#### 2. Iterative extrapolation

The second strategy notices that a strict linearization, pursued in this paper, neglects essential physical properties of gravitational waves propagating on a curved background. For example, this linearization replaces the Laplace-de Rham operator, Eq. (58), for the (curvature-generating;
cf. Sec. III B 1) trace-free gravitomagnetic part $^H\Pi^a$ by a simple Laplacian in local coordinates,

$$ (\Delta^{dR} \Pi^a)_i = (\Delta^H \Pi^a)_i + {}^H \Pi^a_k \mathcal{R}^k_i \approx -\Delta_0 {}^H \Pi^a_i^{(1)}, \quad (90) $$

where the first-order Ricci tensor (33) does not contribute at first order here. (At higher orders, the resulting curvature correction would couple the gravitomagnetic part to the trace and the gravitoelectric part of the perturbations.) It is clear that a propagation of gravitational waves in a curved space will in general imply a nontrivial Laplace–de Rham operator. A possible route toward the understanding of these nonlinear properties consists in keeping the Lagrangian time evolution at first order in the linearized master Eq. (48), but considering the full Laplace–de Rham operator of the Cartan coframes in the curved space section, i.e., by replacing the local Laplacian in (48) by the Laplace–de Rham operator. An iterative procedure to solve the resulting equation is suggested: we calculate the perturbed Cartan deformation for all parts according to the linearized scheme presented in this paper, then evaluate the Laplace–de Rham operator for this perturbed deformation (including the coupling to the trace and gravitoelectric parts), and iterate the equation. An intermediate step consists in just considering the rough Laplacian, thus, avoiding the coupling to the trace and the trace-free gravitoelectric solutions due to the curvature correction. This intermediate step already renders the master equation for the perturbations (48) nonlinear in local coordinates due to the presence of the covariant spatial derivatives. Such an iteration may converge to a nonlinear solution after a sufficient number of steps (see, however, [44], and related iteration procedures for the Newtonian problem [45,46]).

## ACKNOWLEDGMENTS

This work is part of a project that, in the final stage, has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement ERC advanced Grant No. 740021-ARTHUS, PI: T. B.). A part of this project was funded by the National Science Centre, Poland, under Grant No. 2014/13/B/ST9/0084. F. A. R. acknowledges support by the École Doctorale PHAST Lyon. F. A. R. and T. B. acknowledge student exchange support from the French-Bavarian Cooperation Center, BFHZ Munich http://www.bayern-france.org/. A. W. acknowledges support from the German research organization DFG, Grants No. WI 4501/1-1 and No. WI 4501/2-1. We thank Alexandre Alles for his helpful participation in the beginning stage of this paper. Particular thanks go to Léo Brunswic, Mauro Carfora, and Pierre Mourier for many valuable remarks and mathematical insights. We thank them and also Martin Kerscher, Jan Ostrowski, Slava Mukhanov, Herbert Wagner, and David Wiltshire for valuable discussions related to this work and comments.

## APPENDIX A: PROPAGATION EQUATIONS FROM THE MAXWELL-WEYL FORMALISM

In this appendix, we will derive Lagrangian propagation equations for the gravitoelectric and gravitomagnetic spatial parts of the Weyl tensor by employing the Maxwell-Weyl formalism [47,48]. Related to this, the reader may also consult [49–52].

### 1. Gravitoelectric and gravitomagnetic parts of the Weyl tensor

The Riemann curvature tensor can be split into a trace part, the Ricci tensor, that is linked to the matter content of the Universe by the Einstein field equations, and a traceless part, the Weyl tensor. This latter tensor, which encodes the relativistic counterpart of the Newtonian tidal forces and the gravitational waves, can be expressed through the Riemann curvature tensor as follows:

$$ C^{\mu\nu}_{\rho\sigma} = {}^{(4)}R^{\mu\nu}_{\rho\sigma} - 2\delta^{[\mu}_{[\rho} {}^{(4)}R^{\nu]}_{\sigma]} + \frac{1}{3}\delta^{[\mu}_{[\rho}\delta^{\nu]}_{\sigma]} {}^{(4)}R. \quad (A1) $$

Let $u^\mu$ be the 4-velocity of the fluid. The Weyl tensor can be split in an irreducible way into gravitoelectric and gravitomagnetic parts, $E_{\mu\nu}$ and $H_{\mu\nu}$, respectively, defined as follows:

$$ E_{\mu\nu} := C_{\mu\nu\sigma} u^\sigma u^\mu; \quad (A2) $$

$$ H_{\mu\nu} := \frac{1}{2} \sqrt{|{}^{(4)}g|} \epsilon_{\alpha\beta\rho(\mu} C^{\alpha\beta)}_{\nu)\sigma} u^\alpha u^\sigma, \quad (A3) $$

where ${}^{(4)}g$ represents the determinant of the 4-metric tensor ${}^{(4)}g$, and $\epsilon_{\rho\mu\alpha\beta}$ is the four-dimensional Levi-Civita symbol. As these parts are symmetric and traceless, they satisfy by construction the following identities:

$$ \begin{aligned} E^\mu_\mu &= 0; & E_{\mu\nu} &= E_{(\mu\nu)}; & E_{\mu\nu} u^\mu &= 0, \\ H^\mu_\mu &= 0; & H_{\mu\nu} &= H_{(\mu\nu)}; & H_{\mu\nu} u^\mu &= 0. \end{aligned} $$

### 2. Propagation equations for gravitoelectric and gravitomagnetic parts

Combining Bianchi identities with the Einstein field equations we obtain

$$ \nabla^{\kappa} C_{\mu\nu\kappa\lambda} = 8\pi G \left( \nabla_{[\mu} T_{\nu]\lambda} + \frac{1}{3} g_{\lambda[\mu} \nabla_{\nu]} T^{\kappa}_{\kappa} \right). \quad (A4) $$

This equation can be turned into an equation for the gravitoelectric and gravitomagnetic parts, as the Weyl tensor can be expressed as follows:
$$C_{\mu\nu\kappa\lambda} = (g_{\mu\nu\alpha\beta}g_{\kappa\lambda\gamma\delta} - \epsilon_{\mu\nu\alpha\beta}\epsilon_{\kappa\lambda\gamma\delta})u^\alpha u^\gamma E^{\beta\delta} \\ + (\epsilon_{\mu\nu\alpha\beta}g_{\kappa\lambda\gamma\delta} + g_{\mu\nu\alpha\beta}\epsilon_{\kappa\lambda\gamma\delta})u^\alpha u^\gamma H^{\beta\delta}, \quad (A5)$$

with $g_{\mu\nu\alpha\beta} \equiv g_{\mu\alpha}g_{\nu\beta} - g_{\mu\beta}g_{\nu\alpha}$, and $\epsilon_{\mu\nu\kappa\lambda} = \sqrt{-{}^{(4)}g} \epsilon_{\mu\nu\kappa\lambda}$.

In the rest frame of the dust fluid without vorticity, Eq. (A4) yields the following constraints and evolution equations, e.g., [53]:

$$E^k_{i\parallel k} - g_{ik}\epsilon^{kmn}\Theta_{ml}H^l_n = \frac{8\pi G q_{i|j}}{3J}; \quad (A6)$$

$$H^k_{i\parallel k} + g_{ik}\epsilon^{kmn}\Theta_{ml}E^l_n = 0; \quad (A7)$$

$$\dot{E}_{ij} + 2\Theta E_{ij} - 3\Theta_{k(i}E^k_{j)} - \Theta^k{}_l E^l{}_k g_{ij} \\ = g_{m(i} \epsilon^{mkl} H_{j)l\parallel k} - 4\pi G \frac{q_i}{J} \left( \Theta_{ij} - \frac{1}{3} \Theta g_{ij} \right), \quad (A8)$$

$$\dot{H}_{ij} + 2\Theta H_{ij} - 3\Theta_{k(i}H^k_{j)} - \Theta^k{}_l H^l{}_k g_{ij} \\ = -g_{m(i} \epsilon^{mkl} E_{j)l\parallel k}. \quad (A9)$$

Note that, in order to simplify (A9), we used the identity $q_{i|k}\epsilon^{mkl}g_{m(i}g_{j)l} = 0$, which stems from the antisymmetry of the Levi-Civita symbol.

The expansion tensor can be split into its kinematical parts (with the shear $\sigma^i_j$ and the rate of expansion $\Theta$, for vanishing vorticity),

$$\Theta^i_j = \sigma^i_j + \frac{1}{3}\Theta\delta^i_j. \quad (A10)$$

The initial density field $q_i$ can be split into $(q_H)_i + \delta q_i$. At first order, the perturbation of the background density $\delta q_i$ can be expressed through the trace of the nonintegrable generalization of the Newtonian acceleration gradient $W^a$; cf. Eq. (14), the generalization of the Poisson equation, $W = -4\pi G \delta q_i$.

As $E_{ij}$ and $H_{ij}$ are symmetric and traceless tensors, at first order the constraints and evolution equations reduce to the following set:

$$E^k_{i|k} = -\frac{2}{3a^3} W_{|i}; \quad (A11)$$

$$H^k_{i|k} = 0; \quad (A12)$$

$$\dot{E}_{ij} + 3HE_{ij} - \frac{1}{a} \epsilon_{(i}^{kl} H_{j)l|k} = -4\pi G q_H \sigma_{ij}; \quad (A13)$$

$$\dot{H}_{ij} + 3HH_{ij} + \frac{1}{a} \epsilon_{(i}^{kl} E_{j)l|k} = 0, \quad (A14)$$

where $q_H = (q_H)_i/a^3$ denotes the homogeneous part of the density field. In order to obtain a propagation equation for the gravitoelectric spatial part of the Weyl tensor $E_{ij}$, we first multiply Eq. (A13) by $a^4(t)$ and perform a time derivative. Then, we inject Eq. (A14) and simplify the product of the Levi-Civita symbols using the following identity:

$$\epsilon_{kl(i} \epsilon_{j)m} = \delta_{ij}\delta_{km}\delta_{ln} - \delta_{ij}\delta_{lm}\delta_{kn} - \delta_{k(i}\delta_{j)m}\delta_{ln} + \delta_{k(i}\delta_{j)n}\delta_{lm} \\ + \delta_{l(i}\delta_{j)m}\delta_{kn} - \delta_{l(i}\delta_{j)n}\delta_{km}. \quad (A15)$$

Finally, using Eq. (A11), we obtain the following second-order propagation equation:

$$\Box_0 E_{ij} - 7H\dot{E}_{ij} - 4(4\pi G q_H + \Lambda)E_{ij} \\ = -\frac{1}{a^5} D_{ij} W + 4\pi G H q_H \sigma_{ij}, \quad (A16)$$

where $\Box_0 X_{ij} := -\ddot{X}_{ij} + \frac{1}{a^2}\Delta_0 X_{ij}$ denotes the Lagrangian d'Alembertian applied to the tensor field $X_{ij}$. To obtain the previous equation, we used the Friedmann equations to replace $3H^2 + \ddot{a}/a$ by $20\pi G q_H/3 + 4\Lambda/3$ (at order 0 the curvature constant $k$ is assumed null). Furthermore, we used the compact notation $D_{ij} W$ for $W_{|ij} - 1/3\delta_{ij}\Delta_0 W$.

The same logic applied to $H_{ij}$ leads to the following second-order propagation equation:

$$\Box_0 H_{ij} - 7H\dot{H}_{ij} - 4(5\pi G q_H + \Lambda)H_{ij} \\ = -\frac{1}{a} 4\pi G q_H \epsilon_{(i}^{kl} \sigma_{j)l|k}. \quad (A17)$$

[The coefficients in front of the third terms in Eqs. (A16) and (A17) differ, since the first equation features a term $\propto \dot{\sigma}_{ij}$ that can be replaced through $\dot{\sigma}_{ij} = -E_{ij}$; cf. Eq. (111) in [1], which changes the coefficient $5\pi G q_H + \Lambda$ to $4\pi G q_H + \Lambda$.]

## APPENDIX B: ERRATUM TO [11]

In view of the definition $J \equiv \frac{\sqrt{g}}{\sqrt{G}}$ employed in this paper as well as in [10,11], the momentum constraint equations (27) of [11] must contain factors of $\frac{1}{J}$, written correctly in Eq. (6) (in [1] the definition $J \equiv \sqrt{g}$ was used so that the inclusion of the factors is not needed).
[1] 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).

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

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

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

[5] 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 (IOP Press, Amsterdam, 1996).

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

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

[8] 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).

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

[10] 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).

[11] 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).

[12] B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. **119**, 161101 (2017).

[13] C. D. Hill and P. Nurowski, How the green light was given for gravitational wave search, arXiv:1608.08673.

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

[15] 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).

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

[17] 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).

[18] 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).

[19] 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.

[20] F. Al Roumi and T. Buchert, Relativistic structure formation models and gravitoelectromagnetism, in *Proceedings of the Fourteenth Marcel Grossmann Meeting on General Relativity, Rome July 2015*, edited by M. Bianchi, R. T. Jantzen, and R. Ruffini (World Scientific, Singapore, 2017), p. 2345.

[21] W. V. D. Hodge, *The Theory and Applications of Harmonic Integrals* (Cambridge University Press, Cambridge, UK, 1941).

[22] T. Aubin, *Some Nonlinear Problems in Riemannian Geometry* (Springer, New York, 1998).

[23] T. Frankel, *The Geometry of Physics: An Introduction* (Cambridge University Press, Cambridge, UK, 2011).

[24] J. M. Lee, Introduction to Smooth Manifolds, *Graduate Texts in Mathematics* 2nd ed. Vol. 218 (Springer, New York, 2013).

[25] M. Nakahara, *Geometry, Topology and Physics* (Taylor and Francis, Boca Raton, FL, 2003).

[26] W. P. Thurston, *The Geometry and Topology of Three-Manifolds* (Princeton University Press, Princeton, NJ, 1997).

[27] M. Lachieze-Rey and J. P. Luminet, Cosmic topology, Phys. Rep. **254**, 135 (1995).

[28] X. Roy and T. Buchert, Relativistic cosmological perturbation scheme on a general background: Scalar perturbations for irrotational dust, Classical Quantum Gravity **29**, 115004 (2012).

[29] G. Perelman, The entropy formula for the Ricci flow and its geometric applications, arXiv:math/0211159.

[30] G. E. Bredon, *Topology and geometry* (Springer, New York, 1993), Vol. 139.

[31] G. De Rham, Sur la théorie des formes différentielles harmoniques, Ann. Univ. Grenoble **22**, 135 (1946).

[32] W. V. Hodge, *The theory and applications of harmonic integrals*, (CUP Archive, 1989).

[33] C. Voisin, *Hodge Theory and Complex Algebraic Geometry II* (Cambridge University Press, Cambridge, UK, 2003), Vol. 2.

[34] E. Gausmann, R. Lehoucq, J.-P. Luminet, J-P. Uzan, and J. Weeks, Topological lensing in spherical spaces, Classical Quantum Gravity **18**, 5155 (2001).

[35] M. White, D. Scott, and E. Pierpaoli, Boomerang returns unexpectedly, Astrophys. J. **545**, 1 (2000).

[36] R. Aurich and F. Steiner, The cosmic microwave background for a nearly flat compact hyperbolic universe, Mon. Not. R. Astron. Soc. **323**, 1016 (2001).

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

[38] S. Matarrese, S. Mollerach, and M. Bruni, Relativistic second-order perturbations of the Einstein-de Sitter universe, Phys. Rev. D **58**, 043504 (1998).

[39] J. G. Pereira, Gravitational waves: A foundational review, arXiv:1305.0777.

[40] A. Ashtekar and B. Bonga, On the ambiguity in the notion of transverse traceless modes of gravitational waves, Gen. Relativ. Gravit. **49**, 122 (2017).
[41] S. Andrianomena, C. Clarkson, P. Patel, O. Umeh, and J.-P. Uzan, Non-linear relativistic contributions to the cosmological weak-lensing convergence, *J. Cosmol. Astropart. Phys.* **06** (2014) 023.

[42] H. Arcos, M. Krssak, and J. G. Pereira, Exploring higher-order gravitational waves, arXiv:1504.07817.

[43] I. A. Brown, L. Schrempp, and A. Koshore, Accelerating the universe with gravitational waves, *J. Cosmol. Astropart. Phys.* **02** (2010) 011.

[44] D. Christodoulou and B. G. Schmidt, Convergent and asymptotic iteration methods in general relativity, *Commun. Math. Phys.* **68**, 275 (1979).

[45] T. Buchert, An exact Lagrangian integral for the Newtonian gravitational field strength, *Phys. Lett. A* **354**, 8 (2006).

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

[47] G. F. R. Ellis and H. van Elst, Cosmological models (Cargèse lectures 1998), in *Theoretical and Observational Cosmology: Proceedings of the NATO Advanced Study

Institute on Theoretical and Observational Cosmology, Cargèse, France, 1998, edited by M. Lachièze-Rey, NATO Science Series C, Mathematical and Physical Sciences Vol. 541 (Kluwer Academic, Boston, 1999), pp. 1–116.

[48] J. Wainwright and G. F. R. Ellis, *Dynamical Systems in Cosmology* (Cambridge University Press, Cambridge, UK, 1997).

[49] P. Dunsby, B. Bassett, and G. F. R. Ellis, Covariant analysis of gravitational waves in a cosmological context, *Classical Quantum Gravity* **14**, 1215 (1997).

[50] J. Ehlers and T. Buchert, On the Newtonian limit of the Weyl tensor, *Gen. Relativ. Gravit.* **41**, 2153 (2009).

[51] R. Maartens and B. A. Bassett, Gravito-electromagnetism, *Classical Quantum Gravity* **15**, 705 (1998).

[52] R. Maartens, G. F. R. Ellis, and S. Siklos, Local freedom in the gravitational field, *Classical Quantum Gravity* **14**, 1927 (1997).

[53] E. Bertschinger and A. J. S. Hamilton, Lagrangian evolution of the Weyl tensor, *Astrophys. J.* **435**, 1 (1994).