Reduced description method in the kinetic theory of Brownian motion with active fluctuations

Active matter is a field in soft matter physics that studies the properties of aggregates of self-propelled objects that have the ability to take up energy from the environment, store it in an internal depot, and convert internal energy into kinetic energy. Active matter comprises diverse systems spanning from macroscopic (e.g. schools of fish and flocks of birds) to microscopic scales (e.g. migrating cells, motile bacteria, and gels formed through the interaction of nanoscale molecular motors with cytoskeletal filaments within cells). Here we refer the reader to the reviews [1–4] and references therein. In recent years the number of publications devoted to theoretical and experimental studies of the phenomena in active substances has sharply increased, and some recent achievements have been published in the Special Issues of the European Physical Journal Special Topics and the Journal of Statistical Physics [5–9]. The study of nonequilibrium processes in systems of active particles inevitably raises the question of the consistent derivation of the evolution equations for such systems, in particular, the kinetic equations. The derivation of a kinetic theory of active particles is a challenging issue that has attracted significant attention over the past decades. Bertin et al [10–12] has derived the Boltzmann equation for self-propelled point-like particles on a two-dimensional plane with the assumption that the modulus of the velocity vector is fixed and identical for all particles so that only the direction of the vector plays a role in the dynamics. Ihle [13, 14] developed an alternative kinetic approach that is based on the Chapman–Kolmogorov equation for the N-particle probability density. The resulting mean-field kinetic equation has been studied analytically and numerically, and has been extended to the so-called topological interactions [15–17]. Romanczuk et al [18–20] derived and explored the mean-field kinetic equation in two spatial dimensions, starting from the Langevin equation with active friction and active fluctuations, and supplemented with different forces describing the interaction between the particles. In [21], the authors also pursue the Langevin approach to study collective dynamics in the two-dimensional system of active Brownian particles with dissipative interactions.

In this paper we develop a consistent microscopic approach based on the reduced description method of relaxation processes. The reduced description method of Bogolyubov–Peletminskii has a number of notable advantages. In particular, it allows one to formulate approaches for the construction of the kinetic theory of systems of many identical particles, and thus, in essence, provides a dynamic justification of the statistical mechanics of such systems [22, 23]. Bogolyubov suggested a method of reduced description of the evolution of many-particle systems [22], which allowed the construction of a regular procedure for obtaining closed dissipative kinetic equations based on the Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY) chain of reversible equations for many-particle distribution functions. Fundamentals of the reduced description method were formulated in [22] for the classical (nonquantum) systems of many particles. In the case of quantum many-particle systems, the ideas of the Bogolyubov reduced description method were developed in the works by Peletminskii, and the main results were presented in [23]. There are also other approaches for the dynamical justification of statistical mechanics that are different from Bogolyubov's approach; for example, in the works of Prigogine's Brussels school [24], as well as the different formulations of Bogolyubov's ideas, see, e.g. [25–28]. In the present paper we use the reduced description method in the form close to the one of Bogolyubov–Peletminskii [29] to construct the kinetic theory of many-particle systems with active fluctuations and nonlinear friction. For that purpose we need to generalize the canonical Bogolyubov–Peletminskii approach in order to take into account an external stochastic impact and dissipative interactions.

A generalization of the Bogolyubov reduced description method to the case of dissipative many-particle systems in an external stochastic field was first suggested in [30]. In this paper, the authors proposed a formalism for deriving kinetic equations. As a starting point, a stochastic Liouville equation obtained from Hamilton's equations taking dissipation and stochastic perturbations into account was used. The Liouville equation is then averaged over realizations of the stochastic field by an extension of the Furutsu–Novikov formula to the case of a nonGaussian field. As a result, a generalization of the classical Bogolyubov–Born–Green–Kirkwood–Yvon hierarchy is derived. In order to get a kinetic equation for the one-particle distribution function, the authors use a regular breaking procedure of the BBGKY hierarchy by assuming weak interaction between the particles and weak intensity of the field. Within this approximation they get the corresponding Fokker–Planck equation for the system in a nonGaussian stochastic field. Two particular cases assuming either Gaussian statistics of external perturbation or homogeneity of the system are discussed. In that approach, however, the stochastic external forces do not depend on the velocity (or momentum) of the particles. In other words, the formalism developed in [30] can be applied to systems with nonlinear friction, as is the case of active particle systems, but with passive fluctuations of either Gaussian or nonGaussian nature.

As mentioned above, we announce this work as one in which a general microscopic approach is proposed based on the first principles of statistical physics for the construction of the kinetic theory of such complex systems as dissipative media with active fluctuations. These media are proposed to be considered in the manuscript as systems of many identical particles. The structural units of the system are assumed to be point-like, but capable of reacting to an external stochastic effect depending on the magnitude and direction of their own velocity. This is their main individual characteristic, indicating that the point particles have the property of the so-called head–tail asymmetry. In particular, the manifestation of the self-propulsion property of the medium itself can be associated with just this simplest characteristic of the asymmetry of structural units. In this paper, we show that in the studied model of the active medium the self-propulsion property can manifest itself as a collective effect; however, its very manifestation is possible only due to the individual properties of the structural units of the system, i.e. the active particles.

Actually, one of the main goals of this work was to demonstrate both the necessity and the very possibility of developing a microscopic approach for building the kinetics of such complex systems with as accurate formulations as possible. In addition, we wanted to unite the existing phenomenological approaches (if possible, without contradictions) with their generalization to a number of more complex systems with active fluctuations. First of all, this means the correct inclusion of potential and dissipative interactions between the structural units of the systems under study. In parallel, we aimed to take into account the possibility of influencing the system of spatially correlated noise depending on the absolute value and the direction of the velocity of the system's structural unit. At the same time, the formulations should allow for a generalization of the approach to the case of nonGaussian noise. All stated goals could be achieved only within the framework of a microscopic approach based on the first principles of statistical physics. In this sense, we consider reference works [22–27] with a preference for strict approaches and elegant formulations of the monograph [23]. This means that the scheme for constructing the microscopic approach should contain a series of consecutive, strictly controlled steps. In particular, based on the individual equations of motion for the structural units of the system, the Liouville equation must be obtained for the distribution function of all the particles of the system. From the Liouville equation, the transition must be made to the equations of motion for many-particle distribution functions, which will have the form of an infinite chain of kinetic equations, which, in essence, is a generalization of the BBGKY chain to the case of systems with active fluctuations. The next step in the procedure for constructing a kinetic theory is the controlled truncation of this chain in any perturbation theory. In the present article we choose the weak interaction between the structural units of the system and the low intensity of the external stochastic force as the small parameters for such a perturbation theory. In this case, the density of the system is not necessarily small, and the closed kinetic equation for the single-particle distribution function has the form of the Fokker–Planck equation. Recall that, in accordance with [23], it is also possible in perturbation theory to break an infinite chain of equations—an analogue of the BBGKY chain and in the case of low medium density and with arbitrary interaction between particles, but only if this interaction does not lead to the formation of bound states. In this case, we would inevitably come to some analogue of the Boltzmann kinetic equation. However, for the active media under study, the possible (in principle) implementation of this case leads to significant mathematical difficulties and requires separate efforts (see this in regards [31]). Note also that in the transition from the Liouville equation to an infinite chain of equations, an additional intermediate step is required. It is connected with the necessity of averaging the Liouville equation, obtained from the initial equations of motion for structural units, over an external stochastic field. In this case, the work in [30] described above becomes useful, and its technique is easily modified for the systems studied in this article.

Building a consistent microscopic approach requires, first of all, to go beyond the Langevin equations as the initial equations of motion for the structural units of the system (see in this regard [30]). Instead of the Langevin equations, it is necessary to formulate the generalized Hamilton equations, that is, to start from the first principles of classical theoretical mechanics. Within the framework of such formulations it is possible to correctly include the interaction between structural units, both potential and dissipative. In this case, there are no problems with the inclusion of velocity-dependent noise correlated both in time and in coordinates. Note that the proposed theory does not deny the possibility of incorporating the Langevin equations as the initial equations of motion for structural units of the active medium. However, it will not be possible to talk about a rigorous microscopic theory that comes from first principles. In fact, the Langevin equations themselves were originally obtained as phenomenological equations of motion (Newton's equations), in which the force of the medium on a selected (test) particle is divided into a regular friction force and some random force. Thus, the Langevin equations are formulated for obviously noninteracting particles. For this reason, the inclusion of interaction between particles in these equations is again associated with the involvement of phenomenological approaches. Note also that in the 'traditional' Langevin equations, the noise acting on Brownian particles depends only on time and not on the coordinates. The presence of spatially correlated noise in the Langevin equations in the existing methods for deriving the Fokker–Planck kinetic equation would lead to serious mathematical difficulties. Thus, in accordance with the requirements listed above, the task of describing kinetic processes in the systems under study is posed in the most general form possible.

The solution of this problem in the present work is also carried out in the most general and consistent form, adhering to the first principles of statistical physics. First of all, this refers to the derivation of general kinetic equations for a system of active particles. Of course, a number of successive simplifications were used in this paper, connected with the transition from more complex systems, for which more complex equations were derived. In particular, we completely refused in the present work the generalization of the proposed approach in the case of nonGaussian noise. First of all, it was done for the sake of reducing the volume of the article. In addition, such a simplification is possible because, in light of [30], the ways of modifying the approach proposed in this article to the case of nonGaussian noise seem quite understandable and do not contain fundamental difficulties. The paper also carried out a fairly detailed analysis of some more simple specific examples, designed to demonstrate the effectiveness of the proposed method and the kinetic equations obtained in terms of the implementation of the announced goals. Such efficiency is most easily illustrated by analyzing cases in which the results of our theory coincide with or are close to the results of other authors obtained in other approaches. Section 6 of this paper is fully devoted to the solution of the latter problem, in which the results of the analysis of some particular solutions of our equations are given.

Thus, as the result of the consistent implementation of the above tasks, we suggest a generalized formulation of the reduced description method, suitable for describing the kinetics of many-particle dissipative systems with active fluctuations. It is shown in the framework of the microscopic approach developed that it is possible to construct the kinetic theory of active particles in both two- and three-dimensional systems, with the availability of nonlinear friction (dissipative interaction), as well as with local impact of an external active random field. Under the local impact we assume that this field may act differently at different points in space. In other words, the effect of this field on a particle may depend not only on the velocity (or momentum) of that specific particle, but on the point in the coordinate space where the particle is located. The general kinetic equations for such systems are obtained. We also consider special cases in which the obtained kinetic equations give solutions known for the active particles from earlier works [3, 18, 20, 21].

Consider a system consisting of N identical active particles of mass m, each of which is characterized by spatial coordinates ${\mathbf x}_{\alpha}$, $1 \leqslant \alpha \leqslant N$, measured from the center of mass, and momentum ${\mathbf p}_{\alpha}$, $1 \leqslant \alpha \leqslant N$. The interaction between the particles is assumed to consist of two parts—a 'reversible' part described by the Hamiltonian H, and an 'irreversible' one, described by the function $R^{\omega}$, the meaning of which will be explained below.

The Hamiltonian of the system can be written as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_1_} H=H_{0} +V=\sum_{1\leqslant \alpha \leqslant N}\frac{{\mathbf p}_{\alpha }^{2} }{2m} +\sum_{1\leqslant \alpha <\beta \leqslant N}V_{\alpha ,\beta }, \nonumber \end{align} \tag{ 1 }$$

where $V_{\alpha ,\beta }$ is the pair interaction potential,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_2_} V_{\alpha ,\beta } \equiv V\left({\mathbf x}_{\alpha \beta } \right), \nonumber \\ {\mathbf x}_{\alpha \beta } \equiv {\mathbf x}_{\alpha } -{\mathbf x}_{\beta }. \nonumber \end{align} \tag{ 2 }$$

We also assume that the particles of the system are exposed to specific forces that depend on the particle velocity (or momentum) and are characterized by a function R. We assume that the function R can be represented as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_3_} R=R^{{\rm r}} +R^{\omega }, \nonumber \end{align} \tag{ 3 }$$

where $R^{{\rm r}}$ is a regular part of this function.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_4_} \nonumber & R^{{\rm r}} \equiv \sum\limits_{1\leqslant \alpha <\beta \leqslant N}R_{\alpha ,\beta } ,\nonumber \\ \nonumber & R_{\alpha ,\beta } \equiv R\left({\mathbf x}_{\alpha \beta } ,{\mathbf p}_{\alpha \beta } \right),\nonumber \\ & {\mathbf p}_{\alpha \beta } \equiv {\mathbf p}_{\alpha } -{\mathbf p}_{\beta } , \nonumber \end{align} \tag{ 4 }$$

and $R^{\omega }$ is a stochastic part of the function R, which can be written as

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle R^{\omega } \equiv \sum_{1\leqslant \alpha \leqslant N}R^{\omega } \left(x_{\alpha } ,t\right) , \, x_{\alpha } \equiv \left\{{\mathbf x}_{\alpha } ,{\mathbf p}_{\alpha } \right\}.\nonumber \end{align*}$$

The stochastic nature of the function $R^{\omega}$ is formally highlighted by the presence of index $\mathrm{\omega}$.

Note that in the case of nonactive identical particles with the dissipative interaction function, $R^{{\rm r}}$ is treated as a dissipative function (see [32] and [30, 31]). It is usually assumed that the dissipation in the system is related to friction of macroscopic particles. Therefore, in this case the dissipation function $R_{\alpha, \beta}$, following [32], can be chosen as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_5_} R_{\alpha ,\beta } \equiv \frac{1}{2} \tilde{\gamma }\left({\mathbf x}_{\alpha \beta } \right){\mathbf p}_{\alpha \beta }^{2},\qquad \tilde{\gamma }\left({\mathbf x}_{\alpha \beta } \right)>0,\qquad {\mathbf p}_{\alpha \beta } \equiv {\mathbf p}_{\alpha } -{\mathbf p}_{\beta }. \nonumber \end{align} \tag{ 5 }$$

This implies that $\tilde {\gamma} \left ({\mathbf x}_{\alpha \beta} \right) = 0$ if $\left | {\mathbf x}_{\alpha \beta} \right | \mathop {>} \limits_{\sim} r_{0}$, where r₀ is a characteristic range of dissipative forces. In view of the property (5) the friction coefficient is always positive.

However, in the case of active particles the positivity does not always hold [3]. The friction coefficient in the Langevin-type equations for active particles can depend on the velocity and change its sign. Therefore, one cannot use the criteria (5) to determine the properties of dissipation function $R_{\alpha, \beta}$ in the case of active particles.

Following the usual classical theoretical mechanics procedures, and taking into account equations (1)–(4), the generalized Hamilton equations for the system under study can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_6_} \dot{{\mathbf p}}_{\alpha } =-\frac{\partial H}{\partial {\mathbf x}_{\alpha } } -\frac{\partial R}{\partial \mathbf p_{\alpha } }, \qquad \dot{{\mathbf x}}_{\alpha } =\frac{\partial H}{\partial {\mathbf p}_{\alpha } }. \nonumber \end{align} \tag{ 6 }$$

