© The Authors 2025

1. Introduction

The Hubble tension (Riess 2020; Riess et al. 2024a,b; Di Valentino et al. 2021; Dainotti et al. 2023; Scolnic et al. 2024a,b) has attracted attention to the possibly genuine differences in the descriptions of the late and early Universe, including of the structure formation and the large-scale flows (e.g. Capozziello & Lambiase 2022; Bouchè et al. 2022; Bajardi et al. 2022; Capozziello & D’Agostino 2023). The problems of the nature of dark matter and dark energy are still unsolved, and diverse approaches and models are developed to consider the entire scope of the observational challenges (see Capozziello et al. 2024; de Pedreira 2024 and references therein).

The evolution of primordial density fluctuations as the origin of the cosmic structure formation (Peebles 1993) at large scales is described by the hydrodynamical pancake theory of Zeldovich (Zeldovich 1970; Shandarin 1989). It predicts the formation of the cosmic web and the filaments. At smaller scales, scales of clusters and groups of galaxies, the role of local gravitational interaction dominates. Namely, while the hydrodynamic approach assumes little influence of the intrinsic gravity of the involved particles as compared to the influence on the geometry of inhomogeneities 2D caustics; Arnold et al. 1982; Arnold 1982) of the Friedmann flow, the kinetic field approach is assumed to be suitable for relatively fast non-equilibrium phenomena in an external quasi-ordered medium, when the wave fronts overturn and in the intersection of three-dimensional matter flows channels for kinetic processes are formed (Gurzadyan et al. 2022, 2023a,b).

The possible tension of the late (local) and early (global) Universe outlines these differences in the mechanisms of structure formation and of the flows at large and small scales. Namely, the explanation of the Hubble tension as a result of two flows, a local and a global flow, with a non-identical Hubble constant, was provided in Gurzadyan & Stepanian (2021a,b) based on the theorem proved in Gurzadyan (1985). This theorem states the general function that fulfils the equivalency of gravity of a sphere and of a point mass in the following form for the force

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

This function ensures a non-force-free gravitational field inside a shell, and the second term in the left-hand side involves the cosmological constant in weak-field General Relativity and McCrea-Milne non-relativistic cosmological model (McCrea & Milne 1934; Zeldovich 1981). Importantly, Gurzadyan & Stepanian (2018), Gurzadyan (2019), Gurzadyan & Stepanian (2019, 2020) showed that Eq. (1) enables us to describe the observational data of the dynamics of groups and clusters of galaxies. Moreover, observational data support the non-force-free field inside a shell, as predicted by Eq. (1), that is, the influence of the galactic haloes on the properties of the spiral galaxies (Kravtsov 2013).

Then, the late Hubble flow is described by the equation (Gurzadyan & Stepanian 2021a,b)

$$\begin{array}{c}\hfill H_0^2=\frac{8\pi G{\rho }_0}{3}+\frac{\mathrm{\Lambda }{c}^2}{3},\end{array}$$(2)

where H₀ is the local Hubble parameter, and ρ₀ is the local mean density of matter. This equation follows from Eq. (1) and has the structure of the standard Friedmann equation, but it has a different content within the non-relativistic McCrea-Milne model with a cosmological constant.

Below, we continue the Vlasov kinetic analysis (Gurzadyan et al. 2022, 2023a,b) of the structure formation in the small-scale Universe. The essential point in the analysis is the role of the cosmological constant, and hence, that the evolutionary paths of the filament formation are the result of the self-consistent gravitational interaction of the particles along with the repulsive cosmological term. Namely, it appears that zones with a dominating influence of the ascending and descending branches of the modified gravitational potential of Eq. (1) form a topologically different multi-connected matter structure in the Universe. We analysed the of quasi-static processes that occurred in the system defined by the Vlasov-Poisson equations near its state of relative equilibrium. In this case, as shown, it is reasonable to reduce the problem to the analysis of the integral equation of a Hammerstein type. For the latter, the corresponding boundary problem is again based on the theorem in Gurzadyan (1985) because the structure of the of solutions of this integral equation and of the possible singularities can be expressed using known formalisable representations. We investigated the properties of the spectrum of the non-linear integral operator and its relation to the spectrum of the Fredholm operator for the linear potential, which enabled us to reveal the deterministic way of the formation of coherent filaments in local regions in the Universe.

2. The kinetic approach

The starting point for obtaining hydrodynamic models with singularities, that is, caustics and kinetic models without taking into account the influence of gravity, is the introduction of the Friedmann flow. This is characterized by a field of velocities whose values for each pair of points in space are proportional to the distance between these points, v = HΔr. The Friedmann flow in these models was considered as some analogue of an equilibrium state, which is perturbed either by inhomogeneity in the initial data set or by random fluctuations in the density that lead to the formation of a geometrical small inhomogeneity (adiabatically stable) of the spatial distribution of particles of the system on the background of this flow. The smallness is caused by random deviations, which are assumed a priori to be essentially limited in norm in comparison with the basis values of the mean parameters.

The system of Vlasov–Poisson equations for describing of dynamics in a system of N cosmological objects (with masses m_{i = 1, …, N} = m ≡ 1) may be represented as (see Gurzadyan et al. 2022 for details)

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

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathit{x}}^{(3)}\mathrm{\Phi }[F(\mathit{x})]{|}_{t=t_0}=A{S}_{3}G{\displaystyle \int F}(\mathit{x},\mathit{v},t_0)\mathrm{d}\mathit{v}-\frac{c^2\mathrm{\Lambda }}{2},\end{array}$$(4)

$$\begin{array}{c}\hfill S_3\equiv \mathrm{meas}{\mathcal{S}}^2=4\pi ,{\mathcal{S}}^d=\{\mathit{x}\in {\mathbb{R}}^d,|\mathit{x}=1|\},\end{array}$$

where F(x, v, t) is the distribution function of gravitationally interacting particles, A is the normalization factor for the particle density, t₀ is a fixed moment in time, and G is the gravitational constant. The system of objects or particles is considered within the finite domain of configurational space $\mathrm{\Omega }\subseteq {\mathbb{R}}^3$ (diam Ω ≤ ∞), with a C²-smooth boundary ∂Ω.

Equation (4) is the non-linear Poisson equation, which accounts for the cosmological term of Eq. (3). The third term on the right-hand side of the kinetic Equation (3) may be represented as

$$\begin{array}{cc}\hfill \widehat{G}(F;F)& =\mathit{G}(F)\frac{\partial F}{\partial \mathit{v}},\hfill \\ \hfill \mathit{G}(F)& =-{\mathrm{\nabla }}_{\mathit{x}}\mathrm{\Phi }[F(\mathit{x})],\hfill \end{array}$$(5)

$$\begin{array}{c}\hfill \mathrm{\Phi }[F(\mathit{x})]=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 +\frac{\mathrm{\Lambda }{c}^2}{12}{|\mathit{x}|}^2+{\widehat{\mathfrak{B}}}_3(\mathit{x},\mathit{x}\prime ),\end{array}$$

