© The Authors 2023

1. Introduction

The recently emerged Hubble tension can be an indicator of certain, still unnoticed, genuine differences in the descriptions of the late and early Universe (Riess 2020; Riess et al. 2022; Di Valentino et al. 2021; Dainotti et al. 2021, 2023). Namely, the discrepancy between the Hubble constant values obtained from local surveys and global, cosmic microwave background data can imply the need to also consider nonidentical mechanisms of cosmic structure formation at local and global scales (e.g., Capozziello & Lambiase 2022; Bouchè et al. 2022).

The formation of the large-scale matter distribution of the cosmic voids and web, filaments, and clusters of galaxies are considered a result of the evolution of initial density fluctuations (Peebles 1993). Zeldovich pancake theory (Zeldovich 1970; Arnold & Shandarin 1982; Shandarin 1989; Shandarin & Sunyaev 2009) provides deep insights into certain features of cosmic structures.

Among various approaches to addressing the Hubble tension is the prediction of two Hubble flows, local and global ones, which are described by two similar sets of equations but nonidentical Hubble constants. The local flow is based on the following principal concept: Gurzadyan (1985) proved (as a theorem) that the general function for the force satisfying the identity of sphere-point gravity has the form

$$\begin{array}{c}\hfill \mathbf{F}(r)=(-\frac{A}{r^2}+\mathrm{\Lambda }r)\widehat{\mathbf{r}}.\end{array}$$(1)

It is remarkable that this function satisfies the first statement of Newton’s shell theorem (i.e., the sphere-point identity) but not the second part, namely, the force-free field inside a shell (for details, see Gurzadyan 2019; Gurzadyan & Stepanian 2018). Equation (1), as weak field general relativity (GR), allows the dynamics of groups and clusters of galaxies to be described (Gurzadyan & Stepanian 2019, 2020). Observational indications include the influence of halos on the properties of spirals in galaxies (e.g., Kravtsov 2013), which supports the notion of a non-force-free field inside a shell predicted by Eq. (1).

Due to the second term in Eq. (1), the cosmological constant, Λ, emerges in nonrelativistic cosmologies, namely those of McCrea & Milne (1934) and Zeldovich (1981). As a result, one has two equations that look similar but have drastically different contents (Gurzadyan & Stepanian 2021a,b):

$$\begin{array}{cc}\hfill H_{\mathrm{local}}^2& =\frac{8\pi G{\rho }_{\mathrm{local}}}{3}+\frac{\mathrm{\Lambda }{c}^2}{3},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill H_{\mathrm{global}}^2& =\frac{8\pi G{\rho }_{\mathrm{global}}}{3}+\frac{\mathrm{\Lambda }{c}^2}{3}.\hfill \end{array}$$(3)

The first equation defines the nonrelativistic (non-GR) local flow determined by the repulsion of the Λ term in Eq. (1) and, hence, with a local Hubble parameter. The second equation is the Friedmann equation for the Friedmann-Lemaitre-Robertson-Walker metric and flat geometry, with the global Hubble parameter. As shown in Gurzadyan & Stepanian (2021b), an analysis of these equations that takes the difference in the content of the matter mean densities into account – local and global ones (ρ_local and ρ_global) – explains the quantitative discrepancy between the local and global values of the Hubble parameter and thus the Hubble tension.

In this context, our present analysis assumes a principal difference in the treatments of the late (local) and early (global) Universe and hence the need of different techniques to address the formation of structures on local and global scales. Namely, while the evolution of primordial fluctuations is considered to be responsible for the global matter distribution, the role of the self-consistent gravitational interaction can play an important role in the formation of local structures. Hence, we continue the kinetic Vlasov analysis (Gurzadyan et al. 2022) of the formation of the cosmic voids and walls.

We show the emergence of structures, two-dimensional walls and voids, based on an analysis of the properties of solutions of kinetic Vlasov-Poisson system equations, where the fundamental role has the potential of Eq. (1) with the repulsive cosmological constant term. Particularly, using the methods of bifurcation theory and analyzing solutions of integral equations, we show the possibility of multiply connected two-dimensional gravitating structures appearing. The results predict the dependence of void scales on the local physical parameters.

2. Vlasov-Poisson equations for modified gravitational interaction with a cosmological constant

We began with consideration of the mechanism for the emergence of cosmological structures within a nonrelativistic McCrea-Milne model of expansion of the Universe with the modified gravitational interaction of Eq. (1) (Gurzadyan 1985). We considered a sphere of mass M(R) that contains particles of the same mass, m, such that the mass of particles in the ball M = ∫mf(x, v, t) dmdv (including in the case of a uniform distribution, M = 4πρR³/3), and the law of conservation of energy in the form

$$\begin{array}{c}\hfill \frac{1}{2}(\frac{\mathrm{d}R}{\mathrm{d}t}{)}^2-\frac{\mathit{GM}}{R}-\frac{c^2\mathrm{\Lambda }}{12}{R}^2=E=\mathrm{const}.\end{array}$$(4)

For the Hubble parameter, H (here defined as dR/dt = HR), we get the equation

$$\begin{array}{c}\hfill \frac{\mathrm{d}H}{\mathrm{d}t}+H^2+\frac{4\pi G\rho }{3}-\frac{1}{6}{c}^2\mathrm{\Lambda }=0.\end{array}$$(5)

We assumed the values R(t = t₀) = R₀, H(t = t₀) = H₀, and ρ(t = t₀) = ρ₀ at t = t₀. We could thus define the value of the constant $E=E_0=(1/2)H_0^2{R}_0-G·\frac{4\pi }{3}{\rho }_0{R}_0^2-\frac{c^2}{12}\mathrm{\Lambda }{R}_0^2$, and the initial conditions for the differential equation for the evolution of a sphere of radius R(t),

$$\begin{array}{cc}\hfill {\left(\frac{\mathrm{d}R}{\mathrm{d}t}\right)}^2& =\frac{a_1}{R}-a_2{R}^2-a_3,a_1\equiv \frac{8\pi }{3}G{\rho }_0{R}_0^3,a_2\equiv \frac{c^2\mathrm{\Lambda }}{12},\hfill \\ \hfill & a_3\equiv \frac{8\pi }{3}G{R}_0^3({\rho }_0-\frac{3H_0^2}{8\pi G}-\frac{3c^2\mathrm{\Lambda }}{8\pi G{R}_0}).\hfill \end{array}$$(6)

Solutions of Eq. (6) can be expressed in terms of degenerate elliptic Legendre integrals of the third kind (see, for example, Bateman & Erdelyi 1955):

$$\begin{array}{c}\hfill R(t)=I_3^{\ast }(-a_2,0,-a_3/6,a_1/4,0).\end{array}$$

At Λ ≡ 0, one has the classical open (infinitely expanding) and closed Friedmannian models. In this case, the Cauchy data for the considered ordinary differential equations are usually taken as R(t = 0) = 0, $R_t^{\prime }(t=0)=+\infty$. It should be noted that the case of the nonzero cosmological constant Λ does not lead to a fundamental change in the structure of solutions to the Cauchy problem (Eq. (6); see Vedenyapin et al. 2019, 2021).

