Beginnings of Developing Kinetic Scenarios of Plasma Evolution due to Coulomb Collisions

A new logic of reducing the two-time formalism to a highly informative scenario of redistribution of plasma particles in momentum due to Coulomb collisions is reported. Based on objective plasma evolution equations following from a properly reduced full plasma description, it has a more sound foundation than that presented in the previous report on increasing the informativeness of scenarios of the phenomenon. The possibilities of adapting the approach to the further development of more informative scenarios of plasma collisional relaxation and the modelling of transport phenomena are discussed.


Introduction
In physics research, evolving systems are of considerable interest to theorists. Note that even a momentary change in a factor affecting a certain physical system cannot cause an instantaneous transition of the system from an initial state to some final one: the transition from the first to the second state takes some time and, thus, involves the evolution of the system. Meanwhile, the final state itself is often determined by the subtleties of the system evolution scenario. In view of this, the development of evolution scenarios is of particular importance in the above-mentioned studies of evolving systems. In turn, its importance has led to the emergence of the concept of the informativeness of physical theoretical scenarios of system evolution: the longer a scenario adequately depicting the real picture of the macrophysical evolution of a system, the higher the informativeness of the scenario [1][2][3][4][5][6]. Regarding plasma science, we have found that the bulk of the kinetic scenarios of plasma phenomena generated by predecessors have inappropriately low informativeness due to the lack of an adequate understanding of two important aspects. The first of them is the asymptotic nature of the convergence of successive approximations to a plasma scenario. As a matter of fact, modelling plasma kinetics involves reducing the full plasma description that is given by the simultaneous Maxwell [7] and Klimontovich [8,9]-Dupree [10] equations to a simpler kinetic plasma model. In view of this, the researcher inevitably loses the ephemeral possibility of adequately modelling the plasma evolution during an infinite time interval, since any of the possible schemes of this reduction imply a significant reduction in the informational basis of the theory. Further, the reduction inevitably involves the generation of some nonlinear perturbation theory. In any such theory, an increase in the order of consideration entails a factorial increase in the number of terms considered. Consequently, after a certain order of consideration, an improvement in the accuracy of the scenario (which is the only reason for the increase in the order of consideration) will inevitably be replaced by a reduction in accuracy: successive iterations begin to diverge due to a growth in the number of terms. It is this feature of the behaviour of successive approximations of perturbation theory that is termed "asymptotic convergence". Increasing the order of consideration within the framework of the corresponding perturbation theory, the researcher can equally rigorously develop different scenarios of the plasma evolution from a specified initial state by referring to different versions of the leading-order approximation that correspond to this initial state. The point is that, when one takes different lowest-order approximations of the employed nonlinear perturbation theory, its first sequential orders converge to different conditional limits that correspond to different theoretical scenarios of the plasma macrophysical evolution.
The second aspect that was ignored by predecessors is that, with the traditional replacement of real plasmas by plasma ensembles, the theory loses a lot of its informativeness. Indeed, conclusions about the mutual influence of the statistics of the ensemble strongly depend on the composition of the ensemble; therefore, they cannot be treated as objective laws of the physical evolution of the system. Additionally, note that the above two aspects are inseparable of each other: the first implies the second and vice versa [3,4].
Having realised the reasons for the theory non-informativeness, we have formulated principles that help to reduce full plasma descriptions to kinetic models of plasma evolution with as high informativeness as possible. These principles are rejecting the traditional plasma ensemble averaging and the direct time integration of intermediate evolution equations that appear when reducing the full plasma description. Bearing them in mind, we have developed a technique that is suitable for studying phenomena in a weakly turbulent plasma. (In view of the purpose of the corresponding theoretical apparatus, it is most natural to refer to it as a highly informative correlation analysis of plasma kinetics.) Various aspects of this technique were discussed in References [3][4][5][6][11][12][13][14][15][16]. Some stages of its development are reflected in books [1,2]. Our recent contributions to the field are papers on drift in ordinary space and wavenumbers of non-potential plasma waves in an inhomogeneous plasma [6,17] and on the role of forced plasma oscillations that accompany such a drift [18].
Undoubtedly, weakly turbulent plasma phenomena are not the only ones for which it is possible to develop somewhat informative theoretical plasma scenarios. In general, the identification of plasma contexts that have this property and the development of tools for generating informative scenarios in respective physical situations should be an extremely important component of theoretical plasma research. We stress that there is no universal calculation technique: each plasma problem dictates its own logic of developing possibly more informative plasma kinetic scenarios. The simplest illustration of this thesis is the problem of thermodynamic relaxation of a weakly non-ideal nonequilibrium homogeneous plasma. (A plasma is considered to be weakly non-ideal if the number of particles in the Debye sphere N D = n e r 3 D notably exceeds unity: the smaller N D , the less ideal the plasma. Here n e is the density of plasma electrons and r D is the plasma Debye length.) In such a plasma, the particles undergo Coulomb collisions, which lead to plasma thermalisation. For this process, as for the phenomena in a collisionless weakly turbulent plasma, it is possible to develop informative scenarios. At the same time, the machinery for developing highly informative kinetic models of evolving weakly turbulent plasma is useless here.
Informative kinetic scenarios of plasma evolution due to Coulomb collisions can be important, say, for research on nuclear fusion. With this in mind, we have performed a study that aimed at creating suitable tools for developing relevant informative scenarios. A report on this study is currently under consideration for publication in Contributions to Plasma Physics [19]. We note that the leading order of our kinetic scenario of the plasma evolution due to Coulomb collisions of particles in homogeneous plasmas turned out to be consistent with the well-known Lenard [20]-Balescu [21] plasma kinetic equation. It should be stressed that Lenard and Balescu both considered this phenomenon within the plasma ensemble paradigm. The focus on probabilistic plasma ensembles predetermined the sets of ideas underlying the logics of their intermediate calculations. Obviously, a further extension of these logics turned out to be technically unfeasible: the history of publications on modelling the evolution of plasma due to Coulomb collisions does not contain reports on the development of kinetic equations that are superior in accuracy to the Lenard-Balescu equation.
In the above-mentioned research report, our goal was to take the effect of temporal changes in plasma distributions on the plasma evolution scenario into account. This goal was achieved on the basis of our two-time formalism [17] by an auxiliary integration of some formal equation that "governs" time advances of a formal function that is somehow similar to the two-particle distribution function in the concept of plasma kinetics after Bogoluibov [22]-Born-Green [23]-Kirkwood [24]-Yvon [25] (BBGKY kinetics). This integration helped to guess the structure of the real (well defined) function, the two-point correlation function, that governs time advances of one-particle distribution functions. The proposed function turned out to satisfy the natural equation of its evolution. Still, we have some internal dissatisfaction, since our key function of the theory only echoes, in its structure, with the result of the auxiliary calculation. On the one hand, the lack of full correspondence between this function and the result of the formal calculation is quite natural, since the formal function used in the calculation does not reflect the objective physical relations that are pertinent to the plasma under consideration. On the other hand, there should exist a direct derivation of the above key function based on the integration of the natural equation of the plasma evolution, whereas we had not succeeded in deriving it.
In this paper, we report a solution to the latter problem: we present a reasoning approach that is based on a full-scale two-time formalism with integration of the natural evolution equation of the two-point correlation function that strictly justifies the structure of the function.
Our new study will be presented in this paper in accordance with the following plan. In the next section, we briefly outline the principals of the two-time formalism. Section 3 uncovers the logic of consideration in paper [19] and it presents a slightly more clarified version of the two-point correlation function that has been developed according to this logic. Section 4 describes a completely new way of looking at the problem. The results of the study are commented on in the Conclusions.

