© The Authors 2026

1. Introduction

The cosmological tensions that have recently emerged, most notably the Hubble tension (Riess 2020; Riess et al. 2024a,b, 2025; Dainotti & De Simone 2025; Leauthaud & Riess 2025) and the baryon acoustic oscillation signature by DESI DR2 (DESI Collaboration 2025a,b), outline the possible genuine differences in the description of the early and late Universe. The theoretical approaches that address the tensions span a broad spectrum of issues, from the need for new physics to particular models of evolving dark energy and modified gravity, amongst other issues (Capozziello et al. 2024; Alfano et al. 2026; Chaudhary et al. 2025). The origin and the evolution of large-scale low-dimensional structures, such as the long-known void walls and filaments of various scales, are among the goals of the theoretical approaches. The Zeldovich pancake theory (Zeldovich 1970; Arnold & Shandarin 1982; Shandarin 1989) describes the evolution of the primordial density perturbations in hydrodynamical approximation and predicts the formation of the cosmic web on cosmological scales.

The formation mechanisms of the filaments in the early and late Universe can have specific differences that distinguish their features on various scales. The role of the self-consistent interaction in the structure formation in the local Universe was considered in (Gurzadyan et al. 2022, 2023a,b, 2025). The principal aspect in those studies was the consideration of the role of repulsive interaction due to the cosmological constant in the local scales. The cosmological constant in the local scales emerges if based on a theorem (Gurzadyan 1985) that states the most general function for the force satisfying the identity of the gravitational field of a sphere and of a point mass, i.e. when the sphere and the point mass have an identical influence on a test particle. That general function has the form (Gurzadyan 1985)

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

where the second term in the right-hand side marks the cosmological constant term in weak-field General Relativity and McCrea-Milne non-relativistic cosmology (McCrea & Milne 1934; Zeldovich 1981). That second term (i.e., the cosmological term) does not change the O(4) symmetry of the Newtonian field. It is notable that the general function does not contain any other terms besides the second, the cosmological term. Also note that this function does not satisfy a force-free condition inside a spherical shell as distinct from Newtonian gravity. On this point, we mention the observational indications on the influence of galactic halos on the properties of the disks of spiral galaxies (Kravtsov 2013), which supports the presence of a force field inside a shell as predicted by Eq. (1). Then, Eq. (1) enables us to describe the Hubble tension as a result of two flows with non-identical Hubble parameters: a local flow and a global flow. The following two equations describe the two flows (Gurzadyan & Stepanian 2021a,b):

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

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

The first equation follows from Eq. (1) and includes the local mean density, ρ_local, as a parameter; the second equation is the Friedmann equation with the global mean density, ρ_global. The difference in the densities leads to the difference in the relevant Hubble parameters and thus provides an explanation for the Hubble tension. Moreover, as shown in Gurzadyan & Stepanian (2021b), one can derive absolute constraints on the lower and upper values for the local Hubble parameter:

$$\begin{array}{c}\hfill \sqrt{\mathrm{\Lambda }{c}^2/3}\simeq 56.2<H_{\mathrm{local}}<\sqrt{\mathrm{\Lambda }{c}^2}\simeq 97.3({\mathrm{km}/\sec \mathrm{Mpc}}^{-1}),\end{array}$$

which is in agreement with the observational data. This can be considered as an empirical support to the validity of the nonrelativistic description of the local Universe, as outlined in Zeldovich (1981). As shown in Gurzadyan & Stepanian (2018), Gurzadyan (2019), Gurzadyan & Stepanian (2019, 2020), Eq. (1) fits the observational data on the dynamics of galaxy pairs, groups, and clusters. We mention the recent efforts to test the law of Eq. (1) by means of quantum technologies (Fernandez-Melendez et al. 2025).

Although Debye screening is absent in gravitational systems and their problems have to be considered with different methods (e.g., Gurzadyan & Savvidy 1986), plasma theory has a well-developed and efficient mathematical apparatus for analyzing wave motions of various types that are suitable for adaptation to gravitational systems. The methods developed in plasma physics have been applied for certain problems of stellar dynamics (Lynden-Bell 1962; Lynden-Bell et al. 1994), including those associated with Landau damping for small perturbations and with Bernstein–Greene–Kruskal waves (Fridman & Polyachenko 1984; Saslaw 1985; Palmer 1994; Vandervoort 2003; Polyachenko et al. 2021). In Lau & Binney (2021a,b) the authors point out the possibility of using van Kampen wave methods for the large-scale motion of clusters and galaxies.

In this paper, we aim to describe large-scale structures using non-dissipative solutions of the Vlasov–Poisson equations of the van Kampen wave type. The periodicity of the waves is violated when taking into account the repulsive force due to the inclusion of a cosmological term in the consideration, since in this case our system is locally close to weakly inhomogeneous (for a long-range order, it is significantly inhomogeneous). Including the influence of the Λ-term in the modified Poisson equation in the analysis of the behavior of a N-particle system follows from the above-mentioned theorem (Gurzadyan 1985) and Eq. (1). We also consider the possibility of introducing Bernstein–Greene–Kruskal waves (Bernstein et al. 1957; Montgomery 1960) as structural units of cosmological systems. For substantially inhomogeneous systems, we analyze the possibility of a smooth transition if we account for the field of gravitational disturbances to the normal mode method. Thus, aperiodicity emerges as a characteristic feature of the local filaments.

2. Kinetic equations: Linearization

We consider a set of N cosmological objects, “particles” with masses m_{i = 1, …, N} = m ≡ 1, which interact gravitationally. Then the system of Vlasov–Poisson equations for description of its dynamics is represented as

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

$$\begin{array}{cc}\hfill & {\mathrm{\Delta }}_{\mathit{x}}\mathrm{\Phi }[F(\mathit{x})]=4\pi AN\gamma {\displaystyle \int F}(\mathit{x},\mathit{v},t)d\mathit{v}-c^2\mathrm{\Lambda },\hfill \end{array}$$(5)

where F(x, v, t) is the distribution function of gravitationally interacting particles, A is a normalization factor for particle density, and γ is the gravitational constant. The system of particles is situated in a large domain of configurational space Ω ⊂ ℝ_x³ (diam Ω ≡ R_Ω < ∞). The nonlinear Poisson Equation (5) takes the form of an inhomogeneous Liouville–Gelfand equation (Dupaigne 2011) with a local (kinetic) temperature (Vlasov 1966), when using a more general form for the Poisson equation.

Equation (5) is the nonlinear Poisson equation for Newton–type gravitation. The third term on the right hand side of the kinetic Equation (4) may be represented as

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

$$\begin{array}{cc}\hfill \mathrm{\Phi }[F(\mathit{x},t)]& =4\pi AN\gamma {\displaystyle \int \int {\mathfrak{K}}_3}(\mathit{x}-{\mathit{x}}^{\prime })F({\mathit{x}}^{\prime },{\mathit{v}}^{\prime },t)d{\mathit{x}}^{\prime }d{\mathit{v}}^{\prime }\hfill \\ \hfill & +\frac{\mathrm{\Lambda }{c}^2}{6}{|\mathit{x}|}^2+{\widehat{\mathfrak{B}}}_{\partial \mathrm{\Omega }}(\mathit{x},{\mathit{x}}^{\prime }),\hfill \end{array}$$(7)

where 𝔎₃(x − x′) = −|x − x′|⁻¹ (Newtonian interaction kernel) and ${\widehat{\mathfrak{B}}}_{\partial \mathrm{\Omega }}(\mathit{x},\mathit{x}\prime )$ is an operator term that takes into account the influence of the boundary conditions (we take into account the influence of this term by setting the appropriate boundary conditions). Classical Newtonian potential Φ_N(r) =  − γM/r increases monotonically on the interval r ∈ (0, +∞) (Φ_N ∈ ( − ∞, 0)), while the potential of Eq. (1), including a cosmological term Φ_GN(r)≡ − GM/r − c²Λr²/6, increases on the interval r ∈ (0;  r_c] and decreases on the interval r ∈ (r_c;  ∞), where r_c = (3GM/(Λc²))^1/3.

We consider the nonstationary case of dynamics F = F(x, v, t). In previous publications (Gurzadyan et al. 2022, 2023a,b, 2025) we focused on the possibility of transition to the integral form of the equation and the formulation of a boundary value problem of the Dirichlet type (for the gravitational potential) with the aim of determining the Green’s function of the problem for substitution into the kernel of the Hammerstein operator. However, for a nonstationary system of equations for the evolution of a cosmological system of particles in the self-consistent approximation (4)–(5), the main role is played by the formulation of the initial problem for the Vlasov equation. In this case, a direct derivation of the solutions and their study via analytical methods is complicated (Chandrasekhar 1960). In the present work, we restrict ourselves to the study of the properties of the solutions of the linearized version of the Vlasov–Poisson system for gravity for the potential of Eq. (1).

The linearization of the Vlasov equation is quite nontrivial, since its result depends significantly on the type of gravitational field and, due to the self-consistency of the problem, this depends on the distribution function of particles in the system. From a physical point of view, it is natural to single out a stationary homogeneous solution when the distribution function does not depend on the coordinates F = F_M(v; T) or, in a more general case, F = F₀(v), F₀ ∈ C¹ ∩ L²(Ω_v),Ω_v ⊂ ℝ³_v; it corresponds to the point at which the total force acting on the particle is zero, that is, the total potential of the gravitational attractive and repulsive forces is constant. Within the framework of the theorem (Gurzadyan 1985), one can take into account the presence of the previously mentioned local maximum of the two-particle potential. Then, the region near the equilibrium point, denote it x₀, in the interaction channel between two distant external masses of subsystems includes the complete system of gravitating particles (we have a state of unstable equilibrium).

For a broader class of problems, it becomes necessary to consider a more general type of linearization near the equilibrium Maxwell–Boltzmann solution of the stationary Vlasov–Poisson system in the form F_MB ∼ exp( − 𝔈(x, v, t)/T), including the (two-particle) potential, 𝔈 = mv²/2 + Φ(x, t₀), for a fixed instant of time. The dual solutions of the Poisson equation Φ(x, t₀) and the gravitational field strength are expressed through solutions of the Volterra equations of the second kind, and therefore classical dispersion relations for the Vlasov equations cannot be obtained.

It is necessary to note the meaning of the temperature, T (kinetic temperature), in the solutions of the kinetic equation, taking into account the action of the Λ-term in the Poisson equation. Particle density on the right side of the Poisson equation can be expressed in terms of the nonstationary solution of the Vlasov equations. In the simplest case, this solution is identical to unimodal Maxwell distributions; in the general case, one can consider, for example, representing F₀ as a multimodal set of Maxwellians with different amplitudes. However, the physical meaning of the equilibrium (nonuniform stationary states) solution of the Vlasov equation is essentially different from that of the Boltzmann equation. This solution must meet the following requirements: 1)the maximum possible statistical independence, 2)isotropy of the velocity distribution, 3)stationarity of 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 of this 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}$$(8)

and we get a system of ordinary differential equations (ODE)

$$\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}=-T^{-1},\end{array}$$(9)

where T is a constant of separation of the variables; its physical meaning is kinetic temperature in the system of interacting collisionless particles in accordance with Vlasov’s definition (Vlasov 1966, 1978), as collisional equilibrium is globally absent in this system.

Equation (5) for gravitational potential can be written as

$$\begin{array}{cc}\hfill \mathrm{\Delta }\mathrm{\Phi }(\mathit{x})& ={\lambda }^{†}exp(-\mathrm{\Phi }(\mathit{x})/T)-c^2\mathrm{\Lambda },{\lambda }^{†}=4\pi \gamma N{A}_T,\hfill \\ \hfill A_T& \equiv {\rho }_0{\displaystyle \int exp}(-m{v}^2/(2T))v^{2}dv,{\rho }_0=(\frac{m}{2\pi T}{)}^{3/2}.\hfill \end{array}$$(10)

The last equation can be rewritten in the form

$$\begin{array}{cc}\hfill \mathrm{\Delta }W(\mathit{x})& ={\lambda }^{♯}exp(W(\mathit{x})),W(\mathit{x})\equiv -\frac{\mathrm{\Phi }(\mathit{x})}{T}-\frac{c^2\mathrm{\Lambda }{\mathit{x}}^2}{6T},\hfill \\ \hfill {\lambda }^{♯}& =-\frac{{\lambda }^{†}}{T}exp(c^2\mathrm{\Lambda }{\mathit{x}}^2/T).\hfill \end{array}$$(11)

Solutions of the equation ΔW(x) =  − ζexp(W(x)) (ζ ∈ ℝ₊¹) in the 3–dimensional case are radially symmetric (W = W(|x|) by the Gi–Nidas–Nirenberg theorem; Dupaigne 2011) and are unstable with respect to the pre-exponential parameter: their existence and number depend on the value of the parameter ζ. According to Bebernes & Eberly (1989), the solution of the standard Dirichlet problem for it has a structure that can be described as follows. Let ζ_crit = 2 (if the boundary value problem is considered on the reduced interval |x|≡r ∈ [0; 1]); then we have ζ_FK > ζ_crit such that: 1)for ζ = ζ_FK, there is a unique solution (W_FK); 2)for ζ > ζ_FK, there are no solutions; 3)for ζ = ζ_crit, there is a countable infinity of solutions ($W_{\mathrm{crit}}^{(n)}$, n ∈ 𝔙, card(𝔙) = ℵ₀); 4)for ζ ∈ (0, ζ_FK)∖{ζ_crit}, there is a finite number of solutions (W_K^(k), k ∈ {1, 2, …, K},K ≥ 1). Since it is possible to uniquely (for fixed parameters T, N) compare the values of the function −λ^♯(|x|) with the values of the parameter ζ, it can be stated that with an increase in the modulus of the radius vector |x|, three regions of solutions to Equation (11) arise: the region of uniqueness of solutions X₁(x) = {|x|< X^(I)}∪{X^(III)}, the region of multivalued solutions (differing in norm) X₂(x) = {X^(I) < |x|< X^(II)}, and the region of absence of solutions X₃(x) = {|x|> X^(III)}.

Let us consider the linearization of Equation (10) in the neighborhood of the solution W(x) analytic solution to which W can be associated to solution (10) as W + 𝔴(x, t) (∥𝔴∥≪∥W∥ by the chosen norm and, correspondingly, exp(𝔴)≈1 + 𝔴). We obtain the linear Poisson equation

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathfrak{w}(\mathit{x})={\lambda }^{♯}exp(W(\mathit{x}))·\mathfrak{w}(\mathit{x}).\end{array}$$(12)

Obviously, in the neighborhood of x₀ the last equation is simplified, since the gravitational field of a point with an equivalent total mass and cosmological repulsion allows us to set Φ(x₀)≡ − TW(x₀) = Φ⁽⁰⁾(=const). Thus, the equation for the potential perturbation in the above neighborhood O(x₀) takes the form (K⁽⁰⁾ = exp(W(x₀)))

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathfrak{w}={\lambda }^{♯}({\mathit{x}}_0)K^{(0)}·\mathfrak{w}(\mathit{x}).\end{array}$$(13)

The linearization of the Vlasov equation itself is performed (in the simplest case considered) in the neighborhood of the equilibrium function F_M(v) or, in a more general case, F₀(v) with several maxima, which is realized, for example, in the case of co-directional particle beams; so we have

$$\begin{array}{c}\hfill F\to F_{\mathit{MB}}(\mathit{x},\mathit{v})+\stackrel{\sim }{f}(\mathit{x},\mathit{v},t),\end{array}$$(14)

where the perturbation $\stackrel{\sim }{f}(\mathit{x},\mathit{v},t)$ is related to the Poisson Equation (12) with an exponential dependence of the parameter on the spatial variable. Eliminating the quadratic terms a small addition $\stackrel{\sim }{f}$ gives us

$$\begin{array}{c}\hfill \frac{\partial \stackrel{\sim }{f}}{\partial t}+\mathit{v}·{\mathrm{\nabla }}_{\mathit{x}}\stackrel{\sim }{f}-{\mathrm{\nabla }}_{\mathit{v}}{F}_0·{\mathrm{\nabla }}_{\mathit{x}}\varphi [\stackrel{\sim }{f}](\mathit{x},t)=0,-T(W+\mathfrak{w})=\mathrm{\Phi }+\varphi .\end{array}$$(15)

Next, we consider the methodology for studying the linear system of Vlasov-Poisson equations using “normal modes” and the use of the transition to the space of distributions. This allows us to study analogs of the attenuation of Landau waves and longitudinal van Kampen density waves for a system of gravitating particles.

3. Van Kampen modes versus self-consistent gravitational potential with a cosmological constant

Let us consider the invariant properties (independent of solutions) of the linearized Vlasov–Poisson system of Equations (13)–(14). First, we consider the case of a gravitational field strength that corresponds to a local neighborhood of the extremum of the self-consistent potential, taking into account the action of the cosmological term

$$\begin{array}{cc}& \frac{\partial \stackrel{\sim }{f}(\mathit{x},\mathit{v},t)}{\partial t}+\mathit{v}\mathrm{\nabla }\stackrel{\sim }{f}(\mathit{x},\mathit{v},t)={\mathrm{\nabla }}_{\mathit{x}}\varphi (\mathit{x},t)·{\mathrm{\nabla }}_{\mathit{v}}{F}_0(\mathit{v}),\hfill \\ \hfill & {\mathrm{\nabla }}_{\mathit{x}}\varphi (\mathit{x},t)={{\lambda }^{♯}}_0{K}^{(0)}{\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathit{x}\prime }}{\int }_{{\mathrm{\Omega }}_{\mathit{v}}}}{\mathrm{\nabla }}_{\mathit{x}}\frac{\stackrel{\sim }{f}(\mathit{x}\prime ,\mathit{v},t)}{|\mathit{x}-\mathit{x}\prime |}d\mathit{v}d\mathit{x}\prime .\hfill \end{array}$$(16)

We represent $\stackrel{\sim }{f}(\mathit{x},\mathit{v},t)$ via the van Kampen ansatz or “normal modes” (Van Kampen & Felderhof 1967; Holloway & Dorning 1991): 𝔎(v)exp(ikx − iωt) (plane waves are eigenfunctions of the Laplacian from the left-hand side of the Poisson Equation (13)). We are interested in solutions of the system of Equations (13)–(14) in the form of longitudinal waves, therefore in the velocity space we choose axes that are parallel (z) and perpendicular (x, y) to the wave vector, k; then the longitudinal component of the velocity is v_∥ = v_z = e_k ⋅ v (where e_k = k/|k|), and the transverse component is, respectively: u = v − e_kv_∥. In this case, we can introduce distribution functions that only depend on one component of the velocity: $f(\mathit{k},v_{\Vert },t)={\displaystyle \int \stackrel{\sim }{f}}(\mathit{k},\mathit{v},t)\delta (v_{\Vert }-\mathit{k}·\mathit{v}/k)d\mathit{v}={\displaystyle \int \stackrel{\sim }{f}}(\mathit{k},\mathit{v},t)d\mathit{u}$.

We rewrite Equation (14) for these modes, freeing ourselves from the transverse velocity components (and discarding the tilde sign over f):

$$\begin{array}{cc}& (\omega -k{v}_{\Vert }){\displaystyle \int \mathfrak{K}}(\mathit{v})d\mathit{u}=\frac{4\pi }{k^2}\mathit{k}{\lambda }_0^{♯}{K}_{\mathrm{\Lambda }}^{(0)}{\displaystyle \int \frac{\partial F_0}{\partial v_{\Vert }}}d\mathit{u}{\displaystyle \int \mathfrak{K}}(\mathit{v}\prime )d\mathit{v}\prime ,\hfill \\ \hfill & -\widehat{\varphi }(\mathit{k},t)={\lambda }_0^{♯}{K}_{\mathrm{\Lambda }}^{(0)}{\displaystyle \int \int \frac{\mathfrak{K}(\mathit{v}\prime )(\mathit{x}-\mathit{x}\prime )}{{|\mathit{x}-\mathit{x}\prime |}^3}}exp(i\mathit{k}\mathit{x}\prime -i\omega t)d\mathit{v}\prime d\mathit{x}\prime ,\hfill \\ \hfill & {\displaystyle \int \frac{\mathit{x}}{{|\mathit{x}|}^3}}exp(i\mathit{k}\mathit{x})d\mathit{x}=4\pi i\frac{\mathit{k}}{k^2}.\hfill \end{array}$$(17)

We divide both sides of the last equation by (ω − kv_∥) and integrate with respect to the variable v_∥. The integral ∫𝔎(v′)dv′ (an unimportant constant) is canceled out, and we obtain a dispersion relation that is invariant with respect to the form of the solution of the kinetic equation

$$\begin{array}{c}\hfill 1-(\varkappa /k){\displaystyle \int \frac{d{f}_0}{d{v}_{\Vert }}}\frac{d{v}_{\Vert }}{\omega -k{v}_{\Vert }}=0,\varkappa =4\pi {\lambda }_0^{♯}{K}_{\mathrm{\Lambda }}^{(0)}.\end{array}$$(18)

If we do not consider the longitudinal velocity as distinguished, then the general form of the dispersion law has the form: D(k, ω)≡1 − 𝜘(k/k²)∫_L(F₀)_v′(ω − kv)⁻¹dv = 0 (normal modes will correspond to the case Re(ω(k)) ≫ Im(ω(k))).

We are interested in the possibility of obtaining a solution to the Vlasov–Poisson equations that is stable in time and associated with the simplest cosmological structures, i.e. those of low dimensionality. It can be obtained using normal modes in the form

$$\begin{array}{c}\hfill f(z,v_{\Vert },t)={\displaystyle \int \int \mathfrak{O}}(k,\nu )\mathfrak{N}(k,\nu ;v_{\Vert })exp(ikz-ik\nu t){|}_{\nu =\omega /k}dkd\nu ,\end{array}$$(19)

where 𝔒(k, ν) is an (admissible) function that corresponds to certain Cauchy data for the kinetic equation for perturbation f. If the initial condition is represented as f(z, v_∥, t = 0) = ∫g(k, v_∥)exp(ikz)dk, then, obviously, Equation (19) is reduced to the form ∫𝔒(k, ν)𝔑(k, ν; v_∥)dν = g(k, v_∥), and the variable k here acquires the meaning of a parameter.

For what follows, we return to Equation (17) and consider a nonobvious consequence of taking the integral of 𝔎 over the transverse velocities and dividing both parts by (ω − kv_∥). The result here must take into account the possibility of the equation solutions going into the space of generalized functions: as is known, for the functional equation (x − y)μ₁(x) = μ₂(x) (defined on the interval [x₁; x₂] of the real axis) and the point y ∈ (x₁; x₂), the solution must be interpreted as a distribution. This distribution can be written in the following form: ${\mu }_1(x|y)={\mu }_2(x)P.V.\frac{1}{x-y}+{\mu }^{♮}(y)\delta (x-y)$, where the Cauchy principal value in the form of a distribution is defined by the relation $(P.V.\frac{1}{x},\mu )=\underset{ϵ\to 0}{lim}{\displaystyle {\int }_{|x|\ge ϵ}}(\mu (x)/x)dx$), and μ^♮(y) is “the strength of the concentration” of the Dirac function at the point x = y determined from additional conditions imposed on the generalized function μ₁(x|y).

Thus, Equation (17) rewritten as

$$\begin{array}{cc}\hfill (\nu -v_{\Vert })\mathfrak{N}(v_{\Vert })& ={\nu }_{\varkappa }\mathfrak{F}(v_{\Vert }){\displaystyle \int \mathfrak{K}}(v_{\Vert }\prime )d{v}_{\Vert }\prime ,{\displaystyle \int \mathfrak{K}}(v_{\Vert }^{\prime })d{v}_{\Vert }^{\prime }=1,\hfill \\ \hfill \mathfrak{N}(v_{\Vert })& \equiv {\displaystyle \int \mathfrak{K}}(\mathit{v})d\mathit{u},\hfill \\ \hfill \mathfrak{F}(v_{\Vert })& \equiv {\displaystyle \int \frac{\partial F_0(\mathit{v})}{\partial v_{\Vert }}}d\mathit{u},\nu =\frac{\omega }{k},{\nu }_{\varkappa }=\frac{\varkappa }{k^2},\hfill \end{array}$$(20)

after dividing both sides of Equation (20) by (ν − v_∥) and should be written in the sense of distributions

$$\begin{array}{c}\hfill \mathfrak{N}(v_{\Vert })={\nu }_{\varkappa }·P.V.\frac{\mathfrak{F}(v_{\Vert })}{\nu -v_{\Vert }}+{\mu }^{♮}\delta (\nu -v_{\Vert }),(\nu -v_{\Vert })\delta (\nu -v_{\Vert })=0.\end{array}$$(21)

In this case, from the normalization condition in Eq. (20), the intensity value μ^♮ is determined by the condition of its agreement with Formula (21): ${\mu }^{♮}=1-{\nu }_{\varkappa }P.V.{\displaystyle \int (}\mathfrak{F}(v_{\Vert })/(\nu -v_{\Vert }))d{v}_{\Vert }$.

Let us substitute into the equation ∫𝔒(k, ν)𝔑(k, ν; v_∥)dν = g(k, v_∥) the value 𝔑(v_∥) from (21):

$$\begin{array}{c}\hfill \mathfrak{O}(k,v_{\Vert })(1-\pi {\nu }_{\varkappa }^2\widehat{\mathfrak{H}}\mathcal{F}(v_{\Vert }))-\widehat{\mathfrak{H}}\mathfrak{O}(k,v_{\Vert })\pi {\nu }_{\varkappa }^2\mathcal{F}(v_{\Vert })=g(k,v_{\Vert })\end{array}$$(22)

(k is still a parameter). Here, $\widehat{\mathfrak{H}}(\mathrm{\Psi }(x))=(1/\pi )P.V.{\displaystyle \int (}\mathrm{\Psi }(x\prime )/(x-x\prime ))dx\prime$ is the Hilbert transform, which is related to the Fourier transform of the function Ψ(x) = Ψ₊(x)+Ψ₋(x): ${\mathrm{\Psi }}_{+}(x)-{\mathrm{\Psi }}_{-}(x)=i\widehat{\mathfrak{H}}(\mathrm{\Psi }(x))$, where Y₊(x)≡∫₀^∞Y(q)exp(iqx)dq, Y₋(x)≡∫_−∞⁰Y(q)exp(iqx)dq; symbols Y, Y_± are used to denote functions $\mathrm{\Psi },\mathcal{F},\mathfrak{O},g,$ and their decompositions.

The last equation can be rewritten as

$$\begin{array}{cc}\hfill (1+2\pi i{\nu }_{\varkappa }^2& {\mathcal{F}}_{+}(v_{\Vert })){\mathfrak{O}}_{+}(k,v_{\Vert })+(1-2\pi i{\nu }_{\varkappa }^2{\mathcal{F}}_{-}(v_{\Vert })){\mathfrak{O}}_{-}(k,v_{\Vert })\hfill \\ \hfill & =g_{+}(k,v_{\Vert })+g_{-}(k,v_{\Vert })\equiv g(k,v_{\Vert }).\hfill \end{array}$$(23)

The terms on the left-hand side are analytic and have no singularities in the upper (Im(η) > 0) and lower (Im(η) < 0) parts of the complex (η_Re, η_Im)–plane (ℝ ∋ v_∥ → η ∈ ℂ), respectively, and also asymptotically tend to zero in their half-plane. The decomposition of g(η) into two functions with such properties is unique, and therefore (1 ± 2πiν_𝜘²ℱ₊(v_∥))𝔒_±(k, v_∥) = g_±. Therefore, if there is a solution to Eq. (22), then it must coincide with 𝔒 = 𝔒₊ + 𝔒₋, 𝔒_± = g_±/(1 + 2πiν_𝜘²ℱ_±) (the condition for this is 𝔑(v_∥)≠0, which is true, in particular, for the Maxwellian distribution). If we consider on the half-plane Im(η) > 0 a holomorphic and asymptotically close to unity function 𝒵(η) = 1 + 2πiν_𝜘²ℱ₊(v_∥), we can extend it to the half-plane Im(η) < 0: 𝒵(η) = 1 + 4π²iν_𝜘²z𝔑(η)+2πν_𝜘²∫η′𝔑(η′)/(η′−η)dη′. Now we can write out the final form of the solution to the initial value problem with the general solution (19):

$$\begin{array}{cc}\hfill f(z,v_{\Vert },t)& ={(2\pi )}^{-1}{\displaystyle \int \int \int \mathfrak{N}}(k,\nu ;v_{\Vert })exp(ik(z-z\prime )-ik\nu t)\hfill \\ \hfill & (f_{+}(z^{\prime },\nu ,t=0)/\mathcal{Z}(k,\nu )\hfill \\ \hfill & +f_{-}(z^{\prime },\nu ,t=0)/\overline{\mathcal{Z}(k,\nu )})dkdz\prime d\nu ,\hfill \\ \hfill & f_{+}(z,\nu ,0)+f_{-}(z,\nu ,0)=f(z,\nu ,0).\hfill \end{array}$$(24)

For the initial function of the form f(z, v_∥, t = 0) = ∫g(v_∥)exp(ikz)δ(k − k₁)dk (λ = 2π/k₁ = const), the density of particles in the disturbance wave is

$$\begin{array}{cc}\hfill {\varrho }_f(z,t)& =exp(i{k}_{1}z){\displaystyle {\int }_{\mathbb{R}}}exp(-i{k}_1\nu t)(g_{+}(v_{\Vert })/\mathcal{Z}(k_1,\nu )\hfill \\ \hfill & +g_{-}(v_{\Vert })/\overline{\mathcal{Z}(k_1,\nu )}d\nu .\hfill \end{array}$$

In this case, since g₋(ν) is defined through negative frequencies, and $\overline{\mathcal{Z}(\nu )}$ is holomorphic in the lower half-plane and is bounded by unity at infinity, then the integral of $g_{-}/\overline{\mathcal{Z}}$ tends to zero as t > 0. Therefore, ${\varrho }_f(z,t)={\displaystyle \int exp}(i{k}_{1}z-i{k}_1{v}_{\Vert }t)(g_{+}(v_{\Vert })/\overline{\mathcal{Z}}(v_{\Vert }))d{v}_{\Vert }$.

Assuming that 𝒵(ν) can be continued analytically into the strip Im(ν)∈[−|ν_min|;0], and there exists a quantity ν₀ = ν^† − iν_*^† (ν^† ∈ ℝ, ν_*^† ∈ (0, |ν_min|)), we can shift the integration path ∫_ℝ on the left-hand side of the expression for 𝜚_f(z, t) parallel to the real axis down, below the point ν₀: Im(ν) =  − ν_Im, ν_Im ∈ (ν_*^†, |ν_min|). The contribution to the integral from this pole can be obtained by the residue theorem: ${\varrho }_f(z,t)=-2\pi iexp(i{k}_{1}z-i{k}_1{\nu }_{0}t)(g_{+}(\nu )/{\mathcal{Z}}_{\nu }\prime (k_1,\nu ){|}_{\nu ={\nu }_0}$. Since ik₁ν₀t = ik₁ν^†t + ik₁(−iν_*^†)t, the described density wave will be damped with a real damping coefficient β = k₁ν_*^† (β⁻¹ – the wave decay time); that is, in the lower region of the complex plane, Landau damping (Krall & Trivelpiece 1973; Maslov & Fedoryuk 1985) is observed. To determine ν_*^† and ν_Im, we use the expansion of the function Z in the neighborhood of the point ν^†: Z(ν^†)−iν_*^†(dZ/dν)(ν^†) = 0. Thus, if we isolate the real part of the equation (Re(Z)(ν^†) = 0), we determine the condition on the phase velocity ν^†: P.V.∫νf₀(ν)/(ν^† − ν)dν = (2π𝜘/k₁²)⁻¹; if we isolate the condition on the imaginary part, we obtain: πν^†f₀(ν^†) = ν_*^†P.V.∫νF₀(ν)/(ν^† − ν)²dν.

Thus, we obtain a complete description for the density waves of self-gravitating particles moving in one direction, provided that the potential perturbations in the neighborhood of its macro-extrema point (for the equilibrium function F₀(v), coinciding with or being a direct generalization of the Maxwellian) obey the linearized Poisson equation. Van Kampen waves admit a more general form of the ansatz, when normal modes have a more universal form than plane waves (Case 1960). We demonstrate its application to the system of gravitating particles under consideration, which is essential for the two-dimensional geometry of a system with rotation.

Consider the “conjugate” problem to (20) in the following form:

$$\begin{array}{cc}& (\nu -v_{\Vert })\mathfrak{A}(\mathit{k},v_{\Vert };{\omega }^{‡})={\displaystyle \int {\nu }_{\varkappa }}(k,v)\mathfrak{A}(\mathit{k},v;{\omega }^{‡})dv,\hfill \\ \hfill & {\displaystyle \int {\nu }_{\varkappa }}(k,v)\mathfrak{A}(\mathit{k},v;{\omega }^{‡})dv=1,\hfill \\ \hfill & ({\nu }^{‡}-v_{\Vert })\mathfrak{A}(\mathit{k},v_{\Vert };{\omega }^{‡})=1,({\nu }^{‡}-{\omega }^{‡})f(\mathit{k},v_{\Vert };{\omega }^{‡})\mathfrak{A}(\mathit{k},v_{\Vert };{\nu }^{‡})=0,\hfill \end{array}$$(25)

where normal modes are introduced by the relation 𝔄(k, v_∥, t) = 𝔄(k, v_∥; ω)exp(−iωt). If the real eigenvalues ω^‡ are not zeros of the function ν_𝜘(k, v), then the eigenfunctions that correspond to them take the form 𝔄(k, v_∥; ω^‡) = ν^‡(k, ω^‡)δ(ω^‡ − v_∥)+P.V.(1/(ω^‡ − v_∥)); further, we should consider the cases when: 1)ω^‡ are zeros of the function ν_𝜘(k, v_∥), but not ν^‡(k, v_∥); 2)ω^‡ are the zeros of the functions ν_𝜘(k, v_∥) and ν^‡(k, v_∥); 3)ω_j^‡ are the complex zeros of 𝔄(k, v_∥; ω_j^‡). Finally, we obtain

$$\begin{array}{cc}& \mathfrak{A}(\mathit{k},v_{\Vert };{\omega }^{‡})={\displaystyle \sum _j}C(k,j)\mathfrak{A}(\mathit{k},v_{\Vert };{\omega }_j^{‡})+\hfill \\ \hfill & {\displaystyle \int C}(k,j){\omega }^{‡}\mathfrak{A}(\mathit{k},v_{\Vert };{\omega }^{‡})d{\omega }^{‡}.\hfill \end{array}$$

The amplitude of the modes is obtained as the sum over the discrete and continuous spectra of the singularities of the functions ν_𝜘(k, v_∥) and ν^‡(k, v_∥).

Thus, van Kampen waves in the linear approximation for the Poisson equation, with initial conditions that only depend on the particle velocities, can serve as a basis for the quasi-local approximation near the extremum point of the self-consistent potential. In the formulation of the problem of the evolution of cosmological structures, such an approach is applicable for the initial stages of the process of their formation, when the gravitational interaction does not yet have a significant effect on the topological properties of the selected system of particles. It would be interesting to estimate the change in the sizes of proto-structures during the transition to the phase of gravitational interaction dominance from the point of view of an absence of solutions to Equation (11), since this would lead to the proto-structures to a quasi-Jeans-type decay caused by the presence of an additional term – the cosmological term – in the Liouville-Gelfand equation.

4. Aperiodic structures

In addition to van Kampen waves, the Vlasov–Poisson system of equations has wave solutions of a very general type, which can also be associated with cosmological structures. We are talking about one-dimensional Bernstein–Green–Kruskal (BGK) waves (Bernstein et al. 1957; Montgomery 1960; Schwarzmeier et al. 1979). For the simplest one-dimensional case, the Vlasov equation in coordinates (𝔈, x, t) (𝔈 = mv²/2 + mΦ(x) is the energy of a particle in a gravitational field:

$$\begin{array}{cc}\hfill \frac{\partial F(\mathfrak{E},x,t)}{\partial t}& +v(x,\mathfrak{E})\frac{\partial F}{\partial x}+(v(x,\mathfrak{E})/m)(G(x,t)-{\mathrm{\Phi }}^{\prime }(x))\frac{\partial F}{\partial \mathfrak{E}}=0,\hfill \\ \hfill & -\frac{\partial G}{\partial x}=4\pi \gamma N{\displaystyle \int f}(\mathfrak{E},x,t)dv-c^2\mathrm{\Lambda }\hfill \end{array}$$(26)

where the second term on the right-hand side corresponds to the repulsive potential, as before. At equilibrium, f = f₀(𝔈), 𝔈 = −dΦ/dx; if we set F = F₀(𝔈)+f(x, 𝔈, t), G(x, t) =  − Φ′(x)+G₁(x, t), then the linearized Vlasov equation takes the form

$$\begin{array}{c}\hfill (v(x,\mathfrak{E}){)}^{-1}\frac{\partial f}{\partial t}+\frac{\partial f}{\partial x}-\frac{G_1}{m}\frac{d{F}_0}{d\mathfrak{E}}=0\end{array}$$(27)

(the repulsive potential is absent in the equation for perturbations, since its effect is present in the basic macro-potential Φ(x)). We seek a solution to the equation in the form f = ψ(x)exp(−iωt):

$$\begin{array}{cc}& \frac{\partial \psi }{\partial x}-i\omega v^{-1}(x,\mathfrak{E})\psi =G_1/m·\frac{d{F}_0}{d\mathfrak{E}},\hfill \\ \hfill & -i\omega G_1=4\pi \gamma {\displaystyle \int \psi }vdv=4\pi \gamma {\displaystyle {\int }_{{\mathfrak{E}}_0}^{\infty }}\psi (x,\mathfrak{E})d\mathfrak{E},\hfill \end{array}$$(28)

where v⁻¹(x, 𝔈) = (2𝔈 + Φ)^−1/2). If we exclude G₁ from the last two equations, we obtain an equation of the form

$$\begin{array}{c}\hfill \frac{\partial \psi }{\partial x}-i\omega v^{-1}(x,\mathfrak{E})\psi =(4\pi i/m){\omega }^{-1}\frac{d{F}_0}{d\mathfrak{E}}{\displaystyle {\int }_{{\mathfrak{E}}_0}^{\infty }}\psi d\mathfrak{E}.\end{array}$$(29)

If Φ → 0, then the last equation coincides with the eigenvalue equations obtained in the van Kampen method. Therefore, following the previously considered method, we select the “normal” mode with a fixed wave number k = K₁ and the corresponding frequency Ω⁽⁰⁾ related via the dispersion relation (29):

$$\begin{array}{cc}\hfill ik\psi (k;K_1)& -i\omega {\displaystyle {\int }_{\mathbb{R}}}{v}^{-1}(q)\psi (k-q;K_1)dq\hfill \\ \hfill & =i(4\pi F_0^{\prime }\gamma /m)/\omega {\displaystyle {\int }_{{\mathfrak{E}}_0}^{\infty }}\psi (k;K_1)d\mathfrak{E}.\hfill \end{array}$$(30)

This equation can be solved by expanding in powers of the parameter Φ(k)/𝔈₀: Θ_j(k) = ∑_{k = 0, …, ∞}Θ^(j)(k), where Θ(j)(k)∈{v⁻¹(k),ψ(k; K₁),ω}. Putting $v^1(k)=\delta (k)/\sqrt{2\mathfrak{E}}$, in the zeroth approximation we obtain two types of eigenmodes, discrete and continuous ${\psi }^{(0)}(k;K_1)=(K_1-{\omega }^{(0)}/\sqrt{2\mathfrak{E}}{)}^{-1}((4\pi F_0\prime \gamma /m)/{\omega }^{(0)})\delta (k-K_1)$. The criterion for discreteness of the quantities ω⁽⁰⁾ are the conditions (4πF₀′γ/m)⋅((ω⁽⁰⁾)²/(2K₁²)) = 0, or the condition Im(ω⁽⁰⁾)≠0. If (4πF₀′γ/m)⋅((ω⁽⁰⁾)²/(2K₁²)) ≠ 0, the functions ψ⁽⁰⁾(k; K₁) should be considered in the class of distributions, since $(K_1-{\omega }^{(0)}/\sqrt{2\mathfrak{E}}{)}^{-1}=P.V.(1/(K_1-{\omega }^{(0)}/\sqrt{2\mathfrak{E}}))+{({\mu }^{‡})}^{(0)}(K_1,{\omega }^{(0)})\delta (K_1-{\omega }^{(0)}/\sqrt{2\mathfrak{E}})$ (in this case Im(ω⁽⁰⁾)≤0 indicates the asymptotic stability of the complete solution. In a similar way, one can obtain ψ^{(1, 2, …)}(k; K₁).

The main result after constructing the appropriate number of terms in the series for ψ(k; K₁) is the establishment of the density function of the solution of the BGK equations. This expression can be used for comparative calculations of the macro-parameters of cosmological objects (see below).

As can be seen from the form of Equation (29), the initial condition is also taken in the form of a (generalized) Maxwell function, and the methodology of further research makes significant use of this. To what extent is it legitimate in general to use F_MB in the role of the Cauchy conditions for the Vlasov equation for cosmological systems (for the linearized case, f⁽⁰⁾(x, v))? In accordance with the structure of Equation (15), the formal substitution of normal modes (of the form 𝔎₁(v)𝔎₂(x, t), in the simplest case 𝔎₂(z, t) = exp(ikz − iωt)) into this equation at F|_t = 0 = F_MB(x, v) will lead to the appearance of a bilinear dependence on the spatial and temporal variables, which indicates a nonlocal form of interaction of carrier waves, which should be described by an integral relation, which excludes the presence of a local differential dispersion formula. Apparently, the most direct way to study the properties of the linear Vlasov equation for an inhomogeneous field and initial conditions lies through finding the explicit form of the force interaction term (for F₀ → F_MB).

In this case, there are obviously problems when substituting into the equation decomposition solutions of a priori form with independent modulation by coordinates of the extended phase space. Following (Maslov 1994; Maslov & Fedoryuk 1985), we assume that the characteristics of the linear (complete) Vlasov equation coincide with the phase trajectories of the Hamiltonian system dX/dt = V, dV/dt = −dΦ/dX, since one should consider the additional term $\mathcal{T}(\mathrm{\Phi },f)\equiv -{\mathrm{\Phi }}_{\mathit{x}}\prime f_{\mathit{v}}\prime$ on the left-hand side of Equation (16); the spatial changes in the potential of the “main” gravitational field of the system are taken into account. The function Φ satisfies Equation (10) (or (11), if after obtaining the solution we pass from the dependent variable W to Φ). The solution of this dynamical system with the initial conditions X|_t = 0 = x, V|_t = 0 = v is as follows: X(x, v, t), V(x, v, t) (t ∈ R¹). The first integral of the dynamic system is 𝔈 = mv²/2 + Φ(x) (which corresponds to the conservation of energy along the trajectories of the Vlasov equation in the spatially inhomogeneous case, and this is why the term 𝒯(Φ, f) was introduced). For the function f(x, v, t), through the shift along the trajectories from the initial point, we have the IInd type Volterra equation

$$\begin{array}{cc}\hfill f(\mathit{x},\mathit{v},t)& =f^{(0)}(\mathit{X}(\mathit{x},\mathit{v},-t),\mathit{V}(\mathit{x},\mathit{v},-t))\hfill \\ \hfill & +\frac{d{F}_0}{d\mathfrak{E}}{\displaystyle {\int }_0^t}\mathrm{\nabla }\varphi (\mathit{X}(\mathit{x},\mathit{v},\xi -t),\xi )\mathit{V}(\mathit{x},\mathit{v},\xi -t)d\xi ,\hfill \end{array}$$(31)

and, after substituting this expression into the Poisson equation ∇²ϕ = λ^♯exp(W(x)) ⋅ ϕ(x), we have an explicit form for the force term (G → G[Φ]+g[ϕ] when linearized):

$$\begin{array}{cc}\hfill -{\mathfrak{y}}^{-1}\mathrm{\nabla }\varphi (\mathit{x},t)& ={\displaystyle \int \int f^{(0)}}(\mathit{X}(\mathit{x},\mathit{v},-t),\mathit{V}(\mathit{x},\mathit{v},-t))d\mathit{v}d\mathit{x}+\hfill \\ \hfill & +{\displaystyle \int \int {\int }_0^t}\frac{d{F}_0}{d\mathfrak{E}}\mathrm{\nabla }\varphi (\mathit{X}(\mathit{x},\mathit{v},-\xi ),t-\xi )\mathit{V}(\mathit{x},\mathit{v},-\xi )d\xi d\mathit{v}d\mathit{x},\hfill \end{array}$$(32)

where the notation 𝔶 ≡ λ^♯exp( − Φ(x)/T). In accordance with the definition in Formula (11) for the potential value W(x), for the motion in a nonuniform field of a system of gravitating particles, we obtain the influence of two integrand factors at once: dF₀/d𝔈 ⋅ g. This is due to the fact that both 𝔈 and g contain the full Liouville–Gel’fand potential. This significantly complicates the consideration of the question of the uniqueness of the solution, since the values of the potential W(x) in these factors may lie in different regions X_i(x) from Section 2. Apparently, in order to establish the uniqueness of the solution, the behavior of the function v(x)|_{𝔈 = const} should be considered. In addition, the question arises of the physical manifestation of the multi-valuedness of solutions to the Vlasov–Poisson equation in the region X₂(x): since the norms of the solutions W(x) with the same pre-exponential factor differ by finite values, the standard definition of bifurcation of solutions is inapplicable, and smooth solutions of the Vlasov equation that correspond to the minimal norm of the solution must collapse. However, “destruction of the solution” can be expressed in an increase in its norm, for example, due to an increase in the density of particles, which can be a time-dependent process. Consequently, in addition to the wave form of motion, in the simplest case considered in Section 3 using the example of van Kampen waves, there may be processes of local “thickening” of matter over time in a certain region of space (antinodes of a longitudinal wave, in particular), associated with the transition in the region of multivalued solutions of the Liouville-Gel’fand equation to a new norm of its solution.

We point out that the left-hand side of the Vlasov–Poisson equation with an additional term 𝒯(Φ, f) as Φ → const tends can be assumed to be extremely close to the “classical” left-hand side of the linearized Equation (16), however, the right-hand side of the kinetic equation, containing the second term of the right-hand side of the Volterra Equation (31), will retain an unchanged form (𝔈 ≈ v²/2 + Φ(x₀)), and this part only slightly depends on function f. Consequently, we can formally consider the representation of the solution in the form of a normal mode of the above-considered “ansatz” type, divide both parts by (ω − ν) (taking into account the occurrence of the term in the form of a distribution), and repeat all the operations of Section 3. In this regard, van Kampen waves can also be used for the spatially–(weakly)inhomogeneous case. Let us demonstrate this by turning to the one-dimensional case (that corresponds to the previously considered longitudinal waves) for the sake of clarity of the calculations. We integrate both parts (32) over the interval $[0,\stackrel{\sim }{z}]$, rearrange the order of integration, and make a change of the variables η_X = X(x, v, −ξ), η_V = V(x, v, −ξ) in the second term of the right-hand side. Since 𝔈(X, V) = 𝔈(η_X, η_V), dXdV = dη_Xdη_V, this term will take the form of a flow through the surface: ∫∫_σ𝔈(η_X, t − ξ)∂(F₀(η_V²/2 + Φ(η_X))/∂η_V)dη_Xdη_V. The boundary ∂σ is the image of the line η_X on the plane (η_X, η_V) with a shift in time −ξ along the phase trajectories of the dynamic system of the system Ẋ = V, $\dot{V}=-{\mathrm{\Phi }}_X$ (in our case, a small value). We assume that the boundary ∂σ is analytically defined by the relation ${\eta }_V=\beta ({\eta }_X|\xi ,\stackrel{\sim }{z})$ (η_V < β∀(η_X, η_V)∈σ). Then the second term under study will take the form ∫g(η_X, t − ξ)F₀(𝔈[η_X, η_V])dη_X. Therefore, the right-hand side of (32) has the form

$$\begin{array}{cc}\hfill \mathfrak{g}(\stackrel{\sim }{z},t;F_0)& \equiv {\displaystyle {\int }_0^{\stackrel{\sim }{z}}}{g}_0(z,t)dz+{\displaystyle \int {\int }_0^t}g(t-\xi ,{\eta }_X)F_0({\beta }^2({\eta }_X|\xi ,\stackrel{\sim }{z})\hfill \\ \hfill & +\mathrm{\Phi }({\eta }_X))d\xi d{\eta }_X\hfill \end{array}$$(33)

where the tilde sign over the variable z is omitted. If we substitute into the Vlasov equation with this right-hand side (and formal annulment or replacement of the quantity 𝒯(Φ, f) by an approximating term) the normal mode of the van Kampen type 𝔎₁(v), then the left-hand side will take the form (ω − kv)𝔎₁(v)𝔎₂(z, t), and the right-hand side $i\mathfrak{y}{(F_0)}_v\prime \mathfrak{g}(z,t;F_0)\equiv \mathfrak{S}(z,t;v)$. It should be noted that this operation was enabled by the special structure of the Vlasov equation, since the gravitational field strength here is a function closed on itself as solution to the integral equation. Dividing both parts of the resulting equation by (v − ω/k) leads to the need to take into account an additional term, considered as a distribution exiting to the space of generalized functions

$$\begin{array}{cc}\hfill f(k,v;\omega )& =-k^{-1}\mathfrak{S}(z,t;v)·P.V.(1/(\nu -v)){|}_{\nu =\omega /k}\hfill \\ \hfill & +\vartheta (k,\nu )\delta (\nu -v){|}_{\nu =\omega /k},\hfill \end{array}$$

where ϑ(k, ν) is the normalization function (ϑ = 1 + ∫(−k⁻¹𝔖(z, t; v)/(v − ν))dv).

The solution of the initial value problem f(k, v, t) can be represented as an expansion in special solutions f(k, v; ν)exp(−ikνt): f(k, v, t) = ∫𝒰(k, ν)f(k, v; ν)exp(−ikνt)dν; accordingly, the Cauchy condition f⁽⁰⁾(k, v)≡f(k, v, t = 0) = ∫𝒰(k, ν)f(k, v; ν)dν. To determine the coefficients of 𝒰, we obtain a singular integral equation:

$$\begin{array}{cc}& \mathcal{U}(k,v)=-k^{-1}\mathfrak{S}(z,t;v)·P.V.{\displaystyle \int (}\mathcal{U}(k,\nu )/(\nu -v))d\nu +\vartheta (k,v)\mathcal{U}(k,v).\hfill \end{array}$$

Its solution looks like

$$\begin{array}{cc}& \mathcal{U}(k,\nu )=\frac{{\mathcal{G}}_{+}(k,\nu )}{1+2\pi i{\mathcal{H}}_{+}(k,\nu )}-\frac{{\mathcal{G}}_{-}(k,\nu )}{1+2\pi i{\mathcal{H}}_{-}(k,\nu )},\hfill \\ \hfill & {\mathcal{G}}_{+}(k,\nu )-{\mathcal{G}}_{-}(k,\nu )=\mathcal{U}(k,v),\hfill \\ \hfill & {\mathcal{G}}_{+}(k,\nu )+{\mathcal{G}}_{-}(k,\nu )=\frac{1}{\pi }{\displaystyle \int \frac{\mathcal{U}(k,\nu )}{\nu -v}}d\nu ,\hfill \\ \hfill & {\mathcal{H}}_{+}(k,\nu )-{\mathcal{H}}_{-}(k,\nu )=-k^{-1}\mathfrak{S}(z,t;v),\hfill \\ \hfill & {\mathcal{H}}_{+}(k,\nu )+{\mathcal{H}}_{-}(k,\nu )=\frac{1}{\pi i}{\displaystyle \int \frac{-k^{-1}\mathfrak{S}(z,t;v)}{\nu -v}}d\nu .\hfill \end{array}$$

Thus, we have obtained a method for applying van Kampen waves to a formally weakly inhomogeneous system of particles when the gravitational field strength of the complete system changes slowly. Some explanations are required here, which are related to the presence of a cosmological term in the Liouville–Gel’fand equation. The function 𝔖(z, t; v) is defined through relation (33) and contains the factor 𝔶. Recall that in the second term on the right-hand side (33), there is a “full potential” Φ(x), which is a solution to the nonlinear Poisson Equation (10), in which the influence of Λ-repulsion is taken into account due to the cosmological term: ΔΦ(x) = λ^†exp( − Φ(x)/T) − c²Λ. Further, the quantity 𝔶 is defined as 𝔶 ≡ λ^♯exp( − Φ(x)/T), where, in turn, the pre-exponential factor ${\lambda }^{♯}=-\frac{{\lambda }^{†}}{T}exp(c^2\mathrm{\Lambda }{\mathit{x}}^2/T)$, i.e., it also depends significantly on the cosmological term. Thus, the influence of the Λ-term on the dynamics of particles in the system under consideration is critical, and is the important factor that requires modification of standard approaches of van Kampen waves and Landau damping, as a consequence of the expansion of Landau modes in van Kampen waves.

5. Conclusions

The need for a more refined understanding of the possible genuine differences in the features of the early and late Universe is especially sharpened by the Hubble tension, the DESI BAO data, and other challenges. Regarding the global scale, the evolution of primordial density perturbations within various dark sector models is considered to describe the cosmic web, the voids, and large-scale filaments. On the local scale, the role of self-consistent gravitational interaction has become crucial and needs proper techniques to deal with and to reveal the intrinsic features of the filaments on that scale.

In this paper we considered the linearized Vlasov–Poisson equation approach, and applied the profound technique developed to analyse the wave processes in plasma physics. Namely, we showed that the van Kampen waves, associated to Landau damping and phase mixing, mark the appearance of aperiodic solutions to the linearized Vlasov-Poisson equation. Aperiodicity then arises as an intrinsic property of the resulting filamentary structures. Of principal importance is the fact that the aperiodic structures arise when we take into consideration the repulsion of the cosmological term. The cosmological term in the local Universe description appears in view of the theorem on the identity of the gravity of the sphere and point mass within the McCrea-Milne model and weak-field General Relativity. As mentioned, that approach already made it possible to explain naturally the Hubble tension, to describe the dynamics of clusters of galaxies, and to explain the local flows.

The appearance of aperiodicity as an intrinsic feature of the local filaments could become the subject of dedicated analyses of the observational surveys, and thus act as a probe for the role of the cosmological constant in the local Universe. Aperiodicity then has to damp upon the increase in the scale of the filaments.

Acknowledgments

We are thankful to the referee for valuable comments.

References