# On general-relativistic Lagrangian perturbation theory and its non-perturbative generalization

Thomas Buchert, Ismael Delgado Gaspar, Jan Jakub Ostrowski

► To cite this version:

Thomas Buchert, Ismael Delgado Gaspar, Jan Jakub Ostrowski. On general-relativistic Lagrangian perturbation theory and its non-perturbative generalization. Universe, 2022, 8 (11), pp.583. 10.3390/universe8110583. [https://doi.org/10.3390/universe8110583](https://doi.org/10.3390/universe8110583)

HAL Id: hal-03815153

https://hal.science/hal-03815153

Submitted on 19 Mar 2023

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>
<figure><img src="image_2.png" /></figure>

Review

# On General-Relativistic Lagrangian Perturbation Theory and Its Non-Perturbative Generalization

Thomas Buchert <fnref n="1" />,*,† Ⓓ, Ismael Delgado Gaspar <fnref n="1" />,<fnref n="2" />,† Ⓔ and Jan Jakub Ostrowski <fnref n="2" /> Ⓔ

**Abstract:** The Newtonian Lagrangian perturbation theory is a widely used framework to study structure formation in cosmology in the nonlinear regime. We review a general-relativistic formulation of such a perturbation approach, emphasizing results on an already developed extensive formalism including among other aspects: the non-perturbative modeling of Ricci and Weyl curvatures, gravitational waves, and pressure-supported fluids. We discuss subcases of exact solutions related to Szekeres Class II and, as an exact average model, Ricci-flat LTB models. The latter forms the basis of a generalization that we then propose in terms of a scheme that goes beyond the relativistic Lagrangian perturbation theory on a global homogeneous-isotropic background cosmology. This new approximation does not involve a homogeneous reference background and it contains Szekeres class I (and thus general LTB models) as exact subcases. Most importantly, this new approximation allows for the interaction of structure with an evolving "background cosmology", conceived as a spatial average model, and thus includes cosmological backreaction.

**Keywords:** relativistic cosmology; Lagrangian perturbation theory; exact solutions; cosmological backreaction

## 1. Introduction

Relativistic cosmological perturbation theory gave rise to many exciting results regarding the early Universe, CMB (Cosmic Microwave Background) physics, and the late Universe for the formation of large-scale structure. General relativity, unlike other physical theories describing the weak and strong nuclear forces and electrodynamical interactions, does not rely on the notion of *background*, understood as a pre-fixed absolute spacetime. This fact together with the nonlinearity of general relativity introduces conceptual difficulties in the formulation of a cosmological perturbation theory. Difficulties that can be at least partially circumvented in the astrophysical context by a suitable choice of boundary conditions, e.g., asymptotic flatness, but not when we aim at describing the Universe as a whole.

Whether we consider the inflationary epoch, scalar fields, and semi-classical quantum gravity (see, e.g., [1,2]) or the matter-dominated, fully classical Universe, we can usually distinguish two main approaches to relativistic cosmological perturbations, based on which quantity is being perturbed. In both cases, we have to put a reference background by hand and the commonly adopted choice is a member of the family of FLRW (Friedmann-Lemaître-Robertson-Walker) spacetimes. The original work by Lifshitz [3] (see also [4–6], and follow-up papers of this author) introduced perturbations of the FLRW metric coupled to the perturbations of the stress-energy tensor via the linearized Einstein equations. This work employs intrinsic comoving coordinates. The issue of coordinate choices later highlighted the problem of gauge fixing, i.e., choosing the way of mapping the physical

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

check for
updates

**Citation:** Buchert, T.; Delgado Gaspar, I.; Ostrowski, J.J. On General-Relativistic Lagrangian Perturbation Theory and Its Non-Perturbative Generalization. *Universe* **2022**, *8*, 583. https://doi.org/10.3390/universe8110583

Academic Editors: Luigi Tedesco and Giovanni Marozzi

Received: 27 September 2022

Accepted: 27 October 2022

Published: 4 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

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

<footnote n="2">² Department of Fundamental Research, National Centre for Nuclear Research, Pasteura 7, 02-093 Warsaw, Poland</footnote>

<footnote>* Correspondence: buchert@ens-lyon.fr</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
spacetime to the assumed background, which can be thought of as a mapping between two distinct manifolds in the framework of *gauge-invariant perturbation theory* [7–10]. It is worth noticing that the gauge fixing allows to pullback the physical quantities from the physical spacetime, resulting in, e.g., curvature perturbations residing as fields on the background manifold rather than intrinsic geometrical features of the physical manifold. Such a setup is in the spirit of the Eulerian picture in hydrodynamics where we examine the fluid in specific, fixed locations as opposed to the Lagrangian picture where we follow the fluid flow intrinsically [11]. In this context, the background spacetime defines an arbitrary, absolute reference frame. In the seminal paper [12], gauge-invariant variables where defined. However, although the basic physical perturbation variable is the density contrast, none of Bardeen’s variables actually represented this quantity. This issue was partially resolved by John Stewart [13]. A different (although fully compatible) approach was presented in [14], where 1 + 3 threading of spacetime was utilized. It can be considered as a step towards the Lagrangian approach to relativistic cosmological perturbation theory despite the fact that the mapping between two manifolds was still in place. However, in the above-mentioned work, the authors defined the variables related to the congruence of the world-lines of the fluid (spatial projections of density gradient and extrinsic curvature) and the mapping exists between the fluid trajectories not between fixed points in the two manifolds.

While gauge-invariant variables have to be constructed afresh at any order of the perturbations, a genuine general-relativistic treatment does not suffer from this issue when we consider trace-free tensors that vanish on the chosen reference background [13]. We shall be guided by covariant formulations that do not rely on the choice of coordinates, but on the choice of a foliation of spacetime [15,16]. Gauge choices in the two-manifold approaches can be mapped to foliation choices in a covariant setting, see [17] and references therein (e.g., the series of papers by Mario Novello and collaborators). The Lagrangian coordinates, being intrinsic to the fluid congruence, allow us to keep covariance on spatial leaves of a 3 + 1 foliation, while the existence of hypersurfaces allows us to define a background evolution for scalar variables that are conceived as spatial averages of geometric and matter variables [18].

Before moving to a more detailed description of Lagrangian cosmological perturbation
theory (Newtonian and relativistic), we wish to recall the fundamental problem signalized
at the beginning of this introduction: general relativity does neither rely on the notion of
background, nor does it give us a priori any prescription on how to construct one [19]. One
of the aspects of this ever persistent conundrum is the fact that the FLRW background
metric and the spatially homogeneous density obtained by any averaging or smoothing
procedure involving limits (e.g., integrating or taking the limit of some parameter related
to scale or inhomogeneities) are, by construction, not related to each other through the
Einstein equations. The first part of the problem becomes apparent when we try to define
the average of a tensor, which without any additional assumptions, is an ill-defined concept.
Even if we had defined such an operation, it would not commute with taking the inverse
of the metric and calculating the curvature (due to the quadratic terms), and this would
necessarily lead to an additional term in the Einstein equations, the backreaction term. We
note that already volume-averaging and time-evolution are non-commuting operations,
and their difference consists of positive-definite backreaction terms even before invoking
dynamics or a theory of gravity [20,21]. This also holds for operations such as Ricci-
flow smoothing and scaling in a 3-manifold [22] and for other coarse-graining techniques
such as renormalization [23]. The perhaps simplest and to date widely used approach
to the backreaction problem consists in scalar averaging of Einstein’s equations, within
frameworks of 3 + 1 foliations [16,24,25], 2 + 1 + 1 foliations on the lightcone [26], and
their manifestly 4-covariant formulations [27,28], see also [15,29]. Although essential to
the relativistic perturbation theory, the topic of backreaction lies beyond the scope of this
review except for the part where we describe how a concrete model for inhomogeneities
such as the Lagrangian perturbation theory can serve as a closure condition to the otherwise
under-determined system of averaged equations. It was also not our ambition to give an
exhaustive account of all of the many developments in cosmological perturbation theory in general. We rather focus on developments in relativistic Lagrangian perturbation theory, starting however with a brief introduction into Newtonian Lagrangian perturbation theory.

## 1.1. *A Brief History of Lagrangian Perturbation Theory for Newtonian Flows*

The context that led Yakob Zel'dovich to his proposal of extrapolating the Eulerian linear theory of gravitational instability further into the nonlinear regime may be traced back to a book project he conducted together with Anatoly Myshkis, in which inertial motion and the formation of structure was discussed in affectionate spadework [30]. In addition, Zel'dovich might have been aware of related work in plasma physics such as Dawson's discussion of multi-streaming in a cold plasma [31]. In view of this, he must have seen the superiority of a Lagrangian description of fluid motion over the established Eulerian perturbation theory for the purpose of understanding and describing structure formation in the Universe. Upon linearizing the exact solution for Galileian continua, it is evident that a time-dependent rescaling of the dependent and independent variables suffices to match both the Eulerian linear perturbation solution on a Friedmannian background cosmology (see, e.g., [32]) and the powerful nonlinear solution for inertial systems. His early papers [33–36], already explore in detail the consequences of this idea, which— contrary to popular opinions at the time (nonlinear evolution was understood, both analytically and numerically, through spherically symmetric systems)—predicts a highly anisotropic process of structure formation into high-density objects. He called the first collapsing objects *bliny*, the Russian word for small *pancakes*. He had in mind the local picture which was later “americanized” to extended Megaparsec-*pancakes*.

Soon thereafter, confirmations of this idea in the context of the Euler–Newton system for self-gravitating dust (i.e., pressure-less matter) were advanced within Zel'dovich's group. The paper [37] (among others, see, e.g., the considerations on the density probability distribution in [38–40], investigations on self-similar solutions [41], and trajectories [42]) stands out as a convincing analytical justification. The heart of this work is a test of self-consistency of Zel'dovich's approximation, henceforth denoted ZA, conducted through comparison of the density from (i) the field equations and (ii) from the exact density integral (this test is explicitly discussed in [43]). As a side-result, this test also shows that Zel'dovich's approximation provides a special exact solution for plane-symmetric motions on a Friedmannian background. In parallel, Sergei Shandarin performed a numerical simulation, at low resolution at the time (see references in [44]), but confirming the above picture, and reinforcing what was earlier found by Lin, Mestel, and Shu [45] (see also [46,47]).

Further numerical work within Zel'dovich's group [48,49], and especially Sergei Shandarin's move to Lawrence, Kansas, initiated detailed numerical investigations of the *pancake picture* together with Adrian Melott, who developed his expertise in the early 1980's [50]. The "Eastern" *pancake picture* was by then popular for massive neutrino-dominated (*top-down structure formation*) cosmogonies, in which structure on large scales forms first and later fragments into smaller units. It was an equal competitor to the "Western" (*bottom-up structure formation*) hierarchical cosmogony, in which large-scale structure forms out of small-scale inhomogeneities in the course of a merger process. The nowadays still accepted ΛCDM standard model of cosmology performs in between these two scenarios.

Returning to the history of Zel'dovich's analytical model, the collaboration with Vladimir Arnol'd was probably the final highlight for inhomogeneous dust matter models [51], in which the *Lagrange-singularity theory* [52] was employed to furnish a classification scheme of the building blocks of large-scale structure in the Universe. It came right after two further analytical pieces that were missing from the picture: the form of Zel'dovich's approximation for the modeling of structure formation on Friedmannian background cosmologies with curvature [53], and the *derivation* of exact plane-symmetric solutions for self-gravitating dust flows [54], confirming what was suggested by [37]¹. The effect of velocity dispersion was discussed by Shukurov [56]. Thereafter, an extension of Zel'dovich's approximation to phenomenologically access the multi-stream regime was proposed in [57] through
application of the known solution to the 3D Burgers equation: a Laplacian forcing was added to the scaled inertial system. A derivation of this *adhesion approximation* from kinetic theory was provided in [58]; see the review [59] that also investigates deviations from mean field gravity². A further extension was suggested in [65] by including a noise term to lead to KPZ-type equations. For a summary of the different approximations in a cosmological context, including non-perturbative approaches and exact integrals of the field equations, see [66].

Most of these early achievements and a wealth of application work of Zel'dovich's approximation were done before a *Lagrangian theory for self-gravitating flows* and, based on this, a *Lagrangian perturbation approach* was even formulated; the basic mathematical framework (apart from the self-consistency test in [37] and the plane-symmetric investigation in [54]) was still that of inertial motion. Final confirmation and derivation of Zel'dovich's approximation were based on the transformation of the full Euler-Newton system to its Lagrangian form [43,67,68], and, finally, the systematic framework of a Lagrangian perturbation theory, in which Zel'dovich's approximation forms a subclass of the first-order solutions [43,69]; closed form solutions at second [70], third [71], and fourth-order [72] can be reconstructed using the n-th-order perturbation and solution schemes investigated in [73] (see also more recent work on recursion relations [74,75] and related investigations of time-analyticity or smoothness of trajectories [76,77]). The early 1990's can be considered as the birth of Lagrangian perturbation theory, when other groups started to contribute [78-80] (see the early review papers [44,73,81], a brief summary [82], the tutorial [83], the more general review [84], and the more recent review [85]). Most of these investigations employ the matter model dust; velocity dispersion can be treated as an anisotropic pressure that adds a Laplacian forcing as in the adhesion approximation [57,58]. In a Lagrangian perturbation approach, this Laplacian forcing must be transformed to Lagrangian coordinates, which introduces substantial nonlinearity to the problem [86].

We remark that current simulations of structure formation are initialized by the Zel'dovich approximation or its higher-order extensions. What is not well-known is that this approximation and higher-order extensions also perform well in comparison with numerical simulations until the present epoch of structure formation, if the initial spectrum is appropriately truncated to avoid overproduction of structure after shell-crossing [87-91]; for the adhesion approximation [57], see [92]; for recent direct comparisons between higher-order Lagrangian perturbation solutions and numerical simulations, see [93,94]. Shell-crossing itself was also in the focus of interest [70,95,96], a sustained interest up to the present day, see, e.g., the recent papers [97-103]. Despite shell-crossing, Lagrangian schemes are capable of quickly simulating a large number of realizations of Hubble volumes, a task that is still executed with expensive numerical simulations. This remark is especially relevant since the current resolution power of numerical simulations still falls short to access small-scale structure within Hubble volumes, a scale that roughly corresponds to scales where also the analytical schemes fail. The simulation of mock catalogs is an exception, see, e.g., [104-109], and references therein. The corresponding general-relativistic schemes that we shall introduce would in addition allow to map the curvature distribution in mock catalogs, and from the metric, we will be able to implement concrete ray-tracing methods along the lightcone for a precise estimation of lensing effects and for realistic estimates of distances.

Analytical schemes allow for extremely high-resolution studies. After the first study in two spatial dimensions [43], further insights were obtained for three-dimensional special and generic initial data [95,110] including hints for a mechanism of gravitational fragmentation at higher orders in Lagrangian perturbation theory³. Detailed studies of 3D pancakes have been conducted, but did not end up in publications at the time, an example is shown in Figure 1⁴.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** An illustration of the complexity of caustic formation that already arises in a three-dimensional inertial system. Shown are four views of the same evolution stage. Yellow surfaces correspond to caustic surfaces from the degeneracy of the first eigenvalue of the deformation tensor; red surfaces correspond to secondary degeneracies of the second eigenvalue; third degeneracies are barely visible and may occur later. A hierarchy of further nested caustics forms by including gravity. Green dots correspond to *cusp singularities*, while the blue dots correspond to *higher-order singularities*. Pink-grey shaded areas correspond to density contours. Figure by Arno Weiss [112].</figcaption></figure>

During the following decades, many papers referring to Lagrangian approaches have been published, also guided by improvements of numerical technologies. These works initiated a revival of Lagrangian perturbation theory at the beginning of this millennium: the inverse Lagrangian perturbation solutions have been employed for reconstruction techniques, see [113] and references therein to earlier work, Ref. [114] and follow-up papers of these authors, as well as related techniques [115]. Work is pursued, e.g., for convergence issues [99,116], re-expansions [117], renormalization, re-summation, and recursion techniques [74,75,118,119], respectively, redshift space distortions [120,121] and references therein, and statistical biasing relations between the galaxy and matter distributions [122] (a certainly incomplete list).

## 1.2. Relativistic Generalization of Lagrangian Perturbation Theory

For the introduction of a dictionary of Newtonian variables and equations to general relativity, we first recall the governing equations in Newtonian theory and restrict these considerations to a self-gravitating dust continuum (i.e., a pressure-less fluid).

### 1.2.1. Lagrange–Newton System

The description of fluid motion in Newtonian cosmology and, in particular, the problem of perturbations of the Hubble flow have been largely studied in an *Eulerian* coordinate system that is comoving with the Hubble flow, $\mathbf{x} = a(t)\mathbf{q}$, with $a(t)$ a solution in the family of FLRW cosmologies; $\mathbf{x}$ are inertial non-rotating coordinates in Euclidean space. On the contrary, the Lagrangian coordinate system, denoted by $\mathbf{X}$, is a curvilinear, possibly rotating frame in Euclidean space which is defined such as to move with the fluid elements. The Lagrangian description enjoys a number of advantages over the Eulerian one, some of which we shall highlight. This remains true within a perturbation approach. For an early nice introduction into the Lagrangian framework, see [123], for a systematic introduction for self-gravitating Newtonian flows, see [73]. For many complementary investigations of relativistic Lagrangian perturbations see, e.g., [124–129], including gradient
expansions [130,131] and follow-up papers by these authors, as well as the recent papers on gradient expansions [132,133].

Summarizing the Lagrange–Newton system briefly, we note that the one-parameter family of diffeomorphisms, $X \mapsto x = f(X, t)$ allows to reduce Eulerian fields to definitions or provide exact integrals for them, e.g., the velocity field $v$, the acceleration field (equivalent to the Newtonian field strength $g$ for a dust matter model), the density $\rho$, and the vorticity $\omega = \nabla \times v$ can be derived from the knowledge of the family of trajectories $f$:

$$x := f(X, t); \quad v := \dot{f}(X, t); \quad a := \ddot{f}(X, t); \quad \rho(X, t) = \frac{\dot{\rho}}{J(X, t)}; \quad \omega = \frac{\dot{\omega} \cdot \nabla_0 f}{J(X, t)}; \quad \text{etc.,} \quad (1)$$

where $\dot{\rho}$ and $\dot{\omega}$ stand for the initial values of the density and vorticity fields, $\nabla_0$ for the nabla operator with respect to the Lagrangian coordinates, and $J$ is the determinant of the Jacobian matrix ($f^a$),

$$J := \det(f^a)_{(X, t)} = \frac{1}{6} \epsilon_{ijk} \epsilon^{klm} f^i_{|k} f^j_{|l} f^k_{|m}, \quad (2)$$

where we have denoted the spatial derivative with respect to $X_i$ by a vertical slash $|_i$, the partial time derivative by an overdot, $\epsilon_{ijk}$ being the Levi–Cività symbol.⁵ The remaining field equations of the Euler–Newton system, $\nabla \times g = 0; \nabla \cdot g = \Lambda - 4\pi G\rho$, become evolution equations in the Lagrangian picture, comprising the closed Lagrange–Newton system:

$$\delta_{ab} \dot{f}^a_{|i} f^b_{|j} = 0; \quad \frac{1}{2} \epsilon_{abc} \epsilon^{ikl} \dot{f}^a_{|i} f^b_{|k} f^c_{|l} = \Lambda J - 4\pi G \dot{\rho}. \quad (3)$$

We notice that this is a system of equations for the gradient $df^a$, i.e., involving matrices for the derivatives of $f$ that can be easily generalized to tensors as we shall see below.

A Lagrangian perturbation theory results from the ansatz that represents $df^a$ in terms of a homogeneous-isotropic deformation and expands the Lagrange–Newton system to a given order in the deviation gradient field $dP: df = a(t)(X + dP)$.

### 1.2.2. Lagrange–Einstein System

Moving to general relativity, we have to give up the global vector space that is needed to construct vector perturbations and their gradients. Instead, this global vector space reduces to the tangent space at a point on the Riemannian manifold. Thinking of Riemannian space sections, we introduce a spatial metric form $g$ with coefficients $g_{ij}$ in an exact (co-tangential) basis $dX^i \otimes dX^j$. We can write any such metric as a quadratic form of deformation one-forms or co-frame fields, $g = \delta_{ab} \eta^a \otimes \eta^b$, i.e., in terms of coefficients, $g_{ij} = \delta_{ab} \eta^a_i \eta^b_j$. Such a metric form is Ricci-flat, i.e., its Ricci tensor vanishes everywhere, if there exist functions $f^a$, such that the one-forms can be written as exact forms,

$$\eta^a \mapsto df^a. \quad (4)$$

We call (4) the Euclidean restriction. It can be profitably used to either extrapolate Newtonian fields written in Lagrangian form into general-relativistic expressions through inversion of (4), or it transforms any general-relativistic (here spatial) expression into its Newtonian counterpart. This provides a straightforward dictionary that, e.g., in a flow-orthogonal foliation of spacetime of a dust matter model, does not need to send the speed of light to infinity as in standard Newtonian limits [134]. We note that the Eulerian form of the Newtonian equations, written in an inertial coordinate system, depends on these coordinates; to cure this coordinate-dependence, we can transform to rotating coordinates and collect the non-inertial forces in connection coefficients (the Newton–Cartan equations, e.g., [135] and references therein). A connection also arises through the restriction of the spatial relativistic Levi–Cività connection using (4), leading to an integrable Newtonian connection in a Lagrangian frame [55]. For any metric, if we can find a
coordinate transformation $x^i = f^{a \equiv i}(X^j, t)$ that transforms the Euclidean metric coefficients, $\delta_{ij}dx^i dx^j = \delta_{ab}f^a_{\ |i}f^b_{\ |j}dX^i dX^j$, into the metric coefficients $g_{ij} = \delta_{ab}f^a_{\ |i}f^b_{\ |j}$, then these latter are just a rewriting of the flat spatial metric. Given this remark, any perturbation theory that features metric forms of the integrable form does not describe relativistic inhomogeneities; metric coefficients of the above form describe Newtonian (Lagrangian) perturbations on a flat background space. An intrinsic general-relativistic perturbation theory deforms the background geometry; in other words, the perturbations live in a perturbed space, not on a reference background. This remark also shows that perturbation terms can be identified through their integrability in the spirit of (4) to be coordinate artifacts. We shall again encounter these thoughts when we discuss gravitational wave solutions.

The foundational paper I, [55], and the subsequent series of papers, transforms the following Einstein equations for the case of irrotational dust in a 3 + 1 representation, i.e., the 3 irrotationality conditions on the expansion tensor $\Theta_{ij}$ with its trace denoted by $\theta$, the 6 evolution equations for the 6 components of the symmetric spatial metric $g_{ij}$ with the 4 energy and momentum constraint equations:

$$ \Theta_{[ij]} = 0; \qquad (5a) $$

$$ \dot{g}_{ij} = 2g_{ik}\Theta_j^k; \quad \dot{\Theta}_j^i + \theta\Theta_j^i = -\mathcal{R}_j^i + (4\pi G\rho + \Lambda)\delta_j^i; \qquad (5b) $$

$$ \theta^2 - \Theta_j^i \Theta_i^j = -\mathcal{R} + (16\pi G\rho + 2\Lambda); \qquad (5c) $$

$$ \Theta_{j||i}^i - \theta_{|j}^i = 0, \qquad (5d) $$

where $\mathcal{R}_{ij}$ denote the coefficients of the spatial Ricci tensor, with $\mathcal{R}$ its trace; 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-Civita; the dust density is given by $\rho = \partial\rho/\partial J^{-1}$, $J \ge 0$; $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))$.

