© V. G. Gurzadyan et al. 2022

1. Introduction

The evolution of primordial density fluctuations is considered to be responsible for the observed large-scale matter distribution in the Universe (Peebles 1993; Bisnovatij-Kogan 2011). The assumptions on the spectrum of the initial fluctuations and analysis of various phases of the evolution of the perturbations enable one to cover the notable features of the observational surveys.

Zeldovich’s pancake theory (Zeldovich 1970; Arnold et al. 1982; Arnold 1982; Shandarin & Zeldovich 1989; Shandarin & Sunyaev 2009) is currently the basis of studies of the filament formation, including, including a direct comparison with observations. Recently, the Hubble tension (Riess 2020; Di Valentino et al. 2021; Riess et al. 2022), that is the discrepancy between the local and global values of the Hubble parameter, and other observational tensions, appear to indicate the need to reconsider the views regarding the dynamics and structure at least in the local Universe. The modified gravity approaches are actively studied among other options; see Capozziello et al. (2020), Di Valentino et al. (2021), and references therein. One of those approaches is based on a key principle of gravitational interaction being reconsidered, that is to say the equivalency of the gravity of the sphere and of the point mass located in its center. Namely, the theorem proved in Gurzadyan (1985) states that the general function satisfying that principle has the following form:

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

where the first term in the r.h.s. is Newton’s law, while the second term defines the cosmological constant Λ in the McCrea–Milne cosmological scheme (McCrea & Milne 1934) and it corresponds to the cosmological term in the solutions of Einstein equations and weak-field general relativity (GR; Gurzadyan & Stepanian 2018, 2019).

In Zeldovich (1981), the validity of the non-GR description of the local Universe was outlined. Equation (1) enables one to describe the dynamics of groups and clusters of galaxies (Gurzadyan & Stepanian 2018, 2019, 2020; Gurzadyan 2019), that is at distance scales when the second term becomes non-negligible with respect to the Newtonian term.

Equation (1) possesses a remarkable feature, namely, it defines a non-force-free field inside a spherical shell, thus differing from Newton’s law when the shell does not influence its interior. In fact, there are observational indications that the galactic halos do indeed determine certain features of the spirals and disks (Kravtsov 2013), thus supporting the non-force-free shell’s concept. Equation (1) enables one to describe the dynamics of groups and clusters of galaxies, and especially the Hubble tension, as a possible indication of two flows, a local one described by the McCrea–Milne nonrelativistic equation with a cosmological term and the global one, described by Friedmannian equations (Gurzadyan & Stepanian 2021a,b).

We study the role of law Eq. (1) within the Vlasov kinetic formalism (Vlasov 1950, 1967; Vedenyapin et al. 2021), thus outlining the notable differences from previous approaches to filament formation. On the one hand, it is a continuation of Zeldovich’s approach of the validity of the nonrelativistic treatment of the local structure formation. On the other hand, it also considers the role of the self-consistent gravitational interaction in that process, however now with a repulsive term of Eq. (1), thus continuing the studies in Gurzadyan et al. (2022).

The modification of the Newtonian gravitational potential in Eq. (1) (Gurzadyan 1985) leads to an effective negative pressure in the system of interacting particles. So the equilibrium condition for an element of matter, taking the repulsive effect of the Λ term into account, determines the scales of the formation of macro-structures in the Universe.

Namely, at scales smaller than corresponding the one of the condition of balance of gravitational forces and the repulsive Λ term at the Jeans wavelength determined by the inflection point of potential in Eq. (1) (Fig. 1), the collapse of homogeneous matter occurs, and at larger scales fragmentation of the system occurs. From the scale of the inflection in the potential Eq. (1), we obtain the mean dimensions of the voids, that is around 25 Mpc, if oriented on the mean density in the Virgo Supercluster.

Fig. 1. Modified gravitational potential ΦGN(r).

Another representation of the mutual balance of the gravity and repulsive term in Eq. (1) can be the abovementioned two-flow nature of the Hubble flow (Gurzadyan & Stepanian 2021a,b).

The next principal issue is the coherence of the relaxation of perturbations in the system along the selected direction, that is to say at large structures – along two mutually perpendicular directions. Zeldovich’s mechanism of coherence in Lagrangian coordinates (Zeldovich 1970; Arnold et al. 1982) takes the existence of synchronism of anisotropic damping of emerging density fluctuations into account, which is natural within the framework of the theory of long-range interactions. For a system of massive particles, their uniform distribution cannot be considered as a global equilibrium in configuration space (this is implied by the divergence of the potential in the integral form of the Poisson equation). That is why it is necessary to refer to the representation of the aforementioned energy substitution through the “minimum” potential function of the gravitational interactions Φ_min(x) for a many-particle system with an additional condition that singles out its singular point.

In a system of particles obeying the generalized law of Eq. (1), there are several types of singular points in the configuration space (in the general case of its subdomains), the main ones are the following:

• (1)

  libration points in a many-particle system;

• (2)

  inflection points of the generalized gravitational potential of Eq. (1); see Fig. 1.

