# A Nonlinear, Conservative, Entropic Fokker-Planck Model for Multi-Species
Collisions ††thanks:  This manuscript has been authored, in part, by UT-
Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of
Energy (DOE). The U.S. government retains and the publisher, by accepting the
article for publication, acknowledges that the U.S. government retains a
nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce
the published form of this manuscript, or allow others to do so, for U.S.
government purposes. DOE will provide public access to these results of
federally sponsored research in accordance with the DOE Public Access Plan
(https://energy.gov/downloads/doe-public-access-plan).
This work was supported, in part, by the U.S. Department of Energy through the
Los Alamos National Laboratory. Los Alamos National Laboratory is operated by
Triad National Security, LLC, for the National Nuclear Security Administration
of U.S. Department of Energy (Contract No. 89233218CNA000001).
J. Hu’s research is partially supported by AFOSR grant FA9550-21-1-0358, and
DOE grant DE-SC0023164.

Evan Habbershaw, Cory Hauck, Jingwei Hu, and Jeffrey R. Haack  Department of
Mathematics, The University of Tennessee, Knoxville, TN, 37996 Department of
Mathematics, The University of Tennessee, Knoxville, TN, 37996 and Computer
Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge,
TN, 37831 Department of Applied Mathematics, University of Washington,
Seattle, WA, 98195 Computer, Computational, and Statistical Sciences Division,
Los Alamos National Laboratory, Los Alamos, NM 87545

(17 April 2024)

###### Abstract

A multi-species Fokker-Planck model for simulating particle collisions in a
plasma is presented. The model includes various parameters that must be tuned.
Under reasonable assumptions on these parameters, the model satisfies
appropriate conservation laws, dissipates an entropy, and satisfies an
$\mathcal{H}$-Theorem. In addition, the model parameters provide the
additional flexibility that is used to match simultaneously momentum and
temperature relaxation formulas derived from the Boltzmann collision operator
for a binary mixture with Coulomb potentials. A numerical method for solving
the resulting space-homogeneous kinetic equation is presented and two examples
are provided to demonstrate the relaxation of species bulk velocities and
temperatures to their equilibrium values.

Keywords: Multi-species kinetic models, Fokker-Planck, Lenard-Bernstein,
Dougherty, moment equations, entropy, conservation laws, plasma physics

##  1 Introduction

In this paper, we propose and analyze a multi-species version of a kinetic
collision model for plasma physics applications, often referred to as the
Lenard-Bernstein [15] or Dougherty model [5]. The single-species version of
this model (which we refer to as LB) is a nonlinear Fokker-Planck operator
that is used as a computationally cheaper surrogate for the full Landau-
Fokker-Planck (LFP) operator [14]. Roughly speaking, the LB model replaces the
nonlocal advection and diffusion terms of the LFP operator by parameters that
depend only on the bulk velocity and temperature of the underlying kinetic
distribution [16], while maintaining the conservation, entropy dissipation,
and equilibrium properties of the original LFP operator. In this way, the
approximation of LFP by LB is similar in spirit to the approximation of the
Boltzmann collision operator (see e.g, [4]) by the Bhatnagar-Gross-Krook (BGK)
model [1].

Recently, there have been efforts to extend the LB operator to the multi-
species case [20, 7] (see [12] for the well-posedness result of a similar
model) in a manner that maintains relevant (global) conservation properties,
entropy dissipation, and an $\mathcal{H}$-Theorem that uniquely characterizes
equilibrium states. These multi-species models are written as a sum of Fokker-
Planck operators, with coefficients that can be interpreted as tunable
collision frequencies. Each Fokker-Planck term is parameterized by species
bulk velocity and temperature in the case of intra-species collisions or by
mixture values for the bulk velocity and temperature in the case of inter-
species collisions. The choice of mixture values uses the strategy for the
multi-species BGK (M-BGK) model introduced in [9], which preserves the
aforementioned structural properties. However, like [9], the models in [20,
7], can only match the pairwise relaxation rate (derived from the Boltzmann
collision operator for a binary mixture with Coulomb potential) 111This
matching procedure is relevant since for small angle collisions, the Boltzmann
collision operator for a binary mixture with Coulomb potential is approximated
by the LFP operator [8]. of either the bulk velocity or temperature, but not
both at the same time. This limitation is a consequence of the fact that
mixture values for the bulk velocity and temperature are assumed to be
symmetric with respect to a permutation in indices between two species.

A more general M-BGK model was introduced in [13] (see also [2]) that does not
make the symmetry assumption in [9]. Rather, the model in [13], which includes
the one in [9] as a special case, introduces additional parameters in the
mixture values to enable matching of pairwise bulk velocity and temperature
relaxation rates at the same time. In this work, we take a similar approach
and construct a more general multi-species LB (M-LB) model that satisfies the
requisite conservation, entropy dissipation, and equilibrium
($\mathcal{H}$-Theorem) properties, but contains additional parameters for
matching of velocity and temperature relaxation rates at the same time. The
parameters however, have constraints; thus a principled approach for matching
must be carried out. We develop a procedure for doing so, matching the rates
of the M-LB model with those derived from the Boltzmann collision operator
with Coulomb potentials. We then illustrate the application of the matching
procedure with several numerical experiments.

The remainder of the paper is organized as follows. In Section 2, we introduce
the M-LB model. In Section 3, we review conservation laws and derive required
physically relevant conditions on the mixture parameters for the bulk velocity
and temperature. In Section 4, we prove an entropy dissipation law and an
$\mathcal{H}$-Theorem. In Section 5, we describe the procedure for matching
the pairwise relaxation rates to the multi-species Boltzmann equation with
Coulomb potentials. In Section 6, we present an implicit method for solving
the homogeneous M-LB equations, which includes an iterative procedure for
updating the moments when presented with moment-dependent collision rates. In
Section 7, we present some preliminary numerical results for the space
homogeneous problem. In the appendix we present proofs used in the main body
of the paper, and provide comparisons to the recent work in [18].

##  2 Multi-species space homogeneous Lenard-Bernstein equation

Let $f_{i}(t,\boldsymbol{v})$ denote the kinetic distribution of particles of
species $i\in\\{1,\cdots,N\\}$ having mass $m_{i}$. Associated to each $f_{i}$
are the species number density $n_{i}$, mass density $\rho_{i}$, bulk velocity
$\boldsymbol{u}_{i}$, and temperature $T_{i}$, defined by

$\rho_{i}=m_{i}n_{i}=m_{i}{\int_{\mathbb{R}^{d}}}f_{i}\,\,\textnormal{d}\boldsymbol{v},\;\;\;\boldsymbol{u}_{i}=\frac{1}{\rho_{i}}{\int_{\mathbb{R}^{d}}}m_{i}\boldsymbol{v}f_{i}\,\,\textnormal{d}\boldsymbol{v},\;\;\;T_{i}=\frac{m_{i}}{n_{i}d}{\int_{\mathbb{R}^{d}}}\left|\boldsymbol{v}-\boldsymbol{u}_{i}\right|^{2}f_{i}\,\,\textnormal{d}\boldsymbol{v}.$
(1)

The species energy densities, $E_{i}$, are given by

$E_{i}=\frac{1}{2}{\int_{\mathbb{R}^{d}}}m_{i}|\boldsymbol{v}|^{2}f_{i}\,\,\textnormal{d}\boldsymbol{v}=\frac{1}{2}\rho_{i}|\boldsymbol{u}_{i}|^{2}+\frac{d}{2}n_{i}T_{i}.$
(2)

######  Definition 2.1.

For $n>0$, $\boldsymbol{u}\in\mathbb{R}^{d}$, and $\theta>0$, define the
Maxwellian

$M_{n,\boldsymbol{u},\theta}=\frac{n}{(2\pi\theta)^{d/2}}\exp\left(-\frac{|\boldsymbol{v}-\boldsymbol{u}|^{2}}{2\theta}\right).$
(3)

######  Lemma 2.2.

For a general Maxwellian $M_{n,\boldsymbol{u},\theta}$,

$\displaystyle{\int_{\mathbb{R}^{d}}}M_{n,\boldsymbol{u},\theta}\,\,\textnormal{d}\boldsymbol{v}=n,\qquad{\int_{\mathbb{R}^{d}}}\boldsymbol{v}M_{n,\boldsymbol{u},\theta}\,\,\textnormal{d}\boldsymbol{v}=n\boldsymbol{u},$
(4)
$\displaystyle{\int_{\mathbb{R}^{d}}}|\boldsymbol{v}|^{2}M_{n,\boldsymbol{u},\theta}\,\,\textnormal{d}\boldsymbol{v}=n|\boldsymbol{u}|^{2}+dn\theta,\qquad\mbox{and}\qquad\Delta_{\boldsymbol{v}}\log(M_{n,\boldsymbol{u},\theta})=-\frac{d}{\theta}.$
(5)

######  Definition 2.3.

The (space homogeneous) multi-species Lenard Bernstein (M-LB) model is

$\frac{\partial f_{i}}{\partial
t}=\sum_{j=1}^{N}\mathcal{C}_{i,j}(f_{i},f_{j}),\qquad i=1,\cdots,N,$ (6)

where, for $1\leq i,j\leq N$, the pairwise collision operators
$\mathcal{C}_{i,j}$ are given by

$\displaystyle\mathcal{C}_{i,j}(f_{i},f_{j})$
$\displaystyle\coloneqq\lambda_{i,j}\theta_{i,j}\nabla_{\boldsymbol{v}}\cdot\left(\frac{\boldsymbol{v}-\boldsymbol{u}_{i,j}}{\theta_{i,j}}f_{i}+\nabla_{\boldsymbol{v}}f_{i}\right)=\lambda_{i,j}\theta_{i,j}\mathcal{L}(f_{i},M_{i,j}),$
(7a) $\displaystyle\mathcal{L}(f,M)$
$\displaystyle\coloneqq\nabla_{\boldsymbol{v}}\cdot\left(M\nabla_{\boldsymbol{v}}\left(\frac{f}{M}\right)\right)=\nabla_{\boldsymbol{v}}\cdot\left[\nabla_{\boldsymbol{v}}f-f\nabla_{\boldsymbol{v}}(\log(M))\right].$
(7b)

The parameters $\lambda_{i,j}\geq 0$ are the frequencies of collisions between
species $i$ and $j$. The mixture Maxwellians are given by

$M_{i,j}\coloneqq
M_{n_{i},\boldsymbol{u}_{i,j},\theta_{i,j}},\qquad\mbox{where}\qquad\theta_{i,j}\coloneqq
T_{i,j}/m_{i},$ (8)

and the mixture Maxwellian parameters are given by

$\displaystyle\boldsymbol{u}_{i,j}$
$\displaystyle\coloneqq\alpha_{i,j}\boldsymbol{u}_{i}+(1-\alpha_{i,j})\boldsymbol{u}_{j},$
(9a) $\displaystyle T_{i,j}$
$\displaystyle\coloneqq\beta_{i,j}T_{i}+(1-\beta_{i,j})T_{j}+\frac{1}{d}\gamma_{i,j}|\boldsymbol{u}_{i}-\boldsymbol{u}_{j}|^{2},$
(9b)

where the parameters $\alpha_{i,j}$, $\beta_{i,j}$, and $\gamma_{i,j}$ are to
be determined.

######  Remark 2.4.

When $i=j$, we often drop the second index. Thus, for example, $M_{i,i}=M_{i}$
is the unique Maxwellian associated to $f_{i}$ with bulk velocity
$\boldsymbol{u}_{i,i}=\boldsymbol{u}_{i}$ and temperature $T_{i,i}=T_{i}$.

##  3 Conservation laws and mixture parameter constraints

Throughout the paper, we make the following assumption, which allows us to use
the chain rule for time derivatives, apply integration by parts in velocity,
and, when doing the latter, set any boundary terms that arise to zero. While
more precise conditions can be formulated, we do not do so here.

######  Assumption 3.1.

For each $i\in\\{1,\cdots,N\\}$, $f_{i}>0$ is a smooth function of
$\boldsymbol{v}$ and continuously differentiable in $t$. For each $t$, $f_{i}$
decays exponentially to zero as $|\boldsymbol{v}|\to\infty$.

The M-LB model satisfies conservation of species mass:
$m_{i}{\int_{\mathbb{R}^{d}}}\mathcal{C}_{i,j}(f_{i},f_{j})\,\,\textnormal{d}\boldsymbol{v}=0$.
We assume further that it satisfies the pairwise conservation of momentum and
energy.

######  Assumption 3.2.

The collision operators in (6) satisfy the following invariance properties,
which correspond to pairwise momentum conservation and pairwise energy
conservation, respectively. Specifically, for any $i,j$,

$\displaystyle
m_{i}\langle\\!\;\boldsymbol{v}\mathcal{C}_{i,j}(f_{i},f_{j})\\!\;\rangle+m_{j}\langle\\!\;\boldsymbol{v}\mathcal{C}_{j,i}(f_{j},f_{i})\\!\;\rangle$
$\displaystyle=\mathbf{0},$ (10a) $\displaystyle
m_{i}\left\langle\frac{|\boldsymbol{v}|^{2}}{2}\mathcal{C}_{i,j}(f_{i},f_{j})\right\rangle+m_{j}\left\langle\frac{|\boldsymbol{v}|^{2}}{2}\mathcal{C}_{j,i}(f_{j},f_{i})\right\rangle$
$\displaystyle=0,$ (10b)

where brackets are used as a shorthand for velocity integration, i.e.,
$\langle\\!\;\boldsymbol{g}\\!\;\rangle={\int_{\mathbb{R}^{d}}}{\boldsymbol{g}}({\boldsymbol{v}})\,\textnormal{d}{\boldsymbol{v}}$.

Evaluation of these conservation laws relies on the following lemma.

######  Lemma 3.3.

For every $i,j\in\\{1,\cdots,N\\}$,

$\displaystyle\langle\\!\;m_{i}\mathcal{L}(f_{i},M_{i,j})\\!\;\rangle$
$\displaystyle=0,$ (11a)
$\displaystyle\langle\\!\;m_{i}\boldsymbol{v}\mathcal{L}(f_{i},M_{i,j})\\!\;\rangle$
$\displaystyle=-\frac{\rho_{i}(\boldsymbol{u}_{i}-\boldsymbol{u}_{i,j})}{\theta_{i,j}},$
(11b) $\displaystyle\left\langle
m_{i}\frac{|\boldsymbol{v}|^{2}}{2}\mathcal{L}(f_{i},M_{i,j})\right\rangle$
$\displaystyle=-\frac{\rho_{i}\boldsymbol{u}_{i}\cdot(\boldsymbol{u}_{i}-\boldsymbol{u}_{i,j})}{\theta_{i,j}}-\frac{d\rho_{i}(\theta_{i}-\theta_{i,j})}{\theta_{i,j}}.$
(11c)

###### Proof.

The proof is a straightforward, but tedious computation, using integration by
parts several times. ∎

######  Proposition 3.4.

Given the definitions of $\boldsymbol{u}_{i,j}$ and $T_{i,j}$ in (9), the
pairwise conservation of momentum and energy in Assumption 3.2 holds if and
only if the following parameter relationships hold:

$\displaystyle\alpha_{j,i}$
$\displaystyle=1-\frac{\rho_{i}\lambda_{i,j}}{\rho_{j}\lambda_{j,i}}(1-\alpha_{i,j}),$
(12a) $\displaystyle\beta_{j,i}$
$\displaystyle=1-\frac{n_{i}\lambda_{i,j}}{n_{j}\lambda_{j,i}}(1-\beta_{i,j}),$
(12b) $\displaystyle\gamma_{j,i}$
$\displaystyle=\frac{1}{n_{j}\lambda_{j,i}}(\rho_{i}\lambda_{i,j}(1-\alpha_{i,j})-n_{i}\lambda_{i,j}\gamma_{i,j}).$
(12c)

###### Proof.

We show that Assumption 3.2 implies (12). The converse follows by simply
reversing the order of computations. According to (11b), the conservation of
momentum (10a) implies that

$\rho_{i}\lambda_{i,j}\boldsymbol{u}_{i,j}+\rho_{j}\lambda_{j,i}\boldsymbol{u}_{j,i}=\rho_{i}\lambda_{i,j}\boldsymbol{u}_{i}+\rho_{j}\lambda_{j,i}\boldsymbol{u}_{j}.$
(13)

Rearranging this equation to isolate $\boldsymbol{u}_{j,i}$ and applying (9a)
yields (12b). By (11c), the conservation of energy (10b) implies that

$n_{i}\lambda_{i,j}T_{i,j}+n_{j}\lambda_{j,i}T_{j,i}=\frac{1}{d}\rho_{i}\lambda_{i,j}(1-\alpha_{i,j})|\boldsymbol{u}_{i}-\boldsymbol{u}_{j}|^{2}+n_{i}\lambda_{i,j}T_{i}+n_{j}\lambda_{j,i}T_{j}.$
(14)

Rearranging this equation to isolate $T_{j,i}$ and using (9b) yields (12b) and
(12c). ∎

######  Remark 3.5.

The relations in (12) can be rearranged as follows, in order to highlight the
symmetries in the model:

$\displaystyle\delta_{i,j}\coloneqq\rho_{i}\lambda_{i,j}(1-\alpha_{i,j})$
$\displaystyle=\rho_{j}\lambda_{j,i}(1-\alpha_{j,i})\eqqcolon\delta_{j,i},$
(15a) $\displaystyle n_{i}\lambda_{i,j}(1-\beta_{i,j})$
$\displaystyle=n_{j}\lambda_{j,i}(1-\beta_{j,i}),$ (15b)
$\displaystyle\delta_{i,j}-2n_{i}\lambda_{i,j}\gamma_{i,j}$
$\displaystyle=2n_{j}\lambda_{j,i}\gamma_{j,i}-\delta_{j,i}.$ (15c)

######  Remark 3.6.

If we constrain the mixture velocity terms to enforce the symmetry
$\boldsymbol{u}_{i,j}=\boldsymbol{u}_{j,i}$, then
$\alpha_{i,j}+\alpha_{j,i}=1$, and

$\boldsymbol{u}_{i,j}=\alpha_{i,j}\boldsymbol{u}_{i}+\alpha_{j,i}\boldsymbol{u}_{j},\qquad\mbox{where}\qquad\alpha_{i,j}=\frac{\rho_{i}\lambda_{i,j}}{\rho_{i}\lambda_{i,j}+\rho_{j}\lambda_{j,i}}.$
(16)

This is the same expression for the mixture velocity as in [9]. If, in
addition, we require $T_{i,j}=T_{j,i}$, then $\beta_{i,j}+\beta_{j,i}=1$,
$\gamma_{i,j}=\gamma_{j,i}$, and

$T_{i,j}=\beta_{i,j}T_{i}+\beta_{j,i}T_{j}+\frac{1}{d}\frac{\delta_{i,j}|\boldsymbol{u}_{i}-\boldsymbol{u}_{j}|^{2}}{\lambda_{i,j}n_{i}+\lambda_{j,i}n_{j}},$
(17)

where

$\beta_{i,j}=\frac{n_{i}\lambda_{i,j}}{n_{i}\lambda_{i,j}+n_{j}\lambda_{j,i}}\qquad\mbox{and}\qquad\delta_{i,j}=\frac{\rho_{i}\rho_{j}\lambda_{i,j}\lambda_{j,i}}{\rho_{i}\lambda_{i,j}+\rho_{j}\lambda_{j,i}}.$
(18)

In this case

$\displaystyle T_{i,j}$
$\displaystyle=\frac{n_{i}\lambda_{i,j}T_{i}+n_{j}\lambda_{j,i}T_{j}}{n_{i}\lambda_{i,j}+n_{j}\lambda_{j,i}}+\frac{1}{d}\frac{\rho_{i}\rho_{j}\lambda_{i,j}\lambda_{j,i}}{(\rho_{i}\lambda_{i,j}+\rho_{j}\lambda_{i,j})(n_{i}\lambda_{i,j}+n_{j}\lambda_{j,i})}|\boldsymbol{u}_{i}-\boldsymbol{u}_{j}|^{2},$
(19)

which is the mixture temperature in [9].

The following moment ODE system, describing the time evolution of the momenta,
and energies can be derived (see Appendix B) from the model equation (6), by
integrating against $m_{i}(\boldsymbol{v},|\boldsymbol{v}|/2)^{\top}$, and
using (11) and the definitions in (9):

$\displaystyle\frac{\partial(\rho_{i}\boldsymbol{u}_{i})}{\partial t}$
$\displaystyle=\sum_{j}\rho_{i}\lambda_{i,j}(1-\alpha_{i,j})(\boldsymbol{u}_{j}-\boldsymbol{u}_{i}),$
(20a) $\displaystyle\frac{\partial E_{i}}{\partial t}$
$\displaystyle=d\sum_{j}n_{i}\lambda_{i,j}(1-\beta_{i,j})(T_{j}-T_{i})$ (20b)
$\displaystyle\quad+\\!\sum_{j}\bigg{[}(\delta_{i,j}-n_{i}\lambda_{i,j}\gamma_{i,j})\boldsymbol{u}_{i}\cdot(\boldsymbol{u}_{j}-\boldsymbol{u}_{i})-(\delta_{j,i}-n_{j}\lambda_{j,i}\gamma_{j,i})\boldsymbol{u}_{j}\cdot(\boldsymbol{u}_{i}-\boldsymbol{u}_{j})\bigg{]}.$

Using (11) and (20a), the following ODE describing the temperature can be
derived:

$\frac{d}{2}\frac{\partial(n_{i}T_{i})}{\partial
t}=d\sum_{j}n_{i}\lambda_{i,j}(1-\beta_{i,j})(T_{j}-T_{i})+\sum_{j}n_{i}\lambda_{i,j}\gamma_{i,j}|\boldsymbol{u}_{i}-\boldsymbol{u}_{j}|^{2}.$
(21)

######  Lemma 3.7.

The ODE system (20) satisfies the conservation of total momentum and total
energy:

$\sum_{i}\frac{\partial(\rho_{i}\boldsymbol{u}_{i})}{\partial
t}=\mathbf{0}\qquad\mbox{and}\qquad\sum_{i}\frac{\partial E_{i}}{\partial
t}=0,$ (22)

and the quantities

$\boldsymbol{u}^{\infty}=\frac{\sum_{i}\rho_{i}\boldsymbol{u}_{i}}{\sum_{i}\rho_{i}}\in\mathbb{R}^{d}\qquad\mbox{and}\qquad
T^{\infty}=\frac{\sum_{i}n_{i}T_{i}}{\sum_{i}n_{i}}+\frac{\sum_{i}\rho_{i}\left(|\boldsymbol{u}_{i}|^{2}-|\boldsymbol{u}^{\infty}|^{2}\right)}{d\sum_{i}n_{i}}>0$
(23)

are independent of time.

###### Proof.

To verify the conservation of total momentum, sum (20a) over all species
$i\in\\{1,\cdots,N\\}$ and use the symmetry in (15a) to obtain

$\sum_{i}\frac{\partial(\rho_{i}\boldsymbol{u}_{i})}{\partial
t}=\sum_{i}\sum_{j}\rho_{i}\lambda_{i,j}(1-\alpha_{i,j})(\boldsymbol{u}_{j}-\boldsymbol{u}_{i})=\mathbf{0}.$
(24)

Next, to verify the conservation of total energy, sum (20b) over
$i\in\\{1,\cdots,N\\}$ and use (15b). Finally, the time invariance of
$\boldsymbol{u}^{\infty}$ and $T^{\infty}$ follow directly from (22), and the
positivity of $T_{\infty}$ follows as in [10, Proposition 4.1]. ∎

######  Remark 3.8.

Following arguments similar, to that of [10, Theorem 3.2] one can show that
the temperatures $T_{i}$, $i\in\\{1,\cdots,N\\}$, are bounded below and that
$min_{1\leq i\leq N}{T_{i}}$ is non-decreasing. Using this fact, the global
existence and uniqueness of solutions to the moment system and the convergence
for each $i\in\\{1,\cdots,N\\}$ of
$\boldsymbol{u}_{i}\to\boldsymbol{u}^{\infty}$ and $T_{i}\to T^{\infty}$ as
$t\to\infty$ can be established, in a way similar to [10].

##  4 Entropy dissipation and $\mathcal{H}$-Theorem

In this section, we state and prove an entropy dissipation result and an
$\mathcal{H}$-Theorem, under the following physically motivated assumption.

######  Assumption 4.1.

Motivated by physical principles, we assume that for $1\leq i,j\leq N$,

$0\leq\alpha_{i,j}\,,\,\beta_{i,j}\leq 1\qquad\mbox{and}\qquad\gamma_{i,j}\geq
0.$ (25)

######  Definition 4.2.

The total entropy functional for a particle system with kinetic distribution
$\boldsymbol{f}=[f_{1},\cdots f_{N}]^{\top}$ is given by

$\mathcal{H}[\boldsymbol{f}]=\sum_{i=1}^{N}{\int_{\mathbb{R}^{d}}}(f_{i}\log(f_{i})-f_{i})\,\,\textnormal{d}\boldsymbol{v}=\sum_{i=1}^{N}\langle\\!\;f_{i}\log(f_{i})-f_{i}\\!\;\rangle.$
(26)

######  Theorem 4.3.

Let $\boldsymbol{f}=[f_{1},\cdots f_{N}]^{\top}$ satisfy the M-LB system
(6)-(7). Then under the assumptions in (25),
$\frac{\partial\mathcal{H}}{\partial t}\leq 0$, with equality if and only if
each $f_{i}$ is a Maxwellian with common equilibrium velocity
$\boldsymbol{u}\in\mathbb{R}^{d}$ and equilibrium temperature $T>0$; that is,

$\displaystyle f_{i}(\boldsymbol{v})$
$\displaystyle=n_{i}\left(\frac{m_{i}}{2\pi
T^{\infty}}\right)^{d/2}\exp\left(-\frac{m_{i}|\boldsymbol{v}-\boldsymbol{u}^{\infty}|^{2}}{2T^{\infty}}\right),\qquad
1\leq i\leq N,$ (27)

where $\boldsymbol{u}^{\infty}$ and $T^{\infty}$ are defined in (23).

###### Proof.

By the definition of $\mathcal{H}$ in (26) and the chain rule,

$\begin{split}\frac{\partial\mathcal{H}}{\partial
t}&=\sum_{i=1}^{N}\langle\\!\;\log(f_{i})\frac{\partial f_{i}}{\partial
t}\\!\;\rangle\\\
&=\sum_{i=1}^{N}\langle\\!\;\log(f_{i})\mathcal{C}_{i,i}(f_{i},f_{i})\\!\;\rangle+\sum_{i=1}^{N}\sum_{j=i+1}^{N}\langle\\!\;\log(f_{i})\mathcal{C}_{i,j}(f_{i},f_{j})+\log(f_{j})\mathcal{C}_{j,i}(f_{j},f_{i})\\!\;\rangle.\end{split}$
(28)

Thus it is sufficient to show that

$\langle\\!\;\log(f_{i})\mathcal{C}_{i,i}(f_{i},f_{i})\\!\;\rangle\leq 0,$
(29)

with equality if and only if $f_{i}$ takes the form in (27), and

$\langle\\!\;\log(f_{i})\mathcal{C}_{i,j}(f_{i},f_{j})+\log(f_{j})\mathcal{C}_{j,i}(f_{j},f_{i})\\!\;\rangle\leq
0,$ (30)

with equality if and only if $f_{i}$ and $f_{j}$ take the Maxwellian form in
(27) with a common bulk velocity $\boldsymbol{u}$ and temperature $T>0$. We
establish these conditions in the remainder of the section. ∎

######  Lemma 4.4.

For the intra-species collisions,

$\langle\\!\;\log(f_{i})\mathcal{C}_{i,i}(f_{i},f_{i})\\!\;\rangle\leq 0,$
(31)

with equality if and only if $f_{i}=M_{i}$.

###### Proof.

Since

$\log(M_{i})=\log\left(\frac{n_{i}}{(2\pi\theta_{i})^{d/2}}\right)-\frac{|\boldsymbol{v}-\boldsymbol{u}_{i}|^{2}}{2\theta_{i}},$
(32)

it follows from Lemma 3.3 that

$\left\langle\log(M_{i})\mathcal{C}(f_{i},f_{i})\right\rangle=0.$ (33)

Hence, after integrating by parts,

$\begin{split}\left\langle\log(f_{i})\mathcal{C}_{i,i}(f_{i},f_{i})\right\rangle&=\left\langle\log\left(\frac{f_{i}}{M_{i}}\right)\mathcal{C}_{i,i}(f_{i},f_{i})\right\rangle\\\
=\lambda_{i,i}\theta_{i}&\left\langle\log\left(\frac{f_{i}}{M_{i}}\right)\nabla_{\boldsymbol{v}}\cdot\left(M_{i}\nabla_{\boldsymbol{v}}\left(\frac{f}{M_{i}}\right)\right)\right\rangle=-\lambda_{i,i}\theta_{i}\left\langle\frac{M_{i}^{2}}{f}\left|\nabla_{\boldsymbol{v}}\left(\frac{f_{i}}{M_{i}}\right)\right|^{2}\right\rangle.\end{split}$
(34)

Thus the inequality in (31) holds, with equality if and only if $f_{i}/M_{i}$
is constant with respect to $\boldsymbol{v}$. Since
$\langle\\!\;f_{i}\\!\;\rangle=\langle\\!\;M_{i}\\!\;\rangle$, this constant
must be one, and the result follows. ∎

We now prove (30). Several intermediate steps are required.

######  Lemma 4.5.

For a general Maxwellian $M_{n,\boldsymbol{u},\theta}$,

$\langle\\!\;\log(f)\mathcal{L}(f,M_{n,\boldsymbol{u},\theta})\\!\;\rangle=-\langle\\!\;f\Delta_{\boldsymbol{v}}\left(\log(f)-\log(M_{n,\boldsymbol{u},\theta})\right)\\!\;\rangle,$
(35)

which implies that

$\langle\\!\;\log(f_{i})\mathcal{L}(f_{i},M_{i,j})\\!\;\rangle-\langle\\!\;\log(f_{i})\mathcal{L}(f_{i},M_{i,i})\\!\;\rangle=\langle\\!\;f_{i}\Delta_{\boldsymbol{v}}\left(\log(M_{i,i})-\log(M_{i,j})\right)\\!\;\rangle.$
(36)

###### Proof.

Using integration by parts,

$\begin{split}\langle\\!\;\log(f)\mathcal{L}(f,M)\\!\;\rangle&=\langle\\!\;\log(f)\nabla_{\boldsymbol{v}}\cdot\left[\nabla_{\boldsymbol{v}}f-f\nabla_{\boldsymbol{v}}(\log(M))\right]\\!\;\rangle\\\
&=-\langle\\!\;\nabla_{\boldsymbol{v}}\log(f)\cdot\nabla_{\boldsymbol{v}}f\\!\;\rangle+\langle\\!\;f\nabla_{\boldsymbol{v}}\log(f)\cdot\nabla_{\boldsymbol{v}}\log(M)\\!\;\rangle\\\
&=-\langle\\!\;\nabla_{\boldsymbol{v}}\log(f)\cdot\nabla_{\boldsymbol{v}}f\\!\;\rangle+\langle\\!\;\nabla_{\boldsymbol{v}}f\cdot\nabla_{\boldsymbol{v}}\log(M)\\!\;\rangle\\\
&=\langle\\!\;\nabla_{\boldsymbol{v}}f\cdot\nabla_{\boldsymbol{v}}\left(\log(M)-\log(f)\right)\\!\;\rangle\\\
&=-\langle\\!\;f\Delta_{\boldsymbol{v}}\left(\log(f)-\log(M)\right)\\!\;\rangle.\end{split}$
(37)

Identity (36) follows immediately. ∎

The key point of Lemma 4.5 is that the entropy dissipation due to
$\mathcal{L}(f_{i},M_{i,j})$ can be expressed in terms of moments of $f_{i}$
and $f_{j}$. This fact is used to write the inter-species entropy dissipation
in terms of the intra-species entropy dissipation and some moment-dependent
corrections.

######  Lemma 4.6.

$\displaystyle\langle\\!\;\log(f_{i})\mathcal{C}_{i,j}(f_{i},f_{j})\\!\;\rangle$
$\displaystyle=\frac{\lambda_{i,j}\theta_{i,j}}{\lambda_{i,i}\theta_{i,i}}\langle\\!\;\log(f_{i})\mathcal{C}_{i,i}(f_{i},f_{i})\\!\;\rangle+dn_{i}\lambda_{i,j}\frac{T_{i}-T_{i,j}}{T_{i}}.$
(38)

###### Proof.

Adding and subtracting
$\lambda_{i,j}\theta_{i,j}\langle\\!\;\log(f_{i})\mathcal{L}(f_{i},M_{i,j})\\!\;\rangle$
on the left-hand side of (38) and then applying Lemma 4.5 gives

$\displaystyle\langle\\!\;\log(f_{i})\mathcal{C}_{i,j}(f_{i},f_{j})\\!\;\rangle$
$\displaystyle=\lambda_{i,j}\theta_{i,j}\langle\\!\;\log(f_{i})\mathcal{L}(f_{i},M_{i,j})\\!\;\rangle$
$\displaystyle=\lambda_{i,j}\theta_{i,j}\\!\left\langle\log(f_{i})\mathcal{L}(f_{i},M_{i,i})\right\rangle\\!+\\!\lambda_{i,j}\theta_{i,j}\\!\left\langle\log(f_{i})\mathcal{L}(f_{i},M_{i,j})\right\rangle\\!-\\!\lambda_{i,j}\theta_{i,j}\\!\left\langle\log(f_{i})\mathcal{L}(f_{i},M_{i,i})\right\rangle$
$\displaystyle=\frac{\lambda_{i,j}\theta_{i,j}}{\lambda_{i,i}\theta_{i,i}}\langle\\!\;\log(f_{i})\mathcal{C}_{i,i}(f_{i},f_{i})\\!\;\rangle+\lambda_{i,j}\theta_{i,j}\left\langle
f_{i}\Delta_{\boldsymbol{v}}\left[\log(M_{i,i})-\log(M_{i,j})\right]\right\rangle$
$\displaystyle=\frac{\lambda_{i,j}\theta_{i,j}}{\lambda_{i,i}\theta_{i,i}}\langle\\!\;\log(f_{i})\mathcal{C}_{i,i}(f_{i},f_{i})\\!\;\rangle+dn_{i}\lambda_{i,j}\frac{T_{i}-T_{i,j}}{T_{i}},$
(39)

where in the final line above we have used (5) and (8). ∎

######  Lemma 4.7.

The cross terms in (30) satisfy

$\displaystyle\langle\\!\;\log(f_{i})\mathcal{C}_{i,j}(f_{i},f_{j})\\!\;\rangle+\langle\\!\;\log(f_{j})\mathcal{C}_{j,i}(f_{j},f_{i})\\!\;\rangle\leq
0,$ (40)

with equality if and only if $f_{i}$ and $f_{j}$ are Maxwellians with a common
bulk velocity and temperature.

###### Proof.

Using (9b), (15b), and Lemma 4.6,

$\begin{split}\langle\\!\;\log(f_{i})\mathcal{C}_{i,j}(f_{i},f_{j})\\!\;\rangle+\langle\\!\;\log(f_{j})\mathcal{C}_{j,i}(f_{j},f_{i})\\!\;\rangle&=\frac{\lambda_{i,j}\theta_{i,j}}{\lambda_{i,i}\theta_{i,i}}\langle\\!\;\log(f_{i})\mathcal{C}_{i,i}(f_{i},f_{i})\\!\;\rangle+\frac{\lambda_{j,i}\theta_{j,i}}{\lambda_{j,j}\theta_{j,j}}\langle\\!\;\log(f_{j})\mathcal{C}_{j,j}(f_{j},f_{j})\\!\;\rangle\\\
-dn_{i}\lambda_{i,j}(1-\beta_{i,j})&\frac{(T_{i}-T_{j})^{2}}{T_{i}T_{j}}-\left(n_{i}\lambda_{i,j}\gamma_{i,j}T_{j}+n_{j}\lambda_{j,i}\gamma_{j,i}T_{i}\right)\frac{|\boldsymbol{u}_{i}-\boldsymbol{u}_{j}|^{2}}{T_{i}T_{j}}.\end{split}$
(41)

By Lemma 4.4, the right hand side of (41) is non-positive. Moreover, it is
equal to zero if and only if $f_{i}$ and $f_{j}$ are Maxwellians with
$\boldsymbol{u}_{i}=\boldsymbol{u}_{j}$ and $T_{i}=T_{j}$. ∎

######  Remark 4.8.

While the condition on $\alpha_{i,j}\in[0,1]$ is not explicitly used in the
proofs above, the condition that $\gamma_{i,j}\geq 0$, which is used, may fail
to hold if the bound on $\alpha_{i,j}$ does not hold. In particular, (15a) and
(15c) implies that

$\rho_{i}\lambda_{i,j}(1-\alpha_{i,j})=n_{i}\lambda_{i,j}\gamma_{i,j}+n_{j}\lambda_{j,i}\gamma_{j,i}.$
(42)

Thus if $\alpha_{i,j}>1$, either $\gamma_{i,j}$ or $\gamma_{j,i}$ will fail to
be positive.

##  5 Matching relaxation rates

In this section, momentum and energy relaxation rates for the current model
are matched to the relaxation rates for the Boltzmann collision operator with
Coulomb potential (Boltz-C). To derive relaxation rates between species, it is
standard to consider a two species mixture. In this setting [7, 17],

$\displaystyle\left.\frac{\partial}{\partial
t}\left(\rho_{i}\boldsymbol{u}_{i}-\rho_{j}\boldsymbol{u}_{j}\right)\right|_{\textnormal{(Boltz-C)}}$
$\displaystyle=\xi_{i,j}(m_{i}+m_{j})(\boldsymbol{u}_{j}-\boldsymbol{u}_{i}),$
(43a) $\displaystyle\left.\frac{\partial}{\partial
t}\frac{d}{2}\left(n_{i}T_{i}-n_{j}T_{j}\right)\right|_{\textnormal{(Boltz-C)}}$
$\displaystyle=\xi_{i,j}\left[d(T_{j}-T_{i})+\frac{m_{j}-m_{i}}{2}|\boldsymbol{u}_{i}-\boldsymbol{u}_{j}|^{2}\right],$
(43b)

where $\xi_{i,j}$ depends on the vacuum permittivity $\epsilon_{0}$, the
Coulomb logarithm $\log(\Lambda_{i,j})$ and the species charges $q_{i}$ and
$q_{j}$: 222The term $\xi_{i,j}\coloneqq\alpha_{E}$, where $\alpha_{E}$ is
taken from Equation (2.5) of [7].

$\xi_{i,j}=\frac{2}{3(2\pi)^{3/2}\epsilon_{0}^{2}}\frac{|\log\Lambda_{i,j}|(q_{i}q_{j})^{2}n_{i}n_{j}}{m_{i}m_{j}\left(\frac{T_{i}}{m_{i}}+\frac{T_{j}}{m_{j}}\right)^{3/2}}.$
(44)

We will match these rates with the relaxation rates of the M-LB model:

$\displaystyle\left.\frac{\partial}{\partial
t}\left(\rho_{i}\boldsymbol{u}_{i}-\rho_{j}\boldsymbol{u}_{j}\right)\right|_{\textnormal{M-LB}}$
$\displaystyle=2\rho_{i}\lambda_{i,j}(1-\alpha_{i,j})(\boldsymbol{u}_{j}-\boldsymbol{u}_{i}),$
(45a) $\displaystyle\left.\frac{\partial}{\partial
t}\frac{d}{2}\left(n_{i}T_{i}-n_{j}T_{j}\right)\right|_{\textnormal{M-LB}}$
$\displaystyle=2dn_{i}\lambda_{i,j}(1-\beta_{i,j})(T_{j}-T_{i})+(2n_{i}\lambda_{i,j}\gamma_{i,j}-\delta_{i,j})|\boldsymbol{u}_{i}-\boldsymbol{u}_{j}|^{2}.$
(45b)

The derivation of these relaxation rates is given in Appendix B; specifically
(45a) is verified in Lemma B.1 and (45b) is verified in Lemma B.2.

Equating the right-hand sides of (45) and (43) gives the following matching
conditions for the momentum:

$2\rho_{i}\lambda_{i,j}(1-\alpha_{i,j})=\xi_{i,j}(m_{i}+m_{j}),$ (46)

and for the temperature:

$2dn_{i}\lambda_{i,j}(1-\beta_{i,j})=d\xi_{i,j}\qquad\mbox{and}\qquad
2n_{i}\lambda_{i,j}\gamma_{i,j}-\delta_{i,j}=\xi_{i,j}\frac{m_{j}-m_{i}}{2}.$
(47)

######  Theorem 5.1.

For $1\leq i,j\leq N$, the relaxation rates in (45) match the relaxation rates
in (43), with the following expressions for $\alpha_{i,j}$, $\beta_{i,j}$, and
$\gamma_{i,j}$:

$\displaystyle\alpha_{i,j}$
$\displaystyle=1-\frac{1}{2}\frac{m_{i}+m_{j}}{m_{i}}\frac{\xi_{i,j}}{n_{i}\lambda_{i,j}},$
(48a) $\displaystyle\beta_{i,j}$
$\displaystyle=1-\frac{1}{2}\frac{\xi_{i,j}}{n_{i}\lambda_{i,j}}=1-\frac{m_{i}}{m_{i}+m_{j}}(1-\alpha_{i,j}),$
(48b) $\displaystyle\gamma_{i,j}$
$\displaystyle=\frac{1}{2}\frac{1}{n_{i}\lambda_{i,j}}\left[\xi_{i,j}\frac{m_{j}-m_{i}}{2}+\delta_{i,j}\right]=\frac{m_{i}m_{j}}{m_{i}+m_{j}}(1-\alpha_{i,j}).$
(48c)

###### Proof.

The equations in (48) follow immediately by solving (46) and (47) for
$\alpha_{i,j}$, $\beta_{i,j}$, and $\gamma_{i,j}$. ∎

######  Remark 5.2.

The matching conditions between the coefficients in (45) and the coefficients
in (43) are consistent with the relationship between $\alpha_{i,j}$ and
$\alpha_{j,i}$ given in Proposition 3.4. Thus it is sufficient to perform the
matching only for $i<j$. More specifically, if for all $1\leq i<j\leq N$, it
holds that

$2\rho_{i}\lambda_{i,j}(1-\alpha_{i,j})=\xi_{i,j}(m_{i}+m_{j}),$ (49)

then (15a) and the symmetry of $\xi_{i,j}$ imply that

$2\rho_{j}\lambda_{j,i}(1-\alpha_{j,i})=2\rho_{i}\lambda_{i,j}(1-\alpha_{i,j})=\xi_{i,j}(m_{i}+m_{j})=\xi_{j,i}(m_{j}+m_{i}).$
(50)

That is, the matching condition holds for the momentum relaxation with $i>j$.
Similarly, (15b) and the symmetry of $\xi_{i,j}$ imply that

$d\xi_{i,j}=2dn_{i}\lambda_{i,j}(1-\beta_{i,j})=2dn_{j}\lambda_{j,i}(1-\beta_{j,i})=d\xi_{j,i},$
(51)

while (15c), and the symmetry of $\xi_{i,j}$ and $\delta_{i,j}$ imply that

$\xi_{i,j}\frac{m_{j}-m_{i}}{2}=2n_{i}\lambda_{i,j}\gamma_{i,j}-\delta_{i,j}=-(2n_{j}\lambda_{j,i}\gamma_{j,i}-\delta_{j,i})=-\xi_{j,i}\frac{m_{i}-m_{j}}{2}.$
(52)

Based on (48), it is clear that for $1\leq i,j\leq N$, the parameters
$\beta_{i,j}$ and $\gamma_{i,j}$ are completely determined by $\alpha_{i,j}$.
Moreover, we have the following

######  Lemma 5.3.

Suppose that for $1\leq i,j\leq N$, $0\leq\alpha_{i,j}\leq 1$. Then the
remaining conditions in Assumption 4.1 hold. That is

$0\leq\beta_{i,j}\leq 1\qquad\mbox{and}\qquad\gamma_{i,j}\geq 0,\qquad 1\leq
i,j\leq N.$ (53)

###### Proof.

The proof follows immediately from (48b) and (48c) and the assumed bounds on
$\alpha_{i,j}$. ∎

###  5.1 Choosing collision frequencies

What remains to complete the model is to find frequencies $\lambda_{i,j}$ such
that parameters $\alpha_{i,j}$ defined in (48a) satisfy the relation (15a) and
the condition $0\leq\alpha_{i,j}\leq 1$. We propose collision frequencies of
the form

$\lambda_{i,j}=\widehat{\lambda}_{i,j}(T_{i},T_{j})=\frac{\kappa_{i,j}\xi_{i,j}}{n_{i}},$
(54)

where $\xi_{i,j}$ is given in (44) and $\kappa_{i,j}\geq 0$ is a tunable
parameter. Inserting (54) into (48a) gives

$\alpha_{i,j}=1-\frac{1}{2}\frac{m_{i}+m_{j}}{m_{i}}\frac{1}{\kappa_{i,j}},$
(55)

which implies that $\alpha_{i,j}\leq 1$, and

$\alpha_{i,j}\geq
0\quad\iff\quad\kappa_{i,j}\geq\frac{m_{i}+m_{j}}{2m_{i}}\eqqcolon\mu_{i,j}.$
(56)

Below we show how to enforce the desired conditions on $\alpha_{i,j}$ based on
the choice of $\kappa_{i,j}$.

######  Theorem 5.4.

Suppose that the collision frequencies are given by (54), and
$\kappa_{i,j}\geq\mu_{i,j}$ for each $1\leq i,j\leq N$. Then for $1\leq
i,j\leq N$, $0\leq\alpha_{i,j}\leq 1$, and for $1\leq i<j\leq N$, the pairs
$(\alpha_{i,j},\alpha_{j,i})$, $(\beta_{i,j},\beta_{j,i})$, and
$(\gamma_{i,j},\gamma_{j,i})$ given in (48) satisfy (15).

###### Proof.

The bound $0\leq\alpha_{i,j}\leq 1$ is clear from (55) and (56). Moreover
rearranging (55) gives

$1-\alpha_{i,j}=\frac{\mu_{i,j}}{\kappa_{i,j}}.$ (57)

Using (54) and (57), the left-hand side of (15a) is

$\rho_{i}\lambda_{i,j}(1-\alpha_{i,j})=\rho_{i}\frac{\kappa_{i,j}\xi_{i,j}}{n_{i}}\frac{\mu_{i,j}}{\kappa_{i,j}}=\xi_{i,j}\frac{m_{i}+m_{j}}{2},$
(58)

and the right-hand side is

$\rho_{j}\lambda_{j,i}(1-\alpha_{j,i})=\rho_{j}\frac{\kappa_{j,i}\xi_{j,i}}{n_{j}}\frac{\mu_{j,i}}{\kappa_{j,i}}=\xi_{j,i}\frac{m_{j}+m_{i}}{2}.$
(59)

Since $\xi_{i,j}=\xi_{j,i}$, the symmetry in (15a) holds.

Next, using (48b), (54), and (57),

$1-\beta_{i,j}=\frac{1}{2\kappa_{i,j}},$ (60)

so that the left-hand side of (15b) is

$n_{i}\lambda_{i,j}(1-\beta_{i,j})=n_{i}\frac{\kappa_{i,j}\xi_{i,j}}{n_{i}}\frac{1}{2\kappa_{i,j}}=\frac{\xi_{i,j}}{2},$
(61)

and the right-hand side is

$n_{i}\lambda_{j,i}(1-\beta_{j,i})=n_{j}\frac{\kappa_{j,i}\xi_{j,i}}{n_{j}}\frac{1}{2\kappa_{j,i}}=\frac{\xi_{j,i}}{2}.$
(62)

Since $\xi_{i,j}=\xi_{j,i}$, the symmetry in (15b) holds.

Finally, using (48c), (54), and (57),

$\gamma_{i,j}=\frac{m_{j}}{2\kappa_{i,j}},$ (63)

so that the left-hand side of (15c) is

$\delta_{i,j}-2n_{i}\lambda_{i,j}\gamma_{i,j}=m_{i}n_{i}\frac{\kappa_{i,j}\xi_{i,j}}{n_{i}}\frac{\mu_{i,j}}{\kappa_{i,j}}-2n_{i}\frac{\kappa_{i,j}\xi_{i,j}}{n_{i}}\frac{m_{j}}{2\kappa_{i,j}}=\xi_{i,j}\frac{m_{i}-m_{j}}{2},$
(64)

and the right-hand side is

$2n_{j}\lambda_{j,i}\gamma_{j,i}-\delta_{j,i}=2n_{j}\frac{\kappa_{j,i}\xi_{j,i}}{n_{j}}\frac{m_{i}}{2\kappa_{j,i}}-m_{j}n_{j}\frac{\kappa_{j,i}\xi_{j,i}}{n_{j}}\frac{\mu_{j,i}}{\kappa_{j,i}}=\xi_{j,i}\frac{m_{i}-m_{j}}{2}.$
(65)

Since $\xi_{i,j}=\xi_{j,i}$, the symmetry in (15c) holds. ∎

######  Corollary 5.4.1.

For every $1\leq i,j\leq N$, let $\alpha_{i,j}$, $\beta_{i,j}$, and
$\gamma_{i,j}$ be given by (48), and let $\lambda_{i,j}$ be given by (54),
where $\kappa_{i,j}$ satisfies (56). Then (i) $\alpha_{i,j}$, $\beta_{i,j}$,
and $\gamma_{i,j}$ satisfy Assumption 4.1. Furthermore, any solution
$\boldsymbol{f}=[f_{i},\cdots,f_{N}]^{\top}$ of the M-LB model (6)-(7) (ii)
satisfies the pairwise-wise conservation laws (Assumption 3.2), (iii)
satisfies an entropy dissipation law and an $\mathcal{H}$-Theorem (Theorem
4.3), and (iv) matches the pairwise momentum and temperature relaxation rates
derived from the Boltzmann collision operator with Coulomb potential (Theorem
5.1).

###### Proof.

Claim (i) is a restatement of Lemma 5.3, along with the bound in (56). Claim
(ii) follows from Theorem 5.4 and the equivalence of Assumption 3.2 and (12)
(see Proposition 3.4). Claims (iii) and (iv) are summaries of the theorems
referenced in the statement of the corollary. ∎

##  6 Numerics

For simulation purposes, we propose an implicit Euler update of (6) which
takes the form (cf. (7))

$\frac{f_{i}^{\mathfrak{n}+1}-f_{i}^{\mathfrak{n}}}{\Delta
t}=\sum_{j}\lambda_{i,j}^{\mathfrak{n}+1}\theta_{i,j}^{\mathfrak{n}+1}\mathcal{L}_{i,j}^{\mathfrak{n}+1},\;\;\textnormal{where}\;\;\mathcal{L}_{i,j}^{\mathfrak{n}+1}=\nabla_{\boldsymbol{v}}\cdot\left(M_{i,j}^{\mathfrak{n}+1}\nabla_{\boldsymbol{v}}\left(\frac{f_{i}^{\mathfrak{n}+1}}{M_{i,j}^{\mathfrak{n}+1}}\right)\right).$
(66)

Because $\lambda_{i,j}^{\mathfrak{n}+1}$, $\theta_{i,j}^{\mathfrak{n}+1}$, and
the mixture Maxwellians, $M_{i,j}^{\mathfrak{n}+1}$ in
$\mathcal{L}^{\mathfrak{n}+1}_{i,j}$ are defined in terms of the moment
quantities $\rho_{i}^{\mathfrak{n}+1}$, $\boldsymbol{u}_{i}^{\mathfrak{n}+1}$,
and $T_{i}^{\mathfrak{n}+1}$, it is useful to first derive and solve a system
of equations for these moments. What remains afterward is the inversion of a
linear tridiagonal system, the coefficients of which are expressed in terms of
the updated moments.

###  6.1 Implicit moment update

To derive a set of moment equations for $\rho^{\mathfrak{n}}_{i}$,
$\boldsymbol{u}^{\mathfrak{n}}_{i}$, and $T_{i}^{\mathfrak{n}}$, we multiply
(66) by $m_{i}(1,\boldsymbol{v},|\boldsymbol{v}|^{2})^{\top}$ and integrate
over $\boldsymbol{v}\in\mathbb{R}^{d}$. Then, applying Lemma 3.3 and the
definitions of $\boldsymbol{u}_{i,j}$ and $T_{i,j}$ in (9), gives the
following Backward Euler update

$\displaystyle\frac{\rho_{i}^{\mathfrak{n}+1}-\rho_{i}^{\mathfrak{n}}}{\Delta
t}$ $\displaystyle=0,$ (67a)
$\displaystyle\rho_{i}\frac{\boldsymbol{u}_{i}^{\mathfrak{n}+1}-\boldsymbol{u}_{i}^{\mathfrak{n}}}{\Delta
t}$
$\displaystyle=\sum_{j}\rho_{i}\lambda_{i,j}^{\mathfrak{n}+1}(1-\alpha_{i,j})(\boldsymbol{u}_{j}^{\mathfrak{n}+1}-\boldsymbol{u}_{i}^{\mathfrak{n}+1}),$
(67b)
$\displaystyle\begin{split}\rho_{i}\frac{s_{i}^{\mathfrak{n}+1}-s_{i}^{\mathfrak{n}}}{\Delta
t}+dn_{i}\frac{T_{i}^{\mathfrak{n}+1}-T_{i}^{\mathfrak{n}}}{\Delta
t}&=2d\sum_{j}n_{i}\lambda_{i,j}^{\mathfrak{n}+1}(1-\beta_{i,j})(T_{j}^{\mathfrak{n}+1}-T_{i}^{\mathfrak{n}+1})\\\
+2\sum_{j}n_{i}\lambda_{i,j}^{\mathfrak{n}+1}(1-\beta_{i,j}&)\big{[}m_{j}(s_{j}^{\mathfrak{n}+1}-\boldsymbol{u}_{i}^{\mathfrak{n}+1}\cdot\boldsymbol{u}_{j}^{\mathfrak{n}+1})-m_{i}(s_{i}^{\mathfrak{n}+1}-\boldsymbol{u}_{i}^{\mathfrak{n}+1}\cdot\boldsymbol{u}_{j}^{\mathfrak{n}+1})\big{]},\end{split}$
(67c)

where $s_{i}^{\mathfrak{n}}=|\boldsymbol{u}_{i}^{\mathfrak{n}}|^{2}$, and we
drop the superscripts on $n_{i}$ and $\rho_{i}$ due to (67a). The system given
here is consistent with the backward Euler step of the ODE system (20), when
the matching conditions in (48) hold. Summing (67b) and (67c) over
$i\in\\{1,\cdots,N\\}$ and using the symmetries (15a) and (15b) implies that
the discrete versions of the conservation of total momentum and total energy
hold:

$\sum_{i}\rho_{i}\boldsymbol{u}_{i}^{\mathfrak{n}+1}=\sum_{i}\rho_{i}\boldsymbol{u}_{i}^{\mathfrak{n}}\qquad\mbox{and}\qquad\sum_{i}E_{i}^{\mathfrak{n}+1}=\sum_{i}E_{i}^{\mathfrak{n}}.$
(68)

The form of the implicit update in (67) is almost identical to the moment
system in [11] derived for a multi-species BGK model. The only difference is
the factor of two in the energy equation and the appearance of
$\boldsymbol{u}_{i}^{\mathfrak{n}+1}\cdot\boldsymbol{u}_{j}^{\mathfrak{n}+1}$
instead of $|\boldsymbol{u}_{i,j}^{\mathfrak{n}+1}|^{2}$. For this reason, we
adapt the nonlinear Gauss-Seidel type (GST) iterative method of [11] to solve
(67). Given iteration index $\ell$, the iterative update for the momentum and
energy equations within each time step takes the form

$\displaystyle\rho_{i}\frac{\boldsymbol{u}_{i}^{\mathfrak{n}+1,\ell+1}-\boldsymbol{u}_{i}^{\mathfrak{n}}}{\Delta
t}=\sum_{j}\rho_{i}\lambda_{i,j}^{\mathfrak{n}+1,\ell}(1-\alpha_{i,j})(\boldsymbol{u}_{j}^{\mathfrak{n}+1,\ell+1}-\boldsymbol{u}_{i}^{\mathfrak{n}+1,\ell+1})$
(69a)
$\displaystyle\rho_{i}\frac{s_{i}^{\mathfrak{n}+1,\ell+1}\\!-s_{i}^{\mathfrak{n}}}{\Delta
t}+dn_{i}\frac{T_{i}^{\mathfrak{n}+1,\ell+1}\\!-T_{i}^{\mathfrak{n}}}{\Delta
t}=2d\sum_{j}n_{i}\lambda_{i,j}^{\mathfrak{n}+1,\ell}(1-\beta_{i,j})(T_{j}^{\mathfrak{n}+1,\ell+1}\\!-T_{i}^{\mathfrak{n}+1,\ell+1})$
(69b)
$\displaystyle+2\sum_{j}n_{i}\lambda_{i,j}^{\mathfrak{n}+1,\ell}(1-\beta_{i,j})\bigg{[}m_{j}(s_{j}^{\mathfrak{n}+1,\ell+1}-\boldsymbol{u}_{i}^{\mathfrak{n}+1,\ell+1}\cdot\boldsymbol{u}_{j}^{\mathfrak{n}+1,\ell+1})-m_{i}(s_{i}^{\mathfrak{n}+1,\ell+1}-\boldsymbol{u}_{i}^{\mathfrak{n}+1,\ell+1}\cdot\boldsymbol{u}_{j}^{\mathfrak{n}+1,\ell+1})\bigg{]},$

where

$\lambda_{i,j}^{\mathfrak{n}+1,\ell}=\widehat{\lambda}_{i,j}(T_{i}^{\mathfrak{n}+1,\ell},T_{j}^{\mathfrak{n}+1,\ell}),$
(70)

and $\widehat{\lambda}_{i,j}$ is given by (54) and (44). Once
$\boldsymbol{u}^{\mathfrak{n}+1}_{i}$ and $T^{\mathfrak{n}+1}_{i}$ have been
evaluated using the GST method (69), they can be used to update the mixture
terms $\boldsymbol{u}^{\mathfrak{n}+1}_{i,j}$ and $T^{\mathfrak{n}+1}_{i,j}$
via (9), as well as the Maxwellian terms $M_{i,j}^{\mathfrak{n}+1}$ via (8).

###  6.2 Collision term update

After $\boldsymbol{u}_{i}^{\mathfrak{n}+1}$ and $T_{i}^{\mathfrak{n}+1}$ are
determined, (66) is a linear equation. For simplicity, we focus on the case
$d=1$, but the discretization below can be extended to a general dimension
$d$. We truncate the velocity space at some sufficiently large value,
$V=[-v_{\max},v_{\max}]$, and consider a cell centered velocity grid with
$N_{v}$ grid cells of uniform size $h_{v}=2v_{\max}/N_{v}$, with
$v_{k}=-v_{\max}+(k-1/2)h_{v}$ defined for $k$ at integer and half integer
values, so that $v_{1/2}=-v_{\max}$ and $v_{N_{v}+1/2}=v_{\max}$.

Approximating velocity derivatives in $\mathcal{L}_{i,j}$ by central
differences gives, for $k\in\\{1,\cdots,N_{v}\\}$,

$\begin{split}\left[\mathcal{L}_{i,j}\right]_{k}&=\left[\,\textnormal{d}_{v}\left(M_{i,j}\,\textnormal{d}_{v}\left(\frac{f_{i}}{M_{i,j}}\right)\right)\right]_{k}\\\
&\approx\frac{[M_{i,j}]_{k+\frac{1}{2}}[\,\textnormal{d}_{v}(f_{i}/M_{i,j})]_{k+\frac{1}{2}}-[M_{i,j}]_{k-\frac{1}{2}}[\,\textnormal{d}_{v}(f_{i}/M_{i,j})]_{k-\frac{1}{2}}}{\Delta
v}\\\
&\approx\frac{[M_{i,j}]_{k+\frac{1}{2}}\left[\left[\frac{f_{i}}{M_{i,j}}\right]_{k+1}-\left[\frac{f_{i}}{M_{i,j}}\right]_{k}\right]-[M_{i,j}]_{k-\frac{1}{2}}\left[\left[\frac{f_{i}}{M_{i,j}}\right]_{k}-\left[\frac{f_{i}}{M_{i,j}}\right]_{k-1}\right]}{\Delta
v^{2}}\\\
&=\frac{a_{k-1}[f_{i}]_{k-1}-b_{k}[f_{i}]_{k}+c_{k+1}[f_{i}]_{k+1}}{\Delta
v^{2}},\end{split}$ (71)

where

$a_{k-1}=\frac{[M_{i,j}]_{k-\frac{1}{2}}}{[M_{i,j}]_{k-1}},\qquad
b_{k}=\frac{[M_{i,j}]_{k-\frac{1}{2}}}{[M_{i,j}]_{k}}+\frac{[M_{i,j}]_{k+\frac{1}{2}}}{[M_{i,j}]_{k}},\qquad
c_{k+1}=\frac{[M_{i,j}]_{k+\frac{1}{2}}}{[M_{i,j}]_{k+1}},$ (72)

and the dependence of $a_{k}$, $b_{k}$, and $c_{k}$ on $i$ and $j$ has been
suppressed to simplify notation. If we set
$[M_{i,j}]_{k+\frac{1}{2}}=\frac{1}{2}([M_{i,j}]_{k}+[M_{i,j}]_{k+1})$, then

$a_{k-1}=\frac{[M_{i,j}]_{k-1}+[M_{i,j}]_{k}}{2[M_{i,j}]_{k-1}},\quad
b_{k}=\frac{[M_{i,j}]_{k-1}+2[M_{i,j}]_{k}+[M_{i,j}]_{k+1}}{2[M_{i,j}]_{k}},\quad
c_{k+1}=\frac{[M_{i,j}]_{k}+[M_{i,j}]_{k+1}}{2[M_{i,j}]_{k+1}}.$ (73)

To obtain no-flux boundary conditions (for velocity space), we set
$[M_{i,j}]_{\frac{1}{2}}=0=[M_{i,j}]_{N_{v}+\frac{1}{2}}$, which implies that

$a_{0}=\frac{[M_{i,j}]_{\frac{1}{2}}}{[M_{i,j}]_{0}}=0,\quad\\!b_{1}=\frac{[M_{i,j}]_{1+\frac{1}{2}}}{[M_{i,j}]_{1}},\quad\\!b_{N_{v}}=\frac{[M_{i,j}]_{N_{v}-\frac{1}{2}}}{[M_{i,j}]_{N_{v}}},\quad\\!c_{N_{v}+1}=\frac{[M_{i,j}]_{N_{v}+\frac{1}{2}}}{[M_{i,j}]_{N_{v}+1}}=0.$
(74)

Rearranging (66) and using (71) gives

$[f_{i}^{\mathfrak{n}+1}]_{k}\approx[f_{i}^{\mathfrak{n}}]_{k}+\frac{\Delta
t}{\Delta
v^{2}}\sum_{j}\lambda_{i,j}^{\mathfrak{n}+1}\theta_{i,j}^{\mathfrak{n}+1}\left(a_{k-1}^{\mathfrak{n}+1}[f_{i}]_{k-1}^{\mathfrak{n}+1}-b_{k}^{\mathfrak{n}+1}[f_{i}]_{k}^{\mathfrak{n}+1}+c_{k+1}^{\mathfrak{n}+1}[f_{i}]_{k+1}^{\mathfrak{n}+1}\right),$
(75)

or in vectorized form,

$\left(I+\frac{\Delta t}{\Delta
v^{2}}\sum_{j}\lambda_{i,j}^{\mathfrak{n}+1}\theta_{i,j}^{\mathfrak{n}+1}G_{i,j}^{\mathfrak{n}+1}\right)\vec{\boldsymbol{f}}^{\mathfrak{n}+1}_{i}=\vec{\boldsymbol{f}}^{\mathfrak{n}}_{i},$
(76)

where

$\vec{\boldsymbol{f}}^{\mathfrak{n}+1}_{i}=\bigg{(}[f_{i}^{\mathfrak{n}+1}]_{k=1},\cdots,[f_{i}^{\mathfrak{n}+1}]_{k=N_{v}}\bigg{)}^{\top},$
(77)

and

$G_{i,j}^{\mathfrak{n}+1}=\begin{pmatrix}b_{1}^{\mathfrak{n}+1}&-c_{2}^{\mathfrak{n}+1}&0&\cdots&&\cdots&0\\\
-a_{1}^{\mathfrak{n}+1}&b_{2}^{\mathfrak{n}+1}&-c_{3}^{\mathfrak{n}+1}&0&\cdots&\cdots&0\\\
0&-a_{2}^{\mathfrak{n}+1}&b_{3}^{\mathfrak{n}+1}&-c_{4}^{\mathfrak{n}+1}&0&\cdots&0\\\
\vdots&&\ddots&\ddots&\ddots&&\vdots\\\ \\\ \\\
0&\cdots&&0&-a_{N_{v}-2}^{\mathfrak{n}+1}&b_{N_{v}-1}^{\mathfrak{n}+1}&-c_{N_{v}}^{\mathfrak{n}+1}\\\
0&\cdots&&&0&-a_{N_{v}-1}^{\mathfrak{n}+1}&b_{N_{v}}^{\mathfrak{n}+1}\end{pmatrix}.$
(78)

The values of the Maxwellians $M_{i,j}^{\mathfrak{n}+1}$ and hence
coefficients of $G_{i,j}^{\mathfrak{n}+1}$ are explicitly computed with the
moment values determined by the GST algorithm in Section 6.1. After this, the
implicit update of the nonlinear collision operator in (66) requires only the
linear inversion of (76).

##  7 Numerical examples

To provide some numerical examples, we appeal to [3], which considers a multi-
species LFP collision operator. We use simple constants rather than physically
correct values for the physical parameters:

$n_{1}=n_{2}=1,\qquad
q_{i}=q_{j}=1,\qquad|\log\Lambda_{i,j}|=1,\qquad\epsilon_{0}=1.$ (79)

For both tests, the initial distributions are set to Maxwellian initial
conditions with

$\displaystyle f_{1}(0,v)$
$\displaystyle=M_{1\,,\,0.5\,,\,0.25}^{(1)}\qquad\mbox{and}\qquad
f_{2}(0,v)=M_{1\,,\,-0.25\,,\,0.125}^{(1)},$ (80)

where $M_{n,\boldsymbol{u},\theta}^{(d)}$ is the $d$-dimensional Maxwellian
defined in (3). To ensure that $\alpha_{i,j}$, $\beta_{i,j}$, and
$\gamma_{i,j}$ satisfy the desired bounds, we enforce (56), i.e.
$\kappa_{i,j}\geq\frac{m_{i}+m_{j}}{2m_{i}}\eqqcolon\mu_{i,j}$ for all $i,j$.

The computational domain for the velocity space is truncated to $[-4,4]$ (as
in [3]), and is discretized using a cell-centered grid over $N_{v}=80$ uniform
cells. The results below were generated using $\Delta t=0.2$; smaller time
steps gave qualitatively similar results.

For the first test, we set $m_{1}=m_{2}=1$. Thus $\mu_{i,j}\leq 1$, and the
simulation is run with $\kappa_{i,j}=2$ for all $i,j$. The results are given
in Figure 1. For the second test, we let $m_{1}=2$ and $m_{2}=1$, and set
$\kappa_{i,j}=2\max\\{\mu_{i,j}\\}=2\mu_{2,1}=3$. The results are given in
Figure 2. In both examples, the long time behavior of the moments is as
expected. That is, both species relax toward the steady-state bulk velocity
and temperature given by $\boldsymbol{u}^{\infty}$ and $T^{\infty}$ in (23).

(a) Species velocities relax to steady state, $\boldsymbol{u}^{\infty}$.

(b) Species temperatures relax to steady state, $T^{\infty}$.

Figure 1: Case: Equal Masses: The velocities and temperatures relax to steady
state values.

(a) Species velocities relax to steady state, $\boldsymbol{u}^{\infty}$.

(b) Species temperatures relax to steady state, $T^{\infty}$.

Figure 2: Case: Mass Ratio$=2$. The velocities and temperatures relax to
steady state values.

##  8 Conclusions

In this paper, we have derived a multi-species Lenard-Bernstein (M-LB) model
that satisfies (i) conservation of species mass and global conservation of
momentum and energy; (ii) dissipates entropy and satisfies an H-Theorem; and
(iii) matches the pairwise momentum and temperature relaxation rates of the
Boltzmann collision operator with a Coulomb potential. We have also
demonstrated a simple numerical algorithm for implementing the (M-LB) model
with implicit time stepping. In future work, we hope to embed this algorithm
into a full kinetic simulation capability, such as [6] or [19].

##  Appendix A Comparison of models

During the finalization of this manuscript, the authors became aware of very
recent and independent work on the M-LB model [18]. This work does not
consider the problem of matching relaxation rates or perform any numerical
experiments, but it does establish conditions on the coefficients of the
mixture temperature and velocity that are sufficient for the existence of an
$\mathcal{H}$-Theorem. Below we show that the conditions made in [18] are
sufficient to imply the conditions made in Assumption 4.1, which we believe
provide a more transparent physical interpretation, and allow us to establish
the entropy dissipation law and $\mathcal{H}$-Theorem under slightly weaker
conditions. We provide more details below, using the notation of [18] when
possible, but relying on the notation used here when there is a conflict. For
clarity we consider, as in [18], that $i\in\\{1,2\\}$, although this is not a
necessary restriction for either work.

Let

$c_{1,2}=\frac{\lambda_{1,2}}{n_{2}},\qquad
c_{2,1}=\frac{\lambda_{2,1}}{n_{1}},\qquad\mbox{and}\qquad\varepsilon=\frac{c_{1,2}}{c_{2,1}}=\frac{n_{1}}{n_{2}}\frac{\lambda_{1,2}}{\lambda_{1,2}}\geq
0.$ (81)

In terms of $\varepsilon$, the conditions in (12) become

$\displaystyle\alpha_{2,1}$
$\displaystyle=1-\frac{\rho_{1}\lambda_{1,2}}{\rho_{2}\lambda_{2,1}}(1-\alpha_{2,1}),$
(82a) $\displaystyle\beta_{2,1}$
$\displaystyle=1-\frac{n_{1}\lambda_{1,2}}{n_{2}\lambda_{2,1}}(1-\beta_{1,2}),$
(82b) $\displaystyle\gamma_{2,1}$
$\displaystyle=m_{1}\varepsilon(1-\alpha_{1,2})-\varepsilon\gamma_{1,2}.$
(82c)

######  Assumption A.1.

For $i=1,2$, the coefficients $\alpha_{i,j}$, $\beta_{i,j}$, and
$\gamma_{i,j}$ satisfy the following conditions:

$\displaystyle\varepsilon\leq 1,$ (83a)
$\displaystyle\varepsilon\leq\frac{m_{2}}{m_{1}},$ (83b) $\displaystyle
0\leq\gamma_{1,2}\leq m_{1}(1-\alpha_{1,2}),$ (83c)
$\displaystyle\frac{\varepsilon}{1+\varepsilon}\leq\alpha_{1,2}\leq 1,$ (83d)
$\displaystyle\frac{\varepsilon}{1+\varepsilon}\leq\beta_{1,2}\leq 1.$ (83e)

######  Remark A.2.

We list the conditions above in the order they appear in [18]. The conditions
in (83a) and (83b) are given in [18, Eq. 13]. The condition in (83c) is given
in [18, Eq. 20]. The conditions in (83d) and (83e) are given in [18, Eq. 23].
A third condition for entropy dissipation and an H-Theorem is given in [18,
Eq. 23], but we do not need it here.

######  Theorem A.3.

Suppose the pairwise conservation laws in Assumption 3.2 hold. Then for
$i=1,2$, the conditions in Assumption A.1 imply the conditions in Assumption
4.1.

###### Proof.

The conditions in (83d) and (83e) clearly imply that
$0\leq\alpha_{1,2},\beta_{1,2}\leq 1$. The condition in (83a) is equivalent to

$\frac{n_{1}}{n_{2}}\frac{\lambda_{1,2}}{\lambda_{2,1}}\leq 1.$ (84)

According to Proposition 3.4, the pairwise conservation of momentum and energy
in Assumption 3.2 are equivalent to (82). Thus according to (82b), (84), and
the fact that $\beta_{1,2}\in[0,1]$, it holds that $\beta_{2,1}\in[0,1]$.
Similarly, the condition in (83a) is equivalent to

$\frac{\rho_{1}}{\rho_{2}}\frac{\lambda_{1,2}}{\lambda_{2,1}}\leq 1.$ (85)

Acccording to (82a), (85), and the fact that $\alpha_{1,2}\in[0,1]$, it holds
that $\alpha_{2,1}\in[0,1]$. Finally, (83c) states directly that
$\gamma_{1,2}\geq 0$. Moreoever, when combined with (82c), the upper bound in
(83c) implies that $\gamma_{2,1}\geq 0$. ∎

##  Appendix B Model relaxation rates

In this appendix, we prove (45), in Lemma B.1 and Lemma B.2.

######  Lemma B.1.

Momentum relaxation rates between species $i$ and $j$ are given by

$\left.\frac{\partial}{\partial
t}\left(\rho_{i}\boldsymbol{u}_{i}-\rho_{j}\boldsymbol{u}_{j}\right)\right|_{\textnormal{M-LB}}=2\rho_{i}\lambda_{i,j}(1-\alpha_{i,j})(\boldsymbol{u}_{j}-\boldsymbol{u}_{i}).$
(86)

###### Proof.

Integrate (6) against $m_{i}\boldsymbol{v}$ and use (11a) to obtain the
following system of equations that describes the time evolution of the
momentum moment system:

$\displaystyle{\int_{\mathbb{R}^{d}}}\frac{\partial f_{i}}{\partial
t}m_{i}\boldsymbol{v}\,\,\textnormal{d}\boldsymbol{v}$
$\displaystyle=m_{i}\sum_{j}{\int_{\mathbb{R}^{d}}}\mathcal{C}_{i,j}(f_{i},f_{j})\boldsymbol{v}\,\,\textnormal{d}\boldsymbol{v}=\sum_{j}m_{i}\lambda_{i,j}\theta_{i,j}{\int_{\mathbb{R}^{d}}}\boldsymbol{v}\mathcal{L}(f_{i},M_{i,j})\,\,\textnormal{d}\boldsymbol{v}$
$\displaystyle\iff\frac{\partial}{\partial t}(\rho_{i}\boldsymbol{u}_{i})$
$\displaystyle=\sum_{j}\rho_{i}\lambda_{i,j}(\boldsymbol{u}_{i,j}-\boldsymbol{u}_{i}).$
(87)

To derive the relaxation rates between species $i$ and $j$, neglect the
species in $\\{1,\cdots,N\\}\setminus\\{i,j\\}$, take the difference of
relaxation rates in (87), and use (9a) and (15a) to obtain (86). ∎

######  Lemma B.2.

Thermal relaxation rates between species $i$ and $j$ are given by

$\left.\frac{\partial}{\partial
t}\frac{d}{2}\left(n_{i}T_{i}-n_{j}T_{j}\right)\right|_{\textnormal{M-LB}}=2dn_{i}\lambda_{i,j}(1-\beta_{i,j})(T_{j}-T_{i})+(2n_{i}\lambda_{i,j}\gamma_{i,j}-\delta_{i,j})|\boldsymbol{u}_{i}-\boldsymbol{u}_{j}|^{2}.$
(88)

###### Proof.

Integrating (6) against $m_{i}|\boldsymbol{v}|^{2}/2$, and recalling the
definition of the energy $E_{i}$, in (2), gives:

$\displaystyle{\int_{\mathbb{R}^{d}}}\frac{\partial f_{i}}{\partial
t}m_{i}\frac{|\boldsymbol{v}|^{2}}{2}\,\,\textnormal{d}\boldsymbol{v}$
$\displaystyle=\sum_{j}m_{i}\lambda_{i,j}\theta_{i,j}{\int_{\mathbb{R}^{d}}}\frac{|\boldsymbol{v}|^{2}}{2}\mathcal{L}(f_{i},M_{i,j})\,\,\textnormal{d}\boldsymbol{v}$
(89) $\displaystyle\iff\frac{\partial}{\partial
t}\left(\frac{1}{2}\rho_{i}|\boldsymbol{u}_{i}|^{2}+\frac{d}{2}n_{i}T_{i}\right)$
$\displaystyle=\sum_{j}\rho_{i}\lambda_{i,j}\boldsymbol{u}_{i}\cdot(\boldsymbol{u}_{i,j}-\boldsymbol{u}_{i})+d\sum_{j}n_{i}\lambda_{i,j}(T_{i,j}-T_{i}).$

Using (87), the first term is

$\frac{\partial}{\partial
t}\left(\frac{1}{2}\rho_{i}|\boldsymbol{u}_{i}|^{2}\right)=\boldsymbol{u}_{i}\cdot\frac{\partial}{\partial
t}(\rho_{i}\boldsymbol{u}_{i})=\sum_{j}\rho_{i}\lambda_{i,j}\boldsymbol{u}_{i}\cdot(\boldsymbol{u}_{i,j}-\boldsymbol{u}_{i}).$
(90)

With this, it is clear that the terms in (89) involving the velocity cancel.
Thus, (89) becomes

$\frac{\partial}{\partial
t}\left(\frac{d}{2}n_{i}T_{i}\right)=d\sum_{j}n_{i}\lambda_{i,j}(T_{i,j}-T_{i}),$
(91)

a system of equations that describes the time evolution of the temperature
moment system. As above, we neglect species in
$\\{1,\cdots,N\\}\setminus\\{i,j\\}$, take differences of relaxation rates in
(91), and use (9b) and (15c) to obtain (88). ∎

######  Remark B.3.

In this time-continuous case, the relaxation moment equation obtained for the
energy relaxation is simpler than for the BGK model of [9, 10], since, as
shown in (90), the terms involving the bulk velocity cancel directly, whereas,
these terms must remain as part of the energy equation in [10].

## References

  * [1] P.L. Bhatnagar, E.P. Gross, and M. Krook.  A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems.  Phys. Rev., 94:511–525, 1954.
  * [2] Alexander V Bobylev, Marzia Bisi, Maria Groppi, Giampiero Spiga, and Irina F Potapenko.  A general consistent BGK model for gas mixtures.  Kinetic & Related Models, 11(6), 2018.
  * [3] José A. Carrillo, Jingwei Hu, and Samuel Q. Van Fleet.  A particle method for the multispecies Landau equation.  arXiv.org Preprint, 2023.
  * [4] Carlo Cercignani, Reinhard Illner, and Mario Pulvirenti.  The mathematical theory of dilute gases, volume 106.  Springer Science & Business Media, 2013.
  * [5] J.P. Dougherty.  Model Fokker-Planck equation for a plasma and its solution.  The Physics of Fluids, 7(11):1788–1799, 1964.
  * [6] Eirik Endeve and Cory D Hauck.  Conservative DG method for the micro-macro decomposition of the Vlasov–Poisson–Lenard–Bernstein model.  Journal of Computational Physics, 462:111227, 2022.
  * [7] Manaure Francisquez, James Juno, Ammar Hakim, Gregory W. Hammett, and Darin R. Ernst.  Improved multispecies Dougherty collisions.  Journal of Plasma Physics, 88(3):1–36, 2022.
  * [8] Irene M. Gamba and Jeffrey R. Haack.  A conservative spectral method for the Boltzmann equation with anisotropic scattering and the grazing collisions limit.  Journal of Computational Physics, 270:40–57, 2014.
  * [9] J.R. Haack, C.D. Hauck, and M.S. Murillo.  A conservative, entropic multispecies BGK model.  J. Stat Phys, 168:826–856, 2017.
  * [10] Evan Habbershaw, Ryan Glasby, Jeffrey Haack, Cory Hauck, and Steven Wise.  Asymptotic relaxation of moment equations for a multi-species, homogeneous BGK model.  SIAM J. Appl. Math. (in review), 2023.  Preprint at https://arxiv.org/abs/2310.12885.
  * [11] Evan Habbershaw, Cory Hauck, and Steven Wise.  Implicit update of the moment equations for a multi-species, homogeneous BGK model.  arXiv.org Preprint, 2024.
  * [12] Jingwei Hu, Ansgar Jüngel, and Nicola Zamponi.  Global weak solutions for a nonlocal multispecies Fokker-Planck-Landau system.  Kinetic and Related Models, online first, 2024.
  * [13] Christian Klingenberg, Marlies Pirner, and Gabriella Puppo.  A consistent kinetic model for a two-component mixture with an application to plasma.  Kinetic & Related Models, 10(2), 2017.
  * [14] Lev Davidovich Landau.  The kinetic equation in the case of Coulomb interaction.  Technical report, General Dynamics/Astronautics San Diego Calif, 1958.
  * [15] Andrew Lenard and Ira B. Bernstein.  Plasma oscillations with diffusion in velocity space.  Physical Review, 112(5):1456–1459, 1958.
  * [16] Richard L. Liboff.  Kinetic theory: classical, quantum, and relativistic descriptions.  Springer Science & Business Media, 2003.
  * [17] T.F. Morse.  Energy and momentum exchange between nonequipartition gases.  Phys. Fluids, 6(10):1420–1427, 1963.
  * [18] Marlies Pirner.  A consistent kinetic Fokker-Planck model for gas mixtures.  arXiv preprint arXiv:2403.18478, 2024.
  * [19] Stefan Schnake, Coleman Kendrick, Eirik Endeve, Miroslav Stoyanov, Steven Hahn, Cory D Hauck, David L Green, Phil Snyder, and John Canik.  Sparse-grid discontinuous Galerkin methods for the Vlasov-Poisson-Lenard-Bernstein model.  arXiv preprint arXiv:2402.06493, 2024.
  * [20] Philipp Ulbl, Dominik Michels, and Frank Jenko.  Implementation and verification of a conservative, multi-species, gyro-averaged, full-f, Lenard-Bernstein/Dougherty collision operator in the gyrokinetic code GENE-X.  Contributions to Plasma Physics, 62(5-6):e202100180, 2022.