# Beyond relativistic Lagrangian perturbation theory.
I. An exact-solution controlled model for structure formation

Ismael Delgado Gaspar<EMAIL_ADDRESS>Department of
Fundamental Research, National Centre for Nuclear Research, Pasteura 7, 02–093
Warsaw, Poland  Univ Lyon, Ens de Lyon, Univ Lyon1, CNRS, Centre de Recherche
Astrophysique de Lyon UMR5574, F-69007, Lyon, France  Thomas Buchert
<EMAIL_ADDRESS>Univ Lyon, Ens de Lyon, Univ Lyon1, CNRS, Centre de
Recherche Astrophysique de Lyon UMR5574, F-69007, Lyon, France  Jan J.
Ostrowski<EMAIL_ADDRESS>Department of Fundamental Research,
National Centre for Nuclear Research, Pasteura 7, 02–093 Warsaw, Poland

###### Abstract

We develop a new nonlinear method to model structure formation in general
relativity from a generalization of the relativistic Lagrangian perturbation
schemes, controlled by Szekeres (and LTB) exact solutions. The overall
approach can be interpreted as the evolution of a deformation field on an
inhomogeneous reference model, obeying locally Friedmann-like equations. In
the special case of locally one-dimensional deformations, the new model
contains the entire Szekeres family of exact solutions. As thus formulated,
this approach paraphrases the Newtonian and relativistic Zel’dovich
approximations, having a large potential for applications in contexts where
relativistic degrees of freedom are relevant. Numerical simulations are
implemented to illustrate the capabilities and accuracy of the model.

###### pacs:

04.20.-q, 04.20.Jb, 98.80.-k, 04.25.Nx

##  I Introduction

The problem of cosmological structure formation has been approached by a
variety of methods. The exact, inhomogeneous solutions of Einstein’s equations
provide valuable hints on relativistic effects, absent in Newtonian theory,
however, their applications remain limited as the actual cosmic web is much
more complex than the scope covered by the exact solutions. Alternatively, the
standard cosmological perturbation theory can be used in an attempt to trace
the evolution of the density inhomogeneities in the linear regime, the obvious
limit being the overdensities or second derivatives of metric perturbations
reaching relative values of around unity. These limitations could in principle
be overcome with the use of relativistic numerical simulations (see
2020Adameketal  for a recent comparison of performance among the most popular
codes currently being developed). There are, however, several drawbacks of
these numerical methods with a hard to estimate influence on the outcomes,
e.g., global weak field assumption, noncovariance of conformal decomposition,
all-time fixed toroidal topology or current computational capacity
limitations, to name a few. In light of recently revealed, as well as long
known but persistent, tensions haunting modern cosmology, it is very important
to improve on the analytical methods providing an essential counterpart to the
numerical efforts. In this spirit, the so-called silent universe model
(silent1 , silent2 ) was proposed and is believed to be able to describe
structure formation in the nonlinear regime, from a wide variety of initial
data, by an exact method. The silent universe model is based on the specific
restriction on the $1+3$ decomposition of Einstein’s equations for dust, i.e.
the gravitomagnetic part of the projected Weyl tensor is put to zero.
Unfortunately, this model turned out to be insufficient to access the
nonlinear regime (because of the absence of rotation and the requirement for
the shear to be diagonalizable and have two identical eigenvalues) and,
although as desired, it contains several known inhomogeneous exact solutions
as subcases, the span of admissible initial data is very restricted.

In this context, the recent investigation rza6  consolidated and generalized
an earlier insight by Kasai Kasai1995  on the correspondence between Szekeres
class II solutions and the first-order Lagrangian perturbation solutions in
relativistic cosmology. The relativistic Lagrangian perturbation schemes have
been developed in the series of papers rza1 , rza2 , rza3 , rza4 , rza5  and
rza6 . One aspect of this generalization is furnished by the result that an
extrapolated version of the first-order solution scheme, in the spirit of the
original proposal by Zel’dovich in Newtonian cosmology ZA , allows to extend
this correspondence to nonlinear functional expressions of the first-order
scheme, not only for the density but also for the bilinear metric form, the
extrinsic and intrinsic curvatures including their tracefree parts, and other
variables. The exact body of the functionally extrapolated perturbations is
obtained by setting the second and third principal scalar invariants of the
deformation field to zero (for details the reader is referred to rza6 ). In
the relativistic case, this corresponds to Szekeres class II solutions, while
in the Newtonian limit this corresponds to a class of 3D solutions without
symmetry obtained in buchertgoetz  (with empty background), in Buchert1989AA
for backgrounds with zero cosmological constant and in bildhaueretal
including a cosmological constant.

A further insight concerns the way we write the Szekeres class II solutions:
in the so-called Goode-Wainwright parametrization GW1 ; GW2  this class can be
written in the form of deviations off a global FLRW (Friedmann-Lemaître-
Robertson-Walker) background solution. Looking at the spatial average
properties of the inhomogeneous deviations, we find admissible initial data
for deviations that average out on this background solution. This can be
realized, as in Newtonian cosmological simulations, by setting periodic
boundary conditions on the deviation fields on some scale that is commonly
associated with a ‘scale of homogeneity’. The resulting architecture of such a
(relativistic) simulation has the topological structure of a flat $3$-torus
rza6 ; NajeraSussman , very similar to Newtonian simulations, implying
integrability of the deformation fields and an on average zero intrinsic
scalar curvature (for the notion of integrability of deformations, see rza4
and BuchertFocus ). The results on the correspondence with Szekeres class II
solutions carry over to average properties known for constructions with
$3$-torus topology, in Newtonian theory and in general relativity
BuchertEhlers ; buchert:newton . This can be summarized by the property of
zero cosmological backreaction on the scale set by the size of the toroidal
space (for Szekeres class II solutions, see rza6 ). In general, cosmological
backreaction can be nonzero CliftonSussman2019 , and the topology of spatial
sections enjoys rich possibilities in general relativity. For constant-
curvature models, there are relations to the topology of spatial sections, but
in general situations, the scalar curvature function is not tightly
constrained; it has a nonzero average and it can also change its sign during
the evolution buchertcarfora:curvature ; GBC . A different evolution of the
scalar curvature compared with the FLRW evolution of the constant curvature
also furnishes an explanation of the vividly discussed Hubble tension
hubbletension ; Bolejko:2018SilentUniv2  in the standard cosmological model.

In the present paper, we propose to generalize Lagrangian perturbation
schemes. In the example of the first-order scheme we exploit the structure of
Szekeres class I solutions that can be written as deviations off a ‘local
background’ obeying a Friedmann-type equation. These deviations can still be
modeled by the functional expressions developed in relativistic Lagrangian
perturbation theory. A ‘global background’ can be defined through a spatial
averaging operation of the full solution that leads to nonzero cosmological
backreaction, i.e. it allows not only for the impact of the background
evolution on the evolution of inhomogeneities as in Newtonian cosmologies,
relativistic Lagrangian schemes and standard quasi-Newtonian perturbation
theories, but also for the impact of inhomogeneities on the evolution of this
global background that is conceived as the average model. A numerical
implementation of this new model will open the door to answering many of the
questions raised in the context of inhomogeneous relativistic cosmology, and
it may provide a more general architectural setting for relativistic numerical
simulations.

We proceed as follows. In Section II we recall the Einstein equations in terms
of the Lagrangian coframe fields and the definition of the Relativistic
Lagrangian Zel’dovich Approximation (RZA). Section III is devoted to reviewing
the Szekeres models in the Goode and Wainwright parametrization and their
relation with RZA. In Section IV, we develop the proposed generalization of
relativistic Lagrangian perturbation schemes to include both classes of the
Szekeres solutions. Then, its most important subcases, namely the locally one-
dimensional solutions, LTB models and RZA, are discussed in Section V. In
Section VI, we present a family of simple models and implement numerical
simulations aimed at illustrating the capability of the approach to model
realistic cosmological structures. Finally, in Section VII, we put all the
elements of our analysis together, discuss their physical motivation and
conclude with a summary of the paper and final remarks.

The text is complemented by five appendices. Appendix A provides additional
details about the Goode and Wainwright parametrization. Then, in Appendix B,
we have a closer look at Szekeres class I solutions and discuss the
relationship between some of their common parametrizations. The model
equations are obtained in Appendix C. Appendix D contains useful expressions
to compute the kinematical backreaction term and its evolution. Finally, in
Appendix E, we look at the family of models examined in Section VI but from
the perspective of exact solutions.

##  II The relativistic Lagrangian formulation

In the $3+1$ relativistic Lagrangian framework, the Einstein equations’
dynamical freedom is completely encoded in the coframe functions
$\eta^{a}_{\>\>i}$ rza1 .111Indices $i,j,k,\cdots=1,2,3\cdots$ denote
coordinate indices, while indices $a,b,c\cdots=1,2,3\cdots$ are introduced as
counters of components, e.g., of vectors or differential forms. In this paper
we use units where the gravitational constant and the speed of light are set
to $G=c=1$. Restricting our attention to an irrotational dust source and
considering a fluid-flow orthogonal foliation of the spacetime leads to a set
of nine purely spatial coframes, in terms of which the line-element reads:

$^{(4)}\mathbf{g}=-\mathbf{d}t\otimes\mathbf{d}t+^{(3)}\\!\mathbf{g}\quad\hbox{with}\quad^{(3)}\mathbf{g}=G_{ab}\bm{\eta}^{a}\otimes\bm{\eta}^{b}\
.$ (1)

Here, the initial metric coefficients are encoded in the Gram’s matrix:

$G_{ij}(\mathbf{X})\equiv g_{ij}(t_{\textrm{\tiny{ini}}},\mathbf{X})\ ,$ (2)

where $t_{\textrm{\tiny{ini}}}$ denotes some arbitrary initial time.

###  II.1 The Lagrange-Einstein system

The Einstein equations are rewritten as a system of $9$ evolution equations
and $4$ constraint equations rza1 ; rza4 ,

$\displaystyle\hskip 70.0001ptG_{ab}\,\dot{\eta}^{a}_{[i}\eta^{b}_{\ j]}=0\ ;$
(3a)
$\displaystyle\frac{1}{2J}\epsilon_{abc}\epsilon^{ikl}\left(\dot{\eta}^{a}_{\
j}\eta^{b}_{\ k}\eta^{c}_{\ l}\right)^{\hbox{\bigastfont.}}=-{\cal{R}}^{i}_{\
j}+\left(4\pi\varrho+\Lambda\right)\delta^{i}_{\ j}\ ;$ (3b)
$\displaystyle\frac{1}{2J}\epsilon_{abc}\epsilon^{mjk}\dot{\eta}^{a}_{\
m}\dot{\eta}^{b}_{\ j}\eta^{c}_{\
k}=-\frac{{\cal{R}}}{2}+\left(8\pi\varrho+\Lambda\right)\ ;$ (3c)
$\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}\ .$ (3d)

We call this system the Lagrange-Einstein system. In the equations above, as
throughout the text, the overdot represents the partial time derivative (the
covariant time-derivative in the flow-orthogonal foliation), while the single
and the double vertical slashes stand for the partial and covariant spatial
derivatives, respectively. ${\cal{R}}_{ij}$ denotes the spatial Ricci tensor,
${\cal{R}}={\cal{R}}^{i}_{\>\>i}$ the spatial scalar curvature, and the
determinant of the $3\times 3$ coframe matrix is defined as:

$J\equiv\det(\eta^{a}_{\ i})\ .$ (4)

The exact density field follows from the integration of the continuity
equation,

$\varrho=\varrho_{\textrm{\tiny{ini}}}J^{-1}\ ,\quad\hbox{with}\ \quad
J=\sqrt{g}/\sqrt{G}\ ,$ (5)

and $\varrho_{\textrm{\tiny{ini}}}=\varrho(t_{\textrm{\tiny{ini}}})$.222In
this paper we refer to the initial quantities (evaluated at
$t_{\textrm{\tiny{ini}}}$) with the subscript ‘ini’. The expansion tensor can
be computed considering its relation to the extrinsic curvature, which in the
assumed foliation of the spacetime takes the following form:

$\displaystyle\Theta_{ij}$ $\displaystyle=$
$\displaystyle-\mathcal{K}_{ij}=\frac{1}{2}\dot{g}_{ij}\ ;$ (6)
$\displaystyle\Theta^{i}_{\>\>j}=e^{i}_{\>\>a}\,\dot{\eta}^{a}_{\>\>j}\ ,$
$\displaystyle\text{with}\quad
e^{i}_{\>\>a}=\frac{1}{2J}\epsilon_{abc}\epsilon^{ikl}\eta^{b}_{\
k}\eta^{c}_{\ l}\ .$ (7)

For irrotational dust, the expansion tensor consists of only the expansion
scalar (its trace part, $\Theta=\dot{J}/J$) and the shear tensor (its
tracefree symmetric part, $\sigma_{ij}$),

$\Theta^{i}_{\>\>j}=\sigma^{i}_{\>\>j}\,+\,\frac{1}{3}\Theta\,\delta^{i}_{\>\>j}\
.$ (8)

The 3D intrinsic curvature and the gravitoelectric and gravitomagnetic parts
of the spatially projected Weyl curvature tensor reduce to the following
expressions in terms of the coframes rza4 :

$\displaystyle-{\cal{R}}^{i}_{\ j}$ $\displaystyle=$
$\displaystyle\frac{1}{2J}\epsilon_{abc}\epsilon^{ikl}\left(\dot{\eta}^{a}_{\
j}\eta^{b}_{\ k}\eta^{c}_{\
l}\right)^{\hbox{\bigastfont.}}-\left(4\pi\varrho+\Lambda\right)\delta^{i}_{\
j}\ ;\qquad$ (9a) $\displaystyle-\frac{{\cal{R}}}{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2J}\epsilon_{abc}\epsilon^{mjk}\dot{\eta}^{a}_{\
m}\dot{\eta}^{b}_{\ j}\eta^{c}_{\ k}-\left(8\pi\varrho+\Lambda\right)\ ;$ (9c)
$\displaystyle\frac{1}{2J}\epsilon_{abc}\epsilon^{ik\ell}\ddot{\eta}^{a}_{\
i}\eta^{b}_{\ k}\eta^{c}_{\ \ell}=\Lambda-4\pi{\varrho}\ ;$
$\displaystyle-E^{i}_{\>\>j}$ $\displaystyle=$
$\displaystyle\frac{1}{2J}\epsilon_{abc}\epsilon^{ikl}_{\>\>}\ddot{\eta}^{a}_{\>\>j}\eta^{b}_{\>\>k}\eta^{c}_{\>\>l}+\frac{1}{3}\Big{(}4\pi\varrho-\Lambda\Big{)}\delta^{i}_{\>\>j}\
;$ (9d) $\displaystyle-H^{i}_{\>\>j}$ $\displaystyle=$
$\displaystyle\frac{1}{J}G_{ab}\epsilon^{ikl}_{\>\>}\big{(}\dot{\eta}^{a}_{\>\>j\parallel
l}\eta^{b}_{\>\>k}+\dot{\eta}^{a}_{\>\>j}\eta^{b}_{\>\>k\parallel l}\big{)}\
.$ (9e)

###  II.2 Relativistic Zel’dovich Approximation (RZA)

The relativistic Lagrangian perturbation theory consists of perturbing the
trivial coframe set associated with a homogeneous and isotropic spacetime
$\eta^{a}_{\>\>i}=a(t)\delta^{a}_{\>\>i}$, where $a(t)$ is the solution of the
Friedmann equations. Then, RZA emerges as the first-order perturbation of the
deformation field, $P^{a}_{~{}i}$:

$\bm{\eta}^{a}=\eta^{a}_{~{}i}\mathbf{d}X^{i}=a(t)\left(\delta^{a}_{~{}i}+P^{a}_{~{}i}\right)\mathbf{d}X^{i}\
.$ (10)

In this approximation, the line-element (1) takes the following quadratic
bilinear form:

$g_{ij}=a^{2}(t)\left[G_{ij}+G_{ab}\left(\delta^{a}_{~{}i}P^{b}_{~{}j}+\delta^{b}_{~{}j}P^{a}_{~{}i}+P^{a}_{~{}i}P^{b}_{~{}j}\right)\right]\
.$ (11)

Consequently, any relevant field (the Ricci and Weyl curvature tensors, rate
of expansion, shear, etc.) is functionally evaluated in terms of the coframe
fields (10).

##  III Szekeres exact solutions

The Szekeres models are the most general exact cosmological solution of
Einstein’s equations kras1 ; BolejkoInhCosModels .333 This statement requires
some justification. In the present analysis we follow BolejkoInhCosModels  and
consider as cosmological models those exact solutions of Einstein’s equations
with at least a nontrivial subclass of FLRW solutions as a limit. There are
other solutions that can be regarded as more general than Szekeres, for
example, the Lemaître model Lemaitre:1933gd  (although it is spherically
symmetric the source is an inhomogeneous perfect fluid) and other models
containing heat-flow Baysal1992Heatflow ; Goode1986HeatFlow ;
LimaMaia1985HeatFlow ; NajeraSussman , viscosity TomimuraMotta1990aViscosity ;
MottaTomimura1990bViscosity  or electromagnetic fields LimaNobre1989EMfields ;
TomimuraWaga1987EMfields . However, in cosmological applications, the field
source is usually simplified to a dust fluid and the cosmological constant;
then, Szekeres models arise as the most general (exact) cosmological class of
solutions. These models lack symmetries, albeit their ‘quasisymmetries’ impose
significant restrictions on the spacetime’s inhomogeneity and anisotropy. The
source of the original solution consists of irrotational dust Sz75 , which was
later generalized by Szafron szafron1977inhomogeneous  to include a
homogeneous perfect fluid, furnishing the inclusion of the cosmological
constant (see CommentSzeHScaFields  for a discussion of the Szekeres source).
For a detailed review of the properties of these solutions see kras1 ; kras2 ;
BKHC2009 . The general line-element is commonly expressed as (see
wainwright1977characterization ; szafronCollins79 ; BarnesRowlingson1989SzeInv
for coordinate-independent definitions of the Szekeres solution):