Key Equations of the Two-Time Formalism
Stating the need to refrain from ensemble averaging, the theorist is faced with the problem of considering the full plasma description, the integration of which involves the simultaneous integration of the equations of motion for all individual charged plasma particles. On the one hand, such an integration presupposes the knowledge of the initial data on the position and momentum of each charged particle, which are never known. On the other hand, even if it was known, the task of simultaneously integrating infinitely many equations of motion is technically infeasible. This inevitably forces one to construct a Vlasov type statistic of particle distribution f α ( r, p, t) [26] from the Klimontovich's distribution N α = ∑ n δ 3 ( r − r n (t))δ 3 ( p − p n (t)) [8,9]: it was the construction of such statistics that was a consequence of the use of plasma ensemble averaging by our predecessors. (In the Klimontovich distribution, which we call the microdistribution, the subscript n numbers the particles of species α, and the functions r n (t) and p n (t) describe the trajectories of individual particles.) In the corresponding functional role, averaging over a plasma ensemble can only be replaced by contextually oriented averaging over a sufficiently large neighbourhood of the current point ( r, p) of the R ⊗ P phase space. (Here the sign ⊗ denotes the product of two spaces. In what follows, it will also denote the direct tensor product). In this paper, we consider the effect of Coulomb collisions in a macroscopically homogeneous plasma. For this, it is natural to use as the above neighbourhoods uniform six-dimensional parallelepipeds with sufficiently small momentum dimensions and sufficiently large spatial dimensions. As a result, the bulk of the averaged microdistribution f α turns in the sense of mathematics into a well-defined statistic that is independent of r; its transformations in time follow the equation (Here, and in what follows, we will use the notation of tensor analysis, although no freedom to change the reference frame is implied: co-and contravariant indices are distinguished for an unambiguous interpretation of formulae). The equation is written under the usual assumption that the potential part only of microstructural electric fields in the plasma is important in Coulomb collisions. The angle brackets denote averaging over the abovementioned parallelepiped-shaped neighbourhood of the product of the microstructural part of the microdistribution, δN α ( r, p, t) = N α ( r, p, t) − f α ( r, p, t), and the corresponding microstructural part E( r, t) of the electric field. The right-hand side of Equation (1) contains particular data on the two-point correlation function, δN α ( r, p, t) E( r , t ) . Here, we have intentionally split the time and spatial variables of the two objects under the averaging sign to simplify the calculation process. The notation indicates that the arguments r and p of δN α run through the neighbourhood in the phase space that was chosen to construct f α at the current point ( r, p), and the spatial arguments of E( r , t ) vary synchronously with the spatial argument of δN α : the difference between these arguments is fixed in the averaging. (In this regard, our concept of the twopoint correlation function differs substantially from the concepts of two-point functions that can be encountered in many papers that rely on traditional approaches. Particularly, we can point to the concept of phasestrophy (otherwise the two-point phase space density correlation [27]) mentioned in Reference [28]. Although the concept of the latter object involves some averaging in phase space (called averaging over the spatial coordinate of the centre of mass), it was defined as the technical implementation of averaging over an ensemble of plasmas. The reference to this object is motivated by recalling exotic objects, such as convective cells [29] and phase space density granulations [30][31][32]. Meanwhile, had these or any other exotic objects of the traditional plasma theory been of any importance for plasma physical manifestations, then they would have contributed to our two-point correlation function, and their cumulative effect on plasma manifestations would have been adequately described by the equations of two-time formalism provided that the corresponding iterative procedure shows some asymptotic convergence.) Hence the system of spatial arguments in our approach is somewhat redundant. (For more detailed comments on this issue, see References [6,17].) However, this system is very convenient due to its internal symmetry. The second equation of the two-time formalism is either the evolution equation for the two-point correlation function or the evolution equation for the two-time correlation function Φ Φ Φ( r, t, r , t ) (since the magnetic microfields generated by plasma particles can be neglected in the problem under consideration, the last function can be defined as Φ Φ Φ( r, t, r , t ) = E( r, t) ⊗ E( r , t ) ). For modelling the kinetics of plasma collisional relaxation, the use of the first equation is more reasonable. (The other equation, the evolution equation for the two-time correlation function, has been widely used in our previous studies of phenomena in weakly turbulent plasma.) A number of theoretical objects are required to represent this equation. The first is the bare Green function of particles of a given species α, 0 G α ( r, p, t, r , p , t ). It satisfies the causality principle: at t < t , the function is identically zero. In the domain t > t , the function under study evolves according to the equation with the initial data The second object is the electromagnetic Green function F F F , which conceptually corresponds to the well-known delayed potentials. This function is a definite integro-differential operator that gives an expression for the electromagnetic field (EMF) tensor in terms of the charge density. That is, in the absence of external electromagnetic radiation, the EMF tensor in the plasma is The explicit form of F F F is not needed for our considerations. An appropriately accurate equation for the evolution of the two-point correlation function δN α ( r, p, t) E( r , t ) was developed in [33], where its full-scale presentation (and derivation) was given in graphic form. We do not need to repeat this for the sake of the current study. We only note that, in Reference [33], this equation is written in Figure 15 after the graphical notation that was introduced in the paper. The equation governs the advances of δN α ( r, p, t) E( r , t ) with an increase in its entry time t. According to the causality principle, the equation is only consistent for t t . That is, had we known the initial value δN α ( r, p, t ) E( r , t ) , we could have integrated the equation to calculate the function at larger entry times t > t . The above-mentioned equation from [33] has sufficient accuracy for the proper modeling of the plasma kinetics in the next to the leading order of the expansion in N −1 D . The researcher can easily develop more precise analogs of this equation: the principles of corresponding derivations are set out in Reference [34]. Increasing the order of the corresponding derivations makes sense up to the order ∝ N D : higher orders diverge due to the asymptotic nature of the convergence of the theory.

