2. Derivation of hydrodynamic equations tracing the microscopic motion of particles

The self-consistent field or mean-field approximation is well established in plasma physics. However, the derivation of the hydrodynamic and kinetic equations including the explicit operator of averaging on the physically infinitesimal volume suggested by Drofa & Kuz'menkov (Reference Drofa and Kuz'menkov1996) and Kuz'menkov (Reference Kuz'menkov1991) enables a deeper look on this approximation.

Our goal in this paper is to consider the exact dynamics of the arbitrary system of the charged particles and represent this dynamical evolution in terms of functions suitable for the collective processes instead of parameters characterizing each particle in the system like the coordinates and momentums.

From a physical point of view, there is the question on the value of the $\varDelta$-neighbourhood. The problem of estimation of the physically infinitesimal volume is addressed in the literature. For example, Klimontovich (Reference Klimontovich1986) in his book (see also Klimontovich Reference Klimontovich1967) gives this estimation for two regimes: the rarefied gas of neutral atoms and the rarefied plasmas. This estimation shows that the characteristic length for plasmas is of the order of the Debye length $r_{D}=(T/(4{\rm \pi} \sum _{s}n_{0s}q_{s}^{2}))^{1/2}$. However, this clear statement contains a contradiction for the hydrodynamics. Nevertheless, this estimation is confirmed below for the kinetic model. The physically infinitesimal volume has some non-zero value from the microscopic point of view, but it is zero volume on the macroscopic scale. Therefore, it should be expressed via the microscopic parameters related to the motion of individual particles. Let us specify that some authors refer to the notion ‘microscopic’ as to the kinetic description, while ‘macroscopic’ is reserved for the hydrodynamics. In contrast, we call ‘microscopic’ the scale, where we consider the motion of the individual particles, while both the kinetic and hydrodynamic models are considered as ‘macroscopic’. Hence, the macroscopic scale is the scale, where the number density and the distribution function presented below are continuous functions (as continuous as they can be in a world of discrete particles and atoms). The calculation of the Debye length is made on the macroscopic scale. The expression of the Debye length is found via the macroscopic parameters as well (so any attempt to find microscopic derivation does not change our conclusion). Hence, it defines some macroscopic volume which does not correspond to the macroscopically zero volume of the physically infinitesimal volume. Nevertheless, it is essential to have a large number of particles in the physically infinitesimal volume to insure the continuity of hydrodynamic functions. For the low-density plasmas, we expect to have the following relations of the characteristic lengths: $r_{B}\ll a\equiv (n)^{-1/3}\ll \varDelta ^{1/3}\ll r_{D}$, where $r_{B}=\hbar ^{2}/m_{e}e^{2}$ is the Bohr radius. For the high-density plasmas, we have modified the relations of the characteristic lengths: $r_{B}\sim a\ll \varDelta ^{1/3}\ll r_{D}$.

The number density of particles is traditionally used as one of the functions describing the collective effects in the systems of many particles. Moreover, the number density can be introduced as the exact distribution of particles in the coordinate physical space:

(2.1)\begin{equation} n_{m}(\boldsymbol{r},t)=\sum_{i=1}^{N}\delta(\boldsymbol{r}-\boldsymbol{r}_{i}(t)), \end{equation}

where $\boldsymbol {r}_{i}(t)$ is the coordinate of the $i$th particle, and subindex $m$ in $n_{m}$ refers to the fact that we consider the microscopic number density. We do not know the value of coordinates of particles $\boldsymbol {r}_{i}(t)$. However, we do not need to know this information. We need to know the equations of motion of each particle. The equations of motion of particles allow us to derive equations for the evolution of the collective variables. So we will discuss the properties of the system in terms of the collective motion with no further references to the coordinates.

We consider the system of classic particles. We consider the elastic interactions, so we model the dynamics of all particles as structureless objects,

(2.2)\begin{equation} \dot{\boldsymbol{p}}_{i}(t)=\boldsymbol{F}(\boldsymbol{r}_{i}(t), \boldsymbol{v}_{i}(t),t), \end{equation}

where $\boldsymbol {p}_{i}(t)$ is the momentum, which is the function of time, and $\boldsymbol {F}_{i}=\boldsymbol {F}(\boldsymbol {r}_{i}(t),\boldsymbol {v}_{i}(t),t)$ is the force acting on the $i$th particle being in point $\boldsymbol {r}_{i}(t)$. All interactions between objects happens via the fields (usually the electromagnetic field, while gravitational and nuclear fields usually exert no influence on atomic or plasma effects). The force is the projection of the corresponding force field on the trajectory of the $i$th particle $\boldsymbol {F}_{i}=q_{i}\boldsymbol {E}_{i}(\boldsymbol {r}_{i}(t)) +(1/c)q_{i}\boldsymbol {v}_{i}(t)\times \boldsymbol {B}_{i}(\boldsymbol {r}_{i}(t))$ with

(2.3)\begin{equation} \boldsymbol{E}(\boldsymbol{r}_{i}(t),t)=\int {\rm d} \boldsymbol{r}\boldsymbol{E}(\boldsymbol{r},t)\delta(\boldsymbol{r}- \boldsymbol{r}_{i}(t)) \end{equation}

and

(2.4)\begin{equation} \boldsymbol{B}(\boldsymbol{r}_{i}(t),t)=\int {\rm d} \boldsymbol{r}\boldsymbol{B}(\boldsymbol{r},t) \delta(\boldsymbol{r}-\boldsymbol{r}_{i}(t)). \end{equation}

The force $\boldsymbol {F}(\boldsymbol {r}_{i}(t),t)$ is the superposition of interactions with other particles in the system.

Actually, we do not need to know the form of interaction to derive an equation for evolution of function (2.1). Its time derivative gives the following relation:

(2.5)\begin{equation} \partial_{t}n_{m}+\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{j}_{m}=0, \end{equation}

where

(2.6)\begin{equation} \boldsymbol{j}_{m}(\boldsymbol{r},t)=\sum_{i=1}^{N} \boldsymbol{v}_{i}(t)\delta(\boldsymbol{r}-\boldsymbol{r}_{i}(t)) \end{equation}

is the microscopic current, $\boldsymbol {v}_{i}(t)=\dot {\boldsymbol {r}}_{i}(t)$ is the velocity of the $i$th particle, and we also assume that each particle is stable, that is, it does not decay on other particles during evolution. The processes of creation/annihilation or ionization/recombination are not included in our analysis.

The number density (2.1) is the collective variable. However, it is constructed on the microscopic scale. Therefore, the function shows the exact position of each particle in some point of space.

We can make the transition to the macroscopic scale. To this end, we need to introduce the scale giving the macroscopically infinitesimal volume. We use the notation $\varDelta$ for this volume.

In each moment of time $t$, we consider each point of space $\boldsymbol {r}$. We construct the $\varDelta$-neighbourhood of each point of space $\boldsymbol {r}$ and calculate the number density on the chosen scale,

(2.7)\begin{equation} n(\boldsymbol{r},t)=\frac{N(\boldsymbol{r},t)}{\varDelta}. \end{equation}

However, we do not know the number of particles $N(\boldsymbol {r},t)$ in the $\varDelta$-neighbourhood of any point. So this definition is not useful for the further calculation.

Next, we make the proper generalization of definition (2.7) (Drofa & Kuz'menkov Reference Drofa and Kuz'menkov1996)

(2.8)\begin{equation} n_{e}(\boldsymbol{r},t)=\frac{1}{\varDelta}\int_{\varDelta}{\rm d} \boldsymbol{\xi}\sum_{i=1}^{N/2}\delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)). \end{equation}

The $\varDelta$-neighbourhood presented in this formula is illustrated in figure 1. The hydrodynamic model of plasmas requires the introduction of the number density for each species. Each number density and other hydrodynamic functions evolve under the action of all species in the system. We specify that we consider the quasi-neutral plasmas of two species: the electron–ion plasmas. We use the following numeration of particles: $i\in [1,N/2]$ for the electrons and $i\in [N/2+1,N]$ for the ions (below, in the relativistic regime, we consider protons only). We illustrate the derivation following the evolution of the number density of electrons. The explicit contribution of ions is shown in the terms describing the interaction. Similar definition can be introduced for the ions/protons,

(2.9)\begin{equation} n_{p}(\boldsymbol{r},t)= \frac{1}{\varDelta}\int_{\varDelta}{\rm d} \boldsymbol{\xi}\sum_{i=1+N/2}^{N} \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)). \end{equation}

Figure 1.

Illustration of the $\varDelta$-neighbourhood. Vector $\boldsymbol {\xi }$ scanning the $\varDelta$-neighbourhood is illustrated.

Equation (2.8), for electrons, can be interpreted as the action of the operator of averaging (Drofa & Kuz'menkov Reference Drofa and Kuz'menkov1996; Kuz'menkov & Andreev Reference Kuz'menkov and Andreev2012)

(2.10)\begin{equation} \langle \cdots\rangle\equiv\frac{1}{\varDelta}\int_{\varDelta}{\rm d} \boldsymbol{\xi} \sum_{i=1}^{N} \cdots\delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)), \end{equation}

which calculates the number of particles in the $\varDelta$-neighbourhood checking for the presence of each particle in the chosen neighbourhood scanning the neighbourhood by means vector $\boldsymbol {\xi }$. Hence, (2.8) can be rewritten as $n_{e}(\boldsymbol {r},t)= \langle s_{i}\rangle$, where $s_{i}=1$ is the dimensionless indicator. Otherwise, it can be rewritten by mass $n_{e}(\boldsymbol {r},t)= \langle m_{i}\rangle /m_{e}$, where $m_{i}/m_{e}=1$ since we make summation over electrons $i\in [1,N/2]$. Operator (2.10) can be replaced by symbol $\langle \cdots \rangle$ to express the equations in shorter form.

This method is suggested by Kuz'menkov L.S., it appears as the generalization of method suggested by Klimontovich (Reference Klimontovich1967) (see also Weinberg Reference Weinberg1972). However, as it is mentioned in § 1, there is some different physical insight. Therefore, we introduce the space average on the physically infinitesimal volume in contrast to the average on the ensemble of physical systems.

Similarly to (2.1), (2.5) and (2.6), we can derive the continuity equation for the number density (2.8). To get the derivation, we differentiate (2.8) with respect to time and obtain the continuity equation

(2.11)\begin{equation} \partial_{t}n_{e}+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{j}_{e}=0, \end{equation}

where the current $\boldsymbol {j}$ has the following definition:

(2.12)\begin{equation} \boldsymbol{j}_{e}(\boldsymbol{r},t)= \frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\sum_{i=1}^{N/2} \boldsymbol{v}_{i}(t) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)). \end{equation}

We can also introduce the velocity field $\boldsymbol {v}_{e}=\boldsymbol {j}_{e}/n_{e}$. The average velocity is found as the arithmetic mean for all particles being in the neighbourhood. For the current of the single species, we find that the momentum density $\boldsymbol {P}$ is proportional to the current $\boldsymbol {P}_{e}=m_{e} \boldsymbol {j}_{e}$ (for non-relativistic systems only).

To continue the derivation of hydrodynamics, we consider the time evolution of the current $\boldsymbol {j}_{e}$. In this case, we need the expression for the acceleration with the explicit form of interaction

(2.13)\begin{equation} \dot{\boldsymbol{v}}_{i}(t)=\frac{1}{m_{i}} \left(q_{i}\boldsymbol{E}_{\rm ext}(\boldsymbol{r}_{i}(t),t) +\frac{1}{c}q_{i}[\boldsymbol{v}_{i}(t), \boldsymbol{B}_{\rm ext}(\boldsymbol{r}_{i}(t),t)] -\sum_{j=1, j\neq i}^{N}q_{i}q_{j}\nabla_{i}G_{ij}\right), \end{equation}

where $G_{ij}=1/r_{ij}$ is the Green function of the Coulomb interaction, $r_{ij}=|\boldsymbol {r}_{i}(t)-\boldsymbol {r}_{j}(t)|$.

In this section, we consider the non-relativistic plasmas. Therefore, we consider interaction in the quasi-static limit, which is presented by the Coulomb interaction.

The action of the time derivative on the current (2.12) leads to the action of the time derivative on the product of two functions under the integral $\boldsymbol {v}_{i}(t) \delta (\boldsymbol {r}+ \boldsymbol {\xi }-\boldsymbol {r}_{i}(t))$. Hence, the result is the superposition of two terms under the integral. One term $\dot {\boldsymbol {v}}_{i}(t) \delta (\boldsymbol {r}+ \boldsymbol {\xi }-\boldsymbol {r}_{i}(t))$ requires the acceleration (2.13), and the second term $-\boldsymbol {v}_{i}(t) (\boldsymbol {v}_{i}(t)\boldsymbol {\cdot }\boldsymbol {\nabla }) \delta (\boldsymbol {r}+ \boldsymbol {\xi }-\boldsymbol {r}_{i}(t))$ leads to the flux of the momentum.

It gives the general structure of the Euler equation or the momentum balance evolution equation

(2.14)\begin{equation} \partial_{t}j_{e}^{a}+\partial_{b}\varPi_{e}^{ab}= \frac{1}{m_{e}}(\varPhi^{a}_{e,{\rm ext}}+\varPhi_{e}^{a}) , \end{equation}

where

(2.15)\begin{equation} \varPi_{e}^{ab}(\boldsymbol{r},t)= \frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi} \sum_{i=1}^{N/2} v_{i}^{a}(t) v_{i}^{b}(t) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \end{equation}

is the momentum flux,

(2.16)\begin{equation} \boldsymbol{\varPhi}_{e,{\rm ext}}=\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\sum_{i=1}^{N/2} \left(q_{i}\boldsymbol{E}_{\rm ext}(\boldsymbol{r}_{i}(t),t) +\frac{1}{c}q_{i}[\boldsymbol{v}_{i}(t), \boldsymbol{B}_{\rm ext}(\boldsymbol{r}_{i}(t),t)]\right) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \end{equation}

is the density of the force caused by the action of the external fields,

(2.17)\begin{equation} \boldsymbol{\varPhi}_{e}={-}\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\sum_{i=1}^{N/2} \sum_{j=1,j\neq i}^{N}q_{i}q_{j}(\nabla_{i}G_{ij}) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \end{equation}

is the density of the force caused by the interparticle interaction, $a$ and other Latin indexes (from the beginning of the alphabet) correspond to the vector indexes in the Euclidian space, and the Einstein rule of the summation on the repeating indexes is assumed: $\boldsymbol {j} \boldsymbol {\cdot }\boldsymbol {E}=\sum _{a} j_{a} E^{a}=j_{a} E^{a}=j^{a} E^{a}$.

2.1. Multipole moments of physically infinitesimal volume in the external force field

We consider the external force field as two parts, one related to the electric field

(2.18)\begin{equation} \boldsymbol{\varPhi}_{\rm ext,el}=q_{s} \frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\sum_{i=1}^{N/2} \boldsymbol{E}_{\rm ext}(\boldsymbol{r}_{i}(t),t) \delta_{i}, \end{equation}

and the second related to the magnetic field

(2.19)\begin{equation} \boldsymbol{\varPhi}_{{\rm ext},m} =\frac{q_{s}}{c}\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\sum_{i=1}^{N/2} [\boldsymbol{v}_{i}(t), \boldsymbol{B}_{\rm ext}(\boldsymbol{r}_{i}(t),t)] \delta_{i}, \end{equation}

where $\delta _{i}\equiv \delta (\boldsymbol {r}+ \boldsymbol {\xi }-\boldsymbol {r}_{i}(t))$ is the short notation.

We start our discussion with the electric part of the external force field. We use the delta-function under the integral to represent the argument of the electric field

(2.20)\begin{equation} \boldsymbol{\varPhi}_{\rm ext,el}=q_{s}\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\sum_{i=1}^{N/2} \boldsymbol{E}_{\rm ext}(\boldsymbol{r}+ \boldsymbol{\xi},t) \delta_{i}, \end{equation}

where we cannot place the electric field outside of the integral. However, if the electric field changes slowly over the physically infinitesimal volume (over the $\varDelta$-neighbourhood), we can expand function $\boldsymbol {E}_{\rm ext}(\boldsymbol {r}+ \boldsymbol {\xi },t)$ on the vector $\boldsymbol {\xi }$ scanning the $\varDelta$-neighbourhood. The partial sum formed by the first three terms of a Taylor series is considered:

