###### 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 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 exact average model,
Ricci-flat LTB models. This 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.

8 1 583 Academic Editors: Luigi Tedesco and Giovanni Marozzi 27 September 2022
27 October 2022 4 November 2022 https://doi.org/10.3390/universe8110583 On
general-relativistic Lagrangian perturbation theory and its non-perturbative
generalization Relativistic Lagrangian perturbations Thomas Buchert 1,‡*,
Ismael Delgado Gaspar 2,1,‡ and Jan J. Ostrowski 2 Thomas Buchert, Ismael
Delgado Gaspar and Jan J. Ostrowski Buchert, T.; Delgado Gaspar, I.;
Ostrowski, J.J. Correspondence: TB These authors contributed equally to this
work.

##  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. Brustein95 , Gasperini1997 ) 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 lifshitz46  (see also lifshitzkhalatnikov , hawking , tomita  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 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 mukhanovetal
; kodamasasaki ; MalikWands2009 ; Durrer . 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 as intrinsic geometrical features of the physical
manifold. Such a set-up 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
Ehlers61 . In this context, the background spacetime defines an arbitrary,
absolute reference frame. In the seminal paper Bardeen1980 , gauge-invariant
variables where defined. However, although the basic physical perturbation
variable is the density contrast, none of the Bardeen’s variables actually
represented this quantity. This issue was partially resolved by John Stewart
stewart . A different (although fully compatible) approach was presented in
EllisBruni1989 , 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 define 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 stewart . 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 foliations ; buchert:generalfluid . Gauge choices in
the two-manifold approaches can be mapped to foliation choices in a covariant
setting, see vitentietal  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 generalbackground .

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 combi . 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 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 buchert:average ; EllisBuchert2005 . This also holds for operations
such as Ricci-flow smoothing and scaling in a 3-manifold buchertcarfora1  and
for other coarse-graining techniques such as renormalization
carfora:renormalization . 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 buchert:dust ;
buchert:fluid ; buchert:generalfluid , $2+1+1$ foliations on the lightcone
buchert:lightcone , and their manifestly 4-covariant formulations Gasperini2 ;
Gasperini3 , see also foliations ; covariance . 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 zeldovichmyshkis:book . Also, Zel’dovich might have
been aware of related work in plasma physics such as Dawson’s discussion of
multi-streaming in a cold plasma Dawson1959 . 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.
peebles1:book ) and the powerful nonlinear solution for inertial systems. His
early papers zeldovich70a ; zeldovich70b ; zeldovich78 , zeldovichshandarin ,
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, later, was
‘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 doroshkevichetal73  (among others, see
e.g. the considerations on the density probability distribution in
doroshkevich70 , grinsteinwise8789 , bartelmannschneider92 , investigations on
self-similar solutions fillmoregoldreich  and trajectories bartelmann .)
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
buchert89 ). 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
shandarinzeldovich89 ), but confirming the above picture, and reinforcing what
was earlier found by Lin, Mestel and Shu linmestelshu  (see also goodmanbinney
, ZAaccuracy ).