Thus, the force ${{\mathbf F}}_{\alpha, \beta}$ acting on a particle $\alpha$ from the particle $\beta$ consists of two terms:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_7_} {\mathbf F}_{\alpha ,\beta } = {\mathbf F}_{\alpha ,\beta }^{{\rm p}} + {\mathbf F}_{\alpha ,\beta }^{{\rm r}}. \nonumber \end{align} \tag{ 7 }$$

Namely, the force ${{\mathbf F}}_{\alpha, \beta}^{{\rm p}}$, connected with the presence of a potential pair interaction between the particles, and the force ${\mathbf F}_{\alpha, \beta}^{{\rm r}}$, connected with the presence of a dissipative interaction between the particles (in the sense outlined above):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_8_} {\mathbf F}_{\alpha ,\beta }^{{\rm p}} \equiv -\frac{\partial V_{\alpha ,\beta } }{\partial \mathbf x_{\alpha } }, \qquad {\mathbf F}_{\alpha ,\beta }^{{\rm r}} \equiv -\frac{\partial R_{\alpha ,\beta } }{\partial {\mathbf p}_{\alpha } }. \nonumber \end{align} \tag{ 8 }$$

In addition, the $\alpha$th particle is influenced by external random force ${\mathbf Y}_{\alpha}^{\omega}$ (see equation (1), which depends on the momentum of the particle, wherein

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_9_} &\qquad-\frac{\partial R\left(x_{\alpha } ,t\right)}{\partial p_{i\alpha } } \equiv Y_{i}^{\omega } \left(x_{\alpha } ,t\right)=X_{\parallel i}^{\omega } \left( {\mathbf x}_{\alpha }, {\mathbf p}_{\alpha } ,t\right) +X_{\bot i}^{\omega } \left( {\mathbf x}_{\alpha }, {\mathbf p}_{\alpha } ,t\right), \nonumber \\ &\nonumber X_{\parallel i}^{\omega } \left( {\mathbf x}_{\alpha }, {\mathbf p}_{\alpha } ,t\right) \equiv X^{\omega } \left( {\mathbf x}_{\alpha }, t\right) {e} _{\alpha i}^{\rm h}, \quad X_{\bot i}^{\omega } \left( {\mathbf x}_{\alpha }, {\mathbf p}_{\alpha } ,t\right) \equiv X_j^{\omega } \left( {\mathbf x}_{\alpha }, t\right)\left( \delta_{ij} - {e}_{\alpha i}^{\rm h} {e}_{\alpha j}^{\rm h} \right), \quad {e}_{\alpha i}^{\rm h} \equiv \frac{p_{\alpha i} }{\left| {\mathbf p}_{\alpha } \right|}. \nonumber \end{align} \tag{ 9 }$$

The latter expression requires some comments. We emphasize, first of all, that the stochastic force $Y_{i}^{\omega} \left (x_{\alpha}, t \right)$ in equation (9) is written in a form that is not related to the choice of a particular coordinate system. This notation simply reflects the fact that the stochastic force acts differently along and across the direction of a particle velocity. Figure 1 shows a schematic visualization of the assumption written in the form (9). Expression (9) looks like a natural extension of the stochastic force $Y_{i}^{\omega} \left ( {\mathbf x}, t \right)$ typical for the Langevin equation in the case of an ordinary Brownian particle:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \dot{{\mathbf x}}= {\mathbf v}, \qquad \dot{v}_{i} =-\gamma v_{i} +Y_{i}^{\omega } \left( {\mathbf x},t\right).\nonumber \end{align*}$$

Zoom In Zoom Out Reset image size

In fact, the value $Y_{i}^{\omega} \left ( {\mathbf x} , t \right)$ in the last equation can always be identically rewritten as

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle Y_{i}^{\omega } =\left( {\mathbf e}^{\rm h} {\mathbf Y}^{\omega } \right)e_{i}^{\rm h} +Y_{j}^{\omega } \left(\delta_{ij} -e_{i}^{\rm h} e_{j}^{\rm h} \right),\nonumber \end{align*}$$

where $e_i^{\rm h }$ is an arbitrary unit vector, for example, $e_{i}^{\rm h} = p_i/| {\mathbf p}|$. Replacing here the scalar product ${\mathbf e}^{\rm h} {\mathbf Y}^{\omega}$ with $X^{\omega} \left ( {\mathbf x}_{\alpha}, t \right)$, $Y_{j}^{\omega}$ by $X_{j}^{\omega} \left ( {\mathbf x}_{\alpha}, t \right)$, and assuming that $e_{\alpha i}^{\rm h} = p_{\alpha i} / | {\mathbf p}_\alpha |$, we arrive at equation (9). It should be remembered, however, that in equation (9) the values $X^{\omega} \left ( {\mathbf x}_{\alpha}, t \right)$, $X_{j}^{\omega} \left ( {\mathbf x}_{\alpha}, t \right)$ do not relate to each other in general. If necessary, in the three-dimensional case the vector $X_{j}^{\omega} \left ( {\mathbf x}_{\alpha}, t \right)$ can be considered as two-component in a plane perpendicular to the vector $e_{\alpha j}^{\rm h}$. The presence of components $X_{j}^{\omega} \left ( {\mathbf x}_{\alpha}, t \right)$ along $e_{\alpha j}^{\rm h}$ will not affect the description of processes and phenomena in such systems in any case because of the factor $\left (\delta_{ij} -e_{\alpha i}^{\rm h} e_{\alpha j}^{\rm h} \right)$ in the right-hand side of equation (9).

It follows from the above that the stochastic effects on the system under consideration in the form of equation (9) can be regarded as a generalization of stochastic forces used in the theory of two-dimensional systems of active particles, i.e. active fluctuations; see, e.g. [3, 18]. First, equation (9) allows for the possibility of local influence of stochastic forces on the system. Second, equation (9) can be applied to both two- and three-dimensional systems as well. To see this it is sufficient to consider expression (9) as two-dimensional and nonlocal:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_10_} {\mathbf e}_{\alpha }^{\rm h} \equiv {\mathbf e}_{\rm h} , X^{\omega } (t)\equiv \sqrt{D_{{\rm v}} } \xi_{{\rm v}} (t), \qquad \mathbf X^{\omega } \left( {\mathbf x}_{\alpha } ,t\right)\equiv {\mathbf e}_{\phi } \sqrt{D_{\phi } } \xi_{\phi } (t), \qquad {\mathbf e}_{\rm h} {\mathbf e}_{\phi } =0 , \nonumber \end{align} \tag{ 10 }$$

where ${\mathbf e}_{\rm h}$ is a unit vector along the direction of motion of a particle, ${\mathbf e}_{\phi}$ is a unit vector along the azimuthal angle $\phi$, and $D_{\phi}$ and $D_{{\rm v}}$ are angular and velocity noise intensities, respectively [3]. Note that in two-dimensional systems, as is known, the isolated directions can appear in the movement of active particles. The special importance of the direction of particle motion, set by the particle's orientation vector (head–tail axis) is due to the existence of the propulsion mechanism. Thus, due to the head–tail asymmetry in the steady state of the many active particles system it is possible to fix the reference system in a natural way by a special choice of the vectors ${\mathbf e}_{\rm h}$, ${\mathbf e}_{\phi}$. The existence of this asymmetry is reflected in the many-particle system characteristics, such as a one-particle distribution function. As shown below, the existence of the effects of head–tail asymmetry is also possible in three dimensions, even in the case of linear friction (see section 5 of this paper). We emphasize that the source of stochastic effects can be generalized to three dimensions in a form other than that in equation (9). Similar to the two-dimensional case, one can use the spherical coordinates. However, in this paper it is easier to employ the Cartesian coordinates.

Let us also note the following. The time derivative of the total energy of the system in accordance with equations (1) and (6) is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_11_} \frac{{\rm d}H}{{\rm d}t} =-\sum_{1\leqslant \alpha \leqslant N}\frac{{\mathbf p}_{\alpha } }{m} \frac{\partial R}{\partial \mathbf p_{\alpha } }. \nonumber \end{align} \tag{ 11 }$$

If we assume that the system has a dissipation due to friction of macroscopic particles, and a regular part of the dissipation function $R^{{\rm r}}$ is given by equation (5). Then, equation (11) together with equation (9) give

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_12_} \frac{{\rm d}H}{{\rm d}t} =-\frac{2}{m} \sum_{1\leqslant \alpha <\beta \leqslant N}R_{\alpha ,\beta } + \sum_{1\leqslant \alpha \leqslant N}\frac{p_{i\alpha } }{m} \left[X^{\omega } \left( {\mathbf x}_{\alpha } ,t\right)e_{\alpha i}^{\rm h} \right. \left. + X_{j}^{\omega } \left( {\mathbf x}_{\alpha } ,t\right)\left(\delta_{ij} -e_{\alpha i}^{\rm h} e_{\alpha j}^{\rm h} \right)\right], \nonumber \end{align} \tag{ 12 }$$

or

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm d}H}{{\rm d}t} =-\frac{2}{m} \sum_{1\leqslant \alpha <\beta \leqslant N}R_{\alpha ,\beta } +\sum_{1\leqslant \alpha \leqslant N}\frac{\left| {\mathbf p}_{\alpha } \right|}{m} X^{\omega } \left( {\mathbf x}_{\alpha } ,t\right) .\nonumber \end{align*}$$

Taking into account that $\tilde {\gamma} \left ( {\mathbf x}_{\alpha \beta} \right)> 0$ (equation (5)), competition between the dissipation due to friction and the pumping of energy from the stochastic field is possible in such a system.

A further task is to obtain the Liouville equation. To this end, and for the convenience of further calculations, we represent equation (6) in the following form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_13_} \dot{x}_{\alpha } (t)=h_{\alpha }^{\omega } \left(x_{1} (t),...,x_{N} (t)\right), \qquad 1\leqslant \alpha \leqslant N, \nonumber \end{align} \tag{ 13 }$$

where we introduce the notation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_14_} x_{a} (t)\equiv \left({\mathbf x}_{a} (t),\; {\mathbf p}_{a} (t)\right). \nonumber \end{align} \tag{ 14 }$$

Thus, equation (13) together with equations (14) and (5) read

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_15_} \dot{{\mathbf x}}_{\alpha } (t)= \mathbf h_{{\mathbf x}\alpha }^{\omega } \left(x(t)\right), \qquad \dot{{\mathbf p}}_{\alpha } (t)= \mathbf h_{{\mathbf p}\alpha }^{\omega } \left(x(t)\right), \nonumber \end{align} \tag{ 15 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_16_} \mathbf h_{{\mathbf x}\alpha }^{\omega } \left(x(t)\right)=\frac{\partial H}{\partial {\mathbf p}_{\alpha } }, \qquad \mathbf h_{{\mathbf p}\alpha }^{\omega } \left(x(t)\right)=-\frac{\partial H}{\partial {\mathbf x}_{\alpha } } -\frac{\partial R}{\partial {\mathbf p}_{\alpha } }. \nonumber \end{align} \tag{ 16 }$$

The coordinates and momenta of the $\alpha$th particle at time t (see equation (14)), are determined by the coordinates and momenta $\newcommand{\e}{{\rm e}} x_{0} \equiv \left (x_{1} \left (0 \right), ..., x_{N} \left (0 \right) \right)$ of all the particles at the initial time t  =  0:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_17_} x_{\alpha }^{\omega } (t)=X_{\alpha }^{\omega } \left(t,x_{0} \right)\equiv \left( \mathbf X_{\alpha }^{\omega } \left(t,x_{0} \right),\; \mathbf P_{\alpha }^{\omega } \left(t,x_{0} \right)\right), \nonumber \end{align} \tag{ 17 }$$

where the functions $\mathbf X_{\alpha}^{\omega} \left (t, x_{0} \right)$ and $\mathbf P_{\alpha}^{\omega} \left (t, x_{0} \right)$ satisfy the generalized Hamilton equation (5) (or equations (13)–(16)). We introduce the probability density $D \left (x_{1} \left (0 \right),\right.$ $\dots, x_{N} \left (0 \right);$ $\left. 0 \right)$ of the initial conditions $\newcommand{\e}{{\rm e}} x_{0} \equiv \left (x_{1} \left (0 \right), ..., x_{N} \left (0 \right) \right)$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_18_} \int {\rm d}x_{1} (0)...{\rm d}x_{N} (0)D\left(x_{1} (0),...,x_{N} (0);0\right)\equiv \int {\rm d}x_{0} D\left(x_{0} ;0\right) =1 . \nonumber \end{align} \tag{ 18 }$$

Then, at time t the probability density $\newcommand{\e}{{\rm e}} D^{\omega} \left (x_{1}, ..., x_{N}; t \right) \equiv D^{\omega} \left (x; t \right)$, $\newcommand{\e}{{\rm e}} x \equiv \left (x_{1}, ..., x_{N} \right)$, (N-particle distribution function) is defined by the expression

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_19_} D^{\omega } \left(x_{1} ,...,x_{N} ;t\right)=\int {\rm d}x_{0} D\left(x_{0} ;0\right) \mathop{\prod }\limits_{1\leqslant \alpha \leqslant N} \delta \left(x_{\alpha } -X_{\alpha }^{\omega } \left(t,x_{0} \right)\right). \nonumber \end{align} \tag{ 19 }$$

In [33] a detailed procedure for obtaining the Liouville equation for many-particle systems in an external stochastic field neglecting the interaction between the particles is described. In [30] (see also [33]), a similar procedure is used to obtain a generalized Liouville equation for dissipative many-particle systems in the absence of an external stochastic field. The N-particle distribution function obeys the continuity equation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_20_} \frac{\partial D^{\omega } }{\partial t} +\sum_{1\leqslant \alpha \leqslant N}\frac{\partial }{\partial x_{\alpha } } \left(D^{\omega } h_{\alpha }^{\omega } \right)=0, \nonumber \end{align} \tag{ 20 }$$

where the function $h_{\alpha }^{\omega } \left(x(t)\right)$ is given by expressions (16) and (17). This is the Liouville equation generalized to the case of active particles with pair interactions under the influence of external stochastic fields depending on the velocities of the particles. With equations (13)–(16) equation (20) can be written as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_21_} \frac{\partial D^{\omega } }{\partial t} +\sum_{1\leqslant \alpha \leqslant N}\frac{\partial }{\partial {\mathbf x}_{\alpha } } \left(D^{\omega } \frac{\partial H}{\partial {\mathbf p}_{\alpha } } \right) +\sum_{1\leqslant \alpha \leqslant N}\frac{\partial }{\partial {\mathbf p}_{\alpha } } \left(D^{\omega } \left(-\frac{\partial H}{\partial {\mathbf x}_{\alpha } } -\frac{\partial R}{\partial \mathbf p_{\alpha } } \right)\right)=0. \nonumber \end{align} \tag{ 21 }$$