${\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)\
,$ (12)

where $\alpha=\alpha(t,\bm{r})$, $\beta=\beta(t,\bm{r})$, and
$\bm{r}=\left(x,y,z\right)$ are comoving coordinates. The integration of the
Einstein equations is performed by splitting the model into two classes,
namely class I, when $\beta_{,z}\neq 0$, and class II, when $\beta_{,z}=0$.
However, once the solution is obtained, the class II can be formulated as a
limit of class I hellaby1996null .

In the Goode and Wainwright (GW) parametrization GW1  the line-element (12) is
rewritten as:

${\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)\
,$ (13)

where the conformal scale factor, $S(t,z)$, obeys a Friedmann-like equation:

$\dot{{\cal{S}}}^{2}=-k_{0}+\frac{2\mu}{{\cal{S}}}+\frac{\Lambda}{3}{\cal{S}}^{2}\
.$ (14)

Here, $k_{0}$ is a constant taking the values $0,\pm 1$, and $\mu(z)>0$ is an
arbitrary scalar function. Furthermore,

${\cal{G}}={\cal{A}}(\bm{r})-{\cal{F}}(t,z)={\cal{A}}(\bm{r})-\beta_{+}f_{+}-\beta_{-}f_{-}\
,$ (15)

where $f_{\pm}$ are the growing and decaying solutions of the following
equation:

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

The energy-density is given by:

$8\pi\varrho(t,\bm{r})=\frac{6\mu}{{\cal{S}}^{3}}\left(1+\frac{{\cal{F}}}{{\cal{G}}}\right)\
,$ (17)

which can be rewritten as:

$4\pi\varrho(t,\bm{r})=\frac{3\mu{\cal{A}}}{{\cal{S}}^{3}{\cal{G}}}\equiv\frac{M(\bm{r})}{{\cal{S}}^{3}{\cal{G}}}\
,$ (18)

assigning to the term $3\mu{\cal{A}}$ the meaning of a conserved rest mass,
$M$. The functions $f_{\pm}$ can be regarded as the growing and decaying
deviation ‘modes’ at a ‘fictitious background’ (or reference model) with
initial local density $3\mu=4\pi\rho$.

The form of the functions ${\cal{A}}$, $e^{2\nu}$, $\mathcal{W}$,
$\beta_{+}(z)$, $\beta_{-}(z)$, $f_{+}$ and $f_{-}$ vary depending on the
class, and they satisfy:

Class I

${\cal{S}}={\cal{S}}(t,z)$, $\mu=\mu(z)$ and $f_{\pm}=f_{\pm}(t,z)$.

Class II

${\cal{S}}={\cal{S}}(t)\equiv a(t)$, $\mu=\text{const.}$ and
$f_{\pm}=f_{\pm}(t)$.

For better readability of the text, we display their full expressions in
Appendix A.

The quasispherical branch of the class I solutions has been the most widely
used subclass for modeling structure formation, describing nontrivial networks
of cosmic structures Bolejko2006Struformation ; IshakNwankwo2008 ;
Bolejko:2009GERG ; BolCelerier2010 ; BolSuss2011 ; Nwankwo_2011 ; KrasBol2011
; Buckley2013 ; IshakPRL2013 ; Vrba:2014 ; Koksbang2015 ; Koksbang2015II ;
sussman2015multiple ; sussman2016coarse ; Koksbang2017  or even the formation
of primordial black holes HaradaPBHExactSol ; gaspar2018black ; Coley:2019ylo
. Due to its importance, the physical and mathematical properties of this
subclass have also been explored in depth BonnorTomimura1976evolution ;
Bonnor1977 ; HellabyKrasi2002 ; Hellaby:2007hq ; Hellaby:2007hq ;
WaltersHellaby2012 ; NonsphericalSzePerts ; GIraHellaby2017 ; Hellaby2017 .
The other subclasses have found applications in cosmology and astrophysics as
well ishak2012growth ; IshakPeel2012LargeScale ; SzeLamddanot0Meures ; Meures-
Bruni_mnras , although they have received much less attention.

###  III.1 Szekeres class I models with a normalized scale factor

In contrast to class II, the class I solution’s conformal scale factor is a
function of the spatial coordinates,
${\cal{S}}(t_{\textrm{\tiny{ini}}},z)={\cal{S}}_{\textrm{\tiny{ini}}}(z)$. To
have a normalized initial function in a given Szekeres model, we proceed as
follows. First, introduce the function
$\chi=S^{-1}(t_{\textrm{\tiny{ini}}},\bm{r})>0$; then, the conformal scale
factor and the spatial metric are redefined:

$A(t,\bm{r})\equiv{\cal{S}}(t,z)\chi\ ;$ (19)
$g_{ij}=A^{2}~{}\mbox{Diag}\left[\left({\cal{A}}-{\cal{F}}\right)^{2}\left(\frac{{\cal{W}}}{\chi}\right)^{2},\left(\frac{e^{\nu}}{\chi}\right)^{2},\left(\frac{e^{\nu}}{\chi}\right)^{2}\right]\
,$ (20)

where metric functions are displayed in Appendix A.

In terms of the rescaled scale factor, the Friedmann-like equation reads:

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

where:

$\hat{k}(\bm{r})=k_{0}\chi^{2}\
,\quad\frac{4\pi}{3}\hat{\varrho}_{b}(\bm{r})=\mu(z)\chi^{3}\ .$ (22)

Note that the growing and decaying functions still refer to the local
reference background (14), and obey the same equation as in (16):

$\ddot{{\cal{F}}}+2\frac{\dot{A}}{A}\dot{{\cal{F}}}-\frac{4\pi\hat{\varrho}_{b}}{A^{3}}{\cal{F}}=0\
.$ (23)

###  III.2 Relationship between RZA and exact solutions

In rza6  we exploited such a splitting into ‘background’ and the growing and
decaying ‘deviation modes’ to connect the Szekeres class II solution to RZA.
The main results of interest for the present article can be summarized as
follows:

  * $(i)$

RZA contains the class II of the Szekeres exact solutions as a particular case
with the identifications:

$\displaystyle P^{a}_{\>\>i}$ $\displaystyle=$
$\displaystyle(-\widetilde{{\cal{F}}}/\widetilde{{\cal{A}}})\,\delta^{a}_{\>\>3}\delta^{3}_{\>\>i}\
;$ (24a) $\displaystyle G_{ab}$ $\displaystyle=$
$\displaystyle\hbox{Diag}\left[\widetilde{{\cal{A}}}^{2},e^{2\nu},e^{2\nu}\right]\
.$ (24b)

In these equations,
$\widetilde{{\cal{A}}}\equiv{\cal{A}}-{\cal{F}}_{\textrm{\tiny{ini}}}$,
$\widetilde{{\cal{F}}}\equiv{\cal{F}}-{\cal{F}}_{\textrm{\tiny{ini}}}$. Also,
$a(t)$ is the scale factor of the global FLRW background with
$4\pi\varrho_{b}/3=\mu=\mbox{const.}$

  * $(ii)$

The class I can be interpreted as a (constrained) superposition of
nonintersecting world lines, satisfying RZA’s evolution equations. As for the
previous class, the associated deformation field and Gram’s matrix are given
by the identifications:

$\displaystyle\hat{P}^{a}_{\>\>i}$ $\displaystyle=$
$\displaystyle(-\widetilde{{\cal{F}}}/\widetilde{{\cal{A}}})\,\delta^{a}_{\>\>3}\delta^{3}_{\>\>i}\
;$ (25a) $\displaystyle G_{ab}$ $\displaystyle=$
$\displaystyle\hbox{Diag}\left[\widetilde{{\cal{A}}}^{2}\left(\frac{{\cal{W}}}{\chi}\right)^{2},\left(\frac{e^{\nu}}{\chi}\right)^{2},\left(\frac{e^{\nu}}{\chi}\right)^{2}\right]\
,$ (25b)

with the proviso that these expressions need to be evaluated for each world
line independently. Given an arbitrary world line with fixed comoving
coordinates $\bm{r}_{i}=\left(x_{i},y_{i},z_{i}\right)$, $A(t)\equiv
A(t,z_{i})$ can be regarded as the scale function of an associated FLRW ‘local
background’ or ‘reference model’ with initial constant density
$\varrho_{b}(t_{\textrm{\tiny{ini}}})=\hat{\varrho}_{b}(z_{i})$ and curvature
$k_{0}=\hat{k}(z)$. Thus, $f_{\pm}(t)\equiv f_{\pm}(t,z)|_{z=z_{i}}$ are its
associated growing and decaying modes.

##  IV Generalized Relativistic Zel’dovich Approximation (GRZA)

Motivated by the mathematical connection between the Szekeres class II
solutions and RZA, we present a new nonperturbative approach that generalizes
RZA to contain Szekeres class I. The approach retains the mathematical
structure of RZA but without pre-assuming a global background. Instead, we
consider a space-dependent conformal scale factor obeying a Friedmann-like
equation (as in Eq. (22) for Szekeres class I).

Restricting the analysis to an irrotational dust source, the coframe set (10)
finds its generalization in the following expressions:

$\displaystyle\bm{\eta}^{a}=\eta^{a}_{~{}i}\mathbf{d}X^{i}=A\left(\delta^{a}_{~{}i}+\hat{P}^{a}_{\>\>i}\right)\mathbf{d}X^{i}\
;$ (26a) $\displaystyle A=A\left(t,\bm{r}\right)\
,\quad\hat{P}^{a}_{\>\>i}=\hat{P}^{a}_{\>\>i}\left(t,\bm{r}\right)\ ,$ (26b)

where, as discussed above, $A$ satisfies a Friedmanian equation for a
reference model with curvature and matter density parameters
$\hat{k}=\hat{k}\left(\bm{r}\right)$ and $\hat{\varrho}_{b}(\bm{r})$,
respectively:

$\displaystyle 2\ddot{A}/A+\dot{A}^{2}/A^{2}-\Lambda+\hat{k}/A^{2}=0\ ;$ (27a)
$\displaystyle\left(\frac{\dot{A}}{A}\right)^{2}=-\frac{\hat{k}}{A^{2}}+\frac{8\pi}{3}\frac{\hat{\varrho}_{b}}{A^{3}}+\frac{\Lambda}{3}\
.$ (27b)

The line-element keeps the bilinear quadratic mathematical structure for the
deformation field, as in Eq. (1):

$\displaystyle g_{ij}$ $\displaystyle=$ $\displaystyle
G_{ab}\eta^{a}_{~{}i}\eta^{b}_{~{}j}$ (28a) $\displaystyle=$ $\displaystyle
A^{2}\left[G_{ij}+G_{ab}\left(\delta^{a}_{~{}i}\hat{P}^{b}_{~{}j}+\delta^{b}_{~{}j}\hat{P}^{a}_{~{}i}+\hat{P}^{a}_{~{}i}\hat{P}^{b}_{~{}j}\right)\right]\
.\qquad$ (28b)

Since at the initial time $A$ is normalized and we assume (without loss of
generality) that the initial deformation field vanishes
($\hat{P}^{a}_{~{}i}(t_{\textrm{\tiny{ini}}})=0$), the Gram’s matrix is
defined via the initial spatial metric:

$G_{ij}(\bm{r}):=g_{ij}(t_{i},\bm{r})\ .$ (29)

The subsequent approach consists of (i) obtaining (Lagrange-)linear evolution
equations for the deformation field and, then, (ii) injecting the formal
solution into the exact nonlinear functional expressions. This scheme retains
the original Zel’dovich’s extrapolation idea and generalizes RZA to include
the whole family of Szekeres models within its locally one-dimensional
deformation field limit. In light of this, we call the resulting approach
Generalized Relativistic Zel’dovich Approximation (GRZA).

###  IV.1 Functional evaluation

First, let us evaluate all relevant fields as exact functionals of the local
deformation and the conformal scale factor. The exact determinant of the
spatial coframe coefficients is given by:

$\displaystyle
J=A^{3}\mathfrak{J}=A^{3}\left(\mathfrak{J}_{0}+\mathfrak{J}_{1}+\mathfrak{J}_{2}+\mathfrak{J}_{3}\right)\
,$ (30)

where we have introduced the peculiar-determinant $\mathfrak{J}$, and the
quantities $\mathfrak{J}_{n}\equiv\mbox{}^{(n)\\!}\mathfrak{J}^{k}_{\>\>k}$
are defined through:

$\displaystyle\mbox{}^{(0)\\!}\mathfrak{J}^{i}_{\>\>j}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\epsilon_{abc}\epsilon^{ikl}_{\>\>}\delta^{a}_{\>\>j}\delta^{b}_{\>\>k}\delta^{c}_{\>\>l}\
;$ (31a) $\displaystyle\mbox{}^{(1)\\!}\mathfrak{J}^{i}_{\>\>j}$
$\displaystyle=$
$\displaystyle\frac{1}{6}\epsilon_{abc}\epsilon^{ikl}_{\>\>}\hat{P}^{a}_{\>\>j}\delta^{b}_{\>\>k}\delta^{c}_{\>\>l}+\frac{1}{3}\epsilon_{abc}\epsilon^{ikl}_{\>\>}\delta^{a}_{\>\>j}\hat{P}^{b}_{\>\>k}\delta^{c}_{\>\>l}\
;$ (31b) $\displaystyle\mbox{}^{(2)\\!}\mathfrak{J}^{i}_{\>\>j}$
$\displaystyle=$
$\displaystyle\frac{1}{3}\epsilon_{abc}\epsilon^{ikl}_{\>\>}\hat{P}^{a}_{\>\>j}\hat{P}^{b}_{\>\>k}\delta^{c}_{\>\>l}+\frac{1}{6}\epsilon_{abc}\epsilon^{ikl}_{\>\>}\delta^{a}_{\>\>j}\hat{P}^{b}_{\>\>k}\hat{P}^{c}_{\>\>l}\
;\qquad$ (31c) $\displaystyle\mbox{}^{(3)\\!}\mathfrak{J}^{i}_{\>\>j}$
$\displaystyle=$
$\displaystyle\frac{1}{6}\epsilon_{abc}\epsilon^{ikl}_{\>\>}\hat{P}^{a}_{\>\>j}\hat{P}^{b}_{\>\>k}\hat{P}^{c}_{\>\>l}\
.$ (31d)

From its definition, we can see that $\mathfrak{J}_{0}=1$, as expected.
Henceforth, the exact inhomogeneous density field follows from injecting (30)
and (31) into (5).

Next, we express the expansion tensor as a functional of the deformation
field. Writing out Eq. (7) yields:

$\Theta^{i}_{\>\>j}=3\hat{H}\frac{\sum_{n}\mbox{}^{(n)\\!}\mathfrak{J}^{i}_{\>\>j}}{\mathfrak{J}}+\frac{1}{\mathfrak{J}}\left(\mbox{}^{(1)\\!}\overline{I}^{i}_{\>\>j}+\mbox{}^{(2)\\!}\overline{I}^{i}_{\>\>j}+\mbox{}^{(3)\\!}\overline{I}^{i}_{\>\>j}\right)\
.$ (32)

To keep the above expression short, we have introduced the shorthand
notations:

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle\dot{A}/A\ ;$ (33a)
$\displaystyle\mbox{}^{(1)\\!}\overline{I}^{i}_{\>\>j}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\epsilon_{abc}\epsilon^{ikl}_{\>\>}\dot{\hat{P}}^{a}_{\>\>j}\delta^{b}_{\>\>k}\delta^{c}_{\>\>l}\
;$ (33b) $\displaystyle\mbox{}^{(2)\\!}\overline{I}^{i}_{\>\>j}$
$\displaystyle=$
$\displaystyle\epsilon_{abc}\epsilon^{ikl}_{\>\>}\dot{\hat{P}}^{a}_{\>\>j}\hat{P}^{b}_{\>\>k}\delta^{c}_{\>\>l}\
;$ (33c) $\displaystyle\mbox{}^{(3)\\!}\overline{I}^{i}_{\>\>j}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\epsilon_{abc}\epsilon^{ikl}_{\>\>}\dot{\hat{P}}^{a}_{\>\>j}\hat{P}^{b}_{\>\>k}\hat{P}^{c}_{\>\>l}\
.$ (33d)

Taking the trace of the expansion tensor and using a similar notation as
before, $\overline{I}_{n}\equiv\mbox{}^{(n)\\!}\overline{I}^{k}_{\>\>k}$, we
obtain the functional for the expansion scalar:

$\Theta=3\hat{H}+\frac{1}{\mathfrak{J}}\left(\overline{I}_{1}+\overline{I}_{2}+\overline{I}_{3}\right)\
.$ (34)