Further numerical work within Zel’dovich’s group doroshkevichetal80 ,
klypinshandarin , 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 80’s
melottshandarin89 . 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 $\Lambda$CDM standard model of cosmology performs in between
these two scenarii.

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 arnoldshandarinzeldovich , in which the Lagrange-singularity
theory arnold:book  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 shandarin80 , and the
derivation of exact plane-symmetric solutions for self-gravitating dust flows
zentsovachernin , confirming what was suggested by doroshkevichetal73 111For
further discussion of and references to the singularity problem see Ref. rza1
, Section VB.. The effect of velocity dispersion was discussed by Shukurov
shukurov . Thereafter, an extension of Zel’dovich’s approximation to
phenomenologically access the multi-stream regime was proposed in gurbatov89
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 adhesion ; see
the review adhesive  that also investigates deviations from mean field
gravity.222In 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
bertschinger:phasespace . Since the phase space is six-dimensional this,
however, provides limits of resolution. For recent simulations see Abel1 ;
Abel2  and the review Hahn:review . For a discussion of eigenvalues of the
deformation tensor, see sharvari3 . A further extension was suggested in
extending  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
nonperturbative .

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 doroshkevichetal73  and the plane-symmetric investigation
in zentsovachernin ) 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 buchertgoetz , buchert89 ,
bildhaueretal , and, finally, the systematic framework of a Lagrangian
perturbation theory, in which Zel’dovich’s approximation forms a subclass of
the first-order solutions buchert89 , buchert92 ; closed form solutions at
second buchertehlers , third buchert94  and fourth order rampfbuchert  can be
reconstructed using the n-th order perturbation and solution schemes
investigated in ehlersbuchert  (see also more recent work on recursion
relations matsubara2 ; rampf:recursion  and related investigations of time-
analyticity or smoothness of trajectories Zheligovsky2014 , Rampfetal2015 ).
The early 90’s can be considered as the birth of Lagrangian perturbation
theory, when other groups started to contribute moutarde , bouchetetal92 ;
bouchetetal95  (see the early review papers shandarinzeldovich89 , sahnicoles
, ehlersbuchert , a brief summary buchert93L , the tutorial buchert_varenna ,
the more general review bertschinger:review , and the more recent review
bernardeau:review ). 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 gurbatov89 ; adhesion .
In a Lagrangian perturbation approach this Laplacian forcing must be
transformed to Lagrangian coordinates, which introduces substantial
nonlinearity to the problem adlerbuchert .

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
bsw .

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 tza ; melottpellmanshandarin
,buchertmelottweiss ; melottbuchertweiss ; weissetal96 ; for the adhesion
approximation gurbatov89 , see weinberggunn ; for recent direct comparisons
between higher-order Lagrangian perturbation solutions and numerical
simulations, see Michaux2021 ; Schmidt2021 . Shell-crossing itself was also in
the focus of interest buchertehlers ; buchert_performance ; yanoetal , a
sustained interest up to the present day, see e.g. the recent papers
rampffrisch ; colombi:shellcrossing1 ; colombi:shellcrossing2 ;
rampf:shellcrossing1 ; rampf:shellcrossing2 ; rampf:vlasov ; tyger . 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
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. monaco1 ; monacorev ; monaco2 , kitaura1 ;
kitaura2  and ZAnumerics  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 buchert89 , further insights were
obtained for three-dimensional special and generic initial data
buchertbartelmann ; buchert_performance  including hints for a mechanism of
gravitational fragmentation at higher orders in Lagrangian perturbation
theory.333 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. 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.444Recently, using newer
numerical methods, hiddingetal  provided a detailed generalization of the 2D
pancakes of arnoldshandarinzeldovich .

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 millenium: the inverse Lagrangian perturbation solutions
have been employed for reconstruction techniques, see recon  and references
therein to earlier work, courtois:reconstruction  and follow-up papers of
these authors, as well as related techniques frisch:reconstruction . Work is
pursued, e.g. for convergence issues sharvari1 ; rampf:shellcrossing1 , re-
expansions sharvari2 , renormalization, re-summation and recursion techniques
crocce , matsubara1 , matsubara2 ; rampf:recursion , respectively, redshift
space distorsions hivon , white  and references therein, and statistical
biasing relations between the galaxy and matter distributions biasing  (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, ${\bf
x}=a(t){\bf q}$, with $a(t)$ a solution in the family of FLRW cosmologies;
${\bf x}$ are inertial non-rotating coordinates in Euclidean space. On the
contrary, the Lagrangian coordinate system, denoted by $\bf 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 serrin , for a systematic introduction for
self-gravitating Newtonian flows, see ehlersbuchert . For many complementary
investigations of relativistic Lagrangian perturbations see, e.g. ornella1 ;
ornella2 , terranova , kasai93 ; russ1 ; russ2 , including gradient expansions
gradientexpansion1 ; gradientexpansion2  and follow-up papers by these
authors, as well as the recent papers on gradient expansions Rigopoulos2012
and Rampf2013 .

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

$\mathbf{x}:=\mathbf{f}(\mathbf{X},t)\;;\;\mathbf{v}:=\dot{\mathbf{f}}(\mathbf{X},t)\;;\;\mathbf{a}:=\ddot{\mathbf{f}}(\mathbf{X},t)\;;\;\varrho(\mathbf{X},t)=\frac{\mathring{\varrho}}{J(\mathbf{X},t)}\;;\;\bm{\omega}=\frac{\mathring{\bm{\omega}}\cdot\bm{\nabla}_{\bf
0}\mathbf{f}}{J(\mathbf{X},t)}\;\,;\;etc.\,,$ (1)

where $\mathring{\varrho}$ and $\mathring{\bm{\omega}}$ stand for the initial
values of the density and vorticity fields, $\bm{\nabla}_{\bf 0}$ for the
nabla operator with respect to the Lagrangian coordinates, and $J$ is the
determinant of the Jacobian matrix $(f^{a}_{\;|j})$,

$J:=\det(f^{a}_{\;|j})(\mathbf{X},t)=\frac{1}{6}\epsilon_{ijk}\epsilon^{klm}f^{i}_{\
|k}f^{j}_{\ |l}f^{k}_{\ |m}\>\>,$ (2)

where we have denoted by a vertical slash |i the spatial derivative with
respect to $X_{i}$, by an overdot the partial time derivative,
$\epsilon_{ijk}$ being the Levi-Cività symbol.555We have distinguished counter
indices $a,b,c,\cdots$, here counting vector components, from coordinate
indices $i,j,k,\cdots$. 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. The remaining field equations of the
Euler-Newton system, $\bm{\nabla}\times\mathbf{g}=\mathbf{0}\
;\bm{\nabla}\cdot\mathbf{g}=\Lambda-4\pi G\varrho$, become evolution equations
in the Lagrangian picture, comprising the closed Lagrange-Newton system:

$\displaystyle\delta_{ab}\ddot{f}^{a}_{\;\,|[i}\,{f}^{b}_{\;\,|j]}=0\quad;\quad\frac{1}{2}{\epsilon}_{abc}{\epsilon}^{ikl}\ddot{f}^{a}_{\;\,|i}\,f^{b}_{\;\,|k}\,f^{c}_{\;\,|l}=\Lambda
J-4\pi G\mathring{\varrho}\;.$ (3)

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

A Lagrangian perturbation theory results from the ansatz that represents ${\bf
d}f^{a}$ in terms of a homogeneous-isotropic deformation and expands the
Lagrange-Newton system to a given order in the deviation gradient field ${\bf
d}{\bf P}$: ${\bf d}{\bf f}=a(t)({\bf X}+{\bf d}{\bf P})$.

####  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 $\bf g$ with coefficients $g_{ij}$ in an exact (co-tangential) basis
${\bf d}X^{i}\otimes{\bf d}X^{j}$. We can write any such metric as a quadratic
form of deformation one-forms or co-frame fields, ${\bf
g}=\delta_{ab}\,\bm{\eta}^{a}\otimes\bm{\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,

$\bm{\eta}^{a}\mapsto\bm{\mathrm{d}}f^{a}\ .$ (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 bm . 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. QuentinNC  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 rza1 . 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}{\mathrm{d}}x^{i}{\mathrm{d}}x^{j}=\delta_{ab}f^{a}_{\;\,|i}f^{b}_{\;\,|j}{\mathrm{d}}X^{i}{\mathrm{d}}X^{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, rza1 , 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:

$\displaystyle\Theta_{[ij]}=0\;;$ (5a)
$\displaystyle{\dot{g}}_{ij}=2\,g_{ik}\Theta^{k}_{\;j}\quad;\quad{\dot{\Theta}}^{i}_{\;j}+\theta\Theta^{i}_{\;j}=-{\cal
R}^{i}_{\;j}+(\,4\pi G\varrho+\Lambda\,)\,\delta^{i}_{\;j}\;;$ (5b)
$\displaystyle\theta^{2}-\Theta^{i}_{\;j}\,\Theta^{j}_{\;i}=-{\cal R}+(16\pi
G\varrho+2\Lambda)\;;$ (5c)
$\displaystyle\Theta^{i}_{\;j||i}-\theta_{|j}=0\;,$ (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-Cività; the dust density is given by
$\varrho=\mathring{\varrho}J^{-1}$, $J\geq 0$; $J$ defines the coefficient
function of the $3-$volume form, $J\equiv{\sqrt{g}}/{\sqrt{G}}$, with
${\sqrt{g}}\ {\bf 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_{\rm 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:666Starting 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 rza2 ;
rza3 ; for the Euclidean restriction (integrable limit) of Gram’s matrix see
[appendix A]rza3 ), and for the integrable limit of the gravito-magnetic part
see [Sect. III.A.2]rza1 . Note also that Equations (6a) imply
$G_{ab}\,\ddot{\eta}^{a}_{[i}\eta^{b}_{\ j]}=0\ $.

$\displaystyle G_{ab}\,\dot{\eta}^{a}_{[i}\eta^{b}_{\ j]}=0\;;$ (6a)
$\displaystyle\frac{1}{2J}\epsilon_{abc}\epsilon^{ikl}\left(\dot{\eta}^{a}_{\
j}\eta^{b}_{\ k}\eta^{c}_{\ l}\right)\dot{}=-{\mathcal{R}}^{i}_{\
j}+\left(4\pi G\varrho+\Lambda\right)\delta^{i}_{\ j}\;;$ (6b)
$\displaystyle\frac{1}{J}\epsilon_{abc}\epsilon^{mjk}\dot{\eta}^{a}_{\
m}\dot{\eta}^{b}_{\ j}\eta^{c}_{\ k}=-{\mathcal{R}}+\left(16\pi
G\varrho+2\Lambda\right)\;;$ (6c)
$\displaystyle\left(\tfrac{1}{J}\epsilon_{abc}\epsilon^{ikl}\dot{\eta}^{a}_{\
j}\eta^{b}_{\ k}\eta^{c}_{\
l}\right)_{||i}=\left(\tfrac{1}{J}\epsilon_{abc}\epsilon^{ikl}\dot{\eta}^{a}_{\
i}\eta^{b}_{\ k}\eta^{c}_{\ l}\right)_{|j}\;.$ (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}\,\ddot{\eta}^{a}_{[i}\eta^{b}_{\
j]}=0\quad;\quad\frac{1}{2J}\epsilon_{abc}\epsilon^{ik\ell}\ddot{\eta}^{a}_{\
i}\eta^{b}_{\ k}\eta^{c}_{\ \ell}=\Lambda-4\pi G{\varrho}\;.$ (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
rza1 .

We only briefly mention Paper III rza3  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 ehlersbuchert . We discuss more in detail the
results of Paper II rza2  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_{\cal{D}}:\propto
V_{\cal{D}}^{1/3}$, suitably normalized, that appear to be sourced by an
effective perfect fluid energy-momentum tensor buchert:dust ; buchert:fluid ,
reviews are buchert:review ; buchert:focus ; buchertrasanen :

$\displaystyle 3\frac{{\ddot{a}}_{{\cal{D}}}}{a_{{\cal{D}}}}=-4\pi
G(\varrho_{\rm eff}^{{\cal{D}}}+3{p}_{\rm eff}^{{\cal{D}}})+\Lambda\;;\;$
$\displaystyle
3H_{{\cal{D}}}^{2}-\frac{3k_{{\cal{D}}}}{a_{{\cal{D}}}^{2}}=8\pi G\varrho_{\rm
eff}^{{\cal{D}}}+\Lambda\;;\;$ $\displaystyle{\dot{\varrho}}_{\rm
eff}^{{\cal{D}}}+3H_{{\cal{D}}}(\varrho_{\rm eff}^{{\cal{D}}}+{p}_{\rm
eff}^{{\cal{D}}})=0\,,$ (8a) where the effective densities are defined as
$\displaystyle\varrho_{{\rm
eff}}^{{{\cal{D}}}}:=\langle{\varrho}\rangle_{\cal{D}}+\varrho_{\Phi}\;\;\;;\;\;\varrho_{\Phi}:=-\frac{1}{16\pi
G}{{\cal{Q}}}_{{{\cal{D}}}}-\frac{1}{16\pi G}{\cal{W}}_{\cal{D}}\ ;$
$\displaystyle{p}_{{\rm
eff}}^{{{\cal{D}}}}:=p_{\Phi}\;\;\;;\;\;\,p_{\Phi}:=-\frac{1}{16\pi
G}{{\cal{Q}}}_{{{\cal{D}}}}+\frac{1}{48\pi G}{\cal{W}}_{\cal{D}}\ .$ (8b)

In this form the effective equations suggest themselves to interpret the extra
fluctuating sources in terms of a scalar field (the morphon field morphon ),
which refers to the inhomogeneities in the geometric functionals
${\cal{W}}_{\cal{D}}$ and ${\cal{Q}}_{\cal{D}}$.777This 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 morphoninflation , or
models mimicking dark matter quentinDM . The variable
${\cal{W}}_{\cal{D}}:=\langle{\cal
R}\rangle_{\cal{D}}-6k_{{\cal{D}}}/a_{{\cal{D}}}^{2}$ describes the deviations
of the average Ricci scalar curvature from a Friedmannian (eventually scale-
dependent) constant curvature $k_{\cal{D}}$ on the domain ${\cal{D}}$, and the
variable ${\cal{Q}}_{\cal{D}}$ is the backreaction functional, made up of
averaged principal scalar invariants of the expansion tensor:

${\cal Q}_{\cal D}\;=\;2\langle{\mathrm{II}}\rangle_{{\cal
D}}-\frac{2}{3}\langle{\mathrm{I}}\rangle_{{\cal D}}^{2}\ .$ (9)

Let us look at two examples. First, we look at the LTB (Lemaître-Tolman-Bondi)
solution Lemaitre:1933gd ; tolman1934effect ; kras1 ; kras2 . 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) rzaltb .

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 rza2 :

$\left\langle\mathrm{I}(\Theta^{i}_{\;j})\right\rangle_{\rm
LTB}=\frac{4\pi}{V_{LTB}}\int_{0}^{r_{\cal
D}}\frac{\partial_{r}\left(\dot{R}R^{2}\right)}{\sqrt{1+2E}}{\hbox{d}}r\;;$
(10) $\left\langle\mathrm{II}(\Theta^{i}_{\;j})\right\rangle_{\rm
LTB}=\frac{4\pi}{V_{LTB}}\int_{0}^{r_{\cal
D}}\frac{\partial_{r}\left(\dot{R}^{2}R\right)}{\sqrt{1+2E}}{\hbox{d}}r\;;$
(11) $\left\langle\mathrm{III}(\Theta^{i}_{\;j})\right\rangle_{\rm
LTB}=\frac{4\pi}{3V_{LTB}}\int_{0}^{r_{\cal
D}}\frac{\partial_{r}\left(\dot{R}^{3}\right)}{\sqrt{1+2E}}{\hbox{d}}r\;,$
(12)

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

$V_{LTB}=\frac{4\pi}{3}\int_{0}^{r_{\cal{D}}}\frac{\partial_{r}\left(R^{3}\right)}{\sqrt{1+2E}}{\hbox{d}}r\;.$
(13)

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

$\left\langle\mathrm{II}(\Theta^{i}_{\;j})\right\rangle_{\rm
LTB}=\displaystyle\frac{1}{3}\left\langle\mathrm{I}(\Theta^{i}_{\;j})\right\rangle_{\rm
LTB}^{2}\;\;;\;\;\left\langle\mathrm{III}(\Theta^{i}_{\;j})\right\rangle_{\rm
LTB}=\displaystyle\frac{1}{27}\left\langle\mathrm{I}(\Theta^{i}_{\;j})\right\rangle_{\rm
LTB}^{3}\,\;.$ (14)

This proves that, for a spatially Ricci-flat LTB model, the backreaction term
vanishes.

Second, we calculated in rza2  the backreaction model provided by the
Relativistic Zel’dovich approximation (RZA), i.e. the first-order solution of
the Lagrange-Einstein system:

${}^{\rm RZA}{\cal Q}_{\cal D}\;=$
$\displaystyle\displaystyle\frac{\dot{\xi}^{2}\left(\gamma_{1}+\xi\gamma_{2}+\xi^{2}\gamma_{3}\right)}{\left(1+\xi\langle{\rm
I}_{\rm\bf i}\rangle_{{\cal C_{D}}}+\xi^{2}\langle{\rm II}_{\rm\bf
i}\rangle_{{\cal C_{D}}}+\xi^{3}\langle{\rm III}_{\rm\bf i}\rangle_{{\cal
C_{D}}}\right)^{2}}\;\;,$ (15a) where we have defined the set of initial data
featuring the initial principal scalar invariants of the expansion tensor,
their initial average $\langle...\rangle_{\cal
C_{D}}:=(1/V_{{\cal{D}}}(t_{i}))\int_{{\cal{D}}(t_{i})}...{\rm
J}(X^{i},t_{i}){\mathrm{d}}X^{i}$, with ${\rm J}(X^{i},t_{i})=1$, (the first
term is the initial backreaction term): $\displaystyle\gamma_{1}:=2\langle{\rm
II}_{\rm\bf i}\rangle_{{\cal C_{D}}}-\frac{2}{3}\langle{\rm I}_{\rm\bf
i}\rangle_{{\cal C_{D}}}^{2}={\cal{Q}}^{\rm initial}_{{\cal C}_{D}}\;;$
$\displaystyle\gamma_{2}:=6\langle{\rm III}_{\rm\bf i}\rangle_{{\cal
C_{D}}}-\frac{2}{3}\langle{\rm II}_{\rm\bf i}\rangle_{{\cal C_{D}}}\langle{\rm
I}_{\rm\bf i}\rangle_{{\cal C_{D}}}\;;$ $\displaystyle\gamma_{3}:=2\langle{\rm
I}_{\rm\bf i}\rangle_{{\cal C_{D}}}\langle{\rm III}_{\rm\bf i}\rangle_{{\cal
C_{D}}}-\frac{2}{3}\langle{\rm II}_{\rm\bf i}\rangle_{{\cal C_{D}}}^{2}\;.$
(15b)

It is interesting to point out the leading, large-scale behaviour of this
functional: taking as an example an Einstein-de Sitter background (with
vanishing curvature and vanishing $\Lambda$), the leading term is ${}^{\rm
RZA}{\cal Q}_{\cal D}\approx(1/a){\cal{Q}}^{\rm initial}_{{\cal C}_{D}}$. It
corresponds to an approximation of the exact scaling solution of the averaged
system, ${\cal Q}_{\cal D}=(1/a_{\cal{D}}){\cal{Q}}^{\rm initial}_{{\cal
C}_{D}}$, leading to the same scaling behaviour for the average spatial scalar
curvature $\langle{\cal R}\rangle_{\cal{D}}$, via the following coupling
equation (following from the last equation of Eq. (8a)):

$0=\frac{1}{a_{\cal{D}}^{6}}\partial_{t}\left({\cal{Q}}_{\cal{D}}a_{\cal{D}}^{6}\right)+\frac{1}{a_{\cal{D}}^{2}}\partial_{t}\left(\langle{\cal
R}\rangle_{\cal{D}}a_{\cal{D}}^{2}\right)\ .$ (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 asta1 ; asta2 .

The backreaction functional (15) was implemented numerically in the GPL-
licensed code inhomog . The inhomog library performs a numerical integration
of the averaged Raychaudhuri 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 buchert:focus .) A corresponding Newtonian derivation for the ZA
model is provided in bks  and it arises from (15a) and (15) 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.

In Newtonian cosmology, the backreaction functional
${\cal{Q}}_{{\mathbb{T}}^{3}}$ also vanishes for the deviations off 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 buchert94 . 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
buchert:average .

Paper IV rza4  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 tight 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 here give 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:

$\displaystyle E_{ki|k}$ $\displaystyle=$
$\displaystyle-\tfrac{2}{3a^{3}}\,W_{\>|i}\>;$ (17a) $\displaystyle H_{ki|k}$
$\displaystyle=$ $\displaystyle 0\>;$ (17b)
$\displaystyle\dot{E}_{ij}+3H{E}_{ij}-\tfrac{1}{a}\epsilon_{(i}{}^{\,\,kl}H_{j)l|k}$
$\displaystyle=$ $\displaystyle-4\pi G{\varrho}_{H}\,\sigma_{ij}\>;\quad$
(17c)
$\displaystyle\dot{H}_{ij}+3H{H}_{ij}+\tfrac{1}{a}\epsilon_{(i}{}^{\,\,kl}E_{j)l|k}$
$\displaystyle=$ $\displaystyle 0\>,$ (17d)

where ${\varrho}_{H}=({\varrho}_{H})_{\rm\bf i}\,/\,a^{3}$ denotes the
homogeneous part of the density field, $W=-4\pi G\varrho_{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:

$\begin{split}\Box_{0}E_{ij}-7H\dot{E}_{ij}-4(4\pi
G{\varrho}_{H}+\Lambda)E_{ij}=-\frac{1}{a^{5}}\mathcal{D}_{ij}W+4\pi
GH{\varrho}_{H}\,\sigma_{ij}\ ;\end{split}$ (18a)
$\begin{split}\Box_{0}{H}_{ij}-7H\dot{H}_{ij}-4(5\pi
G{\varrho}_{H}+\Lambda)H_{ij}=-\frac{1}{a}4\pi
G{\varrho}_{H}\epsilon_{(i}{}^{\>kl}\sigma_{j)l|k}\ ,\end{split}$ (18b)

where $\Box_{0}{X}_{ij}:=-\ddot{X}_{ij}+\tfrac{1}{a^{2}}\Delta_{0}{X}_{ij}$
denotes the Lagrangian d’Alembertian applied to the tensor field
${X}_{ij}$.888The 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}$, cf. Equation (111) in rza1 , which changes the
coefficient $5\pi G{\varrho}_{H}+\Lambda$ to $4\pi G{\varrho}_{H}+\Lambda$.

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 rza5  in the series takes pressure into account. Contrary to the
Newtonian treatment of a perfect fluid source adlerbuchert , the general-
relativistic approach needs to account for a foliation change, invoking a
lapse function in the simplest case buchert:fluid . This aspect is also
interesting compared with the approach of standard perturbation theory in
cosmology: pressure-supported fluids like 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 rza6  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 kasai95  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 rza6  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 buchertcarfora2 : changes are due to the coupling
of structure formation to the average scalar curvature according to the
general conservation law Eq. (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 parametrization

The family of solutions found by Szekeres Sz75  and later generalized by
Szafron szafron1977inhomogeneous  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 form999Henceforth, we
use units where $c=G=1$. kras1 ; kras2 :

${\hbox{d}}s^{2}=-{\hbox{d}}t^{2}+e^{2\alpha}{\hbox{d}}z^{2}+e^{2\beta}\left({\hbox{d}}x^{2}+{\hbox{d}}y^{2}\right),$
(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
SzaSpatialConfFlat1978  – see CommentSzePress  for a discussion on the
limitations imposed by this property.

In these coordinates, the solution naturally splits into two classes: class I
when $\beta_{,z}\neq 0$ and class II when $\beta_{,z}=0$. In general, neither
class has symmetries SzNoSymm1 , 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
hellaby1996null ; HellabyKrasi2002 ; Bolejko2006Struformation ;
bolejko2007evolution ; WaltersHellaby2012 ; Buckley2013 ; Vrba:2014 ;
Koksbang2015 ; Koksbang2015II ; Koksbang2017 ; Coley:2019ylo ; Hellaby:2007hq
; IshakNwankwo2008 ; Bolejko:2009GERG ; BolCelerier2010 ; BolSuss2011 ;
SussBol2012 ; sussman2015multiple ; sussman2016coarse ; gaspar2018black .

The original solution has been reparametrized multiple times. Among these
parametrizations, the Goode and Wainwright’s GW1  is the one that provides the
most direct connection with RZA. In this representation, the line-element (19)
is rewritten as kras1 ; kras2

${\hbox{d}}s^{2}=-{\hbox{d}}t^{2}+{\cal{S}}^{2}\left({\cal{G}}^{2}{\cal{W}}^{2}{\hbox{d}}z^{2}+e^{2\nu}\left({\hbox{d}}x^{2}+{\hbox{d}}y^{2}\right)\right)\
,$ (20)

where the scale factor ${\cal{S}}$ obeys a Friedmann-like equation:

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

with $k_{0}=0,\pm 1$ and $\mu>0$ being an arbitrary function. Also,

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

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

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

The functions $\mathcal{A}(\bm{r})$, $e^{2\nu(\bm{r})}$, $\mathcal{W}(z)$,
$\beta_{+}(z)$ and $\beta_{-}(z)$ differ for each class, see their specific
form in GW1 ; kras2  or rza6 , 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:

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

  * •

Class II:

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

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

$4\pi\varrho(t,\bm{r})=\frac{3\mu}{{\cal{S}}^{3}}\left(1+\frac{{\cal{F}}}{{\cal{G}}}\right)\
.$ (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)\left(G_{ij}+{h}_{ij}\right)\
,$ (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}\left(\delta^{a}_{\ i}P^{b}_{\
j}+\delta^{b}_{~{}j}P^{a}_{\ i}+P^{a}_{~{}i}P^{b}_{\ j}\right)\ .$ (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\varrho_{b}(t_{i})$), Eq. (21) is precisely the Friedmann equation.
Below, we will relabel ${\cal{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}=\hbox{Diag}\left[\widetilde{{\cal{A}}}^{2},e^{2\nu},e^{2\nu}\right]\quad\hbox{with}\quad\widetilde{{\cal{A}}}:={\cal{A}}-{\cal{F}}_{i}\quad\left(\hbox{and}\quad{\cal{F}}_{i}:={\cal{F}}(t_{i})\right)\
.$ (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{\cal{A}}\widetilde{{\cal{F}}}+\widetilde{{\cal{F}}}^{2}={G}_{ab}\left(\delta^{a}_{~{}3}P^{b}_{~{}3}+\delta^{b}_{~{}3}P^{a}_{~{}3}+P^{a}_{~{}3}P^{b}_{~{}3}\right)=2\widetilde{{\cal{A}}}^{2}P^{3}_{\
3}+\widetilde{{\cal{A}}}^{2}\left(P^{3}_{\ 3}\right)^{2}\ ,$ (28)

where we have defined $\widetilde{{\cal{F}}}={\cal{F}}-{\cal{F}}_{i}$. Thus,
$P^{3}_{\ 3}$ is the only non-trivial component of the deformation field:

$P^{3}_{\ 3}=P=-\widetilde{{\cal{F}}}/\widetilde{{\cal{A}}}\ ,$ (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 Eq. (23) for Szekeres models. Also, all
Szekeres fields are reproduced by the functional expressions of RZA. This
analysis can be summarized in the following remark.

######  Remark.

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
buchert89 .

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 here 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 rza6 .

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 Meures-Bruni_mnras ; Meures-Bruni_prd .

####  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 (${\cal{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}^{a}_{\ i}$ is reinterpreted as a local
deformation field on a Friedmann-like local ‘reference model’ described by Eq.
(21).101010 To avoid confusion, we will use $\hat{P}^{a}_{\ i}$ to refer to
the deformation field with respect a generic reference model and reserve
$P^{a}_{\ i}$ for the deformation with respect to the FLRW background.
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, Eq. (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 BruniPantano1994 ; vanElst:1995 ;
Bolejko:2017SilentUniv1 . Let us consider an arbitrary world line labeled by
the comoving coordinates $\bm{r}_{\ast}$. In the evolution of this fluid line,
the scale factor

$\hat{a}(t):={\cal{S}}(t,z)|_{z=\bm{r}_{\ast}}$ (30)

obeys the Friedmann equations for a model with initial density
$4\pi\varrho_{b}(t_{i})=3\mu(z)|_{z=\bm{r}_{\ast}}$: 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=\bm{r}_{\ast}}$ are interpreted as the growing and decaying
modes on the ‘reference model’.

######  Remark.

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$: ${\cal{S}}(t_{i})={\cal{S}}_{i}(z)$. Then, we define
$A:={\cal{S}}/{\cal{S}}_{i}$, which clearly satisfies $A(t_{i})=1$. Proceeding
along these lines, Eq. (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{\varrho}(z):=3\mu/S_{i}^{3}$:

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

and ${\cal{F}}$ remains as the superposition of the growing and decaying modes
on the ‘reference model’ defined by (31):

$\ddot{{\cal{F}}}+2\frac{\dot{A}}{A}\dot{{\cal{F}}}-\frac{4\pi\hat{\varrho}}{A^{3}}{\cal{F}}=0\
.$ (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
formation111111In this section, we highlight some of the salient results of a
detailed investigation GRZA1 .

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
Bolejko:2017SilentUniv1 .

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,

$\bm{\eta}^{a}=\eta^{a}_{\ i}\mathbf{d}X^{i}=A\left(\delta^{a}_{\
i}+\hat{P}^{a}_{\ i}\right)\mathbf{d}X^{i}\ ,\quad A=A\left(t,\bm{r}\right)\
,\quad\hat{P}^{a}_{\ i}=\hat{P}^{a}_{\ i}\left(t,\bm{r}\right)\ .$ (33)

This ‘reference model’ satisfies Eq. (31) with a general curvature
($\hat{k}=\hat{k}\left(\bm{r}\right)$) and density parameter
($\hat{\varrho}(\bm{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}^{a}_{\ i}$ 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:

$\hbox{RZA}\rightarrow\hbox{GRZA}:\qquad
k_{0}\rightarrow\hat{k}(\bm{r})\quad\land\quad\varrho_{b}(t_{i})\rightarrow\hat{\varrho}(\bm{r})\
.$ (34)

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

GRZA Sze-I RZA Sze-II LTB 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. Also, RZA and the class II of Szekeres (Sze-II) emerge in the
limit when $\hat{k}(\bm{r})\rightarrow k_{0}=\hbox{const.}$ and
$\hat{\varrho}(\bm{r})\rightarrow\varrho_{b}(t_{i})=\hbox{const.}$

####  2.3.1 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,\bm{r})=\frac{\hat{{\cal{F}}}-\hat{{\cal{F}}}_{i}}{{\cal{A}}-\hat{{\cal{F}}}_{i}}\
,\quad\hat{{\cal{F}}}={\cal{F}}+\delta{\cal{F}}\
,\quad\delta{\cal{F}}=\delta\beta_{+}f_{+}+\delta\beta_{-}f_{-}\ .$ (35)

In the equation above, ${\cal{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 Fig. 4.)

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

(a) $\alpha=0$

(b) $\alpha=8.9\times 10^{-4}$

(c) $\alpha=2.2\times 10^{-3}$

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)z\cos\phi$ and $Y=A(t_{0},z)z\sin\phi$. 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 Fig. 3.

##  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 deformation gradient ${\bf d}f^{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 ${\bf d}f^{a}$ by general co-frame fields $\bm{\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 space-time, 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 space-time,
including tilted foliations (for a comprehensive investigation and the
definition of the Lagrangian picture in the general situation, see
buchert:generalfluid  and foliations ).

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
parametrization of the Szekeres solutions, since this allowed 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 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 ${\cal{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
buchertcarfora2 , and the exact volume-partitioning formulas of a multiscale
dynamics multiscale . 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
generalbackground .

TB: conceptualization, writing, methodology, supervision, funding acquisition;
IDG: writing, methodology, software, visualization; JJO: writing, funding
acquisition, software supervision. All authors have read and agreed to the
published version of the manuscript.

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:
TB). IDG and JJO acknowledge support from the National Science Centre (NCN,
Poland) under the Sonata-15 research grant UMO-2019/35/D/ST9/00342. IDG
acknowledges hospitality and support by ERC ARTHUS during a working visit in
Lyon.

The inhomog code referred to in the context of realizing the backreaction
functional in RZA is openly available inhomog . 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.

###### Acknowledgements.

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. The authors declare no conflict of
interest. References

## References

  * (1) Adler, S. and Buchert, T. (1999). Lagrangian theory of structure formation in pressure-supported cosmological fluids. Astron. Astrophys. 343, 317. [arXiv:astro-ph/9806320]
  * (2) Afshordi, N., Mohayaee, R. and Bertschinger, E. (2009). Hierarchical phase space structure of dark matter haloes: Tidal debris, caustics, and dark matter annihilation. Phys. Rev. D 79, 083526. [arXiv:0811.1582]
  * (3) Al Roumi, F., Buchert, T. and Wiegand, A. (2017). Lagrangian theory of structure formation in relativistic cosmology. IV. Lagrangian approach to gravitational waves. Phys. Rev. D 96, 123538. [arXiv:1711.01597]
  * (4) Alles, A., Buchert, T., Al Roumi, F. and Wiegand, A. (2015). Lagrangian theory of structure formation in relativistic cosmology. III. Gravitoelectric perturbation and solution schemes at any order. Phys. Rev. D 92, 023512. [arXiv:1503.02566]
  * (5) Angulo, R.E. and Hahn, O. (2022). Large-scale dark matter simulations. Living Reviews in Comp. Astrophys. 8, 1. [arXiv:2112.05165]
  * (6) Arnol’d, V.I., Gusein–Sade, S.M. and Varchenko, A.N. (1985). Singularities of Differentiable Maps Vol. I, Birkhäuser, Boston; Orig. Nauka, Moscow (1982).
  * (7) Arnol’d, V.I., Shandarin, S.F. and Zel’dovich Ya. B. (1982). The large-scale structure of the Universe. I. General properties. One- and two-dimensional models. Geophys. Astrophys. Fluid Dyn. 20, 111.
  * (8) Bardeen, J.M. (1980). Gauge-invariant cosmological perturbations. Phys. Rev. D 22, 1882.
  * (9) Bartelmann, M. (2014). Trajectories of point particles in cosmology and the Zel’dovich approximation. Phys. Rev. D 91, 083524. [arXiv:arXiv:1411.0805]
  * (10) Bartelmann, M. and Schneider, P. (1992). A large-scale structure model for gravitational lensing. Astron. Astrophys. 259, 413.
  * (11) Bernardeau, F., Colombi, S., Gaztanaga, E. and Scoccimarro, R. (2002). Large-scale structure of the Universe and cosmological perturbation theory. Phys. Rep. 367, 1. [arXiv:astro-ph/0112551]
  * (12) Bertschinger, E. (1995). Cosmological Dynamics. (Les Houches Summer School Lectures on Cosmology, August 1993); Massachusetts Institute of Technology: Cambridge, MA, USA. [arXiv:astro-ph/9503125]
  * (13) Bildhauer, S., Buchert, T. and Kasai, M. (1992). Solutions in Newtonian cosmology – The pancake theory with cosmological constant. Astron. Astrophys. 263, 23.
  * (14) Bolejko, K. (2006). Structure formation in quasispherical Szekeres model. Phys. Rev. D 73, 123508. [arXiv:astro-ph/0604490]
  * (15) Bolejko, K. (2007). Evolution of cosmic structures in different environments in the quasispherical Szekeres model. Phys. Rev. D 75, 043508. [arXiv:astro-ph/0610292]
  * (16) Bolejko, K. (2009). The Szekeres Swiss Cheese model and the CMB observations. Gen. Rel. Grav. 41, 1737. [arXiv:0804.1846]
  * (17) Bolejko, K. (2017). Relativistic numerical cosmology with Silent Universes. Class. Quant. Grav. 35, 024003. [arXiv:1708.09143]
  * (18) Bolejko, K. and Célérier, M.N. (2010). Szekeres Swiss-Cheese model and supernova observations. Phys. Rev. D 82, 103510. [arXiv:1005.2584]
  * (19) Bolejko K. and Sussman, R.A. (2011). Cosmic spherical void via coarse-graining and averaging non-spherical structures. Phys. Lett. B 697, 265. [arXiv:1008.3420]
  * (20) Bonnor, W.B., Sulaiman, A.H. and Tomimura, N. (1977). Szekeres’s Space-Times Have No Killing Vectors. Gen. Rel. Grav. 8, 549.
  * (21) Bouchet, F.R., Juszkiewicz, R., Colombi, S. and Pellat, R. (1992). Weakly nonlinear gravitational instability for arbitrary Omega. The Astrophys. J. 394, L5.
  * (22) Bouchet, F.R., Colombi, S., Hivon, E. and Juszkiewicz, R. (1995). Perturbative Lagrangian approach to gravitational instability. Astron. Astrophys. 296, 575.
  * (23) Brenier, Y., Frisch, U., Henon, M., Loeper, G., Matarrese, S., Mohayaee, R. and Sobolevskii, A. (2003). Reconstruction of the early Universe as a convex optimization problem. M.N.R.A.S. 346, 501. [arXiv:astro-ph/0304214]
  * (24) Bruni, M., Matarrese, S. and Pantano, O. (1995). Dynamics of silent universes. Astrophys. J.445, 958. [arXiv:astro-ph/9406068]
  * (25) Brustein, R., Gasperini, M., Giovannini, M., Mukhanov, V. F. and Veneziano, G. (1995), Metric perturbations in dilaton-driven inflation. Phys. Rev. D 51, 6744. [arXiv:hep-th/9501066]
  * (26) Buchert, T. (1989). A class of solutions in Newtonian cosmology and the pancake theory. Astron. Astrophys. 223, 9.
  * (27) Buchert, T. (1992). Lagrangian theory of gravitational instability of Friedmann–Lemaître cosmologies and the ’Zel’dovich approximation’. M.N.R.A.S. 254, 729.
  * (28) Buchert, T. (1993). Lagrangian perturbation theory – A key-model for large-scale structure. Astron. Astrophys. 267, L51.
  * (29) Buchert, T. (1994). Lagrangian theory of gravitational instability of Friedmann–Lemaître cosmologies – a generic third–order model for nonlinear clustering. M.N.R.A.S. 267, 811. [arXiv:astro-ph/9309055]
  * (30) Buchert, T. (1996). Lagrangian perturbation approach to the formation of large-scale structure. in Proc. International School Enrico Fermi Course CXXXII Dark Matter in the Universe, Varenna 1995, ed. S. Bonometto, J.R. Primack and A. Provenzale, IOS Press Amsterdam, pp.543-64. [arXiv:astro-ph/9509005]
  * (31) Buchert, T. (2000). On average properties of inhomogeneous fluids in general relativity I: Dust cosmologies. Gen. Rel. Grav. 32 105. [arXiv:gr-qc/9906015]
  * (32) Buchert, T. (2001). On average properties of inhomogeneous fluids in general relativity II: Perfect fluid cosmologies. Gen. Rel. Grav. 33 1381. [arXiv:gr-qc/0102049]
  * (33) Buchert, T. (2006). The non-perturbative regime of cosmic structure formation. Astron. Astrophys. 454, 415.
  * (34) Buchert, T. (2008). Dark Energy from structure: a status report. Gen. Rel. Grav. 40 467. [arXiv:0707.2153]
  * (35) Buchert, T. (2011). Toward physical cosmology: focus on inhomogeneous geometry and its non-perturbative effects. Class. Quant. Grav. 28 164007. [arXiv:1103.2016]
  * (36) Buchert, T. and Bartelmann, M. (1991). High-spatial resolution in three dimensions: a challenge for large-scale structure formation models. Astron. Astrophys. 251, 389.
  * (37) Buchert, T. and Carfora, M. (2002). Regional averaging and scaling in relativistic cosmology. Class. Quant. Grav. 19, 6109 [arXiv:gr-qc/0210037]
  * (38) Buchert, T. and Carfora, M. (2008). On the curvature of the present-day Universe. Class. Quant. Grav. 25, 195001 [arXiv:0803.1401]
  * (39) Buchert, T. and Domínguez, A. (1998). Modeling multi-stream flow in collisionless matter: approximations for large-scale structure beyond shell-crossing. Astron. Astrophys. 335, 395. [arXiv:astro-ph/9702139]
  * (40) Buchert, T. and Domínguez, A. (2005). Adhesive Gravitational Clustering. Astron. Astrophys. 438, 443. [arXiv:astro-ph/0502318]
  * (41) Buchert, T., Domínguez, A. and Pérez-Mercader, J. (1999). Extending the scope of models for large-scale structure formation in the Universe. Astron. Astrophys. 349, 343. [arXiv:astro-ph/9709218]
  * (42) Buchert, T. and Ehlers, J. (1993). Lagrangian theory of gravitational instability of Friedmann–Lemaître cosmologies – second-order approach: an improved model for nonlinear clustering. M.N.R.A.S. 264, 375.
  * (43) Buchert, T. and Ehlers, J. (1997). Averaging inhomogeneous Newtonian cosmologies. Astron. Astrophys. 320, 1. [arXiv:astro-ph/9510056]
  * (44) Buchert, T., van Elst, H. and Heinesen, A. (2022). The averaging problem on the past null cone in inhomogeneous dust cosmologies. Gen. Rel. Grav., submitted. [arXiv:2202.10798]
  * (45) Buchert, T. and Götz, G. (1987). A class of solutions for self-gravitating dust in Newtonian gravity. J. Math. Phys. 28, 2714.
  * (46) Buchert, T., Karakatsanis, G., Klaffl, R. and Schiller, P. (1997). The performance of Lagrangian perturbation schemes at high resolution. Astron. Astrophys. 318, 1.
  * (47) Buchert, T., Kerscher, M. and Sicka C. (2000). Backreaction of inhomogeneities on the expansion: the evolution of cosmological parameters. Phys. Rev. D 62, 043525. [arXiv:astro-ph/9912347]
  * (48) Buchert, T., Larena, J. and Alimi, J.M. (2006). Correspondence between kinematical backreaction and scalar field cosmologies – the ‘morphon field’ Class. Quant. Grav. 23 6379. [arXiv:gr-qc/0606020]
  * (49) Buchert, T. and Mädler, T. (2019). Editorial Note to: On the Newtonian Limit of Einstein’s Theory of Gravitation (by Jürgen Ehlers). Gen. Rel. Grav. 51, 162. [arXiv:1910.12106]
  * (50) Buchert, T., Melott, A.L. and Weiß, A.G. (1994). Testing higher-order Lagrangian perturbation theory against numerical simulations – 1. pancake models. Astron. Astrophys. 288, 349.
  * (51) Buchert, T., Mourier, P. and Roy, X. (2018). On cosmological backreaction and its dependence on spacetime foliation. Class. Quant. Grav. 35 24LT02. [arXiv:1805.10455]
  * (52) Buchert, T., Mourier, P. and Roy, X. (2020). On average properties of inhomogeneous fluids in general relativity III: general fluid cosmologies. Gen. Rel. Grav. 52 27. [arXiv:1912.04213]
  * (53) Buchert, T., Nayet, C. and Wiegand, A. (2013). Lagrangian theory of structure formation in relativistic cosmology. II. Average properties of a generic evolution model. Phys. Rev. D 87 123503. [arXiv:1303.6193]
  * (54) Buchert, T. and Obadia N. (2011). Effective inhomogeneous inflation: curvature inhomogeneities of the Einstein vacuum. Class. Quant. Grav. 28, 162002. [arXiv:1010.4512]
  * (55) Buchert, T. and Ostermann, M. (2012). Lagrangian theory of structure formation in relativistic cosmology. I. Lagrangian framework and definition of a nonperturbative approximation. Phys. Rev. D 86 023520. [arXiv:1203.6263]
  * (56) Buchert, T. and Räsänen, S. (2012). Backreaction in Late-Time Cosmology. Annu. Rev. Nucl. Part. Sci. 62 57. [arXiv:1112.5335]
  * (57) Buchert, T., Shandarin, S.F. and Weiß, A.G. (1995). Building blocks of large-scale structure in the Universe. Three-dimensional singularities and metamorphoses. unpublished.
  * (58) Buckley, R.G. and Schlegel, E.M. (2013). CMB dipoles and other low-order multipoles in the quasispherical Szekeres model. Phys. Rev. D 87, 023524. [arXiv:1907.08684]
  * (59) Carfora, M. and Piotrkowska, K. (1995). A Renormalization Group Approach to Relativistic Cosmology. Phys. Rev. D 52, 4393. [arXiv:gr-qc/9502021]
  * (60) Coles, P., Melott, A.L. and Shandarin, S.F. (1993). Testing approximations for nonlinear gravitational clustering. M.N.R.A.S. 260, 765.
  * (61) Coley, A.A., Layden, N. and McNutt, D. (2019) An invariant characterization of the quasi-spherical Szekeres dust models. Gen. Rel. Grav.51, 164 [arXiv:arXiv:2103.05305]
  * (62) Combi, L. (2021). Spacetime is material. In Contemporary Materialism. Its Ontology and Epistemology, Eds.: Romero, E.G., Pérez Jara, J. and Camprubí, L., Synthese Library, Springer Berlin/Heidelberg, Germany. [arXiv:2108.01712]
  * (63) Crocce, M. and Scoccimarro, R. (2006). Renormalized cosmological perturbation theory. Phys. Rev. D 73, 063519. [arXiv:astro-ph/0509418]
  * (64) Dawson, J.M. (1959). Nonlinear Electron Oscillations in a Cold Plasma. Phys. Rev. 113, 383.
  * (65) Delgado Gaspar, I. and Buchert, T. (2021). Lagrangian theory of structure formation in relativistic cosmology. VI. Comparison with Szekeres exact solutions. Phys. Rev. D 103, 023513. [arXiv:2009.06339]
  * (66) Delgado Gaspar, I., Sussman, R.A., McNutt, D.D. and Coley, A.A. (2021). Comment on ‘Szekeres universes with homogeneous scalar fields’. Eur. Phys. J. C 81, 310. [arXiv:2103.06862]
  * (67) Delgado Gaspar, I., Hidalgo, J.C., Sussman, R.A. and Quiros, I. (2018). Black hole formation from the gravitational collapse of a nonspherical network of structures. Phys. Rev. D 97, 104029. [arXiv:1802.09123]
  * (68) Delgado Gaspar, I., Buchert, T. and Ostrowski, J.J. (2022). Beyond relativistic Lagrangian perturbation theory. I. An exact-solution controlled model for structure formation. Phys. Rev. D, submitted. [arXiv:2210.04004]
  * (69) Desgrange, C., Heinesen, A. and Buchert, T. (2019). Dynamical spatial curvature as a fit to type Ia supernovae _Int. J. of Mod. Phys. D_ 28, 1950143 [arXiv:1902.07915]
  * (70) Desjacques, V., Jeong, D. and Schmidt, F. (2018). Large-scale galaxy bias. Phys. Rep. 733, 1–193. [arXiv:1611.09787]
  * (71) Doroshkevich A.G. (1970). The space structure of perturbations and the origin of rotation of galaxies in the theory of fluctuation. Astrofizika 6, 581 (1973); Translated in Astrophysics 6, 320.
  * (72) Doroshkevich, A.G., Ryabenkii, V.S. and Shandarin, S.F. (1973). Non-linear theory of development of potential perturbations. Astrofizika 9, 257 (1975); Translated in Astrophysics 9, 144.
  * (73) Doroshkevich, A.G., Kotok, E.V., Novikov, I.D., Polyudov, A.N., Shandarin, S.F. and Sigov, Yu.S. (1980). Two-dimensional simulation of the gravitational system dynamics and formation of the large-scale structure of the universe. M.N.R.A.S. 192, 321.
  * (74) Doumler, T., Hoffman, Y., Courtois, H. and Gottlöber, S. (2013). Reconstructing cosmological initial conditions from galaxy peculiar velocities. I. Reverse Zel’dovich Approximation. M.N.R.A.S. 430, 888. [arXiv:1212.2806]
  * (75) Durrer, R. (1994). Gauge Invariant Cosmological Perturbation Theory: A General Study and It’s Application to the Texture Scenario of Structure Formation. Fundamentals of Cosmic Physics 15, 209. [arXiv:astro-ph/9311041]
  * (76) Ehlers, J. (1961). Beiträge zur relativistischen Mechanik kontinuierlicher Medien. Akad. Wiss. Lit. (Mainz); Abh. Math.-Nat. Kl. No. 11, 793 (1961); translation: Contributions to the relativistic mechanics of continuous media. Gen. Rel. Grav. 25, 1225 (1993).
  * (77) Ehlers, J. and Buchert, T. (1997). Newtonian cosmology in Lagrangian formulation: foundations and perturbation theory. Gen. Rel. Grav. 29, 733. [arXiv:astro-ph/9609036]
  * (78) Ellis, G.F.R and Bruni, M. (1989). Covariant and gauge-invariant approach to cosmological density fluctuations. Phys. Rev. D 40, 1804.
  * (79) Ellis, G.F.R. and Buchert, T. (2005). The universe seen at different scales. Phys. Lett. A 347, 38. [arXiv:gr-qc/0506106]
  * (80) Fillmore, J.A. and Goldreich, P. (1984). Self-similar gravitational collapse in an expanding Universe. The Astrophys. J. 281, 1.
  * (81) Gasperini, M. (1997). Tensor perturbations in high-curvature string backgrounds. Phys. Rev. D 56, 4815. [arXiv:gr-qc/9704045]
  * (82) Gasperini, M., Marozzi, G. and Veneziano, G. (2010). A covariant and gauge invariant formulation of the cosmological ‘backreaction’ JCAP 02, 009. [arXiv:0912.3244]
  * (83) Gasperini, M., Marozzi, G., Nugier, F. and Veneziano, G. (2011). Light-cone averaging in cosmology: formalism and applications. JCAP 07, 008. [arXiv:1104.1167]
  * (84) Goode, S.W. and Wainwright, J. (1982). Singularities and evolution of the Szekeres cosmological models. Phys. Rev. D 26, 3315.
  * (85) Goodman, J. and Binney, J (1983). Gravitational collapse of pressureless inhomogeneous spheroids. M.N.R.A.S. 203, 265.
  * (86) Grinstein, B. and Wise, M.B. (1987). On the validity of the Zel’dovich Approximation. The Astrophys. J. 320, 448; Erratum (1989). The Astrophys. J. 337, 579.
  * (87) Gurbatov, S.N., Saichev, A.I. and Shandarin, S.F. (1989). The large-scale structure of the universe in the frame of the model equation of non-linear diffusion. M.N.R.A.S. 236, 385.
  * (88) Hahn, O., Abel, T. and Kaehler, R. (2012). Tracing the Dark Matter sheet in phase space. M.N.R.A.S. 427, 61. [arXiv:1111.3944]
  * (89) Hahn, O., Abel, T. and Kaehler, R. (2013). A new approach to simulating collisionless dark matter fluids. M.N.R.A.S. 434, 1171. [arXiv:1210.6652]
  * (90) Hawking, S.W. (1966). Perturbations of an Expanding Universe. The Astrophys. J. 145, 544.
  * (91) Heinesen, A. and Buchert, T. (2020). Solving the curvature and Hubble parameter inconsistencies through structure formation-induced curvature. _Class. Quant. Grav._ 37, 164001. (Focus issue on the Hubble constant tension); Erratum (2020). Class. Quant. Grav. 37, 229601. [arXiv:2002.10831]
  * (92) Heinesen, A., Mourier, P. and Buchert, T. (2019). On the covariance of scalar averaging and backreaction in relativistic inhomogeneous cosmology. _Class. Quant. Grav._ 36 075001. [arXiv:1811.01374]
  * (93) Hellaby, C. (1996). The null and KS limits of the Szekeres model. Class. Quant. Grav. 13, 2537.
  * (94) Hellaby C. and Krasiński, A. (2002). You can’t get through Szekeres wormholes: Or, regularity, topology and causality in quasispherical Szekeres models. Phys. Rev. D 66, 084011. [arXiv:gr-qc/0206052]
  * (95) Hellaby C. and Krasiński, A. (2008). Physical and geometrical interpretation of the $\epsilon\leq 0$ Szekeres models. Phys. Rev. D 77, 023529 (2008) [arXiv:0710.2171]
  * (96) Hidding, J., Shandarin, S.F. and van de Weygaert, R. (2014). The Zel’dovich approximation: key to understanding cosmic web complexity. M.N.R.A.S. 437, 3442.
  * (97) Hivon, E., Bouchet, F.R., Colombi, S. and Juszkiewicz, R. (1995). Redshift distortions of clustering: a Lagrangian approach. Astron. Astrophys. 298, 643. [arXiv:astro-ph/9407049]
  * (98) Ishak, M., Richardson, J., Garred, D., Whittington, D., Nwankwo, A., and Sussman, R.A. (2011). Dark energy or apparent acceleration due to a relativistic cosmological model more complex than the Friedmann-Lemaître-Robertson Walker model? Phys. Rev. D 84, 089902. [arXiv:0708.2943]
  * (99) Kasai, M. (1993). Inhomogeneous cosmological models which are homogeneous and isotropic on average. Phys. Rev. D 47, 3214.
  * (100) Kasai, M. (1995). Tetrad-based perturbative approach to inhomogeneous universes: a general relativistic version of the Zel’dovich approximation. Phys. Rev. D 52, 5605.
  * (101) Kitaura, F.-S. and Hess, S. (2013). Cosmological structure formation with augmented Lagrangian perturbation theory. M.N.R.A.S. 435, L78. [arXiv:1212.3514]
  * (102) Kitaura, F.-S., Balaguera-Antolínez, A., Sinigaglia, F. and Pellejero-Ibàñez, M. (2022). The cosmic web connection to the dark matter halo distribution through gravity. M.N.R.A.S. 512, 2245. [arXiv:2005.11598]
  * (103) Klypin, A.A. and Shandarin, S.F. (1983). Three-dimensional numerical model of the formation of large-scale structure in the Universe. M.N.R.A.S. 204, 891.
  * (104) Kodama, H. and Sasaki, M. (1984). Cosmological Perturbation Theory. Progress of Theoretical Physics Supp. 78, 1.
  * (105) Kokron, N., Chen, S.F., White, M., DeRose, J. and Maus, M. (2022). Accurate predictions from small boxes: variance suppression via the Zel’dovich approximation. JCAP 09, 059. [arXiv:2205.15327]
  * (106) Koksbang, S.M. and Hannestad, S. (2015). Methods for studying the accuracy of light propagation in N-body simulations. Phys. Rev. D 91, 043508. [arXiv:1501.01413]
  * (107) Koksbang, S.M. and Hannestad, S. (2015). Studying the precision of ray tracing techniques with Szekeres models. Phys. Rev. D 92, 023532. [arXiv:1506.09127]
  * (108) Koksbang, S.M. (2017). Light propagation in Swiss-cheese models of random close-packed Szekeres structures: Effects of anisotropy and comparisons with perturbative results. Phys. Rev. D 95, 063532. [arXiv:1703.03572]
  * (109) Krasiński A. (1997). Inhomogeneous Cosmological Models. Cambridge University Press: Cambridge, UK. ISBN 0521481805
  * (110) Lemaître, G. (1933). The expanding universe. Annales Soc. Sci. Bruxelles A 53, 51.
  * (111) Li, Y.Z., Mourier, P., Buchert, T. and Wiltshire, D.L. (2018). Lagrangian theory of structure formation in relativistic cosmology. V. Irrotational fluids. Phys. Rev. D 98, 043507. [arXiv:1806.05467]
  * (112) Lifshitz, E.M. (1946). On the gravitational stability of the expanding universe. ZhETF 16, 587; J. Phys. USSR Acad. Sci. 10, 116 (1946); Republication: Gen. Rel. Grav. 49, 18 (2017).
  * (113) Lifshitz, E.M. and Khalatnikov, I.M. (1963). Investigations in relativistic cosmology. Usp. Fiz. Nauk 80, 391; Sov. Phys. Usp. 6, 522 (1964); Adv. Phys. 12, 185 (1963).
  * (114) Lin, C.C., Mestel, L. and Shu, F.H. (1965). The Gravitational Collapse of a Uniform Spheroid. The Astrophys. J. 142, 1431.
  * (115) Malik, K.A. and Wands, T. (2009). Cosmological perturbations. Phys. Rep. 475, 1. [arXiv:0809.4944]
  * (116) Matarrese, S, Pantano, O. and Saez, D. (1994). General relativistic dynamics of irrotational dust: Cosmological implications. Phys. Rev. Lett. 72, 320. [arXiv:astro-ph/9310036]
  * (117) Matarrese, S, Pantano, O. and Saez, D. (1994). A Relativistic Approach to Gravitational Instability in the Expanding Universe – Second-Order Lagrangian Solutions. M.N.R.A.S. 271, 513. [arXiv:astro-ph/9403032]
  * (118) Matarrese, S. and Terranova, D. (1996). Post-Newtonian cosmological dynamics in Lagrangian coordinates. M.N.R.A.S. 283, 400. [arXiv:astro-ph/9511093]
  * (119) Matsubara, T. (2008). Resumming cosmological perturbations via the Lagrangian picture: One-loop results in real space and in redshift space. Phys. Rev. D 77, 063530. [arXiv:0711.2521]
  * (120) Matsubara, T. (2015). Recursive solutions of Lagrangian perturbation theory. Phys. Rev. D 92, 023534. [arXiv:1505.01481]
  * (121) Melott, A.L., Buchert, T. and Weiß, A.G. (1995). Testing higher-order Lagrangian perturbation theory against numerical simulations – 2. hierarchical models. Astron. Astrophys. 294, 345.
  * (122) Melott, A.L and Shandarin, S.F. (1989). Gravitational Instability with High Resolution. The Astrophys. J. 343, 26.
  * (123) Melott, A.L., Pellman, T.F. and Shandarin, S.F. (1994). Optimizing the Zel’dovich approximation. M.N.R.A.S. 269, 626. [arXiv:astro-ph/9312044]
  * (124) Meures, N. and and Bruni, M. (2011). Exact nonlinear inhomogeneities in $\Lambda$CDM cosmology. Phys. Rev. D 83, 123519. [arXiv:arXiv:1103.0501]
  * (125) Meures, N. and and Bruni, M. (2012). Redshift and distances in a $\Lambda$CDM cosmology with non-linear inhomogeneities. M.N.R.A.S. 419, 3 [arXiv:arXiv:1107.4433]
  * (126) Michaux, M., Hahn, O., Rampf, C. and Angulo, R.E. (2021). Accurate initial conditions for cosmological N-body simulations: minimizing truncation and discreteness errors. M.N.R.A.S. 500, 663. [arXiv:2008.09588]
  * (127) Monaco, P. (2016). Approximate methods for the generation of dark matter halo catalogs in the age of precision cosmology. Galaxies 4, 53. [arXiv:1605.07752]
  * (128) Monaco, P., Sefusatti, E., Borgani, S., Crocce, M., Fosalba, P., Sheth, R. K. and Theuns, T. (2013). An accurate tool for the fast generation of dark matter halo catalogues. M.N.R.A.S. 433, 2389. [arXiv:1305.1505]
  * (129) Munari, E., Monaco, P., Sefusatti, E., Castorina, E., Mohammad, F.G., Anselmi, S. and Borgani, S. (2017). Improving fast generation of halo catalogs with higher-order Lagrangian perturbation theory. M.N.R.A.S. 465, 4658. [arXiv:1605.04788]
  * (130) Morita, M., Nakamura, K. and Kasai, M. (1998). Relativistic Zel’dovich approximation in a spherically symmetric model. Phys. Rev. D 57, 6094. [arXiv:astro-ph/9711026]
  * (131) Moutarde, F., Alimi, J.-M., Bouchet, F.R., Pellat, R. and Ramani, A. (1991). Precollapse scale invariance in gravitational instability. The Astrophys. J. 382, 377.
  * (132) Mukhanov, V.F., Feldman, H.A. and Brandenberger, R.H. (1992). Theory of cosmological perturbations. Phys. Rep. 215, 203.
  * (133) Nadkarni-Ghosh, S. and Chernoff, D.F. (2011). Extending the domain of validity of the Lagrangian approximation. M.N.R.A.S. 410, 1454. [arXiv:1005.1217]
  * (134) Nadkarni-Ghosh, S. and Chernoff, D.F. (2013). Modelling non-linear evolution using Lagrangian Perturbation Theory (LPT) re-expansions. M.N.R.A.S. 431, 799. [arXiv:1211.5777]
  * (135) Nadkarni-Ghosh, S. and Singhal, A. (2016). Phase Space dynamics of triaxial collapse: Joint density-velocity evolution. M.N.R.A.S. 457, 2773. [arXiv:1407.1945]
  * (136) Parry, J., Salopek, D.S. and Stewart, J.M. (1994). Solving the Hamilton-Jacobi equation for general relativity. Phys. Rev. D 49, 2872.
  * (137) Peebles, P.J.E. (1980). The Large-Scale Structure of the Universe. Princeton University Press: Princeton, JS, USA. ISBN 9780691209838
  * (138) Percival, W.J. and White, M. (2009). Testing cosmological structure formation using redshift-space distortions. M.N.R.A.S. 393, 29712. [arXiv:0808.0003]
  * (139) Plebański J. and Krasiński A. (2006). An introduction to general relativity and cosmology. Cambridge University Press: Cambridge, UK. ISBN 0-521-85623-X
  * (140) Rampf, C. (2012). The recursion relation in Lagrangian perturbation theory. JCAP 12, 004. [arXiv:1205.5274]
  * (141) Rampf, C. (2021). Cosmological Vlasov-Poisson equations for dark matter. Rev. Mod. Plasma Phys. 5, 10. [arXiv:2110.06265]
  * (142) Rampf, C. and Buchert, T. (2012). Lagrangian perturbations and the matter bispectrum I: fourth-order model for non-linear clustering. JCAP 06, 021. [arXiv:1203.4260]
  * (143) Rampf, C. and Frisch, U. (2017). Shell-crossing in quasi-one-dimensional flow. M.N.R.A.S. 471, 671. [arXiv:1705.08456]
  * (144) Rampf, C. and Hahn, O. (2021). Shell-crossing in a $\Lambda$CDM Universe. M.N.R.A.S. 501, L71. [arXiv:2010.12584]
  * (145) Rampf, C. and Rigopoulos, G. (2013). Initial conditions for cold dark matter particles and general relativity Phys. Rev. D 87, 123525. [arXiv:1305.0010]
  * (146) Rampf, C., Frisch, U. and Hahn, O. (2021). Unveiling the singular dynamics in the cosmic large-scale structure. M.N.R.A.S. 505, L90.
  * (147) Ray, S.S., Frisch, U., Nazarenko, S. and Matsumoto, T. (2011). Resonance phenomenon for the Galerkin-truncated Burgers and Euler equations. Phys. Rev. E 84, 016301. [arXiv:1011.1826]
  * (148) Rampf, C., Villone, B. and Frisch, U. (2015). How smooth are particle trajectories in a $\Lambda$CDM Universe? M.N.R.A.S. 452, 1421. [arXiv:1504.00032]
  * (149) Rigopoulos, G. and Valkenburg, W. (2012). Nonlinear approximation for perturbations in $\Lambda$CDM. Phys. Rev. D 86, 043523. [arXiv:1203.2796]
  * (150) Roukema, B.F., Ostrowski, J.J. and Borkowska, J. (2017). inhomog: Biscale kinematical backreaction analytical evolution. Astrophysics Source Code Library ascl:1711.002.
  * (151) Roy, X. and Buchert, T. (2012). Relativistic cosmological perturbation scheme on a general background: scalar perturbations for irrotational dust. Class. Quant. Grav. 29, 115004. [arXiv:1202.5766]
  * (152) Russ, H., Morita, M., Kasai, M. and Börner, G. (1996). Zel’dovich-type approximation for an inhomogeneous universe in general relativity: Second-order solutions. Phys. Rev. D 53, 6881. [arXiv:astro-ph/9512071]
  * (153) Russ, H., Soffel, M.H., Kasai, M. and Börner, G. (1997). Age of the universe: Influence of the inhomogeneities on the global expansion factor. Phys. Rev. D 56, 2044. [arXiv:astro-ph/9612218]
  * (154) Saga, S., Taruya, A. and Colombi, S. (2018). Lagrangian Cosmological Perturbation Theory at Shell Crossing. Phys. Rev. Lett. 121, 241302. [arXiv:1805.08787]
  * (155) Saga, S., Taruya, A. and Colombi, S. (2022). Cold dark matter protohalo structure around collapse: Lagrangian cosmological perturbation theory versus Vlasov simulations. Astron. Astrophys. 664, A3. [arXiv:2111.08836]
  * (156) Sahni, V. and Coles, P. (1995). Approximation methods for nonlinear gravitational clustering. Phys. Rep. 262, 1. [arXiv:astro-ph/9505005]
  * (157) Salopek, D.S. and Stewart, J.M. (1992). Hamilton-Jacobi theory for general relativity with matter fields. Class. Quant. Grav. 9, 1943.
  * (158) Schmidt, F. (2021). An n-th order Lagrangian forward model for large-scale structure. JCAP 04, 033. [arXiv:2012.09837]
  * (159) Serrin, I. (1959). Mathematical Principles of Classical Fluid Mechanics. In: Encyclopedia of Physics VIII.1, 125–263, Springer: Berlin/Heidelberg, Germany.
  * (160) Shandarin, S.F. (1980). Evolution of Perturbations in Friedmann Models of the Universe. Astrofizika 16, 769; Translated in Astrophysics 16, 439 (1981).
  * (161) Shandarin, S.F. and Zel’dovich, Ya.B. (1989). The large-scale structure of the universe: Turbulence, intermittency, structures in a self-gravitating medium. Rev. Mod. Phys. 61, 185.
  * (162) Shukurov, A.M. (1981). Nonlinear growth of density perturbations in an expanding collisionless medium. Astrofizika 17, 469; Astrophysics 17, 263 (1982).
  * (163) Spero, A. and Szafron, D. A. (1978). Spatial conformal flatness in homogeneous and inhomogeneous cosmologies. J. Math. Phys. 19, 1536.
  * (164) Stewart, J.M., (1990). Perturbations of Friedmann-Robertson-Walker cosmological models. Class. Quant. Grav. 7, 1169.
  * (165) Sussman, R.A. and Bolejko K., (2012). A Novel approach to the dynamics of Szekeres dust models. Class. Quant. Grav. 29, 065018. [arXiv:1109.1178]
  * (166) Sussman, R.A. and Delgado Gaspar, I. (2015). Multiple nonspherical structures from the extrema of Szekeres scalars. Phys. Rev. D 92, 083533. [arXiv:1508.03127]
  * (167) Sussman, R.A., Delgado Gaspar, I. and Hidalgo, J.C. (2016). Coarse-grained description of cosmic structure from Szekeres models. JCAP 03, 012. [arXiv:1507.02306]
  * (168) Susperregi, M. and Buchert, T. (1997). Cosmic density and velocity fields in Lagrangian perturbation theory. Astron. Astrophys. 323, 295.
  * (169) Szafron, D.A. (1977). Inhomogeneous cosmologies: New exact solutions and their evolution. J. Math. Phys. 18, 1673.
  * (170) Szekeres, P. (1975). Quasispherical Gravitational Collapse. Phys. Rev. D 12, 2941.
  * (171) Tolman, R.C. (1934). Effect of inhomogeneity on cosmological models. Proceedings of the National Academy of Sciences 20, 169.
  * (172) Tomita, K. (1967). Non-Linear Theory of Gravitational Instability in the Expanding Universe. Prog. Theor. Phys. 730, 831.
  * (173) van Elst, H. and Ellis, G.F.R. (1996) The Covariant approach to LRS perfect fluid space-time geometries. Class. Quant. Grav. 13, 1099. [arXiv:gr-qc/9510044]
  * (174) Vigneron, Q. (2021). The 1+3-Newton-Cartan system and Newton-Cartan cosmology. Phys. Rev. D 103, 064064. [arXiv:2012.10213]
  * (175) Vigneron, Q. and Buchert, T. (2019). Dark Matter from Backreaction? Collapse models on galaxy cluster scales. Class. Quant. Grav. 36, 175006. [arXiv:1902.08441]
  * (176) Vitenti, S.D.P., Falciano, F.T. and Pinto-Neto, N. (2014). Covariant Bardeen perturbation formalism. Phys. Rev. D 89,103538. [arXiv:1311.6730]
  * (177) Vrba, D. and Svítek, O. (2014). Modelling inhomogeneity in Szekeres spacetime. Gen. Rel. Grav. 46, 1808. [arXiv:1412.7203]
  * (178) Walters A. and Hellaby, C. (2012). Constructing realistic Szekeres models from initial and final data. Phys. Rev. D 1212, 001. [arXiv:1211.2110]
  * (179) Weinberg, D.H. and Gunn, J.E. (1990). Large-scale Structure and the Adhesion Approximation. M.N.R.A.S. 247, 260.
  * (180) Weiß, A.G., Gottlöber, S. and Buchert, T. (1996). Optimizing higher-order Lagrangian perturbation theory for standard CDM and BSI models. M.N.R.A.S. 278 953.
  * (181) Wiegand, A. and Buchert, T. (2010). Multiscale cosmology and structure-emerging Dark Energy: A plausibility analysis. Phys. Rev. D 82, 023523. [arXiv:1002.3912]
  * (182) Yano, T., Koyama, H., Buchert, T. and Gouda, N. (2004). Universality in the distribution of caustics in the expanding Universe. The Astrophys. J. Suppl. 151 185.
  * (183) Yoshisato, A., Matsubara, T. and Morikawa, M. (1998). Accuracy of Nonlinear Approximations in Spheroidal Collapse: Why Are Zeldovich-Type Approximations Good? The Astrophys. J. 498, 48.
  * (184) Zel’dovich, Ya.B. (1970a). Gravitational instability: an approximate theory for large density perturbations. Astron. Astrophys. 5, 84.
  * (185) Zel’dovich, Ya.B. (1970b). Fragmentation of a homogeneous medium under the action of gravitation. Astrofizika 6, 319; Translated in Astrophysics 6, 164 (1973).
  * (186) Zel’dovich, Ya.B. (1978). The Large Scale Structure of the Universe. Eds.: Longair, M.S. and Einasto, J. In: IAU Symp. 79, 409. Springer: Berlin/Heidelberg, Germany. ISBN 978-94-009-9843-8
  * (187) Zel’dovich, Ya.B. and Myshkis, A.D. (1973). Elements of Mathematical Physics, Medium of Non-interacting Particles. Nauka, Moscow.
  * (188) Zel’dovich, Ya.B. and Shandarin, S.F. (1982). Maximum density in heavy-neutrino clouds. Pis’ma Astron. Zh. 8, 259; Translated in Sov. Astron. Lett. 8, 139 (1982).
  * (189) Zentsova, A.S. and Chernin, A.D. (1980). Evolution of entropy perturbations in the post-recombination epoch II – nonlinear stage. Astrofizika 16, 169; Translated in Astrophysics 16, 108.
  * (190) Zheligovsky, V. and Frisch, U. (2014). Time-analyticity of Lagrangian particle trajectories in ideal fluid flow. J. Fluid Mech. 749, 404. [arXiv:1312.6320]