# Extended relativistic kinetic model composed of the scalar and two vector
distribution functions: Application to the spin-electron-acoustic waves

Pavel A. Andreev<EMAIL_ADDRESS>Department of General Physics,
Faculty of physics, Lomonosov Moscow State University, Moscow, Russian
Federation, 119991.

###### Abstract

Detailed deterministic derivation of the kinetic equations for the
relativistic plasmas is given. Focus is made on the dynamic of one-coordinate
distribution functions of various tensor dimensions, but the closed set of
kinetic equations is constructed of three functions: the scalar distribution
function, the vector distribution function of dipole moment, and the vector
distribution function of velocity (or the dipole moment in the momentum
space). All two-coordinate distributions functions are discussed as well. They
are presented together with their limits existing in the self-consistent field
approximation. The dynamics of the small amplitude spin-electron-acoustic
waves in the dense degenerate plasmas is are studied within the kinetic model.
This work presents the deterministic approach to the derivation and
interpretation of the kinetic equations. So, no probability is introduced
during the transition from the level of individual particles to the collective
functions. Problem of thermalization is not considered, but we can see that
the structure of kinetic equations is kept for the system before and after
thermalization. Hence, the kinetic equations can be used to approach this
item.

spin-electron-acoustic waves, microscopic model, kinetic theory, mean-field
approximation, relativistic plasmas

##  I Introduction

Kinetic theory is one of fundamental theories of the macroscopic phenomena in
various physical systems Klimontovich book , Israel AP 79 , Romatschke IJMPE
10 , Kovtun JP A 12 , Denicol PRD 12 , Heinz ARNPS 13 , Florkowski RPP 18 ,
Denicol PRD 19 . Particularly, is plays essential role in the plasma physics
Hakim book Rel Stat Phys . Kinetic theory requires proper microscopic
derivation, where the evolution of the distribution function is traced from
the evolution of individual particles. The distribution function shows the
relative positions of all particles in the six dimensional phase space and it
can be consistently defined in the deterministic way. So, no probability
theory is applied for the derivation or interpretation of the kinetic
equations or the distribution function.

Such view on the kinetic theory is the extension of the classical
hydrodynamics based on the tracing of the microscopic dynamics of individual
particles Drofa TMP 96 , Kuzmenkov CM 15 .

Classical mechanics gives us dynamics of particles in the 3N-dimensional
configurational space. The model of motion of particles is composed of the
individual trajectories. Each of them, formally, takes place in its own three
dimensional space. this picture of motion is not obvious if we consider the
Newtonian form of mechanics, where we can use our own imagination to place all
particles in the single three-dimensional space. However, the Lagrangian and
Hamiltonian forms of mechanics show the multidimensional structure of the
theory more clearly.

Moreover, the construction of the mathematical model out of trajectories
(which are one dimensional mathematical objects) distinguish the classical
mechanics from the hydrodynamics or the electrodynamics, which are the field
theories.

The derivation of the hydrodynamics from the mechanics requires to represent
the mechanics in the form of the field theory. So, the mechanics would have
the same mathematical structure as the electrodynamics, which is essentially
important for the dynamics of charged particles.

Refs. Drofa TMP 96  and Kuzmenkov CM 15  basically present the field form of
the classical mechanics, where the dynamics of particles takes form of the
dynamics of material fields. Inevitably, the field form of classical mechanics
appears in the form of hydrodynamic equations. These are equations of
nonequilibrium hydrodynamics which requires further reduction for the
particular forms of the collective motion.

Particularly, it leads to the some generalized form of nonequilibrium
thermodynamics. And choosing regime of thermal equilibrium allows to get the
first low of thermodynamics from the energy evolution equation. The
possibility to consider different regimes, like regimes close to the thermal
equilibrium, and regimes out of equilibrium, shows that this approach can be
used for the study of the thermalization process in the physical systems
considering the dynamics of the systems in terms of the collective variables
(material fields). It would require proper truncation of the hydrodynamics out
equilibrium considering higher rank material fields to include the mechanism
of relaxation in the model. However, it is an open problem which can be
addressed within this formulation.

Formulation of the mechanics in terms of the three-dimensional material fields
leads to the hydrodynamics form of equations. While we consider the relative
positions and relative momentums of all particles we go to the distribution of
particles in the six dimensional space (the phase space of coordinate and
momentum). This form of the classical mechanics is considered in this paper.

If we consider the microscopic evolution of point-like particles we get the
following distribution of the particle number in the coordinate space Weinberg
Gr 72

$n_{m}(\textbf{r},t)=\sum_{i=1}^{N}\delta(\textbf{r}-\textbf{r}_{i}(t)),$ (1)

where $\textbf{r}_{i}(t)$ is the coordinate of $i$-th particle. This approach
is suggested by Yu.L. Klimontovich Klimontovich book , Weinberg Gr 72 . While
we use notation $n_{m}(\textbf{r},t)$ on the left-hand side of equation (1) we
should write $n_{m}(\textbf{r},\textbf{r}_{1}(t),...,\textbf{r}_{N}(t))$. So,
the time dependence of function $n_{m}$ is directly constructed out of the
time evolution of coordinates of particles. The parametrization of the
physical space is introduced here r. So this form gives us an intermediate
step from the mechanics of trajectories to the field form of classical
mechanics.

Considering the microscopic distribution (1) different authors applied some
form of averaging (see for instance Klimontovich book , deGrootFoundations1972
). Sometimes, authors use some unspecified distributions $\langle
n_{m}\rangle$ (see for instance Klimontovich book ). Overwise authors use
specific form of averaging via some distribution function (see for instance
deGrootFoundations1972 , Kuzelev UFN 99 ). Anyway, authors try to impose
probabilistic interpretation of the procedure of averaging. Let us mention
that the aim of averaging is to make the transition to the macroscopic scale
by averaging on the physically infinitesimal volume, but usually it is being
replaced by the averaging on some distribution function.

Systematic, nonstatistical, definition of the averaging on the physically
infinitesimal volume is suggested by L.S. Kuz’menkov Drofa TMP 96  (see also
Andreev PIERS 2012 , Andreev 2021 05 )

$n_{e}(\textbf{r},t)=\frac{1}{\Delta}\int_{\Delta}d\mbox{\boldmath$\xi$}\sum_{i=1}^{N/2}\delta(\textbf{r}+\mbox{\boldmath$\xi$}-\textbf{r}_{i}(t)).$
(2)

Equation (2) contains vector $\xi$ which scanning the physically infinitesimal
volume (it is illustrated in equation 1). While the physically infinitesimal
volume appears as the $\Delta$-vicinity of each point of space. It is assumed
that $N$ is the full number of particles in the system. It is composed of two
species: electrons with numbers $i\in[1,N/2]$ and protons $i\in[N/2+1,N]$. We
need to trace the evolution of each species in the system. With no
restrictions, we illustrate our derivation on the subsystem of electrons, so
subindex $e$ is dropped for all functions describing electrons in equations
below $n_{e}\equiv n$.

Figure 1:  The illustration of the $\Delta$-vicinity in the coordinate space.

Similar background is suggested for the kinetic theory. So, the microscopic
distribution function is composed of delta functions in the physical
coordinate space and the momentum space

$f_{m}(\textbf{r},\textbf{p},t)=\sum_{i=1}^{N}\delta(\textbf{r}-\textbf{r}_{i}(t))\delta(\textbf{p}-\textbf{p}_{i}(t)).$
(3)

Transition to the macroscopic level is made similarly to equation (2). We
introduce the physically infinitesimal volume in six-dimensional phase space

$f(\textbf{r},\textbf{p},t)=\frac{1}{\Delta}\frac{1}{\Delta_{p}}\times$
$\times\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\sum_{i=1}^{N/2}\delta(\textbf{r}+\mbox{\boldmath$\xi$}-\textbf{r}_{i}(t))\delta(\textbf{p}+\mbox{\boldmath$\eta$}-\textbf{p}_{i}(t)),$
(4)

where $d\mbox{\boldmath$\eta$}$ is the element of volume in the momentum
space, with $\int_{\Delta_{p}}d\mbox{\boldmath$\eta$}$ integral over
$\Delta_{p}$-vicinity in the momentum space, $\Delta\equiv\Delta_{r}$ is the
delta vicinity in the coordinate space. Recent development of this approach
can be found in Refs. Andreev Ph 15 , Zubkov JPhCS 19 , Zubkov 2006 , Zubkov
VNSU 22 .

This discussion is about the field form of the classic mechanics. However,
similar insight is imposed on the quantum mechanics. So, the formulation of
the many-particle quantum mechanics in the form of evolution of the material
fields is developed Maksimov QHM 99 , MaksimovTMP 2001 , Andreev LP 19 ,
Andreev LP 21 fermions , Andreev LPL 21 , Andreev JPP 21 , Andreev Ch 21 ,
Andreev PoF 21 .

The discussion presented above is focused on the fundamental background of the
hydrodynamics and kinetics, since the theoretical part of this paper is based
on the development of the kinetic theory. However, the application of the
kinetic theory to the high-density magnetized spin polarized degenerate
electron gas is also considered.

Particularly, we focus on the spin-electron-acoustic waves (SEAWs). Existence
of the SEAWs in the partially spin-polarized degenerate plasmas was
theoretically suggested in 2015 Andreev PRE 15 . They appear as the
longitudinal waves of density and electric field, which are propagate parallel
to the anisotropy direction in the magnetized electron gas Andreev PRE 15 ,
Andreev PP 16 SSE kin .

If we have no spin polarization we find single longitudinal wave in the
electron gas under assumption of the motionless ions. It is the Langmuir wave,
which corresponds to the collective oscillations of electrons relatively
motionless ions. If we include the spin-polarization we also find the SEAWs.
The physical picture of dynamics is following. It corresponds to relative
motion of electrons with different spin projections. The spin polarization
manifests itself in different partial concentrations of electrons with the
particular spin polarization. Hence, the relative oscillations of electrons
with different amplitudes of concentration give the oscillation of small
proportion of electrons relatively motionless ions as well Andreev PRE 15 ,
Andreev PP 16 SSE kin . The SEAWs show similarity to the spin-plasmons
considered in two-dimensional structures Ryan PRB 91 , Agarwal PRL 11 ,
Agarwal PRB 14 , Schober EPL 20 , Amico JPD 19 , Yang Ph 21 , Kreil CPP 18 ,
Kreil PRB 15 , Huang PRB 16 , Afanasiev PRB 22 . Recently, the SEAWs have been
considered in the degenerate electron gas Andreev 2112 , while the background
of applied model is given in Refs. Andreev 2021 05 , Andreev 2021 09 , Andreev
2021 11 .

Increase of the concentration of electrons $n_{0e}$ up to values, where the
Fermi energy is proportional to the rest energy of electron
$\varepsilon_{Fe}=(3\pi^{2}n_{0e})^{2/3}\hbar^{2}/2m_{e}$$\sim m_{e}c^{2}$,
gives noticeable change in the dispersion dependence of the SEAWs. It also
corresponds to the change of the dispersion dependence of the Langmuir waves.
It shows the decrease of the relative frequency of the Langmuir wave and SEAW.
Moreover, the relativistic regime shows the growth of the amplitude of the
concentration of the SEAWs relatively the amplitude of the Langmuir waves for
the chosen value of the electric field in these waves.

This paper is organized as follows. In Sec. II the method of derivation of the
kinetic equation from the microscopic motion is demonstrated and the general
structure of the kinetic equation is derived. In Sec. III we present an
approximate kinetic model, where the contribution of physically infinitesimal
volume is considered in minimal regime called the monopole regime. In Sec. IV
the selfconsistent field approximation of monopole regime is considered. In
Sec. V the multipole approximation in the relativistic kinetic equation is
presented as the generalization of the monopole regime considered in Sec. III.
In Sec. VI the selfconsistent field approximation of the multipole
approximation is presented. In Sec. VII the kinetic equation for the evolution
of the dipole moment vector distribution function is obtained. In Sec. VIII
the kinetic equation for the distribution function of velocity is derived. In
Sec. IX closed set of three kinetic equations is discussed. In Sec. X the
spin-electron-acoustic waves in the spin polarized electron gas of high
density are considered. In Sec. XI a brief summary of obtained results is
presented.

##  II Derivation of the Vlasov equation tracing the microscopic motion of
particles

We need to consider the evolution of the distribution function. It requires to
present equation of motion of each particle in full details

$\dot{\textbf{p}}_{i}(t)=\textbf{F}(\textbf{r}_{i}(t),t),$ (5)

where
$\textbf{p}_{i}(t)=m_{i}\textbf{v}_{i}(t)/\sqrt{1-\textbf{v}_{i}^{2}(t)/c^{2}}$
is the relativistic momentum of $i$-th particle, and
$\textbf{F}(\textbf{r}_{i}(t),t)$ is the force acting on $i$-th particle from
the electromagnetic field created by other particles in the system. The
equation of motion for each particle appears as the evolution of the momentum
under action of the Lorentz force

$\textbf{F}(\textbf{r}_{i}(t),t)=\biggl{(}q_{i}\textbf{E}(\textbf{r}_{i}(t),t)+\frac{1}{c}q_{i}[\textbf{v}_{i}(t),\textbf{B}(\textbf{r}_{i}(t),t)]\biggr{)},$
(6)

where $\textbf{E}_{i}=\textbf{E}_{i,ext}+\textbf{E}_{i,int}$,
$\textbf{B}_{i}=\textbf{B}_{i,ext}+\textbf{B}_{i,int}$,
$\textbf{E}_{i}=\textbf{E}(\textbf{r}_{i}(t),t)$, and
$\textbf{B}_{i}=\textbf{B}(\textbf{r}_{i}(t),t)$, while the external fields
are included along with the field of interaction of particles. The electric
$\textbf{E}_{i,int}$ and magnetic $\textbf{B}_{i,int}$ fields caused by
particles surrounding the $i$-th particle are
$\textbf{E}_{i,int}=-\nabla_{i}\varphi(\textbf{r}_{i}(t),t)-\frac{1}{c}\partial_{t}\textbf{A}(\textbf{r}_{i}(t),t)$
and $\textbf{B}_{i,int}=\nabla_{i}\times\textbf{A}(\textbf{r}_{i}(t),t)$ with
Landau v2

$\varphi(\textbf{r}_{i}(t),t)=\sum_{j\neq
i}q_{j}\int\frac{\delta(t-t^{\prime}-\frac{1}{c}\mid\textbf{r}_{i}(t)-\textbf{r}_{j}(t^{\prime})\mid)}{\mid\textbf{r}_{i}(t)-\textbf{r}_{j}(t^{\prime})\mid}dt^{\prime},$
(7)

and

$\textbf{A}(\textbf{r}_{i}(t),t)=\sum_{j\neq
i}q_{j}\int\frac{\delta(t-t^{\prime}-\frac{1}{c}\mid\textbf{r}_{i}(t)-\textbf{r}_{j}(t^{\prime})\mid)}{\mid\textbf{r}_{i}(t)-\textbf{r}_{j}(t^{\prime})\mid}\frac{\textbf{v}_{j}(t^{\prime})}{c}dt^{\prime}.$
(8)

Equations (7) and (8) allow to introduce the Green function of the retarding
electromagnetic interaction Landau v2

$\tilde{G}_{ij}=\frac{\delta(t-t^{\prime}-\frac{1}{c}\mid\textbf{r}_{i}(t)-\textbf{r}_{j}(t^{\prime})\mid)}{\mid\textbf{r}_{i}(t)-\textbf{r}_{j}(t^{\prime})\mid}.$
(9)

