© The Authors 2023

1. Introduction

The observational indications for the accelerated expansion of the Universe and the positive cosmological constant have essential influence on the theoretical studies of the various phases of the cosmological evolution. The Hubble tension (Verde et al. 2019; Riess 2020; Riess & Breuval 2022; Di Valentino et al. 2021; Dainotti et al. 2021, 2023), as the plausible tension between the late and early Universe, has also stipulated the consideration of theories and models (e.g., Capozziello & Lambiase 2022; Bouchà et al. 2022; Bajardi et al. 2022) addressing a broad spectrum of relevant issues, including the formation of cosmic structures.

Zeldovich (Zeldovich 1970; Shandarin & Zeldovich 1989; Shandarin & Sunyaev 2009) pancake theory predicts the formation of cosmic structures based on the evolution of initial density perturbations (Peebles 1993), with the essential use of the Lagrangian singularity theory. The pancake theory originated from cosmology thus linked the “hydrodynamics” of the Universe to symplectic geometry (Arnold et al. 1982; Arnold 1982). The observational facts on the non-zero cosmological constant or the possible tension of the late (local) and early (global) Universe would obviously also have an influence on structure formation and pancake theory.

The concept of two Hubble flows, the local and global one, with a non-identical Hubble constant, is among the suggested approaches to deal with the Hubble tension and the dark sector (Gurzadyan & Stepanian 2021a,b). That approach is based on the theorem proved in Gurzadyan (1985), positing that the general function for the force satisfying the sphere-point mass gravity’s identity takes the following form:

$$\begin{array}{c}\hfill F=-\frac{\mathit{GMm}}{r^2}+\frac{\mathrm{\Lambda }{c}^{2}mr}{3}.\end{array}$$(1)

This function does not satisfy the second statement of Newton’s shell theorem, namely, in this case the gravitational field inside a shell is no more force-free. This formula, that is, taking as the second term the cosmological constant in weak-field general relativity and McCrea–Milne non-relativistic cosmology McCrea & Milne (1934; see also Zeldovich 1981) enables the description of the dynamics of groups and clusters of galaxies, as well as other phenomena (Gurzadyan & Stepanian 2018, 2019, 2020; Gurzadyan 2019). We also note that there are observational indications (Kravtsov 2013) for the influence of the halo on the spiral galaxy structure, thus supporting the non-force-free field inside a shell predicted by Eq. (1).

Our present study is a continuation of previous works (Gurzadyan et al. 2022, 2023), based on the use of the Vlasov kinetic technique applied to gravitating systems and revealing the paths of the formation of the filamentary cosmic structures. The question of the ordered relaxation of perturbations is rather subtle in the system along the selected direction, such as for filaments, along two mutually perpendicular directions. It is clear that the local structure formation essentially differs from the hydrodynamics at large scales, especially regarding the averaging of a self-consistent gravitational field in a system of particles. We show that the deterministic way for emergence of large-scale structures, coherent in the sense of having a translational invariance of substructures, can be crucial in the Local Universe. This approach differs from “order out of chaos” scenario of the pancake theory at global cosmological scales but is complementary to it at local scales. The role of the repulsive second term of Eq. (1) is crucial and also determines the scales of structure formation and sizes of the voids, thus linking the predictions to observations (Böhringer et al. 2021, and references therein).

2. Integral form of the Poisson–Liouville–Gelfand equation for the system of gravitating particles

The system of Vlasov–Poisson equations for describing the cosmological dynamics in a system of N particles of identical masses m_i = m ≡ 1 can be represented in the following form. We note that we implicitly assume that the system is considered in the domain Ω ⊆ ℝ² space with a smooth boundary ∂Ω, where diam Ω ≤ ∞, namely, the size of the region can be continued up to infinity. This system is expressed as follows:

$$\begin{array}{cc}& \frac{\partial F(\mathit{x},\mathit{v},t)}{\partial t}+{\mathrm{div}}_{\mathit{x}}(\mathit{v}F)+\widehat{G}(F;F)=0,\hfill \\ \hfill & \widehat{G}(F;F)\equiv -({\mathrm{div}}_{\mathit{v}}F)({\mathrm{\nabla }}_{\mathit{x}}\mathrm{\Phi }(F)),\hfill \end{array}$$(2)

$$\begin{array}{cc}& {\mathrm{\Delta }}_{\mathit{x}}^{(3)}\mathrm{\Phi }(F){|}_{t=t_{\ast },\forall t_{\ast }\in \mathbb{T}}=A{S}_{3}G{\displaystyle \int }F(\mathit{x},\mathit{v},t_{\ast })\mathrm{d}\mathit{v}-\frac{c^2\mathrm{\Lambda }}{2},\hfill \\ \hfill & S_3\equiv \mathrm{meas}{\mathcal{S}}^2=4\pi ,{\mathcal{S}}^d=\{\mathit{x}\in {\mathbb{R}}^d,|\mathit{x}=1|\},\hfill \end{array}$$(3)

where F(x, v, t) is the distribution function of gravitationally interacting particles, A is a normalization factor for particle density, and t_* is a fixed instant of time. Equation (3) is the Poisson equation for the potential of Eq. (1), taking into account the repulsive cosmological term.

The third term in the r.h.s. of the kinetic equation (Eq. (2)) can be represented as a “source-like” form

$$\begin{array}{cc}\hfill \widehat{G}(F;F)& =\mathit{G}(F)\frac{\partial F}{\partial \mathit{v}},\mathit{G}(F)=-{\mathrm{\nabla }}_{\mathit{x}}\mathrm{\Phi }(F),\hfill \\ \hfill \mathrm{\Phi }(F)& =A{S}_{3}G{\displaystyle \int \int }{\mathfrak{K}}_3(\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 \\ \hfill & +\frac{\mathrm{\Lambda }{c}^2}{12}{|\mathit{x}|}^2+{\widehat{\mathfrak{B}}}_3(\mathit{x},{\mathit{x}}^{\prime }),{\mathfrak{K}}_3(\mathit{x}-{\mathit{x}}^{\prime })=-|\mathit{x}-{\mathit{x}}^{\prime }{|}^{-1},\hfill \end{array}$$(4)

where ${\widehat{\mathfrak{B}}}_3(\mathit{x},{\mathit{x}}^{\prime })$ is an operator term that takes into account border influence domain (for d = 2, the situation is similar, but 𝔎₂(x, x′) = ln |x − x′|, S₂ = 2). The Newtonian potential Φ_N(r) = − Gm/r increases on an interval of r ∈ (0, +∞) (Φ_N ∈ ( − ∞, 0)), while the potential of Eq. (1):

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathit{GN}}(r)\equiv -GM/r-\frac{1}{2}{c}^2\mathrm{\Lambda }{r}^2,\end{array}$$

has a maximum of

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathit{GN}}^{(\max )}=-\frac{1}{2}G(3M{c}^{2/3}){\mathrm{\Lambda }}^{1/3},\end{array}$$

at the critical radius

$$\begin{array}{c}\hfill r_{\mathrm{crit}}=(3GM/\mathrm{\Lambda }{c}^2{)}^{1/3},\end{array}$$(5)

namely, the case where the Λ-term starts to dominate in the l.h.s. of Eq. (1).