Then, the exact functional for the shear tensor can be computed from (8), (32)
and (34).

The gravitoelectric part of the spatially projected Weyl tensor reads:

$\displaystyle E^{i}_{\>\>j}$ $\displaystyle=$
$\displaystyle-3\frac{\ddot{A}}{A}\frac{\sum_{n}\mbox{}^{(n)\\!}\mathfrak{J}^{i}_{\>\>j}}{\mathfrak{J}}-2\hat{H}\frac{\sum_{n}\mbox{}^{(n)\\!}\overline{I}^{i}_{\>\>j}}{\mathfrak{J}}-\frac{\sum_{n}\mbox{}^{(n)\\!}\widehat{I}^{i}_{\>\>j}}{\mathfrak{J}}$
(35)
$\displaystyle\qquad\qquad\qquad\quad-\frac{1}{3}\left(4\pi\varrho-\Lambda\right)\delta^{i}_{\>\>j}\
.$

Above, the quantities $\mbox{}^{(n)\\!}\widehat{I}^{i}_{\>\>j}$ are defined as
$\mbox{}^{(n)\\!}\overline{I}^{i}_{\>\>j}$ in (33), but replacing
$\dot{\hat{P}}^{a}_{\>\>j}$ by $\ddot{\hat{P}}^{a}_{\>\>j}$. The
gravitomagnetic part is given by

$\displaystyle
H^{i}_{\>\>j}=-\frac{\epsilon^{ikl}_{\>\>}}{A\mathfrak{J}}\Bigg{\\{}\left(\frac{\dot{A}_{|l}}{A}+\hat{H}\frac{A_{|l}}{A}\right)h_{jk}$
$\displaystyle\qquad\qquad+\hat{H}\left(G_{ja}\hat{P}^{a}_{\>\>k||l}+G_{ka}\hat{P}^{a}_{\>\>j||l}\right)$
$\displaystyle\qquad\qquad+2\frac{A_{|l}}{A}\left(G_{ka}\dot{\hat{P}}^{a}_{\>\>j}+G_{ab}\hat{P}^{b}_{\>\>k}\dot{\hat{P}}^{a}_{\>\>j}\right)+G_{ka}\dot{\hat{P}}^{a}_{\>\>j||l}\qquad$
$\displaystyle\qquad\qquad+G_{ab}\left[H\left(\hat{P}^{a}_{\>\>j}\hat{P}^{b}_{\>\>k}\right)_{||l}+\left(\dot{\hat{P}}^{a}_{\>\>j}\hat{P}^{b}_{\>\>k}\right)_{||l}\right]\Bigg{\\}}\
.$ (36)

In the first line of the above equation, $h_{jk}$ denotes the conformal
spatial metric $h_{jk}=g_{ij}/A^{2}$, which coincides with the Gram’s matrix
at the initial time (recall that $A_{\textrm{\tiny{ini}}}=1$).

Notice that (36) is a complex expression that in general does not identically
vanish, except in special cases (like in the Szekeres limit). This suggests a
nonsilent character of the GRZA solutions (unless special initial data are
imposed) and paves the way for applications to the physics of gravitational
waves, see rza4  for a study of the Lagrangian approach to gravitational
waves.

Finally, the functional for the spatial Riemann curvature tensor can be
obtained from its well-known expression in terms of partial derivatives of the
metric tensor.

###  IV.2 First-order perturbation scheme

To determine the first-order scheme, we consider the linearized Lagrange-
Einstein system at a local (Lagrangian) background solution. Linearizing only
in the deformation field, i.e., neglecting the terms of second- or higher-
order in $\hat{P}^{a}_{\>\>i}$, we obtain in line with rza4 :

$\displaystyle\left(\frac{\dot{A}}{A}\right)^{2}=-\frac{\hat{k}}{A^{2}}+\frac{8\pi}{3}\frac{\hat{\varrho}_{b}}{A^{3}}+\frac{\Lambda}{3}\
;$ (37a) $\displaystyle\qquad\qquad\dot{\hat{P}}_{[ij]}=0\ ;\quad$ (37b)
$\displaystyle\ddot{\hat{P}}^{i}_{\>\>j}+3\hat{H}\dot{\hat{P}}^{i}_{\>\>j}=-\mathfrak{J}\left[{\cal{R}}{}^{i}_{\>\>j}-\frac{1}{4}\left({\cal{R}}+2\frac{\hat{k}}{A^{2}}\right)\delta^{i}_{\>\>j}\right]\
;\quad$ (37c)
$\displaystyle\hat{H}\dot{\hat{P}}+\frac{4\pi\hat{\varrho}_{b}}{A^{3}}\hat{P}=-\frac{\mathfrak{J}}{4}\left({\cal{R}}-6\frac{\hat{k}}{A^{2}}\right)-\frac{\hat{W}}{A^{3}}\
;$ (37d)
$\displaystyle\frac{\dot{A}_{|j}}{A}-\hat{H}\frac{A_{|j}}{A}+\frac{1}{\mathfrak{J}}\dot{\hat{P}}^{i}_{\>\>[i|j]}+\zeta_{j}=0\
,$ (37e)

where $\mathfrak{J}\equiv J\,A^{-3}$, $\hat{H}$ is the generalized Hubble
function, and $\hat{W}$ the trace of the relativistic analog of the Newtonian
peculiar-acceleration gradient, defined together with $\zeta_{j}$ as follows:

$\hat{W}\equiv-4\pi\hat{\varrho}_{b}\,\hat{\delta}_{\textrm{\tiny{ini}}}\
,\;\;\zeta_{j}\equiv\frac{1}{2\mathfrak{J}}\left(\Gamma^{i}_{\>\>lj}\dot{\hat{P}}^{l}_{\>\>i}-\Gamma^{i}_{\>\>il}\dot{\hat{P}}^{l}_{\>\>j}\right)\
.$ (38)

In the equation above, the Christoffel symbols should be of zeroth-order in
the deformation field; otherwise, $\zeta_{j}$ would be of second or higher
order. (See Appendix C.2 for the explicit expression.) We have opted for
keeping the peculiar determinant in equations (37)-(38). The importance of
this term will be examined elsewhere. Below we will apply a strict
linearization and assume that the curvature terms on the right-hand side of
(37c) and (37d) are small and of the order of $\hat{P}$.

Following the steps of rza4 , we decompose the deformation field into its
trace, tracefree symmetric, and antisymmetric parts:

$\hat{P}_{ij}=\hat{P}_{(ij)}+\hat{P}_{[ij]}=\frac{1}{3}\hat{P}G_{ij}+\hat{\Pi}_{ij}+\hat{{\mathfrak{P}}}_{ij}\
,$ (39a) with the definitions:
$\hat{\Pi}_{ij}\equiv\hat{P}_{(ij)}-\frac{1}{3}\hat{P}G_{ij}\
;\quad\hat{{\mathfrak{P}}}_{ij}\equiv\hat{P}_{[ij]}\ .$ (39b)

The linearized system reads:

$\displaystyle\dot{\hat{{\mathfrak{P}}}}_{ij}=0\ ;$ (40a)
$\displaystyle\ddot{\hat{P}}+3H\dot{\hat{P}}=-\frac{1}{4}\left({\cal{R}}-6\frac{\hat{k}}{A^{2}}\right)\
;$ (40b)
$\displaystyle\ddot{\hat{\Pi}}_{ij}+3\hat{H}\dot{\hat{\Pi}}_{ij}=-\frac{\tau_{ij}}{A^{2}}\
;$ (40c)
$\displaystyle\hat{H}\dot{\hat{P}}+\frac{4\pi\hat{\varrho}_{b}}{A^{3}}\hat{P}=-\frac{1}{4}\left({\cal{R}}-6\frac{\hat{k}}{A^{2}}\right)-\frac{\hat{W}}{A^{3}}\
;\;$ (40d)
$\displaystyle\frac{\dot{A}_{|j}}{A}-\hat{H}\frac{A_{|j}}{A}+\frac{1}{3}\dot{\hat{P}}_{|j}-\frac{1}{2}\dot{\hat{\Pi}}^{i}_{\>\>j|i}=-\zeta_{j}\
,$ (40e)

where $\tau_{ij}$ denotes the tracefree part of the spatial Ricci curvature,

$\tau_{ij}\equiv{\cal{R}}_{ij}-\frac{1}{3}{\cal{R}}\,g_{ij}^{\scriptscriptstyle{(0)}}\
,\quad\hbox{with}\quad g_{ij}^{\scriptscriptstyle{(0)}}=A^{2}G_{ij}\ .$ (41)

Keeping only the terms linear in the deformation, Eq. (38) simplifies to:

$\displaystyle\zeta_{j}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{A_{|j}}{A}\dot{\hat{P}}-3\frac{A_{|l}}{A}\dot{\hat{P}}^{l}_{\>\>j}\right)\qquad$
(42)
$\displaystyle\quad\qquad-\frac{1}{2}\left(\Gamma^{i}_{\>\>il}(\mathbf{G})\dot{\hat{P}}^{l}_{\>\>j}-\Gamma^{l}_{\>\>ij}(\mathbf{G})\dot{\hat{P}}^{i}_{\>\>l}\right)\
.$

In the above equation, $\Gamma^{i}_{\>\>jk}(\mathbf{G})$ denotes the
Christoffel symbols for the initial spatial metric, see Eq. (98).

Substituting (40d) into (40b), we obtain the master equation for the trace
part:

$\ddot{\hat{P}}+2\hat{H}\dot{\hat{P}}-\frac{4\pi\hat{\varrho}_{b}}{A^{3}}\hat{P}=\frac{\hat{W}}{A^{3}}\
,$ (43)

while the traceless symmetric part obeys (40c). Both coincide in form with
their original RZA counterparts (see Eqs. (41) and (24) of rza4 ).

###  IV.3 General initial data setup

To prescribe the initial data, we will proceed as in the standard RZA and
identify the time-derivatives of the deformation field with coefficients of
the (nonintegrable) generalizations of the Newtonian peculiar-velocity
gradient, $\hat{U}^{a}_{\>\>i}$, and peculiar-acceleration gradient,
$\hat{W}^{a}_{\>\>i}$, respectively:

$\displaystyle\hat{P}^{a}_{\>\>i}(t_{\textrm{\tiny{ini}}})$ $\displaystyle=$
$\displaystyle 0\ ;$ (44a)
$\displaystyle\dot{\hat{P}}^{a}_{\>\>i}(t_{\textrm{\tiny{ini}}})$
$\displaystyle=$ $\displaystyle\hat{U}^{a}_{\>\>i}\ ;$ (44b)
$\displaystyle\ddot{\hat{P}}^{a}_{\>\>i}(t_{\textrm{\tiny{ini}}})$
$\displaystyle=$
$\displaystyle\hat{W}^{a}_{\>\>i}-2\hat{H}_{\textrm{\tiny{ini}}}\hat{U}^{a}_{\>\>i}\
.$ (44c)

The initial expansion rate and shear can be computed from the trace of the
generalization of the Newtonian peculiar-velocity gradient, $\hat{U}$:

$\displaystyle{\cal{H}}_{\textrm{\tiny{ini}}}$ $\displaystyle=$
$\displaystyle\frac{\Theta_{\textrm{\tiny{ini}}}}{3}=\hat{H}_{\textrm{\tiny{ini}}}+\frac{\hat{U}}{3}\
;$ (45) $\displaystyle\sigma^{i}_{\>\>j}(t_{\textrm{\tiny{ini}}})$
$\displaystyle=$
$\displaystyle\hat{U}^{i}_{\>\>j}-\frac{\hat{U}}{3}\delta^{i}_{\>\>j}\ ,$ (46)

where
$\hat{H}_{\textrm{\tiny{ini}}}=\dot{A}_{\textrm{\tiny{ini}}}/A_{\textrm{\tiny{ini}}}=\dot{A}_{\textrm{\tiny{ini}}}$
since $A_{\textrm{\tiny{ini}}}=1$. By construction
$\hat{H}_{\textrm{\tiny{ini}}}$ satisfies a Friedmann-like equation,

$\hat{H}_{\textrm{\tiny{ini}}}^{2}=\dot{A}_{\textrm{\tiny{ini}}}^{2}=-\hat{k}(\bm{r})+\frac{8\pi}{3}\hat{\varrho}_{b}(\bm{r})+\frac{\Lambda}{3}\
.$ (47)

Here, $\hat{\varrho}_{b}$ is an auxiliary function interpreted as the local
background density; then, the physical density can be split as follows:

$4\pi\varrho_{\textrm{\tiny{ini}}}=4\pi\hat{\varrho}_{b}+4\pi\hat{\varrho}_{b}\hat{\delta}_{\textrm{\tiny{ini}}}=4\pi\hat{\varrho}_{b}-\hat{W}\
.$ (48)

The initial scalar curvature takes the following form:

${\cal{R}}_{\textrm{\tiny{ini}}}=6\hat{k}-4\left(\hat{W}+\hat{H}_{\textrm{\tiny{ini}}}\hat{U}\right)-{\rm
II}\left(\hat{U}^{a}_{\>\>i}\right)\ ,$ (49)

which obeys the exact energy constraint on the initial time slice. The
quantity ${\rm II}\left(\hat{U}^{a}_{\>\>i}\right)$ represents the second
principal scalar invariant of the generalized Newtonian peculiar-velocity
gradient field. See its definition in Eq. (112).

The initial gravitoelectric and -magnetic parts of the (spatially projected)
Weyl tensor are given by

$\displaystyle E^{i}_{\>\>j}(t_{\textrm{\tiny{ini}}})$
$\displaystyle=-\hat{W}^{i}_{\>\>j}+\frac{\hat{W}}{3}\delta^{i}_{\>\>j}\ ;$
(50) $\displaystyle H^{i}_{\>\>j}(t_{\textrm{\tiny{ini}}})$
$\displaystyle=\\!-\epsilon^{ikl}_{\>\>}\left(\dot{A}_{\textrm{\tiny{ini}}|l}G_{jk}+G_{ka}\hat{U}^{a}_{\>\>j||l}\right)\
.$ (51)

It should be noticed that the splitting of the initial density distribution,
$\varrho_{\textrm{\tiny{ini}}}$, into the associated background density and
its contrast is thus far arbitrary in our model, setting an extra degree of
freedom in the choice of $\hat{\varrho}_{b}$. This property is inherited from
the Szekeres solution, where $\mu(z)$ (i.e., $\hat{\varrho}_{b}$) is an
arbitrary free function.

#### Setting the local background density $\hat{\varrho}_{b}$

We choose a physically motivated form of the local background density, namely
a constant function with value $\varrho_{b}(t_{\textrm{\tiny{ini}}})$, the
mass density of an FLRW solution at $t_{\textrm{\tiny{ini}}}$:

$\hat{\varrho}_{b}=\varrho_{b}(t_{\textrm{\tiny{ini}}})=\mbox{const.}$ (52)

This condition removes the indeterminacy in the local background density and
uniquely defines the set of initial conditions. The physical motivation behind
this choice is the almost-homogeneity that is expected at early cosmic stages.

Note that Eq. (52) does not imply that $\mu$ is constant in the exact Szekeres
limit, which would suppress the growth of structure. In fact,
$\mu=(4\pi\big{/}3)\hat{\varrho}_{b}(t_{\textrm{\tiny{ini}}})\chi^{-3}\neq\mbox{const.}$
since $\chi(\bm{r})\neq\mbox{const.}$444 In the exact solution of class I,
${\cal{S}}_{,z}=0$ necessarily renders the model an FLRW spacetime.

##  V Special subcases: locally one-dimensional deformations

So far, we have kept the approach as general as possible. However, we can
still introduce some additional assumptions to obtain simpler and more
transparent physical models. Let us consider the case of a locally one-
dimensional deformation field, namely

$\hat{P}^{a}_{\>\>i}(t,\bm{r})=\hat{P}(t,\bm{r})\delta^{a}_{\>\>z}\delta^{z}_{\>\>i}\
.$ (53)

Hence, the linearized evolution equations turn out to be exact and the
functional expressions introduced in Section IV.1 become linear,

$\displaystyle\Theta$ $\displaystyle=$ $\displaystyle
3\hat{H}+\frac{\dot{\hat{P}}}{\mathfrak{J}}\ ;$ (54a)
$\displaystyle\sigma^{i}_{\>\>j}$ $\displaystyle=$
$\displaystyle\frac{1}{\mathfrak{J}}\left(\dot{\hat{P}}^{i}_{\>\>j}-\frac{\dot{\hat{P}}}{3}\delta^{i}_{\>\>j}\right)=\Sigma\times\hbox{Diag}\left(1,1,-2\right)\
;$ (54c)
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\hbox{with}\quad\Sigma:=-\frac{1}{3}\frac{\dot{\hat{\hat{P}}}}{\mathfrak{J}}\
;$ $\displaystyle{\cal{R}}$ $\displaystyle=$ $\displaystyle
6\frac{\hat{k}}{A^{2}}-\frac{4}{\mathfrak{J}}\left(3\hat{H}\dot{\hat{P}}+\ddot{\hat{P}}\right)\
;$ (54d) $\displaystyle E^{i}_{\>\>j}$ $\displaystyle=$
$\displaystyle-\frac{1}{\mathfrak{J}}\left(\ddot{\hat{P}}^{i}_{\>\>j}+2\hat{H}\dot{P}^{i}_{\>\>j}\right)+\frac{1}{3\mathfrak{J}}\left(\ddot{\hat{P}}+2\hat{H}\dot{\hat{P}}\right)\delta^{i}_{\>\>j}.\quad$
(54e)