In what follows we will use the Liouville equation (21) transformed with equations (13)–(16), (8), and (9) to the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__1_22_} \frac{\partial D^{\omega } }{\partial t} +\sum_{1\leqslant \alpha \leqslant N}\frac{{\mathbf p}_{\alpha } }{m} \frac{\partial D^{\omega } }{\partial {\mathbf x}_{\alpha } } +\sum_{1\leqslant \alpha <\beta \leqslant N}\frac{\partial }{\partial {\mathbf p}_{\alpha } } D^{\omega } {{\mathbf F}}_{\alpha ,\beta } +\sum_{1\leqslant \alpha \leqslant N}\frac{\partial }{\partial \mathbf p_{\alpha } } D^{\omega } \mathbf Y_{\alpha }^{\omega } =0, \nonumber \end{align} \tag{ 22 }$$

where ${{\mathbf F}}_{\alpha, \beta}$ and $\mathbf Y_{\alpha}^{\omega}$ are determined by equations (7)–(9). Equation (22) is an example of the evolution equation with multiplicative noise. Now, the goal is to average this equation over realizations of the external random force $\mathbf Y_{\alpha}^{\omega}$.

We introduce the N-particle distribution function $D \left (x_{1}, ..., x_{N}; t \right)$, which is the distribution function $D ^ {\omega} \left (x_{1}, ..., x_{N}; t \right)$ (see equation (19) averaged over the random external field ${\mathbf Y} ^ {\omega} \left ( {\mathbf x}, t \right)$ with probability density $W [ {\mathbf Y} ^ {\omega}]$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_1_} D\left(x_{1} ,...,x_{N} ;t\right)\equiv \left\langle D^{\omega } \left(x_{1} ,...,x_{N} ;t\right)\right\rangle_{\omega } , \qquad \left\langle ...\right\rangle_{\omega } \equiv \int D {\mathbf Y}^{\omega } \left( {\mathbf x},t\right)W[ {\mathbf Y}^{\omega } ]... \nonumber \end{align} \tag{ 23 }$$

By using the averaging operation (23) for equation (22), we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_2_} \frac{\partial D}{\partial t} +\sum_{1\leqslant \alpha \leqslant N}\frac{{\mathbf p}_{\alpha } }{m} \frac{\partial D}{\partial {\mathbf x}_{\alpha } } +\sum_{1\leqslant \alpha <\beta \leqslant N}\frac{\partial }{\partial {\mathbf p}_{\alpha } } D {\mathbf F}_{\alpha ,\beta } +\sum_{1\leqslant \alpha \leqslant N}\frac{\partial }{\partial {\mathbf p}_{\alpha } } \left\langle D^{\omega } \mathbf Y_{\alpha }^{\omega } \right\rangle_{\omega } =0. \nonumber \end{align} \tag{ 24 }$$

To have a closed evolution equation for the introduced distribution function, it is necessary to express the value of $\left \langle D ^ {\omega} {\mathbf Y}_{\alpha} ^ {\omega} \right \rangle_{\omega}$ through $D \left (x_{1}, ..., x_{N}; t \right)$. We use the so-called Furutsu–Novikov formula [34, 35], which was proved for the case of Gaussian distributions of the external random field. For nonGaussian random fields, the Furutsu–Novikov formula is generalized in [29] (see also [30]). In this article, we will not recount the latter proof, and refer the reader to the works cited above. We use the result of the proof in [29] in the case of a Gaussian distribution of multiplicative noise. Thus, we get

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_3_} &\nonumber \left\langle Y_{i}^{\omega } \left(x_{\alpha } ,t\right)D^{\omega } [ {\mathbf Y}^{\omega } ]\right\rangle_{\omega } =Y_{i} \left(x_{\alpha } ,t\right)\left\langle D^{\omega } [ {\mathbf Y}^{\omega } ]\right\rangle_{\omega } \nonumber \\ &\qquad\qquad+\int {\rm d}x' \int_{-\infty }^{\infty }{\rm d}t'y_{ij} \left(x_{\alpha } ,x',t-t'\right) \left\langle \frac{\delta D^{\omega } [ {\mathbf Y}^{\omega } ]}{\delta Y_{j} \left(x',t'\right)} \right\rangle_{\omega }, \nonumber \end{align} \tag{ 25 }$$

where $\newcommand{\e}{{\rm e}} Y_{i} \left(x_{\alpha } ,t\right)\equiv \left\langle Y_{i}^{\omega } \left(x_{\alpha } ,t\right)\right\rangle_{\omega }$, $\newcommand{\e}{{\rm e}} x_{\alpha } \equiv \left\{{\mathbf x}_{\alpha } , {\mathbf p}_{\alpha } \right\}$ and $y_{ij} \left(x_{\alpha } ,x',t-t'\right)$ is a pair correlation function of the external Gaussian noise ($\newcommand{\e}{{\rm e}} x'\equiv \left\{{\mathbf x}', {\mathbf p}'\right\}$):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_4_} y_{ij} \left(x_{\alpha } ,x',t-t'\right)=\left\langle Y_{i}^{\omega } \left(x_{\alpha } ,t\right)Y_{j}^{\omega } \left(x',t'\right)\right\rangle_{\omega } -\left\langle Y_{i}^{\omega } \left(x_{\alpha } ,t\right)\right\rangle_{\omega } \left\langle Y_{j}^{\omega } \left(x',t'\right)\right\rangle_{\omega } . \nonumber \end{align} \tag{ 26 }$$

In what follows we use $\newcommand{\e}{{\rm e}} Y_{i} \left (x \right) \equiv 0$. Now, let us consider

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_5_} I_{i} \equiv \int {\rm d}x'\int_{-\infty }^{\infty }{\rm d}t' y_{ij} \left(x_{\alpha } ,x',t-t'\right)\left\langle \frac{\delta D^{\omega } [ {\mathbf Y}^{\omega } ]}{\delta Y_{j}^{\omega } \left(x',t'\right)} \right\rangle_{\omega } \nonumber \end{align} \tag{ 27 }$$

in more detail. We assume that the pair correlation function $y_{ij} \left (x_{\alpha}, x ', t-t' \right)$ is different from zero in the interval $\left | t-t '\right | \leqslant \tau_{0}$. We also assume that when $t \sim t '$, the pair correlation function $y_{ij} \left (x_{\alpha}, x', t-t '\right)$ has a sharp maximum. Then the functional derivative $\frac {\delta D ^ {\omega} [ {\mathbf Y} ^ {\omega}]} {\delta Y_{j} ^ {\omega} \left (x ', t' \right)}$ is to be evaluated only at $t \approx t '$. Moreover, as shown in [29, 34–36], an exact expression for this derivative can be obtained only when $t \approx t '$.

In fact, the variational derivative $\frac {\delta D ^ {\omega} [ {\mathbf Y} ^ {\omega}]} {\delta Y_{j} ^ {\omega} \left (x ', t '\right)}$ at $t \approx t'$ undergoes a jump:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_6_} \frac{\delta D^{\omega } [ {\mathbf Y}^{\omega } ]}{\delta \delta Y_{j}^{\omega } \left(x',t'\right)} \ne 0,\quad t'\leqslant t, \quad \frac{\delta D^{\omega } [ {\mathbf Y}^{\omega } ]}{\delta Y_{j}^{\omega } \left(x',t'\right)} =0,\quad t'>t. \nonumber \end{align} \tag{ 28 }$$

The latter circumstance is due to the fact that according to the equation (22), the value of $D ^ {\omega} \left (t \right)$ cannot depend on the field $Y_{j} ^ {\omega} \left ( {\mathbf x} ', t '\right)$ taken at a time later than t. According to equation (27) the integration over $t '$ in the formula (27) is held in the range of $- \infty$ to t, instead of $- \infty$ to $+ \infty$.

Differentiating equation (22) by $Y {j} ^ {\omega} \left ( {\mathbf x} ', t' \right)$ and noting that according to equation (27) the derivative $\frac {\partial} {\partial t} \frac {\delta D ^ {\omega} [ {\mathbf Y} ^ {\omega}]} {\delta Y_{j} ^ {\omega} \left (x ', t' \right)}$ must have a $\delta$-like shape in time (while the value $\frac {\delta D ^ {\omega} [ {\mathbf Y} ^ {\omega}]} {\delta Y_{j} ^ {\omega} \left (x ', t' \right)}$ does not), the following expression for the functional derivative is obtained (see [29]):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_7_} &\nonumber \frac{\delta D^{\omega } [ {\mathbf Y}^{\omega } ]}{\delta Y_{j}^{\omega } \left(x',t'\right)} \approx -\vartheta \left(t-t'\right)\sum_{1\leqslant \beta \leqslant N}\delta \left(x'-x_{\beta } \right)\frac{\partial D^{\omega } [ {\mathbf Y}^{\omega } ]}{\partial p_{\beta j} }, \nonumber \\ &\delta \left(x'-x_{\beta } \right)\equiv \delta \left( {\mathbf x}'- {\mathbf x}_{\beta } \right)\delta \left( {\mathbf p}'- {\mathbf p}_{\beta } \right), \nonumber \end{align} \tag{ 29 }$$

where $\vartheta \left(t-t'\right)$ is the Heaviside function. This formula allows us to represent I_i, equation (27), in the following form (see equations (23) and (24):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_8_} I_{i} =\int_{-\infty }^{t}{\rm d}t' \sum_{1\leqslant \beta \leqslant N}y_{ij} \left(x,x_{\beta } ;t-t'\right)\frac{\partial D}{\partial p_{\beta j} }. \nonumber \end{align} \tag{ 30 }$$

Thus, the averaged Liouville equation, generalized to the case of systems of many particles with active interaction reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_9_} \nonumber \frac{\partial D}{\partial t} +\sum_{1\leqslant \alpha \leqslant N}\frac{{\mathbf p}_{\alpha } }{m} \frac{\partial D}{\partial {\mathbf x}_{\alpha } } +\sum_{1\leqslant \alpha <\beta \leqslant N}\frac{\partial D {\mathbf F}_{\alpha ,\beta } }{\partial {\mathbf p}_{\alpha } } +\sum_{1\leqslant \alpha \leqslant N}\frac{\partial D {\mathbf Y}\left(x_{\alpha } ,t\right)}{\partial {\mathbf p}_{\alpha } } \nonumber \\ - \int_{-\infty }^{t}{\rm d}t' \sum_{1\leqslant \alpha ,\beta \leqslant N}\frac{\partial }{\partial p_{\alpha i} } y_{ij} \left(x_{\alpha } ,x_{\beta } ;t-t'\right)\frac{\partial D}{\partial p_{\beta j} } =0. \nonumber \end{align} \tag{ 31 }$$

Taking into account that the pair correlation function $y_{ij} \left (x_{\alpha}, x_{\beta}; t-t '\right)$ has a sharp maximum at $t \approx t'$, and also assuming that this function is an even function of difference $t-t'$,

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle y_{ij} \left (x_{\alpha}, x_{\beta}; t-t '\right) = y_{ij} \left (x_{\alpha}, x_{\beta}; t'-t \right), \nonumber \end{align*}$$

then equation (31) gets an even simpler form:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_11_} \nonumber \frac{\partial D}{\partial t} +\sum_{1\leqslant \alpha \leqslant N}\frac{\mathbf {\mathbf p}_{\alpha } }{m} \frac{\partial D}{\partial {\mathbf x}_{\alpha } } +\sum_{1\leqslant \alpha <\beta \leqslant N}\frac{\partial D {\mathbf F}_{\alpha ,\beta } }{\partial \mathbf p_{\alpha } } \nonumber \\ - \frac{1}{2} \sum_{1\leqslant \alpha ,\beta \leqslant N}\frac{\partial }{\partial p_{\alpha i} } y_{ij} \left(x_{\alpha } ,x_{\beta } \right)\frac{\partial D}{\partial p_{\beta j} } =0. \nonumber \end{align} \tag{ 32 }$$

Here, we introduce the following notation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_12_} y_{ij} \left(x_{\alpha } ,x_{\beta } \right)\equiv \int_{-\infty }^{\infty }{\rm d}\tau y_{ij} \left(x_{\alpha } ,x_{\beta } ;\tau \right). \nonumber \end{align} \tag{ 33 }$$

Equation (33) can be put in another form suitable for further calculations:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_13_} &\nonumber \frac{\partial D}{\partial t} +\sum_{1\leqslant \alpha \leqslant N}\frac{{\mathbf p}_{\alpha } }{m} \frac{\partial D}{\partial \mathbf x_{\alpha } } +\sum_{1\leqslant \alpha <\beta \leqslant N}\frac{\partial D {\mathbf F}_{\alpha ,\beta } }{\partial {\mathbf p}_{\alpha } } \hfill \nonumber \\ &\qquad-\frac{1}{2} \sum_{1\leqslant \alpha \leqslant N}\frac{\partial }{\partial p_{\alpha i} } y_{ij} \left(x_{\alpha } ,x_{\alpha } \right)\frac{\partial D}{\partial p_{\alpha j} } - \sum_{1\leqslant \alpha <\beta \leqslant N}\frac{\partial }{\partial p_{\alpha i} } y_{ij} \left(x_{\alpha } ,x_{\beta } \right)\frac{\partial D}{\partial p_{\beta j} } =0. \nonumber \end{align} \tag{ 34 }$$

Note that, in fact, the developed technique allows us to obtain a generalized Liouville equation also in the case of a nonGaussian random field whenever these distributions have moments of any order (see [29]). In the present paper, however, we restrict ourselves to a Gaussian external random field.

For further calculations, we specify the explicit form of the pair correlation function $y_{ij} \left (x_{\alpha}, x_{\beta} \right)$. Using equations (9) and (26) we arrive at the following expression for $y_{ij} \left (x_{\alpha}, x_{\beta} \right)$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_14_} &\nonumber y_{ij} \left(x_{\alpha } ,x_{\beta } \right)=e_{i\alpha }^{\rm h} e_{j\beta }^{\rm h} g\left( {\mathbf x}_{\alpha } , {\mathbf x}_{\beta } \right)+\left(\delta_{i \ell} -e_{i\alpha }^{\rm h} e_{\ell\alpha }^{\rm h} \right)\left(\delta_{j\ell} -e_{j\beta }^{\rm h} e_{\ell\beta }^{\rm h} \right)h\left( {\mathbf x}_{\alpha } , {\mathbf x}_{\beta } \right), \nonumber \\ &e_{\alpha i}^{\rm h} =\frac{p_{i\alpha } }{\left| {\mathbf p}_{\alpha } \right|}. \nonumber \end{align} \tag{ 35 }$$

