# Second-order dust perturbations of the non-flat FLRW model with the positive
cosmological constant.

Szymon Sikora
Astronomical Observatory, Jagiellonian University, Orla 171, 30-244 Kraków,
Poland

###### Abstract

In this paper, a specific solution to the second-order cosmological
perturbation theory is given. Perturbations are performed around any FLRW
spacetime filled with dust and with a positive cosmological constant. In
particular, with a possibly non-vanishing spatial curvature. The adopted
symmetry condition allows us to simplify the equations, leaving us with a
great deal of freedom to choose the density distribution. In the result, we
get a relatively simple metric of an inhomogeneous cosmological model, which
will give a perfect tool for studying the influence of the local
inhomogeneities onto the cosmological observables.

##  1 Introduction

It is not possible to imagine cosmology without the Friedmann-Lemaître-
Robertson-Walker (FLRW) model describing a universe which is spatially
homogeneous and isotropic. However, apart from a significant success in
explaining variety of the cosmological observables by the FLRW model, there
are some issues which remain to be clarified.

In recent years, the precision of the measurements of the Lemaître-Hubble
constant $H_{0}$ has increased. Unfortunately, the estimations of its value
using different methods have become inconsistent. The CMB observations provide
the value $H_{0}=67.37\pm 0.54\,\mathrm{km/s/Mpc}$ 2020A&A…641A…6P lower than
the estimates based on standard candles $H_{0}=74.03\pm
1.42\,\mathrm{km/s/Mpc}$ 2019ApJ…876…85R . This discrepancy is known as the
Hubble tension problem. For more details one can see a review
2021CQGra..38o3001D and the references therein. Among many possible solutions,
the influence onto the local $H_{0}$ measurements by the inhomogeneity around
the observer is reliable and based on standard physics. The idea that the
matter inhomogeneities on a different scales could affect the cosmological
observations is not new and appear also in earlier works, e. g.
2008PhRvD..78h3531L . In the context of the Hubble tension, a directional
measurements of the Lemaître-Hubble constant 2016JCAP…04..036B are important.
A theoretical analysis of this phenomenon within the Szekeres solution can be
found in 2016JCAP…06..035B . We can mention also a recent attempt in
explaining quantitatively Hubble parameter inconsistencies within the
inhomogeneous cosmology 2020CQGra..37p4001H .

Another issue is the problem of the averaging. It is believed that the
homogeneous density of the FLRW model is an effect of the averaging of the
matter inhomogeneities over a large scales. However, calculation of the
averages in a curved spacetime is highly non-trivial. It could happen that the
averaged expansion of the inhomogenous universe differs from the expansion of
the corresponding FLRW model which density equals the averaged density of the
inhomogeneous model. This effect is known as the cosmological backreaction.
There exist several averaging schemes which lead to a different conclusions
about the importance of the backreaction effect. Within the Green-Wald
averaging scheme the authors conclude that the backreaction effect is
negligible 2011PhRvD..83h4020G ; 2016CQGra..33l5027G , while in the Buchert
framework some papers suggest that the backreaction could be important
2000GReGr..32..105B ; 2015CQGra..32u5021B ; 2018MNRAS.473L..46B ;
2018CQGra..35xLT02B or even it could explain the accelerated expansion of the
universe without the need of the cosmological constant 2013JCAP…10..043R .
Recently, other approaches to the averaging problem are investigated. An
interesting example is the analysis of some observables and averages over the
past light cone 2021LMaPh.111…53C ; 2022arXiv220210798B ; 2011JCAP…07..008G .

To study the effects caused by the inhomogeneities one can analyze exact
solutions to the Einstein equations, which are generalizations of the FLRW
spacetime. The most important of such models are the Lemaître-Tolman-Bondi
spacetime 1947MNRAS.107..410B and the Szekeres solution 1975CMaPh..41…55S . It
is also possible to immerse several substructures described by these solutions
into the homogeneous background within the Swiss-cheese framework
2008JCAP…06..021B ; 2013PhRvD..87l3526F ; 2013JCAP…12..051L ;
2015JCAP…10..057L ; 2021PhRvD.104d3505K . Although the Einstein equations are
very complicated and it is hard to discover new solutions, there is some
progress in this field, e. g. 2021CQGra..38a5016N found a new class of
interesting models. Another possibility is to study the evolution of the
inhomogeneities within the numerical relativity 2017PhRvD..95f4028M ;
2016JCAP…07..053A ; 2014CQGra..31w4006A ; 2022PhLB..82636911E or by using the
perturbative methods like the post-Newtonian formalism 2015PhRvD..91j3532S ;
2016PhRvD..93h9903S ; 2017JCAP…07..028S .

In this paper we focus on the cosmological perturbation theory. Many years
after original papers by Lifshitz 1946ZhETF..16..587L ; 2017GReGr..49…18L and
Bardeen 1980PhRvD..22.1882B there is still progress in this area. In
particular, there are attemts to analyze the perturbations beyond the linear
order 2014JCAP…07..037N ; 2017JCAP…10..027G ; 2017PhRvD..96j3522W . In
2009JCAP…11..012C the third-order perturbations are considered. The matter
inhomogeneities in the form of the ensemble of the point masses are studied up
to second order in 2017ApJ…845..153B . There are also investigations of
particular gauge conditions 2006PhRvD..73d4021H and a complementary approach
regarding gauge-independent variables 2005AcPPB..36.2133G .

The perturbative equations in the second order become quite complicated. For
that reason, instead of studying the perturbation theory in general, we are
looking for a specific solution for which the resulting metric tensor is as
simple as possible. In 2017PhRvD..95f3517S we found interesting solution
within the linear perturbation theory, where the matter inhomogeneities form
an infinite periodic cubic lattice on the spatially flat Einstein-de Sitter
background. A periodic lattice is a natural choice for a numerical relativity
and it is present also in the other works, e. g. in the post-Newtonian
formalism 2017JCAP…07..028S . In 2019PhRvD..99h3521S we generalize the
presented solution up to the third-order perturbations and analyze the light
propagation through this inhomogeneous model. In both papers
2017PhRvD..95f3517S ; 2019PhRvD..99h3521S a decaying mode was considered. The
subsequent paper 2021EPJC…81..208S was dedicated to the analysis of the
growing mode and the fourth-order perturbations were taken into account. Since
a high order perturbations were present, a relatively large amplitude of the
inhomogeneities was admissible. The papers 2017PhRvD..95f3517S ;
2019PhRvD..99h3521S ; 2021EPJC…81..208S considered a spatially flat background
without the cosmological constant. The main purpose of the present paper is to
generalize these results such that the considered background could be any FLRW
cosmological model filled with dust, with possibly non-vanishing spatial
curvature, and with a positive cosmological constant. For that purpose, we
consider a two-parameter perturbation theory, where the curvature parameter
$k$ is treated as the second perturbative variable. The idea of a two-
parameter perturbations is not new and appear in different contexts, e. g.
2017PhRvD..95d3503G ; 2017PhRvD..96j3508G ; 2021JCAP…08..048G .

The paper is organized as follows. In Sec. 2 we present the main idea of the
currect work. The specific solution to the perturbation theory in the first
and second order is presented in Sec. 3 and Sec. 4 respectively. Then, the
conclusions are given.

##  2 The model construction.

###  2.1 The metric.

Let us consider Cartesian-like coordinates $(t,x,y,z)$, in which the line
element reads:

$\mathrm{d}s^{2}=-\mathrm{d}t^{2}+a(t)^{2}\,c_{i\>\\!j}(\lambda,k,x^{\mu})\,\mathrm{d}x^{i}\,\mathrm{d}x^{j}\,.$
(1)

For the metric written in this way, $(t,x,y,z)$ are the Gaussian normal
coordinates. It is always possible to transform the metric into the form (1)
locally, but it is not obvious whether one could do it for the whole space-
time. The vector field $(U^{\mu})=(1,0,0,0)$ is tangent to the time-like
geodesic at each point. Due to the attractive nature of gravity, the geodesics
could be focused and eventually cross with each other. In that case, the
Gaussian normal coordinates are not suitable because the point where two
geodesics meet would have two distinct coordinate values. However, this is not
the case for the model presented in this article. As we shall see in Sect.
4.2, the expansion scalar $\theta=\nabla_{\mu}\,U^{\mu}$ is always positive
here. Since space is expanding everywhere, geodesics which start at different
points on some hypersurface of constant time never cross in the future. In
effect, the Gaussian normal coordinates are properly defined in the whole
spacetime. Due to the positive values of $\theta$, the model does not intend
to describe the formation of structures, but provides a tool for handling the
overall differences between cosmic voids and regions populated with galaxy
clusters.

