Construction of the cosmological model with periodically distributed inhomogeneities with growing amplitude

We construct an approximate solution to the cosmological perturbation theory around Einstein-de Sitter background up to the fourth-order perturbations. This could be done with the help of the specific symmetry condition imposed on the metric, from which follows, that the model density forms an infinite, cubic lattice. We show that the perturbative solution obtained this way can be interpreted as the exact solution to the Einstein equations for a dust-like energy-momentum tensor. In our model, it seems that physical quantities averaged over a large scales overlap with the respective Einstein-de Sitter prediction, while local observables could differ significantly from their background counterparts. As an example, we analyze in details a behaviour of the local and the global measurements of the Hubble constant, which is important in the context of a current Hubble tension problem.


Introduction
The studies within the cosmological perturbation theory up to second order yield that the influence of inhomogeneities on the Hubble diagram is visible but small and can reach at most one percent level in the parameters estimation [1,2,3,4,5,6,7]. Similar results are reached by ray tracing into Newtonian N-body numerical simulations [8,9,10], relativistic N-body numerical simulations [11,12,13,14] or relativistic hydrodynamical numerical simulations of an inhomogeneous dust model [15,16,17,18,19,20]. It is suggested that the emergence of spatial curvature during structure formation could be at least partially responsible for this effect [21,22,23,24,25].
The negligible impact of inhomogeneities on the global evolution and observational parameters of cosmological models is supported by general considerations concerning the backreaction effect [26,27,28,29,30,31] and perturbative analysis of weak gravitational lensing [32]. However, these approaches are often criticized as incomplete or inconclusive because of restrictive assumptions made [33,34,35,36]. There are also works which do not confirm the irrelevance of inhomogeneities despite similar methods and techniques used in studies [37,38,39]. Moreover, it is argued that the phenomena of virialization of clusters and volume dominance of voids significantly affect the Hubble diagram and may even explain dark energy [40,41,42,43,44].
Investigations to light propagation in inhomogeneous cosmologies have brought the development of various constructions for the models.
They include a Swiss-cheese arrangement of the Lemaître-Tolman or the Szekeres holes into the Robertson-Walker background space-time [45,46,47,48,49,50] or a regular lattice of the Schwarzschild spheres approximately glued by the Lindquist-Wheeler technique [51,52,53,54,55,56,57,58]. With the help of numerical simulations, it is also possible to study evolving configurations consisting of interacting black holes [59,60,61,62,63,64,65,66,67,68]. An exceptional example of model with point masses periodically distributed on a cubic lattice, presented in [69,70], can be built as a perturbative solution to the Einstein equations for a dust. The post-Newtonian formalism is another framework utilizing a perturbative approach in which a model is hierarchically constructed from autonomous cells whose matching conditions determine an overall dynamics [71,72,73,74,75,76].
Recently, we have proposed a cosmological model containing a continuous periodic distribution of inhomogeneities which are characterized by the equation of state of a dust [77,78]. The space-time of the model is constructed perturbatively and the matter is treated within the fluid hydrodynamic formalism. Here, we present a substantial extension of this model.

