2.1. LBO-G

In the same vein as was done for the BGK in Greene (Reference Greene1973), one may enforce exact momentum and energy conservation, and use Boltzmann relaxation rates (Morse Reference Morse1963) to obtain the cross-species primitive moments appropriate for Dougherty electron–ion collisions. Conservation of momentum and energy in cross-species collisions

(2.1) \begin{align} \int^\infty_{-\infty}v_i\sum_s m_s C[f_s]\,\text{d}^{d_v} v &= 0, \end{align}

(2.2) \begin{align}\int^\infty_{-\infty}\frac{1}{2}v^2\sum_s m_s C[f_s]\,\text{d}^{d_v} v &= 0 \end{align}

($C[f_s]$ the right-hand side of  (1.1) with $r\neq s$) is obeyed pairwise and yields the relations

(2.3) \begin{gather} \begin{gathered}\sum_s {m_s}n_s\nu_{sr}\Delta u_{{s}r,i} = 0,\\ \sum_s {m_s}n_s\nu_{sr}\left({d_v}\Delta v_{t,{s}r}^2+{ \boldsymbol{u} }_s\boldsymbol{\cdot}\Delta{ \boldsymbol{u} }_{{s}r}\right) = 0, \end{gathered} \end{gather}

with the sum running only over two species ($r$ labels the species other than $s$), $\Delta { \boldsymbol {u} }_{{s}r }={ \boldsymbol {u} }_s-{ \boldsymbol {u} }_{{s}r }$ and $\Delta v _{t,{s}r }^2=v_{t,s}^2-v_{t,{s}r }^2$. This system of $2({d_v}+1)$ unknowns can be closed in a number of ways; a particularly simple way is by employing the momentum and thermal relaxation rates of the full Coulomb collision operator (see (15) and (16) in Morse Reference Morse1963). For small-angle collisions these rates are

(2.4) \begin{align} \left.\begin{aligned}\frac{\partial }{\partial t}{m_s}n_{s}u_{s,i}\Big|_{\mathrm{FPO}} &= \frac{\alpha_{E}}{2}\left(m_s+m_r\right)\left(u_{r,i}-u_{s,i}\right), \\ \frac{\partial }{\partial t}\frac{{d_v}}{2}{m_s}n_sv_{t,s}^2\Big|_{\mathrm{FPO}} &= \frac{\alpha_{E}}{2}\left[{d_v}\left({m_r}v_{t,r}^2-m_sv_{t,s}^2\right)+m_r\left({ \boldsymbol{u} }_r-{ \boldsymbol{u} }_s\right)^2\right]. \end{aligned} \right. \end{align}

The parameter $\alpha _{E}$ is inversely proportional to the energy and momentum relaxation times

(2.5)\begin{equation} \alpha_{E} = \frac{2n_sn_r({q_s}q_r)^2\log\varLambda_{{s}r}}{3(2{\rm \pi})^{3/2}\epsilon_0^2m_sm_r(v_{t,s}^2+v_{t,r}^2)^{3/2}}. \end{equation}

The right-hand side of (2.4) originates from the Boltzmann collision integral for Coulomb interactions, truncated at the Debye length, under the premise that $f_s$ are close to Maxwellian. The validity of the relations given below in systems where the plasma may significantly differ from Maxwellian is thus limited.

One can compute the LBO momentum and thermal relaxation rates similar to those for the FPO in (2.4) simply by taking velocity moments of (1.1). These rates areFootnote ¹

(2.6)\begin{equation} \left.\begin{gathered} \frac{\partial }{\partial t}{m_s}n_{s}u_{s,i}\Big|_{\mathrm{LBO}} = {m_s}n_s \nu_{sr}(u_{{s}r,i}-u_{s,i}), \\ \frac{\partial }{\partial t}\frac{{d_v}}{2}{m_s}n_sv_{t,s}^2\Big|_{\mathrm{LBO}} = {d_v}{m_s}n_s\nu_{sr}(v_{t,{s}r}^2-v_{t,s}^2). \end{gathered}\right. \end{equation}

Equating (2.4) and (2.6) does not fully determine $u_{{s}r }$ and $v_{t,{s}r }^2$ because of the as-of-yet arbitrary $\nu _{sr}$. The next step in the Greene methodology is thus to adopt a relationship between the collision frequency in the model operator and $\alpha _{E}$, which for the BGK operator Greene took to be $\nu _{sr}=\alpha _{E}(m_s+m_r)/[(1+\beta )n_{s}m _s]$ with the arbitrary parameter $\beta >-1$. For the LBO-G we will instead use

(2.7)\begin{equation} \nu_{sr} = \frac{\alpha_{E}(m_s+m_r)}{\delta_s(1+\beta){m_s}n_s}, \end{equation}

with $\delta _s=2m_{r}n _r\nu _{rs}/({m_s}n _s\nu _{sr}+m_{r}n _r\nu _{rs})$; it turns out that $\delta _s$ and $\beta$ only appear as $\delta _s(1+\beta )$ so their independent values do not need to be determined separately. We picked this relationship between $\alpha _{E}$ and $\nu _{sr}$ for three reasons. First, we anticipate potential difficulties guaranteeing positivity of $v_{t,{s}r }^2$, although we will see shortly that such problems do not arise with the Dougherty operator for many systems of interest. Second, the formulation presented here avoids the assumption ${m_s}n _s\nu _{sr}=m_{r}n _r\nu _{rs}$ used in earlier work (Greene Reference Greene1973). Lastly, this definition of $\nu _{sr}$ produces relations that more easily enforce exact conservation in their discrete form.

Equipped with (2.7) we can equate (2.4) and (2.6), and together with (2.3) a linear system in $u_{{s}r,i}$, $v_{t,{s}r }^2$, $u_{{r}s,i}$ and $v_{t,rs}^2$ ensues. The solution of this linear problem is

(2.8) \begin{align} u_{{s}r,i} &= u_{s,i} + \delta_s\frac{1+\beta}{2}\left(u_{r,i}-u_{s,i}\right), \end{align}

(2.9) \begin{align}v_{t,{s}r}^2 &= v_{t,s}^2+\frac{\delta_s}{2}\dfrac{1+\beta}{1+\dfrac{m_s}{m_r}}\left[v_{t,r}^2-\dfrac{m_s}{m_r}v_{t,s}^2+\frac{\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}{{d_v}}\right]. \end{align}

One attractive property of these cross-species primitive moments is that, contrary to their BGK counterparts, they do not suffer from the pathology of negative $v_{t,ie}^2$ at supersonic values of the relative drift $|{ \boldsymbol {u} }_s-{ \boldsymbol {u} }_r|$. Positivity of (2.9) does require, however, that

(2.10)\begin{equation} \frac{\delta_s}{2}\frac{1+\beta}{1+\dfrac{m_r}{m_s}}\left[1-\frac{T_r}{T_s}-\frac{\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}{{{d_v}}c_{s,{s}r}^2}\right] < 1, \end{equation}

where $c_{s,{s}r }=\sqrt {T_s/m_r}$. This is true for any choice of $\delta _s$ and $\beta$ provided $\delta _s(1+\beta )<2$, even as the relative drift increases.

Despite such improvements on previous similar multispecies operators, the unspecified $\beta$ parameter poses a clear disadvantage. Dougherty & Watson (Reference Dougherty and Watson1967) had already pointed out that an additional condition is needed to determine all unknowns, and therefore avoid the appearance of any free parameters. As discussed by Haack et al. (Reference Haack, Hauck and Murillo2017), this free parameter can modify the transport coefficients in the associated fluid models. For the BGK operator, Morse (Reference Morse1964) eliminated the need for $\beta$ by assuming $n_s\nu _{sr}=n_r\nu _{rs}$ and requiring that the ratio of the relaxation rate for the momentum difference between the two species to that of the temperature difference be the same for both the FPO and the model operator. However, the resulting multispecies BGK operator does not satisfy the $H$-theorem, discouraging us from pursuing that approach. A possible added constraint that may do away with such parameter is the isotropization rate due to interspecies collisions; imposing such condition is, however, beyond the scope of this manuscript.

2.2. LBO-EM and LBO-ET

Following the path charted by Haack et al. (Reference Haack, Hauck and Murillo2017) for the BGK, one could require that $u_{{s}r }=u_{rs}$ and $m_sv_{t,{s}r }^2=m_{r}v _{t,rs}^2$. Then the momentum conservation constraint in (2.3) results in

(2.11)\begin{equation} u_{{s}r,i} = \frac{{m_s}n_s{\nu_{sr}}u_{s,i}+m_{r}n_r{\nu_{rs}}u_{r,i}}{{m_s}n_s\nu_{sr}+m_{r}n_r\nu_{rs}}, \end{equation}

while energy conservation assuming $m_sv_{t,{s}r }^2=m_{r}v _{t,rs}^2$ yields

(2.12)\begin{align} \left(n_s\nu_{sr}+n_r\nu_{rs}\right)m_sv_{t,{s}r}^2 & = {m_s}n_s\nu_{sr}v_{t,s}^2+m_{r}n_r{\nu_{rs}}v_{t,r}^2 \nonumber\\ & \quad+\frac{{m_s}n_s\nu_{sr}m_{r}n_r\nu_{rs}}{{m_s}n_s\nu_{sr}+m_{r}n_r\nu_{rs}}\frac{\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}{{d_v}}. \end{align}

The next step demands that the momentum relaxation rates are the same for both the LBO and the FPO. Setting the momentum relaxation rates equal to each other

(2.13) \begin{align} \left.\begin{aligned} \frac{\partial }{\partial t}\left({m_s}n_{s}u_{s,i}-m_r{n_r}u_{r,i}\right)\Big|_{\mathrm{FPO}} &= \frac{\partial }{\partial t}\left({m_s}n_{s}u_{s,i}-m_r{n_r}u_{r,i}\right)\Big|_{\mathrm{LBO}},\\ \alpha_{E}\left(m_s+m_r\right)\left(u_{r,i}-u_{s,i}\right) &= {m_s}n_s \nu_{sr}\left(u_{{s}r,i}-u_{s,i}\right) - m_{r}n_r \nu_{rs}\left(u_{{r}s,i}-u_{r,i}\right), \end{aligned}\right.\end{align}