Regarding the first type, we can definitively state that in the absence of a high degree of symmetry in the system, finding such points is practically not feasible (exceptions are only Lagrange points on Keplerian orbits), and their stability requires the development of special methods, see Schaeffer (2004), for example. Particularly, in Gurzadyan et al. (2022), the solutions of the nonlinear Liouville–Gelfand equation (LGE; Gelfand 1963) are studied to describe the inter-particle potential in the stationary version of the boundary value problem. This shows the need for special tools and a technique to reveal the features of the structures that emerged as a result of the Vlasov treatment of the problem.

Below, we concentrate on analyzing the dispersion relations for Vlasov–Poisson equations, corresponding to conditionally equilibrium solutions, defining pseudo-ordered structures of low dimensions (walls), semi-periodic structures and inter-structure regions (voids).

2. The initial perturbation problem for the linearized integral Vlasov–Poisson equations

The system of Vlasov–Poisson equations for describing cosmological dynamics in a system of N particles of equal masses m_i = m (stars, galaxies, ...) can be represented 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 \end{array}$$(2)

$$\begin{array}{cc}& \widehat{G}(F;F)\equiv -{\mathrm{div}}_{\mathit{v}}(m^{-1}({\mathrm{\nabla }}_{\mathit{x}}(\mathrm{\Phi }(F))F),\hfill \\ & {\mathrm{\Delta }}_{\mathit{x}}^{(d)}\mathrm{\Phi }(F){|}_{t=t_{\ast },\forall t_{\ast }\in \mathbb{T}}=4\pi Gm{\displaystyle \int }F(\mathit{x},\mathit{v},t_{\ast })\mathrm{d}\mathit{v},\hfill \end{array}$$(3)

where F(x, v, t) = NF₁(x, v, t). The third term on the right-hand side of this kinetic equation can be represented as a d–dimensional “source–like” form, where

$$\begin{array}{cc}\hfill \widehat{G}(F;F)& =m^{-1}\mathit{G}(F)\frac{\partial F}{\partial \mathit{v}},\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill \mathit{G}(F)& =-{\mathrm{\nabla }}_{\mathit{x}}{\displaystyle \int \int }{\mathfrak{Y}}_d^{(\vartheta )}(\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 {\mathfrak{Y}}_3^{(\vartheta )}(\mathit{x}-{\mathit{x}}^{\prime })& =\frac{\mathit{Gm}}{|\mathit{x}-{\mathit{x}}^{\prime }|}+\frac{\vartheta }{6}{c}^2{\mathrm{\Lambda }}^2{|\mathit{x}-{\mathit{x}}^{\prime }|}^2,\hfill \end{array}$$(5)

where ϑ ∈ {0;  1} and ϑ = 1 corresponds to taking into account of repulsive Λ > 0.

The Newtonian potential Φ_N(r) = − γm/r increases on the interval r ∈ (0, +∞) (Φ_N ∈ ( − ∞, 0)), while the potential of the second term of Eq. (1) (see Fig. 1),

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathit{GN}}(r)\equiv -{\mathfrak{Y}}_3^{(1)}(r),\end{array}$$(6)

has maximum

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

at

$$\begin{array}{c}\hfill r_c=(3Gm/\mathrm{\Lambda }{c}^2{)}^{1/3}.\end{array}$$(8)

As shown in Gurzadyan & Stepanian (2021a), for example for the parameters of the Virgo Supercluster, r_c yields around 13 Mpc.

The Poisson Eq. (3) takes the form of the Liouville–Gelfand equation (Gelfand 1963) if we take the density particle distribution determined by Maxwell–Boltzmann distribution functions F₀(x, v, t_*)∝Ω(ε[v(t_*);Φ_min(x)]|_ext) in terms of velocities.

In this case, in the energy neighborhood of the local extremum

$$\begin{array}{c}\hfill \mathrm{d}\mathrm{\Omega }/\mathrm{d}\epsilon {|}_{{\mathrm{\Phi }}_{\min }\to {\overline{\mathrm{\Phi }}}_{\min }(r_c)}\approx 0,\end{array}$$(9)

we get the relative equilibrium solution for system Eqs. (2) and (3); it is characterized, in addition to the distinguished energy level of particles on the regular branch of the potential (in this case, the force term in the Vlasov equation is either canceled or is the turning point on the phase plane for the evolution of the system), also the average temperature of the system T and potential gauge or density ρ₀ = Aexp(−Φ(0)/T) spatial distribution in the point are the following

$$\begin{array}{cc}\hfill F_0(\mathit{x},\mathit{v}){|}_{t=t_{\ast }}& =Bexp(-\frac{m{\mathit{v}}^2}{2T}-\frac{\mathrm{\Phi }}{T}+a_j{\mathcal{I}}_j){|}_{\mathrm{\Phi }={\mathrm{\Phi }}_{\min }},\hfill \\ \hfill B& \equiv N{\rho }_0(\frac{m}{2\pi T}{)}^{3/2}.\hfill \end{array}$$(10)

In addition to the energy, the arguments of the function are solutions to the first equation of the Ω system of other integrals’ motion of ℐ_j, but no changes in this case develop, and one can consider ℐ_j ≡ 0.

Passing to the density via the integration of the distribution function over velocities, we obtain the Poisson Eq. (3) in the form

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}^{(3)}\mathrm{\Phi }& =4\pi GBmexp(-\frac{\mathrm{\Phi }}{T}){\displaystyle \int \int \int }exp(-\frac{m{\mathit{v}}^2}{2T})\mathrm{d}\mathit{v}\hfill \\ & =4\pi GmN{\rho }_{0}exp(-\frac{\mathrm{\Phi }}{T}).\hfill \end{array}$$(11)

