Extended hydrodynamical models for plasmas

We propose an extended hydrodynamical model for plasmas, based on the moments of the electron distribution function which satisfies the Fokker–Planck–Landau (FPL) transport equation. The equations for the moments can be obtained by multiplying the FPL equation by the corresponding weight functions and integrating over the velocity space. The moments are decomposed in their convective and non–convective parts and closure relations for the fluxes and production terms can be obtained by using the maximum entropy distribution function, which depends on Lagrangian multipliers. These latter can be expressed in terms of the state variables by imposing the constraints that the maximum entropy distribution function reproduces the moments chosen as state variables. In particular, we will concentrate on the 13-moment system. As a first application, we treat the case of the relaxation towards equilibrium of a homogeneous plasma with a temperature anisotropy, showing that the results are in good agreement with those obtained by means of the Kogan solution of the kinetic equation.


Introduction
The problem of collisions in plasmas has been widely studied since the first decades of the last century, starting with Landau's pioneering work in 1936 [9]. The binary collisions between charged particles in plasma physics are caused by the long-range Coulomb interactions and have been modeled through different operators, with different physical characteristics and mathematical structures [9,10,20]. Coulombic interactions in the absence of wave-particle resonance are well described by the Landau integral operator [9], which is a nonlinear, integro-differential, Fokker-Planck type operator that satisfies the H-theorem for entropy growth [7].
It is generally believed that in a weakly collisional plasma, such as the solar wind, collisions are too weak to produce any significant effect on the plasma dynamics. However, in order to consider particle heating and the consequent entropy growth, the collisional approach appears to be necessary [11,14,16]. In fact, because of the H-theorem, thermalization is uniquely due to collisions, which produce heating in general thermodynamic sense.
The main difficulty of a kinetic approach is that the numerical solution of the (FPL) equation requires a huge computational cost. A possible strategy is that of replacing the Landau collisional operator with the Dougherty operator which is a nonlinear differential operator of the Fokker-Planck type which requires a significantly lighter computational effort [5]. A comparison between the Landau and the Dougherty collision operator can be found in [16][17][18].
In this paper, we propose to adopt the alternative approach used by Grad in gas dynamics [6], which approximates the kinetic model by a 13-moment system, with closure relations based on the Maximum Entropy Principle (MEP) [13]. This principle has been successfully used in other fields of physics and applied mathematics such as radiative fluid dynamics [12] and charge and energy transport in semiconductors [1][2][3]15]. The main advantage of this approach is that its solution requires a much lower computational cost with respect to the numerical approximation of the kinetic equation, notwithstanding the use of the Landau operator.
The main objective of this paper is to derive the general moment system from the FPL equation coupled with the Maxwell equations and to underline a general procedure for obtaining closure relations for the extra moments and the moments of the collisional operator by applying the Maximum Entropy Principle. Furthermore, we specialize the above scheme to the 13-moment system, and in this case, we explicitly compute the closure relations. To test the validity of the model, we analytically solve the case of the relaxation towards equilibrium of a homogeneous plasma with a temperature anisotropy [16] and compare the results with those given by the special solution to the kinetic equation proposed by Kogan [8].
The paper is organized as follows: In Sect. 2, we introduce the FPL equation coupled with the Maxwell equations and we find the evolution equations for moments of electron distribution function f with weight functions which are tensorial products of the microscopic velocity. In Sect. 3, we define the internal moments, since their evolution equations are more suitable to be closed on the basis of the MEP. In Sect. 4, we consider the case of 13 moments, which have an immediate physical meaning, by explicitly computing the closure relations for the extra-fluxes and the production terms by expanding the MEP distribution function at the first order in anisotropy. Finally, in Sect. 5, we exploit a numerical test to check up to what extent the small anisotropy approximation gives good results, comparing them with those given by the Kogan solution [8], in the case of the relaxation towards equilibrium of a homogeneous plasma with a temperature anisotropy.

