2 Moment-hierarchy formulation of EPW

We briefly describe the Boltzmann–Poisson system, our starting point for the description of EPW, for an unmagnetized plasma, and derive a moment expansion of the distribution function that allows its numerical solution. The Boltzmann equation for the evolution of the electron distribution function $f$ is given by

(2.1) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}t}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f+\frac{e}{m}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\boldsymbol{v}}={\hat{C}}(f). & & \displaystyle\end{eqnarray}$$

In (2.1), $e$ is the elementary charge, $m$ the electron mass, $\unicode[STIX]{x1D719}$ the electrostatic potential and ${\hat{C}}(f)$ the nonlinear Coulomb (also called Landau) collision operator for electrons

(2.2) $$\begin{eqnarray}\displaystyle {\hat{C}}(f)=\frac{v_{th}^{3}}{n}\mathop{\sum }_{b}\unicode[STIX]{x1D708}_{b}\unicode[STIX]{x2202}_{\boldsymbol{v}}\boldsymbol{\cdot }\left[\frac{m}{m_{b}}(\unicode[STIX]{x2202}_{\boldsymbol{v}}H_{b})f-\unicode[STIX]{x2202}_{\boldsymbol{v}}(\unicode[STIX]{x2202}_{\boldsymbol{v}}G_{b})\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{v}}f\right], & & \displaystyle\end{eqnarray}$$

where $v_{th}=\sqrt{2T/m}$ is the electron thermal velocity with $T$ the electron temperature, and $\unicode[STIX]{x1D708}_{b}$ the characteristic collision frequency between electrons and species $b$ ( $b=e,i$ for electrons and ions, respectively), defined by

(2.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D708}_{e}=\frac{8\sqrt{\unicode[STIX]{x03C0}}}{3\sqrt{m_{e}}}\frac{n\unicode[STIX]{x1D706}e^{4}}{T^{3/2}}, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D706}$ the Coulomb logarithm and $\unicode[STIX]{x1D708}_{i}=\sqrt{2}\unicode[STIX]{x1D708}_{e}$ . $H_{b}=2\int f_{b}(\boldsymbol{v}^{\prime })/|\boldsymbol{v}-\boldsymbol{v}^{\prime }|\,\text{d}\boldsymbol{v}^{\prime }$ and $G_{b}=\int f_{b}(\boldsymbol{v}^{\prime })|\boldsymbol{v}-\boldsymbol{v}^{\prime }|\,\text{d}\boldsymbol{v}^{\prime }$ are the Rosenbluth potentials (Helander & Sigmar Reference Helander and Sigmar2005). We relate the electrostatic potential $\unicode[STIX]{x1D719}$ to $f$ using Poisson’s equation

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}=4\unicode[STIX]{x03C0}e\left(\int f\,\text{d}\boldsymbol{v}-n_{0}\right), & & \displaystyle\end{eqnarray}$$

where the ions are assumed to provide a fixed homogeneous, neutralizing background, with density $n_{0}$ and a Maxwell–Boltzmann equilibrium with the same temperature as the electrons. An atomic number $Z=1$ is considered. Equation (2.1) is linearized by expressing $f$ as $f=f_{M}(1+\unicode[STIX]{x1D6FF}f)$ with $\unicode[STIX]{x1D6FF}f\ll 1$ and $f_{M}$ an isotropic Maxwell–Boltzmann equilibrium distribution with constant density $n_{0}$ and temperature $T_{0}$ , yielding

(2.5) $$\begin{eqnarray}\displaystyle f_{M}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FF}f}{\unicode[STIX]{x2202}t}+f_{M}\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}f+\frac{e}{m}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}\cdot \frac{\unicode[STIX]{x2202}f_{M}}{\unicode[STIX]{x2202}\boldsymbol{v}}=C(f_{M}\unicode[STIX]{x1D6FF}f), & & \displaystyle\end{eqnarray}$$

where we used the fact that ${\hat{C}}(f_{M})=0$ , and $C(f_{M}\unicode[STIX]{x1D6FF}f)$ is the linearized version of the Coulomb collision operator in (2.2) whose expression can be found in Helander & Sigmar (Reference Helander and Sigmar2005). The Boltzmann equation, equation (2.5), is coupled to the Poisson equation $\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}=4\unicode[STIX]{x03C0}e\unicode[STIX]{x1D6FF}n$ with $\unicode[STIX]{x1D6FF}n=\int f_{M}\unicode[STIX]{x1D6FF}f\,\text{d}\boldsymbol{v}$ the perturbed electron density. We now rewrite (2.5) in terms of the Fourier-transformed distribution function $\unicode[STIX]{x1D6FF}f_{k}=\int \unicode[STIX]{x1D6FF}f\exp (-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x})\,\text{d}\boldsymbol{x}$ as

(2.6) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FF}f_{k}}{\unicode[STIX]{x2202}t}+\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}\unicode[STIX]{x1D6FF}f_{k}+\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}\frac{4\unicode[STIX]{x03C0}e^{2}\unicode[STIX]{x1D6FF}n_{k}}{k^{2}T_{0}}=\frac{C(f_{M}\unicode[STIX]{x1D6FF}f_{k})}{f_{M}}, & & \displaystyle\end{eqnarray}$$

where we used the Fourier-transformed Poisson equation $-k^{2}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}_{k}=4\unicode[STIX]{x03C0}e\unicode[STIX]{x1D6FF}n_{k}$ , with $\unicode[STIX]{x1D6FF}n_{k}=\int \unicode[STIX]{x1D6FF}n\exp (-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x})\,\text{d}\boldsymbol{x}$ and $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}_{k}=\int \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}\exp (-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x})\,\text{d}\boldsymbol{x}$ .

Similarly to previous studies on the collisional damping of EPW (Brantov et al. Reference Brantov, Bychenkov and Rozmus2012; Banks et al. Reference Banks, Brunner, Berger and Tran2016), a three-dimensional cylindrical $(v_{\bot },\unicode[STIX]{x1D711},v_{z})$ velocity coordinate system is used, therefore decomposing the velocity vector $\boldsymbol{v}$ as

(2.7) $$\begin{eqnarray}\displaystyle \boldsymbol{v}=v_{z}\boldsymbol{e}_{z}+v_{\bot }(\cos \unicode[STIX]{x1D711}\boldsymbol{e}_{x}+\sin \unicode[STIX]{x1D711}\boldsymbol{e}_{y}), & & \displaystyle\end{eqnarray}$$

where $(\boldsymbol{e}_{x},\boldsymbol{e}_{y},\boldsymbol{e}_{z})$ are Cartesian unit vectors with $z$ the direction of the wave vector $\boldsymbol{k}=k\boldsymbol{e}_{z}$ . In this work, we consider the azimuthally isotropic, i.e. $\unicode[STIX]{x1D711}$ independent, subset of solutions $\unicode[STIX]{x1D6FF}f_{k}=\unicode[STIX]{x1D6FF}f_{k}(v_{\bot },v_{z})$ of the Boltzmann equation, equation (2.6), which further allows us to reduce our study to the one of a two-dimensional model. This can be done by applying the averaging operator $\langle \cdots \rangle$ defined by

(2.8) $$\begin{eqnarray}\displaystyle \langle g\rangle (v_{\bot },v_{z})=\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}g(v_{\bot },\unicode[STIX]{x1D711},v_{z})\,\text{d}\unicode[STIX]{x1D711}, & & \displaystyle\end{eqnarray}$$

to the Boltzmann equation, equation (2.5), yielding

(2.9) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\langle \unicode[STIX]{x1D6FF}f_{k}\rangle }{\unicode[STIX]{x2202}t}+\text{i}kv_{z}\langle \unicode[STIX]{x1D6FF}f_{k}\rangle +\text{i}kv_{z}\frac{4\unicode[STIX]{x03C0}e^{2}\unicode[STIX]{x1D6FF}n_{k}}{k^{2}T_{0}}=\frac{\langle C(f_{M}\unicode[STIX]{x1D6FF}f_{k})\rangle }{f_{M}}. & & \displaystyle\end{eqnarray}$$

As the linearized Coulomb collision operator satisfies $\langle C(f_{M}\unicode[STIX]{x1D6FF}f_{k})\rangle =C(f_{M}\langle \unicode[STIX]{x1D6FF}f_{k}\rangle )$ , equation (2.9) can be used to obtain the subset of azimuthally symmetric solutions in velocity space $\langle \unicode[STIX]{x1D6FF}f_{k}\rangle$ , which are decoupled from the azimuthally asymmetric solutions $\unicode[STIX]{x1D6FF}\tilde{f}_{k}=\unicode[STIX]{x1D6FF}f_{k}-\langle \unicode[STIX]{x1D6FF}f_{k}\rangle$ . Finally, we rewrite (2.9) by normalizing time to $kv_{th0}$ with $v_{th0}=\sqrt{2T_{0}/m}$ the electron thermal velocity, $\unicode[STIX]{x1D6FF}n_{k}$ to $n_{0}$ , and $v_{z}$ to $v_{th}$ , yielding

(2.10) $$\begin{eqnarray}\displaystyle \text{i}\frac{\unicode[STIX]{x2202}\langle \unicode[STIX]{x1D6FF}f_{k}\rangle }{\unicode[STIX]{x2202}t}-v_{z}\langle \unicode[STIX]{x1D6FF}f_{k}\rangle -\frac{v_{z}\unicode[STIX]{x1D6FF}n_{k}}{\unicode[STIX]{x1D6FC}_{D}}=\text{i}\frac{\langle C(f_{M}\unicode[STIX]{x1D6FF}f_{k})\rangle }{f_{M}}, & & \displaystyle\end{eqnarray}$$

where we defined $\unicode[STIX]{x1D6FC}_{D}=k^{2}\unicode[STIX]{x1D706}_{D}^{2}$ with $\unicode[STIX]{x1D706}_{D}=\sqrt{T_{0}/(4\unicode[STIX]{x03C0}e^{2}n_{0})}$ the Debye length, and where the collision frequency coefficient present in $C(f_{M}\unicode[STIX]{x1D6FF}f_{k})$ is now in units of $kv_{th}$ . Furthermore, we define $\unicode[STIX]{x1D708}$ as the electron–ion collision frequency normalized to $kv_{th}$ , namely