Model construction
Consider the following space-time metric in the Cartesianlike coordinates (t, x, y, z): where the six metric functions cij = cij (λ, t, x, y, z) depend on the space-time coordinates and some auxiliary parameter λ. We adopt geometrized units for which c = 1 and G = 1. Since the Christoffel symbols Γ µ 00 = 0 vanish for all µ, the vector U µ = (1, 0, 0, 0) is always tangent to some time-like geodesic. The worldlines of dust particles are geodesics, therefore for the universe filled with dust, the coordinates we use are comoving.
The task we want to address in this article is the following. Suppose that the Einstein equations are satisfied exactly. How to find the metric functions cij so that the energy-momentum tensor Tµν = Gµν /(8π) is a dust-like energy-momentum tensor and the density distribution has inhomogeneities with the amplitude growing during the time evolution?
By the dust-like energy-momentum tensor we mean such Tµν for which all of the elements, except the density ρ = T00, are negligible compared to ρ. A variety of exact solutions to the Einstein equations, which tend to describe the Universe in the matter-dominated era, assume the dust energy-momentum tensor for simplicity. In that case, all of the Tµν elements other than the density are exactly zero. However, many physical processes could contribute to these energy-momentum tensor elements. For example, proper motions of the galaxies act as the pressure, while this contribution is very small. This justifies the dust-like assumption for Tµν in the matter-dominated era if only smallness of the elements other than the density is properly controlled.
The task of finding metric functions cij for which the energy-momentum tensor has the properties described above is quite complicated in general. To handle it in the specific case we consider two simplifying assumptions: (i) The metric functions cij are invariant over every permutation of the spatial variables (x, y, z).
(ii) The metric functions cij are analytic in λ and the parameter λ is proportional to the amplitude of the inhomogeneities. Additionally, we assume that for λ = 0 we recover the Einstein-de Sitter (EdS) spacetime, which is the spatially flat universe homogeneously filled with dust. This means that the metric functions in this limit reads: cij (0, t, x, y, z) = a(t) 2 for i = j and zero otherwise. The scale factor in this case is a(t) = C t 2/3 , where C is a constant.
The assumption (ii) enables us to consider the usual perturbation theory around the EdS background. We may take Taylor expansion of the metric and resulting energy-momentum tensor around λ = 0: and analyze k-order metric elements g (k) ij and k-order µν order by order. The form of the metric (1) implies that perturbations are performed in the synchronous comoving gauge in each order. The scalar perturbations in the linear order admit two modes: the decaying and the growing one. In our prevuous paper [78] we analyze the decaying mode consistent with the assumption (i) in the orders higher than the linear order. In the current paper we concentrate on the growing mode. We will start with the linear perturbations and then we correct the solution in the consecutive orders. At the end of this procedure, we get the desired exact dust-like solution, finding of which is the main goal of the present paper.

Perturbation theory in the linear order.
We assume the following ansatz for the spatial part of the metric in the k-order, which is consistent with the symmetry assumption (i): where the following abbreviated notation is used: In each order we introduced three functions A (k) , B (k) , F (k) and three coefficients α k , β k and φ k , which should be specified in the process of the model construction. The symmetry condition (i) implicates that the energy-momentum tensor in each order k ≥ 1 has the four types of components: The structural form of the components within each type is invariant over every permutation of spatial variables.
Let us start with the linear order perturbations. If we specify the powers: then the elements T (1)0 i and T (1)i j | i =j are equal to zero and all of the terms of T (1)0 0 together with all of the terms of T (1)i j |i=j has a simple power law dependence on time. For simplicity, we put the function F (1) = 0. Then, we may cancel out the pressure-like terms T (1)i j |i=j = 0 demanding that the function B (1) satisfies differential equation: After that, the first order density ρ (1) ≡ −T (1)0 0 is: This way we obtain exact dust solution to the cosmological perturbation theory in the linear order, in which the density distribution in space is given by the arbitrary function A (1) . Let us analyze in details the model for one exemplary function A (1) given below.
In the beginning, we comment on the dimensions of the physical quantities important for model construction. We are working in the geometrized units c = 1 and G = 1. In this system of units, all of the physical quantities are dimensionless or their dimension can be expressed as some power of the unit of length. We choose the megaparsec as basic unit of length. Then, the age of the EdS universe t (EdS) 0 = 9.32 Gyr, compatibile with the Hubble constant value H0 = 70 km/s/Mpc, reads t (EdS) 0 = 2855.57 Mpc. There is a natural convention of normalizing the scale factor to unity at the age of the universe a(t (EdS) 0 ) = 1. From that follows the value of the constant C = 4.97 × 10 −3 . The density of the EdS model ρ (0) = 1/(6πt 2 ) evaluated at the universe age defines the critical density ρcr = ρ (0) (t (EdS) 0 ), which value is ρcr = 6.51 × 10 −9 Mpc −2 . The critical density introduces a natural density scale. We define the quantity Ω := ρ/ρcr as a density measured in the critical density units. Now, we choose the function A (1) as: where s0 = 1, s1 = 0.5, B0 = π/25 and B1 = π/5. If we fix the lambda parameter value as λ = 4.42 × 10 −4 , then the maximal density at t (EdS) 0 is Ω = 1.2. Density distribution illustrated in Fig. 1 forms periodic, cubic lattice. The linear size of the elementary cell at t (EdS) 0 is around 50 Mpc. Because the metric is nontrivial, the length of a segment depends on its orientation and position in space. The linear size of the elementary cell measured near to the overdense region gives the value slightly lower than 50 Mpc, while the result of the same measurement performed close to the underdense region could be slightly higher than 50 Mpc. For a scales much larger than the size of the elementary cell the model universe becomes homogeneous and isotropic in common sense, and FLRW space-time arises as a natural candidate for an effective average model. Although the density distribution profile is restricted by the model symmetry, quite complicated distributions are allowed. In the given example, one can identify large overdensity and underdensity regions of a size comparable to the typical size of superclusters of galaxies and smaller substructures with a size around a few megaparsecs.

