Hydrodynamic and kinetic representation of the microscopic classic dynamics at the transition on the macroscopic scale

Abstract

An open problem of the derivation of the relativistic Vlasov equation for systems of charged particles moving with velocities up to the speed of light and creating the electromagnetic field in accordance with the full set of the Maxwell equations is considered. Moreover, the method of derivation is illustrated on the non-relativistic kinetic model. Independent derivation of the relativistic hydrodynamics is also demonstrated. The key role of these derivations of the hydrodynamic and kinetic equations includes the explicit operator of averaging on the physically infinitesimal volume suggested by L.S. Kuzmenkov.

1. Introduction

The self-consistent field approximation, or the mean-field approximation, is well established in plasma physics. The Vlasov equation and corresponding hydrodynamics (the non-relativistic five moments, the non-relativistic thirteen moments and the relativistic five moments approximations) are considered in the literature to study the collective plasma phenomena developing on the time scale, where the contribution of ‘collisions’ is negligible.

The derivation of the Vlasov equation (Vlasov 1938; Weinberg 1972; Akhiezer 1975; Landau & Lifshitz 1980; Aleksandrov, Bogdankevich & Rukhadze 1984) based on the atomic structure of matter (composition of plasmas as the number of individual electrons and ions) was developed by N.N. Bogoliubov in 1946 (let us make references to the representation of this work in textbooks (Akhiezer 1975; Landau & Lifshitz 1980)). This derivation is made for the non-relativistic motion of charged particles interacting via the Coulomb interaction. So, the equations of field are reduced to the Poisson equation. Further generalization of this approach is made (Zaslavskii 1966; Pavlotskii 1973; Orlov & Pavlotsky 1989), where the weakly relativistic effects are considered. It includes the account of the interaction of moving charges via the magnetic field created in accordance with the Biot–Savart law. Hence, the problem of a kinetic description of particles that create an electromagnetic field in accordance with the full Maxwell equations is an open problem of existing kinetic theory. Some important steps in the direction of solving this problem are made in the works of Y.L. Klimontovich and L.S. Kuzmenkov. Y.L. Klimontovich suggested the microscopic form of the number density (and other hydrodynamic functions) and the distribution function, which are constructed of the Dirac delta functions (see for instance Klimontovich 1986). Y.L. Klimontovich also proposed a kinetic equation for the microscopic distribution function. However, the complete derivation requires the transition to the macroscopic scale. Our method of the explicit transition on the macroscopic scale via the integral operator is suggested by L.S. Kuzmenkov in a set of works (Drofa & Kuz'menkov 1996; Kuz'menkov & Andreev 2012; Kuz'menkov 2015) (see also an earlier work on kinetic equations Kuz'menkov 1991).

There are several definitions of the average value in classical statistical physics (Klimontovich 1986). They are the time average and the phase average. There is also the ergodic hypothesis, which states that these averages give the same result. The phase average can give different limit regimes depending on the number of controlled parameters. If we know the initial conditions of all particles in each system of the ensemble, we get a full dynamic description of each system. In this regime, we do not need to use the ensemble of physical systems at all. We can consider the dynamical evolution of the single system.

This regime, to some extent, resembles the analysis given in this paper. However, in spite of the similarity, there is a fundamental difference. If we consider any physical system, it does not ‘know’ which level of knowledge of this system we have. Our goal is to consider the evolution of this single system. This evolution is governed by the interaction between particles keeping the cause and effect consequence. Analysis is made for the structureless particles as the electrons and protons (we do not need to consider such microscopic particles as quarks on the chosen energy scale). If there is some structure which can be modified on the chosen energy scale, like the ionization of an atom or ion, excitation of electron inside of ion, etc., we have additional dependence on the structure of particles. It would make the model more complex, so we consider the ‘structureless’ electrons and protons. Hence, we introduce some microscopic functions. For instance, the number density if we want to study the evolution in the coordinate physical space or the distribution function if we want to study the evolution in the six-dimensional coordinate-momentum phase space.