The radius, r_crit, determines the mutual contribution of the gravitational attraction and repulsive Λ-term. That is to say, at relevant scales, it defines the local Hubble flow given by the equation:

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

where H_local and ρ_local are the local values of the Hubble constant and the matter density, respectively (Gurzadyan & Stepanian 2021b). The critical radius, r_crit, also defines the scaling between the semi-periodic coherent structures (i.e., the voids) formed as a result of action of attraction and repulsion (Gurzadyan et al. 2022, 2023).

We go on to consider the stationary dynamics, thus we set F = F(x, v). The subsequent analysis mainly concerns the second equation of the system Eqs. (2) and (3), which is the Poisson equation relative to the potential with no explicit time dependence. That is why varying the Hilbert–Einstein–Poisson action (Vedenyapin et al. 2021, 2011) includes a separate variation of fields (for a fixed particle distribution) and the distribution function. Thus, the considered below approach is applicable for adiabatic processes at quasi-equilibrium (weakly varying) particle distribution functions. In this case, for the distribution function, we can use the energy substitution (Vedenyapin et al. 2011), namely: $F(\mathit{x},\mathit{v})=f(\epsilon )\in C_{+}^1({ℝ}^1)$, where ε = v²/2 + Φ(x).

Thus, the particle density on the right side of the Poisson equation can be expressed in terms of the integral of the known equilibrium solution of Vlasov equations, having the most transparent physical meaning, namely, that of Maxwell–Boltzmann distributions f = f₀(ε) = ANexp(−ε/θ),

$$\begin{array}{cc}\hfill \mathrm{\Delta }\mathrm{\Phi }(\mathit{x})& =AN{G}_3{S}_3^2({\displaystyle {\int }_{y\in [0,\infty ]}exp(-y^2/(2\theta ))y^2\mathrm{d}y})·\hfill \\ \hfill & exp(-\mathrm{\Phi }/\theta )-\frac{c^2\mathrm{\Lambda }}{2},A,\theta ,R_{\mathrm{\Omega }}\in {\mathbb{R}}^1,\hfill \end{array}$$(6)

where R_Ω is the radius of the region Ω as a ball in accordance with the Gidas-Ni-Nirenberg theorem (Gidas et al. 1979; Dupaigne 2011). Here, θ has the meaning of dynamic kinetic temperature of the system of interacting particles. The meaning of the last value is by no means obvious, as thermodynamic equilibrium is globally absent here.

Introducing the (kinetic) temperature according to the general definition (Vlasov 1961, 1978), θ ≡ −(∂Φ/∂x_j)/ $(\partial \rho /\partial x_j){|}_{j=\overline{1,d}}$, we get its dependence on the local properties of the self-consistent potential. We note that the temperature, generally speaking, can be considered as a tensor quantity, but for simplicity, we confine ourselves to the scalar temperature in the above formula, implying summation over the repeated indices.

Thus, the Poisson equation (Eq. (3)) takes the form of an inhomogeneous Liouville–Gelfand (LG) equation (Gelfand 1963; Dupaigne 2011) with local (generalized) temperature changing the sign depending on the value of the derivative of the potential at a given point. As mentioned above, for two-particle problem – in particular, for a formal pair of a central coalescence of the main part of the particles and the conditional “far” particle (generalized Milne–McCrea model), a dominant one can become the repulsive force due to the presence of a quadratic term in the radius-vector containing cosmological parameter of Eq. (1). While the generalized indefinite thermodynamics of a system of gravitating particles becomes similar to that for the Onsager vortices in the classical hydrodynamics (Fimin & Chechetkin 2020) in the context of the formation of quasi-ordered large-scale coherent structures). The existence of solutions of the LG equation for large system sizes supports (unlike the case θ > ;0 for a non-positive Newton’s gravitational attractive potential) the existence of solutions to the Vlasov equation. In this case, the gaps in the solution due to the fact that for ${\mathrm{\Phi }}_r^{\prime }\gtrless 0$ solutions can exist on non-intersecting intervals (here variable r = |x|) on both sides of the maximum of the self-consistent potential, corresponding to Φ_GN for the two-particle potential). This can be shown using the parametric Young’s inequality (Pokhozhaev 2010). To do this, we multiply both sides of non-homogeneous LG equation to the expression q₁(Φ)q₂(x), where q₁ ∈ C¹(ℝ¹) is such that $Q_1\equiv -{(1+\exp (\mathrm{\Phi }/\theta ))}^{-1}={\int }_0^{\mathrm{\Phi }}{q}_1(z)\text{d}z$, and the test function q₂(x)∈C²(ℝ^d) can be represented as:

$$\begin{array}{cc}\hfill q_2(\mathit{x})& =q_2(\mathit{y}),\mathit{y}=\mathit{x}/R(0\le R\le R_{\mathrm{\Omega }}/2),\hfill \\ \hfill q_2(\mathit{y})& =1\mathrm{for}|\mathit{y}|\le 1,q_2(\mathit{y})=0\mathrm{for}|\mathit{y}|\ge 2,\hfill \\ \hfill q_0& \equiv {\displaystyle \int }(|\mathrm{\Delta }{q}_2{|}^2/q_2)\mathrm{d}\mathit{y}<\infty .\hfill \end{array}$$

Then, we integrate this and using the Young inequality obtain the expression:

$$\begin{array}{cc}& {\displaystyle {\int }_{\mathrm{\Omega }}{|\mathrm{\nabla }\mathrm{\Phi }|}^2\frac{exp(\mathrm{\Phi }/\theta )q_2(\mathit{x})(1-exp(\mathrm{\Phi }/\theta ))}{(1+exp(\mathrm{\Phi }/\theta ){)}^3{\theta }^2}\mathrm{d}\mathit{x}}\hfill \\ \hfill & \le -{\displaystyle {\int }_{\mathrm{\Omega }}\frac{\lambda }{2\theta (1+exp(\mathrm{\Phi }/\theta ){)}^2}\mathrm{d}\mathit{x}}\hfill \\ \hfill & +\frac{q_0\theta }{2R\lambda }-{\displaystyle {\int }_{\mathrm{\Omega }}\frac{c^2\mathrm{\Lambda }}{2}{q}_1(\mathrm{\Phi })q_2(\mathit{x})\mathrm{d}\mathit{x}\le -\frac{\lambda S_3}{2\theta }{R}^3+\frac{q_0\theta }{2\lambda R}+\frac{2c^2\mathrm{\Lambda }}{3\theta }{R}^3,}\hfill \\ \hfill & q_0\equiv {\displaystyle \int }({|\mathrm{\Delta }{q}_2|}^2/q_2)\mathrm{d}\mathit{y}<\infty ,\hfill \end{array}$$(7)

where

$$\begin{array}{c}\hfill \lambda \equiv AN{\gamma }_3{S}_3^{2}J(\theta ),J(\theta )\equiv {\displaystyle {\int }_0^{v_{\max }}exp(-v^2/(2\theta ))v^2\mathrm{d}v.}\end{array}$$

The positivity condition on the left-hand-side gives us the parameter link: for θ ≷ 0 the compatibility condition is c²Λ ≷ 3πλ.