and using (2.11) for $u_{{s}r }$ one obtains the relationship

(2.14)\begin{equation} \alpha_{E}\left(m_s+m_r\right) = \frac{2{m_s}n_s\nu_{sr}m_{r}n_r\nu_{rs}}{{m_s}n_s\nu_{sr}+m_{r}n_r\nu_{rs}}. \end{equation}

On the other hand, equivalence between thermal relaxation rates

(2.15) \begin{align} \frac{\partial}{\partial t}\frac{{d_v}}{2}\left({m_s}n_sv_{t,s}^2-m_r{n_r}v_{t,r}^2\right)\Big|_{\mathrm{FPO}} &= \frac{\partial }{\partial t}\frac{{d_v}}{2}\left({m_s}n_sv_{t,s}^2-m_r{n_r}v_{t,r}^2\right)\Big|_{\mathrm{LBO}},\nonumber\\ \alpha_{E}\left[{d_v}\left({m_r}v_{t,r}^2-m_sv_{t,s}^2\right)+\frac{m_r-m_s}{2}\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2\right] & = {d_v}\left[{m_s}n_s\nu_{sr}\left(v_{t,{s}r}^2-v_{t,s}^2\right)\right. \nonumber\\ & \left.\quad-m_{r}n_r\nu_{rs}\left(v_{t,rs}^2-v_{t,r}^2\right)\right], \end{align}

with the $v_{t,{s}r }^2$ from (2.12) implies that

(2.16)\begin{align} & \alpha_{E}\left[{m_r}v_{t,r}^2-m_sv_{t,s}^2+\frac{m_r-m_s}{2{d_v}}\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2\right]\nonumber\\ & \quad = \frac{2n_s\nu_{sr}n_r\nu_{rs}}{n_s\nu_{sr}+n_r\nu_{rs}}\left({m_r}v_{t,r}^2-m_sv_{t,s}^2\right) +\frac{n_s\nu_{sr}-n_r\nu_{rs}}{n_s\nu_{sr}+n_r\nu_{rs}} \frac{{m_s}n_s\nu_{sr}m_{r}n_r\nu_{rs}}{{m_s}n_s\nu_{sr}+m_{r}n_r\nu_{rs}}\frac{\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}{{d_v}}. \end{align}

Although they may look strongly nonlinear, one can solve (2.14) and (2.16) in order to obtain an expression for $\nu _{sr}$.Footnote ² The result is

(2.17)\begin{equation} \nu_{sr} = \frac{\alpha_{E}\left(m_s+m_r\right)}{n_{s}m_s}\cdot \frac{\dfrac{m_s-m_r}{2m_sm_r}{d_v}\left({m_r}v_{t,r}^2-m_sv_{t,s}^2\right)+ \left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}{\dfrac{1}{m_r}{d_v}\left({m_r}v_{t,r}^2-m_sv_{t,s}^2\right)+\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}, \end{equation}

but we can immediately notice that this can lead to negative collision frequencies in some parameter regimes; for example, the electron–ionfrequency $\nu _{ei}$ with zero relative drift. This indicates that enforcing the equality of momentum and thermal relaxation rates while using the assumptions $u_{{s}r,i}=u_{{r}s,i}$ and $m_sv_{t,{s}r }^2=m_{r}v _{t,rs}^2$ leads to unphysical behaviour. We nevertheless present two slight variations in the following subsections, as was also recently done by Ulbl et al. (Reference Ulbl, Michels and Jenko2021), in order to provide a point of reference for the LBO-G and comparing against Haack et al. (Reference Haack, Hauck and Murillo2017).

2.2.1. LBO-EM

Instead of trying to match both the momentum and thermal relaxation rates, we could satisfy ourselves with only attaining the same momentum relaxation rate. We can do this by employing (2.14), which we obtained from setting LBO and FPO momentum relaxation rates equal to each other, and further assuming that

(2.18)\begin{equation} {m_s}n_s\nu_{sr} = m_{r}n_r\nu_{rs}. \end{equation}

These two equations together set the collision frequency in our model to

(2.19)\begin{equation} \nu_{sr}^M = \alpha_{E}\frac{m_s+m_r}{{m_s}n_s} = \frac{2\left(m_s+m_r\right)({q_s}q_r)^2n_r\log\varLambda_{{s}r}}{3(2{\rm \pi})^{3/2}\epsilon_0^2m_s^2m_r\left(v_{t,s}^2+v_{t,r}^2\right)^{3/2}}. \end{equation}

This choice of collision frequency reduces the cross-primitive moments to

\begin{align} \left.\begin{aligned}u_{{s}r,i} &= \frac{u_{s,i}+u_{r,i}}{2}, \\ v_{t,{s}r}^2 &= \frac{1}{1+\dfrac{m_s}{m_r}}\left[v_{t,s}^2+v_{t,r}^2 + \frac{\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}{2{d_v}}\right]. \end{aligned}\right. \end{align}

We call (1.1) with collision frequency and cross-primitive moments in (2.19)–(2.20) the LBO-EM. Compared with the equations that led to LBO-G, (2.19) suggests that LBO-EM is LBO-G in the limit of $\beta =0$ and $\delta _s=1$. In this case the cross-species flow velocity in LBO-G (2.8) does reduce to that in LBO-EM, but the LBO-G cross-species thermal velocity in this limit does not equal that in (2.20). Interestingly, for vanishing relative drifts, $\beta =1$ leads to an agreement between $v_{t,{s}r }^2$ for LBO-G and LBO-EM, but leads to $u_{{s}r,i}=u_{r,i}$, which disagrees with LBO-EM's $u_{{s}r,i}$. Therefore, as with BGK, LBO-EM is not a special case of LBO-G.

2.2.2. LBO-ET

Alternatively, we could choose to approximately match the thermal relaxation rate of the FPO. Focusing on the temperature difference term in (2.4), we see that the relaxation rate due to temperature differences alone is the same for both species. We could choose to mimic this behaviour, and examining (2.6) the conclusion would be that we have to require

(2.21)\begin{equation} n_s\nu_{sr} = n_r\nu_{rs}. \end{equation}

This assumption renders (2.11)–(2.12) into

(2.22)\begin{gather} u_{{s}r,i} = \frac{{m_s}u_{s,i}+{m_r}u_{r,i}}{m_s+m_r}, \end{gather}

(2.23)\begin{gather}v_{t,{s}r}^2 = \frac{1}{2}\left[v_{t,s}^2+\frac{m_r}{m_s}v_{t,r}^2 + \frac{m_r}{m_s+m_r}\frac{\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}{{d_v}}\right]. \end{gather}

Although we took up relation (2.21) we have not specified the collision frequency precisely yet. We can do so by returning to the thermal relaxation rate equivalence (2.15) and inserting the cross-species temperature in (2.23). The result isFootnote ³

(2.24) \begin{align} \begin{aligned}& \alpha_{E}\left[{d_v}\left({m_r}v_{t,r}^2-m_sv_{t,s}^2\right)+\frac{m_r-m_s}{2}\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2 \right] \\ & \quad={d_v}\left[\frac{n_s\nu_{sr}}{2}\left({m_r}v_{t,r}^2-m_sv_{t,s}^2 + \frac{m_sm_r}{m_s+m_r} \frac{\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}{{d_v}}\right)\right.\\ & \qquad\left.- \frac{n_r\nu_{rs}}{2}\left(m_sv_{t,s}^2-{m_r}v_{t,r}^2 + \frac{m_sm_r}{m_s+m_r}\frac{\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}{{d_v}}\right)\right],\\ & \alpha_{E}\left[\left({m_r}v_{t,r}^2-m_sv_{t,s}^2\right)+\frac{m_r-m_s}{2}\frac{\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}{{d_v}}\right]\\ & \quad=\frac{n_s\nu_{sr}+n_r\nu_{rs}}{2}\left({m_r}v_{t,r}^2-m_sv_{t,s}^2 \right) + \frac{1}{2}\frac{m_sm_r}{m_s+m_r}\left(n_s\nu_{sr}-n_r\nu_{rs}\right)\frac{\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}{{d_v}}. \end{aligned}\end{align}

These thermal relaxation rates cannot agree exactly because, under the assumption $n_s\nu _{sr}=n_r\nu _{rs}$, the relative drift term vanishes for this LBO, but we could at least match the response due to the temperature difference, leading to the collision frequency for this model

(2.25)\begin{equation} \nu_{sr}^T = \frac{\alpha_{E}}{n_s} = \frac{2n_r({q_s}q_r)^2\log\varLambda_{{s}r}}{3(2{\rm \pi})^{3/2}\epsilon_0^2m_sm_r\left(v_{t,s}^2+v_{t,r}^2\right)^{3/2}}. \end{equation}

The operator (1.1) with this collision frequency and the cross-primitive moments in (2.22)–(2.23) is referred to as the LBO-ET. As with the previous operator, setting $\beta =0$ in the LBO-G leads to the same cross-species flow velocity $u_{{s}r,i}$, but then the thermal speeds $v_{t,{s}r }$ do not agree. More importantly, we restate that the LBO-ET did not exactly match the FPO thermal relaxation rate because of the difference in response to relative drifts. For plasmas where the relative drifts are small relative to temperature differences (not an uncommon situation), these rates agree exactly.

2.3. The $H$-theorem

The full FPO does not decrease entropy, i.e. it satisfies an $H$-theorem, and as a model FPO it is desirable that this formulation of multispecies Dougherty collisions retains such a property. The original paper on the multispecies Dougherty operator demonstrated a non-decreasing entropy only to second order after linearization (Dougherty & Watson Reference Dougherty and Watson1967), hinting at the possibility that the above full-$f$ equivalent operator could possess an $H$-theorem. It is in fact possible, however, to show that the Dougherty models for multispecies full-$f$ collisions presented here do have an $H$-theorem, even for species with unequal temperatures. A more detailed proof of this statement is given in Appendix B, and we give an outline of the argument here.