The system (5) is transformed into a set of 13 equations for the three spatial co-frame fields, i.e., 9 evolution equations for the 9 components of the spatial co-frames and 4 equations stemming from the 4 constraint equations in the classical metric representation (5). We call the resulting system the *Lagrange-Einstein system*:⁶

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

$$ \frac{1}{2J} \epsilon_{abc} \epsilon^{ikl} (\dot{\eta}_j^a \eta_k^b \eta_l^c) = -\mathcal{R}_j^i + (4\pi G\rho + \Lambda)\delta_j^i; \qquad (6b) $$

$$ \frac{1}{J} \epsilon_{abc} \epsilon^{mjk} \dot{\eta}_j^a \eta_k^b \eta_l^c = -\mathcal{R} + (16\pi G\rho + 2\Lambda); \qquad (6c) $$

$$ \left( \frac{1}{J} \epsilon_{abc} \epsilon^{ikl} \dot{\eta}_j^a \eta_k^b \eta_l^c \right)_{||i} = \left( \frac{1}{J} \epsilon_{abc} \epsilon^{ikl} \dot{\eta}_i^a \eta_k^b \eta_l^c \right)_{|j}. \qquad (6d) $$

In this system, we left the Ricci tensor and its trace as well as the covariant derivative implicit, since writing them out in terms of the co-frames becomes non-transparent for our discussion.

We can recover the previously presented Lagrange-Newton system through elimination of the Ricci curvature. We keep the first 3 equations (6a) that also hold for the second time-derivative of the co-frame fields, and Raychaudhuri's equation (derived from the trace of the equation of motion (6b) combined with the energy constraint (6c)):

$$ G_{ab} \dot{\eta}_{[i}^a \eta_{j]}^b = 0; \quad \frac{1}{2J} \epsilon_{abc} \epsilon^{ikl} \dot{\eta}_i^a \eta_k^b \eta^c_l = \Lambda - 4\pi G\rho. \qquad (7) $$
The dictionary (4) works and reproduces what we call the gravito-electric part of the Lagrange–Einstein system consisting of 4 equations. Its Euclidean restriction is a closed system (the Lagrange–Newton system (3)) for the three vector components $f^a$ of the Lagrangian deformation. The 9 components of co-frame fields are solutions of this system complemented by the gravito-magnetic part of the Lagrange–Einstein system. Perturbation solutions are tensorial and feature a richer structure within the full system. We can trace this back to the fact that after elimination of the trace of the Ricci tensor, we are left with Equation (6b) for the trace-free part of the Ricci tensor. An evolution equation for trace-free parts is not needed to find Newtonian solutions, e.g., the Newtonian tidal tensor can be represented in terms of solutions of (3).