The metric functions $c_{i\>\\!j}$ depend on the coordinates and two
parameters with the following meaning. For $\lambda$ equal to zero, the model
should reduce to the FLRW spacetime with arbitrary curvature parameter $k$:

$\lim\limits_{\lambda\rightarrow
0}c_{i\>\\!j}=\frac{\delta_{i\>\\!j}}{\left(1+\frac{1}{4}k\,(x^{2}+y^{2}+z^{2})\right)^{2}}\equiv\widetilde{c}_{i\>\\!j}\,.$
(2)

This means that $a(t)$ is the scale factor of the background. The energy-
momentum tensor of matter that fills spacetime is
$T^{\mu}\>\\!{}_{\nu}=\rho\,U^{\mu}\,U_{\nu}$, where $(U^{\mu})=(1,0,0,0)$.
For $\lambda\neq 0$, the density $\rho$ could be inhomogeneous and the
parameter $\lambda$ is proportional to the amplitude of the inhomogeneities.

###  2.2 The perturbations.

It is difficult to find the functions $c_{i\>\\!j}$ that meet the conditions
given above. To deal with this problem, we will consider some simplifications.
Let us assume that we can express metric functions as a power series:

$c_{i\>\\!j}(\lambda,k,x^{\mu})=\sum\limits_{l=0}^{N}\sum\limits_{m=0}^{l}c_{i\>\\!j}^{(l-m,m)}(x^{\mu})\,\lambda^{l-m}\,k^{m}\,,$
(3)

where the functions $c_{i\>\\!j}^{(l,m)}$ in each order, labeled with the
indices $(l,m)$, depend on the coordinates only. Then, the metric is
differentiable in parameters $(\lambda,k)$, and it is possible to calculate
the two-dimensional Taylor expansion of the Einstein tensor around
$(\lambda,k)=(0,0)$. I denote the two-dimensional Taylor series of any
function $f(\lambda,k)$, truncated at the $N$-th order by:

$\mathcal{T}_{N}[f(\lambda,k)]=\sum\limits_{n=0}^{N}\frac{1}{n!}\sum\limits_{m=0}^{n}\binom{n}{m}\frac{\partial^{n}f}{\partial\lambda^{n-m}\partial
k^{m}}\Bigg{|}_{\begin{subarray}{c}\lambda=0\\\
k=0\end{subarray}}\lambda^{n-m}\,k^{m}\,.$ (4)

The Einstein equations hold approximately:

$\mathcal{T}_{N}[G^{\mu}\>\\!{}_{\nu}(\lambda,k)]=8\pi\,\rho\,U^{\mu}\,U_{\nu}-\Lambda\,\delta^{\mu}_{\nu}+8\pi\Delta
T^{\mu}\>\\!{}_{\nu}\,,$ (5)

since the truncation of the Taylor series of Einstein tensor elements
introduces some uncertainty $\Delta T^{\mu}\>\\!{}_{\nu}$. We demand that the
four-velocity $(U^{\mu})=(1,0,0,0)$ is unperturbed; the density
$\rho(t,\vec{x})$ is a known specified function and the uncertainty $|\Delta
T^{\mu}\>\\!{}_{\nu}|\ll 1$ is small. For exact solutions to the Einstein
equations, one would expect $G^{\mu}\>\\!{}_{\nu}$ instead of
$\mathcal{T}_{N}[G^{\mu}\>\\!{}_{\nu}]$ on the left-hand side of (5) and
$\Delta T^{\mu}\>\\!{}_{\nu}=0$. However, in the real Universe the pressure of
the intra-cluster gas and the proper motions of the galaxy cluster members
introduce departures from the energy-momentum tensor of the dust. In the
previous paper 2021EPJC…81..208S , we estimated the upper bound of the
uncertainty $|\Delta T^{\mu}\>\\!{}_{\nu}|\approx 10^{-6}\,\rho$ consistent
with a typical value of the velocity dispersion of galaxies observed in galaxy
clusters. In the matter-dominated era, it is then justified to describe the
matter content of the universe with the approximate energy-momentum tensor
$T^{\mu}\>\\!{}_{\nu}=\rho\,U^{\mu}\,U_{\nu}+\Delta T^{\mu}\>\\!{}_{\nu}$
satisfying this empirical constraint.

###  2.3 The background.

The second simplification concerns the background. The form of the metric
functions (3) is not compatible with the limit (2). For that reason, we
consider a weaker condition:

$\lim\limits_{\lambda\rightarrow
0}c_{i\>\\!j}=\mathcal{T}_{N}[\widetilde{c}_{i\>\\!j}]\,.$ (6)

Because the parameter $k$ is small, the difference between
$\widetilde{c}_{i\>\\!j}$ and $\mathcal{T}_{N}[\widetilde{c}_{i\>\\!j}]$
remains small either, if only the coordinate distance
$r=\sqrt{x^{2}+y^{2}+z^{2}}$ is not too large. For that reason, the model
presented in this paper can only describe the local neighborhood of the chosen
observer. The form of the metric functions $c_{i\>\\!j}^{(0,m)}$ follows from
the condition (6). The first few of these functions are:

$\displaystyle c_{i\>\\!j}^{(0,0)}=\delta_{i\>\\!j},\quad
c_{i\>\\!j}^{(0,1)}=-\frac{1}{2}r^{2}\,\delta_{i\>\\!j},\quad
c_{i\>\\!j}^{(0,2)}=\frac{3}{16}r^{4}\,\delta_{i\>\\!j}$ (7) $\displaystyle
c_{i\>\\!j}^{(0,3)}=-\frac{1}{16}r^{6}\,\delta_{i\>\\!j},\quad
c_{i\>\\!j}^{(0,4)}=\frac{5}{256}r^{8}\,\delta_{i\>\\!j}\,.$

When $\lambda=0$, the metric functions other than $c_{i\>\\!j}^{(0,m)}$ are
equal to zero, and the Einstein tensor becomes a function of the curvature
parameter only $G^{\mu}\>\\!{}_{\nu}(0,k)$. One can check by inspection that:

$\displaystyle\mathcal{T}_{4}[G^{0}\>\\!{}_{0}(0,k)]=-\frac{3(\dot{a}^{2}+k)}{a^{2}}\,,$
(8)
$\displaystyle\mathcal{T}_{4}[G^{i}\>\\!{}_{j}(0,k)]\Big{|}_{i=j}=-\frac{2a\ddot{a}+\dot{a}^{2}+k}{a^{2}}\,,$

where the overdot denotes differentiation with respect to time $t$. After
inserting these formulas into (5), one can see that the scale factor $a(t)$
satisfies the usual Friedmann equations for a universe filled with dust, with
an arbitrary curvature parameter $k$ and the positive cosmological constant
$\Lambda$. It is convenient to introduce the Lemaître-Hubble parameter
$H(t)=\dot{a}(t)/a(t)$, the Lemaître-Hubble constant $H_{0}=H(t_{0})$ in which
$t_{0}$ is the age of the universe, and the cosmological parameters
$\Omega_{m}=8\pi\rho(t_{0})/(3H_{0}^{2})$ and
$\Omega_{\Lambda}=\Lambda/(3H_{0}^{2})$. Then, the Friedmann equations reduce
to one ordinary differential equation for the scale factor:

$\frac{\mathrm{d}a}{\mathrm{d}t}=H_{0}\,\sqrt{\frac{\Omega_{m}}{a(t)}+\Omega_{\Lambda}\,a(t)^{2}-\frac{k}{H_{0}^{2}}}\,,$
(9)

where the curvature parameter is
$k=H_{0}^{2}\,\left(\Omega_{m}+\Omega_{\Lambda}-1\right)$. Solutions to an
equation of this type are known in terms of the Weierstrass elliptic
$\wp$-function (e. g. 1987AJ…..94.1373D ). However, for the purpose of this
article, we solve this equation backward in time numerically, with the help of
the fourth-order Runge-Kutta method and the initial condition $a(t_{0})=1$.