Preliminary Explanations
Following Equation (1), we need the feature δN α ( r, p, t) E( r, t) to determine the rate of change of the distribution function. In order to obtain it, we will develop a twopoint correlation function δN α ( r, p, t) E( r , t ) as a natural solution to the corresponding evolution equation. In fact, we will concentrate our efforts on developing an approximation of the function that evolves following only the linear version of this equation, i.e., the one that is based on the terms containing a single wavy line in Figure 15 of paper [33]. This time, our goal will be to properly account for the terms that are due to the next to the next to the leading order of δN α ( r, p, t) E( r , t ) in its expansion in the ratio of the inverse plasma frequency to the characteristic time scale of the plasma evolution.
Note that in a quiet plasma that evolves exclusively due to Coulomb collisions, the ratio of inverse plasma frequency to the characteristic time scale of the plasma evolution is of the order 1/N D . Correspondingly, we will develop a solution to the above linearised evolution equation up to corrections scaled as the leading-order function multiplied by 1/N 2 D . The nonlinear corrections due to the terms that are represented by three wavy lines from the more precise version of the equation depicted in the above-mentioned figure have the same order. (The reader can develop the corresponding version using the ideas outlined in Reference [34]. An appropriate starting point for the prospective development of such a version is Figure 4 in that paper.) These corrections can be developed independently, and their proper versions will be obtained via iterations on the basis of the above-mentioned approximation of the two-point correlation function (we mean the versions of the nonlinear corrections that are most consistent with the motivation for the development of a plasma kinetic scenario with maximum informativeness). We will not discuss the corresponding iterations in the present paper.
Our approximate equation for the evolution of the two-point correlation function takes the analytical form Here, a key property of two-time functions is that they have rather sharp dependencies on the difference between their entry and exit time variables, whereas their dependencies on the half-sum of these variables are rather smooth. That is, the characteristic time scales of the first dependencies are of the order of the inverse electron plasma frequency ω −1 pe , whereas the time scales of the dependencies on the half-sum of the entry and exit time variables are of the order of the lowest of the inverse collision frequencies of plasma particles, which is the inverse frequency of electron-electron Coulomb collisions ν −1 ee (the smallness of ν ee as compared to ω pe is equivalent to N −1 D 1). A second property of these functions is that they only depend on the difference between the entry and exit spatial variables because of the macroscopic homogeneity of the plasma. Therefore, we will use the Fourier transform in space according to the conventions This yields The two-time correlation function satisfies the Poisson equation: Thus, Equation (5) is conceptually an equality that contains some integral operator acting on the two-point correlation function δN α E γ k ( p, t, t ). Using this equation, one can express the two-point correlation function The equation cannot be integrated for t < t : recall the causality principle. Meanwhile, the function δN α E γ k ( p, t, t ) is the spatial Fourier transform of a statistic that is well defined for any sequence of its entry and exit temporal arguments. Moreover, the more accurate evolution equation that is presented in Figure 15 from Reference [33] shows that the advances in time of the two-point correlation function δN α E γ k ( p, t, t ) depend not only on the data of the function with an entry time exceeding the exit time, but also on the data of the function with the reverse order of these times. That is, we should develop an approximation of the two-time correlation function δN α E γ k ( p, t, t ) for an arbitrary sequence of its entry and exit time arguments. The idea of this approximation was put forward in Reference [19] based on auxiliary integration of a formal equation defining the time advances of the formal function δN α δN α k ( p, t, p , t ). We stress that this function is not a statistic. Its spatial original (inverse Fourier transform) is represented as the average where the difference p − p is fixed similarly to the difference r − r while r and p run through the neighbourhood selected to define the distribution function f α at current point r, p. This is an extremely rugged function containing spikes that are scattered in the phase space, depending on the details of the microdistribution.
The "evolution equation" itself is written as Using the Laplace transform and its inversion, we have expressed δN α δN α k ( p, t, p , t 0 ) for an entry time t that satisfies the constraint ω pe −1 t − t 0 (ν ee ) −1 through the data δN α δN α k ( p, t 0 , p , t 0 ).
(Note: the effective expression of this data is 3 . It is rather rough: it does not provide an understanding of whether the plasma is drastically dynamic or whether it has reached a rather quiet state and evolves slowly further on. This "roughness" of the data is just the reason for the persistence of some uncertainty as a result of the formal calculation). Subsequently, we have used the Poisson equation to form the initial data of the two-point . After that, we have applied the Laplace transform to the Equation (5) and obtained some expression for the two-point The time delays t − t 0 and t − t 0 satisfying the above constraints are quite sufficient for the transition of the plasma to a relatively quiet state, even with an exotic initial arrangement of its particles in space (when the plasma is in a notably nonequilibrium state) and, at the same time, these time delays are too small for some significant changes in the distributions of plasma particles f α to occur. This ultimately yielded some formal expression for the two-point correlation function δN α E γ k ( p, t, t ) that depended on three time variables: the entry time t, the exit time t , and the variable t 0 , which stands for a somewhat arbitrary previous moment. Here, we will not present the expression itself. Instead, we present the result of a more accurate calculation following the same principles. By presenting the one, we pursue two goals. The first is to obtain a clearer understanding of the structure of the two-point correlation function, and the second is to further illustrate the consistency of the new logic of developing the two-point correlation function, which is discussed in the next section.
The new expression of the two-point correlation function is Here, is the linear dielectric function of the plasma, and we have also introduced notation The function G α kω (t , t 0 ) is well defined on the axis of real ω, but we intend to manipulate with its analytical continuation at least to a narrow strip in the plane of complex ω around the line Im ω = 0. The variable ω on the right-hand side of relation (8) is also implied to be a complex one. The line of integration over this variable goes from the left to right and it passes over all singularities of the integrand (note that the only important singularity of the integrand is the one due to the resonance ω = k · v ).
Note that the variable t enters the expression δN α E γ k ( p, t, t ) exclusively within the argument of the exponent: all other components of the expression depend on t and t 0 . The temporal variable t 0 stands for a somewhat arbitrary previous moment. Its presence in the expression is just due to the fact that corresponding calculation is purely formal. We have concluded that one should set t 0 = t in the case t > t and t 0 = t in the case t < t to obtain a reliable approximation of the two-point correlation function. The discussion shown in the next section will provide additional support for the specified value assignments to the variable t 0 .
Note that the above-mentioned formal calculation leading to relation (8) is a tedious task. For those who wish to repeat it, in Appendix A we list the basic mathematical relations that we used in this calculation.