### 1.3. Summary of Results on Relativistic Lagrangian Perturbations

We now discuss the main contents of the series of papers following Paper I [55].

We only briefly mention Paper III [136] in the series that gives the n-th-order perturbation and solution schemes, where emphasis is on the gravito-electric part inspired by the corresponding perturbation and solution schemes in the Lagrange–Newton framework [73]. We discuss in more detail the results of Paper II [137] in the series that invokes spatial averaging of the first-order solutions and provides an explicit form of the backreaction functional. We give a brief account of the basic notions of spatial averaging and the main results of Paper II in what follows.

Spatially averaging Einstein's equations, we can write the result in terms of Friedmann-like equations for a volume scale factor $a_D \propto V_D^{1/3}$, suitably normalized, that appear to be sourced by an effective perfect fluid energy-momentum tensor [24,25]; reviews are [138–140]:

$$ 
\begin{aligned}
3 \frac{\ddot{a}_D}{a_D} &= -4\pi G(\rho_{\text{eff}}^D + 3p_{\text{eff}}^D) + \Lambda; \\
3H_D^2 - \frac{3k_D}{a_D^2} &= 8\pi G\rho_{\text{eff}}^D + \Lambda; \\
\dot{\rho}_{\text{eff}}^D + 3H_D(\rho_{\text{eff}}^D + p_{\text{eff}}^D) &= 0,
\end{aligned} 
\quad (8a) $$

where the effective densities are defined as

$$ 
\begin{aligned}
\rho_{\text{eff}}^D &:= \langle\rho\rangle_D + \rho_{\Phi} \quad ; \quad \rho_{\Phi} := -\frac{1}{16\pi G} Q_D - \frac{1}{16\pi G} W_D; \\
p_{\text{eff}}^D &:= p_{\Phi} \quad ; \quad p_{\Phi} := -\frac{1}{16\pi G} Q_D + \frac{1}{48\pi G} W_D.
\end{aligned} 
\quad (8b) $$

In this form, the effective equations suggest themselves to interpret the extra fluctuating sources in terms of a scalar field (the morphon field [141]), which refers to the inhomogeneities in the geometric functionals $W_D$ and $Q_D$.⁷ The variable $W_D := \langle R \rangle_D - 6k_D/a_D^2$ describes the deviations of the average Ricci scalar curvature from a Friedmannian (eventually scale-dependent) constant curvature $k_D$ on the domain $D$, and the variable $Q_D$ is the backreaction functional, made up of averaged principal scalar invariants of the expansion tensor:

$$ Q_D = 2\langle\text{II}\rangle_D - \frac{2}{3}\langle\text{I}\rangle_D^2. \quad (9) $$

Let us look at two examples. First, we look at the LTB (Lemaître–Tolman–Bondi) solution [144–147]. Being an exact solution, the LTB models have been used as a reference to show that Zel'dovich-type approximations are more accurate than the standard linear and second-order perturbation theory (although the modeling of the velocity field is not as accurate as that for the density field) [148].
Using the relation between the expansion tensor and the metric tensor in the coordinate form $\Theta^i_j := \frac{1}{2}g^{ik}\dot{g}_{kj}$, the averaged scalar invariants of the expansion tensor on a simply-connected LTB-domain can be calculated [137]:

$$ \langle \mathrm{I}(\Theta^i_j) \rangle_{\mathrm{LTB}} = \frac{4\pi}{V_{\mathrm{LTB}}} \int_0^{r_D} \frac{\partial_r (\dot{R}R^2)}{\sqrt{1+2E}} dr; \qquad (10) $$

$$ \langle \mathrm{II}(\Theta^i_j) \rangle_{\mathrm{LTB}} = \frac{4\pi}{V_{\mathrm{LTB}}} \int_0^{r_D} \frac{\partial_r (\dot{R}^2 R)}{\sqrt{1+2E}} dr; \qquad (11) $$

$$ \langle \mathrm{III}(\Theta^i_j) \rangle_{\mathrm{LTB}} = \frac{4\pi}{3V_{\mathrm{LTB}}} \int_0^{r_D} \frac{\partial_r (\dot{R}^3)}{\sqrt{1+2E}} dr, \qquad (12) $$

where R and E are functions of the coordinates r, t, and where the volume is given by

$$ V_{\mathrm{LTB}} = \frac{4\pi}{3} \int_0^{r_D} \frac{\partial_r (R^3)}{\sqrt{1+2E}} dr. \qquad (13) $$

These integrals can be straightforwardly solved for $E(r) = E_0 = \text{const.}$ yielding for the averaged invariants:

$$ \langle \mathrm{II}(\Theta^i_j) \rangle_{\mathrm{LTB}} = \frac{1}{3} \langle \mathrm{I}(\Theta^i_j) \rangle_{\mathrm{LTB}}^2 ; \quad \langle \mathrm{III}(\Theta^i_j) \rangle_{\mathrm{LTB}} = \frac{1}{27} \langle \mathrm{I}(\Theta^i_j) \rangle_{\mathrm{LTB}}^3 . \qquad (14) $$

This proves that, for a spatially Ricci-flat LTB model, the backreaction term vanishes. Second, we calculated in [137] the backreaction model provided by the Relativistic Zel'dovich approximation (RZA), i.e., the first-order solution of the Lagrange–Einstein system:

$$ \mathrm{RZA} Q_D = \frac{\xi^2 (\gamma_1 + \xi\gamma_2 + \xi^2\gamma_3)}{(1 + \xi\langle\mathrm{I}_i\rangle_{C_D} + \xi^2\langle\mathrm{II}_i\rangle_{C_D} + \xi^3\langle\mathrm{III}_i\rangle_{C_D})^2}, \qquad (15a) $$

where we have defined the set of initial data featuring the initial principal scalar invariants of the expansion tensor, their initial average $\langle \dots \rangle_{C_D} := (1/V_D(t_i)) \int_{D(t_i)} \dots J(X^i, t_i) dX^i$, with $J(X^i, t_i) = 1$ (the first term is the initial backreaction term):

$$ \begin{aligned} \gamma_1 &:= 2\langle\mathrm{II}_i\rangle_{C_D} - \frac{2}{3}\langle\mathrm{I}_i\rangle_{C_D}^2 = Q_{C_D}^{\mathrm{initial}}; \\ \gamma_2 &:= 6\langle\mathrm{III}_i\rangle_{C_D} - \frac{2}{3}\langle\mathrm{II}_i\rangle_{C_D}\langle\mathrm{I}_i\rangle_{C_D}; \\ \gamma_3 &:= 2\langle\mathrm{I}_i\rangle_{C_D}\langle\mathrm{III}_i\rangle_{C_D} - \frac{2}{3}\langle\mathrm{II}_i\rangle_{C_D}^2. \end{aligned} \qquad (15b) $$

It is interesting to point out the leading, large-scale behavior of this functional: taking as an example an Einstein–de Sitter background (with vanishing curvature and vanishing $\Lambda$), the leading term is $\mathrm{RZA} Q_D \approx (1/a) Q_{C_D}^{\mathrm{initial}}$. It corresponds to an approximation of the exact scaling solution of the averaged system, $Q_D = (1/a_D) Q_{C_D}^{\mathrm{initial}}$, leading to the same scaling behavior for the average spatial scalar curvature $\langle \mathcal{R} \rangle_D$, via the following coupling equation (following from the last equation of Equation (8a)):

$$ 0 = \frac{1}{a_D^6} \partial_t (Q_D a_D^6) + \frac{1}{a_D^2} \partial_t (\langle \mathcal{R} \rangle_D a_D^2). \qquad (16) $$

This shows that the averaged curvature evolution can be very different from the evolution of a constant-curvature model, and this result can be profitably used to provide a fit to supernova data without dark energy [149,150].

The backreaction functional (15) was implemented numerically in the GPL-licensed code [151]. The INHOMOG library performs a numerical integration of the averaged Ray-
chaudhuri equation for the domain-dependent volume scale factor using the RZA analytical formula for kinematical backreaction.

As a further result, the RZA backreaction model also predicts vanishing backreaction in the case of the averaged Ricci-flat LTB model. (For more details, see [139].) A corresponding Newtonian derivation for the ZA model is provided in [152] and it arises from (15a) and (15b) through the dictionary (4). The result that only a restricted class of averaged LTB models is contained in the RZA backreaction functional forms one of the motivations that led us to generalize the relativistic Lagrangian perturbation theory or the RZA model, to which we come later.

**Remark 1.** In Newtonian cosmology, the backreaction functional $Q_{\mathbb{T}^3}$ also vanishes for the deviations of a Friedmannian background cosmology by employing periodic boundary conditions, i.e., confining the deviations to a flat torus $\mathbb{T}^3$. This is mainly a result of the flatness of space and the fact that the principle scalar invariants of the expansion tensor are integrable, i.e., they can be written as divergences of vector fields as demonstrated in [71]. The backreaction functional can then be written as integrals over flux terms that have to vanish on a toroidal space $\mathbb{T}^3$. The proof can be found in [20].

Paper IV [153] in the series gives a comprehensive investigation of the linearization of the tensorial perturbations, involving decompositions into trace and trace-free parts, where the transverse trace-free part is identified with Lagrangian gravitational waves. This paper also invokes Hodge theory to provide detailed insight into the integrable parts of the solution. It also provides a detailed dictionary of the variables if compared with standard perturbation theory, i.e., it shows which additional steps of linearization are involved to recover the Eulerian-linearized gravitational waves solutions, and it explains the relation of the Hodge decomposition to the standard Scalar–Vector–Tensor decomposition. This paper also makes clear the role of the momentum constraints (6d) for the interpretation of the trace-free part of the Ricci tensor in (6b): these latter are used to split the trace-free part of the Ricci tensor into gravito-electric and gravito-magnetic parts, where the former are tied to the sources (e.g., perturbations evolve according to the trace-solutions), and the latter to freely propagating waves in the linearized system.

As an example, we give here the evolution equations for the gravito-electric and gravito-magnetic spatially projected parts of the Weyl tensor. 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_{ki|k} = -\frac{2}{3a^3} W_{|i}; \qquad (17a)$$

$$H_{ki|k} = 0; \qquad (17b)$$

$$\dot{E}_{ij} + 3H E_{ij} - \frac{1}{a} \epsilon_{(i}^{kl} H_{j)l|k} = -4\pi G \rho_H \sigma_{ij}; \qquad (17c)$$

$$\dot{H}_{ij} + 3H H_{ij} + \frac{1}{a} \epsilon_{(i}^{kl} E_{j)l|k} = 0, \qquad (17d)$$

where $\rho_H = (\rho_H)_i / a^3$ denotes the homogeneous part of the density field, $W = -4\pi G \rho_H \delta$ is related to the density perturbation $\delta$, and $\sigma_{ij}$ are the components of the shear tensor. We obtained the following second-order propagation equations:

$$\Box_0 E_{ij} - 7H \dot{E}_{ij} - 4(4\pi G \rho_H + \Lambda)E_{ij} = -\frac{1}{a^5} D_{ij} W + 4\pi G H \rho_H \sigma_{ij}; \qquad (18a)$$

$$\Box_0 H_{ij} - 7H \dot{H}_{ij} - 4(5\pi G \rho_H + \Lambda)H_{ij} = -\frac{1}{a} 4\pi G \rho_H \epsilon_{(i}^{kl} \sigma_{j)l|k}, \qquad (18b)$$

where $\Box_0 X_{ij} := -\dot{X}_{ij} + \frac{1}{a^2} \Delta_0 X_{ij}$ denotes the Lagrangian d'Alembertian applied to the tensor field $X_{ij}^8$.

This work also illustrates the inherently nonlinear character of a Lagrangian linearization, in contrast to a linearized propagation of gravitational waves written in global
background coordinates, as was mentioned in the Introduction: gravitational waves are perturbations of space and they propagate within the perturbed space in a Lagrangian-linear approach.

Paper V [154] in the series takes pressure into account. Contrary to the Newtonian treatment of a perfect fluid source [86], the general-relativistic approach needs to account for a foliation change, invoking a lapse function in the simplest case [25]. This aspect is also interesting compared with the approach of standard perturbation theory in cosmology: pressure-supported fluids such as a radiation fluid are treated within the same foliation of spacetime, while an intrinsic description of a dust and a radiation fluid, respectively, physically live on different foliations. Since substantial information is drawn from Cosmic Microwave Background radiation properties—such as the ratio of peaks in the spectrum of baryon acoustic oscillations—and compared with late-time structure formation in a dust matter model, such issues might alter those interpretations and it would be worthwhile to dive deeper into this difference.

## 2. Non-Perturbative Generalization of Relativistic Lagrangian Perturbation Theory