The gravitomagnetic part of the Weyl tensor can be obtained by inserting (53)
into (36).

The presence of shear, Ricci, and spatially projected gravitoelectric Weyl
tensors with two equal eigenvalues constitutes a characteristic feature of
these solutions. The importance of the present toy model is two-fold: first,
it is the simplest realization of GRZA and its particular case, RZA, and
second, it contains the Szekeres models and the subfamily of exact solutions
resulting from them.

###  V.1 Exact subcases and RZA limit

The architecture of GRZA allows us to reproduce the class I of Szekeres
solutions when the Gram’s matrix and the deformation field take the form
specified in Eq. (25), furnishing a subset of the admissible initial data for
the locally one-dimensional models. Under these assumptions, the constraints
are identically satisfied, while Eq. (27) and (43) are the only nontrivial
dynamical equations. Equation (27) is nothing more than Eq. (21), or (14); and
(43) reduces to (23), or (16).

The well-known LTB models emerge as the spherically symmetric limit of the
quasispherical Szekeres models of class I. This occurs when the unspecified
Szekeres functions in (25) ($\nu$, ${\cal{A}}$ and ${\cal{W}}$) are set as in
(70) with $\epsilon=1$ and $c_{0},c_{1},c_{2}=\hbox{const.}$ 555The whole
class of LTB models is obtained from the Szekeres solution when the dipole
vanishes, i.e., when $\mathcal{X}=\mathcal{Y}=\mathcal{Z}=0$ in Eq. (78) or
equivalently, $S$, $P$ and $Q$ are constant in Eq. (79).

RZA results from GRZA in the limit when the curvature and density parameters
of the reference model are constant:

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

Then, the conformal scale factor is a function of time only,
$A(t,\bm{r})\rightarrow a(t)$. In this limit, Szekeres models of class II
follow from a simple specialization of the initial data, Eq. (24) and (71).
Figure 1 provides a visual representation of the relationship between these
solutions; see Fig. 2.1 and 2.4 in kras1  for a diagrammatic representation of
the exact solutions that can be deduced as a limit of the Szekeres spacetimes.

Figure 1: Special subcases of GRZA. Other exact limits of Szekeres class I
(Sze-I) and Szekeres class II (Sze-II) are not represented in the figure; see
kras1  for their classification.

##  VI A numerical example

To illustrate the potential of GRZA, we explore the numerical solutions of a
family of locally one-dimensional models, containing the quasispherical
Szekeres solution as a particular case.666A jupyter notebook containing a
pedagogical introduction to the numerical example shown in this section is
available at
https://github.com/idgaspar/Beyond_relativistic_Lagrangian_perturbation_theory.
A key aspect of these models is that they retain the spatial foliation of the
hypersurfaces by nonconcentric 2-spheres, a characteristic feature of the
quasispherical Szekeres models, but generalize Szekeres in the way the
2-spheres are arranged. Such a nontrivial foliation leads to more general
networks of structures than the standard Szekeres dipolar mass distribution.

####  VI.0.1 Szekeres-like initial data

The initial data is specified at the last scattering time. First, we set the
free parameters of the set of local Friedmannian reference models. For
simplicity, we take for the density field a single global value (while we keep
the dependence on the radial coordinate of the curvature parameter) of a
conventional FLRW universe at the initial time
($\varrho_{b}(t_{\textrm{\tiny{ini}}}):=\varrho_{bi}$):

$\hat{\varrho}_{b}(\bm{r})=\varrho_{bi}\ .$ (56)

To work with dimensionless quantities, we take the Hubble parameter at the
initial time ($H_{\ast}=H(t_{ls})$) as the timescale and $\ell=2\times
10^{-2}$ Mpc as the characteristic length of the inhomogeneous region. Next,
we introduce the following parameters and functions,777With a slight abuse of
notation, in this section we let $\mu$ denote the dimensionless density.

$\tau=H_{\ast}t\ ,\quad\xi=\frac{z}{\ell}\
,\quad\mu=\frac{4\pi}{3}\frac{\varrho}{H_{\ast}^{2}}\
,\quad\varkappa=\frac{\hat{k}}{H_{\ast}^{2}}\ .$ (57)

In these variables, $\mu_{bi}$ denotes the density value of the equivalent
FLRW model at the initial time, and the remaining arbitrary free function of
the Friedmann-like reference model is given by

$\displaystyle\hat{\varkappa}(\xi)$ $\displaystyle=$
$\displaystyle\hat{\varkappa}_{a}(\xi)\xi^{2}\ ,$ (58)
$\displaystyle\hat{\varkappa}_{a}(\xi)$ $\displaystyle=$
$\displaystyle\hat{\varkappa}_{1}+\hat{\varkappa}_{2}\left(1+\tanh
4\left(\xi-1/2\right)\right)\ ,$ (59)

with $\hat{\varkappa}_{1}=10^{-5}$ and $\hat{\varkappa}_{2}=37\times 10^{-5}$.
For our analysis, we consider these functions in the interval $0\leq\xi\leq
1$, $-\infty<x,y<+\infty$. In this model, $\xi$ should be interpreted as a
radial comoving coordinate, while $x$ and $y$ are the angular coordinates; see
Section VI.0.2 and Appendix E below.

The Gram’s matrix (initial metric) is taken as

$G_{ij}=\mbox{Diag}\left[\left({\cal{A}}-\hat{{\cal{F}}}_{\textrm{\tiny{ini}}}\right)^{2}\left(\frac{{\cal{W}}}{\chi}\right)^{2},\xi^{2}\mathcal{E}^{-2},\xi^{2}\mathcal{E}^{-2}\right]\
,$ (60)

where $\chi=\sqrt{\hat{\varkappa}}$ and the metric functions $\mathcal{E}$,
${\cal{A}}$ and ${\cal{W}}$ are chosen as in Szekeres models, Eq. (84a)-(84b),
with $S(\xi)=1$, $Q(\xi)=0$ and $P(\xi)$ is a third order polynomial function
satisfying $P(0)=P^{\prime}(0)=P^{\prime}(\xi_{1})=0$ and $P(\xi_{1})=-0.6$.
Moreover, $\hat{{\cal{F}}}_{\textrm{\tiny{ini}}}$ is the initial value of the
function

$\hat{{\cal{F}}}={\cal{F}}+\delta{\cal{F}}\
,\quad\delta{\cal{F}}=\delta\beta_{+}f_{+}+\delta\beta_{-}f_{-}\ .$ (61)

Here, ${\cal{F}}=\beta_{+}f_{+}+\beta_{-}f_{-}$ is the Szekeres function
introduced in Section III, and $\delta\beta_{\pm}(\bm{\xi})$ are arbitrary
functions of the spatial coordinates, $\bm{\xi}=\left(\xi,x,y\right)$. In this
simulation, we take

$\delta\beta_{+}(\bm{\xi})=\alpha\sin^{3}(\gamma\pi y)\sin^{4}(\pi\xi)\ ,$
(62)

with $\gamma=3$, and vary $\alpha$. Also, in order to have a model where the
voids compensate for the mass excess in the overdensities (i.e., a mass-
compensated array of structures), $\delta\beta_{-}(\bm{\xi})$ satisfies

$\delta\beta_{-}(\bm{\xi})=-\frac{\delta\beta_{+}f_{+}(t_{\textrm{\tiny{ini}}},\xi)}{f_{-}(t_{\textrm{\tiny{ini}}},\xi)}\
.$ (63)

Then, the initial energy density is set to

$\mu_{\textrm{\tiny{ini}}}=\mu_{bi}\left(\frac{{\cal{A}}}{{\cal{A}}-\hat{{\cal{F}}}_{\textrm{\tiny{ini}}}}\right)\
,$ (64)

which determines the generalized Newtonian peculiar-acceleration divergence,

$\hat{W}=-3\mu_{bi}\frac{\hat{{\cal{F}}}_{\textrm{\tiny{ini}}}}{{\cal{A}}-\hat{{\cal{F}}}_{\textrm{\tiny{ini}}}}\
.$ (65)

The whole set is completed by the initial values of the deformation,
$P(t_{\textrm{\tiny{ini}}})=0$, and the generalized Newtonian peculiar-
velocity divergence,

$\hat{U}=-\frac{\dot{\hat{{\cal{F}}}}_{\textrm{\tiny{ini}}}}{{\cal{A}}-\hat{{\cal{F}}}_{\textrm{\tiny{ini}}}}\
.$ (66)

####  VI.0.2 Physical interpretation

This family of models corresponds to a solution with a deformation field of
the form (53) with

$\hat{P}(t,\bm{\xi})=-\frac{\hat{{\cal{F}}}-\hat{{\cal{F}}}_{\textrm{\tiny{ini}}}}{{\cal{A}}-\hat{{\cal{F}}}_{\textrm{\tiny{ini}}}}\
,$ (67)

containing the exact Szekeres solution in the case
$\delta\beta_{+}(\bm{\xi})=0$. As the amplitude of $\delta\beta_{+}(\bm{\xi})$
increases, the model deviates from Szekeres. At the same time, its coupling to
the growing mode leads to significantly different evolutions, even for small
values of $\delta\beta_{+}(\bm{\xi})$, compared to the exact solutions. The
function $\delta\beta_{+}(\bm{\xi})$, its derivatives and the Szekeres dipole
functions vanish at the boundary of the domain considered in this simulation,
$\xi=\xi_{1}$, indicating that the solution is an exact LTB model at this
shell. This property allows matching the inhomogeneous region simulated here
to an external Szekeres, LTB, FLRW, or Schwarzschild solution.

Writing out the line-element (28), one can check that the metric of the
surfaces $\left\\{t=\hbox{const.},\xi=\hbox{const.}\right\\}$ coincides with
its equivalent in the exact solution, and then it can be transformed into that
of a 2-sphere with radius $A(t,\xi)\chi(\xi)$ (by the transformation (76)).
The deviation from the exact solution comes from the metric coefficient
$g_{\xi\xi}$, which is generalized to resemble the usual approximation where
the growing (decaying) mode is coupled to a nontrivial function of the spatial
coordinates. In this way, the models retain the characteristic foliation by
nonconcentric 2-spheres of Szekeres; the term $\delta\beta_{+}f_{+}$ leads to
more general arrangements of the 2-spheres. This feature will be examined in
detail in Section VI.1.

####  VI.0.3 Estimation of the accuracy

The GRZA model, although containing exact solutions as special subcases, is in
general an approximation. The Hamiltonian constraint, or more specifically its
violation, can be used to estimate the error generated by this approach.
Taking into account that the fluid flow $4-$velocity is the vector normal to
our spatial hypersurfaces, we define

$\mathbb{H}:=\left(G_{\mu\nu}-8\pi T_{\mu\nu}\right)u^{\mu}u^{\nu}=G_{00}-8\pi
T_{00}\ .$ (68)

Then, the function $\mathbb{H}$ vanishes when the solution is exact, and it is
nontrivial for approximate models. In general, we can tolerate small
violations of the constraints as long as they remain below a certain threshold
value. Nonnegligible values of $\mathbb{H}(\bm{\xi},t)$ can be attributed to
either an increase of the numerical error or to the natural failure of the
approximate approach as the solution goes deep into the nonlinear regime.

We exploit the relationship between the present models and the exact solutions
to work out a closed-form expression for $\mathbb{H}$. The final result is Eq.
(117), displayed in Appendix E (where we rewrite this section in the language
of exact solutions). Note that the error does not grow linearly in time but
exhibits a nontrivial dependence on the ‘perturbation’, reinforcing the
nonlinear nature of the current approach.

We note in this context that a self-consistency test, proposed by Doroshkevich
et al. for the Newtonian Zel’dovich approximation in Doroshkevichetal73 , see
also Buchert1989AA , would return a vanishing error, since the comparison of
the density field (i) as a solution of the field equation and (ii) as a
solution of the continuity equation, leads to an exact agreement in the case
considered here: this error depends on second and third principal scalar
invariants of the generalized deformation gradient that both vanish in the
class of models considered here.

(a) $\alpha=0$

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

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

(d) $\alpha=2.9\times 10^{-3}$

(e) $\alpha=3.8\times 10^{-3}$

(f) $\alpha=4.8\times 10^{-3}$

Figure 2: Equatorial projection of the density contrast
$\delta=\varrho/\varrho_{b}-1$ at the present cosmic time, $t_{0}$.
Coordinates are defined as $X=A(t_{0},z)z\cos\phi$ and
$Y=A(t_{0},z)z\sin\phi$, where $A(t_{0},z)z$ is the areal radius. The plots
correspond to different values of $\alpha$. For $\alpha=0$, the solution is
exact, exhibiting the distinctive dipolar matter distribution of Szekeres
models.  Figure 3:  Estimation of the error from the violation of the
Hamiltonian constraint at the present cosmic time. Equation (68) is evaluated
along the radial rays defined by $\theta=\pi/2$ and
$\phi=\left[\phi_{\ast}^{\tiny{\hbox{(1)}}},\phi_{\ast}^{\tiny{\hbox{(2)}}},0,\pi/4,\pi/2,\pi,3\pi/2\right]$
for the solution with $\alpha=4.8\times 10^{-3}$.
$\phi_{\ast}^{\tiny{\hbox{(1)}}}$ and $\phi_{\ast}^{\tiny{\hbox{(2)}}}$ denote
the azimuthal angles of the two largest overdensities in the order seen by the
reader in Fig. 2(f).

####  VI.0.4 Numerical results

The GRZA model equations are solved in the comoving domain specified above,
${\cal D}_{i}:0\leq\xi\leq 1$. The inhomogeneous solution is contrasted by a
theoretical FLRW model with the standard cosmological constant and energy
density values. Our choice of the initial data ensures that the mass in the
initial domain is the same as it would be in an FLRW model
($\varrho_{b}(t_{\textrm{\tiny{ini}}})={\left\langle\varrho_{\textrm{\tiny{ini}}}\right\rangle_{{\cal
D}_{H\bf i}}}$) – this mass is subsequently conserved through the continuity
equation. Furthermore, the initial curvature parameter of this FLRW model is
taken as the average initial curvature,
$k_{b}={\left\langle{\cal{R}}_{\textrm{\tiny{ini}}}\right\rangle_{{\cal
D}_{H\bf i}}}/6$. In more realistic cosmological simulations, such a
background will emerge from the averaging on a sufficiently large compact
domain. However, the scales considered here correspond to the size of a
typical supercluster, and the background is set by hand. For merely
pedagogical reasons, we have chosen to set $\hat{k}$ as a strictly positive
function and assumed a background with positive curvature as well.888 In the
GW parametrization, the conformal scale factor has a vertical asymptote at the
points where $\hat{k}(z)=0$. Although the rescaling of the solution removes
this singular behavior, we took a well-behaved model to illustrate how the
initial data of Szekeres can be set in Hellaby’s formulation, then
reparametrized in the language of GW, rescaled, rewritten in the language of
GRZA, and finally generalized to obtain the approximate solutions discussed in
this section. This is the approach followed in Appendix E.

(a) LTB

(b) Szekeres

(c) GRZA: Numerical example

Figure 4: Foliation by 2-spheres plotted in terms of the radial proper length.
Time slices of LTB and Szekeres solutions and of our numerical example model
are leaves of a foliation by topological 2-spheres. Only in the LTB case, the
2-spheres are geometric and concentric. Overdense structures (or density
maxima) arise at the loci of points where the proper distance between shells
is minimal and becomes zero at shell-crossing. In the Szekeres and more
general GRZA cases, these degeneracies also arise along the angular variables
but have well-localized spatial positions at points of degeneracy of the
deformation field, $\det(\eta^{a}_{\ i})=0$. The regions marked in the figure
are prone to forming shell-crossings.

Figure 2 shows snapshots of the density contrast at the present cosmic time
for the different values of $\alpha$, including $\alpha=0$, the exact Szekeres
case. In all of these simulations, the error remains small. Figure 3 shows
$\mathbb{H}$ as a function of $z$ for different values of the polar angle for
the case $\alpha=4.5\times 10^{-3}$. While for the other cases the error
remains smaller, being negligible (purely numerical) in the exact case
($\alpha=0$).

###  VI.1 Overdense structures and shell-crossings

In the presented models, as in other dust-sourced models, structures emerge
from the inhomogeneous expansion and shear of the cosmic fluid. Portions of
the fluid expanding faster than the global average distribution with
subdominant shear become voids (underdensities), whereas a collapse, i.e. an
expansion slower than the global average distribution or a dominant shear
leads to clustering (overdensities) BolejkoInhCosModels . In the Szekeres
solutions their complex spatial matter distribution can be traced back to
overdense spatial patterns in angular directions of a foliation of the
constant time slices by nonconcentric 2-spheres, see Figure 4.

To gain a better insight into this property let us consider first the LTB
models, where the 2-spheres are concentric and structures correspond to
spherical dust shells. Szekeres models can be thought of as a nonspherical
generalization of LTB via the dipole function (see Appendix B.2), where the
(topological) 2-spheres are no longer orthogonal to the radial rays when they
are plotted in terms of proper length sussman2015multiple . Similarly, the
models developed in this section constitute a generalization of Szekeres
through the function $\delta\beta_{+}$, furnishing a more general foliation
into topological 2-spheres (with “less quasisymmetries”), and in turn, leading
to a more generic network of structures. Figure 4 shows the spatial foliation
by 2-spheres in terms of the proper length; we also highlight the regions
where the proper distance between adjacent shells is minimal, giving rise to
local density maxima, i.e., overdensities. Eventually, this distance will go
to zero resulting in shell-crossing singularities kras2 .