The number density is the deterministic scalar field. Being the field means that this is the distribution of the physical parameter (characteristic) in the physical space arithmetized by $\boldsymbol {r}$. The change of this distribution in time corresponds to the motion of particles via the dependence on coordinates of particles $\boldsymbol {r}_{i}(t)$ on time $t$. This function describes the exact evolution of the system, but this evolution is given via the collective variables describing the whole system. The hydrodynamics itself can be used instead of the mechanical laws of motion for each particle, but the number of required hydrodynamic functions (the number density, the three projections of the velocity field, etc.) should be equal to the number of degrees of freedom of the system.

Large-scale dynamics, including the hydrodynamic and kinetic approaches, of interacting particles is considered in the book by Spohn (1991) relying on the probabilistic methods of description of the many-particle systems.

A direct derivation of both Debye shielding and some collective effects are demonstrated by Elskens, Escande & Doveil (2014), where the many-body dynamics is considered with no transition on the macroscopic scale or usage of average functions. A unified N-body description of Debye shielding and Landau damping is given by Escande, Doveil & Elskens (2016), where the action on the single particle of its own field is included. Multiple interpretations of Landau damping are given in the literature. It includes the mechanical N-body description of this damping (Escande et al. 2018), which enables its rigorous and simple calculation. It also gives its interpretation as the synchronization of almost resonant passing particles. This concept shows that the Debye shielding results from a cooperative dynamical self-organization process (Escande et al. 2018). The N-body approach incorporates spontaneous emission naturally, whose compound effect with Landau damping drives a thermalization of Langmuir waves (Escande et al. 2018). Comparison of the time evolution of the one-particle density of the N-particle Vlasov–Maxwell system with the Vlasov–Poisson equation is presented by Chen et al. (2020). It demonstrates the closeness of both time evolutions for number of particles $N$ and speed of light $c$ being large enough (Chen et al. 2020). The convergence in any time interval of a point-particle approximation of the Vlasov equation by particles initially equally separated for a force in $1/|x|^{\alpha }$ is demonstrated by Hauray & Jabin (2007), with $\alpha \leq 1$. The mean-field convergence for systems of points evolving along the gradient flow of their interaction energy for the Coulomb interaction is established by Serfaty (2020). A review of some classical results on the mean-field limit and propagation of chaos for systems of many particles, leading to the Vlasov or macroscopic equation, is presented by Jabin (2014).

A probabilistic proof of the mean-field limit, and propagation of chaos N-particle systems and convergence of the empirical distributions to solutions of the Vlasov–Poisson system are given by Lazarovici & Pickl (2017). Some mathematical properties of probabilistic justification of the mean-field approach is given by Dobrushin (1979). Dobrushin (1979) and Lazarovici & Pickl (2017) deal with the methods of description of the random fields, while we try to construct the macroscopic field following the cause-effect consequence.

The kinetic equations of Vlasov theory, in the weak formulation, are rigorously shown to govern the infinite number of particles limit of the Newtonian dynamics (Kiessling 2014). The proof is based on the Liouville equation, more precisely, the first member of the pertinent BBGKY hierarchy, and does not invoke the Hewitt–Savage theorem, nor any regularization of the interactions (Kiessling 2014). A non-trivial electromagnetic vacuum is considered by Elskens & Kiessling (2020) to model the radiation reaction and include it in the relativistic Vlasov–Maxwell equations via N-particle dynamics. Some mathematical concepts regarding the mean-field limit of an N-particle system towards a regularized variant of the relativistic Vlasov–Maxwell system are presented by Golse (2012). Some mathematical aspects related to the microscopic derivation of Vlasov equations with singular potentials are discussed by Golse (2022) and Grass (2021), where the transition to the mean-field in the quantum theory is discussed as well. These works do not include the concept of construction of the material fields based of motion of classical particles suggested by Drofa & Kuz'menkov (1996) and developed in this paper.