The McCrea-Milne model with the modified potential of Eq. (1) is stable with respect to a change in the initial data, which is due to the self-similarity of its solutions and the theorem in Gurzadyan (1985) on the general form of the potential for a spherical region. Thus, to study the structure formation in the local Universe, one can use the Vlasov-Poisson equations without taking the GR effects into account. Obviously, the presence of an additional repulsive Λ term on the right side of the Poisson equation is not due to the GR metric but to Eq. (1).

We are interested in analyzing the properties of deviations from the state of local equilibrium of systems governed by Vlasov-Poisson equations. As the uniform distribution does not imply equilibrium in the general case in a system of massive particles with an additional oscillatory interaction of geometric genesis (for the total super-harmonic potential), one can study the evolution of matter in the vicinity of particle clusters as a result of the evolution of cosmological fluctuations. Then, one can introduce the concept of relative equilibrium, where the distribution of particles of a given equilibrium state satisfies the principle of maximum entropy and, for example, have the property of averaging over a set of hydrodynamic macro-parameters (see Vedenyapin et al. 2019, 2021).

The system of Vlasov-Poisson equations that describes the dynamics of a many-particle system for a d-dimensional case (d = 2, 3) for particles of equal masses, m, within the representation of a self-consistent gravitational field is written in the following form (Vlasov 1961, 1978),

$$\begin{array}{cc}& \frac{\partial F(\mathit{x},\mathit{v},t)}{\partial t}+{\mathrm{div}}_{\mathit{x}}(\mathit{v}F)-\frac{\partial F}{\partial \mathit{v}}{\mathrm{\nabla }}_{\mathit{x}}\mathrm{\Phi }(F)=0,\hfill \\ \hfill & \mathrm{\Phi }(F)\equiv {\displaystyle {\int }_{{\omega }_{{\mathit{x}}^{\mathbf{\prime }}}}{\int }_{{\omega }_{{\mathit{x}}^{\mathbf{\prime }}}}{\mathfrak{Y}}_d(|\mathit{x}-{\mathit{x}}^{\prime }|)F({\mathit{x}}^{\prime },{\mathit{v}}^{\prime },t_{\ast })\mathrm{d}{\mathit{x}}^{\prime }\mathrm{d}{\mathit{v}}^{\prime },}\hfill \end{array}$$(7)

where ${\mathfrak{Y}}_d(\mathit{x}-{\mathit{x}}^{\prime })\sim Y_d^{(1)}(|\mathit{x}-{\mathit{x}}^{\prime }|)+Y_d^{(2)}·{|\mathit{x}-{\mathit{x}}^{\prime }|}^2$ (the function corresponding to the attraction between particles) is $Y_d^{(1)}(\mathit{x}-{\mathit{x}}^{\prime })\in L_1^{\mathit{loc}}\bigcap C^2(\overline{\omega })$ and is integrable (and has a proper smoothness outside zero) in some bounded domain ${\omega }_d={\omega }_{\mathit{x}}\subset {\mathbb{R}}^d$ function. The $Y_d^{(2)}\in {ℝ}^1$ is a function that depends on the cosmological constant, Λ. We have assumed that, from the point of view of an external observer, the cosmological expansion affects the substructures of the quasi-two-dimensional physical structures in the same way as in the three-dimensional case.

It can be shown (see Sect. 3) that, for an arbitrary finite time $t_{\ast }\in {\mathbb{R}}_{+}^1$, the modified Poisson equation is

$$\begin{array}{cc}& {\mathrm{\Delta }}_{\mathit{x}}^{(d)}\mathrm{\Phi }(F){|}_{t=t_{\ast },\forall t_{\ast }\in {\mathbb{R}}_{+}^1}=\frac{2{\pi }^{d/2}}{\mathrm{\Gamma }(d/2)}G\rho (\mathit{x})-C_1,\hfill \\ \hfill & \rho (\mathit{x})\sim m{\displaystyle \int }F(\mathit{x},\mathit{v},t_{\ast })\mathrm{d}\mathit{v}(G_3\equiv G,C_1=C_1(\mathrm{\Lambda })).\hfill \end{array}$$(8)

Equation (8) takes the form of the so-called (inhomogeneous) Liouville–Gelfand equation (LGE; Gelfand 1963) if we take the density distribution of particles equal to that in the case of the Maxwell–Boltzmann integration distribution function. Hereinafter we do not take into account the dependence of the equilibrium distribution function on integrals of motion other than energy: F = F₀(x, v, t_*)∝Ω(ε[v(t_*);Φ(x)]|_ext) in terms of velocities (a regular branch of Φ is chosen). In this case, in the energy neighborhood of the local extremum, one has

$$\begin{array}{c}\hfill \mathrm{d}\mathrm{\Omega }/\mathrm{d}\epsilon {|}_{r\to r_{\mathrm{c}}}\approx 0,r_{\mathrm{c}}=(Gm/(3\mathrm{\Lambda }{c}^2){)}^{1/3}.\end{array}$$(9)

In this case, the force term in the Vlasov equation is either vanishing or possesses a turning point on the phase plane of the evolution of the system, and we get a relative equilibrium (or quasi-stationary) state as a solution of the Vlasov-Poisson equations system. It is characterized by the average quasi-equilibrium statistical temperature of the system, T, and the potential gauge or density, ρ₀ = Aexp(−Φ(0)/T), or the spatial distribution in a point: x = 0, as well as by the selected energy level of particles on the regular branch of the potential. For the sake of simplicity, we provide Eq. 10, an explicit form of the LGE for the two-dimensional case:

$$\begin{array}{cc}\hfill -{\mathrm{\Delta }}^{(2)}U& =\lambda exp(U)-c^2\mathrm{\Lambda }/T,U=-\mathrm{\Phi }/T,\hfill \\ \hfill \lambda & =2GmN{\rho }_0/T,T\equiv -\frac{\partial \mathrm{\Phi }/\partial x_j}{\partial \rho /\partial x_j}{|}_{j=1,2}.\hfill \end{array}$$(10)

From a mathematical point of view, the potentials (solutions of the LGE) correspond to the continual matter density distribution (in the $\omega \subset {\mathbb{R}}^d$ calculation area with given conditions on the ∂ω boundary). Obviously, this distribution is definitely idealized from a physical point of view. However, based on Theorem 1.3 in Esposito et al. (2005) and Brezis & Merle (1991), it can be argued that, when the internal parameters of the function λ(ρ₀, T) are changed in a limited area of space, there is a system of discrete masses, ${\mathit{q}}_i{|}_{i=\overline{1,k}}$ (defined in terms of singular points of Robin functions), whose action on a distant particle is equivalent to the given continuous distribution of matter:

$$\begin{array}{c}\hfill \frac{{\displaystyle \lambda exp(U_{\lambda })}}{{\displaystyle {\displaystyle {\int }_{\omega }}exp(U_{\lambda })\mathrm{d}\mathit{x}}}\rightharpoonup 8\pi {\displaystyle \sum _{i=1}^k}{\delta }_{{\mathit{q}}_i}.\end{array}$$(11)

Since there is scale invariance here, we can see that the replacement of discrete particles on a continuum distribution is valid both for a local single cluster (for example, understood as an elementary substructure) and for the system of two-dimensional walls. The gravitational potential in the neighborhood of the extremum point is close to constant, and the influence of the self-consistent fields on particle motion is absent; therefore, the inertial motion of particles leads to relatively long-duration formations whose local topology is not distorted.

We have assumed that the variation in the action over the gravitational field (resulting in the Poisson equation) and variation over particles (leading to the equations of motion along geodesics, and after the introduction of the distribution function particles to the Liouville equation) are independent operations. Thus, we can separate in the first approximation the particle distribution evolution equation and equations for self-consistent gravitational fields, assuming that particles move in a given field and that changes in the particle distribution function occur in accordance with the solution of the Liouville equations (Vedenyapin et al. 2011). Thus, we can study the equation for the potential separately, although it is essentially nonlinear. To study it, one needs to use methods of branching theory for solutions of nonlinear equations. However, important conclusions can be drawn even from the linearized version of the equation for the gravitational potential, which we demonstrate below.

3. Equation for the gravitational potential in integral representation

We considered a system of N point bodies with equal masses, m, interacting via the modified gravitational potential of Eq. (1). The Hamiltonian of the system has the form

$$\begin{array}{c}\hfill {({\mathfrak{H}}_N)}_{ϵ}={\displaystyle \sum _j^N(\frac{{\mathit{p}}_j^2}{2m}+{\sum _i}^{\prime }{\mathrm{\Phi }}_{ϵ}^{(2)}(|{\mathit{x}}_i-{\mathit{x}}_j|)+m{\mathrm{\Phi }}_{\mathrm{ext}}({\mathit{x}}_j)+B(\partial \omega ,{\mathit{x}}_j)),}\end{array}$$(12)

where Φ_ext(x_j) is the potential of the external gravitational field at a point (x_j), and B(∂ω,  bfx_j) is the potential energy of the interaction of particles with boundaries, which includes the reflection of particles from the boundary (periodicity of the boundaries) and the impact due to a change in the local temperature at the boundary of near-equilibrium particles inside the ω region (see Kiessling 1989, 1993). The index ϵ of the two-particle potential ${\mathrm{\Phi }}^{(2)}(r(\mathit{x},{\mathit{x}}^{\prime }))\sim c_1/|\mathit{x}-{\mathbf{x}}^{\prime }|+c_2{|\mathit{x}-{\mathit{x}}^{\prime }|}^2$ implies the regularization of a weakly polar function with a small argument; in other words, under the integral we replaced (e.g., in the three-dimensional case) the function Φ⁽²⁾(r) with the function