Here, we introduce the notations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_15_} & \nonumber g\left( {\mathbf x}_{\alpha } , {\mathbf x}_{\beta } \right)\equiv \int_{-\infty }^{\infty }{\rm d}t \left\langle X^{\omega } \left( {\mathbf x}_{\alpha } ,t\right) X^{\omega } \left( {\mathbf x}_{\beta } ,t'\right)\right\rangle_{\omega }, \nonumber \\ &\qquad\qquad\delta_{\ell k} { h} \left( {\mathbf x}_{\alpha } , {\mathbf x}_{\beta } \right)\equiv \int_{-\infty }^{\infty }{\rm d}t \left\langle X_{\ell}^{\omega } \left( {\mathbf x}_{\alpha } ,t\right)X_{k}^{\omega } \left( {\mathbf x}_{\beta } ,t'\right)\right\rangle_{\omega }. \nonumber \end{align} \tag{ 36 }$$

When obtaining expressions (35) we assumed that the stochastic force $Y_{i} ^ {\omega} \left (x, t \right)$ has the following properties:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__2_16_} \left\langle X^{\omega } \left( {\mathbf x}_{\alpha } ,t\right)X_{i}^{\omega } \left( {\mathbf x}_{\beta } ,t'\right)\right\rangle =0, \qquad \left\langle X^{\omega } \left( {\mathbf x}_{\alpha } ,t\right)\right\rangle =0, \qquad \left\langle X_{i}^{\omega } \left( {\mathbf x}_{\beta } ,t'\right)\right\rangle =0 . \nonumber \end{align} \tag{ 37 }$$

The last two formulas in equation (37) are the result of the requirement $\newcommand{\e}{{\rm e}} Y_{i} \left (x, t \right) \equiv \left \langle Y_{i} ^ {\omega} \left (x, t \right) \right \rangle_{\omega} = 0$ (see equation (9)).

Along with the probability density $D(x_1,...,x_N,t)$ we can introduce the probability of finding one or more particles in the given elements of phase space, regardless of the positions of the remaining particles (see also [22, 23]). These probabilities can be obtained by integrating the function D over all variables except those that relate to the particles under consideration:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__3_1_} f_{S} \left(x_{1} ,...,x_{S} ;t\right)={{\mathcal V}}^{S} \int {\rm d}x_{S+1} ...\int {\rm d}x_{N} D\left(x_{1} ,...,x_{N} ;t\right), \quad x_{\alpha } \equiv \left( {\mathbf x}_{\alpha } , {\mathbf p}_{\alpha } \right), \nonumber \end{align} \tag{ 38 }$$

where $D \left (x_{1}, ..., x_{N}; t \right)$ satisfies equation (34) and ${{\mathcal V}}$ is the system volume. Following the procedure described in [37, 38], after some transformations we arrive at the following equation for the S-particle distribution function $f_{S} \left (x_{1}, ..., x_{S}; t \right)$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__3_2_} \nonumber \frac{\partial f_{S} }{\partial t} +&\sum_{1\leqslant \alpha \leqslant S}\frac{{\mathbf p}_{\alpha } }{m} \frac{\partial f_{S} }{\partial x_{\alpha } } -\frac{1}{2} \sum_{1\leqslant \alpha \leqslant S}\frac{\partial }{\partial p_{\alpha i} } y_{ij} \left(x_{\alpha } ,x_{\alpha } \right)\frac{\partial f_{S} }{\partial p_{\alpha j} } \qquad \qquad \qquad \nonumber \\ \nonumber & +\sum_{1\leqslant \alpha <\beta \leqslant S}\frac{\partial f_{S} {\mathbf F}_{\alpha ,\beta } }{\partial \mathbf p_{\alpha } } -\sum_{1\leqslant \alpha <\beta \leqslant S}\frac{\partial }{\partial p_{\alpha i} } y_{ij} \left(x_{\alpha } ,x_{\beta } \right)\frac{\partial f_{S} }{\partial p_{\beta j} } \nonumber \\ &\nonumber =-\frac{1}{\mathfrak{v}} \sum_{1\leqslant \alpha \leqslant S}\frac{\partial }{\partial \mathbf p_{\alpha } } \int {\rm d}x_{S+1} f_{S+1} {{\mathbf F}}_{\alpha ,S+1} \nonumber \\ &+\frac{1}{\mathfrak{v}} \sum_{1\leqslant \alpha \leqslant S}\frac{\partial }{\partial p_{\alpha i} } \int {\rm d}x_{S+1} y_{ij} \left(x_{\alpha } ,x_{S+1} \right)\frac{\partial f_{S+1} }{\partial p_{S+1j} }, \qquad {\mathfrak{v}} \equiv \frac{{{\mathcal V}}}{N}, \nonumber \end{align} \tag{ 39 }$$

where the quantities ${\mathbf F}_{\alpha, \beta}$ and $y_{ij} \left (x_{\alpha}, x_{\beta} \right)$ are still given by equations (8), (26), (35) and (36). It is easy to see that the equation for the S-particle distribution function includes an S  +  1-particle distribution function. Thus, we obtain an infinite chain of kinetic equation (39). These chains are a generalization of the well-known chain of Bogolyubov–Born–Green–Kirkwood–Yvon equations in the case of identical active interacting particles under the influence of external stochastic fields. It is necessary to make the following remark. According to the definition (38), the distribution functions of a higher order contain all the information contained in the functions of lower order [22]. This leads to the fact that with the increase in the order S, the distribution functions $f_{S} \left (x_{1}, ..., x_{S}; t \right)$ become increasingly complex. According to equation (39), it is necessary to consider the distribution functions up to S  =  N, and we conclude that the resulting chain of equation (39) is equivalent to the Liouville equation (34). In other words, the most complete description of the studied systems is equally complex both within the framework of the full the distribution function $D \left (x_{1}, ..., x_{N}; t \right)$ and within the one of the many-particle distribution function $f_{S} \left (x_{1}, ..., x_{S}; t \right)$.

A significant simplification in description of the state of the system occurs in two cases: when the interaction between the particles is small, and when the number density of particles is small and the interaction is arbitrary (but is such that it does not lead to the formation of bound states) [38]. This simplification in the description is the consequence of the difference in the evolutionary behavior of many- and single-particle distribution functions. In fact, at an early stage of evolution, when the time t is small compared to the characteristic time of chaotization $\tau_{0}$, the multiparticle distribution functions $f_{S} \left (x_{1}, .. ., x_{S}; t \right)$ change rapidly over time, in contrast to the single-particle distribution function $f_{1} \left (x, t \right)$. The single-particle distribution function experiences significant changes in time at times much longer than the relaxation time of the system $\tau_{{\rm r}}$, and $\tau_{{\rm r}} \gg \tau_{0}$. The order of magnitude of $\tau_{0}$ is determined by the duration of one collision. It can also be estimated as a value of the order of $r_{{\rm c}} / {\bar {v}}$, where $\bar {v}$ is the average velocity of the particles and $r_{{\rm c}}$ is the correlation radius, that is, the distance from which the multiparticle distribution functions break up into a single-particle product. Usually, the order of magnitude of the correlation radius coincides with the radius of action of the forces. This is why $\tau_{0}$ can be considered of the order of the duration of one collision (see also [22, 23] in this regard). Time $\tau_{{\rm r}}$ in order of magnitude should be the same as the time of the establishment of the statistical equilibrium state in the system (for more details see [23]). Such a difference in the evolutionary behavior of the single-particle and many-particle distribution functions formed the basis of Bogolyubov's ideas about a hierarchy of the system relaxation times [22].

According to the idea of a hierarchy of relaxation times, the evolution of the many-particle system can be divided into several stages. Each subsequent stage of evolution differs from the previous by a simplification in the description of the evolution of a system of many particles. The simplest scenario for the evolution of systems of many particles is as follows. When $\tau_{0} \ll t \ll \tau_{{\rm r}}$ there takes place a kinetic stage of evolution of the system where the system behavior can be described by a single-particle distribution function. This description of the system evolution is much easier than that using the multiparticle distribution functions. Further simplification of the description of many-particle systems occurs when $t \gg \tau_{{\rm r}}$ (the hydrodynamic stage of evolution of the system) and when the behavior of the system can be described by the hydrodynamic description parameters such as the particle number density, the average velocity, and the temperature of the medium. Such a gradual simplification of the system description of the reduced description method are based on [22, 23].

In this paper, the method of reduced description of nonequilibrium processes will be used for the derivation of the kinetic equations describing the evolution of systems of interacting active particles in an external random field. The initial equations will be the chain equation (39). The mathematical formulation of the idea of a hierarchy of relaxation times of the system is a time-functional dependence of many-particle distribution functions $f_{S} \left (x_{1}, ..., x_{S}; t \right)$ only through a dependence on time of the parameters of the reduced description at the appropriate stage of evolution. In particular, at the kinetic stage of the evolution, the many-particle distribution functions depend on time only through the one-particle distribution function $f_{1} (x ', t)$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__3_3_} f_{S} \left(x_{1} ,...,x_{S} ;t\right)=f_{S} \left(x_{1} ,...,x_{S} ;f_{1} (x',t)\right). \nonumber \end{align} \tag{ 40 }$$

In addition to the functional hypothesis (40), the reduced description method is also based on the principle of spatial correlation weakening. In the language of multiparticle distribution functions, this principle can be summarized as follows [23]. Let S of the particles be divided into two subgroups of particles containing $S'$ and $S''$ particles, respectively, and $S = S' + S''$. If the distance R between these subgroups of particles increases infinitely, $R \to \infty$, then due to the weakening of correlations between particles the S-particle distribution function decomposes into the product of the distribution functions related to each particle subgroup:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__3_4_} f_{S} \left(x_{1} ,...,x_{S} ;t\right)\mathop{\to }\limits_{R\to \infty } f_{S'} \left(x'_{1} ,...,x'_{S} ;t\right)\,f_{S''} \left(x''_{1} ,...,x''_{S} ;t\right). \nonumber \end{align} \tag{ 41 }$$

In equation (41), the 'prime' sign is used to indicate the coordinates and momenta of the particles of the subgroup $S '$, and the 'double prime' sign is used to indicate the coordinates and momenta of the particles of the second subgroup. It should be noted, however, that the principle of spatial correlation weakening equation (41) refers to the many-particle distribution functions for which the thermodynamic limit is made: $N \to \infty$, $\mathcal{V} \to \infty$, and $\left (N / \mathcal{V} \right) = const$ [38].

According to equation (40), the time derivative of ${\partial f_{S}} / {\partial t}$ in equation (39) when $S \ne 1$ must be understood as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__3_5_} \frac{\partial }{\partial t} f_{S} \left(x_{1} ,...,x_{S} ;f_{1} (x,t)\right)=\int {\rm d}x'\frac{\delta f_{S} \left(x_{1} ,...,x_{S} ;f_{1} (x,t)\right)}{\delta f_{1} (x',t)} \frac{\partial f_{1} (x',t)}{\partial t}, \nonumber \end{align} \tag{ 42 }$$

where ${\delta f_{S} \left (\,f_{1} (x, t) \right)} / {\delta f_{1} (x ', t)}$ is the functional derivative. According to equation (39), the single-particle distribution function itself $f_{1} (x ', t)$ must satisfy the following equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__3_6_} \frac{\partial f_{1} }{\partial t} +\frac{{\mathbf p}_{1} }{m} \frac{\partial f_{1} }{\partial {\mathbf x}_{1} } -\frac{1}{2} \frac{\partial }{\partial p_{1i} } y_{ij} \left(x_{1} ,x_{1} \right)\frac{\partial f_{1} }{\partial p_{1j} } =\frac{1}{\mathfrak{v}} L\left(x_{1} ;f_{1} \right), \nonumber \end{align} \tag{ 43 }$$

where as before $\mathfrak{v} = {\mathcal{V}} / {N}$, and $L \left (x_{1}; f_{1} \right)$ is the generalized collision integral defined by the formula

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__3_7_} \nonumber L\left(x_{1} ;f_{1} \right)\equiv -\frac{\partial }{\partial {\mathbf p}_{1} } \int {\rm d}x_{2} f_{2} \left(x_{1} ,x_{2} ;f_{1} \right) {\mathbf F}_{1,2} \qquad \qquad \qquad \nonumber \\ \hfill +\frac{\partial }{\partial p_{\alpha i} } \int {\rm d}x_{2} y_{ij} \left(x_{1} ,x_{2} \right)\frac{\partial }{\partial p_{2j} } f_{2} \left(x_{1} ,x_{2} ;f_{1} \right). \nonumber \end{align} \tag{ 44 }$$

To make the equation (43) closed, one must obtain the collision integral (44) as a functional of the particle distribution function, and, therefore, it is necessary to cut off an infinite chain of equation (39). Clearly, this can only be done within a certain approximation. In particular, in the system of nonactive particles, such a cutoff may be implemented in the two cases mentioned above, i.e. when the interaction between the particles is small, and when the particle density is low and the interaction is arbitrary but does not lead to the formation of bound states [22]. Similar situations can be implemented in the case of a system of identical active particles with interaction; this is also discussed in this paper. We will demonstrate this in the case of weak interaction of all kinds between the active particles and the external noise of low intensity. In other words, we assume that the forces ${{\mathbf F}}_{\alpha, \beta}$ and the correlation functions of an external random field are small.

First, however, we make some remarks. The functional relation (40) does not necessarily imply an expansion of $f_{S} \left (x_{1}, ..., x_{S}; f_{1} (x ', t) \right)$ in the functional perturbation series by the one-particle distribution function. This expansion must be realized only in one of the above-mentioned cases of chain breaking, namely when the particle density is low. We recall that thus arises the famous question of the possible divergences in higher orders of perturbation theory by a small particle number density and about a renormalization of this theory (see, e.g. [39–41]). In the case of perturbation theory by the weak interaction between the particles, these issues do not appear, as is easily seen from the subsequent calculations (see also [23, 29]).

Here we follow the methodology suggested in [23]. Using equations (42) and (43), a chain of equation (39) can be written as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__4_1_} -\int {\rm d}x\frac{\delta f_{S} \left(\,f_{1} \right)}{\delta f_{1} (x,t)} \frac{{\mathbf p}}{m} \frac{\partial f_{1} (x,t)}{\partial {\mathbf x}} +\sum_{1\leqslant \alpha \leqslant S}\frac{{\mathbf p}_{\alpha } }{m} \frac{\partial f_{S} \left(\,f_{1} \right)}{\partial {\mathbf x}_{\alpha } } =\frac{1}{\mathfrak{v}} K_{S} \left(\,f_{1} \right), \nonumber \end{align} \tag{ 45 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__4_2_} \nonumber K_{S} \left(\,f_{1} \right)\equiv & -\mathfrak{v}\sum_{1\leqslant \alpha <\beta \leqslant S}\frac{\partial f_{S} {\mathbf F}_{\alpha ,\beta } }{\partial {\mathbf p}_{\alpha } } +\frac{1}{2} \sum_{1\leqslant \alpha \leqslant S}\frac{\partial }{\partial p_{\alpha i} } y_{ij} \left(x_{\alpha } ,x_{\alpha } \right)\frac{\partial f_{S} }{\partial p_{\alpha j} } \nonumber \\ \nonumber &+\sum_{1\leqslant \alpha <\beta \leqslant S}\frac{\partial }{\partial p_{\alpha i} } y_{ij} \left(x_{\alpha } ,x_{\beta } \right)\frac{\partial f_{S} }{\partial p_{\beta j} } - \sum_{1\leqslant \alpha \leqslant S}\frac{\partial }{\partial {\mathbf p}_{\alpha } } \int {\rm d}x_{S+1} f_{S+1} {{\mathbf F}}_{\alpha ,S+1} \nonumber \\ \nonumber & +\sum_{1\leqslant \alpha \leqslant S}\frac{\partial }{\partial p_{\alpha i} } \int {\rm d}x_{S+1} y_{ij} \left(x_{\alpha } ,x_{S+1} \right)\frac{\partial f_{S+1} }{\partial p_{S+1j} } \nonumber \\ \hfill &-\int {\rm d}x_{1} \frac{\delta f_{S} \left(\,f_{1} \right)}{\delta f_{1} (x_{1} ,t)} \left\{L\left(x_{1} ;f_{1} \right)+\frac{1}{2} \frac{\partial }{\partial p_{1i} } y_{ij} \left(x_{1} ,x_{1} \right)\frac{\partial f_{1} (x_{1} ,t)}{\partial p_{1j} } \right\}. \nonumber \end{align} \tag{ 46 }$$

The chain of equations (45) and (46) must be supplemented by the 'initial conditions'. To this end, following [22, 23], we introduce an auxiliary parameter $\tau$; it has the dimension of time but does not necessarily represent physical time. We next consider the many-particle distribution function $f_ {S} \left ( {\mathbf x}_ {1} - \frac {{\mathbf p}_ {1}} {m} \tau, {\mathbf p}_ {1}, ..., \mathbf x_ {S} - \frac {{\mathbf p}_ {S}} {m} \tau, {\mathbf p}_ {S}; f_ {1} \right)$. According to equation (41), this function must satisfy the asymptotic relation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__4_3_} f_{S} \left( {\mathbf x}_{1} -\frac{{\mathbf p}_{1} }{m} \tau , {\mathbf p}_{1} ,..., {\mathbf x}_{S} -\frac{{\mathbf p}_{S} }{m} \tau , {\mathbf p}_{S} ;f_{1} \right)\mathop{\to }\limits_{\tau \to \infty } \prod_{1\leqslant \alpha \leqslant S}f_{1} \left( {\mathbf x}_{\alpha } -\frac{{\mathbf p}_{\alpha } }{m} \tau , {\mathbf p}_{\alpha } \right). \nonumber \end{align} \tag{ 47 }$$

If we define further the shift operator $\hat {\Lambda} _ {S} ^ {0}$ in the coordinate space with the formula

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__4_4_} {\rm i}\hat{\Lambda }_{S}^{0} \equiv \sum_{1\leqslant \alpha \leqslant S}\frac{{\mathbf p}_{\alpha } }{m} \frac{\partial }{\partial {\mathbf x}_{\alpha } }, \nonumber \end{align} \tag{ 48 }$$

