Exact collisional moments for plasma fluid theories

TL;DR

This paper derives exact expressions for collisional moments in plasma fluid theories using Hermite-polynomial moments, enabling precise modeling of plasma collisions without restrictions on species properties.

Contribution

It provides a novel exact formula for collisional moments in plasma physics, applicable to both Chapman-Enskog and Grad's moment methods, with no restrictions on plasma parameters.

Findings

Exact collisional moments expressed via Hermite-polynomial moments.

Derived collisional ten-moment equations with conservation laws.

Applicable to various plasma transport and fluid models.

Abstract

The velocity-space moments of the often troublesome nonlinear Landau collision operator are expressed exactly in terms of multi-index Hermite-polynomial moments of the distribution functions. The collisional moments are shown to be generated by derivatives of two well-known functions, namely the Rosenbluth-MacDonald-Judd-Trubnikov potentials for a Gaussian distribution. The resulting formula has a nonlinear dependency on the relative mean flow of the colliding species normalised to the root-mean-square of the corresponding thermal velocities, and a bilinear dependency on densities and higher-order velocity moments of the distribution functions, with no restriction on temperature, flow or mass ratio of the species. The result can be applied to both the classic transport theory of plasmas, that relies on the Chapman-Enskog method, as well as to deriving collisional fluid equations that follow Grad's moment approach. As an illustrative example, we provide the collisional ten-moment equations with exact conservation laws for momentum- and energy-transfer rate.

Equations

$$\frac{d{f}_s}{dt}=\sum _{s^{\prime }}{C}_{s{s}^{\prime }}[f_s,f_{s^{\prime }}].$$

$${\bar{\mathbf{G}}}_{(k)}(\mathbf{x}-\mathbf{\mu };{\sigma }^2)=\frac{1}{{\mathcal{N}}_{{\sigma }^2}(\mathbf{x}-\mathbf{\mu })}(-{\nabla }_{\mathbf{x}}{)}^{(k)}{\mathcal{N}}_{{\sigma }^2}(\mathbf{x}-\mathbf{\mu }),$$

$${\bar{\mathbf{H}}}_{(k)}(\mathbf{x}-\mathbf{\mu };{\sigma }^2)=\frac{1}{{\mathcal{N}}_{{\sigma }^2}(\mathbf{x}-\mathbf{\mu })}(-{\sigma }^2{\nabla }_{\mathbf{x}}{)}^{(k)}{\mathcal{N}}_{{\sigma }^2}(\mathbf{x}-\mathbf{\mu })={\sigma }^{2k}{\bar{\mathbf{G}}}_{(k)}(\mathbf{x}-\mathbf{\mu };{\sigma }^2),$$

$${\mathcal{N}}_{{\sigma }^2}(\mathbf{x}-\mathbf{\mu })\equiv \frac{e^{-(\mathbf{x}-\mathbf{\mu }{)}^2/2{\sigma }^2}}{(2\pi {)}^{3/2}{\sigma }^3}.$$

$${\nabla }^{(k)}\equiv \underset{k\text{ terms}}{\underbrace{\nabla \otimes \cdots \otimes \nabla }},{\mathbf{x}}^{(k)}\equiv \underset{k\text{ terms}}{\underbrace{\mathbf{x}\otimes \cdots \otimes \mathbf{x}}},$$

$${\bar{\mathbf{G}}}_{(k)}(\mathbf{x})={\bar{\mathbf{G}}}_{(k)}^{k_1\dots k_k}(\mathbf{x})\text{where}{k}_i\in \{1,2,3\}.$$

$$\displaystyle\int_{\mathbb{R}^{3}}d\bm{y}\bm{\bar{H}}_{(i)}(\bm{y};\sigma^{2})\bm{\bar{G}}_{(j)}(\bm{y};\sigma^{2})\mathcal{N}_{\sigma^{2}}(\bm{y})\overset{(\ref{eq:grad_gaussian},\ref{eq:weierstrass_hermite})}{=}\nabla_{\bm{x}}^{(j)}\bm{x}^{(i)}\Big{|}_{\bm{x}=\bm{0}}=\delta^{(i)}_{[(j)]},$$

$$\frac{f_s(\mathbf{v})}{n_s}={\mathcal{N}}_{{\sigma }_s^2}(\mathbf{v}-{\mathbf{V}}_s)\sum _{i=0}^{\infty }\frac{1}{i!}{\mathbf{c}}_{s(i)}{\bar{\mathbf{G}}}_{(i)}(\mathbf{v}-{\mathbf{V}}_s;{\sigma }_s^2)\stackrel{(\text{ref}eq:covarian{t}_{g}en)}{=}\sum _{i=0}^{\infty }\frac{{\mathbf{c}}_{s(i)}}{i!}{\nabla }_{{\mathbf{V}}_s}^{(i)}{\mathcal{N}}_{{\sigma }_s^2}(\mathbf{v}-{\mathbf{V}}_s),$$

$${\mathbf{c}}_{s(j)}\equiv {\int }_{{\mathbb{R}}^3}d\mathbf{v}\frac{f_s(\mathbf{v})}{n_s}{\bar{\mathbf{H}}}_{(j)}(\mathbf{v}-{\mathbf{V}}_s;{\sigma }_s^2).$$

$$\mathbf{J}\cdot {\nabla }_{\mathbf{v}}{\bar{\mathbf{H}}}_{(k)}(\mathbf{v}-\mathbf{V};{\sigma }^2)=\frac{1}{(k-1)!}{J}^{[k_1}{\bar{\mathbf{H}}}_{(k-1)}^{k_2\cdots k_k]}(\mathbf{v}-\mathbf{V};{\sigma }^2)=k\text{ Sym}[\mathbf{J}{\bar{\mathbf{H}}}_{(k-1)}(\mathbf{v}-\mathbf{V};{\sigma }^2)]$$

$$C_{s{s}^{\prime }}[f_s,f_{s^{\prime }}](\mathbf{v})\equiv -\frac{c_{s{s}^{\prime }}}{m_s}{\nabla }_{\mathbf{v}}\cdot {\mathbf{J}}_{s{s}^{\prime }}[f_s,f_{s^{\prime }}](\mathbf{v}).$$

$${\mathbf{J}}_{s{s}^{\prime }}[f_s,f_{s^{\prime }}](\mathbf{v})\equiv {\mu }_{s{s}^{\prime }}({\nabla }_{\mathbf{v}}{\varphi }_{s^{\prime }})f_s-m_s^{-1}{\nabla }_{\mathbf{v}}\cdot \left[({\nabla }_{\mathbf{v}}{\nabla }_{\mathbf{v}}{\psi }_{s^{\prime }})f_s\right],$$

$${\varphi }_s(\mathbf{v})$$

$${\mathbf{C}}_{s{s}^{\prime }(k+1)}$$

$$\equiv {\bar{c}}_{s{s}^{\prime }}(k+1)\text{Sym}\left[{\mu }_{s{s}^{\prime }}{\mathbf{R}}_{s{s}^{\prime }(k+1)}+\frac{k}{m_s}{\mathbf{D}}_{s{s}^{\prime }(k+1)}\right],$$

$${\mathbf{R}}_{s{s}^{\prime }(k+1)}$$

$${\mathbf{D}}_{s{s}^{\prime }(k+1)}$$

$$\frac{1}{n_{s^{\prime }}}{\nabla }_{\mathbf{v}}{\varphi }_{s^{\prime }}(\mathbf{v})$$

$$\frac{1}{n_{s^{\prime }}}{\nabla }_{\mathbf{v}}{\nabla }_{\mathbf{v}}{\psi }_{s^{\prime }}(\mathbf{v})$$

$$\frac{1}{n_s}{f}_s(\mathbf{v}){\bar{\mathbf{H}}}_{(k)}(\mathbf{v}-{\mathbf{V}}_s;{\sigma }_s^2)$$

$${\bar{\mathbf{a}}}_{(i)(j)}^{(l)}$$

$${\iint }_{{\mathbb{R}}^3}d\mathbf{v}d{\mathbf{v}}^{\prime }\frac{{\mathcal{N}}_{{\sigma }_{s^{\prime }}^2}({\mathbf{v}}^{\prime }-{\mathbf{V}}_{s^{\prime }}){\mathcal{N}}_{{\sigma }_s^2}(\mathbf{v}-{\mathbf{V}}_s)}{|\mathbf{v}-{\mathbf{v}}^{\prime }|}$$

$$\frac{1}{2}{\iint }_{{\mathbb{R}}^3}d\mathbf{v}d{\mathbf{v}}^{\prime }{\mathcal{N}}_{{\sigma }_{s^{\prime }}^2}({\mathbf{v}}^{\prime }-{\mathbf{V}}_{s^{\prime }}){\mathcal{N}}_{{\sigma }_s^2}(\mathbf{v}-{\mathbf{V}}_s)|\mathbf{v}-{\mathbf{v}}^{\prime }|$$

$$\Phi (z)$$

$$\Psi (z)$$

$${\mathbf{R}}_{s{s}^{\prime }(k+1)}$$

$${\mathbf{D}}_{s{s}^{\prime }(k+1)}$$

$$\Phi ({\Delta }_{ei})$$

$$\chi (\mathbf{w})=\sum _{n=2}^{\infty }{a}_n{L}_n^{1/2}\left(\frac{w^2}{2{\sigma }^2}\right)+\mathbf{w}\cdot \sum _{n=1}^{\infty }{\mathbf{b}}_n{L}_n^{3/2}\left(\frac{w^2}{2{\sigma }^2}\right)+\left(\mathbf{w}\mathbf{w}-w^2\frac{\mathbf{I}}{3}\right):\sum _{n=0}^{\infty }{\mathbf{d}}_n{L}_n^{5/2}\left(\frac{w^2}{2{\sigma }^2}\right),$$

$$\chi (\mathbf{w})=$$

$${\mathbf{c}}_{(2n)}=(-1{)}^n\frac{(2n)!}{2^{n}n!}{\sigma }^{2n}\text{Sym}\left[a_n{\mathbf{\delta }}_{2n}-2n{\sigma }^2({\mathbf{d}}_{n-1}{\mathbf{\delta }}_{2n-2}-\frac{1}{3}\text{tr}[{\mathbf{d}}_{n-1}]{\mathbf{\delta }}_{2n})\right]$$

$${\mathbf{c}}_{(2n+1)}=(-1{)}^n\frac{(2n+1)!}{2^{n}n!}{\sigma }^{2n+2}\text{Sym}\left[{\mathbf{b}}_n{\mathbf{\delta }}_{2n}\right].$$

$$f_s(\mathbf{v})$$

$$=n_s{\left(\frac{m_s}{2\pi T_s}\right)}^{3/2}{e}^{-\frac{m_s}{2T_s}(\mathbf{v}-{\mathbf{V}}_s{)}^2}\left[1+\frac{m_s}{2T_s}(\mathbf{v}-{\mathbf{V}}_s)\cdot \left(\frac{{\mathbf{P}}_s}{p_s}-\mathbf{I}\right)\cdot (\mathbf{v}-{\mathbf{V}}_s)\right].$$

Full text

collision operator; Hermite polynomials; fluid moments

I Introduction

Fluid models have been widely employed in many fields of science, ranging from astronomy and physics to biology and chemistry. The fundamental principle, and motivation, behind fluid models is to provide an effective macroscopic representation of the collective behaviour arising from a large number of microscopic events. Thus, the main advantage of fluid models is a reduction in complexity, while still capturing the essential characteristics of the macroscopic system.

Typically, the fluid equations are derived from a parent kinetic model Ichimaru (1973); Villani (2002) in which the particle dynamics is governed by the equation

[TABLE]

In the above formula, $f_{s}(t,\bm{x},\bm{v})$ denotes the phase-space distribution function of species $s$, and $d/dt=\partial/\partial t+\dot{\bm{x}}\cdot\nabla_{\bm{x}}+\dot{\bm{v}}\cdot\nabla_{\bm{v}}$ is the free-streaming Vlasov operator. On the RHS, $C_{ss^{\prime}}[f_{s},f_{s^{\prime}}]$ is the (bilinear) collision operator between the particle species $s$ and $s^{\prime}$, thereby embodying the transition from many-body dynamics to the dynamical evolution of a single-particle distribution function.

The difficulty in constructing fluid models from kinetic theory arises from the presence of the collision operator on the RHS of (1). There exist two primary approaches for addressing this issue: 1) the Chapman-Enskog procedure Chapman and Cowling (1970) is a perturbation theory that relies upon a small-parameter expansion in the Knudsen number of the kinetic equation, and 2) Grad’s procedure Grad (1949a) is a Galerkin projection based on the expansion of the distribution functions in terms of orthogonal polynomials. The starting point in both cases is a Maxwellian distribution function, corresponding to the null space of the collision operator, while further refinements require evaluation of velocity moments of the collision operator.

Recently, the fluid moments of the nonlinear Landau collision operator were provided in a systematic and programmable way Hirvijoki et al. (2016), and were demonstrated to be generated by the gradients of three scalar valued integrals. However, closed-form expressions for the integrals were not found. In the present work, we improve upon previous findings Hirvijoki et al. (2016), and provide this time a closed, analytic form. The derivation exploits some remarkable properties of the Hermite polynomials and the Maxwellian distribution function under convolution. The result is expressed as gradients of two well-known functions, namely the Rosenbluth-MacDonald-Judd-Trubnikov potentials (Rosenbluth et al., 1957; Trubnikov, 1958) for a Gaussian distribution, taken with respect to a dimensionless variable denoting the relative mean flow of the colliding species normalised to root-mean-square of the corresponding thermal velocities. By this fact, the result is manifestly Galilean invariant, which is a property that is usually not preserved under various approximations of the Landau collision operator. The formula is also bilinear with respect to species densities and the so-called second- or higher-order Hermite moments of the distribution functions. The procedure is valid regardless of the mass ratio, temperature or flow difference between species and, since the Hermite polynomials form a complete basis, it is exact. In other words, the knowledge of the Hermite moments of any distribution function is sufficient to provide the collisional moments of the nonlinear Landau operator exactly.

Our representation is most convenient for the hierarchy of moment equations obtained via Grad’s approach. Generating extended collisional fluid equations for plasmas to arbitrary order is then expected to be straight-forward with the help of computer algebra systems. Given the equivalence between Laguerre and contracted multi-index Hermite polynomials – discussed in detail in appendix C – a linearised version of our general formula can also be used within a Chapman-Enskog approach to recover Braginskii’s transport coefficients, and possibly to extend the calculation to arbitrary species, flows and temperature differences.

The paper begins in Sec. II with a thorough discussion on the similarities and differences between the Chapman-Enskog and Grad’s approach. In particular, it is motivated why truncated distribution functions are required in both cases. Sec. III recaps Grad’s Hermite expansion for square-integrable functions and provides important definitions and identities for further computations. After presenting the necessary tools, the collisional moments are given explicit expressions in Sec. IV. To consider how the result could be applied to the Chapman-Enskog theory, a correspondence between the Laguerre expansion – typically used to solve the so-called correction equations – and Grad’s Hermite expansion of the distribution function is established in Sec. V. To illustrate how the resulting formulae can be applied in the context of Grad’s moment approach, the collisional ten-moment fluid equations are derived in Sec. VI and the nonlinear expressions for the momentum- and energy-transfer rate are proven to exactly satisfy the conservation laws. Sec. VII concludes the work.

II Approaches to fluid theories

Chapman-Enskog:

With the Chapman-Enskog approach, near-Maxwellian corrections are derived through a hierarchy of linear asymptotic equations, order by order. At each order, an integral equation must be solved that involves the Vlasov operator acting on the distribution function of the previous order on the LHS and the collision operator acting linearly on the current order on the RHS. Spatial gradients of the Maxwellian variables (temperature, density and flow) thus become associated with the collisional moments of higher-order corrections to the distribution function. The parameters in front of the spatial gradients are collected as transport coefficients; for example at first order, the stress tensor is proportional to the strain tensor (or velocity gradient) via viscosity and the heat flux to temperature gradient via conductivity. For neutral fluids, the Euler equations, which originate from purely Maxwellian behaviour in velocity space, become at first-order the Navier-Stokes equations and at higher-order the Burnett and super-Burnett equations.

For magnetised plasmas, the Chapman-Enskog procedure is a more complex multi-parameter perturbation theory due to the dominance of the Lorentz force, the long-range Coulomb interactions and the presence of multiple species. The first-order solution is known as the Braginskii equations Braginskii (1965) and is commonly used to describe (classical) transport in magnetised plasmas where the mean-free-path is comparable to the Larmor radius but shorter than characteristic gradient scale lengths. Different orderings and second-order solutions have been attempted Mikhailovskii and Tsypin (1971); Catto and Simakov (2004). In order to make the problem tractable, the Landau collision operator Landau (1937) between different species must be approximated and/or split into simpler forms, especially in the case where the species flows and temperatures are different Ramos (2007). In this regard, the Maxwellian solution and thus the starting point for the perturbative treatment of the kinetic equation is valid only in the limit of vanishing mass ratio between electron and ion species, in which case the electron-ion collision operator is approximated by the Lorentz collision operator.

Solving the integral equation at the next order is far from trivial. By virtue of the self-adjoint properties of the linearised collision operator, it is typically reformulated as a variational problem Robinson and Bernstein (1962), where a functional of the solution is maximised by varying the coefficients of a polynomial expansion of the distribution function. Depending on the choice of polynomial basis (Sonine, Laguerre or Legendre), the transport coefficients obtained converge rapidly to their physical values as a function of the number of terms in the series Hinton and Hazeltine (1976). The Chapman-Enskog procedure, applied to plasmas or neutral fluids, is thus formally solved using a truncated polynomial expansion of the distribution function.

Grad:

Grad’s method relies upon solving the kinetic equation indirectly by projecting it onto a set of orthogonal polynomials, namely the multi-index Hermite polynomials. The Hermite basis naturally arises from the Gaussian measure, i.e. the null-space solution of the collision operator. The procedure yields a weak solution to the kinetic equation, the same way finite-element methods Brenner and Scott (2007) are constructed, or observables are obtained as expectation values of the Schrödinger equation in quantum mechanics (Weinberg, 2012). The result is an infinite hierarchy of dynamical equations for the projection coefficients. The projection coefficients correspond to moments of the distribution function and are evolved in time by these equations from known initial conditions. The evolution of all moments is equivalent to the evolution of the distribution function, in the sense that the complete set of moment equations is a spectral representation of the kinetic equation. In its infinite form, this method remains a microscopic description of the fluid at hand (Balescu, 1988, Chapter 4.5).

A closed set of fluid equations is obtained by setting, above a given order, all expansion coefficients to zero, thereby restricting the solution of the distribution function to a finite dimensional vector space. The justification for a specific truncation depends on the studied flow and is done a posteriori by verifying certain realisability conditions (Levermore, 1996, Eq.(3.24)). We note that these realisability conditions also apply to the Chapman-Enskog approach. Orthogonality among Hermite polynomials guarantees that all the fluid contributions of a given order have been properly captured. The projection of the collision operator does not yield spatial gradients, and therefore the collisional moments are not immediately related to orthodox transport coefficients, as in the Chapman-Enskog approach.

In the limit of vanishing mean-free-path, one can apply a Chapman-Enskog analysis to Grad’s hierarchy of moment equations and show under certain assumptions on timescales and flow that they faithfully reduce to the Navier-Stokes equations Grad (1949a). The effect of the truncation is mainly to underestimate the coefficient of viscosity, as proven for generic collision operators in neutral fluids by Levermore (Levermore, 1996, Eq.(5.30)). Levermore also shows that there is a natural way of ordering corrections to transport coefficients involving higher-order moments; a BGK analysis of the collision operator in the limit of zero mean-free-path reveals that the associated relaxation rates are increasingly stronger (Levermore, 1996, Eq.(6.28)). Grad’s moment equations, although limited in the accuracy of asymptotic transport coefficients, have the advantage over Chapman-Enskog transport equations of retaining more kinetic features and being valid beyond Maxwellian behaviour, as discussed in (Grad, 1949a, Appendix 4), thereby embedding classical transport theory as a special case of its realisable flows. Grad’s method is thus often applied to the investigation of the nonlinear properties of fluids and the formation of shocks in rarefied gases.

In plasma physics, Grad’s moment equations are not used nearly as much as Braginskii’s, even though the trial solution to the first order Chapman-Enskog correction equation written with Laguerre polynomials corresponds identically to a contracted multi-index Hermite polynomial expansion. This means that the trial solution with only the first ($N=1$) Laguerre polynomial coincides with Grad’s 13-moment expansion and the trial solution with $N=2$ truncation matches Grad’s 21-moment and so on. This correspondence was noticed, e.g., by Balescu (Balescu, 1988, chapter 4) and is revisited in section V for the expansion coefficients in terms of which the collisional moments of the Landau operator are expressed.

III Hermite expansion of square-integrable functions

Inspired by the seminal work by Grad on the asymptotic theory of the Boltzmann equation Grad (1949a), we consider the spectral expansion of square-integrable functions in terms of multi-index Hermite polynomials. While various definitions exist, we will employ, from the review by Holmquist Holmquist (1996), the so-called covariant Hermite polynomials

[TABLE]

as well as the so-called contravariant Hermite polynomials

[TABLE]

both generated by the three-dimensional Normal distribution

[TABLE]

Regarding our notation for outer products of vectors (such as consecutive application of the gradient operator), it is implied throughout this document that

[TABLE]

whereas for tensors (such as the Hermite polynomials), the notation refers to the rank of the multi-indexing according to

[TABLE]

One important property of the Hermite polynomials is that they are orthogonal to each other with respect to the Gaussian measure,

[TABLE]

where $[(j)]=[j_{1}\cdots j_{j}]=\sum_{\pi}\pi(j_{1})\cdots\pi(j_{j})$ is the sum over all permutations of $j$ indices and $\delta^{(i)}_{(j)}=\delta^{i_{1}}_{j_{1}}\dots\delta^{i_{i}}_{j_{i}}$. This orthogonality property allows for a convenient expansion of any square-integrable function, such as the distribution function $f_{s}$ for species $s$, according to

[TABLE]

where the (symmetric) expansion coefficients are the so-called Hermite-moments of the distribution function

[TABLE]

As per the Einstein summation convention on repeated indices $(i)$, the tensors $\bm{\bar{G}}_{(i)}(\bm{x};\sigma^{2})$ in the expansion are fully contracted with the tensors $\bm{c}_{(i)}$.

Written in the form (8), the distribution function of each species is automatically normalised to the species’ density and thus $c_{s(0)}=1$. The mean velocity, $\bm{V}_{s}$, and the variance of each species (half thermal velocity squared), $\sigma_{s}^{2}=\tfrac{1}{2}v_{\text{th},s}^{2}=T_{s}/m_{s}$, is contained in the Gaussian envelope, so that $\bm{c}_{s(1)}=\bm{0}$ and $\bm{c}_{s(2)}=(\bm{P}_{s}-p_{s}\bm{I})/n_{s}m_{s}=\sigma_{s}^{2}(\bm{P}_{s}/p_{s}-\bm{I})$ represents the trace-less pressure tensor measuring the degree of anisotropy and off-diagonal features, where $m_{s}$, $n_{s}$ and $T_{s}$ are the species’ mass, density and temperature respectively, $p_{s}=n_{s}T_{s}=\frac{1}{3}\text{tr}\bm{P}_{s}$ is the isotropic pressure and $\bm{P}_{s}\equiv m_{s}\int_{\mathbb{R}^{3}}d\bm{v}^{3}\;(\bm{v}-\bm{V}_{s})(\bm{v}-\bm{V}_{s})f_{s}$ is the pressure tensor. Third-order tensorial moments of the distribution function are captured by the tensor $\bm{c}_{s(3)}=\int_{\mathbb{R}^{3}}d\bm{v}^{3}\;(\bm{v}-\bm{V}_{s})(\bm{v}-\bm{V}_{s})(\bm{v}-\bm{V}_{s})f_{s}/n_{s}$ which, when maximally contracted, represents the heat-flux, $q^{i}_{s}=m_{s}n_{s}c_{s(3)}^{ikk}=m_{s}\int_{\mathbb{R}^{3}}d\bm{v}^{3}\;|\bm{v}-\bm{V}_{s}|^{2}(\bm{v}-\bm{V}_{s})f_{s}$.

The temporal and spatial dependence of the kinetic distribution $f_{s}(\bm{v})$ has been omitted for convenience, but it is understood that the coefficients $n_{s}$, $\bm{V}_{s}$, $T_{s}$ and $\bm{c}_{s(i)}$ (i.e. the fluid variables) vary with respect to time and spatial coordinates. It is also important to note that the definition for the projection coefficients (9) is independent of how the distribution function is presented.

Before computing the Hermite-moments of the collision operator, one more identity is needed, namely the directional derivative of a given Hermite polynomial. For any vector $\bm{J}$, one has (Holmquist, 1996, Eq.(6.2))

[TABLE]

where $\text{Sym}\bm{A}^{a_{1}\cdots a_{k}}=\bm{A}^{[a_{1}\cdots a_{k}]}/k!$ is the symmetrization of a tensor. This result appears in Grad’s original work for dimensionless Hermite polynomials (Grad, 1949b, Eq.(17)).

IV Landau Collision operator and Hermite-moments

In warm plasmas, collisions are dominated by continuous small-angle Coulomb scattering. The appropriate operator to describe the collective effect of these events was derived by Landau Landau (1937), and can be expressed as a velocity-space divergence of a collisional velocity-space flux defined by

[TABLE]

Here, $c_{ss^{\prime}}=\ln\Lambda(e_{s}e_{s^{\prime}})^{2}/4\pi\varepsilon_{0}^{2}$, $\ln\Lambda$ denotes the Coulomb Logarithm, and $e_{s}$ is the species charge. The collisional velocity-space flux $\bm{J}_{ss^{\prime}}$ can be represented in its original integral form, or through the so-called Rosenbluth-MacDonald-Judd-Trubnikov potential functions Trubnikov (1958); Rosenbluth et al. (1957) according to

[TABLE]

where $\mu_{ss^{\prime}}=1/m_{s}+1/m_{s^{\prime}}$ and the potential functions, $\phi_{s}(\bm{v})$, and $\psi_{s}(\bm{v})$, are defined through

[TABLE]

We observe that $\nabla_{\bm{v}}\cdot\nabla_{\bm{v}}\psi_{s}=\phi_{s}$ and $\nabla_{\bm{v}}\cdot\nabla_{\bm{v}}\phi_{s}=-4\pi f_{s}$.

To derive collisional fluid equations based on Grad’s expansion (8), we consider the Hermite-moments of the collision operator

[TABLE]

where $\bar{c}_{ss^{\prime}}=n_{s}n_{s^{\prime}}c_{ss^{\prime}}$ and the integral has been split, after first integrating by parts and then using identity (10), into drag- and diffusion-related terms

[TABLE]

One may already notice that $C_{ss^{\prime}(0)}=0$ and $\bm{D}_{ss^{\prime}(1)}=\bm{0}$ such that $\bm{C}_{ss^{\prime}(1)}=\bar{c}_{ss^{\prime}}\mu_{ss^{\prime}}\bm{R}_{ss^{\prime}(1)}\equiv\bm{F}_{ss^{\prime}}$ represents the collisional momentum transfer rate between species $s$ and $s^{\prime}$.

To proceed, the velocity gradients of the potential functions are manipulated in order to extract derivatives with respect to the mean velocity $\bm{V}_{s^{\prime}}$ according to

[TABLE]

Next, the products of two Hermite polynomials in $f_{s}(\bm{v})\bm{\bar{H}}_{(k)}(\bm{v}-\bm{V}_{s};\sigma_{s}^{2})$ are expressed as a series of single Hermite polynomials (linearisation) using several identities derived in Appendix A so to extract derivatives with respect to the mean velocity $\bm{V}_{s}$ according to

[TABLE]

The expression for the so-called linearisation coefficient $\bm{\bar{a}}^{(l)}_{(i)(j)}$, derived explicitly in Appendix B, is

[TABLE]

Since the gradients with respect to $\bm{V}_{s}$ and $\bm{V}_{s^{\prime}}$ can be brought out of the integrals (15) and (16), one is only left with the task of evaluating the following integrals

[TABLE]

where $\bm{U}_{ss^{\prime}}=\bm{V}_{s}-\bm{V}_{s^{\prime}}$, $\Sigma_{ss^{\prime}}^{2}=\sigma_{s}^{2}+\sigma_{s^{\prime}}^{2}=\frac{1}{2}(v^{2}_{\text{th},s}+v^{2}_{\text{th},s^{\prime}})$, and the functions $\Phi(z)$ and $\Psi(z)$ are nothing but the Rosenbluth-MacDonald-Judd-Trubnikov potentials for a Normal distribution ${\cal N}_{1/2}(z)$, defined according to

[TABLE]

Intermediate steps to yield (21) and (22) rely on the fact that the convolution of two Normal distribution results in a Normal distribution with the sum of the variances, as seen from equation (62).

The final task is to define a dimensionless parameter $\bm{\Delta}_{ss^{\prime}}=-\bm{\Delta}_{s^{\prime}s}=\bm{U}_{ss^{\prime}}/\sqrt{2}\Sigma_{ss^{\prime}}$, and to transform gradients with respect to $\bm{V}_{s}$ and $\bm{V}_{s^{\prime}}$ into gradients with respect to $\bm{\Delta}_{ss^{\prime}}$. In physical terms, $\bm{\Delta}_{ss^{\prime}}$ is the ratio of the relative mean flow and (total) thermal velocities. As a result, the drag- and diffusion-related term are expressed as

[TABLE]

where $\Delta_{ss^{\prime}}=|\bm{\Delta}_{ss^{\prime}}|$. The collisional moments of the nonlinear Landau operator can thus be computed exactly, given the knowledge of the projection coefficients (9) of the species distribution functions. Since the Hermite polynomials form a complete basis, the result is independent of how the distribution function is presented. It is directly applicable to Grad’s expansion but can be accommodated to other polynomials (such as Laguerre) in the context of the Chapman-Enskog procedure, as discussed in section V.

The convergence of the bilinear series depends on the ratio between the Hermite-moments $\bm{c}_{s(i)}$ and the $i$-th power of the (total) variance $\Sigma_{ss^{\prime}}=\sqrt{v_{th,s}^{2}+v_{th,s^{\prime}}^{2}}$, multiplied by the $i$-th gradient of the special functions $\Phi$ and $\Psi$. The latter are functions of the dimensionless parameter $\Delta_{ss^{\prime}}$ which is a small parameter in most physical cases Lingam et al. (2016); for ion-electron plasmas, $\Delta_{ei}=\sqrt{m_{e}J^{2}/2ne^{2}p_{e}}\sim(d_{e}/L)\beta_{e}^{-1/2}\ll 1$ on dimensional grounds (where $d_{e}/L$ is the normalised electron skin depth and $\beta_{e}$ the electron plasma beta). A Taylor expansion of the special functions around zero,

[TABLE]

confirms that applying any number of derivatives on these alternating fast-decaying series does not give rise to singularities nor does it affect the convergence of the collisional moments. The collisional momentum transfer rate is provided at lowest order and studied in more detail in the companion paper Lingam et al. (2016).

V Chapman-Enskog compatible Hermite expansion

In order to solve the linear integral equation in the Chapman-Enskog approach (or the so-called Spitzer problem), the total distribution function is expressed as $f/n={\cal N}_{\sigma^{2}}(\bm{w})[1+\chi(\bm{w})]$, where the random velocity $\bm{w}=\bm{v}-\bm{V}$ is adopted for convenience. Considering the tensorial and vectorial invariance of the first-order correction equations, the term $\chi(\bm{w})$ is generically of the form

[TABLE]

where most notably, the first vector valued expansion coefficient corresponds to heat-flux $\bm{b}_{1}=-m\bm{q}/5pT$ and the first tensor valued expansion coefficient corresponds to viscosity $\mathbf{d}_{0}=m(\bm{P}-p\bm{I})/2pT$, and is therefore traceless $\text{tr}(\mathbf{d}_{0})=0$. The linear integral equation in the Chapman-Enskog theory is then efficiently converted into a linear algebraic equations for the scalar, vector, and tensor coefficients $a_{n}$, $\bm{b}_{n}$, and $\mathbf{d}_{n}$ Robinson and Bernstein (1962); Braginskii (1965). Using the results from Appendix C, an equivalent expression for $\chi(\bm{w})$ can be given in terms of the irreducible Hermite polynomials of scalar, vector and two-rank tensor kind as

[TABLE]

where $N_{n}=(-1)^{n}2^{n}\sigma^{2n}n!$ is the conversion factor between the Laguerre and Hermite basis. The coefficients $\bm{c}_{(i)}$ for a Hermite expansion of the distribution function compatible with the Chapman-Enskog correction equations are thus given, for $n\geq 1$, by

[TABLE]

where the coefficient $a_{1}=0$ is reminded to be zero and the short-hand notation, $\bm{\delta}_{2n}=\delta_{i_{1}i_{2}}\delta_{i_{3}i_{4}}\cdots\delta_{i_{2n-1}i_{2n}}$, for the pair-wise contraction operator is used. Essentially, the result emerges from expressing spherically based tensor objects in Cartesian coordinates (Weinberg, 2012, Chapter 4). These coefficients can be used in (25) and (26) to express the collisional moments of the Landau operator, where only the linear terms would be retained to be consistent with the Chapman-Enskog procedure.

VI Collisional ten-moment equations

The ten-moment equations represent the simplest possible projection of the kinetic equation beyond Maxwellian behaviour. In effect, they are derived by truncating the distribution functions after the second order Hermite polynomials. This way, a closed set of collisional fluid equations is obtained in which the stress tensor features as a dynamical moment on an equal footing to density, flow and temperature. In the limit of vanishing mean-free-path, only the concept of viscosity emerges from the ten-moment model. Conductivity comes in pair with heat flux, for example in Grad’s 13-moment equations and higher-moment theories (Balescu, 1988, Chapter 4). Nevertheless, the simple ten-moment model is useful to illustrate how to apply our formulae for the collisional moments of the Landau operator, as well as to demonstrate their conservation properties. We thus proceed to consider distribution functions of the form

[TABLE]

for which there are 10 variables for each species, namely the density $n_{s}$ (1 scalar), mean velocity $\bm{V}_{s}$ (3-component vector) and pressure tensor $\bm{P}_{s}$ (6-component symmetric matrix, including $p_{s}=n_{s}T_{s}$ as its trace). The extended fluid equations are then determined via the system

[TABLE]

VI.1 Moments of the Vlasov operator

Considering that the Vlasov operator for plasmas, in the absence of gravitational forces, is given by

[TABLE]

we may write the fluid equations explicitly. Defining the mass density $\rho_{s}\equiv m_{s}n_{s}$ and the momentum vector $\bm{K}_{s}\equiv\rho_{s}\bm{V}_{s}$, equation (33) corresponds to the familiar continuity equation

[TABLE]

since $C_{ss^{\prime}(0)}=0$. Defining the stress tensor $\bm{\Pi}_{s}\equiv\bm{P}_{s}+\rho_{s}\bm{V}_{s}\bm{V}_{s}$, equation (34) becomes the momentum equation

[TABLE]

and equation (35) provides the evolution equation for the stress tensor

[TABLE]

Rather than using the number density, mean velocity, and pressure tensor, the mass density, the momentum vector, and the stress tensor were respectively introduced for the sake of expressing the moment equations in a divergence form; in numerical implementations of these equations, the divergence form allows for the use of conservative discretisation methods.

VI.2 Moments of the collision operator

The collisional contributions to the ten-moment equations require the determination of the coefficients $\bm{R}_{ss^{\prime}(1)}$, $\bm{R}_{ss^{\prime}(2)}$, and $\bm{D}_{ss^{\prime}(2)}$. By virtue of (25) and (26) one readily finds

[TABLE]

where the scalar differential operator ${\cal O}_{ss^{\prime}}$ is given by

[TABLE]

and depends on the species $s$ and $s^{\prime}$ via $\bm{\Delta}_{ss^{\prime}}$ and the normalised viscosity tensor,

[TABLE]

The operator ${\cal O}_{ss^{\prime}}$ is observed to be symmetric with respect to exchanging the species’ indices. Thus, it is seen that $\bm{R}_{ss^{\prime}(1)}$ is antisymmetric while $\bm{D}_{ss^{\prime}(2)}$ is symmetric under this operation. This property remains true when higher-moments are included.

VI.3 Conservation laws of the collisional moments

The Landau collision operator conserves particle densities, total kinetic momentum, and total kinetic energy. Since our procedure to compute the collisional moments is exact, all fluid equations derived by applying (25) and (26) automatically satisfy the same conservation properties, regardless of the order of truncation. This statement is proven explicitly for the ten-moment equations, although generalizations to higher moment fluid theories are straightforward.

The conservation of particle densities is trivial because $C_{ss^{\prime}(0)}=0$. The collisional momentum-transfer rate, given by $\bm{F}_{ss^{\prime}}=\bar{c}_{ss^{\prime}}\mu_{ss^{\prime}}\bm{R}_{ss^{\prime}(1)}$, is anti-symmetric with respect to changing the species indices, i.e., $\bm{F}_{ss^{\prime}}=-\bm{F}_{ss^{\prime}}$. This follows from the symmetry of the operator ${\cal O}_{ss^{\prime}}$ and the anti-symmetry of $\nabla_{\bm{\Delta}_{ss^{\prime}}}=-\nabla_{\bm{\Delta}_{s^{\prime}s}}$, thereby establishing the conservation of total kinetic momentum.

The collisional energy-exchange rate to species $s$ from species $s^{\prime}$ is defined as

[TABLE]

The total energy-exchange rate can be expressed, thanks to the symmetry of $\bm{D}_{ss^{\prime}(2)}$ and the anti-symmetry of $\bm{R}_{ss^{\prime}(1)}$ with respect to species indices, as

[TABLE]

It is observed, after applying the identity $\nabla^{(n)}(\bm{x}\cdot\nabla)=(\bm{x}\cdot\nabla)\nabla^{(n)}+n\nabla^{(n)}$, that

[TABLE]

and thus the total energy-exchange rate vanishes identically by virtue of

[TABLE]

The last step depends on the properties of the functions $\Phi$ and $\Psi$, namely $\nabla\cdot\nabla\Psi=\Phi$ and $\nabla\cdot\nabla\Phi=-4\pi\mathcal{N}_{1/2}(\bm{x})=-2\text{erf}^{\prime}$ and $\bm{x}\cdot\nabla\Phi=x\Phi^{\prime}=\text{erf}^{\prime}-\Phi$.

VII Conclusion

The formalism originally introduced by Grad to derive collisional fluid theories has been applied to Coulomb interactions in warm plasmas, with the sole (implicit) assumption that the Landau collision operator is valid. As our main result, analytic expressions for the collisional moments were obtained, that are represented fully in terms of derivatives of the two Rosenbluth-MacDonald-Judd-Trubnikov potentials for a Normal distribution with respect to the relative mean flow normalised to the root-mean-square of the thermal velocities. The formulae are manifestly Galilean invariant and highlight the nonlinear dependency of the collisional moments on local equilibrium Maxwellian variables, as well as the bilinearity with respect to higher-order moments of the distribution functions. Thanks to the correspondence between Laguerre and irreducible Hermite polynomials, a linearised version of our collisional moments can be used within the Chapman-Enskog approach to derive transport equations for plasmas and extend their validity to broader classes of flows, temperature differences and arbitrary species mass ratio. Within Grad’s approach, one can truncate the expansion of the distribution function at a given polynomial degree, and arrive at a fully closed, self-consistent system of extended fluid equations. As an example, the ten-moment equations were presented, and the collisional momentum- and energy-transfer rate were demonstrated to preserve exact conservation properties, which is an often overlooked requirement.

Employing our formulae to include collisional effects in fluid codes will not only improve the physical accuracy of the modelling but will also help alleviate numerical difficulties such as the build-up of sharp gradients and formation of fine structures in a consistent and controllable way. Another important consequence of our methodology is that one can provide an exact dependency of the effective electrical resistivity on all available fluid variables, not only on density and temperature. The physical consequences of such resistivity tensor, and its relevance in space and astrophysical plasmas (and fast magnetic reconnection) is explored in the companion paper Lingam et al. (2016).

Expressions for higher-order moments based on our formalism, although more complex, are naturally programmable. The formulation could possibly be widened to encompass also anisotropic distribution functions and other non-Maxwellian features by generalising to multivariate Gaussians and Hermite polynomials. Given the similarities between the Coulomb and gravitational forces Binney and Tremaine (1987), we suggest that our methodology could also be applied to the latter. Lastly, given the importance of the Hermite polynomials in the quantum harmonic oscillator (Weinberg, 2012), it is plausible that some of the mathematical identities derived here could be utilised to study this system.

Acknowledgements.

The authors would like to acknowledge useful discussions with A.Bhattacharjee, J.A.Krommes and G.W.Hammett. The insightful comments from the anonymous referee were greatly appreciated. The authors were supported by the Department of Energy Contract No. DE-AC02-09CH11466 during the course of this work.

Appendix A Hermite polynomials: definitions and properties

Many properties of the physicists’ Hermite polynomials, $\bm{H}_{(k)}(\bm{x})=\bm{\bar{G}}_{(k)}(\bm{x};\frac{1}{2})$ (in our notation), were originally derived by Grad Grad (1949b) and successfully used in his work on the collisional moments of hard spheres Grad (1949a). The Landau collision operator in plasmas being a much more complicated convolution, our previous work Hirvijoki et al. (2016) relied on both the physicists’ and probabilists’ Hermite polynomials, $\bm{He}_{(k)}(\bm{x})=\bm{H}_{(k)}(\bm{x};1)=\bm{\bar{G}}_{(k)}(\bm{x};1)$, in order to transfer the Gaussian convolution onto the Rosenbluth potentials. The collisional moments were then generated by the gradients of three scalar valued integrals, which proved to be quite tedious to evaluate. The general approach of this paper is based on a wider class of Hermite polynomials described by Holmquist Holmquist (1996) and fully exploits their properties under convolution. In this section, some fundamental identities from Holmquist (1996) are reviewed and extended.

Contravariant Hermite polynomials:

Owing to their definition in equation (3) or (Holmquist, 1996, eq.2.1), the contravariant Hermite polynomials can be obtained equivalently from

[TABLE]

The first few Hermite polynomials are listed

[TABLE]

The contravariant Hermite polynomials scale into each other upon multiplication by a scalar $c$ (Holmquist, 1996, eq.3.7) as

[TABLE]

Covariant polynomials:

The covariant Hermite polynomials of equation (2) or (Holmquist, 1996, eq.3.8),

[TABLE]

are conveniently used to express gradients of the Normal distribution with respect to its mean as

[TABLE]

The first few covariant Hermite polynomials are listed

[TABLE]

The scaling properties of covariant Hermite polynomials are opposite to equation (51)

[TABLE]

Gauss-Weierstrass transform:

Probably the most useful property of Hermite polynomials is that their Gaussian mean (or Gaussian convolution) is a Hermite polynomial with a different variance (Holmquist, 1996, eq.9.5)

[TABLE]

provided that $\sigma^{2}\geq\rho^{2}$. For the covariant Hermite polynomials, it reads

[TABLE]

The special case where $\sigma^{2}=\rho^{2}$, $\bm{x}\rightarrow\bm{y}$, $\bm{\mu}=\bm{x}$ and $\bm{\nu}=\bm{0}$ in (56) leads to the following expression for the Gauss-Weierstrass transform of a Hermite polynomial

[TABLE]

which is basically the inverse of (49), proving that the inverse Gauss-Weierstrass transform is $\mathcal{G}^{-1}_{\sigma^{2}}\equiv e^{-\frac{1}{2}\sigma^{2}\nabla^{2}}$, as far as these expansions are concerned.

Convolution identities:

It is often useful to split the product of two Gaussians in order to isolate one of the variables ($\bm{y}$ in the case that follows). For this, the quadratic exponent is arranged as

[TABLE]

where

[TABLE]

so that the product of two Gaussians becomes

[TABLE]

With this decomposition, one immediately shows that the convolution of two Gaussians is a Gaussian with the sum of the variances

[TABLE]

More generally, the Gaussian convolution of a Normal distribution times a Hermite polynomial is found to be

[TABLE]

The latter formula is particularly elegant for the case where $\rho^{2}=\tau^{2}$,

[TABLE]

This identity was used in our previous work Hirvijoki et al. (2016) to present the distribution function as a Gaussian convolution.

Linearisation:

The product of two Hermite polynomials is conveniently linearised as

[TABLE]

where the linearisation coefficient is found to be (see appendix B for a detailed proof)

[TABLE]

In particular, we evidently have

[TABLE]

We also note that for $i>0$

[TABLE]

i.e. it amounts to

[TABLE]

which can also be obtained by using the second relation in (66). The above identity comes in handy for expressing (with $i>0$)

[TABLE]

Appendix B Linearisation coefficient for products of Hermite polynomials

The linearisation coefficient for the tensor product of two Hermite polynomials is computed explicitly

[TABLE]

Appendix C Correspondence between irreducible Hermite and Laguerre polynomials

It is useful to remark that maximally contracted multi-index Hermite polynomials are proportional to a family of Laguerre polynomials. This connection was noticed by Balescu (Balescu, 1988, appendix G) and is essentially related to changing from Cartesian to spherical coordinates. A formal proof of Balescu’s construction is provided here.

Laguerre polynomials:

The generalised Laguerre polynomials are obtained via the following Rodrigues formula

[TABLE]

where the generalised binomial expansion,

[TABLE]

and the Leibniz rule, $d^{n}(fg)=\sum_{i=0}^{n}{n\choose i}(d^{i}f)(d^{n-i}g)$, are used to carry out the last step. It can be shown that they satisfy the following property under differentiation

[TABLE]

as well as the so-called three-point rules

[TABLE]

Irreducible Hermite polynomials:

A family of maximally contracted even-order Hermite polynomials, or so-called irreducible Hermite polynomials, is derived via the relation (49) by exponentiation of the Laplace operator $\nabla^{2}=\nabla\cdot\nabla$,

[TABLE]

where $x=|\bm{x}|$. Expressing the Laplacian in terms of the normalised radial variable $y=x^{2}/2\sigma^{2}$, the above definition becomes

[TABLE]

where the Laplacian operator $\Delta=\frac{1}{y^{1/2}}\frac{d}{dy}y^{3/2}\frac{d}{dy}$ is defined with respect to the variable $y$. The last equality in equation (79) follows from the fact that $\Delta$ acts as a ladder operator and decreases the powers of $y$ by one, i.e. $\Delta(y^{n})=n\left(n+\frac{1}{2}\right)y^{n-1}$, such that only terms up to $n$ matter in the sum and the order of indices can be reversed to $k=n-i$.

From (78), (79) and (72), the following correspondence is made between the even-order irreducible Hermite polynomials and the $1/2$-Laguerre basis,

[TABLE]

where the conversion coefficient, $N_{n}=(-1)^{n}2^{n}\sigma^{2n}n!$, is identified.

A family of maximally contracted odd-order Hermite vectors is derived from $h_{|2n|}$ by applying the recursion formula (49),

[TABLE]

A family of second-rank even-order Hermite tensors are generated by repeating the recursion on $\bm{h}_{|2n|+(1)}$,

[TABLE]

It is easily shown from (77) that the contraction of the expression above correctly corresponds to $h_{|2n+2|}(x;\sigma^{2})$. Hence, noting that $N_{n+1}=-2\sigma^{2}(n+1)N_{n}$, the following traceless object is formed by considering

[TABLE]

The hierarchy of higher-order tensorial Hermite polynomials extends beyond this point by successively applying the recursion formula. At each level, the traceless subparts can be identified with the corresponding Laguerre polynomial thanks to the three-point rules. A one-to-one correspondence between the Hermite and Laguerre basis is thus established, as in Balescu (Balescu, 1988, appendix G).