$$\begin{array}{cc}& {\mathrm{\Phi }}_{ϵ}^{(2)}(r(\mathit{x},{\mathit{x}}^{\prime }))=\hfill \\ \hfill & \{\begin{array}{cc}{(c_{1}r(\mathit{x},{\mathit{x}}^{\prime }))}^{-1}+c_2{r}^2(\mathit{x},{\mathit{x}}^{\prime }),\hfill & r(\mathit{x},{\mathit{x}}^{\prime })\ge ϵ;\hfill \\ {(2ϵ)}^{-1}·(3-{|r(\mathit{x},{\mathit{x}}^{\prime })|}^2/{ϵ}^2)+c_2{r}^2(\mathit{x},{\mathit{x}}^{\prime }),\hfill & r(\mathit{x},{\mathit{x}}^{\prime })<ϵ.\hfill \end{array}\hfill \end{array}$$(13)

In this case, the volume integral potential with a summable bounded density of particles is a continuous function up to the boundaries of the system. In the absence of significant gradients of macro-characteristics, the particles of the system can be described using the formalism of the canonical ensemble, characterized by the corresponding (near-)equilibrium density in the phase 6N-dimensional space:

$$\begin{array}{cc}& {\varrho }_N={(N!{\mathrm{\Omega }}_{ϵ})}^{-1}exp(-{\mathfrak{H}}_N/T),\hfill \end{array}$$(14)

$$\begin{array}{cc}& {\mathrm{\Omega }}_{ϵ}={(N!)}^{-1}{\displaystyle \int exp(-{({\mathfrak{H}}_N)}_{ϵ}/T)\prod _j\mathrm{d}{\mathit{x}}_j\mathrm{d}{\mathit{p}}_j.}\hfill \end{array}$$(15)

The canonical probability measure of particles in the configuration subspace has the form

$$\begin{array}{cc}\hfill {\mu }_N(T)& ={\mathrm{\Theta }}_{ϵ}^{-1}exp(-\frac{1}{2}{\displaystyle \sum _{j<i}{\mathrm{\Phi }}_{ϵ}(|{\mathit{x}}_j-{\mathit{x}}_i|)/T)\mathrm{d}{\omega }^N,}\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill \mathrm{d}{\omega }^N& \equiv {\displaystyle \prod _k^{N}exp(-m\mathrm{\Phi }({\mathit{x}}_k)/T)\mathrm{d}{\mathit{x}}_k,}\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill {\mathrm{\Theta }}_{ϵ}& ={\displaystyle \int {exp}_{\omega }(-\frac{1}{2}\sum _{j,i;j\ne i}{\mathrm{\Phi }}_{ϵ}(|{\mathit{x}}_j-{\mathit{x}}_i|)/T)\mathrm{d}{\omega }^N,}\hfill \end{array}$$(18)

and the divergence of the configuration integral Θ_ϵ → ∞ at ϵ → 0 yields |x_j − x_i|→0. Separating two particles from the N-particle ensemble (x_N ≡ x, x_N−1 ≡ x′), we can consider the density of the probability measure ${\mathfrak{P}}_N(\mathit{x},{\mathit{x}}^{\prime },{\omega }_{N-2};ϵ)$:

$$\begin{array}{cc}& {\mathfrak{P}}_N(\mathit{x},{\mathit{x}}^{\prime },{\omega }_{N-2};ϵ)=\hfill \\ \hfill & {\mathrm{\Theta }}_{ϵ}^{-1}(N,T)exp(-{\mathrm{\Phi }}_{ϵ}(|\mathit{x}-{\mathit{x}}^{\prime }|))\mathcal{P}(\mathit{x},{\mathit{x}}^{\prime },{\omega }_{N-2};ϵ),\hfill \end{array}$$(19)

where 𝒫(…) is a positive function that includes the interaction potentials of other particles, and ${\mathfrak{P}}_N(\mathit{x},{\mathit{x}}^{\prime },{\omega }_{N-2};ϵ\to 0)\to 0$ (x ≠ x′). So, $w^{\ast }-{lim}_{ϵ\to 0}{\mathfrak{P}}_N(\mathit{x},{\mathit{x}}^{\prime },{\omega }_{N-2};ϵ)={\mathfrak{P}}_{N-2}(\mathit{r})\delta (\mathit{x}-{\mathit{x}}^{\prime })$. Then, we can consider the canonical ensemble in the mean field limit for de-singularized potentials. According to the Hewitt–Savage representation theorem on the permutation-invariant probability measure (Hewitt & Savage 1955), any such (admissible) measure, μ, of the configuration subspace of the phase space of the many-particle system under consideration can be represented as an integral over d𝜚=ρ(x)d³x, ρ(x) and can be considered in the classical density form. The free energy functional for our system in this case has the form

$$\begin{array}{cc}\hfill \mathcal{F}(\rho )& =\frac{1}{2}{\displaystyle {\int }_{\omega }{\int }_{{\omega }^{\prime }}\rho (\mathit{x})\rho ({\mathit{x}}^{\prime })N{\mathrm{\Phi }}_{ϵ}(|\mathit{x}-{\mathit{x}}^{\prime }|)\mathrm{d}\mathit{x}\mathrm{d}{\mathit{x}}^{\prime }}\hfill \\ \hfill & +{\displaystyle {\int }_{\omega }\frac{1}{2}\rho (\mathit{x})\zeta (\mathit{x})\mathrm{d}\mathit{x}+{\int }_{\omega }\rho (\mathit{x})ln(\rho (\mathit{x}))\mathrm{d}\mathit{x},}\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill \zeta (\mathit{x})& =m{\mathrm{\Phi }}_{\mathrm{ext}}(\mathit{x})+B(\partial \omega ,\mathit{x}).\hfill \end{array}$$(21)

The minimization of ℱ(ρ) occurs with the functions

$$\begin{array}{cc}\hfill \rho (\mathit{x})& ={\mathfrak{R}}^{-1}·exp(T^{-1}({\displaystyle {\int }_{{\omega }^{\prime }}\rho ({\mathit{x}}^{\prime })N{\mathrm{\Phi }}_{ϵ}(|\mathit{x}-{\mathit{x}}^{\prime }|)\mathrm{d}{\mathit{x}}^{\prime }+\zeta (\mathbf{x}))),}\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill \mathfrak{R}& \equiv {\displaystyle {\int }_{{\omega }^{\prime }}exp(T^{-1}({\int }_{{\omega }^{\prime }}\rho ({\mathit{x}}^{\prime })N{\mathrm{\Phi }}_{ϵ}(|{\mathit{x}}^{″}-{\mathit{x}}^{\prime }|)\mathrm{d}{\mathit{x}}^{\prime }+\zeta ({\mathit{x}}^{″})))\mathrm{d}{\mathit{x}}^{″}.}\hfill \end{array}$$(23)

Since the most probable state of a thermodynamic system corresponds to the minimum free energy, the solution of Eq. (23) corresponds to the equilibrium density distribution in the multi-particle system.

If we recall the original Vlasov-Poisson equation, Eq. (2), the simplest stationary solution of Eq. (23) is the function $f(\mathit{x},\mathit{v})=\widehat{f}(\chi )$, where χ = mv²/2 + Φ(x) (energy substitution) and $\widehat{f}(\chi )\in C^2(\mathbb{R})$. The Poisson equation when choosing $\widehat{f}(\chi )=c_{1}exp(-\chi /T)$ (the Maxwell-Boltzmann equilibrium distribution known from collisional theory) can be represented as

$$\begin{array}{cc}& {\mathrm{\Delta }}^{(d)}\mathrm{\Phi }=G{\varkappa }_d(T)exp(-\mathrm{\Phi }/T)-c^2\mathrm{\Lambda },\hfill \end{array}$$(24)

$$\begin{array}{cc}& {\varkappa }_d(T)=c_0{|s_d|}^2{\displaystyle {\int }_{{\mathbb{R}}_{+}^1}exp(-{\eta }^2/(2T)){\eta }^{d-1}\mathrm{d}\eta ,}\hfill \end{array}$$(25)

where |s_d| is the area of the d-dimensional sphere of unit radius. In this case, the potential Φ(x) refers to the averaged self-consistent field determined by integrating into formula (22) the local densities of particles with a nucleus in the form of an inter-particle potential, Φ_ϵ(x − x′). Accordingly, the density of particles, as noted in Sect. 2, can be considered as a continuous argument function, and with it we can pass to the weak topology of isolated points. The equivalence of temperatures for the canonical and micro-canonical descriptions, generally speaking, is not too obvious, but significant problems can only arise for strongly nonequilibrium systems with negative thermodynamic temperatures. The influence of system boundaries (Eq. (22)), in turn, is extremely important. We will deal with this issue in the future, using a more appropriate formalism.

The form of the LGE for the gravitational potential in cosmological systems is essentially local (as it is a differential equation), and it is difficult to use it to describe global processes in a self-consistent field in order to analyze the evolution of the entire system of particles (substructures). When studying nonlocal effects, it is expedient to consider the integral version of the equation for the potential. However, a natural limitation for the use of this form in analytical and numerical calculations is the obvious need to take into account the boundaries of the region that contains the system of particles with regards to the value of the potential at a given point (i.e., the explicit form and effects of the term ζ(x) in Eq. (22)). Of course, this is due to the fact that, for a system of particles interacting according to the law of Eq. (1), the zero Dirichlet condition, Φ(x|_∂ω) = 0, has to apply to the boundary of the region under study. This condition is subtle, for example, for the general, though unrealistic, case of an asymmetric many-particle system.

We considered the following physical problem regarding the structure of the local (late) Universe where there are high density regions, consequences of cosmological fluctuations at earlier phases: can the gravitational interaction between these ω_j regions (each of which contains, e.g., a hierarchy of substructures of various sizes) or in the vicinity of each region form a (quasi-)stationary ordered structure? The essential difference compared to a charged plasma with Debye screening, which allows pseudo-homogeneous equilibrium particle distributions, is that the cosmological substructures evolving from primary perturbations of density macroscale objects are spherically symmetric (d = 3) or radially symmetric (for d = 2), according to the Gidas–Ni–Nirenberg theorem (Gidas et al. 1979). In accordance with the structure of the inter-particle potential, ${\mathfrak{Y}}_d$, in the neighborhood spheres with increased density, a radially symmetric layer arises, in which the attraction to the center of the sphere interplays with a semi-infinite spherical layer with a prevailing repulsion. For a system of two spheres that contain matter of increased density, there is competition between those two centers, leading either to the destruction of one of them or, in connection with the cosmological expansion, to the emergence of repulsion zones between them.

Two approaches can be considered when searching for the gravitational potential in the vicinity of a region that contains a system of interacting particles.

The first is a solution of the equation for the potential (inhomogeneous LGE) under given Dirichlet conditions on some a priori defined boundary of the region, ∂ω_j; this boundary should include only one local overdensity. In this case, one can consider both the internal and external Dirichlet problems, which correspond to attempts at establishing the values of the potential (and, hence, the density of particles) inside and outside the ∂ω_j boundary shell, respectively. Influence potentials of neighboring regions can be neglected under certain conditions. The data at the boundaries must be consistent with the resulting solution (i.e., they must be a continuity of the solution).

A second, similar situation can also be considered for the Neumann problem; however, it is possible to define the boundary of the ∂ω_j domain due to the presence of a maximum of the interaction potential (at this point, the derivative of the potential vanishes). However, in such a formulation, a question naturally arises regarding the legitimacy of describing the dynamics of the region of the particles only as the zone of attraction of the potential.

The physical aspects of this problem obviously determine the way additional conditions for the nonlinear potential equation are set. We restricted ourselves to using the Dirichlet conditions (due to the shell theorem on the equivalence of the gravitational field of a sphere and a point at its center). Thus, we are interested in whether inside (or outside) the fixed region, ω, for a given value of the potential at the boundary, there are secondary solutions of the LGE, Eq. (24), that possess a (quasi-)periodicity property.

The solution of the Dirichlet problem for the Poisson equation, according to Gilbarg & Trudinger (1983) and Sauvigny (2006), can be represented through the Green’s elliptic operator and the value of the potential at the borders. Explicitly, this solution is

$$\begin{array}{cc}\hfill \mathrm{\Phi }(\mathit{x})& =-{\displaystyle {\int }_{\omega }\mathcal{G}(\mathit{x},{\mathit{x}}^{\prime })(G\rho ({\mathit{x}}^{\prime })-\frac{c^2\mathrm{\Lambda }}{4\pi })\mathrm{d}{\mathit{x}}^{\prime }}\hfill \\ \hfill & -\frac{1}{4\pi }{\displaystyle {\int }_{\partial {\omega }^{\prime }}\mathrm{\Phi }{{|}_{\partial \omega }\mathit{n}·\frac{\partial }{\partial {\mathit{x}}^{\prime }}\mathcal{G}(\mathit{x},{\mathit{x}}^{\prime })|}_{\partial \omega }\mathrm{d}S\prime ,}\hfill \end{array}$$(26)

where 𝒢(x, x′) is Green’s function of the Dirichlet problem for the inhomogeneous Poisson equation. If we choose a ball of radius ℛ ($\omega =\{\stackrel{\sim }{\mathit{x}};0\le |\stackrel{\sim }{\mathit{x}}|\le \mathcal{R}\}$), then it is possible, if the distribution of particles is close to the isotropic one, to assume the values of the potential at its boundary, in accordance with the theorem in Gurzadyan (1985), as

$$\begin{array}{c}\hfill {\mathrm{\Phi }|}_{\partial \omega }=-GM/\mathcal{R}-c^2\mathrm{\Lambda }{\mathcal{R}}^2/6,(M=Nm).\end{array}$$(27)

Replacing ρ(x) = (4π)⁻¹𝜘_d(T)Gexp( − Φ/T), and defining, according to the Nowakowski (2001), the Green’s function for the considered Dirichlet domain,

$$\begin{array}{cc}& \mathcal{G}(\mathit{x},{\mathit{x}}^{\prime })\equiv 4\pi {\displaystyle \sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }\frac{Y_{\ell m}^{\ast }(\theta \prime ,\phi \prime )Y_{\ell m}(\theta ,\phi )}{2\ell +1}\frac{x_{<}^{\ell }{x}_{>}^{\ell }}{{\mathcal{R}}^{2\ell +1}},}\hfill \\ \hfill & x_{<}=\min (|\mathit{x}|,|{\mathit{x}}^{\prime }|),x_{>}=\max (|\mathit{x}|,|{\mathit{x}}^{\prime }|),\hfill \end{array}$$(28)