It is important to emphasize that the models examined here (and displayed in
Fig. 2) inevitably lead to shell-crossings, simply because any inhomogeneous
physical source that would allow to counteract gravity (such as velocity
dispersion, pressure or vorticity) is absent in the basic system of equations
(3) studied. Although these shell-crossing singularities will emerge in the
future, i.e., $t>t_{0}$, their spatial locations coincide with regions where
the three overdensities in our example model are located. Once singularities
in the density field emerge, the simulation would ask for regularization of
these singularities. We remark here that, in any realization of dust models
and their shell-crossings at finite resolution, these singularities are
smoothed out, so that the distribution of caustics in the density field is
indicative for the network of high-density structures, see the methods of
realization and illustrations at high resolution in Newtonian models
Buchert1989AA ; buchertbartelmann ; highresolution , and for further
discussions on shell-crossing singularities in the literature and their
morphological classification in generic Lagrangian-to-Eulerian maps, see
[sect.VB]rza1 , as well as the recent review LagTheorUniver  for related
discussions and a list of recent work on shell-crossing singularities in
Newtonian models.999Most of the literature on shell-crossing singularities
draws conclusions from Newtonian models. While indicative for shell-crossings
in the general-relativistic context, care has to be taken as to whether these
studies are applicable. Singularities in the density field have deeper
implications in general relativity, since these will in general come with
singularities in the geometry, see the discussion in [sect.VB]rza1 . Questions
of extendability of solutions beyond caustics are addressed in e.g. joshi ;
nolan .

In the mathematical literature, the singularity morphologies are classified
topologically. The geometry of caustics emanates from local germs of
classified stable caustic patterns, while their global geometry depends on the
large-scale distribution of the density field. Locally, the degeneracies of
directional scale factors can be used to locate these singularities in the
density field (corresponding to degeneracies of the determinant of the coframe
deformation field, $\det(\eta^{a}_{\ i})=0$, cf. Eq. (5)). They can be
identified as the eigendirections of the expansion tensor GW1 :
$\dot{\l}_{i}/\l_{i}=\Theta_{ii}$ (no summation implied). Then, for the models
examined in this section we simply have:

$\l_{3}=A(1+P)\ ,\qquad\l_{1}=\l_{2}=A\ .$ (69)

In general, shell-crossings arise due to the vanishing of the term $1+P$
(i.e., $\l_{3}\to 0$) while $\l_{1}=\l_{2}>\l_{3}>0$; resulting in what in
cosmology is called a ‘pancake singularity’, i.e. a 2-surface of infinite
density in the case of a 3-space.

##  VII General discussion and concluding remarks

This section brings together the steps and procedure followed to develop GRZA
and discusses their physical interpretation. Our analysis is put in
perspective with well-established features of exact solutions and the
Lagrangian perturbation theory. Finally, we summarize the main points of this
article and provide final remarks.

###  VII.1 Physical motivation

The relativistic Lagrangian perturbation schemes and exact solutions of
Einstein’s equations provide complementary approaches to analytically model
the large-scale structure formation in the Universe. On the one hand, the
Lagrangian description is known for its high accuracy; its first-order
solution (known as the Zel’dovich approximation in Newtonian theory for
initial data where peculiar-velocity and peculiar-acceleration are assumed to
be parallel, so-called ‘slaving’) remains valid beyond the linear regime
LagTheorUniver . On the other hand, exact solutions are a fundamental tool to
understand inhomogeneous cosmology, describing the structure formation in
idealized scenarios with symmetries or quasisymmetries but otherwise not tied
to any particular regime kras1 ; kras2 .

Motivated by our recent results on the correspondence between the relativistic
Lagrangian theory of structure formation and the Szekeres family of exact
solutions (the first-order Lagrangian perturbation solutions (RZA) contain the
entire class II of Szekeres solutions) rza6 , in this paper, we proposed a
generalization of the Lagrangian theory, controlled by Szekeres models of
class I. This furnishes a new method for structure formation in relativistic
cosmology: the generalized Relativistic Zel’dovich Approximation (GRZA).

The GRZA architecture relies on a combination of standard elements from the
Lagrangian theory and exact solutions. Firstly, we retain the essence of RZA
and consider the spatial coframe coefficients as the only variables of the
$3+1$ Lagrange-Einstein system (3) for irrotational dust in the comoving frame
rza1 . Next, the coframes are split into the expansion of a trivial background
model and a deformation field. In this splitting, we introduce a space-
dependent scale factor, which obeys a Friedmann-like evolution equation with
in general space-dependent density and curvature parameters, Eq. (26)-(27).
Still, the overall model consists of a deformation field evolving on this
inhomogeneous Friedmann-like reference model. As in RZA, the evolution
equations are linearized for the deformation, Eq. (37) (or (40)), while the
functional evaluation of the physical quantities is kept exact, Section IV.1.

The inclusion of such an inhomogeneous Friedmann-like reference model is
motivated by the structure of Szekeres class I solutions. However, it can be
related to a characteristic feature of the silent solutions of Einstein’s
equations, where some Bianchi and local rotational symmetric solutions, LTB
and Szekeres models are governed by a common set of evolution equations. Their
discrimination is reflected in the constraint equations
Bolejko:2018SilentUniv1 . In this way, our space-dependent scale factor allows
for fulfilling the constraints associated with more general solutions,
including the Szekeres class I exact case.

To put it differently, GRZA is obtained as a generalization of RZA, where the
constant curvature and density of the Friedmann background are extended to be
space-dependent functions (leading to a space-dependent scale factor):
$k_{0}\rightarrow\hat{k}(\bm{r})$ and
$\varrho_{b}(t_{i})\rightarrow\hat{\varrho}(\bm{r})$. This procedure
constitutes the most direct extension of RZA that allows reproducing the
entire Szekeres class as an exact limit. When the solution is not exact we can
control its overall error by estimating the violation of the constraints.

###  VII.2 Backreaction in cosmological models

Let us now turn to the backreaction problem in cosmology, i.e, how small-scale
structures can back-react and change the Universe’s average expansion rate on
large scales. We here refer to the scalar volume-averaging scheme
Buchert2000GERG , reviewed in buchert:review ; BuchertFocus ; buchertrasanen .
As inherited from Szekeres models, GRZA solutions present a nontrivial global
backreaction term, which can be computed by inserting (8), (32) and (34) into
(111). However, proper physical models should be compatible with the current
observational probes, and the backreaction calculated from classes of exact
solutions might be sensitive to these constraints. This is the case of the LTB
void models used to fit the observations without dark energy, which have been
ruled out by the constraints imposed by the kinematic Sunyaev Zel’dovich
effect SussmanVoidsNoBack ; Zibin:2011 ; Zibin:2011ma ; Bull:2011wi .

The situation will be different when considering a generic model: with GRZA we
keep control of approximations made since its construction contains the exact
subcases while going beyond the restricted space of initial data. Foremost we
emphasize the role of higher scalar invariants contained in GRZA but required
to strictly vanish in the exact subcases. Especially for a discussion of
cosmological backreaction, a nonvanishing second invariant implies the
possibility of a sign change in the backreaction functional: a positive
backreaction contributes to acceleration. This generic feature of backreaction
is otherwise suppressed, but it is a consequence of the volume-dominance of
void regions. Expressed in geometric terms, given a large generic simulation
that includes large-volume voids, the average scalar curvature obtains
negative contributions mimicking a dark energy-like effect. These remarks show
that the implementation of GRZA for generic initial data is key to the future
development of the proposed simulations. They will help to quantify
cosmological backreaction globally.

In this context, we think of two questions to be answered. First, does a
generic GRZA simulation change the conclusions drawn from LTB models when
compared with observational constraints, i.e. what are the quantitative
changes by including more degrees of freedom? Second, how backreaction or the
emergence of an average negative curvature can solve the tension in the
estimation of $H_{0}$, given the controlled quantification proposed in this
paper. How does this compare with silent universe simulations
Bolejko:2018SilentUniv1 ; Bolejko:2018SilentUniv2 ;
BolejkoOstrowski2019SilentUniv3 ?

These and further questions will be addressed in a forthcoming paper.

###  VII.3 Source of the gravitational field

The present study has been restricted to models sourced by an irrotational
dust fluid, dictated by the basic system of equations (3), and aiming at
applications in the theory of large-scale structure formation in relativistic
cosmology. Both, RZA and GRZA schemes can be extended to include more general
sources (isotropic and anisotropic) pressure, velocity dispersion and heat-
flow terms, and further sources in the energy-momentum tensor. In these more
general cases of a single-fluid model, other nontrivial components should be
included in the metric (28) and a more general four-velocity should be assumed
to assure the compatibility of the overall system of equations for the
deformation field, e.g. if a nonzero vorticity of the fluid is to be included.

In concrete terms, the basic system (3) can be generalized by including lapse
and shift degrees of freedom, as well as a four-velocity that is tilted with
respect to the normals to the spatial hypersurface foliation, or if
hypersurfaces do not exist, a threading of spacetime can be implemented by
using threading lapse and shift functions, see e.g. BMR  for the general
foliation and threading formalisms. Writing such a general foliated or
threaded geometry of spacetime in terms of coframe fields will determine the
added geometrical degrees of freedom through the energy-momentum tensor via
Einstein’s equations, e.g. a pressure gradient will require a nonconstant
lapse function, a rotational flow or a nondiagonal energy-momentum tensor will
require shift functions and a tilted four-velocity. For the RZA model an
extension of the general-relativistic Lagrangian perturbation theory to
irrotational perfect fluid models with a barotropic equation of state has been
investigated rza5  using an inhomogeneous lapse function buchert:fluid .

###  VII.4 Summary and concluding remarks

We have presented a new nonlinear method to model large-scale structure
formation in the Universe (GRZA). In its development, we merged elements of
the Lagrangian perturbation theory and exact solutions; consequently, GRZA
contains Szekeres (then LTB) models and RZA as particular limits. This new
approach retains the essence of the Lagrangian perturbation theory: the
dynamics is described in terms of a deformation field evolving on a
cosmological background. But, it generalizes the global FLRW background to an
inhomogeneous Friedmann-like reference model.

The evolution of the GRZA models is formally governed by the same system of
equations as the one defining RZA; in this sense, their difference is
reflected by the constraints. This similarity in the evolution allows to
translate the approach into the Newtonian language, which can be exploited to
examine relativistic corrections to the current N-body numerical simulations.
Through their nonperturbative nature, GRZA models are well-suited to explore
the impact of the inhomogeneities on the prediction of $H_{0}$ or the
Universe’ large-scale average evolution (i.e., to address the Hubble tension
and the backreaction problem), to mention a few examples.

To illustrate the method’s potential, we provided a family of numerical
examples developed as a generalization of the quasispherical Szekeres models.
To a certain extent, the path followed in this examination remains pedagogical
and explorative. More general sets of initial data will open the door to
realistic simulations, complementing the current relativistic numerical
approaches. We highlight the simplified initial conditions employed to run
simulations based on the silent universes approach bolejko2007evolution ;
Bolejko:2018SilentUniv1 ; Bolejko:2018SilentUniv2 ;
BolejkoOstrowski2019SilentUniv3  as a ‘must-explore’ set.

The fact that GRZA contains the entire family of Szekeres (and thus LTB)
models is a remarkable theoretical feature. However, this does not necessarily
justify its use over RZA, since GRZA has a more complex mathematical
structure. Physical applications will ultimately reveal the quality of the
approach and whether or not its use is justified.

Acknowledgments: This work is part of a project that has received funding from
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (grant agreement ERC advanced grant
740021-ARTHUS, PI: TB). IDG and JJO acknowledge the support of the National
Science Centre (NCN, Poland) under the Sonata-15 research grant
UMO-2019/35/D/ST9/00342. We would like to thank Roberto Sussman for valuable
remarks and discussion, and the anonymous referee for valuable suggestions to
improve the paper. IDG and JJO also acknowledge hospitality during their
working visits at CRAL-ENS, Lyon.

##  Appendix A The Goode and Wainwright parametrization of Szekeres models

We start by giving the relevant functions for class I ($\beta_{,z}\neq 0$):

$\displaystyle{\cal{S}}={\cal{S}}(t,z)\
,\quad\hbox{with}\quad{\cal{S}}_{,z}\neq 0\ ,$ (70a) $\displaystyle
f_{\pm}=f_{\pm}(t,z)\ ,\quad t_{B}=t_{B}(z)\ ,\quad\mu=\mu(z)\ ,$ (70b) and
$\displaystyle e^{\nu}=f(z)\mathcal{E}^{-1}\ ;$ (70c)
$\displaystyle\mathcal{E}=\Big{[}c_{0}(z)(x^{2}+y^{2})+2c_{1}(z)x+2c_{2}(z)y+c_{3}(z)\Big{]}^{-1};$
(70d) $\displaystyle c_{0}c_{3}-c_{1}^{2}-c_{2}^{2}=\epsilon/4\
,\,\quad\quad\epsilon=0\ ,\pm 1\ ;$ (70e)
$\displaystyle{\cal{A}}=f\nu_{,z}-k_{0}\beta_{+}\
,\quad\;{\cal{W}}^{2}=(\epsilon-k_{0}f^{2})^{-1}\ ;\;$ (70f)
$\displaystyle\beta_{+}=-k_{0}f\mu_{,z}/(3\mu)\ ,\,\quad\beta_{-}=f\,t_{B,z}\
.$ (70g)

The arbitrary function $t_{B}(z)$ is the inhomogeneous ‘big bang time’.

Correspondingly, for class II ($\beta_{,z}=0$):