It is worth mentioning that the limit $\lim\limits_{\lambda\rightarrow
0}g_{\mu\>\\!\nu}$ provides the FLRW model as a background spacetime for which
the amplitude of inhomogeneities vanishes while the spatially flat background
$\lim\limits_{(\lambda,k)\rightarrow(0,0)}g_{\mu\>\\!\nu}$ is auxiliary and
has no physical meaning. Since $a(t)$ satisfies Friedmann’s equations with an
arbitrarily $k$ and a positive cosmological constant, the quantities
$3\dot{a}^{2}/a^{2}$ and $-(2a\ddot{a}+\dot{a}^{2})/a^{2}$ do not represent
the density and the pressure of any physical fluid.

The main task of this article is to present a method to find the remaining
metric functions $c_{i\>\\!j}^{(l,m)}$ for $l\geq 1$, so that the approximate
Einstein equations (5) up to second order hold within an acceptable
uncertainty $|\Delta T^{\mu}\>\\!{}_{\nu}|$.

##  3 Perturbations in the linear order.

When $N=1$, the summation (3) consists of elements $c_{i\>\\!j}^{(1,0)}$ and
$c_{i\>\\!j}^{(0,1)}$. The second one is given by (7), so it remains only
$c_{i\>\\!j}^{(1,0)}$. Let us introduce only scalar perturbations in the
linear order. We restrict the metric elements to the following form:

$c_{i\>\\!j}^{(1,0)}=\mathscr{A}_{10}(t)\frac{\partial^{2}}{\partial
x^{i}\partial x^{j}}A_{10}(x,y,z)+B_{10}(x,y,z)\,\delta_{i\>\\!j}\,.$ (10)

To complete the construction of the model in the first order, one should
determine the two functions of the spatial variables $A_{10}$, $B_{10}$, and
the one function of time $\mathscr{A}_{10}$. For convenience, we denote the
contributions to the Einstein tensor components from the consecutive orders in
a way similar to (3):

$\mathcal{T}_{N}[G^{\mu}\>\\!{}_{\nu}(\lambda,k)]=\sum\limits_{l=0}^{N}\sum\limits_{m=0}^{l}[G^{\mu}\>\\!{}_{\nu}]^{(l-m,m)}\,\lambda^{l-m}\,k^{m}\,.$
(11)

The terms $[G^{\mu}\>\\!{}_{\nu}]^{(0,0)}$ and
$[G^{\mu}\>\\!{}_{\nu}]^{(0,1)}$ are included in (8). In the following, we
analyze the remaining part $[G^{\mu}\>\\!{}_{\nu}]^{(1,0)}$.

The assumed form of the metric (10) guarantees that the elements
$[G^{i}\>\\!{}_{0}]^{(1,0)}$ are equal to zero. The other two types of
Einstein tensor elements $[G^{i}\>\\!{}_{j}]^{(1,0)}\big{|}_{i=j}$ and
$[G^{i}\>\\!{}_{j}]^{(1,0)}\big{|}_{i\neq j}$ contain non-zero terms with
different time dependence. However, one can unify the time dependence of all
the terms if only the following differential equation holds:

$a^{2}(t)\,\ddot{\mathscr{A}}_{10}(t)+3\,a(t)\,\dot{a}(t)\,\dot{\mathscr{A}}_{10}(t)=\alpha\,,$
(12)

where $\alpha$ is a constant. Since the scale factor $a(t)$ is known, this
becomes a well-posed equation for the metric function $\mathscr{A}_{10}(t)$.
With the auxiliary function $f_{1}$, such that
$f_{1}(t)=\dot{\mathscr{A}}_{10}(t)$, this equation reduces to the
inhomogeneous first-order linear ODE, which can be solved with the variation
of constant method. Specifying the initial conditions at some initial time
$t_{i}$, one can write the solution as:

$\displaystyle\mathscr{A}_{10}(t)=\int\limits_{t_{i}}^{t}f_{1}(t^{\prime})\,\mathrm{d}t^{\prime}+\mathcal{C}_{1},\quad
P_{1}(t)=3\int\limits_{t_{i}}^{t}\frac{\dot{a}(t^{\prime})}{a(t^{\prime})}\,\mathrm{d}t^{\prime}$
(13)
$\displaystyle\mathrm{and}\>\>f_{1}(t)=\left(\mathcal{C}_{2}+\alpha\int\limits_{t_{i}}^{t}\frac{\mathrm{e}^{P_{1}(t^{\prime})}}{a^{2}(t^{\prime})}\,\mathrm{d}t^{\prime}\right)\,\mathrm{e}^{-P_{1}(t)}\,.$

For simplicity, we may ignore the integration constants $\mathcal{C}_{1}=0$
and $\mathcal{C}_{2}=0$. This means that
$\mathscr{A}_{10}(t_{i})=f_{1}(t_{i})=0$. If $t_{i}$ is the time of the Big
Bang $t_{i}=0$, then this condition is justifiable.

In the previous paper 2021EPJC…81..208S , we investigated perturbations around
the Einstein-de Sitter (EdS) background, with $\Omega_{m}=1$,
$\Omega_{\Lambda}=0$ and $k=0$. The EdS scale factor is
$a^{(EdS)}(t)=\mathcal{C}\,t^{2/3}$. The value of the constant
$\mathcal{C}=4.85\times 10^{-3}$ can be found using the megaparsec as a unit
of length, considering normalization $a^{(EdS)}(t_{0})=1$ and calculating the
age of the EdS universe $t_{0}$ for the Lemaître-Hubble constant value
$H_{0}=67.37\,\mathrm{[km/s/Mpc]}$ taken from the Planck data 2020A&A…641A…6P
. In the current notation, the respective metric function from the work
2021EPJC…81..208S is $\mathscr{A}_{10}(t)=t^{2/3}$. It is easy to check that
this form of $\mathscr{A}_{10}(t)$ follows from (13), when
$\alpha=\frac{10}{9}\mathcal{C}^{2}=2.61\times 10^{-5}$. For generic values of
the background parameters $(\Omega_{m},\Omega_{\Lambda})$, the scale factor
could no longer have a simple power-law dependence. However, taking the
initial time slightly above the Big Bang $t_{i}=10^{-3}$, the integrals (13)
can be easily evaluated numerically. Figure 1 shows the resulting metric
function $\mathscr{A}_{10}(t)$ for various cosmological parameters
$\Omega_{m}$, $\Omega_{\Lambda}$, and fixed $\alpha=2.61\times 10^{-5}$.

Figure 1:  The metric function $\mathscr{A}_{10}(t)$ calculated for various
$\Omega_{m}$, $\Omega_{\Lambda}$ parameters. The curves are plotted in the
range $t\in(t_{i},t_{0}(\Omega_{m},\Omega_{\Lambda}))$, where $t_{i}$ is the
initial time described in the text and $t_{0}(\Omega_{m},\Omega_{\Lambda})$ is
the universe age of each cosmological model evaluated for
$H_{0}=67.37\,\mathrm{[km/s/Mpc]}$.

When the temporal metric function $\mathscr{A}_{10}(t)$ satisfies (13), the
Einstein tensor elements $[G^{i}\>\\!{}_{j}]^{(1,0)}\big{|}_{i=j}$ and
$[G^{i}\>\\!{}_{j}]^{(1,0)}\big{|}_{i\neq j}$ are equal to zero if only the
additional condition on the spatial metric functions is satisfied:

$B_{10}(x,y,z)=\alpha\,A_{10}(x,y,z)\,.$ (14)

In this way, one can construct a specific dust solution to the linear
perturbation theory.

The energy density $\rho^{(1,0)}=-[G^{0}\>\\!{}_{0}]^{(1,0)}/(8\pi)$ becomes:

$\rho^{(1,0)}(t,x,y,z)=\frac{a(t)\,\dot{a}(t)\,\dot{\mathscr{A}}_{10}(t)-\alpha}{8\pi\,a^{2}(t)}\,\triangle
A_{10}(x,y,z)\,.$ (15)

The symbol $\triangle$ denotes the Laplace operator in Cartesian coordinates.
The specification of the cosmological parameters $\Omega_{m}$ and
$\Omega_{\Lambda}$ constrains the time evolution of the first-order density,
while its spatial distribution depends on the arbitrary function $A_{10}$.