where ${\mathfrak{K}}_3(\mathit{x}-\mathit{x}\prime )=-{|\mathit{x}-\mathit{x}\prime |}^{-1}$, ${\widehat{\mathfrak{B}}}_3(\mathit{x},\mathit{x}\prime )$ is an operator term that takes into account the influence of the boundary conditions. The classical Newtonian potential Φ_N(r) =  − Gm/r increases monotonically in the interval r ∈ (0, +∞) (Φ_N ∈ ( − ∞, 0)), while the generalised (with cosmological term) Newton gravity potential

$$\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

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

where r_c = (Gm/(3Λc²))^1/3 (it increases in the interval r ∈ (0;  r_c] and decreases in the interval r ∈ (r_c;  ∞)).

We consider the stationary case of dynamics, F = F(x, v). However, the further analysis mainly concerns the second equation of the system, which is the Poisson equation relative to the potential, and no explicit time dependence is observed in it. When it varies, the Hilbert–Einstein–Maxwell action (Vedenyapin et al. 2021) can therefore be a separate variation over the fields (for a fixed particle distribution) and a variation over the distribution functions (with fixed fields). The approach we considered is therefore applicable for adiabatic processes at a quasi-equilibrium (weakly varying) particle distribution functions. In this case, we can use the energy substitution for a unique variable of the distribution function Vedenyapin et al. (2011): $F(\mathit{x},\mathit{v})=f(\epsilon )\in C_{+}^1({\mathbb{R}}^1)$, where ε = mv²/2 + Φ(x). Thus, the particle density in the right side of the Poisson equation can be expressed in terms of the equilibrium solution of the Vlasov equations. This solution is identical to Maxwell–Boltzmann distributions f = f₀(ε) = ANexp(−ε/θ). However, the physical meaning of the equilibrium solution of the Vlasov equation is essentially different from that of the Boltzmann equation. This solution must meet the following requirements: (1)a maximum possible statistical independence, (2)an isotropy of the velocity distribution, and (3)stationarity of the distribution in the form $F(\mathit{x},\mathit{v})=\rho (\mathit{x}){\displaystyle \prod _{i=1,2,3}}\mathfrak{f}(v_i^2)$. The substitution expression into the Vlasov equation gives

$$\begin{array}{c}\hfill {\displaystyle \sum _i}(v_i\frac{\partial ln(\rho )}{\partial x_i}-\frac{\partial \mathrm{\Phi }}{m\partial x_i}\frac{\partial \mathfrak{f}(v_i^2)}{\mathfrak{f}(v_i^2)\partial v_i})F=0,\end{array}$$(6)

and we obtain a system of ODEs,

$$\begin{array}{c}\hfill \frac{\partial (ln\rho )/\partial x_i}{-\partial \mathrm{\Phi }/\partial x_i}=\frac{\partial ln(\mathfrak{f}(v_i^2)/\partial v_i)}{m{v}_i}=-{\theta }^{-1},\end{array}$$(7)

where θ is a constant of separation of variables, and its physical meaning is kinetic temperature in the system of interacting collisionless particles (in accordance with Vlasov’s definition (Vlasov 1961, 1978), thermodynamic/collisional equilibrium is globally absent in this system).

Equation (4) for the gravitational potential can be written as

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Phi }(\mathit{x})=ANG{S}_3^2({\displaystyle {\int }_{y\in [0,\infty ]}}exp(-y^2/(2\theta ))y^2\mathrm{d}y)·\end{array}$$

$$\begin{array}{c}\hfill exp(-\mathrm{\Phi }/\theta )-\frac{c^2\mathrm{\Lambda }}{2},A,\theta ,R_{\mathrm{\Omega }}\in {\mathbb{R}}^1,\end{array}$$(8)

where R_Ω is the radius of the region Ω and is accepted in the form of a ball in configurational 3–space (as the simplest physically realizable case).

The Poisson equation (4) takes the form of the inhomogeneous equation of Liouville–Gelfand (LG; Dupaigne 2011), in which the local (generalized) temperature changes sign depending on the value of the derivative of the potential at a given point: As mentioned above, for the two-particle problem (in particular, for a formal pair in the form of a centre coalescence of the main part of the particles and the conditional extremely distant particle) the repulsive force can dominate due to the presence of a quadratic term ∼|x|². Generalised to the indefinite thermodynamics of a system of gravitating particles, this becomes similar to that for the Onsager vortices in classical hydrodynamics (Fimin & Chechetkin 2020), and the existence of solutions of the LG equation for large system sizes ensures the existence of solutions to the Vlasov Equation (3). This can be shown using the parametric Young inequality (Pokhozhaev 2010). It was shown (Gurzadyan et al. 2022) that for the conditions c²Λ ≷ 3πλ^† (θ ≷ 0), the system of Vlasov–Poisson equations has solutions of the type of distribution functions that admit the energy substitution, and the potential of the gravitational field, which has the property of convexity (in the general case, for an arbitrary R_Ω ≤ ∞, in contrast to the case of the attraction potential, for which there is a limit R_Ω < (C₀θ²/(λ^†/S₃²)²)^1/4). We used the notation ${\lambda }^{†}\equiv ANG{S}_3^2\mathcal{J}(\theta )$, $\mathcal{J}(\theta )\equiv {\displaystyle {\int }_0^{v_{\max }}}exp(-v^2/(2\theta ))v^2\mathrm{d}v$.

As already noted in Gurzadyan et al. (2022), in the formulation of the Dirichlet problem for the Poisson Equation (4) (or (8)) with a constant right-hand side on the boundary of the Ω region according to a McCrea-Milne averaging gravitational field outside the compact subdomain Ω₀ that contains a system of particles (that is situated in the region Ω (meas Ω₀ ≪ meas Ω)), we can assume the boundary condition on the ∂Ω as given in accordance with the theorem in Gurzadyan (1985).

The solution of the Dirichlet problem may be obtained with the help of an integral representation of the equation for a gravitational (double–layer) potential (Nowakowski 2001). The equation for the potential with a Maxwell–Boltzmann particle density, corresponding to an internal Dirichlet problem in a bounded domain Ω (under boundary conditions corresponding to the Milne–McCrea model) has the following form:

$$\begin{array}{c}\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^{\prime },\end{array}$$(9)

$$\begin{array}{c}\hfill \mathcal{G}(\mathit{x},\mathit{x}\prime )\equiv 4\pi {\displaystyle \sum _{\ell =0}^{\infty }}{\displaystyle \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}},\end{array}$$

$$\begin{array}{c}\hfill x_{<}=\min (|\mathit{x}|,|\mathit{x}\prime |),x_{>}=\max (|\mathit{x}|,|\mathit{x}\prime |),\end{array}$$

$$\begin{array}{c}\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\prime =-\frac{\mathit{GNm}}{R_{\mathrm{\Omega }}}-\frac{c^2\mathrm{\Lambda }{R_{\mathrm{\Omega }}}^2}{12},\end{array}$$

$$\begin{array}{c}\hfill {\lambda }_I={\lambda }^{†}/S_3.\end{array}$$