the condition (47) may be rewritten as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__4_5_} {\rm e}^{{\rm i}\tau \hat{\Lambda }_{S}^{0} } f_{S} \left(\tau \right)\mathop{\to }\limits_{\tau \to \infty } \prod_{1\leqslant \alpha \leqslant S}f_{1} \left(x_{\alpha } \right), \nonumber \end{align} \tag{ 49 }$$

where $\newcommand{\e}{{\rm e}} \exp \left({\rm i}\tau \hat{\Lambda }_{S}^{0} \right)$ is a so-called 'free evolution operator' and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__4_6_} f_{S} \left(\tau \right)\equiv f_{S} \left(x_{1} ,...,x_{S} ;{\rm e}^{-{\rm i}\tau \hat{\Lambda }_{1}^{0} } f_{1} \left(x'\right)\right) =f_{S} \left(x_{1} ,...,x_{S} ;f_{1} \left( {\mathbf x}'-\frac{{\mathbf p}'}{m} \tau , {\mathbf p}'\right)\right). \nonumber \end{align} \tag{ 50 }$$

Now, equation (45) can be written in the following way:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__4_7_} \frac{\partial }{\partial \tau } {\rm e}^{{\rm i}\tau \hat{\Lambda }_{S}^{0} } f_{S} \left(\tau \right)=\frac{1}{\mathfrak{v}} {\rm e}^{{\rm i}\tau \hat{\Lambda }_{S}^{0} } K_{S} \left(\tau \right), \nonumber \end{align} \tag{ 51 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__4_8_} K_{S} \left(\tau \right) \equiv K_{S} \left(x_{1} ,...,x_{S} ;{\rm e}^{-{\rm i}\tau \hat{\Lambda }_{1}^{0} } f_{1} \left(x'\right)\right) =K_{S} \left(x_{1} ,...,x_{S} ;f_{1} \left( {\mathbf x}'-\frac{{\mathbf p}'}{m} \tau , {\mathbf p}'\right)\right). \nonumber \end{align} \tag{ 52 }$$

By integrating equation (51) over $\tau$ within the limits from $-\infty$ to 0 and using the asymptotic conditions (49), we get

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__4_9_} f_{S} \left(x_{1} ,...,x_{S} ;f_{1} \left(x'\right)\right)=\prod_{1\leqslant \alpha \leqslant S}f_{1} \left(x_{\alpha } \right) +\frac{1}{\mathfrak{v}} \int_{-\infty }^{0}{\rm d}\tau {\rm e}^{{\rm i}\tau \hat{\Lambda }_{S}^{0} } K_{S} \left(\tau \right). \nonumber \end{align} \tag{ 53 }$$

The ratio (53) allows us to develop a perturbation theory over the two small parameters: the weak interaction and the intensity of stochastic effects. Under such assumptions, the value $K_ {S} \left (\tau \right)$ (see equation (47) can be considered small; therefore, in the main approximation we have

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle f_{S} \left(x_{1} ,...,x_{S} ;f_{1} \left(x'\right)\right)=\prod_{1\leqslant \alpha \leqslant S}f_{1} \left(x_{\alpha } \right) ,\nonumber \end{align*}$$

which implies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__4_10_} f_{2} \left(x_{1} ,x_{2} \right)=f_{1} \left(x_{1} \right)\,f_{1} \left(x_{2} \right). \nonumber \end{align} \tag{ 54 }$$

By substituting equation (54) into equation (44) and using equations (7) and (8), we obtain the following closed kinetic equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{4.11a} \nonumber \frac{\partial f_{1} \left(x_{1} ,t\right)}{\partial t} +\frac{{\mathbf p}_{1} }{m} \frac{\partial f_{1} \left(x_{1} ,t\right)}{\partial \mathbf x_{1} } -\frac{1}{2} \frac{\partial }{\partial p_{1i} } y_{ij} \left(x_{1} ,x_{1} \right)\frac{\partial f_{1} \left(x_{1} ,t\right)}{\partial p_{1j} } \nonumber \\ \nonumber = \frac{1}{\mathfrak{v}} \frac{\partial }{\partial {\mathbf p}_{1} } f_{1} \left(x_{1} \right)\int {\rm d}x_{2} f_{1} \left(x_{2} \right)\left(\frac{\partial V_{1,2} }{\partial {\mathbf x}_{1} } +\frac{\partial R_{1,2}}{\partial {\mathbf p}_{1} } \right)\nonumber \\ +\frac{1}{\mathfrak{v}} \frac{\partial }{\partial p_{1i} } f_{1} \left(x_{1} \right)\int {\rm d}x_{2} y_{ij} \left(x_{1} ,x_{2} \right)\frac{\partial f_{1} \left(x_{2} \right)}{\partial p_{2j} }, \nonumber \end{align} \tag{ 55 }$$

where the values $V_ {1,2}$ and $R_ {1,2} ^ {}$ are given by equations (2)–(4), and the correlation function $y_ {ij} \left (x_ {1 }, x_ {2} \right)$ is still given by equation (35). Equation (55) and can be rewritten in a slightly different form:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{4.11b} &\nonumber \frac{\partial f_{1} \left(x_{1} ,t\right)}{\partial t} +\frac{{\mathbf p}_{1} }{m} \frac{\partial f_{1} \left(x_{1} ,t\right)}{\partial {\mathbf x}_{1} } -\frac{\partial U\left( {\mathbf x}_{1} ,t\right)}{\partial {\mathbf x}_{1} } \frac{\partial f_{1} \left(x_{1} ,t\right)}{\partial {\mathbf p}_{1} } \nonumber \\ &\qquad\qquad\nonumber = \frac{1}{2} \frac{\partial }{\partial p_{1i} } y_{ij} \left(x_{1} ,x_{1} \right)\frac{\partial f_{1} \left(x_{1} ,t\right)}{\partial p_{1j} } \nonumber \\ &\qquad\qquad\hfill +\frac{1}{\mathfrak{v}} \frac{\partial }{\partial p_{1i} } f_{1} \left(x_{1} ,t\right)\int {\rm d}x_{2} f_{1} \left(x_{2} ,t\right)\left[\frac{\partial R_{1,2} }{\partial p_{1i} } -\frac{\partial y_{ij} \left(x_{1} ,x_{2} \right)}{\partial p_{2j} } \right], \nonumber \end{align} \tag{ 56 }$$

or

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{4.11c} &\nonumber \frac{\partial f_{1} \left(x_{1} ,t\right)}{\partial t} +\frac{{\mathbf p}_{1} }{m} \frac{\partial f_{1} \left(x_{1} ,t\right)}{\partial \mathbf x_{1} } -\frac{\partial U\left( {\mathbf x}_{1} ,t\right)}{\partial \mathbf x_{1} } \frac{\partial f_{1} \left(x_{1} ,t\right)}{\partial {\mathbf p}_{1} } \nonumber \\ &\qquad\qquad\nonumber = \frac{1}{2} \frac{\partial }{\partial p_{1i} } y_{ij} \left(x_{1} ,x_{1} \right)\frac{\partial f_{1} \left(x_{1} ,t\right)}{\partial p_{1j} } \nonumber \\ &\qquad\qquad\hfill +\frac{1}{\mathfrak{v}} \frac{\partial }{\partial p_{1i} } f_{1} \left(x_{1} ,t\right)\int {\rm d}x_{2} \left[R_{1,2} \delta_{i,j} +y_{ij} \left(x_{1} ,x_{2} \right)\right]\frac{\partial f_{1} \left(x_{2} ,t\right)}{\partial p_{2j} } . \nonumber \end{align} \tag{ 57 }$$

We may consider an average field $U\left( {\mathbf x}_{1} ,t\right)$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__4_12_} U\left( {\mathbf x}_{1} ,t\right)=\frac{1}{\mathfrak{v}} \int {\rm d} {\mathbf x}_{2} V\left( {\mathbf x}_{1} - {\mathbf x}_{2} \right)\int {\rm d} {\mathbf p}_{2} f_{1} \left(x_{2} ,t\right),\quad f_{1} \left(x_{2} \right)\equiv f_{1} \left( {\mathbf x}_{2} , {\mathbf p}_{2} \right) , \nonumber \end{align} \tag{ 58 }$$

where $V\left( \mathbf x_{1} - \mathbf x_{2} \right)$ is the interaction potential, see equation (2). Equations (55)–(57) are the kinetic equations for the active particles with pair interactions (potential and 'dissipative' ones) between the particles under the influence of active space-dependent fluctuations. We emphasize that all the equations (55)–(57) are obtained without using the explicit form of the potential interaction $\newcommand{\e}{{\rm e}} V_ {1,2} \equiv V \left ( {\mathbf x}_ {1} - {\mathbf x}_ {2} \right)$, dissipation function R12, and the correlation function $y_ {ij} \left (x_ {1}, x_ {2} \right)$.

Note that the random force (9), typical for active fluctuations and having a local effect on the particles, as it is seen from equations (56) and (57), is determined by the pair correlation function $y_ {ij} \left (x_ {1}, x_ {2} \right)$.

We also note the following circumstance. The basis for obtaining equations (55)–(58) is formula (54), which reflects the main approximation of the perturbation theory developed in the paper for small interactions between particles, including dissipative interactions, as well as for a low external noise intensity. Naturally, formula (54), indicating that there are no correlations between particles, is already understood at the intuitive level. After all, it is clear that, in this form, the two-particle distribution function allows one to truncate the BBGKY chain and obtain a closed kinetic equation for the one-particle distribution function. However, we cannot immediately assume that the correlations between the particles can be neglected due to the small interactions and the weakness of the stochastic effects and thereby completely ignore the material of section 5 up to equations (55)–(58). First, it is not obvious. One needs to know what may be neglected. The purpose of the present study, as mentioned above, is the development of a microscopic approach based on first principles of statistical physics to construct a kinetic theory of dissipative media with active fluctuations, like done in [22–27] but with a preference for strict approaches and formulations of the monograph [23]. This means that the scheme for constructing such a microscopic approach should contain a number of consecutive, strictly controlled steps, including steps with respect to truncating the BBGKY chain. Secondly, the scheme given in section 5 allows us to, if necessary, amend the right-hand side of formula (54), first-order terms for interactions, and weak stochastic processes. Thus, we can obtain expressions that determine the correlations between the particles of the medium in the first approximation of perturbation theory. It is important to note that in equations (55)–(57) quadratic terms (see [31]) and cross terms appear in the interactions. However, we considered that within the framework of this paper it is better to refrain from considering such effects. First of all, this decision is made in order to not clutter the article with material that is unnecessary from the point of view of the article's main aims. Nevertheless, it should be recognized that the discussed terms may be responsible for very interesting effects; for example, effects related to the above-mentioned effect of the appearance of additional interaction between structural units of the environment. Currently, the authors are studying the issues raised here. As it turns out, such research requires the involvement of completely nontrivial numerical methods along with analytical methods. For this reason, the analysis of the issue requires additional time.

We make a few more comments at the end of this section. We have already emphasized above the reasons why we refused in this article to generalize the developed approach to the case of nonGaussian noise. Here, in connection with this generalization, we draw attention to the fact that when comparing the materials of this work and [30], the mathematical similarity of the derivation procedures for chains of equations (analogs of BBGKY chains) for many-particle distribution functions and kinetic equations for single-particle distribution functions is noteworthy. The similarity of the noted procedures indicates the possibility of 'organic' generalization of the approach of the present work in the case of nonGaussian noise. Such a task will not be too difficult, but will greatly clutter the calculations and significantly worsen the visibility of the results. Connected to this, the external mathematical similarity of some parts of this article and [30] in the particular case of the Gaussian stochastic effect on the system is also noteworthy. In particular, we are talking about the similarity of the mathematical structure of equations (34) and (57) of this work and, accordingly, equations (40) and (63) from [30]. Such a similarity should be considered formal if we take into account the fact that this article and [30] deal with completely different systems. In fact, equations (34) and (57) of this article describe dissipative systems with nonlinear (in general) friction and active fluctuations. This is reflected in the dependence of the pair correlation functions included in (34) and (57), not only from the coordinates, but also from the momenta (or velocities), which, as noted above, is characteristic of such systems. Equations (40) and (63) from [30] describe the evolution of dissipative systems with nonlinear (in general) friction and a source of stochastic noise, depending only on coordinates. A similar circumstance characterizes the system studied in [30] as a dissipative system with passive fluctuations.

Thus, we have developed a microscopic approach to the description of the kinetics of the active media that are allowed to come to common kinetic equations (55)–(57) in the approximation of the average or mean field; see (7), (8), and (58). The general nature of the obtained kinetic equations must be confirmed by demonstrating the facts that in some cases these equations lead to results obtained by other authors in other approaches or for other models of active media.

Here we demonstrate that the kinetic equations (55)–(57) involve some known special cases for systems of active particles. For this, we selected examples that would give results known from earlier works by other authors, using primarily phenomenological approaches to solving such problems. In this case, it is necessary to introduce a number of successive simplifications into the obtained general equations (55)–(57). In particular, it is necessary to proceed to the consideration of spatially homogeneous cases. In addition, we recall that we allowed a further simplification of the nature of noise. Namely, from a more general form of Gaussian noise (26) to time-delta-correlated noise (see (33), and (35)–(37)). Thus, there was a shift from equation (31) (a more complex form) to the simpler Liouville equation (34). For this reason, the path to equation (34) may seem somewhat long. After all, it would be possible to come to equation (34) more directly by immediately assuming the delta-correlation of noise, but then we would move away from controllability of the considered approximations and would not have a more general equation (31). Therefore, let us start by considering a spatially homogeneous state. Then, a single-particle distribution function $f_ {1} \left ( {\mathbf x}, {\mathbf p}, t \right)$ does not depend on the coordinates:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_1_} f_{1} \left( {\mathbf x}, {\mathbf p},t\right)\equiv f_{1} \left( {\mathbf p},t\right). \nonumber \end{align} \tag{ 59 }$$

We should note that the spatially correlated noise impacting the system (see equation (9)) does not necessarily prevent the existence of the spatially homogeneous states (equation (59)). The latter are possible only in the case of a zero mean of the external random force acting on the system, and the pair correlation functions $y_{ij} \left(x_{1}, x_{2} \right)$ (see equations (35)–(37)) should depend only on the difference of the coordinates of the two particles:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_2_} y_{ij} \left(x_{1} ,x_{2} \right)=g\left( {\mathbf x}_{1} - {\mathbf x}_{2} \right)e_{1i}^{\rm h} e_{2j}^{\rm h} +h\left( {\mathbf x}_{1} - {\mathbf x}_{2} \right)\left(\delta_{i \ell} -e_{1i}^{\rm h} e_{1l}^{\rm h} \right)\left(\delta_{j \ell} -e_{2j}^{\rm h} e_{2l}^{\rm h} \right). \nonumber \end{align} \tag{ 60 }$$

We recall that the assumption of a zero average of the random field was made in the present paper beginning from equation (32). We also emphasize that both specific kinds of functions $g\left( {\mathbf x}_{1} - {\mathbf x}_{2} \right)$ and $h\left( {\mathbf x}_{1} - {\mathbf x}_{2} \right)$, and their dependence on the coordinates difference, may be different due to uncorrelated noise components along and across the particle velocity (see equations (35)–(37)). As an example, this distinction is illustrated in figure 2 for the simple case when the pair correlation functions $g\left( {\mathbf x}_{1} - {\mathbf x}_{2} \right)$ and $h\left( {\mathbf x}_{1} - {\mathbf x}_{2} \right)$ decrease with the distance $\newcommand{\e}{{\rm e}} r \equiv |\mathbf{x_1} - \mathbf{x_2}|$ according to the Gaussian law but have different values of variance and amplitude. However, it is easy to see that the calculations that lead to the results of the present and subsequent sections of this work do not require the specification of the explicit form of the functions mentioned. We recall that according to equation (4), all restrictions on the general properties of functions R1,2 are contained in the expression

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_3_} R_{1,2} \equiv R\left( {\mathbf x}_{1} - {\mathbf x}_{2} , {\mathbf p}_{1} - {\mathbf p}_{2} \right), \nonumber \end{align} \tag{ 61 }$$