The Fokker-Planck-Landau equation and the moments equations
From a kinetic point of view, the dynamics of electrons in a plasma is described by means of their distribution function f (x, v, t). This function represents the density of electrons at time t in an elementary volume dxdv, around position x and velocity v. The distribution f satisfies the Fokker-Planck-Landau (FPL) transport equation coupled with the self-consistent Maxwell equations: Here, the Einstein convention about the sum over repeated indices is adopted, q is the charge of the carriers (with sign), v is the carrier velocity, E is the electric field, φ is the electric potential, B is the magnetic field, is the dielectric constant, μ is the permeability, n is the carrier number density (concentration) and u is the charge average velocity. The operator C[ f ] gives the time rate of change of f due to collisions, the prime denotes evaluation at v , as in f = f (x, v , t). Thus, (1) is a nonlinear integro-differential equation in seven independent variables.
An important property of the collision operator C[ f ] is the existence of physically relevant collisional invariants: This property implies conservation of the carrier number density, the average momentum and the average energy.
Another property of C[ f ] is the existence of an H-theorem, which can be stated in the following form: with T the temperature, m the particle mass and y = 2 m . This property establishes the existence of an entropy, i.e., h( f ) = − R 3 f (log f /y − 1)dv, and of equilibrium distributions, i.e., shifted Maxwellians, of the system. Macroscopic models can be obtained from the FPL equation by introducing moments of the distribution function which correspond to suitable choices of weight functions. One of the most common choices is that of taking tensor products of the microscopic velocity as weight functions, which gives the following moments The first few moments have an immediate physical interpretation, that is, M is the number density, m M i is the momentum density, 1 2 m M ii is the energy density, m M i j is the momentum flux, and 1 2 m M i j j is the energy flux. The evolution equations for these moments can be easily obtained by multiplying equation (1) by v i 1 v i 2 · · · v i k and integrating over R 3 , with a suitable vanishing condition for f as |v| → ∞. So doing, one gets with k ≥ 0 and production terms, and i jk the three-dimensional Levi-Civita symbol. Parentheses in suffices in (10) denote symmetrization with respect to the indices enclosed by them. In this way, there are only four independent variables, but the system is infinite. One needs to truncate the system at a finite order N , but the resulting system is not closed, since there appear extra variables which are the production terms and the flux M ri 1 i 2 ...i N . Constitutive equations are needed for these extra variables, which can be obtained by exploiting the MEP.

Internal moments and the Maximum Entropy Principle closure
For the system of the moments with 0 ≤ k ≤ N , the MEP states that the closure relations for the extra moments and for the moments of C[ f ] can be evaluated by using the distribution function which maximizes the entropy under the constraint that the leading moments (9) are preserved [13]. Since, as has been proved in [19], the entropy does not depend on the mean velocity u i = M i M , it is convenient to introduce the so-called internal or non-convective moments which are defined as follows: Also in this case, the first moments have a particular interpretation, that isM = M,M i = 0, mM ii = 2ne = 3 p, with e the specific internal energy and p the pressure, mM i j = P i j , with P i j the pressure tensor, and mM i j j = 2q i with q i the heat flux.
It can be shown [19] that the internal moments are related to the old ones through the relation ⎠ is the vector of all the moments taken into account, andM is the vector of the corresponding internal moments, where the bottom indices are column indices and the upper ones are row indices. The above-written relation can be retrieved directly from the definitions of M i 1 ...i k andM i 1 ...i k . However (11) relates all tensorial state quantities to their evaluation in the local rest frame, if the evolution evolution equations of these quantities are Galilean invariant [19], as is the case for the moments of the distribution function. The evolution equations for the internal moments can be obtained by substituting (11) into (10), which gives where C = X(u)Ĉ, with obvious meaning of C andĈ. At this point, in order to obtain the necessary constitutive relations, the first step is to find the extremum of the entropy density The solution is given by The physical interpretation of the Lagrange multipliersλ follows from the generalized Gibbs relation which holds as a consequence of the MEP [13].
Inserting (14) into the constraints (13) and solving the resulting equations with respect to the Lagrange multipliers, in principle it is possible to express the ME distribution as function of the momentsM. Eventually, substituting the so obtained f ME into the integrals defining the extra variables, the needed constitutive relations could be found. Actually, in general, inversion can be done only numerically; however if the system is not too far from local equilibrium that it is slightly anisotropic, the ME distribution (14) can be expanded around the local equilibrium values of the Lagrange multipliers, so simplifying the inversion. At the first order, the expanded f ME is strictly related to the Grad distribution which is widely used in extended thermodynamics. In the next section, we will present this procedure in the particular case when thirteen moments are used, the moments with an immediate physical interpretation.