In essence, the above is the explicit form of the equation for the potential introduced in expression (5), where the Green function $\mathcal{G}(\mathit{x},\mathit{x}\prime )$ for the inner boundary value problem in the domain Ω (in this case, due to symmetry of the latter, we have ${\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}(\propto 1/R_{\mathrm{\Omega }})$). We introduce a new variable U(x)≡(Φ(x)−C₀′−C₁)/θ + α|x|², α ≡ c²Λ/(12θ), and the above equation can be written as the uniform Hammerstein integral equation,

$$\begin{array}{c}\hfill U(\mathit{x})={\lambda }_{\theta }\widehat{\mathfrak{G}}(U),\widehat{\mathfrak{G}}(U)\equiv {\displaystyle {\int }_{\mathrm{\Omega }\prime }}\mathcal{K}(|\mathit{x}-\mathit{y}|)\mathrm{\Psi }(\mathit{y},U(\mathit{y}))\mathrm{d}\mathit{y},\end{array}$$(10)

$$\begin{array}{c}\hfill {\lambda }_{\theta }\equiv \frac{{\lambda }^{†}}{\theta S_3}exp((-C_0-C_1)/\theta ),\end{array}$$

$$\begin{array}{c}\hfill \mathcal{K}(|\mathit{x}-\mathit{y}|)={-|\mathit{x}-\mathit{y}|}^{-1},\mathrm{\Psi }(\mathit{y},U(\mathit{y}))\equiv -exp(-\alpha {\mathit{y}}^2-U(\mathit{y})).\end{array}$$

For θ > 0, we obtain λ_θ > 0, the mapping $\widehat{\mathfrak{G}}(U)$ is a compact (in L₂(Ω)) non-linear operator, since the conditions of the Nemytskij–Vainberg theorem (Krasnosel’sky 1964) are satisfied for it, (${\displaystyle {\int }_{\mathrm{\Omega }}}{\displaystyle {\int }_{\mathrm{\Omega }}}{\mathcal{K}}^2(|\mathit{x}-\mathit{y}|)\mathrm{d}\mathit{x}\mathrm{d}\mathit{y}={\mathcal{K}}_{\mathrm{\Omega }}<\infty$, $\mathrm{\Psi }(\mathit{y},U(\mathit{y}))\in C(\mathrm{\Omega }\otimes \mathbb{R})$ and |Ψ(y, U(y))|≤g(y)+C_Ψ ⋅ |U|, g ∈ L₂(Ω), g(y),C_Ψ > 0).

In the expression Ψ^† ≡ ∫₀^UΨ(y, U)dU = exp(−αy²)⋅(exp(−U)−1) it is clear that Ψ^† ≤ τ₁|U|+τ₂ (x ∈ Ω), where τ_1, 2 > 0, τ₁ < 1/λ_𝒦, λ_𝒦 is the maximum eigenvalue of the integral equation kernel 𝒦. Then, in accordance with theorem 2.8 in Zabreiko et al. (1975), the Hammerstein equation has at least one solution U₀(x). We consider the question of the uniqueness of a solution or the presence of many solutions in the next paragraph.

3. Solutions and cosmological consequences

We assumed that for λ = (λ_θ)₀ Eq. (10) has a non-trivial solution U = U₀. We considered the Fredholm determinant 𝒟(λ) (Zemyan 2012) for an integral kernel of the linearised (in the vicinity O(U₀) of the basic solution U₀) Hammerstein equation $\stackrel{\sim }{\mathcal{K}}(\mathit{x},\mathit{y})\equiv \mathcal{K}(|\mathit{x}-\mathit{y}|)·\partial \mathrm{\Psi }(\mathit{y},U_0(\mathit{y}))/\partial U_0(\mathit{y})$. After linearisation (an application of the Frechet derivative), we obtained the linear Fredholm self-adjoint compact operator with a discrete spectrum (on the real axis; Krasnosel’sky 1964). If $\mathcal{D}({({\lambda }_{\theta })}_0)\ne 0$ (i. e. (λ_θ)₀ is not the characteristic value of the kernel $\stackrel{\sim }{\mathcal{K}}(\mathit{x},\mathit{y})$), then in the vicinity O((λ_θ)₀) Eq. (10) has a unique analytic (by powers of (λ_θ − (λ_θ)₀)^j, j = 1, 2, …) solution U(x|λ_θ) (Corduneanu 1991), for which lim_{λθ → (λθ)0}U(x|λ_θ) = U₀(x) (∀x ∈ Ω).

We considered the solution of (10) in the vicinity O(U₀)×O(λ_θ)₀, and introduce the perturbed characteristic value λ_θ = (λ_θ)₀ + ξ and the perturbed solution U = U₀(x)+ζ(x). Following the method in Akhmedov (1957), we substituted these expressions in Eq. (10), written in the form

$$\begin{array}{cc}\hfill -\zeta (\mathit{x})& ={({\lambda }_{\theta })}_0{\displaystyle {\int }_{\mathrm{\Omega }}}\frac{({\omega }_1(\mathit{y})\zeta +{\omega }_2(\mathit{y}){\zeta }^2+\dots )\mathrm{d}\mathit{y}}{|\mathit{x}-\mathit{y}|}\hfill \\ \hfill & +\xi {\displaystyle {\int }_{\mathrm{\Omega }}}\frac{({\omega }_0+{\omega }_1(\mathit{y})\zeta +{\omega }_2(\mathit{y}){\zeta }^2+\dots )\mathrm{d}\mathit{y}}{|\mathit{x}-\mathit{y}|},\hfill \end{array}$$(11)

and taking into account the Taylor expansion Ψ(x, U) = ∑_{j = 0, 1, …}ω_j(x)ζ^j(x) (here ζ(x) = ∑_k = 1ξ^kζ_k(x)),

$$\begin{array}{cc}\hfill {\omega }_0(\mathit{x})& =\mathrm{\Psi }(\mathit{x},U_0)=-exp(-\alpha {\mathit{x}}^2-U_0(\mathit{x})),\hfill \\ \hfill {\omega }_1(\mathit{x})& =\frac{\partial \mathrm{\Psi }(\mathit{x},U)}{1!\partial U}{|}_{U=U_0}=exp(-\alpha {\mathit{x}}^2-U_0(\mathit{x})),\dots ,\hfill \end{array}$$(12)

we obtained a system of recurrent linear and multi-linear equations for the variables ζ_i(x)

$$\begin{array}{cc}\hfill -{\zeta }_1(\mathit{x})& ={({\lambda }_{\theta })}_0{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}{\omega }_1(\mathit{y}){\zeta }_1(\mathit{y})\mathrm{d}\mathit{y}\hfill \\ \hfill & +{({\lambda }_{\theta })}_0{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}{\omega }_0(\mathit{y})\mathrm{d}\mathit{y},\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill -{\zeta }_2(\mathit{x})& ={({\lambda }_{\theta })}_0{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}{\omega }_1(\mathit{y}){\zeta }_2(\mathit{y})\mathrm{d}\mathit{y}\hfill \\ \hfill & +{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}({({\lambda }_{\theta })}_0{\omega }_2(\mathit{y}){\zeta }_1^2(\mathit{y})+{\omega }_1(\mathit{y}){\zeta }_1(\mathit{y}))\mathrm{d}\mathit{y},\dots \hfill \end{array}$$(14)

Since, by assumption, 𝒟(λ₀) ≠ 0, then for a linear non-uniform Fredholm II type equation (13) there exists a resolvent R(x, y; λ₀), and we can write

$$\begin{array}{cc}& {\zeta }_k(\mathit{x})=-{({\lambda }_{\theta })}_0{\displaystyle {\int }_{\mathrm{\Omega }}}R(\mathit{x},\mathit{y};{\lambda }_0){\mathcal{H}}_k({\zeta }_1(\mathit{y}),{\zeta }_1(\mathit{y}),\dots ,\hfill \\ \hfill & {\zeta }_{k-1}(\mathit{y}))\mathrm{d}\mathit{y}+{\mathcal{H}}_k({\zeta }_1(\mathit{x}),{\zeta }_1(\mathit{x}),\dots ,{\zeta }_{k-1}(\mathit{x}))\hfill \end{array}$$(15)

(here the operator function ℋ_k is the sum of all integrals including ζ₁, ζ₂, … up to (k − 1)–th power: ${\zeta }_k=-{({\lambda }_{\theta })}_0{\displaystyle \int {|\mathit{x}-\mathit{y}|}^{-1}}{\zeta }_1(\mathit{y})\mathrm{d}\mathit{y}+{\mathcal{H}}_k(\mathit{x})$).

Thus, we consistently and uniquely defined all functions ζ_k. Consequently, we constructed a power series for ζ(x) = ζ(ξⁱ, ζ_i(x)) ≡ ξⁱζ_i (summation over repeating indices), and it only remains to prove that the constructed series converges (for all values from some convergence circle) for the deviation parameter ξ (this completely establishes the equivalence of the expansion ζ(ξⁱ, ζ_i(x)) and the function ζ(x)).

Based on the properties of the Lyapunov–Schmidt integral operator with a weakly polar kernel (lemma 8.1 in Krasnosel’sky 1964), and explicit forms of ω_k ∝ exp(−αx² − U)/k!, we can state |∫_Ω|x − y|⁻¹ω_k(y)dy|< Z₁(=const) (k = 0, 1, …, ∀x ∈ Ω). Next, by definition, |R(x, y; λ₀)| < Z₂(=const) (∀x, y ∈ Ω). Then, we can write the major series for series ζ(ξⁱ, ζ_i(x)). We introduced an algebraic equation ζ^† = Z₃(ξ + ξζ^† + (ξ + (λ_θ)₀)⋅((ζ^†)² + (ζ^†)³ + …)), Z₃ ≡ Z₁(1 + Z₂|(λ_θ)₀|). We substituted ζ(x) = ξ^j𝜘_j|_{j = 1, 2, …} into the last equation and then compared the coefficients at different powers of the deviation variable ξ,

$$\begin{array}{cc}\hfill {\varkappa }_1& =Z_3,{\varkappa }_2=Z_3({\varkappa }_1+|{({\lambda }_{\theta })}_0|·{\varkappa }_1^2),\hfill \\ \hfill {\varkappa }_3& =Z_3({\varkappa }_2+{\varkappa }_1^2+2|{({\lambda }_{\theta })}_0|·{\varkappa }_1{\varkappa }_2+|{({\lambda }_{\theta })}_0|{\varkappa }_1^3),\dots \hfill \end{array}$$(16)

Consequently, |ζ_k(x)| < 𝜘_k (x ∈ Ω), and the convergence region of the series ζ = Z₃(ξ + ξζ + …) is equivalent to the convergence region of the series ζ(x) = ξ^jζ_j(x)|_{j = 1, 2, …}. To define the convergence region of ζ = ξⁱ𝜘_i, we ought to investigate the implicit function ζ = Z₃(ξ + ξζ + …) as a function of $\zeta =\stackrel{\sim }{\zeta }(\xi )$

$$\begin{array}{c}\hfill \stackrel{\sim }{\zeta }(\xi )=Z_3\xi +(Z_3\xi -1){\zeta }^{†}+Z_3(\xi +|{({\lambda }_{\theta })}_0|)·({({\zeta }^{†})}^2+{({\zeta }^{†})}^3+\dots )=0.\end{array}$$(17)

By the implicit function theorem (Deimling 2010), since $\partial \stackrel{\sim }{\zeta }/\partial {\zeta }^{†}=-1$ for ξ = ζ^† = 0, then there exists a convergence circle with a positive radius for the series ζ^† = Z₃(ξ + ξζ^† + (ξ + |(λ_θ)₀|)((ζ^†)² + (ζ^†)³ + …). Consequently, there exists a function ζ(x), and the function

$$\begin{array}{cc}\hfill U(\mathit{x})& =U_0(\mathit{x})+\zeta (\mathit{x})=U_0(\mathit{x})\hfill \\ \hfill & +({\lambda }_{\theta }-{({\lambda }_{\theta })}_0){\zeta }_1(\mathit{x})+({\lambda }_{\theta }-{({\lambda }_{\theta })}_0{)}^2{\zeta }_2(\mathit{x})+\dots ,\hfill \end{array}$$(18)

which is a holomorphic solution of the Hammerstein equation (10) in the vicinity O((λ_θ)₀) (herewith lim_{(λθ)→(λθ)0}U = U₀).

We recall some of the properties of the linear compact operator $\widehat{L}\varphi =\lambda {\displaystyle {\int }_{\mathrm{\Omega }}}\mathcal{K}(|\mathit{x}-\mathit{y}|)\varphi (\mathit{y})\mathrm{d}\mathit{y}$. The kernel of this operator is symmetric, of weak polar type; consequently, the operator $\widehat{L}$ belongs to the Hilbert–Schmidt type of operators. The existence of at least one characteristic function has been known since Kellogg (Kellogg 1922). Moreover, there exists a sequence of characteristic values and corresponding eigenfunctions of the investigated kernel of the linear integral operator (theorem 146 in Titchmarsh 1957). In accordance with Kalmenov & Suragan (2011), we considered

$$\begin{array}{c}\hfill {\lambda }_{\ell ,j}=R_{\mathrm{\Omega }}^{-2}·({\phi }_j^{(\ell +1/2)}{)}^2,\ell \ge 0,j\ge 1,\end{array}$$(19)

where φ_j^(ℓ+1/2) are the roots of the transcendental equation

$$\begin{array}{cc}& (2\ell +1)J_{\ell +1/2}({\phi }_j^{(\ell +1/2)})+\frac{{\phi }_j^{(\ell +1/2)}}{2}(J_{\ell -1/2}({\phi }_j^{(\ell +1/2)})\hfill \\ \hfill & -J_{\ell +3/2}({\phi }_j^{(\ell +1/2)}))=0,\hfill \end{array}$$(20)

where $J_{\nu }{(\dots )|}_{\nu \in \mathbb{R}}$ refers to the Bessel function of fractional order. The eigenfunctions corresponding to each eigenvalue λ_ℓ, j can be represented in spherical coordinates in the form ${\varphi }_{\ell ,j,m}(r,\theta ,\chi )=J_{\ell +1/2}(\sqrt{{\lambda }_{\ell ,j}}r)Y_{\ell }^m(\theta ,\chi )$ (|m|≤ℓ), Y_ℓ^m(θ, χ) = P_ℓ^m(cos(θ))cos(mχ). It should be noted that we considered the Hammerstein equation (10) for an arbitrarily shaped region Ω, but for a spherical region, we assumed ℓ = m = 0 ($Y_0^0=\sqrt{1/(4\pi )}$). The kernel $\stackrel{\sim }{\mathcal{K}}(\mathit{x},\mathit{y})=-{|\mathit{x}-\mathit{y}|}^{-1}{\omega }_1(\mathit{y},U_0(\mathbf{y}))$ of the Fredholm equation for the potential belongs to the Schmidt class, and it can be transformed into symmetrical form,

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathcal{K}}(\mathit{x},\mathit{y})=-\sqrt{{\omega }_1(\mathit{y})/{\omega }_1(\mathit{x})}{|\mathit{x}-\mathit{y}|}^{-1}\sqrt{{\omega }_1(\mathit{y}){\omega }_1(\mathit{x})}.\end{array}$$(21)

In this case, the resolvent for kernel $\stackrel{\sim }{\mathcal{K}}(\mathit{x},\mathit{y})$ is equal to $-\sqrt{{\omega }_1(\mathit{y})/{\omega }_1(\mathit{x})}{R}_1(\mathit{x},\mathit{y};\lambda )$, where R₁ is a resolvent for the symmetric kernel ${-|\mathit{x}-\mathit{y}|}^{-1}\sqrt{{\omega }_1(\mathit{y}){\omega }_1(\mathit{x})}$. In fact, ω₁(y) is a weight for the Newtonian kernel. We denote ${\lambda }_1\equiv {({\lambda }_{\theta })}_0={\lambda }_{0,1}(\stackrel{\sim }{\mathcal{K}})$ as the characteristic number to which the eigenfunction corresponds ${\varphi }_1\sim J_{1/2}(\sqrt{{\lambda }_1}r)$. To simplify the further calculations, we assumed ω₁ ≡ 1 (this can always be accomplished by multiplying of the left- and right-hand sides of Eq. (11) by $\sqrt{{\omega }_1}$ and redefining $\sqrt{{\omega }_1}\zeta \to \zeta \prime$, ${\omega }_0/\sqrt{{\omega }_1}\to {\omega }_0\prime$ etc.). We searched for solutions of Eq. (11) in the form of a series ζ = ξⁱζ_i|_{i = 1, 2, …}. After substituting the last series, we obtained a system of recurrent integral equations,

$$\begin{array}{cc}\hfill -{\zeta }_1(\mathit{x})& ={\lambda }_1{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}{\zeta }_1(\mathit{y})\mathrm{d}\mathit{y}+{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}{\omega }_0\mathrm{d}\mathit{y},\hfill \\ \hfill -{\zeta }_2(\mathit{x})& ={\lambda }_1{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}{\zeta }_2(\mathit{y})\mathrm{d}\mathit{y}\hfill \\ \hfill & +{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}({\lambda }_1{\omega }_2{\zeta }_2^2(\mathit{y})+{\zeta }_1(\mathit{y}))\mathrm{d}\mathit{y},...\hfill \end{array}$$(22)

According to the Fredholm alternative, the condition of the existence of the solution for the first equation of the system is the orthogonality of the characteristic function ϕ₁ and of the second term in the right-hand side of Eq. (19): $\widehat{\mathfrak{H}}({\varphi }_1)\equiv {\displaystyle {\int }_{\mathrm{\Omega }}}{\omega }_0(\mathit{y}){\varphi }_1(\mathit{y})\mathrm{d}\mathit{y}=0$. This case is physically unrealizable (for ϕ₁, ϕ₂, …), which can be checked directly. This means the absence (in the neighborhood of the characteristic value λ_k, k ≥ 1) of the analytic solution of the non-linear equation for the gravitational potential. Therefore, we turn to case $\widehat{\mathfrak{H}}({\varphi }_1)\ne 0$. Then, Eq. (14) is unsolvable (the condition of the Fredholm theorem is absent), and consequently, the analytic series ζ = ξⁱζ_i (for integer indices and powers) does not exist. However, we can consider the representation ζ(x) as a Puiseux series: ζ = ξ^k/2ζ_k|_{k = 1, 2, …}. We denote ξ^1/2 ≡ ν, and then, equation (11) takes the form

$$\begin{array}{cc}\hfill -\zeta (\mathit{x})& ={\nu }^2{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}({\omega }_0+\zeta (\mathit{y})+{\omega }_2{\zeta }^2(\mathit{y})+\dots )\mathrm{d}\mathit{y}\hfill \\ \hfill & +{\lambda }_1{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}(\zeta (\mathit{y})+{\omega }_2{\zeta }^2(\mathit{y})+\dots )\mathrm{d}\mathit{y},\hfill \end{array}$$(23)

and as mentioned above, the Puiseux series takes the form ζ = ν^kζ_k|_{k = 1, 2, …}. Substituting this series into equation (23), we obtain an infinite system of equations,

$$\begin{array}{cc}\hfill -{\zeta }_1(\mathit{x})& ={\lambda }_1{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}{\zeta }_1(\mathit{y})\mathrm{d}\mathit{y},\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill -{\zeta }_2(\mathit{x})& ={\lambda }_1{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}{\zeta }_2(\mathit{y})\mathrm{d}\mathit{y}\hfill \\ \hfill & +{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}({\lambda }_1{\omega }_2{\zeta }_1^2(\mathit{y})+{\zeta }_1(\mathit{y}))\mathrm{d}\mathit{y},\dots ,\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill -{\zeta }_n(\mathit{x})& ={\lambda }_1{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}{\zeta }_n(\mathit{y})\mathrm{d}\mathit{y}\hfill \\ \hfill & +{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}(2{\lambda }_1{\omega }_2{\zeta }_1(\mathit{y}){\zeta }_{n-1}(\mathit{y})+\mathfrak{Q}({\zeta }_1,\dots ,{\zeta }_{n-1}))\mathrm{d}\mathit{y}.\hfill \end{array}$$(26)

Then, we obtain ${\zeta }_1(\mathit{x})={\mathfrak{E}}_1·{\varphi }_1(\mathit{x})$, ${\mathfrak{E}}_1=\mathrm{const}.$ The condition (by the Fredholm alternative) for the existence of a solution of equation (25) can be written as

$$\begin{array}{c}\hfill {\displaystyle {\int }_{\mathrm{\Omega }}}{\displaystyle {\int }_{\mathrm{\Omega }}}{|\mathit{x}-\mathit{y}|}^{-1}({\lambda }_1{\omega }_2{\zeta }_1^2(\mathit{y})+{\omega }_0){\varphi }_1(\mathit{x})d\mathit{x}d\mathit{y}=0,\end{array}$$(27)

or, after integration by the variable x: ∫(λ₁ω₂ζ₁²(y)+ω₀)ϕ₁(y)dy = 0. When we substitute the obtained above expression ${\zeta }_1={\mathfrak{E}}_1·{\varphi }_1$ in this formula, then we have for the constant ${\mathfrak{E}}_1$ an explicit expression,