Determining the Structure of the Two-Point Correlation Function from First Principles
We assume that the initial data of the two-point correlation function δN α E γ k ( p, t , t ) is known. Using the Laplace transform technique, we can integrate Equation (5) and obtain its solution for the time domain ω −1 pe t − t ν −1 ee . In this domain, the contributions to the inverse Laplace transform of the poles due to the zeros of the linear dielectric function ε kω (t ) are exponentially damped, which allows for obtaining an easily structured result. However, we will not follow the corresponding calculation program straightforwardly. Our calculations will include two stages. We will first develop an expression for the twotime correlation function Φ k Φ k Φ k (t, t ), and then derive a two-point correlation function that generates the last expression.

Two-Time Correlation Function
In view of the Poisson Equation (6), the two-time correlation function is nothing else than the integral of contributions from groups of particles that differ in momentum, and each contribution can be calculated well while ignoring the zeros of the dielectric function. Accordingly, it is possible to change the sequence of integration in momentum and frequency after the formation of the integral based on inverse Laplace transforms for differing groups of particles. That is, we can form the Laplace transform of Φ k Φ k Φ k (t, t ) and ignore the exponentially damped contributions due to the zeros of the dielectric function when inverting this Laplace transform.
The Laplace transforms of two-time functions will be introduced according to the rule and the inverse Laplace transforms, according to Here, C denotes the line in the complex plane ω that goes from left to right and passes over all singularities of the integrand. After the Laplace transform, Equation (5) takes the form Note: the solenoidal components of microfields in the plasma are assumed to be negligible, at least in terms of their contribution to the two-time correlation function. Therefore, the following equality is valid: Similarly, This makes it possible to simplify notation and, hence, the calculus. Namely, we denote This yields the following form of Equation (11) The symbol ∇ p denotes the derivative over the momentum p. Solving the equation for F α kω ( p, t ) and substituting the result into the Poisson equation, we obtain The solution to this new equation can be constructed by iterations. We have: Obviosuly, all of the components of h contain the summation over α with the weight 4πe α and the p -integration of F α k ( p, t , t ). Bearing this in mind, we only control the dependence of these components on ω: An equivalent form of the right-hand side here is Performing the inversion of the Laplace transform, we find that this form leads to the following structure Hence, for the time domain ω −1 pe t − t ν −1 ee , we have developed: Note: regardless of whether the plasma was initially drastically dynamic or not, we assume that, by the time t = t , it has acquired a quiet state with a rather large time scale of further evolution, which is of the order of the inverse electron-electron collision frequency ν −1 ee . Therefore, expression (13) should be reliable not only at t − t ω −1 pe , but also at smaller times t, up to t = t . Similarly, at times −ν −1 ee t − t < 0, the two-time correlation function should comply with the formula Naturally, when the time variable t passes through t = t , no discontinuity in the twotime correlation function should be observed. That is, with the growth of t, the two-time correlation function following Equation (14) is transformed after t = t into the one that corresponds to Equation (13). Accordingly, a compound form that fits the entire time where Q α is some function that only depends on the momentum p and time, and we have introduced the variable t 0 to denote the last dependence. This variable stands for the argument t for t > t and the argument t for t < t . The meaning of the function Q α can be stated, as follows. Setting t = t , we can directly calculate the two-time correlation function Φ Φ Φ( r, t, r , t ) for the case where the displacement | r − r | is small when compared to the mean interparticle distance in the plasma. At corresponding distances, the influence of neighbours on the field around a plasma particle can be neglected, and the function becomes the sum of the equal terms standing for individual charged particles: We can now modify this result to take plasma evolution into account on a time scale that is much shorter than the time of particle free flight through the mean interparticle distance. This time depends on electrons and, since the typical velocity of bulk electrons is of the order of the electron thermal velocity v Te , we obtain Using this formula, we can obtain the data of the spatial Fourier transform of the function for large wavenumbers at small time delays (k (n 0 ) 1/3 , 0 < t − t 1/(kv Te )) [35]: Naturally, the right-hand side here should satisfy Equation (15). Hence, there should be (We comment that the dielectric function ε k( k · v) is indistinguishable from unity for the above large wave vectors.)