We consider the high-temperature plasmas , where electrons (and may be other
species) have relativistic temperatures comparable with the rest energy of the
electron $m_{e}c^{2}$ (of the particle of corresponding species). Obtaining
such huge temperatures of the species would lead to high degree of ionization
of the atomic objects. Moreover, the two-particle interactions (collision-like
processes) of the high-energy electrons with the ions would change the
electron configurations of ions or provide additional ionization. These
dynamical processes should contribute in the model. To avoid these
complexifications we consider the hydrogen plasmas, where we have electrons
and protons only.

To find the equation for evolution of the distribution function (4) we
consider the time derivative of this function to obtain the following equation

$\partial_{t}f(\textbf{r},\textbf{p},t)+\nabla\cdot\frac{1}{\Delta}\frac{1}{\Delta_{p}}\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\sum_{i=1}^{N/2}\dot{\textbf{r}}_{i}(t)\delta(\textbf{r}+\mbox{\boldmath$\xi$}-\textbf{r}_{i}(t))\delta(\textbf{p}+\mbox{\boldmath$\eta$}-\textbf{p}_{i}(t))$
$+\nabla_{\textbf{p}}\cdot\frac{1}{\Delta}\frac{1}{\Delta_{p}}\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\sum_{i=1}^{N/2}\dot{\textbf{p}}_{i}(t)\delta(\textbf{r}+\mbox{\boldmath$\xi$}-\textbf{r}_{i}(t))\delta(\textbf{p}+\mbox{\boldmath$\eta$}-\textbf{p}_{i}(t))=0.$
(10)

The second (third) term in this equation appears as the result of action of
the time derivative on the delta-function containing the coordinate (the
momentum).

The second term in equation (10) contains the following function

$\langle\textbf{v}_{i}(t)\rangle\equiv\frac{1}{\Delta}\frac{1}{\Delta_{p}}\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\sum_{i=1}^{N/2}\textbf{v}_{i}(t)\delta_{\textbf{r}i}\delta_{\textbf{p}i},$
(11)

where $\textbf{v}_{i}(t)=d\textbf{r}_{i}(t)/dt$,
$\delta_{\textbf{r}i}\equiv\delta(\textbf{r}+\mbox{\boldmath$\xi$}-\textbf{r}_{i}(t))$,
and
$\delta_{\textbf{p}i}\equiv\delta(\textbf{p}+\mbox{\boldmath$\eta$}-\textbf{p}_{i}(t))$.
For the further representation of this function we apply relativistic
expression for the velocity of particle via its momentum
$\textbf{v}_{i}(t)=\textbf{p}_{i}(t)c/\sqrt{\textbf{p}_{i}^{2}(t)+m_{i}^{2}c^{2}}$,
where additional replacement of the momentum $\textbf{p}_{i}(t)$ on
$\textbf{p}+\mbox{\boldmath$\eta$}$ can be applied. It gives the
representation of function $\langle\textbf{v}_{i}(t)\rangle$:

$\langle\textbf{v}_{i}(t)\rangle=\frac{1}{\Delta}\frac{1}{\Delta_{p}}\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\sum_{i=1}^{N/2}\frac{(\textbf{p}+\mbox{\boldmath$\eta$})c}{\sqrt{(\textbf{p}+\mbox{\boldmath$\eta$})^{2}+m_{i}^{2}c^{2}}}\delta_{\textbf{r}i}\delta_{\textbf{p}i}.$
(12)

Here, our goal is the extraction of traditional for the physical kinetics
terms $\textbf{v}\cdot f$. Hence, expression (12) is represented in required
form

$\langle\textbf{v}_{i}(t)\rangle=\textbf{v}\cdot
f(\textbf{r},\textbf{p},t)+\textbf{F}(\textbf{r},\textbf{p},t),$ (13)

where

$\textbf{F}(\textbf{r},\textbf{p},t)\equiv\langle\Delta\textbf{v}_{i}(t)\rangle$
(14)

with

$\langle\Delta\textbf{v}_{i}\rangle\equiv\frac{1}{\Delta}\frac{1}{\Delta_{p}}\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\sum_{i=1}^{N/2}\biggl{(}\frac{\mbox{\boldmath$\eta$}c}{\sqrt{\textbf{p}^{2}+m_{i}^{2}c^{2}}}$
$-\frac{\textbf{p}(\textbf{p}\cdot\mbox{\boldmath$\eta$})c}{(\sqrt{\textbf{p}^{2}+m_{i}^{2}c^{2}})^{3}}\biggr{)}\delta_{\textbf{r}i}\delta_{\textbf{p}i},$
(15)

with $\textbf{p}=m_{s}\textbf{v}/\sqrt{1-\textbf{v}^{2}/c^{2}}$ and
$\textbf{v}=\textbf{p}c/\sqrt{\textbf{p}^{2}+m_{s}^{2}c^{2}}$.

We include representation (13) in equation (10). We also include equations of
motion of the individual particles (5)-(8) in the last term in equation (10):

$\partial_{t}f(\textbf{r},\textbf{p},t)+(\textbf{v}\cdot\nabla)f(\textbf{r},\textbf{p},t)+\nabla\cdot\textbf{F}(\textbf{r},\textbf{p},t)+\frac{q_{s}}{m_{s}}\frac{1}{\Delta}\frac{1}{\Delta_{p}}\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\sum_{i=1}^{N/2}\biggl{[}\textbf{E}_{ext}(\textbf{r}+\mbox{\boldmath$\xi$},t)$
$+\frac{1}{c}[(\textbf{v}+\Delta\textbf{v}_{i})\times\textbf{B}_{ext}(\textbf{r}+\mbox{\boldmath$\xi$},t)]+q_{s^{\prime}}\sum_{j=1,j\neq
i}^{N}\int
dt^{\prime}\biggl{(}-\nabla_{\textbf{r}}-\frac{1}{c}\frac{\textbf{v}_{j}(t^{\prime})}{c}\partial_{t}+\frac{1}{c^{2}}[\textbf{v}_{i}(t)\times[\nabla_{\textbf{r}}\times\textbf{v}_{j}(t^{\prime})]]\biggr{)}G(\textbf{r}+\mbox{\boldmath$\xi$}-\textbf{r}_{j}(t))\biggr{]}\delta_{\textbf{r}i}\cdot\nabla_{\textbf{p}}\delta_{\textbf{p}i}=0,$
(16)

where we use symbol $\hat{\partial_{t}}$ instead of the time derivative due to
the change of the structure of argument of the Green function $G$. We replaced
$\textbf{r}_{i}(t)$ by $\textbf{r}+\mbox{\boldmath$\xi$}$ in the arguments of
the external electric field, the external magnetic field, and the Green
function using the delta-function $\delta_{\textbf{r}i}$. Initially the action
of the time derivative has the following form

$\partial_{t}G(t,t^{\prime},\textbf{r}_{i}(t),\textbf{r}_{j}(t^{\prime}))=\frac{\delta^{\prime}}{\mid\textbf{r}_{i}(t)-\textbf{r}_{j}(t^{\prime})\mid},$
(17)

where $\delta^{\prime}$ is the derivative of the delta function on its
argument. It is also can be rewritten in the following form

$\partial_{t}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})=\frac{\delta^{\prime}}{\mid\textbf{r}+\mbox{\boldmath$\xi$}-\textbf{r}^{\prime}-\mbox{\boldmath$\xi$}^{\prime}\mid},$
(18)

For the introduction of $\mbox{\boldmath$\xi$}^{\prime}$ see the equations
below.

Next, we use representation (13) in terms describing the interaction:

$\partial_{t}f+\textbf{v}\cdot\nabla
f+\nabla\cdot\textbf{F}+\frac{q_{s}}{m_{s}}\frac{1}{\Delta}\frac{1}{\Delta_{p}}\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\sum_{i=1}^{N/2}\biggl{(}\textbf{E}(\textbf{r}+\mbox{\boldmath$\xi$},t)+\frac{1}{c}[(\textbf{v}+\Delta\textbf{v}_{i})\times\textbf{B}(\textbf{r}+\mbox{\boldmath$\xi$},t)]\biggr{)}\delta_{\textbf{r}i}\cdot\nabla_{\textbf{p}}\delta_{\textbf{p}i}$
$+\frac{q_{s}}{m_{s}}q_{s^{\prime}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int_{\Delta,\Delta_{p}}\frac{d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}d\mbox{\boldmath$\xi$}^{\prime}d\mbox{\boldmath$\eta$}^{\prime}}{\Delta^{2}\Delta_{p}^{2}}\sum_{i=1}^{N/2}\sum_{j=1,j\neq
i}^{N}\int
dt^{\prime}\biggl{(}-\nabla_{\textbf{r}}-\frac{1}{c}\frac{(\textbf{v}^{\prime}+\Delta\textbf{v}_{j}^{\prime})}{c}\partial_{t}$
$-\frac{1}{c^{2}}[(\textbf{v}+\Delta\textbf{v}_{i})\times[(\textbf{v}^{\prime}+\Delta\textbf{v}_{j}^{\prime})\times\nabla_{\textbf{r}}]]\biggr{)}G(\textbf{r}+\mbox{\boldmath$\xi$}-\textbf{r}^{\prime}-\mbox{\boldmath$\xi$}^{\prime})\delta_{\textbf{r}i}\cdot\nabla_{\textbf{p}}\delta_{\textbf{p}i}\delta_{\textbf{r}^{\prime}j}\delta_{\textbf{p}^{\prime}j}=0,$
(19)

where
$\delta_{\textbf{r}^{\prime}j}\equiv\delta(\textbf{r}^{\prime}+\mbox{\boldmath$\xi$}^{\prime}-\textbf{r}_{j}(t))$,
and
$\delta_{\textbf{p}^{\prime}j}\equiv\delta(\textbf{p}^{\prime}+\mbox{\boldmath$\eta$}^{\prime}-\textbf{p}_{j}(t))$.

Presented kinetic equation (19) contains the deviation of coordinates $\xi$
from the center of the $\Delta$-vicinity r. It also contain the deviation of
the velocity $\Delta\textbf{v}_{i}$ from the value corresponding to the center
of the $\Delta_{\textbf{p}}$-vicinity in the momentum space
$\textbf{v}=\textbf{p}c/\sqrt{\textbf{p}^{2}+m_{s}^{2}c^{2}}$. This is the
intermediate general representation of the kinetic equation which is
considered below under some additional assumptions.

##  III Monopole approximation of the kinetic equation for the scalar
distribution function

To consider the monopole approximation of the kinetic equation we need to
neglect contribution of $\xi$, $\mbox{\boldmath$\xi$}^{\prime}$, $\eta$
(including $\Delta\textbf{v}_{i}$), and $\mbox{\boldmath$\eta$}^{\prime}$
(including $\Delta\textbf{v}_{j}^{\prime}$) in dynamical functions. However,
on this step, we present the kinetic equation, where the monopole
approximation is considered for the space variables only. Hence, we neglect
contribution of $\xi$, and $\mbox{\boldmath$\xi$}^{\prime}$, but we keep the
contribution of $\eta$:

$\partial_{t}f+\textbf{v}\cdot\nabla
f+\nabla\cdot\textbf{F}+\frac{q_{s}}{m_{s}}\biggl{(}\textbf{E}(\textbf{r},t)+\frac{1}{c}[\textbf{v}\times\textbf{B}(\textbf{r},t)]\biggr{)}\cdot\nabla_{\textbf{p}}f$
$+\frac{q_{s}}{m_{s}c}(\nabla_{\textbf{p}}\cdot[\textbf{F}(\textbf{r},\textbf{p},t)\times\textbf{B}(\textbf{r},t)])$
$+\frac{q_{s}}{m_{s}}q_{s^{\prime}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int_{\Delta,\Delta_{p}}\frac{d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}d\mbox{\boldmath$\xi$}^{\prime}d\mbox{\boldmath$\eta$}^{\prime}}{\Delta^{2}\Delta_{p}^{2}}\sum_{i=1}^{N/2}\sum_{j=1,j\neq
i}^{N}\int dt^{\prime}\biggl{(}-\nabla_{\textbf{r}}$
$-\frac{\textbf{v}^{\prime}}{c^{2}}\partial_{t}+\frac{1}{c^{2}}[\textbf{v}\times[\nabla\times\textbf{v}^{\prime}]]\biggr{)}G(\textbf{r},\textbf{r}^{\prime})\delta_{\textbf{r}i}\cdot\nabla_{\textbf{p}}\delta_{\textbf{p}i}\delta_{\textbf{r}^{\prime}j}\delta_{\textbf{p}^{\prime}j}=0.$
(20)

Transition from equation (19) to equation (20) includes the transformation
like
$\textbf{E}(\textbf{r}+\mbox{\boldmath$\xi$},t)\approx\textbf{E}(\textbf{r},t)$
and

$\frac{1}{\Delta}\frac{1}{\Delta_{p}}\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\sum_{i=1}^{N/2}\textbf{E}(\textbf{r},t)\delta_{\textbf{r}i}\cdot\nabla_{\textbf{p}}\delta_{\textbf{p}i}=\textbf{E}(\textbf{r},t)\cdot\nabla_{\textbf{p}}f,$
(21)

where the electric field $\textbf{E}(\textbf{r},t)$ is placed outside of the
integral, while the rest is the derivative of the distribution function on the
momentum.

Equation (20) allows to introduce the two-particle (or in other terms two-
coordinate) distribution functions. To the best of my knowledge, majority of
papers and books, including my own works, use notion ”two-particle
distribution functions”. It can be confusing to some extend since the model
describes the many-particle systems. Hence, the distribution function and the
two-particle distribution function actually describe the dynamics of whole
system. Sometimes, the two-particle distribution functions are called ”two-
coordinate distribution functions”. It looks more logical since it contains
two coordinates r and $\textbf{r}^{\prime}$ (or more exactly two coordinates
in the phase space $\\{\textbf{r},\textbf{p}\\}$ and
$\\{\textbf{r}^{\prime},\textbf{p}^{\prime}\\}$).

Let us present the following kinetic equation in the monopole approximation
both on coordinate and momentum (we do it due to the technical reasons, since
we do not want to introduce several two-coordinate distribution functions, but
we show the complete picture below):

$\partial_{t}f+\textbf{v}\cdot\nabla
f+\nabla\cdot\textbf{F}+\frac{q_{s}}{m_{s}}\biggl{(}\textbf{E}+\frac{1}{c}[\textbf{v}\times\textbf{B}]\biggr{)}\cdot\nabla_{\textbf{p}}f$
$+\frac{q_{s}}{m_{s}c}(\nabla_{\textbf{p}}\cdot[\textbf{F}\times\textbf{B}])+\frac{q_{s}}{m_{s}}q_{s^{\prime}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\biggl{(}-\nabla_{\textbf{r}}$
$-\frac{\textbf{v}^{\prime}}{c^{2}}\partial_{t}+\frac{1}{c^{2}}[\textbf{v}\times[\nabla\times\textbf{v}^{\prime}]]\biggr{)}G(\textbf{r},\textbf{r}^{\prime})\cdot\nabla_{\textbf{p}}f_{2}=0,$
(22)

where $q_{s^{\prime}}f_{2}=q_{e}f_{2,ee}+q_{i}f_{2,ei}$ and

$f_{2,ee}(\textbf{r},\textbf{r}^{\prime},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})=\int_{\Delta,\Delta_{p}}\frac{d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}d\mbox{\boldmath$\xi$}^{\prime}d\mbox{\boldmath$\eta$}^{\prime}}{\Delta^{2}\Delta_{p}^{2}}\times$
$\times\sum_{i=1}^{N/2}\sum_{j=1,j\neq
i}^{N/2}\delta_{\textbf{r}i}\delta_{\textbf{p}i}\delta_{\textbf{r}^{\prime}j}\delta_{\textbf{p}^{\prime}j}\equiv\langle\langle
1\rangle\rangle$ (23)