$$\begin{array}{c}\hfill {\mathfrak{E}}_1=\pm \sqrt{-\frac{{\displaystyle {\int }_{\mathrm{\Omega }}}{\omega }_0{\varphi }_1(\mathit{y})\mathrm{d}\mathit{y}}{{\displaystyle {\int }_{\mathrm{\Omega }}}{\lambda }_1{\omega }_2{\zeta }_1^3(\mathit{y})\mathrm{d}\mathit{y}}}.\end{array}$$(28)

Equation (25) may be written as

$$\begin{array}{cc}\hfill -{\zeta }_2(\mathit{x})& ={\mathcal{P}}_2(\mathit{x})+{\mathfrak{E}}_1{\varphi }_1(\mathit{x}),{\mathcal{P}}_2(\mathit{x})\hfill \\ \hfill & \equiv {\displaystyle \sum _{j=2,\dots }}\frac{{\varphi }_j(\mathit{x})}{{\lambda }_j-{\lambda }_1}{\displaystyle {\int }_{\mathrm{\Omega }}}({\lambda }_1{\omega }_2{\zeta }_1^2(\mathit{y})+{\omega }_0){\varphi }_j(\mathit{y})\mathrm{d}\mathit{y}.\hfill \end{array}$$(29)

In the general case n ≥ 3, we obtain

$$\begin{array}{cc}\hfill -{\zeta }_3(\mathit{x})& ={\mathcal{P}}_2(\mathit{x})+{\mathfrak{E}}_3{\varphi }_1(\mathit{x}),\hfill \\ \hfill {\mathcal{P}}_2(\mathit{x})& \equiv {\displaystyle \sum _{j=2,\dots }}\frac{{\varphi }_j}{{\lambda }_j-{\lambda }_1}{\displaystyle {\int }_{\mathrm{\Omega }}}(2{\lambda }_1{\omega }_2{\zeta }_1(\mathit{y}){\zeta }_2(\mathit{y})\hfill \\ \hfill & +\mathfrak{Q}({\zeta }_1,{\zeta }_2)){\varphi }_j(\mathit{y})\mathrm{d}\mathit{y},\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill -{\zeta }_n(\mathit{x})& ={\mathcal{P}}_n(\mathit{x})+{\mathfrak{E}}_n{\varphi }_1(\mathit{x}),\hfill \\ \hfill {\mathcal{P}}_n(\mathit{x})& \equiv {\displaystyle \sum _{j=2,\dots }}\frac{{\varphi }_j}{{\lambda }_j-{\lambda }_1}{\displaystyle {\int }_{\mathrm{\Omega }}}(2{\lambda }_1{\omega }_2{\zeta }_1(\mathit{y}){\zeta }_{n-1}(\mathit{y})\hfill \\ \hfill & +\mathfrak{Q}({\zeta }_1,{\zeta }_2,\dots ,{\zeta }_{n-1})){\varphi }_j(\mathit{y})\mathrm{d}\mathit{y},\hfill \\ \hfill \frac{{\mathfrak{E}}_n}{{\mathfrak{E}}_1}{\displaystyle {\int }_{\mathrm{\Omega }}}{\omega }_0{\varphi }_1(\mathit{y})\mathrm{d}\mathit{y}& ={\displaystyle {\int }_{\mathrm{\Omega }}}(2{\lambda }_1{\omega }_2{\zeta }_1(\mathit{y}){\mathcal{P}}_n(\mathit{y})\hfill \\ \hfill & +\mathfrak{Q}({\zeta }_1,{\zeta }_2,\dots ,{\zeta }_{n+1})){\varphi }_1(\mathit{y})\mathrm{d}\mathit{y}.\hfill \end{array}$$(31)

Here, ${\mathcal{P}}_n(\mathit{x})$ can be associated with the Ruston pseudo-resolvent (Ruston 1948). Consequently, we can write the series ζ = ξ^k/2ζ_k|_{k = 1, 2, ...}, and this series formally satisfies Eq. (11). We ought to demonstrate the convergence of this series in the neighbourhood O(ξ). We introduce (1)the constant ${\mathfrak{X}}_0$, defined by the condition $|{\mathfrak{E}}_1{\varphi }_1(\mathit{x})|=|{\zeta }_1(\mathit{x})|<{\mathfrak{X}}_0$ (∀x ∈ Ω), and (2)the function $S(z)=(|{\lambda }_1|+{\nu }^2)\mathfrak{M}{z}^2/(1-z/{\rho }^{‡})+{\nu }^2({\omega }_0^{(m)}+z)$, where |ω₀|< ω₀^(m), ρ^‡ ∈ (0,  ρ_max^‡), ρ_max^‡ is the convergence radius of the series ω₂(x)+ω₃(x)z + ω₄(x)z² + …, $|{\omega }_2(\mathit{x})|<\mathfrak{M}$. We substitute in the definition S(z) the decomposition $z=\nu {\mathfrak{X}}_0+{\nu }^2({\mathfrak{X}}_1+{\mathfrak{Y}}_1)+{\nu }^3({\mathfrak{X}}_1+{\mathfrak{Y}}_1)+\dots$ and formally obtain the series S(z) = ν²S₂ + ν³S₃ + …. The value S_n(z) is a major function for expression $2{\lambda }_1{\omega }_2{\zeta }_1(\mathit{x}){\zeta }_{n-1}(\mathit{x})+\mathfrak{Q}({\zeta }_1,{\zeta }_2,\dots ,{\zeta }_n)$, and the condition ${\mathfrak{X}}_n+{\mathfrak{Y}}_n$ is a majorant function for ζ_n + 1(x), n = 1, 2, …

We denote (1)$(1/{\mathfrak{E}}_1){\displaystyle {\int }_{\mathrm{\Omega }}}{\omega }_0{\varphi }_1(\mathit{x})\mathrm{d}\mathit{x}=\mathfrak{N}(=\mathrm{const})$, (2)${\mathfrak{Y}}_n=S_{n+1}{\rho }_{\mathit{max}}^{‡}$, and (3)$(\mathfrak{N}+2|{\lambda }_1|\mathfrak{M}{\mathfrak{X}}_0{a}^2){\mathfrak{X}}_n=S_{n+1}{a}^2$ (a > |ϕ₁(x)|). Consequently, $|{\mathcal{P}}_n(\mathit{x})|<{\mathfrak{Y}}_n$, ${\mathfrak{E}}_{n+1}<{\mathfrak{X}}_n$. We introduce the functions $\mathfrak{X}\equiv {\nu }^2{\mathfrak{X}}_1+{\nu }^3{\mathfrak{X}}_2+\dots$, $\mathfrak{Y}\equiv {\nu }^2{\mathfrak{Y}}_1+{\nu }^3{\mathfrak{Y}}_2+\dots$, and then, the above given definition of variables ${\mathfrak{X}}_1,\dots ,{\mathfrak{X}}_n,\dots$ and ${\mathfrak{Y}}_1,\dots ,{\mathfrak{Y}}_n,\dots$ is equivalent to solving the system of equations

$$\begin{array}{cc}\hfill \mathfrak{Y}& ={\rho }_{\max }^{‡}((|{\lambda }_1|+{\nu }^2)\frac{\mathfrak{M}{(\nu {\mathfrak{X}}_0+\mathfrak{X}+\mathfrak{Y})}^2}{1-(\nu {\mathfrak{X}}_0+\mathfrak{X}+\mathfrak{Y})/{\rho }^{‡}}\hfill \\ \hfill & +{\nu }^2({\omega }_0^{(m)}+\nu {\mathfrak{X}}_0+\mathfrak{Y}+\mathfrak{X})),\hfill \end{array}$$(32)

$$\begin{array}{cc}& (|\mathfrak{N}|+2|{\lambda }_1|\mathfrak{M}{\mathfrak{X}}_0{a}^2)\nu \mathfrak{X}=a^2(({\nu }^2+|{\lambda }_1|)\frac{\mathfrak{M}{(\nu {\mathfrak{X}}_0+\mathfrak{X}+\mathfrak{Y})}^2}{1-(\nu {\mathfrak{X}}_0+\mathfrak{X}+\mathfrak{Y})/{\rho }^{‡}}\hfill \\ \hfill & +{\nu }^2({\rho }_{\mathit{max}}^{‡}+\nu {\mathfrak{X}}_0+\mathfrak{X}+\mathfrak{Y})-{\nu }^2({\rho }_{\mathit{max}}^{‡}+\mathfrak{M}|{\lambda }_1|{\mathfrak{X}}_0^2)).\hfill \end{array}$$(33)

We replace the variables $\mathfrak{X}=\nu {\mathfrak{X}}^{†}$, $\mathfrak{Y}=\nu {\mathfrak{Y}}^{†}$. Then, the system (32)–(33) takes the form

$$\begin{array}{cc}\hfill {\mathrm{\Theta }}_1& ={\mathfrak{Y}}^{†}-\nu {\mathfrak{Y}}^{†}\frac{{\mathfrak{X}}_0+{\mathfrak{Y}}^{†}+{\mathfrak{X}}^{†}}{{\rho }^{‡}}-{\rho }_{\mathit{max}}^{‡}(({\nu }^2+|{\lambda }_1|)\mathfrak{M}\nu ({\mathfrak{X}}_0+{\mathfrak{Y}}^{†}+{\mathfrak{X}}^{†})\hfill \\ \hfill & +\nu ({\rho }_{\mathit{max}}^{‡}+\nu ({\mathfrak{X}}_0+{\mathfrak{Y}}^{†}+{\mathfrak{X}}^{†}))(1-\nu \frac{{\mathfrak{X}}_0+{\mathfrak{Y}}^{†}+{\mathfrak{X}}^{†}}{{\rho }^{‡}}))=0,\hfill \end{array}$$(34)

$$\begin{array}{cc}\hfill {\mathrm{\Theta }}_2& =(|\mathfrak{N}|+2|{\lambda }_1|{\mathfrak{X}}_0\mathfrak{M}{a}^2){\mathfrak{X}}^{†}-(|\mathfrak{N}|+2|{\lambda }_1|{\mathfrak{X}}_0\mathfrak{M}{a}^2){\mathfrak{X}}^{†}\nu \frac{{\mathfrak{X}}_0+{\mathfrak{Y}}^{†}+{\mathfrak{X}}^{†}}{{\rho }^{‡}}\hfill \\ \hfill & -(|{\lambda }_1|+{\nu }^2)\mathfrak{M}{({\mathfrak{X}}_0+{\mathfrak{Y}}^{†}+{\mathfrak{X}}^{†})}^2{a}^2-(\nu ({\mathfrak{X}}_0+{\mathfrak{Y}}^{†}+{\mathfrak{X}}^{†})\hfill \\ \hfill & -\mathfrak{M}|{\lambda }_1|{\mathfrak{X}}_0^2)(1-\nu \frac{{\mathfrak{X}}_0+{\mathfrak{Y}}^{†}+{\mathfrak{X}}^{†}}{{\rho }^{‡}})=0.\hfill \end{array}$$(35)

The Jacobi determinant of the latter system of equations is $\mathrm{\Delta }=D({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)/D({\mathfrak{X}}^{†},{\mathfrak{Y}}^{†})=-|\mathfrak{N}|<0$ for ${\mathfrak{X}}^{†}={\mathfrak{Y}}^{†}=\nu =0$. Consequently, the series of defined variables $\mathfrak{X}$ and $\mathfrak{Y}$ converges in the vicinity of the point ν = 0. The series $\nu {\mathfrak{X}}_0+{\nu }^2({\mathfrak{X}}_1+{\mathfrak{Y}}_1)+{\nu }^3({\mathfrak{X}}_2+{\mathfrak{Y}}_2)+\dots$ has a convergence circle (with a centre in the point ν = 0). This indicates that the Puiseux series ζ = ξ^k/2ζ_k|_{k = 1, 2, …} converges in the vicinity of the point ξ = 0.

We can conclude that the Hammerstein equation for the gravitational potential in the vicinity O(λ₁)∋λ (where λ₁ is one of the characteristic values of the linear Fredholm equation for the Newton potential) has two non-holomorphic solutions of the form

$$\begin{array}{cc}\hfill U(\mathit{x})& =U_0(\mathit{x})+{({\lambda }_{\theta }-{({\lambda }_{\theta })}_0)}^{1/2}{\zeta }_1(\mathit{x})\hfill \\ \hfill & +({\lambda }_{\theta }-{({\lambda }_{\theta })}_0){\zeta }_2(\mathit{x})+{({\lambda }_{\theta }-{({\lambda }_{\theta })}_0)}^{3/2}{\zeta }_3(\mathit{x})+\dots \hfill \end{array}$$(36)

We return to Eq. (10) and to the assumption that ((λ_θ)₀, U₀) is a characteristic value and the eigenvalue of the Hammerstein operator $\widehat{\mathfrak{G}}(U)$. The above analysis of the structure of the solutions of the Hammerstein equation used the fact that the kinetic temperature in the system was assumed to be positive (θ > 0). For negative temperatures (θ < 0), the properties of the solutions of the equation for the potential change slightly, although the method of their constructive continuation on the parameter plane remains the same. Since the presence of negative temperatures in the system is due to the effect of anti-gravity (due to the cosmological term in the integral equation), we assumed the Newtonian kernel to be continuous and a bounded function (although we can also use the Hilbert-Schmidt theory), and, on the basis of theorem 2 in Jingxian & Bendong (1997), we conclude that Eq. (10) has at least two analytic solutions outside the points of the spectrum of the Fréchet derivative of the Hammerstein operator ${({\lambda }_{\theta <0})}_k(\widehat{\mathfrak{G}}\prime )$. In the vicinity of the characteristic points corresponding to the solutions of the non-linear equation and simultaneous coincidence with the special points of the resolvent of the linearized equation, the solutions are continued algebraically and can be written out in the form of Puiseux series on half-integer degrees of the deviation from the solution points.

In physical terms, the expression (36) can be written as

$$\begin{array}{c}\hfill \mathrm{\Phi }(\mathit{x})=(C_1^{‡}+c^2\mathrm{\Lambda }/12){|\mathit{x}|}^2+U_0(\mathit{x}){\displaystyle \sum _i}(ANG\mathcal{J}(\theta )exp(C_1^{‡}/\theta )-{({\lambda }_{\theta })}_0{)}^{i/2},\end{array}$$(37)

$$\begin{array}{cc}\hfill \mathcal{J}(\theta >0)& =-exp(-v_{\max }^2/(2\theta ))v_{\max }\theta +{\theta }^{3/2}\sqrt{\pi /2}\mathrm{erf}(v_{\max }/\sqrt{2\theta }),\hfill \\ \hfill \mathcal{J}(\theta <0)& =exp(v_{\max }^2/(2\theta ))v_{\max }\theta -\theta \sqrt{-\pi \theta /2}\mathrm{erf}(v_{\max }/\sqrt{-2\theta }),\hfill \end{array}$$

where U₀(x) is a solution of the Hammerstein equation for the potential (existing due to the properties of the corresponding integral operator), λ₁(=λ_θ)₀) is the characteristic value of kernel ω₁/|x − y| (for the Fredholm equation) in the vicinity of the point λ_θ ∈ ((λ_θ)₀ − ϵ, (λ_θ)₀ + ϵ). If sgn(∫_Ωω₀ϕ₁(y)dy) ≠ sgn(∫_Ωλ₁ω₂ϕ₁³(y)dy), then for λ > λ₁ there exist two solutions, and for λ < λ₁ there are no solutions. For sgn(∫_Ωω₀ϕ₁(y)dy) = sgn(∫_Ωλ₁ω₂ϕ₁³(y)dy), then for λ < λ₁ there exist two solutions, and for λ > λ₁ there are no solutions.