we obtain a nonlinear integral equation for the potential Φ(x) (for d = 3):

$$\begin{array}{cc}\hfill \mathrm{\Phi }(\mathit{x})& =-{\varkappa }_d(T)G{\displaystyle {\int }_{{\omega }^{\prime }}(\frac{1}{|\mathit{x}-{\mathit{x}}^{\prime }|}-\mathcal{G}(\mathit{x},{\mathit{x}}^{\prime }))}\hfill \\ \hfill & exp(-\mathrm{\Phi }({\mathit{x}}^{\prime })/T)\mathrm{d}{\mathit{x}}^{\prime }-\frac{c^2\mathrm{\Lambda }}{6}{\mathit{x}}^2+C_0,\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill C_0& =-\frac{\mathit{GM}}{\mathcal{R}}-\frac{c^2\mathrm{\Lambda }{\mathcal{R}}^2}{6}.\hfill \end{array}$$(30)

After obvious transformations, this can be rewritten as an inhomogeneous Hammerstein-type equation for the dimensionless potential, U_H:

$$\begin{array}{cc}\hfill U_{\mathrm{H}}(\mathit{x})& =\widehat{\mathfrak{W}}(U_{\mathrm{H}}(\mathit{x})),\widehat{\mathfrak{W}}(U_{\mathrm{H}}(\mathit{x}))\hfill \\ \hfill & \equiv {\lambda }_{\mathrm{H}}(T){\displaystyle {\int }_{{\omega }^{\prime }}\underset{\mathcal{K}(\mathit{x},{\mathit{x}}^{\prime })}{\underbrace{(\frac{1}{|\mathit{x}-{\mathit{x}}^{\prime }|}-\mathcal{G}(\mathit{x},{\mathit{x}}^{\prime }))}}exp(-U_{\mathrm{H}}({\mathit{x}}^{\prime }))\mathrm{d}{\mathit{x}}^{\prime }+}\hfill \\ \hfill & +{\alpha (\mathrm{\Lambda },T)|\mathit{x}|}^2;\hfill \end{array}$$(31)

$$\begin{array}{cc}\hfill {\lambda }_{\mathrm{H}}(T)& \equiv \frac{-{\varkappa }_3(T)G_3}{T}exp(-\frac{C_0}{T}),\alpha (\mathrm{\Lambda },T)=-\frac{c^2\mathrm{\Lambda }}{6T},\hfill \\ \hfill & U_{\mathrm{H}}=\frac{\mathrm{\Phi }-C_0}{T}.\hfill \end{array}$$(32)

It should be noted that, for a spherically symmetric density distribution,