is the two-coordinate distribution function. The compact notation is also
introduced for the two-coordinate distribution function constructed of double
brackets $\langle\langle 1\rangle\rangle$. We use this notation below for
other two-coordinate distribution functions.

##  IV The selfconsistent field approximation in classic monopole kinetics for
the equation of evolution of scalar distribution function

In order to consider the selfconsistent field (the mean-field) approximation
we need to split the two-coordinate distribution function
$f_{2,ee}(\textbf{r},\textbf{r}^{\prime},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})$
into the product of two one-coordinate distribution functions
$f_{2,ee}(\textbf{r},\textbf{r}^{\prime},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})=f(\textbf{r},\textbf{p},t)\cdot
f(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})$. Hence, the term
containing interaction in equation (22) reappears in the following form

$\frac{q_{s}}{m_{s}}q_{s^{\prime}}\nabla_{\textbf{p}}f(\textbf{r},\textbf{p},t)\cdot\biggl{(}-\nabla_{\textbf{r}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(\textbf{r},\textbf{r}^{\prime})f(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})$
$-\frac{1}{c^{2}}\int d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\textbf{v}^{\prime}\partial_{t}G(\textbf{r},\textbf{r}^{\prime})f(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})$
$+\frac{1}{c^{2}}[\textbf{v}\times[\nabla\times\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\textbf{v}^{\prime}]]G(\textbf{r},\textbf{r}^{\prime})f(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})=0.$
(24)

It allows to introduce the scalar and vector potentials of the self-consistent
electromagnetic field

$\varphi_{int}(\textbf{r},t)=\int d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(\textbf{r},\textbf{r}^{\prime})f(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})$
(25)

and

$\textbf{A}_{int}(\textbf{r},t)=\frac{1}{c}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\textbf{v}^{\prime}G(\textbf{r},\textbf{r}^{\prime})f(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}).$
(26)

Some discussion on the selfconsistent field approximation for the kinetics
based on the deterministic microscopic motion of particles is gives in Ref.
Andreev 2204  for the nonrelativistic kinetics and hydrodynamics.

So, we find the well-known form of the Vlasov kinetic equation for the system
of relativistic particles Vlasov JETP 38

$\partial_{t}f+\textbf{v}\cdot\nabla
f+q_{s}\biggl{(}\textbf{E}+\frac{1}{c}\textbf{v}\times\textbf{B}\biggr{)}\cdot\frac{\partial
f}{\partial\textbf{p}}=0,$ (27)

where $\textbf{E}=\textbf{E}_{ext}+\textbf{E}_{int}$,
$\textbf{B}=\textbf{B}_{ext}+\textbf{B}_{int}$, the electric field is caused
by the distribution of charges in the coordinate space:
$\textbf{E}_{int}=-\nabla\varphi-\partial_{t}\textbf{A}_{int}/c$,
$\textbf{B}_{int}=\nabla\times\textbf{A}_{int}$. Equation (27) appears to be
coupled with the electromagnetic field:
$\nabla\times\textbf{E}_{int}=-\partial_{t}\textbf{B}_{int}/c$,
$\nabla\cdot\textbf{B}_{int}=0$,

$\nabla\times\textbf{B}_{int}=\partial_{t}\textbf{E}_{int}/c+(4\pi/c)\sum_{s}q_{s}\int\textbf{v}f_{s}(\textbf{r},\textbf{p},t)d\textbf{p},$
(28)

and

$\nabla\cdot\textbf{E}_{int}=4\pi\sum_{s}q_{s}\int
f_{s}(\textbf{r},\textbf{p},t)d\textbf{p}.$ (29)

##  V Multipole approximation in the relativistic kinetic equation for the
scalar distribution function (the Vlasov equation)

###  V.1 The interaction of particles with the external electromagnetic field

The multipole expansion of the electric field
$\textbf{E}_{ext}(\textbf{r}+\mbox{\boldmath$\xi$},t)$, the magnetic field
$\textbf{B}_{ext}(\textbf{r}+\mbox{\boldmath$\xi$},t)$, and the Green function
$G(\textbf{r}-\textbf{r}^{\prime}+\mbox{\boldmath$\xi$}-\mbox{\boldmath$\xi$}^{\prime},t-t^{\prime})$
leads to appearance of new distribution functions of different tensor ranks.
Let us introduce these functions before we find they appearance from the
expansion of equation (19) on $\xi$ and $\mbox{\boldmath$\xi$}^{\prime}$ up to
the second order on these vectors.

Firstly, we introduce the vector distribution function of the electric dipole
moment (divided by the charge $q_{s}$) or the displacement of particles
relatively the center of the $\Delta$-vicinity:

$d^{a}(\textbf{r},\textbf{p},t)=\langle\xi^{a}\rangle=\int_{\Delta,\Delta_{p}}\frac{d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}}{\Delta\Delta_{p}}\sum_{i=1}^{N}\xi^{a}\delta_{\textbf{r}i}\delta_{\textbf{p}i}.$
(30)

We also include in our analysis the distribution function of the electric
quadrupole moment (divided by the charge $q_{s}$)

$Q^{ab}(\textbf{r},\textbf{p},t)=\langle\xi^{a}\xi^{b}\rangle=\int_{\Delta,\Delta_{p}}\frac{d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}}{\Delta\Delta_{p}}\sum_{i=1}^{N}\xi^{a}\xi^{b}\delta_{\textbf{r}i}\delta_{\textbf{p}i}.$
(31)

Presented functions contain vector $\xi$ scanning the $\Delta$-vicinity.
Functions (30) and (31) appear as two examples of the infinite set of
distribution functions containing different degrees of vector $\xi$ (the
product of different number of projections to construct an element of tensor
of corresponding rank).

As it is demonstrated by equations (11)-(14) we can find the velocity of
particle under integral defining the distribution function. We can also find
the product of several projections of the velocity (of the same particle, in
order to get one-coordinate distribution function) or the product of the
projection of the velocity on the projections of vector $\xi$.

One of such distribution functions appearing in our derivation is the
distribution function of flux of the electric dipole moment (divided by the
charge $q_{s}$)

$J_{D}^{ab}(\textbf{r},\textbf{p},t)=\langle
v_{i}^{a}(t)\xi^{b}\rangle=\frac{1}{\Delta^{2}}\frac{1}{\Delta_{p}^{2}}\times$
$\times\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}d\mbox{\boldmath$\xi$}^{\prime}d\mbox{\boldmath$\eta$}^{\prime}\sum_{i,j=1,j\neq
i}^{N}v_{i}^{a}(t)\xi^{b}\delta_{\textbf{r}i}\delta_{\textbf{p}i}\delta_{\textbf{r}^{\prime}j}\delta_{\textbf{p}^{\prime}j}.$
(32)

It is more useful to extract the velocity corresponding to the center of
$\Delta_{\textbf{p}}$-vicinity in the momentum space, like we make it for
$\langle v_{i}^{a}(t)\rangle$ in equations (11)-(14). Therefore, function (32)
leads to the following distribution function:

$j_{D}^{ab}(\textbf{r},\textbf{p},t)=\langle\Delta
v_{i}^{a}(t)\xi^{b}\rangle=\frac{1}{\Delta^{2}}\frac{1}{\Delta_{p}^{2}}\times$
$\times\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}d\mbox{\boldmath$\xi$}^{\prime}d\mbox{\boldmath$\eta$}^{\prime}\sum_{i,j=1,j\neq
i}^{N}(\Delta
v_{i}^{a}(t))\xi^{b}\delta_{\textbf{r}i}\delta_{\textbf{p}i}\delta_{\textbf{r}^{\prime}j}\delta_{\textbf{p}^{\prime}j}.$
(33)

Let us introduce one more distribution function, which is the third rank
tensor. It is the distribution function of flux of the electric quadrupole
moment (divided by the charge $q_{s}$)

$J_{Q}^{abc}(\textbf{r},\textbf{p},t)=\langle
v_{i}^{a}(t)\xi^{b}\xi^{c}\rangle=\frac{1}{\Delta^{2}}\frac{1}{\Delta_{p}^{2}}\times$
$\times\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}d\mbox{\boldmath$\xi$}^{\prime}d\mbox{\boldmath$\eta$}^{\prime}\sum_{i,j=1,j\neq
i}^{N}v_{i}^{a}(t)\xi^{b}\xi^{c}\delta_{\textbf{r}i}\delta_{\textbf{p}i}\delta_{\textbf{r}^{\prime}j}\delta_{\textbf{p}^{\prime}j}.$
(34)

For this distribution function, we can also introduce the reduced flux of the
electric quadrupole moment on the velocities in the local frame $\Delta
v_{i}^{a}(t)$: $j_{Q}^{abc}(\textbf{r},\textbf{p},t)=\langle\Delta
v_{i}^{a}(t)\xi^{b}\xi^{c}\rangle$.

Functions (30)-(34) appear in the Vlasov-like kinetic equation for the scalar
distribution function. Firstly, these functions appear at the analysis of the
action of the external electromagnetic field on the system of charged
particles. Here, we need to consider the fourth term in equation (19). We
apply the expansion of the external electromagnetic field on the deviations of
coordinate $\xi$ from the center of the $\Delta$-vicinity of point r. This
expansion requires slow change of the fields E and B over the diameter
$\sqrt[3]{\Delta}$ of the $\Delta$-vicinity. Particularly, if we consider the
electromagnetic wave we assume that the wavelength $\lambda$ is larger than
the diameter $\sqrt[3]{\Delta}$: $\lambda\geq\sqrt[3]{\Delta}$.

For instance, if we consider diameter $\sqrt[3]{\Delta}\approx 0.1$ $\mu$m we
get estimation on the minimal value of the frequency of the electromagnetic
wave of order of $\omega_{max}\sim 10^{16}$ s-1. So, we can consider
electromagnetic waves with frequencies up to the frequencies of the visible
light, but the ultraviolet radiation is our of the range of applicability of
the model. If we deal with the dense medium we can choose the delta-vicinity
of the smaller value. So, the range of the ultraviolet radiation can be
included in our analysis.

During derivation of the kinetic equations we consider the expansion up to the
second order on $\xi^{a}$:
$\textbf{E}(\textbf{r}+\mbox{\boldmath$\xi$},t)=\textbf{E}(\textbf{r},t)$$+\xi_{a}\partial^{a}_{\textbf{r}}\textbf{E}(\textbf{r},t)$
$+(1/2)\xi_{a}\xi_{b}\partial^{a}_{\textbf{r}}\partial^{b}_{\textbf{r}}\textbf{E}(\textbf{r},t)+...$,
and
$\textbf{B}(\textbf{r}+\mbox{\boldmath$\xi$},t)=\textbf{B}(\textbf{r},t)$$+\xi_{a}\partial^{a}_{\textbf{r}}\textbf{B}(\textbf{r},t)$
$+(1/2)\xi_{a}\xi_{b}\partial^{a}_{\textbf{r}}\partial^{b}_{\textbf{r}}\textbf{B}(\textbf{r},t)+...$.
Hence, the fourth term in equation (19) can be rewritten as

$\Omega_{0}=\frac{q_{s}}{m_{s}}\langle(\textbf{E}_{ext}(\textbf{r}+\mbox{\boldmath$\xi$},t)$
$+\frac{1}{c}[\textbf{v}\times\textbf{B}_{ext}(\textbf{r}+\mbox{\boldmath$\xi$},t)]+\frac{1}{c}[\Delta\textbf{v}_{i}\times\textbf{B}_{ext}(\textbf{r}+\mbox{\boldmath$\xi$},t)])\nabla_{\textbf{p}}\rangle$
(35)

before the expansion. After described expansion we obtain

$\Omega_{0}=\frac{q_{s}}{m_{s}}\biggl{\\{}\textbf{E}_{ext}(\textbf{r},t)\nabla_{\textbf{p}}f(\textbf{r},\textbf{p},t)+(\partial^{a}_{\textbf{r}}\textbf{E}_{ext})\nabla_{\textbf{p}}d^{a}(\textbf{r},\textbf{p},t)$
$+\frac{1}{2}(\partial^{a}_{\textbf{r}}\partial^{b}_{\textbf{r}}\textbf{E}_{ext})\nabla_{\textbf{p}}Q^{ab}(\textbf{r},\textbf{p},t)+\frac{1}{c}[\textbf{v}\times\textbf{B}_{ext}(\textbf{r},t)]\nabla_{\textbf{p}}f(\textbf{r},\textbf{p},t)$
$+\frac{1}{c}[\textbf{v}\times(\partial^{a}_{\textbf{r}}\textbf{B}_{ext})]\nabla_{\textbf{p}}d^{a}(\textbf{r},\textbf{p},t)$
$+\frac{1}{c}[\textbf{v}\times(\partial^{a}_{\textbf{r}}\partial^{b}_{\textbf{r}}\textbf{B}_{ext})]\nabla_{\textbf{p}}Q^{ab}(\textbf{r},\textbf{p},t)$
$+\frac{1}{c}\varepsilon^{abc}\biggl{[}B^{b}_{ext}\nabla_{\textbf{p}}^{c}F^{a}(\textbf{r},\textbf{p},t)+(\partial^{d}_{\textbf{r}}B^{b}_{ext})\nabla_{\textbf{p}}^{c}j^{ad}(\textbf{r},\textbf{p},t)$
$+(\partial^{d}_{\textbf{r}}\partial^{f}_{\textbf{r}}B^{b}_{ext})\nabla_{\textbf{p}}^{c}j^{adf}(\textbf{r},\textbf{p},t)\biggr{]}\biggr{\\}}.$
(36)

###  V.2 The interparticle interaction

Let us introduce function $\Omega$ which presents a term in the kinetic
equation describing the interparticle interaction (the last term in equation
(19)). We split this term on three parts
$\Omega=\Omega_{1}+\Omega_{2}+\Omega_{3}$ which are
$\Omega_{1}=-(q_{s}/m_{s})\langle\nabla_{i}\varphi_{i,int}\cdot\nabla_{p,i}\rangle$,
$\Omega_{2}=-(q_{s}/m_{s}c)\langle\partial_{t}\textbf{A}_{i,int}\cdot\nabla_{p,i}\rangle$,
$\Omega_{3}=(q_{s}/m_{s}c)\varepsilon^{abc}\varepsilon^{cfg}\langle
v_{i}^{b}(t)\partial_{i}^{f}A_{i,int}^{g}\cdot\nabla_{p,i}\rangle$, where we
use short notation for the average on the physically infinitesimal volume,
like

$\langle\nabla_{i}\varphi_{i,int}\cdot\nabla_{p,i}\rangle=\int_{\Delta_{\textbf{r}}\Delta_{\textbf{p}}}\frac{d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}}{\Delta_{\textbf{r}}\Delta_{\textbf{p}}}\sum_{i}\nabla_{i}\varphi_{i,int}\delta_{\textbf{r}i}\cdot\nabla_{p,i}\delta_{\textbf{p}i}.$
(37)

Here, functions $\varphi_{i,int}$ and $\textbf{A}_{i,int}$ are given by
equations (7) and (8).

Let us consider the multipole expansion for each of function $\Omega_{j}$
(j=1,2,3). We start this presentation with $\Omega_{1}$ which can be
represented in the following form

$\Omega_{1}=-\int\frac{d\textbf{r}^{\prime}d\textbf{p}^{\prime}}{(\Delta_{\textbf{r}}\Delta_{\textbf{p}})^{2}}\int
dt^{\prime}\int_{\Delta_{\textbf{r}}\Delta_{\textbf{p}}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\int_{\Delta_{\textbf{r}^{\prime}}\Delta_{\textbf{p}^{\prime}}}d\mbox{\boldmath$\xi$}^{\prime}d\mbox{\boldmath$\eta$}^{\prime}\times$
$\times\sum_{i,j,j\neq
i}\frac{q_{i}q_{j}}{m_{s}}\nabla_{\textbf{r}}\biggl{(}G(t,t^{\prime},\textbf{r}+\mbox{\boldmath$\xi$},\textbf{r}^{\prime}+\mbox{\boldmath$\xi$}^{\prime})\biggr{)}\delta_{\textbf{r}i}\cdot\nabla_{p,i}\delta_{\textbf{p}i}\delta_{\textbf{r}^{\prime}j}\delta_{\textbf{p}^{\prime}j}.$
(38)

We make expansion of the Green functions on $\xi$ and
$\mbox{\boldmath$\xi$}^{\prime}$. We also expand the velocity of particle
$v_{i}^{b}$ on the deviation from the average velocity $\eta$. For the Green
function
$G(t,t^{\prime},\textbf{r}+\mbox{\boldmath$\xi$},\textbf{r}^{\prime}+\mbox{\boldmath$\xi$}^{\prime})$
we have $G=G_{0}+\xi^{b}G_{1,b}+(1/2)\xi^{b}\xi^{c}G_{2,bc}$, where
$G_{0}=\delta(t-t^{\prime}-\mid\textbf{r}-\textbf{r}^{\prime}\mid/c)/\mid\textbf{r}-\textbf{r}^{\prime}\mid$,
$G_{1,b}=\partial_{\textbf{r}}^{b}G_{0}$, and
$G_{2,bc}=\partial_{\textbf{r}}^{b}\partial_{\textbf{r}}^{c}G_{0}$.

Function $\Omega_{1}$ can be represented in terms of two-coordinate
distribution functions:

$\Omega_{1}=-\frac{q_{s}q_{s^{\prime}}}{m_{s}}\nabla_{\textbf{p}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\biggl{[}(\nabla_{\textbf{r}}G_{0})f_{2,ss^{\prime}}(\textbf{r},\textbf{r}^{\prime},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})$
$+(\nabla_{\textbf{r}}G_{1}^{\beta})\biggl{(}d_{2,ss^{\prime}}^{b}(\textbf{r},\textbf{r}^{\prime},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})-d_{2,ss^{\prime}}^{b}(\textbf{r}^{\prime},\textbf{r},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})\biggr{)}$
$+\frac{1}{2}(\nabla_{\textbf{r}}G_{2}^{bc})\biggl{(}q_{2,ss^{\prime}}^{bc}(\textbf{r},\textbf{r}^{\prime},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})+q_{2,ss^{\prime}}^{bc}(\textbf{r}^{\prime},\textbf{r},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})$
$-D_{2,ss^{\prime}}^{bc}(\textbf{r},\textbf{r}^{\prime},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})-D_{2,ss^{\prime}}^{bc}(\textbf{r}^{\prime},\textbf{r},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})\biggr{)}\biggr{]},$
(39)

where we use the following notations for the two-coordinate distribution
functions

$d_{2,ss^{\prime}}^{b}(\textbf{r},\textbf{r}^{\prime},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})\equiv\langle\langle\xi^{b}\rangle\rangle\rightarrow
d^{b}_{s}(\textbf{r},\textbf{p},t)\cdot
f_{s^{\prime}}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (40)
$q_{2,ss^{\prime}}^{bc}(\textbf{r},\textbf{r}^{\prime},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})\equiv\langle\langle\xi^{b}\xi^{c}\rangle\rangle\rightarrow
q_{s}^{bc}(\textbf{r},\textbf{p},t)\cdot
f_{s^{\prime}}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (41)

and

$D_{2,ss^{\prime}}^{bc}(\textbf{r},\textbf{r}^{\prime},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})\equiv\langle\langle\xi^{b}\xi^{\prime
c}\rangle\rangle\rightarrow d^{b}_{s}(\textbf{r},\textbf{p},t)\cdot
d_{s^{\prime}}^{c}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}).$ (42)

Potential part of the electric force $q_{s}\textbf{E}$ of the interparticle
interaction is presented in terms of the two-coordinate distribution function.
Next, we make same representation for the part of the electric force
$q_{s}\textbf{E}$ expressed via the vector potential. Therefore, let us
demonstrate the explicit form of $\Omega_{2}$:

$\Omega_{2}=-\int\frac{d\textbf{r}^{\prime}d\textbf{p}^{\prime}}{(\Delta_{\textbf{r}}\Delta_{\textbf{p}})^{2}}\int
dt^{\prime}\int_{\Delta_{\textbf{r}}\Delta_{\textbf{p}}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\int_{\Delta_{\textbf{r}^{\prime}}\Delta_{\textbf{p}^{\prime}}}d\mbox{\boldmath$\xi$}^{\prime}d\mbox{\boldmath$\eta$}^{\prime}\times$
$\times\sum_{i,j,j\neq
i}\frac{q_{i}q_{j}}{m_{s}c^{2}}\partial_{t}\biggl{(}\textbf{v}_{j}(t^{\prime})G(t,t^{\prime},\textbf{r}+\mbox{\boldmath$\xi$},\textbf{r}^{\prime}+\mbox{\boldmath$\xi$}^{\prime})\biggr{)}\delta_{\textbf{r}i}\cdot\nabla_{p,i}\delta_{\textbf{p}i}\delta_{\textbf{r}^{\prime}j}\delta_{\textbf{p}^{\prime}j}.$
(43)

Function $\Omega_{2}$ can be approximately represented in terms of two-
coordinate distribution functions after expansion on $\xi$ and
$\mbox{\boldmath$\xi$}^{\prime}$:

$\Omega_{2}=-\frac{1}{c^{2}}\frac{q_{s}q_{s^{\prime}}}{m_{s}}\partial_{\textbf{p}}^{a}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\biggl{[}(\partial_{t}G_{0})\langle\langle
v_{i}^{a}(t)\rangle\rangle$
$+(\partial_{t}G_{1}^{b})\biggl{(}\langle\langle\xi^{b}v_{j}^{a}(t^{\prime})\rangle\rangle-\langle\langle\xi^{\prime
b}v_{j}^{a}(t^{\prime})\rangle\rangle\biggr{)}+\frac{1}{2}(\partial_{t}G_{2}^{bc})\times$
$\times\biggl{(}\langle\langle\xi^{b}\xi^{c}v_{j}^{a}(t^{\prime})\rangle\rangle+\langle\langle\xi^{\prime
b}\xi^{\prime
c}v_{j}^{a}(t^{\prime})\rangle\rangle-2\langle\langle\xi^{b}\xi^{\prime
c}v_{j}^{a}(t^{\prime})\rangle\rangle\biggr{)}\biggr{]},$ (44)

where we use the following notations for the two-coordinate distribution
functions

$\langle\langle v_{i}^{a}(t)\rangle\rangle\rightarrow
j^{a}_{s}(\textbf{r},\textbf{p},t)\cdot
f_{s^{\prime}}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (45)
$\langle\langle v_{j}^{a}(t^{\prime})\rangle\rangle\rightarrow
f_{s}(\textbf{r},\textbf{p},t)\cdot
j_{s^{\prime}}^{a}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (46)
$\langle\langle\xi^{b}v_{i}^{a}(t)\rangle\rangle\rightarrow
J_{D,s}^{ab}(\textbf{r},\textbf{p},t)\cdot
f_{s^{\prime}}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (47)
$\langle\langle\xi^{b}v_{j}^{a}(t^{\prime})\rangle\rangle\rightarrow
d_{s}^{b}(\textbf{r},\textbf{p},t)\cdot
j^{a}_{s^{\prime}}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (48)
$\langle\langle\xi^{\prime b}v_{j}^{a}(t^{\prime})\rangle\rangle\rightarrow
f_{s}(\textbf{r},\textbf{p},t)\cdot
J_{D,s^{\prime}}^{ab}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$
(49)
$\langle\langle\xi^{b}\xi^{c}v_{j}^{a}(t^{\prime})\rangle\rangle\rightarrow
q^{bc}_{s}(\textbf{r},\textbf{p},t)\cdot
j_{s^{\prime}}^{a}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (50)
$\langle\langle\xi^{\prime b}\xi^{\prime
c}v_{j}^{a}(t^{\prime})\rangle\rangle\rightarrow
f_{s}(\textbf{r},\textbf{p},t)\cdot
J_{Q,s^{\prime}}^{abc}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$
(51)

and

$\langle\langle\xi^{b}\xi^{\prime
c}v_{j}^{a}(t^{\prime})\rangle\rangle\rightarrow
d_{s}^{b}(\textbf{r},\textbf{p},t)\cdot
J_{D,s^{\prime}}^{ac}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}).$
(52)

Finally, we present the explicit form of the magnetic part of the
interparticle interaction:

$\Omega_{3}=\frac{1}{c^{2}}\varepsilon^{abc}\varepsilon^{cfg}\frac{1}{(\Delta_{\textbf{r}}\Delta_{\textbf{p}})^{2}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int dt^{\prime}$
$\int_{\Delta_{\textbf{r}}\Delta_{\textbf{p}}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\int_{\Delta_{\textbf{r}^{\prime}}\Delta_{\textbf{p}^{\prime}}}d\mbox{\boldmath$\xi$}^{\prime}d\mbox{\boldmath$\eta$}^{\prime}\sum_{i,j,j\neq
i}\frac{q_{i}q_{j}}{m_{s}}v_{i}^{b}(t)v_{j}^{g}(t^{\prime})\times$
$\times\partial_{\textbf{r}}^{f}\biggl{(}G(\textbf{r}+\mbox{\boldmath$\xi$},\textbf{r}^{\prime}+\mbox{\boldmath$\xi$}^{\prime})\biggr{)}\delta_{\textbf{r}i}\cdot\partial_{p,i}^{a}\delta_{\textbf{p}i}\delta_{\textbf{r}^{\prime}j}\delta_{\textbf{p}^{\prime}j}$
(53)

Here, we make expansion on $\xi$ and $\mbox{\boldmath$\xi$}^{\prime}$ and make
interpretation of found terms via the two-coordinate distribution functions:

$\Omega_{3}=\frac{1}{c^{2}}\frac{q_{s}q_{s^{\prime}}}{m_{s}}\varepsilon^{abc}\varepsilon^{cfg}\partial_{\textbf{p}}^{a}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\biggl{[}(\partial_{\textbf{r}}^{f}G_{0})\langle\langle
v_{i}^{b}(t)v_{j}^{g}(t^{\prime})\rangle\rangle$
$+(\partial_{\textbf{r}}^{f}G_{1}^{k})\biggl{(}\langle\langle\xi^{k}v_{i}^{b}(t)v_{j}^{g}(t^{\prime})\rangle\rangle-\langle\langle\xi^{\prime
k}v_{i}^{b}(t)v_{j}^{g}(t^{\prime})\rangle\rangle\biggr{)}$
$+\frac{1}{2}(\partial_{\textbf{r}}^{f}G_{2}^{kl})\biggl{(}\langle\langle\xi^{k}\xi^{l}v_{i}^{b}(t)v_{j}^{g}(t^{\prime})\rangle\rangle$
$+\langle\langle\xi^{\prime k}\xi^{\prime
l}v_{i}^{b}(t)v_{j}^{g}(t^{\prime})\rangle\rangle-2\langle\langle\xi^{k}\xi^{\prime
l}v_{i}^{b}(t)v_{j}^{g}(t^{\prime})\rangle\rangle\biggr{)}\biggr{]},$ (54)

where we use no specific notations for the two-coordinate distribution
functions since we have two many two-coordinate distribution functions. We
have no reason to use special symbol for each of them. We also have following
limits for the self-consistent field approximation:

$\langle\langle v_{i}^{b}(t)v_{j}^{g}(t^{\prime})\rangle\rangle\rightarrow
j^{b}(\textbf{r},\textbf{p},t)\cdot
j^{g}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (55)
$\langle\langle\xi^{a}v_{i}^{b}(t)v_{j}^{g}(t^{\prime})\rangle\rangle\rightarrow
J_{D}^{ba}(\textbf{r},\textbf{p},t)\cdot
j^{g}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (56)
$\langle\langle\xi^{\prime
a}v_{i}^{b}(t)v_{j}^{g}(t^{\prime})\rangle\rangle\rightarrow
j^{b}(\textbf{r},\textbf{p},t)\cdot
J_{D}^{ga}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (57)
$\langle\langle\xi^{s}\xi^{l}v_{i}^{b}(t)v_{j}^{g}(t^{\prime})\rangle\rangle\rightarrow
J^{bsl}_{Q}(\textbf{r},\textbf{p},t)\cdot
j^{g}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (58)
$\langle\langle\xi^{\prime s}\xi^{\prime
l}v_{i}^{b}(t)v_{j}^{g}(t^{\prime})\rangle\rangle\rightarrow
j^{b}(\textbf{r},\textbf{p},t)\cdot
J^{gsl}_{Q}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (59)

and

$\langle\langle\xi^{s}\xi^{\prime
l}v_{i}^{b}(t)v_{j}^{g}(t^{\prime})\rangle\rangle\rightarrow
J^{bs}_{D}(\textbf{r},\textbf{p},t)\cdot
J_{D}^{gl}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}).$ (60)

To conclude this section, we point out that we made the generalization of
equation (22) up to the account of the qudrupole moment distribution function
in the multipole expansion while equation (22) includes the distribution of
charger (and currents as well) in the monopole approximation. Presented
generalization is demonstrated by four terms labeled as $\Omega_{0}$,
$\Omega_{1}$, $\Omega_{2}$, and $\Omega_{3}$.

##  VI Selfconsistent field approximation in the multipole approximation for
the scalar distribution function evolution equation

We consider the dynamics of charged particles. Therefore, the mean-field or
the self-consistent field approximation gives the major contribution in the
dynamics of system. Hence, we consider all two-coordinate distribution
functions as the product of the corresponding one-coordinate distribution
functions. Therefore, we combine equations (36), (39), (44), (54) represented
in the selfconsistent field approximation and find the following equation for
the scalar distribution function:

$\partial_{t}f_{s}+\textbf{v}\cdot\nabla
f_{s}+\nabla\cdot\textbf{F}_{s}+q_{s}\biggl{(}\textbf{E}+\frac{1}{c}\textbf{v}\times\textbf{B}\biggr{)}\cdot\frac{\partial
f_{s}}{\partial\textbf{p}}$
$+q_{s}\biggl{(}\partial^{b}\textbf{E}+\frac{1}{c}\textbf{v}\times\partial^{b}\textbf{B}\biggr{)}\cdot\frac{\partial
d_{s}^{b}}{\partial\textbf{p}}+q_{s}\biggl{(}\partial^{b}\partial^{c}\textbf{E}$
$+\frac{q_{s}}{c}\textbf{v}\times\partial^{b}\partial^{c}\textbf{B}\biggr{)}\cdot\frac{\partial
Q_{s}}{\partial\textbf{p}}+\frac{q_{s}}{c}\varepsilon^{abc}\biggl{(}B^{c}\cdot\partial_{a,\textbf{p}}F_{s}^{b}$
$+\partial^{d}B^{c}\cdot\partial_{a,\textbf{p}}j_{D,s}^{bd}+\partial^{d}\partial^{f}B^{c}\cdot\partial_{a,\textbf{p}}j_{Q,s}^{bdf}\biggr{)}=0,$
(61)

where the velocity is extracted from functions $j_{s}^{a}$, $j_{D,s}^{ab}$,
and $j_{Q,s}^{abc}$. We also have
$\textbf{E}=\textbf{E}_{ext}+\textbf{E}_{int}$, and
$\textbf{B}=\textbf{B}_{ext}+\textbf{B}_{int}$.

The integral terms presenting the sources of the electromagnetic field are
written as $\textbf{E}_{int}$ and $\textbf{B}_{int}$, while these functions
obey the Maxwell equations in accordance with the explicitly found integral
expressions for the electromagnetic field
$\\{\textbf{E}_{int},\textbf{B}_{int}\\}$ via the one-coordinate distribution
functions:

$\begin{array}[]{cc}\nabla\cdot\textbf{B}_{int}=0,&\nabla\times\textbf{E}_{int}=-\frac{1}{c}\partial_{t}\textbf{B}_{int},\end{array}$
(62)
$(\nabla\times\textbf{B}_{int})^{a}=\frac{1}{c}\partial_{t}E_{int}^{a}+\frac{4\pi}{c}\sum_{s}q_{s}\biggl{(}\int
j_{s}^{a}(\textbf{r},\textbf{p},t)d\textbf{p}$ $+\partial^{b}\int
J_{s,D}^{ab}(\textbf{r},\textbf{p},t)d\textbf{p}+\frac{1}{2}\partial^{b}\partial^{c}\int
J_{s,Q}^{abc}(\textbf{r},\textbf{p},t)d\textbf{p}+...\biggr{)},$ (63)
$\nabla\cdot\textbf{E}_{int}=4\pi\sum_{s}q_{s}\biggl{(}\int
f_{s}(\textbf{r},\textbf{p},t)d\textbf{p}$ $+\partial^{b}\int
d_{s}^{b}(\textbf{r},\textbf{p},t)d\textbf{p}+\frac{1}{2}\partial^{b}\partial^{c}\int
Q_{s}^{bc}(\textbf{r},\textbf{p},t)d\textbf{p}+...\biggr{)}.$ (64)

The selfconsistent field approximation shows that all introduced distribution
functions appear as the sources of the electromagnetic field in the Maxwell
equations.

##  VII Equation for the evolution of the dipole moment distribution function

Let us repeat the definition of the distribution function of dipole moment

$d^{a}(\textbf{r},\textbf{p},t)=\frac{1}{\Delta}\frac{1}{\Delta_{p}}\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\sum_{i=1}^{N}\xi^{a}\delta_{\textbf{r}i}\delta_{\textbf{p}i}$
(65)

If we want to consider the evolution of the distribution function of dipole
moment we can consider the time derivative of function $\langle\xi^{a}\rangle$
(65). We can consider alternative function $\langle
r^{a}(t)\rangle$$=r^{a}f(\textbf{r},\textbf{p},t)+\langle\xi^{a}\rangle$.

The time derivative of function (65) has the following structure

$\partial_{t}d^{a}+\partial^{b}_{\textbf{r}}\langle\xi^{a}v_{i,b}(t)\rangle+\partial^{b}_{\textbf{p}}\langle\xi^{a}\dot{p}_{i,b}(t)\rangle=0.$
(66)

The second term in this kinetic equation can be represented via functions (32)
or (33). While the last term in this equation requires longer discussion
similar to the analysis of interaction in the Vlasov equation presented above.

Substituting the time derivative of the momentum $\dot{p}_{i}^{b}$ into
equation (66) (equation of motion of each particle is given by equation
(5))-(8) we obtain

$\partial_{t}d^{a}+\partial^{b}_{\textbf{r}}J_{D}^{ba}+q_{s}[\langle\xi^{a}E^{b}_{ext}(\textbf{r}+\mbox{\boldmath$\xi$},t)\partial_{b,\textbf{p}}\rangle+\frac{1}{c}\varepsilon^{bcd}\langle\xi^{a}v_{i}^{c}(t)B^{d}_{ext}(\textbf{r}+\mbox{\boldmath$\xi$},t)\partial_{b,\textbf{p}}\rangle]$
$-q_{s}q_{s^{\prime}}\int d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\langle\langle\xi^{a}(\partial_{\textbf{r}}^{b}G(t,t^{\prime},\textbf{r}+\mbox{\boldmath$\xi$},\textbf{r}^{\prime}+\mbox{\boldmath$\xi$}^{\prime}))\partial_{b,\textbf{p}}\rangle\rangle-\frac{q_{s}q_{s^{\prime}}}{c^{2}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\langle\langle\xi^{a}v_{j}^{b}(t^{\prime})(\partial_{t}G(t,t^{\prime},\textbf{r}+\mbox{\boldmath$\xi$},\textbf{r}^{\prime}+\mbox{\boldmath$\xi$}^{\prime}))\partial_{b,\textbf{p}}\rangle\rangle$
$+\frac{q_{s}q_{s^{\prime}}}{c^{2}}\varepsilon^{bcd}\varepsilon^{dfg}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\langle\langle\xi^{a}v_{i}^{c}(t)v_{j}^{g}(t^{\prime})(\partial_{\textbf{r}}^{f}G(t,t^{\prime},\textbf{r}+\mbox{\boldmath$\xi$},\textbf{r}^{\prime}+\mbox{\boldmath$\xi$}^{\prime}))\partial_{b,\textbf{p}}\rangle\rangle=0.$
(67)

Next, we make expansion on $\xi$ and $\mbox{\boldmath$\xi$}^{\prime}$ in order
to get the multipole expansion. We keep terms up to the electric quadrupole
moment and find:

$\partial_{t}d^{a}+\partial^{b}_{\textbf{r}}J_{D}^{ba}+q_{s}\biggl{[}E^{b}_{ext}(\textbf{r},t)\partial_{b,\textbf{p}}\langle\xi^{a}\rangle+(\partial^{c}_{\textbf{r}}E^{b}_{ext}(\textbf{r},t))\partial_{b,\textbf{p}}\langle\xi^{a}\xi^{c}\rangle+\frac{1}{c}\varepsilon^{bcd}B^{d}_{ext}(\textbf{r},t)\partial_{b,\textbf{p}}\langle\xi^{a}v_{i}^{c}(t)\rangle$
$+\frac{1}{c}\varepsilon^{bcd}(\partial^{f}_{\textbf{r}}B^{d}_{ext}(\textbf{r},t))\partial_{b,\textbf{p}}\langle\xi^{a}\xi^{f}v_{i}^{c}(t)\rangle\biggr{]}-q_{s}q_{s^{\prime}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}(\partial_{\textbf{r}}^{b}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime}))\partial_{b,\textbf{p}}\langle\langle\xi^{a}\rangle\rangle$
$-q_{s}q_{s^{\prime}}\int d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}(\partial_{\textbf{r}}^{b}\partial_{\textbf{r}}^{c}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime}))\partial_{b,\textbf{p}}\langle\langle\xi^{a}\xi^{c}\rangle\rangle+q_{s}q_{s^{\prime}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}(\partial_{\textbf{r}}^{b}\partial_{\textbf{r}}^{c}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime}))\partial_{b,\textbf{p}}\langle\langle\xi^{a}\xi^{\prime
c}\rangle\rangle$ $-\frac{q_{s}q_{s^{\prime}}}{c^{2}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}(\partial_{t}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime}))\partial_{b,\textbf{p}}\langle\langle\xi^{a}v_{j}^{b}(t^{\prime})\rangle\rangle-\frac{q_{s}q_{s^{\prime}}}{c^{2}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}(\partial_{t}\partial_{\textbf{r}}^{c}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime}))\partial_{b,\textbf{p}}\langle\langle\xi^{a}\xi^{c}v_{j}^{b}(t^{\prime})\rangle\rangle$
$+\frac{q_{s}q_{s^{\prime}}}{c^{2}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}(\partial_{t}\partial_{\textbf{r}}^{c}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime}))\partial_{b,\textbf{p}}\langle\langle\xi^{a}\xi^{\prime
c}v_{j}^{b}(t^{\prime})\rangle\rangle+\frac{q_{s}q_{s^{\prime}}}{c^{2}}\varepsilon^{bcd}\varepsilon^{dfg}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}(\partial_{\textbf{r}}^{f}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime}))\times$
$\times\partial_{b,\textbf{p}}\langle\langle\xi^{a}v_{i}^{c}(t)v_{j}^{g}(t^{\prime})\rangle\rangle+\frac{q_{s}q_{s^{\prime}}}{c^{2}}\varepsilon^{bcd}\varepsilon^{dfg}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}(\partial_{\textbf{r}}^{f}\partial_{\textbf{r}}^{l}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime}))\partial_{b,\textbf{p}}\langle\langle\xi^{a}\xi^{l}v_{i}^{c}(t)v_{j}^{g}(t^{\prime})\rangle\rangle$
$-\frac{q_{s}q_{s^{\prime}}}{c^{2}}\varepsilon^{bcd}\varepsilon^{dfg}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}(\partial_{\textbf{r}}^{f}\partial_{\textbf{r}}^{l}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime}))\partial_{b,\textbf{p}}\langle\langle\xi^{a}\xi^{\prime
l}v_{i}^{c}(t)v_{j}^{g}(t^{\prime})\rangle\rangle=0$ (68)

The two-coordinate distribution functions are presented as the double brackets
of the corresponding values with no usage of specific notations since part of
them are introduced above, while part of them get no specific notation. Here,
we make transition to the self-consistent field approximation. To this end, we
split the two-coordinate distribution functions on the product of one-
coordinate distribution functions (we do it in the same way as it is done
above for the equation of evolution of the scalar distribution function)

$\partial_{t}d^{a}+\partial^{b}_{\textbf{r}}J_{D}^{ba}+q_{s}(\partial_{b,\textbf{p}}d^{a})\biggl{[}E^{b}_{ext}(\textbf{r},t)-q_{s}\biggl{(}\partial^{b}_{\textbf{r}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})f(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})$
$+\frac{1}{c^{2}}\int d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\partial_{t}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})\langle
v_{j}^{b}(t^{\prime})\rangle(\textbf{r}^{\prime},\textbf{p}^{\prime})\biggr{)}+q_{s}\partial^{c}_{\textbf{r}}\biggl{(}\partial^{b}_{\textbf{r}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})d^{c}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})$
$+\frac{1}{c^{2}}\int d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\partial_{t}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})J_{D}^{bc}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})\biggr{)}\biggr{]}+q_{s}(\partial_{b,\textbf{p}}q^{ac})\biggl{[}\partial_{c,\textbf{r}}E^{b}_{ext}(\textbf{r},t)$
$-q_{s}\partial_{c,\textbf{r}}\biggl{(}\partial^{b}_{\textbf{r}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})f(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})+\frac{1}{c^{2}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\partial_{t}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})\langle
v_{j}^{b}(t^{\prime})\rangle(\textbf{r}^{\prime},\textbf{p}^{\prime})\biggr{)}\biggr{]}+\frac{q_{s}}{c}\varepsilon^{bcd}(\partial_{b,\textbf{p}}J_{D}^{ca})\times$
$\times\biggl{[}B_{ext}^{d}+\frac{q_{s^{\prime}}}{c}\varepsilon^{dfg}\partial_{f,\textbf{r}}\biggl{(}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})\langle
v_{j}^{g}(t^{\prime})\rangle(\textbf{r}^{\prime},\textbf{p}^{\prime})-\partial_{l,\textbf{r}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})J_{D}^{gl}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})\biggr{)}\biggr{]}$
$+\frac{q_{s}}{c}\varepsilon^{bcd}(\partial_{b,\textbf{p}}J_{Q}^{caf})\biggl{[}\partial_{f,\textbf{r}}B_{ext}^{d}+\frac{q_{s^{\prime}}}{c}\varepsilon^{dfg}\partial_{f,\textbf{r}}\partial_{l,\textbf{r}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})\langle
v_{j}^{g}(t^{\prime})\rangle(\textbf{r}^{\prime},\textbf{p}^{\prime})\biggr{]}=0.$
(69)

On this stage we can introduce the averaged scalar and vector potentials of
the electromagnetic field and corresponding representation of the integral
terms within $\textbf{E}_{int}$ and $\textbf{B}_{int}$. This electromagnetic
field obeys the Maxwell equations (62), (63), (64).

Finally, equation (69) is represented in terms of $\textbf{E}_{int}$ and
$\textbf{B}_{int}$:

$\partial_{t}d^{a}+\partial_{\textbf{r}}^{b}J_{D}^{ba}+q_{s}(\partial_{\textbf{p}}^{b}d^{a})E^{b}+q_{s}(\partial_{\textbf{p}}^{b}q^{ac})\partial_{\textbf{r}}^{c}E^{b}$
$+\frac{q_{s}}{c}\varepsilon^{bcd}(\partial_{\textbf{p}}^{b}J^{ca}_{D})B^{d}+\frac{q_{s}}{c}\varepsilon^{bcd}(\partial_{\textbf{p}}^{b}J^{caf}_{Q})\partial_{\textbf{r}}^{f}B^{d}=0.$
(70)

##  VIII Equation for the distribution function of velocity

Let us repeat the definition of the considering vector distribution function
(11)

$\textbf{j}(\textbf{r},\textbf{p},t)\equiv\langle\textbf{v}_{i}(t)\rangle\equiv\frac{1}{\Delta}\frac{1}{\Delta_{p}}\int_{\Delta,\Delta_{p}}d\mbox{\boldmath$\xi$}d\mbox{\boldmath$\eta$}\sum_{i=1}^{N/2}\textbf{v}_{i}(t)\delta_{\textbf{r}i}\delta_{\textbf{p}i},$
(71)

If we consider the evolution of function (71) we would have four terms in the
initial form of the kinetic equation

$\partial_{t}\textbf{j}(\textbf{r},\textbf{p},t)=\partial_{t}\langle\textbf{v}_{i}(t)\rangle$
$=\langle\dot{\textbf{v}}_{i}(t)\rangle-\partial_{r}^{b}\langle\textbf{v}_{i}(t)v_{i}^{b}(t)\rangle-\partial_{p}^{b}\langle\textbf{v}_{i}(t)\dot{p}_{i}^{b}(t)\rangle.$
(72)

However, we can represent the velocity of each particle $\textbf{v}_{i}(t)$
via parameters p and $\eta$ which do not depend on time (see equation (12))

$\partial_{t}\textbf{j}(\textbf{r},\textbf{p},t)=\partial_{t}\langle\textbf{v}+\Delta\textbf{v}\rangle$
$=-\partial_{r}^{b}\langle\textbf{v}_{i}(t)v_{i}^{b}(t)\rangle-\partial_{p}^{b}\langle\textbf{v}_{i}(t)\dot{p}_{i}^{b}(t)\rangle.$
(73)

Equation (72) can be simplified to equation (72) using the kinetic equation
for the scalar distribution function $f(\textbf{r},\textbf{p},t)$.

We consider equation (73) in more details

$\partial_{t}j^{a}(\textbf{r},\textbf{p},t)+\partial_{r}^{b}\langle
v_{i}^{a}(t)v_{i}^{b}(t)\rangle$ $+q_{s}\partial_{p}^{b}\langle
v_{i}^{a}E_{i}^{b}(\textbf{r}+\mbox{\boldmath$\xi$})+\varepsilon^{bcd}v_{i}^{a}v_{i}^{c}B_{i}^{d}(\textbf{r}+\mbox{\boldmath$\xi$})\rangle$
$-q_{s}q_{s^{\prime}}\int d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\langle\langle
v_{i}^{a}(t)\partial_{r}^{b}G(t,t^{\prime},\textbf{r}+\mbox{\boldmath$\xi$},\textbf{r}^{\prime}+\mbox{\boldmath$\xi$}^{\prime})\partial_{p}^{b}\rangle\rangle$
$-\frac{q_{s}q_{s^{\prime}}}{c^{2}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int dt^{\prime}\langle\langle
v_{i}^{a}(t)v_{j}^{b}(t^{\prime})\partial_{t}G(t,t^{\prime},\textbf{r}+\mbox{\boldmath$\xi$},\textbf{r}^{\prime}+\mbox{\boldmath$\xi$}^{\prime})\partial_{p}^{b}\rangle\rangle$
$+\frac{q_{s}q_{s^{\prime}}}{c^{2}}\varepsilon^{bcd}\varepsilon^{dfg}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int dt^{\prime}\times$
$\times\langle\langle
v_{i}^{a}(t)v_{i}^{c}(t)v_{j}^{g}(t^{\prime})\partial_{r}^{f}G(t,t^{\prime},\textbf{r}+\mbox{\boldmath$\xi$},\textbf{r}^{\prime}+\mbox{\boldmath$\xi$}^{\prime})\partial_{p}^{b}\rangle\rangle=0.$
(74)

Next, we represent equation (74) via the multipole expansion

$\partial_{t}j^{a}(\textbf{r},\textbf{p},t)+\partial_{r}^{b}\langle
v_{i}^{a}(t)v_{i}^{b}(t)\rangle+q_{s}\biggl{(}E^{b}(\textbf{r},t)\partial_{p}^{b}\langle
v_{i}^{a}\rangle$ $+(\partial_{r}^{c}E^{b})\partial_{p}^{b}\langle
v_{i}^{a}\xi^{c}\rangle+\varepsilon^{bcd}B^{d}\partial_{p}^{b}\langle
v_{i}^{a}v_{i}^{c}\rangle+\varepsilon^{bcd}(\partial^{f}B^{d})\partial_{p}^{b}\langle
v_{i}^{a}v_{i}^{c}\xi^{f}\rangle\biggr{)}$ $-q_{s}q_{s^{\prime}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\partial_{r}^{b}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})\partial_{p}^{b}\langle\langle
v_{i}^{a}(t)\rangle\rangle$ $-q_{s}q_{s^{\prime}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\partial_{r}^{b}\partial_{r}^{c}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})\partial_{p}^{b}\langle\langle\xi^{c}v_{i}^{a}(t)\rangle\rangle$
$+q_{s}q_{s^{\prime}}\int d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\partial_{r}^{b}\partial_{r}^{c}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})\partial_{p}^{b}\langle\langle\xi^{\prime
c}v_{i}^{a}(t)\rangle\rangle$ $-\frac{q_{s}q_{s^{\prime}}}{c^{2}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}\partial_{t}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})\partial_{p}^{b}\langle\langle
v_{i}^{a}(t)v_{j}^{b}(t^{\prime})\rangle\rangle$
$+\frac{q_{s}q_{s^{\prime}}}{c^{2}}\varepsilon^{bcd}\varepsilon^{dfg}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int dt^{\prime}\times$
$\times\partial_{r}^{f}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})\partial_{p}^{b}\langle\langle
v_{i}^{a}(t)v_{i}^{c}(t)v_{j}^{g}(t^{\prime})\rangle\rangle$
$+\frac{q_{s}q_{s^{\prime}}}{c^{2}}\varepsilon^{bcd}\varepsilon^{dfg}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int dt^{\prime}\times$
$\times\partial_{r}^{f}\partial_{r}^{h}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})\partial_{p}^{b}\langle\langle
v_{i}^{a}(t)v_{i}^{c}(t)v_{j}^{g}(t^{\prime})\xi^{h}\rangle\rangle=0.$ (75)

Finally, we present equation (75) in the self-consistent field approximation

$\partial_{t}j^{a}(\textbf{r},\textbf{p},t)+\partial_{r}^{b}\langle
v_{i}^{a}(t)v_{i}^{b}(t)\rangle+q_{s}\biggl{(}E^{b}_{ext}(\textbf{r},t)\partial_{p}^{b}\langle
v_{i}^{a}\rangle$ $+(\partial_{r}^{c}E^{b}_{ext})\partial_{p}^{b}\langle
v_{i}^{a}\xi^{c}\rangle+\varepsilon^{bcd}(B^{d}_{ext}\partial_{p}^{b}\langle
v_{i}^{a}v_{i}^{c}\rangle+(\partial^{f}B^{d})\partial_{p}^{b}\langle
v_{i}^{a}v_{i}^{c}\xi^{f}\rangle)\biggr{)}$
$-q_{s}q_{s^{\prime}}(\partial_{p}^{b}j^{a})\partial_{r}^{b}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})f(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})$
$-q_{s}q_{s^{\prime}}(\partial_{p}^{b}J_{D}^{ac})\partial_{r}^{b}\partial_{r}^{c}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})f(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})$
$+q_{s}q_{s^{\prime}}(\partial_{p}^{b}j^{a})\partial_{r}^{b}\partial_{r}^{c}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})d^{c}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})$
$-\frac{q_{s}q_{s^{\prime}}}{c^{2}}(\partial_{p}^{b}j^{a})\partial_{t}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})j^{c}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})$
$+\frac{q_{s}q_{s^{\prime}}}{c^{2}}\varepsilon^{bcd}\varepsilon^{dfg}(\partial_{p}^{b}\langle
v_{i}^{a}(t)v_{i}^{c}(t)\rangle)\times$ $\times\partial_{r}^{f}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})j^{g}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})$
$+\frac{q_{s}q_{s^{\prime}}}{c^{2}}\varepsilon^{bcd}\varepsilon^{dfg}(\partial_{p}^{b}\langle
v_{i}^{a}(t)v_{i}^{c}(t)\xi^{h}\rangle)\times$
$\times\partial_{r}^{f}\partial_{r}^{h}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})j^{g}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime})=0,$
(76)

where

$\langle\langle
v_{i}^{a}(t)\rangle\rangle(\textbf{r},\textbf{r}^{\prime},\textbf{p},\textbf{p}^{\prime},t,t^{\prime})=j^{a}(\textbf{r},\textbf{p},t)\cdot
f(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (77)
$\langle\langle\xi^{c}v_{i}^{a}(t)\rangle\rangle=J_{D}^{ac}(\textbf{r},\textbf{p},t)\cdot
f(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (78)
$\langle\langle\xi^{\prime
c}v_{i}^{a}(t)\rangle\rangle=j^{a}(\textbf{r},\textbf{p},t)\cdot
d^{c}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (79)
$\langle\langle
v_{i}^{a}(t)v_{j}^{b}(t^{\prime})\rangle\rangle=j^{a}(\textbf{r},\textbf{p},t)\cdot
j^{b}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (80)
$\langle\langle
v_{i}^{a}(t)v_{i}^{c}(t)v_{j}^{g}(t^{\prime})\rangle\rangle=\langle
v_{i}^{a}(t)v_{i}^{c}(t)\rangle\cdot
j^{g}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$ (81)

and

$\langle\langle
v_{i}^{a}(t)v_{i}^{c}(t)v_{j}^{g}(t^{\prime})\xi^{h}\rangle\rangle=\langle
v_{i}^{a}(t)v_{i}^{c}(t)\xi^{h}\rangle\cdot
j^{g}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}).$ (82)

Finally, we introduce electric $\textbf{E}_{int}$ and magnetic
$\textbf{B}_{int}$ fields caused by the charges of the system and present the
kinetic equation for the velocity distribution function in compact form:

$\partial_{t}j^{a}+\partial_{\textbf{r}}^{b}\langle
v_{i}^{a}(t)v_{i}^{b}(t)\rangle+q_{s}E^{b}\partial_{\textbf{r}}^{b}j^{a}+q_{s}(\partial_{\textbf{r}}^{c}E^{b})\partial_{\textbf{p}}^{b}J_{D}^{ac}$
$+\frac{q_{s}}{c}\varepsilon^{bcd}[B^{d}\partial_{\textbf{p}}^{b}\langle
v_{i}^{a}(t)v_{i}^{c}(t)\rangle+(\partial^{f}B^{d})\partial_{\textbf{p}}^{b}\langle
v_{i}^{a}(t)v_{i}^{c}(t)\xi^{f}\rangle]=0.$ (83)

##  IX On the extended closed set of kinetic equations containing scalar and
vector distribution functions

Equations for evolution of $f(\textbf{r},\textbf{p},t)$,
$d^{a}(\textbf{r},\textbf{p},t)$ and $j^{a}(\textbf{r},\textbf{p},t)$ (see
equations (61), (70), (83)) give incomplete set of equations. This set
requires additional assumptions to make final truncation.

Let us repeat equations (61), (70), (83) together to make additional analysis
of these equations

$\partial_{t}f+\textbf{v}\cdot\nabla
f+\nabla\cdot\textbf{F}+q_{s}\biggl{(}\textbf{E}+\frac{1}{c}\textbf{v}\times\textbf{B}\biggr{)}\cdot\frac{\partial
f}{\partial\textbf{p}}$
$+q_{s}\biggl{(}\partial^{b}\textbf{E}+\frac{1}{c}\textbf{v}\times\partial^{b}\textbf{B}\biggr{)}\cdot\frac{\partial
d_{s}^{b}}{\partial\textbf{p}}+q_{s}\biggl{(}\partial^{b}\partial^{c}\textbf{E}$
$+\frac{q_{s}}{c}\textbf{v}\times\partial^{b}\partial^{c}\textbf{B}\biggr{)}\cdot\frac{\partial
Q^{bc}}{\partial\textbf{p}}+\frac{q_{s}}{c}\varepsilon^{abc}\biggl{(}B^{c}\cdot\partial_{a,\textbf{p}}F^{b}$
$+\partial^{d}B^{c}\cdot\partial_{a,\textbf{p}}j_{D}^{bd}+\partial^{d}\partial^{f}B^{c}\cdot\partial_{a,\textbf{p}}j_{Q}^{bdf}\biggr{)}=0,$
(84)

and

$\partial_{t}d^{a}+\partial_{\textbf{r}}^{b}J_{D}^{ba}+q_{s}(\partial_{\textbf{p}}^{b}d^{a})E^{b}+q_{s}(\partial_{\textbf{p}}^{b}Q^{ac})\partial_{\textbf{r}}^{c}E^{b}$
$+\frac{q_{s}}{c}\varepsilon^{bcd}(\partial_{\textbf{p}}^{b}J^{ca}_{D})B^{d}+\frac{q_{s}}{c}\varepsilon^{bcd}(\partial_{\textbf{p}}^{b}J^{caf}_{Q})\partial_{\textbf{r}}^{f}B^{d}=0,$
(85)

and

$\partial_{t}j^{a}+\partial_{\textbf{r}}^{b}\langle
v_{i}^{a}(t)v_{i}^{b}(t)\rangle+q_{s}E^{b}\partial_{\textbf{r}}^{b}j^{a}+q_{s}(\partial_{\textbf{r}}^{c}E^{b})\partial_{\textbf{p}}^{b}J_{D}^{ac}$
$+\frac{q_{s}}{c}\varepsilon^{bcd}[B^{d}\partial_{\textbf{p}}^{b}\langle
v_{i}^{a}(t)v_{i}^{c}(t)\rangle+(\partial^{f}B^{d})\partial_{\textbf{p}}^{b}\langle
v_{i}^{a}(t)v_{i}^{c}(t)\xi^{f}\rangle]=0.$ (86)

Here we have equations for three distribution functions $f$, $d^{a}$, and
$j^{a}=v^{a}f+F^{a}$, but we have a number of other distribution functions
like $Q^{bc}$, $J_{D}^{bc}$, $J_{Q}^{bcd}$, $\langle
v_{i}^{a}(t)v_{i}^{c}(t)\rangle$, and $\langle
v_{i}^{a}(t)v_{i}^{c}(t)\xi^{f}\rangle$. We need to obtain the closed set of
equations at this step. Before we make this truncation, we need to provide
discussion of physical picture hidden in these functions.

Figure 2:  The illustration of the $\Delta_{\textbf{p}}$-vicinity in the
momentum space with few numbered particles being in the vicinity. Here, p is
the center of the vicinity, $\eta$ is the vector scanning the vicinity,
$\textbf{p}_{i}(t)$ are the momentums of particles, $\bar{\textbf{p}}$ is the
average value of the momentum of particles being in the vicinity, while vector
$\bar{\textbf{p}}$ differs from the center of the vicinity p.

###  IX.1 Kinetic theory of fluctuations

Considering transition from the kinetic theory (the field form of the
classical mechanics in the six dimensional space) to the hydrodynamic theory
(the field form of the classical mechanics in the three dimensional space) we
get the concentration of particles $n(\textbf{r},t)=\int
f(\textbf{r},\textbf{p},t)d\textbf{p}$, the current of particles
$\textbf{j}(\textbf{r},t)=n(\textbf{r},t)\textbf{v}(\textbf{r},t)$
$=\int\textbf{v}f(\textbf{r},\textbf{p},t)d\textbf{p}$. Difference
$\textbf{v}-\textbf{v}(\textbf{r},t)$ is the measure of chaotic motion, which
can be associated with the thermal motion. This interpretation getting closer
to the thermodynamic temperature if the system get though the process of
thermalization. For the quantum systems we also need to include other than
thermal mechanisms for the symmetric distribution of particles over the
quantum states in the momentum space. Major contribution is given by the Pauli
blocking existing in systems of fermions.

The average velocity $\int\textbf{v}f(\textbf{r},\textbf{p},t)d\textbf{p}$
gives the velocity field. The deviation from the average velocity
$\textbf{v}-\textbf{v}(\textbf{r},t)$ leads to other hydrodynamic functions,
like the temperature scalar field, the pressure tensor field, etc. These
functions are also smooth functions which are defined on certain scale. It
leaves as the question: how fluctuations can appear in our model? Since the
fluctuations is the natural phenomenon following from the individual motion of
particles governed be the microscopic equations.

Basically, the chaotic motion of particles gives small variation of
macroscopic functions due to the continuous unequal exchange of particles
between nearest vicinities. Hence, the concentration $n(\textbf{r},t)=\int
f(\textbf{r},\textbf{p},t)d\textbf{p}$ (2), the distribution function
$f(\textbf{r},\textbf{p},t)$ (4), etc show small variation over time and
coordinate. These fluctuations can be considered via the correlations like
space correlations $\overline{n(\textbf{r},t)n(\textbf{r}^{\prime},t)}$ or the
time correlations $\widetilde{n(\textbf{r},t)n(\textbf{r},t^{\prime})}$, where
we introduce additional average on the macroscopic space scale
$\bar{A}=(1/\Delta V)\int_{\Delta V}AdV$, or the average over the interval of
time $\Delta T$: $\tilde{A}=(1/\Delta T)\int_{\Delta T}Adt$.

The distribution function is constructed on a certain scale in the coordinate
and momentum space. If we want to trace the deviation from average
(particularly in the momentum space) we need to focus our attention on the
deviation of momentum of all particles in $\Delta$-vicinity from the middle
point of the vicinity $\eta$. Parameter $\Delta\textbf{v}_{i}(t)$ is also
associated with $\eta$. If we imagine a $\Delta$-vicinity with 5 particles (it
is an imaginary example, for the illustration, see Fig. 2, real vicinity
contain a large number of particles) we can count the average momentum
$\bar{\textbf{p}}=(1/5)\sum_{i=1}^{5}\textbf{p}_{i}(t)\neq\textbf{p}$ and we
get a deviation from p (the moment of the center of the vicinity). However,
the distribution function $f(\textbf{r},\textbf{p},t)$ count all these
particles as the particles with momentum p. In other words
$\sum_{i=1}^{5}\mbox{\boldmath$\eta$}_{i}\neq 0$ (in general $\eta$ is a
parameter scanning the vicinity, but here we introduce particular
$\mbox{\boldmath$\eta$}_{i}$ as
$\mbox{\boldmath$\eta$}_{i}(t)=\textbf{p}_{i}(t)-\textbf{p}$ the deviation of
momentum of particle being in the vicinity from the momentum of the center of
vicinity). This example shows that physical nature of fluctuations can be
presented via functions $\langle\eta^{a}\rangle$ and
$\langle\eta^{a}\eta^{b}\rangle$ (or $\langle\Delta v_{i}^{a}(t)\rangle$ and
$\langle\Delta v_{i}^{a}(t)\Delta v_{i}^{b}(t)\rangle$). It explains necessity
of the account of these functions in our model. We consider equation for
evolution of function $\langle\Delta v_{i}^{a}(t)\rangle$, while function
$\langle\Delta v_{i}^{a}(t)\Delta v_{i}^{b}(t)\rangle$ is assumed to be
approximately expressed via functions $f(\textbf{r},\textbf{p},t)$ and
$\langle\Delta v_{i}^{a}(t)\rangle$. Here, we introduce equation of state in
order to make truncation of the set of kinetic equations (in addition to the
self-consistent field approximation).

In the chosen moment of time $t$ we have
$\langle\Delta\textbf{v}_{i}\rangle\neq 0$. However, the fluctuations happen
on some time scale $\tau$, which can be associated with the time of
propagation of average particle via the $\Delta$-vicinity:
$\tau=\sqrt[3]{\Delta}/v_{0}$, where the average velocity of the
$\Delta_{p}$-vicinity is $v_{0}=pc/\sqrt{p^{2}+m_{s}^{2}c^{2}}$. We can go
further and make averaging of the kinetic equation over this time interval
$\tilde{a}=(1/\tau)\int_{0}^{\tau}adt$. Here, we find
$\widetilde{\langle\Delta\textbf{v}_{i}\rangle}=0$, but
$\widetilde{\langle\Delta v_{i}^{a}\Delta v_{i}^{b}\rangle}\neq 0$. As an
estimation we can choose $\Delta\textbf{v}_{i,max}=\sqrt[3]{\Delta_{v}}$
$=\sqrt[3]{\Delta_{p}}c/\sqrt{\sqrt[3]{\Delta_{v}}^{2}+m_{s}^{2}c^{2}}$.

###  IX.2 Truncation method

Let us consider simplified form of the kinetic equation for the vector
distribution function of velocity (76). Function
$j^{a}(\textbf{r},\textbf{p},t)$ can be represented via two other functions
introduced above $j^{a}(\textbf{r},\textbf{p},t)=v^{a}\cdot
f(\textbf{r},\textbf{p},t)$ $+F^{a}(\textbf{r},\textbf{p},t)$. Equation (76)
includes function $\langle v_{i}^{a}(t)v_{i}^{b}(t)\rangle$
$=\langle(v^{a}+\Delta v_{i}^{a})(v^{b}+\Delta v_{i}^{b})\rangle$
$=v^{a}v^{b}\langle 1\rangle+v^{a}\langle\Delta
v_{i}^{b}\rangle+v^{b}\langle\Delta v_{i}^{a}\rangle$ $+\langle\Delta
v_{i}^{a}\Delta v_{i}^{b}\rangle$ $=v^{a}v^{b}f(\textbf{r},\textbf{p},t)$
$+v^{a}F^{b}(\textbf{r},\textbf{p},t)$ $+v^{b}F^{a}(\textbf{r},\textbf{p},t)$
$+\langle\Delta v_{i}^{a}\Delta v_{i}^{b}\rangle$. Here we extract the
velocity from functions $j_{s}^{a}$, $j_{D,s}^{ab}$, $\langle
v_{i}^{a}(t)v_{i}^{b}(t)\rangle$, and $j_{Q,s}^{abc}$ to show the Lorentz
force in more familiar way. It helps us to make the truncation as well.
Finally, we represent $\langle v_{i}^{a}(t)v_{i}^{b}(t)\xi^{c}\rangle$
$=v^{a}v^{b}\langle\xi^{c}\rangle+v^{a}\langle\Delta
v_{i}^{b}\xi^{c}\rangle+v^{b}\langle\Delta v_{i}^{a}\xi^{c}\rangle$
$+\langle\Delta v_{i}^{a}\Delta v_{i}^{b}\xi^{c}\rangle$.

Our model includes the multipole moments in the coordinate and momentum space,
like $d^{a}$ and $F^{a}$. The multipole moments are defined on the scale of
the $\Delta$-vicinity. Therefore, they appear to be small. On the macroscopic
scale this value tends to zero. Hence, we can introduce the small parameter
$\varepsilon$ to indicate the relative value of these functions. We find
$d^{a}\sim\varepsilon f$, $F^{a}\sim\varepsilon f$,
$Q^{ab}\sim\varepsilon^{2}f$, $j^{ab}\sim\varepsilon^{2}f$, etc. To make the
truncation we need to estimate the high-rank tensor distribution functions
included in the model
$Q^{ab}=\langle\xi^{a}\xi^{b}\rangle=\delta^{ab}\sqrt[3]{\Delta^{2}}f(\textbf{r},\textbf{p},t)$
$\langle\Delta v^{a}\Delta
v^{b}\rangle=\delta^{ab}\sqrt[3]{\Delta_{v}^{2}}f(\textbf{r},\textbf{p},t)$,
and
$j_{D}^{ab}=\delta^{ab}\sqrt[3]{\Delta\Delta_{v}}f(\textbf{r},\textbf{p},t)$.
The third rank tensors are neglected. Function $\langle\Delta v^{a}\Delta
v^{b}\rangle$ is symmetric, so it is reasonable to be proportional to
$\delta^{ab}$. However, function $j_{D}^{ab}$ is not symmetric, but we assume
that deviations of same projections of the coordinate and momentum are
correlated. This correlation is included in the model via nonzero value of
$j_{D}^{ab}$, which can be traced in calculations of any specific problem.
Further comparison with experiment can show validity of our assumption. If no
correlation is found this function can be considered equal to zero in the
future modification of the model.

We start the discussion of equations (84), (85), (86) with the first of them
(84). We can extract terms proportional to the zero degree of the small
parameter $\varepsilon$:

$\partial_{t}f+\textbf{v}\cdot\nabla
f+q_{s}\biggl{(}\textbf{E}+\frac{1}{c}\textbf{v}\times\textbf{B}\biggr{)}\cdot\frac{\partial
f}{\partial\textbf{p}},$ (87)

which correspond to the well-known Vlasov equation. Next, we present terms
proportional to the first degree of the small parameter $\varepsilon$:

$q_{s}\biggl{(}\partial^{b}\textbf{E}+\frac{1}{c}\textbf{v}\times\partial^{b}\textbf{B}\biggr{)}\cdot\frac{\partial
d_{s}^{b}}{\partial\textbf{p}}+\nabla\cdot\textbf{F}+\frac{q_{s}}{c}\varepsilon^{abc}B^{c}\cdot\partial_{a,\textbf{p}}F^{b}.$
(88)

Other terms in equation (84) correspond to the second and third degree of the
small parameter $\varepsilon$. They are neglected since we include the first
correction to the Vlasov equation.

Next we consider equation (85). Lowest order on parameter $\varepsilon$ is
equal to one in this equation. We consider the lowest order and one correction
to it. Let us show term existing in the lowest order:

$\partial_{t}d^{a}+(\textbf{v}\cdot\nabla)d^{a}+q_{s}(\partial_{\textbf{p}}^{b}d^{a})E^{b}+\frac{q_{s}}{c}\varepsilon^{bcd}v^{c}(\partial_{\textbf{p}}^{b}d^{a})B^{d},$
(89)

and in the next order

$\partial_{\textbf{r}}^{b}j_{D}^{ba}+\frac{q_{s}}{c}\varepsilon^{bcd}(\partial_{\textbf{p}}^{b}j^{ca}_{D})B^{d}$
$+q_{s}(\partial_{\textbf{p}}^{b}Q^{ac})\partial_{\textbf{r}}^{c}E^{b}+\frac{q_{s}}{c}\varepsilon^{bcd}v^{c}(\partial_{\textbf{p}}^{b}Q^{af})\partial_{\textbf{r}}^{f}B^{d}.$
(90)

Let us consider equation (86). Function $j^{a}$ splits on two terms
$v^{a}f+F^{a}$. Obviously, equation (86) contains contribution of terms (87)
and (88) multiplied by $v^{a}$. Hence, all these terms are equation to zero in
accordance with the equation of evolution of the scalar distribution function.
After this simplification we find that the lowest order of terms in equation
(86) on the parameter $\varepsilon$ is equal to one. It contains the following
terms:

$\partial_{t}F^{a}+(\textbf{v}\cdot\nabla)F^{a}+q_{s}E^{b}\partial_{\textbf{p}}^{b}F^{a}+\frac{q_{s}}{c}\varepsilon^{bcd}v^{c}B^{d}\partial_{\textbf{p}}^{b}F^{a}$
(91)

In the next order on $\varepsilon$ we obtain

$\partial_{\textbf{r}}^{b}\langle\Delta v^{a}\Delta
v^{b}\rangle+q_{s}(\partial_{\textbf{r}}^{c}E^{b})\partial_{p}^{b}j_{D}^{ac}$
$+\frac{q_{s}}{c}\varepsilon^{bcd}B^{d}\partial_{\textbf{p}}^{b}\langle\Delta
v^{a}\Delta
v^{c}\rangle+\frac{q_{s}}{c}\varepsilon^{bcd}v^{c}(\partial_{\textbf{r}}^{f}B^{d})\partial_{\textbf{p}}^{b}j_{D}^{af}.$
(92)

###  IX.3 Final set of kinetic equations

Let us combine together final parts of equations (84), (85), (86) with the
equations of state

$\partial_{t}f+\textbf{v}\cdot\nabla
f+q_{s}\biggl{(}\textbf{E}+\frac{1}{c}\textbf{v}\times\textbf{B}\biggr{)}\cdot\frac{\partial
f}{\partial\textbf{p}}$
$+q_{s}\biggl{(}\partial^{b}\textbf{E}+\frac{1}{c}\textbf{v}\times\partial^{b}\textbf{B}\biggr{)}\cdot\frac{\partial
d_{s}^{b}}{\partial\textbf{p}}+\nabla\cdot\textbf{F}+\frac{q_{s}}{c}\varepsilon^{abc}B^{c}\cdot\partial_{a,\textbf{p}}F^{b}=0,$
(93)
$\partial_{t}d^{a}+(\textbf{v}\cdot\nabla)d^{a}+q_{s}(\partial_{\textbf{p}}^{b}d^{a})E^{b}+\frac{q_{s}}{c}\varepsilon^{bcd}v^{c}(\partial_{\textbf{p}}^{b}d^{a})B^{d}$
$+\sqrt[3]{\Delta\Delta_{v}}\partial_{\textbf{r}}^{a}f-\frac{q_{s}}{c}\sqrt[3]{\Delta\Delta_{v}}\varepsilon^{abd}(\partial_{\textbf{p}}^{b}f)B^{d}$
$+q_{s}\sqrt[3]{\Delta^{2}}(\partial_{\textbf{p}}^{b}f)\partial_{\textbf{r}}^{a}E^{b}+\frac{q_{s}}{c}\sqrt[3]{\Delta^{2}}\varepsilon^{bcd}v^{c}(\partial_{\textbf{p}}^{b}f)\partial_{\textbf{r}}^{a}B^{d}=0,$
(94)

and

$\partial_{t}F^{a}+(\textbf{v}\cdot\nabla)F^{a}+q_{s}E^{b}\partial_{\textbf{p}}^{b}F^{a}+\frac{q_{s}}{c}\varepsilon^{bcd}v^{c}B^{d}\partial_{\textbf{p}}^{b}F^{a}$
$+\sqrt[3]{\Delta_{v}^{2}}\partial_{\textbf{r}}^{a}f+q_{s}\sqrt[3]{\Delta\Delta_{v}}(\partial_{\textbf{r}}^{a}E^{b})\partial_{\textbf{p}}^{b}f$
$-\frac{q_{s}}{c}\sqrt[3]{\Delta_{v}^{2}}\varepsilon^{abd}B^{d}\partial_{\textbf{p}}^{b}f+\frac{q_{s}}{c}\sqrt[3]{\Delta\Delta_{v}}\varepsilon^{bcd}v^{c}(\partial_{\textbf{r}}^{a}B^{d})\partial_{\textbf{p}}^{b}f=0.$
(95)

The electromagnetic field in equations (93), (94), (95) obeys the following
Maxwell equations

$\begin{array}[]{cc}\nabla\cdot\textbf{B}_{int}=0,&\nabla\times\textbf{E}_{int}=-\frac{1}{c}\partial_{t}\textbf{B}_{int},\end{array}$
(96) $(\nabla\times\textbf{B}_{int})^{a}=\frac{1}{c}\partial_{t}E_{int}^{a}+$
$+\frac{4\pi}{c}\sum_{s}q_{s}\biggl{(}\int
v^{a}f_{s}(\textbf{r},\textbf{p},t)d\textbf{p}+\partial^{b}\int
v^{a}d_{s}^{b}(\textbf{r},\textbf{p},t)d\textbf{p}$ $+\int
F_{s}^{a}(\textbf{r},\textbf{p},t)d\textbf{p}+\sqrt[3]{\Delta\Delta_{v}}\partial^{a}\int
f(\textbf{r},\textbf{p},t)d\textbf{p}\biggr{)},$ (97)
$\nabla\cdot\textbf{E}_{int}=4\pi\sum_{s}q_{s}\biggl{(}\int
f_{s}(\textbf{r},\textbf{p},t)d\textbf{p}$ $+\partial^{b}\int
d_{s}^{b}(\textbf{r},\textbf{p},t)d\textbf{p}+\frac{1}{2}\sqrt[3]{\Delta^{2}}\partial^{b}\partial^{b}\int
f(\textbf{r},\textbf{p},t)d\textbf{p}+...\biggr{)}.$ (98)

which appear as corresponding modification of equations (64), (62), and (63).

##  X The spin-electron-acoustic waves propagation in the spin polarized
electron gas of high density

The degenerate macroscopically motionless (being in the equilibrium state)
electron gas are described via the Vlasov kinetic equation for each spin
projection of electrons Andreev PP 16 SSE kin

$\partial_{t}f_{s}+\textbf{v}\cdot\nabla
f_{s}+q_{s}\biggl{(}\textbf{E}+\frac{1}{c}\textbf{v}\times\textbf{B}\biggr{)}\cdot\frac{\partial
f_{s}}{\partial\textbf{p}}=0.$ (99)

Subindex $s$ corresponds to the electrons with the spin-up $s=\uparrow$ or
spin-down $s=\downarrow$. We consider small perturbations of the equilibrium
state $f_{s}=f_{0s}+\delta f_{s}$, where $f_{0s}$ is the equilibrium
distribution function, and $\delta f_{s}$ is the perturbation of the
distribution function, which is chosen as the plane wave in the coordinate
space $\delta f_{s}=F_{s}e^{-\imath\omega t+\imath k_{z}z}$, with the constant
amplitude $F_{s}$. It lead to the linearized kinetic equation

$-\imath(\omega-k_{z}v_{z})\delta
f_{s}+q_{s}\delta\textbf{E}\cdot\frac{\partial f_{0s}}{\partial\textbf{p}}=0,$
(100)

where $\delta\textbf{B}=0$ for the longitudinal waves.

Equilibrium distribution functions for each subspecies of electrons is chosen
as the Fermi step $f_{0s}=\theta(p_{Fs}-p)/(2\pi\hbar)^{3}$, where
$p_{Fs}=(6\pi n_{0s})^{1/3}\hbar$ is the radius of the Fermi sphere in the
momentum space for species $s$.

We consider the longitudinal waves, hence the perturbation of the electric
field is parallel to the direction of the wave propagation
$\textbf{k}\parallel\delta\textbf{E}$.

Equation (100) give expression for the perturbation of the distribution
function of electrons

$\delta f_{s}=-\imath q_{s}(\textbf{v}\cdot\delta\textbf{E})\frac{\partial
f_{0s}}{\partial\varepsilon}\frac{1}{\omega-k_{z}v_{z}},$ (101)

where

$\frac{\partial f_{0s}}{\partial\varepsilon}=\frac{\partial f_{0s}}{\partial
p}\frac{\partial
p}{\partial\varepsilon}=-\frac{1}{(2\pi\hbar)^{3}}\delta(p-p_{Fs})\frac{pc^{2}}{\varepsilon}.$
(102)

Expression (101) allows us to calculate the perturbations of concentration of
the degenerate electrons

$\delta n_{s}=\int\delta f_{s}(\textbf{r},\textbf{p},t)d\textbf{p}$
$=-\frac{2\pi\imath q_{s}p_{Fs}^{2}\delta
E_{z}}{(2\pi\hbar)^{3}k_{z}v_{Fs}}\biggl{[}2+\frac{\omega}{k_{z}v_{Fs}}ln\biggl{(}\frac{\omega-
k_{z}v_{Fs}}{\omega+k_{z}v_{Fs}}\biggr{)}\biggr{]}.$ (103)

It leads to the dispersion equation

$1+\frac{3}{2}\frac{\omega_{Le}^{2}}{k_{z}^{2}c^{2}}\sum_{s=u,d}\frac{n_{0s}}{n_{0e}}\gamma_{Fs}\frac{m_{s}^{2}c^{2}}{p_{Fs}^{2}}\biggl{[}2+\frac{\omega}{k_{z}v_{Fs}}ln\biggl{(}\frac{\omega-
k_{z}v_{Fs}}{\omega+k_{z}v_{Fs}}\biggr{)}\biggr{]}=0.$ (104)

Presented dispersion equation is found for the motionless ions, but we have
two species of particles: spin-up electrons and spin-down electrons. Both
species of electrons are degenerate, but they have different equilibrium
concentrations. So we have nonzero equilibrium spin polarization. The spin
polarization is caused by the external magnetic field. However, it is not
included in equation (100), since is does not affect the final result (101)
for the waves propagating parallel to the magnetic field.

Equation (104) gives two solutions: the Langmuir wave and the SEAW. The SEAW
has been studied in different regimes and using different method. It is
considered within nonrelativistic hydrodynamics Andreev PRE 15 ,
nonrelativistic kinetics Andreev PP 16 SSE kin , and relativistic
hydrodynamics Andreev 2112 . Here, we present relativistic kinetic theory of
the SEAW. Numerical analysis of equation (104) shows that the dispersion
dependencies of the Langmuir wave and the SEAW getting close to each other at
the large wave vectors, which leads to an instability in the short-wavelength
limit. However, it can be subject of another work specified on the
instabilities in degenerate plasmas related to the SEAWs, while this paper is
focused on derivation of the extended kinetic model.

##  XI Conclusion

Detailed derivation of kinetic equations for the relativistic plasmas has been
presented. It concludes in the set of three kinetic equations (for each
species of particles): the Vlasov kinetic equation for the scalar distribution
function, the equation of evolution of the vector distribution function of the
electric dipole moment (the fluctuation of local center of mass of the
particles being in the physically infinitesimal volume), and the equation of
evolution of the vector distribution function of the velocity (the fluctuation
of average momentum of the particles being in the physically infinitesimal
volume). This closed set of equations appears as the result of truncation of
the quasi-infinite set of kinetic equations (number of degrees of freedom is
restricted due to the finite number of particles N), which extends via the
account of two-, three- and more coordinate distribution functions, and higher
multipole moments of one-coordinate (along with two-, three- and more
coordinate) distribution functions.

The method of derivation shows the representation of the microscopic motion of
the individual particles in terms of the macroscopic distribution functions.
It gives representation of deterministic classical mechanics in terms of
evolution collective functions. No probabilistic approaches are used for the
derivation or interpretation of presented equations and found functions.

Equations have been found for the relativistic plasmas meaning both the large
temperatures of the system (large velocities of the chaotic motion) and large
velocities of the ordered motion. The motion of particles with the velocities
close to the speed of light requires the account of full retarding potentials
for the electromagnetic field created by particles. It reflects in the full
Maxwell equations obtained in the self-consistent field approximation. The
Maxwell equations also include the distribution function of the electric
dipole moment, distribution function of the velocity, etc as the sources of
the electromagnetic field.

The model is found for the hot plasmas, but it can be used for the degenerate
electron gas with the Fermi velocity close to the speed of light. Hence, it is
applied to the spin-electron acoustic waves in the extremely dense plasmas.
The SEAWs are considered in the regime of account of the scalar distribution
function governed by the Vlasov equation.

##  XII DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or
analyzed in this study, which is a purely theoretical one.

##  Appendix A Potentials of the electromagnetic field in the multipole
approximation

The application of the full retarding potentials for the description of the
field created by particles (7) and (8) reveals in the full set of Maxwell
equations found in the self-consistent field approximation (64), (62), and
(63). However, firstly it appears (in the self-consistent field approximation)
in the kinetic equation as the integral terms simplifies after splitting of
the two-coordinate distribution functions on the products of the corresponding
one-coordinate distribution functions. These integral terms are shown in the
text, but we want to specify the structure of the scalar and vector potentials
of the electromagnetic field caused by different distribution functions since
the Maxwell equations (64), (62), and (63) show the combined contribution of
all terms.

Therefore, the partial scalar potentials are:

$\varphi_{0}=\int d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})f(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$
(105)

$\varphi_{D}=-\partial^{c}_{\textbf{r}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})d^{c}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$
(106)

and

$\varphi_{Q}=\frac{1}{2}\partial^{c}_{\textbf{r}}\partial^{b}_{\textbf{r}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})q^{bc}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}).$
(107)

The partial vector potentials are

$A^{a}_{0}=\frac{1}{c}\int d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})j^{a}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$
(108) $A_{D}^{a}=-\frac{1}{c}\partial^{b}_{\textbf{r}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})J_{D}^{ab}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}),$
(109)

and

$A_{Q}^{a}=\frac{1}{2}\frac{1}{c}\partial^{c}_{\textbf{r}}\partial^{b}_{\textbf{r}}\int
d\textbf{r}^{\prime}d\textbf{p}^{\prime}\int
dt^{\prime}G(t,t^{\prime},\textbf{r},\textbf{r}^{\prime})J_{Q}^{abc}(\textbf{r}^{\prime},\textbf{p}^{\prime},t^{\prime}).$
(110)

These potentials correspond to the Maxwell equations presented in the text
(62), (64), (63).

## References

  * (1) Yu. L. Klimontovich, _Statistical Physics_ Harwood, New York (1986).
  * (2) W. Israel, J. M. Stewart, Transient relativistic thermodynamics and kinetic theory, Annals Phys. 118, 341 (1979).
  * (3) P. Romatschke, ”New developments in relativistic viscous hydrodynamics”, Int. J. Mod. Phys. E 19, 1 (2010).
  * (4) P. Kovtun, ”Lectures on hydrodynamic fluctuations in relativistic theories” J. Phys. A: Math. Theor. 45, 473001 (2012).
  * (5) G. S. Denicol, H. Niemi, E. Molnar, D. H. Rischke, ”Derivation of transient relativistic fluid dynamics from the Boltzmann equation” Phys. Rev. D 85, 114047 (2012).
  * (6) U. Heinz, R. Snellings, ”Collective Flow and Viscosity in Relativistic Heavy-Ion Collisions”, Annu. Rev. Nucl. Part. Sci. 63, 123 (2013).
  * (7) W. Florkowski, M. P. Heller, M. Spalinski, ”New theories of relativistic hydrodynamics in the LHC era”, Rep. Prog. Phys. 81, 046001 (2018).
  * (8) G. S. Denicol, E. Molnar, H. Niemi, D. H. Rischke, ”Resistive dissipative magnetohydrodynamics from the Boltzmann-Vlasov equation”, Phys. Rev. D 99, 056017 (2019).
  * (9) Remi Hakim, Introduction to Relativistic Statistical Mechanics Classical and Quantum, World Scientific Publishing Co. Pte. Ltd., 2011.
  * (10) M. A. Drofa, L. S. Kuz’menkov, ”Continual approach to multiparticle systems with long-range interaction. Hierarchy of macroscopic fields and physical consequences”, Theoretical and Mathematical Physics 108, 849 (1996).
  * (11) L. S. Kuzmenkov, Theoretical Physics: Classical Mechanics (Nauka, Moscow, 2015) [in Russian].
  * (12) S. Weinberg, Gravitation and Cosmology (John Wiley and Sons, Inc., NewYork, 1972).
  * (13) S.R. de Groot, L.G. Suttorp, ”Foundations of Electrodynamics” (North-Holland 1972).
  * (14) M. V. Kuzelev, A. A. Rukhadze, ”On the quantum description of the linear kinetics of a collisionless plasma”, Phys. Usp. 42, 603–605 (1999).
  * (15) L. S. Kuz’menkov and P. A. Andreev, ”Microscopic Classic Hydrodynamic and Methods of Averaging”, presented in PIERS Proceedings, p. 158, August 19-23, Moscow, Russia 2012.
  * (16) P. A. Andreev, ”On the structure of relativistic hydrodynamics for hot plasmas”, Phys. Scr. 97, 085602 (2022).
  * (17) P. A. Andreev, ”Quantum kinetics of spinning neutral particles: General theory and Spin wave dispersion” Physica A 432, 108 (2015).
  * (18) V. V. Zubkov, ”Distribution functions for continuous medium without probability hypotheses”, J. Phys.: Conf. Ser. 1352, 012067 (2019).
  * (19) V.V. Zubkov, A.V. Zubkova, ”Irreversibility in classical kinetic theory: retardation of interaction and distribution functions”, arXiv:2006.11565.
  * (20) V. V. Zubkov, D. A. Mayfat, Yashkin K.Yu. Tensor field method in the linear response theory // Vestnik NovSU. Issue: Engineering Sciences. 2022. № 3(128). P.21–25. DOI: https://doi.org/10.34680/2076-8052.2022.3(128).21-25
  * (21) L. S. Kuz’menkov, S. G. Maksimov, ”Quantum hydrodynamics of particle systems with Coulomb interaction and quantum Bohm potential,” Theor. Math. Phys. 118, 227 (1999).
  * (22) L. S. Kuz’menkov, S. G. Maksimov, V. V. Fedoseev, ”Microscopic quantum hydrodynamics of systems of fermions: part I”, Theor. Math. Phys. 126, 110 (2001).
  * (23) P. A. Andreev, Hydrodynamic model of a Bose–Einstein condensate with anisotropic short-range interaction and bright solitons in a repulsive Bose–Einstein condensate, Laser Phys. 29, 035502 (2019).
  * (24) P. A. Andreev, ”Extended hydrodynamics of degenerate partially spin polarized fermions with short-range interaction up to the third order by interaction radius approximation”, Laser Phys. 31, 045501 (2021).
  * (25) P. A. Andreev, ”Hydrodynamics of the atomic Bose–Einstein condensate beyond the mean-field approximation”, Laser Phys. Lett. 18, 055501 (2021).
  * (26) P. A. Andreev, ”Hydrodynamics of quantum corrections to the Coulomb interaction via the third rank tensor evolution equation: application to Langmuir waves and spin-electron acoustic waves”, J. Plasma Phys. 87, 905870511 (2021).
  * (27) P. A. Andreev, ”Quantum hydrodynamic theory of quantum fluctuations in dipolar Bose–-Einstein condensate”, Chaos 31, 023120 (2021).
  * (28) P. A. Andreev, I. N. Mosaki, and M. I. Trukhanova, ”Quantum hydrodynamics of the spinor Bose–Einstein condensate at non-zero temperatures”, Phys. Fluids 33, 067108 (2021).
  * (29) P. A. Andreev, ”Separated spin-up and spin-down quantum hydrodynamics of degenerated electrons”, Phys. Rev. E 91, 033111 (2015).
  * (30) P. A. Andreev, ”Spin-electron acoustic waves: The Landau damping and ion contribution in the spectrum”, Phys. Plasmas 23, 062103 (2016).
  * (31) J. C. Ryan, ”Collective excitations in a spin-polarized quasi-two-dimensional electron gas”, Phys. Rev. B 43, 4499 (1991).
  * (32) A. Agarwal, M. Polini, R. Fazio, G. Vignale, ”Persistent Spin Oscillations in a Spin-Orbit-Coupled Superconductor”, Phys. Rev. Lett. 107, 077004 (2011).
  * (33) A. Agarwal, M. Polini, G. Vignale, M. E. Flatte, ”Long-lived spin plasmons in a spin-polarized two-dimensional electron gas”, Phys. Rev. B 90, 155409 (2014).
  * (34) M. Schober, D. Kreil, and H. M. Bohm, ”Dynamic response of partially spin- and valley-polarised two-dimensional electron liquids”, EPL 129, 17001 (2020).
  * (35) I. D’Amico, F. Perez, and C. A. Ullrich, ”Chirality and intrinsic dissipation of spin modes in two-dimensional electron liquids”, J. Phys. D: Appl. Phys. 52, 203001 (2019).
  * (36) C. Yang, Z.-C. Wang, G. Su, ”Spin-polarized plasmon in ferromagnetic metals”, Physica A 575, 126043 (2021).
  * (37) D. Kreil, C. Staudinger, K. Astleithner, H. M. Bohm, ”Resonant and anti-resonant modes of the dilute, spin-inbalanced, two-dimensional electron liquid including correlations”, Contrib. Plasma Phys. 58, 179 (2018).
  * (38) D. Kreil, R. Hobbiger, J. T. Drachta, and H. M. Bohm, ”Excitations in a spin-polarized two-dimensional electron gas”, Physical Review B 92, 205426 (2015).
  * (39) C. Huang, Y. D. Chong, G. Vignale, and M. A. Cazalilla, ”Graphene electrodynamics in the presence of the extrinsic spin Hall effect”, Physical Review B 93, 165429 (2016)
  * (40) A. N. Afanasiev, ”Acoustic plasmons and isotropic short-range interaction in two-component electron liquids”, Physical Review B 106, 224301 (2022).
  * (41) P. A. Andreev, ”Spin-electron-acoustic waves and solitons in high-density degenerate relativistic plasmas”, Phys. Plasmas 29, 122102 (2022).
  * (42) P. A. Andreev, ”Microscopic model for relativistic hydrodynamics of ideal plasmas”, arXiv:2109.14050.
  * (43) P. A. Andreev, ”Relativistic hydrodynamic model with the average reverse gamma factor evolution for the degenerate plasmas: high-density ion-acoustic solitons”, Phys. Plasmas 29, 062109 (2022).
  * (44) L. D. Landau and E. M. Lifshitz, _The Classical Theory of Fields_ (Butterworth-Heinemann, 1975). L. D. Landau, E.M. Lifshitz, ”Field theory” , Vol. 2 of Course of Theoretical Physics.
  * (45) P. A. Andreev, ”Hydrodynamic and kinetic representation of the microscopic dynamics as the transitions on the macroscopic scale of description and meaning of the self-consistent field approximation in these models”, arXiv:2204.11141.
  * (46) A. A. Vlasov, J. Exp. Theor. Phys. 8, 291 (1938); A. A. Vlasov Sov. Phys. Usp. 10, 721 (1968).