The total entropy $\mathcal {S}=-\sum _s\int \textrm {d}^{d_v} v\,f_s\ln f _s$ can be shown to obey

(2.26)\begin{equation} \dot{\mathcal{S}} = \frac{\partial \mathcal{S}}{\partial t} ={-}\sum_s\int\text{d}^{d_v} v~\nu_{sr}\left(\ln f_s+1\right)\nabla_{\boldsymbol{v}}\cdot{{ \boldsymbol{J_{{s}r}} }}, \end{equation}

where the flux ${ \boldsymbol {J_{{s}r }} }$ is the term in square brackets in (1.1). Using the definition of this flux and integration by parts twice together with the fact that $f_s\to 0$ faster than powers of $v_i$ as $v_i\to \pm \infty$ one is led to

(2.27)\begin{equation} \dot{\mathcal{S}} = \sum_s\nu_{sr}\left(-{{d_v}}n_s+v_{t,{s}r}^2\int\text{d}^{d_v} v\,\nabla_{\boldsymbol{v}}{f_s}\boldsymbol{\cdot}\nabla_{\boldsymbol{v}}{\ln f_s}\right). \end{equation}

At this point we can perform a variational minimization of this functional in order to determine if that minimum is below zero (indicating a violation of thermodynamic law). For a given set of primitive moments ($n_s$, $u_{s,i}$, $v_{t,s}$, $u_{{s}r,i}$, $v_{t,{s}r }$) and the virtual displacement $\delta f _s=f_s-f_{s0}$, the response of $\dot {\mathcal {S}}$ is

(2.28)\begin{equation} \delta\dot{\mathcal{S}} = \sum_sv_{t,{s}r}^2\int\text{d}^{d_v} v\,\nabla_{\boldsymbol{v}}{f_{s0}}\boldsymbol{\cdot}\left[\frac{2}{\delta f_s}\nabla_{\boldsymbol{v}}{\delta f_s}-\frac{1}{f_{s0}} \nabla_{\boldsymbol{v}}{f_{s0}}\right]\frac{\delta f_s}{f_{s0}}. \end{equation}

At an extremum in $\dot {\mathcal {S}}$ this function must vanish, and since (2.27) has no upper bound this extremum must be a minimum. We are also interested in virtual displacements that do not alter the moments of each distribution, that is

(2.29)\begin{equation} \int\text{d}^{d_v} v\,{ \boldsymbol{v} }^k\delta f_s =0 \quad \text{for}\ k\in\{0,1,2\}. \end{equation}

Further imposing that the displacement $\delta f _s$ vanishes at infinity, $\delta \dot {\mathcal {S}}=0$ and (2.28)–(2.29) yield the nonlinear inhomogeneous equation for the minimizing distribution $f_{s0}$

(2.30)\begin{equation} \left|\nabla_{\boldsymbol{v}}{\ln f_{s0}}\right|^2+2\nabla^2_{\boldsymbol{v}}{\ln f_{s0}} = h_0^2+2{d_v} h_1+2h_0{ \boldsymbol{h_1} }\boldsymbol{\cdot}{ \boldsymbol{v} }+h_1^2v^2, \end{equation}

with $h_0$, ${ \boldsymbol {h_1} }$ and $h_2$ undetermined linear coefficients. The solution to this equation is $f_{s0}\propto \exp (h_{0,i}v_i+h_1v^2/2)$. Enforcing the condition that it has the same number density $n_s$ and primitive moments ($u_{s,i}$ and $v_{t,s}$) as the original distribution, $f_s$, reveals that the distribution that minimizes the rate of entropy change of this operator is a Maxwellian with $n_s$, $u_{s,i}$ and $v_{t,s}$. The final step is to insert this distribution back into our expression for entropy change, (2.26), and check that the minimum entropy rate of change does not fall below zero. Such a procedure results in

(2.31)\begin{equation} \min\left(\frac{\partial \mathcal{S}}{\partial t}\right) = {d_v} \sum_s \frac{n_s\nu_{sr}}{v_{t,s}^2}\left(v_{t,{s}r}^2-v_{t,s}^2\right), \end{equation}

and since the procedure obtained a single minimum it must be the global minimum.

2.3.1. LBO-G $H$-theorem

Use the definition of $v_{t,{s}r }^2$ for the LBO-G model (2.9) to arrive at

(2.32)\begin{align} \min\left(\frac{\partial \mathcal{S}}{\partial t}\right) & = \frac{\delta_s{m_s}n_s\nu_{sr}}{2v_{t,s}^2}\frac{1+\beta}{m_s+m_r}\left[{d_v}\frac{\left(T_r-T_s\right)^2}{T_{s}T_r} \right.\nonumber\\ & \left.\quad+\left(\frac{m_r}{m_s}\frac{1}{v_{t,s}^2}+\frac{m_s}{m_r}\frac{1}{v_{t,r}^2}\right)\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2\right] \geqslant 0. \end{align}

We are thus led to the conclusion that the LBO-G model of full-$f$ multispecies collisions does not decrease the entropy. This is in contrast to the BGK-G operator, for which the $H$-theorem could not be proven or disproven (Greene Reference Greene1973; Haack et al. Reference Haack, Hauck and Murillo2017).

2.3.2. LBO-EM and LBO-ET $H$-theorem

Using the relationship between collision frequencies for the LBO-EM (2.19) and the corresponding cross-species thermal speeds one obtains

(2.33)\begin{align} \min\left(\frac{\partial \mathcal{S}}{\partial t}\right) = \frac{{{d_v}}n_s\nu_{sr}}{m_s+m_r}\left[ \frac{\left({m_r}v_{t,r}^2-m_sv_{t,s}^2\right)^2}{v_{t,s}^2{m_r}v_{t,r}^2} + \left(\frac{m_r}{v_{t,s}^2}+\dfrac{m_s}{m_r}\frac{m_s+m_r}{v_{t,r}^2}\right)\frac{\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}{2{d_v}}\right] \geqslant 0. \end{align}

Similarly, using the $n_s\nu _{sr}=n_r\nu _{rs}$ assumption of the LBO-ET model and the definition of $v_{t,{s}r }^2$ in (2.23) turns (2.31) into

(2.34)\begin{align} \min\left(\frac{\partial \mathcal{S}}{\partial t}\right) = \frac{{{d_v}}n_s\nu_{sr}}{2}\left[ \frac{{m_r}v_{t,r}^2}{m_sv_{t,s}^2} + \frac{m_sv_{t,s}^2}{{m_r}v_{t,r}^2} + \left(\frac{m_r}{v_{t,s}^2}+\frac{m_s}{v_{t,r}^2}\right)\frac{1}{m_s+m_r}\frac{\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2}{{d_v}}\right] \geqslant 0. \end{align}

Therefore, both LBO-EM and LBO-ET models satisfy the $H$-theorem.

2.4. Gyroaveraged multispecies Dougherty operator

This model operator is also used by modern, long-wavelength full-$f$ gyrokinetic codes (Pan et al. Reference Pan, Told, Shi, Hammett and Jenko2018; Gkeyll 2020) in its gyroaveraged form. Its form, conservative properties and discontinuous Galerkin discretization for self-species collisions have been presented by Francisquez et al. (Reference Francisquez, Bernard, Mandell, Hammett and Hakim2020). The operator can, however, be extended to incorporate cross-species collisions. For that purpose, we write the gyroaveraged operator acting on the guiding centre distribution function $f_s({ \boldsymbol {R} },v_\parallel,\mu )$ as

(2.35)\begin{align} \left(\frac{\partial \mathcal{J} f_s}{\partial t}\right)_c & = \sum_r\nu_{sr} \left\{\frac{\partial }{\partial v_\parallel}\left[\left(v_\parallel{-}u_{{\parallel}{s}r}\right)\mathcal{J}f_s+v_{t,{s}r}^2\frac{\partial \mathcal{J}f_s}{\partial v_\parallel}\right] \right.\nonumber\\ & \left.\quad+\frac{\partial }{\partial \mu}2\mu\left[\mathcal{J}f_s+\frac{m_sv_{t,{s}r}^2}{B}\frac{\partial \mathcal{J}f_s}{\partial \mu}\right]\right\}, \end{align}

where $\mathcal {J}$ represents the Jacobian of the guiding centre coordinates, ${ \boldsymbol {R} }$ is the guiding centre position, $v_\parallel$ is the velocity along the background magnetic field and $\mu$ is the adiabatic moment; see Francisquez et al. (Reference Francisquez, Bernard, Mandell, Hammett and Hakim2020) for more details.

In order to use this multispecies gyroaveraged operator one must then determine the multispecies parallel flow velocities $u_{\parallel {s}r }$ and thermal speed $v_{t,{s}r }$. Our proposal is to use the same LBO-G (2.8), (2.9), LBO-EM (2.20) and LBO-ET (2.22)–(2.23) models with this gyroaveraged operator. The only difference is that the self-species primitive moments are defined by

(2.36)\begin{equation} \left.\begin{gathered} {u_{\parallel{s}}}M_{0,s} = M_{1,s}, \\ {u_{\parallel{s}}}M_{1,s} + {d_v}v_{t,s}^2M_{0,s} = M_{2,s}, \end{gathered}\right. \end{equation}

where ${d_v}=1$ or ${d_v}=3$ depending on whether one is considering $v_\parallel$- or $(v_\parallel,\mu )$-space, respectively. The velocity moments in the gyroaveraged model are

(2.37) \begin{align} \left.\begin{aligned} M_{0,s} &= (2{\rm \pi}/m_s)\int\mathcal{J}f_s({ \boldsymbol{R} },v_\parallel,\mu)\,\mathrm{d}v_\parallel\,\mathrm{d}\mu, \\ M_{1\parallel,s} &= (2{\rm \pi}/m_s)\int v_\parallel\mathcal{J}f_s({ \boldsymbol{R} },v_\parallel,\mu)\,\mathrm{d}v_\parallel\,\mathrm{d}\mu, \\ M_{2,s} &= (2{\rm \pi}/m_s)\int\left(v_\parallel^2+2\mu B/m_s\right)\mathcal{J}f_s({ \boldsymbol{R} },v_\parallel,\mu)\,\mathrm{d}v_\parallel\,\mathrm{d}\mu. \end{aligned}\right. \end{align}

3.1. Discrete momentum conservation

In order to formulate a momentum-conserving discretization based on (3.1), we can set $\psi _\ell =m_s v$ and sum over all cells along velocity space. According to (2.1) this sum has to be equal and opposite to that of the other species it is colliding with. Therefore discrete momentum conservation requires that

(3.3)\begin{align} & \sum_k\int_{x_{j-1/2}}^{x_{j+1/2}}\left(m_s\nu_{sr}\left\{\left(v G_s-v_{t,{s}r}^2\,\hat{f_s}\right)_{v_{k-1/2}}^{v_{k+1/2}} - \int_{v_{k-1/2}}^{v_{k+1/2}} \left(v-u_{{s}r}\right)f_s\,\mathrm{d}v\right\} \right.\nonumber\\ & \quad\left.+ m_r\nu_{rs}\left\{\left(v G_r-v_{t,rs}^2\,\hat{f_r}\right)_{v_{k-1/2}}^{v_{k+1/2}} - \int_{v_{k-1/2}}^{v_{k+1/2}} \left(v-u_{rs}\right)f_r\,\mathrm{d}v\right\}\right){\rm d} x = 0. \end{align}

Carry out the velocity-space integrals and sum over all velocity-space cells. Use the fact that the numerical fluxes $G_s$ are continuous and have opposite signs on either side of a cell boundary, and that $\hat {f_s}$ is continuous across cell boundaries as well. Furthermore, we impose the zero-flux boundary conditions $G_s(v=v_{\max })=G_s(v=v_{\min })=0$ in order to arrive at

(3.4)\begin{align} & \int_{x_{j-1/2}}^{x_{j+1/2}}\left[m_s\nu_{sr}\left(v_{t,{s}r}^2\,f_s\Big|_{v_{\min}}^{v_{\max}} + M_{1,s} - u_{{s}r}M_{0,s}\right) \right.\nonumber\\ & \quad\left.+ m_r\nu_{rs}\left(v_{t,rs}^2\,f_r\Big|_{v_{\min}}^{v_{\max}} + M_{1,r} - {u_{rs}}M_{0,r}\right)\right]{{\rm d} x} = 0, \end{align}

having substituted the discrete form of the velocity moments (3.2). This relation is satisfied if

(3.5)\begin{equation} m_s\nu_{sr}\left(v_{t,{s}r}^2\,f_s\Big|_{v_{\min}}^{v_{\max}} + M_{1,s} - u_{{s}r}M_{0,s}\right) + m_r\nu_{rs}\left(v_{t,rs}^2\,f_r\Big|_{v_{\min}}^{v_{\max}} + M_{1,r} - {u_{rs}}M_{0,r}\right) \doteq 0. \end{equation}

Note that we have used the same $v_{\max }$ and $v_{\min }$ for both species for pedagogical reasons only. In fact since (3.5) is a constraint on position-space fields only, the discrete velocity space of each species can be completely different. For completeness we state that in the ${d_v}$-dimensional case the condition for the scheme in (3.1) to conserve momentum is that the cross-species primitive moments $u_{{s}r }$ and $v_{t,{s}r }$ must satisfy

(3.6)\begin{align} & m_s\nu_{sr}\left(M_{1i,s} - {u_{sr,i}}M_{0,s} + v_{t,{s}r}^2\int f_s\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i\right) \nonumber\\ & \quad+ m_r\nu_{rs}\left(M_{1i,r} - u_{{r}s,i}M_{0,r} + v_{t,rs}^2\int f_r\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i\right) \doteq 0, \end{align}

using $\int \mathrm {d}S_i$ as an integral over the velocity-space boundaries orthogonal to the $i$th velocity dimension, and the repeated index $i$ implies summation.

3.2. Discrete energy conservation

Energy conservation will impose a secondary constraint on how the discrete cross-species primitive moments must be computed. In order to obtain such condition we substitute $\psi _\ell =m_s v^2/2$ in (3.1), and sum over velocity-space cells and species. This action leads to

(3.7)\begin{align} & \sum_k\int_{x_{j-1/2}}^{x_{j+1/2}}\left(m_s\nu_{sr}\left\{\left.\left(\frac{v^2}{2}G_s-vv_{t,{s}r}^2\,\hat{f_s}\right)\right|_{v_{k-1/2}}^{v_{k+1/2}} - \int_{v_{k-1/2}}^{v_{k+1/2}} \left[v\left(v-u_{{s}r}\right)f_s-v_{t,{s}r}^2\,f_s\right]\mathrm{d}v\right\}\right.\nonumber\\ & \quad \left.+m_r\nu_{rs}\left\{\left.\left(\frac{v^2}{2}G_r-vv_{t,rs}^2\,\hat{f_r}\right)\right|_{v_{k-1/2}}^{v_{k+1/2}} - \int_{v_{k-1/2}}^{v_{k+1/2}} \left[v\left(v-u_{rs}\right)f_r-v_{t,rs}^2\,f_r\right]\mathrm{d}v\right\}\right){{\rm d} x} = 0. \end{align}

Once again, we employ continuity of $G_s$ and $\hat {f_s}$ and boundary conditions so that after performing the velocity-space integrals and carrying out the sum over velocity-space cells this relation is transformed into

(3.8)\begin{align} & \int_{x_{j-1/2}}^{x_{j+1/2}}\left(m_s\nu_{sr}\left\{vv_{t,{s}r}^2\,f_s\Big|_{v_{\min}}^{v_{\max}} + \left(M_{2,s}-{u_{sr}}M_{1,s}\right)-v_{t,{s}r}^2M_{0,s}\right\} \right.\nonumber\\ & \left.\quad+m_r\nu_{rs}\left\{ vv_{t,rs}^2\,f_r\Big|_{v_{\min}}^{v_{\max}} + \left(M_{2,r}-u_{rs}M_{1,r}\right)-v_{t,rs}^2M_{0,r}\right\}\right){{\rm d} x} = 0. \end{align}

Therefore, we can guarantee that our DG discretization exactly conserves energy if we enforce

(3.9)\begin{align} & m_s\nu_{sr}\left\{ vv_{t,{s}r}^2\,f_s\Big|_{v_{\min}}^{v_{\max}} + \left(M_{2,s}-{u_{sr}}M_{1,s}\right)-v_{t,{s}r}^2M_{0,s}\right\} \nonumber\\ & \quad+m_r\nu_{rs}\left\{ vv_{t,rs}^2\,f_r\Big|_{v_{\min}}^{v_{\max}} + \left(M_{2,r}-u_{rs}M_{1,r}\right)-v_{t,rs}^2M_{0,r}\right\} \doteq 0, \end{align}

when computing $u_{{s}r }$, $u_{rs}$, $v_{t,{s}r }$ and $v_{t,rs}$. In the case of ${d_v}$ velocity dimensions this constraint becomes

(3.10)\begin{align} & m_s\nu_{sr}\left[M_{2,s}-{u_{sr,i}}M_{1i,s}-v_{t,{s}r}^2\left({d_v}M_{0,s}-\int v_if_s\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i\right)\right]\nonumber\\ & \quad+m_r\nu_{rs}\left[M_{2,r}-{u_{rs}}iM_{1i,r}-v_{t,rs}^2\left({d_v}M_{0,r}-\int v_if_r\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i\right)\right] \doteq 0. \end{align}

We make a final comment that the substitution $\psi _\ell =m_s v^2/2$ is only valid if $v^2$ belongs to the space spanned by the basis, which for piecewise-linear basis functions ($p=1$) it does not. In order for the algorithm to be conservative with $p=1$ additional precautions must be taken, a topic that is deferred to Appendix C.

3.3. Discrete relaxation rates

Together with the relations $u_{{s}r,i}=u_{{r}s,i}$ and $m_sv_{t,{s}r }^2=m_{r}v _{t,rs}^2$, and the definitions of the collision frequency given in §§ 2.2.1–2.2.2, (3.6) and (3.10) are all one needs to compute the cross-primitive moments for the LBO-EM and LBO-ET. The LBO-G, however, needs to further incorporate the equivalence between the momentum and thermal relaxation rates of the FPO (2.4) and those of the LBO (2.6) in the discrete sense. First, we obtain the weak form of (2.13) by projecting it onto the $\psi _\ell$ basis function

(3.11)\begin{align} & m_s\nu_{sr}\left(M_{1i,s} - {u_{sr,i}}M_{0,s} + v_{t,{s}r}^2\int f_s\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i\right) \nonumber\\ & \quad- m_r\nu_{rs}\left(M_{1i,r} - u_{{r}s,i}M_{0,r} + v_{t,rs}^2\int f_r\Big|_{v_{i,\min}}^{v_{i,\max}}\mathrm{d}S_i\right) \doteq \alpha_{E}\left(m_s+m_r\right)\left(u_{s,i}-u_{r,i}\right), \end{align}

where we used the same series of steps that led to (3.6).

The equivalent condition on the thermal relaxation rates necessitates the discrete thermal speed moment of the LBO. We can obtain it by substituting $\psi _\ell =m_s(v-u_{s})^2/2$ into the weak scheme in (3.1) and summing over velocity-space cells, resulting in

(3.12)\begin{align} & \sum_k\int_{K_{j,k}}\frac{m_s}{2}\left(v-u_{s}\right)^2\left(\frac{{\rm d} f_s}{{\rm d} t}\right)_c\,\mathrm{d}x\,\mathrm{d}v \nonumber\\ & \quad = \sum_k\int_{x_{j-1/2}}^{x_{j+1/2}}m_s\nu_{sr}\left\{\left.\left[\frac{1}{2}\left(v-u_{s}\right)^2G_s -\left(v-u_{s}\right)v_{t,{s}r}^2\,\hat{f_s}\right]\right|_{v_{k-1/2}}^{v_{k+1/2}}\right.\nonumber\\ & \qquad \left.- \int_{v_{k-1/2}}^{v_{k+1/2}} \left[\left(v-u_{s}\right)\left(v-u_{{s}r}\right)-v_{t,{s}r}^2\right]f_s\,\mathrm{d}v\right\}{{\rm d} x}. \end{align}

Performing the velocity-space integrals, carrying out the $k$-sum, accounting for the continuity of $G_s$ and $\hat {f_s}$ and using the zero-flux boundary conditions on $G_s$ ushers us to

(3.13)\begin{align} & \sum_k\int_{K_{j,k}}\left(\frac{{\rm d} \left({m_s}n_sv_{t,s}^2/2\right)}{{\rm d} t}\right)_c\,\mathrm{d}x\,\mathrm{d}v \nonumber\\ & \quad = \int_{x_{j-1/2}}^{x_{j+1/2}}m_s\nu_{sr}\left\{\left[- \left(v-u_{s}\right)v_{t,{s}r}^2\,f_s\right]\Big|_{v_{\min}}^{v_{\max}} \right.\nonumber\\ & \left.\qquad- \vphantom{\left[- \left(v-u_{s}\right)v_{t,{s}r}^2\,f_s\right]\Big|_{v_{\min}}^{v_{\max}}}\left[M_{2,s}-{u_{s}}M_{1,s}+u_{{s}r}\left(u_{s}M_{0,s}-M_{1,s}\right)-v_{t,{s}r}^2M_{0,s}\right]\right\}{{\rm d} x}, \end{align}

or in the ${d_v}$-dimensional velocity space

(3.14)\begin{align} & \sum_k\int_{K_{j,k}}\left(\frac{{\rm d} \left({m_s}n_sv_{t,s}^2/2\right)}{{\rm d} t}\right)_c\,\mathrm{d}x\,\mathrm{d}v \nonumber\\ & \quad = \int_{x_{j-1/2}}^{x_{j+1/2}}m_s\nu_{sr}\left\{v_{t,{s}r}^2\int\left[- \left(v_i-u_{s,i}\right)f_s\right]\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i \right.\nonumber\\ & \qquad-\left.\vphantom{\int\left[- \left(v_i-u_{s,i}\right)f_s\right]}\left[M_{2,s}-u_{s,i}M_{1i,s}+u_{{s}r,i}\left(u_{s,i}M_{0,s}-M_{1i,s}\right)-v_{t,{s}r}^2{d_v}M_{0,s}\right]\right\}{{\rm d} x}, \end{align}

where once again the repeated index $i$ implies summation. Equipped with this formula we can write down the discrete equivalence between thermal relaxation rates (2.15) as

(3.15)\begin{align} & m_s\nu_{sr}\left\{v_{t,{s}r}^2\int\left[- \left(v_i-u_{s,i}\right)f_s\right]\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i \right.\nonumber\\ & \left.\qquad- \vphantom{\int\left[- \left(v_i-u_{s,i}\right)f_s\right]\Big|_{v_{i,\min}}^{v_{i,\max}}} \left[M_{2,s}-u_{s,i}M_{1i,s}+u_{{s}r,i}\left(u_{s,i}M_{0,s}-M_{1i,s}\right)-v_{t,{s}r}^2{d_v}M_{0,s}\right]\right\} \nonumber\\ & \qquad- m_r\nu_{rs}\left\{v_{t,rs}^2\int\left[- \left(v_i-u_{r,i}\right)f_r\right]\Big|_{v_{i,\min}}^{v_{i,\max}}\, \mathrm{d}S_i \right.\nonumber\\ & \left.\qquad- \vphantom{\int\left[- \left(v_i-u_{r,i}\right)f_r\right]\Big|_{v_{i,\min}}^{v_{i,\max}}} \left[M_{2,r}-u_{r,i}M_{1i,r}+u_{{r}s,i}\left(u_{r,i}M_{0,r}-M_{1i,r}\right)-v_{t,rs}^2{d_v}M_{0,r}\right]\right\} \nonumber\\ & \quad\doteq \alpha_{E} \left[{d_v}\left({m_r}v_{t,r}^2-m_sv_{t,s}^2\right)+\frac{m_r-m_s}{2}\left({ \boldsymbol{u} }_s-{ \boldsymbol{u} }_r\right)^2\right]. \end{align}

The two discrete relaxation rate equivalences in (3.11) and 3.15 in conjunction with the discrete momentum and energy conservation constraints (3.6) and (3.10) provide the four equations for the calculation of the $u_{{s}r,i}$, $u_{{r}s,i}$, $v_{t,{s}r }$, $v_{t,rs}$ unknowns in the LBO-G, provided a value of $\beta$. Equations (3.11) and (3.15), however, are written in terms of the self-primitive moments (e.g. $u_{s,i}$, $v_{t,s}$) of each species, imbuing such equations with some ambiguity as to whether the calculation of self-primitive moments should include the corrections from velocity-space boundaries or not (Hakim et al. Reference Hakim, Francisquez, Juno and Hammett2020; Francisquez et al. Reference Francisquez, Bernard, Mandell, Hammett and Hakim2020). We therefore opt to instead write those relations in terms of the velocity moments as follows:

(3.17) \begin{align} \begin{aligned} &m_s\nu_{sr}\left(M_{1i,s} - {u_{sr,i}}M_{0,s} + v_{t,{s}r}^2\int f_s\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i\right) \nonumber\\ &\qquad - m_r\nu_{rs}\left(M_{1i,r} - u_{{r}s,i}M_{0,r} + v_{t,rs}^2\int f_r\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i\right) \nonumber\\ &\quad\doteq \frac{\alpha_{E}\left(m_s+m_r\right)}{{M_{0,s}}M_{0,r}}\left({M_{0,r}}M_{1i,s}-{M_{0,s}}M_{1i,r}\right),\nonumber\\ &m_s\nu_{sr}\left\{{-}v_{t,{s}r}^2\int \left(v_i-u_{s,i}\right)f_s\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i \right.\nonumber\\ &\qquad\left. - \vphantom{\int \left(v_i-u_{s,i}\right)f_s\Big|_{v_{i,\min}}^{v_{i,\max}}}\left[M_{2,s}-u_{s,i}M_{1i,s}+u_{{s}r,i}\left(u_{s,i}M_{0,s}-M_{1i,s}\right)-v_{t,{s}r}^2{d_v}M_{0,s}\right]\right\} \nonumber\\ &\qquad - m_r\nu_{rs}\left\{{-}v_{t,rs}^2\int \left(v_i-u_{r,i}\right)f_r\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i \right.\nonumber\\ &\qquad\left. -\vphantom{\int \left(v_i-u_{r,i}\right)f_r\Big|_{v_{i,\min}}^{v_{i,\max}}} \left[M_{2,r}-u_{r,i}M_{1i,r}+u_{{r}s,i}\left(u_{r,i}M_{0,r}-M_{1i,r}\right)-v_{t,rs}^2{d_v}M_{0,r}\right]\right\} \nonumber\\ &\quad\doteq \frac{\alpha_{E}}{{M_{0,s}}M_{0,r}}\left[\vphantom{\frac{m_r-m_s}{2}}m_{r}M_{0,s}\left(M_{2,r}-u_{r,i}M_{1i,r}\right)-m_sM_{0,r}\left(M_{2,s}-u_{s,i}M_{1i,s}\right) \right.\nonumber\\ &\qquad\left.+\frac{m_r-m_s}{2}\left(u_{s,i}-u_{r,i}\right)\left({M_{0,r}}M_{1i,s}-{M_{0,s}}M_{1i,r}\right)\right]. \end{aligned}\end{align}

The division by ${M_{0,s}}M _{0,r}$ on the right side of these equations is to be performed weakly (Hakim et al. Reference Hakim, Francisquez, Juno and Hammett2020) in order to avoid aliasing errors.

3.4. Summary of discrete equations

In summary, in order for the algorithm based on (3.1) to conserve momentum and energy the discrete cross-primitive moments are computed using the conservation constraints

(3.19)\begin{align} \begin{aligned} & m_s\nu_{sr}\left({u_{sr,i}}M_{0,s} - v_{t,{s}r}^2\int f_s\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i\right) \nonumber\\ &\quad + m_r\nu_{rs}\left(u_{{r}s,i}M_{0,r} - v_{t,rs}^2\int f_r\Big|_{v_{i,\min}}^{v_{i,\max}}\, \mathrm{d}S_i\right) \doteq m_s\nu_{sr}M_{1i,s} + m_r{\nu_{rs}}M_{1i,r},\nonumber\\ & m_s\nu_{sr}\left[{u_{sr,i}}M_{1i,s} + v_{t,{s}r}^2\left({d_v}M_{0,s}-\int v_if_s\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i\right)\right]\nonumber\\ &\quad +m_r\nu_{rs}\left[{u_{rs}}iM_{1i,r} + v_{t,rs}^2 \left({d_v}M_{0,r}-\int v_if_r\Big|_{v_{i,\min}}^{v_{i,\max}}\, \mathrm{d}S_i\right)\right] \doteq m_s\nu_{sr}M_{2,s} + m_r\nu_{rs}M_{2,r}. \end{aligned}\end{align}

Additionally, the LBO-EM and LBO-ET use $u_{{s}r,i}=u_{{r}s,i}$ and $m_sv_{t,{s}r }^2=m_{r}v _{t,rs}^2$, respectively, as well as their corresponding collision frequencies (2.19) and (2.25). When the discrete expansions are inserted in (3.17) one is faced with a linear problem of size $({d_v}+1)N_b^x$ that must be solved in every position-space cell ($N_b^x$ is the number of monomials of the basis spanning position space). The LBO-G on the other hand uses the equality between discrete LBO relaxation rates and the FPO relaxation rates

(3.20)\begin{align} & m_s\nu_{sr}\left({u_{sr,i}}M_{0,s} - v_{t,{s}r}^2\int f_s\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i\right) - m_r\nu_{rs}\left(u_{{r}s,i}M_{0,r} - v_{t,rs}^2\int f_r\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i\right)\nonumber\\ & \quad\doteq m_s\nu_{sr}M_{1i,s} - m_r{\nu_{rs}}M_{1i,r} + \frac{\alpha_{E}\left(m_s+m_r\right)}{{M_{0,s}}M_{0,r}}\left({M_{0,s}}M_{1i,r}-{M_{0,r}}M_{1i,s}\right), \end{align}

(3.21)\begin{align} & m_s\nu_{sr}\left\{u_{{s}r,i}\left(M_{1i,s}-u_{s,i}M_{0,s}\right)+v_{t,{s}r}^2\left[{d_v}M_{0,s} - \int\left(v_i-u_{s,i}\right)f_s\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i\right]\right\}\nonumber\\ & \qquad- m_r\nu_{rs}\left\{u_{{r}s,i}\left(M_{1i,r}-u_{r,i}M_{0,r}\right)+v_{t,rs}^2\left[{d_v}M_{0,r} - \int\left(v_i-u_{r,i}\right)f_r\Big|_{v_{i,\min}}^{v_{i,\max}}\,\mathrm{d}S_i\right]\right\} \nonumber\\ & \quad\doteq m_s\nu_{sr}\left(M_{2,s}-u_{s,i}M_{1i,s}\right) - m_r\nu_{rs}\left(M_{2,r}-u_{r,i}M_{1i,r}\right) \nonumber\\ & \qquad+ \frac{\alpha_{E}}{{M_{0,s}}M_{0,r}}\left[\vphantom{\frac{m_r-m_s}{2}} m_{r}M_{0,s}\left(M_{2,r}-u_{r,i}M_{1i,r}\right)-m_sM_{0,r}\left(M_{2,s}-u_{s,i}M_{1i,s}\right) \right.\nonumber\\ & \qquad \left.+\frac{m_r-m_s}{2}\left(u_{s,i}-u_{r,i}\right)\left({M_{0,r}}M_{1i,s}-{M_{0,s}}M_{1i,r}\right)\right], \end{align}

where the relationship between $\alpha _{E}$ and $\nu _{sr}$ is given by (2.7). Therefore for the LBO-G (3.17)–(3.19) signify a $2({d_v}+1)N_b^x$ linear problem that must be solved in every position-space cell.

3.4.1. Discrete equations for the gyroaveraged operator

In § 2.4 we introduced the gyroaveraged cross-species LBO. Its discretization follows that outlined in Francisquez et al. (Reference Francisquez, Bernard, Mandell, Hammett and Hakim2020) and the calculation of the cross-primitive moments $u_{\parallel {s}r }$ and $v_{t,{s}r }$ is similar to that done for the non-gyroaveraged operator. The main differences arise from the fact that moments are defined via $v_\parallel$-$\mu$ integrals (e.g. (2.37)) and that the momentum density $M_{1\parallel,s}$ is a scalar instead of a vectorial quantity. The equations that arise from momentum and energy conservation in the gyroaveraged case are thus

(3.22)\begin{align} & m_s\nu_{sr}\left(M_{1\parallel,s} - {u_{\parallel{s}}}rM_{0,s} + v_{t,{s}r}^2\frac{2{\rm \pi}}{m_s}\int\mathcal{J}\hat{f_s}\Big|_{v_{{\parallel}\min}}^{v_{{\parallel}\max}}\,\mathrm{d}\mu\right) \nonumber\\ & \quad +m_r\nu_{rs}\left(M_{1\parallel,r} - {u_{\parallel{r}{s}}}M_{0,r} + v_{t,rs}^2\frac{2{\rm \pi}}{m_r}\int\mathcal{J}\hat{f_r}\Big|_{v_{{\parallel}\min}}^{v_{{\parallel}\max}}\,\mathrm{d}\mu\right) \doteq 0, \end{align}

(3.23)\begin{align} & m_s\nu_{sr}\left[{-}v_{t,{s}r}^2\left(3M_{0,s}-\frac{2{\rm \pi}}{m_s}\int v_\parallel\mathcal{J}\hat{f_s}\Big|_{v_{{\parallel}\min}}^{v_{{\parallel}\max}}\,\mathrm{d}\mu-\frac{2{\rm \pi}}{m_s}\int 2\mu\mathcal{J}\hat{f_s}\Big|_{\mu_{\min}}^{\mu_{\max}}\,\mathrm{d}v_\parallel\right)\right]\nonumber\\ & \quad +m_r\nu_{rs}\left[\vphantom{\frac{2{\rm \pi}}{m_r}}M_{2,r}-{u_{\parallel{r}{s}}}M_{1\parallel,r} \right.\nonumber\\ & \left.\quad-v_{t,rs}^2\left(3M_{0,r}-\frac{2{\rm \pi}}{m_r}\int v_\parallel\mathcal{J}\hat{f_r}\Big|_{v_{{\parallel}\min}}^{v_{{\parallel}\max}}\, \mathrm{d}\mu-\frac{2{\rm \pi}}{m_r}\int 2\mu\mathcal{J}\hat{f_r}\Big|_{\mu_{\min}}^{\mu_{\max}}\,\mathrm{d}v_\parallel\right)\right]\doteq 0. \end{align}

For the gyroaveraged LBO-EM and LBO-ET (3.20)–(3.21), and the relations $u_{\parallel {s}r }=u_{\parallel{r}{s} }$ and $m_sv_{t,{s}r }^2=m_{r}v _{t,rs}^2$, is all that is needed to compute the cross-primitive moments. This requires a solution to a linear problem of size $2N_b^x$ in each position-space cell. In the case of the LBO-G operator we make use of the discrete moment relaxation equations once again, as in (3.18)–(3.19) but using the gyroaveraged moments instead

(3.24)\begin{align} & m_s\nu_{sr}\left({u_{\parallel{s}}}rM_{0,s} - v_{t,{s}r}^2\frac{2{\rm \pi}}{m_s}\int\mathcal{J}\hat{f_s}\Big|_{v_{{\parallel}\min}}^{v_{{\parallel}\max}}\,\mathrm{d}\mu\right) \nonumber\\ & \qquad - m_r\nu_{rs}\left({u_{\parallel{r}{s}}}M_{0,r} - v_{t,rs}^2\frac{2{\rm \pi}}{m_r}\int\mathcal{J}\hat{f_r}\Big|_{v_{{\parallel}\min}}^{v_{{\parallel}\max}}\,\mathrm{d}\mu\right) \nonumber\\ & \quad\doteq m_s\nu_{sr}M_{1\parallel,s} - m_r{\nu_{rs}}M_{1\parallel,r} + \frac{\alpha_{E}\left(m_s+m_r\right)}{{M_{0,s}}M_{0,r}}\left({M_{0,s}}M_{1\parallel,r}-M_{0,r}M_{1\parallel,s}\right), \end{align}

(3.25)\begin{align} & m_s\nu_{sr}\left\{\vphantom{\frac{2{\rm \pi}}{m_s}}u_{{\parallel}{s}r}\left(M_{1\parallel,s}-{u_{\parallel{s}}}M_{0,s}\right) \right.\nonumber\\ & \qquad\left. +v_{t,{s}r}^2\left[3M_{0,s} - \frac{2{\rm \pi}}{m_s}\int\left(v_\parallel{-}u_{\parallel{s}}\right)\mathcal{J}\hat{f_s}\Big|_{v_{{\parallel}\min}}^{v_{{\parallel}\max}}\,\mathrm{d}\mu - \frac{2{\rm \pi}}{m_s}\int 2\mu\mathcal{J}\hat{f_s}\Big|_{\mu_{\min}}^{\mu_{\max}}\,\mathrm{d}v_\parallel\right]\right\}\nonumber\\ & \qquad -m_r\nu_{rs}\left\{\vphantom{\frac{2{\rm \pi}}{m_s}}u_{\parallel{r}{s}}\left(M_{1\parallel,r}-u_{\parallel{r}}M_{0,r}\right) \right.\nonumber\\ & \qquad \left. +v_{t,rs}^2\left[3M_{0,r} - \frac{2{\rm \pi}}{m_r}\int \left(v_\parallel{-}u_{\parallel{r}}\right)\mathcal{J}\hat{f_r}\Big|_{v_{{\parallel}\min}}^{v_{{\parallel}\max}}\,\mathrm{d}\mu - \frac{2{\rm \pi}}{m_r}\int 2\mu\mathcal{J}\hat{f_r}\Big|_{\mu_{\min}}^{\mu_{\max}}\,\mathrm{d}v_\parallel\right]\right\} \nonumber\\ & \quad\doteq m_s\nu_{sr}\left(M_{2,s}-u_{\parallel{s}}M_{1\parallel,s}\right) - m_r\nu_{rs}\left(M_{2,r}-{u_{\parallel{r}}}M_{1\parallel,r}\right) \nonumber\\ & \qquad+ \frac{\alpha_{E}}{{M_{0,s}}M_{0,r}}\left[\vphantom{\frac{m_r-m_s}{2}}m_{r}M_{0,s} \left(M_{2,r}-{u_{\parallel{r}}}M_{1\parallel,r}\right)-m_sM_{0,r}\left(M_{2,s}-u_{\parallel{s}}M_{1\parallel,s}\right) \right.\nonumber\\ & \qquad \left. +\frac{m_r-m_s}{2}\left(u_{\parallel{s}}-u_{\parallel{r}}\right)\left(M_{0,r}M_{1\parallel,s}-{M_{0,s}}M_{1\parallel,r}\right)\right], \end{align}

where once again weak division by ${M_{0,s}}M _{0,r}$ on the right-hand side of these equations assumed (Hakim et al. Reference Hakim, Francisquez, Juno and Hammett2020) in order to avoid aliasing errors. In (3.20)–(3.23) we assumed a $v_\parallel$-$\mu$ simulation such that ${d_v}=3$. For a $v_\parallel$ simulation the 3 in front of $M_{0,s}$ and $M_{0,r}$ would be simply a 1, $(2{\rm \pi} /m)\int \mathrm {d}\mu$ integrals would vanish and so would the the $2\mu \mathcal {J}\hat {f}$ terms. In either case, the gyroaveraged LBO-G requires inverting a matrix with $4N_b^x\times 4N_b^x$ matrix in each configuration space cell.

4.1. Conservation tests

4.1.1. Vlasov LBO conservation

We check that momentum and energy are indeed conserved by our discrete scheme by initiating two populations of electrons and protons with an arbitrary non-Maxwellian distribution function given by

(4.1)\begin{equation} f_s({ \boldsymbol{v} },t=0) = a\left[1+d\cos\left(k_s v\right)\right]\exp\left[-\frac{\left({ \boldsymbol{v} }-{ \boldsymbol{b} }_s\right)^2}{2\sigma_s^2}\right], \end{equation}

with $a=7\times 10^{19}$, $d=0.5$, $k_s={\rm \pi} /v_{t,s}$, $\sigma _s=v_{t,s}$, ${ \boldsymbol {b} }_e=\{v_{t,e}/2,-v_{t,e}/2,0\}$, ${ \boldsymbol {b} }_i=\{3v_{t,i}/2,3v_{t,i}/2,v_{t,i}/2\}$ with $v_{t,s}=\sqrt {T_{s0}/m_s}$, and $T_{e0}=40$ eV and $T_{i0}=80$ eV. We discretize these distribution functions in phase spaces restricted to $[-1,1]\times [-5v_{t,s},5v_{t,s}]^{{d_v}}$ and meshed with $1\times N _v^{{d_v}}$ cells. We show conservation properties for both piecewise-linear ($p=1$) and piecewise-quadratic ($p=2$) serendipity basis functions (Juno et al. Reference Juno, Hakim, TenBarge, Shi and Dorland2018); higher-order basis functions may also be used but the results do not change. For the same reason we use a single cell in position space; the results in this section are independent of position-space dimensionality although we checked such cases anyway to make sure there are no errors in the implementation.

We time integrate the cross-species collision terms (no self-species collisional or collisionless terms are included here) with constant collision frequency using a strong-stability preserving (SSP) third-order Runge–Kutta method (RK3). As electrons and ions collide with each other their temperatures and flow velocities relax to a common value, a process that is more carefully benchmarked in § 4.3. We also see that whatever anisotropies were present at $t=0$ go away on the $\nu _{sr}^{-1}$ time scale. In figure 1, for example, we illustrate the isotropization of the electrons after several $(\nu _{ei}^M)^{-1}$ periods as they collide with the ions using the LBO-EM, but since $\nu _{ie}$ is smaller by $m_e/m_i$ the ions will take much longer to isotropize as they collide with the electrons.

Figure 1.

Initial and final (after $\nu _{ei}^M t=100$) electron and ion distribution functions as they collide with each other with the LBO-EM model using $1\times 64^2$ cells and a $p=2$ serendipity basis. Initial conditions are given in (4.1). Colour bars are normalized to the extrema at $\nu _{ei}^M t=100$.

We ran this simulation for ${d_v}=\{1,2,3\}$, $p=\{1,2\}$ and using both the LBO-EM and the LBO-G. The ability to conserve the first three volume-integrated velocity moments of the distribution function was quantified in each case by integrating the equations for $N_t$ time steps, and computing the relative error per time step in the volume-integrated particle, momentum and kinetic energy density. The relative error per time step in the number density $M_0$ is given by

(4.2)\begin{equation} E_{r,M_0} = \frac{1}{N_t}\frac{\left\langle M_{0e}+M_{0i} \right\rangle(t=N_t\Delta t)-\left\langle M_{0e}+M_{0i} \right\rangle(t=0)}{\left\langle M_{0e}+M_{0i} \right\rangle(t=0)}, \end{equation}

where $\left \langle \cdot \right \rangle$ indicates a volume average and $N_t=10^4$. The relative error per time step in momentum and kinetic energy conservation takes into account the mass of each species

(4.3)\begin{gather} E_{r,M_{1k}} = \frac{1}{N_t}\frac{\left\langle m_eM_{1e,k}+m_iM_{1i,k} \right\rangle(t=N_t\Delta t)- \left\langle m_eM_{1e,k}+m_iM_{1i,k} \right\rangle(t=0)}{\left\langle m_eM_{1e,k}+m_iM_{1i,k} \right\rangle(t=0)}, \end{gather}

(4.4)\begin{gather}E_{r,M_2} = \frac{1}{N_t}\frac{\left\langle \frac{1}{2}m_eM_{2e}+\frac{1}{2}m_iM_{2i} \right\rangle(t=N_t\Delta t)- \left\langle \frac{1}{2}m_eM_{2e}+\frac{1}{2}m_iM_{2i} \right\rangle(t=0)}{\left\langle \frac{1}{2}m_eM_{2e}+\frac{1}{2}m_iM_{2i} \right\rangle(t=0)}. \end{gather}

The results as a function of velocity-space resolution (i.e. $N_v$) are given in figure 2. The middle row plots have ${d_v}$ lines for each operator and polynomial order $p$ because the relative error per time step in the volume-integrated momentum density is measured along each direction separately. In all cases we see that the errors in momentum conservation per time step remain of the order of machine precision. This is true even for the simulations with piecewise-linear basis functions or very coarse velocity-space meshes. The LBO-ET uses the same algorithm and implementation as LBO-EM but with a different collision frequency, so its conservation errors are similar to those of the LBO-EM shown here.

Figure 2.

Relative error per time step in number density (a–c), momentum density in each direction (d–f) and kinetic energy density (g–i). The middle row plots contain ${d_v}$ lines for each operator and $p$, corresponding to the error per time step in the conservation of momentum along each direction.

These conservation properties do not depend on the large mass disparity between ions and electrons; the algorithm's ability to conserve the velocity moments is also independent of the mass ratio. We provide as an example the ${d_v}=2$ and $p=2$ simulation with the LBO-EM, scanning the number of velocity-space cells in one direction ($N_v$) and using the mass ratios $m_i/m_e=\{300,600,1000,1836\}$. The conservation errors for these simulations are provided in figure 3, once again shown that for all mass ratios and resolutions used, the error per time step in the volume-integrated velocity moments remains of the order of machine precision.

Figure 3.

Relative error per time step in number density (a), momentum density in each direction (b) and kinetic energy density (c) in two-dimensional velocity space (${d_v}=2$).

4.1.2. Gyroaveraged LBO conservation

Similar tests were run with the gyroaveraged version of the LBO operators in order to guarantee that the algorithm remains conservative in that case as well. For these tests we initialize the ion and electron distribution functions with

(4.5)\begin{equation} f_s(v_\parallel,\mu,t=0) = \frac{a\left[1+d\cos\left(k_s v\right)\right]}{\left(2{\rm \pi}\sigma_s^2\right)^{3/2}}\exp\left[-\frac{\left(v_\parallel{-}b_s\right)^2+2\mu Bm_s}{2\sigma_s^2}\right], \end{equation}

and the parameters $B=1.2$ T, $a=7\times 10^{19}$, $d=0.5$, $k_s={\rm \pi} /v_{t,s}$, $\sigma _s=v_{t,s}$, $b_e=5v_{t,i}/4$, $b_i=v_{t,i}$, $v_{t,s}=\sqrt {T_{s 0}/m_s}$, $T_{e0}=40$ eV and $T_{i0}=80$ eV. The phase space $[-1,1]\times [-5v_{t,s},5v_{t,s}]\times [0,m_s(5v_{t,s})^2/(2B)]$ is meshed with $1\times N _v^2$ cells and functions are expanded on piecewise-linear ($p=1$) or piecewise-quadratic ($p=2$) serendipity basis.

We allow the electrons and ions to collide with each other but not with themselves, and we do not apply the collisionless terms either. The cross-species collision terms were integrated in time for $10^4$ time steps using a third-order SSP RK3, and we computed the relative error per time step in the volume-integrated velocity moments as in § 4.1.2. The results in figure 4 demonstrate how the relative error per time step in the conservation of velocity moments stays of order of machine precision for all velocity-space resolutions, and even for $p=1$. Figure 4 gives conservation errors for the LBO-ET and the LBO-G; the LBO-EM has similar conservative properties as the LBO-ET since it only differs by the definition of the collision frequency.

Figure 4.

Relative error per time step in number density (a,b), momentum density (c,d) and kinetic energy density (e, f) as a function of the number of cells along one direction of velocity space ($N_v$) for the gyroaveraged LBO-ET (solid) and LBO-G (dashed). Top row contains tests with $v_\parallel$-space only, while the bottom row contains tests in $v_\parallel -\mu$ space.

4.2. Landau damping of electron Langmuir waves

A seminal test bed for collision operators is the Landau damping of plasma waves across the collisional range. We pursued this analysis to examine the effect that these collision models have on the Landau damping rate of electrostatic electron Langmuir waves. For this purpose we employ the Vlasov–Maxwell solver in Gkeyll (Juno et al. Reference Juno, Hakim, TenBarge, Shi and Dorland2018; Hakim & Juno Reference Hakim and Juno2020) with the self-species collision terms (Hakim et al. Reference Hakim, Francisquez, Juno and Hammett2020) and the multispecies collision models described in this work. The hydrogen ions are fixed in time so the equations solved are

(4.6)\begin{align} \frac{\partial f_e}{\partial t} + { \boldsymbol{v} }\boldsymbol{\cdot}\nabla f_e - { \boldsymbol{E} }\boldsymbol{\cdot}\nabla_{\boldsymbol{v}}{f_e} &= \sum_{r=e,i}\nu_{er} \nabla_{\boldsymbol{v}}\cdot{\left[\left({ \boldsymbol{v} }-{ \boldsymbol{u} }_{er}\right)f_e+v_{ter}^2\nabla_{\boldsymbol{v}}{f_e}\right]}, \end{align}

(4.7)\begin{align}\frac{\partial { \boldsymbol{E} }}{\partial t} &={-}{ \boldsymbol{J} }, \end{align}

where we used normalized units,Footnote ⁴ the current density is given by ${ \boldsymbol {J} }=-{ \boldsymbol {M_{1e}} }$ and we solve (4.7) in a way that keeps the simulation electrostatic.Footnote ⁵ We use four-dimensional simulations with the phase space $[-{\rm \pi} /k,{\rm \pi} /k]\times [-5v_{t,e},5v_{t,e}]^3$ discretized by $16\times 36^3$ cells and $p=2$ basis functions, for the wavenumber $k\lambda _{De}=0.3$, with $\lambda _{De}$ being the electron Debye length. We confirmed that the resolution used is the minimum needed to obtain converged results by scanning the position- and velocity-space resolution as well as the velocity-space extents. The static ions have the normalized density $n_i(x)=1$ while the electrons are initialized with a non-drifting Maxwellian distribution that has the temperature $T_e=1=T_i$ and the density $n_e(x,t=0)=1+\alpha \cos (kx)$, with $\alpha =10^{-4}$. The electric field is initialized in a manner consistent with Poisson's equation ${ \boldsymbol {E} }=- {\boldsymbol {\hat {x}}}\, \alpha \sin (kx)/k$.

As the simulation proceeds we see the amplitude of the electrostatic wave damp, which can be appreciated by examining the volume-integrated field energy over time as shown in figure 5(a). We can quantify the rate at which these waves damp and plot it as a function of collision frequency as is done in figure 5(b). If one were to only use self-species collisions one would obtain the results shown with black crosses, and for that case the equations are sufficiently simple that one can obtain an analytic dispersion relation (Anderson & O'Neil Reference Anderson and O'Neil2007; Francisquez et al. Reference Francisquez, Bernard, Mandell, Hammett and Hakim2020) which agrees well with the numerical results (black circles), providing additional confidence in the Gkeyll implementation. When we introduce electron–ion collisions obtaining analytic growth rates is more difficult. So we instead compare the results obtained with the LBO-G (solid blue), the LBO-EM (dashed purple) and the LBO-ET (dotted green) with previously reported results for the FPO (Jorge et al. Reference Jorge, Ricci, Brunner, Gamba, Konovets, Loureiro, Perrone and Teixeira2019).

Figure 5.

(a) Normalized volume-integrated field energy in simulations of electron Langmuir waves as they Landau damp over time. (b) Damping rates of electron Langmuir waves given as the offset from the collisionless value of $\gamma _{\mathrm {collisionless}}/\omega _{pe}=-1.247\times 10^{-2}$ for the LBO-G with $\beta =0$ and $\beta =-0.65$ (blue circles), the LBO-EM (purple triangles) and the LBO-ET (dotted green) compared with those with a full FPO (Jorge et al. Reference Jorge, Ricci, Brunner, Gamba, Konovets, Loureiro, Perrone and Teixeira2019) (solid red). If using only self-collisions, numerical Gkeyll data (black crosses) agrees with theory (black circles).

The LBO-G simulations were performed using $\nu _{ei}^G=\sqrt {2}\nu _{ee}$ ($=(m_i/m_e)\nu _{ie}$) since this is the relationship assumed in the reference FPO work (Jorge et al. Reference Jorge, Ricci, Brunner, Gamba, Konovets, Loureiro, Perrone and Teixeira2019). Figure 5(b) suggests that the LBO-G can provide a more accurate description of this kinetic phenomenon than, say, using self-species collisions only. There is the caveat, however, that we have not established from first principles what the most suitable choice of the free parameter $\beta$ ought to be at any given collision frequency. We scanned this parameter and show the results for $\beta =0$ and $\beta =-0.65$, the latter bringing the damping rates closer to those of the full FPO. But it is apparent that the wrong choice of $\beta$ can also result in significant deviation from the FPO.

Also shown in figure 5(b) are the damping rates obtained when using the LBO-EM and the LBO-ET. Despite having a different model for $u_{{s}r }$ and $v_{t,{s}r }$, LBO-EM has the same collision frequency ($\nu _{ei}^M=\nu _{ei}^G(\beta =0)$) and gives the same damping rates as the LBO-G with $\beta =0$ (top blue and dashed purple lines). The LBO-ET on the other hand has a collision frequency that is smaller by the mass ratio ($\nu _{ei}^T = m_e\nu _{ei}^M/(m_e+m_i)\simeq (m_e/m_i)\nu _{ei}^M$), and therefore is essentially equivalent to neglecting cross-species collisions for this problem; i.e. solid black and dotted green lines agree. If we were to run the simulation with the LBO-ET but the same value of $\nu _{ei}$ as the LBO-EM then we would simply obtain the same results as if we had used the LBO-EM.

4.3. Velocity and temperature relaxation

As a final benchmark of the multispecies LBO algorithms and solvers we employ the gyroaveraged LBO to model the relaxation of a deuterium plasma to thermal equilibrium. We employ identical conditions, as best as we can tell, to those used in a benchmark of the FPO in the XGC code (Hager et al. Reference Hager, Yoon, Ku, D'Azevedo, Worley and Chang2016). The same test was recently performed with a finite-volume implementation of the LBO-ET and LBO-EM in GENE-X Ulbl et al. (Reference Ulbl, Michels and Jenko2021). This means that the initial distribution functions are described by the bi-Maxwellians

(4.8)\begin{equation} f_s\left(v_\parallel,\mu,t=0\right) = \frac{n_0}{\alpha\left(2{\rm \pi} v_{t,s}^2\right)^{3/2}}\exp\left[-\frac{\left(v_\parallel{-}u_{\parallel{s}}\right)^2+2\mu B/(m_s\alpha)}{2v_{t,s}^2}\right], \end{equation}

where $B=1$ T, $v_{t,s}=\sqrt {T_{s 0}/m_s}$ $\alpha =1.3$, $u_{\parallel i}=50(m_e/m_i)v_{t,i}$, $u_{\parallel e}=0.5\sqrt {m_e/m_i}v_{t,s}$, $T_{i0}=200$ eV and $T_{e0}=300$ eV. Note that these reference temperatures are slightly different than the true initial temperatures $T_i(t=0)=240$ eV and $T_e(t=0)=360$ eV given by $T_s=(2T_{\perp s} + T_{\parallel s})/3$. For this test we once again neglect the collisionless terms and use a collision frequency that depends on time, i.e. $\nu _{sr}=\nu _{sr}(T_s(t),T_r(t))$. The phase space $[-2,2]\times [-5v_{t,s},5v_{t,s}]\times [0,m_s(5v_{t,s})^2/(2B)]$ is meshed with $1\times 16^2$ cells and dynamic fields are expanded in a piecewise-linear ($p=1$) basis. This resolution and velocity-space extents were confirmed as sufficient by convergence tests.

Figure 6 provides the time evolution of the parallel and perpendicular temperatures for the LBO-G ($\beta =0$), LBO-EM and LBO-ET operators compared with the previously reported FPO resultsFootnote ⁶ (Hager et al. Reference Hager, Yoon, Ku, D'Azevedo, Worley and Chang2016). The first event is the isotropization of the electrons followed by the isotropization of the ions, happening on the $\nu _{ss}^{-1}$ time scale. We used

(4.9)\begin{equation} \nu_{ss} = \frac{1}{\sqrt{2}}\frac{q_{s}^4n_s\log\varLambda_{ss}}{3\left(2{\rm \pi}\right)^{3/2}\epsilon_0^2m_s^2v_{t,s}^3}, \end{equation}

for the like-species scattering rate. Note that the LBO-G and LBO-EM exhibit a delayed isotropization time compared with the FPO's, an observation that Pezzi, Valentini & Veltri (Reference Pezzi, Valentini and Veltri2015) had also made while comparing the self-species Dougherty operator with the FPO (in Landau form). Later, on the $\nu _{ie}^{-1}$ time scale, we see the electrons and ions come into thermal equilibrium with each other, a process that is better described by both the LBO-G and the LBO-ET operators since, after all, the LBO-EM made no attempt at matching the FPO thermal relaxation rates. The time axis on these plots has been normalized to the isotropization rate (Huba Reference Huba2013)

(4.10)\begin{equation} \nu_T^i = \frac{2\sqrt{\rm \pi}q_i^4n_i\log\varLambda_{ii}}{\left(4{\rm \pi}\epsilon_0\right)^2m_i^2v_{t,i}^3}A^{{-}2}\left[{-}3+\left(A+3\right)a_F\right], \end{equation}

where $v_{t,i}$ and $\log \varLambda _{ii}$ use the initial ion temperature (240 eV), and $A=T_{\perp i}/T_{\parallel i}-1$ and $a_F=A^{-1/2}\tan ^{-1}(A^{1/2})$ if $A>0$ or $a_F=(-A)^{1/2}\tanh ^{-1}(-A)^{1/2}$ if $A<0$.

Figure 6.

Isotropization and relaxation of temperatures as deuterium ions and electrons collide with themselves and each other, compared with the results from the full FPO (Hager et al. Reference Hager, Yoon, Ku, D'Azevedo, Worley and Chang2016). (a) LBO-G with $\beta =0$. (b) LBO-EM. (c) LBO-ET.

We can also examine the velocity evolution as the plasma approaches an equilibrium as is done in figure 7. The first thing we notice is that the LBO-ET (green dash-dot with circles) grossly overestimates the time scale on which the electron flow relaxes to the ion flow, which happens because the ions are so much more massive. By definition the LBO-ET did not attempt to match the momentum relaxation rate, and we see the result of that here. On the other hand, the LBO-G (solid blue) and LBO-EM (dashed orange with crosses) models do a better job of approximating the FPO results for the slowing down of electrons, since their formulation included matching the FPO's momentum relaxation rates. There is still a discrepancy, e.g. between solid blue and dashed orange lines, although we point out that the LBO-G and LBO-EM would appear to match the analytic result based on the flow relaxation frequency given by the friction force at large mass ratio (Hinton & Hazeltine Reference Hinton and Hazeltine1976) (see figure 4 of Hager et al. Reference Hager, Yoon, Ku, D'Azevedo, Worley and Chang2016).

Figure 7.

Relaxation of parallel flow speed as deuterium ions and electrons collide with themselves and each other. Electron flow speeds are compared with the results from the full FPO (dotted black, Hager et al. (Reference Hager, Yoon, Ku, D'Azevedo, Worley and Chang2016)). Flow speeds when using the LBO-G ($\beta =0$) are given in solid dark blue, results for the LBO-EM in dashed orange with crosses and green dash-dot lines with circles represent results obtained with the LBO-ET.