(2.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D708}=\frac{\unicode[STIX]{x1D708}_{i}}{kv_{th}}. & & \displaystyle\end{eqnarray}$$

We note that the number of particles in a Debye sphere $N$ can be cast in terms of $\unicode[STIX]{x1D708}$ and $\unicode[STIX]{x1D6FC}_{D}$ as

(2.12) $$\begin{eqnarray}\displaystyle N=\frac{1}{9}\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D708}\sqrt{\unicode[STIX]{x1D6FC}_{D}}}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D706}$ is the Coulomb logarithm.

Following Jorge et al. (Reference Jorge, Ricci and Loureiro2017, Reference Jorge, Ricci and Loureiro2018), we solve the linearized kinetic equation, equation (2.10), at arbitrary collisionalities by expanding the perturbed distribution function $\langle \unicode[STIX]{x1D6FF}f_{k}\rangle$ into an orthogonal Hermite–Laguerre polynomial basis of the form

(2.13) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D6FF}f_{k}\rangle =\mathop{\sum }_{p,j=0}^{\infty }\frac{N^{pj}}{\sqrt{2^{p}p!}}H_{p}(v_{z})L_{j}(v_{\bot }^{2}), & & \displaystyle\end{eqnarray}$$

where $H_{p}$ are physicists’ Hermite polynomials of order $p$ , defined by the Rodrigues’ formula

(2.14) $$\begin{eqnarray}\displaystyle H_{p}(x)=(-1)^{p}e^{x^{2}}\frac{\text{d}^{p}}{\text{d}x^{p}}e^{-x^{2}}, & & \displaystyle\end{eqnarray}$$

and normalized via the orthogonality condition based on the scalar product with weight factor $e^{-x^{2}}$

(2.15) $$\begin{eqnarray}\displaystyle \int _{-\infty }^{\infty }\text{d}xH_{p}(x)H_{p^{\prime }}(x)e^{-x^{2}}=2^{p}p!\sqrt{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D6FF}_{pp^{\prime }}, & & \displaystyle\end{eqnarray}$$

and $L_{j}$ the Laguerre polynomials of order $j$ , defined by the corresponding Rodrigues’ formula

(2.16) $$\begin{eqnarray}\displaystyle L_{j}(x)=\frac{e^{x}}{j!}\frac{\text{d}^{\,j}}{\text{d}x^{\,j}}(e^{-x}x^{\,j}), & & \displaystyle\end{eqnarray}$$

and orthonormal with respect to the weight $e^{-x}$

(2.17) $$\begin{eqnarray}\displaystyle \int _{0}^{\infty }\text{d}xL_{j}(x)L_{j^{\prime }}(x)e^{-x}=\unicode[STIX]{x1D6FF}_{jj^{\prime }}. & & \displaystyle\end{eqnarray}$$

Due to the orthogonality of the Hermite–Laguerre basis, the coefficients $N^{pj}=N^{pj}(k,t)$ of the expansion in (2.13) can be computed via the expression

(2.18) $$\begin{eqnarray}\displaystyle N^{pj}=\int \frac{H_{p}(v_{z})L_{j}(v_{\bot }^{2})\langle \unicode[STIX]{x1D6FF}f_{k}\rangle }{\sqrt{2^{p}p!}}\frac{e^{-v_{z}^{2}-v_{\bot }^{2}}}{\sqrt{\unicode[STIX]{x03C0}}}\,\text{d}v_{z}\,\text{d}v_{\bot }^{2}. & & \displaystyle\end{eqnarray}$$

With respect to a grid treatment using collocation points based on orthogonal polynomials (Belli & Candy Reference Belli and Candy2008; Landreman & Ernst Reference Landreman and Ernst2013), we note that the Hermite–Laguerre decomposition used in this work allows us to find the Rosenbluth potentials analytically as linear combinations of the $N^{pj}$ moments of the distribution function, to find linear and nonlinear closures at an arbitrary number of moments and mean-free path, and to make much more evident the role of the parallel and perpendicular phase-mixing processes.

By projecting the Boltzmann equation, equation (2.10), onto a Hermite–Laguerre basis, a system of differential equations (henceforth called moment hierarchy) for the coefficients $N^{pj}$ is obtained

(2.19) $$\begin{eqnarray}\displaystyle \text{i}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}N^{pj}=\sqrt{\frac{p+1}{2}}N^{p+1j}+\sqrt{\frac{p}{2}}N^{p-1j}+\frac{N^{00}}{\unicode[STIX]{x1D6FC}_{D}}\frac{\unicode[STIX]{x1D6FF}_{p,1}\unicode[STIX]{x1D6FF}_{j,0}}{\sqrt{2}}+\text{i}C^{pj}, & & \displaystyle\end{eqnarray}$$

with $C^{pj}$ the projection of the linearized collision operator onto a Hermite–Laguerre basis

(2.20) $$\begin{eqnarray}\displaystyle C^{pj}=\int \frac{H_{p}(v_{z})L_{j}(v_{\bot }^{2})\langle C(f_{M}\unicode[STIX]{x1D6FF}f_{k})\rangle }{\sqrt{2^{p}p!}}\,\text{d}v_{z}\,\text{d}v_{\bot }^{2}. & & \displaystyle\end{eqnarray}$$

The Coulomb collisional moments $C^{pj}$ are derived leveraging the result in Ji & Held (Reference Ji and Held2006), where the collision operator is projected onto a tensorial basis of the form $\boldsymbol{p}^{ls}=\boldsymbol{P}^{l}(\boldsymbol{c})L_{s}^{l+1/2}(c^{2})$ , with the irreducible tensorial Hermite polynomials defined by the recurrence $\boldsymbol{P}^{l+1}(\boldsymbol{c})=\boldsymbol{c}\boldsymbol{P}^{l}(\boldsymbol{c})-c^{2}\unicode[STIX]{x2202}_{\boldsymbol{c}}\boldsymbol{P}^{l}(\boldsymbol{c})/(2l+1)$ , being $\boldsymbol{P}^{0}(\boldsymbol{c})=1$ with $\boldsymbol{c}=\boldsymbol{v}/v_{th}$ , and the associated Laguerre polynomials $L_{s}^{l+1/2}(x)=\sum _{m=0}^{s}L_{sm}^{l}x^{m}$ with $L_{sm}^{l}=[(-1)^{m}(l+s+1/2)!]/[(s-m)!(l+m+1/2)!m!]$ . Indeed, expanding the distribution function as $\unicode[STIX]{x1D6FF}f_{a}=\sum _{l,s}\boldsymbol{M}_{a}^{ls}\boldsymbol{\cdot }\boldsymbol{p}^{ls}/\unicode[STIX]{x1D70E}_{s}^{l}$ , and with $\unicode[STIX]{x1D70E}_{s}^{l}=l!(l+s+1/2)!/(2^{l}(l+1/2)!s!)$ , Ji and Held showed that the linearized collision operator can be written as

(2.21) $$\begin{eqnarray}\displaystyle C(f_{M}\unicode[STIX]{x1D6FF}f)=f_{M}\mathop{\sum }_{b}\mathop{\sum }_{l,s=0}^{\infty }\frac{\boldsymbol{P}^{l}(\hat{v})}{\unicode[STIX]{x1D70E}_{s}^{l}}\boldsymbol{\cdot }(\boldsymbol{M}_{e}^{ls}\unicode[STIX]{x1D708}_{eb}^{ls,0}+\boldsymbol{M}_{b}^{ls}\unicode[STIX]{x1D708}_{eb}^{0,ls}), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D708}_{eb}^{ls,0}(v)$ and $\unicode[STIX]{x1D708}_{eb}^{0,ls}(v)$ are linear combinations of the error function and its derivatives (for their expression, see Ji & Held (Reference Ji and Held2006)), and represent the test-particle and field-particle (back-reaction) parts of the linearized collision operator, respectively. We remark that a similar expansion in Legendre-associated Laguerre polynomials was used in Brantov et al. (Reference Brantov, Bychenkov and Rozmus2012) in order to derive a simplified dispersion relation applicable to the study of EPW, ion-acoustic waves and entropy modes.

In order to evaluate $C^{pj}$ , we Fourier transform in space and average the operator $C(f_{M}\unicode[STIX]{x1D6FF}f)$ in (2.21) according to (2.8), using the averaging identity $\langle \boldsymbol{P}^{l}(\boldsymbol{c})\rangle =c^{l}P_{l}(v_{z}/v)\boldsymbol{P}^{l}(\hat{e}_{z})$ , with $P_{l}(x)=\unicode[STIX]{x2202}_{x}^{l}(x^{2}-1)^{l}/(2^{l}l!)$ the Legendre polynomials. This yields

(2.22) $$\begin{eqnarray}\displaystyle \langle C(f_{M}\unicode[STIX]{x1D6FF}f_{k})\rangle =\mathop{\sum }_{b}\mathop{\sum }_{l,s=0}^{\infty }(C_{eb}^{ls,0}+C_{eb}^{0,ls}), & & \displaystyle\end{eqnarray}$$

where $C_{eb}^{ls,0}$ is the averaged test-particle operator

(2.23) $$\begin{eqnarray}\displaystyle C_{eb}^{ls,0}=f_{M}\unicode[STIX]{x1D708}_{eb}^{ls,0}(v)P_{l}\left(\frac{v_{z}}{v}\right)M_{e}^{ls}\frac{2^{l}(l!)^{2}}{(2l!)}, & & \displaystyle\end{eqnarray}$$

and $C_{eb}^{0,ls}$ the field-particle (back-reaction) operator

(2.24) $$\begin{eqnarray}\displaystyle C_{eb}^{0,ls}=f_{M}\unicode[STIX]{x1D708}_{eb}^{0,ls}(v)P_{l}\left(\frac{v_{z}}{v}\right)M_{b}^{ls}\frac{2^{l}(l!)^{2}}{(2l!)}, & & \displaystyle\end{eqnarray}$$

with the fluid moments defined as

(2.25) $$\begin{eqnarray}\displaystyle M_{b}^{ls}=\mathop{\sum }_{p=0}^{l+2s}\mathop{\sum }_{j=0}^{s+l/2}T_{ls}^{pj}N_{b}^{pj}\sqrt{\frac{2^{p}p!}{\unicode[STIX]{x1D70E}_{s}^{l}}}. & & \displaystyle\end{eqnarray}$$