$\displaystyle{\cal{S}}={\cal{S}}(t)\ ,\quad f_{\pm}=f_{\pm}(t)\ ,\quad
t_{B},\ \mu=\hbox{const.}\ ;$ (71a) $\displaystyle k_{0}=0,\pm 1\
,\quad{\cal{W}}=1\ ;$ (71b) $\displaystyle
e^{\nu}=\left[1+\frac{k_{0}}{4}(x^{2}+y^{2})\right]^{-1}\ ,$ (71c) and
${\cal{A}}=\left\\{\begin{array}[]{l
l}e^{\nu}\Big{[}c_{0}(z)\left(1-\frac{k_{0}}{4}(x^{2}+y^{2})\right)+c_{1}(z)x\\\
\hskip 25.6073pt+c_{2}(z)y\Big{]}-k_{0}\beta_{+}\ ,\;\hbox{for}\;k_{0}=\pm 1\
;\\\ c_{0}(z)+c_{1}(z)x+c_{2}(z)y\\\ \hskip
25.6073pt-\beta_{+}(z)(x^{2}+y^{2})/2\ ,\;\hbox{for}\;k_{0}=0\
;\end{array}\right.$ (71d)

$c_{i}(z)$ and $\beta_{\pm}(z)$ arbitrary and $\mu$, $t_{B}=\text{const.}$

##  Appendix B A deeper insight into Szekeres models of class I

Szekeres models of class I are traditionally examined in an LTB-like
parametrization introduced by Hellaby hellaby1996null . In this
representation, the line-element is given by

${\hbox{d}}s^{2}=-{\hbox{d}}t^{2}+\frac{\left(\Phi_{,z}-\Phi\mathcal{E}_{,z}/\mathcal{E}\right)^{2}}{\epsilon-k(z)}{\hbox{d}}z^{2}+\frac{\Phi^{2}}{\mathcal{E}^{2}}\left({\hbox{d}}x^{2}+{\hbox{d}}y^{2}\right)\
,$ (72)

with

$\mathcal{E}=\frac{S(z)}{2}\left(\left(\frac{x-P(z)}{S(z)}\right)^{2}+\left(\frac{y-Q(z)}{S(z)}\right)^{2}+\epsilon\right)\
,$ (73)

where the latter expression is just a reparametrization of (70d). Here,
$\Phi=\Phi(t,z)$, $\epsilon=\pm 1,0$ (for quasispherical, quasihyperbolic and
quasiplane models, respectively), and $k(z)$, $P(z)$, $Q(z)$ and $S(z)$ are
arbitrary functions.

The function $\Phi$ obeys a Friedmann-like equation,

$\dot{\Phi}^{2}=-k(z)+\frac{2M(z)}{\Phi}+\frac{1}{3}\Lambda\Phi^{2}\ ,$ (74a)
which can be integrated to give
$t-t_{B}(z)=\int_{0}^{\Phi}\frac{{\hbox{d}}\tilde{\Phi}}{\sqrt{-k+2M/\tilde{\Phi}+\frac{1}{3}\Lambda\tilde{\Phi}^{2}}}\
.$ (74b)

$M(z)$ and $t_{B}(z)$ are arbitrary integration functions, identified as the
effective gravitational mass and the time of the initial singularity (the ‘big
bang time’).

In these variables the density reads:

$4\pi\varrho=\frac{M_{,z}-3M\mathcal{E}_{,z}/\mathcal{E}}{\Phi^{2}\left(\Phi_{,z}-\Phi\mathcal{E}_{,z}/\mathcal{E}\right)}\
,$ (75)

while the expressions for the expansion, shear and curvature tensors can be
seen in Section 2.3.2 of Ref. BKHC2009 .

###  B.1 Coordinate system

The coordinate system of the standard representation of Szekeres models has
the physical interpretation of a stereographic projection of a sphere
($\epsilon=1$), plane ($\epsilon=0$), or hyperboloid ($\epsilon=-1$) kras1 ;
kras2 ; Bolejko2006Struformation . In the quasispherical case, the
transformation

$\displaystyle x$ $\displaystyle=$ $\displaystyle
P+S\cot\left(\vartheta/2\right)\cos(\phi)\ ;$ (76a) $\displaystyle y$
$\displaystyle=$ $\displaystyle Q+S\cot\left(\vartheta/2\right)\sin(\phi)\ ,$
(76b)

turns the metric of the surfaces
$\left\\{t=\mbox{const.},z=\mbox{const.}\right\\}$ into that of a 2-sphere,

$\frac{\Phi^{2}}{{\cal{E}}^{2}}\left({\hbox{d}}x^{2}+{\hbox{d}}y^{2}\right)=\Phi^{2}\left({\hbox{d}}\vartheta^{2}+\sin^{2}\vartheta{\hbox{d}}\phi^{2}\right)\
.$ (77)

###  B.2 Modeling of the cosmic network

It is often convenient to employ spherical coordinates to examine the dynamics
of the quasispherical models. Performing the stereographic transformation of
coordinates (76) and $z=r$, one finds that all the information about the
departure from spherical symmetry is contained in the ‘Szekeres dipole’
function,

$\mathfrak{W}:=\frac{r{\cal{E}}^{\prime}}{{\cal{E}}}=-\mathcal{X}\,\sin\vartheta\,\cos\phi-\mathcal{Y}\,\sin\vartheta\,\sin\phi-\mathcal{Z}\,\cos\vartheta\
,\quad$ (78)

where we have defined $\mathcal{X}$, $\mathcal{Y}$ and $\mathcal{Z}$ as

$\mathcal{X}=\frac{rP^{\prime}}{S};\qquad\mathcal{Y}=\frac{rQ^{\prime}}{S};\qquad\mathcal{Z}=\frac{rS^{\prime}}{S}\
,$ (79)

which in terms of the GW free functions take the form,

$\displaystyle\mathcal{X}=2z\left(\frac{c_{1}c_{0}^{\prime}}{c_{0}}-c_{1}^{\prime}\right)\
,\;\mathcal{Y}$ $\displaystyle=$ $\displaystyle
2z\left(\frac{c_{2}c_{0}^{\prime}}{c_{0}}-c_{2}^{\prime}\right)\ ,$
$\displaystyle\mathcal{Z}$ $\displaystyle=$
$\displaystyle-z\frac{c_{0}^{\prime}}{c_{0}}\ .$ (80)

For more details see Ref. SussBol2012 , where Sussman and Bolejko introduced a
description of this subclass employing weighed proper volume averages of the
local scalars (so-called q-scalars) and their fluctuations. The advantages of
this notation were exploited in Ref. sussman2015multiple  to obtain sufficient
conditions for the existence of extrema of the Szekeres scalars. In
particular, the functions $\mathcal{X}$, $\mathcal{Y}$ and $\mathcal{Z}$ (and
in turn $S$, $P$ and $Q$) set the angular position of the extrema in the curve
defined by $\mathfrak{W}_{,\vartheta}=\mathfrak{W}_{,\phi}=0$:

${\cal{B}}_{\pm}(r)=[r,\,\vartheta_{\pm}(r),\,\phi_{\pm}(r)],$ (81)

where $\vartheta_{\pm}$ and $\phi_{\pm}$ denote two antipodal positions for
each value of $r$,

$\displaystyle\phi_{-}$ $\displaystyle=$
$\displaystyle\arctan\left(\frac{\mathcal{Y}}{\mathcal{X}}\right),\;\;\vartheta_{-}=\arccos\left(\frac{\mathcal{Z}}{\sqrt{\mathcal{X}^{2}+\mathcal{Y}^{2}+\mathcal{Z}^{2}}}\right)\
,$ $\displaystyle\phi_{+}$ $\displaystyle=$
$\displaystyle\pi+\phi_{-},\qquad\vartheta_{+}=\pi-\vartheta_{-}\ .$ (82b)

The extrema of all the scalars lie along the curve ${\cal{B}}_{+}$. And we
have overdensities (/voids) in a given interval of $r$, if $M_{,r}>0$
(/$M_{,r}<0$). Section VI.0.4 considers $S=1$ and $Q=0$; thus, the extremum of
the exact case is located on the equatorial plane and along the x-axis.

The limit $\mathfrak{W}=0$ defines a unique spherically symmetric LTB model,
where the 2-spheres foliating the time slices become concentric (see Section
VI.1). This idea was exploited in sussman2015multiple  to introduce the term
term ‘LTB seed model’; then, every Szekeres model can be developed from an
‘LTB seed model’ by specifying the dipole parameters.

###  B.3 Relationship between the LTB-like and GW parametrizations

To obtain the GW parametrization, firstly, we introduce the constant
$k_{0}=|k(z)|/k(z)$, and define the metric function $f(z)$ as

$k_{0}f^{2}(z)=k(z)\quad\mbox{when}\quad k_{0}\neq 0\ .$ (83)

Otherwise, $f$ is arbitrary.

The remaining functions relate to those in the LTB-like representation as
follows IshakPeel2012LargeScale :101010When $k=0$, the scale factor is defined
by ${\cal{S}}=\Phi/\mu^{1/3}\implies\mu=1$ kras1 ; kras2 .

$\displaystyle c_{0}=1/(2S)\ ,\;c_{1}=-P/(2S)\ ,\;c_{2}=-Q/(2S)\ ;\qquad$
(84a) $\displaystyle c_{3}=(P^{2}+Q^{2}+\epsilon\,S^{2})/(2S)\ ;\qquad$ (84b)
$\displaystyle S=\Phi/f\ ,\;\;\mu=M/f^{3}\
;\;\;{\cal{W}}=\left(\epsilon-k\right)^{-\frac{1}{2}}\ ;\qquad$ (84c)
$\displaystyle{\cal{G}}=f\left(\frac{\Phi_{,z}}{\Phi}-\frac{{\cal{E}}_{,z}}{{\cal{E}}}\right)\
.\qquad\;$ (84d)

To complete the analysis, we recall some results of rza6 . Focusing on the
case $k\neq 0$, Eq. (84d) gives

$\displaystyle{\cal{G}}$ $\displaystyle=$ $\displaystyle
f\left(\frac{\Phi_{,z}}{\Phi}-\frac{{\cal{E}}_{,z}}{{\cal{E}}}\right)=\left(f\frac{{\cal{S}}_{,z}}{{\cal{S}}}+f_{,z}-f\frac{{\cal{E}}_{,z}}{{\cal{E}}}\right)$
(85) $\displaystyle=$ $\displaystyle
f\frac{{\cal{S}}_{,z}}{{\cal{S}}}+f\nu_{,z}\ ,$

with $e^{\nu}:=f/{\cal{E}}$. Under the parametrization (84), Eq. (74)
transforms into (14), which has the following formal integral (similar to
(74b)):

$t-t_{B}(z)=\int^{{\cal{S}}}_{0}\frac{{\hbox{d}}\hat{{\cal{S}}}}{\left(-k_{0}+\frac{2\mu}{\hat{{\cal{S}}}}+\frac{\Lambda}{3}\hat{{\cal{S}}}^{2}\right)^{1/2}}\
,$ (86a) where, as before, $t_{B}(z)$ is the ‘big bang time’ and
$\dot{{\cal{S}}}=\left(-k_{0}+\frac{2\mu}{{\cal{S}}}+\frac{\Lambda}{3}{\cal{S}}^{2}\right)^{1/2}\
.$ (86b) Differentiating (86a) (and substituting (86b)) yields
$\frac{{\cal{S}}_{,z}}{\dot{{\cal{S}}}}-\mu_{,z}\int^{{\cal{S}}}_{0}\frac{{\hbox{d}}\hat{{\cal{S}}}}{\hat{{\cal{S}}}(\dot{\hat{{\cal{S}}}})^{3}}=-t_{B,z}\
.$ (86c) The integral above can be rewritten as
$\int^{{\cal{S}}}_{0}\frac{{\hbox{d}}\hat{{\cal{S}}}}{\hat{{\cal{S}}}(\dot{\hat{{\cal{S}}}})^{3}}=\frac{1}{3\mu}\left\\{k_{0}\int^{{\cal{S}}}_{0}\frac{{\hbox{d}}\hat{{\cal{S}}}}{(\dot{\hat{{\cal{S}}}})^{3}}+\left(\frac{\dot{{\cal{S}}}}{{\cal{S}}}\right)^{-1}\right\\}\
.$ (86d) Then, from (86c) and (86d), we obtain:
$\frac{{\cal{S}}_{,z}}{{\cal{S}}}=-t_{B,z}\left(\frac{\dot{{\cal{S}}}}{{\cal{S}}}\right)+k_{0}\frac{\mu_{,z}}{3\mu}\left(\frac{\dot{{\cal{S}}}}{{\cal{S}}}\int^{{\cal{S}}}_{0}\frac{{\hbox{d}}\hat{{\cal{S}}}}{(\dot{\hat{{\cal{S}}}})^{3}}\right)+\frac{\mu_{,z}}{3\mu}\
,$ (86e) where we can identify the growing and decaying modes on the local
background described by Eq. (14) PeacockBookCosPhys :
$f_{-}=\frac{\dot{{\cal{S}}}}{{\cal{S}}}\ ,\qquad
f_{+}=\frac{\dot{{\cal{S}}}}{{\cal{S}}}\int^{{\cal{S}}}_{0}\frac{{\hbox{d}}\hat{{\cal{S}}}}{(\dot{\hat{{\cal{S}}}})^{3}}\
.$ (86f) In addition to Eq. (16), these functions satisfy
${\cal{S}}\,\dot{{\cal{S}}}\,\dot{f}_{\pm}-\left(k_{0}-\frac{3\mu}{{\cal{S}}}\right)\,f_{\pm}=\alpha_{\pm}\
,$ (86g)

with $\alpha_{+}=1$ and $\alpha_{-}=0$.

Finally, the metric coefficient ${\cal{G}}$ can be cast into

${\cal{G}}=f\frac{{\cal{S}}_{,z}}{{\cal{S}}}+f\nu_{,z}={\cal{A}}-\left(\beta_{+}f_{+}+\beta_{-}f_{-}\right)\
,$ (87)

where

$\displaystyle\beta_{+}\equiv-f\frac{k_{0}\mu_{,z}}{3\mu}\
,\quad\beta_{-}\equiv ft_{B,z}\ ;$ (88a)
$\displaystyle{\cal{A}}=f\nu_{,z}-k_{0}\beta_{+}=f_{,z}-f\frac{{\cal{E}}_{,z}}{{\cal{E}}}-k_{0}\beta_{+}\
.$ (88b)

Note that our definition of $\beta_{-}$ differs from the Goode and Wainwright
choice by a multiplicative function constant in time, $6\mu$.

##  Appendix C GRZA model equations

Let us obtain the evolution equations for the deformation field. By direct
substitution, we find that

$\displaystyle\frac{1}{2J}\epsilon_{abc}\epsilon^{ikl}\left(\dot{\eta}^{a}_{\
j}\eta^{b}_{\ k}\eta^{c}_{\
l}\right)^{\hbox{\bigastfont.}}=\frac{1}{\mathfrak{J}}\left(\ddot{\hat{P}}^{i}_{\>\>j}+3\hat{H}\dot{\hat{P}}^{i}_{\>\>j}\right)\qquad\qquad$
$\displaystyle\qquad\qquad\qquad+\left(\frac{\hat{H}\dot{\hat{P}}}{\mathfrak{J}}+2\hat{H}^{2}+\frac{\ddot{A}}{A}\right)\delta^{i}_{\>\>j}+\mathcal{O}(\hat{P}^{2})\
;$ (89)
$\displaystyle\frac{1}{2J}\epsilon_{abc}\epsilon^{ik\ell}\ddot{\eta}^{a}_{\
i}\eta^{b}_{\ k}\eta^{c}_{\
\ell}=\frac{1}{\mathfrak{J}}\left(\ddot{\hat{P}}+2\hat{H}\dot{\hat{P}}\right)+3\frac{\ddot{A}}{A}+\mathcal{O}(\hat{P}^{2})\
.$ (90)

Using (9a) and (9c) and keeping only the linear terms we obtain:

$\displaystyle-{\cal{R}}^{i}_{\
j}=\frac{1}{\mathfrak{J}}\left(\ddot{\hat{P}}^{i}_{\>\>j}+3\hat{H}\dot{\hat{P}}^{i}_{\>\>j}\right)$
$\displaystyle\quad+\left[\frac{1}{\mathfrak{J}}\left(\ddot{\hat{P}}+3\hat{H}\dot{\hat{P}}\right)+4\frac{\ddot{A}}{A}+2\left(\frac{\dot{A}}{A}\right)^{2}-2\Lambda\right]\delta^{i}_{\>\>j}$
$\displaystyle\quad=\frac{1}{\mathfrak{J}}\left(\ddot{\hat{P}}^{i}_{\>\>j}+3\hat{H}\dot{\hat{P}}^{i}_{\>\>j}\right)+\left[\frac{1}{\mathfrak{J}}\left(\ddot{\hat{P}}+3\hat{H}\dot{\hat{P}}\right)-2\frac{\hat{k}}{A^{2}}\right]\delta^{i}_{\>\>j}\
.$ (91)

Taking the trace leads to the evolution equation for $\hat{P}$:

$\ddot{\hat{P}}+3\hat{H}\dot{\hat{P}}=-\frac{\mathfrak{J}}{4}\left({\cal{R}}-6\frac{\hat{k}}{A^{2}}\right)\simeq-\frac{1}{4}\left({\cal{R}}-6\frac{\hat{k}}{A^{2}}\right)\
.$ (92)

Note that ${\cal{R}}-6\frac{\hat{k}}{A^{2}}$ is of the order of $\hat{P}$.
From now on (as in the main text), we will also assume that the Ricci tensor
and $\hat{k}/A^{2}$ are of the order of the deformation field as well.
Substituting the above relation into Eq. (91) yields (37c). Multiplying this
equation by the metric tensor $\mathbf{g}^{\hbox{\tiny{(1)}}}$, contracting
indices and neglecting the nonlinear terms, we get

$\left(\ddot{\hat{P}}_{ij}+3\hat{H}\dot{\hat{P}}_{ij}\right)=-\left[\frac{{\cal{R}}_{ij}}{A^{2}}-\frac{1}{4}\left({\cal{R}}+2\frac{\hat{k}}{A^{2}}\right)G_{ij}\right]\
,$ (93)

Next, we decompose the deformation field into the trace, tracefree symmetric,
and antisymmetric parts, Eq. (39a). Since $\dot{{\mathfrak{P}}}_{ij}=0$,
$\dot{\hat{P}}_{ij}=\frac{1}{3}\dot{\hat{P}}G_{ij}+\dot{\hat{\Pi}}_{ij}$ and
$\ddot{\hat{P}}_{ij}=\frac{1}{3}\ddot{\hat{P}}G_{ij}+\ddot{\hat{\Pi}}_{ij}$,
the evolution for the tracefree part simplifies to (40c). For Szekeres models,
this equation is equivalent to (92).

###  C.1 Hamiltonian constraint

Let us turn now to the Hamiltonian constraint; the left-hand term of Eq. (3c)
can be expanded to give

$\frac{1}{2J}\epsilon_{abc}\epsilon^{mjk}\dot{\eta}^{a}_{\ m}\dot{\eta}^{b}_{\
j}\eta^{c}_{\
k}=3\left(\frac{\dot{A}}{A}\right)^{2}+2\hat{H}\frac{\dot{\hat{P}}}{\mathfrak{J}}+\mathcal{O}(\hat{P}^{2})\
.\qquad$ (94)

The density can be rewritten as

$\varrho=\frac{\varrho_{\textrm{\tiny{ini}}}}{J}=\frac{\hat{\varrho}_{b}}{A^{3}\mathfrak{J}}+\frac{\hat{\varrho}_{b}\hat{\delta}}{A^{3}\mathfrak{J}}=\frac{\hat{\varrho}_{b}}{A^{3}\mathfrak{J}}-\frac{\hat{W}}{4\pi
A^{3}\mathfrak{J}}\ ,$ (95)

with $\hat{W}=-4\pi\hat{\varrho}_{b}\,\hat{\delta}$. Finally, Eq. (37d)
results from substituting the previous expressions into (3c) and using (22).

###  C.2 Connection coefficients

Since in this paper we only use the zero-order spatial Christoffel symbols, we
will not compute the higher order quantities. Injecting the zero-order spatial
metric ($g^{\scriptscriptstyle{(0)}}_{ij}=A^{2}G_{ij}$) and its inverse
($g^{\scriptscriptstyle{(0)}ij}=A^{-2}G^{ij}$), into the definition of the
Christoffel symbols,we obtain:

$\displaystyle{}^{\hbox{\tiny{(0)}}}\Gamma^{i}_{\
kl}=\\!{}^{\hbox{\tiny{(0)}}}\overline{\Gamma}^{i}_{\
kl}+\\!{}^{\hbox{\tiny{(0)}}}\widetilde{\Gamma}^{i}_{\ kl},\qquad\hbox{with}$
(96) $\displaystyle{}^{\hbox{\tiny{(0)}}}\overline{\Gamma}^{i}_{\
kl}=\frac{A_{|l}}{A}\delta^{i}_{\>\>k}+\frac{A_{|k}}{A}\delta^{i}_{\>\>l}-\frac{A_{|m}}{A}G^{im}G_{kl}\
;$ (97) $\displaystyle{}^{\hbox{\tiny{(0)}}}\widetilde{\Gamma}^{i}_{\
kl}=\frac{G^{im}}{2}\left(G_{mk|l}+G_{lm|k}-G_{kl|m}\right)\ .$ (98)

Note that ${}^{\hbox{\tiny{(0)}}}\widetilde{\Gamma}^{i}_{\ kl}$ are the
Christoffel symbols for $G_{ij}$:

${}^{\hbox{\tiny{(0)}}}\widetilde{\Gamma}^{i}_{\
kl}=\Gamma^{i}_{\>\>jk}(\mathbf{G})\ .$ (99)

In the context of large-scale structure formation, it is convenient to assume
that the initial metric is close to the flat FLRW line-element in Cartesian
coordinates:

$G_{ij}\simeq\delta_{ij}+\delta G_{ij}\ ,$ (100)

where the quantities $\delta G_{ij}$ are small and of the order of magnitude
of the deformation field. Note that, in general, the Szekeres metric does not
admit a splitting of the form (100). This is due to the fact that, in Szekeres
models, FLRW emerges in an unfamiliar coordinate system. In this work, we have
adopted a general Gram’s matrix to represent the initial metric, which allows
for recovering the complete Szekeres family under a proper selection of the
initial data.

###  C.3 Momentum constraints

Let us consider the momentum constraints, Eq. (3d), and note that

$\displaystyle\lambda^{i}_{\>\>j}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{J}\epsilon_{abc}\epsilon^{ikl}\dot{\eta}^{a}_{\
j}\eta^{b}_{\ k}\eta^{c}_{\
l}=2\left(\hat{H}\delta^{i}_{\>\>j}+\dot{\hat{P}}^{i}_{\>\>j}/\mathfrak{J}\right)+\mathcal{O}(\hat{P}^{2})\
;$ $\displaystyle\hat{\lambda}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{J}\epsilon_{abc}\epsilon^{ikl}\dot{\eta}^{a}_{\
i}\eta^{b}_{\ k}\eta^{c}_{\
l}=6\hat{H}+2\dot{\hat{P}}/\mathfrak{J}+\mathcal{O}(\hat{P}^{2})\ .$ (102)

Then, let us compute each term of the covariant derivative
($\lambda^{i}_{\>\>j||i}=\lambda^{i}_{\>\>j|i}+\Gamma^{i}_{\>\>il}\lambda^{l}_{\>\>j}-\Gamma^{l}_{\>\>ij}\lambda^{i}_{\>\>l}$)
on the left-hand side of Eqs. (3d),

$\displaystyle\lambda^{i}_{\>\>j|i}$ $\displaystyle=$ $\displaystyle
2\frac{\dot{A}_{|j}}{A}-2\hat{H}\frac{A_{|j}}{A}+2\frac{\dot{\hat{P}}^{i}_{\>\>j|i}}{\mathfrak{J}}+\mathcal{O}(\hat{P}^{2})\
;$ (103) $\displaystyle\Gamma^{i}_{\>\>il}\lambda^{l}_{\>\>j}$
$\displaystyle-$ $\displaystyle\Gamma^{l}_{\>\>ij}\lambda^{i}_{\>\>l}$ (104)
$\displaystyle=$ $\displaystyle
2\hat{H}\left(\Gamma^{i}_{\>\>il}\delta^{l}_{\>\>j}-\Gamma^{l}_{\>\>ij}\delta^{i}_{\>\>l}\right)+\frac{2}{\mathfrak{J}}\left(\Gamma^{i}_{\>\>il}\dot{\hat{P}}^{l}_{\>\>j}-\Gamma^{l}_{\>\>ij}\dot{\hat{P}}^{i}_{\>\>l}\right)$
$\displaystyle=$
$\displaystyle\frac{2}{\mathfrak{J}}\left(\\!{}^{\hbox{\tiny{(0)}}}\Gamma^{i}_{\>\>il}\dot{\hat{P}}^{l}_{\>\>j}-\\!{}^{\hbox{\tiny{(0)}}}\Gamma^{l}_{\>\>ij}\dot{\hat{P}}^{i}_{\>\>l}\right)+\mathcal{O}(\hat{P}^{2})\
.\quad$

Putting all these results together, we obtain:

$\displaystyle\left(\tfrac{1}{J}\epsilon_{abc}\epsilon^{ikl}\dot{\eta}^{a}_{\
j}\eta^{b}_{\ k}\eta^{c}_{\
l}\right)_{||i}=2\frac{\dot{A}_{|j}}{A}-2\hat{H}\frac{A_{|j}}{A}$
$\displaystyle\qquad\qquad\qquad\qquad+2\frac{\dot{\hat{P}}^{i}_{\>\>j|i}}{\mathfrak{J}}-4\zeta_{j}+\mathcal{O}(\hat{P}^{2})\
;$ (105)
$\displaystyle\left(\tfrac{1}{J}\epsilon_{abc}\epsilon^{ikl}\dot{\eta}^{a}_{\
i}\eta^{b}_{\ k}\eta^{c}_{\
l}\right)_{|j}=6\frac{\dot{A}_{|j}}{A}-6\hat{H}\frac{A_{|j}}{A}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+2\frac{\dot{\hat{P}}^{i}_{\>\>i|j}}{\mathfrak{J}}+\mathcal{O}(\hat{P}^{2})\
,$ (106)

where
$\zeta_{j}\equiv-\frac{1}{2\mathfrak{J}}\left(\\!{}^{\hbox{\tiny{(0)}}}\Gamma^{i}_{\>\>il}\dot{\hat{P}}^{l}_{\>\>j}-\\!{}^{\hbox{\tiny{(0)}}}\Gamma^{l}_{\>\>ij}\dot{\hat{P}}^{i}_{\>\>l}\right)$;
then, we obtain (37e), where $\zeta_{j}$ can be written out to give Eq. (42).

Using (39a), the term $\dot{\hat{P}}^{i}_{\>\>[i|j]}$ can be rewritten as
follows:

$\displaystyle\dot{\hat{P}}^{i}_{\>\>[i|j]}$
$\displaystyle=\frac{1}{2}\left(\dot{\hat{P}}^{i}_{\>\>i|j}-\dot{\hat{P}}^{i}_{\>\>j|i}\right)=\frac{1}{2}\left(\dot{\hat{P}}_{|j}-\dot{\hat{P}}^{i}_{\>\>j|i}\right)$
$\displaystyle=\frac{1}{2}\left[\dot{\hat{P}}_{|j}-\left(\frac{1}{3}\dot{\hat{P}}\delta^{i}_{\>\>j}+\dot{\hat{\Pi}}^{i}_{\>\>j}\right)_{|i}\right]=\frac{1}{3}\dot{\hat{P}}_{|j}-\frac{1}{2}\dot{\hat{\Pi}}^{i}_{\>\>j|i}\
.$

Finally, we obtain the linear expression for the momentum constraints, Eq.
(40e). Under the assumption (100), $\zeta_{j}$ simplifies to

$\zeta_{j}=\frac{1}{2}\left(\frac{A_{|j}}{A}\dot{\hat{P}}-3\frac{A_{|l}}{A}\dot{\hat{P}}^{l}_{\>\>j}\right)\;.$
(109)

##  Appendix D Backreaction and emergence of spatial curvature

Consider a compact domain ${{\cal D}}$, the volume scale factor of which being
defined as follows:

$a_{{\cal D}}(t)=\left(V_{{\cal D}}/V_{{\cal D}_{i}}\right)^{\frac{1}{3}}\
;\quad V_{{\cal D}_{i}}=V_{{\cal D}}(t_{i})\ .$ (110)

Then, $H_{\cal D}=\dot{a}_{\cal D}/a_{\cal D}$ represents the volume Hubble
expansion rate.

In a general spacetime, the volume expansion and volume acceleration
generalize the kinematic laws of a homogeneous and isotropic FLRW spacetime.
In this framework, emerging curvature is a result of its coupling to the
kinematical backreaction term (involving the first two principal scalar
invariants of the expansion tensor) via Eq. (113d) below Buchert2000GERG ;
BuchertFocus :

$\displaystyle{\cal Q}_{{\cal D}}$ $\displaystyle\equiv$ $\displaystyle
2\left\langle{\rm II}\right\rangle_{\cal D}-\frac{2}{3}\left\langle{\rm
I}\right\rangle_{\cal D}^{2}$ (111a) $\displaystyle=$
$\displaystyle\frac{2}{3}\left(\left\langle\Theta^{2}\right\rangle_{\cal
D}-\left\langle\Theta\right\rangle_{\cal
D}^{2}\right)-2\left\langle\sigma^{2}\right\rangle_{\cal D}\ .$ (111b)

In these equations,
$\sigma^{2}\equiv\frac{1}{2}\sigma^{a}_{\>\>b}\sigma^{b}_{\>\>a}$ while ${\rm
I}$ and ${\rm II}$ (and ${\rm III}$) represent the principal scalars invariant
of the expansion tensor, ${\rm I}=\Theta^{i}_{\>\>i}=\Theta$, ${\rm
II}=\frac{1}{2}\left(\Theta^{2}-\Theta^{i}_{\ j}\Theta^{j}_{\
i}\right)=\frac{1}{3}\Theta^{2}-\sigma^{2}$, which are defined as follows for
an arbitrary matrix $B^{i}_{\>\>j}$

$\displaystyle{\rm I}(B^{i}_{\>\>j})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\epsilon_{abc}\epsilon^{ikl}_{\>\>}B^{a}_{\>\>i}\delta^{b}_{\>\>j}\delta^{c}_{\>\>k}\
;$ (112a) $\displaystyle{\rm II}(B^{i}_{\>\>j})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\epsilon_{abc}\epsilon^{ikl}_{\>\>}B^{a}_{\>\>i}B^{b}_{\>\>j}\delta^{c}_{\>\>k}\
;$ (112b) $\displaystyle{\rm III}(B^{i}_{\>\>j})$ $\displaystyle=$
$\displaystyle\frac{1}{6}\epsilon_{abc}\epsilon^{ikl}_{\>\>}B^{a}_{\>\>i}B^{b}_{\>\>j}B^{c}_{\>\>k}\
.$ (112c)

Introducing the deviation from a (scale-dependent) constant curvature model as
a new backreaction variable, ${\cal{W}}_{\cal
D}:=\left\langle{\cal{R}}\right\rangle_{\cal D}-6k_{{\cal D}_{\rm i}}/a_{\cal
D}^{2}$, the volume-averaged set of equations reads Buchert2000GERG :

$\displaystyle\left(\frac{\dot{a}_{\cal D}}{a_{\cal
D}}\right)^{2}-\frac{\Lambda}{3}-\frac{8\pi}{3}\left\langle\varrho\right\rangle_{\cal
D}+\frac{k_{{\cal D}_{\rm i}}}{a_{\cal D}^{2}}=-\frac{1}{6}({\cal Q}_{\cal
D}+{\cal{W}}_{\cal D})\,;\qquad$ (113a) $\displaystyle\frac{\ddot{a}_{\cal
D}}{a_{\cal
D}}-\frac{\Lambda}{3}+\frac{4\pi}{3}\left\langle\varrho\right\rangle_{\cal
D}=\frac{1}{3}{\cal Q}_{{\cal D}}\,;\qquad$ (113b)
$\displaystyle\left\langle\varrho\right\rangle^{\hbox{\bigastfont.}}_{\cal
D}+3\frac{\dot{a}_{\cal D}}{a_{\cal D}}\left\langle\varrho\right\rangle_{\cal
D}=0\,;\qquad$ (113c) $\displaystyle\frac{1}{a_{\cal D}^{6}}\left({\cal
Q}_{{\cal D}}a_{\cal D}^{6}\right)^{\hbox{\bigastfont.}}+\frac{1}{a_{\cal
D}^{2}}\left({\cal{W}}_{\cal D}a_{\cal
D}^{2}\right)^{\hbox{\bigastfont.}}=0,\qquad$ (113d)

where (113b), (113c) and (113c) correspond to the averaged energy constraint,
Raychaudhuri’s and rest mass conservation equations, respectively. Equation
(113d) is not independent; it emerges as a condition of integrability in order
to obtain (113b) through the time-derivative of (113a).

##  Appendix E Numerical example in the language of exact solutions

The initial data of the family of models examined in Section VI is set up in
the Hellaby’s parametrization, where the physical interpretation of the
variables is clearer than in the GW’s. The initial areal distance is taken
equal to the ‘radial’ comoving coordinate

$\Phi(t_{\textrm{\tiny{ini}}},z)=z\ .$ (114)

The solution of interest is contained in the domain ${\cal D}:0\leq z\leq
z_{1},-\infty<x,y<+\infty$. Here, $z_{1}$ sets the length scale of the initial
configuration.

The free function $M(z)$ is chosen as

$M(z)=\frac{4\pi}{3}\varrho_{bi}z^{3}\ ,$ (115)

where $\varrho_{bi}$ is the density of the homogeneous universe at the initial
time (i.e., the last scattering time) and $\hat{k}(z)$ is specified as in Eq.
(58).

Then, we introduce the dimensionless variables (57) using the Hubble parameter
at the last scattering time as the timescale and the characteristic length
$\ell=2\times 10^{-2}$. (This characteristic length expands to supercluster
scales at present time.) This solution is rewritten in the GW language using
the results of Appendix B.3, and then rescaled as sketched in Section III.1.
In terms of the Szekeres functions and the perturbation, the line-element
reads:

$\displaystyle{\hbox{d}}s^{2}=-{\hbox{d}}t^{2}+A^{2}\bigg{[}\left({\cal{A}}-\hat{{\cal{F}}}\right)^{2}\left(\frac{{\cal{W}}}{\chi}\right)^{2}{\hbox{d}}\xi^{2}$
$\displaystyle+\xi^{2}\mathcal{E}^{-2}\left({\hbox{d}}x^{2}+{\hbox{d}}y^{2}\right)\bigg{]}\
,$ (116)

where we have used that $\mathcal{E}$ was defined in (70d).

Inserting Eq. (116) into (68), taking back the rescaling of the solution and
using the Szekeres exact relations, we obtain the following expression for the
violation of the Hamiltonian constraint:

$\displaystyle\mathbb{H}(\bm{r},t):=\frac{2}{A^{2}\xi^{2}{\cal{W}}^{2}}\bigg{[}-\frac{1}{2}\left(1+\frac{{\cal{G}}}{\hat{{\cal{G}}}}\right)-\frac{1}{\hat{{\cal{G}}}}f\nu_{,\xi}$
$\displaystyle+\epsilon
f_{,\xi}+\frac{f{\cal{G}}_{,\xi}}{\hat{{\cal{G}}}^{2}}\bigg{]}\left(1-\frac{{\cal{G}}}{\hat{{\cal{G}}}}\right)$
$\displaystyle-\frac{\mathcal{E}^{2}}{A^{2}\xi^{2}\hat{{\cal{G}}}}\bigg{[}\left(k_{0}+f_{+}\right)\left(\delta\beta_{+,xx}+\delta\beta_{+,yy}\right)$
$\displaystyle+f_{-}\left(\delta\beta_{-,xx}+\delta\beta_{-,yy}\right)\bigg{]}$
$\displaystyle+\frac{2f}{A^{2}\xi^{2}{\cal{W}}^{2}\hat{{\cal{G}}}^{2}}\left(\frac{{\cal{G}}}{\hat{{\cal{G}}}}\right)\bigg{[}\delta\beta_{+}(k_{0}+f_{+})+\delta\beta_{-}f_{-}\bigg{]}_{,\xi}.\quad$
(117)

## References

  * (1) J. Adamek, C. Barrera-Hinojosa, M. Bruni, B. Li, H. Macpherson and J.B. Mertens, Numerical solutions to Einstein’s equations in a shearing-dust universe: a code comparison. Class. Quant. Grav. 37, 154001 (2020) [arXiv:2003.08014]
  * (2) A. Alles, T. Buchert, F. Al Roumi and A. Wiegand, Lagrangian theory of structure formation in relativistic cosmology. III. Gravitoelectric perturbation and solution schemes at any order. Phys. Rev. D 92, 023512 (2015) [arXiv:1503.02566]
  * (3) F. Al Roumi, T. Buchert and A. Wiegand, Lagrangian theory of structure formation in relativistic cosmology. IV. Lagrangian approach to gravitational waves. Phys. Rev. D 96, 123538 (2017) [arXiv:1711.01597]
  * (4) A. Barnes and R.R. Rowlingson, Irrotational perfect fluids with a purely electric Weyl tensor. Class. Quant. Grav. 6, 949 (1989)
  * (5) H. Baysal, Inhomogeneous cosmologic models with heat flow. Astrophys. Space Sci. 196, 345 (1992)
  * (6) S. Bildhauer, T. Buchert and M. Kasai, Solutions in Newtonian cosmology – the pancake theory with cosmological constant. Astron. Astrophys. 263, 23 (1992)
  * (7) K. Bolejko, Structure formation in quasispherical Szekeres model. Phys. Rev. D 73, 123508 (2006) [arXiv:astro-ph/0604490]
  * (8) K. Bolejko, Evolution of cosmic structures in different environments in the quasispherical Szekeres model. Phys. Rev. D 75, 043508 (2007) [arXiv:astro-ph/0610292]
  * (9) K. Bolejko, The Szekeres Swiss Cheese model and the CMB observations. Gen. Rel. Grav. 41, 1737 (2009) [arXiv:0804.1846]
  * (10) K. Bolejko, Relativistic numerical cosmology with Silent Universes. Class. Quant. Grav. 35, 024003 (2017) [arXiv:1708.09143]
  * (11) K. Bolejko, Emerging spatial curvature can resolve the tension between high-redshift CMB and low-redshift distance ladder measurements of the Hubble constant. Phys. Rev. D 97, 103529 (2018) [arXiv:1712.02967]
  * (12) K. Bolejko and M.N. Célérier, Szekeres Swiss-Cheese model and supernova observations. Phys. Rev. D 82, 103510 (2010) [arXiv:1005.2584]
  * (13) K. Bolejko, M.N. Célérier and A. Krasiński, Inhomogeneous cosmological models: exact solutions and their applications. Class. Quant. Grav. 28, 164002 (2011) [arXiv:1102.1449]
  * (14) K. Bolejko, A. Krasiński, C. Hellaby and M.N. Célérier, Structures in the Universe by exact methods: formation, evolution, interactions. Cambridge University Press (2009). ISBN, 9780521769143
  * (15) K. Bolejko and J.J. Ostrowski, Environment dependence of the growth of the most massive objects in the Universe. Phys. Rev. D 99, 124036 (2019) [arXiv:1805.11047]
  * (16) K. Bolejko and R.A. Sussman, Cosmic spherical void via coarse-graining and averaging non-spherical structures. Phys. Lett. B 697, 265 (2011) [arXiv:1008.3420]
  * (17) W.B. Bonnor, A.H. Sulaiman and N. Tomimura, Szekeres’s Space-Times Have No Killing Vectors. Gen. Rel. Grav. 8, 549 (1977)
  * (18) W.B. Bonnor and N. Tomimura, Evolution of Szekeres’s cosmological models. Mon. Not. Roy. Astron. Soc. 175, 85 (1976)
  * (19) M. Bruni, S. Matarrese and O. Pantano, Dynamics of Silent Universes. Astroph. J 445, 958 (1995) [arXiv:astro-ph/9406068]
  * (20) L. Brunswic and T. Buchert, Gauss-Bonnet-Chern approach to the averaged Universe. Class. Quant. Grav. 37, 215022 (2020) [arXiv:2002.08336]
  * (21) T. Buchert and G. Götz, A class of solutions for self-gravitating dust in Newtonian gravity. J. Math. Phys. 28, 2714 (1987)
  * (22) T. Buchert, A class of solutions in Newtonian cosmology and the pancake theory. Astron. Astrophys. 223, 9 (1989)
  * (23) T. Buchert, On Average Properties of Inhomogeneous Fluids in General Relativity: Dust Cosmologies. Gen. Rel. Grav. 32, 105 (2000) [arXiv:gr-qc/9906015]
  * (24) T. Buchert, On average properties of inhomogeneous fluids in general relativity II: Perfect fluid cosmologies. Gen. Rel. Grav. 33, 1381 (2001) [arXiv:gr-qc/0102049]
  * (25) T. Buchert, Dark Energy from structure: a status report. Gen. Rel. Grav. 40, 467 (2008) [arXiv:0707.2153]
  * (26) T. Buchert, Toward physical cosmology: focus on inhomogeneous geometry and its non-perturbative effects. Class. Quant. Grav. 28, 164007 (2011) [arXiv:1103.2016]
  * (27) T. Buchert, On backreaction in Newtonian cosmology. Mon. Not. Roy. Astron. Soc. 473, L46 (2018) [arXiv:1704.00703]
  * (28) T. Buchert and M. Bartelmann, High-spatial resolution in three dimensions: a challenge for large-scale structure formation models. Astron. Astrophys. 251, 389 (1991)
  * (29) T. Buchert and M. Carfora, On the curvature of the present-day Universe. Class. Quant. Grav. 25, 195001 (2008) [arXiv:0803.1401]
  * (30) T. Buchert, I. Delgado Gaspar and J.J. Ostrowski, On general-relativistic Lagrangian perturbation theory and its non-perturbative generalization. Universe 8, 583 (2022) [arXiv:2209.13417]
  * (31) T. Buchert and J. Ehlers, Averaging inhomogeneous Newtonian cosmologies. Astron. Astrophys. 320, 1 (1997) [arXiv:astro-ph/9510056]
  * (32) T. Buchert, G. Karakatsanis, R. Klaffl and P. Schiller, The performance of Lagrangian perturbation schemes at high resolution. Astron. Astrophys. 318, 1 (1997) [arXiv:astro-ph/9510133]
  * (33) T. Buchert, P. Mourier and X. Roy, On average properties of inhomogeneous fluids in general relativity III: general fluid cosmologies. Gen. Rel. Grav. 52, 27 (2020) [arXiv:1912.04213]
  * (34) T. Buchert, C. Nayet and A. Wiegand, Lagrangian theory of structure formation in relativistic cosmology. II. Average properties of a generic evolution model. Phys. Rev. D 87, 123503 (2013) [arXiv:1303.6193]
  * (35) T. Buchert and M. Ostermann, Lagrangian theory of structure formation in relativistic cosmology. I. Lagrangian framework and definition of a nonperturbative approximation. Phys. Rev. D 86, 023520 (2012) [arXiv:1203.6263]
  * (36) T. Buchert and S. Räsänen, Backreaction in Late-Time Cosmology. Annu. Rev. Nucl. Part. Sci. 62, 57 (2012) [arXiv:1112.5335]
  * (37) R.G. Buckley and E.M. Schlegel, CMB dipoles and other low-order multipoles in the quasispherical Szekeres model. Phys. Rev. D 87, 023524 (2013) [arXiv:1907.08684]
  * (38) P. Bull, T. Clifton and P.G. Ferreira, The Kinematic Sunyaev-Zel’dovich effect as a test of general radial inhomogeneity in Lemaître-Tolman-Bondi cosmology. Phys. Rev. D 85, 024002 (2012) [arXiv:1108.2222]
  * (39) T. Clifton and R.A. Sussman, Cosmological backreaction in spherical and plane symmetric dust-filled space-times. Class. Quant. Grav. 36, 205004 (2019) [arXiv:1904.02557]
  * (40) A.A. Coley, N. Layden and D. McNutt, An invariant characterization of the quasi-spherical Szekeres dust models. Gen. Rel. Grav. 51, 164 (2019) [arXiv:2103.05305]
  * (41) I. Delgado Gaspar and T. Buchert, Lagrangian theory of structure formation in relativistic cosmology. VI. Comparison with Szekeres exact solutions. Phys. Rev. D 103, 023513 (2021) [arXiv:2009.06339]
  * (42) I. Delgado Gaspar, J.C. Hidalgo, R.A. Sussman and I. Quiros, Black hole formation from the gravitational collapse of a nonspherical network of structures. Phys. Rev. D 97, 104029 (2018) [arXiv:1802.09123]
  * (43) I. Delgado Gaspar, R.A. Sussman, D.D. Mcnutt and A.A. Coley, Comment on Szekeres universes with homogeneous scalar fields. Eur. Phys. J. C. 81, 310 (2021) [arXiv:2103.06862]
  * (44) A.G. Doroshkevich, V.S. Ryabenkii and S.F. Shandarin, Non-linear theory of development of potential perturbations. Astrofizika 9, 257 (1973); Astrophysics 9, 144 (1973)
  * (45) H. van Elst, C. Uggla, W.M Lesame, G.F.R Ellis and R. Maartens Integrability of irrotational silent cosmological models, Class. Quant. Grav. 14, 1151 (1997) [arXiv:gr-qc/9611002]
  * (46) I. Georg and C. Hellaby, Symmetry and Equivalence in Szekeres Models. Phys. Rev. D 95, 124016 (2017) [arXiv:1702.05347]
  * (47) S.W. Goode, Spatially inhomogeneous cosmologies with heat flow. Class. Quant. Grav. 3, 1247 (1986)
  * (48) S.W. Goode and J. Wainwright, Singularities and evolution of the Szekeres cosmological models. Phys. Rev. D 26, 3315 (1982)
  * (49) S.W. Goode and J. Wainwright, Friedmann-like singularities in Szekeres’ cosmological models. Mon. Not. Roy. Astron. Soc. 198, 83 (1982)
  * (50) T. Harada and S. Jhingan, Spherical and nonspherical models of primordial black hole formation: exact solutions. Prog. Theor. Exp. Phys. 9, 093E04 (2016)
  * (51) A. Heinesen and T. Buchert, Solving the curvature and Hubble parameter inconsistencies through structure formation-induced curvature. Class. Quant. Grav. 37, 164001 (2020) (Focus issue on the Hubble constant tension); Erratum (2020): Class. Quant. Grav. 37, 229601 [arXiv:2002.10831]
  * (52) C. Hellaby, The null and KS limits of the Szekeres model. Class. Quant. Grav. 13, 2537 (1996)
  * (53) C. Hellaby, Frame Rotation in the Szekeres Spacetimes. Class. Quant. Grav. 34, 145006 (2017) [arXiv:1706.00622]
  * (54) C. Hellaby and A. Krasiński, You can’t get through Szekeres wormholes: Or, regularity, topology and causality in quasispherical Szekeres models. Phys. Rev. D 66, 084011 (2002) [arXiv:gr-qc/0206052]
  * (55) C. Hellaby and A. Krasiński, Physical and geometrical interpretation of the $\epsilon\leq 0$ Szekeres models. Phys. Rev. D 77, 023529 (2008) [arXiv:0710.2171]
  * (56) M. Ishak and A. Peel, The growth of structure in the Szekeres inhomogeneous cosmological models and the matter-dominated era. Phys. Rev. D 85, 083502 (2012) [arXiv:1104.2590]
  * (57) M. Ishak, A. Peel and M. A. Troxel, Stringent Restriction from the Growth of Large-Scale Structure on Apparent Acceleration in Inhomogeneous Cosmological Models. Phys. Rev. Lett. 111, 251302 (2013) [arXiv:1307.0723]
  * (58) M. Ishak, J. Richardson, D. Whittington and D. Garred, Dark Energy or Apparent Acceleration Due to a Relativistic Cosmological Model More Complex than FLRW? Phys. Rev. D 84, 089902 (2011) [arXiv:0708.2943]
  * (59) P.S. Joshi, Gravitational collapse: The story so far. Pramana 55, 529 (2000) [arXiv:gr-qc/0006101]
  * (60) M. Kasai, Tetrad-based perturbative approach to inhomogeneous universes: a general relativistic version of the Zel’dovich approximation. Phys. Rev. D 52, 5605 (1995)
  * (61) S.M. Koksbang and S. Hannestad, Methods for studying the accuracy of light propagation in N-body simulations. Phys. Rev. D 91, 043508 (2015) [arXiv:1501.01413]
  * (62) S.M. Koksbang and S. Hannestad, Studying the precision of ray tracing techniques with Szekeres models. Phys. Rev. D 92, 023532 (2015) [arXiv:1506.09127]
  * (63) S.M. Koksbang, 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 (2017) [arXiv:1703.03572]
  * (64) A. Krasiński, Inhomogeneous Cosmological Models. Cambridge University Press (1997). ISBN 0521481805
  * (65) A. Krasiński and K. Bolejko, Redshift propagation equations in the $\beta^{\prime}\neq 0$ Szekeres models. Phys. Rev. D 83, 083503 (2011) [arXiv:1007.2083]
  * (66) G. Lemaître, The expanding universe. Annales Soc. Sci. Bruxelles A 53, 51-85 (1933); reprinted: Gen. Rel. Grav. 29, 641(1997)
  * (67) Y.Z. Li, P. Mourier, T. Buchert and D.L. Wiltshire, Lagrangian theory of structure formation in relativistic cosmology. V. Irrotational fluids. Phys. Rev. D 98, 043507 (2018) [arXiv:1806.05467]
  * (68) J.A.S. Lima and M.R. de Garcia Maia, Chaotic cosmology: Inhomogeneous models with heat flow. Phys. Lett. A 110, 366 (1985)
  * (69) J.A.S. Lima and M.A.S. Nobre, Inhomogeneous two-fluid cosmologies with electromagnetic field. Class. Quant. Grav. 7, 399 (1990)
  * (70) N. Meures and M. Bruni, Exact non-linear inhomogeneities in $\Lambda$CDM cosmology. Phys. Rev. D 83, 123519 (2011) [arXiv:1103.0501]
  * (71) N. Meures and M. Bruni, Redshift and distances in a $\Lambda$CDM cosmology with non-linear inhomogeneities. Mon. Not. Roy. Astron. Soc. 3, 1937 (2012) [arXiv:arXiv:1107.4433]
  * (72) D.C. Motta and N. Tomimura, Some solutions of Einstein’s equations with viscous fluids. Astrophys. Space Sci. 165, 237 (1990)
  * (73) S. Nájera and R.A. Sussman, Pancakes as opposed to Swiss cheese. Class. Quant. Grav. 38, 015016 (2021) [arXiv:2010.04027]
  * (74) B.C. Nolan, Dynamical extensions for shell-crossing singularities. Class. Quant. Grav. 20, 575 (2003) [arXiv:gr-qc/0301028]
  * (75) A. Nwankwo, M. Ishak and J. Thompson, Luminosity distance and redshift in the Szekeres inhomogeneous cosmological models. JCAP 1105, 028 (2011) [arXiv:1005.2989]
  * (76) J.A. Peacock: Cosmological physics. Cambridge University Press (1998). ISBN, 9780511804533
  * (77) A. Peel, M. Ishak and M.A. Troxel, Large-scale growth evolution in the Szekeres inhomogeneous cosmological models with comparison to growth data. Phys. Rev. D 86, 123508 (2012) [arXiv:1212.2298]
  * (78) J. Plebański and A. Krasiński, An introduction to general relativity and cosmology. Cambridge University Press (2006). ISBN 0-521-85623-X
  * (79) R.A. Sussman, Back-reaction and effective acceleration in generic LTB dust models. Class. Quant. Grav. 28, 235002 (2011) [arXiv:1102.2663]
  * (80) R.A. Sussman and K. Bolejko, A Novel approach to the dynamics of Szekeres dust models. Class. Quant. Grav. 29, 065018 (2012) [arXiv:1109.1178]
  * (81) R.A. Sussman and I. Delgado Gaspar, Multiple nonspherical structures from the extrema of Szekeres scalars. Phys. Rev. D 92, 083533 (2015) [arXiv:1508.03127]
  * (82) R.A. Sussman, I. Delgado Gaspar and J.C. Hidalgo, Coarse-grained description of cosmic structure from Szekeres models. JCAP 1603, 012 (2016) [arXiv:1507.02306]
  * (83) R.A. Sussman, J.C. Hidalgo, I. Delgado Gaspar and G. German, Nonspherical Szekeres models in the language of cosmological perturbations. Phys. Rev. D 95, 064033 (2017) [arXiv:1701.00819]
  * (84) D.A. Szafron, Inhomogeneous cosmologies: New exact solutions and their evolution. J. Math. Phys. 18, 1673 (1977)
  * (85) D.A. Szafron and C.B. Collins, A new approach to inhomogeneous cosmologies: intrinsic symmetries. II. Conformally flat slices and an invariant classification. J. Math. Phys. 20, 2354 (1979)
  * (86) P. Szekeres, Quasispherical Gravitational Collapse. Phys. Rev. D 12, 2941 (1975)
  * (87) N. Tomimura and D.C. Motta, A new class of Szekeres inhomogeneous cosmologies with imperfect fluids. Astrophys. Space Sci. 165, 243 (1990)
  * (88) N. Tomimura and I. Waga, Inhomogeneous Cosmological Models: New Exact Solutions with Dust, Isotropic Radiation and Electromagnetic Field. Astrophys. J. 317, 52 (1987)
  * (89) D. Vrba and O. Svítek, Modelling inhomogeneity in Szekeres spacetime. Gen. Rel. Grav. 46, 1808 (2014) [arXiv:1412.7203]
  * (90) J. Wainwright, Characterization of the Szekeres inhomogeneous cosmologies as algebraically special space-times. J. Math. Phys. 18, 672 (1977)
  * (91) A. Walters and C. Hellaby, Constructing realistic Szekeres models from initial and final data. Phys. Rev. D 1212, 001 (2012) [arXiv:1211.2110]
  * (92) Ya.B. Zel’dovich, Gravitational instability: an approximate theory for large density perturbations. Astron. Astrophys. 5, 84 (1970); reprinted: Astron. Astrophys. 500, 13 (2009)
  * (93) J.P. Zibin, Can decaying modes save void models for acceleration? Phys. Rev. D 84, 123508 (2011) [arXiv:1108.3068]
  * (94) J.P. Zibin and A. Moss, Linear kinetic Sunyaev-Zel’dovich effect and void models for acceleration. Class. Quant. Grav. 28, 164005 (2011) [arXiv:1105.0909]