The 13-moment model
In this section, we will explicitly write the moment system, including the closure relations, when the first thirteen moments with an immediate physical interpretation are taken as state variables. We recall that these moments are as follows: n =M, u i , P i j = mM i j , and q i = 1 2 mM ikk , which respectively correspond to the weight functions {1, v i , mc i c j , 1 2 mc 2 c i }, with c 2 = c · c. The corresponding moment system, written in a conservative form, can be easily retrieved from (12) and reads ∂n ∂t ∂ ∂t mnu i u j + P i j + ∂ ∂ x r mnu i u j u r + P i j u r + P ri u j + P r j u i + mM ri j ∂ ∂t +P ik u k u r + P rk u i u k + 1 2 P ri u 2 + q i u r + q r u i + mM rik u k + 1 2 mM rikk = mĈ ik u k + 1 2 mĈ ikk + q m E i 1 2 mnu 2 + ne + E k (mnu i u k + P ik ) + B r 2 lr j 1 2 mM i jl +3u ( j P il) + 1 2 mu i u j u l + lri q l + 2u j P jl + u l ne + 1 2 mu 2 u l .
Closure relations are needed for the extra-fluxesM ri j ,M rikk and the collision termsĈ i j ,Ĉ ikk : c r c i c j f ME dc, Here, the Maximum Entropy distribution function f ME is given by: and the Lagrangian multipliers can be obtained by inverting the constraints The resulting 13-moment system (15)-(18) is a symmetrizable hyperbolic system, coupled with a Poisson equation.

Inversion of the constraint relations
In order to invert the constraints, we consider the trace-free part of the tensor P i j , defined as P i j = P i j − 1 3 P kk δ i j and, recalling that P kk = 2ne, we write P i j = P i j − 2 3 neδ i j . We notice that P i j and ne are the moments associated to the weights mc i c i and m 1 2 c 2 . Then, the Maximum Entropy distribution function becomes withλ e = 2 3λ P kk . We introduce a distinction between isotropic moments, n, ne, and anisotropic moments, u i , P i j , q i , by assuming small anisotropic Lagrange multipliers,λ u i ,λ P i j ,λ q i . Then, at first order in anisotropy, the ME function reads with Substituting the f ME (25) into the constraint relations, it is immediate to obtain the following system 2 mλ e c 2 dc , Defining n i := c i /c and recalling the two properties where S 2 is the unit sphere of R 3 , d is the element of the solid angle and is the Gamma function; it is possible to solve the previous system with respect to the Lagrange multipliers to get These results are typical of the first-order approximation of extended thermodynamics of monatomic gases. It is well-known in the literature that this causes a reduction of the hyperbolicity region of the system, which improves if the second order approximation is used [4]. As said, the hypothesis of small anisotropy could be removed by numerically inverting the constraints which would require additional computational effort. We plan to do this in a future paper.

Closure for the extra-fluxes and the collision terms
Having analytically solved the constraints under the small anisotropy assumption, one can find the constitutive relations by inserting the approximate expression into the integral expressions of the extra variables. For the extra-fluxes, one easily obtainŝ The calculation of the production terms is more involved, due to the expression of the collision operator. We first notice that the kernel of the collision operator satisfies the following Moreover, at the first order, we have Indicating by ψ(c) a generic weight function, the corresponding production term iŝ where we have exploited the behaviour at infinity of f ME and changed the role of c and c for one half of the integral. Using Property 1 and (27), one has As a consequence the generic production term can be written as follows: This expression can be further simplified, if we distinguish between ψ(c) being an even or an odd function. In fact, indicating the even and odd functions by ψ e and ψ o , respectively, we havê In the thirteen moment model, the production terms to compute are as follows: C, C i ,Ĉ i j , andĈ ill . It is immediate to see that the first two vanish together withĈ ll , the latter as consequence of Property 1. This is consistent with the conservation of number of particles, momentum and energy in the collisions between particles. Therefore, the collision terms which remain to be computed areĈ i j andĈ ill , the first involves an even function, c i c j , the second an odd one, c i c l c l . After some algebra and making the substitution c = c−c , one findsĈ ME (c )c 3 c 2 δ k(i n j) − n i n j n k n l d d dc dc , ME (c )c 3 c 2 c + 2c n · n 2 δ ik − c 2 + 8c c n · n + 16 n · n 2 c 2 − 4c 2 n i n k +4c c + 2c n · n n (i n k) d d dc dc , where c = |c |, c = |c |. If in the integration with respect to c one uses as reference frame one which has the z-axis parallel to n , the above-written integrals can be explicitly computed to get As can be seen, after lengthy and cumbersome computations, the closure relations have very simple explicit expressions.
The final system is the following: ∂ ∂t +P ik u k u r + P rk u k u i + 1 2 P ri u 2 + 7 5 (q i u r + q r u i ) + 2 5 q k u k δ ri Equation (32) can be replaced by the equations ∂ ∂t where