Perturbation theory in the second
order.
The function B (1) which is a solution of the Eqn. 5 for the above proposition of the arbitray function A (1) is: In the above formula there appears two small constants: (C/B0) 2 = 1.56 × 10 −3 and (C/B1) 2 = 6.25 × 10 −5 , where C 2 = 2.47 × 10 −5 . The metric function B (1) is then much smaller than the metric function A (1) . However, in the second-order energy-momentum tensor elements appear some terms containing derivatives of the function B (1) divided by C 2 . These terms are the same order of magnitude as the terms with the function A (1) alone. For that reason, one cannot neglect the metric function B (1) in the first-order metric g (1) , if the second-order energymomentum tensor is considered. Nevertheless, the existence of these small constants enables one to identify the leading terms in the second-order energy-momentum tensor and to neglect the terms which are much smaller in comparison to the leading terms.
From now on, we denote by the symbol ≈ the approximation of some expression, in which all of the terms proportional to C 2 are neglected. Because some subexpressions could have a different time dependence, the validity of this approximation should be checked in each epoch of time. Lets write the leading terms of T (2)x x.
In the beginning, time dependence of different terms was not the same. However, it can be simplified to a single power-law t −2/3 , when we fix the following values of the powers: If we wish that the terms containing the second-order metric functions B (2) and F (2) are the same order of magnitude as the terms containing functions A (k) , then functions B (2) and F (2) should be proportional to C 2 . We will verify this assumption at the end of the presented procedure. In that case, T (2)0 i ≈ 0 and T (2)i j ≈ 0 for i = j. The corresponding elements T (2)y y and T (2)z z one would obtain by performing permutation of the spatial variables in the formula Eq. 9.
On the right-hand side of Eq. 9 it is possible to separate terms depending on two variables. These terms are equal to zero when the following differential equation is satisfied: The remaining terms depending on one variable can be canceled out by the following condition: If the above differential equations are satisfied, then all of the elements T (2)i j ≈ 0 for i = j. This is guaranteed by the symmetry condition imposed on the metric.
The conditions Eq. 11 and Eq. 12 enable to simplify the form of the second-order density ρ (2) = −T (2)0 0. As the result, one gets: Specification of the arbitrary function A (2) determines the second-order density. In our example, we choose this function such that perturbed density at the second order has the same spatial distribution as the first-order density perturbation: We introduced here the parameter K, which will control the growth rate of the density contrast. In this case, the function A (2) takes the form: The right-hand sides of equations Eq. 11 and Eq. 12 depend only on the first-order metric functions and the arbitrary function A (2) . In the case of our example, these functions are composed of the trigonometric functions and it is very simple to find solutions to Eq. 11 and Eq. 12. Taking the constants of integration equal to zero for simplicity we obtained: and As it was expected, the functions F (2) and B (2) are proportional to C 2 , so the method of construction of the metric functions is self-consistent. This way we end up with a dust-like solution up to the second order of the perturbation theory.

Third and fourth order perturbations.
It is possible to apply the procedure given in the previous subsection in the consecutive orders of the perturbation theory. First, we fix the values of the powers: which appear in the time evolution part of the metric functions. This way we simplify the time dependence of the energy-momentum tensor subexpressions to a single power law. Next, we assume that metric functions B (k) and F (k) are proportional to C 2 and neglect the terms of the energy-momentum tensor which are small in comparison to the leading-order terms. Since T (k)0 i and T (k)i j | i =j contains only the functions B (k) and F (k) and the constant C 2 is small, we may expect that T (k)0 i ≈ 0 and T (k)i j ≈ 0 for i = j. In the remaining elements T (k)i j |i=j one can identify terms depending on two variables, which can be cancel out when the function F (k) satisfies appropriate differential equation similar to Eqn. 11. Then, in the formula for T (k)i j |i=j remain some terms depending on the one variable only. They can be set to zero, by demanding that the function B (k) satisfies some differential equation similar to Eqn. 12.
The differential equations for the metric functions F (k) and B (k) enable one to simplify the formula for the korder density ρ (k) . In effect, ρ (k) depends only on the metric functions A (l) , for l ≤ k. This way, specification of the arbitrary metric function A (k) determines the spatial profile of the k-order density. In our example, we analyze the case for which the spatial distribution of the density in each order is the same. This means that: for k ≥ 2. The time dependence t (2k−6)/3 of the k-order density is a consequence of the specific values of the powers Eq. 18.
The right-hand sides of the differential equations for F (k) and B (k) depend on the metric functions F (l) and B (l) known from the previous orders l < k and the metric functions A (m) , for m ≤ k. After the arbitrary function A (k) is fixed, one can solve these differential equations and obtain the resulting F (k) and B (k) . It is easy to verify that these functions are proportional to C 2 , so the method is self-consistent.
We apply the presented method up to fourth order of the perturbation theory. The resulting metric functions are made of simple trigonometric and monomial functions, but the expressions become quite large and we decided not to display them.

Exact solution
By the method presented in the previous subsections we construct a dust-like solution to the fourth-order perturbation theory. Now, we would like to use the solution we obtained so far as a guess of the metric functions cij(λ, t, x, y, z) of some exact solution to the Einstein equations. More precisely, we assume that the space-time metric of the model universe is a fourth-order polynomial in the parameter λ: where the metric functions in each order g (k) were constructed by the method described in the subsections 2.1-2.3. In the next section we verify that the energymomentum tensor Tµν = Gµν/(8π) refers to some dustlike matter distribution and we analyze in details the model properties.

Model properties 3.1 Density distribution
Let us begin with the analysis of the distribution of the model density ρ = −T 0 0. In the model construction method presented in the previous section, in each order k ≥ 1 one can specify the shape of the k-order density distribution ρ (k) by means of arbitrary function A (k) . We fix the A (1) function as in the definition Eq. 7. The functions A (k) for 1 < k ≤ 4 has been specified so that the density ρ (k) has the same spatial distribution as ρ (1) . The formula for ρ (k) is given by the Eqn. 19. The density in each order ρ (k) has a different power in a time dependence. To some extent, we can control the growth rate of the inhomogeneities by manipulating the contribution of each order ρ (k) to the total density. In our example, the parameter K describes the contribution of the high-order densities ρ (k) | 1<k≤4 in comparison to the linear order contribution ρ (1) , given by the Eq. 6.
From now on, we fix the parameter λ = 4.42 × 10 −4 . For this value of λ the maximal density at the EdS universe age t (EdS) 0 up to the first order is where Ω (k) = ρ (k) /ρcr. For K = 0 we expect that the firstorder density gives a dominant contribution to the total density. Indeed, the total density expressed in the units of the critical density Ω = ρ/ρcr evaluated at the maximum of the overdensity region xO = (12.5, 12.5, 12.5) and at the t (EdS) 0 is Ω = 1.1999. The shape of the isodensity surfaces remains practically unchanged in comparison to the isodensity surfaces of the perturbation theory up to the first order. Therefore, density of the full model ρ is approximated very well by the density of the perturbation theory up to the first-order perturbations. Also, one can think about Figure 1 as it describes the density distribution in space of the density of the full model ρ. We examined also two specific values K = 0.1 and K = −0.1, for which the density of the full model is very close to the density of the fourth-order perturbation theory.
The density has a spatial distribution, which forms an infinite, cubic lattice. For a scales much larger than the elementary cell the model becomes homogeneous and isotropic in common sense. It is then reasonable to approximate our inhomogeneous universe by the FLRW solution with some average density distribution ρ if only scales much larger than the elementary cell are considered. We will use a natural definition of the average of some physical quantity f over the elementary cell D, at the hypersurface of the constant time t: It should be pointed out that the volume element is not trivial and depend on the position in space. The isosurfaces of the square root of the spatial part of the metric determinant are shown in Fig. 2. Comparison of this figure with Fig. 1 shows that underdensity regions have larger values of the volume element than overdensities. In Figure 3 we present by the solid blue curve the model average density over the elementary cell Ω D (t) as a function of time, calculated for the case K = 0. For two other values of K = 0.1 and K = −0.1 we obtain the same result. For comparison, by a red, dashed curve we plot on the same figure the density of the Einstein-de Sitter model, which was used as a background space-time g (0) for the model construction. It is evident that both curves overlap, so the density of the EdS model is indeed an averaged density of the full inhomogeneous model considered here.
The notion of the average density enables one to define the density contrast as: In Figure 4 we present the density contrast as a function of time for three values of K. For K = 0 density contrast grows with time exactly as the growing mode of the first-order perturbation theory, where δ ∝ t 2/3 . For the parameter K = 0.1 density contrast grows faster than t 2/3 , while for the value K = −0.1 it grows slower. It is interesting to notice, that the density contrast of the exact solution to the Einstein equations could differ from the prediction of the first-order perturbation theory. Important difference appears when second and higher-order terms contribute significantly to the total density.

Is the energy-momentum tensor dust-like?
We developed our model within the framework of the perturbation theory up to the fourth order. Then, we treated the resulting fourth-order polynomial in λ parameter (Eqn. 20) as a metric of the full model, for which the Einstein equations are satisfied exactly. In effect, we have to deal with the contributions to the energy-momentum tensor coming from fifth and higher orders, which we do not control. We have to check whether the resulting energymomentum tensor of the full model Tµν = Gµν /(8π) resembles the properties of the energy-momentum tensor from the fourth-order perturbation theory.
In the previous subsection, we have checked that the density of the full model practically does not change in comparison with the EdS model density perturbed up to the fourth order. Now, we analyze the values of the other elements of the energy-momentum tensor. Because of the symmetry condition imposed on the metric, there are four types of the T µ ν components. In Figure 5 we present the absolute value of the remaining three types of the energy-momentum elements relative to the density.  The blue curve corresponds to the elements T i j |i=j representing the pressure measured by the observer which is in rest in the comoving reference frame (t, x, y, z). To verify whether the values of these elements are small enough we have to consider some interpretation of the source of this pressure.
If the model energy-momentum tensor describes galaxies which are members of rich galaxy clusters and have its proper motions, then the energy-momentum tensor could be interpreted in the framework of the Jeans theory of a collisionless system of particles. Within this model framework, the stress-energy tensor (the spatial part of the energy-momentum tensor) is directly related to the velocity dispersion tensor of these particles σij: Since the density is positive, the elements of T i j should be positive also. Unfortunately, the resulting pressure T i j |i=j is negative in some regions. However, this problem can be solved in the following way. We increase a bit the pressure by adding the small positive term P t −2/3 to the right-hand side of the second-order formula Eqn. 9, where P = 1.006 × 10 −4 . Then, after recalculation of the metric functions, one can check that in the case K = 0 the order of magnitude of the energy-momentum tensor elements remain unchanged in comparison to the case P = 0, but the pressure T i j |i=j is always positive within the elementary cell. For other values of K the corresponding P should be different. In Figure 6 we present one of the elements of the velocity dispersion tensor σxx, which is related to the model energy-momentum tensor element T x x by the formula Eqn. 23, for the case K = 0 and at the time t (EdS) 0 . This figure shows that the values of the pressure of order 10 −6 in the geometrized units correspond to the velocity dispersion of order of 1000 km/s. This value of the velocity dispersion is consistent with observations of galaxy clusters. The nondiagonal elements T i j | i =j are very close to zero. The order of magnitude of these elements is 10 −10 . This suggests that the distribution of the velocities should be isotropic. As can be seen in Fig. 6, this is not the case of the resulting σij. However, the order of magnitude of the spatial part of the resulting energy-momentum tensor T i j is consistent with the values of the velocity dispersion found in the galaxy clusters.
The energy flux terms T 0 i are one order of magnitude lower than the pressure-like terms T i j |i=j , therefore we conclude, that the resulting energy-momentum tensor of the full inhomogeneous solution considered here is really dust-like.

Curvature of space.
The EdS model is the background space-time for the perturbative construction scheme presented in Section 2. The EdS universe is spatially flat. In this subsection, we will analyze the behavior of the spatial curvature of the full model.
In the Einstein-de Sitter model understood as a spe- cial case of the Friedmann-Lemaître cosmological model, additionally, there vanish the curvature scalar of hypersurfaces orthogonal to the fluid flow R and the isotropic pressure p. By the Stewart-Walker lemma [79], in the perturbed Einstein-de Sitter model at the first order, these two scalar fields are gauge-invariant. In contrast, the perturbation of the matter energy density ρ is not gauge-invariant but one can consider its spatial gradient Xν = P ν α ∇αρ as a suitable perturbative quantity. Geometric, kinematic and dynamic gauge-invariant quantities which characterize properties of the perturbed spacetime, energy-momentum field and its flow are mutually related by the Ellis-Bruni equations [80]. In special cases, when this set of equations become closed, it provides analytic solution for the behavior of inhomogeneities. When we restrict our considerations only to perturbations of the scalar type then the fluid flow is necessarily irrotational and the magnetic part of the Weyl tensor vanishes [81]. If we further assume that the perturbed flow is geodesic and the fluid is nonconductive and inviscid then we arrive at the following set of equations for perturbative quantities where P µν = U µU ν +gµν is the projection tensor and θ is the expansion scalar in the background. It appears that the imposed assumptions do not allow for nonzero pressure perturbations. Furthermore, they confine only the temporal evolution of perturbations leaving their spatial variability free. Curvature perturbations decrease with the expansion of the model but solving for density perturbations reveals two separate modes, one decaying and one growing. These modes differ in their physical nature, since the growing mode is govern entirely by the curvature perturbations of hypersurfaces. In particular, imposing zero curvature perturbations on hypersurfaces eliminates the growing mode of density perturbations. After these general considerations, let us go back to the specific exact solution presented in Sec 2.4. The scalar curvature of the hypersurface of a constant time R is conventionally related to the quantity Ω k = −R/(6 H(t) 2 ), where H(t) is the Hubble parameter. We used the Hubble parameter H(t) of the background EdS universe. In Figure 7 we show a dependence of the Ω k on the position within the elementary cell, at the time t (EdS) 0 . The overdense regions have negative values of Ω k , so the scalar curvature is positive in these regions. Within the underdense regions the situation is opposite. The Ω k parameter is positive and the scalar curvature R is negative there. is seen that the scalar curvature decreases with time and space tends to flatten during the time evolution. The orange curve corresponds to the model with the value K = 0.1, while the light-green curve refers to K = −0.1. We note that for the case K = 0.1, for which the growth rate of the inhomogeneities is higher than for the case K = −0.1, the scalar curvature R decreases more slowly than for the case with K = −0.1. This shows, that the growth rate of the inhomogeneities is related to the scalar curvature R behavior. Finally, we calculate the average over the elementary cell of the parameter Ω k . The results for three different values of K and for different instants of time are plotted on Fig. 9. In the paper [23], the authors based on their silent universe model suggest that the scalar curvature of space could emerge with time. In our case, we observe that the average Ω k parameter grows slightly with time, however, values of Ω k D are very small. The growth rate of Ω k D depends on the K parameter, but in each case, the values of the Ω k D are smaller than 10 −2 . This means that on average the space is almost flat.

Local measurments of the Hubble constant.
At the end of this paragraph let us analyze expansion of our inhomogeneous universe. In Figure 10 we of this figure we deduce that underdense regions expand faster than overdense regions, although on average the model universe expand exactly as the Einstein-de Sitter homegeneous case. Therefore, local measurments of the Hubble constant could differ from the EdS Hubble parameter H(t) evaluated at the universe age t (EdS) 0 , while the measurements of the Hubble constant on basis of some obervables related to scales much larger than the elementary cell should be consistent with the EdS prediction.
To simulate how observer living in our inhomogeneous universe would perform local measurement of the Hubble constant, we made the following numerical experiment. For a given observer position x0 = (x0, y0, z0) we generate ten random directions (θ, φ) with probability distribution uniform on the unit sphere. For each direction we generate randomly ten points belonging to the line γ(l) = (t (EdS) 0 , x0+l sin θ cos φ, y0 +l sin θ sin φ, z0+l cos θ). This way we simulate one hundred sources distributed randomly in the close neighborhood of the observer. To generate the Hubble diagram we have to calculate for each source the physical distance to the observer d and its time derivativeḋ. The physical distance we obtain by numerical integration: Since the metric elements depend explicitly on time, to get the respectiveḋ we need to take the time derivative of the integration kernel and perform numerical integration again.  The resulting Hubble diagram generated for two different observer's positions is plotted in Figure 11. The blue points correspond to the observer located at the maximum density xO, whereas the light-blue points are generated for the observer located at the minimum density xU . For the points in the range of distances d ∈ (5, 20) Mpc we perform the linear fit to get the resulting local Hubble constant. For the observer located at the overdensity xO we obtain H0 = 69.12 km/s/Mpc, while the observer located at the underdensity measures H0 = 71.11 km/s/Mpc. It is clear that in both cases the resulting value differs significantly from the EdS Hubble constant H0 = 70 km/s/Mpc, which was fixed at the beggining of the perturbative approach described in Sec. 2. The difference we obtained here is slightly lower than the current observational difference between the Hubble constant estimation from CMB H0 = 67.37 ± 0.54 km/s/Mpc [82] and from the local measurements H0 = 74.03 ± 1.42 km/s/Mpc [83]. However, the order of magnitude of the difference we get here is comparable to the current observational difference. It is then reasonable to expect that inhomogeneities could play an important role concerning local Hubble constant measurements.
In our previous paper [78] we haven't noticed a dependence of the local measurement of H0 on the position within the elementary cell. However, in the previous model we considered smaller scale of the inhomogeneities, of the order of 3 Mpc. Currently, the scale of the inhomogeneous region is close to 25 Mpc with the additional substructures of the scale around 3 Mpc. It is quite clear that the local measurements of the Hubble constant should depend on the scale of the inhomogeneities under consideration.

Conlusions.
In the current paper, we constructed an example of the dust-like exact solution to the Einstein equations representing an inhomogeneous cosmological model with growing amplitude of the inhomogeneities. By the term dust-like we mean that an observer which is in rest in the comoving reference frame measures nonzero energymomentum tensor elements T µ ν other than the energy density −T 0 0, but which are negligible in comparison to −T 0 0. In the interpretation of the Jeans theory of the collisionless system of particles as a source of the pressurelike terms of the energy-momentum tensor, we checked that the values of the resulting energy-momentum tensor elements correspond to the particles velocity dispersion around 1000 km/s, which is reasonable value for the velocity dispersion of the galaxy cluster members.
The model construction method is based on the perturbation theory around the Einstein-de Sitter background. By using the additional simplifying symmetry condition and identifying leading terms of the energy-momentum tensor elements we are able to obtain the solution to the perturbation theory up to the fourth-order perturbations. Then, we consider the resulting fourth-order solution as a guess of the metric of some exact solution to the Einstein equations. We checked that this solution remains dust-like and analyze its basic properties.
Many of the current discussions concerning the possible influence of the inhomogeneities on the cosmological observations focus on the problem of backreaction. In these approaches, one asks whether the existence of the inhomogeneities affects properties of the average space-time.
Recently, most of the researchers conclude that the effect of backreaction is possible but it is rather negligible. In the presented model there is no backreaction effect at all. If we consider Eqn. 20 as a definition of the Green-Wald family of metrics gµν(λ) [84], then the effective energymomentum tensor tµν is zero since our model is based on the ordinary perturbation theory. Moreover, the average over the whole space of the density ρ D (t), the average curvature parameter Ω k D and the average expansion θ D overlap with the respective quantities of the background EdS model. At the same time, local physical quantities could differ significantly from the EdS background. We note nonnegligible differences considering the local volume element detgij d 3 x, local curvature of space and local expansion parameter θ. In effect, we show that some local measurements could differ from the expectations of the EdS background model. As the example, we verify that the local measurements of the Hubble constant could differ from the EdS value. Such an effect could possibly explain the current Hubble tension problem.
We want to stress, that our model is not a complete description of the real Universe. It is rather a step toward understanding the role of the large-scale inhomogeneities in the interpretation of the cosmological observables. The present paper is a large improvement of our previous work [78,77], but many issues could be done better in the future. First of all, there is a much more challenging task of considering perturbations around the nonflat background with nonzero cosmological constant. One can also look for more complicated density distributions beyond our symmetry restrictions. We point also, that in the present model the density, the scalar curvature of space and the expansion parameter are decreasing functions of time in all spatial positions. Therefore, the presented model could be interpreted as a descrip-tion for a large scale inhomogeneous regions behaviour but does not provide the framework for the formation of the individual structures for which we expect some kind of collapse and very high values of the density contrast. Nevertheless, the model gives an explicit example of an important influence of the inhomogeneities on the local observations of the Hubble constant, while at the same time the observables related to the scales much larger than the elementary cell overlap with the prediction of the background homogeneous model.