The relativistic plasmas are actively studied for the classical and quantum regimes (Hakim & Sivak 1982; Kuz'menkov 1991; Hakim et al. 1992; Shatashvili, Javakhishvili & Kaya 1997; Shatashvili & Rao 1999; Hazeltine & Mahajan 2002; Mahajan & Hazeltine 2002; Melrose 2008; Melrose & Weise 2009; Romatschke 2010; Asenjo et al. 2011; Bret & Haas 2011; Hakim 2011; Mahajan & Yoshida 2011; Mendonca 2011; Asenjo et al. 2012; Melrose & Weise 2012; Zhu & Ji 2012; Comisso & Asenjo 2014; Ivanov, Andreev & Kuz'menkov 2014,  2015; Ruiz & Dodin 2015; Shatashvili, Mahajan & Berezhiani 2020). Therefore, the detailed analysis of the derivation of the classic hydrodynamics and kinetics is important for the better understanding of these physical processes. Recently, some steps in the direction of analysis of the relativistic hydrodynamics have been made. However, the contribution of the multipole moments of the physically infinitesimal volume is ignored in previous papers on this subject (Andreev 2022a, 2023a, 2023b) (see also Refs. Andreev 2022b, 2022c, 2023c for the adaptation of this approach for the degenerate plasmas). The multipole moments and equations for their evolution are partially described by Drofa & Kuz'menkov (1996) for the non-relativistic regime, where some nonlinear phenomena are also discussed. The systematic description of the appearance of the multipole moments in the hydrodynamic and kinetic equations is presented in this paper.

This paper is organized as follows. In § 2, the non-relativistic systems of charged particles are considered on the microscopic scale and corresponding hydrodynamic equations are obtained to make the transition on the language suitable for the description of the collective phenomena. In § 3, the self-consistent field approximation is considered for the non-relativistic hydrodynamic equations. In § 4, the non-relativistic kinetic Vlasov equation for the systems of charged particles is derived. In § 5, the self-consistent field approximation is considered for the non-relativistic kinetic equations. In § 6, the analysis of the self-consistent field approximation is made for the relativistic hydrodynamic model with the average reverse gamma factor evolution. In § 7, the relativistic kinetic Vlasov equation is derived and the self-consistent field approximation is considered for the relativistic kinetic Vlasov equation. In § 8, a brief summary of the obtained results is presented.

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 (1996) and Kuz'menkov (1991) 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 (1986) in his book (see also Klimontovich 1967) 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 1996)

(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 1996; Kuz'menkov & Andreev 2012)

(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 (1967) (see also Weinberg 1972). 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 1980).

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 1975; Landau & Lifshitz 1980; Aleksandrov et al. 1984).

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 (1996).

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 (1996), 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 (2011) and Li, Hwang & Sarma (2011), 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 1986):

(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 1986):

(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 (2022a, 2023b). 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 1990; Spohn 2004) 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 (2022a, 2023b). 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 (2022a, 2023b). 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 (2022a) 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 2022a). 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 (2022a, 2023b):

(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 2022a, 2023b):

(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 2022a). Equation (6.41) shows that we need to include the third equation of state for the function $\tilde {M}$. It is found in Andreev (2022a). 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 2022a).

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.

8. Conclusion

A problem of derivation of the Vlasov kinetic equation for the full relativistic regime has been considered, where the electromagnetic field created by each particle satisfies the full set of Maxwell equations. It has been included along with the high temperatures of plasmas. A similar problem has been considered for the relativistic hydrodynamics. Both problems have been addressed in this paper and the method of solving these problems has been demonstrated explicitly, including some essential technical details.