Since the functions $\mathscr{A}_{10}$ and $B_{10}$ are proportional to
$\alpha$, the metric elements $c_{i\>\\!j}^{(1,0)}$ are proportional to
$\alpha$. From (15) it follows that the first-order density $\rho^{(1,0)}$ is
also proportional to $\alpha$. This means that the value of the constant
$\alpha$ does not matter. The value of $\alpha$ can be arbitrarily changed,
and then the amplitude $\lambda$ should be redefined so that the first-order
metric $\lambda\,c_{i\>\\!j}^{(1,0)}$ and the first-order density
$\lambda\rho^{(1,0)}$ remain unaffected. For numerical reasons, it is
convenient to keep the constant $\alpha$, because it enables one to normalize
the function $\mathscr{A}_{10}$.

Figure 2:  The time dependence of the first-order contribution to the density.
Some specific values of the cosmological parameters $\Omega_{m}$ and
$\Omega_{\Lambda}$ are considered.

Figure 2 shows the time dependence of the first-order density
$\rho^{(1,0)}(t)=(a(t)\,\dot{a}(t)\,\dot{\mathscr{A}}_{10}(t)-\alpha)/(8\pi\,a^{2}(t))$.
For every value of the background parameters $\Omega_{m}$ and
$\Omega_{\Lambda}$, the first-order density is a decreasing function of time.
The homogeneous density of the FLRW background model $\rho^{(0)}(t)$ is also a
decreasing function. If $\rho^{(0)}$ decreases faster than $\rho^{(1,0)}$, the
density contrast $\delta=\lambda\,\rho^{(1,0)}/\rho^{(0)}$ may increase with
time.

Figure 3:  The density contrast $\delta\rho(t)$ for different background FLRW
models and the fixed value of the amplitude $\lambda=9.68\times 10^{-4}$.
Figure 4:  The density contrast $\delta\rho(t)$ evaluated for different values
of the amplitude $\lambda$. For the set of the background parameters
$(\Omega_{m}\Omega_{\Lambda})\in\\{(1,0),(0.3,0.7),(0.3,0),(0.7,0.5),(0.05,0)\\}$
the respective values of the amplitude are $\lambda\in\\{9.68\times
10^{-4},3.71\times 10^{-4},3.19\times 10^{-4},7.45\times 10^{-4},5.68\times
10^{-5}\\}$.

In Fig. 3, the density contrast is plotted for different cosmological
parameters $\Omega_{m}$, $\Omega_{\Lambda}$, and the same value of the
amplitude $\lambda=9.68\times 10^{-4}$. In Fig. 4, we show the density
contrast for different values of $\lambda$. For each background spacetime, we
chose the respective value of $\lambda$ such that the density contrast at the
universe age is equal to 0.1. The black and green curves correspond to models
that have the same $\Omega_{m}$ but differ by the value of the cosmological
constant. One can observe that the presence of the cosmological constant
causes a smaller grouth-rate of the density contrast.

##  4 Perturbations in the second order.

For $N=2$, in the summation (3) appear the elements $c_{i\>\\!j}^{(2,0)}$,
$c_{i\>\\!j}^{(1,1)}$, and $c_{i\>\\!j}^{(0,2)}$. The latter is defined by
(7), so the elements $c_{i\>\\!j}^{(2,0)}$ and $c_{i\>\\!j}^{(1,1)}$ remain to
be determined. From (8) it follows that $[G^{\mu}\>\\!{}_{\nu}]^{(0,2)}=0$. To
achieve the dust solution in second-order perturbation theory, one should
analyze contributions to the Einstein tensor elements
$[G^{\mu}\>\\!{}_{\nu}]^{(2,0)}$ and $[G^{\mu}\>\\!{}_{\nu}]^{(1,1)}$.
Equations in the second order become quite tough. To simplify them, let us
introduce the following restriction.

As explained in the previous section, the spatial distribution of the density
in the first order is given by the arbitrary function $A_{10}$. Let us assume
that this function is separable:

$A_{10}(x,y,z)=C_{10}(x)+C_{10}(y)+C_{10}(z)\,,$ (16)

where $C_{10}$ is an arbitrary function of one variable. To fulfill (14) one
should also take:

$B_{10}(x,y,z)=D_{10}(x)+D_{10}(y)+D_{10}(z)\,,$ (17)

with $D_{10}(w)=\alpha\,C_{10}(w)$ and $w=x,y,z$.

Figure 5:  The isodensity surfaces of the first-order density $\rho^{(1,0)}$,
for which the function $C_{10}$ is given by (18), with parameters $M=1$,
$P_{1}=-625/\pi^{2}$, $Q_{1}=\pi/25$ and $R_{1}=0$.

If one gives the arbitrary function $C_{10}$ as a linear combination of the
sine functions:

$C_{10}(w)=\sum\limits_{s=1}^{M}P_{s}\,\sin(Q_{s}\,w+R_{s})\,,$ (18)

then the density $\rho^{(1,0)}$ forms an infinite cubic lattice. We plot an
exemplary density distribution of this type in Figure 5 for one mode only.
Although condition (16) constrains the space of all possible density
distributions, quite complicated distributions are still admissible. The
subsequent modes present in (18) could describe different substructures
distributed over the elementary cell. Moreover, such a model is locally
inhomogeneous, but for scales much larger than the size of the elementary
cell, it becomes homogeneous and isotropic in common sense.

The simplification (16) eliminates the mixed derivatives that appear in the
Einstein tensor elements. After that, elements $[G^{i}\>\\!{}_{0}]^{(2,0)}$
are proportional to $\alpha\,\dot{\mathscr{A}_{10}}(t)$, elements
$[G^{i}\>\\!{}_{j}]^{(2,0)}\big{|}_{i\neq j}$ are proportional to
$\alpha^{2}\,a(t)^{-2}$, while diagonal elements
$[G^{i}\>\\!{}_{j}]^{(2,0)}\big{|}_{i=j}$ comprise three types of terms. The
terms of the first type are proportional to
$\alpha\,\mathscr{A}_{10}(t)\,a(t)^{-2}$, the terms of the second type are
proportional to $\dot{\mathscr{A}_{10}}(t)^{2}$, and the terms of the third
type are proportional to $\alpha^{2}\,a(t)^{-2}$.

Figure 6:  The time dependence of various components of the Einstein tensor
elements, calculated for $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$ and
$\alpha=2.61\times 10^{-5}$. For comparison, the time-dependent metric
functions $\mathscr{D}_{20}(t)$, $\mathscr{F}_{20}(t)$, introduced further in
the text, and their derivatives.

In Fig. 6 these time-dependent functions are plotted for the scale factor
$a(t)$ and the function $\mathscr{A}_{10}(t)$ calculated for the parameters
$\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$ and the already used value
$\alpha=2.61\times 10^{-5}$. For other values of the cosmological parameters
$\Omega_{m}$ and $\Omega_{\Lambda}$ considered in this paper, this figure
looks the same. One can observe that for every time $t$ between the Big Bang
and the age of the universe, the functions $\alpha\,\dot{\mathscr{A}_{10}}(t)$
and $\alpha^{2}\,a(t)^{-2}$ are at least five orders of magnitude smaller than
$\alpha\,\mathscr{A}_{10}(t)\,a(t)^{-2}$ and $\dot{\mathscr{A}_{10}}(t)^{2}$.
The energy-momentum tensor elements consist of the time-dependent functions
mentioned above multiplied by some terms that contain spatial metric functions
and their derivatives. If the frequency of the sine functions in (18) is not
too high, then the derivatives of the spatial metric functions are relatively
small, and all the $[G^{\mu}\>\\!{}_{\nu}]^{(2,0)}$ terms proportional to
$\alpha\,\dot{\mathscr{A}_{10}}(t)$ and $\alpha^{2}\,a(t)^{-2}$ could be
neglected.

The remaining two types of terms present in
$[G^{i}\>\\!{}_{j}]^{(2,0)}\big{|}_{i=j}$ can be eliminated if the metric
functions in the order $(2,0)$ are the following:

$\displaystyle
c_{i\>\\!j}^{(2,0)}=\left(\mathscr{A}_{20}(t)C_{20}(x^{i})+\mathscr{D}_{20}(t)\sum\limits_{s=1}^{3}D_{20}(x^{s})\right)\delta_{i\>\\!j}+$
(19)
$\displaystyle+\mathscr{F}_{20}(t)\,F_{20}(x^{i},x^{j})\,(1-\delta_{i\>\\!j})\,.$

The two of the time-dependent functions are given explicitly:

$\mathscr{D}_{20}(t)=a^{2}(t)\,\dot{\mathscr{A}}_{10}(t)^{2}\,,$ (20)
$\mathscr{F}_{20}(t)=\frac{1}{4}a^{2}(t)\,\dot{\mathscr{A}}_{10}(t)^{2}+\alpha\,\mathscr{A}_{10}(t)\,.$
(21)

The third one should satisfy the differential equation:

$a^{2}(t)\,\ddot{\mathscr{A}}_{20}(t)+3\,a(t)\,\dot{a}(t)\,\dot{\mathscr{A}}_{20}(t)=2\,\alpha\,\mathscr{A}_{10}(t)\,.$
(22)

This equation is similar to (12). Since the function $\mathscr{A}_{10}$ is
known from the first order, it can be integrated numerically with the help of
the same method.

The spatial metric functions should satisfy the following equations:

$D^{\prime\prime}_{20}(w)=-\frac{1}{2}\,C^{\prime\prime}_{10}(w)^{2}\,,$ (23)
$\frac{\partial^{2}}{\partial w\partial
v}F_{20}(v,w)=-C^{\prime\prime}_{10}(v)\,C^{\prime\prime}_{10}(w)\,,$ (24)
$C_{20}(w)=-\frac{1}{4}\,C^{\prime}_{10}(w)\,C^{\prime\prime\prime}_{10}(w)\,,$
(25)

where the prime denotes the derivative with respect to one of the spatial
variables $u,v\in\\{x,y,z\\}$. For the function $C_{10}$ given by (18), it is
easy to find the explicit solution to these equations:

$\displaystyle
D_{20}(w)\\!=\\!-\frac{1}{8}\sum\limits_{s=1}^{M}P_{s}^{2}Q_{s}^{4}\\!\left(w^{2}+\frac{\cos(2Q_{s}w\\!+\\!2R_{s})}{2Q_{s}^{2}}\right)+$
(26)
$\displaystyle+\frac{1}{4}\sum\limits_{s=1}^{M}\sum\limits_{\begin{subarray}{c}r=1\\\
r\neq
s\end{subarray}}^{M}P_{s}P_{r}Q_{s}^{2}Q_{r}^{2}\\!\left[\frac{\cos((Q_{s}\\!-\\!Q_{r})w+R_{s}\\!-\\!R_{r})}{(Q_{s}\\!-\\!Q_{r})^{2}}-\right.$
$\displaystyle\left.-\frac{\cos((Q_{s}\\!+\\!Q_{r})w+R_{s}\\!+\\!R_{r})}{(Q_{s}\\!+\\!Q_{r})^{2}}\right]\,,$
$F_{20}(v,w)\\!=\\!-\\!\sum\limits_{s=1}^{M}\sum\limits_{r=1}^{M}\\!\\!P_{s}P_{r}Q_{s}Q_{r}\\!\cos(Q_{s}v+R_{s})\cos(Q_{r}w+R_{r})\,,$
(27)
$C_{20}(w)\\!=\\!\frac{1}{4}\sum\limits_{s=1}^{M}\sum\limits_{r=1}^{M}P_{s}P_{r}Q_{s}Q_{r}^{3}\cos(Q_{s}v+R_{s})\cos(Q_{r}w+R_{r})\,.$
(28)

Integration constants are ignored for simplicity.

For the metric functions described above, the leading-order terms proportional
to $\alpha\,\mathscr{A}_{10}(t)\,a(t)^{-2}$ and
$\dot{\mathscr{A}_{10}}(t)^{2}$ are eliminated from
$[G^{i}\>\\!{}_{j}]^{(2,0)}\big{|}_{i=j}$. Instead, some terms proportional to
the derivatives of the metric functions $\dot{\mathscr{D}}_{20}$,
$\dot{\mathscr{F}}_{20}$, and $\ddot{\mathscr{F}}_{20}$ appear in the other
elements of $[G^{\mu}\>\\!{}_{\nu}]^{(2,0)}$. The function
$\dot{\mathscr{D}}_{20}$ changes sign, but as could be seen in Fig. 6, its
absolute value $|\dot{\mathscr{D}}_{20}|$ remains very small compared to the
leading-order terms. Similarly, the derivatives $\dot{\mathscr{F}}_{20}$ and
$\ddot{\mathscr{F}}_{20}$ are also small and could be neglected. In effect,
all the Einstein tensor elements $[G^{\mu}\>\\!{}_{\nu}]^{(2,0)}\approx 0$,
apart from the second-order contribution to the energy density
$\rho^{(2,0)}=-[G^{0}\>\\!{}_{0}]^{(2,0)}/(8\pi)$. With the help of equations
(20-25) the energy density could be simplified to the following form:

$\rho^{(2,0)}(t,x,y,z)=H(t,x)+H(t,y)+H(t,z)\,,$ (29)

where

$\displaystyle
H(t,w)=\left[\frac{\alpha\mathscr{A}_{10}(t)}{2a^{2}(t)}-\frac{\dot{\mathscr{A}}_{20}(t)\dot{a}(t)}{4a(t)}\right]C^{\prime}_{10}(w)C^{\prime\prime\prime}_{10}(w)+$
(30)
$\displaystyle+\left[\frac{\dot{\mathscr{A}}_{10}(t)^{2}}{2}+\frac{\alpha\mathscr{A}_{10}(t)}{a^{2}(t)}-\frac{\mathscr{A}_{10}(t)\dot{\mathscr{A}_{10}(t)}\dot{a}(t)}{a(t)}\right]C^{\prime\prime}_{10}(w)^{2}-$
$\displaystyle-\frac{\alpha\dot{\mathscr{A}}_{10}(t)\dot{a}(t)}{a(t)}C_{10}(w)C^{\prime\prime}_{10}(w)+\frac{3\alpha\dot{a}(t)}{a(t)}\frac{\mathrm{d}}{\mathrm{d}t}\left[a^{2}\dot{\mathscr{A}}_{10}^{2}\right]C_{10}(w)$

Each term of $H(t,w)$ represents a decreasing function of time. The last two
are the smallest.

Similarly, for the order $(1,1)$, one can choose the same proposition for the
structure of the metric functions:

$\displaystyle
c_{i\>\\!j}^{(1,1)}=\left(\mathscr{A}_{11}(t)C_{11}(x^{i})+\mathscr{D}_{11}(t)\sum\limits_{s=1}^{3}D_{11}(x^{s})\right)\delta_{i\>\\!j}+$
(31)
$\displaystyle+\mathscr{F}_{11}(t)\,F_{11}(x^{i},x^{j})\,(1-\delta_{i\>\\!j})\,.$

The Einstein tensor elements $[G^{i}\>\\!{}_{j}]^{(1,1)}\big{|}_{i=j}=0$ and
$[G^{i}\>\\!{}_{0}]^{(1,1)}=0$ vanish if the metric functions satisfy the
following relations. The two of the time-dependent functions are equal to the
first-order function $\mathscr{A}_{10}$:

$\mathscr{D}_{11}(t)=\mathscr{F}_{11}(t)=\mathscr{A}_{10}(t)\,,$ (32)

while the third one should satisfy the differential equation:

$a^{2}(t)\,\ddot{\mathscr{A}}_{11}(t)+3\,a(t)\,\dot{a}(t)\,\dot{\mathscr{A}}_{11}(t)=\mathscr{A}_{10}(t)\,.$
(33)

Again, this equation is of the same type as (12) and (22); therefore, it can
be integrated numerically in a similar manner. The spatial metric functions
are given by:

$D_{11}(w)=C_{10}(w)-\frac{1}{2}\,w\,C^{\prime}_{10}(w)\,,$ (34)
$F_{11}(v,w)=\frac{1}{2}\left(v\,C^{\prime}_{10}(w)+w\,C^{\prime}_{10}(v)\right)\,,$
(35) $C_{11}(w)=C^{\prime\prime}_{10}(w)\,.$ (36)

After that, the only non-zero components of the Einstein tensor in order
$(1,1)$ are $[G^{0}\>\\!{}_{0}]^{(1,1)}$ and
$[G^{i}\>\\!{}_{j}]^{(1,1)}\big{|}_{i\neq j}$. However, the components
$[G^{i}\>\\!{}_{j}]^{(1,1)}\big{|}_{i\neq j}$ are proportional to
$\alpha/a^{2}(t)$ and can be neglected. The density of order $(1,1)$
$\rho^{(1,1)}=-[G^{0}\>\\!{}_{0}]^{(1,1)}/(8\pi)$ can be simplified to the
following form:

$\rho^{(1,1)}(t,x,y,z)=K(t,x,r)+K(t,y,r)+K(t,z,r)\,,$ (37)

where

$\displaystyle
K(t,w,r)=\left[\frac{\dot{\mathscr{A}}_{10}(t)\dot{a}(t)}{2a(t)}\,r^{2}-\frac{\mathscr{A}_{10}(t)}{a^{2}(t)}\right]C^{\prime\prime}_{10}(w)+$
(38)
$\displaystyle+\ddot{\mathscr{A}}_{10}(t)\left[\frac{1}{2}w\,C^{\prime}_{10}(w)-C_{10}(w)\right]+\frac{\dot{\mathscr{A}}_{11}(t)\dot{a}(t)}{a(t)}C_{10}(w)\,.$

Once again, this point shows that the model presented here is adequate only
for the description of the local neighborhood of the selected observer, where
the coordinate distance $r^{2}=x^{2}+y^{2}+z^{2}$ is small.

###  4.1 The convergence

To construct an approximate solution to the Einstein equations, we used
perturbative methods and neglected some terms identified as small compared to
the leading-order terms. It is then necessary to verify whether the solution
converges to a universe filled with dust. If the model converges properly,
then the uncertainty of the energy-momentum tensor $\Delta
T^{\mu}\>\\!{}_{\nu}$ should decrease when the contributions to the metric
tensor from the consecutive orders are taken into account.

In Fig. 7 we plot the components $\Delta T^{\mu}\>\\!{}_{\nu}$ in units of
energy density $\rho$ for the Taylor expansion of the Einstein tensor
$\mathcal{T}_{N}[G^{\mu}\>\\!{}_{\nu}(\lambda,k)]$ evaluated up to the third
order $N=3$. We take the arbitrary function $A_{10}$ as the one mode sine
function for which the density distribution is plotted in Fig. 5. The value of
the amplitude $\lambda$ is such that the first-order density $\rho^{(1,0)}$
evaluated at the maximum of the overdensity equals $0.025$ in the critical
density units. In each plot, we show the maximum of the absolute value
$|\,\Delta T^{\mu}\>\\!{}_{\nu}\,|$ among the values at one hundred randomly
chosen points within the elementary cell. The blue curves represent pressure-
like terms $\Delta T^{i}\>\\!{}_{j}\big{|}_{i=j}$, the orange curves
correspond to $\Delta T^{i}\>\\!{}_{0}$, and the green curves describe $\Delta
T^{i}\>\\!{}_{j}\big{|}_{i\neq j}$. The results for the metric truncated at
linear order are plotted by dashed lines. When the metric functions are
considered up to the second order, the corresponding results are plotted by
solid lines.

In each case, pressure-like terms give the greatest contribution to the
discrepancy $\Delta T^{\mu}\>\\!{}_{\nu}$ between the energy-momentum tensor
of the model and the energy-momentum tensor of dust. When the second-order
metric is taken into account, the absolute value of the pressure-like terms is
smaller than in the case of linear-order perturbations. In the case of the
spatially flat universe, with $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, and
$k=0$, the second-order results are improved by more than two orders of
magnitude compared to the first-order predictions. The value of the pressure-
like terms $\Delta T^{i}\>\\!{}_{j}\big{|}_{i=j}/\rho$ below $10^{-5}$ is
acceptable, since the pressure generated by the proper motions of members of a
typical galaxy cluster is of the order of $10^{-6}\,\rho$. For a non-flat
background, the value of the pressure-like terms is higher. We observe the
values of $\Delta T^{i}\>\\!{}_{j}\big{|}_{i=j}/\rho$ around $10^{-4}$ and
$10^{-3}$ for the parameters $(\Omega_{m},\Omega_{\Lambda})=(0.3,0.6)$ and
$(\Omega_{m},\Omega_{\Lambda})=(0.3,0.0)$, respectively. In these cases, to
achieve better accuracy, third-order corrections to the metric tensor are
needed.

In each plot, the terms $\Delta T^{i}\>\\!{}_{0}$ and $\Delta
T^{i}\>\\!{}_{j}\big{|}_{i\neq j}$ are improved when the second-order metric
is taken into account or their values do not change too much. Generally, we
may conclude that the convergence of the presented model is satisfactory.

Figure 7:  The maximum of $|\,\Delta T^{\mu}\>\\!{}_{\nu}\,(t)|/\rho(t)$ over
the elementary cell. The _top_ graph corresponds to
$(\Omega_{m},\Omega_{\Lambda})=(0.3,0.7)$, the _middle_ is done for
$(\Omega_{m},\Omega_{\Lambda})=(0.3,0.6)$, and the _bottom_ \-
$(\Omega_{m},\Omega_{\Lambda})=(0.3,0.0)$. The _blue_ curves describe $\Delta
T^{i}\>\\!{}_{j}\big{|}_{i=j}$, the _orange_ curves - $\Delta
T^{i}\>\\!{}_{0}$, and the _green_ curves - $\Delta
T^{i}\>\\!{}_{j}\big{|}_{i\neq j}$. _Dashed_ lines are generated for the
linear order metric, while the _solid_ lines correspond to the metric up to
the second order.

###  4.2 The gauge

In Fig. 8 we show the expansion scalar $\theta=\nabla_{\mu}\,U^{\mu}$ for an
exemplary model. We take the background parameters
$(\Omega_{m},\Omega_{\Lambda})=(0.3,0.7)$, fix the arbitrary function $A_{10}$
as the one mode sine function defining the density distribution presented in
Fig. 5, and take the value of the amplitude $\lambda$ so that the first-order
density $\rho^{(1,0)}$ at the maximum of the overdensity is $0.025$ in the
critical density units. Then, we randomly generate one hundred points within
the elementary cell. The time dependence of $\theta$ at the points whose
density is higher than the background density $\Omega_{m}=0.3$ is plotted in
Fig. 8 by the blue curves. The relation $\theta(t)$ at the points where the
model density is lower than $0.3$ is plotted by the light blue curves. By the
red dashed curve, we show the time dependence of the expansion scalar of the
background. It is seen that in the considered inhomogeneous model the space is
expanding everywhere, and the rate of expansion is slightly changed compared
to the background. The overdense regions expand a bit slower than the
underdense areas. Because the bundle of time-like geodesics is not collapsing,
the world-lines of initially separated observers would not cross with each
other, and the Gaussian normal coordinates (1) cover the whole space-time.

Figure 8:  The expansion scalar $\theta$ as a function of time for an
exemplary model on the background $(\Omega_{m},\Omega_{\Lambda})=(0.3,0.7)$.
The expansion scalar evaluated in the overdensity regions is plotted by the
_blue_ curves, while $\theta$ in the underdensity regions are drawn by the
_light blue_ curves. The _dashed red_ curve represents the expansion scalar in
the background space-time. The axes range in the _bottom_ panel helps us to
visualize more details.

The perturbation theory is carried out within the synchronous comoving gauge
when the form of the metric tensor (1) is fixed, so that $g_{0\>\\!0}=-1$,
$g_{0\>\\!i}=0$, and the velocity field $U^{\mu}$ is unperturbed. These
conditions do not completely determine the gauge. There is the following gauge
freedom left:

$x^{i}\mapsto\chi^{i}(x^{j})\,,$ (39)

where $\chi^{i}$, for $i=1,2,3$, are three arbitrary functions of the spatial
variables only. The gauge transformation (39) leaves the form of the metric
tensor (1) unchanged and the velocity field unperturbed. It is necessary to
check whether the considered perturbations are not fictitious. The phenomena
which are not physical can be removed by a gauge transformation.

Let us specify the following transformation:

$\chi^{i}(x^{j})=x^{i}+\lambda\,[\chi^{(1,0)}]^{i}(x^{j})+k\,\,[\chi^{(0,1)}]^{i}(x^{j})\,,$
(40)

where $\chi^{(1,0)}$ and $\chi^{(0,1)}$ can be chosen arbitrarily. The model
density is transformed under this gauge transformation in the following way:

$\rho(t,\chi^{i}(x^{j}))\approx\rho(t,x^{j})+\lambda\,\frac{\partial\rho}{\partial
x^{m}}\,[\chi^{(1,0)}]^{m}(x^{j})+\\\ +k\,\frac{\partial\rho}{\partial
x^{m}}\,[\chi^{(0,1)}]^{m}(x^{j})\,.$ (41)

If the derivative $\partial\rho/\partial x^{j}$ is non-zero, then spurious
density fluctuations can be introduced in each order of perturbation theory.
The background spacetime is homogeneous, so its density $\rho_{0}$ does not
depend on the spatial coordinates $x^{j}$. Then it is not possible to generate
the first-order density perturbation $\rho^{(1,0)}$ by a gauge transformation.
The first-order density is not fictitious.

Suppose that we have the following density:

$\rho(t,x^{j})=\rho_{0}(t)+\lambda\,\rho^{(1,0)}(t,x^{j})\,.$ (42)

Since it depends on the spatial variables, the gauge transformation (41)
introduces non-physical density in the second order:

$\widetilde{\rho}^{(2,0)}(t,x^{j})=\frac{\partial\rho^{(1,0)}}{\partial
x^{m}}\,[\chi^{(1,0)}]^{m}(x^{j})\,,$ (43)

and

$\widetilde{\rho}^{(1,1)}(t,x^{j})=\frac{\partial\rho^{(1,0)}}{\partial
x^{m}}\,[\chi^{(0,1)}]^{m}(x^{j})\,.$ (44)

However, since the gauge functions $\chi^{i}$ depend only on spatial
variables, the time dependence of $\widetilde{\rho}^{(2,0)}$ and
$\widetilde{\rho}^{(1,1)}$ is inherited from the first-order density
$\rho^{(1,0)}$. The time dependence of $\widetilde{\rho}^{(2,0)}$ differs from
the time dependence of every term in (30). As a result, the second-order
density $\rho^{(2,0)}$ is also physical and cannot be produced from
$\rho^{(1,0)}$ by a gauge transformation. The same reasoning shows that the
density $\rho^{(1,1)}$, given by the formula (37), is also physical. In this
way, analysis of the time dependence of the density distribution enables one
to remove the gauge freedom.

##  5 Conclusions.

In this work we present a specific solution to the second-order cosmological
perturbation theory around a generic FLRW model filled with dust. The non-zero
spatial curvature of the background and a positive cosmological constant are
allowed. We simplified the metric tensor so that only one arbitrary function
$C_{10}$ ramains and defines the spatial distribution of the density, while
the other metric functions are determined from the model construction. If
$C_{10}$ is a periodic function, for example a linear combination of the sine
functions, then the density forms an infinite cubic lattice. The distribution
of the substructures within the elementary cell depends on the choice of
$C_{10}$. Although the density contrast could be an increasing function of
time, the space of the model universe is expanding at each point. For that
reason, the model does not intend to describe the formation of individual
structures but gives a tool for handling differences between large areas
occupied by the galaxy clusters and the cosmic voids. Since the metric tensor
is explicitly given and has a relatively simple form, the model can be used to
examine quantitatively how these differences affect cosmological observables.