The basis transformation coefficients $T_{ls}^{pj}$ in (2.25) are defined by

(2.26) $$\begin{eqnarray}\displaystyle v^{l}P_{l}\left(\frac{v_{z}}{v}\right)L_{s}^{l+1/2}(v^{2})=\mathop{\sum }_{p=0}^{l+2s}\mathop{\sum }_{j=0}^{s+l/2}T_{ls}^{pj}H_{p}(v_{z})L_{j}(v_{\bot }^{2}), & & \displaystyle\end{eqnarray}$$

and their closed form expression can be found in Jorge et al. (Reference Jorge, Ricci and Loureiro2017). While an expansion in tensorial Hermite polynomials $\boldsymbol{P}^{l}(\boldsymbol{c})$ allows us to conveniently express the linearized collision operator in terms of $M_{b}^{ls}$ moments, the basis transformation of (2.26) is needed to cast the velocity dependence of $\langle C(f_{M}\unicode[STIX]{x1D6FF}f_{k})\rangle$ in a Hermite–Laguerre polynomial basis, and to calculate its velocity moments.

We now use the expression of $\langle C(f_{M}\unicode[STIX]{x1D6FF}f_{k})\rangle$ contained in (2.22) and inject it in (2.18). By defining the fluid moments $A_{eb}^{lts}$ as

(2.27) $$\begin{eqnarray}\displaystyle A_{eb}^{lts}=\int v^{l}L_{t}^{l+1/2}(v^{2})f_{M}\unicode[STIX]{x1D708}_{b}^{ls,0}(v)\,\text{d}v, & & \displaystyle\end{eqnarray}$$

and $B_{eb}^{lts}$ as

(2.28) $$\begin{eqnarray}\displaystyle B_{eb}^{lts}=\int v^{l}L_{t}^{l+1/2}(v^{2})f_{M}\unicode[STIX]{x1D708}_{b}^{0,ls}(v)\,\text{d}v, & & \displaystyle\end{eqnarray}$$

the resulting collision operator moments $C^{pj}$ can be written as

(2.29) $$\begin{eqnarray}\displaystyle C^{pj}=\mathop{\sum }_{b}\mathop{\sum }_{s=0}^{\infty }\mathop{\sum }_{l=0}^{p+2j}\mathop{\sum }_{t=0}^{j+\lfloor p/2\rfloor }\frac{(T^{-1})_{pj}^{lt}2^{l}(l!)^{2}}{(2l)!\unicode[STIX]{x1D70E}_{s}^{l}\sqrt{2^{p}p!}}\frac{\unicode[STIX]{x1D708}_{b}}{(2l+1)}(M_{e}^{ls}A_{eb}^{lts}+M_{b}^{ls}B_{eb}^{lts}), & & \displaystyle\end{eqnarray}$$

where we introduce the inverse transformation coefficients

(2.30) $$\begin{eqnarray}\displaystyle (T^{-1})_{pj}^{lt}=T_{lt}^{pj}\frac{\sqrt{\unicode[STIX]{x03C0}}2^{p}p!(l+1/2)t!}{(t+l+1/2)!}. & & \displaystyle\end{eqnarray}$$

The analytical expressions for $A_{eb}^{lts}$ and $B_{eb}^{lts}$ suitable for numerical implementation are given in Ji & Held (Reference Ji and Held2006). The moments of the collision operator, $C^{pj}$ , correspond to the ones derived in Jorge et al. (Reference Jorge, Ricci and Loureiro2018) for the study of drift waves, and can also be obtained by linearizing the electron collisional moments presented in Jorge et al. (Reference Jorge, Ricci and Loureiro2017).

Besides the Coulomb collision operator, the Hermite–Laguerre expansion described above can be advantageously applied to describe other collision operators. We consider here the Lenard–Bernstein (Lenard & Bernstein Reference Lenard and Bernstein1958), the Dougherty (Dougherty Reference Dougherty1964) and the electron–ion collision operators that are used for comparison with the full Coulomb one. The Lenard–Bernstein and Dougherty operators are implemented in a number of advanced kinetic codes (Grandgirard et al. Reference Grandgirard, Abiteboul, Bigot, Cartier-Michaud, Crouseilles, Dif-Pradalier, Ehrlacher, Esteve, Garbet and Ghendrih2016; Loureiro et al. Reference Loureiro, Dorland, Fazendeiro, Kanekar, Mallet, Vilelas and Zocco2016; Nakata et al. Reference Nakata, Honda, Yoshida, Urano, Nunami, Maeyama, Watanabe and Sugama2016; Pan et al. Reference Pan, Told, Shi, Hammett and Jenko2018), and are frequently used to introduce collisional effects in weakly collisional plasmas (Zocco & Schekochihin Reference Zocco and Schekochihin2011; Zocco et al. Reference Zocco, Loureiro, Dickinson, Numata and Roach2015; Shi et al. Reference Shi, Hammett, Stoltzfus-Dueck and Hakim2017; Mandell, Dorland & Landreman Reference Mandell, Dorland and Landreman2018). Therefore, a comparison between the Coulomb and the Lenard–Bernstein and Dougherty operators, even in simplified systems such as the case of EPW, is important to determine the accuracy and validity of these operators. The Lenard–Bernstein collision operator $C_{LB}(f)$ , first derived in 1891 by Lord Rayleigh (Wax Reference Wax1954) is of the Fokker–Planck type. It conserves particle number and satisfies the H-theorem, and it can be written as (Lenard & Bernstein Reference Lenard and Bernstein1958)

(2.31) $$\begin{eqnarray}\displaystyle C_{LB}(f)=\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{v}}\boldsymbol{\cdot }\left(\boldsymbol{v}f+\frac{v_{th}^{2}}{2}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\boldsymbol{v}}\right). & & \displaystyle\end{eqnarray}$$

This operator can be derived from the Fokker–Planck equation, equation (2.2), by assuming $(m/m_{b})\unicode[STIX]{x2202}_{\boldsymbol{v}}H_{b}=\boldsymbol{v}$ and $\unicode[STIX]{x2202}_{\boldsymbol{v}}\unicode[STIX]{x2202}_{\boldsymbol{v}}G_{b}=-\unicode[STIX]{x1D644}v_{th}^{2}/2$ with $\unicode[STIX]{x1D644}$ the identity matrix. By projecting the Lenard–Bernstein operator onto a Hermite–Laguerre basis according to (2.20), one obtains

(2.32) $$\begin{eqnarray}\displaystyle C_{LB}^{pj}=-\unicode[STIX]{x1D708}(p+2j)N^{pj}. & & \displaystyle\end{eqnarray}$$

Equation (2.32) can then be used in the moment-hierarchy equation (2.19), yielding

(2.33) $$\begin{eqnarray}\displaystyle \text{i}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}N^{pj}=\sqrt{\frac{p+1}{2}}N^{p+1j}+\sqrt{\frac{p}{2}}N^{p-1j}+\frac{N^{00}}{\unicode[STIX]{x1D6FC}_{D}}\frac{\unicode[STIX]{x1D6FF}_{p,1}\unicode[STIX]{x1D6FF}_{j,0}}{\sqrt{2}}-\text{i}\unicode[STIX]{x1D708}(p+2j)N^{pj}. & & \displaystyle\end{eqnarray}$$

The linearized Dougherty collision operator $C_{D}(f)$ , on the other hand, adds the necessary field-particle collisional terms to the Lenard–Bernstein operator in order to provide momentum and energy conservation properties. Namely, it sets $(m/m_{b})\unicode[STIX]{x2202}_{\boldsymbol{v}}H_{b}=\boldsymbol{v}-\boldsymbol{u}$ , with $\boldsymbol{u}=\int \boldsymbol{v}f\,\text{d}v_{z}\,\text{d}v_{\bot }^{2}\,\text{d}\unicode[STIX]{x1D711}/n_{0}$ , and $\unicode[STIX]{x2202}_{\boldsymbol{v}}\unicode[STIX]{x2202}_{\boldsymbol{v}}G_{b}=-\unicode[STIX]{x1D644}T/m_{a}$ with $T=\int m(\boldsymbol{v}-\boldsymbol{u})^{2}f\,\text{d}v_{z}\,\text{d}v_{\bot }^{2}\,\text{d}\unicode[STIX]{x1D711}/(3n_{0})=(\sqrt{2}N^{20}-2N^{01})/3$ . The Hermite–Laguerre moments of the linearized Dougherty collision operator $C_{D}^{pj}$ are given by

(2.34) $$\begin{eqnarray}\displaystyle C_{D}^{pj}=-\unicode[STIX]{x1D708}[(p+2j)N^{pj}-N^{10}\unicode[STIX]{x1D6FF}_{p1}\unicode[STIX]{x1D6FF}_{j0}+T(\sqrt{2}\unicode[STIX]{x1D6FF}_{p0}\unicode[STIX]{x1D6FF}_{j1}-2\unicode[STIX]{x1D6FF}_{p2}\unicode[STIX]{x1D6FF}_{j0})], & & \displaystyle\end{eqnarray}$$

yielding the moment-hierarchy equation

(2.35) $$\begin{eqnarray}\displaystyle \text{i}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}N^{pj} & = & \displaystyle \sqrt{\frac{p+1}{2}}N^{p+1j}+\sqrt{\frac{p}{2}}N^{p-1j}+\frac{N^{00}}{\unicode[STIX]{x1D6FC}_{D}}\frac{\unicode[STIX]{x1D6FF}_{p,1}\unicode[STIX]{x1D6FF}_{j,0}}{\sqrt{2}}\nonumber\\ \displaystyle & & \displaystyle -\,\text{i}\unicode[STIX]{x1D708}[(p+2j)N^{pj}-N^{10}\unicode[STIX]{x1D6FF}_{p1}\unicode[STIX]{x1D6FF}_{j0}+T(\sqrt{2}\unicode[STIX]{x1D6FF}_{p0}\unicode[STIX]{x1D6FF}_{j1}-2\unicode[STIX]{x1D6FF}_{p2}\unicode[STIX]{x1D6FF}_{j0})].\end{eqnarray}$$

We note that the moment hierarchies with the Lenard–Bernstein or the Dougherty collision operator, equations (2.33) and (2.35), respectively, do not couple different Laguerre moments and, therefore, one can focus on obtaining the coefficients $N^{p0}$ for solving the moment hierarchy.