Therefore, if there are appropriate relations between the parameters of the problem, the system of Vlasov–Poisson equations has solutions of the type of distribution functions that admit the energy substitution. The potentials of the gravitational field, which have the property of convexity (cf. the difference with tidal fields; Gurzadyan & Ozernoi 1981) in the general case, for an arbitrary R_Ω ≤ ∞ (in contrast to the case of the attraction potential), there is a limitation R_Ω < (q₀θ²/(4πλ²))^1/4.

Now we go on to examine in detail the properties of the gravitational potential and the influence of the cosmological term, as well as its influence at a given point of boundary conditions for region of space Ω. As already noted in Gurzadyan et al. (2022, 2023), in the formulation of the Dirichlet problem for the Poisson equation with a constant right-hand side on the boundary of the Ω region, according to averaging gravitational field outside the compact subdomain Ω₀ containing system of particles and located inside domain Ω (meas Ω₀ ≪ meas Ω), we can assume that the data on the ∂Ω boundary is given by in accordance with the theorem in Gurzadyan (1985). Then we consider what happens in the context of other similar equations.

We should first consider the properties of the boundary value problem for Eq. (6) in the following:

$$\begin{array}{c}\hfill \mathrm{\Delta }\varphi +\lambda exp(\varphi )=\frac{c^2\mathrm{\Lambda }}{2\theta },\varphi \equiv -\frac{\mathrm{\Phi }}{\theta },\varphi (\mathit{x}){|}_{\mathit{x}\in \partial \mathrm{\Omega }}={\varphi }_b(R_{\mathrm{\Omega }}).\end{array}$$(8)

If we assume a priori the quantity ϕ to be small in norm, we obtain from this a linear equation; in principle, a linearization near any norm-finite solution ϕ₀ homogeneous equation gives a similar result, although with the change λ → λexp(ϕ₀); for simplicity, we further assume ϕ₀ = 0): Δϕ + λϕ = c²Λ/(2θ). We note that the boundary condition at r → R_Ω must take the form: ∼k₁(θ)/R_Ω + k₂(θ)).

Thus, we have the Dirichlet problem for the inhomogeneous Helmholtz equation, the solution of which in a three-dimensional radially-symmetric case takes the form

$$\begin{array}{c}\hfill {\varphi }_3(\mathit{x})\sim r^{-1}(cos(\sqrt{\lambda }r)+sin(\sqrt{\lambda }r)){|}_{r=|\mathit{x}|}+c^2\mathrm{\Lambda }/(2\theta \lambda ),\end{array}$$(9)

which is analogous to the consideration of the inhomogeneous Poisson equation (Gurzadyan 1985) and an additional condition, ϕ(0) < ∞, can be satisfied by a discrete or locally finite density of mass carriers ρ, namely, $\exists \omega :\omega =\{0\le r_{-}\le R_{\min }\ll R_{\mathrm{\Omega }}\}:\rho (r_{-})=0$). We note that the Dirichlet condition is consistent with the fundamental system of solutions of the considered equation. For the stability of the solution of the boundary value problem, and to set the corresponding condition for the case λ ≡ 0 (for a non-homogeneous Poisson equation), the Dirichlet condition on ∂Ω is expressed as:

$$\begin{array}{c}\hfill {\varphi }_3(|\mathit{x}|\to R_{\mathrm{\Omega }})\sim R_{\mathrm{\Omega }}^{-1}+c^2\mathrm{\Lambda }{R}_{\mathrm{\Omega }}^2/(2\theta ).\end{array}$$(10)

Obviously, we can pose the asymptotic conditions at infinity for a linearized model reducing it to the Helmholtz equation, in accordance with the known Sommerfeld radiation conditions (Vladimirov 1971). At the same time, for the Poisson equation with the right-hand side, this is impossible, because lim_{r → ∞}|ϕ(r)| → ∞. However, in a formal sense, for the last equation, it is possible to associate solutions of exact hydrodynamic consequence of the system of Vlasov–Poisson equations with Milne–McCrea model. Hence, the latter model should be seen as a solution the Dirichlet problem in a ball with a large but finite radius R_Ω; moreover, the possibility of setting periodic boundary conditions from the point of view of physical problem statement is highly doubtful. The Helmholtz equation can be solved on a semi-infinite interval (moreover, periodic conditions may well be set), and it corresponds to non-relativistic cosmological model with pseudo-periodic decrease of density amplitude. We note, in fact, that these linear models are hierarchical steps for a more general models with a conditionally equilibrium distribution of particles in space.

For a more detailed analysis of the situation, it is expedient to pass to the integral representation of the equation for gravitational potential. The equation for the potential with Maxwell–Boltzmann particle distribution, corresponding to internal Dirichlet problem in a bounded domain Ω (under boundary conditions corresponding to the McCrea–Milne model), takes the following form:

$$\begin{array}{cc}\hfill \mathrm{\Phi }(\mathit{x})& ={\lambda }_I{\displaystyle {\int }_{\mathrm{\Omega }\prime }\mathcal{K}(\mathit{x},{\mathit{x}}^{\prime })exp(-\mathrm{\Phi }({\mathit{x}}^{\prime })/\theta )\mathrm{d}{\mathit{x}}^{\prime }-\frac{c^2\mathrm{\Lambda }}{12}{\mathit{x}}^2+C_0,}\hfill \\ \hfill \mathcal{G}(\mathit{x},{\mathit{x}}^{\prime })& \equiv 4\pi {\displaystyle \sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }\frac{Y_{\ell m}^{\ast }({\vartheta }^{\prime },{\phi }^{\prime })Y_{\ell m}(\vartheta ,\phi )}{2\ell +1}\frac{x_{<}^{\ell }{x}_{>}^{\ell }}{R_{\mathrm{\Omega }}^{2\ell +1}},}\hfill \\ \hfill & x_{<}=\min (|\mathit{x}|,|{\mathit{x}}^{\prime }|),\hfill \\ \hfill & x_{>}=\max (|\mathit{x}|,|{\mathit{x}}^{\prime }|),\hfill \\ \hfill \mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)& \equiv \mathcal{G}(\mathit{x},{\mathit{x}}^{\prime })-\frac{1}{|\mathit{x}-{\mathit{x}}^{\prime }|},C_0=-\frac{{\gamma }_{3}Nm}{R_{\mathrm{\Omega }}}-\frac{c^2\mathrm{\Lambda }{R_{\mathrm{\Omega }}}^2}{12},\hfill \\ \hfill & {\lambda }_I=\lambda /S_3.\hfill \end{array}$$(11)

In essence, the above is the explicit form of the equation for the potential introduced in the expression (Eq. (3)), where 𝒢(x, x′) is the Green function for the inner boundary value problem in the domain Ω; in this case, due to symmetry of the latter, we have:

$$\begin{array}{c}\hfill {\displaystyle {\int }_{{\mathrm{\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}$$(12)

We go on to introduce a new variable U(x)≡(Φ(x)−C₀)/θ, for a symmetric problem, where C₁ will enter the constant C₀, the above equation can be written as the Hammerstein integral equation of the form:

$$\begin{array}{cc}\hfill U(\mathit{x})& =\stackrel{\sim }{\lambda }\widehat{\mathfrak{G}}(U)+\alpha (\theta ,\mathrm{\Lambda }){\mathit{x}}^2,\hfill \\ \hfill \widehat{\mathfrak{G}}(U)& \equiv {\displaystyle {\int }_{{\mathrm{\Omega }}^{\prime }}\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)exp(-U({\mathit{x}}^{\prime }))\mathrm{d}{\mathit{x}}^{\prime },}\hfill \\ \hfill \stackrel{\sim }{\lambda }& \equiv \stackrel{\sim }{{\lambda }_{\mathit{II}}}(\theta )\equiv \frac{{\lambda }_I}{\theta }exp(-C_0/\theta ),\alpha (\theta ,\mathrm{\Lambda })=-\frac{\mathrm{\Lambda }{c}^2}{12\theta }.\hfill \end{array}$$(13)

To clarify the basic properties of the solution of the potential equation, first, we assume that the value of the potential differs little from the constant value U₀, namely,U(x) = U₀ − ϕ(x), |ϕ|≪U₀; from a physical point of view, we consider the region near the inflection point of the potential, where for multi-particle systems there should be an equilibrium plateau forces of attraction and repulsion (unless, of course, the size of the Ω system is large enough). A linearization of the equation near (quasi)constant solution U₀ leads to the homogeneous Fredholm second kind equation

$$\begin{array}{cc}& \varphi (\mathit{x})+{\lambda }_{\mathit{II}}{\displaystyle {\int }_{{\mathrm{\Omega }}^{\prime }}\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)\varphi ({\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }=0,}\hfill \\ \hfill & {\lambda }_{\mathit{II}}=\stackrel{\sim }{{\lambda }_{\mathit{II}}}(\theta )exp(-U_0).\hfill \end{array}$$(14)

In addition to the trivial solution, this equation has a set of periodic (in the 3D space) solutions. We introduce eigenfunctions of symmetrical kernel 𝒦(|x − x′|), and consider representation of solution in the form of Fourier series

$$\begin{array}{c}\hfill \varphi (\mathit{x})={\displaystyle \sum _{\mathit{n}}{\zeta }_{\mathit{n}}exp(i\mathit{n}\mathit{x}),\mathit{n}=\mathrm{colon}{(n_i)}_{i=1,2,3},}\end{array}$$(15)

where n_i = 2π/T₃, T₃ is a spatial period. For the 1D case and Ω = {r < R_Ω}, we have:

$$\begin{array}{cc}\hfill {\varphi }_n(r)& =\sqrt{2/R_{\mathrm{\Omega }}}sin(\pi nr/R_{\mathrm{\Omega }}),\hfill \\ \hfill {\lambda }_n^{-1}& =4\pi n^{-2}(1-cos(n{R}_{\mathrm{\Omega }})),\underset{R_{\mathrm{\Omega }}\to \infty }{lim}{\lambda }_n^{-1}=4\pi /n^2.\hfill \end{array}$$(16)

We substitute this expression into the linearized equation and obtain:

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

The criterion for the existence of periodic solutions is:

$$\begin{array}{cc}& {\lambda }_{\mathit{II}}\le {\lambda }_{\mathit{II}}^{\mathrm{crit}},1/{\lambda }_{\mathit{II}}=-S_3{\displaystyle {\int }_0^{R_{\mathrm{\Omega }}}\mathcal{K}(y)·(sin(ny)/n)y\mathrm{d}y,}\hfill \\ \hfill & 1/{\lambda }_{\mathit{II}}^{\mathrm{crit}}\equiv -S_3{\displaystyle {\int }_0^{R_{\mathrm{\Omega }}}\mathcal{K}(y)y^2\mathrm{d}y,}\hfill \end{array}$$(18)

the latter critical value corresponds to the T₃ → ∞ case. For ${\lambda }_{\mathit{II}}<{\lambda }_{\mathit{II}}^{\mathit{crit}}$, the solutions do not become periodic. We can represent the condition of spatial periodicity of solutions in terms of “critical temperature”, namely:

$$\begin{array}{c}\hfill S_3{\lambda }_I{\theta }_{\mathrm{crit}}^{-1}exp(-C_0{\theta }_{\mathrm{crit}}^{-1}){\displaystyle {\int }_0^{R_{\mathrm{\Omega }}}|\mathcal{K}(y)|y^2\mathrm{d}y=1,}\end{array}$$(19)

whose nontrivial solution is unique and is expressed in terms of transcendental Lambert W-function. It should be noted that with a formally negative kinetic temperature criterion is satisfied for the prevailing repulsive forces in the system, that is, over very long distances. Thus, in principle, two types of periodicity of large structures exist – for the potential proper attraction at positive temperatures and for the repulsion potential for negative temperatures.

The linear equation (Eq. (14)) above does not contain any inhomogeneity. However, it does carry information about the Dirichlet conditions and contains Green’s function associated with the original Poisson equation. As the next step, we can consider inhomogeneous linear equation taking into account quadratic repulsion α(θ, Λ)x².

For this purpose, it is necessary to consider a linearization of Eq. (13) with selection of the inhomogeneous term: U(x) = U₀ − α|x|² + ϕ(x). In this case, we obtain a linear equation:

$$\begin{array}{cc}& \varphi (\mathit{x})={\lambda }^{†}{\displaystyle {\int }_{{\mathrm{\Omega }}^{\prime }}\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)\varphi ({\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }+\widehat{\beta }(\mathit{x}),}\hfill \\ \hfill & \widehat{\beta }(\mathit{x})\equiv -{\lambda }^{†}{\displaystyle {\int }_{{\mathrm{\Omega }}^{\prime }}\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)\alpha |{\mathit{x}}^{\prime }{|}^2\mathrm{d}{\mathit{x}}^{\prime }-\alpha {|\mathit{x}|}^2,{\lambda }^{†}\equiv -\stackrel{\sim }{\lambda }exp(-U_0).}\hfill \end{array}$$(20)

Obviously, we then get an inhomogeneous Fredholm-type equation with weak polar kernel 𝒦(|x − x′|); this is in accordance with (Vasil’eva & Tikhonov 2002), who introduced a new function $g(\mathit{x})\equiv \varphi (\mathit{x})-\widehat{\beta }(\mathit{x})$. By construction it is sourcewise representable in terms of the kernel 𝒦. Therefore, according to the Hilbert–Schmidt theorem (Tricomi 1957), the function g(x) can be expanded into a Fourier series in terms of the eigenfunctions of the mentioned kernel:

$$\begin{array}{cc}\hfill g(\mathit{x})& ={\displaystyle \sum _{\mathit{n}}{g}_{\mathit{n}}exp(-\mathrm{i}\mathit{n}\mathit{x}),}\hfill \\ \hfill g_{\mathit{n}}& =⟨\widehat{\beta }(\mathit{x}),{\varphi }_{\mathit{n}}(\mathit{x})⟩\equiv {\displaystyle {\int }_{\mathrm{\Omega }}\widehat{\beta }(\mathit{x})exp(-\mathrm{i}\mathit{n}\mathit{x}).}\hfill \end{array}$$(21)

Since

$$\begin{array}{cc}\hfill g_{\mathit{n}}& =\frac{1}{{\lambda }_{\mathit{n}}}⟨\varphi (\mathit{x}),\hfill \\ \hfill {\varphi }_{\mathit{n}}(\mathit{x})⟩& =\frac{{\lambda }^{†}}{{\lambda }_{\mathit{n}}}⟨\varphi (\mathit{x})-\widehat{\beta }(\mathit{x}),{\varphi }_{\mathit{n}}(\mathit{x})⟩,\hfill \end{array}$$(22)

then if the solution ϕ(x) exists, then (for λ^† ≠ λ_n)

$$\begin{array}{c}\hfill \varphi (\mathit{x})={\displaystyle \sum _{\mathit{n}}{g}_{\mathit{n}}{\varphi }_{\mathit{n}}(\mathit{x})-\widehat{\beta }(\mathit{x})=-\widehat{\beta }(\mathit{x})+{\lambda }^{†}\sum _{\mathit{n}}\frac{⟨-\widehat{\beta }(\mathit{x}),{\varphi }_{\mathit{n}}⟩{\varphi }_{\mathit{n}}(\mathit{x})}{{\lambda }_{\mathit{n}}-{\lambda }^{†}}.}\end{array}$$(23)

The solution in the form of a series (23) was obtained under the assumption that it exists. Let’s check that expression (Eq. (23)) is indeed a solution, that is, it satisfies Eq. (20), with the right side $-\widehat{\beta }(\mathit{x})$; to do this, we substitute into this inhomogeneous equation (Eq. (23)):

$$\begin{array}{cc}\hfill {\lambda }^{†}{\widehat{\mathfrak{G}}}^{\prime }{|}_{U=U_0}\varphi & ={\lambda }^{†}{\widehat{\mathfrak{G}}}^{\prime }{|}_{U=U_0}(-\widehat{\beta }(\mathit{x})+{\lambda }^{†}{\displaystyle \sum _{\mathit{n}}⟨-\widehat{\beta }(\mathit{x}),{\varphi }_{\mathit{n}}⟩{\varphi }_{\mathit{n}}/({\lambda }_{\mathit{n}}-{\lambda }^{†}))}\hfill \\ \hfill & ={\lambda }^{†}{\displaystyle \sum _{\mathit{n}}\frac{-\widehat{\beta }(\mathit{x}),{\varphi }_{\mathit{n}}⟩{\varphi }_{\mathit{n}}}{{\lambda }_{\mathit{n}}}+{({\lambda }^{†})}^2\sum _{\mathit{n}}\frac{⟨-\widehat{\beta }(\mathit{x}),{\varphi }_{\mathit{n}}⟩{\varphi }_{\mathit{n}}}{{\lambda }_{\mathit{n}}({\lambda }_{\mathit{n}}-{\lambda }^{†})}}\hfill \\ \hfill & ={\lambda }^{†}{\displaystyle \sum _{\mathit{n}}\frac{⟨-\widehat{\beta }(\mathit{x}),{\varphi }_{\mathit{n}}⟩{\varphi }_{\mathit{n}}}{{\lambda }_{\mathit{n}}-{\lambda }^{†}}=\varphi +\widehat{\beta }(\mathit{x}).}\hfill \end{array}$$(24)

We represent ϕ(x) (solution of the inhomogeneous equation) using (Eq. (24)) as:

$$\begin{array}{cc}\hfill \varphi (\mathit{x})& =-\widehat{\beta }(\mathit{x})+{\lambda }^{†}{\displaystyle \sum _{\mathit{n}}⟨-\widehat{\beta }(\mathit{x}),{\varphi }_{\mathit{n}}⟩{\varphi }_{\mathit{n}}{\lambda }_{\mathit{n}}^{-1}}\hfill \\ \hfill & +{({\lambda }^{†})}^2{\displaystyle \sum _{\mathit{n}}⟨-\widehat{\beta }(\mathit{x}),{\varphi }_{\mathit{n}}⟩{\varphi }_{\mathit{n}}{\lambda }_{\mathit{n}}^{-1}{({\lambda }_{\mathit{n}}-{\lambda }^{†})}^{-1}=-\widehat{\beta }(\mathit{x})}\hfill \\ \hfill & +{\lambda }^{†}{\displaystyle {\int }_{\mathrm{\Omega }}\underset{\mathrm{\Gamma }(\mathit{x},{\mathit{x}}^{\prime },{\lambda }^{†})}{\underbrace{(\mathcal{K}(\mathit{x},{\mathit{x}}^{\prime })+{\lambda }^{†}{\sum }_{\mathit{n}}{\varphi }_{\mathit{n}}(\mathit{x}){\varphi }_{\mathit{n}}({\mathit{x}}^{\prime }){\lambda }_{\mathit{n}}^{-1}{({\lambda }_{\mathit{n}}-{\lambda }^{†})}^{-1})}}}\hfill \\ \hfill & (-\widehat{\beta }(\mathit{x}))\mathrm{d}{\mathit{x}}^{\prime }.\hfill \end{array}$$(25)

Thus, the solution of the inhomogeneous Fredholm equation for λ^† ≠ λ_n expressed using the resolvent Γ(x, x′,λ^†) for the integral kernel −|x − x′|⁻¹ (this representation of the solution, generally speaking, is valid also in the case of deviation from the spherical symmetry of the problem). For the case λ^† = λ_n, there is no solution to the inhomogeneous equation, since $⟨\widehat{\beta }(\mathit{x}),sin(\mathit{n}\mathit{x})⟩\ne 0$ (according to Fredholm’s 3rd theorem).

For an integral kernel that differs from the classical Newtonian potential, a representation of the solution of the Fredholm equation in terms of a resolvent with sinusoidal eigenfunctions is impossible. We go on to demonstrate the technique for constructing the resolving resolvent for semi-symmetric integral Schmidt-type kernels.

We now consider the connection between the criterion for the implementation of periodic solutions and the previously introduced Helmholtz equation. The equation can be represented as the following expansion

$$\begin{array}{c}\hfill {\displaystyle \int }\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)\varphi ({\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }={\displaystyle \sum _{j=0}}{\xi }_{2j}{\mathrm{d}}^{2j}\varphi (\mathit{x})/\mathrm{d}{\mathit{x}}^{2j},\end{array}$$(26)

where

$$\begin{array}{c}\hfill {\xi }_{2j}={\displaystyle \int }\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)|\mathit{x}-{\mathit{x}}^{\prime }{|}^{2j}\mathrm{d}{\mathit{x}}^{\prime },\end{array}$$(27)