Paper VI [155] in the series takes its motivation from the remarks of the spatially averaged LTB model above, but also from an early insight of Masumi Kasai [156] that the relativistic form of Zel'dovich's approximation contains some of the exact Szekeres solutions in a subcase. We face the property that Lagrangian perturbation solutions at first order (RZA) provide an exact description in the maximally anisotropic case, while, on average, the opposite behavior is also contained in a subclass, namely the averaged LTB model. The paper VI [155] headed at making the link to the exact class of Szekeres solutions precise with the main results that (i) Szekeres Class II solutions are contained in a subclass of RZA and they furnish exact deviations from a global FLRW background cosmology for all times, and that (ii) the spatial average of these exact deviations has to vanish. The latter result reveals a crucial restriction of RZA: vanishing averages of deviations from a global FLRW background imply that these models cannot describe global backreaction, i.e., the inhomogeneous deviations evolve according to the chosen background cosmology, but the background cosmology is not affected by the inhomogeneities. This construction implies a conservation law for the average scalar curvature in conformity with the background curvature evolution, which is not the case in general [157]: changes are due to the coupling of structure formation to the average scalar curvature according to the general conservation law Equation (16). This remark is equally relevant for general-relativistic numerical simulations that are globally constructed in the same way: the evolution of the global scalar curvature can be used to test as to whether the simulation allows for global deviations from the background curvature model, or whether its architecture is, by construction, limited to conserve the imposed background curvature.

In view of these remarks, it is natural to ask whether we could generalize the relativistic first-order solutions and their averages in such a way that, locally, also Szekeres Class I solutions are contained in a subclass and hence also LTB models with curvature. Before we come to the proposed generalization, we shall look into the family of Szekeres solutions.

### 2.1. Szekeres Solutions in the Goode and Wainwright Parameterization

The family of solutions found by Szekeres [158] and later generalized by Szafron [159] constitutes the most general exact solution to Einstein's equations that can be regarded as an inhomogeneous cosmological model. The general line-element can be written in the form⁹ [144,146]:

$$ds^2 = -dt^2 + e^{2\alpha} dz^2 + e^{2\beta} (dx^2 + dy^2), \qquad (19)$$

where the metric functions $\alpha$ and $\beta$, in general, depend on all the comoving coordinates. The field source of the Szekeres solution is a perfect fluid with inhomogeneous density but an only time-dependent pressure $p = p(t)$, which for physical models reduces to
an irrotational but inhomogeneous dust fluid with the cosmological constant as the dark energy model [160]—see [161] for a discussion on the limitations imposed by this property.

In these coordinates, the solution naturally splits into two classes: class I when $β_{,z} \neq 0$ and class II when $β_{,z} = 0$. In general, neither class has symmetries [162], and other well-known exact solutions are contained as particular cases, e.g., the LTB model follows from the class I, and some subclasses of the Bianchi models emerge from the homogeneous limit of class II. In general, both classes have a well-defined FLRW limit. Szekeres models have been used to examine several astrophysical and cosmological problems, being the first class the one that has caught the most attention [163–182].

The original solution has been reparameterized multiple times. Among these parameterizations, the Goode and Wainwright parameterization [183] is the one that provides the most direct connection with RZA. In this representation, the line-element (19) is rewritten as [144,146]

$$ds^2 = -dt^2 + S^2(G^2W^2dz^2 + e^{2\nu}(dx^2 + dy^2)), \qquad (20)$$

where the scale factor $S$ obeys a Friedmann-like equation:

$$\left(\frac{\dot{S}}{S}\right)^2 = -\frac{k_0}{S^2} + \frac{2\mu}{S^3} + \frac{\Lambda}{3}, \qquad (21)$$

with $k_0 = 0, \pm 1$ and $\mu > 0$ being an arbitrary function. In addition,

$$G = A(r) - F(t, z) \quad \text{with} \quad F := \beta_{+}f_{+} + \beta_{-}f_{-}. \qquad (22)$$

One of the most enlightening features of this parameterization is that $f_{\pm}$ (and then $F$) satisfy the same equation that leads to the growing and decaying modes in cosmological perturbation theory:

$$\ddot{F} + 2\frac{\dot{S}}{S}\dot{F} - \frac{3\mu}{S^3}F = 0. \qquad (23)$$

The functions $A(r)$, $e^{2\nu(r)}$, $W(z)$, $\beta_+(z)$ and $\beta_-(z)$ differ for each class, see their specific form in [146,183] or [155], the latter for a presentation using the notation used in this paper. Here, we restrict ourselves to highlight the most important characteristics of each class:

### Class I:

$S = S(t, z)$, $\mu = \mu(z)$, $f_{\pm} = f_{\pm}(t, z)$, $\beta_{+} = -k_0 f_{\frac{\mu}{3\mu}}$ and $\beta_{-} = f_{T,z}$ where $T(z)$ is the time of the initial singularity and $f(z)$ is an arbitrary function.

### Class II:

$S = S(t)$, $f_{\pm} = f_{\pm}(t)$, $\mu = \text{const.}$, $W = 1$.

Regardless of the class, the energy density can be integrated to the following compact expression:

$$4\pi\rho(t, r) = \frac{3\mu}{S^3}\left(1 + \frac{F}{G}\right). \qquad (24)$$

## 2.2. Szekeres Models and Lagrangian Perturbation Theory

The starting point to relate the Lagrangian formalism and Szekeres exact solutions is the comparison of their line-elements. We first identify the initial metric ($G_{ij}$) and the exact deviation ($h_{ij}$):

$$g_{ij} \equiv a^2(t)\gamma_{ij} \equiv a^2(t)(G_{ij} + h_{ij}), \qquad (25)$$

and then find the Szekeres deformation field by solving the equations from the equality of the Szekeres metric and the bilinear functional of RZA:

$$h_{ij} \equiv G_{ab} (\delta_i^a P_j^b + \delta_j^b P_i^a + P_i^a P_j^b). \qquad (26)$$
Finally, we compare the evolution equation(s) for the thus defined Szekeres deformation field with the model equations of RZA.

At this point, it is convenient to split the analysis into the two different classes and treat them separately.

### 2.2.1. Class II Models as the Exact Limit of RZA

For this class, the parameter $\mu$ is a constant, and the conformal scale factor is a function of the time only. If we set $\mu$ equal to the energy density of the homogeneous universe model at the initial time ($3\mu = 4\pi\rho_b(t_i)$), Equation (21) is precisely the Friedmann equation. Below, we will relabel $S = a(t)$ and assume that $a(t)$ is initially normalized: $a(t_i) = 1$. The initial metric is via Gram's matrix determined by the direct evaluation of (20) at $t = t_i$:

$$G_{ab} = \text{Diag}[\tilde{\mathcal{A}}^2, e^{2\nu}, e^{2\nu}] \quad \text{with} \quad \tilde{\mathcal{A}} := \mathcal{A} - \mathcal{F}_i \quad (\text{and} \quad \mathcal{F}_i := \mathcal{F}(t_i)). \qquad (27)$$

Comparing the line-elements (20) and (26), we obtain that the only non-vanishing component of $h_{ij}$ is $h_{33}$ (or $h_{zz}$), which is given by

$$h_{33} = -2\mathcal{A}\tilde{\mathcal{F}} + \tilde{\mathcal{F}}^2 = G_{ab}(\delta_3^a P_3^b + \delta_3^b P_3^a + P_3^a P_3^b) = 2\tilde{\mathcal{A}}^2 P_3^3 + \tilde{\mathcal{A}}^2 (P_3^3)^2, \qquad (28)$$

where we have defined $\tilde{\mathcal{F}} = \mathcal{F} - \mathcal{F}_i$. Thus, $P_3^3$ is the only non-trivial component of the deformation field:

$$P_3^3 = P = -\tilde{\mathcal{F}}/\tilde{\mathcal{A}}, \qquad (29)$$

furnishing a locally one-dimensional solution.

It can also be verified by direct substitution that $P$ obeys the evolution equation for the trace part of the deformation field—the only non-trivial one for Szekeres and equivalent to Equation (23) for Szekeres models. In addition, all Szekeres fields are reproduced by the functional expressions of RZA. This analysis can be summarized in the following remark.

**Remark 2.** The class II of Szekeres models is contained in RZA within its locally one-dimensional case, and the dynamics is exclusively governed by the first principal scalar invariant of the deformation field. Applying the dictionary (4), we find the exact 3D solutions without symmetry in Newtonian cosmology [43].

Hence, Szekeres class II emerges from a proper set of the initial data, for which the constraints are satisfied on the initial hypersurface. This result illustrates the nonlinear nature of RZA. We note, in particular, that the metric is a bilinear functional of the linearized perturbation, hence quadratic, and the density is given by the exact integral of the continuity equation with the second and third principal scalar invariants vanishing. As for the averages, since the background model is global for this class, the deviation fields lead to a vanishing backreaction functional for periodic boundary conditions [155].

Although these solutions (Szekeres II and its Newtonian analog) are valuable for understanding the theory, their applications as cosmological models are limited by the difficulty of finding initial data without singularities. In Szekeres II models with a flat associated background, the presence of a growing mode necessarily leads to singularities on the initial hypersurface. Still, these models have found cosmological applications considering only the regular part of the manifold [184,185].

### 2.2.2. Class I Models and Lagrangian Perturbation Theory

The relationship with the Lagrangian perturbation theory is more involved for class I than for class II models. Since class I lacks a global background, all relevant dynamical functions ($S$ and $f_\pm$) exhibit a nontrivial dependence on the spatial coordinates, making it impossible to separate the space- and time-dependencies.
Still, we consider the co-frame set as the fundamental dynamical field and retain the mathematical structure of the Lagrangian perturbation theory. In this way, the field $\hat{P}_i^a$ is reinterpreted as a local deformation field on a Friedmann-like local “reference model” described by Equation (21)$^{10}$. Following similar arguments to those used for class II, it can be shown that such a deformation field has only one non-trivial component, which has the same functional expression as its equivalent for class II, Equation (29).

To dig more deeply into the connection of RZA and the exact solutions, we exploit the silent property of Szekeres, under which each world line evolves “independently” of the others [186–188]. Let us consider an arbitrary world line labeled by the comoving coordinates $r_*$.

$$ \hat{a}(t) := S(t, z)|_{z=r_*} \qquad (30) $$

obeys the Friedmann equations for a model with initial density $4\pi q_b(t_i) = 3\mu(z)|_{z=r_*}$: the “reference model”. Then, the local deformation field introduced in the above paragraph satisfies the RZA evolution equation for its trace part (the only non-trivial one). As for class II, this evolution equation reduces to (23), and $f_\pm(t, z)|_{z=r_*}$ are interpreted as the growing and decaying modes on the “reference model”.

**Remark 3.** The dynamics of the Szekeres models of class I can be interpreted as a constrained superposition of world lines satisfying the RZA model equations sourced by different “backgrounds”, i.e., a local Friedmann-like reference model. All the local quantities characterizing the world lines can also be obtained from their RZA functional expressions. Averaging over these “local evolutions” will provide a new global average model that interacts with the inhomogeneities, evolves differently compared with a Friedmannian model and, hence, includes cosmological backreaction.

This analysis paves the way for a strategy to generalize RZA so that it can reproduce the entire family of Szekeres solutions (and thus the LTB models) for a proper set of the initial data.

### 2.2.3. Rescaled Class I Solution

To examine the relation with RZA, it is convenient to normalize the initial conformal scale factor. While for class II this can be easily done by a re-scaling of the comoving coordinates, for class I, the initial scale factor is a function of z: $S(t_i) = S_i(z)$. Then, we define $A := S/S_i$, which clearly satisfies $A(t_i) = 1$. Proceeding along these lines, Equation (21) keeps the form of a Friedmann equation but in terms of space-dependent curvature and density parameters, $\hat{k}(z) := k_0/S_i^2$ and $4\pi\hat{\rho}(z) := 3\mu/S_i^3$:

$$ \left(\frac{\dot{A}}{A}\right)^2 = -\frac{\hat{k}}{A^2} + \frac{8\pi}{3}\frac{\hat{\rho}}{A^3} + \frac{\Lambda}{3}, \qquad (31) $$

and $\mathcal{F}$ remains as the superposition of the growing and decaying modes on the “reference model” defined by (31):

$$ \ddot{\mathcal{F}} + 2\frac{\dot{A}}{A}\mathcal{F} - \frac{4\pi\hat{\rho}}{A^3}\mathcal{F} = 0. \qquad (32) $$

Other minor changes to the arbitrary functions are also needed to leave the solution unaltered.

### 2.3. Generalized RZA: An Exact-Solution-Controlled Model of Structure Formation$^{11}$

Although RZA cannot reproduce the class I of Szekeres models, both solutions are governed by the same set of evolution equations. Their difference purely relies on the constraints. This is a characteristic feature of the silent class of exact solutions to Einstein’s equations, where some local rotational symmetric and Bianchi solutions and the LTB models have a common set of evolution equations. It is precisely the constraints that make them different [186].
Motivated by these results, we propose a generalization of RZA with as minimal modifications as possible to contain Szekeres class I as a particular case, while still retaining the mathematical structure of the Lagrangian perturbation theory. Proceeding along these lines, the co-frames are kept as the fundamental field of the approach. Next, we introduce a Friedmann-like “reference model” as a generalization of the homogeneous background and the deformation field on it,