A numerical test
From now on, time is scaled to the inverse plasma frequency ω p , lengths to the Debye length λ D and velocities to the thermal speed v th , the electron density n = 1, and all other physical quantities will be scaled using these characteristic parameters. In this section, we exploit a numerical test to check to what extent the small anisotropy approximation gives good results, by a comparison with the results obtained by the kinetic model [8], in the case of the relaxation towards equilibrium of a homogeneous field-free plasma with a temperature anisotropy. Therefore, at a microscopic level, the electron state can be described by a bi-Maxwellian distribution function: Here, the subscript || indicates the z direction, while x and y are the perpendicular (⊥) directions. Let us define the temperature anisotropy, at time t, as A(t) = T ⊥ (t)/T || (t), with A 0 = A(0). By assuming that the distribution function remains a bi-Maxwellian during the process of collisional relaxation, and recalling that the total temperature T = 2T ⊥ +T || 3 remains constant in time, the evolution of the distribution is determined by that of T ⊥ .
The evolution equation of T ⊥ was found by Kogan [8] and reads ν K being a thermalization frequency given by: whereÃ := A − 1 := T ⊥ /T || − 1 and The solution of (37) can be determined numerically.
To obtain the corresponding evolution equation for T ⊥ given by the model presented in this paper, first of all we have to expand the bi-Maxwellian in small anisotropy, which gives where 2T || +T ⊥ T || T ⊥ is the trace of the inverse temperature tensor, and the expansion is made with respect to otherwise. After that, comparing this expansion with the MEP distribution (26), taken with q = 0 and P <11> = P <22> , we find Now, we have to write the moments system in the case where no spatial dependence and no fields are present.
The resulting system is as follows: ∂ ∂t ∂ ∂t We get that n is constant, as well as u and e. If the velocity and the heat flux are posed equal to zero, we remain with the equation: where ν ME = 3 5 nU 0 6 πe 3 1 2 .
Therefore, by solving this equation and using (38), we find In Figs. 1 and 2 and in the left side of Fig. 3, we compare the results for the two models, in the cases A 0 = 0.03, 0.9, 1.5, 15, 30. We can see that, for A 0 > 1, the results of the macroscopic (ME) model are in excellent agreement with the kinetic ones, with a maximum relative error which is equal to around 1% for A 0 = 30, see Fig. 4. Worse are the results for A 0 < 1 with a maximum relative error which rapidly increases from around 3.4% for A 0 = 0.5 to around 47% for A 0 = 0.03, see Fig. 4.
This different behaviour can be explained by looking at the trend of the right-hand side of dT ⊥ dt as a function of T ⊥ for the two models, which is shown in the left side of Fig. 5. As can be seen, the two production functions are much nearer for T ⊥ > T . It can be noticed that the ME model works well beyond the limit of small anisotropy, which, in the present case, is related to the smallness of the trace-free part of the inverse temperature tensor which is proportional to T ⊥ −T || T ⊥ T || . In fact, as can be seen in the right side of Fig. 5, the neighbourhood of T ⊥ = T where |T ⊥ −T || | T ⊥ T || is small is much narrower than that where the results of the ME model are in good agreement with those of the kinetic model.

Conclusions
In this paper, we have described a general procedure to obtain a physics-based closure to the moment equations derived from the FPL transport equation for plasmas, with an arbitrary number of moments. In particular, we have explicitly written the 13-moment model for the intrinsic moments and we have tested it in the case of a homogeneous field-free plasma with an anisotropic temperature, obtaining results which are in good agreement with those given by the FPL equation. This agreement gives hints that the proposed macroscopic model, which requires a significantly lighter computational effort with respect to the kinetic model, makes self-consistent simulations of plasmas in presence of collisions affordable, even in the three-dimensional space geometry. This will be the object of a future paper.