represents the even moments of interaction energy. Then the previously presented homogeneous Fredholm equation of the second kind can be reduced (in the case j = 0, 1) to the Helmholtz equation:

$$\begin{array}{c}\hfill \mathrm{\Delta }\varphi +{\lambda }_{\mathrm{H}}\varphi =0,{\lambda }_{\mathrm{H}}=(1+{\lambda }_{\mathit{II}}{\xi }_0)/({\lambda }_{\mathit{II}}{\xi }_2).\end{array}$$(28)

The solutions of this equation are spherical waves decreasing in amplitude, as follows:

$$\begin{array}{c}\hfill \sim r^{-1}{\displaystyle \sum _{j=1}^2}{\mathfrak{a}}_{j}exp(i\sqrt{{\lambda }_{\mathrm{H}}}r),\end{array}$$(29)

with a spatial period of:

$$\begin{array}{c}\hfill T=2\pi {{\lambda }_{\mathrm{H}}}^{-1/2}=2\pi \sqrt{{\lambda }_{\mathit{II}}{\xi }_2/(1+{\lambda }_{\mathit{II}}{\xi }_0)}.\end{array}$$(30)

Moreover, the admissible size of the region with the Dirichlet boundary conditions can be unlimited; in this case, assuming the Sommerfeld conditions, we can obtain the only solution to the problem for a gravitational source (including the right-hand-side of the Helmholtz equation) with the conditions at infinity as follows:

$$\begin{array}{c}\hfill \varphi (r)\sim {\displaystyle \int }{c}^2\mathrm{\Lambda }/(r\theta )exp(i\sqrt{{\lambda }_{\mathrm{H}}}r)\mathrm{d}r.\end{array}$$(31)

However, if we accept that it is possible to set spatial periodic conditions corresponding to the existence of massive multi-particle systems, then the solution of such a problem should be considered in the form of superpositions of divergent and convergent spherical waves from adjacent sources. At the same time, an interesting problem of interpretation in the physical space of the results of interference of these waves, including quasi-two-dimensional formations with increased density, is arising due to the implementation of conditions for occurrence of wave beats.

A similar situation arises in the plane case (d = 2) when the solutions of the Helmholtz equation are representable in terms of zeroth-order Hankel functions (Sveshnikov et al. 2004):

$$\begin{array}{c}\hfill {\varphi }_2={\displaystyle \sum _{j=1}^2{\mathfrak{b}}_j{H}_0^{(j)},}\end{array}$$(32)

with a change of semi-periodicity vs. monotonic decrease.

A nonlinear Hammerstein equation (Eq. (13)) with symmetric kernel and quadratic nonhomogeneous term satisfies the existence conditions for the (λ₀, U₀(x)) solutions (O’Regan & Meehan 1998). We can then consider the possibility of constructing a locally unique solutions of (Eq. (13)) with a nearly non-trivial solution of the inhomogeneous linear equation obtained above (on a “plateau” corresponding to the cosmological self-consistent gravitational potential. To do this, we apply the Bratu perturbation method, introducing new variables δλ, δU (λ = λ₀ + δλ, U(x) = U₀(x)+δU(x) and we rewrite the equation (Eq. (13)) accordingly:

$$\begin{array}{cc}\hfill \delta U(\mathit{x})& =\delta \lambda {\displaystyle {\int }_{{\mathrm{\Omega }}^{\prime }}\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)\sum _{j\ge 0}\frac{{(\delta U)}^j}{j!}{(exp(-U))}^{(j)}{|}_{U=U_0}\mathrm{d}{\mathit{x}}^{\prime }}\hfill \\ \hfill & +{\lambda }_0{\displaystyle {\int }_{{\mathrm{\Omega }}^{\prime }}\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)\sum _{\ell \ge 1}\frac{{(\delta U)}^{\ell }}{\ell !}{(exp(-U))}^{(\ell )}{|}_{U=U_0}\mathrm{d}{\mathit{x}}^{\prime },}\hfill \\ \hfill U_0& ={\lambda }_0{\displaystyle {\int }_{{\mathrm{\Omega }}^{\prime }}\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)exp(-U_0)+\alpha {\mathit{x}}^2.}\hfill \end{array}$$(33)

Representing the solution as a series in powers of a small parameter δU(x) = ∑_j ≥ 1h_j(δλ)^j and substituting into the equation, we obtain a (infinite) system of equations for determining the coefficients:

$$\begin{array}{cc}\hfill h_1(\mathit{x})& ={\lambda }_0{\displaystyle \int }\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)(-exp(U_0({\mathit{x}}^{\prime })))h_1({\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }\hfill \\ \hfill & +{\displaystyle \int }\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)exp(U_0({\mathit{x}}^{\prime }))\mathrm{d}{\mathit{x}}^{\prime },\dots ,\hfill \end{array}$$(34)

$$\begin{array}{cc}\hfill h_j(\mathit{x})& ={\lambda }_0{\displaystyle \int }\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)(-exp(U_0({\mathit{x}}^{\prime })))h_j({\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }+{\mu }_j(\mathit{x}),\hfill \\ \hfill & |{\mu }_j(\mathit{x})|<P_1·{\max }_{{\mu }_j}(|{\lambda }_0|({\mu }_{j,2}/2!+{\mu }_{j,3}/3!)+\dots +{\mu }_{j,j}/j!)\hfill \\ \hfill & +({\mu }_{j-1,1}+{\mu }_{j-1,2}/2!+\dots +{\mu }_{j-1,j-1}/(j-1)!)),\hfill \\ \hfill & |h_j(\mathit{x})|<|{\mu }_j(\mathit{x})|(1+|{\lambda }_0|P_2)<P_3,\hfill \end{array}$$(35)

$$\begin{array}{cc}\hfill & |h_j(\mathit{x})|<|{\mu }_j(\mathit{x})|(1+|{\lambda }_0|P_2)<P_3,\hfill \end{array}$$(35)

$$\begin{array}{cc}\hfill {\mu }_{j,s}& =\sum {\mu }_{i_1}\dots {\mu }_{i_s},i_1+\dots +i_s=j,i_{1,2,\dots }\in \{1,2,\dots ,j\}.\hfill \end{array}$$(37)

The above representation of δU(x) as a regular series in the neighborhood of U₀(x) due to the presence of majorants of the coefficients only. However, if the value of the Fredholm determinant for the linearized equation D(λ)→0 (at the point (λ₀, U₀)), then the second branch of the Hammerstein equation solution becomes relevant.

Thus, we get an analytical representation of the solution the Hammerstein equation for the potential, which is unique by construction and whose existence can be guaranteed in some neighborhood of equilibrium solution, U₀. This means the resistance of the system to external influences (the presence of external gravitational fields, varying within a certain spatial region) and final changes in the most self-consistent system field. The solutions of the linearized and nonlinear equations for the potential, define of a quasi-periodic nature of the structures, when the difference in the amplitudes of the local maxima of the potential depends on the relation α/λ_II. If we accept that extrema density distribution of matter correspond to nodes and anti-nodes on solutions of equation (Eq. (13)), it becomes obvious that one-dimensional (1D) and two-dimensional (2D) structures, observable in cosmological scales can be explained as completely deterministic objects, without the involvement of stochastic dynamics.

3. Linearized equations of general form and construction of a new branch of a nonlinear equation

Above, we consider the case of Fredholm equations linearized near the point of constant potential. In doing so, we assumed that the kinetic temperature parameter has a fixed value. Clearly, there is a question of clarification of the equilibrium function and of the subsequent behavior of the system of gravitating particles.

We assume that the Hammerstein equation (Eq. (13)) has a local equilibrium solution of U = U₀(x) and we consider a small deviation from it, φ(x) = U − U₀. For this, the difference of two equations of the form (13) is usually considered as:

$$\begin{array}{c}\hfill (U_0+\phi )-U_0=\lambda (\widehat{\mathfrak{G}}(U_0+\phi )-\widehat{\mathfrak{G}}(U_0)),\end{array}$$(38)