$$ \eta^a = \eta_i^a \mathbf{d}X^i = A(\delta_i^a + \hat{P}_i^a) \mathbf{d}X^i, \quad A = A(t, \mathbf{r}), \quad \hat{P}_i^a = \hat{P}_i^a(t, \mathbf{r}). \qquad (33) $$

This “reference model” satisfies Equation (31) with a general curvature ($\hat{k} = \hat{k}(\mathbf{r})$) and density parameter ($\hat{\rho}(\mathbf{r})$).

The subsequent approach consists of inserting this co-frame set into the Lagrange-Einstein system (6) and linearizing the evolution equations in the deformation field. As in RZA, this formal solution for $\hat{P}_i^a$ is injected into the nonlinear functional expressions for the physical fields, e.g., density, 3D Ricci scalar, expansion scalar, and the (bilinear) metric. We call this approach Generalized RZA (GRZA).

It is instructive to rephrase the previous analysis in different terms. GRZA results from adding two extra degrees of freedom to RZA; namely, the constant curvature and density of the FLRW model are generalized by space-dependent functions:

$$ \text{RZA} \rightarrow \text{GRZA}: \quad k_0 \rightarrow \hat{k}(\mathbf{r}) \wedge \rho_b(t_i) \rightarrow \hat{\rho}(\mathbf{r}). \qquad (34) $$

This architecture allows reproducing the class I of the Szekeres solutions and its well-known particular cases, e.g., the LTB models. In addition, RZA and the class II of Szekeres emerge in the limit when the curvature and density parameters of the reference model become constant, see Figure 2.

<figure><img src="image_6.png" /><figcaption>Figure 2. Special subcases of GRZA. GRZA contains class I of the Szekeres (Sze-I) solutions and all its subcases, such as the LTB models. In addition, RZA and the class II of Szekeres (Sze-II) emerge in the limit when $\hat{k}(\mathbf{r}) \rightarrow k_0 = \text{const.}$ and $\hat{\rho}(\mathbf{r}) \rightarrow \rho_b(t_i) = \text{const.}$</figcaption></figure>

## An Example of GRZA Numerical Simulation

GRZA has enormous potential to model the large-scale structure formation in general relativity without suffering from the constraints of a global background architecture. As a proof of principle, we show a simple numerical example obtained from a set of Szekeres-like initial data at the last scattering time. This family of models is characterized by a locally one-dimensional deformation field, where the Szekeres exact expression is modified to include a fluctuation proportional to the growing and decaying modes:

$$ \hat{P}(t, \mathbf{r}) = \frac{\hat{\mathcal{F}} - \hat{\mathcal{F}}_i}{\mathcal{A} - \hat{\mathcal{F}}_i}, \quad \hat{\mathcal{F}} = \mathcal{F} + \delta\mathcal{F}, \quad \delta\mathcal{F} = \delta\beta_+ f_+ + \delta\beta_- f_- . \qquad (35) $$

In the equation above, $\mathcal{F} = \beta_+ f_+ + \beta_- f_-$ is the Szekeres function introduced in (22).

As thus built, these models retain the characteristic foliation of the hypersurfaces by non-concentric two-spheres of the quasi-spherical models, but generalize the Szekeres solution in the way these 2-spheres are arranged. Figures 3 and 4 show the present-day
projection of the density contrast in the equatorial plane for different solutions with different deviations. The function $\delta\beta_+$ is assumed to have the form $\beta_+ \propto \alpha \sin^3(\gamma\pi y) \sin^4(\pi z)$; hence, the case $\alpha = 0$ corresponds to the exact solution. While the models with $\alpha \neq 0$ are not exact, the violation of the constraints remains small throughout the entire evolution, indicating that they are valid approximate solutions of general relativity. (For these simulations, the largest value of the violation of the constraints at the present cosmic time is of the order of $10^{-3}$ for the case depicted in Figure 4).

Figure 3 illustrates the breaking of the Szekeres “dipolar symmetry”. As can be seen in the sequence of Figure 3a-c, the characteristic Szekeres pancake-shaped overdensity splits into three structures (at the present cosmic time) as the value of $\alpha$ increases. Then, Figure 4 shows the density contrast for a larger value of $\alpha$, leading to a network of structures undergoing nonlinear collapse. To achieve a better visualization, in this graph, we have used a different color scale than the one used in Figure 3—the scale is shown in the right bar of each graph.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Density contrast in the equatorial plane at the present cosmic time. The case $\alpha = 0$ corresponds to an exact Szekeres model; then, the Szekeres pancake-shaped overdensity splits into three structures. The coordinates are defined as $X = A(t_0, z)\cos\phi$ and $Y = A(t_0, z)\sin\phi$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Density contrast in the equatorial plane at the present cosmic time for $\alpha = 5.2 \times 10^{-3}$. The coordinates are defined as in Figure 3.</figcaption></figure>

**3. Summary, Discussion and Outlook**

After a historical review on the development of the Lagrangian perturbation theory in Newtonian cosmology, we recalled the Lagrange–Newton system (3) governing the defor-
mation gradient $df^a$, which is the starting point of a Lagrangian perturbation expansion on a global background cosmology. Since the deformation gradient is the only dynamical variable, we are entitled to calculate any field variable through its exact definition for a given order of the perturbative expansion.

Through the introduction of a simple dictionary, which replaces the exact one-form fields $df^a$ by general co-frame fields $\eta^a$, we obtain a general-relativistic form of the Lagrangian deformation and its perturbative expansions. We recalled that the same logic applies to the Einstein equations, exemplified for an irrotational dust continuum in a flow-orthogonal foliation of spacetime, namely that the Einstein equations can be represented through the Lagrange-Einstein system (6) for the co-frame fields as the only dynamical variable. Consequently, we are entitled to calculate any variable in terms of exact functionals of the co-frame fields, e.g., the metric as a bilinear quadratic form of the co-frames, the expansion and Ricci tensors, etc., valid for any order of the perturbative deformation that we consider. We remark that the Lagrange-Einstein system can be generalized for general energy-momentum tensors and general foliations of spacetime, including tilted foliations (for a comprehensive investigation and the definition of the Lagrangian picture in the general situation, see [15,16]).

We reviewed the main contents of a series of papers on the general-relativistic form of the Lagrangian perturbation theory on a global background cosmology. In particular, we also considered average properties of the Lagrangian first-order solutions that gave rise to a compact backreaction functional capturing the deviations from the assumed background cosmology. We listed a few known theorems: backreaction vanishes for Newtonian cosmological models that employ periodic boundary conditions. The same applies for spherically symmetric models on a flat space section. These results carry over to corresponding theorems in general relativity, i.e., for simulations that impose the torus architecture on flat space sections, and for flat LTB models. In a next step, we looked at the class of Szekeres solutions, and we have shown that Szekeres Class II is contained in a subclass of the first-order relativistic Lagrangian perturbation solutions (RZA model). Integral constraints can be imposed to comply with the torus architecture of relativistic simulations that are conceived as deviations from a global background cosmology.

We then motivated a non-perturbative generalization of the Lagrangian perturbation theory with the aim to also include Szekeres Class I and general LTB models. This could be achieved by using the Goode-Wainwright parameterization of the Szekeres solutions, since this allowed us to represent the solutions in terms of a Friedmann-like equation for the “background” and deviations thereof that also comply with the form of the perturbation equations known in Lagrangian perturbation theory. We presented the results that (i) Szekeres Class II solutions with a homogeneous scale factor correspond to the general-relativistic form of the Lagrangian first-order perturbation solutions, while (ii) Szekeres Class I solutions with an inhomogeneous scale factor still allow us to employ the Lagrangian formalism, but only locally.

We gave a numerical example of a Szekeres Class I model on the scale of a collapsing structure in comparison with initial data that were slightly deformed within the class of the corresponding Lagrangian perturbation solution. These deformations we controlled via the condition of smallness of violation of the energy constraint.

To proceed further towards a generic situation for the formation of large-scale structure in the Universe, we have to implement a local background and a corresponding deviation dynamics that is sourced by this local background evolution for a large-scale simulation. As this construction is local, we have to average the resulting dynamics to obtain the physical average cosmology. By construction, the resulting generalization of the Lagrangian model (GRZA) contains Szekeres Class I solutions in a subclass, but it also provides an exact-solution-controlled approximation that features an evolving average distribution. The result is a model that is able to quantify the global cosmological backreaction. It thus realizes the expectation that fluctuations “talk” to the “global background”, conceived as the average model.
Further investigations are needed to numerically realize this new approximation. We point out that an alternative implementation is possible: instead of locally splitting the dynamics into a “local background” and deviations thereof, we may directly use the averaged variables such as the backreaction functional. This approach is quasi-local in the sense that average properties are implemented on a domain of the size of the resolution scale that we choose. The limit where the averaging domain $\mathcal{D}$ tends to zero agrees with the local dynamics, so that we may apply averages to sufficiently small domains, then summing up the quasi-local results by using the exact volume-partitioning formula for the backreaction functional [157], and the exact volume-partitioning formulas of a multiscale dynamics [190]. Control on this procedure is provided by the spatially averaged constraints of general relativity.

A further alternative that is independent of the here proposed strategy is iterative rather than perturbative in nature: we can set up the deviation dynamics on a large-scale or “global background model” that is identified with the exact average dynamics. This approach has been outlined in [18].

**Author Contributions:** T.B.: conceptualization, writing, methodology, supervision, funding acquisition; I.D.G.: writing, methodology, software, visualization; J.J.O.: writing, funding acquisition, software supervision. All authors have read and agreed to the published version of the manuscript.

**Funding:** 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 programme (grant agreement ERC advanced grant 740021-ARTHUS, PI: T.B.). I.D.G. and J.J.O. acknowledge support from the National Science Centre (NCN, Poland) under the Sonata-15 research grant UMO-2019/35/D/ST9/00342. I.D.G. acknowledges hospitality and support by ERC ARTHUS during a working visit in Lyon.

**Data Availability Statement:** The INHOMOG code referred to in the context of realizing the backreaction functional in RZA is openly available [151]. A pedagogical jupyter notebook to realize the example for GRZA in comparison with Szekeres Class I solutions is available at https://github.com/idgaspar/Beyond_relativistic_Lagrangian_perturbation_theory (accessed on 26 September 2022).

**Acknowledgments:** We would like to thank the guest editors Giovanni Marozzi and Luigi Tedesco for inviting us to contribute to this special issue. We would like to emphasize the fruitful inspirations that we got from Maurizio Gasperini's work on covariant averaging formalisms that he investigated with Gabriele Veneziano and collaborators (not being the main subject of this review). TB wishes to thank Oliver Hahn and Cornelius Rampf for valuable discussions during a visit to Vienna Observatory, Austria. We are thankful to a referee for insightful comments and suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Notes

1. For further discussion of and references to the singularity problem, see [55], Section VB.

2. In this context, it is important to remark that the best numerical methods are simulations in phase space that provide a controlled and detailed access of the formation of structure via a projection into space [60]. Since the phase space is six-dimensional, this, however, provides limits of resolution. For recent simulations, see [61,62] and the review [63]. For a discussion of eigenvalues of the deformation tensor, see [64].

3. In order to achieve even higher resolution, it is possible to interpolate initial data, since they are smooth and the analytical scheme maps any (calculated or interpolated) point to the nonlinear stage.

4. Recently, using newer numerical methods, Ref. [111] provided a detailed generalization of the 2D pancakes of [51].

5. We have distinguished counter-indices $a, b, c, \dots$, here counting vector components, from coordinate indices $i, j, k, \dots$. We shall later see that this is convenient when a coordinate basis does not exist; here, the vector components are also coordinate indices in the Eulerian basis.

6. Starting with Paper II in the series, we employ the metric form coefficients $g_{ij} = G_{ab}\eta^a_i\eta^b_j$, using Gram's matrix $G_{ab}$ instead of $\delta_{ab}$ for orthonormal co-frames. It has the advantage that the initial deformation can be set to be undeformed and therefore provides a closer correspondence to Newtonian solutions. The initial metric is then not encoded in initial deformations but in Gram's matrix (for details, see [136,137]); for the Euclidean restriction (integrable limit) of Gram's matrix, see Appendix A [136]), and for the integrable limit of the gravito-magnetic part, see Sect. III.A.2 [55]. Note also that Equations (6a) imply $G_{ab} \eta^a_i \eta^b_j = 0$.
7. This rewriting of the backreaction terms in the form of a minimally coupled effective scalar field can be profitably used to establish quintessence models without dark energy, or inflationary models without fundamental scalar field [142], or models mimicking dark matter [143].

8. The coefficients in front of the third terms in Equations (18a) and (18b) differ since the first equation features a term $\propto \dot{\sigma}_{ij}$ that can be replaced through $\dot{\sigma}_{ij} = -E_{ij}$, see Equation (111) in [55], which changes the coefficient $5\pi G\rho_H + \Lambda$ to $4\pi G\rho_H + \Lambda$.