Thus, we showed that the solution of the system of Vlasov-Poisson equations in the case of using the energy substitution in the stationary Vlasov equation and choosing the quasi-Maxwell distribution (depending on the kinetic temperature) as a basis solution of the kinetic equation, leads to the non-linear Hammerstein integral equation for the potential, including repulsion due to the cosmological term in Eq. (1). The Hammerstein equation in essence is an equation for the self-consistent field in the system of massive many-particles on cosmological scales. This leads to the simultaneous realisation in the cosmological system with a self-consistent field of two types of ordering: A significant increase in the matter density and its decrease (at the one-dimensional motion of matter in the channel, walls appear with a high value of the gravitational potential, and voids, which are practically empty spaces between walls). Oscillations of the eigenfunctions even in the simplest case of spherical symmetry of the Green function can form analogues of interference patterns, which leads to the production of structure formation in a homogeneous medium far from the line that connects the massive objects. Since the eigenfunctions of the linearised version of the equation for the potential are proportional to the Bessel function (with a rapidly decreasing weight function), the loss of periodicity of the solutions in physical space becomes obvious. Since for non-linear integral operators it is possible to establish the the continuum character of the spectrum, we can observe the effect of secondary solutions (on the hyperplane of parameters near the basic solution of the equation) that coexist with the initial solutions, which leads to the construction of the second type of periodicity: The translation of the solution isotropically along all directions. Thus, the composition of solutions with a potential of two types in the presence of additional translational transgression due to the structure of the spectrum results in a distribution of the field and matter in space that is peculiar to the cosmic web. The scaling of the semi-periodic structures corresponds to the scale of the involved constant, that is, the cosmological constant in agreement with observational data (Gurzadyan et al. 2023a,b). Incidentally, the change in the kinetic temperature at the zero-point transition can lead to a dipole-type structures involving a repeller.

The study performed above of the properties of solutions of the non-linear inhomogeneous integral Hammerstein equation obtained from the Dirichlet boundary value problem allowed us without the need to invoke the assumption of a correlation of fluctuations to describe in a completely deterministic manner the existence in the vicinity of equilibrium points of many-particle systems (at different scales) defined by secondary solutions with significantly different properties. We also proved in this analysis that these solutions exist non-locally. The presence of analytical and non-holomorphic branches of solutions of the equation for the potential given by Eq. (10) from a physical point of view in distinguished directions with sharp gradients of the density parameter of the medium has to be associated with the walls as the internal spaces between the voids. It should be noted that the cosmological term in the Hammerstein equation leads to the necessity of introducing of the concept of a modified kinetic temperature as defined by Vlasov (Vlasov 1961, 1978). In contrast to the linear analysis based on the properties of the Fredholm equations of the second kind, which, as was previously established by the authors (Gurzadyan et al. 2022, 2023a,b), should lead to periodic and mixed-periodic cosmological structures, the method proposed in this study is more refined and suitable for numerical modelling. It can be used in certain cases as an alternative to statistical computations.

4. Conclusions

We performed an analysis of the filament formation in the late Universe based on the Vlasov kinetic technique and using the modified law of gravitational interaction with a cosmological constant, Eq. (1). While the primordial density fluctuations were considered to lead to an early structure formation via the Zeldovich pancake theory, that is, via stochastic dynamics, the Vlasov self-consistent field mechanism can provide structure features at late stages of the evolution.

We continued the rigorous analysis of the Vlasov-Poisson equations of Gurzadyan et al. (2022, 2023a,b). The equations reduced to the Hammerstein integral equation were shown to predict semi-periodic structures, that is, those peculiar to the cosmic web. The role of the repulsive term, that is, of the cosmological constant, is a key role because it is also defines the scaling of the formed semi-periodic structures, that is, of the voids. In this context, the considered mechanism of the formation of the structures in the local Universe complements the role of the cosmological constant term in Eq. (1) in the appearance of two galactic flows. The local flow is determined by the in McCrea-Milne non-relativistic model, and the global flow is determined by the Friedmann equations, while the discrepancy of the Hubble constants of the flows thus naturally explains the Hubble tension.

Thus, the Hubble tension (Riess et al. 2024a,b; Scolnic et al. 2024a,b) and other observational indications of the possible genuine difference in the properties of the early and late Universe can serve as tests for the mechanism of the structure formation on the local scale, as discussed above. In particular, the measurement by the Dark Energy Spectroscopic Instrument (DESI) collaboration (Said et al. 2024) of the local value of the Hubble constant and the distance to the Coma cluster using the fundamental plane of galaxies confirmed this tension. An additional informative window for tracing the cosmological evolution arises from the released DESI first-year data (Adame et al. 2024a,b,c) on the baryon acoustic oscillation using several tracers. The data revealed additional tensions with the ΛCDM model. These observational data support the approach of our study, that is, the need of developing different, more refined descriptions of the properties of the late Universe. In particular, this approach, in which the scaling radius r_c = (Gm/(3Λc²))^1/3 defined the role of the cosmological constant term in Eq. (1), enabled us to fit the dynamics of the galactic flow in the vicinity of the Virgo supercluster or the Laniakea supercluster, for instance (Gurzadyan & Stepanian 2021a,b).

The revelation of the dominating role of the self-consistent gravitational field in the formation of the cosmic structures on local cosmological scales can be considered the main result of our rigorous analysis. The refined comparison of the predictions of the considered mechanism with the observational surveys on the voids and of the filaments, for example, to study the distributions of their scales versus redshift, can therefore be of particular interest.

Acknowledgments

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

References