where we obtain the (homogeneous) Fredholm equation of the second kind, which takes the following form:

$$\begin{array}{c}\hfill \phi (\mathit{x})+\lambda {\displaystyle {\int }_{\mathrm{\Omega }}\underset{\mathcal{L}(\mathit{x},{\mathit{x}}^{\prime })}{\underbrace{\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)exp(-U_0({\mathit{x}}^{\prime }))}}\phi ({\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }=0.}\end{array}$$(39)

We consider the last equation, including (as before) the inhomogeneous term αx². This corresponds to the replacement of the nonlinear integral operator with its Frechet derivative (near the solution U₀). Generally speaking, to interpret the statement of the problem of inclusion in the consideration in a linear approximation, we can either assume: (1) the equilibrium state U₀(x) corresponds to homogeneous nonlinear equation and the presence of an inhomogeneous right-hand-side a priori is due to solutions-deviations from U₀; (2) or it is possible in a nonlinear Hammerstein equation to change the variable U^† = U − αx² and at the same time to modify the integral kernel appropriately

$$\begin{array}{c}\hfill \mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)\to {\mathcal{K}}^{†}(|\mathit{x}-{\mathit{x}}^{\prime }|)\equiv \mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)exp(-\alpha {\mathit{x}}^2).\end{array}$$(40)

Then Eq. (33) can be tranformed to a linear equation for U^†(x), whose kernel will contain information about the inhomogeneity of the original Hammerstein equation.

The discussion that follows mainly refers to the first variant of accounting for inhomogeneity, but it can be extended in a completely trivial way to the second option (if we are further taking U₀ → U₀ − αx²).

The kernel of the integral ℒ(x, x′) belongs to the class of Schmidt kernels and can be written in the following form:

$$\begin{array}{cc}\hfill \mathcal{L}(\mathit{x},{\mathit{x}}^{\prime })& =\mathrm{{\rm Y}}({\mathit{x}}^{\prime },\mathit{x})\mathcal{K}(|\mathit{x}-{\mathit{x}}^{\prime }|)\sqrt{exp(-U_0(\mathit{x}))·exp(-U_0({\mathit{x}}^{\prime }))}\hfill \\ \hfill & =\mathrm{{\rm Y}}({\mathit{x}}^{\prime },\mathit{x}){\mathcal{K}}^{†}(\mathit{x},{\mathit{x}}^{\prime }),\mathrm{{\rm Y}}({\mathit{x}}^{\prime },\mathit{x})\equiv \sqrt{\frac{exp(-U_0({\mathit{x}}^{\prime }))}{exp(-U_0(\mathit{x}))}}.\hfill \end{array}$$(41)

We use Γ(x, x′,λ) to denote the resolving kernel for the (weakly polar) integral kernel 𝒦(x, x′), whose expression for it was obtained in Sect. 2. Thus the resolvent satisfies the functional equation:

$$\begin{array}{c}\hfill \mathrm{\Gamma }(\mathit{x},{\mathit{x}}^{\prime },\lambda )=\mathcal{K}(\mathit{x},{\mathit{x}}^{\prime })+\lambda {\displaystyle {\int }_{\mathrm{\Omega }}(\mathit{x},{\mathit{x}}^{″})\mathrm{\Gamma }({\mathit{x}}^{″},{\mathit{x}}^{\prime },\lambda )\mathrm{d}{\mathit{x}}^{″}.}\end{array}$$(42)

If we set:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}^{†}(\mathit{x},{\mathit{x}}^{\prime },\lambda )={(\mathrm{{\rm Y}}({\mathit{x}}^{\prime },\mathit{x}))}^{-1}\mathrm{\Gamma }(\mathit{x},{\mathit{x}}^{\prime },\lambda ),\end{array}$$(43)

then we have:

$$\begin{array}{cc}& {\mathrm{\Gamma }}^{†}(\mathit{x},{\mathit{x}}^{\prime },\lambda )=\mathcal{L}(\mathit{x},{\mathit{x}}^{\prime })+\lambda {\displaystyle {\int }_{\mathrm{\Omega }}\mathcal{L}(\mathit{x},{\mathit{x}}^{″}){\mathrm{\Gamma }}^{†}({\mathit{x}}^{″},{\mathit{x}}^{\prime },\lambda )\mathrm{d}{\mathit{x}}^{″},}\hfill \end{array}$$(44)

$$\begin{array}{cc}& {\mathcal{K}}^{†}(\mathit{x},{\mathit{x}}^{″}){\mathrm{\Gamma }}^{†}({\mathit{x}}^{″},{\mathit{x}}^{\prime },\lambda )=\mathrm{{\rm Y}}({\mathit{x}}^{\prime },\mathit{x})\mathcal{K}(\mathit{x},{\mathit{x}}^{″})\mathrm{\Gamma }({\mathit{x}}^{″},{\mathit{x}}^{\prime },\lambda ).\hfill \end{array}$$(45)

The resolvent kernel corresponding to the integral kernel 𝒦(x, x′)exp(−U₀(x′)) will be the product Υ(x′,x)Γ^†(x, x′,λ), where Γ^†(x, x′,λ) is in resolution for the integral kernel 𝒦^†(x, x′).

By construction, for the ℒ kernel, we use sinusoidal oscillations with space-dependent variable coefficients as eigenfunctions, which means that when constructing a solution to a compositional solution, we can observe a superposition of a function close to a two-particle potential with repulsion of Eq. (1) and oscillations with a weakly varying amplitude.

It should be noted that based on the structure of solution (14), as well as the corresponding solution of the linearized equation for 𝒦 → ℒ, we can see that the presence of distinguished directions in expansions of the anisotropy of the vectors,

$$\begin{array}{c}\hfill {(0,0,j)}_{j=1,\dots ,n},{(0,j_1,j_2)}_{j_1=1,\dots ,J_1;j_2=1,\dots ,J_2}),\end{array}$$(46)

in the background of a homogeneous cosmological expansion and leads to the formation 1D and 2D large structures. We note that this mechanism drastically differs from those of Zeldovich pancakes, since it is completely deterministic and has nothing to do with stochastic perturbations.

The form of solving the linearized equation obtained above can be used for construction of analytic branches in the neighborhood of the above solution, and of a possible new branch of the Hammerstein equation (D(λ)→0). If λ₀ is not a kernel eigenvalue ℒ(x, x′), then in the neighborhood of λ₀ of the original Hammerstein equation has a unique holomorphic (with respect to λ − λ₀) solution U₀(x, λ), tending toward U₀(x) as λ → λ₀.

The original equation can be written as:

$$\begin{array}{cc}\hfill U_0+\varphi (\mathit{x})& =({\lambda }_0+\delta \lambda ){\displaystyle \int }{\mathcal{L}}_{+}(\mathit{x},\mathit{y})exp(U_0)(1-\varphi ({\mathit{x}}^{\prime })\hfill \\ \hfill & +(1/2!)\varphi ({\mathit{x}}^{\prime })-\dots )\mathrm{d}{\mathit{x}}^{\prime },{\mathcal{L}}_{+}=\mathcal{L}exp(-\alpha {\mathit{x}}^2).\hfill \end{array}$$(47)