(2.21)\begin{equation} \boldsymbol{E}_{\rm ext}(\boldsymbol{r}+ \boldsymbol{\xi},t) \approx \boldsymbol{E}_{\rm ext}(\boldsymbol{r},t) +( \boldsymbol{\xi}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{E}_{\rm ext}(\boldsymbol{r},t) +\tfrac{1}{2}( \boldsymbol{\xi}\boldsymbol{\cdot}\boldsymbol{\nabla})^{2} \boldsymbol{E}_{\rm ext}(\boldsymbol{r},t)+\cdots . \end{equation}

The substitution of this expression in the electric part of the external force field (2.20) gives a corresponding expression of the electric part of the external force field:

(2.22)\begin{equation} \boldsymbol{\varPhi}_{\rm ext,el}=q_{s} n_{s}\boldsymbol{E}_{\rm ext}(\boldsymbol{r},t) +q_{s} (\boldsymbol{d}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{E}_{\rm ext}(\boldsymbol{r},t) +q_{s}Q^{ab}\partial_{a}\partial_{b}\boldsymbol{E}_{\rm ext}(\boldsymbol{r},t) +\cdots, \end{equation}

where

(2.23)\begin{equation} \boldsymbol{d}\equiv \boldsymbol{d}(\boldsymbol{r},t) =\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\sum_{i=1}^{N/2} \boldsymbol{\xi}\delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \end{equation}

is the electric dipole moment of the $\varDelta$-neighbourhood divided by the charge $q_{s}$, and

(2.24)\begin{equation} Q^{ab}(\boldsymbol{r},t)= \frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\sum_{i=1}^{N/2} \xi^{a}\xi^{b}\delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \end{equation}

is the electric quadrupole moment of the $\varDelta$-neighbourhood divided by the charge $q_{s}$.

2.2. Multipole moments of the Lorentz force field

Let us consider the magnetic part of the force field (2.26), where we expand the magnetic field on the vector $\boldsymbol {\xi }$ scanning the $\varDelta$-neighbourhood:

(2.25)\begin{equation} \boldsymbol{B}_{\rm ext}(\boldsymbol{r}_{i}(t),t)=\boldsymbol{B}_{\rm ext}(\boldsymbol{r}+ \boldsymbol{\xi},t) \approx \boldsymbol{B}_{\rm ext}(\boldsymbol{r},t) +( \boldsymbol{\xi}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{B}_{\rm ext}(\boldsymbol{r},t) +\tfrac{1}{2}( \boldsymbol{\xi}\boldsymbol{\cdot}\boldsymbol{\nabla})^{2} \boldsymbol{B}_{\rm ext}(\boldsymbol{r},t)+\cdots . \end{equation}

Therefore, (2.26) can be rewritten as

(2.26)\begin{equation} \varPhi_{{\rm ext},m}^{a} =\frac{q_{s}}{c}\varepsilon^{abc}(j^{b} B_{\rm ext}^{c}(\boldsymbol{r},t) + J_{D}^{bd}\partial^{d}B_{\rm ext}^{c}(\boldsymbol{r},t) +J_{Q}^{bdf} \partial^{d}\partial^{f}B_{\rm ext}^{c}(\boldsymbol{r},t))+\cdots, \end{equation}

where we also use the Levi–Civita symbol $\varepsilon ^{abc}$ for the vector product in the tensor notation, where

(2.27)\begin{equation} J_{D}^{ab}(\boldsymbol{r},t) =\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\sum_{i=1}^{N/2} v_{i}^{a}(t)\xi^{b}\delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \end{equation}

is the flux of the electric dipole moment of the $\varDelta$-neighbourhood divided by the charge $q_{s}$, and

(2.28)\begin{equation} J_{Q}^{abc}(\boldsymbol{r},t)=\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\sum_{i=1}^{N/2} v_{i}^{a}(t)\xi^{b}\xi^{c}\delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \end{equation}

is the flux of the electric quadrupole moment of the $\varDelta$-neighbourhood divided by the charge $q_{s}$.

3. Self-consistent field approximation in non-relativistic hydrodynamics

In the previous section, we presented the derivation of the general form of the Euler equation. Moreover, we considered the multipole expansion of the external force field. Here, we consider the interparticle interaction. We have two goals to achieve in this section. The first goal is the analysis of the self-consistent field approximation to understand its properties in the deterministic derivation of the hydrodynamic equations. The second goal is the multipole expansion of the interparticle interaction force field.

We repeat (2.17) to discuss it in more detail. Moreover, it is useful to specify the number of species in the system. To get the most simple presentation, we chose the electron–proton plasmas (or completely ionized hydrogen plasmas). Let us numerate electrons as the particles with numbers $i\in [1,N/2]$ and ions as the particles with numbers $i\in [N/2+1,N]$. Let us also point out that a set of hydrodynamic equations is obtained for each species. We focus on the dynamics of electrons. Therefore, the force field acting on electrons is composed of the electron–electron interaction, and the force field created by ions and acting on the electrons. The force field created by electrons, which acts on the electrons (the self-action of the electron material field), has the following form:

(3.1)\begin{equation} \boldsymbol{\varPhi}_{e-e}={-}\frac{q_{e}^{2}}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\sum_{i,j=1,j\neq i}^{N/2} \nabla_{i}G (|\boldsymbol{r}_{i}(t)-\boldsymbol{r}_{j}(t)|) \delta_{i}. \end{equation}

The action of ions on the electrons can be written in the following form:

(3.2)\begin{equation} \boldsymbol{\varPhi}_{e-i}={-}q_{e}q_{i}\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi} \sum_{i=1}^{N/2}\sum_{j=N/2+1}^{N} \nabla_{i}G (|\boldsymbol{r}_{i}(t)-\boldsymbol{r}_{j}(t)|) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)). \end{equation}

We continue our analysis for the expression (3.1). This force field is not symmetric relative to the $i$th and $j$th particles. Introducing the integral over the whole space, we include the delta function containing the coordinate of the $j$th particle

(3.3)\begin{equation} \boldsymbol{\varPhi}_{e-e}={-}\frac{q_{e}^{2}}{\varDelta}\int {\rm d} \boldsymbol{r}'\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\sum_{i,j=1,j\neq i}^{N/2} \nabla_{i}G (|\boldsymbol{r}_{i}(t)-\boldsymbol{r}'|) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{r}'-\boldsymbol{r}_{j}(t)), \end{equation}

where we consider the integral over the whole space on $\boldsymbol {r}'$. Here we see that the delta functions containing the $i$th and $j$th particles have a different structure of arguments. We need to continue the symmetrization of the force field. To this end, we use the following mathematical relation: if we have two functions with the relation $f(\boldsymbol {r})=(1/\varDelta ) \int _{\varDelta }g(\boldsymbol {r}+ \boldsymbol {\xi }) \, {\rm d} \boldsymbol {\xi }$, we find that their integrals over the whole space are equal to each other $\int {\rm d} \boldsymbol {r} f(\boldsymbol {r})=\int {\rm d} \boldsymbol {r} g(\boldsymbol {r})$. Consequently, we can represent integral $\int {\rm d} \boldsymbol {r}'\nabla _{i}G (|\boldsymbol {r}_{i}(t)-\boldsymbol {r}'|) \delta (\boldsymbol {r}'-\boldsymbol {r}_{j}(t))$ as the following structure: $(1/\varDelta )\int {\rm d} \boldsymbol {r}' \int _{\varDelta }\, {\rm d} \boldsymbol {\xi } \nabla _{i}G (|\boldsymbol {r}_{i}(t)-\boldsymbol {r}'- \boldsymbol {\xi }|) \delta (\boldsymbol {r}'+ \boldsymbol {\xi }-\boldsymbol {r}_{j}(t))$. It leads to the symmetric form of

(3.4)\begin{align} \boldsymbol{\varPhi}_{e-e}& ={-}\frac{q_{e}^{2}}{\varDelta^{2}} \int {\rm d} \boldsymbol{r}'\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}' \sum_{i,j=1,j\neq i}^{N/2}\nonumber\\ & \quad \times\nabla_{i}G (|\boldsymbol{r}_{i}(t)-\boldsymbol{r}_{j}(t)|) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{r}'+ \boldsymbol{\xi}'-\boldsymbol{r}_{j}(t)). \end{align}

3.1. Monopole approximation and the self-consistent field approximation

We can use the delta function to express the coordinates of particles in the Green function on $\boldsymbol {r}+ \boldsymbol {\xi }$ and $\boldsymbol {r}'+ \boldsymbol {\xi }'$:

(3.5)\begin{align} \boldsymbol{\varPhi}_{e-e}& ={-}\frac{q_{e}^{2}}{\varDelta^{2}} \int {\rm d} \boldsymbol{r}'\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}' \sum_{i,j=1,j\neq i}^{N/2}\nonumber\\ & \quad \times\nabla_{\boldsymbol{r}}G (|\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}'- \boldsymbol{\xi}'|) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{r}'+ \boldsymbol{\xi}'-\boldsymbol{r}_{j}(t)). \end{align}

Both (3.4) and (3.5) show that we cannot introduce the two-particle number density at this stage of the derivation. This is because we cannot extract the Green function of the Coulomb interaction $G$ from under integrals over $\boldsymbol {\xi }$ and $\boldsymbol {\xi }'$.

First, we consider the multipole expansion of the Green function $G (|\boldsymbol {r}+ \boldsymbol {\xi }-\boldsymbol {r}'- \boldsymbol {\xi }'|)$ assuming that it slowly changes on the scale of the $\varDelta$-neighbourhood. Let us specify the assumptions for this expansion. We have two $\varDelta$-neighbourhoods in (3.5); one is around point $\boldsymbol {r}$ and the second one around point $\boldsymbol {r}'$. We can consider two regimes. The first regime corresponds to points $\boldsymbol {r}$ and $\boldsymbol {r}'$ far from each other, so these neighbourhoods do not overlap each other. The second regime corresponds to points $\boldsymbol {r}$ and $\boldsymbol {r}'$ relatively close to each other, so these neighbourhoods overlap. Let us consider the Green function $G (|\boldsymbol {r}-\boldsymbol {r}'+ \boldsymbol {\xi }- \boldsymbol {\xi }'|)$ in the first regime, where $|\boldsymbol {r}-\boldsymbol {r}'| >\varDelta ^{1/3}$. We have $|\boldsymbol {\xi }|$, $|\boldsymbol {\xi }|'\leq \varDelta ^{1/3}$ and $|\boldsymbol {\xi }- \boldsymbol {\xi }'| <|\boldsymbol {\xi }|+|\boldsymbol {\xi }'| \leq 2\varDelta ^{1/3}$. The Green function is basically proportional to $1/R$. In the first regime, we have small variation of $R\in ((ar_{De})^{1/2},r_{De})$ on the scale of $|\boldsymbol {\xi }- \boldsymbol {\xi }'| \sim \varDelta ^{1/3}=(ar_{De})^{1/2}$. In this case, we have a relatively small value of $\boldsymbol {\xi }- \boldsymbol {\xi }'$ in comparison with $\boldsymbol {r}-\boldsymbol {r}'$. It justifies the expansion of the Green function on $\boldsymbol {\xi }- \boldsymbol {\xi }'$ by accounting for a few of the first terms in the expansion. In the second regime, we find $|\boldsymbol {r}-\boldsymbol {r}'| <|\boldsymbol {r}|+|\boldsymbol {r}'| \leq 2\varDelta ^{1/3}$ and $|\boldsymbol {\xi }- \boldsymbol {\xi }'| <2\varDelta ^{1/3}$. So, we can make the formal expansion on $\boldsymbol {\xi }- \boldsymbol {\xi }'$, but we need to keep an infinite number of terms of the expansion. However, we can specify a limit of the second regime. This limit appears at $|\boldsymbol {r}-\boldsymbol {r}'|\ll \varDelta ^{1/3}$. In this case, there is large overlap of two neighbourhoods. Interaction of particles in the area of the overlapping cancel each other due to Newton's third law. However, there is the interaction of particles in the area of the overlap with the particles in areas beyond the overlap, as well as the interaction of particles in areas beyond the overlap. The described overlapping decreases the role of interaction for the close points $\boldsymbol {r}$ and $\boldsymbol {r}'$.

For simplicity, in this subsection, we consider the zero-order expansion (the monopole limit). So, we have $G (|\boldsymbol {r}+ \boldsymbol {\xi }-\boldsymbol {r}'- \boldsymbol {\xi }|) \approx G (|\boldsymbol {r}-\boldsymbol {r}'|)$. Hence, the Green function can be placed outside of the integral on the $\varDelta$-neighbourhood:

(3.6)\begin{equation} \boldsymbol{\varPhi}_{e-e}={-}q_{e}^{2} \int {\rm d} \boldsymbol{r}' \nabla_{\boldsymbol{r}}G (|\boldsymbol{r}-\boldsymbol{r}'|) \boldsymbol{\cdot} n_{2,ee}(\boldsymbol{r},\boldsymbol{r}',t), \end{equation}

where

(3.7)\begin{equation} n_{2,ee}(\boldsymbol{r},\boldsymbol{r}',t)= \frac{1}{\varDelta^{2}} \int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}' \sum_{i,j=1,j\neq i}^{N/2} \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{r}'+ \boldsymbol{\xi}'-\boldsymbol{r}_{j}(t)) \end{equation}

is the two-particle number density.

The Debye radius $r_{D}$ is the distance where the Coulomb field of the charge is screened. The screening is a macroscopic effect which requires a macroscopic number of particles in the Debye sphere. The average interparticle distance is $a\equiv n^{-1/3}$ and we have $r_{D}\gg a$.

If we consider the neutral particles, we have a strong decrease of the potential of interaction. Hence, the considerable changes in the state of motion of neutral particles are interpreted as the collisions since it happens at the small aiming parameter. The neutral particles move as the free particles between collisions.

In plasmas, we have the long-range interaction. The interaction is screened, but it happens on the scale of the Debye radius $r_{D}$. However, the interaction of the charged particles inside the Debye sphere is not interpreted as collisions, at least, not all of these interactions have an interpretation as these collisions, but only a small part of them.

There is a mechanism of chaotic interactions which transits the system to the equilibrium state and leads to the increase of the entropy. It corresponds to the interaction at small interparticle distances. It can be interpreted as the scattering with a small aiming parameter. The interaction at large distances do not lead to the relaxation (see also Landau & Lifshitz Reference Landau and Lifshitz1980).

Most of the time, particles are at distances of the order of the average interparticle distance $a=n^{-1/3}$, but particles with large velocities can come close to relatively small distances $\Delta r\ll a$ overcoming the Coulomb repulsion. Particles mostly are located at distances corresponding to the average interparticle distances from their neighbours. The particles move in the average collective field. However, incident convergence to distances $\Delta r\ll a$ gives strong scattering. The described scales allow to introduce a corresponding scaling of the electromagnetic field acting on the $i$th particle located at point $\boldsymbol {r}_{i}(t')=\boldsymbol {r}'$. For instance, let us present the decomposition of the electromagnetic field vector $\boldsymbol {E}_{i}(\boldsymbol {r}_{i}(t'),t')= \boldsymbol {E}_{>\varDelta }(\boldsymbol {r}',t')+ \boldsymbol {e}_{<\varDelta }(\boldsymbol {r}',t')$. Here, vector $\boldsymbol {E}_{>\varDelta }$ is the electric field created by the particles beyond the $\varDelta$-neighbourhood surrounding the $i$th particle (see figure 2). Vector $\boldsymbol {e}_{<\varDelta }$ is the electric field created by the particles inside the $\varDelta$-neighbourhood surrounding the $i$th particle.

Figure 2.

Illustration of the $\varDelta$-neighbourhood around an arbitrary $i$th particle.

However, if we consider the $\varDelta$-neighbourhood around the arbitrary point of space $\boldsymbol {r}$ and consider the evolution of particles inside (getting in or out) the $\varDelta$-vicinity, we have a distribution of particles in different points of the neighbourhood (not in its centre), as is demonstrated in figure 3. So, we have the following picture for the arbitrary particle in the $\varDelta$-neighbourhood.

Figure 3.

Illustration of the $\varDelta$-neighbourhood around an arbitrary point of space $\boldsymbol {r}$, while some particles are around the centre.

Ratio $(a/r_{D})^{1/2}$ is the small parameter $\epsilon =(a/r_{D})^{1/2}$. This parameter allows us to introduce the intermediate scale $\varDelta ^{1/3}$. It leads to the following explicit expression for the radius of the $\varDelta$-neighbourhood: $\varDelta ^{1/3}=\sqrt {a \cdot r_{D}}$. There is no specific value of $\varDelta$, but we estimate the minimal value of $\varDelta$. So we find a spatial scale, which is large from a microscopic point of view associated with the average interparticle distance, and small from the macroscopic point of view associated with the Debye length.

Figure 4 shows the $i$th picture which belongs to the $\varDelta _{\boldsymbol {r}}$-neighbourhood of point $\boldsymbol {r}$ and to the $\varDelta _{\boldsymbol {r}'}$-neighbourhood of point $\boldsymbol {r}'$. Hence, the particle $i$ is under action of ‘collisions’ from the particle $j\in \varDelta _{\boldsymbol {r}'}$ (but beyond $\varDelta _{\boldsymbol {r}}$, see figure 5). Figure 5 shows particles $j_{k}\in \varDelta _{\boldsymbol {r}'}$, but they do not belong to $\varDelta _{\boldsymbol {r}}$, which ‘collide’ with particle $i$ to change the momentum of particles in the neighbourhood $\varDelta _{\boldsymbol {r}}$. Figure 6 shows the $i$th and $j$th particles which simultaneously belong to $\varDelta _{\boldsymbol {r}}$ and $\varDelta _{\boldsymbol {r}'}$. Their interaction does not change the momentum of all particles in $\varDelta _{\boldsymbol {r}}$ due to Newton's third law. To neglect the contribution of the collisions of the $i$th and $j$th particles illustrated in figure 6 in the evolution of the particles in the $\varDelta _{\boldsymbol {r}}$-neighbourhood, we need to keep $\boldsymbol {r}'$ at distances larger than $2\varDelta ^{1/3}$ from point $\boldsymbol {r}$. Therefore, the vicinities do not cross each other. Particles $i_{k}\in \varDelta _{\boldsymbol {r}}$ interact with particles $j_{k}\in \varDelta _{\boldsymbol {r}'}$ up to distances $|\boldsymbol {r}-\boldsymbol {r}'|\sim r_{De}$. The interaction can be neglected completely for the larger distances, as it follows from the well-known Debye screening effect (see Akhiezer Reference Akhiezer1975; Landau & Lifshitz Reference Landau and Lifshitz1980; Aleksandrov et al. Reference Aleksandrov, Bogdankevich and Rukhadze1984).

Figure 4.

Interaction is represented via the two-particle functions, which includes the consideration of delta vicinities of two arbitrary points $\boldsymbol {r}$ and $\boldsymbol {r}'$. The regime of overlapping delta vicinities of points $\boldsymbol {r}$ and $\boldsymbol {r}'$ is illustrated. Positions of the $i$th particles belonging to both vicinities are illustrated as well.

Figure 5.

Two groups of interacting particles. One group is illustrated via the single particle $i$ belonging to both neighbourhoods. The second group of particles is illustrated by $j_{1}$, $j_{2}$ and $j_{3}$, which belong to the $\varDelta$-neighbourhood of point $\boldsymbol {r}'$.

Figure 6.

Illustration of two interacting particles simultaneously being parts of two $\varDelta$-neighbourhoods.

We use the self-consistent field approximation in (3.6), so we assume $n_{2}(\boldsymbol {r},\boldsymbol {r}',t)=n(\boldsymbol {r},t)n(\boldsymbol {r}',t)$, and find

(3.8)\begin{equation} \boldsymbol{\varPhi}_{e-e}={-}q_{e}^{2}n(\boldsymbol{r},t) \nabla_{\boldsymbol{r}}\int {\rm d} \boldsymbol{r}' G (|\boldsymbol{r}-\boldsymbol{r}'|) n(\boldsymbol{r}',t). \end{equation}

Equation (3.8) allows us to introduce the electrostatic potential of the electric field created by electrons as $\varphi _{e}(\boldsymbol {r},t)=q_{e}\int {\rm d} \boldsymbol {r}' G (|\boldsymbol {r}-\boldsymbol {r}'|) n(\boldsymbol {r}',t)$, and the corresponding electric field $\boldsymbol {E}_{e}=-\nabla _{\boldsymbol {r}}\varphi _{e}(\boldsymbol {r},t)$. The obtained electric field satisfies the following equations: $\boldsymbol {\nabla } \times \boldsymbol {E}_{e}=0$, and $\boldsymbol {\nabla } \boldsymbol {\cdot }\boldsymbol {E}_{e}=4{\rm \pi} q_{e}n_{e}$. The complete electric field is the superposition of the electric fields created by all species of the system $\boldsymbol {E}=\sum _{s=e,i}\boldsymbol {E}_{s}$, which obeys the electrostatic Maxwell equations $\boldsymbol {\nabla } \times \boldsymbol {E}=0$ and

(3.9)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{E}=4{\rm \pi} \sum_{s} q_{s}n_{s}. \end{equation}

At this stage, we can present the intermediate form of the Euler equation (2.14):

(3.10)\begin{equation} \partial_{t}j_{s}^{a}+\partial_{b}\varPi_{s}^{ab}= \frac{q_{s}}{m_{s}}(n_{s}(E_{\rm ext}^{a}+E^{a})+ \varepsilon^{abc}j_{s}^{b}B_{\rm ext}^{c}). \end{equation}

3.2. Multipole expansion

Here, we consider the multipole expansion for the interaction field. Formally, we can consider the expansion of the Green function $G (|\boldsymbol {r}+ \boldsymbol {\xi }-\boldsymbol {r}'- \boldsymbol {\xi }'|)$ if it slowly changes on the scale of the $\varDelta$-neighbourhood. We need to consider (3.5) in more detail. First, we present the expansion of the Green function

(3.11)\begin{align} & G (|\boldsymbol{r}-\boldsymbol{r}'+ \boldsymbol{\xi}- \boldsymbol{\xi}'|) \approx G (|\boldsymbol{r}-\boldsymbol{r}'|) + (\xi^{a}-\xi'^{a})\partial_{a} G (|\boldsymbol{r}-\boldsymbol{r}'|)\nonumber\\ & \quad +\tfrac{1}{2}(\xi^{a}-\xi'^{a})(\xi^{b}-\xi'^{b}) \partial_{a} \partial_{b} G (|\boldsymbol{r}-\boldsymbol{r}'|). +\cdots . \end{align}

The direct substitution in (3.5) gives the following expression:

(3.12)\begin{align} \varPhi^{a}_{e-e}& ={-}\frac{q_{e}^{2}}{\varDelta^{2}} \int {\rm d} \boldsymbol{r}' \partial_{\boldsymbol{r}}^{a}G (|\boldsymbol{r}-\boldsymbol{r}'|) \int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}' \sum_{i,j=1,j\neq i}^{N/2} \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{r}'+ \boldsymbol{\xi}'-\boldsymbol{r}_{j}(t))\nonumber\\ & \quad -\frac{q_{e}^{2}}{\varDelta^{2}} \int {\rm d} \boldsymbol{r}' \partial_{\boldsymbol{r}}^{a}\partial_{\boldsymbol{r}}^{b}G (|\boldsymbol{r}-\boldsymbol{r}'|) \int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\nonumber\\& \quad \int_{\varDelta}\, {\rm d} \boldsymbol{\xi}' \sum_{i,j=1,j\neq i}^{N/2} (\xi^{a}-\xi'^{a}) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{r}'+ \boldsymbol{\xi}'-\boldsymbol{r}_{j}(t))\nonumber\\ & \quad-\frac{q_{e}^{2}}{\varDelta^{2}} \int {\rm d} \boldsymbol{r}' \partial_{\boldsymbol{r}}^{a}\partial_{\boldsymbol{r}}^{b}\partial_{\boldsymbol{r}}^{c}G (|\boldsymbol{r}-\boldsymbol{r}'|) \int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\nonumber\\& \quad \int_{\varDelta}\, {\rm d} \boldsymbol{\xi}' \sum_{i,j=1,j\neq i}^{N/2} \frac{1}{2}(\xi^{a}\xi^{b}-\xi'^{a}\xi^{b}-\xi^{a}\xi'^{b}+\xi'^{a}\xi'^{b}) \delta_{\boldsymbol{r}i}\delta_{\boldsymbol{r}'j} +\cdots, \end{align}

where $\delta _{\boldsymbol {r}i}=\delta (\boldsymbol {r}+ \boldsymbol {\xi }-\boldsymbol {r}_{i}(t))$ and $\delta _{\boldsymbol {r}'j}=\delta (\boldsymbol {r}'+ \boldsymbol {\xi }'-\boldsymbol {r}_{j}(t))$.

The obtained expression can be rewritten via the corresponding two-particle functions

(3.13)\begin{align} \varPhi^{a}_{e-e}& ={-}q_{e}^{2} \int {\rm d} \boldsymbol{r}' \partial_{\boldsymbol{r}}^{a}G (|\boldsymbol{r}-\boldsymbol{r}'|) \boldsymbol{\cdot} n_{2}(\boldsymbol{r},\boldsymbol{r}',t)\nonumber\\& \quad -q_{e}^{2} \int {\rm d} \boldsymbol{r}' \partial_{\boldsymbol{r}}^{a}\partial_{\boldsymbol{r}}^{b}G (|\boldsymbol{r}-\boldsymbol{r}'|) \boldsymbol{\cdot} (d_{2}^{b}(\boldsymbol{r},\boldsymbol{r}',t)-d_{2}^{b}(\boldsymbol{r}', \boldsymbol{r},t))\nonumber\\ & \quad -\frac{1}{2}q_{e}^{2} \int {\rm d} \boldsymbol{r}'\partial_{\boldsymbol{r}}^{a} \partial_{\boldsymbol{r}}^{b}\partial_{\boldsymbol{r}}^{c}G (|\boldsymbol{r}-\boldsymbol{r}'|)\boldsymbol{\cdot}(Q_{2}^{bc}(\boldsymbol{r},\boldsymbol{r}',t) \nonumber\\ & \quad +Q_{2}^{bc}(\boldsymbol{r}',\boldsymbol{r},t) -D_{2}^{bc}(\boldsymbol{r},\boldsymbol{r}',t)-D_{2}^{cb} (\boldsymbol{r},\boldsymbol{r}',t)), \end{align}

where we introduce three two-particle functions

(3.14)\begin{equation} d_{2}^{b}(\boldsymbol{r},\boldsymbol{r}',t)= \frac{1}{\varDelta^{2}}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}' \sum_{i=1,j\neq i}^{N/2}\xi^{b} \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{r}'+ \boldsymbol{\xi}'-\boldsymbol{r}_{j}(t)), \end{equation}

where the permutation of its arguments leads to

(3.15)\begin{equation} d_{2}^{b}(\boldsymbol{r}',\boldsymbol{r},t)= \frac{1}{\varDelta^{2}}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}' \sum_{i=1,j\neq i}^{N/2}\xi'^{b} \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{r}'+ \boldsymbol{\xi}'-\boldsymbol{r}_{j}(t)), \end{equation}

which are the two forms of the two-particle function of number density-polarization (3.14) or polarization-number density (3.15),

(3.16)\begin{equation} Q_{2}^{bc}(\boldsymbol{r},\boldsymbol{r}',t)= \frac{1}{\varDelta^{2}}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}' \sum_{i=1,j\neq i}^{N/2}\xi^{b}\xi^{c} \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{r}'+ \boldsymbol{\xi}'-\boldsymbol{r}_{j}(t)), \end{equation}

which is the two-particle function of number density-quadrupole moment, and

(3.17)\begin{equation} D_{2}^{bc}(\boldsymbol{r},\boldsymbol{r}',t)= \frac{1}{\varDelta^{2}}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}' \sum_{i=1,j\neq i}^{N/2}\xi^{b}\xi'^{c} \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{r}'+ \boldsymbol{\xi}'-\boldsymbol{r}_{j}(t)), \end{equation}

which is the two-particle polarization–polarization function.

3.3. Self-consistent field approximation in the multipole regime

Let us repeat the conclusion about the meaning of the self-consistent field approximation. Our discussion presented above leads to conclusion that the $\varDelta$-neighbourhood has the minimal radius of the order of $\sqrt {a \cdot r_{D}}$, where $a=n^{-1/3}$ is the average interparticle distance and $r_{D}$ is the Debye radius. The self-consistent field approximation corresponds to the regime of interaction of particles being in the non-overlapping $\varDelta$-neighbourhoods. This condition allows to split the two-particle hydrodynamic functions to the product of corresponding one-particle hydrodynamic functions.

We use the self-consistent field approximation in (3.13), so we assume $n_{2}(\boldsymbol {r},\boldsymbol {r}',t)=n(\boldsymbol {r},t)n(\boldsymbol {r}',t)$, $d_{2}^{b}(\boldsymbol {r},\boldsymbol {r}',t)=d^{b}(\boldsymbol {r},t) n(\boldsymbol {r}',t)$, $d_{2}^{b}(\boldsymbol {r}',\boldsymbol {r},t)=n(\boldsymbol {r},t) d^{b}(\boldsymbol {r}',t)$, $Q_{2}^{bc}(\boldsymbol {r},\boldsymbol {r}',t)=Q^{bc}(\boldsymbol {r},t)n(\boldsymbol {r}',t)$, $Q_{2}^{bc}(\boldsymbol {r}',\boldsymbol {r},t)=n(\boldsymbol {r},t) Q^{bc}(\boldsymbol {r}',t)$, and $D_{2}^{bc}(\boldsymbol {r},\boldsymbol {r}',t)=d^{b}(\boldsymbol {r},t)d^{c}(\boldsymbol {r}',t)$. Consequently, (3.13) transforms into

(3.18)\begin{align} \varPhi^{a}_{e-e}& ={-}q_{e}^{2}\left[ n(\boldsymbol{r},t)\partial^{a}\left(\int {\rm d} \boldsymbol{r}' G (|\boldsymbol{r}-\boldsymbol{r}'|) n(\boldsymbol{r}',t) -\partial^{b}\int {\rm d} \boldsymbol{r}' G (|\boldsymbol{r}-\boldsymbol{r}'|) d^{b}(\boldsymbol{r}',t)\right.\right.\nonumber\\ & \quad \left.+\,\frac{1}{2}\partial^{b}\partial^{c}\int {\rm d} \boldsymbol{r}' G (|\boldsymbol{r}-\boldsymbol{r}'|) n(\boldsymbol{r}',t)+\cdots\right) +d^{b}(\boldsymbol{r},t)\partial^{a}\partial^{b}\left(\int {\rm d} \boldsymbol{r}' G (|\boldsymbol{r}-\boldsymbol{r}'|) n(\boldsymbol{r}',t)\right.\nonumber\\ & \quad\left.-\,\partial^{c}\int {\rm d} \boldsymbol{r}' G (|\boldsymbol{r}-\boldsymbol{r}'|) d^{c}(\boldsymbol{r}',t)+\cdots\right)\nonumber\\ & \quad\left.+\,Q^{bc}(\boldsymbol{r},t)\partial^{a}\partial^{b}\left(\int {\rm d} \boldsymbol{r}' G (|\boldsymbol{r}-\boldsymbol{r}'|) n(\boldsymbol{r}',t) +\cdots\right) +\cdots\right]. \end{align}

Equation (3.18) allows to introduce the multipole expansion of the electrostatic potential

(3.19)\begin{align} \varphi_{e}& =q_{e}\left(\int {\rm d} \boldsymbol{r}' G (\boldsymbol{r},\boldsymbol{r}') n(\boldsymbol{r}',t) -\partial^{b}\int {\rm d} \boldsymbol{r}' G (\boldsymbol{r},\boldsymbol{r}') d^{b}(\boldsymbol{r}',t) \right. \nonumber\\ & \left. \quad +\, \frac{1}{2}\partial^{b}\partial^{c}\int {\rm d} \boldsymbol{r}' G (\boldsymbol{r},\boldsymbol{r}') n(\boldsymbol{r}',t)+\cdots\right). \end{align}

Therefore, (3.18) can be represented in the structural form

(3.20)\begin{equation} \boldsymbol{\varPhi}_{e-e}=q_{e}n_{e}\boldsymbol{E}_{e} +(\boldsymbol{d}(\boldsymbol{r},t)\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{E}_{e}+ \tfrac{1}{2}Q^{bc}(\boldsymbol{r},t)\partial^{b}\partial^{c}\boldsymbol{E}_{e}+\cdots, \end{equation}

where $\boldsymbol {E}_{e}=-\boldsymbol {\nabla }\varphi _{e}(\boldsymbol {r},t)$, and we also can introduce full electric field $\boldsymbol {E}=\sum _{s=e,i}\boldsymbol {E}_{s}$, which satisfies the following quasi-static Maxwell equations: $\boldsymbol {\nabla } \times \boldsymbol {E}=0$ and

(3.21)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{E}=4{\rm \pi} \sum_{s} q_{s}\left(n_{s} +(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{d}(\boldsymbol{r},t)) +\frac{1}{2}\partial^{b}\partial^{c}Q^{bc}(\boldsymbol{r},t) +\cdots \right). \end{equation}

The full set of hydrodynamic equations requires the equations for the evolution of functions $\boldsymbol {d}(\boldsymbol {r},t)$, $Q^{bc}(\boldsymbol {r},t)$, etc. We do not present or discuss these equations. Some information can be found from Drofa & Kuz'menkov (Reference Drofa and Kuz'menkov1996).

While we present systematic derivation of the contribution of the electric dipole moment in the set of hydrodynamic equations, we do not consider the evolution of the electric dipole moment density. However, it is considered by Drofa & Kuz'menkov (Reference Drofa and Kuz'menkov1996), where a closed set of hydrodynamic equations is developed, and some applications to the waves and solitons are considered. Study of the evolution of the electric dipole moment density based on the quantum methods is suggested by Andreev, Kuz'menkov & Trukhanova (Reference Andreev, Kuz'menkov and Trukhanova2011) and Li, Hwang & Sarma (Reference Li, Hwang and Sarma2011), where some applications to waves in the low-dimensional systems are studied.

4. Derivation of the Vlasov equation tracing the microscopic motion of particles

To derive the kinetic equations, we need to introduce the distribution function in the six-dimensional coordinate-momentum space. We start with the microscopic definition for the system of the point-like particles (Klimontovich Reference Klimontovich1986):

(4.1)\begin{equation} f(\boldsymbol{r},\boldsymbol{p},t)=\sum_{i=1}^{N/2} \delta(\boldsymbol{r}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{p}-\boldsymbol{p}_{i}(t)), \end{equation}

where we have the coordinate of the $i$th particle $\boldsymbol {r}_{i}(t)$ and the momentum of the $i$th particle $\boldsymbol {p}_{i}(t)=m_{i}\dot {\boldsymbol {r}}_{i}(t)$ (the non-relativistic regime). We use the notation $f(\boldsymbol {r},\boldsymbol {p},t)$ for the microscopic distribution function, same as the notation for the macroscopic function below. However, we can give a more detailed representation of arguments of the distribution function (4.1) via tracing the time dependence as follows: $f(\boldsymbol {r},\boldsymbol {p},\boldsymbol {r}_{i}(t),\boldsymbol {p}_{i}(t))$ or $f(\boldsymbol {r},\boldsymbol {p},\{\boldsymbol {r}_{1}(t),\boldsymbol {p}_{1}(t),\ldots, \boldsymbol {r}_{N/2}(t),\boldsymbol {p}_{N/2}(t)\})$. Let us repeat that the consideration of the kinetic model of plasmas requires the introduction of the distribution function for each species. Each distribution function evolves under the action of all species in the system. We specify that we consider the quasi-neutral plasmas of two species: the electron–ion plasmas. We use the following numeration of particles: $i\in [1,N/2]$ for the electrons and $i\in [N/2+1,N]$ for the ions. It can be easily assigned to the arbitrary set of species in plasmas including neutral particles. We illustrate the derivation following the evolution of the distribution function of electrons. The explicit contribution of ions is shown in the terms describing the interaction.

Next, we make the transition to the physically infinitesimal volume. Moreover, we consider the physically infinitesimal areas both for the coordinate space and the momentum space. It is constructed in the similar way as the number density (2.8) presented above,

(4.2)\begin{equation} f(\boldsymbol{r},\boldsymbol{p},t)= \frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \sum_{i=1}^{N/2} \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{p}+\boldsymbol{\eta}-\boldsymbol{p}_{i}(t)), \end{equation}

where ${\rm d}\boldsymbol {\eta }$ is the element of volume in the momentum space, $\int _{\varDelta,\varDelta _{p}} \, {\rm d} \boldsymbol {\xi } \,{\rm d}\boldsymbol {\eta } =\int _{\varDelta }\, {\rm d} \boldsymbol {\xi } \int _{\varDelta _{p}}\,{\rm d}\boldsymbol {\eta }$, with $\int _{\varDelta _{p}}\,{\rm d}\boldsymbol {\eta }$ integral over the $\varDelta _{p}$-neighbourhood in the momentum space, $\varDelta \equiv \varDelta _{r}$ is the $\varDelta$-neighbourhood in the coordinate space. As the notion, it is the same $\varDelta$-neighbourhood which is used above in the derivation of hydrodynamics. However, its value in kinetics should be larger to get the continuous distribution function. Momentum space related to definition (4.2) is shown in Fig. 7.

Figure 7.

Consideration of physical kinetics requires analysis of the six-dimensional phase space. We need to construct the $\varDelta$-neighbourhood in the six-dimensional space. It is basically the $\varDelta$-neighbourhoods in the coordinate space and in the momentum space. The $\varDelta$-neighbourhood in the coordinate space is illustrated in figure 1. The $\varDelta$-neighbourhood in the momentum space is illustrated here. Corresponding notation including illustration of vector $\boldsymbol {\eta }$ scanning the vicinity are shown in this figure.

We need to derive the equations for the evolution of the distribution functions. Hence, we need to use the equations of motion of the particles. Here, we consider the non-relativistic regime for the systems of charged particles. Consequently, we can use (2.13). First, we consider the kinetic equation for the microscopic distribution function (4.1). We take the derivative on time of function (4.1) and obtain the equation for its evolution:

(4.3)\begin{equation} \partial_{t}f + \boldsymbol{\nabla}\boldsymbol{\cdot}\sum_{i=1}^{N/2} \boldsymbol{v}_{i}(t)\delta(\boldsymbol{r}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{p}-\boldsymbol{p}_{i}(t)) +\sum_{i=1}^{N/2} \delta(\boldsymbol{r}-\boldsymbol{r}_{i}(t)) (\dot{\boldsymbol{p}}_{i}(t)\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}}) \delta(\boldsymbol{p}-\boldsymbol{p}_{i}(t))=0. \end{equation}

We can replace $\nabla _{\boldsymbol {p}}$ to put it in front of the whole term, but it would interfere with the following manipulations of this equation. We use the equations of motion for each particle. Using the delta functions, we replace $\boldsymbol {r}_{i}(t)\rightarrow \boldsymbol {r}$ and $\boldsymbol {p}_{i}(t)\rightarrow \boldsymbol {p} (\boldsymbol {v}_{i}(t)\rightarrow \boldsymbol {v}=\boldsymbol {p}/m)$ in the expression for the external field acting on the $i$th particle.

This calculation gives the untruncated microscopic kinetic equation (Klimontovich Reference Klimontovich1986):

(4.4)\begin{align} & \partial_{t}f(\boldsymbol{r},\boldsymbol{p},t) +(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla})f(\boldsymbol{r},\boldsymbol{p},t) +\frac{q_{s}}{m_{s}}((\boldsymbol{E}_{\rm ext}(\boldsymbol{r},t)+\boldsymbol{v}\times \boldsymbol{B}_{\rm ext}(\boldsymbol{r},t)/c) \boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}})f(\boldsymbol{r},\boldsymbol{p},t)\nonumber\\ & \quad-\frac{q_{s}q_{s'}}{m_{s}}\nabla_{\boldsymbol{p}} \boldsymbol{\cdot}\sum_{i=1}^{N/2}\sum_{j=1,j\neq i}^{N} (\nabla_{i}G_{ij})\delta(\boldsymbol{r}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{p}-\boldsymbol{p}_{i}(t))=0. \end{align}

The last term in (4.4) describes the interaction. The presence of the Green function of the Coulomb interaction leads to the fact that it cannot be expressed via the distribution function, but it requires the introduction of the two-particle distribution function.

Let us derive the evolution of the distribution function based on the exact microscopic motion, but the motion is considered on the scale of the $\varDelta$-neighbourhood (4.2). We calculate the time derivative of function (4.2) and find the following intermediate equation:

(4.5)\begin{align} & \partial_{t}f(\boldsymbol{r},\boldsymbol{p},t)+ \boldsymbol{\nabla}\boldsymbol{\cdot}\frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \sum_{i=1}^{N/2} \dot{\boldsymbol{r}}_{i}(t) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{p}+\boldsymbol{\eta}-\boldsymbol{p}_{i}(t))\nonumber\\ & \quad +\nabla_{\boldsymbol{p}}\boldsymbol{\cdot}\frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \sum_{i=1}^{N/2} \dot{\boldsymbol{p}}_{i}(t) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{p}+\boldsymbol{\eta}-\boldsymbol{p}_{i}(t))=0. \end{align}

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

Let us make some simple transformations of (4.5). We use the delta-function depending on the momentum to replace the velocity of the $i$th particle $\boldsymbol {v}_{i}(t)=\dot {\boldsymbol {r}}_{i}(t)$ by $(\boldsymbol {p}+\boldsymbol {\eta })/m_{s}$ ($m_{i}=m_{s}$ for all particles of the species under consideration) in the second term of (4.5). In the third term of (4.5), we use the equation of motion for the $i$th particle (2.13). Next, we replace the coordinates $\boldsymbol {r}_{i}(t)$ (velocities $\boldsymbol {v}_{i}(t)$) in the force acting on the $i$th particle by $\boldsymbol {r}+ \boldsymbol {\xi }$ (by $(\boldsymbol {p}+\boldsymbol {\eta })/m_{s}$). Hence, we obtain the following representation of (4.5):

(4.6)\begin{align} & \partial_{t}f(\boldsymbol{r},\boldsymbol{p},t)+ \frac{1}{m_{s}}\boldsymbol{\nabla}\boldsymbol{\cdot}\frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \sum_{i=1}^{N/2} (\boldsymbol{p}+\boldsymbol{\eta}) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{p}+\boldsymbol{\eta}-\boldsymbol{p}_{i}(t))\nonumber\\ & \quad+\frac{q_{s}}{m_{s}}\frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \sum_{i=1}^{N/2} \left(\vphantom{\left.-\,q_{s'}\nabla_{\boldsymbol{r}}\sum_{j=1,j\neq i}^{N}G(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{j}(t))\right)}\boldsymbol{E}(\boldsymbol{r}+ \boldsymbol{\xi},t) +\frac{1}{m_{s}c}[(\boldsymbol{p}+\boldsymbol{\eta})\times \boldsymbol{B}(\boldsymbol{r}+ \boldsymbol{\xi},t)]\right.\nonumber\\ & \quad \left.-\,q_{s'}\nabla_{\boldsymbol{r}}\sum_{j=1,j\neq i}^{N}G(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{j}(t))\right) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}}\delta(\boldsymbol{p}+\boldsymbol{\eta}-\boldsymbol{p}_{i}(t))=0. \end{align}

The second term in (4.6) splits into two terms. The first of them has the well-known form $\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {\nabla } f_{e}$. However, the second part of this term has a rather unusual structure $({1}/{m_{s}})\boldsymbol {\nabla }\boldsymbol {\cdot }({1}/{\varDelta })({1}/{\varDelta _{p}}) \int _{\varDelta,\varDelta _{p}} \, {\rm d} \boldsymbol {\xi } \,{\rm d}\boldsymbol {\eta } \sum _{i=1}^{N/2} \boldsymbol {\eta } \delta _{\boldsymbol {r}i} \delta _{\boldsymbol {p}i}\equiv ({1}/{m_{s}})\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {f}(\boldsymbol {r},\boldsymbol {p},t)$ related to the deviation of the average momentum of particles in the $\varDelta _{p}$-neighbourhood from the value $\boldsymbol {p}$ being the centre of the neighbourhood, where $\delta _{\boldsymbol {r}i}\equiv \delta (\boldsymbol {r}+ \boldsymbol {\xi }-\boldsymbol {r}_{i}(t))$ and $\delta _{\boldsymbol {p}i}\equiv \delta (\boldsymbol {p}+ \boldsymbol {\eta }-\boldsymbol {p}_{i}(t))$. Therefore, the second term in (4.6) contains the vector distribution function, which has the following explicit form:

(4.7)\begin{equation} \boldsymbol{f}(\boldsymbol{r},\boldsymbol{p},t)= \frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \sum_{i=1}^{N/2} \boldsymbol{\eta} \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{p}+\boldsymbol{\eta}-\boldsymbol{p}_{i}(t)). \end{equation}

The third term in (4.6) contains three terms, where the second of them presents the Lorentz force. The Lorentz force $({1}/{m_{s}c})[(\boldsymbol {p}+\boldsymbol {\eta })\times \boldsymbol {B}(\boldsymbol {r}+ \boldsymbol {\xi },t)]$ splits into two terms due to the deviation of the average momentum $\boldsymbol {p}+\boldsymbol {\eta }$ from $\boldsymbol {p}$. Therefore, we have the contribution of the vector distribution function $\boldsymbol {f}(\boldsymbol {r},\boldsymbol {p},t)$ mentioned above.

In addition to the presence of the vector distribution function $\boldsymbol {f}(\boldsymbol {r},\boldsymbol {p},t)$ in kinetic equation (4.6), we see the necessity to make the multipole expansion of the electric field, the magnetic field and the Green function of the electron–electron interaction. Before we make the expansion of the Green function, we need to give the symmetric form to the term containing this function. Technical steps are the same as we use for the transformation of the hydrodynamic equations above. Equation (4.6) can be rewritten in the following form:

(4.8)\begin{align} & \partial_{t}f+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f+\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{f} +\frac{q_{s}}{m_{s}}\frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta}\nonumber\\ & \quad \sum_{i=1}^{N/2} \left(\boldsymbol{E}(\boldsymbol{r}+ \boldsymbol{\xi},t) +\frac{1}{m_{s}c}[(\boldsymbol{p}+\boldsymbol{\eta})\times \boldsymbol{B}(\boldsymbol{r}+ \boldsymbol{\xi},t)]\right) \delta_{\boldsymbol{r}i} \boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}}\delta_{\boldsymbol{p}i}\nonumber\\ & \quad -\frac{q_{s}}{m_{s}}q_{s'}\nabla_{\boldsymbol{p}}\boldsymbol{\cdot} \frac{1}{\varDelta^{2}}\frac{1}{\varDelta_{p}^{2}} \int {\rm d} \boldsymbol{r}' \, {\rm d} \boldsymbol{p}' \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi}' \,{\rm d}\boldsymbol{\eta}'\nonumber\\ & \quad \sum_{i=1}^{N/2}\sum_{j=1, j\neq i}^{N} \nabla_{\boldsymbol{r}}G(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}'- \boldsymbol{\xi}') \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i} \delta_{\boldsymbol{r}'j} \delta_{\boldsymbol{p}'j} =0, \end{align}

where $\delta _{\boldsymbol {r}'j}\equiv \delta (\boldsymbol {r}'+ \boldsymbol {\xi }'-\boldsymbol {r}_{j}(t))$ and $\delta _{\boldsymbol {p}'j}\equiv \delta (\boldsymbol {p}'+ \boldsymbol {\eta }'-\boldsymbol {p}_{j}(t))$. The fourth and fifth terms describe the action of the external fields and the field of other particles on the $i$th particle (with the further summation on $i$ over all particles). In this form, we cannot include the distribution function in these terms due to the presence of the electric field, the magnetic field and the Green function under the integral.

4.1. Monopole approximation of the kinetic equation

The dependence of the electric and magnetic fields on $\boldsymbol {\xi }$ and dependence of the Green function on $\boldsymbol {\xi }- \boldsymbol {\xi }'$ do not allow to replace these functions outside of integrals on $d \boldsymbol {\xi }$. We need to expand these functions on $\boldsymbol {\xi }$ to introduce the distribution function in these terms and obtain the closed model. If functions $\boldsymbol {E}$, $\boldsymbol {B}$ slowly change on scale of the $\varDelta$-neighbourhood, we can expand these functions on $\boldsymbol {\xi }$. It requires that $\varDelta ^{1/3}\partial _{r}E^{a}\ll E^{a}$ and $\varDelta ^{1/3}\partial ^{b}B^{a}\ll B^{a}$. Or, more accurately, we need $\varDelta ^{(n+1)/3}(\partial _{r})^{n+1}E^{a}\ll \varDelta ^{n/3}(\partial _{r})^{n}E^{a}$ and $\varDelta ^{(n+1)/3}(\partial _{r})^{n+1}B^{a}\ll \varDelta ^{n/3}(\partial _{r})^{n}B^{a}$. Conditions for the Green function of the Coulomb interaction $G$ are discussed after (3.5). In this section, we consider the first terms in these expansions. We call it the monopole approximation.

Monopole approximation of (4.8) in the coordinate space (on $\boldsymbol {\xi }$) of the kinetic equation has the following form:

(4.9)\begin{align} & \partial_{t}f +\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f +\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{f} +\frac{q_{s}}{m_{s}}\left(\boldsymbol{E}(\boldsymbol{r},t) +\frac{1}{c}[\boldsymbol{v}\times \boldsymbol{B}(\boldsymbol{r},t)]\right)\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}}f\nonumber\\ & \quad +\frac{q_{s}}{m_{s}^{2}c}(\nabla_{\boldsymbol{p}} \boldsymbol{\cdot}[\boldsymbol{B}(\boldsymbol{r},t)\times \boldsymbol{f}(\boldsymbol{r},\boldsymbol{p},t)])\nonumber\\ & \quad-\frac{q_{s}}{m_{s}}q_{s'}\nabla_{\boldsymbol{p}}\boldsymbol{\cdot} \int {\rm d} \boldsymbol{r}' \, {\rm d} \boldsymbol{p}' (\nabla_{\boldsymbol{r}}G(\boldsymbol{r}-\boldsymbol{r}')) f_{2}(\boldsymbol{r},\boldsymbol{r}',\boldsymbol{p},\boldsymbol{p}',t) =0, \end{align}

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

(4.10)\begin{equation} f_{2,ee}(\boldsymbol{r},\boldsymbol{r}',\boldsymbol{p},\boldsymbol{p}',t) =\frac{1}{\varDelta^{2}}\frac{1}{\varDelta_{p}^{2}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \, {\rm d} \boldsymbol{\xi}' \,{\rm d}\boldsymbol{\eta}' \sum_{i,j=1, j\neq i}^{N/2} \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i} \delta_{\boldsymbol{r}'j}\delta_{\boldsymbol{p}'j} \end{equation}

is the two-particle electron–electron distribution function.

The third and fifth terms should be dropped if we consider the monopole approximation in the momentum space. Neglecting $\boldsymbol {\eta }$ in comparison with the momentum $\boldsymbol {p}$ corresponds to the neglecting of the vector distribution function in (4.9).

5. Self-consistent field approximation in non-relativistic kinetics

5.1. Self-consistent field approximation in the monopole approximation

Formally, the coordinate part of the $\varDelta$-neighbourhood of the point in six-dimensional phase space is introduced in the same way as it is made for the hydrodynamics. Above, we give an estimation of the radius of the $\varDelta _{r}$-neighbourhood as $\varDelta _{r}^{1/3}=\sqrt {a r_{D}}$, where $a=n^{1/3}$ is the average interparticle distance and $r_{D}$ is the Debye radius. This physical estimation reflects the mathematical requirement for the construction of the continuous functions (like the number density, velocity field, etc.) on the macroscopic scale. The $\varDelta$-neighbourhood of each point $\boldsymbol {r}$ should contain the large number of particles $N(\boldsymbol {r},t)\gg 1$, so the change of this number to one or few particles gives the small change of the hydrodynamic functions. We also have the same requirement for the distribution function $f(\boldsymbol {r},\boldsymbol {p},t)$. However, there are stronger restrictions on the number of particles in the six-dimensional $\varDelta$-neighbourhood $N(\boldsymbol {r},\boldsymbol {p},t)$. Hence, if we have fixed interval of momentum ($\varDelta _{p}$-neighbourhood of point $\boldsymbol {p}$ in the momentum space), we should have $f(\boldsymbol {r},\boldsymbol {p},t)$ continuous in the coordinate space. This property should be satisfied for all points in the momentum space. Hence, for each point in the momentum space, we should have in the coordinate space the same number of particles as for the hydrodynamic description $\tilde {N}$. However, the further summation over all momentum space gives the full number of particles in the coordinate $\varDelta _{r}$-neighbourhood up to $\tilde {N}^{2}$. For the fixed number density, it increases the scale of the $\varDelta _{r}$-neighbourhood in the coordinate space up to $\varDelta _{r\in ph}^{1/3}=r_{D}$, where the subindex $r\in ph$ means that it is the coordinate $\varDelta _{r}$-neighbourhood being part of the six-dimensional $\varDelta$-neighbourhood.

We drop the third and fifth terms in (4.9) to get the full monopole approximations in the coordinate space and the momentum space. The monopole approximation is a necessary step to obtain the two-particle distribution function. To get the self-consistent field approximation, we need to split the two-particle distribution function as a product of the single-particle distribution functions. It corresponds to the interaction of particles from two non-overlapping $\varDelta$-neighbourhoods. Finally, we obtain the Vlasov kinetic equation in the quasi-electrostatic approximation

(5.1)\begin{equation} \partial_{t}f+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f +q_{s}\left(\boldsymbol{E}_{\rm ext}+\boldsymbol{E}+ \frac{1}{c}\boldsymbol{v}\times\boldsymbol{B}_{\rm ext}\right) \boldsymbol{\cdot}\frac{\partial f}{\partial \boldsymbol{p}}=0, \end{equation}

where the electric field is caused by the distribution of charges in the coordinate space: $\boldsymbol {\nabla }\times \boldsymbol {E}=0$ and

(5.2)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{E}=4{\rm \pi} \sum_{s} q_{s}\int f_{s}(\boldsymbol{r},\boldsymbol{p},t)\, {\rm d} \boldsymbol{p}. \end{equation}

There is the question of justification of neglecting several terms at the transition from (4.9) to (5.1). Reasonable justification cannot be made on this step. It would require application of equations for new functions. Systematic derivation of equations for new functions will be made in further work. So generalizations of the Vlasov equation and estimation of its applicability will be made there.

The structure of the equation for the two species electron–ion regime is discussed. To simplify the presentation, we explicitly show the evolution of electrons under the interaction with electrons. An account of ions can be made in the way described above.

5.2. Multipole approximation of the kinetic equation

We use the expansions of the electric field (2.21), the magnetic field (2.25), and the Green function of the Coulomb interaction (3.11) on $\boldsymbol {\xi }$ and $\boldsymbol {\xi }- \boldsymbol {\xi }'$ in the kinetic equation up to the second order on $\boldsymbol {\xi }$. Equation (4.9) appears in the zeroth-order multipole expansion in the coordinate space. Here, we consider the multipole expansion of (4.8) in the coordinate space.

We consider terms up to the second order on $\boldsymbol {\xi }$ or $\boldsymbol {\xi }- \boldsymbol {\xi }'$, and find the following equations presented in terms of several one-particle distribution functions:

(5.3) \begin{align} & \partial_{t}f_{s}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f_{s}+\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{f}_{s} +q_{s}\left(\boldsymbol{E}_{\rm ext}+\frac{1}{c}\boldsymbol{v}\times\boldsymbol{B}_{\rm ext}\right)\boldsymbol{\cdot}\frac{\partial f_{s}}{\partial \boldsymbol{p}}\nonumber\\ & \quad+q_{s}\left(\partial^{b}\boldsymbol{E}_{\rm ext}+\frac{1}{c}\boldsymbol{v}\times\partial^{b}\boldsymbol{B}_{\rm ext}\right) \boldsymbol{\cdot}\frac{\partial d_{s}^{b}(\boldsymbol{r},\boldsymbol{p},t)}{\partial \boldsymbol{p}}\nonumber\\ & \quad+q_{e}\left(\partial^{b}\partial^{c}\boldsymbol{E}_{\rm ext} +\frac{1}{c}\boldsymbol{v}\times\partial^{b}\partial^{c}\boldsymbol{B}_{\rm ext}\right) \boldsymbol{\cdot}\frac{\partial Q_{s}(\boldsymbol{r},\boldsymbol{p},t)}{\partial \boldsymbol{p}}\nonumber\\ & \quad+q_{s}\frac{1}{c}\varepsilon^{abc}B_{\rm ext}^{c}\partial_{a,\boldsymbol{p}} f_{s}^{b} +q_{s}\frac{1}{c}\varepsilon^{abc}\partial^{d}B_{\rm ext}^{c} \partial_{a,\boldsymbol{p}}J_{D,s}^{bd} \nonumber\\ & \quad+q_{s}\frac{1}{c}\varepsilon^{abc}\partial^{d}\partial^{f}B_{\rm ext}^{c} \partial_{a,\boldsymbol{p}}J_{Q,s}^{bdf} -\frac{q_{s}}{m_{s}}\nabla_{\boldsymbol{p}}\boldsymbol{\cdot} \int {\rm d} \boldsymbol{r}' \, {\rm d} \boldsymbol{p}' (\nabla_{\boldsymbol{r}}G(\boldsymbol{r}-\boldsymbol{r}')) f_{2}(\boldsymbol{r},\boldsymbol{r}',\boldsymbol{p},\boldsymbol{p}',t)\nonumber\\ & \quad-\frac{q_{s}}{m_{s}}\nabla_{\boldsymbol{p}}\boldsymbol{\cdot} \int {\rm d} \boldsymbol{r}' \, {\rm d} \boldsymbol{p}' (\nabla_{\boldsymbol{r}}\partial^{b}G(\boldsymbol{r}-\boldsymbol{r}')) [d_{2}^{b}(\boldsymbol{r},\boldsymbol{r}',\boldsymbol{p}, \boldsymbol{p}',t)-d_{2}^{b}(\boldsymbol{r}',\boldsymbol{r}, \boldsymbol{p}',\boldsymbol{p},t)]\nonumber\\ & \quad-\frac{q_{s}}{m_{s}}\nabla_{\boldsymbol{p}}\boldsymbol{\cdot} \int {\rm d} \boldsymbol{r}' \, {\rm d} \boldsymbol{p}' (\nabla_{\boldsymbol{r}}\partial^{b}\partial^{c}G(\boldsymbol{r}-\boldsymbol{r}')) [Q_{2}^{bc}(\boldsymbol{r},\boldsymbol{r}',\boldsymbol{p},\boldsymbol{p}',t)\nonumber\\ & \quad+Q_{2}^{bc}(\boldsymbol{r}',\boldsymbol{r},\boldsymbol{p}',\boldsymbol{p},t) -D_{2}^{bc}(\boldsymbol{r},\boldsymbol{r}',\boldsymbol{p},\boldsymbol{p}',t) -D_{2}^{cb}(\boldsymbol{r},\boldsymbol{r}',\boldsymbol{p},\boldsymbol{p}',t)] =0, \end{align}

where

(5.4)\begin{equation} d^{a}(\boldsymbol{r},\boldsymbol{p},t)=\frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \sum_{i=1}^{N} \xi^{a} \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i} \end{equation}

is the distribution function of dipole moment divided by the charge $q_{s}$,

(5.5)\begin{equation} Q^{ab}(\boldsymbol{r},\boldsymbol{p},t)=\frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \int_{\varDelta,\varDelta_{p}} \xi^{a}\xi^{b} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \sum_{i=1}^{N} \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i} \end{equation}

is the distribution function of quadrupole moment divided by the charge $q_{s}$,

(5.6)\begin{equation} d_{2}^{a}(\boldsymbol{r},\boldsymbol{r}',\boldsymbol{p},\boldsymbol{p}',t) =\frac{1}{\varDelta^{2}}\frac{1}{\varDelta_{p}^{2}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi}' \,{\rm d}\boldsymbol{\eta}' \sum_{i,j=1, j\neq i}^{N} \xi^{a} \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i} \delta_{\boldsymbol{r}'j} \delta_{\boldsymbol{p}'j} \end{equation}

is the two-particle dipole–charge distribution function divided by the charge $q_{s}$,

(5.7)\begin{equation} Q_{2}^{ab}(\boldsymbol{r},\boldsymbol{r}',\boldsymbol{p},\boldsymbol{p}',t) =\frac{1}{\varDelta^{2}}\frac{1}{\varDelta_{p}^{2}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi}' \,{\rm d}\boldsymbol{\eta}' \sum_{i,j=1, j\neq i}^{N} \xi^{a}\xi^{b} \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i} \delta_{\boldsymbol{r}'j} \delta_{\boldsymbol{p}'j} \end{equation}

is the two-particle quadrupole–charge distribution function divided by the charge $q_{s}$,

(5.8)\begin{equation} D_{2}^{ab}(\boldsymbol{r},\boldsymbol{r}',\boldsymbol{p},\boldsymbol{p}',t) =\frac{1}{\varDelta^{2}}\frac{1}{\varDelta_{p}^{2}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi}' \,{\rm d}\boldsymbol{\eta}' \sum_{i,j=1, j\neq i}^{N} \xi^{a} \xi'^{b} \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i} \delta_{\boldsymbol{r}'j} \delta_{\boldsymbol{p}'j}, \end{equation}

is the two-particle dipole–dipole distribution function divided by the charge $q_{s}$,

(5.9)\begin{equation} J_{D}^{ab}(\boldsymbol{r},\boldsymbol{p},t) =\frac{1}{\varDelta^{2}}\frac{1}{\varDelta_{p}^{2}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi}' \,{\rm d}\boldsymbol{\eta}' \sum_{i,j=1, j\neq i}^{N} \eta^{a} \xi^{b} \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i} \delta_{\boldsymbol{r}'j} \delta_{\boldsymbol{p}'j} \end{equation}

is the distribution function of flux of dipole moment divided by the charge $q_{s}$, and

(5.10)\begin{equation} J_{Q}^{abc}(\boldsymbol{r},\boldsymbol{p},t) =\frac{1}{\varDelta^{2}}\frac{1}{\varDelta_{p}^{2}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi}' \,{\rm d}\boldsymbol{\eta}' \sum_{i,j=1, j\neq i}^{N} \eta^{a} \xi^{b} \xi^{c} \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i} \delta_{\boldsymbol{r}'j} \delta_{\boldsymbol{p}'j} \end{equation}

is the distribution function of flux of the electric quadrupole moment divided by the charge $q_{s}$.

5.3. Self-consistent field approximation in the multipole approximation of kinetic equation

Physically, the self-consistent field approximation in the multipole regime is the same as the self-consistent field approximation in the monopole regime described above. Technically, we have splitting of the two-particle distribution functions $d_{2}^{a}$, $Q_{2}^{ab}$ and $D_{2}^{ab}$ into the corresponding one-particle distribution functions. It gives a simplification of kinetic equation (5.3):

(5.11)\begin{align} & \partial_{t}f_{s}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f_{s} +\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{f}_{s} +q_{s}\left(\boldsymbol{E}_{\rm ext}+\boldsymbol{E}+\frac{1}{c}\boldsymbol{v}\times\boldsymbol{B}_{\rm ext}\right) \boldsymbol{\cdot}\frac{\partial f_{s}}{\partial \boldsymbol{p}}\nonumber\\ & \quad+q_{s}\left(\partial^{b}(\boldsymbol{E}_{\rm ext}+\boldsymbol{E})+\frac{1}{c}\boldsymbol{v}\times\partial^{b}\boldsymbol{B}_{\rm ext}\right) \boldsymbol{\cdot}\frac{\partial d_{s}^{b}}{\partial \boldsymbol{p}}\nonumber\\ & \quad+q_{e}\left(\partial^{b}\partial^{c}(\boldsymbol{E}_{\rm ext}+\boldsymbol{E}) +\frac{1}{c}\boldsymbol{v}\times\partial^{b}\partial^{c}\boldsymbol{B}_{\rm ext}\right) \boldsymbol{\cdot}\frac{\partial Q_{s}}{\partial \boldsymbol{p}}\nonumber\\ & \quad+\frac{q_{s}}{c}\varepsilon^{abc}B_{\rm ext}^{c}\partial_{a,\boldsymbol{p}} f_{s}^{b} +\frac{q_{s}}{c}\varepsilon^{abc}\partial^{d}B_{\rm ext}^{c} \partial_{a,\boldsymbol{p}}J_{D,s}^{bd} +\frac{q_{s}}{c}\varepsilon^{abc}\partial^{d}\partial^{f}B_{\rm ext}^{c} \partial_{a,\boldsymbol{p}}J_{Q,s}^{bdf} =0. \end{align}

Kinetic equation (5.11) contains the self-consistent electric field, which satisfies the following equations: $\boldsymbol {\nabla }\times \boldsymbol {E}=0$ and

(5.12)\begin{align} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{E}& =4{\rm \pi} \sum_{s} q_{s}\left(\int f_{s}(\boldsymbol{r},\boldsymbol{p},t)\, {\rm d} \boldsymbol{p} +\partial^{b}\int {\rm d} _{s}^{b}(\boldsymbol{r},\boldsymbol{p},t)\, {\rm d} \boldsymbol{p}\right. \nonumber\\ & \left. \quad +\,\frac{1}{2} \partial^{b}\partial^{c}\int Q_{s}^{bc}(\boldsymbol{r},\boldsymbol{p},t)\, {\rm d} \boldsymbol{p} +\cdots\right), \end{align}

where the Coulomb interaction leads to the contribution of the charge dynamics along with the dynamics of the kinetic multipole distribution functions.

Complete analysis of (5.11) requires the kinetic equations for the additional vector and tensor distribution functions. These equations can be derived by the method demonstrated in this paper. However, we do not present these equations here. This paper is focused on the method of derivation of hydrodynamic and kinetic equations, and on the ways of further generalizations of these models including the account of hydrodynamic or kinetic multipole functions. Examples of the closed set of hydrodynamic or kinetic equations consistently including these effects are the subject of further work.

6. Relativistic hydrodynamic model with the average reverse gamma factor evolution and the self-consistent field approximation in the relativistic hydrodynamics

Non-relativistic hydrodynamics and kinetics are considered above. Basic definitions are given above as well. The self-consistent field approximation is discussed in terms of the suggested model. Our next goal is the generalization of this method for the systems of relativistic particles. First, we consider the relativistic hydrodynamics. We choose the form of relativistic hydrodynamics in the form of the hydrodynamic model with the average reverse gamma factor evolution recently suggested by Andreev (Reference Andreev2022a, Reference Andreev2023b). The hydrodynamic model with the average reverse gamma factor evolution is obtained in the monopole approximation. Moreover, the interaction between particles and the self-consistent field approximation is not considered explicitly in the cited papers.

The equation of motion for each particle appears as the evolution of the momentum under action of the Lorentz force:

(6.1)\begin{equation} \dot{\boldsymbol{p}}_{i}(t)=q_{i} \left(\boldsymbol{E}(\boldsymbol{r}_{i}(t),t) +\frac{1}{c}[\boldsymbol{v}_{i}(t), \boldsymbol{B}(\boldsymbol{r}_{i}(t),t)] \right), \end{equation}

where $\boldsymbol {p}_{i}(t)=m_{i}\boldsymbol {v}_{i}(t)/\sqrt {1- \boldsymbol {v}_{i}^{2}(t)/c^{2}}$ is the relativistic momentum, $\boldsymbol {E}_{i}=\boldsymbol {E}_{i,{\rm ext}}+\boldsymbol {E}_{i,{\rm int}}$, $\boldsymbol {B}_{i}=\boldsymbol {B}_{i,{\rm ext}}+\boldsymbol {B}_{i,{\rm int}}$, $\boldsymbol {E}_{i}=\boldsymbol {E}(\boldsymbol {r}_{i}(t),t)$ and $\boldsymbol {B}_{i}=\boldsymbol {B}(\boldsymbol {r}_{i}(t),t)$. The electric $\boldsymbol {E}_{i,{\rm int}}$ and magnetic $\boldsymbol {B}_{i,{\rm int}}$ fields caused by particles surrounding the $i$th particle are $\boldsymbol {E}_{i,{\rm int}}=-\nabla _{i}\varphi (\boldsymbol {r}_{i}(t),t)- ({1}/{c})\partial _{t}\boldsymbol {A}(\boldsymbol {r}_{i}(t),t)$ and $\boldsymbol {B}_{i,{\rm int}}=\nabla _{i}\times \boldsymbol {A}(\boldsymbol {r}_{i}(t),t)$ with

(6.2)\begin{equation} \varphi(\boldsymbol{r}_{i}(t),t)=\sum_{j\neq i}q_{j}\int \frac{\delta \left(t-t'-\dfrac{1}{c}|\boldsymbol{r}_{i}(t)- \boldsymbol{r}_{j}(t')|\right)}{|\boldsymbol{r}_{i}(t)-\boldsymbol{r}_{j}(t')|}\, {\rm d} t' \end{equation}

and

(6.3)\begin{equation} \boldsymbol{A}(\boldsymbol{r}_{i}(t),t)=\sum_{j\neq i}q_{j}\int \frac{\delta\left(t-t'-\dfrac{1}{c}|\boldsymbol{r}_{i}(t)- \boldsymbol{r}_{j}(t')|\right)}{|\boldsymbol{r}_{i}(t)-\boldsymbol{r}_{j}(t')|} \frac{\boldsymbol{v}_{j}(t')}{c}\, {\rm d} t'. \end{equation}

Here, we use the Green function of the retarding electromagnetic interaction

(6.4)\begin{equation} \tilde{G}_{ij}= \frac{\delta\left(t-t'-\dfrac{1}{c}|\boldsymbol{r}_{i}(t)- \boldsymbol{r}_{j}(t')|\right)}{|\boldsymbol{r}_{i}(t)- \boldsymbol{r}_{j}(t')|}. \end{equation}

Here, we show the scalar and vector potentials of the electromagnetic field acting on $i$th particles and, therefore, created by the surrounding particles. These potentials (6.2) and (6.3) can be represented via corresponding Maxwell equations. The radiation reaction (see Rohrlich Reference Rohrlich1990; Spohn Reference Spohn2004) is not considered in this paper. It is left to be considered in future works. We use the Lorenz gauge $\partial _{t}\varphi /c+\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {A}=0$ in the sections related to the relativistic plasmas.

In relativistic plasmas, the electromagnetic field has dynamical degrees of freedom. While the system of $N$ particles has 3$N$ degrees of freedom, the electromagnetic field has an infinite number of degrees of freedom. However, this continuum of degrees of freedom can be considered via the positions and velocities of particles in the previous moments of time. Hence, the continuum of degrees of freedom of the field is represented by the continuum of time in accordance with (6.2) and (6.3).

Relativistic hydrodynamics is derived by the method described above by Andreev (Reference Andreev2022a, Reference Andreev2023b). However, the details of the self-consistent field approximation is not considered in these papers, while equations are obtained there in this approximation.

Let us consider the definition of number density $n$ in (2.8) for the relativistic regime. It has same form (2.8) which is represented in the shorter form $n=\langle m_{i}\rangle /m$. We contract the $\varDelta$-neighbourhoods in the arbitrary inertial frame. If we consider the transition to another inertial frame which moves relatively the first frame with the constant velocity $\boldsymbol {V}=\{V,0,0\}$, we use the global Lorentz transformation. All $\varDelta$-neighbourhoods change their form due to the contraction of distance in the direction of motion of the second inertial frame: $\Delta x'=\sqrt {1-V^{2}/c^{2}}\Delta x$. So, the $\varDelta$-vicinities in the second frame are not spherical any more. Their volume also changes. However, the number of particles in each neighbourhood does not change. So, we have a formal change of the number density at the change of the frame. Anyway, direct transition of the $\varDelta$-neighbourhoods to another frame change their essential property. They are constructed in the first frame as the motionless $\varDelta$-neighbourhoods of the points of space. While they move in the second frame, they are vicinities of moving points. Hence, the proper formulation of hydrodynamics in the second frame requires the reconstruction of the $\varDelta$-neighbourhoods around a motionless point of space. Therefore, we are focused on the study of the relativistic effects in the fixed frame. We do not make reference to the ‘rest frame’ since it exists for the relatively simple collective motion of plasmas, but it does not exist in the general case.

Let us consider the evolution of number density (2.8) in the inertial frame. We find that it expresses itself via the average velocity of particles $\boldsymbol {j}=\langle \dot {\boldsymbol {r}}_{i}\rangle =\langle \boldsymbol {v}_{i}\rangle \equiv n\boldsymbol {v}$. It manifests itself via the continuity equation:

(6.5)\begin{equation} \partial_{t}n+\boldsymbol{\nabla}\boldsymbol{\cdot}(n\boldsymbol{v})=0. \end{equation}

If we want to continue the set of hydrodynamic equations, we need to consider the evolution of current $\boldsymbol {j}$. Since the current is the average velocity $\langle \boldsymbol {v}_{i}\rangle$, we need to consider the accelerations of all particles. Therefore, we need to rewrite the equations of motion of each particle (6.1) for the velocity change instead of the momentum change:

(6.6)\begin{equation} \dot{\boldsymbol{v}}_{i} =\frac{e_{i}}{m_{i}}\sqrt{1-\frac{\boldsymbol{v}_{i}^{2}}{c^{2}}} \left[\boldsymbol{E}_{i}+\frac{1}{c}[\boldsymbol{v}_{i}\times\boldsymbol{B}_{i}] -\frac{1}{c^{2}}\boldsymbol{v}_{i}(\boldsymbol{v}_{i}\boldsymbol{\cdot}\boldsymbol{E}_{i})\right]. \end{equation}

Some further details for the derivation of the general structure of the hydrodynamic model with the average reverse gamma factor evolution presented below can be found from Andreev (Reference Andreev2022a, Reference Andreev2023b). However, we present the part related to the interaction.

Here, we consider the evolution of the average velocity $\boldsymbol {j}=\langle \boldsymbol {v}_{i}\rangle$:

(6.7)\begin{equation} \partial_{t}j^{a} =\partial_{t}\left[\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi} \sum_{i=1}^{N} v_{i}^{a}(t) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t))\right] ={-}\partial_{b}\varPi^{ab} +\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi} \sum_{i=1}^{N} \dot{v}_{i}^{a}(t) \delta_{i}, \end{equation}

where

(6.8)\begin{equation} \varPi^{ab}= \frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi} \sum_{i=1}^{N} v_{i}^{a}(t)v_{i}^{b}(t) \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \end{equation}

and $\delta _{i}=\delta (\boldsymbol {r}+ \boldsymbol {\xi }-\boldsymbol {r}_{i}(t))$.

Use of (6.6) in (6.7) gives the following representation of $\partial _{t}j^{a}$:

(6.9)\begin{align} \partial_{t}j^{a} & ={-}\partial_{b}\varPi^{ab} +\frac{q_{s}}{m_{s}}\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi} \sum_{i=1}^{N}\frac{\delta_{i}}{\gamma_{i}(t)} \left(E^{a}_{i,{\rm ext}} +\frac{1}{c}\varepsilon^{abc}v_{i}^{b}B_{i,{\rm ext}}^{c} -\frac{1}{c^{2}}v_{i}^{a} v_{i}^{b} E^{b}_{i,{\rm ext}}\right)\nonumber\\ & \quad +\frac{q_{s}}{m_{s}}\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi} \sum_{i=1}^{N}\frac{\delta_{i}}{\gamma_{i}(t)} \left(-\partial^{a}_{i}\varphi_{i,{\rm int}}-\frac{1}{c}\partial_{t}A_{i,{\rm int}}^{a}\right.\nonumber\\ & \quad \left.+\,\frac{1}{c}\varepsilon_{abc}v_{i}^{b}\varepsilon^{cdf}\partial_{i}^{d}A_{i,{\rm int}}^{f} -\frac{1}{c^{2}}v_{i}^{a} v_{i}^{b}\left(-\partial^{b}_{i}\varphi_{i,{\rm int}}-\frac{1}{c} \partial_{t}A_{i,{\rm int}}^{b}\right)\right). \end{align}

For the external field, we can make the following transformation using the delta function $\delta _{i}$: $E^{a}_{i,{\rm ext}}=E^{a}_{\rm ext}(\boldsymbol {r}_{i}(t),t)=E^{a}_{\rm ext}(\boldsymbol {r}+ \boldsymbol {\xi },t)$ and $B^{a}_{i,{\rm ext}}=B^{a}_{\rm ext}(\boldsymbol {r}_{i}(t),t)=B^{a}_{\rm ext}(\boldsymbol {r}+ \boldsymbol {\xi },t)$. So, the external electromagnetic field can be expanded on $\boldsymbol {\xi }$ if the external electromagnetic field has small changes over the $\varDelta$ volume. Hence, we find

(6.10)\begin{align} F^{a}_{\rm ext}& =\frac{q_{s}}{m_{s}}\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi} \sum_{i=1}^{N}\frac{\delta_{i}}{\gamma_{i}(t)}\left(E^{a}_{i,{\rm ext}} +\frac{1}{c}\varepsilon^{abc}v_{i}^{b}B_{i,{\rm ext}}^{c} -\frac{1}{c^{2}}v_{i}^{a} v_{i}^{b} E^{b}_{i,{\rm ext}}\right)\nonumber\\ & =\frac{q_{s}}{m_{s}}\left[ \varGamma E^{a} +\frac{1}{c}\varepsilon^{abc}\varTheta^{b}B^{c}-\frac{1}{c^{2}}\varXi^{ab}E^{b} +\varGamma_{D}^{b}\partial_{b}E^{a} +\frac{1}{c}\varepsilon^{abc}\varTheta_{D}^{bd}\partial_{d}B^{c}\right.\nonumber\\ & \quad -\left.\frac{1}{c^{2}}\varXi_{D}^{abc}\partial_{c}E^{b} +\varGamma_{Q}^{bc}\partial_{b}\partial_{c}E^{a}+\frac{1}{c}\varepsilon^{abc}\varTheta_{Q}^{bdf}\partial_{d}\partial_{f}B^{c} -\frac{1}{c^{2}}\varXi_{Q}^{abcd}\partial_{c}\partial_{d}E^{b}+\cdots\right], \end{align}

where the subindex $D$ refers to dipolar and the subindex $Q$ refers to quadrupolar. The monopolar terms found in Andreev (Reference Andreev2022a) contain the following functions: $\varGamma =\langle \gamma _{i}^{-1}(t)\rangle$, $\varTheta ^{b}=\langle \gamma _{i}^{-1}(t)v_{i}^{b}\rangle$ and $\varXi ^{ab}=\langle \gamma _{i}^{-1}(t)v_{i}^{a}v_{i}^{b}\rangle$. The dipolar terms contain the following functions: $\varGamma _{D}^{b}=\langle \gamma _{i}^{-1}(t)\xi ^{b}\rangle$, $\varTheta _{D}^{bd}=\langle \gamma _{i}^{-1}(t)v_{i}^{b}\xi ^{d}\rangle$ and $\varXi _{D}^{abc}=\langle \gamma _{i}^{-1}(t)v_{i}^{a}v_{i}^{b}\xi ^{c}\rangle$. The quadrupolar terms contain the following functions $\varGamma _{Q}^{bc}=\langle \gamma _{i}^{-1}(t)\xi ^{b}\xi ^{b}\rangle$, $\varTheta _{Q}^{bdf}=\langle \gamma _{i}^{-1}(t)v_{i}^{b}\xi ^{d}\xi ^{f}\rangle$ and $\varXi _{Q}^{abcd}=\langle \gamma _{i}^{-1}(t)v_{i}^{a}v_{i}^{b}\xi ^{c}\xi ^{d}\rangle$.

To consider the force field of interaction, we start with a single term (it is constructed of the two first terms of (6.9))

(6.11)\begin{equation} F^{a}_{{\rm int},1}=\frac{q_{s}}{m_{s}}\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi} \sum_{i=1}^{N}\frac{\delta_{i}}{\gamma_{i}(t)}E_{i,{\rm int}}^{a}, \end{equation}