which follows from the Galilean invariance of the system in the absence of external influences. Moreover, since the function R1,2 is a scalar quantity, its dependence on the differences ${\mathbf x}_ {1} - {\mathbf x}_ {2}$, ${\mathbf p}_ {1} - {\mathbf p}_ {2}$ should be characterized by the following expression:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_4_} R\left( {\mathbf x}, {\mathbf p}\right)\equiv R\left( {\mathbf x}^{2} , {\mathbf p}^{2} , \mathbf{xp}\right). \nonumber \end{align} \tag{ 62 }$$

Zoom In Zoom Out Reset image size

According to equations (59)–(62), equation (55) transforms into

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_5_} & \nonumber \frac{\partial f_{1} \left( {\mathbf p}_{1} ,t\right)}{\partial t} - \frac{1}{2} \frac{\partial }{\partial p_{1i} } \left[g(0)e_{1i}^{\rm h} e_{1j}^{\rm h} +h(0)\left(\delta_{ij} -e_{1i}^{\rm h} e_{1j}^{\rm h} \right)\right]\frac{\partial f_{1} \left( {\mathbf p}_{1} ,t\right)}{\partial p_{1j} } \nonumber \\ &\qquad\qquad\nonumber =\frac{\partial }{\partial p_{1i} } f_{1} \left( {\mathbf p}_{1} ,t\right)\frac{\partial }{\partial p_{1i} } \int {\rm d} {\mathbf p}_{2} f_{1} \left( {\mathbf p}_{2} ,t\right)\bar{R}\left(\left( {\mathbf p}_{1} - {\mathbf p}_{2} \right)^{2} \right)\nonumber \\ &\quad \qquad\qquad\nonumber +\frac{\partial }{\partial p_{1i} } f_{1} \left( {\mathbf p}_{1} ,t\right)\int {\rm d} {\mathbf p}_{2} \left[\bar{g}e_{1i}^{\rm h} e_{2j}^{\rm h} +\bar{h}\left(\delta_{i \ell} -e_{1i}^{\rm h} e_{1l}^{\rm h} \right)\left(\delta_{j \ell} -e_{2j}^{\rm h} e_{2l}^{\rm h} \right)\right]\nonumber \\ &\qquad\qquad\quad \times \frac{\partial f_{1} \left( {\mathbf p}_{2} ,t\right)}{\partial p_{2j} } , \nonumber \end{align} \tag{ 63 }$$

where we introduce

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_6_} &\nonumber \bar{R}\left(\left( {\mathbf p}_{1} - {\mathbf p}_{2} \right)^{2} \right) \equiv \frac{1}{\mathfrak{v}} \int {\rm d} {\mathbf x} R\left( {\mathbf x}^{2} ,\left( {\mathbf p}_{1} - {\mathbf p}_{2} \right)^{2} , {\mathbf x}\left( {\mathbf p}_{1} - {\mathbf p}_{2} \right)\right),\nonumber \\ &\bar{g}\equiv \frac{1}{\mathfrak{v}} \int {\rm d} {\mathbf x} g\left( {\mathbf x}\right) ,\nonumber \\ &\nonumber \bar{h}\equiv \frac{1}{\mathfrak{v}} \int {\rm d} {\mathbf x} h\left( {\mathbf x}\right) . \nonumber \end{align} \tag{ 64 }$$

6.1. Brownian particles with active fluctuations. Space-independent noise

Here we study solutions of the kinetic equation (63), which are isotropic in the momentum space.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_7_} f_{1} \left( {\mathbf p},t\right)\equiv f_{1} \left(\,p,t\right). \nonumber \end{align} \tag{ 65 }$$

By taking equations (64) and (65), it is possible to reduce equation (63) to the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_8_} \nonumber \frac{\partial f_{1} \left(\,p_{1} ,t\right)}{\partial t} =&\frac{\partial }{\partial p_{1i} } e_{1i}^{\rm h} \left\{\,f_{1} \left(\,p_{1} ,t\right)\gamma \left(\,p_{1} ,t\right)\,p_{1} +\frac{1}{2} g(0)\frac{\partial f_{1} \left(\,p_{1} ,t\right)}{\partial p_{1} } \right.\nonumber \\ &+ \left.\bar{g}f_{1} \left(\,p_{1} ,t\right)\int {\rm d} {\mathbf p}_{2} \frac{\partial f_{1} \left(\,p_{2} ,t\right)}{\partial p_{2} } \right\}, \nonumber \end{align} \tag{ 66 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_9_} &\nonumber \gamma \left(\,p,t\right)\equiv 2\frac{\partial \bar{R}\left(\,p^{2} ,t\right)}{\partial p^{2} },\nonumber \\ &\bar{R}\left(\,p^{2} ,t\right)\equiv \int {\rm d} {\mathbf p}_{2} f_{1} \left(\,p_{2} ,t\right)\bar{R}\left(\left( {\mathbf p}- {\mathbf p}_{2} \right)^{2} \right). \nonumber \end{align} \tag{ 67 }$$

The resulting equation (66) is the kinetic equation for active particles with time-dependent nonlinear friction (friction factor $\gamma \left (\,p, t \right)$). This equation can be regarded as a three-dimensional generalization of the kinetic equation for quasi-Brownian particles with active fluctuations, dissipative interaction, and a space-dependent external stochastic field. This fact may be proven if we make some simplifications of equation (66).

First of all, note that the term 'Brownian particles with active fluctuations' is commonly understood as a system of particles under the influence of a space-independent stochastic field given by equations (9) and (10) [3, 21]. Consequently, to prove the above assumption, we should consider the linear friction case in equation (66) and refuse the dependence of the external noise on the coordinates. In the case of linear friction, the friction coefficient $\gamma \left (\,p \right)$ does not depend on the momentum, $\newcommand{\e}{{\rm e}} \gamma \left (\,p \right) \equiv \gamma$, and, according to equations (5) and (64), the value of $\gamma$ in this case is given by (see [30]):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_10_} \gamma =\frac{1}{\mathfrak{v}^{2} } \int {\rm d} {\mathbf x}\tilde{\gamma }\left( {\mathbf x}\right) ,\qquad \int {\rm d} {\mathbf p} f_{1} \left(\,p,t\right)=\frac{1}{\mathfrak{v}}. \nonumber \end{align} \tag{ 68 }$$

However, the consequences of the noise space-independence in equation (66) are rather hard to see immediately. For this we need to repeat the whole procedure of the kinetic equation derivation until equations (55)–(57)—assuming that the values $X^\omega(\mathbf{x}, t)$, $X^\omega_j(\mathbf{x}, t)$ in equation (9) are independent of the coordinates—and conditions (35)–(37) are fulfilled. It turns out that the result of this procedure is equivalent to equating the value $\bar g$ in equation (66) to zero. Thus, we come to the following equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_11_} &\nonumber \frac{\partial f_{1} \left(\,p_{1} ,t\right)}{\partial t} =\frac{\partial }{\partial p_{1i} } e_{1i}^{\rm h} \left\{\gamma p_{1} f_{1} \left(\,p_{1} ,t\right)+D_{\rm p} \frac{\partial f_{1} \left(\,p_{1} ,t\right)}{\partial p_{1} } \right\}, \nonumber \\ &g(0)\equiv 2D_{\rm p}=2m^2D_{\rm v}. \nonumber \end{align} \tag{ 69 }$$

The quantity $\newcommand{\e}{{\rm e}} g(0) \equiv 2D_{\rm p}$ is still defined by equations (36) and (60), keeping in mind the fact that the noise characteristic $X^\omega(\mathbf{x},t)$ does not depend on the coordinate in this case. If in equation (69) we substitute the particles' momentum distribution function f ₁(p , t) by the distribution function in the velocity, $f_1(v,t)$, $\mathbf{p} = m \mathbf{v}$, then the equation takes the form usual for the case of quasi-Brownian particles with active fluctuations (see, e.g. [3, 20, 21]). At the same time, the second formula in equation (69) connects the intensity of the 'momentum' noise $D_{\rm p}$ introduced here with the intensity of the 'velocity' noise $D_{\rm v}$ (see equations (9) and (10)). This implies that equation (66) may be regarded as a kinetic equation for quasi-Brownian particles with active fluctuations, which is generalized for the case of a three-dimensional system with dissipative interaction and a nonlocal external stochastic field.

The stationary solution $f_{\infty } \left(\,p\right)=\mathop{\lim }\nolimits_{t\to \infty } f_{1} \left(\,p,t\right)$ of equation (69) has a Boltzmann form,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_12_} f_{\infty } \left(\,p\right)=A\exp\left(-\frac{\gamma }{2D_{p} } p^{2} \right), \nonumber \end{align} \tag{ 70 }$$

which is different in two-dimensional and three-dimensional cases only by the value of the normalizing constant A (see equation (68)):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_13_} &\nonumber A=\frac{\gamma }{2\pi D_{p} \mathfrak{v}} \quad &\mathrm{for~ the~\quad two{\text -}dimensional ~case}, \nonumber \\ &A=\frac{1}{\mathfrak{v}} \left(\frac{2\pi D_{p} }{\gamma } \right)^{-3/2} \quad &\mathrm{for~ the \quad three{\text -}dimensional ~case.} \nonumber \end{align} \tag{ 71 }$$

Taking into account normalization ((68) and (71)) in the two-dimensional case, the formula (70) for the stationary distribution function of the active particles coincides with the corresponding expression in [3].

6.2. Brownian particles with active local fluctuations

We now investigate the spatially homogeneous stationary states of the system under study in the case of a spatially inhomogeneous external impact. As we already noted (see text explaining formulas (59) and (60), as well as figure 2), spatially correlated noise in certain cases does not prevent the existence of spatially homogeneous states in the system. Let us now consider the stationary solution $f_{\infty } \left(\,p\right)=\mathop{\lim }\nolimits_{t\to \infty } f_{1} \left(\,p,t\right)$ of equation (66), which is more general than equation (69). Now equation (66) in the limit $t\to \infty$ can be written as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_14_} \frac{\partial }{{\partial {p_i}}}e_i^{\rm h}\left\{{{f_\infty }\left(\,p \right)\gamma \left(\,p \right)\,p + {D_{\rm p}}\frac{{\partial {f_\infty }\left(\,p \right)}}{{\partial p}} + \tilde g{f_\infty }\left(\,p \right)} \right\} = 0, \nonumber \end{align} \tag{ 72 }$$

where we introduce

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_16_} g(0)\equiv 2D_{\rm p} \, , \qquad \tilde{g}\equiv \bar{g}\int {\rm d} {\mathbf p}_{2} \frac{\partial f_{\infty } \left(\,p_{2} \right)}{\partial p_{2} } , \nonumber \end{align} \tag{ 73 }$$