Two-Point Correlation Function
In this subsection, our goal is to invent a procedure that allows us to reverse the relation that is given by the Poisson Equation (6), i.e., to develop a unique two-point correlation function δN α E γ k ( p, t, t ) that, in each situation, generates the given two-time correlation function. For the class of evolving plasmas of our interest, i.e., those that have reached a relatively quiet state at the moment under consideration, the latter problem has a logical solution.
We take a finalized expression for the two-time correlation function, (Here, we have replaced the integration variable p and the type of particles under the sign of summation α with p and α , respectively). Its entry time variable should be the entry time variable for the desired two-point correlation function δN α E γ k ( p, t, t ). We separate the structure in the pre-exponent of the integrand that depends on t 0 and t . In addition, in this structure, we will explicitly distinguish the dependence on k · v . The latter presupposes the use of the notation (10). This allows us to write We now expand the remaining part of the pre-exponent in powers of t − t up to the next to the next leading order: Integrating the right-hand side by parts, we reduce it to the form The last transformation in this chain of equations was carried out, as follows. In the integrand, we first separated the part of the pre-exponent that contains the second-order time derivative of ε kω −1 . The fact that the term with this derivative enters the expression without any differentiation over ω indicates that the previous calculation was correct. (This is in our order of consideration. Desiring to generate higher orders, one should make sure that in the above-mentioned pre-exponent in the term with the highest-order time derivative of the function ε kω −1 , this function is undifferentiated over ω.) Subsequently, we separated those terms with the first-order time derivative of ε kω −1 in which this derivative was not differentiated with respect to ω (the terms with both time-and ω-derivatives of ε kω −1 may well be present in the multipliers of these terms). Finally, terms remained that contained ε kω −1 without any differentiations. Now consider the last steps. In the developed formula, we make the following substitutions: 1 ε kω (t ) ∂ ∂t platform will form a natural basis for modelling plasma transport phenomena for classical high-temperature plasma. Thus, a basis is created for the development of more informative kinetic scenarios of the plasma evolution due to Coulomb collisions as compared to the previously formulated scenarios of this evolution.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Mathematical Tricks Used in the Formal Calculation Following the Logics of a Preceding Study
The main difficulty in the formal calculation of the two-point correlation function δN α E γ k ( p, t, t ) is to express it in terms of its "initial data" δN α E γ k ( p, t 0 , t ) where ω −1 pe t − t 0 ν −1 ee . Using the Laplace transform to integrate the corresponding evolution equation, we find two terms in the inverse Laplace transform. The first term is , and the second is the integral over ω, in which the integrand contains the integral over the dummy momentum p . The discussion below will focus solely on this integral. A significant contributions to this integral come from the poles originating from the multipliers 1/ ω − k · v and 1/ ω − k · v n . In the intended order of expansion in powers of the ratio of the inverse plasma frequency to the characteristic time of the plasma evolution, the exponent n takes values n = 1 − 3. That is, under the sign of integration over ω, we find structures of the type To express the contributions to the inverse Laplace transform due to these structures, we substitute them into the integrand according to the following presentations: Correctly approximating the inverse Laplace transform (with integration over ω by parts, if necessary), we obtain an integrand in the integral over dummy p that has no singularities at k · v = k · v . This makes it possible to treat the integral over k · v as a principal value integral. The initial data δN α E γ k ( p , t 0 , t ) in the p -integral is proportional to exp −i k · v (t 0 − t ) . The integration over ω provides either the exponent exp −i k · v (t − t 0 ) or the exponent exp −i k · v (t − t 0 ) , depending on the summand. Accordingly, under the sign of integration over k · v , we find terms with smooth functions of this variable that are multiplied by either exp i k · v − k · v (t − t 0 ) or exp −i k · v − k · v (t − t ) . Due to the natural structure of the distribution function, the p -integrals of these terms converge absolutely. In terms of the first type, the exponent is rapidly oscillating, which makes it possible to rearrange their integrands via the use of the following simplifying substitutions: