Nikolay N. Fimin
Keldysh Institute of Applied Mathematics of Russian Academy of Sciences, Moscow, Russian Federation.
DOI: 10.4236/jamp.2026.145091   PDF    HTML   XML   24 Downloads   127 Views   Citations

Abstract

This paper is devoted to analysis of the evolution of cosmological structures using the invariant properties of the Vlasov-Poisson system of equations. It is shown that the energy substitution in the Vlasov equation leads to a class of undamped and damped wave oscillations of the van Kampen, Landau and Bernstein type, internally related to each other. These waves in space have very universal properties, in many ways identical to those considered in plasma theory. Large-scale astrophysical systems can be considered as consequences of the implementation of density waves, and since their density is anisotropic, this can affect the measurement of various astronomical parameters.

Keywords

Van Kampen Normal Modes, Dispersion Relation, Vlasov-Poisson Equation, Liouville-Gelfand Equation

Share and Cite:

1. Introduction

The emergence and evolution of large-scale low-dimensional structures (such as the long-known void walls and filaments in Laniakea-type superclusters, as well as recently discovered megascale arc objects) are the subject of close study not only from the point of view of observers recording the time spectrum of their states (in the Earth’s reference frame), but also represents an extremely effective testing ground for testing the modeling of various variants of inhomogeneous modifications of the Friedmann concept of the expansion of the Universe. It is obvious that at present it is impossible to state with complete certainty that we know for sure all the mechanisms of formation and realization of a high degree of coherence of the majority of large-scale structures. In addition to the approach that studies the formation of caustic features in the macromotions of matter at the early stages after de Sitter inflation, mechanisms characteristic of later times are currently being studied, such as, for example, self-assembly due to fluctuations in a preferred direction in a system of gravitating particles, or the formation of a large structure as a local (possibly multiply connected) topological object possessing the property of a “state of relative equilibrium” with extrema of certain (entropy, free energy) thermodynamic potentials; both approaches turn out to be internally deeply connected, although separated by the scale of scales.

Plasma theory has a well-developed and effective mathematical apparatus for analyzing wave motions of various types, suitable for adaptation for gravitational systems. Using it together with the methods of the theory of kinetic equations with a self-consistent field allows us to identify not entirely obvious properties of the solutions of these equations (such as collisionless attenuation), which can help us discover new physical phenomena or explain the nature of obscure observable phenomena.

The above-mentioned formalism has been successfully applied in astrophysics before, although in a rather limited set of problems. Here, among others, we should mention the works of D. Lynden-Bell [1] [2] (they considered the method of applying Landau damping for small perturbations of equilibrium in spherical star clusters and drew an analogy with the Bernstein-Greene-Kruskal waves in plasma theory), P. Vandervoort [3] (for the problem of stationary oscillations of galaxies, the van Kampen wave method in “action-angle” coordinates was proposed and implemented, and the possibility of applying the theory of BGK waves in the model approximation to the study of the properties of the mentioned problem was investigated), V.L. and E.V. Polyachenko [4] [5] (stability of many astrophysical problems in various geometries was studied; it turns out that in the unstable regime, the Landau-damped waves can be represented as a superposition of van Kampen modes plus a discrete damped mode in dynamically stable spherical stellar systems), W.C. Saslaw [6] (the main approaches to modeling clusters of astrophysical objects of various scales were analyzed, including the use of the kinetic approach taking into account collisionless damping of oscillations of limited amplitude), P.L. Palmer [7] (a theory of stability of star clusters and galaxies was constructed based on the theory of eigenfunctions of the perturbed part of the gravitational potential operator, which is equivalent to taking into account Landau damping); specially, it is necessary to highlight the works [8] [9], where the authors directly point to the possibility of using van Kampen wave methods for large-scale movements of clusters and galaxies.

In this paper, an attempt is made to describe large-scale structures using non-dissipation solutions of the Vlasov–Poisson equations of the van Kampen wave type. The periodicity of the waves should be violated when taking into account the repulsive forces due to the inclusion of a cosmological term in the considerations, since in this case our system is locally close to weakly inhomogeneous (for a long-range order—significantly inhomogeneous), which is associated with the inclusion in the analysis of the behavior of a many-particle megasystem of the influence of the cosmological term, which is included in the modified Poisson equation as a source of antigravity [10]. The possibility of introducing Bernstein-Greene-Kruskal waves as structural units of cosmological systems as an alternative to van Kampen waves is considered. For substantially inhomogeneous systems, the possibility of a smooth transition from the integral accounting of the field of gravitational disturbances to the normal mode method is substantiated.

The structure of the article is as follows: in paragraph 2, the system of Vlasov–Poisson equations for gravitationally interacting particles (with taking into account the action of the cosmological term), and also makes a basic assumption (“energy substitution”) about the form of the potential, corresponding to the basic unperturbed distribution function, which is a solution to the Vlasov equation; paragraph 3 establishes the possibility liearization of the Vlasov–Poisson system in a self-gravitating system of particles and the construction of an explicit form of normal modes is carried out for a kinetic equation taking the form of an integral Fredholm equation of the third kind; paragraph 4 discusses the methodology the second type of solution (Bernstein) of the linearized Vlasov equation, obtained by energy substitution, as well as the possibility of generalizing the formalism of van Kampen normal modes to the weakly inhomogeneous case of particle density was demonstrated; paragraph 5 is a conclusion that describes the main results of the work.

2. Linearization of Vlasov-Poisson Equations with Cosmological Term

We will consider the set of $N$ cosmological objects (“particles” with masses $m_{i=1,\cdots ,N}=m\equiv 1$ ), interacting with each other gravitationally. The system of Vlasov–Poisson equations for describing its dynamics may be represented as

$\frac{\partial F\left(x,v,t\right)}{\partial t}+{\nabla }_x\left(vF\right)+\widehat{\mathcal{Q}}\left(F;F\right)=0,$ (1)

${\text{Δ}}_x\Upsilon \left[F\left(x\right)\right]=4\pi AN\gamma {\displaystyle \int F\left(x,v,t\right)\text{d}v}-c^2\Lambda ,$ (2)

where $F\left(x,v,t\right)$ is the distribution function of gravitationally interacting particles, $A$ is a normalization factor for particle density, $\gamma$ is the gravitational constant. The system of particles is situated in large domain of configurational space $\Omega \subset {ℝ}_x^3$ ($\text{diam}\Omega \equiv R_{\Omega }<\infty$ ). The nonlinear Poisson Equation (2) takes the form of an inhomogeneous Liouville-Gelfand equation [11] with local temperature [12].

The third term on the right hand side of the kinetic Equation (1) may be represented as

$\widehat{\mathcal{Q}}\left(F;F\right)=-{\nabla }_x\Upsilon \left[F\left(x\right)\right]\cdot \frac{\partial F}{\partial v},$ (3)

$\Upsilon \left[F\left(x,t\right)\right]=4\pi AN\gamma {\displaystyle \iint \mathfrak{U}\left(x-x^{\prime }\right)F\left(x^{\prime },v^{\prime },t\right)\text{d}{x}^{\prime }\text{d}{v}^{\prime }}+\frac{\Lambda c^2}{6}{\left|x\right|}^2+\widehat{\mathfrak{O}}\left(x,x^{\prime }\right),$(4)

where: $\mathfrak{U}\left(x-x^{\prime }\right)=-{\left|x-x^{\prime }\right|}^{-1}$ (Newtonian interaction kernel), $\widehat{\mathfrak{O}}\left(x,x^{\prime }\right)$ is an operator term that takes into account the influence of the boundary conditions (we will take into account the influence of this term by setting the appropriate boundary conditions). Classical Newtonian potential ${\Upsilon }_N\left(r\right)=-\gamma M/r$ increases monotonically on the interval $r\in \left(0,+\infty \right)$ (${\Upsilon }_N\in \left(-\infty ,0\right)$ ), while the generalized (including a cosmological term) Newton interparticle gravity potential ${\Upsilon }_{GN}\left(r\right)\equiv -\gamma M/r-c^2\Lambda r^2/6$ , increases on the interval $r\in \left(0;r_c\right]$ and decreases on the interval $r\in \left(r_c;\infty \right)$ , where $r_c={\left(3\gamma M/\left(\Lambda c^2\right)\right)}^{1/3}$ .

We’ll consider the nonstationary case of dynamics: $F=F\left(x,v,t\right)$ . For a non-stationary system of equations for the evolution of a cosmological system of particles in the self-consistent approximation (1)-(2), the main role is played by the formulation of the initial problem for the Vlasov equation; in this case, the complete problem for the kinetic equation with a self-consistent gravitational field becomes mixed. In this case, direct derivation of solutions and their study by analytical methods are difficult. In the present work, we restrict ourselves to studying the properties of solutions of the linearized version of the Vlasov–Poisson system for gravity (including both gravity and antigravity for particles cosmological system). We will rely, in particular, on the results of the works [13] [14], in which the validity of using the “energy substitution” in the Vlasov-Poisson equation for systems of particles with a periodic density distribution was established.

The linearization of the Vlasov equation is quite non-trivial, since its result depends significantly on the type of gravitational field (and this type, due to the self-consistency of the problem, 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\left(v;T\right)$ or, in a more general case, $F=F_0\left(v\right)$ , $F_0\in C^1\cap L^2\left({\Omega }_v\right)$ , ${\Omega }_v\subset {ℝ}_v^3$ ).

For a certain class of (much broader) 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}~\text{exp}\left(-\mathfrak{E}\left(x,v,t\right)/T\right)$ , including the (2-particle) potential of this introduced force: $\mathfrak{E}=m{v}^2/2+\Upsilon \left(x,t_0\right)$ (for a fixed instant of the current time); the dual solution of the Poisson equation $\Upsilon \left(x,t_0\right)$ 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).

Equation (2) for gravitational potential can be written as

$\text{Δ}\Upsilon \left(x\right)={\varsigma }^{†}\exp \left(-\Upsilon \left(x\right)/T\right)-c^2\Lambda ,{\varsigma }^{†}=4\pi \gamma N{A}_T,$ (5)

$A_T\equiv {\rho }_0{\displaystyle \int \exp \left(-m{v}^2/\left(2T\right)\right)v^2\text{d}v},{\rho }_0={\left(\frac{m}{2\pi T}\right)}^{3/2}.$

The last equation can be rewritten in the form

$\Delta \mathfrak{M}\left(x\right)={\varsigma }^{♯}\exp \left(\mathfrak{M}\left(x\right)\right),\mathfrak{M}\left(x\right)\equiv -\frac{\Upsilon \left(x\right)}{T}-\frac{c^2\Lambda x^2}{6T},$ (6)

${\varsigma }^{♯}=-\frac{{\varsigma }^{†}}{T}\exp \left(c^2\Lambda x^2/T\right).$

Solutions of the equation $\Delta \mathfrak{M}\left(x\right)=-\zeta \exp \left(\mathfrak{M}\left(x\right)\right)$ ($\zeta \in {ℝ}_{+}^1$ ) in the 3-dimensional case are radially symmetric ($\mathfrak{M}=\mathfrak{M}\left(\left|x\right|\right)$ in the accordance with the Gi-Nidas-Nirenberg theorem [11]) and are unstable with respect to the pre-exponential parameter: their existence and number depend on the value of the parameter $\zeta$ . According to [15], the solution of the standard Dirichlet problem for it has a structure that can be described as follows. Let ${\zeta }_{crit}=2$ (if the boundary value problem is considered on the reduced interval $\left|x\right|\equiv r\in \left[0;1\right]$ ); then there exists ${\zeta }_{FK}>{\zeta }_{crit}$ such that: 1) for $\zeta ={\zeta }_{FK}$ , there is a unique solution (${\mathfrak{M}}_{FK}$ ); 2) for $\zeta >{\zeta }_{FK}$ , there are no solutions; 3) for $\zeta ={\zeta }_{crit}$ , there is a countable infinity of solutions (${\mathfrak{M}}_{crit}^{\left(n\right)}$ , $n\in \mathfrak{V}$ , $card\left(\mathfrak{V}\right)={\aleph }_0$ ); 4) for $\zeta \in \left(0,{\zeta }_{FK}\right)\backslash \left\{{\zeta }_{crit}\right\}$ , there is a finite number of solutions (${\mathfrak{M}}_K^{\left(k\right)}$ , $k\in \left\{1,2,\cdots ,K\right\}$ , $K\ge 1$ ). Since it is possible to uniquely (for fixed parameters $T,N$ ) compare the values of the function $-{\varsigma }^{♯}\left(\left|x\right|\right)$ with the values of the parameter $\zeta$ , it can be stated that with an increase in the modulus of the radius vector $\left|x\right|$ , three regions of solutions to the Equation (6) arise: the region of uniqueness of solutions $X_1\left(x\right)=\left\{\left|x\right|<X^{\left(I\right)}\right\}\cup \left\{X^{\left(III\right)}\right\}$ , the region of multivalued solutions $X_2\left(x\right)=\left\{X^{\left(I\right)}<\left|x\right|<X^{\left(II\right)}\right\}$ , the region of absence of solutions $X_3\left(x\right)=\left\{\left|x\right|>X^{\left(III\right)}\right\}$ .

Let us consider the linearization of the Equation (5) in the neighborhood of the solution $W\left(x\right)$ analytic solution with which $\mathfrak{M}$ ) can be associated by representing the solution (5) as $\mathfrak{M}+\mathfrak{w}\left(x,t\right)$ ($\Vert \mathfrak{w}\Vert \ll \Vert \mathfrak{M}\Vert$ by the chosen norm, and, correspondingly, $\text{exp}\left(\mathfrak{w}\right)\approx 1+\mathfrak{w}$ ). We obtain the linear Poisson equation

$\Delta \mathfrak{w}\left(x\right)={\varsigma }^{♯}\exp \left(\mathfrak{M}\left(x\right)\right)\cdot \mathfrak{w}\left(x\right).$ (7)

Obviously, in the neighborhood of m. $x_0$ the last equation is simplified, since the gravitational field of a point with equivalent total mass and cosmological repulsion allow us to set $\Upsilon \left(x_0\right)\equiv -T\mathfrak{M}\left(x_0\right)={\Upsilon }^{\left(0\right)}\left(=\text{const}\right)$ ; thus, the equation for the potential perturbation in the above neighborhood $O\left(x_0\right)$ takes the form ($K^{\left(0\right)}=\text{exp}\left(\mathfrak{M}\left(x_0\right)\right)$ ):

$\text{Δ}\mathfrak{w}={\varsigma }^{♯}\left(x_0\right)K^{\left(0\right)}\cdot \mathfrak{w}\left(x\right).$ (8)

The linearization of the Vlasov equation itself is performed (in the simplest case under consideration) in the neighborhood of the equilibrium function $F_M\left(v\right)$ or, in a more general case, $F_0\left(v\right)$ with several maxima, which is realized, for example, in the case of codirectional particle beams (for a more general linearization we have:

$F \to F_{MB}\left(x,v\right)+\tilde{f}\left(x,v,t\right),$ (9)

where the perturbation $\tilde{f}\left(x,v,t\right)$ is related to the Poisson Equation (7) with an exponential dependence of the parameter on the spatial variable). Eliminating quadratic by $\tilde{f}$ terms gives us

$\frac{\partial \tilde{f}}{\partial t}+v\cdot {\nabla }_x\tilde{f}-{\nabla }_v{F}_0\cdot {\nabla }_x{\Upsilon }^{♭}\left[\tilde{f}\right]\left(x,t\right)=0,-T\left(\mathfrak{M}+\mathfrak{w}\right)=\Upsilon +{\Upsilon }^{♭}.$(10)

Note that further construction of solutions to the system of Vlasov-Poisson equations will be carried out near the equilibrium $F_M\left(v\right)$ (the formalism proposed by the author does not reveal any significant difference with “multi-peak” distributions that depend only on speed; Apparently, nontrivial conclusions for multibeam systems in a self-gravitating medium will be associated with local the emergence of Landau damping, but this issue is beyond the scope of this work). Weakly inhomogeneous distribution $F_{MB}\left(x,v\right)$ will be discussed in paragraph 4.

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

3. Linearized Vlasov Equation and Van Kampen Modes and Wave Motion for a Self-Consistent Gravitational Potential

Let us consider the invariant properties (independent of solutions) of the linearized Vlasov-Poisson system of Equations (8)-(9), first for the case of the gravitational field strength corresponding to a local neighborhood of the extremum of the self-consistent potential, taking into account the action of the cosmological term:

$\frac{\partial \tilde{f}\left(x,v,t\right)}{\partial t}+v\nabla \tilde{f}\left(x,v,t\right)={\nabla }_x{\Upsilon }^{♭}\left(x,t\right)\cdot {\nabla }_v{F}_0\left(v\right),$ (11)

${\nabla }_x{\Upsilon }^{♭}\left(x,t\right)={\varsigma }_0^{♯}{K}^{\left(0\right)}{\displaystyle {\int }_{{\Omega }_{x^{\prime }}}{\displaystyle {\int }_{{\Omega }_v}{\nabla }_x\frac{\tilde{f}\left(x^{\prime },v,t\right)}{\left|x-x^{\prime }\right|}\text{d}v\text{d}{x}^{\prime }}}.$

We represent $\tilde{f}\left(x,v,t\right)$ via the van Kampen ansatz or “normal modes” [16]-[18]: $\mathfrak{K}\left(v\right)\text{exp}\left(ikx-i\omega t\right)$ (plane waves are eigenfunctions of the Laplacian from the left-hand side of the Poisson Equation (8)). We will be interested in solutions–perturbations of the system of Equations (8)-(9) in the form of longitudinal waves, therefore we choose in the velocity space axes parallel ($z$ ) and perpendicular ($x,y$ ) to the wave vector $k$ ; then the longitudinal component of the velocity is $v_{\parallel }=v_z=e_k\cdot v$ (where $e_k=k/\left|k\right|$ ), the transverse component, respectively: $u=v-e_k{v}_{\parallel }$ . In this case, we can introduce distribution functions that depend only on one component of the velocity:

$f\left(k,v_{\parallel },t\right)={\displaystyle \int \tilde{f}\left(k,v,t\right)\delta \left(v_{\parallel }-k\cdot v/k\right)\text{d}v}={\displaystyle \int \tilde{f}\left(k,v,t\right)\text{d}u}$.

Let us rewrite the Equation (9) for these modes, freeing ourselves from the transverse velocity components (and discarding the tilde sign over $f$ ):

$\left(\omega -k{v}_{\parallel }\right){\displaystyle \int \mathfrak{K}\left(v\right)\text{d}u}=\frac{4\pi }{k^2}k{\varsigma }_0^{♯}{K}_{\Lambda }^{\left(0\right)}{\displaystyle \int \frac{\partial F_0}{\partial v_{\parallel }}\text{d}u}{\displaystyle \int \mathfrak{K}\left(v^{\prime }\right)\text{d}{v}^{\prime }},$ (12)

$-\widehat{{\Upsilon }^{♭}}\left(k,t\right)={\varsigma }_0^{♯}{K}_{\Lambda }^{\left(0\right)}{\displaystyle \iint \frac{\mathfrak{K}\left(v^{\prime }\right)\left(x-x^{\prime }\right)}{{\left|x-x^{\prime }\right|}^3}\text{exp}\left(ik{x}^{\prime }-i\omega t\right)\text{d}{v}^{\prime }\text{d}{x}^{\prime }},$

${\displaystyle \int \frac{x}{{\left|x\right|}^3}\text{exp}\left(ikx\right)\text{d}x}=4\pi i\frac{k}{k^2}.$

We divide both sides of the last equation by $\left(\omega -k{v}_{\parallel }\right)$ and integrate with respect to the variable $v_{\parallel }$ . The integral ${\displaystyle \int \mathfrak{K}\left(v^{\prime }\right)\text{d}{v}^{\prime }}$ (an unimportant constant) cancels out, and we obtain a dispersion relation that is invariant with respect to the form of the solution of the kinetic equation:

$1-\left(\varkappa /k\right){\displaystyle \int \frac{\text{d}{f}_0}{\text{d}{v}_{\parallel }}\frac{\text{d}{v}_{\parallel }}{\omega -k{v}_{\parallel }}\text{d}{v}_{\parallel }}=0,\varkappa =4\pi {\varsigma }_0^{♯}{K}_{\Lambda }^{\left(0\right)}.$ (13)

If we do not consider the longitudinal velocity as distinguished, then the general form of the dispersion law has the form: $D\left(k,\omega \right)\equiv 1-\varkappa \left(k/k^2\right){\displaystyle {\int }_L{\left(F_0\right)}^{\prime }{}_v{\left(\omega -kv\right)}^{-1}\text{d}v}=0$ (normal modes will correspond to the case $\mathrm{Re}\left(\omega \left(k\right)\right)\gg \mathrm{Im}\left(\omega \left(k\right)\right)$ ).

We will be interested in the possibility of obtaining a solution of the Vlasov–Poisson equations that is stable in time and associated with the simplest cosmological structures (low dimensionality). It can be obtained using normal modes in the form

$f\left(z,v_{\parallel },t\right)={\displaystyle \iint {\mathfrak{O}\left(k,\nu \right)\mathfrak{N}\left(k,\nu ;v_{\parallel }\right)\text{exp}\left(ikz-ik\nu t\right)|}_{\nu =\omega /k}\text{d}k\text{d}\nu },$ (14)

where $\mathfrak{O}\left(k,\nu \right)$ is some (admissible) function (which corresponds to certain Cauchy data for the kinetic equation for the perturbation $f$ ). If the initial condition is represented as $f\left(z,v_{\parallel },t=0\right)={\displaystyle \int g\left(k,v_{\parallel }\right)\text{exp}\left(ikz\right)\text{d}k}$ , then, obviously, equation (14) is reduced to the form ${\displaystyle \int \mathfrak{O}\left(k,\nu \right)\mathfrak{N}\left(k,\nu ;v_{\parallel }\right)\text{d}\nu }=g\left(k,v_{\parallel }\right)$ , and the variable $k$ here acquires the meaning of a parameter.

For what follows, we return to Equation (12) and consider a non-obvious consequence of taking the integral of $\mathfrak{K}$ over the transverse velocities and dividing both parts by $\left(\omega -k{v}_{\parallel }\right)$ . 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 $\left(x-y\right){\mu }_1\left(x\right)={\mu }_2\left(x\right)$ (defined on the interval $\left[x_1;x_2\right]$ of the real axis) and the point $y\in \left(x_1;x_2\right)$ , the solution must be interpreted as a distribution. This distribution can be written in the following form: ${\mu }_1\left(x|y\right)={\mu }_2\left(x\right)P.V.\frac{1}{x-y}+{\mu }^{♮}\left(y\right)\delta \left(x-y\right)$ (where the Cauchy principal value in the form of a distribution is defined by the relation $\left(P.V.\frac{1}{x},\mu \right)={\text{lim}}_{ϵ\to 0}{\displaystyle {\int }_{\left|x\right|\ge ϵ}\left(\mu \left(x\right)/x\right)\text{d}x}$ ), and ${\mu }^{♮}\left(y\right)$ is “the strength of the concentration” of the Dirac function at the point $x=y$ determined from additional conditions imposed on the generalized function ${\mu }_1\left(x|y\right)$ .

Thus, Equation (12), rewritten as

$\left(\nu -v_{\parallel }\right)\mathfrak{N}\left(v_{\parallel }\right)={\nu }_{\varkappa }\mathfrak{F}\left(v_{\parallel }\right){\displaystyle \int \mathfrak{K}\left({v^{\prime }}_{\parallel }\right)\text{d}{v^{\prime }}_{\parallel }},{\displaystyle \int \mathfrak{K}\left({v^{\prime }}_{\parallel }\right)\text{d}{v^{\prime }}_{\parallel }}=1,$ (15)

$\mathfrak{N}\left(v_{\parallel }\right)\equiv {\displaystyle \int \mathfrak{K}\left(v\right)\text{d}u},\mathfrak{F}\left(v_{\parallel }\right)\equiv {\displaystyle \int \frac{\partial F_0\left(v\right)}{\partial v_{\parallel }}\text{d}u},\nu =\frac{\omega }{k},{\nu }_{\varkappa }=\frac{\varkappa }{k^2},$

after dividing both sides of Equation (15) by $\left(\nu -v_{\parallel }\right)$ should be written in the sense of distributions:

$\mathfrak{N}\left(v_{\parallel }\right)={\nu }_{\varkappa }\cdot P.V.\frac{\mathfrak{F}\left(v_{\parallel }\right)}{\nu -v_{\parallel }}+{\mu }^{♮}\delta \left(\nu -v_{\parallel }\right),\left(\nu -v_{\parallel }\right)\delta \left(\nu -v_{\parallel }\right)=0,$ (16)

in this case, from the normalization condition in (15), the intensity value ${\mu }^{♮}$ is determined by the condition of its agreement with the formula (16): ${\mu }^{♮}=1-{\nu }_{\varkappa }P.V.{\displaystyle \int \left(\mathfrak{F}\left(v_{\parallel }\right)/\left(\nu -v_{\parallel }\right)\right)\text{d}{v}_{\parallel }}$ .

Let’s substitute into the equation ${\displaystyle \int \mathfrak{O}\left(k,\nu \right)\mathfrak{N}\left(k,\nu ;v_{\parallel }\right)\text{d}\nu }=g\left(k,v_{\parallel }\right)$ the value $\mathfrak{N}\left(v_{\parallel }\right)$ from (16):

$\mathfrak{O}\left(k,v_{\parallel }\right)\left(1-\pi {\nu }_{\varkappa }^2\widehat{ℌ}ℱ\left(v_{\parallel }\right)\right)-\widehat{ℌ}\mathfrak{O}\left(k,v_{\parallel }\right)\pi {\nu }_{\varkappa }^2ℱ\left(v_{\parallel }\right)=g\left(k,v_{\parallel }\right)$ (17)

($k$ is still a parameter). Here $\widehat{ℌ}\left(\Psi \left(x\right)\right)=\left(1/\pi \right)P.V.{\displaystyle \int \left(\Psi \left(x^{\prime }\right)/\left(x-x^{\prime }\right)\right)\text{d}{x}^{\prime }}$ is the Hilbert transform, which is related to the Fourier transform of the function $\Psi \left(x\right)={\Psi }_{+}\left(x\right)+{\Psi }_{-}\left(x\right)$ : ${\Psi }_{+}\left(x\right)-{\Psi }_{-}\left(x\right)=i\widehat{ℌ}\left(\Psi \left(x\right)\right)$ , where $Y_{+}\left(x\right)\equiv {\displaystyle {\int }_0^{\infty }Y\left(q\right)\exp \left(iqx\right)\text{d}q}$ , $Y_{-}\left(x\right)\equiv {\displaystyle {\int }_{-\infty }^{0}Y\left(q\right)\text{exp}\left(iqx\right)\text{d}q}$ (symbols $Y,Y_{\pm }$ are used to denote the functions $\Psi ,ℱ,\mathfrak{O},g$ and their decompositions).

The last equation can be rewritten as follows:

$\begin{array}{l}\left(1+2\pi i{\nu }_{\varkappa }^2{ℱ}_{+}\left(v_{\parallel }\right)\right){\mathfrak{O}}_{+}\left(k,v_{\parallel }\right)+\left(1-2\pi i{\nu }_{\varkappa }^2{ℱ}_{-}\left(v_{\parallel }\right)\right){\mathfrak{O}}_{-}\left(k,v_{\parallel }\right)\\ =g_{+}\left(k,v_{\parallel }\right)+g_{-}\left(k,v_{\parallel }\right)\equiv g\left(k,v_{\parallel }\right).\end{array}$ (18)

The terms on the left-hand side are analytic and have no singularities in the upper ($Im\left(\eta \right)>0$ ) and lower ($Im\left(\eta \right)<0$ ) parts of the complex $\left({\eta }_{Re},{\eta }_{Im}\right)$ –plane ($ℝ\ni v_{\parallel }\to \eta \in ℂ$ ), respectively, and also asymptotically tend to zero in their half-plane. The decomposition of $g\left(\eta \right)$ into two functions with such properties is unique, and therefore $\left(1\pm 2\pi i{\nu }_{\varkappa }^2{ℱ}_{+}\left(v_{\parallel }\right)\right){\mathfrak{O}}_{\pm }\left(k,v_{\parallel }\right)=g_{\pm }$ . Therefore, if there is a solution (17), then it must coincide with $\mathfrak{O}={\mathfrak{O}}_{+}+{\mathfrak{O}}_{-}$ , ${\mathfrak{O}}_{\pm }=g_{\pm }/\left(1+2\pi i{\nu }_{\varkappa }^2{ℱ}_{\pm }\right)$ (the condition for this is $\mathfrak{N}\left(v_{\parallel }\right)\ne 0$ , which is true, in particular, for the Maxwellian distribution). Considering on the half-plane $Im\left(\eta \right)>0$ a holomorphic and asymptotically close to unity function $\mathcal{Z}\left(\eta \right)=1+2\pi i{\nu }_{\varkappa }^2{ℱ}_{+}\left(v_{\parallel }\right)$ , we can extend it to the half-plane $Im\left(\eta \right)<0$ : $\mathcal{Z}\left(\eta \right)=1+4{\pi }^{2}i{\nu }_{\varkappa }^{2}z\mathfrak{N}\left(\eta \right)+2\pi {\nu }_{\varkappa }^2{\displaystyle \int {\eta }^{\prime }\mathfrak{N}\left({\eta }^{\prime }\right)/\left({\eta }^{\prime }-\eta \right)\text{d}{\eta }^{\prime }}$ . Now we can write out the final form of the solution of the initial value problem with the general solution (14):

$\begin{array}{l}f\left(z,v_{\parallel },t\right)={\left(2\pi \right)}^{-1}\int \int \int \mathfrak{N}\left(k,\nu ;v_{\parallel }\right)\text{exp}\left(ik\left(z-z^{\prime }\right)-ik\nu t\right)\\ \times \left(f_{+}\left(z^{\prime },\nu ,t=0\right)/\mathcal{Z}\left(k,\nu \right)+f_{-}\left(z^{\prime },\nu ,t=0\right)/\overline{\mathcal{Z}\left(k,\nu \right)}\right)\text{d}k\text{d}{z}^{\prime }\text{d}\nu ,\\ f_{+}\left(z,\nu ,0\right)+f_{-}\left(z,\nu ,0\right)=f\left(z,\nu ,0\right).\end{array}$(19)

For the initial function of the form $f\left(z,v_{\parallel },t=0\right)={\displaystyle \int g\left(v_{\parallel }\right)\text{exp}\left(ikz\right)\delta \left(k-k_1\right)\text{d}k}$ ($\varsigma =2\pi /k_1=\text{const}$ ) the density of particles in the disturbance wave is:

${\varrho }_f\left(z,t\right)=\text{exp}\left(i{k}_{1}z\right){\displaystyle {\int }_{ℝ}\text{exp}\left(-i{k}_1\nu t\right)\left(g_{+}\left(v_{\parallel }\right)/\mathcal{Z}\left(k_1,\nu \right)+g_{-}\left(v_{\parallel }\right)/\overline{\mathcal{Z}\left(k_1,\nu \right)}\right)\text{d}\nu }.$

In this case, since $g_{-}\left(\nu \right)$ is defined through negative frequencies, and $\overline{\mathcal{Z}\left(\nu \right)}$ 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\left(z,t\right)={\displaystyle \int \text{exp}\left(i{k}_{1}z-i{k}_1{v}_{\parallel }t\right)\left(g_{+}\left(v_{\parallel }\right)/\overline{\mathcal{Z}}\left(v_{\parallel }\right)\right)\text{d}{v}_{\parallel }}$.

Assuming that $\mathcal{Z}\left(\nu \right)$ can be continued analytically into the strip $Im\left(\nu \right)\in \left[-\left|{\nu }_{\min }\right|;0\right]$ , and there exists a quantity ${\nu }_0={\nu }^{†}-i{\nu }_{\text{*}}^{†}$ (${\nu }^{†}\in ℝ$ ,${\nu }_{*}^{†}\in \left(0,\left|{\nu }_{\min }\right|\right)$ ), we can shift the integration path ${{\displaystyle \int }}_{ℝ}$ on the left-hand side of the expression for ${\varrho }_f\left(z,t\right)$ parallel to the real axis down, below the point ${\nu }_0$ : $Im\left(\nu \right)=-{\nu }_{Im}$ , ${\nu }_{Im}\in \left({\nu }_{\text{*}}^{†},\left|{\nu }_{\min }\right|\right)$ . The contribution to the integral from this pole can be obtained by the residue theorem: ${\varrho }_f\left(z,t\right)=-2\pi i\exp \left(i{k}_{1}z-i{k}_1{\nu }_{0}t\right){\left(g_{+}\left(\nu \right)/{{\mathcal{Z}}^{\prime }}_{\nu }\left(k_1,\nu \right)\right)|}_{\nu ={\nu }_0}$ . Since $i{k}_1{\nu }_{0}t=i{k}_1{\nu }^{†}t+i{k}_1\left(-i{\nu }_{\text{*}}^{†}\right)t$ , the described density wave will be damped with a real damping coefficient $\beta =k_1{\nu }_{\text{*}}^{†}$ (${\beta }^{-1}$ —the wave decay time), that is, in the lower region of the complex plane, Landau damping [19] [20] is observed. To determine ${\nu }_{\text{*}}^{†}$ and ${\nu }_{Im}$ , we use the expansion of the function $Z$ in the neighborhood of the point ${\nu }^{†}$ : $Z\left({\nu }^{†}\right)-i{\nu }_{\text{*}}^{†}\left(\text{d}Z/\text{d}\nu \right)\left({\nu }^{†}\right)=0$ . Thus, isolating the real part of the equation ($Re\left(Z\right)\left({\nu }^{†}\right)=0$ ), we determine the condition on the phase velocity ${\nu }^{†}$ : $P.V.{\displaystyle \int \nu f_0\left(\nu \right)/\left({\nu }^{†}-\nu \right)\text{d}\nu }={\left(2\pi \varkappa /k_1^2\right)}^{-1}$ ; isolating the condition on the imaginary part, we obtain: $\pi {\nu }^{†}{f}_0\left({\nu }^{†}\right)={\nu }_{\text{*}}^{†}P.V.{\displaystyle \int \nu F_0\left(\nu \right)/{\left({\nu }^{†}-\nu \right)}^2\text{d}\nu }$ .

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 macroextremum point (for the equilibrium function $F_0\left(v\right)$ , 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; we will demonstrate its application to the system of gravitating particles under consideration (this is essential for the 2-dimensional geometry of a system with rotation).

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

$\left(\nu -v_{\parallel }\right)\mathfrak{A}\left(k,v_{\parallel };{\omega }^{‡}\right)={\displaystyle \int {\nu }_{\varkappa }\left(k,v\right)\mathfrak{A}\left(k,v;{\omega }^{‡}\right)\text{d}v},{\displaystyle \int {\nu }_{\varkappa }\left(k,v\right)\mathfrak{A}\left(k,v;{\omega }^{‡}\right)\text{d}v}=1,$ (20)

$\left({\nu }^{‡}-v_{\parallel }\right)\mathfrak{A}\left(k,v_{\parallel };{\omega }^{‡}\right)=1,\left({\nu }^{‡}-{\omega }^{‡}\right)f\left(k,v_{\parallel };{\omega }^{‡}\right)\mathfrak{A}\left(k,v_{\parallel };{\nu }^{‡}\right)=0,$

where normal modes are introduced by the relation $\mathfrak{A}\left(k,v_{\parallel },t\right)=\mathfrak{A}\left(k,v_{\parallel };\omega \right)\text{exp}\left(-i\omega t\right)$ . If the real eigenvalues ${\omega }^{‡}$ are not zeros of the function ${\nu }_{\varkappa }\left(k,v\right)$ , then the eigenfunctions corresponding to them take the form $\mathfrak{A}\left(k,v_{\parallel };{\omega }^{‡}\right)={\nu }^{‡}\left(k,{\omega }^{‡}\right)\delta \left({\omega }^{‡}-v_{\parallel }\right)+P.V.\left(1/\left({\omega }^{‡}-v_{\parallel }\right)\right)$ ; further, we should consider the cases when: 1) ${\omega }^{‡}$ are zeros of the function ${\nu }_{\varkappa }\left(k,v_{\parallel }\right)$ , but not ${\nu }^{‡}\left(k,v_{\parallel }\right)$ ; 2) ${\omega }^{‡}$ are the zeros of the functions ${\nu }_{\varkappa }\left(k,v_{\parallel }\right)$ and ${\nu }^{‡}\left(k,v_{\parallel }\right)$ ; 3) ${\omega }_j^{‡}$ are the complex zeros of $\mathfrak{A}\left(k,v_{\parallel };{\omega }_j^{‡}\right)$ . Finally, we obtain

$\mathfrak{A}\left(k,v_{\parallel };{\omega }^{‡}\right)={\displaystyle \sum _{j}C\left(k,j\right)\mathfrak{A}\left(k,v_{\parallel };{\omega }_j^{‡}\right)}+{\displaystyle \int C\left(k,j\right){\omega }^{‡}\mathfrak{A}\left(k,v_{\parallel };{\omega }^{‡}\right)\text{d}{\omega }^{‡}}.$

The amplitude of the modes is obtained as the sum over the discrete and continuous spectra of the singularities of the functions ${\nu }_{\varkappa }\left(k,v_{\parallel }\right)$ and ${\nu }^{‡}\left(k,v_{\parallel }\right)$ .

Thus, van Kampen waves in the linear approximation for the Poisson equation, with initial conditions that depend only 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 seems interesting to estimate the change the sizes of protostructures during the transition to the phase of gravitational interaction dominance from the point of view of the absence of solutions to Equation (6), since this would lead to the protostructures to a quasi-Jeans type decay (caused by the presence of an additional term, including the cosmological term, in the Liouville-Gelfand equation).

4. Van Kampen Waves in the Case of Non-Uniform Structure of the Initial Field and Its strength. Solutions of the Vlasov Equation of the Bernstein Wave Type

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 [21] [22]. For the simplest 1-dimensional case, the Vlasov equation in coordinates $\left(\mathfrak{E},x,t\right)$ ($\mathfrak{E}=m{v}^2/2+m\Upsilon \left(x\right)$ is the energy of a particle in a gravitational field):

$\frac{\partial F\left(\mathfrak{E},x,t\right)}{\partial t}+v\left(x,\mathfrak{E}\right)\frac{\partial F}{\partial x}+\left(v\left(x,\mathfrak{E}\right)/m\right)\left(G\left(x,t\right)-{\Upsilon }^{\prime }\left(x\right)\right)\frac{\partial F}{\partial \mathfrak{E}}=0,$ (21)

$-\frac{\partial G}{\partial x}=4\pi \gamma N{\displaystyle \int f\left(\mathfrak{E},x,t\right)\text{d}v}-c^2\Lambda$

(the second term on the right-hand side corresponds to the repulsive potential, as before). At equilibrium $f=f_0\left(\mathfrak{E}\right)$ , $\mathfrak{E}=-\text{d}\Upsilon /\text{d}x$ ; if we set $F=F_0\left(\mathfrak{E}\right)+f\left(x,\mathfrak{E},t\right)$ , $G\left(x,t\right)=-{\Upsilon }^{\prime }\left(x\right)+G_1\left(x,t\right)$ , then the linearized Vlasov equation takes the form:

${\left(v\left(x,\mathfrak{E}\right)\right)}^{-1}\frac{\partial f}{\partial t}+\frac{\partial f}{\partial x}-\frac{G_1}{m}\frac{\text{d}{F}_0}{\text{d}\mathfrak{E}}=0.$ (22)

(the repulsive potential is absent in the equation for perturbations, since its effect is present in the basic macropotential $\Upsilon \left(x\right)$ ). We will seek a solution to the equation in the form $f=\psi \left(x\right)\text{exp}\left(-i\omega t\right)$ :

$\frac{\partial \psi }{\partial x}-i\omega v^{-1}\left(x,\mathfrak{E}\right)\psi =G_1/m\cdot \frac{\text{d}{F}_0}{\text{d}\mathfrak{E}},$ (23)

$-i\omega G_1=4\pi \gamma {\displaystyle \int \psi v\text{d}v}=4\pi \gamma {\displaystyle {\int }_{{\mathfrak{E}}_0}^{\infty }\psi \left(x,\mathfrak{E}\right)\text{d}\mathfrak{E}},$

where $v^{-1}\left(x,\mathfrak{E}\right)={\left(2\mathfrak{E}+\Upsilon \right)}^{-1/2}$ . If we exclude $G_1$ from the last two equations, we obtain an equation of the form

$\frac{\partial \psi }{\partial x}-i\omega v^{-1}\left(x,\mathfrak{E}\right)\psi =\left(4\pi i/m\right){\omega }^{-1}\frac{\text{d}{F}_0}{\text{d}\mathfrak{E}}{\displaystyle {\int }_{{\mathfrak{E}}_0}^{\infty }\psi \text{d}\mathfrak{E}}.$ (24)

If $\Upsilon \to 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_1$ and the corresponding frequency ${\Omega }^{\left(0\right)}$ (they are related via the dispersion relation (24)):

$ik\psi \left(k;K_1\right)-i\omega {\displaystyle {\int }_{ℝ}{v}^{-1}\left(q\right)\psi \left(k-q;K_1\right)\text{d}q}=i\left(4\pi {F^{\prime }}_0\gamma /m\right)/\omega {\displaystyle {\int }_{{\mathfrak{E}}_0}^{\infty }\psi \left(k;K_1\right)\text{d}\mathfrak{E}}.$ (25)

This equation can be solved by expanding in powers of the parameter $\Upsilon \left(k\right)/{\mathfrak{E}}_0$ : ${\Theta }_j\left(k\right)={\displaystyle {\sum }_{k=0,\cdots ,\infty }{\Theta }^{\left(j\right)}\left(k\right)}$ , where $\Theta \left(j\right)\left(k\right)\in \left\{v^{-1}\left(k\right),\psi \left(k;K_1\right),\omega \right\}$ . Putting $v^1\left(k\right)=\delta \left(k\right)/\sqrt{2\mathfrak{E}}$ , we obtain in the zeroth approximation two types of eigenmodes, discrete and continuous: ${\psi }^{\left(0\right)}\left(k;K_1\right)={\left(K_1-{\omega }^{\left(0\right)}/\sqrt{2\mathfrak{E}}\right)}^{-1}\left(\left(4\pi {F^{\prime }}_0\gamma /m\right)/{\omega }^{\left(0\right)}\right)\delta \left(k-K_1\right)$ (the criterion for discreteness of the quantities ${\omega }^{\left(0\right)}$ are the conditions $\left(4\pi {F^{\prime }}_0\gamma /m\right)\cdot \left({\left({\omega }^{\left(0\right)}\right)}^2/\left(2K_1^2\right)\right)=0$ , or the condition $Im\left({\omega }^{\left(0\right)}\right)\ne 0$ . If $\left(4\pi {F^{\prime }}_0\gamma /m\right)\cdot \left({\left({\omega }^{\left(0\right)}\right)}^2/\left(2K_1^2\right)\right)\ne 0$ , the functions ${\psi }^{\left(0\right)}\left(k;K_1\right)$ should be considered in the class of distributions, since ${\left(K_1-{\omega }^{\left(0\right)}/\sqrt{2\mathfrak{E}}\right)}^{-1}=P.V.\left(1/\left(K_1-{\omega }^{\left(0\right)}/\sqrt{2\mathfrak{E}}\right)\right)+{\left({\mu }^{‡}\right)}^{\left(0\right)}\left(K_1,{\omega }^{\left(0\right)}\right)\delta \left(K_1-{\omega }^{\left(0\right)}/\sqrt{2\mathfrak{E}}\right)$ (in this case $Im\left({\omega }^{\left(0\right)}\right)\le 0$ indicates the asymptotic stability of the complete solution). In a similar way, one can obtain ${\psi }^{\left(1,2,\cdots \right)}\left(k;K_1\right)$ .

The main result after constructing the appropriate number of terms in the series for $\psi \left(k;K_1\right)$ is the establishment of the density function of the solution of the BGK equations [23] [24]. This expression can be used for comparative calculations of the macroparameters of cosmological objects.

As can be seen from the form of the Equation (24), 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—accordingly, $f^{\left(0\right)}\left(x,v\right)$ )? In accordance with the structure of the Equation (10), the formal substitution of normal modes (of the form ${\mathfrak{U}}_1\left(v\right){\mathfrak{U}}_2\left(x,t\right)$ , in the simplest case ${\mathfrak{U}}_2\left(z,t\right)=\text{exp}\left(ikz-i\omega t\right)$ ) into this equation at ${F|}_{t=0}=F_{MB}\left(x,v\right)$ will lead to the appearance of a bilinear dependence on the spatial and temporal variables, which indicates a non-local 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_0\to 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 [20], we assume that the characteristics of the linear (complete) Vlasov equation coincide with the phase trajectories of the dynamic Hamiltonian system $\text{d}X/\text{d}t=V$ , $\text{d}V/\text{d}t=-\text{d}\Upsilon /\text{d}X$ , since one should consider the additional term $\mathcal{T}\left(\Upsilon ,f\right)\equiv -{{\Upsilon }^{\prime }}_x{f^{\prime }}_v$ on the left-hand side of Equation (11) (the spatial changes in the potential of the “main” gravitational field of the system are taken into account); $\Upsilon$ satisfies Equation (5) (or (6), if after obtaining the solution we pass from the dependent variable $\mathfrak{M}$ to $\Upsilon$ ). The solution of this dynamic system with initial conditions ${X|}_{t=0}=x$ , ${V|}_{t=0}=v$ is as follows: $X\left(x,v,t\right)$ , $V\left(x,v,t\right)$ ($t\in R^1$ ). The first integral of the dynamic system: $\mathfrak{E}=m{v}^2/2+\Upsilon \left(x\right)$ (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 $\mathcal{T}\left(\Upsilon ,f\right)$ was introduced). For the function $f\left(x,v,t\right)$ , through the shift along the trajectories from the initial point, we have the Volterra equation of the $II$ nd kind:

$\begin{array}{c}f\left(x,v,t\right)=f^{\left(0\right)}\left(X\left(x,v,-t\right),V\left(x,v,-t\right)\right)\\ +\frac{\text{d}{F}_0}{\text{d}\mathfrak{E}}{\displaystyle {\int }_0^t\nabla {\Upsilon }^{♭}\left(X\left(x,v,\xi -t\right),\xi \right)V\left(x,v,\xi -t\right)\text{d}\xi },\end{array}$ (26)

and, after substituting this expression into the Poisson equation ${\nabla }^2{\Upsilon }^{♭}={\varsigma }^{♯}\exp \left(\mathfrak{M}\left(x\right)\right)\cdot {\Upsilon }^{♭}\left(x\right)$ , we have an explicit form for the force term ($G\to G\left[\Upsilon \right]+g\left[\Upsilon \right]$ when linearized):

$\begin{array}{c}-{\mathfrak{y}}^{-1}\nabla {\Upsilon }^{♭}\left(x,t\right)={\displaystyle \iint f^{\left(0\right)}\left(X\left(x,v,-t\right),V\left(x,v,-t\right)\right)\text{d}v\text{d}x}\\ +{\displaystyle \iint {\displaystyle {\int }_0^t\frac{\text{d}{F}_0}{\text{d}\mathfrak{E}}\nabla {\Upsilon }^{♭}\left(X\left(x,v,-\xi \right),t-\xi \right)V\left(x,v,-\xi \right)\text{d}\xi \text{d}v\text{d}x}},\end{array}$ (27)

where we use notation $\mathfrak{y}\equiv {\varsigma }^{♯}\exp \left(-\Upsilon \left(x\right)/T\right)$ . In accordance with the definition in formula (6 of the potential value $\mathfrak{M}\left(x\right)$ we obtain for the motion in a non-uniform field of a system of gravitating particles the influence of two integrand factors at once: $\text{d}{F}_0/\text{d}\mathfrak{E}\cdot g$ . This is due to the fact that both $\mathfrak{E}$ 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 $\mathfrak{M}\left(x\right)$ in these factors may lie in different regions $X_i\left(x\right)$ from item 2 (apparently, in order to establish the uniqueness of the solution, the behavior of the function ${v\left(x\right)|}_{\mathfrak{E}=\text{const}}$ should be considered). In addition, the question arises of the physical manifestation of the multivaluedness of solutions to the Vlasov–Poisson equation in the region $X_2\left(x\right)$ : since the norms of the solutions $\mathfrak{M}\left(x\right)$ 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 corresponding 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 p. 3 using the example of van Kampen waves, there may be processes of local “thickening” over time in a certain region of space (antinodes of a longitudinal wave, in particular) of matter, associated with the transition in the region of multivalued solutions of the Liouville-Gelfand 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 $\mathcal{T}\left(\Upsilon ,f\right)$ as $\Upsilon \to \text{const}$ tends can be assumed to be extremely close to the “classical” left-hand side of the linearized Equation (11), however, the right-hand side of the kinetic equation, containing the second term of the right-hand side of the Volterra Equation (26), will retain an unchanged form ($\mathfrak{E}\approx v^2/2+\Upsilon \left(x_0\right)$ ), and this part depends only slightly on the 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 $\left(\omega -\nu \right)$ (taking into account the occurrence of the term in the form of a distribution), and repeat all the operations of p. 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 (corresponding to the previously considered longitudinal waves) for the sake of clarity of the calculations. We integrate both parts (27) over the interval $\left[0,\tilde{z}\right]$ , rearrange the order of integration and make a change of variables ${\eta }_X=X\left(x,v,-\xi \right)$ , ${\eta }_V=V\left(x,v,-\xi \right)$ in the second term of the right-hand side. Since $\mathfrak{E}\left(X,V\right)=\mathfrak{E}\left({\eta }_X,{\eta }_V\right)$ , $\text{d}X\text{d}V=\text{d}{\eta }_X\text{d}{\eta }_V$ , this term will take the form of a flow through the surface: ${\displaystyle \int {\displaystyle {\int }_{\sigma }\mathfrak{E}\left({\eta }_X,t-\xi \right)\partial \left(F_0\left({\eta }_V^2/2+\Upsilon \left({\eta }_X\right)\right)/\partial {\eta }_V\right)\text{d}{\eta }_X\text{d}{\eta }_V}}$ . The boundary $\partial \sigma$ is the image of the line ${\eta }_X$ on the plane $\left({\eta }_X,{\eta }_V\right)$ with a shift in time $-\xi$ along the phase trajectories of the dynamic system of the system $\dot{X}=V$ , $\dot{V}=-{\Upsilon }_X$ (in our case, a small value). We will assume that the boundary $\partial \sigma$ is analytically defined by the relation ${\eta }_V=\beta \left({\eta }_X|\xi ,\tilde{z}\right)$ (${\eta }_V<\beta$ $\forall \left({\eta }_X,{\eta }_V\right)\in \sigma$ ). Then the second term under study will take the form ${\displaystyle \int g\left({\eta }_X,t-\xi \right)F_0\left(\mathfrak{E}\left[{\eta }_X,{\eta }_V\right]\right)\text{d}{\eta }_X}$ . Therefore, the right-hand side of (27) has the form:

$\mathfrak{g}\left(\tilde{z},t;F_0\right)\equiv {\displaystyle {\int }_0^{\tilde{z}}{g}_0\left(z,t\right)\text{d}z}+{\displaystyle \int {\displaystyle {\int }_0^{t}g\left(t-\xi ,{\eta }_X\right)F_0\left({\beta }^2\left({\eta }_X|\xi ,\tilde{z}\right)+\Upsilon \left({\eta }_X\right)\right)\text{d}\xi \text{d}{\eta }_X}}$(28)

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

$f\left(k,v;\omega \right)=-k^{-1}\mathfrak{S}\left(z,t;v\right)\cdot P.V.{\left(1/\left(\nu -v\right)\right)|}_{\nu =\omega /k}+{\vartheta \left(k,\nu \right)\delta \left(\nu -v\right)|}_{\nu =\omega /k},$

where $\vartheta \left(k,\nu \right)$ —normalization function ($\vartheta =1+{\displaystyle \int \left(-k^{-1}\mathfrak{S}\left(z,t;v\right)/\left(v-\nu \right)\right)\text{d}v}$ ).

The solution of the initial value problem $f\left(k,v,t\right)$ can be represented as an expansion in special solutions $f\left(k,v;\nu \right)\text{exp}\left(-ik\nu t\right)$ : $f\left(k,v,t\right)={\displaystyle \int \mathcal{U}\left(k,\nu \right)f\left(k,v;\nu \right)\text{exp}\left(-ik\nu t\right)\text{d}\nu }$ ; accordingly, the Cauchy condition $f^{\left(0\right)}\left(k,v\right)\equiv f\left(k,v,t=0\right)={\displaystyle \int \mathcal{U}\left(k,\nu \right)f\left(k,v;\nu \right)\text{d}\nu }$ . To determine the coefficients of $\mathcal{U}$ , we obtain a singular integral equation:

$\mathcal{U}\left(k,v\right)=-k^{-1}\mathfrak{S}\left(z,t;v\right)\cdot P.V.{\displaystyle \int \left(\mathcal{U}\left(k,\nu \right)/\left(\nu -v\right)\right)\text{d}\nu }+\vartheta \left(k,v\right)\mathcal{U}\left(k,v\right).$

His solution looks like:

$\mathcal{U}\left(k,\nu \right)=\frac{{\mathcal{G}}_{+}\left(k,\nu \right)}{1+2\pi i{\mathscr{H}}_{+}\left(k,\nu \right)}-\frac{{\mathcal{G}}_{-}\left(k,\nu \right)}{1+2\pi i{\mathscr{H}}_{-}\left(k,\nu \right)},$

${\mathcal{G}}_{+}\left(k,\nu \right)-{\mathcal{G}}_{-}\left(k,\nu \right)=\mathcal{U}\left(k,v\right),{\mathcal{G}}_{+}\left(k,\nu \right)+{\mathcal{G}}_{-}\left(k,\nu \right)=\frac{1}{\pi }{\displaystyle \int \frac{\mathcal{U}\left(k,\nu \right)}{\nu -v}\text{d}\nu },$

${\mathscr{H}}_{+}\left(k,\nu \right)-{\mathscr{H}}_{-}\left(k,\nu \right)=-k^{-1}\mathfrak{S}\left(z,t;v\right),$

${\mathscr{H}}_{+}\left(k,\nu \right)+{\mathscr{H}}_{-}\left(k,\nu \right)=\frac{1}{\pi i}{\displaystyle \int \frac{-k^{-1}\mathfrak{S}\left(z,t;v\right)}{\nu -v}\text{d}\nu }.$

Thus, we have obtained a method for applying van Kampen waves to a formally weakly inhomogeneous system of particles (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 $\mathfrak{S}\left(z,t;v\right)$ is defined through the relation (28) and contains the factor $\mathfrak{y}$ . Recall that in the second term on the right-hand side (28) there is a “general potential” $\Upsilon \left(x\right)$ , which is a solution to the nonlinear Poisson Equation (5), in which the influence of antigravity is taken into account (the term with the cosmological term): $\Delta \Upsilon \left(x\right)={\varsigma }^{†}\exp \left(-\Upsilon \left(x\right)/T\right)-c^2\Lambda$ . Further, the quantity $\mathfrak{y}$ is defined as $\mathfrak{y}\equiv {\varsigma }^{♯}\exp \left(-\Upsilon \left(x\right)/T\right)$ , where, in turn, the pre-exponential factor ${\varsigma }^{♯}=-\frac{{\varsigma }^{†}}{T}\exp \left(c^2\Lambda x^2/T\right)$ , 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 most important factor that requires modification of standard approaches such as van Kampen waves and Landau damping (as a consequence of the expansion of Landau modes in van Kampen waves).

5. Conclusion

The application of plasma theory concepts, including Landau damping, van Kampen and Case waves, and the normal mode method, to describe phenomena and processes in astrophysical conditions are of considerable interest both in terms of searching for manifestations of known physical aspects of plasma oscillations, resonances, and nonlinear interactions of waves by observation, and for identifying new patterns in known experimental material—the interpretation of observations may well be not entirely legitimate, obscured by statistical noise. Therefore, the development and application of approaches in cosmology for which a powerful mathematical apparatus has already been developed has enormous practical meaning. In this paper, a method for applying wave processes in plasma associated with invariant properties of the linearized Vlasov-Poisson equation is proposed and implemented. It is established that for gravitational interaction, including antigravity (due to the inclusion of the cosmological term in the considerations), the structures arising in the process of wave motion have very nontrivial dynamic properties associated with solutions of the Liouville-Gelfand equation. Motion in a quasi-homogeneous gravitational field is very similar in terms of description methods to an electromagnetic field in plasma, however, when taking into account the inhomogeneity of the initial field and small gravitational perturbations, the situation changes fundamentally. Direct analogues of van Kampen waves can be constructed only locally, in the case of an extremum of global interaction in a system of separated masses (this is due to the fact that in the vicinity of the extremum of the global field the potential is close to a constant value, which means the absence of external and self-consistent forces; due to this, the distribution function should be Maxellian equilibrium, which strictly corresponds the applicability limit of the unmodified normal mode method).