$\bar{g}$ is given by equation (64), and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_17_} \gamma \left(\,p_{1} \right)=2\mathop{\lim }\limits_{t\to \infty } \frac{\partial \bar{R}\left(\,p^{2} ,t\right)}{\partial p^{2} } =2\frac{\partial }{\partial p^{2} } \int {\rm d} {\mathbf p}_{2} f_{\infty } \left(\,p_{2} \right)\bar{R}\left(\left( {\mathbf p}- {\mathbf p}_{2} \right)^{2} \right) \nonumber \end{align} \tag{ 74 }$$

defines the nonlinear friction forces. The solution of equation (72) is defined by the relation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_15_} f_{\infty } \left(\,p\right)\sim \exp \left\{-\frac{1}{D_{\rm p} } \int_{}^{\,p}{\rm d}p' \left(\gamma \left(\,p'\right)\,p'+\tilde{g}\right)\right\}. \nonumber \end{align} \tag{ 75 }$$

The expression similar to equation (75) is given in [18] for the case of two-dimensional active systems with a nonlinear friction coefficient. We emphasize that in this paper, as well as in the general analysis of the above expression, the friction coefficient $\gamma(\,p)$ is considered to be an arbitrary function of the momentum (in fact, of the velocity $v$, $\gamma(v)$, since the authors of this work performed calculations in the space of velocities not momenta.)

The general analysis of this expression in [18] with an arbitrary dependence of the coefficient of friction on velocity shows that it is the nonlinearity of this coefficient that is one of the necessary conditions for the occurrence of self-propulsion in the system. Indeed, as seen from analysis of the case with nonlinear friction, for one velocity range the dissipation forces may be negative while within the other velocity range the forces may be positive. Thus, the particles with the velocities within the first interval slow down; consequently, this interval corresponds to regular friction. On the contrary, the particles with velocities within the second interval accelerate, that is, they experience propulsion. As a result, a symmetric distribution with two peaks at $\pm v_0$ is established. We here note that the stationary distribution cannot be isotropic in the velocity (or momentum) space due to the existence of a selected direction $\mathbf{e}_{0i}^h = \mathbf{v}_{0i}/v_0$ in the stationary state. Such distribution corresponds to the general solution of the stationary kinetic equation derived in [18] on the basis of the two-dimensional Langevin equation for active particles with a nonlinear friction coefficient. We also note that in [18] the case of the special dependence of the friction coefficient $\gamma(v)$ on the velocity is analyzed in detail, leading to a linear dependence on the velocity of the friction force itself, but with a shift in the velocity space. It is shown that such a case of linear friction does not prevent the appearance of self-movement in the system. As far as we know, this result was obtained for the first time in the mentioned work. Finally, we note that the special case of linear friction considered in [18] corresponds (both for two-dimensional and three-dimensional systems) to the case $\newcommand{\e}{{\rm e}} \gamma(v)\equiv \gamma$ of linear friction in formula (75) of this work. For more details, see below.

Thus, as will be shown below, it follows from the solution (75) of equation (72) that the stationary distribution function with two maxima (self-propelled particles) can also be realized in the case of linear friction, namely, when $\newcommand{\e}{{\rm e}} \gamma \left (\,p \right) \equiv \gamma> 0$ (see equation (68)). This is due to the local impact of stochastic forces with active fluctuations. This leads to the appearance of the factor $\tilde{g}$ in equation (75), which is defined by the pair correlation functions of the external force and the distribution function itself. This is why the solution of equation (72) is a complex self-consistent task. In fact, the general solution in the case of linear friction, as follows from equation (75), is given by the formula

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_18_} &\nonumber f_{\infty } \left(\,p\right) \sim {} \exp \left\{-\frac{\gamma }{2D_{\rm p} } \left(\,p+\frac{\tilde{g}}{\gamma } \right)^{2} \right\},\nonumber \\ &\nonumber \tilde{g}\equiv {} \bar{g}\int {\rm d} {\mathbf p}_{2} \frac{\partial f_{\infty } \left(\,p_{2} \right)}{\partial p_{2} },\nonumber \\ &\bar{g}\equiv {} \frac{1}{\mathfrak{v}} \int {\rm d} {\mathbf x} g\left( {\mathbf x}\right). \nonumber \end{align} \tag{ 76 }$$

The display of the head–tail asymmetry is related to the sign of $\tilde {g}$. Namely, since $\gamma> 0$, the positivity of this value, $\tilde {g}> 0$, must comply with a purely dissipative case. When $\tilde {g} <0$, there are values of momenta for which the inequality $\gamma p + \tilde {g} <0$ is true. Particles moving with the corresponding momentum exhibit self-propulsion. Naturally, the very possibility of the existence of a stationary state presupposes the establishment of a stationary distribution combining both described cases. One of the possible options for achieving such a stationary state in spatially homogeneous systems is the establishment of the regime of two differently directed motions of its structural units with an average momentum p ₀ along some selected direction $\newcommand{\e}{{\rm e}} \mathbf{e}_0^h \equiv \mathbf{p_0}/p_0$. It is this version of the stationary distribution that would correspond to the solution described in [18] for the two-dimensional case. However, the appearance of the selected direction in the momentum space indicates that the stationary state of the system ceases to be isotropic in the momentum space. For this reason, the formal construction of the corresponding stationary distribution on the basis of solution (76) of equation (72) would be incorrect. The reason is that equation (72) was obtained from equation (66) precisely under the assumption of isotropy of the momentum space.

Therefore, in describing the stationary distribution discussed above, it is necessary to start from the most general equation (63). We recall that this equation characterizes the evolution of spatially homogeneous states of the system. In view of the foregoing, the stationary distribution under the steady-state regime of multidirectional movements along one direction should be described by the function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:new77} {f_\infty }\left( {\bf{p}} \right) \equiv {f_\infty }\left( {{p_\parallel }} \right), \quad {p_\parallel } \equiv {\bf{pe}}_0^{\rm h} = p{{\bf{e}}^{\rm h}}{\bf{e}}_0^{\rm h}, \quad {{\bf{e}}^{\rm h}} \equiv {{\bf{p}}/p}, \quad {\bf{e}}_0^{\rm h} \equiv {{{{\bf{p}}_0}} / {{p_0}}}. \nonumber \end{align} \tag{ 77 }$$

We emphasize that the stationary mode defined by formulas (77)—the possibility of which is predicted by formulas (75) and (76)—is one of the multiple modes possible for the system, not the only one. In this sense, it is one of the particular stationary solutions of the general equation (63) under the assumption (77). An equation for which $f_\infty(\,p_\parallel)$ holds follows from equation (63) within a limit $t \to \infty$ and respects relations (77). It has a form similar to equation (72) if in the latter we assume $\newcommand{\e}{{\rm e}} \gamma \left(\,p \right) \equiv \gamma > 0$ and the absolute values of momentum we replace with $p_\parallel$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:new78} \left( {\gamma {p_\parallel } + \tilde g} \right){f_\infty }\left( {{p_\parallel }} \right) + {D_{\rm p}}\frac{{\partial {f_\infty }\left( {{p_\parallel }} \right)}}{{\partial {p_\parallel }}} = 0, \qquad {D_{{{\rm p}_\parallel }}} \equiv D_{\rm p} \equiv \frac{1}{2}g\left( 0 \right). \nonumber \end{align} \tag{ 78 }$$

For this reason, an analysis of the solutions of this equation can be carried out by analogy of the analysis of formulas (75) and (76) given above (see also [18]). Here, the momentum absolute values are to be replaced by the quantity $p_\parallel$, and we assume $\newcommand{\e}{{\rm e}} \gamma \left(\,p \right) \equiv \gamma > 0$.

Thus, taking into account the foregoing considerations for the case of a steady-state regime of multidirectional motions along $\newcommand{\e}{{\rm e}} {\bf{e}}_0^{\rm h} \equiv {{{{\bf{p}}_0}} /{{p_0}}}$, the stationary distribution function of active particles has the form (see [3, 18, 21])

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_19_} \nonumber {f_\infty }\left( {{p_\parallel }} \right) = C\left\{{\vartheta \left( {{{\bf{e}}^{\rm h}}{\bf{e}}_0^{\rm h}} \right)\exp \left[ {- \frac{\gamma }{{2{D_{{{\rm p}_\parallel }}}}}{{\left( {{p_\parallel } - {p_0}} \right)}^2}} \right]} \right.\nonumber \\ \quad \quad \quad \quad + \left. {\vartheta \left( {- {{\bf{e}}^{\rm h}}{\bf{e}}_0^{\rm h}} \right)\exp \left[ {- \frac{\gamma }{{2{D_{{{\rm p}_\parallel }}}}}{{\left( {{p_\parallel } + {p_0}} \right)}^2}} \right]} \right\}, \nonumber \end{align} \tag{ 79 }$$

where $\vartheta \left( x \right)$ is the Heaviside step function, and C is the normalization constant. Momentum p ₀ in equation (79), characterizing the location of maxima of the distribution function symmetric with respect to the point p   =  0, is determined by $\tilde {g}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_20_} p_{0} ={\left|\tilde{g}\right|\mathord{\left/ {\vphantom {\left|\tilde{g}\right| \gamma }} \right. } \gamma } . \nonumber \end{align} \tag{ 80 }$$

The value of $\tilde {g}$ itself, according to equations (73) and (76), depends on the derivative of the unknown momentum distribution function. Thus, the definition (76) with the explicit form of the distribution function (79) should be considered as an equation that connects $\tilde {g}$ to the normalization constant C:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_21_} \nonumber \tilde g = - C\frac{{\gamma \bar g}}{{2{D_{{{\rm p}_\parallel }}}}}\int {{\rm d}{\bf{p}}\left\{{\vartheta \left( {{{\bf{e}}^{\rm h}}{\bf{e}}_0^{\rm h}} \right)\left( {{p_\parallel } - {p_0}} \right)\exp \left[ {- \frac{\gamma }{{2{D_{{{\rm p}_\parallel }}}}}{{\left( {{p_\parallel } - {p_0}} \right)}^2}} \right] }\right.}\nonumber \\ +{\left.{\vartheta \left( {- {{\bf{e}}^{\rm h}}{\bf{e}}_0^{\rm h}} \right)\left( {{p_\parallel } + {p_0}} \right)\exp \left[ {- \frac{\gamma }{{2{D_{{{\rm p}_\parallel }}}}}{{\left( {{p_\parallel } + {p_0}} \right)}^2}} \right]} \right\}}. \nonumber \end{align} \tag{ 81 }$$

In turn, the constant C is determined from the normalization condition (see equation (68)

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \int {\rm d} {\mathbf p} f_{\infty } \left(\,p_\parallel\right)=\frac{1}{\mathfrak{v}} ,\nonumber \end{align*}$$

which can be rewritten after combining with equation (79) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_22_} \nonumber \frac{1}{\mathfrak{v}} = C\int {{\rm d}{\bf{p}}} \left\{{\vartheta \left( {{{\bf{e}}^{\rm h}}{\bf{e}}_0^{\rm h}} \right)\exp \left[ {- \frac{\gamma }{{2{D_{{{\rm p}_\parallel }}}}}{{\left( {{p_\parallel } - {p_0}} \right)}^2}} \right] } \right. \nonumber \\ + \left. {\vartheta \left( {- {{\bf{e}}^{\rm h}}{\bf{e}}_0^{\rm h}} \right)\exp \left[ {- \frac{\gamma }{{2{D_{{p_\parallel }}}}}{{\left( {{p_\parallel } + {p_0}} \right)}^2}} \right]} \right\}. \nonumber \end{align} \tag{ 82 }$$

The latter expression is also an equation relating the constant C and the unknown quantity $\tilde {g}$. Thus, equations (81) and (82) represent a system of two equations to determine two unknown quantities, C and $\tilde {g}$, in terms of parameters characterizing the system, namely, the friction coefficient $\gamma$, the number density of particles $1 / \mathfrak{v}$, and the parameters of noise with active fluctuations, i.e. the pair correlation function $\bar {g}$ and $g \left (0 \right) = 2D_ {\rm p}$ (see equations (73) and (76)). We note that, in contrast to the magnitude of the momentum p ₀, the direction $\newcommand{\e}{{\rm e}} {\bf{e}}_0^h \equiv {{{{\bf{p}}_0}} /{{p_0}}}$ of steady motion (the self-propelled motion) in the framework of a spatially homogeneous model of a medium with active particles cannot be determined. To determine it, it is necessary to break the spatial isotropy of the system, for example, by introducing anisotropic interaction between particles, even an infinitely small interaction. This will also be mentioned in Conclusion. We also note that because of the integration with respect to the total volume in momentum space, equations (81) and (82) have different forms for two- and three-dimensional cases. Furthermore, the calculation of the integrals over $\mathbf{p}$ in equations (81) and (82) must take into account the anisotropy of the integrands in momentum space. The procedure for calculating these integrals in a spherical coordinate system in momentum space is greatly simplified if the polar axis is assumed to be directed along the vector $\mathbf{e}_0^{\rm h}$, see (77).

We first consider the two-dimensional case. Equations (81) and (82) take the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_23_} &\nonumber \tilde{g}=-2\pi ^{3/2} \sqrt{\frac{2D_{\rm p} }{\gamma } } \bar{g}C, \nonumber \\ &\frac{1}{\mathfrak{v}} =C\left\{4\pi \frac{D_{\rm p} }{\gamma } +2\pi ^{3/2} p_{0} \sqrt{\frac{2D_{\rm p} }{\gamma } } \mathop{\mathrm{erf}}\left(\,p_{0} \sqrt{\frac{\gamma }{2D_{\rm p} } } \right)\right\}, \quad p_{0} \equiv \frac{\left|\tilde{g}\right|}{\gamma } , \nonumber \end{align} \tag{ 83 }$$

where $\mathop{\mathrm{erf}}(x)$ is the error integral:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_24_} \mathop{\mathrm{erf}}(x)\equiv \frac{2}{\sqrt{\pi } } \int_{0}^{x}{\rm d}y \exp \left(-y^{2} \right). \nonumber \end{align} \tag{ 84 }$$

In general, equations (83) are complex transcendental equations that can be solved numerically. However, in the two extreme cases, namely those of small and large values of the argument $p_ {0} \sqrt {{\gamma \mathord {\left / {\vphantom {\gamma 2D_ {\rm p}}} \right. } 2D_ {\rm p}}} = \sqrt {{\tilde {g} ^ {2} \mathord {\left / {\vphantom {\tilde {g} ^ {2} 2 \gamma D_{\rm p}}} \right. } 2 \gamma D_ {\rm p}}}$ of the error integral (84), these equations can be solved analytically. In the case $p_ {0} \sqrt {{\gamma \mathord {\left / {\vphantom {\gamma 2D_ {\rm p}}} \right. } 2D_ {\rm p}}} = \sqrt {{\tilde {g} ^ {2} \mathord {\left / {\vphantom {\tilde {g} ^ {2} 2 \gamma D_{\rm p}}} \right. } 2 \gamma D_ {\rm p}}} \ll 1$, these solutions are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_25_} \nonumber & C\approx \frac{\gamma }{4\pi D_{\rm p} \mathfrak{v}}, \nonumber \\ \nonumber & \tilde{g}=-\pi ^{1/2} \frac{\bar{g}}{\mathfrak{v}} \sqrt{\frac{\gamma }{2D_{\rm p} } } , \nonumber \\ & p_{0} \approx \pi ^{1/2} \frac{1}{\mathfrak{v}} \sqrt{\frac{\bar{g}^{2} }{2\gamma D_{\rm p} } } . \nonumber \end{align} \tag{ 85 }$$