9. Henceforth, we use units where $c = G = 1$.

10. To avoid confusion, we use $\hat{P}_i^a$ to refer to the deformation field with respect a generic reference model and reserve $P_i^a$ for the deformation with respect to the FLRW background.

11. In this section, we highlight some of the salient results of a detailed investigation [189].

References

1. Brustein, R.; Gasperini, M.; Giovannini, M.; Mukhanov, V.F.; Veneziano, G. Metric perturbations in dilaton-driven inflation. *Phys. Rev. D* **1995**, *51*, 6744. [CrossRef] [PubMed]

2. Gasperini, M. Tensor perturbations in high-curvature string backgrounds. *Phys. Rev. D* **1997**, *56*, 4815. [CrossRef]

3. Lifshitz, E.M. On the gravitational stability of the expanding universe. *J. Phys.* **1946**, *10*, 116–129; reprinted in *Gen. Rel. Grav.* **2017**, *49*, 18. [CrossRef]

4. Hawking, S.W. Perturbations of an Expanding Universe. *Astrophys. J.* **1966**, *145*, 544. [CrossRef]

5. Lifshitz, E.M.; Khalatnikov, I.M. Investigations in relativistic cosmology. *Adv. Phys.* **1963**, *12*, 185–249; reprinted in *Sov. Phys. Usp. 1964*, *6*, 522. [CrossRef]

6. Tomita, K. Non-Linear Theory of Gravitational Instability in the Expanding Universe. *Prog. Theor. Phys.* **1967**, *37*, 831–846. [CrossRef]

7. Durrer, R. Gauge Invariant Cosmological Perturbation Theory: A General Study and Its Application to the Texture Scenario of Structure Formation. *Fundam. Cosm. Phys.* **1994**, *15*, 209.

8. Kodama, H.; Sasaki, M. Cosmological Perturbation Theory. *Prog. Theor. Phys. Suppl.* **1984**, *78*, 1. [CrossRef]

9. Malik, K.A.; Wands, T. Cosmological perturbations. *Phys. Rep.* **2009**, *475*, 1. [CrossRef]

10. Mukhanov, V.F.; Feldman, H.A.; Brandenberger, R.H. Theory of cosmological perturbations. *Phys. Rep.* **1992**, *215*, 203. [CrossRef]

11. Ehlers, J. Beiträge zur relativistischen Mechanik kontinuierlicher Medien. *Abh. Math.-Naturwissenschaftlichen Kl.* **1961**, 793–836; Translated: Contributions to the relativistic mechanics of continuous media. *Gen. Rel. Grav.* **1993**, *25*, 1225–1266. [CrossRef]

12. Bardeen, J.M. Gauge-invariant cosmological perturbations. *Phys. Rev. D* **1980**, *22*, 1882. [CrossRef]

13. Stewart, J.M. Perturbations of Friedmann-Robertson-Walker cosmological models. *Class. Quant. Grav.* **1990**, *7*, 1169. [CrossRef]

14. Ellis, G.F.R.; Bruni, M. Covariant and gauge-invariant approach to cosmological density fluctuations. *Phys. Rev. D* **1989**, *40*, 1804. [CrossRef] [PubMed]

15. Buchert, T.; Mourier, P.; Roy, X. On cosmological backreaction and its dependence on spacetime foliation. *Class. Quant. Grav.* **2018**, *35*, 24LT02. [CrossRef]

16. Buchert, T.; Mourier, P.; Roy, X. On average properties of inhomogeneous fluids in general relativity III: General fluid cosmologies. *Gen. Rel. Grav.* **2020**, *52*, 27. [CrossRef]

17. Vitenti, S.D.P.; Falciano, F.T.; Pinto-Neto, N. Covariant Bardeen perturbation formalism. *Phys. Rev. D* **2014**, *89*, 103538. [CrossRef]

18. Roy, X.; Buchert, T. Relativistic cosmological perturbation scheme on a general background: Scalar perturbations for irrotational dust. *Class. Quant. Grav.* **2012**, *29*, 115004. [CrossRef]

19. Combi, L. Spacetime is material. In *Contemporary Materialism. Its Ontology and Epistemology*; Romero, E.G., Pérez Jara, J., Camprubí, L., Eds.; Synthese Library, Springer: Berlin/Heidelberg, Germany, 2021.

20. Buchert, T.; Ehlers, J. Averaging inhomogeneous Newtonian cosmologies. *Astron. Astrophys.* **1997**, *320*, 1.

21. Ellis, G.F.R.; Buchert, T. The universe seen at different scales. *Phys. Lett. A* **2005**, *347*, 38. [CrossRef]

22. Buchert, T.; Carfora, M. Regional averaging and scaling in relativistic cosmology. *Class. Quant. Grav.* **2002**, *19*, 6109. [CrossRef]

23. Carfora, M.; Piotrkowska, K. A Renormalization Group Approach to Relativistic Cosmology. *Phys. Rev. D* **1995**, *52*, 4393. [CrossRef] [PubMed]

24. Buchert, T. On average properties of inhomogeneous fluids in general relativity I: Dust cosmologies. *Gen. Rel. Grav.* **2000**, *32* 105. [CrossRef]

25. Buchert, T. On average properties of inhomogeneous fluids in general relativity II: Perfect fluid cosmologies. *Gen. Rel. Grav.* **2001**, *33*, 1381. [CrossRef]

26. Buchert, T.; van Elst, H.; Heinesen, A. The averaging problem on the past null cone in inhomogeneous dust cosmologies. *Gen. Rel. Grav.* **2022**, submitted.

27. Gasperini, M.; Marozzi, G.; Veneziano, G. A covariant and gauge invariant formulation of the cosmological 'backreaction'. *J. Cosmol. Astropart. Phys.* **2010**, *2*, 009. [CrossRef]

28. Gasperini, M.; Marozzi, G.; Nugier, F.; Veneziano, G. Light-cone averaging in cosmology: Formalism and applications. *J. Cosmol. Astropart. Phys.* **2011**, *7*, 008. [CrossRef]

29. Heinesen, A.; Mourier, P.; Buchert, T. On the covariance of scalar averaging and backreaction in relativistic inhomogeneous cosmology. *Class. Quant. Grav.* **2019**, *36*, 075001. [CrossRef]
30. Zel'dovich, Y.B.; Myshkis, A.D. *Elements of Mathematical Physics, Medium of Non-Interacting Particles*; Nauka: Moscow, Russia, 1973.

31. Dawson, J.M. *Nonlinear Electron Oscillations* in a Cold Plasma. *Phys. Rev.* **1959**, *113*, 383. [CrossRef]

32. Peebles, P.J.E. *The Large-Scale Structure of the Universe*; Princeton University Press: Princeton, JS, USA, 1980; ISBN 9780691209838.

33. Zel'dovich, Y.B. *Gravitational instability: An approximate theory for large density perturbations*. *Astron. Astrophys.* **1970**, *5*, 84.

34. Zel'dovich, Y.B. *Fragmentation of a homogeneous medium under the action of gravitation*. *Astrofizika* **1970**, *6*, 319; Translated in *Astrophysics* **1973**, *6*, 164–174. [CrossRef]

35. Zel'dovich, Y.B. *The Large Scale Structure of the Universe*. In *International Astronomical Union Symposia*; Longair, M.S., Einasto, J., Eds.; Springer: Berlin/Heidelberg, Germany, 1978; Volume 79, p. 409. ISBN 978-94-009-9843-8.

36. Zel'dovich, Y.B.; Shandarin, S.F. Maximum density in heavy-neutrino clouds. *Sov. Astron. Lett.* **1982**, *8*, 139; Translated in *Pis'ma Astron. Zh.* **1982**, *8*, 259–262.

37. Doroshkevich, A.G.; Ryabenkii, V.S.; Shandarin, S.F. *Non-linear theory of development of potential perturbations*. *Astrofizika* **1975**, *9*, 257; Translated in *Astrophysics* **1973**, *9*, 144.

38. Bartelmann, M.; Schneider, P. *A large-scale structure model for gravitational lensing*. *Astron. Astrophys.* **1992**, *259*, 413.

39. Doroshkevich, A.G. The space structure of perturbations and the origin of rotation in galaxies in the theory of fluctuation. *Astrofizika* **1973**, *6*, 581–600; Translated in *Astrophysics* **1970**, *6*, 320–330. [CrossRef]

40. Grinstein, B.; Wise, M.B. On the validity of the Zel'dovich Approximation. *Astrophys. J.* **1987**, *320*, 448; Erratum in *Astrophys. J.* **1989**, *337*, 579. [CrossRef]

41. Fillmore, J.A.; Goldreich, P. *Self-similar gravitational collapse in an expanding Universe*. *Astrophys. J.* **1984**, *281*, 1. [CrossRef]

42. Bartelmann, M. *Trajectories of point particles in cosmology and the Zel'dovich approximation*. *Phys. Rev. D* **2014**, *91*, 083524. [CrossRef]

43. Buchert, T. *A class of solutions in Newtonian cosmology and the pancake theory*. *Astron. Astrophys.* **1989**, *223*, 9.

44. Shandarin, S.F.; Zel'dovich, Y.B. The large-scale structure of the universe: Turbulence, intermittency, structures in a self-gravitating medium. *Rev. Mod. Phys.* **1989**, *61*, 185. [CrossRef]

45. Lin, C.C.; Mestel, L.; Shu, F.H. *The Gravitational Collapse of a Uniform Spheroid*. *Astrophys. J.* **1965**, *142*, 1431. [CrossRef]

46. Goodman, J.; Binney, J. *Gravitational collapse of pressureless inhomogeneous spheroids*. *Mon. Not. R. Astron. Soc.* **1983**, *203*, 265. [CrossRef]

47. Yoshisato, A.; Matsubara, T.; Morikawa, M. *Accuracy of Nonlinear Approximations in Spheroidal Collapse: Why Are Zel'dovich-Type Approximations Good?* *Astrophys. J.* **1998**, *498*, 48. [CrossRef]

48. Doroshkevich, A.G.; Kotok, E.V.; Novikov, I.D.; Polyudov, A.N.; Shandarin, S.F.; Sigov, Y.S. *Two-dimensional simulation of the gravitational system dynamics and formation of the large-scale structure of the universe*. *Mon. Not. R. Astron. Soc.* **1980**, *192*, 321. [CrossRef]

49. Klypin, A.A.; Shandarin, S.F. *Three-dimensional numerical model of the formation of large-scale structure in the Universe*. *Mon. Not. R. Astron. Soc.* **1983**, *204*, 891. [CrossRef]

50. Melott, A.L.; Shandarin, S.F. *Gravitational Instability with High Resolution*. *Astrophys. J.* **1989**, *343*, 26. [CrossRef]

51. Arnol'd, V.I.; Shandarin, S.F.; Zel'dovich, Y.B. The large-scale structure of the Universe. I. General properties. One- and two-dimensional models. *Geophys. Astrophys. Fluid Dyn.* **1982**, *20*, 111. [CrossRef]

52. Arnol'd, V.I.; Gusein-Sade, S.M.; Varchenko, A.N. * singularities of Differentiable Maps*; Birkhäuser: Boston, MA, USA, 1985; USSR Academy of Sciences: Moscow, Russia, 1982; Volume I.

53. Shandarin, S.F. *Evolution of Perturbations in Friedmann Models of the Universe*. *Astrofizika* **1980**, *16*, 769–779; Translated in *Astrophysics* **1980**, *16*, 439–445. [CrossRef]

54. Zentsova, A.S.; Chernin, A.D. *Evolution of entropy perturbations in the post-recombination epoch II—Nonlinear stage*. *Astrofizika* **1980**, *16*, 169–177; Translated in *Astrophysics* **1980**, *16*, 108–113. [CrossRef]

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

56. Shukurov, A.M. *Nonlinear growth of density perturbations in an expanding collisionless medium*. *Astrofizika* **1981**, *17*, 469; Translated in *Astrophysics* **1982**, *17*, 263. [CrossRef]

57. Gurbatov, S.N.; Saichev, A.I.; Shandarin, S.F. The large-scale structure of the universe in the frame of the model equation of non-linear diffusion. *Mon. Not. R. Astron. Soc.* **1989**, *236*, 385. [CrossRef]

58. Buchert, T.; Domínguez, A. *Modeling multi-stream flow in collisionless matter: Approximations for large-scale structure beyond shell-crossing*. *Astron. Astrophys.* **1998**, *335*, 395.

59. Buchert, T.; Domínguez, A. *Adhesive Gravitational Clustering*. *Astron. Astrophys.* **2005**, *438*, 443. [CrossRef]

60. Afshordi, N.; Mohayaee, R.; Bertschinger, E. *Hierarchical phase space structure of dark matter haloes: Tidal debris, caustics, and dark matter annihilation*. *Phys. Rev. D* **2009**, *79*, 083526. [CrossRef]