We represent the function ϕ(x) as a series:

$$\begin{array}{c}\hfill \varphi (\mathit{x})=(\delta \lambda ){\varphi }_1(\mathit{x})+{(\delta \lambda )}^2{\varphi }_2(\mathit{x})+\dots ,\end{array}$$(48)

where

$$\begin{array}{c}\hfill {\varphi }_m(\mathit{x})={\psi }_m({\varphi }_1,\dots ,{\varphi }_{m-1};\mathit{x})+{\lambda }_0{\displaystyle \int }\mathrm{\Gamma }(\mathit{x},{\mathit{x}}^{\prime },{\lambda }_0){\psi }_m({\varphi }_1,\dots ,{\varphi }_{m-1};{\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }.\end{array}$$(49)

The branching condition for the equilibrium solution U₀(x) is the presence of λ₀ an eigenvalue of the kernel ℒ₊, as well as the following condition:

$$\begin{array}{c}\hfill {\displaystyle {\int }_{\omega }exp(U_0({\mathit{x}}^{\prime })){\varphi }_1({\mathit{x}}^{\prime })\mathrm{d}{\mathit{x}}^{\prime }\ne 0.}\end{array}$$(50)

In this case, we need to consider the Puiseux series in powers of (δλ)^1/2

$$\begin{array}{c}\hfill \varphi (\mathit{x})={(\delta \lambda )}^{1/2}{\varphi }_1(\mathit{x})+{(\delta \lambda )}^{2/2}{\varphi }_2(\mathit{x})+\dots \end{array}$$(51)

Then, the considered Hammerstein equation has: (1) two different solutions corresponding to one λ ≥ λ₀ and none solution corresponding to λ < λ₀, provided that ⟨exp(−U₀),ϕ₁⟩ and $\langle {\lambda }_0\exp (-U_0),{\varphi }_1^3\rangle$ have different signs; (2) two solutions for λ < λ₀ and none for λ > ;λ₀, if ⟨exp(−U₀),ϕ₁⟩ and $\langle {\lambda }_0\exp (-U_0),{\varphi }_1^3\rangle$ have the same signs (in the system we are considering, conditions must arise for the implementation of the first option).

Thus, we carried on with the following steps. Although the linear equation (Eq. (14)) does not contain an inhomogeneity, it carries information about the Dirichlet conditions and contains Green’s function (Nowakowski 2001) associated with the original Poisson equation. So, we considered the inhomogeneous linear equation taking into account quadratic repulsion α(θ, Λ)x². The use of an integral equation of the second kind for solving the internal Dirichlet problem is due to the presence of inhomogeneous boundary conditions associated with this additional term growing with increasing distance, as follows from the theorem (Gurzadyan 1985) on the general form of the gravitational potential of a sphere and a point. Therefore, we use the representation of the solution of the internal problem in the form of a double layer potential (Nowakowski 2001), as considered in detail in (Gurzadyan et al. 2023), as per Eqs. (26)–(29). The integral representation of the solution of the Dirichlet problem leads us to an equation Eq. (11), then via Hammerstein equation, to its linearized form of the Fredholm equation.

For the homogeneous Fredholm equation of the form given in Eq. (14) there is a countable system of fundamental functions in the general form of spherical harmonics (Atkinson 1997). These functions correspond to the fundamental values λ(=λ_n), which are (specifically, degenerate, generally speaking, when the non-sphericity of the problem is taken into account) roots of the algebraic equation D(λ) = 0, where $D(\lambda )=1+\sum _{j=1}^{\infty }{(-1)}^j{A}_j/j!$ is the Fredholm determinant of the kernel 𝒦 (Zemyan 2012); coefficients A_j are expressed in terms of the integral of the determinants of a symmetric matrix with elements 𝒦(x_k, x_ℓ)|_{k, ℓ = 1, …, j}. In accordance with the first Fredholm theorem (Wazwaz 2011), the solution of an inhomogeneous second kind equation is unique (of class C ∪ L₂) for values of the parameter λ that do not coincide with the set of roots D(λ); moreover (for D(λ)≠0), in accordance with the Fredholm alternative, the corresponding homogeneous equation has only a trivial solution. Then, the construction of a solution to the inhomogeneous Fredholm equation is obtained using a (resolvent) meromorphic function D(x, x′;λ)/D(λ), where D(x, x′;λ) is the first minor of the Fredholm determinant.

4. Conclusions

We consider the emergence of matter structures in the Local Universe based on the Vlasov kinetic technique. We show a deterministic mechanism of structure formation that is due to the self-consistent gravitational field of particles on local scale, as distinct to the stochastic mechanism due to the evolution of the primordial density perturbations within pancake theory on global scale. Hence, within this approach the local (late) and global (early) Universe reveal different mechanisms of structure formation. The crucial aspect of our analysis is the consideration of the modified gravitational interaction given by Eq. (1) with a repulsive cosmological term, which is based on the theorem (Gurzadyan 1985) of the identity of the sphere and point mass gravitational fields (however, with the non-force-free field within a shell).

Our analysis includes the transformation of the Poisson equation to Liouville–Gelfand equation, then a hierarchy of Dirichlet problems for differential equations of Helmholtz type. The transition from differential equations for gravitation to integral equations of the Hammerstein type is considered through bulk density. The resulting inhomogeneous Fredholm integral equation of the second kind was analysed using the eigenfunctions of the weakly polar kernel of the homogeneous equation. The solutions were found through the resolvent predicting the formation of 1D and 2D filamentary semi-periodic structures as a result of self-consistent gravitational interaction involving the cosmological constant. Methodically, the novelty of our analysis is in the use of inhomogeneous Fredholm integral equation of the second kind for solving the internal Dirichlet problem, with inhomogeneous boundary conditions and with an additional term growing with distance following from sphere-point identity theorem.

Observationally, the characteristic scale of the considered structure formation mechanism is given by the critical radius r_crit = (3GM/Λc²)^1/3, defining the mutual role of the attracting and repulsive terms in Eq. (1). This implies that the sizes of the voids can vary depending on the local region (e.g., Gurzadyan et al. 2014; Samsonyan et al. 2020, 2021) and local galactic flows can occur. To illustrate, let us consider the case of the Virgo cluster: the mass is 1.2 × 10¹⁵ M_⊙ within the radius 2.2 Mpc, the center of the cluster is located in 16.5 ± 0.1 Mpc from us (Fouque 2001; Mei et al. 2007; Kashibadze et al. 2020), then r_crit = 11.80 Mpc. This implies, that since the r_crit exceeds the distance of the Local Group from the center of the Virgo cluster the Λ-term is able to describe their observed mutual repel, as a local H-flow, in entire accordance with the reported value of the Hubble constant. Thus, the cosmological constant in the potential determines not only the appearance of the structures but also defines their scale, for instance, the sizes of the voids depending on the mean density of the local region and hence differing for various regions.

The principal conclusion of this work is that the stochastic mechanism of structure formation can be aptly replaced by deterministic one in local regions of the Universe. Along with the Hubble tension, this can be another indication for intrinsic differences in the early and late Universe.

Acknowledgments

We are thankful to the referee for valuable comments. V.M.C. is acknowledging the Russian Science Foundation grant 20-11-20165.

References