With equation (85), the inequality above can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_26_} \frac{\bar{g}}{D_{\rm p} \mathfrak{v}} \ll 1. \nonumber \end{align} \tag{ 86 }$$

We have already mentioned that the presence of self-propelled properties depends on the sign of $\tilde {g}$. Since $\bar{g}$ is positive (see equation (76)), the quantity $\tilde {g}$ in equation (85) is negative, i.e. $\tilde {g} < 0$. Consequently, a state like a distribution function with two bell-like peaks and parameters defined by equation (85) can be realized.

Now consider the opposite case $p_{0} \sqrt{\gamma / 2D_{\rm p} } =\sqrt{\tilde{g}^{2} / 2\gamma D_{\rm p} } \gg 1$. Then, in the main order we get from equation (83):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_27_} &\nonumber C\approx \frac{1}{4\pi ^{3/2} \mathfrak{v}} \frac{\gamma }{D_{\rm p} } \sqrt{\frac{2\gamma D_{\rm p} }{\tilde{g}^{2} } } , \nonumber \\ &\nonumber \tilde{g}=-2\pi ^{3/2} \sqrt{\frac{2D_{\rm p} }{\gamma } } \bar{g}C, \nonumber \\ &p_{0} =\sqrt{\frac{\bar{g}}{\gamma \mathfrak{v}} } . \nonumber \end{align} \tag{ 87 }$$

We see that also in this case the quantity $\tilde g$ is negative, i.e. $\tilde g < 0$, and the particle distribution function is defined by equations (79) and (87). We also add that the inequality $p_ {0} \sqrt {{\gamma \mathord {\left / {\vphantom {\gamma 2D_ {\rm p}}} \right. } 2D_ {\rm p}}} = \sqrt {{\tilde {g} ^ {2} \mathord {\left / {\vphantom {\tilde {g} ^ {2} 2 \gamma D_ {\rm p}}} \right. } 2 \gamma D_ {\rm p}}} \gg 1$, can be transformed into the relation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_28_} \frac{\bar{g}}{D_{\rm p} \mathfrak{v}} \gg 1 \nonumber \end{align} \tag{ 88 }$$

with the use of equation (87). This result is the opposite of equation (86).

Now let us return to equations (81) and (82) and study their solution in the case of a three-dimensional system of active particles with linear friction and active space-dependent (i.e. local in space) fluctuations. In this case, equations (81) and (82) are transformed to a form substantially different from equation (83):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_29_} &\nonumber C=\frac{1}{2v} \left(\frac{\pi \gamma }{2D_{\rm p} } \right)^{3/2} , \nonumber \\ &\nonumber \tilde{g}+2\frac{\left|\tilde{g}\right|}{\gamma } \frac{\pi ^{3} \bar{g}}{\mathfrak{v}} \frac{\gamma }{D_{\rm p} } \mathop{\mathrm{erf}}\left(\sqrt{\frac{\tilde{g}^{2} }{2\gamma D_{\rm p} } } \right)=-8\frac{\pi \bar{g}D_{\rm p} }{\gamma \mathfrak{v}} \left(\frac{\pi \gamma }{2D_{\rm p} } \right)^{3/2} , \nonumber \\ &p_{0} \equiv \frac{\left|\tilde{g}\right|}{\gamma } , \nonumber \end{align} \tag{ 89 }$$

where $\mathop{\mathrm{erf}}(x)$ is still given by equation (84). As in the previous case of the two-dimensional system of active particles, the second equation in equation (89) in general form can be solved only numerically. However, in the two extreme cases examined above, this equation can be solved analytically. Namely, in the case of small values of the parameter $p_{0} \sqrt{{\gamma \mathord{\left/ {\vphantom {\gamma 2D_{\rm p} }} \right. } 2D_{\rm p} } } =\sqrt{{\tilde{g}^{2} \mathord{\left/ {\vphantom {\tilde{g}^{2} 2\gamma D_{\rm p} }} \right. } 2\gamma D_{\rm p} } } \ll 1$ the solution of equation (89) is determined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_30_} \nonumber \tilde{g}\approx -4\pi ^{2} \frac{\bar{g}}{\mathfrak{v}} \left(\frac{\pi \gamma }{2D_{\rm p} } \right)^{1/2}, \nonumber \\ p_{0} \approx 4\pi ^{2} \frac{\bar{g}}{\mathfrak{v}D_{\rm p} } \left(\frac{\pi D_{\rm p} }{2\gamma } \right)^{1/2}. \nonumber \end{align} \tag{ 90 }$$

Furthermore, as is easily seen directly, the ratio $p_{0} \sqrt{{\gamma \mathord{\left/ {\vphantom {\gamma 2D_{\rm p} }} \right. } 2D_{\rm p} } } =\sqrt{{\tilde{g}^{2} \mathord{\left/ {\vphantom {\tilde{g}^{2} 2\gamma D_{\rm p} }} \right. } 2\gamma D_{\rm p} } } \ll 1$ can be reduced to the form equation (86), which is valid for the two-dimensional case.

At large values of the parameter $p_{0} \sqrt{\gamma / 2D_{\rm p} } =\sqrt{\tilde{g}^{2} / 2\gamma D_{\rm p} }$, $p_{0} \sqrt{\gamma / 2D_{\rm p} } =\sqrt{\tilde{g}^{2} / 2\gamma D_{\rm p}} \gg 1$ the solution of equation (89) is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_31_} \tilde{g}\approx -8\frac{\pi \bar{g}}{\mathfrak{v}} \left(\frac{\pi \gamma }{2D_{\rm p} } \right)^{3/2} \frac{D_{\rm p} }{\gamma } \left(1-2\pi ^{3} \frac{\bar{g}}{\mathfrak{v}D_{\rm p} } \right)^{-1} , \nonumber \end{align} \tag{ 91 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle p_{0} \approx 4\pi ^{2} \frac{\bar{g}}{\mathfrak{v}D_{\rm p} } \sqrt{\frac{\pi D_{\rm p} }{2\gamma } } \left|1-2\pi ^{3} \frac{\bar{g}}{\mathfrak{v}D_{\rm p} } \right|^{-1}. \nonumber \end{align} \tag{ 92 }$$

By analyzing equation (91) we find that at $\bar{g}>0$ in the domain of large values of the parameter, the negative $\tilde{g}$ is only possible if $2\pi ^{3} \frac{\bar{g}}{\mathfrak{v}D_{\rm p} } \mathop{<}\limits_{\sim } 1$. In this case the expression for p ₀ can be simplified:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_32_} p_{0} \approx \sqrt{\frac{2D_{\rm p} }{\pi \gamma } } \left|1-2\pi ^{3} \frac{\bar{g}}{\mathfrak{v}D_{\rm p} } \right|^{-1} , \nonumber \end{align} \tag{ 93 }$$

and the condition $\sqrt{{\tilde{g}^{2} \mathord{\left/ {\vphantom {\tilde{g}^{2} 2\gamma D_{\rm p} }} \right. } 2\gamma D_{\rm p} } } \gg 1$ can be replaced with an equivalent one,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{GrindEQ__5_33_} 0<\frac{\bar{g}}{\mathfrak{v}D_{\rm p} } \mathop{<}\limits_{\sim } \frac{1}{2\pi ^{3} }. \nonumber \end{align} \tag{ 94 }$$

Thus, we conclude that self-propelled particles can be realized in the case of large values of the parameter $\sqrt{{\tilde{g}^{2} \mathord{\left/ {\vphantom {\tilde{g}^{2} 2\gamma D_{\rm p} }} \right. } 2\gamma D_{\rm p} } }$ and in the three-dimensional system; however, the criteria (94) in this motion differ significantly from that in the two-dimensional case (see equation (88)).

It should be noted that if $\bar {g} = 0$ then the characteristic momentum p ₀ is always zero, in both two- and three-dimensional cases; see equation (80) and the original equation (72). This should have been expected since this case corresponds to the degeneration of the two bell-like peaks single-particle distribution function into a symmetrical one (with respect to p   =  0) of the Gaussian type, with parameters that coincide with those of equations (70) and (71).

It is necessary to emphasize once again that the example we considered in detail above is one of the special cases of the stationary solutions of equation (63). It is clear that the stationary states of the system, which satisfy the stationary variant of equation (63), can also have a more complex form while remaining spatially homogeneous. We have chosen this type of solution (see (76)–(94)) precisely because it reproduces formulas that demonstrate the self-propelled properties of the systems under consideration, which are similar to the self-propelled results, for example, of [3, 18, 21]. It should be noted here that the self-propulsion in active media in this paper is described as a collective effect and a consequence of the influence of spatially correlated noise. But at the same time, for the manifestation of the self-propulsion effect in the system, the particles (even if they are considered to be points) should have the property of head–tail asymmetry. That is, one of the reasons for the possible manifestation of the properties of the self-propelled system are the individual characteristics of the particles—the structural units of the system. In our model, this was ensured by the dependence of the influence of noise on the magnitude and direction of the velocity of the particle itself. Thus, with the external similarity of the results of the present work and, for example, the work [18], the reasons leading to these results may be different. Recall that in [18] one particle is considered, which is logical for an approach that uses the Langevin equations as the initial equations for the individual motion of particles. The driving force arising in such a system is not a collective effect. However, it should still be borne in mind that a phenomenological approach was used in [18], where a nonlinear friction force in the medium was assumed, and a very specific dependence of the friction coefficient on speed was postulated.

In concluding this section, it is necessary to reiterate the following. In a number of places above, we specifically separated the cases of two-dimensional and three-dimensional systems. From the procedure for constructing the kinetic theory of active media in the microscopic approach described in the paper, it can be seen that there is no conceptual difference between the derivation of equations in two- and three-dimensional cases. However, the situation looks this way only after the equations of motion (6)–(9) are written out. These equations are formulated to replace the two-dimensional Langevin equations, which usually underlie the construction of the kinetic theory of active media. The generalization of the initial equations of motion for structural units, as could be seen, is completely nontrivial. In fact, as was noted in Introduction, historically, when constructing the kinetic theory of active media, authors, as a rule, proceeded from the two-dimensional Langevin type equation (10). In some cases, the interactions between the structural units of the system were included in the aforementioned Langevin equations; see, for example, [1–6, 19–21]. However, building a consistent microscopic approach to deriving the kinetic equations of such media requires going beyond the scope of the Langevin equations as the initial equations of motion for the structural units of the system. Instead of the Langevin equations, it is necessary to formulate the generalized Hamilton equations, that is, to start from the first principles of classical theoretical mechanics. For this reason, in this paper we proceeded from the generalized Hamilton equations, which naturally allow us to include both interactions between particles (including dissipative interactions) and a three-dimensional view of the stochastic effect on the system, depending on the speed and coordinates of structural units. Furthermore, in order to compare the results that follow from our theory in certain particular cases with the results of other authors (for two-dimensional systems), we separated in some cases the results for three-dimensional and two-dimensional systems. It should also be noted that in the framework of the model we are considering, the conditions (or criteria) for the existence of self-propulsion in the system may differ for two-dimensional and three-dimensional systems (see, for example, (88) and (94)). At the same time, the very possibility of the existence of such self-propulsion modes in both two-dimensional and three-dimensional cases is not prohibited.

In this paper we propose a microscopic approach to construct a kinetic theory of many-particle systems with dissipative and potential interactions in the presence of active fluctuations. The approach is based on a generalization of the Bogolyubov–Peletminskii reduced description method applied to the systems of many active particles. The microscopic approach allows us to construct the kinetic theory of active particles in the case of two-dimensional and three-dimensional systems, as well as in the presence of a nonlinear friction and local active fluctuations. We obtain general kinetic equations for these systems in the case of a weak interaction between the particles (both potential and dissipative) and low-intensity active fluctuations. We define some particular cases in which the derived kinetic equations have solutions that match with the results for the systems of active particles known from earlier works by other authors. It is also shown that one of the consequences of the local nature of the active fluctuations is a manifestation of a head–tail asymmetry and a self-propulsion, which is typical for systems of active particles, even in the case of a linear friction (see equations (83)–(94)).

We recall that in this context formulas (83)–(94) describe only two special limiting cases of two-dimensional and three-dimensional systems with self-propelled particles. The expressions obtained for a two-dimensional case coincide with those derived in [3, 21]. However, in our approach the nature of the phenomenon of self-propulsion is associated with a local (individual) exposure to the particles of the external stochastic field with active fluctuations, see equation (9). Besides, the parameters of the self-propelled motion are self-consistently expressed by the internal characteristics of a many-body system, namely the number density of particles in the system, the parameters of the dissipative function, and the characteristics of the external influence, that is, the pair correlation function of the active random field. Note that the stationary direction of the head–tail asymmetry within the spatially homogeneous model (see equations (65) and (66)) cannot be determined. To define it, we should introduce an interaction, even an arbitrarily small one; however, this violates the spatial isotropy of the problem. In this sense, the situation expressed by equations (79)–(88) is similar to the situation with a phase transition to the magnetic ordering in ferromagnets, see, e.g. [42]. As is known, the value of the total magnetic momentum in a ferromagnet in the main approximation is determined by the isotropic exchange interaction. The direction of the magnetization is at the same time given by nonisotropic weak relativistic interactions. The situation with the appearance of anisotropy in the momentum space as a result of a phase transition in a normal Fermi-liquid with a weakly anisotropic interaction was also considered in [43].

Here we note that the kinetic equations (55)–(57) are general in the sense that they describe quite a number of different many-particle systems, both two-dimensional and three-dimensional ones, with active local fluctuations, homogeneous and inhomogeneous space, and different variations in the nonlinear friction. The study of various particular cases of solutions of kinetic equations (55)–(57) is beyond the scope of this paper. As outlined above, the main objective of this work is the development of a microscopic approach to derive the general kinetic equations for active particles with nonlinear friction under the influence of active fluctuations, including a generalization to the case of three-dimensional systems.

We also note that the microscopic approach suggested here for the construction of the kinetic theory of many-particle systems with dissipative interaction and active fluctuations allows for further generalization, e.g. to the case of simultaneous presence of both active and passive fluctuations. Furthermore, a nonGaussian stochastic nature of external effects can be taken into account; see also in this regard [18].

The authors would like to thank Prof A Chechkin for fruitful discussions and valuable comments on the paper.

OS acknowledges the support by the Spanish Ministry of Economy and Competitiveness MINECO through BCAM Severo Ochoa excellence accreditation SEV-2013-0323 and SEV-2017-0718.