$$\begin{array}{c}\hfill {\displaystyle {\int }_{{\omega }^{\prime }}\mathcal{G}(\mathit{x},{\mathit{x}}^{\prime })\rho ({\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }\to C_1(=\mathrm{const}).}\end{array}$$(33)

4. Formation of voids

Equation (31) contains information about the behavior of the considered cosmological, non-GR dynamics of a many-particle system with the gravitational interaction of Eq. (1). First of all, we should point out the obvious absence of a solution for this equation for d = 3 in the entire area of study (and in the whole space, except for two-dimensional surface isograves) in the form of a constant potential, in contrast to the classical Newtonian case of only attractive gravity. So, in this case, the global uniform distribution of matter, as a condition for matching the modified Poisson equation and the Vlasov equation, can only be considered as a certain approximation; a more comprehensive analysis is needed for more general cases.

We considered a linearized (near the point U_H(0) = U₀) version of Eq. (31) in the homogeneous and inhomogeneous cases:

$$\begin{array}{c}\hfill U^{†}={\widehat{\mathfrak{W}}}_{0^{\prime }}(U^{†}),U^{‡}={\widehat{\mathfrak{W}}}_{0^{\prime }}(U^{‡})+\alpha {|\mathit{x}|}^2.\end{array}$$(34)

Operator $\widehat{I}-{\widehat{\mathfrak{W}}}_{0\prime }$, where ${\widehat{\mathfrak{W}}}_{0\prime }\equiv \widehat{\mathfrak{W}}\prime [U_0]$ is a Frechet derivative of the Hammerstein integral operator on the right side of the first equation, belongs to the class of zero-index Noetherian operators with a weak singularity. Since Green’s function is symmetric to its coordinates, for matter distributions that slightly deviate from spherically symmetry, ${\widehat{\mathfrak{W}}}_{0\prime }$ can also be considered self-adjoint. For a small deviation in the amplitude of δU^† from the selected solution, one can apply the well-known mathematical apparatus for the analysis of Fredholm operators to the homogeneous linear equation $\delta U^{†}-\lambda {\int }_{{\omega }^{\prime }}\mathcal{K}(\mathit{x}-{\mathit{x}}^{\prime })\exp (-U_0^{†})\delta U^{†}\text{d}{\mathit{x}}^{\prime }=0$. To simplify the calculations without losing generality, we can take $U_0^{†}=0$ and look for periodic solutions of the last equation in the form of an expansion in terms of eigenfunctions 𝔟_j = c_ωexp(iqx) of the kernel 𝒦 in $\omega \subseteq {\mathbb{R}}^3$) in the form δU^† = ∑_j𝔞_j(c_ω)exp(iqx) (q_{ℓ = 1, 2, 3} = 2π/d; the cubic case can be generalized to d_ℓ ≠ d_k). Substituting this expression into the first population equation, Eq. (34), we have

$$\begin{array}{c}\hfill 1=\lambda {\displaystyle {\int }_{{\omega }^{\prime }}\mathcal{K}(\mathit{x}-{\mathit{x}}^{\prime })exp(-i\mathit{q}(\mathit{x}-{\mathit{x}}^{\prime }))\mathrm{d}|\mathit{x}-{\mathit{x}}^{\prime }|.}\end{array}$$(35)

Here we introduce a critical value of the parameter λ = λ_c, corresponding to the case d → ∞, q_ℓ ≡ 0. Using it, we can write the criterion for the existence of periodic solutions for the linearized integral homogeneous Poisson equation:

$$\begin{array}{cc}\hfill \lambda & =({\displaystyle {\int }_{\omega }\mathcal{K}(r)\frac{sin(qr)}{\mathit{qr}}{r}^2\mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi {)}^{-1}\ge {\lambda }_{\mathrm{c}}}\hfill \\ \hfill & \equiv ({\displaystyle {\int }_{\omega }\mathcal{K}(r)r^{2}sin(\theta )\mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi {)}^{-1}.}\hfill \end{array}$$(36)

Obviously, this criterion is suitable only for the case Λ ≡ 0, and only for the two-body interaction approximation. This implies that the root, q, of Eq. (35) is unique and, thus, that the potential distribution will be purely periodic. As there are several incommensurable roots of q_s, the distribution of the potential will belong to the class of almost-periodic functions. The accounting of collective interactions of N particles leads to a (homogeneous) equation of the form

$$\begin{array}{cc}& \delta U^{†}(\mathit{x})={\displaystyle \sum _{k=1,\dots ,N}{\int }_{{\omega }_1}\dots {\int }_{{\omega }_k}{\mathcal{K}}_k(\mathit{x},{\mathit{x}}_1,\dots ,{\mathit{x}}_k)}\hfill \\ \hfill & (exp(-\delta U^{†}({\mathit{x}}_1))\dots exp(-\delta U^{†}({\mathit{x}}_k)))\mathrm{d}{\mathit{x}}_1\dots \mathrm{d}{\mathit{x}}_k.\hfill \end{array}$$(37)

After linearizing this equation, we get (see Vlasov 1961)

$$\begin{array}{cc}\hfill \delta U^{†}(\mathit{x})& ={\displaystyle \sum _k{\lambda }_k{\int }_{{\omega }_1}\dots {\int }_{{\omega }_k}{\mathcal{K}}_k(\mathit{x},{bfx}_1,\dots ,{\mathit{x}}_k)}\hfill \\ \hfill & {\displaystyle \sum _{\ell }^k\delta U^{†}({\mathit{x}}_{\ell })\prod _{s=1}^k\mathrm{d}{\mathit{x}}_s.}\hfill \end{array}$$(38)

The corresponding criterion, criterion (4), for the occurrence of three-dimensional periodic solutions can be represented in the following form:

$$\begin{array}{cc}& {\displaystyle {\int }_{\omega }\sum _k\sum _s^k{\lambda }_k\frac{sin(qr)}{\mathit{qr}}({\int }_{{\omega }_1}\dots {\int }_{{\omega }_{s-1}}{\int }_{{\omega }_{s+1}}\dots {\int }_{{\omega }_k}}\hfill \\ \hfill & {\mathcal{K}}_s(\mathit{x},{\mathit{x}}_1,\dots ,{\mathit{x}}_k)\frac{{\prod }_{n=1}^k\mathrm{d}{\mathit{x}}_n}{\mathrm{d}{\mathit{x}}_s})r^{2}sin(\theta )\mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi =1.\hfill \end{array}$$(39)

Therefore, the accounting of the cluster interactions in a many-particle system is consistent with the accounting of two-particle interactions: the appearance of a periodic structure potential and of the density of matter in a linear homogeneous approximation occurs abruptly upon reaching a certain, quasi-equilibrium temperature in the system. The periods of the structures are determined from conditions (35) and (39).

For an inhomogeneous linear equation (the second half of Eq. (34)), it is possible, using the Hilbert–Schmidt theorem (Tricomi 1957), to obtain an explicit form of the (unique) solutions in the form of a resolvent series, as uniformly and absolutely convergent Fourier series in the eigenfunctions of the kernel 𝒦, for λ ≠ λ_ℓ (ℓ = 1, 2, …):

$$\begin{array}{c}\hfill \delta U^{‡}(\mathit{x})=\lambda {\displaystyle {\int }_{\omega }\sum _{\ell =1}^{\infty }\frac{{\mathfrak{b}}_{\ell }(\mathit{x}){\mathfrak{b}}_{\ell }({\mathit{x}}^{\prime })}{{\lambda }_{\ell }-\lambda }\alpha |{\mathit{x}}^{\prime }{|}^2\mathrm{d}{\mathit{x}}^{\prime }+\alpha {|\mathit{x}|}^2,}\end{array}$$(40)