The Hermite–Laguerre expansion procedure can also be applied to the electron–ion operator $C_{ei}(f)$ that is often used for EPW studies (Epperlein et al. Reference Epperlein, Short and Simon1992; Banks et al. Reference Banks, Brunner, Berger and Tran2016), that is

(2.36) $$\begin{eqnarray}\displaystyle C_{ei}(f)=\frac{\unicode[STIX]{x1D708}}{2v^{3}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\left[(1-\unicode[STIX]{x1D709}^{2})\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\right], & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D709}=v_{z}/v$ . This operator describes the pitch-angle scattering of electrons due to collisions with ions. By projecting (2.36) into a Hermite–Laguerre basis, we obtain (Jorge et al. Reference Jorge, Ricci and Loureiro2017)

(2.37) $$\begin{eqnarray}\displaystyle C_{ei}^{pj}=-\frac{\unicode[STIX]{x1D708}}{8\unicode[STIX]{x03C0}^{3/2}}\mathop{\sum }_{l=0}^{p+2j}\mathop{\sum }_{f=0}^{j+\lfloor p/2\rfloor }\frac{(T^{-1})_{pj}^{lf}}{\sqrt{2^{p}p!}}\mathop{\sum }_{s=0}^{\infty }A_{ei}^{lf,s}M_{e}^{ls}, & & \displaystyle\end{eqnarray}$$

where the $A_{ei}^{lf,s}$ coefficients are given by

(2.38) $$\begin{eqnarray}\displaystyle A_{ei}^{lf,s}=\frac{l(l+1)}{l+1/2}\frac{2^{l}(l!)^{2}}{(2l)!}\mathop{\sum }_{m=0}^{f}\mathop{\sum }_{n=0}^{s}\frac{L_{fm}^{l}L_{sn}^{l}}{\sqrt{\unicode[STIX]{x1D70E}_{s}^{l}}}(l+m+n-1)!. & & \displaystyle\end{eqnarray}$$

As an aside, we note that previous studies on EPW have shown that the solutions of the linearized Boltzmann equation are, in fact, sensitive to the discretization method used. For example, it was shown that finite-difference methods, when applied to the problem of EPW, produce a number of numerical, non-physical modes with a rather small damping rate that do not lie in the vicinity of the collisionless solutions, even for weak collisionalities and a very high resolution (Bratanov et al. Reference Bratanov, Jenko, Hatch and Brunner2013). On the other hand, a discretization scheme based on a Hermite–Laguerre polynomial decomposition yields a large number of roots that lie in the vicinity of the collisionless solution. Since previous EPW studies using the electron–ion collision operator have been performed using a discretization of the distribution function into a set of Legendre polynomials, i.e. (Epperlein et al. Reference Epperlein, Short and Simon1992; Brantov et al. Reference Brantov, Bychenkov and Rozmus2012; Banks et al. Reference Banks, Brunner, Berger and Tran2016)

(2.39) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D6FF}f_{k}\rangle =\mathop{\sum }_{l=0}^{\infty }a_{l}(v)P_{l}(\unicode[STIX]{x1D709}), & & \displaystyle\end{eqnarray}$$

as a test of our approach, we compare in § 6 our results with the Legendre decomposition in (2.39). By projecting the Boltzmann equation, equation (2.10), with an electron–ion collision operator into a Legendre basis, equation (2.39), the following moment-hierarchy equation is obtained

(2.40) $$\begin{eqnarray}\displaystyle \frac{\text{i}}{v}\frac{\text{d}a_{l}(v)}{\text{d}t}=\frac{l}{2l-1}a_{l-1}(v)+\frac{l+1}{2l+3}a_{l+1}(v)+\frac{\unicode[STIX]{x1D6FF}_{l,1}}{\unicode[STIX]{x1D6FC}_{D}}\int f_{M}v^{2}a_{0}(v)\,\text{d}v-\frac{\text{i}\unicode[STIX]{x1D708}}{v^{4}}l(l+1)a_{l}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

A relation between the Hermite–Laguerre $N^{pj}$ and Legendre moments $a_{l}$ can be found by comparing (2.13) and (2.39), yielding

(2.41) $$\begin{eqnarray}\displaystyle a_{l}(v)=\mathop{\sum }_{p=0}^{\infty }\mathop{\sum }_{j=0}^{\infty }\mathop{\sum }_{s=0}^{p+2j}\mathop{\sum }_{t=0}^{j+\lfloor p/2\rfloor }(T^{-1})_{pj}^{st}\frac{N^{pj}}{\sqrt{2^{p}p!}}v^{s}L_{t}^{s+1/2}(v^{2})\unicode[STIX]{x1D6FF}_{ls}, & & \displaystyle\end{eqnarray}$$

and

(2.42) $$\begin{eqnarray}\displaystyle N^{pj}=\mathop{\sum }_{s=0}^{p+2j}\mathop{\sum }_{t=0}^{j+\lfloor p/2\rfloor }\frac{(T^{-1})_{pj}^{st}}{\sqrt{2^{p}p!}(2l+1)}\int a_{s}(v)v^{s+2}L_{t}^{s+1/2}(v^{2})\,\text{d}v. & & \displaystyle\end{eqnarray}$$

3 Collisionless dispersion relation

As a first step in the analysis of EPW, and for comparison with the results in the presence of collisions, we derive the EPW dispersion relation in the collisionless limit. We first Fourier transform in time the collisionless limit of the moment-hierarchy equation, equation (2.19), by imposing $\unicode[STIX]{x1D6FF}f_{k}\sim \text{e}^{(\unicode[STIX]{x1D6FE}+\text{i}\unicode[STIX]{x1D714})t}$ , obtaining

(3.1) $$\begin{eqnarray}\displaystyle \text{i}(\unicode[STIX]{x1D6FE}+\text{i}\unicode[STIX]{x1D714})N^{pj}=\sqrt{\frac{p+1}{2}}N^{p+1j}+\sqrt{\frac{p}{2}}N^{p-1j}+\frac{N^{00}}{\unicode[STIX]{x1D6FC}_{D}}\frac{\unicode[STIX]{x1D6FF}_{p,1}\unicode[STIX]{x1D6FF}_{j,0}}{\sqrt{2}}. & & \displaystyle\end{eqnarray}$$

A closed form solution of the collisionless moment hierarchy in (3.1) can be obtained by dividing the Boltzmann equation, equation (2.10), by the resonant $\text{i}\unicode[STIX]{x1D6FE}-v_{z}$ factor, multiplying it by the Hermite–Laguerre polynomial basis functions and, finally, integrating it over velocity space. This yields

(3.2) $$\begin{eqnarray}\displaystyle N^{pj}=-\frac{N^{00}}{\unicode[STIX]{x1D6FC}_{D}}\left[-\text{i}(\unicode[STIX]{x1D6FE}+\text{i}\unicode[STIX]{x1D714})\frac{(-1)^{p}}{\sqrt{2^{p}p!}}Z^{(p)}(\unicode[STIX]{x1D714}-\text{i}\unicode[STIX]{x1D6FE})+\unicode[STIX]{x1D6FF}_{p,0}\right]\unicode[STIX]{x1D6FF}_{j,0}, & & \displaystyle\end{eqnarray}$$

where $Z^{(p)}$ is the $p$ th derivative of the plasma dispersion function $Z^{(0)}$ , defined by

(3.3) $$\begin{eqnarray}\displaystyle Z^{(p)}(u)=\frac{(-1)^{p}}{\sqrt{\unicode[STIX]{x03C0}}}\int _{-\infty }^{\infty }\frac{H_{p}(x)e^{-x^{2}}}{x-u}\,\text{d}x. & & \displaystyle\end{eqnarray}$$

By setting $(p,j)=(0,0)$ in (3.2), the collisionless dispersion relation is found

(3.4) $$\begin{eqnarray}\displaystyle D=1+\unicode[STIX]{x1D6FC}_{D}-\text{i}(\unicode[STIX]{x1D6FE}+\text{i}\unicode[STIX]{x1D714})Z(\unicode[STIX]{x1D714}-\text{i}\unicode[STIX]{x1D6FE})=0. & & \displaystyle\end{eqnarray}$$

Alternatively, equation (3.4) can be derived from the collisionless limit of the Boltzmann equation, equation (2.10), upon division by the factor $\text{i}\unicode[STIX]{x1D6FE}-\unicode[STIX]{x1D714}-v_{z}$ and integration with respect to $v_{z}$ . The numerical solution of (3.4) is obtained by discretizing $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D714}$ into a two-dimensional $[\unicode[STIX]{x1D714},\unicode[STIX]{x1D6FE}]$ grid, evaluating $D$ on the grid, and storing the values where $[\text{Re}(D),\text{Im}(D)]$ vanishes. To evaluate $Z$ , we make use of the identity $Z(x)=\text{i}\sqrt{\unicode[STIX]{x03C0}}e^{-x^{2}}\text{erfc}(-ix)$ with $\text{erfc}(x)=1-\text{erf}(x)$ and $\text{erf}(x)$ the error function, and use the algorithm developed in Gautschi (Reference Gautschi1970) to numerically compute $\text{erf}(x)$ for complex arguments.

Figure 1.

Collisionless frequency (blue line) and damping rate (red line) of the least-damped solution of the collisionless dispersion relation, equation (3.4), as a function of $\unicode[STIX]{x1D6FC}_{D}$ , for $0.05<\unicode[STIX]{x1D6FC}_{D}<0.4$ .

The $\unicode[STIX]{x1D6FC}_{D}$ dependence of the least-damped solution of (3.4) is shown in figure 1, where both its damping rate $\unicode[STIX]{x1D6FE}$ and frequency $\unicode[STIX]{x1D714}$ are seen to be monotonic functions of $\unicode[STIX]{x1D6FC}_{D}$ , which is in agreement with previous EPW studies (Banks et al. Reference Banks, Brunner, Berger, Arrighi and Tran2017). In the following, without loss of generality and similarly to previous studies of collisional damping of EPW (Banks et al. Reference Banks, Brunner, Berger and Tran2016, Reference Banks, Brunner, Berger, Arrighi and Tran2017), we select the value of $\unicode[STIX]{x1D6FC}_{D}=0.09$ when fixed $\unicode[STIX]{x1D6FC}_{D}$ studies are performed, which corresponds to $k\unicode[STIX]{x1D706}_{D}=\sqrt{\unicode[STIX]{x1D6FC}_{D}}=0.3$ . While this value of $\unicode[STIX]{x1D6FC}_{D}$ is typical for EPW driven by stimulated Raman scattering (Brunner & Valeo Reference Brunner and Valeo2004; Winjum et al. Reference Winjum, Berger, Chapman, Banks and Brunner2013), we add that the particular choice of $\unicode[STIX]{x1D6FC}_{D}$ has no quantitative impact on the conclusions we draw.

4 Temporal evolution of EPW

In this section, the moment-hierarchy equation, equation (2.19), is solved numerically as a time-evolution problem using an implicit variable-step variable-order solver that employs backward differentiation formulas (Shampine Reference Shampine2002) with a maximum time step of $10^{-6}kv_{th}$ . For the numerical solution, the moment hierarchy is truncated at a maximum Hermite–Laguerre index $(P,J)$ by setting

(4.1) $$\begin{eqnarray}\displaystyle N^{pj}=0,\quad \text{for}~(p,j)>(P,J). & & \displaystyle\end{eqnarray}$$

We consider as initial condition $N^{pj}(t=0)=\unicode[STIX]{x1D6FF}_{p0}\unicode[STIX]{x1D6FF}_{j0}$ , such that the perturbed density and electrostatic potential are initially excited, while higher moments of the distribution function are set to zero. The initial perturbation of the velocity distribution function is then taken to be a Maxwellian in velocity space. The temporal evolution of $N^{00}$ (and therefore of $\unicode[STIX]{x1D719}$ ) is shown in figure 2 for different collisionalities and $\unicode[STIX]{x1D6FC}_{D}$ values. In this section, we focus on the oscillating initial phase of figure 2, where the EPW dominate the dynamics. We fit the amplitude of $N^{00}$ to an exponentially damped sinusoidal wave with real frequency $\unicode[STIX]{x1D714}$ and damping rate $\unicode[STIX]{x1D6FE}$ , taking into account a minimum of three oscillation periods. The later phase, where a purely damped behaviour is observed at higher collisionalities due to the presence of an entropy mode, is investigated in § 5.

Figure 2.

Time evolution of the absolute value of $N^{00}$ using a truncation with $(P,J)=(18,2)$ , evaluated using the full linearized Coulomb collision operator. Different values of $\unicode[STIX]{x1D708}$ and $\unicode[STIX]{x1D6FC}_{D}$ are shown.

A convergence study with the truncation indices $(P,J)$ is shown in figure 3. Convergence is observed for $(P,J)=(18,2)$ in the range of collisionalities and $\unicode[STIX]{x1D6FC}_{D}$ investigated (a variation of less than 3 % is observed between damping rates evaluated with a truncation at $(P,J)=(18,2)$ and a truncation at higher values of $P$ and $J$ ). The values of oscillation frequency $\unicode[STIX]{x1D714}(\unicode[STIX]{x1D6FC}_{D},\unicode[STIX]{x1D708})$ and damping rate $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D6FC}_{D},\unicode[STIX]{x1D708})$ obtained as a fit of the initial damping phase are shown in figure 4 for $0.075<\unicode[STIX]{x1D6FC}_{D}<0.2$ and $0.015<\unicode[STIX]{x1D708}<0.7$ using the Coulomb collision operator, where a truncation at $(P,J)=(18,2)$ is used. Such values of $(P,J)$ are in line with the estimate found in Jorge et al. (Reference Jorge, Ricci and Loureiro2018) which, by underestimating the effect of the collisional term $C^{pj}$ in the moment-hierarchy equation leads to $P\sim 4/\sqrt{2}\unicode[STIX]{x1D708}$ and $J\sim 2$ . For $\unicode[STIX]{x1D708}\sim 0.1$ , this yields $(P,J)\sim (30,2)$ . The largest deviation of the damping rate from the collisionless case is seen to occur for large values of collisionality and small $\unicode[STIX]{x1D6FC}_{D}$ . This is expected, as for large $\unicode[STIX]{x1D708}$ and small $\unicode[STIX]{x1D6FC}_{D}$ the collisional fluid limit is retrieved. The dependence on $\unicode[STIX]{x1D6FC}_{D}$ may be attributed to the decreasing magnitude of the Landau damping rate for decreasing $\unicode[STIX]{x1D6FC}_{D}$ (see figure 1). This makes the ratio between the collisional and the collisionless damping rates increasingly larger. Finally, we remark that the presence of several competing eigenmodes in the initial transients of the temporal evolution of $N^{00}$ contribute to the presence of a transition at $\unicode[STIX]{x1D6FC}_{D}\sim 0.1$ visible in figure 4.

Figure 3.

Comparison of the normalized damping rate $\unicode[STIX]{x1D6FE}$ for $\unicode[STIX]{x1D6FC}_{D}=0.09$ as a function of the normalized collision frequency considering a truncation at different values of $(P,J)$ (solid lines) and using the full linearized Coulomb collision operator. The collisionless least-damped Landau solution is shown for comparison (dashed blue line). All frequencies are normalized to $kv_{th}$ .

Figure 4.

Damping rate $\unicode[STIX]{x1D6FE}$ (a) and oscillation frequency $\unicode[STIX]{x1D714}$ (b) of the electron-plasma wave obtained from the moment-hierarchy equation, equation (2.19), as a function of $\unicode[STIX]{x1D6FC}_{D}$ and $\unicode[STIX]{x1D708}$ for $(P,J)=(18,2)$ . The full linearized Coulomb collision operator is considered.

Following a long standing tradition established by Jackson (Jackson Reference Jackson1960), and used in consequent EPW studies (Opher, Morales & Leboeuf Reference Opher, Morales and Leboeuf2002), we also display the frequency and the damping rate of EPW normalized to the plasma frequency $\unicode[STIX]{x1D714}_{pe}=v_{th}/\unicode[STIX]{x1D706}_{D}$ in figure 5 as a function of $k\unicode[STIX]{x1D706}_{D}$ for fixed $\unicode[STIX]{x1D708}_{ei}/\unicode[STIX]{x1D714}_{pe}$ . We show the results for the collisionless $\unicode[STIX]{x1D708}_{ei}/\unicode[STIX]{x1D714}_{pe}=0$ and the collisional $\unicode[STIX]{x1D708}_{ei}/\unicode[STIX]{x1D714}_{pe}=0.1$ case. Consistently with figure 4, it is observed that $\unicode[STIX]{x1D714}/\unicode[STIX]{x1D714}_{pe}$ has a negligible relative variation with the collision frequency when compared with the relative variation of the damping rate. Furthermore, as observed in Opher et al. (Reference Opher, Morales and Leboeuf2002), the damping rate of EPW with finite collision frequency is seen to approach the collisionless result for high values of $k\unicode[STIX]{x1D706}_{D}$ . A further comparison of the collisional component of the damping rate $\unicode[STIX]{x1D6FE}_{\text{coll}}$ , obtained with different collision models, is shown in figure 6, where $\unicode[STIX]{x1D6FE}_{\text{coll}}=\unicode[STIX]{x1D6FE}-\unicode[STIX]{x1D6FE}_{\text{collisionless}}$ with $\unicode[STIX]{x1D6FE}$ the total damping rate and $\unicode[STIX]{x1D6FE}_{\text{collisionless}}$ the collisionless Landau damping. Results considering Lenard–Bernstein, Dougherty, electron–ion and the full Coulomb collision operators are shown. We note that when the Lenard–Bernstein and the Dougherty operators are considered, only self-collisions are taken into account, and with the electron–ion operator only unlike-particle collisions are included. In general, the Coulomb operator yields a damping rate smaller than the Lenard–Bernstein and larger than the Dougherty one, with deviations of up to 50 % between different operators. The use of an electron–ion collision operator is preferable since it yields damping rates and frequencies similar, just slightly lower, than the Coulomb operator. The collisional damping rate, $\unicode[STIX]{x1D6FE}_{fluid}=-0.532\unicode[STIX]{x1D708}/2$ , obtained from a fluid description using the Braginskii equations (Banks et al. Reference Banks, Brunner, Berger, Arrighi and Tran2017), is also shown for comparison. We remark that the results in figure 6 for the collisional component of the damping rate of the fluid, purely electron–ion collisions, and the full Coulomb operator are in close agreement with the findings of Banks et al. (Reference Banks, Brunner, Berger, Arrighi and Tran2017).

Figure 5.

Oscillation frequency $\unicode[STIX]{x1D714}_{ph}=\unicode[STIX]{x1D714}kv_{th}$ and damping rate $\unicode[STIX]{x1D6FE}_{ph}=\unicode[STIX]{x1D6FE}kv_{th}$ of EPW in physical units, normalized to the plasma frequency $\unicode[STIX]{x1D714}_{pe}$ for the collisionless $\unicode[STIX]{x1D708}_{ei}/\unicode[STIX]{x1D714}_{pe}=0$ and collisional $\unicode[STIX]{x1D708}_{ei}/\unicode[STIX]{x1D714}_{pe}=0.1$ case.

Figure 6.

Difference between the collisionless damping rate $\unicode[STIX]{x1D6FE}$ with the one resulting from the moment-hierarchy equation, equation (2.19), with the full Coulomb, Lenard–Bernstein, Dougherty and electron–ion collision operators at $\unicode[STIX]{x1D6FC}_{D}=0.09$ and $(P,J)=(18,2)$ . The collisional damping rate $\unicode[STIX]{x1D6FE}_{fluid}=-0.266\unicode[STIX]{x1D708}$ obtained from a fluid description is also shown for comparison.

For $\unicode[STIX]{x1D708}\ll 1$ , it is seen that the damping rates of all solutions approach the collisionless limit regardless of the collision operator used. When $\unicode[STIX]{x1D708}$ is increased, figure 6 shows that the differences between the collision operators still persist. This allows us to draw arguments for the difference between the different collision operators by using a low number of moments. The lowest-order particle conservation $C^{00}=0$ and collisional friction $C^{10}=-\unicode[STIX]{x1D708}N^{10}$ moments $C^{pj}$ are the same between the Coulomb, Lenard–Bernstein and electron–ion collision operators, while the Dougherty operator has $C_{D}^{10}=0$ . This effectively reduces the damping rate evaluated with the Dougherty operator with respect to the Coulomb case, as seen in figure 6.

On the other hand, only the Coulomb, the Dougherty and the electron–ion collision operators are energy conserving, i.e. satisfying $(1/2)\int m(v_{z}^{2}+v_{\bot }^{2})C(f)\,\text{d}\boldsymbol{v}=0$ or, equivalently, $C^{20}=\sqrt{2}C^{01}$ , while the Lenard–Bernstein operator does not conserve energy. In fact, the energy moments are given by $C_{LB}^{01}=-\unicode[STIX]{x1D708}2N^{01}$ and $C_{LB}^{20}=-\unicode[STIX]{x1D708}2N^{20}$ , yielding $C_{LB}^{01}=C_{LB}^{20}$ . However, despite the additional conservation properties, the agreement of the Dougherty operator is rather poor, as seen in figure 6. We conclude therefore that the presence of additional momentum and energy conserving terms in the Dougherty operator with respect to the Lenard–Bernstein operator does not yield a damping rate closer to the Coulomb one. This was also pointed out in Jorge et al. (Reference Jorge, Ricci and Loureiro2018), where a similar framework was used to derive the growth rate of the drift-wave instability.

5 Entropy mode

We focus on the latter stage of the time evolution of $N^{00}$ shown in figures 2 and 7, where a purely damped behaviour is found at high collisionalities. In order to enhance the role of the zero-frequency mode, we consider the collision frequency $\unicode[STIX]{x1D708}=5$ , while decreasing the role of Landau damping by setting $\unicode[STIX]{x1D6FC}_{D}=0.01$ (the purely damped mode is not affected by the value of $\unicode[STIX]{x1D6FC}_{D}$ , if $\unicode[STIX]{x1D6FC}_{D}\ll 1$ ). Indeed, the transition to a purely damped behaviour is seen to occur at times that decrease with the collision frequency (see figure 7 (black)). The resulting time traces of $|N^{00}|$ using a full Coulomb collision operator are shown in figure 7(a) for $(P,J)=(18,0)$ , $(P,J)=(18,2)$ and $(P,J)=(18,4)$ , while time traces using the full Coulomb, electron–ion, Lenard–Bernstein and the Dougherty collision operators with $(P,J)=(18,2)$ are shown in figure 7(b). We observe that for the Coulomb and electron–ion cases, there is a transition to a purely damped mode at $t\simeq 7$ only when perpendicular velocity dynamics is introduced, $J\geqslant 2$ . At the same time, while for the Coulomb operator, the purely damped mode that sets the late time evolution of the system in figure 7 has a damping rate $\unicode[STIX]{x1D6FE}\simeq -0.202$ , the electron–ion collision operator yields a damping rate one order of magnitude smaller, $\unicode[STIX]{x1D6FE}\simeq -0.024$ . This purely damped decay is not present when the Lenard–Bernstein or the Dougherty operators are considered. We therefore conclude that in order to obtain the correct long-term behaviour of the Boltzmann equation, Coulomb self-collisions must be included in the description.

Figure 7.

Time evolution of the absolute value of $N^{00}$ with $\unicode[STIX]{x1D708}=5$ and $\unicode[STIX]{x1D6FC}_{D}=0.01$ . (a) Convergence study with a full Coulomb operator and $(P,J)=(18,0)$ (green), $(P,J)=(18,2)$ (red) and $(P,J)=(18,4)$ (blue). The time evolution of $N^{00}$ with $\unicode[STIX]{x1D708}=3.5$ and $(P,J)=(18,4)$ is also shown for comparison (black). (b) Truncation at $(P,J)=(18,2)$ using the full Coulomb collision operator (red), electron–ion collisions only (green), the Lenard–Bernstein (orange) and the Dougherty (blue) collision operators.

The long-term behaviour observed in figure 7 is due to the presence of the entropy mode (Banks et al. Reference Banks, Brunner, Berger and Tran2016). An analytical framework to model the entropy mode can be derived by noting that this is a purely damped mode with a damping rate much smaller than both the plasma and the collision frequencies. Previous studies on the collisional damping of EPW show that such a mode results from the effect of pitch-angle scattering at high collisionality (Banks et al. Reference Banks, Brunner, Berger and Tran2016). Considering that in the high-collisionality limit only the lowest-order terms in the expansion of $\langle \unicode[STIX]{x1D6FF}f_{k}\rangle$ in (2.13) play a role and that, according to figure 7, a finite perpendicular velocity-space resolution is an essential element for the entropy mode, we consider in the moment-hierarchy equation, equation (2.19), the six lowest-order Hermite–Laguerre expansion coefficients, namely $N^{00},N^{10},N^{20},N^{30},N^{01}$ and $N^{11}$ with the ordering $N^{30}\sim N^{11}\sim \unicode[STIX]{x1D716}N^{20}\sim \unicode[STIX]{x1D716}N^{01}\sim \unicode[STIX]{x1D716}\unicode[STIX]{x1D719}$ , with $\unicode[STIX]{x1D716}$ the small expansion parameter

(5.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}\sim \frac{1}{\unicode[STIX]{x1D708}}\sim \unicode[STIX]{x1D6FE}, & & \displaystyle\end{eqnarray}$$

so that the particle mean-free path $\unicode[STIX]{x1D706}_{mfp}=v_{th}/\unicode[STIX]{x1D708}_{ei}$ is small compared with typical wavelengths of the perturbed quantities, i.e. $k\unicode[STIX]{x1D706}_{mfp}\ll 1$ . Higher-order moments are considered to be $O(\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D719})$ . Charge neutrality is kept up to second order, i.e.

(5.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}_{D}\sim \unicode[STIX]{x1D716}^{2}. & & \displaystyle\end{eqnarray}$$

Using Poisson’s equation, we find that density perturbations $N^{00}$ are negligible when compared with electrostatic fluctuations, namely

(5.3) $$\begin{eqnarray}\displaystyle \frac{N^{00}}{\unicode[STIX]{x1D719}}=-\unicode[STIX]{x1D6FC}_{D}\ll 1. & & \displaystyle\end{eqnarray}$$

The moment-hierarchy equation, equation (2.19), at $(p,j)=(0,0)$ , shows that $N^{10}/N^{00}\sim \unicode[STIX]{x1D6FE}$ . Together with the estimate in (5.3), this yields

(5.4) $$\begin{eqnarray}\displaystyle \frac{N^{10}}{\unicode[STIX]{x1D719}}\sim \unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FC}_{D}\ll 1. & & \displaystyle\end{eqnarray}$$

Using the moment-hierarchy equation (2.19) and neglecting second-order terms in the parallel $(p,j)=(2,0)$ and perpendicular $(p,j)=(0,1)$ temperature equations, we find that

(5.5) $$\begin{eqnarray}\displaystyle \text{i}\unicode[STIX]{x1D6FE}N^{20}\simeq \sqrt{\frac{3}{2}}N^{30}-\text{i}\unicode[STIX]{x1D708}(0.45N^{01}+0.64N^{20}), & & \displaystyle\end{eqnarray}$$

and

(5.6) $$\begin{eqnarray}\displaystyle \text{i}\unicode[STIX]{x1D6FE}N^{01}\simeq \frac{N^{11}}{\sqrt{2}}-\text{i}\unicode[STIX]{x1D708}(0.32N^{01}+0.45N^{20}), & & \displaystyle\end{eqnarray}$$

respectively. The same procedure in the $(p,j)=(3,0)$ and $(p,j)=(1,1)$ moment equations yields

(5.7) $$\begin{eqnarray}\displaystyle 0\simeq \sqrt{\frac{3}{2}}N^{20}-\text{i}\unicode[STIX]{x1D708}(0.15N^{11}+1.03N^{30}), & & \displaystyle\end{eqnarray}$$

and

(5.8) $$\begin{eqnarray}\displaystyle 0\simeq \frac{N^{01}}{\sqrt{2}}-\text{i}\unicode[STIX]{x1D708}(1.09N^{11}+0.15N^{30}), & & \displaystyle\end{eqnarray}$$

respectively.

As a consequence, the truncated moment-hierarchy equations, equations (5.5)–(5.8), yield the following dispersion relation

(5.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}^{2}+1.96\unicode[STIX]{x1D6FE}\left(\frac{1}{\unicode[STIX]{x1D708}}+0.49\unicode[STIX]{x1D708}\right)+\frac{0.69}{\unicode[STIX]{x1D708}^{2}}+1.4\times 10^{-5}\unicode[STIX]{x1D708}^{2}+0.88\simeq 0 & & \displaystyle\end{eqnarray}$$

that up to second order in $\unicode[STIX]{x1D716}$ yields the solutions $\unicode[STIX]{x1D6FE}\simeq -0.96\unicode[STIX]{x1D708}-1.04/\unicode[STIX]{x1D708}$ and $\unicode[STIX]{x1D6FE}\simeq -0.92/\unicode[STIX]{x1D708}$ . The least damped solution, which is the one consistent with the ordering $\unicode[STIX]{x1D6FE}\sim 1/\unicode[STIX]{x1D708}$ in (5.1), when applied to the $\unicode[STIX]{x1D708}=5$ case of figure 7, leads to $\unicode[STIX]{x1D6FE}\simeq -0.18$ , which has a relative difference of 11 % with respect to the $\unicode[STIX]{x1D6FE}=-0.202$ value obtained numerically.

We note that, with the same ordering above, a purely damped solution can also be obtained from the one-dimensional linearized Braginskii equations (Braginskii Reference Braginskii1965). In this limit, in fact, the following linearized electron temperature equation is found

(5.10) $$\begin{eqnarray}\displaystyle n_{0}\frac{3}{2}\frac{\unicode[STIX]{x2202}T_{e}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}_{z}(-\unicode[STIX]{x1D712}_{\Vert }^{e}\unicode[STIX]{x1D735}_{z}T_{e})\simeq 0, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D712}_{\Vert }^{e}=3.2n_{0}T_{e}/(m_{e}\unicode[STIX]{x1D708}kv_{th})$ , with the Joule heating term proportional to $m_{e}/m_{i}$ neglected. Equation (5.10) yields the electron Braginskii entropy mode $\unicode[STIX]{x1D6FE}\simeq -1.1/\unicode[STIX]{x1D708}$ , a value that is close to the estimate above based on the truncated moment-hierarchy equation.

Finally, we remark that a purely damped mode is only observed in the temporal evolution of $N^{00}$ for values of $\unicode[STIX]{x1D708}\gtrsim 1$ , while for $\unicode[STIX]{x1D708}\lesssim 1$ a transition from damped oscillations to a purely damped behaviour is not seen to occur for the range of values of $\unicode[STIX]{x1D6FC}_{D}$ considered here even at later times. The value of $\unicode[STIX]{x1D708}$ where a transition from collisional Landau damping to a purely damped entropy mode occurs after an initial transient is visible in the time evolution of $N^{00}$ can be estimated by balancing the damping rate of the collisional damping of EPW with the damping rate of entropy modes. Estimating the former as $\unicode[STIX]{x1D6FE}\simeq -0.03-0.26\unicode[STIX]{x1D708}$ from figure 6, and the latter as $\unicode[STIX]{x1D6FE}\simeq -0.92/\unicode[STIX]{x1D708}$ , the collision frequency at which the transition occurs is therefore estimated to be $\unicode[STIX]{x1D708}\simeq 1.8$ , in agreement with the numerical results.

6 Eigenvalue spectrum

We now compute the eigenmode spectra of EPW, and highlight the differences between the spectra of the full Coulomb, electron–ion and Lenard–Bernstein operators using a Hermite–Laguerre decomposition, and the electron–ion operator using a Legendre polynomial decomposition. We note that subdominant and stable modes can be nonlinearly excited to finite amplitude (Terry et al. Reference Terry, Baver and Gupta2006; Hatch et al. Reference Hatch, Terry, Jenko, Merz and Nevins2011, Reference Hatch, Jenko, Navarro, Bratanov, Terry and Pueschel2016; Pueschel et al. Reference Pueschel, Faber, Citrin, Hegna, Terry and Hatch2016) and have a major role in nonlinear energy dissipation and turbulence saturation, affecting structure formation, as well as heat and particle transport. We note that both the Lenard–Bernstein and the Dougherty collision operators are seen to yield similar eigenmode spectra. Therefore, we do not consider the Dougherty operator for this analysis. To compute the Hermite–Laguerre EPW eigenmode spectrum, the moment-hierarchy equation, equation (2.19), is truncated at a maximum index $(P,J)$ , Fourier transformed in time and the resulting eigenvalue problem solved numerically, yielding the spectrum of solutions at arbitrary collisionality. In matrix form, this yields

(6.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63C}N=(\unicode[STIX]{x1D714}+\text{i}\unicode[STIX]{x1D6FE})N, & & \displaystyle\end{eqnarray}$$

where $N=[N^{00}N^{01}\cdots N^{0J}N^{10}N^{11}\cdots N^{pj}]$ is the moment vector and $\unicode[STIX]{x1D63C}$ the $(P+1)(J+1)\times (P+1)(J+1)$ matrix of moment-hierarchy coefficients of elements $\unicode[STIX]{x1D608}_{n}^{m}$ with $m$ and $n$ the row and column, respectively

(6.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D608}_{p^{\prime }J+j^{\prime }}^{pJ+j}=\sqrt{\frac{p+1}{2}}\unicode[STIX]{x1D6FF}_{p+1,p^{\prime }}\unicode[STIX]{x1D6FF}_{j,j^{\prime }}+\sqrt{\frac{p}{2}}\unicode[STIX]{x1D6FF}_{p-1,p^{\prime }}\unicode[STIX]{x1D6FF}_{j,j^{\prime }}+\frac{\unicode[STIX]{x1D6FF}_{p,1}\unicode[STIX]{x1D6FF}_{j,0}}{\sqrt{2}}\frac{\unicode[STIX]{x1D6FF}_{p^{\prime },0}\unicode[STIX]{x1D6FF}_{j^{\prime },0}}{\unicode[STIX]{x1D6FC}_{D}}+\text{i}C_{p^{\prime }J+j^{\prime }}^{pJ+j}, & & \displaystyle\end{eqnarray}$$

which can be written in matrix form as

(6.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63C}=\left[\begin{array}{@{}cccccc@{}}0 & 0 & \cdots \, & 1/\sqrt{2} & 0 & \cdots \\ 0 & \text{i}C_{01}^{01} & \cdots \, & \text{i}C_{10}^{01} & 1/\sqrt{2}+\text{i}C_{11}^{01} & \cdots \\ \vdots & \vdots & \cdots \, & \vdots & \vdots & \cdots \\ (1+1/\unicode[STIX]{x1D6FC}_{D})/\sqrt{2} & \text{i}C_{01}^{10} & \cdots \, & \text{i}C_{10}^{10} & \text{i}C_{11}^{10} & \cdots \\ 0 & 1\sqrt{2}+C_{01}^{11} & \cdots \, & \text{i}C_{10}^{11} & \text{i}C_{11}^{11} & \cdots \\ \vdots & \vdots & \cdots \, & \vdots & \vdots & \cdots \end{array}\right]. & & \displaystyle\end{eqnarray}$$

In (6.2) and (6.3), we have defined the collisional coefficients $C_{st}^{pj}$ in terms of the collisional moments $C^{pj}$ as $C^{pj}=\sum _{p^{\prime },j^{\prime }}C_{p^{\prime }j^{\prime }}^{pj}N^{p^{\prime }j^{\prime }}$ and $C_{p^{\prime }J+j^{\prime }}^{pJ+j}=C_{p^{\prime }j^{\prime }}^{pj}$ . The spectrum of $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D714}$ is then found by computing the eigenvalues of the matrix $\unicode[STIX]{x1D63C}$ .

The resulting eigenvalue spectrum for the Coulomb collision case is shown in figure 8 for $\unicode[STIX]{x1D708}=0.1$ (a) and $\unicode[STIX]{x1D708}=1$ (b), with $\unicode[STIX]{x1D6FC}_{D}=0.09$ and $(P,J)=(18,2)$ , together with the corresponding collisionless Landau root (red marker), i.e. the least-damped solution of (3.4). The resulting collisional spectrum is discrete, contrary to the continuous collisionless Van Kampen spectrum, as noted in previous studies of weakly collisional plasma systems (Ng, Bhattacharjee & Skiff Reference Ng, Bhattacharjee and Skiff1999; Bratanov et al. Reference Bratanov, Jenko, Hatch and Brunner2013). Figure 8 shows that the damping rate of the Coulomb eigenmodes decreases with the corresponding frequency, which is possibly related to the fact that the collisional drag force decreases with the particle velocity in the Coulomb collision operator. We also note that the least-damped Coulomb eigenvalue in figure 8 is not the one closest to the Landau collisionless solution, as there are modes with higher oscillation frequency $\unicode[STIX]{x1D714}$ that are less damped than the collisionless damping rate. These eigenvalue solutions, however, are related to eigenvectors that mainly involve moments $N^{pj}$ with large values of $p$ and $j$ , and have therefore a negligible contribution to the initial damping of $N^{00}$ and $\unicode[STIX]{x1D719}$ .

Finally, the Coulomb eigenmode spectrum in figure 8 includes modes with vanishing frequency and damping that increases with $\unicode[STIX]{x1D708}$ . These modes correspond therefore to purely damped modes with a damping rate that at low collisionalities can be comparable to the collisionless Landau one. These zero-frequency solutions have also been previously observed in the analysis of linear EPW when pitch-angle scattering effects are included (Epperlein et al. Reference Epperlein, Short and Simon1992; Banks et al. Reference Banks, Brunner, Berger and Tran2016), and correspond to the entropy mode studied in § 5.

As an aside, we note that when the moment-hierarchy equation, equation (2.19), is truncated at a higher $P$ , i.e. using a higher number of Hermite polynomials, the number of eigenmodes with high frequency and small damping rate increases. On the other hand, when the number of Laguerre polynomials, hence $J$ , is increased, the eigenmode spectrum present and increasing number of modes with similar frequencies but increasingly higher damping rates. However, as shown by figure 3, the damping rates $\unicode[STIX]{x1D6FE}$ closest to the collisionless solution have negligible variation when $P$ and $J$ are increased (for $P\geqslant 18$ and $J\geqslant 0$ the variation is smaller than 3 %).

Figure 8.

Complete eigenvalue spectrum of the truncated moment-hierarchy equation with $\unicode[STIX]{x1D6FC}_{D}=0.09$ and $(P,J)=(18,2)$ (yielding $(P+1)(J+1)=57$ eigenvalues), using the full Coulomb collision operator (a,b), the Lenard–Bernstein operator (c,d) and the electron–ion Coulomb operator only (e,f), with $\unicode[STIX]{x1D708}=0.1$ (a,c,e) and $\unicode[STIX]{x1D708}=1$ (b,d,f). The collisionless least damped solution is shown as a red marker, and the red vertical lines are the solutions of (6.6).

The eigenmode spectra using a Lenard–Bernstein collision operator are also shown in figure 8 for $\unicode[STIX]{x1D708}=0.1$ (c) and $\unicode[STIX]{x1D708}=1$ (d), with $\unicode[STIX]{x1D6FC}_{D}=0.09$ and $(P,J)=(18,2)$ . A clear difference is seen between the eigenmode spectra of the Coulomb and Lenard–Bernstein operators. Contrary to the Coulomb case, the damping rate of the EPW modes increases with the frequency $\unicode[STIX]{x1D714}$ when a Lenard–Bernstein collision operator is used. Also, contrary to the Coulomb case, the Lenard–Bernstein root closest to the Landau collisionless root is the least-damped one, as also noted in previous weakly collisional studies of EPW (Bratanov et al. Reference Bratanov, Jenko, Hatch and Brunner2013).

Finally, the eigenvalue spectrum using the electron–ion Coulomb operator introduced in (2.36) is shown in figure 8 for $\unicode[STIX]{x1D708}=0.1$ (e) and $\unicode[STIX]{x1D708}=1$ (f). The spectrum is qualitatively similar to the Coulomb one, with high frequency modes being less damped than modes with smaller oscillation frequency. As for the Coulomb collision operator, such frequency dependence may be due to the dependence of the drag force on the particle velocity. Indeed, the electron–ion collision operator contains a drag force that decreases with the particle velocity, similarly to the Coulomb operator.

We now estimate the frequency $\unicode[STIX]{x1D714}$ of the modes in figure 8 with a damping rate $\unicode[STIX]{x1D6FE}$ different than the ones closest to the collisionless roots, by noting that the values of $\unicode[STIX]{x1D714}$ in figure 8 are seen to be weakly dependent on $\unicode[STIX]{x1D708},\unicode[STIX]{x1D6FC}_{D}$ , and the collision operator for the range of values used. We therefore solve the moment-hierarchy equation, equation (2.19), in the $\unicode[STIX]{x1D719}=0$ limit, which effectively neglects the roots related to EPW. Furthermore, in order to retrieve purely oscillatory solutions, the collisional damping terms $C^{pj}$ in (2.19) are neglected. The time Fourier-transformed moment-hierarchy equation in the $\unicode[STIX]{x1D719}=C^{pj}=0$ limit reads

(6.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}N^{pj}=\sqrt{\frac{p+1}{2}}N^{p+1j}+\sqrt{\frac{p}{2}}N^{p-1j}. & & \displaystyle\end{eqnarray}$$

We recognize in (6.4) the recursion relation for the Hermite polynomials

(6.5) $$\begin{eqnarray}\displaystyle N^{pj}=\frac{H_{p}(\unicode[STIX]{x1D714})}{\sqrt{2^{p}p!}}. & & \displaystyle\end{eqnarray}$$

The roots $\unicode[STIX]{x1D714}$ can be found by applying the truncation condition in (4.1) to the solution in (6.5), yielding

(6.6) $$\begin{eqnarray}\displaystyle H_{P+1}(\unicode[STIX]{x1D714})=0. & & \displaystyle\end{eqnarray}$$

The solutions $\unicode[STIX]{x1D714}$ in (6.6) are purely real, yielding frequencies that closely follow the ones observed in the eigenvalue spectra (red vertical lines in figure 8).

Finally, we present two tests to assess the validity of the results in figure 8, first for the Lenard–Bernstein case and then for the electron–ion case. Focusing on the Lenard–Bernstein spectrum, we derive a polynomial in $\unicode[STIX]{x1D6FE}$ whose roots closely follow the modes in figure 8(c,d) that appear with damping rates larger than the ones of the two least-damped roots. Fourier transforming the Boltzmann equation, equation (2.10), in time and in velocity space similarly to Ng, Bhattacharjee & Skiff (Reference Ng, Bhattacharjee and Skiff2004), with $C(f)$ the Lenard–Bernstein collision operator, the following differential equation for $g(s)=\int _{-\infty }^{\infty }\exp (\text{i}sv_{z}-\unicode[STIX]{x1D6FE}t+\text{i}\unicode[STIX]{x1D714}t)f_{M}\unicode[STIX]{x1D6FF}f_{k}\,\text{d}v_{z}\,\text{d}t$ is obtained

(6.7) $$\begin{eqnarray}\displaystyle g(s)\left(\unicode[STIX]{x1D6FE}+\text{i}\unicode[STIX]{x1D714}+\frac{\unicode[STIX]{x1D708}}{2}s^{2}\right)+(1+\unicode[STIX]{x1D708}s)\frac{\text{d}g(s)}{\text{d}s}=-s\frac{\sqrt{\unicode[STIX]{x03C0}}}{2\unicode[STIX]{x1D6FC}_{D}}e^{-(s^{2}/4)}g(0). & & \displaystyle\end{eqnarray}$$

We solve (6.7) neglecting the coupling with the electrostatic potential $\unicode[STIX]{x1D719}$ by setting $\unicode[STIX]{x1D6FC}_{D}\gg 1$ (or, equivalently, setting $\unicode[STIX]{x1D719}=0$ in the Boltzmann equation), and define $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D708}^{-2}/2$ and $\unicode[STIX]{x1D6E4}=\sqrt{2\unicode[STIX]{x1D706}}(\unicode[STIX]{x1D6FE}+\text{i}\unicode[STIX]{x1D714})-\unicode[STIX]{x1D706}$ , yielding

(6.8) $$\begin{eqnarray}\displaystyle g(s)=g(0)\left(1+\frac{s}{2\unicode[STIX]{x1D706}}\right)^{\unicode[STIX]{x1D6E4}}e^{-(s^{2}/4)+s\unicode[STIX]{x1D706}}. & & \displaystyle\end{eqnarray}$$

Similarly, Fourier transforming the Hermite–Laguerre expansion of $\langle \unicode[STIX]{x1D6FF}f_{k}\rangle$ in velocity space, we obtain

(6.9) $$\begin{eqnarray}\displaystyle g(s)=\mathop{\sum }_{p=0}^{\infty }\frac{\text{i}^{p}N^{p0}}{\sqrt{2^{p+1}p!}}s^{p}e^{-(s^{2}/4)}. & & \displaystyle\end{eqnarray}$$

Equating the two expressions above, we find

(6.10) $$\begin{eqnarray}\displaystyle N^{p0}=N^{00}(-i)^{p}\sqrt{\frac{2^{p+1}}{p!}}\frac{\text{d}^{p}}{\text{d}s^{p}}\left[e^{s\unicode[STIX]{x1D706}}\left(1+\frac{s}{2\unicode[STIX]{x1D706}}\right)^{\unicode[STIX]{x1D6E4}}\right]_{s=0}. & & \displaystyle\end{eqnarray}$$

Therefore, the truncation condition in (4.1) in the $\unicode[STIX]{x1D719}=0$ limit is equivalent to imposing

(6.11) $$\begin{eqnarray}\displaystyle \frac{\text{d}^{P}}{\text{d}s^{P}}\left[e^{s\unicode[STIX]{x1D706}}\left(1+\frac{s}{2\unicode[STIX]{x1D706}}\right)^{\unicode[STIX]{x1D6E4}}\right]_{s=0}=0. & & \displaystyle\end{eqnarray}$$

Finally, we can rewrite (6.11) as a polynomial in $\unicode[STIX]{x1D6FE}+\text{i}\unicode[STIX]{x1D714}$

(6.12) $$\begin{eqnarray}\displaystyle \mathop{\sum }_{t=0}^{P}a_{Pt}(\unicode[STIX]{x1D706})(\unicode[STIX]{x1D6FE}+\text{i}\unicode[STIX]{x1D714})^{t}=0, & & \displaystyle\end{eqnarray}$$

with

(6.13) $$\begin{eqnarray}\displaystyle a_{pt}(\unicode[STIX]{x1D706})=\mathop{\sum }_{l=t}^{p}\mathop{\sum }_{n=l}^{p}\binom{p}{n}\binom{l}{t}s(n,l)\unicode[STIX]{x1D706}^{p+l-n-t/2}(-1)^{l-t}2^{t/2}, & & \displaystyle\end{eqnarray}$$

where $s(n,l)$ are the Stirling numbers of the first kind (Moser & Wyman Reference Moser and Wyman1958; Qi Reference Qi2014). As shown in figure 9, the polynomial expression in (6.12) closely reproduces the eigenvalue spectrum observed in figure 8(c,d). Therefore, although the coupling of the electron distribution function with $\unicode[STIX]{x1D719}$ is crucial to reproduce the EPW roots, additional modes in the eigenmode spectrum are related to solutions decoupled from the electrostatic potential $\unicode[STIX]{x1D719}$ , subject to the truncation condition of the Hermite–Laguerre series, equation (4.1), with frequencies similar to the ones of (6.6).

Figure 9.

Blue dots: roots of the polynomial in (6.12), which corresponds to the solution of the Boltzmann equation with a Lenard–Bernstein collision operator where the distribution function is approximated by a truncated Hermite expansion, for $\unicode[STIX]{x1D706}=50$ (a) and $\unicode[STIX]{x1D706}=0.5$ (b) (corresponding to $\unicode[STIX]{x1D708}=0.1$ and 1, respectively) at $P=20$ . Red vertical lines: solutions of (6.6).

As a second test, to assess the validity of the eigenmode spectrum found with an electron–ion collision operator in figure 8, we solve the Boltzmann equation using a different set of basis functions, namely expanding $\langle \unicode[STIX]{x1D6FF}f_{k}\rangle$ in Legendre polynomials, equation (2.40), and solving the resulting moment-hierarchy equation, equation (2.40), numerically. In this case, the expansion of $\langle \unicode[STIX]{x1D6FF}f_{k}\rangle$ in (2.39) is truncated at $l_{max}=L$ by setting $a_{L+1}=0$ . The velocity $v$ is discretized over an interval $[0,v_{max}]$ with an equidistant mesh made of $n_{v}$ points, and the integral estimated with a composite trapezoidal rule. The resulting spectrum is shown in figure 10. When compared with the Hermite–Laguerre spectrum in figure 8(e,f), the two spectra look qualitatively similar, confirming the validity of the Hermite–Laguerre approach. However, a higher number of small-damped low-frequency solutions is observed when a Legendre decomposition is used. The appearance of small-damped non-physical eigenmodes when using a finite-difference discretization in $v$ was also noted by Bratanov et al. (Reference Bratanov, Jenko, Hatch and Brunner2013), leading to the conclusion that, in general, a Hermite discretization of the distribution function is in fact superior to a finite-difference one. Furthermore, for the values of $\unicode[STIX]{x1D708}=0.02$ and $\unicode[STIX]{x1D6FC}_{D}=0.09$ where the Hermite–Laguerre formulation with $(P+1)(J+1)=19\times 3=57$ polynomials is seen to converge to the collisionless Landau solution in figure 6 with a relative difference of ${\sim}16\,\%$ , when using a Legendre decomposition, a total of $n_{v}\times L\simeq 100$ equations is needed to yield a similar accuracy on $\unicode[STIX]{x1D6FE}$ .

Figure 10.

Eigenvalue spectrum of the truncated moment-hierarchy equation using an electron–ion Coulomb collision operator for $\unicode[STIX]{x1D708}=0.1$ (a) and $\unicode[STIX]{x1D708}=1.0$ (b) with $\unicode[STIX]{x1D6FC}_{D}=0.09$ , $n_{v}=12$ and $L=7$ , with a Legendre decomposition. The collisionless solution is shown with a red marker.