where $E_{i,{\rm int}}^{a}=-\partial ^{a}_{i}\varphi _{i,{\rm int}}-({1}/{c})\partial _{t}A_{i,{\rm int}}^{a}$. Using potentials (6.2) and (6.3), we represent the force via the Green function $G_{ij}$ in (6.4):

(6.12)\begin{equation} F^{a}_{{\rm int},1}={-}\frac{q_{s}q_{j}}{m_{s}}\frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\int {\rm d} t' \sum_{i,j=1, i\neq j}^{N}\frac{\delta_{i}}{\gamma_{i}(t)} \left(\partial_{i}^{a}G_{ij}+\frac{v_{j}^{a}(t')}{c^{2}}\partial_{t}G_{ij}\right). \end{equation}

Next, we include the additional delta function

(6.13)\begin{align} & F^{a}_{{\rm int},1}={-}\frac{q_{s}q_{j}}{m_{s}} \frac{1}{\varDelta}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi} \int {\rm d} t' \int {\rm d} \boldsymbol{r}' \sum_{i,j=1, i\neq j}^{N}\frac{1}{\gamma_{i}(t)}\delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t))\nonumber\\ & \quad \times\delta(\boldsymbol{r}'-\boldsymbol{r}_{j}(t')) \left(\partial_{i}^{a}+\frac{v_{j}^{a}(t')}{c^{2}}\partial_{t}\right) G_{ij}(|\boldsymbol{r}_{i}(t)-\boldsymbol{r}_{j}(t')|,t-t'). \end{align}

We modify the argument of the delta function containing the coordinate of the $j$th particle

(6.14)\begin{align} & F^{a}_{{\rm int},1}={-}\frac{q_{s}q_{j}}{m_{s}} \frac{1}{\varDelta^{2}}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}' \int {\rm d} t' \int {\rm d} \boldsymbol{r}' \sum_{i,j=1, i\neq j}^{N} \frac{1}{\gamma_{i}(t)} \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t))\nonumber\\ & \quad\delta(\boldsymbol{r}'+ \boldsymbol{\xi}'-\boldsymbol{r}_{j}(t')) \left(\partial_{i}^{a}+\frac{v_{j}^{a}(t')}{c^{2}}\partial_{t}\right) G_{ij}(|\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}- \boldsymbol{\xi}'|,t-t'), \end{align}

where we also used the modified argument of $\delta _{j}$ to represent the argument of the Green function.

In the following step, we are ready to expand the Green function on $\boldsymbol {\xi }- \boldsymbol {\xi }'$ assuming that the Green function has small change over the $\varDelta$-neighbourhood scale:

(6.15)\begin{align} F^{a}_{{\rm int},1}& ={-}\frac{q_{s}q_{j}}{m_{s}} \frac{1}{\varDelta^{2}}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}' \int {\rm d} t' \int {\rm d} \boldsymbol{r}' \nonumber\\ & \quad\sum_{i,j=1, i\neq j}^{N}\frac{1}{\gamma_{i}(t)} \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{r}'+ \boldsymbol{\xi}'-\boldsymbol{r}_{j}(t')) \left(\partial_{i}^{a}+\frac{v_{j}^{a}(t')}{c^{2}}\partial_{t}\right)\nonumber\\ & \quad\left(G_{ij}(|\boldsymbol{r}-\boldsymbol{r}|,t-t') +(\xi^{b}-\xi'^{b})\partial_{r}^{b}G_{ij}(|\boldsymbol{r}-\boldsymbol{r}|,t-t')\vphantom{\frac{1}{2}} \right.\nonumber\\ & \quad+\left.\frac{1}{2}(\xi^{b}-\xi'^{b})(\xi^{c}-\xi'^{c})\partial_{r}^{b}\partial_{r}^{c} G_{ij}(|\boldsymbol{r}-\boldsymbol{r}|,t-t')+\cdots\right). \end{align}

Let us to interpret the expression presented for the force field via the two-particle macroscopic functions

(6.16)\begin{align} F^{a}_{{\rm int},1}& ={-}\frac{q_{s}q_{j}}{m_{s}} \int {\rm d} t' \int {\rm d} \boldsymbol{r}' \left[ \partial_{r}^{a}G \cdot \varGamma_{2}+\frac{1}{c^{2}}\partial_{t}G \cdot X_{1}^{a}+\partial_{r}^{a}\partial_{r}^{b}G \cdot X_{2}^{b}+\frac{1}{c^{2}}\partial_{r}^{b}\partial_{t}G \cdot X_{3}^{ab}\right.\nonumber\\ & \quad -\partial_{r}^{a}\partial_{r}^{b}G \cdot X_{4}^{b}-\frac{1}{c^{2}}\partial_{r}^{b}\partial_{t}G \cdot X_{5}^{ab} +\partial_{r}^{a}\partial_{r}^{b}\partial_{r}^{c}G \cdot X_{6}^{bc}+\frac{1}{c^{2}}\partial_{r}^{b}\partial_{r}^{c}\partial_{t}G \cdot X_{7}^{abc}\nonumber\\ & \quad+\left.\partial_{r}^{a}\partial_{r}^{b}\partial_{r}^{c}G \cdot X_{8}^{bc}+\frac{1}{c^{2}}\partial_{r}^{b}\partial_{r}^{c}\partial_{t}G \cdot X_{9}^{abc} -\partial_{r}^{a}\partial_{r}^{b}\partial_{r}^{c}G \cdot X_{10}^{bc}-\frac{1}{c^{2}}\partial_{r}^{b}\partial_{r}^{c}\partial_{t}G \cdot X_{11}^{abc} \right], \end{align}

where $G=G(\boldsymbol {r}-\boldsymbol {r}',t-t')$. We use the following two-particle macroscopic functions. Here we present them together with their limits for the mean-field approximation. Limits required to get the mean-field approximation are discussed above. It is placed in § 3.1. The first part is placed in the paragraph before (3.6), where the relative position of two $\varDelta$-neighbourhoods in the coordinate space is discussed. We also need to refer to the discussion below (3.6) related to figures 3–6. However, in the relativistic regime, we also should keep in mind that we consider the position of the $j$th particle in different moments of time $t'< t$. Equation (6.16) contains the following functions:

(6.17)\begin{align} & \varGamma_{2}(\boldsymbol{r},\boldsymbol{r}',t,t')= \frac{1}{\varDelta^{2}}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}\int_{\varDelta}\, {\rm d} \boldsymbol{\xi}' \sum_{i,j=1, i\neq j}^{N}\frac{1}{\gamma_{i}(t)} \delta(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}_{i}(t)) \delta(\boldsymbol{r}'+ \boldsymbol{\xi}'-\boldsymbol{r}_{j}(t'))\nonumber\\ & \quad \equiv \langle\langle\gamma_{i}^{{-}1}(t)\rangle\rangle \rightarrow \varGamma(\boldsymbol{r},t)\boldsymbol{\cdot} n(\boldsymbol{r}',t'), \end{align}

(6.18)\begin{gather} X_{1}^{a}(\boldsymbol{r},\boldsymbol{r}',t,t') =\langle\langle\gamma_{i}^{{-}1}(t)\boldsymbol{\cdot} v_{j}^{a}(t')\rangle\rangle \rightarrow \varGamma(\boldsymbol{r},t)\boldsymbol{\cdot} j^{a}(\boldsymbol{r}',t'), \end{gather}

(6.19)\begin{gather}X_{2}^{b}(\boldsymbol{r},\boldsymbol{r}',t,t') =\langle\langle\gamma_{i}^{{-}1}(t)\boldsymbol{\cdot} \rangle\rangle \rightarrow \varGamma_{D}^{b}(\boldsymbol{r},t)\boldsymbol{\cdot} n(\boldsymbol{r}',t'), \end{gather}

(6.20)\begin{gather}X_{3}^{ab}(\boldsymbol{r},\boldsymbol{r}',t,t') =\langle\langle\gamma_{i}^{{-}1}(t)\boldsymbol{\cdot} v_{j}^{a}(t')\rangle\rangle \rightarrow \varGamma_{D}^{b}(\boldsymbol{r},t)\boldsymbol{\cdot} j^{a}(\boldsymbol{r}',t'), \end{gather}

(6.21)\begin{gather}X_{4}^{b}(\boldsymbol{r},\boldsymbol{r}',t,t') =\langle\langle\gamma_{i}^{{-}1}(t)\boldsymbol{\cdot} \rangle\rangle \rightarrow \varGamma(\boldsymbol{r},t)\boldsymbol{\cdot} d^{b}(\boldsymbol{r}',t'), \end{gather}

(6.22)\begin{gather}X_{5}^{ab}(\boldsymbol{r},\boldsymbol{r}',t,t') =\langle\langle\gamma_{i}^{{-}1}(t)\boldsymbol{\cdot} v_{j}^{a}(t')\rangle\rangle \rightarrow \varGamma(\boldsymbol{r},t)\boldsymbol{\cdot} J_{D}^{ab}(\boldsymbol{r}',t'), \end{gather}

(6.23)\begin{gather}X_{6}^{bc}(\boldsymbol{r},\boldsymbol{r}',t,t') =\langle\langle\gamma_{i}^{{-}1}(t)\boldsymbol{\cdot} \rangle\rangle \rightarrow \varGamma_{Q}^{bc}(\boldsymbol{r},t)\boldsymbol{\cdot} n(\boldsymbol{r}',t'), \end{gather}

(6.24)\begin{gather}X_{7}^{abc}(\boldsymbol{r},\boldsymbol{r}',t,t') =\langle\langle\gamma_{i}^{{-}1}(t)\boldsymbol{\cdot} v_{j}^{a}(t')\rangle\rangle \rightarrow \varGamma_{Q}^{bc}(\boldsymbol{r},t)\boldsymbol{\cdot} j^{a}(\boldsymbol{r}',t'), \end{gather}

(6.25)\begin{gather}X_{8}^{bc}(\boldsymbol{r},\boldsymbol{r}',t,t') =\langle\langle\gamma_{i}^{{-}1}(t)\boldsymbol{\cdot} \rangle\rangle \rightarrow \varGamma(\boldsymbol{r},t)\boldsymbol{\cdot} Q^{bc}(\boldsymbol{r}',t'), \end{gather}

(6.26)\begin{gather}X_{9}^{abc}(\boldsymbol{r},\boldsymbol{r}',t,t') =\langle\langle\gamma_{i}^{{-}1}(t)\boldsymbol{\cdot} v_{j}^{a}(t')\rangle\rangle \rightarrow \varGamma(\boldsymbol{r},t)\boldsymbol{\cdot} J_{Q}^{abc}(\boldsymbol{r}',t'), \end{gather}

(6.27)\begin{gather}X_{10}^{bc}(\boldsymbol{r},\boldsymbol{r}',t,t') =\langle\langle\gamma_{i}^{{-}1}(t)\boldsymbol{\cdot} \rangle\rangle \rightarrow \varGamma_{D}^{b}(\boldsymbol{r},t)\boldsymbol{\cdot} d^{c}(\boldsymbol{r}',t') \end{gather}

and

(6.28)\begin{equation} X_{11}^{abc}(\boldsymbol{r},\boldsymbol{r}',t,t') =\langle\langle\gamma_{i}^{{-}1}(t)\boldsymbol{\cdot} v_{j}^{a}(t')\rangle\rangle \rightarrow \varGamma_{D}^{b}(\boldsymbol{r},t)\boldsymbol{\cdot} J_{D}^{ac}(\boldsymbol{r}',t'). \end{equation}

The force field (6.16) being represented in the self-consistent field (mean-field) approximation can be reconstructed as the six groups of terms

(6.29)\begin{align} & F^{a}_{{\rm int},1}={-}\frac{q_{s}q_{j}}{m_{s}}\left[ -\varGamma\partial^{a}\left(\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') n(\boldsymbol{r}',t') -\partial^{b}\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') d^{b}(\boldsymbol{r}',t')\right.\right.\nonumber\\ & \quad \left.+\,\frac{1}{2}\partial^{b}\partial^{c}\!\int {\rm d} t' \!\int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') Q^{bc}(\boldsymbol{r}',t') +\cdots\right) -\frac{1}{c}\varGamma\boldsymbol{\cdot}\frac{1}{c}\partial_{t}\left(\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') j^{a}(\boldsymbol{r}',t')\right.\nonumber\\ & \quad-\left.\partial^{b}\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') J_{D}^{ab}(\boldsymbol{r}',t') +\frac{1}{2}\partial^{b}\partial^{c}\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') J_{Q}^{abc}(\boldsymbol{r}',t') +\cdots\right)\nonumber\\ & \quad-\varGamma_{D}^{b}\partial^{a}\partial^{b}\left(\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') n(\boldsymbol{r}',t') -\partial^{b}\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') d^{c}(\boldsymbol{r}',t') +\cdots\right)\nonumber\\ & \quad-\frac{1}{c}\varGamma_{D}^{b}\boldsymbol{\cdot}\frac{1}{c}\partial^{b}\partial_{t}\left(\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') j^{a}(\boldsymbol{r}',t') -\partial^{c}\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') J_{D}^{ac}(\boldsymbol{r}',t') +\cdots\right)\nonumber\\ & \quad-\frac{1}{2}\varGamma_{Q}^{bc}\partial^{a}\partial^{b}\partial^{c}\left(\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') n(\boldsymbol{r}',t') +\cdots\right)\nonumber\\ & \quad-\left.\frac{1}{2}\frac{1}{c}\varGamma_{Q}^{bc}\boldsymbol{\cdot} \frac{1}{c}\partial_{t}\partial^{b}\partial^{c}\left(\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') j^{a}(\boldsymbol{r}',t') +\cdots\right) \right]. \end{align}

These structures allows to introduce the macroscopic scalar and vector potentials

(6.30)\begin{align} \varphi(\boldsymbol{r},t)& =\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') n(\boldsymbol{r}',t')-\partial^{b}\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') d^{b}(\boldsymbol{r}',t')\nonumber\\ & \quad +\frac{1}{2}\partial^{b}\partial^{c}\int {\rm d} t' \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') Q^{bc}(\boldsymbol{r}',t') +\cdots, \end{align}

and

(6.31)\begin{align} A^{a}(\boldsymbol{r},t)& =\int {\rm d} t' \left( \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') j^{a}(\boldsymbol{r}',t')-\partial^{b} \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') J_{D}^{ab}(\boldsymbol{r}',t')\right.\nonumber\\ & \quad \left.+\,\frac{1}{2}\partial^{b}\partial^{c} \int {\rm d} \boldsymbol{r}' G(\boldsymbol{r},\boldsymbol{r}',t,t') J_{Q}^{abc}(\boldsymbol{r}',t') +\cdots\right). \end{align}

These potentials allow to introduce the macroscopic electric and magnetic fields $\boldsymbol {E}=-\boldsymbol {\nabla }\varphi (\boldsymbol {r},t)-(1/c)\partial _{t}\boldsymbol {A}$ and $\boldsymbol {B}=curl \boldsymbol {A}$. Considering part of the force field $F^{a}_{{\rm int},1}$, (6.29) is expressed via the electric field only

(6.32)\begin{equation} F^{a}_{{\rm int},1}={-}\frac{q_{s}q_{j}}{m_{s}} \left(\varGamma E^{a} +\varGamma_{D}^{b}\partial^{b}E^{a}+\frac{1}{2} \varGamma_{Q}^{bc}\partial^{b}\partial^{c}E^{a}\right). \end{equation}

Considering other parts of the force field in (6.9), we can obtain the contribution of the dipolar and quadrupolar effects there. These expressions are also presented by the scalar and vector potentials (6.30) and (6.31). The electromagnetic field is presented with the potentials (6.30) and (6.31), in accordance with the explicit form of the Green function (6.4), which satisfies the Maxwell equations

(6.33)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{B}=0, \end{gather}

(6.34)\begin{gather}\boldsymbol{\nabla}\times \boldsymbol{E}={-}\frac{1}{c}\partial_{t}\boldsymbol{B}, \end{gather}

(6.35)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{E}=\sum_{s}4{\rm \pi} q_{s}\left(n_{s}-\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{d}_{s}+\frac{1}{2}\partial_{a}\partial_{b}Q_{s}^{ab}+\cdots\right) \end{gather}

and

(6.36)\begin{equation} (\boldsymbol{\nabla}\times \boldsymbol{B})^{a}=\frac{1}{c}\partial_{t}E^{a} +\sum_{s}\frac{4{\rm \pi} q_{s}}{c}\left(n_{s}v_{s}^{a}-\partial_{b}J^{ab}_{D,s} +\frac{1}{2}\partial_{b}\partial_{c}J_{Q,s}^{abc} +\cdots\right). \end{equation}

Other group of terms in (6.9) can be presented by the macroscopic scalar and vector potentials in a similar way.

In the monopole approximation of the hydrodynamic model with the average reverse gamma factor evolution, we have the following equation of the velocity field evolution:

(6.37)\begin{align} & n\partial_{t}v^{a}+n(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla})v^{a}+\partial^{a}\tilde{p} =\frac{q}{m}\varGamma E^{a}+\frac{q}{mc}\varepsilon^{abc}(\varGamma v^{b}+t^{b})B^{c}\nonumber\\ & \quad -\frac{q}{mc^{2}}(\varGamma v^{a} v^{b}+v^{a}t^{b}+v^{b}t^{a})E^{b} -\frac{e}{mc^{2}}\tilde{t}E^{a}, \end{align}

where $\varGamma =\langle {1}/{\gamma _{i}}\rangle$, $t^{a}=\langle ({1}/{\gamma _{i}})v_{i}^{a} \rangle -\varGamma v^{a}$, $p^{ab}=\langle v_{i}^{a}v_{i}^{b} \rangle -n v^{a}v^{b}$, $t^{ab}=\langle ({1}/{\gamma _{i}})v_{i}^{a}v_{i}^{b} \rangle -\varGamma v^{a}v^{b}-t^{a}v^{b}- v^{a}t^{b}$, $\gamma _{i}=1/\sqrt {1-\boldsymbol {v}_{i}(t)^{2}/c^{2}}$. Two equations of state should be applied for functions $\tilde {p}$ ($p^{ab}=\tilde {p}\delta ^{ab}$) and $\tilde {t}$ ($t^{ab}=\tilde {t}\delta ^{ab}$) (see Andreev Reference Andreev2022a). We also have the corresponding simplification of the Maxwell equations (6.33)–(6.36): $\boldsymbol {\nabla } \boldsymbol {\cdot }\boldsymbol {B}=0$,

(6.38a,b)\begin{gather} \boldsymbol{\nabla}\times \boldsymbol{E}={-}\frac{1}{c}\partial_{t}\boldsymbol{B}, \quad \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{E}=4{\rm \pi} (en_{i}-en_{e}), \end{gather}

(6.39)\begin{gather} \boldsymbol{\nabla}\times \boldsymbol{B}=\frac{1}{c}\partial_{t}\boldsymbol{E}+\sum_{s}\frac{4{\rm \pi} q_{s}}{c}n_{s}\boldsymbol{v}_{s}. \end{gather}

Next, the hydrodynamic model with the average reverse gamma factor evolution suggests the derivation of the equation for evolution of the average reverse gamma factor $\varGamma$. This derivation is similar to the derivation presented above for the particle current evolution. It also includes the contribution of the multipole moments. We do not show this derivation assuming that the illustration made for the particle current evolution is enough for the purpose of this paper. The equation for the average reverse gamma factor evolution $\varGamma$ has the following form in agreement with Andreev (Reference Andreev2022a, Reference Andreev2023b):

(6.40)\begin{equation} \partial_{t}\varGamma+\partial_{b}(\varGamma v^{b}+t^{b}) ={-}\frac{q}{mc^{2}}n\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{E} \left(1-\frac{1}{c^{2}}\left(\boldsymbol{v}^{2}+\frac{5p}{n}\right)\right). \end{equation}

The fourth and final equation for the evolution of the material field in the hydrodynamic model with the average reverse gamma factor evolution is the evolution of the flux of the average reverse gamma factor (Andreev Reference Andreev2022a, Reference Andreev2023b):

(6.41)\begin{align} & (\partial_{t}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla})t^{a}+\partial_{a}\tilde{t} +(\boldsymbol{t}\boldsymbol{\cdot}\boldsymbol{\nabla}) v^{a}+t^{a} (\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{v}) +\varGamma(\partial_{t}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla})v^{a}\nonumber\\ & \quad=\frac{q}{m}nE^{a}\left[1-\frac{\boldsymbol{v}^{2}}{c^{2}}-\frac{3p}{nc^{2}}\right]+\frac{q}{mc}\varepsilon^{abc}nv^{b}B^{c}\left[1-\frac{\boldsymbol{v}^{2}}{c^{2}}-\frac{5p}{nc^{2}}\right]\nonumber\\ & \qquad-\frac{2q}{mc^{2}}E^{a}p\left[1-\frac{\boldsymbol{v}^{2}}{c^{2}}\right]-\frac{q}{mc^{2}}nv^{a}v^{b}E^{b}\left[1-\frac{\boldsymbol{v}^{2}}{c^{2}}-\frac{9p}{nc^{2}}\right] +\frac{10q}{3mc^{4}}\tilde{M} E^{a}. \end{align}

Function $\tilde {M}$ appears as the simplification for the fourth rank tensor $M^{abcd}=({\tilde {M}}/{3})(\delta ^{ab}\delta ^{cd}+ \delta ^{ac}\delta ^{bd}+\delta ^{ad}\delta ^{bc})$ (see (17) of Andreev Reference Andreev2022a). Equation (6.41) shows that we need to include the third equation of state for the function $\tilde {M}$. It is found in Andreev (Reference Andreev2022a). Functions $\langle \gamma _i \rangle$ and $\langle v_i / \gamma _i \rangle$ are independent functions along with the concentration and velocity field, but we need to calculate their value in the equilibrium state (Andreev Reference Andreev2022a).

7. Self-consistent field approximation in the relativistic Vlasov equation

Above we present the derivation of the kinetic equation in the non-relativistic regime, where the interaction is restricted by the Coulomb interaction. However, the contribution of the multipole moments is considered as well. Here, we consider the kinetic model for the relativistic motion of particles (6.1). Moreover, we assume that the interaction between particles is the full electromagnetic interaction, so the field acting on the $i$th particle is created by surrounding particles in accordance with the microscopic Maxwell equations (see (6.1)–(6.3)). We use the definition (4.2) for the distribution function. Hence, its time derivative has the form of (4.5). We use (6.1) for the time derivative of the momentum of the $i$th particle, which gives the following equation:

(7.1)\begin{align} & \partial_{t}f +\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{F}(\boldsymbol{r},\boldsymbol{p},t) +\frac{q_{s}}{m_{s}}\nabla_{\boldsymbol{p}}\boldsymbol{\cdot} \frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta}\nonumber\\ & \quad\sum_{i=1}^{N/2} \left(\boldsymbol{E}_{\rm ext}(\boldsymbol{r}+ \boldsymbol{\xi},t) +\frac{1}{c}[\boldsymbol{v}_{i}(t)\times \boldsymbol{B}_{\rm ext}(\boldsymbol{r}+ \boldsymbol{\xi},t)]\right) \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i}\nonumber\\ & \quad-\frac{q_{s}}{m_{s}}q_{s'}\nabla_{\boldsymbol{p}}\boldsymbol{\cdot} \frac{1}{\varDelta^{2}}\frac{1}{\varDelta_{p}^{2}} \int {\rm d} t' \int {\rm d} \boldsymbol{r}' \, {\rm d} \boldsymbol{p}' \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \, {\rm d} \boldsymbol{\xi}' \,{\rm d}\boldsymbol{\eta}' \sum_{i=1}^{N/2}\sum_{j=1, j\neq i}^{N} \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i} \delta_{\boldsymbol{r}'j} \delta_{\boldsymbol{p}'j}\nonumber\\ & \quad \times\left( \left(1-\frac{\boldsymbol{v}_{i}(t)\boldsymbol{\cdot}\boldsymbol{v}_{j}(t')}{c^{2}}\right) \nabla_{\boldsymbol{r}}G(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}'- \boldsymbol{\xi}')\right.\nonumber\\ & \quad+\left.\frac{\boldsymbol{v}_{j}(t')}{c^{2}}(\partial_{t}+ (\boldsymbol{v}_{i}(t)\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{r}})G(\boldsymbol{r}+ \boldsymbol{\xi}-\boldsymbol{r}'- \boldsymbol{\xi}')\right) =0, \end{align}

where $\delta _{\boldsymbol {r}'j}\equiv \delta (\boldsymbol {r}'+ \boldsymbol {\xi }'-\boldsymbol {r}_{j}(t))$ and $\delta _{\boldsymbol {p}'j}\equiv \delta (\boldsymbol {p}'+ \boldsymbol {\eta }'-\boldsymbol {p}_{j}(t))$. Equation (7.1) contains function

(7.2)\begin{equation} \boldsymbol{F}(\boldsymbol{r},\boldsymbol{p},t)= \frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \sum_{i=1}^{N/2} \boldsymbol{v}_{i}(t) \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i}. \end{equation}

We can use the relativistic expression for the velocity via the momentum $\boldsymbol {v}_{i}(t)=\boldsymbol {p}_{i}(t)c/\sqrt {\boldsymbol {p}_{i}^{2}(t) +m_{i}^{2}c^{2}}$ with further replacement of the momentum $\boldsymbol {p}_{i}(t)$ on $\boldsymbol {p}+\boldsymbol {\eta }$. Hence, function $\boldsymbol {F}(\boldsymbol {r},\boldsymbol {p},t)$ contains non-polinomic dependence on $\boldsymbol {\eta }$:

(7.3)\begin{equation} \boldsymbol{F}(\boldsymbol{r},\boldsymbol{p},t)=\frac{1}{\varDelta} \frac{1}{\varDelta_{p}}\int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \sum_{i=1}^{N/2}\frac{(\boldsymbol{p}+\boldsymbol{\eta})c}{\sqrt{(\boldsymbol{p}+ \boldsymbol{\eta})^{2}+m_{i}^{2}c^{2}}} \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i}. \end{equation}

The same replacement should be made for the velocities in other terms in (7.1).

Let us consider the expansion on $\boldsymbol {\eta }$ as the small value in comparison with $\boldsymbol {p}$ and $m_{i}c$:

(7.4)\begin{equation} \boldsymbol{F}(\boldsymbol{r},\boldsymbol{p},t)\!=\!\boldsymbol{v}\boldsymbol{\cdot} f(\boldsymbol{r},\boldsymbol{p},t)\!+\!\frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \sum_{i=1}^{N/2} \left(\frac{\boldsymbol{\eta}c}{\sqrt{\boldsymbol{p}^{2}+m_{i}^{2}c^{2}}} -\frac{\boldsymbol{p}(\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\eta})c}{(\sqrt{\boldsymbol{p}^{2}+m_{i}^{2}c^{2}})^{3}}\right) \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i}. \end{equation}

If we consider the monopole regime in the momentum space, we neglect $\boldsymbol {\eta }$ and find $\boldsymbol {F}(\boldsymbol {r},\boldsymbol {p},t)=\boldsymbol {v}\boldsymbol {\cdot } f(\boldsymbol {r},\boldsymbol {p},t)$. Here, we also use $\boldsymbol {p}=m_{s}\boldsymbol {v}/\sqrt {1-\boldsymbol {v}^{2}/c^{2}}$ and $\boldsymbol {v}=\boldsymbol {p}c/\sqrt {\boldsymbol {p}^{2}+m_{s}^{2}c^{2}}$.

We consider this equation in the monopole approximation. Hence, (7.1) simplifies to

(7.5)\begin{align} & \partial_{t}f(\boldsymbol{r},\boldsymbol{p},t)\!+\!\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f \!+\!\frac{q_{s}}{m_{s}}\frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta}\sum_{i=1}^{N/2} \left(\boldsymbol{E}_{\rm ext}(\boldsymbol{r},t) +\frac{1}{c}[\boldsymbol{v}\times \boldsymbol{B}_{\rm ext}(\boldsymbol{r},t)]\right) \delta_{\boldsymbol{r}i} \boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}}\delta_{\boldsymbol{p}i}\nonumber\\ & \quad -\frac{q_{s}}{m_{s}}q_{s'}\boldsymbol{\cdot}\frac{1}{\varDelta^{2}}\frac{1}{\varDelta_{p}^{2}} \int {\rm d} t' \int {\rm d} \boldsymbol{r}' \, {\rm d} \boldsymbol{p}'\times \left( \left(1-\frac{\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{v}'}{c^{2}}\right) \nabla_{\boldsymbol{r}}G(\boldsymbol{r}-\boldsymbol{r}', t-t')\right.\nonumber\\ & \quad+\left.\frac{\boldsymbol{v}'}{c^{2}}(\partial_{t}+(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{r}}))G(\boldsymbol{r}-\boldsymbol{r}', t-t')\right) \nabla_{\boldsymbol{p}}f_{2}(\boldsymbol{r},\boldsymbol{p},t, \boldsymbol{r}',\boldsymbol{p}',t') =0, \end{align}

where the two-particle distribution function is presented in accordance with definition (4.10), but it includes the dependence on $t'$:

(7.6)\begin{equation} f_{2}(\boldsymbol{r},\boldsymbol{p},t, \boldsymbol{r}',\boldsymbol{p}',t')=\int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \, {\rm d} \boldsymbol{\xi}' \,{\rm d}\boldsymbol{\eta}' \sum_{i=1}^{N/2}\sum_{j=1, j\neq i}^{N} \delta_{\boldsymbol{r}i} \delta_{\boldsymbol{p}i} \delta_{\boldsymbol{r}'j} \delta_{\boldsymbol{p}'j}. \end{equation}

Here, we make transition to the electromagnetic field instead of the integral form of the kinetic equation

(7.7)\begin{align} & \partial_{t}f(\boldsymbol{r},\boldsymbol{p},t)+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f +\frac{q_{s}}{m_{s}}\frac{1}{\varDelta}\frac{1}{\varDelta_{p}} \nonumber\\ & \quad \int_{\varDelta,\varDelta_{p}} \, {\rm d} \boldsymbol{\xi} \,{\rm d}\boldsymbol{\eta} \sum_{i=1}^{N/2} \left(\boldsymbol{E}(\boldsymbol{r},t) +\frac{1}{c}[\boldsymbol{v}\times \boldsymbol{B}(\boldsymbol{r},t)]\right) \delta_{\boldsymbol{r}i} \boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}}\delta_{\boldsymbol{p}i}=0, \end{align}

where $\boldsymbol {E}(\boldsymbol {r},t)$ and $\boldsymbol {B}(\boldsymbol {r},t)$ are the full fields composed of the external fields (like $\boldsymbol {E}_{\rm ext}(\boldsymbol {r},t)$) and the fields of interparticle interaction (like $\boldsymbol {E}_{\rm int}(\boldsymbol {r},t)$) $\boldsymbol {E}(\boldsymbol {r},t)=\boldsymbol {E}_{\rm ext}(\boldsymbol {r},t)+\boldsymbol {E}_{\rm int}(\boldsymbol {r},t)$ and $\boldsymbol {B}(\boldsymbol {r},t)=\boldsymbol {B}_{\rm ext}(\boldsymbol {r},t)+\boldsymbol {B}_{\rm int}(\boldsymbol {r},t)$. The field of interaction satisfies the Maxwell equations: $\boldsymbol {\nabla } \boldsymbol {\cdot }\boldsymbol {B}_{\rm int}=0$,

(7.8a,b)\begin{gather} \boldsymbol{\nabla}\times \boldsymbol{E}_{\rm int}={-}\frac{1}{c}\partial_{t}\boldsymbol{B}_{\rm int}, \quad \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{E}_{\rm int}=4{\rm \pi} \sum_{s}q_{s}\int f(\boldsymbol{r},\boldsymbol{p},t)\, {\rm d} \boldsymbol{p}, \end{gather}

(7.9)\begin{gather} \boldsymbol{\nabla}\times \boldsymbol{B}_{\rm int}=\frac{1}{c}\partial_{t}\boldsymbol{E}_{\rm int}+\frac{4{\rm \pi} }{c}\sum_{s}q_{s}\int \boldsymbol{v} f(\boldsymbol{r},\boldsymbol{p},t)\, {\rm d} \boldsymbol{p}. \end{gather}

Finally, the full set of the Vlasov–Maxwell equations is obtained in the self-consistent field approximation for the relativistic plasmas.