61. Hahn, O.; Abel, T.; Kaehler, R. *Tracing the Dark Matter sheet in phase space*. *Mon. Not. R. Astron. Soc.* **2012**, *427*, 61. [CrossRef]

62. Hahn, O.; Abel, T.; Kaehler, R. *A new approach to simulating collisionless dark matter fluids*. *Mon. Not. R. Astron. Soc.* **2013**, *434*, 1171. [CrossRef]

63. Angulo, R.E.; Hahn, O. *Large-scale dark matter simulations*. *Living Rev. Comput. Astrophys.* **2022**, *8*, 1. [CrossRef]

64. Nadkarni-Ghosh, S.; Singhal, A. *Phase Space dynamics of triaxial collapse: Joint density-velocity evolution*. *Mon. Not. R. Astron. Soc.* **2016**, *457*, 2773. [CrossRef]
65. Buchert, T.; Domínguez, A.; Pérez-Mercader, J. Extending the scope of models for large-scale structure formation in the Universe. *Astron. Astrophys.* **1999**, *349*, 343.

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

67. Bildhauer, S.; Buchert, T.; Kasai, M. Solutions in Newtonian cosmology—The pancake theory with cosmological constant. *Astron. Astrophys.* **1992**, *263*, 23.

68. Buchert, T.; Götz, G. A class of solutions for self-gravitating dust in Newtonian gravity. *J. Math. Phys.* **1987**, *28*, 2714. [CrossRef]

69. Buchert, T. Lagrangian theory of gravitational instability of Friedmann–Lemaître cosmologies and the ‘Zel’dovich approximation’. *Mon. Not. R. Astron. Soc.* **1992**, *254*, 729. [CrossRef]

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

71. Buchert, T. Lagrangian theory of gravitational instability of Friedmann–Lemaître cosmologies—A generic third-order model for nonlinear clustering. *Mon. Not. R. Astron. Soc.* **1994**, *267*, 811. [CrossRef]

72. Rampf, C.; Buchert, T. Lagrangian perturbations and the matter bispectrum I: Fourth-order model for non-linear clustering. *J. Cosmol. Astropart. Phys.* **2012**, *06*, 021. [CrossRef]

73. Ehlers, J.; Buchert, T. Newtonian cosmology in Lagrangian formulation: Foundations and perturbation theory. *Gen. Rel. Grav.* **1997**, *29*, 733. [CrossRef]

74. Matsubara, T. Recursive solutions of Lagrangian perturbation theory. *Phys. Rev. D* **2015**, *92*, 023534. [CrossRef]

75. Rampf, C. The recursion relation in Lagrangian perturbation theory. *J. Cosmol. Astropart. Phys.* **2012**, *12*, 004. [CrossRef]

76. Rampf, C.; Villone, B.; Frisch, U. How smooth are particle trajectories in a ΛCDM Universe? *Mon. Not. R. Astron. Soc.* **2015**, *452*, 1421. [CrossRef]

77. Zheligovsky, V.; Frisch, U. Time-analyticity of Lagrangian particle trajectories in ideal fluid flow. *J. Fluid Mech.* **2014**, *749*, 404. [CrossRef]

78. Bouchet, F.R.; Juszkiewicz, R.; Colombi, S.; Pellat, R. Weakly nonlinear gravitational instability for arbitrary Omega. *Astrophys. J.* **1992**, *394*, L5. [CrossRef]

79. Bouchet, F.R.; Colombi, S.; Hivon, E.; Juszkiewicz, R. Perturbative Lagrangian approach to gravitational instability. *Astron. Astrophys.* **1995**, *296*, 575.

80. Moutarde, F.; Alimi, J.-M.; Bouchet, F.R.; Pellat, R.; Ramani, A. Precollapse scale invariance in gravitational instability. *Astrophys. J.* **1991**, *382*, 377. [CrossRef]

81. Sahni, V.; Coles, P. Approximation methods for nonlinear gravitational clustering. *Phys. Rep.* **1995**, *262*, 1. [CrossRef]

82. Buchert, T. Lagrangian perturbation theory—A key-model for large-scale structure. *Astron. Astrophys.* **1993**, *267*, L51.

83. Buchert, T. Lagrangian perturbation approach to the formation of large-scale structure. In *Proceedings of the International School of Physics “Enrico Fermi”*; IOS Press: Amsterdam, The Netherlands, 1995; pp. 543–564. [CrossRef]

84. Bertschinger, E. *Cosmological Dynamics*. (Les Houches Summer School Lectures on Cosmology, August 1993); Massachusetts Institute of Technology: Cambridge, MA, USA, 1995. [CrossRef]

85. Bernardeau, F.; Colombi, S.; Gaztanaga, E.; Scoccimarro, R. Large-scale structure of the Universe and cosmological perturbation theory. *Phys. Rep.* **2002**, *367*, 1–248. [CrossRef]

86. Adler, S.; Buchert, T. Lagrangian theory of structure formation in pressure-supported cosmological fluids. *Astron. Astrophys.* **1999**, *343*, 317. [CrossRef]

87. Buchert, T.; Melott, A.L.; Weiß, A. Testing higher-order Lagrangian perturbation theory against numerical simulations—1. pancake models. *Astron. Astrophys.* **1994**, *288*, 349.

88. Coles, P.; Melott, A.L.; Shandarin, S.F. Testing approximations for nonlinear gravitational clustering. *Mon. Not. R. Astron. Soc.* **1993**, *260*, 765. [CrossRef]

89. Melott, A.L.; Buchert, T.; Weiß, A.G. Testing higher-order Lagrangian perturbation theory against numerical simulations—2. hierarchical models. *Astron. Astrophys.* **1995**, *294*, 345.

90. Melott, A.L.; Pellman, T.F.; Shandarin, S.F. Optimizing the Zel'dovich approximation. *Mon. Not. R. Astron. Soc.* **1994**, *269*, 626. [CrossRef]

91. Weiß, A.G.; Gottlöber, S.; Buchert, T. Optimizing higher-order Lagrangian perturbation theory for standard CDM and BSI models. *Mon. Not. R. Astron. Soc.* **1996**, *278*, 953. [CrossRef]

92. Weinberg, D.H.; Gunn, J.E. Large-scale Structure and the Adhesion Approximation. *Mon. Not. R. Astron. Soc.* **1990**, *247*, 260.

93. Michaux, M.; Hahn, O.; Rampf, C.; Angulo, R.E. Accurate initial conditions for cosmological N-body simulations: Minimizing truncation and discreteness errors. *Mon. Not. R. Astron. Soc.* **2021**, *500*, 663. [CrossRef]

94. Schmidt, F. An n-th order Lagrangian forward model for large-scale structure. *J. Cosmol. Astropart. Phys.* **2021**, *04*, 033. [CrossRef]

95. Buchert, T.; Karakatsanis, G.; Klaffl, R.; Schiller, P. The performance of Lagrangian perturbation schemes at high resolution. *Astron. Astrophys.* **1997**, *318*, 1.

96. Yano, T.; Koyama, H.; Buchert, T.; Gouda, N. Universality in the distribution of caustics in the expanding Universe. *Astrophys. J. Suppl.* **2004**, *151*, 185. [CrossRef]

97. Rampf, C. Cosmological Vlasov-Poisson equations for dark matter. *Rev. Mod. Plasma Phys.* **2021**, *5*, 10. [CrossRef]

98. Rampf, C.; Frisch, U. Shell-crossing in quasi-one-dimensional flow. *Mon. Not. R. Astron. Soc.* **2017**, *471*, 671. [CrossRef]

99. Rampf, C.; Hahn, O. Shell-crossing in a ΛCDM Universe. *Mon. Not. R. Astron. Soc.* **2021**, *501*, L71. [CrossRef]
100. Rampf, C.; Frisch, U.; Hahn, O. Unveiling the singular dynamics in the cosmic large-scale structure. *Mon. Not. R. Astron. Soc.* **2021**, *505*, L90. [CrossRef]

101. Ray, S.S.; Frisch, U.; Nazarenko, S.; Matsumoto, T. Resonance phenomenon for the Galerkin-truncated Burgers and Euler equations. *Phys. Rev. E* **2011**, *84*, 016301. [CrossRef] [PubMed]

102. Saga, S.; Taruya, A.; Colombi, S. Lagrangian Cosmological Perturbation Theory at Shell Crossing. *Phys. Rev. Lett.* **2018**, *121*, 241302. [CrossRef] [PubMed]

103. Saga, S.; Taruya, A.; Colombi, S. Cold dark matter protohalo structure around collapse: Lagrangian cosmological perturbation theory versus Vlasov simulations. *Astron. Astrophys.* **2022**, *664*, A3. [CrossRef]

104. Kitaura, F.-S.; Hess, S. Cosmological structure formation with augmented Lagrangian perturbation theory. *Mon. Not. R. Astron. Soc.* **2013**, *435*, L78. [CrossRef]

105. Kitaura, F.-S.; Balaguera-Antolínez, A.; Sinigaglia, F.; Pellejero-Ibàñez, M. The cosmic web connection to the dark matter halo distribution through gravity. *Mon. Not. R. Astron. Soc.* **2022**, *512*, 2245. [CrossRef]

106. Kokron, N.; Chen, S.F.; White, M.; DeRose, J.; Maus, M. Accurate predictions from small boxes: Variance suppression via the Zel'dovich approximation. *J. Cosmol. Astropart. Phys.* **2022**, *09*, 059. [CrossRef]

107. Monaco, P. Approximate methods for the generation of dark matter halo catalogs in the age of precision cosmology. *Galaxies* **2016**, *4*, 53. [CrossRef]

108. Monaco, P.; Sefusatti, E.; Borgani, S.; Crocce, M.; Fosalba, P.; Sheth, R.K.; Theuns, T. An accurate tool for the fast generation of dark matter halo catalogues. *Mon. Not. R. Astron. Soc.* **2013**, *433*, 2389. [CrossRef]

109. Munari, E.; Monaco, P.; Sefusatti, E.; Castorina, E.; Mohammad, F.G.; Anselmi, S.; Borgani, S. Improving fast generation of halo catalogs with higher-order Lagrangian perturbation theory. *Mon. Not. R. Astron. Soc.* **2017**, *465*, 4658. [CrossRef]

110. Buchert, T.; Bartelmann, M. High-spatial resolution in three dimensions: A challenge for large-scale structure formation models. *Astron. Astrophys.* **1991**, *251*, 389.

111. Hidding, J.; Shandarin, S.F.; van de Weygaert, R. The Zel'dovich approximation: Key to understanding cosmic web complexity. *Mon. Not. R. Astron. Soc.* **2014**, *437*, 3442. [CrossRef]

112. Buchert, T.; Shandarin, S.F.; Weiß, A. Building Blocks of Large-Scale Structure in the Universe. Three-Dimensional Singularities and Metamorphoses. 1995, unpublished.

113. Susperregi, M.; Buchert, T. Cosmic density and velocity fields in Lagrangian perturbation theory. *Astron. Astrophys.* **1997**, *323*, 295.

114. Doumler, T.; Hoffman, Y.; Courtois, H.; Gottlöber, S. Reconstructing cosmological initial conditions from galaxy peculiar velocities. I. Reverse Zel'dovich Approximation. *Mon. Not. R. Astron. Soc.* **2013**, *430*, 888. [CrossRef]

115. Brenier, Y.; Frisch, U.; Henon, M.; Loeper, G.; Matarrese, S.; Mohayaee, R.; Sobolevskii, A. Reconstruction of the early Universe as a convex optimization problem. *Mon. Not. R. Astron. Soc.* **2003**, *346*, 501–524. [CrossRef]

116. Nadkarni-Ghosh, S.; Chernoff, D.F. Extending the domain of validity of the Lagrangian approximation. *Mon. Not. R. Astron. Soc.* **2011**, *410*, 1454. [CrossRef]

117. Nadkarni-Ghosh, S.; Chernoff, D.F. Modelling non-linear evolution using Lagrangian Perturbation Theory (LPT) re-expansions. *Mon. Not. R. Astron. Soc.* **2013**, *431*, 799. [CrossRef]

118. Crocce, M.; Scoccimarro, R. Renormalized cosmological perturbation theory. *Phys. Rev. D* **2006**, *73*, 063519. [CrossRef]

119. Matsubara, T. Resumming cosmological perturbations via the Lagrangian picture: One-loop results in real space and in redshift space. *Phys. Rev. D* **2008**, *77*, 063530. [CrossRef]

120. Hivon, E.; Bouchet, F.R.; Colombi, S.; Juszkiewicz, R. Redshift distortions of clustering: A Lagrangian approach. *Astron. Astrophys.* **1995**, *298*, 643.

121. Percival, W.J.; White, M. Testing cosmological structure formation using redshift-space distortions. *Mon. Not. R. Astron. Soc.* **2009**, *393*, 29712. [CrossRef]

122. Desjacques, V.; Jeong, D.; Schmidt, F. Large-scale galaxy bias. *Phys. Rep.* **2018**, *733*, 1–193. [CrossRef]

123. Serrin, I. Mathematical Principles of Classical Fluid Mechanics. In *Encyclopedia of Physics*; Springer: Berlin/Heidelberg, Germany, 1959; Volume VIII.1, pp. 125–263. [CrossRef]

124. Kasai, M. Inhomogeneous cosmological models which are homogeneous and isotropic on average. *Phys. Rev. D* **1993**, *47*, 3214. [CrossRef]

125. Matarrese, S.; Pantano, O.; Saez, D. General relativistic dynamics of irrotational dust: Cosmological implications. *Phys. Rev. Lett.* **1994**, *72*, 320. [CrossRef]

126. Matarrese, S.; Pantano, O.; Saez, D. A Relativistic Approach to Gravitational Instability in the Expanding Universe—Second-Order Lagrangian Solutions. *Mon. Not. R. Astron. Soc.* **1994**, *271*, 513. [CrossRef]

127. Matarrese, S.; Terranova, D. Post-Newtonian cosmological dynamics in Lagrangian coordinates. *Mon. Not. R. Astron. Soc.* **1996**, *283*, 400. [CrossRef]

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

129. Russ, H.; Soffel, M.H.; Kasai, M.; Börner, G. Age of the universe: Influence of the inhomogeneities on the global expansion factor. *Phys. Rev. D* **1997**, *56*, 2044. [CrossRef]

130. Parry, J.; Salopek, D.S.; Stewart, J.M. Solving the Hamilton-Jacobi equation for general relativity. *Phys. Rev. D* **1994**, *49*, 2872. [CrossRef]
131. Salopek, D.S.; Stewart, J.M. Hamilton-Jacobi theory for general relativity with matter fields. *Class. Quant. Grav.* **1992**, *9*, 1943.
[CrossRef]

132. Rampf, C.; Rigopoulos, G. Initial conditions for cold dark matter particles and general relativity. *Phys. Rev. D* **2013**, *87*, 123525.
[CrossRef]

133. Rigopoulos, G.; Valkenburg, W. Nonlinear approximation for perturbations in $\Lambda$CDM. *Phys. Rev. D* **2012**, *86*, 043523. [CrossRef]

134. Buchert, T.; Mädler, T. Editorial Note to: On the Newtonian Limit of Einstein's Theory of Gravitation (by Jürgen Ehlers). *Gen. Rel. Grav.* **2019**, *51*, 162. [CrossRef]

135. Vigneron, Q. The 1+3-Newton-Cartan system and Newton-Cartan cosmology. *Phys. Rev. D* **2021**, *103*, 064064. [CrossRef]

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

137. Buchert, T.; Nayet, C.; Wiegand, A. Lagrangian theory of structure formation in relativistic cosmology. II. Average properties of a generic evolution model. *Phys. Rev. D* **2013**, *87*, 123503. [CrossRef]

138. Buchert, T. Dark Energy from structure: A status report. *Gen. Rel. Grav.* **2008**, *40*, 467. [CrossRef]

139. Buchert, T. Toward physical cosmology: Focus on inhomogeneous geometry and its non-perturbative effects. *Class. Quant. Grav.* **2011**, *28*, 164007. [CrossRef]

140. Buchert, T.; Räsänen, S. Backreaction in Late-Time Cosmology. *Annu. Rev. Nucl. Part. Sci.* **2012**, *62*, 57. [CrossRef]

141. Buchert, T.; Larena, J.; Alimi, J.M. Correspondence between kinematical backreaction and scalar field cosmologies—The ‘morphon field’. *Class. Quant. Grav.* **2006**, *23*, 6379. [CrossRef]

142. Buchert, T.; Obadia, N. Effective inhomogeneous inflation: Curvature inhomogeneities of the Einstein vacuum. *Class. Quant. Grav.* **2011**, *28*, 162002. [CrossRef]

143. Vigneron, Q.; Buchert, T. Dark Matter from Backreaction? Collapse models on galaxy cluster scales. *Class. Quant. Grav.* **2019**, *36*, 175006. [CrossRef]

144. Krasiński, A. *Inhomogeneous Cosmological Models*; Cambridge University Press: Cambridge, UK, 1997; ISBN 0521481805.

145. Lemaître, G. The expanding universe. *Ann. Soc. Sci. Brux. A* **1933**, *53*, 51. [CrossRef]

146. Plebański, J.; Krasiński, A. *An Introduction to General Relativity and Cosmology*; Cambridge University Press: Cambridge, UK, 2006; ISBN 0-521-85623-X.

147. Tolman, R.C. Effect of inhomogeneity on cosmological models. *Proc. Natl. Acad. Sci. USA* **1934**, *20*, 169. [CrossRef]

148. Morita, M.; Nakamura, K.; Kasai, M. Relativistic Zel'dovich approximation in a spherically symmetric model. *Phys. Rev. D* **1998**, *57*, 6094. [CrossRef]

149. Desgrange, C.; Heinesen, A.; Buchert, T. Dynamical spatial curvature as a fit to type Ia supernovae. *Int. J. Mod. Phys. D* **2019**, *28*, 1950143. [CrossRef]

150. Heinesen, A.; Buchert, T. Solving the curvature and Hubble parameter inconsistencies through structure formation-induced curvature. *Class. Quant. Grav.* **2020**, *37*, 164001. [CrossRef]

151. Roukema, B.F.; Ostrowski, J.J.; Borkowska, J. inhomog: Biscale Kinematical Backreaction Analytical Evolution. Astrophysics Source Code Library. 2017. Available online: https://ui.adsabs.harvard.edu/abs/2017ascl.soft11002R (accessed on 26 September 2022).

152. Buchert, T.; Kerscher, M.; Sicka, C. Backreaction of inhomogeneities on the expansion: The evolution of cosmological parameters. *Phys. Rev. D* **2000**, *62*, 043525. [CrossRef]

153. Al Roumi, F.; Buchert, T.; Wiegand, A. Lagrangian theory of structure formation in relativistic cosmology. IV. Lagrangian approach to gravitational waves. *Phys. Rev. D* **2017**, *96*, 123538. [CrossRef]

154. Li, Y.Z.; Mourier, P.; Buchert, T.; Wiltshire, D.L. Lagrangian theory of structure formation in relativistic cosmology. V. Irrotational fluids. *Phys. Rev. D* **2018**, *98*, 043507. [CrossRef]

155. Delgado Gaspar, I.; Buchert, T. Lagrangian theory of structure formation in relativistic cosmology. VI. Comparison with Szekeres exact solutions. *Phys. Rev. D* **2021**, *103*, 023513. [CrossRef]

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

157. Buchert, T.; Carfora, M. On the curvature of the present-day Universe. *Class. Quant. Grav.* **2008**, *25*, 195001. [CrossRef]

158. Szekeres, P. Quasispherical Gravitational Collapse. *Phys. Rev. D* **1975**, *12*, 2941. [CrossRef]

159. Szafron, D.A. Inhomogeneous cosmologies: New exact solutions and their evolution. *J. Math. Phys.* **1977**, *18*, 1673. [CrossRef]

160. Spero, A.; Szafron, D.A. Spatial conformal flatness in homogeneous and inhomogeneous cosmologies. *J. Math. Phys.* **1978**, *19*, 1536. [CrossRef]

161. Delgado Gaspar, I.; Sussman, R.A.; McNutt, D.D.; Coley, A.A. Comment on 'Szekeres universes with homogeneous scalar fields'. *Eur. Phys. J. C* **2021**, *81*, 310. [CrossRef]

162. Bonnor, W.B.; Sulaiman, A.H.; Tomimura, N. Szekeres's Space-Times Have No Killing Vectors. *Gen. Rel. Grav.* **1977**, *8*, 549. [CrossRef]

163. Bolejko, K. Structure formation in quasispherical Szekeres model. *Phys. Rev. D* **2006**, *73*, 123508. [CrossRef]

164. Bolejko, K. Evolution of cosmic structures in different environments in the quasispherical Szekeres model. *Phys. Rev. D* **2007**, *75*, 043508. [CrossRef]

165. Bolejko, K. The Szekeres Swiss Cheese model and the CMB observations. *Gen. Rel. Grav.* **2009**, *41*, 1737. [CrossRef]
166. Bolejko, K.; Célérier, M.N. Szekeres Swiss-Cheese model and supernova observations. *Phys. Rev. D* **2010**, *82*, 103510. [CrossRef]

167. Bolejko, K.; Sussman, R.A. Cosmic spherical void via coarse-graining and averaging non-spherical structures. *Phys. Lett. B* **2011**, *697*, 265. [CrossRef]

168. Buckley, R.G.; Schlegel, E.M. CMB dipoles and other low-order multipoles in the quasispherical Szekeres model. *Phys. Rev. D* **2013**, *87*, 023524. [CrossRef]

169. Coley, A.A.; Layden, N.; McNutt, D. An invariant characterization of the quasi-spherical Szekeres dust models. *Gen. Rel. Grav.* **2019**, *51*, 164. [CrossRef]

170. Delgado Gaspar, I.; Hidalgo, J.C.; Sussman, R.A.; Quiros, I. Black hole formation from the gravitational collapse of a nonspherical network of structures. *Phys. Rev. D* **2018**, *97*, 104029. [CrossRef]

171. Hellaby, C. The null and KS limits of the Szekeres model. *Class. Quant. Grav.* **1996**, *13*, 2537. [CrossRef]

172. Hellaby, C.; Krasiński, A. You can't get through Szekeres wormholes: Or, regularity, topology and causality in quasispherical Szekeres models. *Phys. Rev. D* **2002**, *66*, 084011. [CrossRef]

173. Hellaby, C.; Krasiński, A. Physical and geometrical interpretation of the $\epsilon \le 0$ Szekeres models. *Phys. Rev. D* **2008**, *77*, 023529. [CrossRef]

174. Ishak, M.; Richardson, J.; Garred, D.; Whittington, D.; Nwankwo, A.; Sussman, R.A. Dark energy or apparent acceleration due to a relativistic cosmological model more complex than the Friedmann-Lemaître-Robertson Walker model? *Phys. Rev. D* **2011**, *84*, 089902. [CrossRef]

175. Koksbang, S.M.; Hannestad, S. Methods for studying the accuracy of light propagation in N-body simulations. *Phys. Rev. D* **2015**, *91*, 043508. [CrossRef]

176. Koksbang, S.M.; Hannestad, S. Studying the precision of ray tracing techniques with Szekeres models. *Phys. Rev. D* **2015**, *92*, 023532. [CrossRef]

177. Koksbang, S.M. Light propagation in Swiss-cheese models of random close-packed Szekeres structures: Effects of anisotropy and comparisons with perturbative results. *Phys. Rev. D* **2017**, *95*, 063532. [CrossRef]

178. Sussman, R.A.; Bolejko, K. A Novel approach to the dynamics of Szekeres dust models. *Class. Quant. Grav.* **2012**, *29*, 065018. [CrossRef]

179. Sussman, R.A.; Delgado Gaspar, I. Multiple nonspherical structures from the extrema of Szekeres scalars. *Phys. Rev. D* **2015**, *92*, 083533. [CrossRef]

180. Sussman, R.A.; Delgado Gaspar, I.; Hidalgo, J.C. Coarse-grained description of cosmic structure from Szekeres models. *J. Cosmol. Astropart. Phys.* **2016**, *03*, 012. [CrossRef]

181. Vrba, D.; Svítek, O. Modelling inhomogeneity in Szekeres spacetime. *Gen. Rel. Grav.* **2014**, *46*, 1808. [CrossRef]

182. Walters, A.; Hellaby, C. Constructing realistic Szekeres models from initial and final data. *Phys. Rev. D* **2012**, *121*, 001. [CrossRef]

183. Goode, S.W.; Wainwright, J. Singularities and evolution of the Szekeres cosmological models. *Phys. Rev. D* **1982**, *26*, 3315. [CrossRef]

184. Meures, N.; Bruni, M. Exact nonlinear inhomogeneities in $\Lambda$CDM cosmology. *Phys. Rev. D* **2011**, *83*, 123519. [CrossRef]

185. Meures, N.; Bruni, M. Redshift and distances in a $\Lambda$CDM cosmology with non-linear inhomogeneities. *Mon. Not. R. Astron. Soc.* **2012**, *419*, 3. [CrossRef]

186. Bolejko, K. Relativistic numerical cosmology with Silent Universes. *Class. Quant. Grav.* **2017**, *35*, 024003. [CrossRef]

187. Bruni, M.; Matarrese, S.; Pantano, O. Dynamics of silent universes. *Astrophys. J.* **1995**, *445*, 958. [CrossRef]

188. van Elst, H.; Ellis, G.F.R. The Covariant approach to LRS perfect fluid space-time geometries. *Class. Quant. Grav.* **1996**, *13*, 1099. [CrossRef]

189. Delgado Gaspar, I.; Buchert, T.; Ostrowski, J.J. Beyond relativistic Lagrangian perturbation theory. I. An exact-solution controlled model for structure formation. *Phys. Rev. D* **2022**, submitted.

190. Wiegand, A.; Buchert, T. Multiscale cosmology and structure-emerging Dark Energy: A plausibility analysis. *Phys. Rev. D* **2010**, *82*, 023523. [CrossRef]