where λ_ℓ are the characteristic numbers of the homogeneous equation (i.e., they correspond to the eigenfunctions 𝔟_ℓ(x)).

Thus, the periodicity of the structure, when Λ repulsion is taken into account, degenerates (in the simplest case) into a composition of sinusoidal functions and growing branches of parabolas, which leads to a smoothing of the periodicity and the dominance of repulsion at a certain distance, in the line of interactions between two formed inhomogeneities in a pseudo-homogeneous distribution of matter. This process can be considered a pair interaction and a mechanism for the formation of voids in an initially uniformly distributed material continuum. The anisotropy of particle momentum directions in the channels between I and II inhomogeneities forms not a cubic lattice (d_1, 2, 3 = d) but quasi-one-dimensional layers (d₃ ≪ d_1, 2). Their velocity profile in these channels is decelerated for initially fast particles and thus is synchronized as a uniform distribution in velocity space, due to repulsion in the far zone of the second component of the macro-system (i.e., the II inhomogeneity at the other end of the channel). It can be assumed that a significant part of the density inhomogeneity transforms into a flat structure corresponding to the first potential minimum and that this situation is mirrored in the channel of pair interaction. The repulsion acts as an external force, and from the II side leads to the quasi-stationarity of the wall formed by the I inhomogeneity and vice versa. The Λ repulsion leads to initially attracting massive inhomogeneities appearing at distances at which the repulsion becomes more significant due to the influence of the cosmological term.

As shown in Gurzadyan & Stepanian (2021a,b), based on an analysis of observational data, the repulsive term in Eq. (1) can dominate the attraction term depending on the mean matter density of certain galaxy clusters. Then, based on the above analysis, one can conclude that, although the considered mechanism of the void formation is common – the mutual interplay of the self-consistent gravitational attraction and of the repulsion due to the cosmological constant term – the voids can possess different mean scales that are determined by the local conditions, the density, and the initial momenta of particle flows. Thus, our proposed method for analyzing the solutions to the Poisson equation is able to explain the presence of voids of various scales.

The solution of the nonlinear Hammerstein equation for the potential differs significantly from the solution of the Fredholm equations. If first we turn to the case Λ ≡ 0 (homogeneity of the equation), we can immediately see that Eq. (31) has a constant solution, $U_{\mathrm{H}}=U_{\mathrm{H}}^{(0)}=\mathrm{const}$. From a physical point of view, this means that we are considering a medium that consists of particles, the interaction between which corresponds to only their mutual attraction, in which there are no primary perturbations. Obviously, this is an extremely unstable system. But even without an external influence leading to the local coalescence of particles, when the parameter λ reaches a certain value, solutions of a new type arise, branching off from the constant.

We denote ${\lambda }_{\mathrm{H}}={\lambda }_{\mathrm{H}}^{(0)}+\eta$, $U_{\mathrm{H}}(\mathit{x})=U_{\mathrm{H}}^{(0)}+u(\mathit{x})$. Then, if we expand the exponent in the integrand expression in a Taylor series, following Bratu (1914), we get