## References

  * (1) Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, R. Battye, K. Benabed, J. P. Bernard, M. Bersanelli, P. Bielewicz, J. J. Bock, J. R. Bond, J. Borrill, F. R. Bouchet, F. Boulanger, M. Bucher, C. Burigana, R. C. Butler, E. Calabrese, J. F. Cardoso, J. Carron, A. Challinor, H. C. Chiang, J. Chluba, L. P. L. Colombo, C. Combet, D. Contreras, B. P. Crill, F. Cuttaia, P. de Bernardis, G. de Zotti, J. Delabrouille, J. M. Delouis, E. Di Valentino, J. M. Diego, O. Doré, M. Douspis, A. Ducout, X. Dupac, S. Dusini, G. Efstathiou, F. Elsner, T. A. Enßlin, H. K. Eriksen, Y. Fantaye, M. Farhang, J. Fergusson, R. Fernandez-Cobos, F. Finelli, F. Forastieri, M. Frailis, A. A. Fraisse, E. Franceschi, A. Frolov, S. Galeotta, S. Galli, K. Ganga, R. T. Génova-Santos, M. Gerbino, T. Ghosh, J. González-Nuevo, K. M. Górski, S. Gratton, A. Gruppuso, J. E. Gudmundsson, J. Hamann, W. Handley, F. K. Hansen, D. Herranz, S. R. Hildebrandt, E. Hivon, Z. Huang, A. H. Jaffe, W. C. Jones, A. Karakci, E. Keihänen, R. Keskitalo, K. Kiiveri, J. Kim, T. S. Kisner, L. Knox, N. Krachmalnicoff, M. Kunz, H. Kurki-Suonio, G. Lagache, J. M. Lamarre, A. Lasenby, M. Lattanzi, C. R. Lawrence, M. Le Jeune, P. Lemos, J. Lesgourgues, F. Levrier, A. Lewis, M. Liguori, P. B. Lilje, M. Lilley, V. Lindholm, M. López-Caniego, P. M. Lubin, Y. Z. Ma, J. F. Macías-Pérez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli, A. Marcos-Caballero, M. Maris, P. G. Martin, M. Martinelli, E. Martínez-González, S. Matarrese, N. Mauri, J. D. McEwen, P. R. Meinhold, A. Melchiorri, A. Mennella, M. Migliaccio, M. Millea, S. Mitra, M. A. Miville-Deschênes, D. Molinari, L. Montier, G. Morgante, A. Moss, P. Natoli, H. U. Nørgaard-Nielsen, L. Pagano, D. Paoletti, B. Partridge, G. Patanchon, H. V. Peiris, F. Perrotta, V. Pettorino, F. Piacentini, L. Polastri, G. Polenta, J. L. Puget, J. P. Rachen, M. Reinecke, M. Remazeilles, A. Renzi, G. Rocha, C. Rosset, G. Roudier, J. A. Rubiño-Martín, B. Ruiz-Granados, L. Salvati, M. Sandri, M. Savelainen, D. Scott, E. P. S. Shellard, C. Sirignano, G. Sirri, L. D. Spencer, R. Sunyaev, A. S. Suur-Uski, J. A. Tauber, D. Tavagnacco, M. Tenti, L. Toffolatti, M. Tomasi, T. Trombetti, L. Valenziano, J. Valiviita, B. Van Tent, L. Vibert, P. Vielva, F. Villa, N. Vittorio, B. D. Wandelt, I. K. Wehus, M. White, S. D. M. White, A. Zacchei, and A. Zonca.  Planck 2018 results. VI. Cosmological parameters.  Astronomy&Astrophysics, 641:A6, September 2020.
  * (2) Adam G. Riess, Stefano Casertano, Wenlong Yuan, Lucas M. Macri, and Dan Scolnic.  Large Magellanic Cloud Cepheid Standards Provide a 1% Foundation for the Determination of the Hubble Constant and Stronger Evidence for Physics beyond $\Lambda$CDM.  The Astrophysical Journal, 876(1):85, May 2019.
  * (3) Eleonora Di Valentino, Olga Mena, Supriya Pan, Luca Visinelli, Weiqiang Yang, Alessandro Melchiorri, David F. Mota, Adam G. Riess, and Joseph Silk.  In the realm of the Hubble tension-a review of solutions.  Classical and Quantum Gravity, 38(15):153001, July 2021.
  * (4) Nan Li and Dominik J. Schwarz.  Scale dependence of cosmological backreaction.  Phys. Rev. D, 78(8):083531, October 2008.
  * (5) Jr. Bengaly, C. A. P.  Constraining the local variance of H0 from directional analyses.  JCAP, 2016(4):036, April 2016.
  * (6) Krzysztof Bolejko, M. Ahsan Nazer, and David L. Wiltshire.  Differential cosmic expansion and the Hubble flow anisotropy.  JCAP, 2016(6):035, June 2016.
  * (7) Asta Heinesen and Thomas Buchert.  Solving the curvature and Hubble parameter inconsistencies through structure formation-induced curvature.  Classical and Quantum Gravity, 37(16):164001, August 2020.
  * (8) Stephen R. Green and Robert M. Wald.  New framework for analyzing the effects of small scale inhomogeneities in cosmology.  Phys. Rev. D, 83(8):084020, April 2011.
  * (9) Stephen R. Green and Robert M. Wald.  A simple, heuristic derivation of our ‘no backreaction’ results.  Classical and Quantum Gravity, 33(12):125027, June 2016.
  * (10) Thomas Buchert.  On Average Properties of Inhomogeneous Fluids in General Relativity: Dust Cosmologies.  General Relativity and Gravitation, 32:105–126, January 2000.
  * (11) T. Buchert, M. Carfora, G. F. R. Ellis, E. W. Kolb, M. A. H. MacCallum, J. J. Ostrowski, S. Räsänen, B. F. Roukema, L. Andersson, A. A. Coley, and D. L. Wiltshire.  Is there proof that backreaction of inhomogeneities is irrelevant in cosmology?  Classical and Quantum Gravity, 32(21):215021, November 2015.
  * (12) Thomas Buchert.  On Backreaction in Newtonian cosmology.  Monthly Notices of the Royal Astronomical Society, 473(1):L46–L49, January 2018.
  * (13) Thomas Buchert, Pierre Mourier, and Xavier Roy.  Cosmological backreaction and its dependence on spacetime foliation.  Classical and Quantum Gravity, 35(24):24LT02, December 2018.
  * (14) Boudewijn F. Roukema, Jan J. Ostrowski, and Thomas Buchert.  Virialisation-induced curvature as a physical explanation for dark energy.  JCAP, 2013(10):043, October 2013.
  * (15) Mauro Carfora and Francesca Familiari.  A comparison theorem for cosmological lightcones.  Letters in Mathematical Physics, 111(2):53, April 2021.
  * (16) Thomas Buchert, Henk van Elst, and Asta Heinesen.  The averaging problem on the past null cone in inhomogeneous dust cosmologies.  arXiv e-prints, page arXiv:2202.10798, February 2022.
  * (17) M. Gasperini, G. Marozzi, F. Nugier, and G. Veneziano.  Light-cone averaging in cosmology: formalism and applications.  JCAP, 2011(7):008, July 2011.
  * (18) H. Bondi.  Spherically symmetrical models in general relativity.  Monthly Notices of the Royal Astronomical Society, 107:410, January 1947.
  * (19) P. Szekeres.  A class of inhomogeneous cosmological models.  Communications in Mathematical Physics, 41(1):55–64, February 1975\.
  * (20) Tirthabir Biswas and Alessio Notari.  ’Swiss-cheese’ inhomogeneous cosmology and the dark energy problem.  JCAP, 2008(6):021, June 2008.
  * (21) Pierre Fleury, Hélène Dupuy, and Jean-Philippe Uzan.  Interpretation of the Hubble diagram in a nonhomogeneous universe.  Phys. Rev. D, 87(12):123526, June 2013.
  * (22) Mikko Lavinto, Syksy Räsänen, and Sebastian J. Szybka.  Average expansion rate and light propagation in a cosmological Tardis spacetime.  JCAP, 2013(12):051, December 2013.
  * (23) Mikko Lavinto and Syksy Räsänen.  CMB seen through random Swiss Cheese.  JCAP, 2015(10):057, October 2015.
  * (24) S. M. Koksbang.  Understanding the Dyer-Roeder approximation as a consequence of local cancellations of projected shear and expansion rate fluctuations.  Phys. Rev. D, 104(4):043505, August 2021.
  * (25) S. Nájera and R. A. Sussman.  Pancakes as opposed to Swiss cheese.  Classical and Quantum Gravity, 38(1):015016, January 2021.
  * (26) Hayley J. Macpherson, Paul D. Lasky, and Daniel J. Price.  Inhomogeneous cosmology with numerical relativity.  Phys. Rev. D, 95(6):064028, March 2017.
  * (27) Julian Adamek, David Daverio, Ruth Durrer, and Martin Kunz.  gevolution: a cosmological N-body code based on General Relativity.  JCAP, 2016(7):053, July 2016.
  * (28) Julian Adamek, Ruth Durrer, and Martin Kunz.  N-body methods for relativistic cosmology.  Classical and Quantum Gravity, 31(23):234006, December 2014.
  * (29) Maxim Eingorn, A. Emrah Yükselci, and Alexander Zhuk.  Screening vs. gevolution: In chase of a perfect cosmological simulation code.  Physics Letters B, 826:136911, March 2022.
  * (30) Viraj A. A. Sanghai and Timothy Clifton.  Post-Newtonian cosmological modelling.  Phys. Rev. D, 91(10):103532, May 2015.
  * (31) Viraj A. A. Sanghai and Timothy Clifton.  Erratum: Post-Newtonian cosmological modelling [Phys. Rev. D 91, 103532 (2015)].  Phys. Rev. D, 93(8):089903, April 2016.
  * (32) Viraj A. A. Sanghai, Pierre Fleury, and Timothy Clifton.  Ray tracing and Hubble diagrams in post-Newtonian cosmology.  JCAP, 2017(7):028, July 2017.
  * (33) E. M. Lifshitz.  On the gravitational stability of the expanding universe.  Zhurnal Eksperimentalnoi i Teoreticheskoi Fiziki, 16:587–602, January 1946.
  * (34) E. Lifshitz.  Republication of: On the gravitational stability of the expanding universe.  General Relativity and Gravitation, 49(2):18, February 2017.
  * (35) James M. Bardeen.  Gauge-invariant cosmological perturbations.  Phys. Rev. D, 22(8):1882–1905, October 1980.
  * (36) Hyerim Noh.  Fully nonlinear and exact perturbations of the Friedmann world model: non-flat background.  JCAP, 2014(7):037, July 2014.
  * (37) Jinn-Ouk Gong, Jai-chan Hwang, Hyerim Noh, David Chan Lon Wu, and Jaiyul Yoo.  Exact non-linear equations for cosmological perturbations.  JCAP, 2017(10):027, October 2017.
  * (38) Bo Wang and Yang Zhang.  Second-order cosmological perturbations. I. Produced by scalar-scalar coupling in synchronous gauge.  Phys. Rev. D, 96(10):103522, November 2017.
  * (39) Adam J. Christopherson and Karim A. Malik.  Practical tools for third order cosmological perturbations.  JCAP, 2009(11):012, November 2009.
  * (40) Ruslan Brilenkov and Maxim Eingorn.  Second-order Cosmological Perturbations Engendered by Point-like Masses.  The Astrophysical Journal, 845(2):153, August 2017.
  * (41) Jai-Chan Hwang and Hyerim Noh.  Second-order perturbations of a zero-pressure cosmological medium: Comoving versus synchronous gauge.  Phys. Rev. D, 73(4):044021, February 2006.
  * (42) Z. A. Golda, A. Woszczyna, and K. Zawada.  Canonical Gauge-Invariant Variables for Scalar Perturbations in Synchronous Coordinates.  Acta Physica Polonica B, 36(6):2133, June 2005.
  * (43) Szymon Sikora and Krzysztof Głód.  Example of an inhomogeneous cosmological model in the context of backreaction.  Phys. Rev. D, 95(6):063517, March 2017.
  * (44) Szymon Sikora and Krzysztof Głód.  Perturbatively constructed cosmological model with periodically distributed dust inhomogeneities.  Phys. Rev. D, 99(8):083521, April 2019.
  * (45) Szymon Sikora and Krzysztof Głód.  Construction of the cosmological model with periodically distributed inhomogeneities with growing amplitude.  European Physical Journal C, 81(3):208, March 2021.
  * (46) Sophia R. Goldberg, Timothy Clifton, and Karim A. Malik.  Cosmology on all scales: A two-parameter perturbation expansion.  Phys. Rev. D, 95(4):043503, February 2017.
  * (47) Sophia R. Goldberg, Christopher S. Gallagher, and Timothy Clifton.  Perturbation theory for cosmologies with nonlinear structure.  Phys. Rev. D, 96(10):103508, November 2017.
  * (48) Christopher S. Gallagher, Timothy Clifton, and Chris Clarkson.  Multi-scale perturbation theory II: Solutions and leading-order bispectrum in the $\Lambda$CDM universe.  JCAP, 2021(8):048, August 2021.
  * (49) Mariusz P. Dabrowski and Jerzy Stelmach.  Astrophysical formulas for Friedman models with cosmological constant and radiation.  The Astronomical Journal, 94:1373–1379, November 1987.