Renaming U = −Φ/T, λ = 4πγmNρ₀T⁻¹, we obtain an equation in the form Δ⁽³⁾U + λexp(U) = 0. In what follows, as the gravitational potential, we use the function Φ(x) and the reduced function U(x). According to Gidas–Ni–Nirenberg (Gidas et al. 1979) all solutions of the last equation are radially symmetric. Therefore, we need to consider the solutions for sustainable branches U_λ(r) of the equation

$$\begin{array}{c}\hfill U_{\mathit{rr}}^{″}+(2/r)U_r^{\prime }+\lambda exp(U)=0\end{array}$$(12)

and corresponding to them in Eq. (6) radial-symmetric functions

$$\begin{array}{cc}& {\mathrm{\Phi }}_{\lambda }(r)=-TU(r),\Vert U\Vert (\lambda )\mathrm{is}\mathcal{J}-\mathrm{valued}\mathrm{function},\hfill \end{array}$$(13)

$$\begin{array}{cc}& \mathcal{J}{|}_{\lambda \in (0,2)\cup (2,3.32)}\in \{1,2,\dots ,{\mathcal{J}}_{\max }\},\mathcal{J}{|}_{\lambda =2}=\infty ,\hfill \\ & \mathcal{J}{|}_{\lambda =3.32}=1\hfill \end{array}$$(14)

(for λ ≳ 3.32 there are no regular solutions of the Liouville-Gelfand equation).

We are interested in a qualitative analysis and description of the properties of solutions of a linearized Vlasov–Liouville–Gelfand system, leading to the possible realization of pseudo–ordered structures. Thus, we assume that the transfer processes proceed slowly enough, so that an unbalanced particle system allows for the use of a quasi–equilibrium microcanonical ensemble.

We consider the case of a d-dimensional dynamical system (d = 2 and 3) and the linearized version of the Vlasov equation; in this case, this means that the dynamics of the system of gravitating particles is considered with respect to the background of its conditional equilibrium state that changes over time in a known way, in other words autonomously, on the full iso–energetic surface of the system.

Thus, we arrive at a linear integro–differential equation, as in (2d + 1)–dimensional phase space ${\mathbb{R}}_{\mathit{x}}^{\mathrm{d}}\oplus {\mathbb{R}}_{\mathit{v}}^{\mathrm{d}}\oplus \mathbb{T}$) reduced to the integral equation (over $\mathbb{T}\subseteq {ℝ}_t^1$). We represent its solution in the following form: F(x, v, t) = F₀(x, v)+f(x, v, t), |F₀|≫|f|. Substituting this expression into the Vlasov equation, we obtain

$$\begin{array}{c}\hfill \frac{\partial f}{\partial t}+\mathit{v}\frac{\partial f}{\partial \mathit{x}}=\frac{\partial F_0(\mathit{x},\mathit{v})}{\partial \mathit{v}}·\frac{\partial }{\partial \mathit{x}}{\displaystyle \int \int }{\mathfrak{Y}}_3^{(\vartheta )}(|\mathit{x}-{\mathit{x}}^{\prime }|)f({\mathit{x}}^{\prime },{\mathit{v}}^{\prime },t)\mathrm{d}{\mathit{v}}^{\prime }\mathrm{d}{\mathit{x}}^{\prime }.\end{array}$$(15)

Using the inversion of the total derivative operator (as semi-groups of linear operators), we arrive at the form of the Vlasov equation in the representation of shifts by trajectories

$$\begin{array}{cc}\hfill f(\mathit{x},\mathit{v},t)=& f(\mathit{x}-\mathit{v}(t-t_0),\mathit{v},t_0)\hfill \\ \hfill & +{\displaystyle {\int }_{t_0}^t}\mathrm{d}{t}_{\prime }^{}\stackrel{\sim }{B}{\mathrm{\nabla }}_{\mathit{v}}exp(-\frac{m{\mathit{v}}^2}{2T}-\frac{1}{T}{\mathrm{\Phi }}_{-}(t-t_{\prime }^{}))·{\widehat{\mathcal{K}}}_d^{(-)}[\vartheta ]\varrho ,\hfill \end{array}$$(16)

$$\begin{array}{cc}& {\mathrm{\Phi }}_{-}(t-t^{\prime })\equiv {\mathrm{\Phi }}_{\min }(\mathit{x}-\mathit{v}(t-t^{\prime })),\hfill \end{array}$$(17)

$$\begin{array}{cc}& {\widehat{\mathcal{K}}}_d^{(-)}[\vartheta ]\varrho \equiv \frac{\partial }{\partial \mathit{x}}{\displaystyle \int }{\mathfrak{Y}}_d^{(\vartheta )}(|\mathit{x}-{\mathit{x}}^{\prime }-\mathit{v}(t-t^{\prime })|)\varrho ({\mathit{x}}^{\prime },t^{\prime })\mathrm{d}{\mathit{x}}^{\prime }.\hfill \end{array}$$(18)