$$\begin{array}{cc}& u(\mathit{x})-{\lambda }_{\mathrm{H}}^{(0)}exp(-U_{\mathrm{H}}^{(0)}){\displaystyle {\int }_{{\omega }^{\prime }}\mathcal{K}(\mathit{x}-{\mathit{x}}^{\prime })u({\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }}\hfill \\ \hfill & =\eta exp(-U_{\mathrm{H}}^{(0)}){\displaystyle {\int }_{{\omega }^{\prime }}\mathcal{K}(\mathit{x}-{\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }-\eta exp(-U_{\mathrm{H}}^{(0)})}\hfill \\ \hfill & {\displaystyle {\int }_{{\omega }^{\prime }}\mathcal{K}(\mathit{x}-{\mathit{x}}^{\prime })u({\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }}\hfill \\ \hfill & +(\eta +{\lambda }_{\mathrm{H}}^{(0)})exp(-U_{\mathrm{H}}^{(0)}){\displaystyle \sum _{j=2}^{\infty }{(-1)}^j{(j!)}^{-1}{\int }_{{\omega }^{\prime }}\mathcal{K}(\mathit{x}-{\mathit{x}}^{\prime })(u({\mathit{x}}^{\prime }){)}^j\mathrm{d}{\mathit{x}}^{\prime }.}\hfill \end{array}$$(41)

Substituting $u(\mathit{x})=\sum _{s=1}^{\infty }{(\eta exp(-U_{\mathrm{H}}^{(0)}))}^{s/2}\mathcal{X}(\mathit{x})$ (here 𝒳(x) are unknown functions to be defined), we obtain a sequence of linking equations with a linear left-hand side of the form

$$\begin{array}{cc}& \widehat{\mathfrak{N}}({\mathcal{D}}_1)={\displaystyle {\int }_{{\omega }^{\prime }}\mathcal{K}(\mathit{x}-{\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }=0,\widehat{\mathfrak{N}}({\mathcal{D}}_2)={\int }_{{\omega }^{\prime }}\mathcal{K}(\mathit{x}-{\mathit{x}}^{\prime })}\hfill \\ \hfill & (1+{\lambda }_{\mathrm{H}}^{(0)}exp(-U_{\mathrm{H}}^{(0)})\frac{{({\mathcal{D}}_1({\mathit{x}}^{\prime }))}^2}{2!})\mathrm{d}{\mathit{x}}^{\prime },\dots ,\hfill \end{array}$$(42)

$$\begin{array}{cc}& \widehat{\mathfrak{N}}({\mathcal{D}}_j)\equiv {\mathcal{D}}_j+{\lambda }^{(0)}exp(-U^{(0)}){\displaystyle \int }\mathcal{K}(\mathit{x}-{\mathit{x}}^{\prime }){\mathcal{D}}_j({\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }\hfill \end{array}$$(43)

(the operator associated with the homogeneous Fredholm equation). We have already obtained the form of the first term of a series of successive approximations: 𝒟₁(x) = c_ω ⋅ sin(qx). The equation $\widehat{\mathfrak{N}}({\mathcal{D}}_1)=0$ has periodic solutions under the criterion λ⁽⁰⁾exp(−U⁽⁰⁾Ω(q)+1 = 0, Ω(q)≡4π∫𝒦(r)(sin(qr)/(qr))r²dr). In accordance with Vlasov (1961), we successively solved the above equations and obtain a solution for u(x):

$$\begin{array}{cc}\hfill u(\mathit{x})& ={\displaystyle \sum _{m=0}(\eta exp(-U^{(0)}){)}^{m+1/2}({\sum _{\ell =1}^m}_{(I)}{C}_{2\ell +1}^{(2m+1)}}\hfill \\ \hfill & sin((2\ell +1)\mathit{q}\mathit{x}){+}_{(III)}{C}_{2m+1}sin(\mathit{q}\mathit{x}))\hfill \\ \hfill & +{\displaystyle \sum _{m=1}^{\infty }(\eta exp(-U^{(0)}){)}^m({\sum _{\ell =1}^m}_{(II)}{C}_{2\ell }^{(2m)}cos(2m\mathit{q}\mathit{x})}\hfill \\ \hfill & {+}_{(III)}{C}_{2m}sin(\mathit{q}\mathit{x})).\hfill \end{array}$$(44)

The resulting series, subject to its convergence, is a periodic function with T = T(|q|). Since all coefficients _(i)C_… can be obtained explicitly by sequential calculation, we are able to determine the numerical values of the amplitudes of harmonics in a given Fourier series. Thus, one can obtain a “fine structure” of the potential, which means that the density of particles in the emerging structures branching off from the cosmological solution with a constant density parameter, η, is related to the difference between the temperature of the medium and the critical temperature corresponding to the critical parameter, λ_c.

As already mentioned, this technique requires knowledge of the “primary” solution to the homogeneous Hammerstein equation for the potential. Using the methods of the perturbation theory of Fredholm operators (e.g., Keller & Langford 1972; Krasnosel’sky 1964), one can construct a solution for the inhomogeneous Hammerstein equation based on the results obtained above. If we denote ${\widehat{I}}_{+}\equiv \widehat{I}-\alpha {|\mathit{x}|}^2$, then ${\widehat{I}}_{+}-{\widehat{\mathfrak{W}}}_{0\prime }$ is a Noetherian operator zero index, dim𝔜 = 1, $\mathfrak{Y}\equiv \mathrm{Ker}({\widehat{I}}_{+}-{\widehat{\mathfrak{W}}}_{0\prime })$. The condition that ${\widehat{I}}_{+}-{\widehat{\mathfrak{W}}}_{0\prime }$ has a pseudo-inverse – that is, a bounded inverse from the co-kernel complement 𝔜^*⊥ in some complement $\overline{\mathfrak{Y}}$ of the kernel 𝔜 – is equivalent to the statement about the existence of a unique element $v\in \overline{\mathfrak{Y}}$, such that $({\widehat{I}}_{+}-{\widehat{\mathfrak{W}}}_{0\prime })v=w$, w ∈ 𝔜^*⊥. In this case, the nonlinear Hammerstein equation has a unique solution: (λ_H(ε),u_H(ε)), where ε is a small parameter associated with the temperature deviation from the critical one. The values λ_H and u_H are the limits of the sequences λ^μ + 1(ε) = F₊(u^μ(ε)), u^μ + 1(ε) = ε(ϕ₀ + εv^μ + 1(ε)), and $v^{\mu +1}\in \overline{\mathfrak{Y}}$ (u⁰(ε) = εϕ₀, v⁰(ε) = 0, where ϕ₀ is an element ${\widehat{I}}_{+}-{\widehat{\mathfrak{W}}}_{0\prime }$ kernels). Thus, we have at our disposal a procedure for obtaining a (unique) solution of an inhomogeneous nonlinear equation for the potential. Near it, we can consider a certain neighborhood where new non-holomorphic branches of its solution can appear. However, it is already clear that the Hammerstein equation under study has a solution that is locally close in its construction to the inhomogeneous Fredholm equation. Further obtaining a solution constructed in a consistent way outside the uniqueness-domain branching solutions (replacing the constant solution for the homogeneous case) is practically equivalent to the homogeneous case.

Thus, it can be argued that the behavior of the potential in a many-particle system is capable of creating conditions for the emergence of two-dimensional structures, which can be associated with the walls of voids. The inclusion in the Poisson equation of an additional repulsive Λ term from Eq. (1) plays a fundamental role here. Moreover, using the formalism developed above, one can obtain the fine structure of the voids themselves, since the solutions of the Hammerstein equation for the potential have a very complex multi-periodic structure, which can be used for comparison with observational data.

5. Conclusions

The emergence of two-dimensional structures such as Zeldovich pancakes is associated with the density perturbations described by classical or (weakly) relativistic hydrodynamics. Recent observational tensions, such as the Hubble constant tension, can require the consideration of nonidentical processes that lead to structure formation and evolution on local and global cosmological scales.

Thus, we have considered a description of the local Universe, taking into account (i) the self-consistent many-particle interaction by means of the Vlasov kinetic formalism and (ii) modified gravity with the cosmological constant term from Eq. (1), that is to say, the repulsion at galaxy cluster scales. This description predicts two Hubble flows, local and global ones, with nonidentical Hubble parameters. Thus, in a certain sense, the Hubble tension can be considered as empirical support of Eq. (1) and of the presented mechanism of structure formation.

Based on these aspects, we have developed a kinetic model for the occurrence of voids separated by two-dimensional surfaces, following from a rigorous analysis of the Poisson equation and its quasi-oscillatory solutions. We show the appearance mechanism of voids of different scales, depending on the local region of the Universe.

Its appearance in the Poisson equation with the modified potential following from the theorem in Gurzadyan (1985) on a general function for the “sphere-point” identity poses a mathematical problem of the study of inhomogeneous Fredholm integral equations.

The principal consequence of the developed mechanism is that the cosmological constant poses a natural scaling for the voids, via the characteristic distance when the second, repulsive Λ term dominates over the first term in Eq. (1):

$$\begin{array}{c}\hfill r_{\mathrm{cr}}={\left(\frac{3GM}{\mathrm{\Lambda }{c}^2}\right)}^{1/3}.\end{array}$$(45)

This distance, r_cr, together with the total mass within the corresponding volume, or, equivalently, with the mean density in that volume, will determine the actual size of the voids (Gurzadyan & Stepanian 2021a).

We can illustrate this quantitatively with the available observational parameters of the Virgo Supercluster and the Laniakea Supercluster. For the Virgo Supercluster, one has M = 1.48 × 10¹⁵ M_⊙ (Einasto et al. 2007; Reid et al. 2019), and the scale at which the repulsive Λ term can lead to a local Hubble flow is in the range 17.3 < r < 18.4 Mpc. For the Laniakea Supercluster, one has the mass M = 10¹⁷ M_⊙ (Tully et al. 2014) and the critical radius when the repulsion overwhelms the gravitational attraction, r_cr ≃ 51 Mpc. These parameters imply diameters of voids from 35 Mpc up to 100 Mpc in such environments.

These predictions can be directly tested in observational surveys, namely, via the search for correlations in the sizes of the voids versus the total masses and mean densities of their local regions (e.g., Ceccarelli et al. 2006; Gurzadyan et al. 2014; Samsonyan et al. 2021a,b). Such studies, along with the dynamics and the flows of galaxy clusters and superclusters, can be instrumental in revealing whether there are additional tensions in descriptions of the late and early Universe.

Acknowledgments

We are thankful to the referee for helpful comments.

References