To simplify the presentation, the paper has been split into several sections. First, the non-relativistic hydrodynamics has been considered. So, the interparticle interaction has been reduced to the Coulomb interaction. It has been derived from the microscopic motion of classic particles. The whole paper has been focused on classical systems with no discussion of the quantum effects. The presented method of derivation includes the transition on the macroscopic scale from the microscopic level of description. Hence, the physically infinitesimal volume has been presented analytically. Dealing with the microscopic description, we consider the finite elements of volume as the macroscopically point-like objects. However, these volumes (which are called $\varDelta$-neighbourhoods) are characterized by the dipole moment density, the quadrupole moment density, etc., in addition to the charge density. The contribution of the multipole moments in the Euler equation and the Poisson equation has been illustrated. The equations for the multipole moments evolution are not presented since our goal in this paper is to give a background for the further derivation of the relativistic hydrodynamic and kinetic models. Hence, the features related to the multipole moments have been mentioned, but no stress has been made on this item. The self-consistent field approximation (the mean-field approximation) has been discussed for the non-relativistic plasmas as well (§ 3).

The second part of this paper has been focused on the derivation of the Vlasov kinetic equation, in the non-relativistic regime, where particles interact via the Coulomb interaction. The derivation includes the transition on the macroscopic scale both in the coordinate and momentum space. Hence, the distribution function for the multipole moments have been found. No additional kinetic equations have been considered for these functions, but their presence in the Vlasov equation and the Poisson equation has been highlighted. Neglecting the multipole expansion of the Coulomb interaction, we still have an additional vector distribution function. It can be characterized as the ‘dipole moment of the neighbourhood in the momentum space’. If we neglect the contribution of all additional functions, we find the kinetic equation with the two-particle distribution function. Further application of the self-consistent field approximation leads to the well-known Vlasov equation in the Coulomb approximation.

The third item is one of two major items in this paper. It shows the derivation of the relativistic hydrodynamics with no analysis of the radiation reaction. Let us repeat that the suggested method of derivation of hydrodynamic and kinetic models shows the method of explicit transition to the macroscopic scale. It has been called the method of averaging, but it has no statistical or probabilistic meaning. It gives the transition of the deterministic behaviour on another scale of space parameters. The suggested method allows to derive the relativistic hydrodynamics in the well-known form of the set of continuity equation and four momentum evolution equation (see for instance Andreev 2023b). So, this model will include the multipole expansion demonstrated in this paper. However, another form of the relativistic hydrodynamic model is chosen for the analysis. It is the relativistic hydrodynamic model with the average reverse gamma factor evolution. The continuity equation and the equation for the evolution of the current of particles, which transforms to the velocity field evolution equation, are derived in full details. The evolution of the average reverse gamma factor and its current demonstrate similar structures, but they have not been demonstrated. The derivation includes two essential elements. First is the relativistic temperatures of the plasmas. Second is the full relativistic interaction between particles, so the electromagnetic field created by each particle satisfies the full set of Maxwell equations. Let us point out that in the non-relativistic systems, the electric and magnetic field are bound to the particles. Therefore, the system is described by the degrees of freedom of the particles. In the relativistic regime, the electromagnetic field has its own infinite number of degrees of freedom. However, the 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 the integral form of the Maxwell equations.

The fourth item is the full relativistic derivation of the Vlasov kinetic equation. This element of the paper is the major part from the fundamental point of view. The derivation itself has been presented via a few equations. From a technical point of view, it is similar to the derivation of the non-relativistic version. However, it is essential that the described method allows to give the derivation at the relativistically large temperatures while the interparticle interaction happens in accordance with the full set of Maxwell equations at the microscopic level. So, the self-consistent macroscopic field also satisfies the full set of Maxwell equations. Moreover, this derivation is rather simple and straightforward. So, it is easy to understand.

The described derivation includes some open problems for the hydrodynamics and kinetics. A major problem is the complete analysis of the multipole moments and their dynamics. Corresponding equations for the evolution of these functions are not discussed. Moreover, a closed set of equations consistently describing these effects, especially in kinetics, is to be found.