Here, 𝜚(x′,t′) ≡ ∫f(x′,v′,t′) dv′ and t₀ certain initial moment, allocated on the time axis, for example, at a perturbation in a system described by retarded potentials, $\stackrel{\sim }{B}=B/m$.

To simplify the calculations, we assume the temperature to be constant in that part of the system in which we are interested; if we introduce a temperature field, then it is also necessary to take the shift of the argument of this field into account when transitioning to Eq. (16) (the presence of a nonconstant temperature field can modify the obtained results).

The solution f(x, v, t) is sought in the form of integrals of Fourier–Laplace type from spatial harmonics’ waves with decreasing amplitude (or waves with gaps in the coefficients). Then it is in the form ${\int }_{{ℝ}_{\mathit{w}}^3}{f}_{\mathit{w}}(\mathit{v},t)\exp (i\mathit{w}x)\text{d}\mathit{w}$, $\mathit{w}={(w_{\mathcal{l}})}_{\mathcal{l}=1,2,3}^T$, ${\mathit{w}}_{\ell }={\mathit{w}}_{\ell }^{(Re)}+i{\mathit{w}}_{\ell }^{(Im)}$ (for the density 𝜚(x, t), we have $\mathrm{\#x03F1;}(\mathit{x},t)={\int }_{{ℝ}_{\mathit{w}}^3}{\varrho }_{\mathit{w}}(t)\exp (i\mathit{w}x)$dw).

The kernel on the right-hand side of the definition of ${\widehat{\mathcal{K}}}_3\varrho$, together with the exponentially decreasing part of density, can become an integrable function (in a limited region of space ℝ³ due to the “infinite mass problem”), obviously, when the appropriate conditions are met for imaginary indicators’ exponents, ${\mathit{w}}_{\ell }^{(Im)}$. The elements f_w and 𝜚_w belong to a set of coefficients of the generalized integral Fourier (Titchmarsh 1962), and the existence conditions for these integrals a priori lack an exponential growth distribution function f(x,  v,  t) (as we see below, one can deduce the quite obvious, albeit somewhat cumbersome, criterion for this fact).

Before moving on to Fourier–images, we must integrate Eq. (16) in terms of velocities, to form an equation for density. On its left-hand side is the Fourier coefficient 𝜚_w(t), and the first term on the right–hand side of the resulting integral equation has the form

$$\begin{array}{c}\hfill {}^{[1]}\mathcal{D}{f}_{\mathit{w}}(t)={\displaystyle \int }{f}_{\mathit{w}}(\mathit{v},t_0)exp(-i\mathit{w}\mathit{v}(t-t_0))\mathrm{d}\mathit{v}.\end{array}$$(19)

The second term on the right-hand side of Eq. (16) can be converted to ${\int }^{t_{t_0}}{}^{[2]}{\mathcal{D}}_{\mathit{w}}(t-t^{\prime })\mathrm{d}{t}^{\prime }$, where

$$\begin{array}{c}{}^{[2]}{\mathcal{D}}_{\mathit{w}}(t-t^{\prime })={\displaystyle {\int }_{{\mathbb{R}}_{\mathit{v}}^3}}\mathrm{d}\mathit{v}[G\stackrel{\sim }{B}{\mathrm{\nabla }}_{\mathit{v}}\hfill \\ exp(-\frac{m{\mathit{v}}^2}{2T}-\frac{{\mathrm{\Phi }}_{-}(t-t^{\prime })}{T}){]}_{1}exp(i\mathit{w}(\mathit{x}-\mathit{v}(t-t^{\prime })))\hfill \\ \times \frac{\partial }{\partial \mathit{x}}[\int \mathrm{d}{\mathit{x}}^{\prime }{|\mathit{x}-\mathit{v}(t-t^{\prime })-{\mathit{x}}^{\prime }|}^{-1}\hfill \\ exp(-i\mathit{w}(\mathit{x}-\mathit{v}(t-t^{\prime })-{\mathit{x}}^{\prime }))\varrho ({\mathit{x}}^{\prime },t^{\prime }){]}_2.\hfill \end{array}$$(20)

It should be noted that, in the vicinity of the point $\mathrm{\Phi }={\overline{\mathrm{\Phi }}}_{\min }$ in the phase plane, the fixed branch of the function of the extrema of the gravitational potential is close to constant or has a locally parabolic structure (which corresponds to the motion of a particle near a state of stability at the generalized Lagrangian point).

The change in the temperature parameter can lead to a change in the multivalued potential branch and the loss of the characteristics of a singular point of the function Φ_min. As a consequence, among other things, the impossibility of a significant difference in the velocity dispersion in structures can, in principle, be explained via the Vlasov–Poisson formalism: their ordering is controlled by the local temperature of the system, the change of which leads to a nonsmooth, jump-like bifurcation of the potential branch, with the inevitable displacement of libration points and, accordingly, a possible restructuring of the existing structure.

Renaming |x − v(t − t′) − x′| = ξ and integrating [...]₂ in brackets, we get

$$\begin{array}{c}\hfill 4\pi {\displaystyle {\int }_0^{\infty }}exp(w^{\mathit{Im}}\xi )\frac{sin(w^{\mathit{Re}}\xi )}{w^{\mathit{Re}}\xi }\xi \mathrm{d}\xi \equiv \mathcal{H}(w^{\mathit{Re}},w^{\mathit{Im}}),\end{array}$$(21)

and finally

$$\begin{array}{c}\hfill {}^{[2]}{\mathcal{D}}_{\mathit{w}}(t-t^{\prime })=\end{array}$$(22)

$$\begin{array}{c}\hfill \mathcal{H}(w^{\mathit{Re}},w^{\mathit{Im}}){\displaystyle {\int }_{{\mathbb{R}}_{\mathit{v}}^3}}d\mathit{v}exp(i\mathit{w}(\mathit{x}-\mathit{v}(t-t^{\prime })))i(\mathit{w}·[...{]}_1).\end{array}$$

Thus, to obtain the Fourier transforms of the density 𝜚_w, we use the Volterra (second kind) integral equation

$$\begin{array}{c}\hfill {\varrho }_{\mathit{w}}(t)={}^{[1]}\mathcal{D}{f}_{\mathit{w}}(t)+{\displaystyle {\int }_{t_0}^t}{}^{[2]}{\mathcal{D}}_{\mathit{w}}(t-t^{\prime }){\varrho }_{\mathit{w}}(t^{\prime })\mathrm{d}{t}^{\prime }.\end{array}$$(23)

The solution for this type of equation can be obtained in terms of the one–sided Laplace transform (by the time variable)

$$\begin{array}{c}\hfill \widehat{\mathcal{L}}{f}_{\mathit{w}}(\omega )={\displaystyle {\int }_0^{\infty }}{f}_{\mathit{w}}(t)exp(-\omega t)\mathrm{d}t,\end{array}$$(24)

$$\begin{array}{c}\hfill f_{\mathit{w}}(t)=Heav(t){(2\pi i)}^{-1}{\displaystyle {\int }_{{\zeta }_0+i\infty }^{{\zeta }_0+i\infty }}\widehat{\mathcal{L}}{f}_{\mathit{w}}(\omega )exp(\omega t)\mathrm{d}\omega ,\end{array}$$

$$\begin{array}{c}\hfill {\varrho }_{\mathit{w}}(t)={(2\pi i)}^{-1}{\displaystyle {\int }_{{\zeta }_0-i\infty }^{{\zeta }_0+i\infty }}exp(\omega t)\frac{\widehat{\mathcal{L}}{(}^{[1]}\mathcal{D}{f}_{\mathit{w}})(\omega )}{1-\widehat{\mathcal{L}}{(}^{[2]}{\mathcal{D}}_{\mathit{w}})(\omega )}\mathrm{d}\omega ,\end{array}$$(25)

$$\begin{array}{c}\hfill \widehat{\mathcal{L}}{(}^{[k]}{\mathcal{D}}_{\mathit{w}})(\omega )={\displaystyle {\int }_0^{\infty }}exp{(-\omega t)}^{[k]}{\mathcal{D}}_{\mathit{w}}(t)\mathrm{d}t.\end{array}$$

The integral is taken in the right half–plane of the complex plane ω = −ζ − i∞, along a line parallel to axes ℜ(ω) = 0 and Heav(t), where is the Heaviside function. The poles of the integrand are determined by the equation $\widehat{\mathcal{L}}{(}^{[2]}{\mathcal{D}}_{\mathit{w}})=1$, which is reduced to the form (for ℜ(ω) > 0)

$$\begin{array}{c}\hfill \gamma \mathcal{H}(w_{\ell }^{\mathit{Re}},w_{\ell }^{\mathit{Im}}){\displaystyle {\int }_{{\mathbb{R}}^3}}\frac{i\mathit{w}{\mathrm{\nabla }}_{\mathit{v}}{F}_0}{\omega +i\mathit{w}\mathit{v}}=-1.\end{array}$$(26)

This is the desired dispersion equation in the general case, which, of course, can only be solved numerically. However, in the limiting cases ω = ω(w), it is fairly easy to determine using asymptotic estimates. In particular, for the small wave-vector modulus, we can expand it in a Taylor series of |w| (up to the first nonvanishing term) the factor in the integrand (ω/(ω + ivw)≈w/ω − ivw²/ω²), and choose the coordinate system with the abscissa axis aligned with the vector w (to get a scalar dependency).

Then

$$\begin{array}{c}\hfill \gamma \mathcal{H}(\mathit{w})\frac{{\mathit{w}}^2}{{\omega }^2}{\displaystyle \int }[B{T}^{-1}(-m{v}_1+{\mathrm{\Phi }}_{\min }^{\prime }·(t-t^{\prime }))\end{array}$$(27)

$$\begin{array}{c}\hfill exp(-\frac{m{v_1}^2}{2T}-\frac{{\mathrm{\Phi }}_{-}(t-t^{\prime })}{T}){]}_3\mathrm{d}{v}_1=1.\end{array}$$

Thus, we obtain an approximate dispersion relation (integrated over time), which can be briefly written as

$$\begin{array}{c}\hfill {\omega }^2=\gamma {\mathit{w}}^2\mathcal{H}(\mathit{w}){\displaystyle \int }[\dots {]}_3\mathrm{d}{v}_1.\end{array}$$(28)

The quantity ℋ(w^Re, w^Im) > 0, when the dynamical system depends on the sign of the expression in brackets […]₃, if […]₃ < 0, then, ω ∈ ℂ (a purely imaginary number), otherwise, ω is real. This means that the perturbation in our system, ∼exp(ωt + iwv), will have an oscillatory character in time if the system is in a conditional equilibrium at high wavelength and when criterion (28) is satisfied for imaginary ω. For real frequencies, the primary perturbation is fading or growing (unstable) depending on the sign ω ≶ 0. In addition, the integrand can change the sign, passing from the oscillatory behavior of the system to the osculatory one; this is affected by the increase in observation time, as well as at the transition to another branch of the multivalued function Φ_min.

For the case of the d = 2 system, the situation is simplified radically. This is due to the nature of the solutions of the Liouville–Gelfand equation (LGE) in the two-dimensional case: the number of solution branches varies depending on the λ parameter (the same as in the three-dimensional case), but in that case it is either two (for λ < λ_crit) or zero – when λ > λ_crit and there are no regular solutions.

For the critical value of the solution parameter of LGE, one can consider the corresponding pair (λ_crit, Φ_crit) as a singular point of the potential; however, the set of changes in the λ parameter do not always include the critical value. It is essential that we know the explicit form of the LGE solutions to therefore be able to move away from the local description of the system: neighborhoods of libration points in the three-dimensional system are spatially separated areas, and the above mathematical apparatus mainly focused on the application of limited systems. Naturally, one can move away from the mechanical equilibrium and turn to the thermodynamic description of the state systems; however, in an explicit form, this is extremely time consuming due to writing the potential in implicit representation and its near–equilibrium thermodynamics. The latter can make the analysis of physical applications of the kinetics of the system difficult at significant deviations from the initial position taken as equilibrium one. However, for the two-dimensional system, the situation is different: as previously mentioned, the gravitational potential of the N–particle system can be represented in a visible form, and we actually always use its uniqueness (the minimal solution of LGE).

The nonlinear Poisson Eq. (6) for d = 2 has the form

$$\begin{array}{c}\hfill {\mathrm{\Delta }}^{(2)}\varphi (\mathit{x})=4\pi G_2{m}^{2}Nexp(-\varphi (\mathit{x})/T),\end{array}$$(29)

$$\begin{array}{c}\hfill \varphi (\mathit{x})=2G_2{m}^{2}N{\displaystyle \int }F({\mathit{x}}^{\prime },{\mathit{v}}^{\prime })ln(|\mathit{x}-{\mathit{x}}^{\prime }|/\mathfrak{L})\mathrm{d}{\mathit{x}}^{\prime }\mathrm{d}{\mathit{v}}^{\prime }.\end{array}$$

In this case, the asymptotic condition is lim_|x|→∞|ϕ(x)−2G₂Nm²ln |x|)| = 0. Here G₂ is the gravitational constant of the two–dimensional model, and 𝔏 is the dimensionless factor, which we omit below. If you apply this to the model representation of our system isolated from external mass flows and preserve the global temperature regime (in the three–dimensional case, we imply a similar situation), then in determining the conditional extremum (N = const₁, ∑ε_k = const₂), the Gibbs entropy S(F) of the system on the plane in closed region $\overline{\mathrm{\Omega }}$ ($\mathit{x}\in \overline{\mathrm{\Omega }}$, $\overline{\mathrm{\Omega }}:|\mathit{x}|⩽R$) corresponding to the maximum of S(F) is the particle distribution function

$$\begin{array}{cc}& F_m(\mathit{x},\mathit{v})=\frac{a_1}{a_2+a_3·{|\mathit{x}|}^2}exp(-a_4{\mathit{v}}^2)\stackrel{R\to \infty }{\to }\hfill \end{array}$$(30)

$$\begin{array}{cc}& \frac{exp(2E/N^2-2)}{{\pi }^2}(exp(2E/N^2-2)+{\mathit{x}}^2{)}^{-1}exp(-{\mathit{v}}^2/N),\hfill \end{array}$$(31)

$$\begin{array}{cc}& E=\frac{N}{2}{\displaystyle \int }(m{\mathit{v}}^2+\varphi (\mathit{x}))F(\mathit{x},\mathit{v})\mathrm{d}\mathit{x}\mathrm{d}\mathit{v},\hfill \end{array}$$

$$\begin{array}{cc}& a_1=\frac{G_2{m}^{2}N/T}{2{\pi }^2}(2-G_2{m}^{2}N/T)\frac{R^2}{2},a_2=(2-{\gamma }_2{m}^{2}N/T)\frac{R^2}{2},\hfill \end{array}$$(32)

$$\begin{array}{cc}& a_3=G_2{m}^{2}N/2T,a_4=\frac{m}{2T}.\hfill \end{array}$$

Such a minimal state is realized (uniquely) only if T > T_c, T_c = G₂m²N/2.

Thus, substituting the above formulas into the system of Eqs. (2) and (3), taking the replacement in (4) of the integral kernels Y₃ → Y₂(x, x′) = 2G₂mln(|x − x′|) into account, we can repeat all the calculations of this section for the two–dimensional case, taking the Maxwell–Boltzmann distribution F_m(x, v) for a conditionally equilibrium state. This allows us to perform explicit differentiation in square brackets […]₁ in the definition of ^[2]𝒟_w(t−t′). This leads to the possibility of explicitly observing the change of sign in brackets […]₃ in the integrand on the left-hand side dispersion relation

$$\begin{array}{c}\hfill -\frac{{\mathrm{\Pi }}_1([a_k],v_1^j,\mathrm{\Delta }t)}{{\mathrm{\Pi }}_2([a_i],v_1,\mathrm{\Delta }t)}{|}_{i,k=1,\dots 4;j=1,\dots 3}\times a_{1}exp(a_4{v}_1^2)\equiv \mathfrak{D},\end{array}$$(33)

where ${\mathrm{\Pi }}_{1,2}([a_k],v_1^j,\mathrm{\Delta }t)$ are polynomials in variable v₁ (depending on parameters a_k,  Δt).

We emphasize a significant difference from the case d = 3: we are now considering the only minimal branch of the inter-particle gravitational potential; therefore, a slight change in the parameter norms does not lead to a qualitative change in the behavior of the system (far from the singular points’ dispersion equation). At the same time, for three–dimensional gravitational potential under conditions of polysemy of the function Φ_min for a fixed λ, in principle, a situation of branch inversion arises (the potential jumps move to a new stable position corresponding to a more energetically favorable state of the system), when the norm of the particle distribution function increases (in nonstationary systems with matter inflow). In addition, when we use the Fourier–Laplace transform in the three-dimensional case, we can see that the values of w^Im are allowed only in the lower half–plane. Otherwise, for the convergence of the integral on the infinite interval in the expression for ℋ(w) (solutions are only equivalent to damped density waves in space), and for two–dimensional systems, it is enough to use the usual Fourier transform. In this case the size of the system can be arbitrarily large and the periodicity of the system is violated only when changing the sign of the integrand in the dispersion relation.

3. The unperturbed problem

In the previous section, we considered the situation when a system of particles, initially in a conditional equilibrium at some point in time, undergoes a perturbation and the explicit expression was obtained for the dispersion relation in the case of long-wave perturbations. It is possible to analyze the unperturbed case in infinite or semi-infinite time intervals. At the same time, however, the conditions for the formation of inhomogeneities in space and the solution of the Vlasov–Poisson system becomes “quasi-discrete.” Namely, randomly distributed, intersecting one-dimensional “channels” do arise, consisting bunches of particles separated by voids, as the intersecting three-dimensional structure of the cosmic web (Gurzadyan et al. 2022).

Let us turn to a particular solution of the Vlasov–Poisson equation (for the three-dimensional case) of the form

$$\begin{array}{c}\hfill f^{+}(\mathit{x},\mathit{v},t)\sim c(k)f_{\mathit{w}}^{†}(\mathit{v})exp(i\omega t-i\mathit{w}\mathit{x}).\end{array}$$(34)

Let us take as a simplifying assumption that, in the expression

$$\begin{array}{c}\hfill F_0(\mathit{x},\mathit{v})=Bexp(-\frac{m{\mathit{v}}^2}{2T}-\frac{1}{T}{\mathrm{\Phi }}_{\min }(\mathit{x})),\end{array}$$(35)

the norm of quantity Φ_min ≈ φ = const, that is when the state of the system is close to relative equilibrium, in which by definition ${\mathrm{\Phi }}_{\text{min}}^{\prime }=0$.

We, as before, follow the analysis methodology proposed in Vlasov (1950). As a result, we have

$$\begin{array}{c}\hfill f_{\mathit{w}}^{†}(\mathit{v})=\frac{G}{m}\mathcal{H}(\mathit{w})\frac{(\mathit{w}{\mathrm{\nabla }}_{\mathit{v}}\mathfrak{M}(\mathit{v}))}{\mathit{w}\mathit{v}-\omega }{\displaystyle \int }{f}_{\mathit{w}}^{†}(\mathit{v})\mathrm{d}{v}_1\mathrm{d}{v}_2\mathrm{d}{v}_3,\end{array}$$(36)

$$\begin{array}{c}\hfill \mathfrak{M}(\mathit{v})=Bexp(-\phi /T)exp(-m{\mathit{v}}^2/2).\end{array}$$

The condition for the existence of a nontrivial solution as of functions f⁺ is obtained after integrating over the velocities of both parts (dispersion relation ω = ω(w))

$$\begin{array}{c}\hfill G{m}^{-1}\mathcal{H}(\mathit{w})Bexp(-\phi /T){\displaystyle \int }\frac{(\mathit{w}{\mathrm{grad}}_{\mathit{v}})}{\mathit{w}\mathit{v}-\omega }\mathrm{d}\mathit{v}=1.\end{array}$$(37)

The values of the components of the wave vector w are complex and the integral is taken in the sense of the principal value symmetrically along the axis v₁, if the w axis is along the axis of the radial coordinate. The sum of two terms when going around the pole corresponds to the retarded and leading potentials.

Hence, we have

$$\begin{array}{c}\hfill f_{\mathit{w}}^{†}(\mathit{v})=c(\mathit{w})G{m}^{-1}\mathcal{H}(\mathit{w})Bexp(-\phi /T)(\mathit{w}{\mathrm{\nabla }}_{\mathit{v}}/(\omega -\mathit{w}\mathit{v})).\end{array}$$(38)

One can obtain an analytical representation of the dispersion relation if we approximate the Maxwellian 𝔐

$$\begin{array}{c}\hfill \underset{|\mathit{v}|\ll \infty }{lim}\mathfrak{M}\to N\sqrt{2T/m}/\pi (2T/m+v_1^2{)}^{-1}\equiv {\mathfrak{M}}^{†}.\end{array}$$(39)

Then

$$\begin{array}{c}\hfill {\displaystyle \int }\mathrm{d}{v}_1\frac{\mathrm{d}{\mathfrak{M}}^{†}/\mathrm{d}{v}_1}{m(v_1-\omega /w)}=-\frac{N}{\pi T}\frac{\pi }{2}\frac{1-{\nu }^2}{{(1+{\nu }^2)}^2},\nu =\omega /w\sqrt{m/(2T)}.\end{array}$$(40)

Therefore, the dispersion relation takes the form

$$\begin{array}{c}\hfill (\frac{4\pi }{m}B{\displaystyle {\int }_0^{\infty }}exp(w^{\mathit{Im}}\xi )\frac{sin(w^{\mathit{Re}}\xi )}{w^{\mathit{Re}}\xi }\xi \mathrm{d}\xi {)}^{-1}=-\frac{1-{\nu }^2}{{(1+{\nu }^2)}^2}.\end{array}$$(41)

Similarly, one can obtain the dispersion relation for any function 𝔐(v₁). So, we got a ω(w) dependency for the system of particles with self-gravity. This provides information about the nature of the dynamics of the system, that is it is oscillatory in time and/or space with decreasing or conserving amplitude being exponentially decreasing and unstable.

The importance of the analysis of dispersion solutions lie in the fact that no random fluctuations have a role in the mechanism of emergence of pseudo-oscillatory (in time and in selected spatial directions) structures of matter. While explicitly expressed spatial periodicity of structures does not occur, the dispersion relations reveal certain ordering (sequencing) of the structures.

Thus, both the emergence of macro-structures, and their distribution is immanently inherent in the technique using the Vlasov-Poisson equations at the potential of Eq. (1).

4. Conclusions

The emergence of two-dimensional structures as Zeldovich pancakes is currently associated with density perturbations, whose growth is described by classical or weakly relativistic hydrodynamics, without considering the role of gravity. Then, certainly the problem is to ensure the emergence of the needed perturbations.

The issue we dealt with in this study is whether a kinetic description – taking the role of self-consistent gravity of the many-particle ensemble into account, that is to say within the formalism of the Vlasov equation – can also lead to similar filamentary structures.

Within the Vlasov formalism (Vlasov 1950, 1967), we analyzed a kinetic model and revealed the occurrence of semi-periodicity of the arisen structures – voids – separated by two-dimensional surfaces (walls). On the one hand, the emergence of two-dimensional structures have clear differences from the hydrodynamical approach, since is not depending on the form of the initial perturbations, but it is due to the self-consistent gravitational interaction only. On the other hand, the emergence of two-dimensional structures (walls) in the three-dimensional N-body problem is not trivial a priori.

The crucial point in our kinetic approach is the use of the potential of Eq. (1) with the repulsive cosmological Λ term. The latter does not include any free parameters (e.g., scalar fields, etc); however, as outlined above, it is based on one of the very first principles, the generalized form of a potential derived in Gurzadyan (1985) to ensure an identical gravity for the sphere and point mass. As previously stated, this potential already has certain empirical bases, enabling one to describe the dynamics of galaxy clusters, and the Hubble tension as a result of local and global flows of galaxies. We now see that the same approach within the Vlasov formalism is able to predict certain features of the cosmic structures, for example of voids.

Of course, the important point for any theory or model is the study of the correspondence – qualitative and whenever possible, quantitative – of the predictions with the observations, as for example the analysis in Bouchè et al. (2022) of the predictions of the nonlocal gravity model with gravitational lensing observations of the CLASH survey (Postman et al. 2012).

In our case, qualitative correspondence can be considered as the prediction of two-dimensional structures (walls) and voids between them, obtained at the rigorous analysis of the three-dimensional system. This nontrivial result is due to the self-consistent interaction being considered which has both attractive and repulsive terms.

As for the quantitative criterion, we have derived the scale (size) of the voids, namely, that scale is determined by the critical radius r_c = (3Gm/Λc²)^1/3, which is balancing the Newtonian and repulsive terms in Eq. (1). This enables a direct comparison with observations, for example the size of voids yields around 25 Mpc when using the data of the Virgo Supercluster, which is in agreement with observational data (e.g., Ceccarelli et al. 2006; Gurzadyan & Kocharyan 2009; Samsonyan et al. 2021; Nadathur et al. 2020). As mentioned in Gurzadyan & Stepanian (2021a), the observational data defining r_c can vary from survey to survey, thus indicating that the numerical scale of the voids can vary depending on the mean matter density of the local region of the Universe. In other words, although the mechanism of the void formation is common – the interplay between the self-consistent contribution of gravitational attraction and the repulsion with the cosmological constant – the voids can have a different mean size determined by the local conditions (density) of the Universe. This can become a vast topic for an analysis involving observational surveys.

Our analysis thus indicates the important role that the self-interaction can have in the emergence of the voids in the large-scale matter distribution, thus supporting further studies in this direction.

Acknowledgments

We thank the referee for valuable comments.

References

All Figures

	Fig. 1. Modified gravitational potential ΦGN(r).
In the text