# Implicit Update of the Moment Equations for a Multi-Species, Homogeneous BGK
Model ††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).

Evan Habbershaw, Cory D. Hauck, and Steven M. Wise  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 Mathematics, The University of Tennessee, Knoxville, TN, 37996

(27 April 2024)

###### Abstract

A simple iterative approach for solving a set of implicit kinetic moment
equations is proposed. This implicit solve is a key component in the IMEX
discretization of the multi-species Bhatnagar-Gross-Krook (M-BGK) model with
nontrivial collision frequencies depending on individual species temperatures.
We prove that under mild time step restrictions, the iterative method
generates a contraction mapping. Numerical simulations are provided to
illustrate results of the IMEX scheme using the implicit moment solver.

Keywords: Multi-species BGK, fixed point iterative solvers, implicit time
stepping strategies, moment equations, kinetic theory, rarefied gas dynamics

MSC Codes: 76P05, 47J25, 35Q70, 82C40, 65L04

##  1 Introduction

Bhatnagar-Gross-Krook (BGK) models (see e.g. [22]) are a computationally
inexpensive alternative to the Boltzmann equation of rarefied gas dynamics.
They are derived by replacing the sophisticated Boltzmann collision operator
with a nonlinear relaxation model. In the single species case, the target of
this relaxation model is a Maxwellian distribution that shares the same mass,
momentum, and energy density as the kinetic distribution. This model was first
introduced in [4], and rigorous mathematical results about the model can be
found in [25, 24].

Like the Boltzmann equation, the single-species BGK model recovers the
compressible system of Euler equations in the limit of infinite collisions
[9]. For simulation purposes, implicit-explicit (IMEX) methods are often
employed that treat the advection term explicitly and the stiff collision term
implicitly. The strategy for the implicit evaluation of the BGK operator is to
recognize that the moments that characterize the Maxwellian target are
conserved by the space-homogeneous problem. As a result, the implicit part of
the relaxation operator is reduced to a simple multiple of the identity. This
trick, introduced in [12], is the basis of many semi-implicit schemes.

More recently, multi-species versions of the BGK equation have been derived
that satisfy the conservation and entropy dissipation properties of the multi-
species Boltzmann equation. Some models still use a single relaxation operator
[1, 8] while others use a sum of operators with Maxwellian-like targets for
each binary interaction [6, 20, 16, 15]. We focus here on the model from [16],
but the approach and analysis can be extended to the generalizations in [6,
20]. In these models, each relaxation operator is equipped with a collision
frequency. Conservation of the number density in the multi-species setting
implies that frequencies that depend only on the number densities will not be
affected by the collision process. However, for general molecular models,
collision frequencies may depend on the species temperature, average
velocities, and even microscopic velocities [5], and removing such
dependencies can significantly alter the behavior of solutions [15]. In
addition, the analysis and simulation of the BGK equations becomes more
difficult when the frequencies depend on time and/or phase-space velocity.

The focus of the present paper is the implicit update of the collision term in
the multi-species setting when more realistic collision frequencies are
employed. In order to follow the strategy used for the single-species model,
one must update an ODE system for the species’ average velocities and
temperatures, thereby enabling the evaluation of the Maxwellian target at the
next time step. As in the single species case, this approach effectively
linearizes the problem, allowing for a simple update of the collision term.
However, unlike the single species case, the moments of each species are not
conserved; rather the individual species densities, total momentum, and total
energy are conserved. As a result, the ODE for the average velocities and
temperatures requires an implicit solve. For collision frequencies that depend
only on species densities, such a solve is still trivial to implement because
the frequencies are not affected by the collisions. However, for more
complicated collision frequencies, a fully implicit solve must take such
variations into account.

We propose here to solve the multi-species BGK collision operator in a fully
implicit manner using a nonlinear Gauss-Seidel-type iterative solver. The
focus of the paper is on the numerical analysis of the solver. In particular,
we demonstrate a contraction mapping under mild time step restrictions.
Importantly, we show that, given a user specified time step, there is a stable
time step that, at worst, is smaller by an $\mathcal{O}(1)$ constant that does
not depend on the stiffness of the collision operator. Thus, in the fluid
regime, the method can still take time steps determined by particle advection
and not particle collisions. This result and the underlying analysis can be
extended to general relaxation systems with state-dependent relaxation times.

The remainder of the paper is organized as follows. In Section 2, we introduce
the model and discuss the IMEX strategy that provides the context for the
implicit solver. In Section 3, we introduce the solver. In Section 4 and
Section 5, we analyze the convergence of the solver, and in Section 6, we
analyze the associated time step restriction. In Section 7, we provide several
numerical demonstrations, and in Section 8 we make conclusions and discuss
future work. The Appendix contains several results used in the main body of
the paper.

##  2 Background

###  2.1 Models and equations

We denote by $f_{i}$ the kinetic distribution of particles of species
$i\in\\{1,\cdots,N\\}$ having mass $m_{i}>0$. More specifically,
$f_{i}(\boldsymbol{x},\boldsymbol{v},t)$ is the density of species $i$
particles at the point $\boldsymbol{x}\in\Omega\subset\mathbb{R}^{d}$, with
microscopic velocity $\boldsymbol{v}\in\mathbb{R}^{d}$, at time $t\geq 0,$
with respect to the measure
$\,\textnormal{d}\boldsymbol{v}\,\,\textnormal{d}\boldsymbol{x}$. Associated
to each $f_{i}$ are the species number density $n_{i}$, bulk velocity
$\boldsymbol{u}_{i}$, and temperature $T_{i}$, defined by

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

The species mass density $\rho_{i}$, momentum density $\boldsymbol{k}_{i}$,
and energy density $E_{i}$ are

$\rho_{i}=m_{i}n_{i},\qquad\boldsymbol{k}_{i}=\rho_{i}\boldsymbol{u}_{i},\qquad\mbox{and}\qquad
E_{i}=\frac{1}{2}\rho_{i}|\boldsymbol{u}_{i}|^{2}+\frac{d}{2}n_{i}T_{i}.$ (2)

BGK models are expressed in terms of Maxwellian distributions that depend on
these moments.

######  Definition 2.1.

Given $n>0$, $\boldsymbol{u}\in\mathbb{R}^{d}$, and $\theta>0$, the Maxwellian
$M_{n,\boldsymbol{u},\theta}$ is a function of the form

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

Straightforward computations show that the moments of
$M_{n,\boldsymbol{u},\theta}$ satisfy

${\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},\qquad{\int_{\mathbb{R}^{d}}}|\boldsymbol{v}|^{2}M_{n,\boldsymbol{u},\theta}\,\,\textnormal{d}\boldsymbol{v}=n|\boldsymbol{u}|^{2}+dn\theta.$
(4)

We consider in this paper the multi-species BGK equation from [16, 17], given
by

$\frac{\partial f_{i}}{\partial
t}+\boldsymbol{v}\cdot\nabla_{\boldsymbol{x}}f_{i}=\frac{1}{\varepsilon}\sum_{j=1}^{N}\lambda_{i,j}(M_{i,j}-f_{i}),\qquad\forall
i\in\\{1,\cdots,N\\},$ (5)

where $\varepsilon>0$ is a dimensionless parameter that comes from a rescaling
of space and time, $\lambda_{i,j}>0$ is the frequency of collisions between
species $i$ and $j$ (independent of $\boldsymbol{v}$), and

$M_{i,j}(\boldsymbol{v})\coloneqq
M_{n_{i},\boldsymbol{u}_{i,j},T_{i,j}/m_{i}}(\boldsymbol{v})$ (6)

is a Maxwellian defined by (3), using the mixture velocities and temperatures

$\displaystyle\boldsymbol{u}_{i,j}$
$\displaystyle=\frac{\rho_{i}\lambda_{i,j}\boldsymbol{u}_{i}+\rho_{j}\lambda_{j,i}\boldsymbol{u}_{j}}{\rho_{i}\lambda_{i,j}+\rho_{j}\lambda_{j,i}},$
(7a) $\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_{j,i}}\frac{|\boldsymbol{u}_{i}-\boldsymbol{u}_{j}|^{2}}{n_{i}\lambda_{i,j}+n_{j}\lambda_{j,i}}.$
(7b)

As in [17], we consider a hard spheres collision model for the collision
frequencies:

$\lambda_{i,j}\coloneqq\lambda_{i,j}^{\textnormal{HS}}(n_{j},T_{i},T_{j})=\frac{32\pi^{2}}{3(2\pi)^{3/2}}\frac{m_{i}m_{j}}{(m_{i}+m_{j})^{2}}(d_{i}+d_{j})^{2}n_{j}\sqrt{\frac{T_{i}}{m_{i}}+\frac{T_{j}}{m_{j}}},$
(8)

where $d_{i}$ and $d_{j}$ are reference diameters for the particles of species
$i$ and $j$. For simplicity, we write $\lambda_{i,j}=c_{i,j}n_{j}\Psi_{i,j},$
where

$\Psi_{i,j}(T_{i},T_{j})=\sqrt{\frac{T_{i}}{m_{i}}+\frac{T_{j}}{m_{j}}}$ (9)

is the only time-dependent term. The challenging features of
$\lambda_{i,j}^{\textnormal{HS}}$ are (i) it depends on species temperatures
that are not conserved by collision dynamics, and (ii) this temperature
dependence is not Lipschitz at zero temperatures.

Going forward, we will abbreviate summation symbols $\sum_{j=1}^{N}$ as
$\sum_{j}$, when the upper and lower index is understood as $j=1$ and $j=N$,
respectively.

###  2.2 IMEX strategy and backward Euler step

Implicit-explicit (IMEX) Runge Kutta methods are a common tool for simulating
kinetic equations. These methods typically treat advection explicitly and
collisions implicitly. A diagonally implicit IMEX Runge-Kutta method for (5)
with $\nu$ stages and time step $\Delta t^{\mathfrak{n}}>0$ takes the form
[26]

$\displaystyle f_{i}^{(s)}$ $\displaystyle=f_{i}^{(s,\star)}+\frac{\Delta
t^{\mathfrak{n}}}{\varepsilon}a_{s,s}\sum_{j}\lambda_{i,j}^{(s)}(M^{(s)}_{i,j}-f^{(s)}_{i}),\qquad
s\in\\{1,\cdots,\nu\\},$ (10a) $\displaystyle f^{\mathfrak{n}+1}_{i}$
$\displaystyle=f_{i}^{\mathfrak{n}}-\Delta
t^{\mathfrak{n}}\sum_{s=1}^{\nu}\widetilde{b}_{s}\boldsymbol{v}\cdot\nabla_{\boldsymbol{x}}f_{i}^{(s)}+\frac{\Delta
t^{\mathfrak{n}}}{\varepsilon}\sum_{s=1}^{\nu}b_{s}\sum_{j}\lambda_{i,j}^{(s)}(M_{i,j}^{(s)}-f_{i}^{(s)}),$
(10b)

where

$f_{i}^{(s,\star)}=f_{i}^{\mathfrak{n}}-\Delta
t^{\mathfrak{n}}\sum_{r=1}^{s-1}\widetilde{a}_{s,r}\boldsymbol{v}\cdot\nabla_{\boldsymbol{x}}f^{(r)}_{i}+\frac{\Delta
t^{\mathfrak{n}}}{\varepsilon}\sum_{r=1}^{s-1}a_{s,r}\sum_{j}\lambda_{i,j}^{(r)}(M^{(r)}_{i,j}-f^{(r)}_{i}),$
(11)

and
$\lambda_{i,j}^{(r)}=\lambda_{i,j}^{\textnormal{HS}}(n_{j}^{(r)},T^{(r)}_{i},T^{(r)}_{j})$.
The coefficients $a_{s,r}$ and $\widetilde{a}_{s,r}$ are determined by order
conditions and other desired properties of the method [23, 3, 7].

A prototype for the update in (10a), and the focus of this paper, is the
Backward Euler update, namely,

$f_{i}^{\mathfrak{n}+1}=f_{i}^{\mathfrak{n}}+\frac{\Delta
t}{\varepsilon}\sum_{j}\lambda_{i,j}^{\mathfrak{n}+1}(M^{\mathfrak{n}+1}_{i,j}-f^{\mathfrak{n}+1}_{i}).$
(12)

Indeed, given $f_{i}^{(s,\star)}$, solving (12) is equivalent to solving (10)
with slightly modified constants. The standard approach to solve the nonlinear
equation (12) is to first solve moment equations for $n^{\mathfrak{n}+1}_{i}$,
$\boldsymbol{u}^{\mathfrak{n}+1}_{i}$, and $T^{\mathfrak{n}+1}_{i}$, which can
then be used to define $M^{\mathfrak{n}+1}_{i,j}$ in an explicit fashion; what
remains is a linear solve for $f_{i}^{\mathfrak{n}+1}$ [26]. In the single
species case, the update is trivial: due to conservation properties of the
collision operator
$(n_{i}^{\mathfrak{n}+1},\boldsymbol{u}_{i}^{\mathfrak{n}+1},T_{i}^{\mathfrak{n}+1})=(n_{i}^{\mathfrak{n}},\boldsymbol{u}_{i}^{\mathfrak{n}},T_{i}^{\mathfrak{n}})$
[12]. In the multi-species case, taking moments of (12) and applying the
definitions of $\boldsymbol{u}_{i,j}$ and $T_{i,j}$ from (7) gives the
following algebraic system:

######  Proposition 2.2 (Backward Euler).

If scheme (12) is applied for the kinetic densities, then, for each
$i\in\\{1,\cdots,N\\}$,

$\displaystyle\frac{n_{i}^{\mathfrak{n}+1}-n_{i}^{\mathfrak{n}}}{\Delta t}$
$\displaystyle=0,$ (13a)
$\displaystyle\rho_{i}\frac{\boldsymbol{u}_{i}^{\mathfrak{n}+1}-\boldsymbol{u}_{i}^{\mathfrak{n}}}{\Delta
t}$
$\displaystyle=\frac{1}{\varepsilon}\sum_{j}\left(A_{i,j}^{\mathfrak{n}+1}-D_{i,j}^{\mathfrak{n}+1}\right)\boldsymbol{u}_{j}^{\mathfrak{n}+1},$
(13b)
$\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}&=\frac{d}{\varepsilon}\sum_{j}\left(B_{i,j}^{\mathfrak{n}+1}-F_{i,j}^{\mathfrak{n}+1}\right)T_{j}^{\mathfrak{n}+1}\\\
+\frac{1}{\varepsilon}\sum_{j}B_{i,j}^{\mathfrak{n}+1}\\!&\,\left[m_{j}(s_{j}^{\mathfrak{n}+1}-S_{i,j}^{\mathfrak{n}+1})-m_{i}(s_{i}^{\mathfrak{n}+1}-S_{i,j}^{\mathfrak{n}+1})\right],\end{split}$
(13c)

where $s_{i}^{\mathfrak{n}}=|\boldsymbol{u}_{i}^{\mathfrak{n}}|^{2}$ and the
${N\times N}$ matrices $A$, $B$, $D$, $F$, and $S$, at time step
$t^{\mathfrak{n}+1}$ are given by

$\displaystyle[A]_{i,j}^{\mathfrak{n}+1}=\frac{\rho_{i}\rho_{j}\lambda_{i,j}^{\mathfrak{n}+1}\lambda_{j,i}^{\mathfrak{n}+1}}{\rho_{i}\lambda_{i,j}^{\mathfrak{n}+1}+\rho_{j}\lambda_{j,i}^{\mathfrak{n}+1}},\qquad[B]_{i,j}^{\mathfrak{n}+1}=\frac{n_{i}n_{j}\lambda_{i,j}^{\mathfrak{n}+1}\lambda_{j,i}^{\mathfrak{n}+1}}{n_{i}\lambda_{i,j}^{\mathfrak{n}+1}+n_{j}\lambda_{j,i}^{\mathfrak{n}+1}},\qquad[S]_{i,j}^{\mathfrak{n}+1}=\left|\boldsymbol{u}_{i,j}^{\mathfrak{n}+1}\right|^{2},$
(14a)
$\displaystyle[D]_{i,j}^{\mathfrak{n}+1}=\Big{(}\sum_{k}A_{i,k}^{\mathfrak{n}+1}\Big{)}\delta_{i,j},\qquad[F]_{i,j}^{\mathfrak{n}+1}=\Big{(}\sum_{k}B_{i,k}^{\mathfrak{n}+1}\Big{)}\delta_{i,j}.$
(14b)

######  Remark 2.3.

According to (13a), the quantities $n_{i}$ and $\rho_{i}$ are independent of
$\mathfrak{n}$. We therefore do not include the superscript on these
variables.

We summarize some useful properties of the backward Euler method below.

######  Proposition 2.4 (Discrete Conservation Laws).

The backward Euler scheme in (13) satisfies discrete versions of the
conservation of total momentum and conservation of total energy. Specifically,

$\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}}.$
(15)

###### Proof.

Sum (13b) and (13c) over all species $i$, and use the symmetry of
$A_{i,j}^{\mathfrak{n}+1}$ and $B_{i,j}^{\mathfrak{n}+1}$ to obtain the
result. For example, summing (13b) on $i$, using the definition of
$D_{i,j}^{\mathfrak{n}+1}$, and symmetry of $A_{i,j}^{\mathfrak{n}+1}$ gives
$\frac{1}{\varepsilon}\sum_{i,j}(A_{i,j}^{\mathfrak{n}+1}-D_{i,j}^{\mathfrak{n}+1})\boldsymbol{u}_{j}^{\mathfrak{n}+1}=\frac{1}{\varepsilon}\sum_{i,j}A_{i,j}^{\mathfrak{n}+1}(\boldsymbol{u}_{j}^{\mathfrak{n}+1}-\boldsymbol{u}_{i}^{\mathfrak{n}+1})=0$.
∎

######  Proposition 2.5.

The quantities

$\boldsymbol{u}^{\infty}=\frac{\sum_{i}\rho_{i}\boldsymbol{u}_{i}^{\mathfrak{n}}}{\sum_{i}\rho_{i}}\qquad\mbox{and}\qquad
T^{\infty}=\frac{\sum_{i}n_{i}T_{i}^{\mathfrak{n}}}{\sum_{i}n_{i}}+\frac{\sum_{i}\rho_{i}\left(|\boldsymbol{u}_{i}^{\mathfrak{n}}|^{2}-|\boldsymbol{u}^{\infty}|^{2}\right)}{d\sum_{i}n_{i}}$
(16)

are independent of the time step index $\mathfrak{n}$.

###### Proof.

The fact that $\boldsymbol{u}^{\infty}$ is independent of $\mathfrak{n}$
follows immediately from the first equation in (15). Meanwhile,

$T^{\infty}=\frac{2\sum_{i}E_{i}^{\mathfrak{n}}-\sum_{i}\rho_{i}|\boldsymbol{u}^{\infty}|^{2}}{d\sum_{i}n_{i}}$
(17)

is independent of $\mathfrak{n}$ according to the second equation in (15) and
the fact that $\boldsymbol{u}^{\infty}$ is also independent of $\mathfrak{n}$.
∎

When collision frequencies are constant or depend only on the number
densities, $\lambda_{i,j}^{\mathfrak{n}+1}=\lambda_{i,j}^{\mathfrak{n}}$, and
solving (13) is straightforward [26]. First, (13b) requires a linear solve for
$\boldsymbol{u}^{\mathfrak{n}+1}_{i}$. The new velocity can then be used to
evaluate implicit source terms in (13c), leaving another linear solve to find
each $T_{i}$. However, for more realistic frequencies that depend on $T_{i}$,
this approach no longer works. Two options remain: (i) lag $\lambda_{i,j}$ at
the previous time step or (ii) solve (13) in a fully nonlinear fashion.

Since the first option may lead to loss of accuracy or instability, we focus
here on the second option. In terms of cost, the linear solve is trivial
compared to the rest of the kinetic simulation. The challenge is to establish
an iterative algorithm that is provably convergent, even when $\Delta
t/\varepsilon$ is large.

##  3 Iterative Solver

###  3.1 Definition and Main Result

We propose the following Gauss-Seidel type (GST) iterative solver to solve the
backward Euler system in (13). We show below that, under a mild time step
restriction which is independent of $\varepsilon$, the GST iteration has a
fixed point that yields the unique solution to (13).

######  Definition 3.1 (GST method).

Let $\boldsymbol{u}_{i}^{\mathfrak{n}+1,0}=\boldsymbol{u}_{i}^{\mathfrak{n}}$
and $T_{i}^{\mathfrak{n}+1,0}=T_{i}^{\mathfrak{n}}$ be given. For
$\ell\in\\{0,1,2,\cdots\\}$, let $\boldsymbol{u}_{i}^{\mathfrak{n}+1,\ell+1}$
and $T_{i}^{\mathfrak{n}+1,\ell+1}$ be given by

$\displaystyle\rho_{i}\frac{\boldsymbol{u}_{i}^{\mathfrak{n}+1,\ell+1}-\boldsymbol{u}_{i}^{\mathfrak{n}}}{\Delta
t}$
$\displaystyle=\frac{1}{\varepsilon}\sum_{j}\left(A_{i,j}^{\mathfrak{n}+1,\ell}-D_{i,j}^{\mathfrak{n}+1,\ell}\right)\boldsymbol{u}_{j}^{\mathfrak{n}+1,\ell+1},$
(18a)
$\displaystyle\begin{split}\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}&=\frac{d}{\varepsilon}\sum_{j}\left(B_{i,j}^{\mathfrak{n}+1,\ell}-F_{i,j}^{\mathfrak{n}+1,\ell}\right)T_{j}^{\mathfrak{n}+1,\ell+1}\\\
&+\frac{1}{\varepsilon}\sum_{j}B_{i,j}^{\mathfrak{n}+1,\ell}\left[m_{j}\left(s_{j}^{\mathfrak{n}+1,\ell+1}-S_{i,j}^{\mathfrak{n}+1,\ell+1}\right)-m_{i}\left(s_{i}^{\mathfrak{n}+1,\ell+1}-S_{i,j}^{\mathfrak{n}+1,\ell+1}\right)\right].\end{split}$
(18b)

To simplify notation, we drop the superscript $\mathfrak{n}+1$, since it is
clear that we are seeking the update at this time level. Thus, instead of
$\boldsymbol{u}_{i}^{\mathfrak{n}+1,\ell}$ we use $\boldsymbol{u}_{i}^{\ell}$,
and, likewise, $T_{i}^{\mathfrak{n}+1,\ell}$ is replaced by $T_{i}^{\ell}$. In
particular, the initialization
$\boldsymbol{u}_{i}^{\mathfrak{n}+1,0}\coloneqq\boldsymbol{u}_{i}^{\mathfrak{n}}$
is replaced by
$\boldsymbol{u}_{i}^{\ell=0}\coloneqq\boldsymbol{u}_{i}^{\mathfrak{n}}$, and
$T_{i}^{\mathfrak{n}+1,0}\coloneqq T_{i}^{\mathfrak{n}}$ is replaced by
$T_{i}^{\ell=0}\coloneqq T_{i}^{\mathfrak{n}}$. Furthermore, to facilitate the
analysis of the GST method, we introduce the auxiliary variables

$\boldsymbol{W}=[\boldsymbol{w}_{1},\cdots,\boldsymbol{w}_{N}]^{\top}=P^{\nicefrac{{1}}{{2}}}\boldsymbol{U}\qquad\mbox{and}\qquad\boldsymbol{\eta}=Q^{\nicefrac{{1}}{{2}}}\boldsymbol{T},$
(19)

where $P=\text{diag}\\{\rho_{k}\\}\in\mathbb{R}^{N\times N}$,
$Q=\text{diag}\\{n_{k}\\}\in\mathbb{R}^{N\times N}$,
$\boldsymbol{U}=[\boldsymbol{u}_{1},\cdots,\boldsymbol{u}_{N}]^{\top}\in\mathbb{R}^{N\times
d}$, and $\boldsymbol{T}=(T_{1},\cdots,T_{N})^{\top}\in\mathbb{R}^{N}$.

######  Proposition 3.2 (GST Method for $\boldsymbol{W}$ and
$\boldsymbol{\eta}$).

The iterates
$\boldsymbol{W}^{\ell}=P^{\nicefrac{{1}}{{2}}}\boldsymbol{U}^{\ell}$ and
$\boldsymbol{\eta}^{\ell}=Q^{\nicefrac{{1}}{{2}}}\boldsymbol{T}^{\ell}$
satisfy

$\displaystyle\left(I+\frac{\Delta
t}{\varepsilon}Z^{\ell}\right)\boldsymbol{W}^{\ell+1}$
$\displaystyle=\boldsymbol{W}^{\mathfrak{n}},$ (20a)
$\displaystyle\left(I+\frac{\Delta
t}{\varepsilon}\widehat{Z}^{\ell}\right)\boldsymbol{\eta}^{\ell+1}$
$\displaystyle=\boldsymbol{\eta}^{\mathfrak{n}}+\frac{1}{d}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}^{(1),\ell+1}+\frac{\Delta
t}{\varepsilon}\frac{1}{d}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}^{(2),\ell+1,\ell},$
(20b)

where the $N\times N$ matrices

$Z^{\ell}=P^{-\nicefrac{{1}}{{2}}}(D^{\ell}-A^{\ell})P^{-\nicefrac{{1}}{{2}}}\qquad\mbox{and}\qquad\widehat{Z}^{\ell}=Q^{-\nicefrac{{1}}{{2}}}(F^{\ell}-B^{\ell})Q^{-\nicefrac{{1}}{{2}}}$
(21)

are symmetric, and the source terms in (20b) are given by

$\displaystyle\boldsymbol{s}^{(1),\ell+1}$
$\displaystyle=\boldsymbol{s}_{\boldsymbol{W}}^{\mathfrak{n}}-\boldsymbol{s}_{\boldsymbol{W}}^{\ell+1},\qquad\text{with
}\qquad\boldsymbol{s}_{\boldsymbol{W}}=\left(|\boldsymbol{w}_{1}|^{2},\cdots,|\boldsymbol{w}_{N}|^{2}\right)^{\top},$
(22a) $\displaystyle\text{and
}\qquad\left[\boldsymbol{s}^{(2),\ell+1,\ell}\right]_{i}$
$\displaystyle=\sum_{j}B_{i,j}^{\ell}\left[m_{j}(s_{j}^{\ell+1}-S_{i,j}^{\ell+1})-m_{i}(s_{i}^{\ell+1}-S_{i,j}^{\ell+1})\right].$
(22b)

######  Definition 3.3.

Let

$c_{i,j}^{A}=\frac{16}{3}\sqrt{\frac{\pi}{2}}\frac{m_{i}m_{j}(d_{i}+d_{j})^{2}}{(m_{i}+m_{j})^{3}}\rho_{i}\rho_{j}\qquad\mbox{and}\qquad
c_{i,j}^{B}=\frac{8}{3}\sqrt{\frac{\pi}{2}}\frac{(d_{i}+d_{j})^{2}}{(m_{i}+m_{j})^{2}}\rho_{i}\rho_{j},$
(23)

so that $A^{\mathfrak{n}}_{i,j}=c_{i,j}^{A}\Psi_{i,j}^{\mathfrak{n}}$ and
$B^{\mathfrak{n}}_{i,j}=c_{i,j}^{B}\Psi_{i,j}^{\mathfrak{n}}$, with $\Psi$ as
in (9). Also define the bounds

$\displaystyle A_{\min}^{\mathfrak{n}}$
$\displaystyle\coloneqq\min_{i,j}A_{i,j}^{\mathfrak{n}}\leq
A_{i,j}^{\mathfrak{n}}\leq\max_{i,j}A_{i,j}^{\mathfrak{n}}\eqqcolon
A_{\max}^{\mathfrak{n}},$ (24a) $\displaystyle B_{\min}^{\mathfrak{n}}$
$\displaystyle\coloneqq\min_{i,j}B_{i,j}^{\mathfrak{n}}\leq
B_{i,j}^{\mathfrak{n}}\leq\max_{i,j}B_{i,j}^{\mathfrak{n}}\eqqcolon
B_{\max}^{\mathfrak{n}}.$ (24b)

Finally, define the following maximum and minimum values:

$\displaystyle m_{\min}\coloneqq\min\\{m_{k}\\},\qquad
n_{\min}\coloneqq\min\\{n_{k}\\},\qquad\rho_{\min}\coloneqq\min\\{\rho_{k}\\},$
(25a) $\displaystyle m_{\max}\coloneqq\max\\{m_{k}\\},\qquad
n_{\max}\coloneqq\max\\{n_{k}\\},\qquad\rho_{\max}\coloneqq\max\\{\rho_{k}\\}.$
(25b)

######  Proposition 3.4.

The eigenvalues of the matrices
$Z^{\ell}=P^{-\nicefrac{{1}}{{2}}}(D^{\ell}-A^{\ell})P^{-\nicefrac{{1}}{{2}}}$
and
$\widehat{Z}^{\ell}=Q^{-\nicefrac{{1}}{{2}}}(F^{\ell}-B^{\ell})Q^{-\nicefrac{{1}}{{2}}}$
satisfy

$\displaystyle 0$ $\displaystyle=z_{0}<z_{\min}^{\mathfrak{n}}\leq
z_{1}^{\ell}\leq\cdots\leq z_{N-1}^{\ell}\leq z_{\max}^{\mathfrak{n}},$ (26a)
$\displaystyle 0$
$\displaystyle=\widehat{z}_{0}<\widehat{z}_{\min}^{\mathfrak{n}}\leq\widehat{z}_{1}^{\ell}\leq\cdots\leq\widehat{z}_{N-1}^{\ell}\leq\widehat{z}_{\max}^{\mathfrak{n}},$
(26b)

where

$z_{\min}^{\mathfrak{n}}=\frac{A_{\min}^{\mathfrak{n}}N}{\rho_{\max}},\qquad
z_{\max}^{\mathfrak{n}}=\frac{A_{\max}^{\mathfrak{n}}(N-1)}{\rho_{\min}},\qquad\widehat{z}_{\min}^{\mathfrak{n}}=\frac{B_{\min}^{\mathfrak{n}}N}{n_{\max}},\qquad\widehat{z}_{\max}^{\mathfrak{n}}=\frac{B_{\max}^{\mathfrak{n}}(N-1)}{{n_{\min}}}.$
(27)

###### Proof.

The proof of this result follows as in [17, Appendix B], and uses (33). ∎

######  Definition 3.5.

For any quantity $K$ that depends on $\ell$, define the difference operator
$\Delta^{\ell}$ by $\Delta^{\ell}K\equiv\Delta^{\ell}(K)=K^{\ell+1}-K^{\ell}$.

The main result of this paper is the following theorem.

######  Theorem 3.6.

The differences $\Delta^{\ell}\boldsymbol{W}$ and
$\Delta^{\ell}\boldsymbol{\eta}$, defined according to (20), satisfy

$\displaystyle\|\Delta^{\ell}\boldsymbol{W}\|_{\textnormal{F}}$
$\displaystyle\leq
C_{\boldsymbol{W}}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2},$
(28a) $\displaystyle\|\Delta^{\ell}\boldsymbol{\eta}\|_{2}$ $\displaystyle\leq
C_{0}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2}+C_{1}^{\mathfrak{n}}\Gamma_{\boldsymbol{X}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2},$
(28b)

where the positive constants $C_{\boldsymbol{W}}^{\mathfrak{n}}$,
$C_{0}^{\mathfrak{n}}$, and $C_{1}^{\mathfrak{n}}$ depend only on
$\left\\{N,m_{k},d_{k},n_{k},\boldsymbol{W}^{\mathfrak{n}},\boldsymbol{\eta}^{\mathfrak{n}}\right\\}$;

$\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}=\frac{\Delta
t\varepsilon}{(\varepsilon+\Delta
tz_{\min}^{\mathfrak{n}})^{2}},\qquad\Gamma_{\boldsymbol{\eta}}^{\mathfrak{n}}=\frac{\Delta
t\varepsilon}{(\varepsilon+\Delta
t\widehat{z}_{\min}^{\mathfrak{n}})^{2}},\qquad\mbox{and}\qquad\Gamma_{\boldsymbol{X}}^{\mathfrak{n}}=\max\\{\Gamma_{\boldsymbol{W}}^{\mathfrak{n}},\Gamma_{\boldsymbol{\eta}}^{\mathfrak{n}}\\};$
(29)

and $z_{\min}^{\mathfrak{n}}$ and $\widehat{z}_{\min}^{\mathfrak{n}}$ are
defined in Proposition 3.4.

As a corollary, we obtain the following sufficient condition on the time step
to ensure the method is convergent. Importantly, this condition does not
depend on $\varepsilon$.

######  Corollary 3.6.1.

Given $\Delta t_{0}>0$, there exists $\Delta t_{1}$ such that (i) $C\Delta
t_{0}\leq\Delta t_{1}\leq\Delta t_{0}$ for a constant $C\in(0,1]$ that is
independent of $\varepsilon$, and (ii) for this $\Delta t_{1}$,

$C_{\boldsymbol{W}}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}+C_{0}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}+C_{1}^{\mathfrak{n}}\Gamma_{\boldsymbol{X}}^{\mathfrak{n}}<1,$
(30)

thus ensuring convergence of the GST method.

Most of the remainder of the paper is dedicated to the proof of Theorem 3.6.
In the next subsection, we collect some useful properties of the method.
Section 4 is dedicated to proving the bound in (28a), and Section 5 is
dedicated to proving the bound in (28b); see Theorems 4.9 and 5.15,
respectively. The proof of Corollary 3.6.1 is contained in Section 6.

###  3.2 Basic properties of the GST method

######  Proposition 3.7.

The iterates of the GST method satisfy the discrete version of the
conservation of total momentum and total energy:

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

###### Proof.

The proof follows as in Proposition 2.4, with $\mathfrak{n}+1$ replaced by
$\ell$. ∎

######  Proposition 3.8.

The velocity and temperature iterates satisfy the bounds

$\displaystyle\min_{j}\\{u_{j,k}^{\mathfrak{n}}\\}\leq\min_{j}\\{u_{j,k}^{\ell+1}\\}\leq\max_{j}\\{u_{j,k}^{\ell+1}\\}\leq\max_{j}\\{u_{j,k}^{\mathfrak{n}}\\},\qquad
k\in\\{1,\cdots,d\\},$ (32) $\displaystyle
T_{i}^{\ell+1}\geq\min_{j}\\{T_{j}^{\mathfrak{n}}\\}\eqqcolon
T_{\min}^{\mathfrak{n}},\qquad i\in\\{1,\cdots,N\\},$ (33)
$\displaystyle\|\boldsymbol{u}_{i}^{\ell+1}\|_{2}\leq\|\boldsymbol{u}_{\max}^{\mathfrak{n}}\|_{2}\eqqcolon
u_{\max}^{\mathfrak{n}},\qquad i\in\\{1,\cdots,N\\},$ (34)

where
$\left[\boldsymbol{u}_{\max}^{\mathfrak{n}}\right]_{k}=\max_{j}\\{|u_{j,k}^{\mathfrak{n}}|\\}$,
for $k\in\\{1,\cdots,d\\}$.

###### Proof.

We first verify the lower bound in (32). The velocity update (18a) is
equivalent to

$u_{i,k}^{\ell+1}=u_{i,k}^{\mathfrak{n}}+\frac{\Delta
t}{\varepsilon}\sum_{j}\lambda_{i,j}^{\ell}\alpha_{j,i}\left(u_{j,k}^{\ell+1}-u_{i,k}^{\ell+1}\right),$
(35)

where
$\alpha_{j,i}=\frac{\rho_{j}\lambda_{j,i}^{\ell}}{\rho_{i}\lambda_{i,j}^{\ell}+\rho_{j}\lambda_{j,i}^{\ell}}$
is independent of $\ell$ (see Lemma C.1). Write
$\underline{u}_{k}^{\ell}=\min_{j}\\{u_{j,k}^{\ell}\\}$. Then, with
$c_{i}^{\ell}=\sum_{j}\lambda_{i,j}^{\ell}\alpha_{j,i}$,

$\displaystyle u_{i,k}^{\ell+1}$
$\displaystyle\geq\underline{u}_{k}^{\mathfrak{n}}+\frac{\Delta
t}{\varepsilon}\sum_{j}\lambda_{i,j}^{\ell}\alpha_{j,i}\left(\underline{u}_{k}^{\ell+1}\\!-u_{i,k}^{\ell+1}\right)=\underline{u}_{k}^{\mathfrak{n}}+\frac{\Delta
t}{\varepsilon}c_{i}^{\ell}\underline{u}_{k}^{\ell+1}-\frac{\Delta
t}{\varepsilon}c_{i}^{\ell}u_{i,k}^{\ell+1}$ (36a)
$\displaystyle\Longrightarrow\quad u_{i,k}^{\ell+1}$
$\displaystyle\geq\frac{\varepsilon}{\varepsilon+\Delta
tc_{i}^{\ell}}\underline{u}_{k}^{\mathfrak{n}}+\frac{\Delta
tc_{i}^{\ell}}{\varepsilon+\Delta tc_{i}^{\ell}}\underline{u}_{k}^{\ell+1}.$
(36b)

Choose $i_{\star}\in\\{1,\cdots,N\\}$ such that
$u_{i_{\star},k}^{\ell+1}=\underline{u}_{k}^{\ell+1}$, and let
$c_{\star}^{\ell}=c_{i_{\star}}^{\ell}$. Then after some simple algebra, (36b)
implies that

$\underline{u}_{k}^{\ell+1}\geq\frac{\varepsilon}{\varepsilon+\Delta
tc_{\star}^{\ell}}\underline{u}_{k}^{\mathfrak{n}}+\frac{\Delta
tc_{\star}^{\ell}}{\varepsilon+\Delta
tc_{\star}^{\ell}}\underline{u}_{k}^{\ell+1}\qquad\Longrightarrow\qquad\underline{u}_{k}^{\ell+1}\geq\underline{u}_{k}^{\mathfrak{n}}.$
(37)

The proof of the upper bound in (32) is similar, and (34) follows easily from
(32).

We next show (33). Rearranging (18b) for the temperature gives

$\displaystyle T_{i}^{\ell+1}$
$\displaystyle=T_{i}^{\mathfrak{n}}+\frac{\Delta
t}{\varepsilon}\sum_{j}\lambda_{i,j}^{\ell}\beta_{j,i}\left(T_{j}^{\ell+1}-T_{i}^{\ell+1}\right)+\frac{m_{i}}{d}(s_{i}^{\mathfrak{n}}-s_{i}^{\ell+1})+\frac{\Delta
t}{\varepsilon
d}\frac{1}{n_{i}}\left[\boldsymbol{s}^{(2),\ell+1,\ell}\right]_{i}$
$\displaystyle\geq T_{i}^{\mathfrak{n}}+\frac{\Delta
t}{\varepsilon}\sum_{j}\lambda_{i,j}^{\ell}\beta_{j,i}\left(T_{j}^{\ell+1}-T_{i}^{\ell+1}\right),$
(38)

where the inequality in (38) follows from Proposition D.3. The rest of the
proof of (33) proceeds in the same way as the velocity lower bound, starting
at (36a). ∎

######  Remark 3.9.

The results of Proposition 3.8 can be verified for the backward Euler step
directly. The proof is a slight modification of the one given above.

##  4 Analysis of Velocity Iterates

Let $\boldsymbol{W}_{\mathcal{N}}^{\ell}$ and
$\boldsymbol{W}_{\mathcal{R}}^{\ell}$ be the null and range space components,
respectively, of $\boldsymbol{W}^{\ell}$ with respect to $Z^{\ell}$. Since
$Z^{\ell}$ is symmetric, its null space and range space are orthogonal, and
$\boldsymbol{W}^{\ell}=\boldsymbol{W}_{\mathcal{N}}^{\ell}+\boldsymbol{W}_{\mathcal{R}}^{\ell}$.
The following property will be used throughout the section.

######  Lemma 4.1 ([17], Section 4).

The null space $\mathcal{N}$ and range space $\mathcal{R}$ of the matrix
$Z^{\ell}$ are independent of $\ell$. The null space has dimension one and is
spanned by the vector $P^{\nicefrac{{1}}{{2}}}\mathbf{1}$.

It will be convenient for the analysis of velocity terms to use the Frobenius
norm

$\|\boldsymbol{W}\|_{\textnormal{F}}^{2}=\sum_{i=1}^{N}\sum_{j=1}^{d}w_{i,j}^{2}=\sum_{i=1}^{N}\|\boldsymbol{w}_{i}\|_{2}^{2}=\sum_{j=1}^{d}\left\|\boldsymbol{W}_{(\cdot),j}\right\|_{2}^{2}.$
(39)

###  4.1 Velocity Iterates on Null Space

######  Lemma 4.2.

The null space iterates are
$\boldsymbol{W}^{\ell}_{\mathcal{N}}=P^{\nicefrac{{1}}{{2}}}\boldsymbol{U}^{\infty}$,
where $\boldsymbol{U}^{\infty}=\mathbf{1}(\boldsymbol{u}^{\infty})^{\top}$. In
particular, they are independent of $\ell$. As a result,
$\|\Delta^{\ell}\boldsymbol{W}\|_{\textnormal{F}}=\|\Delta^{\ell}(\boldsymbol{W}_{\mathcal{R}})\|_{\textnormal{F}}$.

###### Proof.

According to Lemma 4.1,
$\boldsymbol{W}_{\mathcal{N}}^{\ell}=P^{\nicefrac{{1}}{{2}}}\mathbf{1}(\boldsymbol{c}^{\ell})^{\top}$,
for some $\boldsymbol{c}^{\ell}\in\mathbb{R}^{d}$. Since
$\boldsymbol{W}^{\ell}_{\mathcal{N}}\perp\boldsymbol{W}^{\ell}_{\mathcal{R}}$,
it follows that
$\mathbf{1}^{\top}P^{\nicefrac{{1}}{{2}}}\boldsymbol{W}^{\ell}=\mathbf{1}^{\top}P^{\nicefrac{{1}}{{2}}}\boldsymbol{W}^{\ell}_{\mathcal{N}}=\mathbf{1}^{\top}P\mathbf{1}(\boldsymbol{c}^{\ell})^{\top}.$
Hence,

$\boldsymbol{c}^{\ell}=\frac{(\boldsymbol{W}^{\ell})^{\top}P^{\nicefrac{{1}}{{2}}}\mathbf{1}}{\mathbf{1}^{\top}P\mathbf{1}}=\frac{(\boldsymbol{U}^{\ell})^{\top}P\mathbf{1}}{\mathbf{1}^{\top}P\mathbf{1}}=\frac{\sum_{i}\rho_{i}\boldsymbol{u}_{i}^{\ell}}{\sum_{i}\rho_{i}}=\boldsymbol{u}^{\infty},$
(40)

which is independent of $\ell$ by Proposition 3.7. Thus,
$\boldsymbol{W}^{\ell}_{\mathcal{N}}=P^{\nicefrac{{1}}{{2}}}\mathbf{1}(\boldsymbol{u}^{\infty})^{\top}=P^{\nicefrac{{1}}{{2}}}\boldsymbol{U}^{\infty}$
is independent of $\ell$. As a consequence,

$\|\Delta^{\ell}\boldsymbol{W}\|_{\textnormal{F}}^{2}=\|\Delta^{\ell}(\boldsymbol{W}_{\mathcal{R}})\|_{\textnormal{F}}^{2}+\|\Delta^{\ell}(\boldsymbol{W}_{\mathcal{N}})\|_{\textnormal{F}}^{2}=\|\Delta^{\ell}(\boldsymbol{W}_{\mathcal{R}})\|_{\textnormal{F}}^{2}.$
(41)

∎

###  4.2 Velocity Iterates on Range Space and Initial Norm Bounds

######  Definition 4.3.

Let the columns of
$V=[\boldsymbol{v}_{1},\cdots,\boldsymbol{v}_{N-1}]\in\mathbb{R}^{N\times(N-1)}$
form an orthonormal basis of $\mathcal{R}$. For each $\ell$, let
$\boldsymbol{A}^{\ell}=[\boldsymbol{a}^{\ell}_{1},\cdots,\boldsymbol{a}^{\ell}_{d}]\in\mathbb{R}^{(N-1)\times
d}$ be the unique matrix of coefficients such that
$\boldsymbol{W}^{\ell}_{\mathcal{R}}=V\boldsymbol{A}^{\ell}$.

######  Lemma 4.4.

Given $V$ as in 4.3 and any $\boldsymbol{Y}\in\mathbb{R}^{N\times d}$,

$\|V\|_{2}=1,\qquad\|V^{\top}\|_{2}\leq
1,\qquad\|V^{\top}\boldsymbol{Y}\|_{\textnormal{F}}\leq\|\boldsymbol{Y}\|_{\textnormal{F}},$
(42)

where $\|\cdot\|_{2}$ denotes the induced (operator) 2-norm when the argument
is an operator. In particular, setting $\boldsymbol{Y}=\boldsymbol{A}$ and
$\boldsymbol{Y}=\Delta^{\ell}\boldsymbol{A}$, respectively, gives

$\displaystyle\|\boldsymbol{W}^{\ell}_{\mathcal{R}}\|_{\textnormal{F}}=\|V\boldsymbol{A}^{\ell}\|_{\textnormal{F}}=\|\boldsymbol{A}^{\ell}\|_{\textnormal{F}},$
(43a)
$\displaystyle\|\Delta^{\ell}\boldsymbol{W}\|_{\textnormal{F}}=\|\Delta^{\ell}(\boldsymbol{W}_{\mathcal{R}})\|_{\textnormal{F}}=\|V\Delta^{\ell}(\boldsymbol{A})\|_{\textnormal{F}}=\|\Delta^{\ell}\boldsymbol{A}\|_{\textnormal{F}}.$
(43b)

Based on (43), we seek to bound $\|\boldsymbol{A}^{\ell}\|_{\textnormal{F}}$,
and the difference $\|\Delta^{\ell}\boldsymbol{A}\|_{\textnormal{F}}$.

######  Lemma 4.5.

The matrix iterates $\boldsymbol{A}^{\ell+1}$ satisfy

$\boldsymbol{A}^{\ell+1}=M^{\ell}V^{\top}\boldsymbol{W}^{\mathfrak{n}}\qquad\mbox{and}\qquad\Delta^{\ell}\boldsymbol{A}=\Delta^{\ell-1}(M)V^{\top}\boldsymbol{W}^{\mathfrak{n}},$
(44)

where $M^{\ell}=\left(I+\frac{\Delta
t}{\varepsilon}V^{\top}Z^{\ell}V\right)^{-1}$.

###### Proof.

Write
$\boldsymbol{W}^{\ell+1}=\boldsymbol{W}_{\mathcal{R}}^{\ell+1}+\boldsymbol{W}_{\mathcal{N}}^{\ell+1}=V\boldsymbol{A}^{\ell+1}+\boldsymbol{W}_{\mathcal{N}}^{\ell+1}$
in (20a) and multiply on the left by $V^{\top}$:

$V^{\top}\left(I+\frac{\Delta
t}{\varepsilon}Z^{\ell}\right)\left(V\boldsymbol{A}^{\ell+1}+\boldsymbol{W}_{\mathcal{N}}^{\ell+1}\right)=V^{\top}\boldsymbol{W}^{\mathfrak{n}}.$
(45)

Since
$Z^{\ell}\boldsymbol{W}_{\mathcal{N}}^{\ell+1}=\mathbf{0}\in\mathbb{R}^{N\times
d}$ and
$V^{\top}\boldsymbol{W}_{\mathcal{N}}^{\ell+1}=\mathbf{0}\in\mathbb{R}^{(N-1)\times
d}$, the result follows. ∎

######  Lemma 4.6.

Let $\gamma^{\mathfrak{n}}=\frac{\varepsilon}{\varepsilon+\Delta
tz_{\min}^{\mathfrak{n}}}$, where $z_{\min}^{\mathfrak{n}}>0$ is a lower bound
on the positive eigenvalues of $Z^{\mathfrak{n}}$, defined in Proposition 3.4.
Then

$\|M^{\ell}\|_{2}\leq\gamma^{\mathfrak{n}}\qquad\mbox{and}\qquad\|\Delta^{\ell-1}M\|_{2}\leq\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\|\Delta^{\ell-1}{Z}\|_{2},$
(46)

where $\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}=\frac{\Delta
t}{\varepsilon}(\gamma^{\mathfrak{n}})^{2}$.

###### Proof.

The first bound in (46) is verified in Appendix A. For the second, the
identity $X^{-1}-Y^{-1}=X^{-1}(Y-X)Y^{-1}$ with $X^{-1}=M^{\ell}$ and
$Y^{-1}=M^{\ell-1}$ implies

$\Delta^{\ell-1}M=-\frac{\Delta
t}{\varepsilon}M^{\ell}V^{\top}\Delta^{\ell-1}(Z)VM^{\ell-1}.$ (47)

The result follows after taking operator norms. ∎

######  Lemma 4.7.

The matrix iterates $\boldsymbol{A}^{\ell}$ satisfy

$\displaystyle\|\boldsymbol{A}^{\ell+1}\|_{\textnormal{F}}$
$\displaystyle\leq\gamma^{n}\|\boldsymbol{W}^{n}\|_{\textnormal{F}},$ (48a)
$\displaystyle\|\Delta^{\ell}\boldsymbol{A}\|_{\textnormal{F}}$
$\displaystyle\leq\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\|\boldsymbol{W}^{\mathfrak{n}}\|_{\textnormal{F}}\|\Delta^{\ell-1}Z\|_{2},$
(48b)

where $\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}=\frac{\Delta
t}{\varepsilon}(\gamma^{\mathfrak{n}})^{2}$.

###### Proof.

The result follows using Lemmas 4.4, 4.5 and 4.6. ∎

######  Lemma 4.8.

The difference of matrices $Z^{\ell}$ satisfy

$\|\Delta^{\ell-1}Z\|_{2}\leq
C_{Z}^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2},\qquad\mbox{where}\qquad
C_{Z}^{\mathfrak{n}}=\frac{2(N+N^{\nicefrac{{1}}{{2}}})C_{A}^{\mathfrak{n}}}{\rho_{\min}n_{\min}^{\nicefrac{{1}}{{2}}}}.$
(49)

###### Proof.

The proof is contained in Appendix B. ∎

###  4.3 Summary of Velocity Term Analysis

######  Theorem 4.9.

The iterations of the velocity terms satisfy the following bound:

$\|\Delta^{\ell}\boldsymbol{W}\|_{\textnormal{F}}=\|\Delta^{\ell}\boldsymbol{A}\|_{\textnormal{F}}\leq
C_{\boldsymbol{W}}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2},$
(50)

where
$C_{\boldsymbol{W}}^{\mathfrak{n}}=C_{Z}^{\mathfrak{n}}\|\boldsymbol{W}^{\mathfrak{n}}\|_{\textnormal{F}}$
and $\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}=\frac{\Delta
t\varepsilon}{(\varepsilon+\Delta tz_{\min}^{\mathfrak{n}})^{2}}$, with
$C_{Z}^{\mathfrak{n}}$ as in (49) and $z_{\min}^{\mathfrak{n}}>0$ defined in
Proposition 3.4.

###### Proof.

The result follows using (43b), (48b), and (49). ∎

##  5 Analysis of Temperature Iterations

In this section, we analyze the temperature iterations; the results are
summarized by Theorem 5.15. The basic strategy is similar to that of the
previous section, but there are important differences. Denote the null space
and range space of $\widehat{Z}^{\ell}$ by $\widehat{\mathcal{N}}$ and
$\widehat{\mathcal{R}}$, respectively, and let
$\boldsymbol{\eta}^{\ell}=\boldsymbol{\eta}_{\mathcal{N}}^{\ell}+\boldsymbol{\eta}_{\mathcal{R}}^{\ell}$,
where $\boldsymbol{\eta}_{\mathcal{N}}^{\ell}\in\widehat{\mathcal{N}}$ and
$\boldsymbol{\eta}_{\mathcal{R}}^{\ell}\in\widehat{\mathcal{R}}$. As before,
the fact that $\widehat{\mathcal{R}}\perp\widehat{\mathcal{N}}$ implies that
$\|\Delta^{\ell}\boldsymbol{\eta}\|_{2}^{2}=\|\Delta^{\ell}\boldsymbol{\eta}_{\mathcal{N}}\|_{2}^{2}+\|\Delta^{\ell}\boldsymbol{\eta}_{\mathcal{R}}\|_{2}^{2}$.
However, unlike $\boldsymbol{W}^{\ell}_{\mathcal{N}}$, the iterates
$\boldsymbol{\eta}^{\ell}_{\mathcal{N}}$ depend on $\ell$, and thus
$\|\boldsymbol{\eta}^{\ell+1}_{\mathcal{N}}-\boldsymbol{\eta}^{\ell}_{\mathcal{N}}\|_{2}$
must be bounded. In addition, there are source terms in the update for
$\boldsymbol{\eta}^{\ell}$ that depend on velocity iterates and require
technical bounds.

###  5.1 Temperature Iterates on Null Space

######  Lemma 5.1 ([17], Sec. 4).

The null space $\widehat{\mathcal{N}}$ and range space $\widehat{\mathcal{R}}$
of $\widehat{Z}^{\ell}$ are independent of $\ell$. The null space has
dimension one and is spanned by the vector
$Q^{\nicefrac{{1}}{{2}}}\mathbf{1}$.

######  Lemma 5.2.

The null space iterations $\boldsymbol{\eta}^{\ell}_{\mathcal{N}}$ satisfy

$\displaystyle\left\|\Delta^{\ell}(\boldsymbol{\eta}_{\mathcal{N}})\right\|_{2}\leq
C_{0}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2},\qquad\mbox{where}\qquad$
(51) $\displaystyle
C_{0}^{\mathfrak{n}}=\frac{\sqrt{2}b}{d}\rho_{\max}^{\nicefrac{{1}}{{2}}}u_{\max}^{\mathfrak{n}}N^{\nicefrac{{1}}{{2}}}C_{\boldsymbol{W}}^{\mathfrak{n}},\qquad
b=\|Q^{\nicefrac{{1}}{{2}}}\mathbf{1}\|_{2}^{-1},\qquad\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}=\frac{\Delta
t\varepsilon}{(\varepsilon+\Delta tz_{\min}^{\mathfrak{n}})^{2}}.$ (52)

###### Proof.

Write
$\boldsymbol{\eta}_{\mathcal{N}}^{\ell+1}=b_{0}^{\ell+1}\widehat{\boldsymbol{v}}_{0},$
where $\widehat{\boldsymbol{v}}_{0}=bQ^{\nicefrac{{1}}{{2}}}\mathbf{1}$
satisfies $\|\widehat{\boldsymbol{v}}_{0}\|_{2}=1$. Then

$\left\|\Delta^{\ell}(\boldsymbol{\eta}_{\mathcal{N}})\right\|_{2}=\left\|\Delta^{\ell}(b_{0})\widehat{\boldsymbol{v}}_{0}\right\|_{2}=\left|\Delta^{\ell}(b_{0})\right|.$
(53)

Now multiply (20b) by $\widehat{\boldsymbol{v}}_{0}^{\top}$. Since
$\widehat{\mathcal{R}}\perp\widehat{\mathcal{N}}$, the left-hand side gives

$\widehat{\boldsymbol{v}}_{0}^{\top}\left(I+\frac{\Delta
t}{\varepsilon}\widehat{Z}^{\ell}\right)(b_{0}^{\ell+1}\widehat{\boldsymbol{v}}_{0}+\boldsymbol{\eta}_{\mathcal{R}}^{\ell+1})=b_{0}^{\ell+1}\widehat{\boldsymbol{v}}_{0}^{\top}\widehat{\boldsymbol{v}}_{0}=b_{0}^{\ell+1}.$
(54)

Meanwhile, the right-hand side gives

$\left(\widehat{\boldsymbol{v}}_{0},\boldsymbol{\eta}^{\mathfrak{n}}\right)_{2}+\frac{1}{d}\left(\widehat{\boldsymbol{v}}_{0},Q^{-\nicefrac{{1}}{{2}}}\left(\boldsymbol{s}_{\boldsymbol{W}}^{\mathfrak{n}}-\boldsymbol{s}_{\boldsymbol{W}}^{\ell+1}\right)\right)_{2}+\frac{\Delta
t}{\varepsilon}\frac{1}{d}\left(\widehat{\boldsymbol{v}}_{0},Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}^{(2),\ell+1,\ell}\right)_{2}\\\
=\left(\widehat{\boldsymbol{v}}_{0},\boldsymbol{\eta}^{\mathfrak{n}}\right)_{2}+\frac{1}{d}\left(\widehat{\boldsymbol{v}}_{0},Q^{-\nicefrac{{1}}{{2}}}\left(\boldsymbol{s}_{\boldsymbol{W}}^{\mathfrak{n}}-\boldsymbol{s}_{\boldsymbol{W}}^{\ell+1}\right)\right)_{2},$
(55)

where the second source term vanishes due to the symmetry of $B$.
Specifically, using (22b) and the definition of
$\widehat{\boldsymbol{v}}_{0}$,

$\left(\widehat{\boldsymbol{v}}_{0},Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}^{(2),\ell+1,\ell}\right)_{2}=b\sum_{i,j}B_{i,j}^{\ell}\left[m_{j}(s_{j}^{\ell+1}-S_{i,j}^{\ell+1})-m_{i}(s_{i}^{\ell+1}-S_{i,j}^{\ell+1})\right]=0.$
(56)

Combining (54) and (55) gives

$b_{0}^{\ell+1}=\left(\widehat{\boldsymbol{v}}_{0},\boldsymbol{\eta}^{\mathfrak{n}}\right)_{2}+\frac{b}{d}\left(\mathbf{1},\boldsymbol{s}_{\boldsymbol{W}}^{\mathfrak{n}}\right)_{2}-\frac{b}{d}\left(\mathbf{1},\boldsymbol{s}_{\boldsymbol{W}}^{\ell+1}\right)_{2}.$
(57)

Therefore, with the definition of $\boldsymbol{s}_{\boldsymbol{W}}$ in (22a),

$\begin{split}\left|\Delta^{\ell}(b_{0})\right|&=\frac{b}{d}\left|\left(\mathbf{1},\Delta^{\ell}(\boldsymbol{s}_{\boldsymbol{W}})\right)_{2}\right|\leq\frac{b}{d}\sum_{i=1}^{N}\left|\;\left\|\boldsymbol{w}_{i}^{\ell+1}\right\|_{2}^{2}-\left\|\boldsymbol{w}_{i}^{\ell}\right\|_{2}^{2}\;\right|\\\
&\leq\frac{b}{d}\sum_{i=1}^{N}\left\|\Delta^{\ell}(\boldsymbol{w}_{i})\right\|_{2}\left\|\boldsymbol{w}_{i}^{\ell+1}+\boldsymbol{w}_{i}^{\ell}\right\|_{2}\leq\frac{2b}{d}\rho_{\max}^{\nicefrac{{1}}{{2}}}u_{\max}^{\mathfrak{n}}\sum_{i=1}^{N}\left\|\Delta^{\ell}(\boldsymbol{w}_{i})\right\|_{2}\\\
&\leq\frac{2b}{d}\rho_{\max}^{\nicefrac{{1}}{{2}}}u_{\max}^{\mathfrak{n}}\sqrt{N\sum_{i=1}^{N}\left\|\Delta^{\ell}(\boldsymbol{w}_{i})\right\|^{2}_{2}}=\frac{2b}{d}\rho_{\max}^{\nicefrac{{1}}{{2}}}u_{\max}^{\mathfrak{n}}N^{\nicefrac{{1}}{{2}}}\|\Delta^{\ell}\boldsymbol{W}\|_{\textnormal{F}}\\\
&\leq\frac{2b}{d}\rho_{\max}^{\nicefrac{{1}}{{2}}}u_{\max}^{\mathfrak{n}}N^{\nicefrac{{1}}{{2}}}C_{\boldsymbol{W}}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2},\end{split}$
(58)

where the last line follows from Theorem 4.9. ∎

###  5.2 Bounding Source Terms in the Temperature Equation

######  Lemma 5.3.

The vector
$\boldsymbol{s}_{\boldsymbol{W}}^{\ell}=\left(|\boldsymbol{w}_{1}^{\ell}|^{2},\cdots,|\boldsymbol{w}_{N}^{\ell}|^{2}\right)^{\top}$
satisfies the bounds

$\|\boldsymbol{s}_{\boldsymbol{W}}^{\ell}\|_{2}\leq(u_{\max}^{\mathfrak{n}})^{2}\|P\mathbf{1}\|_{2}\qquad\mbox{and}\qquad\|\Delta^{\ell}(\boldsymbol{s}_{\boldsymbol{W}})\|_{2}\leq
C_{\boldsymbol{s}}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2},$
(59)

where
$C_{\boldsymbol{s}}^{\mathfrak{n}}=2\rho_{\max}^{\nicefrac{{1}}{{2}}}u_{\max}^{\mathfrak{n}}C_{\boldsymbol{W}}^{\mathfrak{n}}$
and $\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}=\frac{\Delta
t\varepsilon}{(\varepsilon+\Delta tz_{\min}^{\mathfrak{n}})^{2}}$.

###### Proof.

The first inequality follows from the bound on the velocity terms (34) in
Proposition 3.8:

$\|\boldsymbol{s}_{\boldsymbol{W}}^{\ell}\|_{2}^{2}=\sum_{i}[\boldsymbol{s}_{\boldsymbol{W}}]_{i}^{2}=\sum_{i}\rho_{i}^{2}\|\boldsymbol{u}_{i}^{\ell}\|_{2}^{4}\leq(u_{\max}^{\mathfrak{n}})^{4}\sum_{i}\rho_{i}^{2}.$
(60)

For the second inequality,

$\begin{split}\|\Delta^{\ell}(\boldsymbol{s}_{\boldsymbol{W}})\|_{2}^{2}&=\sum_{i}\left|\left(\Delta^{\ell}(\boldsymbol{w}_{i}),\boldsymbol{w}_{i}^{\ell+1}+\boldsymbol{w}_{i}^{\ell}\right)_{2}\right|^{2}\leq\sum_{i}\|\Delta^{\ell}(\boldsymbol{w}_{i})\|_{2}^{2}\;\|\boldsymbol{w}_{i}^{\ell+1}+\boldsymbol{w}_{i}^{\ell}\|_{2}^{2}\\\
&\leq
4\rho_{\max}(u_{\max}^{\mathfrak{n}})^{2}\|\Delta^{\ell}\boldsymbol{W}\|_{\textnormal{F}}^{2}.\end{split}$
(61)

Applying the bound in Theorem 4.9 and taking square roots gives the result. ∎

######  Lemma 5.4.

$\displaystyle\left\|\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}_{\infty}\right\|_{2}$
$\displaystyle\leq\frac{\|\boldsymbol{W}^{\mathfrak{n}}\|_{\textnormal{F}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}\;\gamma^{\mathfrak{n}},$
(62a) $\displaystyle\left|s_{i}^{\ell+1}-S_{i,j}^{\ell+1}\right|$
$\displaystyle\leq\frac{4\alpha_{j,i}u_{\max}^{\mathfrak{n}}\|\boldsymbol{W}^{\mathfrak{n}}\|_{\textnormal{F}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}\;\gamma^{\mathfrak{n}},$
(62b)

where $\gamma^{\mathfrak{n}}$ is defined in Lemma 4.6.

###### Proof.

We first prove (62a). From (43a) and the bound in (48a),

$\|\boldsymbol{W}_{\mathcal{R}}^{\ell+1}\|_{\textnormal{F}}=\|\boldsymbol{A}^{\ell+1}\|_{\textnormal{F}}\leq\gamma^{\mathfrak{n}}\|\boldsymbol{W}^{\mathfrak{n}}\|_{\textnormal{F}}.$
(63)

Then using Lemma 4.2 and the fact that
$\boldsymbol{W}_{\mathcal{R}}^{\ell+1}=\boldsymbol{W}^{\ell+1}-\boldsymbol{W}_{\mathcal{N}}$,

$\|\boldsymbol{U}^{\ell+1}-\boldsymbol{U}_{\infty}\|_{\textnormal{F}}^{2}=\|P^{-\nicefrac{{1}}{{2}}}(\boldsymbol{W}^{\ell+1}-\boldsymbol{W}_{\mathcal{N}})\|_{\textnormal{F}}^{2}\leq\frac{1}{\rho_{\min}}\|\boldsymbol{W}_{\mathcal{R}}^{\ell+1}\|_{\textnormal{F}}^{2}\overset{\eqref{eq:velocity_norm_0}}{\leq}\frac{1}{\rho_{\min}}(\gamma^{\mathfrak{n}})^{2}\|\boldsymbol{W}^{\mathfrak{n}}\|_{\textnormal{F}}^{2}.$
(64)

Therefore,

$\|\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}_{\infty}\|_{2}\leq\|\boldsymbol{U}^{\ell+1}-\boldsymbol{U}^{\infty}\|_{\textnormal{F}}\leq\gamma^{\mathfrak{n}}\frac{\|\boldsymbol{W}^{\mathfrak{n}}\|_{\textnormal{F}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}.$
(65)

We next prove (62b). Since $\alpha_{i,j}+\alpha_{j,i}=1$,

$s_{i}^{\ell+1}-S_{i,j}^{\ell+1}=\alpha_{j,i}\left(\boldsymbol{u}_{i}^{\ell+1}+\boldsymbol{u}_{i,j}^{\ell+1},\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}_{j}^{\ell+1}\right)_{2}.$
(66)

Thus, using (62a),

$\begin{split}\left|s_{i}^{\ell+1}-S_{i,j}^{\ell+1}\right|&\leq\max\\{\alpha_{i,j}\\}\left[\left\|\boldsymbol{u}_{i}^{\ell+1}\right\|_{2}+\left\|\boldsymbol{u}_{i,j}^{\ell+1}\right\|_{2}\right]\left\|\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}^{\infty}+\boldsymbol{u}^{\infty}-\boldsymbol{u}_{j}^{\ell+1}\right\|_{2}\\\
&\leq
4\max\\{\alpha_{i,j}\\}u_{\max}^{\mathfrak{n}}\left(\frac{\|\boldsymbol{W}^{\mathfrak{n}}\|_{\textnormal{F}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}\right)\gamma^{\mathfrak{n}}.\end{split}$
(67)

∎

######  Lemma 5.5.

For all $\ell\geq 1$,

$\displaystyle\|\boldsymbol{s}^{(2),\ell,\ell-1}\|_{2}\leq
C_{\boldsymbol{s}^{(2)}_{0}}^{\mathfrak{n}}\gamma^{\mathfrak{n}},\quad\mbox{where}\quad
C_{\boldsymbol{s}^{(2)}_{0}}^{\mathfrak{n}}=\frac{8\max\\{\alpha_{i,j}\\}B_{\max}u_{\max}^{\mathfrak{n}}m_{\max}N^{\nicefrac{{1}}{{2}}}(N-1)\|\boldsymbol{W}^{\mathfrak{n}}\|_{\textnormal{F}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}.$
(68)

###### Proof.

$\begin{split}\|&\boldsymbol{s}^{(2),\ell,\ell-1}\|_{2}^{2}=\sum_{i}\bigg{|}\sum_{j\neq
i}B_{i,j}^{\ell}\big{[}m_{j}(s_{j}^{\ell+1}-S_{i,j}^{\ell+1})-m_{i}(s_{i}^{\ell+1}-S_{i,j}^{\ell+1})\big{]}\bigg{|}^{2}\\\
&\leq\sum_{i,j\neq i}\sum_{k\neq
i}B_{i,j}B_{i,k}\bigg{[}m_{i}^{2}\big{|}s_{i}^{\ell+1}-S_{i,j}^{\ell+1}\big{|}\big{|}s_{i}^{\ell+1}-S_{i,k}^{\ell+1}\big{|}+m_{i}m_{j}\big{|}s_{i}^{\ell+1}-S_{i,k}^{\ell+1}\big{|}\big{|}s_{j}^{\ell+1}-S_{i,j}^{\ell+1}\big{|}\\\
&\qquad+m_{i}m_{k}\big{|}s_{i}^{\ell+1}-S_{i,j}^{\ell+1}\big{|}\big{|}s_{k}^{\ell+1}-S_{i,k}^{\ell+1}\big{|}+m_{j}m_{k}\big{|}s_{j}^{\ell+1}-S_{i,j}^{\ell+1}\big{|}\big{|}s_{k}^{\ell+1}-S_{i,k}^{\ell+1}\big{|}\bigg{]}\\\
&\leq
B_{\max}^{2}m_{\max}^{2}\;4\;\left(\frac{4\max\\{\alpha_{i,j}\\}u_{\max}^{\mathfrak{n}}\|\boldsymbol{W}^{\mathfrak{n}}\|_{\textnormal{F}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}\gamma^{\mathfrak{n}}\right)^{2}N(N-1)^{2}.\end{split}$
(69)

Combining constants and taking the square root gives the result. ∎

######  Lemma 5.6.

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

$\displaystyle\left|s_{k}^{\ell+1}-s_{k}^{\ell}\right|$
$\displaystyle\leq\frac{2u_{\max}^{\mathfrak{n}}C_{\boldsymbol{W}}^{\mathfrak{n}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2},$
(70a)
$\displaystyle\text{and}\qquad\left|S_{i,j}^{\ell+1}-S_{i,j}^{\ell}\right|$
$\displaystyle\leq\frac{2u_{\max}^{\mathfrak{n}}C_{\boldsymbol{W}}^{\mathfrak{n}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2}.$
(70b)

###### Proof.

We first prove (70a). As a simple corollary to Theorem 4.9,

$\begin{split}\|\boldsymbol{u}_{k}^{\ell+1}-\boldsymbol{u}_{k}^{\ell}\|_{2}^{2}&\leq\sum_{i=1}^{N}\|\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}_{i}^{\ell}\|_{2}^{2}=\|\boldsymbol{U}^{\ell+1}-\boldsymbol{U}^{\ell}\|_{\textnormal{F}}^{2}=\|P^{-\nicefrac{{1}}{{2}}}\left(\Delta^{\ell}\boldsymbol{W}\right)\|_{\textnormal{F}}^{2}\\\
&\overset{\ref{theorem:VelocityError}}{\leq}\frac{1}{\rho_{\min}}(C_{\boldsymbol{W}}^{\mathfrak{n}})^{2}(\Gamma_{\boldsymbol{W}}^{\mathfrak{n}})^{2}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2}^{2}.\end{split}$
(71)

Thus,

$\left|s_{k}^{\ell+1}-s_{k}^{\ell}\right|=\left|\left(\boldsymbol{u}_{k}^{\ell+1}-\boldsymbol{u}_{k}^{\ell},\boldsymbol{u}_{k}^{\ell+1}+\boldsymbol{u}_{k}^{\ell}\right)_{2}\right|\leq
2u_{\max}^{\mathfrak{n}}\|\boldsymbol{u}_{k}^{\ell+1}-\boldsymbol{u}_{k}^{\ell}\|_{2}\leq\frac{2u_{\max}^{\mathfrak{n}}C_{\boldsymbol{W}}^{\mathfrak{n}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2},$
(72)

as desired. Next, note that

$\begin{split}\left|S_{i,j}^{\ell+1}-S_{i,j}^{\ell}\right|&\;\;=\;\;\left|\left(\boldsymbol{u}_{i,j}^{\ell+1}-\boldsymbol{u}_{i,j}^{\ell},\boldsymbol{u}_{i,j}^{\ell+1}+\boldsymbol{u}_{i,j}^{\ell}\right)_{2}\right|\leq
2u_{\max}^{\mathfrak{n}}\left\|\boldsymbol{u}_{i,j}^{\ell+1}-\boldsymbol{u}_{i,j}^{\ell}\right\|_{2}\\\
&\;\;=\;\;2u_{\max}^{\mathfrak{n}}\left\|\alpha_{i,j}(\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}_{i}^{\ell})+\alpha_{j,i}(\boldsymbol{u}_{j}^{\ell+1}-\boldsymbol{u}_{j}^{\ell})\right\|_{2}\\\
&\;\;\leq\;\;2u_{\max}^{\mathfrak{n}}\left[\alpha_{i,j}\left\|\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}_{i}^{\ell}\right\|_{2}+\alpha_{j,i}\left\|\boldsymbol{u}_{j}^{\ell+1}-\boldsymbol{u}_{j}^{\ell}\right\|_{2}\right]\\\
&\overset{\eqref{eq:u_k_ell1_ell_diff}}{\leq}2u_{\max}^{\mathfrak{n}}(\alpha_{i,j}+\alpha_{j,i})\frac{C_{\boldsymbol{W}}^{\mathfrak{n}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2}\\\
&\;\;=\;\;\frac{2u_{\max}^{\mathfrak{n}}C_{\boldsymbol{W}}^{\mathfrak{n}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2}.\end{split}$
(73)

∎

######  Lemma 5.7.

$\left\|\boldsymbol{s}^{(2),\ell+1,\ell}-\boldsymbol{s}^{(2),\ell,\ell-1}\right\|_{2}\leq
C_{\boldsymbol{s}^{(2)}_{1}}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2}+C_{\boldsymbol{s}^{(2)}_{2}}^{\mathfrak{n}}\gamma^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2},$
(74)

for
$C_{\boldsymbol{s}^{(2)}_{1}}^{\mathfrak{n}}=\frac{6B_{\max}u_{\max}^{\mathfrak{n}}m_{\max}N(N-1)C_{\boldsymbol{W}}^{\mathfrak{n}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}},$
and
$C_{\boldsymbol{s}^{(2)}_{2}}^{\mathfrak{n}}=C_{B}^{\mathfrak{n}}\frac{16\max\\{\alpha_{i,j}\\}u_{\max}^{\mathfrak{n}}m_{\max}(N-1)\|\boldsymbol{W}^{\mathfrak{n}}\|_{\textnormal{F}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}$.

###### Proof.

Let $\xi_{i,j}^{k}=m_{j}(s_{j}^{k}-S_{i,j}^{k})-m_{i}(s_{i}^{k}-S_{i,j}^{k})$,
so that $\left[\boldsymbol{s}^{(2),k_{1},k_{2}}\right]_{i}=\sum_{j\neq
i}B_{i,j}^{k_{2}}\xi_{i,j}^{k_{1}}$. Then using the norm equivalence of the 2-
and 1- norms, adding and subtracting $B_{i,j}^{\ell}\xi_{i,j}^{\ell}$, and
using triangle inequality,

$\begin{split}\left\|\boldsymbol{s}^{(2),\ell+1,\ell}-\boldsymbol{s}^{(2),\ell,\ell-1}\right\|_{2}&\leq\left\|\boldsymbol{s}^{(2),\ell+1,\ell}-\boldsymbol{s}^{(2),\ell,\ell-1}\right\|_{1}=\sum_{i=1}^{N}\bigg{|}\sum\limits_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{N}B_{i,j}^{\ell}\xi_{i,j}^{\ell+1}-B_{i,j}^{\ell-1}\xi_{i,j}^{\ell}\bigg{|}\\\
&\leq\sum_{i,j\neq
i}B_{i,j}^{\ell}\left|\xi_{i,j}^{\ell+1}-\xi_{i,j}^{\ell}\right|+\sum_{i,j\neq
i}\left|B_{i,j}^{\ell}-B_{i,j}^{\ell-1}\right|\left|\xi_{i,j}^{\ell}\right|.\end{split}$
(75)

Using Lemma 5.6, the first term of (75) can be bounded as follows:

$\displaystyle\sum_{i,j\neq i}$ $\displaystyle
B_{i,j}^{\ell}|\xi_{i,j}^{\ell+1}-\xi_{i,j}^{\ell}|\leq
B_{\max}^{\mathfrak{n}}\sum_{i,j\neq
i}\bigg{|}m_{j}(s_{j}^{\ell+1}-S_{i,j}^{\ell+1})-m_{i}(s_{i}^{\ell+1}-S_{i,j}^{\ell+1})-m_{j}(s_{j}^{\ell}-S_{i,j}^{\ell})+m_{i}(s_{i}^{\ell}-S_{i,j}^{\ell})\bigg{|}$
$\displaystyle\leq B_{\max}^{\mathfrak{n}}m_{\max}\sum_{i,j\neq
i}\left[|s_{j}^{\ell+1}-s_{j}^{\ell}|+|s_{i}^{\ell+1}-s_{i}^{\ell}|+|S_{i,j}^{\ell+1}-S_{i,j}^{\ell}|\right]$
$\displaystyle\leq\frac{6B_{\max}^{\mathfrak{n}}u_{\max}^{\mathfrak{n}}m_{\max}N(N-1)C_{\boldsymbol{W}}^{\mathfrak{n}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2}.$
(76)

Using Lemmas B.2 and 5.4, the second term of (75) can be bounded as follows:

$\begin{split}\sum_{i,j\neq
i}\left|B_{i,j}^{\ell}-B_{i,j}^{\ell-1}\right|\left|\xi_{i,j}^{\ell}\right|&\leq
C_{B}^{\mathfrak{n}}\left[\gamma^{\mathfrak{n}}\frac{8\max\\{\alpha_{i,j}\\}u_{\max}^{\mathfrak{n}}m_{\max}\|\boldsymbol{W}^{\mathfrak{n}}\|_{\textnormal{F}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}\right]\sum_{i,j\neq
i}\left(|\Delta^{\ell-1}T_{i}|+|\Delta^{\ell-1}T_{j}|\right)\\\
&=C_{B}^{\mathfrak{n}}\frac{8\max\\{\alpha_{i,j}\\}m_{\max}u_{\max}^{\mathfrak{n}}\|\boldsymbol{W}^{n}\|_{\textnormal{F}}}{\rho_{\min}^{\nicefrac{{1}}{{2}}}}2(N-1)\gamma^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2}.\end{split}$
(77)

The result follows upon putting together the terms in (76) and (77). ∎

###  5.3 Temperature Iterates on Range Space

######  Definition 5.8.

Let the columns of
$\widehat{V}=[\widehat{\boldsymbol{v}}_{1},\cdots,\widehat{\boldsymbol{v}}_{N-1}]\in\mathbb{R}^{N\times(N-1)}$
form an orthonormal basis of $\widehat{\mathcal{R}}$. For each $\ell$, let
$\boldsymbol{b}^{\ell}\in\mathbb{R}^{N-1}$ be the unique vector of
coefficients such that
$\boldsymbol{\eta}_{\mathcal{R}}^{\ell}=\widehat{V}\boldsymbol{b}^{\ell}$.

######  Lemma 5.9.

Given $\widehat{V}\in\mathbb{R}^{N\times(N-1)}$ as in 5.8 and any
$\boldsymbol{y}\in\mathbb{R}^{N}$,

$\|\widehat{V}\|_{2}=1,\qquad\|\widehat{V}^{\top}\|_{2}\leq
1,\qquad\|\widehat{V}^{\top}\boldsymbol{y}\|_{2}\leq\|\boldsymbol{y}\|_{2}.$
(78)

######  Lemma 5.10.

The iterates $\boldsymbol{b}^{\ell+1}$ satisfy

$\boldsymbol{b}^{\ell+1}=\widehat{M}^{\ell}\left[\widehat{V}^{\top}\boldsymbol{\eta}^{\mathfrak{n}}+\frac{1}{d}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\left(\boldsymbol{s}_{\boldsymbol{W}}^{\mathfrak{n}}-\boldsymbol{s}_{\boldsymbol{W}}^{\ell+1}\right)+\frac{\Delta
t}{\varepsilon
d}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}^{(2),\ell+1,\ell}\right],$
(79)

where $\widehat{M}^{\ell}=\left(I+\frac{\Delta
t}{\varepsilon}\widehat{V}^{\top}\widehat{Z}^{\ell}\widehat{V}\right)^{-1}$.
Therefore

$\begin{split}\left\|\Delta^{\ell}\boldsymbol{b}\right\|_{2}&\leq\left\|\Delta^{\ell-1}\widehat{M}\widehat{V}^{\top}\left[\boldsymbol{\eta}^{\mathfrak{n}}+\frac{1}{d}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}_{\boldsymbol{W}}^{\mathfrak{n}}\right]\right\|_{2}\\\
&\quad+\frac{1}{d}\left\|\widehat{M}^{\ell}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}_{\boldsymbol{W}}^{\ell+1}-\widehat{M}^{\ell-1}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}_{\boldsymbol{W}}^{\ell}\right\|_{2}\\\
&\quad+\frac{\Delta t}{\varepsilon
d}\left\|\widehat{M}^{\ell}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}^{(2),\ell+1,\ell}-\widehat{M}^{\ell-1}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}^{(2),\ell,\ell-1}\right\|_{2}.\end{split}$
(80)

###### Proof.

Write
$\boldsymbol{\eta}^{\ell+1}=\boldsymbol{\eta}_{\mathcal{R}}^{\ell+1}+\boldsymbol{\eta}_{\mathcal{N}}^{\ell+1}=\widehat{V}\boldsymbol{b}^{\ell+1}+\boldsymbol{\eta}_{\mathcal{N}}^{\ell+1}$,
multiply (20b) by $\widehat{V}^{\top}$ to cancel the null space component
$\boldsymbol{\eta}_{\mathcal{N}}^{\ell+1}$, then multiply by $\widehat{M}$ to
obtain (79). To obtain (80), take the norm of the difference of consecutive
iterates, and apply the triangle inequality. ∎

The remainder of this subsection is devoted to bounding the terms in (80).

######  Lemma 5.11.

Let $\widehat{\gamma}^{\mathfrak{n}}=\frac{\varepsilon}{\varepsilon+\Delta
t\widehat{z}_{\min}^{\mathfrak{n}}}$, where
$\widehat{z}_{\min}^{\mathfrak{n}}>0$ is a lower bound on the positive
eigenvalues of $\widehat{Z}^{\ell}$, defined in Proposition 3.4. Then the
matrices $\widehat{M}^{\ell}$ and $\Delta^{\ell-1}\widehat{M}$ satisfy the
bounds

$\displaystyle\|\widehat{M}^{\ell}\|_{2}$
$\displaystyle\leq\widehat{\gamma}^{\mathfrak{n}},$ (81a)
$\displaystyle\|\Delta^{\ell-1}(\widehat{M})\|_{2}$
$\displaystyle\leq\Gamma_{\boldsymbol{\eta}}^{\mathfrak{n}}\|\Delta^{\ell-1}\widehat{Z}\|_{2}\leq
C_{\widehat{Z}}^{\mathfrak{n}}\Gamma_{\boldsymbol{\eta}}^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2},$
(81b)

where

$C_{\widehat{Z}}^{\mathfrak{n}}=\frac{2\max_{i,j}\\{c_{i,j}^{B}\\}m_{\max}^{\nicefrac{{1}}{{2}}}N^{\nicefrac{{1}}{{2}}}(N^{\nicefrac{{1}}{{2}}}+1)}{\sqrt{2}m_{\min}n_{\min}^{\nicefrac{{3}}{{2}}}\sqrt{T_{\min}^{\mathfrak{n}}}}\qquad\mbox{and}\qquad\Gamma_{\boldsymbol{\eta}}^{\mathfrak{n}}=\frac{\Delta
t}{\varepsilon}\left(\widehat{\gamma}^{\mathfrak{n}}\right)^{2}.$ (82)

###### Proof.

The bound in (81a) is verified in Appendix A. To prove (81b), use (81a), the
matrix identity $X^{-1}-Y^{-1}=X^{-1}(Y-X)Y^{-1}$, Lemma 5.9, and Appendix B.
∎

######  Lemma 5.12.

The first term of (80) satisfies

$\left\|\Delta^{\ell-1}\widehat{M}\left[\boldsymbol{\eta}^{\mathfrak{n}}+\frac{1}{d}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}_{\boldsymbol{W}}^{\mathfrak{n}}\right]\right\|_{2}\leq
C_{1,1}^{\mathfrak{n}}\Gamma_{\boldsymbol{\eta}}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2},$
(83)

where
$C_{1,1}^{\mathfrak{n}}=\left\|\boldsymbol{\eta}^{\mathfrak{n}}+\frac{1}{d}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}_{\boldsymbol{W}}^{\mathfrak{n}}\right\|_{2}C_{\widehat{Z}}^{\mathfrak{n}}$.

###### Proof.

Use Lemma 5.11. ∎

######  Lemma 5.13.

The second term of (80) satisfies

$\frac{1}{d}\left\|\widehat{M}^{\ell}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}_{\boldsymbol{W}}^{\ell+1}-\widehat{M}^{\ell-1}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}_{\boldsymbol{W}}^{\ell}\right\|_{2}\leq
C_{2,1}^{\mathfrak{n}}\widehat{\gamma}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2}+C_{2,2}^{\mathfrak{n}}\Gamma_{\boldsymbol{\eta}}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2},$
(84)

where
$C_{2,1}^{\mathfrak{n}}=\frac{C_{\boldsymbol{s}}^{\mathfrak{n}}}{dn_{\min}^{\nicefrac{{1}}{{2}}}}$
(with $C_{\boldsymbol{s}}^{\mathfrak{n}}$ defined in Lemma 5.3) and
$C_{2,2}^{\mathfrak{n}}=\frac{(u_{\max}^{\mathfrak{n}})^{2}\|P\mathbf{1}\|_{2}C_{\widehat{Z}}^{\mathfrak{n}}}{dn_{\min}^{\nicefrac{{1}}{{2}}}}$.

###### Proof.

Add and subtract
$\widehat{M}^{\ell}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}_{\boldsymbol{W}}^{\ell}$,
factor, use standard inequalities and apply Lemmas 5.11 and 5.3:

$\begin{split}\frac{1}{d}&\left\|\widehat{M}^{\ell}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}_{\boldsymbol{W}}^{\ell+1}-\widehat{M}^{\ell-1}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}_{\boldsymbol{W}}^{\ell}\right\|_{2}\\\
&\quad\leq\frac{1}{d}\left\|\widehat{M}^{\ell}\right\|\left\|Q^{-\nicefrac{{1}}{{2}}}\right\|_{2}\left\|\Delta^{\ell}(\boldsymbol{s}_{\boldsymbol{W}})\right\|_{2}+\frac{1}{d}\left\|\Delta^{\ell-1}\widehat{M}\right\|_{2}\left\|Q^{-\nicefrac{{1}}{{2}}}\right\|_{2}\left\|\boldsymbol{s}_{\boldsymbol{W}}^{\ell}\right\|_{2}\\\
&\quad\leq\frac{1}{d}\widehat{\gamma}^{\mathfrak{n}}\frac{1}{n_{\min}^{\nicefrac{{1}}{{2}}}}C_{\boldsymbol{s}}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2}+\frac{\|P\mathbf{1}\|_{2}(u_{\max}^{\mathfrak{n}})^{2}}{dn_{\min}^{\nicefrac{{1}}{{2}}}}C_{\widehat{Z}}^{\mathfrak{n}}\Gamma_{\boldsymbol{\eta}}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2}.\end{split}$
(85)

∎

######  Lemma 5.14.

The third term of (80) satisfies

$\begin{split}\frac{\Delta t}{\varepsilon
d}&\left\|\widehat{M}^{\ell}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}^{(2),\ell+1,\ell}-\widehat{M}^{\ell-1}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}^{(2),\ell,\ell-1}\right\|_{2}\\\
&\leq
C_{3,1}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2}+C_{3,2}^{\mathfrak{n}}(\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}\Gamma_{\boldsymbol{\eta}}^{\mathfrak{n}})^{\nicefrac{{1}}{{2}}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2}+C_{3,3}^{\mathfrak{n}}\Gamma_{\boldsymbol{\eta}}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2},\end{split}$
(86)

where

$C_{3,1}^{\mathfrak{n}}=\frac{C_{\boldsymbol{s}_{1}^{(2)}}^{\mathfrak{n}}}{dn_{\min}^{\nicefrac{{1}}{{2}}}}\frac{1-\widehat{\gamma}^{\mathfrak{n}}}{\widehat{z}^{\mathfrak{n}}_{\min}},\qquad
C_{3,2}^{\mathfrak{n}}=\frac{C_{\boldsymbol{s}_{2}^{(2)}}^{\mathfrak{n}}}{dn_{\min}^{\nicefrac{{1}}{{2}}}},\qquad\mbox{and}\qquad
C_{3,3}^{\mathfrak{n}}=\frac{C_{\boldsymbol{s}_{0}^{(2)}}^{\mathfrak{n}}C_{\widehat{Z}}^{\mathfrak{n}}}{dn_{\min}^{\nicefrac{{1}}{{2}}}}\frac{1-\widehat{\gamma}^{\mathfrak{n}}}{\widehat{z}^{\mathfrak{n}}_{\min}}.$
(87)

###### Proof.

Add and subtract
$\widehat{M}^{\ell}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}^{(2),\ell,\ell-1}$,
and factor:

$\begin{split}&\frac{\Delta t}{\varepsilon
d}\left\|\widehat{M}^{\ell}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}^{(2),\ell+1,\ell}-\widehat{M}^{\ell-1}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}^{(2),\ell,\ell-1}\right\|_{2}\\\
&\leq\frac{\Delta t}{\varepsilon
d}\left\|\widehat{M}^{\ell}\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\left(\boldsymbol{s}^{(2),\ell+1,\ell}-\boldsymbol{s}^{(2),\ell,\ell-1}\right)\right\|_{2}+\frac{\Delta
t}{\varepsilon
d}\left\|\Delta^{\ell-1}(\widehat{M})\,\widehat{V}^{\top}Q^{-\nicefrac{{1}}{{2}}}\boldsymbol{s}^{(2),\ell,\ell-1}\right\|_{2}\end{split}$
(88)

The first term of (88) can be bounded using Lemmas 5.7 and 5.11:

$\begin{split}\frac{\Delta t}{\varepsilon
d}\left\|\widehat{M}^{\ell}\right\|_{2}\left\|Q^{-\nicefrac{{1}}{{2}}}\right\|_{2}&\left\|\boldsymbol{s}^{(2),\ell+1,\ell}-\boldsymbol{s}^{(2),\ell,\ell-1}\right\|_{2}\\\
&\leq\frac{1}{dn_{\min}^{\nicefrac{{1}}{{2}}}}\frac{\Delta
t}{\varepsilon}\widehat{\gamma}^{\mathfrak{n}}\left[C_{\boldsymbol{s}^{(2)}_{1}}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2}+C_{\boldsymbol{s}^{(2)}_{2}}^{\mathfrak{n}}\widehat{\gamma}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2}\right]\\\
&=\frac{C_{\boldsymbol{s}^{(2)}_{1}}^{\mathfrak{n}}}{dn_{\min}^{\nicefrac{{1}}{{2}}}}\frac{1-\widehat{\gamma}^{\mathfrak{n}}}{\widehat{z}_{\min}^{\mathfrak{n}}}\Gamma_{\boldsymbol{W}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2}+\frac{C_{\boldsymbol{s}^{(2)}_{2}}^{\mathfrak{n}}}{dn_{\min}^{\nicefrac{{1}}{{2}}}}(\Gamma_{\boldsymbol{W}}\Gamma_{\boldsymbol{\eta}})^{\nicefrac{{1}}{{2}}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2}.\end{split}$
(89)

The second term of (88) can be bounded using Lemmas 5.11 and 5.5:

$\frac{\Delta t}{\varepsilon
d}\left\|Q^{-\nicefrac{{1}}{{2}}}\right\|_{2}\left\|\boldsymbol{s}^{(2),\ell,\ell-1}\right\|_{2}\left\|\widehat{M}^{\ell}-\widehat{M}^{\ell-1}\right\|_{2}\leq\frac{C_{\boldsymbol{s}_{0}^{(2)}}^{\mathfrak{n}}}{dn_{\min}^{\nicefrac{{1}}{{2}}}}\frac{1-\widehat{\gamma}^{\mathfrak{n}}}{\widehat{z}_{\min}^{\mathfrak{n}}}C_{\widehat{Z}}^{\mathfrak{n}}\Gamma_{\boldsymbol{\eta}}^{\mathfrak{n}}\left\|\Delta^{\ell-1}\boldsymbol{\eta}\right\|_{2}.$
(90)

Putting together (89) and (90), the result follows. ∎

###  5.4 Summary of Temperature Term Analysis

The results of the above sections are summarized by the following:

######  Theorem 5.15.

The temperature iterations satisfy the following bounds:

$\|\Delta^{\ell}\boldsymbol{\eta}\|_{2}\leq\bigg{[}C_{0}^{\mathfrak{n}}\Gamma^{\mathfrak{n}}_{\boldsymbol{W}}+C_{1,1}^{\mathfrak{n}}\Gamma^{\mathfrak{n}}_{\boldsymbol{\eta}}+C_{2,1}^{\mathfrak{n}}\Gamma^{\mathfrak{n}}_{\boldsymbol{W}}\widehat{\gamma}^{\mathfrak{n}}+C_{2,2}^{\mathfrak{n}}\Gamma^{\mathfrak{n}}_{\boldsymbol{\eta}}+C_{3,1}^{\mathfrak{n}}\Gamma^{\mathfrak{n}}_{\boldsymbol{W}}+C_{3,2}^{\mathfrak{n}}(\Gamma^{\mathfrak{n}}_{\boldsymbol{W}}\Gamma^{\mathfrak{n}}_{\boldsymbol{\eta}})^{\nicefrac{{1}}{{2}}}+C_{3,3}^{\mathfrak{n}}\Gamma^{\mathfrak{n}}_{\boldsymbol{\eta}}\bigg{]}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2},$
(91)

where the constants, $C_{(\cdot)}$ can be found using the following table:

Constant: | $C_{0}$ | $C_{1,1}$ | $C_{2,1},\;C_{2,2}$ | $C_{3,1},\;C_{3,2},\;C_{3,3}$
---|---|---|---|---
Lemma: | 5.2 | 5.12 | 5.13 | 5.14

###### Proof.

Since
$\|\Delta^{\ell}\boldsymbol{\eta}\|_{2}\leq\|\Delta^{\ell}\boldsymbol{\eta}_{\mathcal{N}}\|_{2}+\|\Delta^{\ell}\boldsymbol{\eta}_{\mathcal{R}}\|_{2}$,
then using Lemmas 5.2, 5.12, 5.13 and 5.14, the result follows. ∎

##  6 Convergence Criterion and Time Step Selection

In this section, we use Theorem 4.9 and 5.15 to establish the time step
restriction in Corollary 3.6.1

###  6.1 Preliminaries

For $\ell\in\mathbb{N}\cup\\{0\\},$ define the tuple
$\boldsymbol{X}^{\ell}\coloneqq\left(\boldsymbol{A}^{\ell},\boldsymbol{b}^{\ell},b_{0}^{\ell}\right)\in\mathbb{R}^{(N-1)\times
d}\times\mathbb{R}^{N-1}\times\mathbb{R}$ and the norm
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\boldsymbol{X}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|=\|\boldsymbol{A}\|_{\textnormal{F}}+\|\boldsymbol{b}\|_{2}+|b_{0}|$.

######  Theorem 6.1.

The Cauchy errors of the iterates $\boldsymbol{X}^{\ell}$ satisfy

$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\Delta^{\ell}\boldsymbol{X}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|\leq
C_{\boldsymbol{X}}^{\mathfrak{n}}\Gamma_{\boldsymbol{X}}^{\mathfrak{n}}(\Delta
t,\varepsilon)\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\Delta^{\ell-1}\boldsymbol{X}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|,$
(92)

where $\Gamma_{\boldsymbol{X}}^{\mathfrak{n}}=\frac{\Delta
t\varepsilon}{(\varepsilon+\Delta
tz^{\mathfrak{n}})^{2}},\,z^{\mathfrak{n}}=\min\\{z_{\min}^{\mathfrak{n}},\widehat{z}_{\min}^{\mathfrak{n}}\\}\eqqcolon
z$, and

$C_{\boldsymbol{X}}^{\mathfrak{n}}=C_{\boldsymbol{W}}^{\mathfrak{n}}+C_{0}^{\mathfrak{n}}+C_{1,1}^{\mathfrak{n}}+C_{2,1}^{\mathfrak{n}}+C_{2,2}^{\mathfrak{n}}+C_{3,1}^{\mathfrak{n}}+C_{3,2}^{\mathfrak{n}}+C_{3,3}^{\mathfrak{n}}.$
(93)

###### Proof.

The proof is simply a summary of the results in Theorems 4.9 and 5.15. ∎

######  Lemma 6.2.

The bound
$C_{\boldsymbol{X}}^{\mathfrak{n}}\Gamma_{\boldsymbol{X}}^{\mathfrak{n}}\leq
r$ is equivalent to $Q_{r}(\Delta t)\geq 0$, where

$Q_{r}(\Delta t)=rz^{2}\Delta t^{2}+(2rz-
C_{\boldsymbol{X}}^{\mathfrak{n}})\varepsilon\Delta t+r\varepsilon^{2}.$ (94)

###### Proof.

The proof is a simple calculation. ∎

######  Lemma 6.3.

If
$C_{\boldsymbol{X}}^{\mathfrak{n}}\Gamma_{\boldsymbol{X}}^{\mathfrak{n}}\leq
r<1,$ then convergence of the method is guaranteed.

###### Proof.

By Theorem 6.1, if
$C_{\boldsymbol{X}}^{\mathfrak{n}}\Gamma_{\boldsymbol{X}}^{\mathfrak{n}}<r<1,$
then

$\lim_{\ell\to\infty}\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\Delta^{\ell}\boldsymbol{X}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|\leq\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\boldsymbol{X}^{1}-\boldsymbol{X}^{0}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|\lim_{\ell\to\infty}r^{\ell}=0.$
(95)

∎

###  6.2 Time Step Selection

Denote the zeros of $Q_{r}$, defined in (94), by

$\Delta
t_{\pm}=\frac{C_{\boldsymbol{X}}^{\mathfrak{n}}-2rz}{2rz^{2}}\varepsilon\pm\frac{C_{\boldsymbol{X}}^{\mathfrak{n}}\sqrt{1-\frac{4rz}{C_{\boldsymbol{X}}^{\mathfrak{n}}}}}{2rz^{2}}\varepsilon.$
(96)

Since the $y$-intercept of $Q_{r}$ is positive, there are three possible
cases: (1) $\Delta t_{\pm}$ are complex conjugates, which implies that,
$Q_{r}\geq 0$ for any $\Delta t>0$; (2) $\Delta t_{\pm}$ are both real-valued
and negative, so that $Q_{r}>0$ for any $\Delta t>0$; (3) $\Delta t_{\pm}$ are
both real-valued and positive, in which case, selecting $\Delta t$ in the set
$S_{0}\coloneqq(0,\Delta t_{-}]\cup[\Delta t_{+},\infty)$ ensures that
$Q_{r}(\Delta t)\geq 0,$ and thus, the method will converge.

######  Theorem 6.4.

Given $\Delta t_{0}>0$, there exists $\Delta t_{1}$ such that (i) $C\Delta
t_{0}\leq\Delta t_{1}\leq\Delta t_{0}$ for a constant $C\in(0,1]$ that is
independent of $\varepsilon$, and (ii) for this $\Delta t_{1}$,

$C_{\boldsymbol{W}}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}+C_{0}^{\mathfrak{n}}\Gamma_{\boldsymbol{W}}^{\mathfrak{n}}+C_{1}^{\mathfrak{n}}\Gamma_{\boldsymbol{X}}^{\mathfrak{n}}<1,$
(97)

thus ensuring convergence of the GST method.

###### Proof.

The only nontrivial case to consider is when $\Delta t_{\pm}$ are both real-
valued and positive. If $\Delta t_{0}\leq\Delta t_{-}$ or $\Delta
t_{+}\leq\Delta t_{0}$, then the result follows trivially with $\Delta
t_{0}=\Delta t_{1}$, i.e., $C=1$. If $\Delta t_{-}<\Delta t_{0}<\Delta t_{+}$,
then setting $\Delta t_{1}=\Delta t_{-}$ gives, via (96),

$\frac{\Delta t_{1}}{\Delta t_{0}}>\frac{\Delta t_{-}}{\Delta
t_{+}}=\frac{C_{\boldsymbol{X}}^{\mathfrak{n}}-2rz-
C_{\boldsymbol{X}}^{\mathfrak{n}}\sqrt{1-\frac{4rz}{C_{\boldsymbol{X}}^{\mathfrak{n}}}}}{C_{\boldsymbol{X}}^{\mathfrak{n}}-2rz+C_{\boldsymbol{X}}^{\mathfrak{n}}\sqrt{1-\frac{4rz}{C_{\boldsymbol{X}}^{\mathfrak{n}}}}}\eqqcolon
C\in(0,1],$ (98)

and $C$ is independent of $\varepsilon$. ∎

##  7 Numerical demonstration

To demonstrate the application of the GST method, we consider a slab geometry,
in which the distribution functions $f_{i}$ are invariant in the $x_{2}$ and
$x_{3}$ direction. In this setting, (5) takes the form

$\frac{\partial f_{i}}{\partial
t}+v_{1}\partial_{x}f_{i}=\frac{1}{\varepsilon}\sum_{j}\lambda_{i,j}(M_{i,j}-f_{i}),\quad\forall
i\in\\{1,\cdots,N\\},$ (99)

where $[\boldsymbol{u}_{i,j}]_{2}=[\boldsymbol{u}_{i,j}]_{3}=0$, and to
simplify the notation, we let $x_{1}=x$, $[\boldsymbol{u}_{i}]_{1}=u_{i}$, and
$[\boldsymbol{u}_{i,j}]_{1}=u_{i,j}$. Rather than solve (99) in the full
three-dimensional velocity space, we apply a multi-species extension of Chu
reduction [11]. For each $i\in\\{1,\cdots,N\\}$, let

$\displaystyle g_{i}(x,v_{1},t)$
$\displaystyle=\int_{\mathbb{R}^{2}}f_{i}(x,\boldsymbol{v},t)\,\,\textnormal{d}v_{2}\,\,\textnormal{d}v_{3},$
(100a) $\displaystyle h_{i}(x,v_{1},t)$
$\displaystyle=\int_{\mathbb{R}^{2}}(v_{2}^{2}+v_{3}^{2})f_{i}(x,\boldsymbol{v},t)\,\,\textnormal{d}v_{2}\,\,\textnormal{d}v_{3}.$
(100b)

Then (99) can be reduced to a system of two equations for $g_{i}$ and $h_{i}$:

$\displaystyle\frac{\partial g_{i}}{\partial t}+v_{1}\frac{\partial
g_{i}}{\partial x}$
$\displaystyle=\frac{1}{\varepsilon}\sum_{j}\lambda_{i,j}(G_{i,j}-g_{i}),$
(101a) $\displaystyle\frac{\partial h_{i}}{\partial t}+v_{1}\frac{\partial
h_{i}}{\partial x}$
$\displaystyle=\frac{1}{\varepsilon}\sum_{j}\lambda_{i,j}(H_{i,j}-h_{i}),$
(101b)

where $G_{ij}(v_{1})=M_{n_{i},u_{i,j},T_{i,j}/m_{i}}^{(1)}(v_{1})$ is a
Maxwellian (refer to 2.1) with $d=1$, $H_{i,j}=2\frac{T_{i,j}}{m_{i}}G_{i,j}$,
and

$n_{i}=\int_{\mathbb{R}}g_{i}\,\,\textnormal{d}v_{1},\qquad
n_{i}u_{i}=\int_{\mathbb{R}}v_{1}g_{i}\,\,\textnormal{d}v_{1},\qquad\frac{3n_{i}}{m_{i}}T_{i}=\int_{\mathbb{R}}(v_{1}-u_{i})^{2}g_{i}\,\,\textnormal{d}v_{1}+\int_{\mathbb{R}}h_{i}\,\,\textnormal{d}v_{1}.$
(102)

We use the IMEX approach outlined in Section 2.2 to solve (101). The Butcher
tableaux for the IMEX method is taken from Section 2.6 of [3]. The implicit
step takes the form of a backward Euler method and the associated moment
updates for $n_{i}$, $u_{i}$, and $T_{i}$ take exactly the form of (13). Thus,
the GST method still forms the key component of the implicit step, and all of
the analysis above carries through without any modification. The spatial
discretization of the advection operator uses a second-order finite volume
scheme with minmod-type slope limiters. The velocity is discretized using
standard discrete velocity models based on evenly spaced points and a mid-
point quadrature rule. Details of these discretizations can be found in [18].
The strategy for selecting the time step is to use the natural hyperbolic
restriction of the advection operators in the equations. In practice, we
observe that the GST method converges without the restrictions set in Section
6.2. However, for general problems, these restrictions provide a rigorous
fall-back whenever the GST method fails to converge at the hyperbolic time
step.

###  7.1 Sod problem with two identical species

To verify the multi-species code, we perform a simulation of the Sod shock
tube problem [27], a Riemann problem for the compressible Euler equations of a
perfect gas. The initial data are
$(n_{\textnormal{L}},\boldsymbol{u}_{\textnormal{L}},T_{\textnormal{L}})=(1,\mathbf{0},1)$
for $x\leq 0$, and
$(n_{\textnormal{R}},\boldsymbol{u}_{\textnormal{R}},T_{\textnormal{R}})=(0.125,\mathbf{0},0.8)$
for $x\geq 0$. An analytical solution for the Riemann problem can be found in
[28]. For small values of $\varepsilon$, both the single species BGK equation
[12] and multi-species BGK (M-BGK) equations [16] with properly initialized
data formally recover the Sod solution. Because the slab problem (99) is
formulated in three velocity dimensions, the adiabatic index for the limiting
Euler equations is $\gamma=5/3$.

In the current test, we treat the left- and right-hand side of the domain as
different particle species, each having mass $m_{i}=1$. The kinetic initial
data is

$f_{1}(x,\boldsymbol{v},0)=\begin{cases}M_{n_{\textnormal{L}},\boldsymbol{u}_{\textnormal{L}},T_{\textnormal{L}}}^{(3)}&x\leq
0\\\
M_{0.001,\boldsymbol{u}_{\textnormal{R}},T_{\textnormal{R}}}^{(3)}&x>0\end{cases}\qquad\mbox{and}\qquad
f_{2}(x,\boldsymbol{v},0)=\begin{cases}M_{0.001,\boldsymbol{u}_{\textnormal{L}},T_{\textnormal{L}}}^{(3)}&x\leq
0\\\
M_{n_{\textnormal{R}},\boldsymbol{u}_{\textnormal{R}},T_{\textnormal{R}}}^{(3)}&x>0\end{cases},$
(103)

where $M_{n,\boldsymbol{u},\theta}^{(d)}$ denotes the $d$-dimensional
Maxwellian, defined in (3). Trace densities are used in (103) instead of voids
in order to avoid negativity in the numerical solutions.

The numerical solution is computed on a discretization of the phase space
$X\times V=[-1,1]\times[-10,10]$, using periodic boundary conditions and a
mesh with $N_{x}=256$ by $N_{v}=192$ cells of uniform size. The time step
selection is driven by the advection in (101): $\Delta
t=0.9\frac{2/Nx}{2\max\\{|v_{l}|\\}}$. Figure 1 shows single species and
multi-species simulations of the model with $\varepsilon=10^{-4}$, run to a
final time of $t_{\textnormal{F}}=0.2$, compared with the Sod solution for the
associated Euler equations. The results show that both BGK and M-BGK models
perform well as numerical solvers for the Sod problem in this highly
collisional regime.

(a) Total number density.

(b) Total velocity.

(c) Total temperature.

Figure 1: The total density, total bulk velocity, and total temperature
moments are compared to the solution of the Euler equations at time $t=0.2$,
with $\varepsilon=10^{-4}$.

###  7.2 Three Species Mixture: Ar-Kr-Xe

The goal of this test is to demonstrate the benefits of the IMEX method when
the problem becomes stiff (i.e., as $\varepsilon\to 0$). The following CFL-
type conditions on the time step are required for the explicit and IMEX
methods (which use a minmod-type slope limiter in the flux computation) to
ensure stability and positivity of the numerical solution:

$\Delta t^{\mathfrak{n}}\leq\begin{cases}\frac{\varepsilon\Delta
x}{2\varepsilon\max\\{|v_{l}|\\}+\Delta
x\Lambda^{\mathfrak{n}}}&\\!\\!\\!\\!\eqqcolon\Delta
t_{\textnormal{EX}}^{\mathfrak{n}}\\\ \frac{\Delta
x}{2\max\\{|v_{l}|\\}}&\\!\\!\\!\\!\eqqcolon\Delta
t_{\textnormal{IMEX}}\end{cases},\qquad\mbox{where}\qquad\Lambda^{\mathfrak{n}}=\max_{i=1}^{N}\left\\{\max_{\kappa}\left\\{\sum_{j=1}^{N}\lambda_{i,j}^{\mathfrak{n}}(x_{\kappa})\right\\}\right\\}.$
(104)

As $\varepsilon\to 0$, the condition on $\Delta
t_{\textnormal{EX}}^{\mathfrak{n}}$ becomes restrictive.

To compare the methods, we simulate an interface problem, with a mixture of
Argon, Krypton, and Xenon; the masses and diameters of these particle species
are taken from [10] and are given below in SI units:

$\begin{pmatrix}m_{\textnormal{Ar}}\\\ m_{\textnormal{Kr}}\\\
m_{\textnormal{Xe}}\end{pmatrix}=\begin{pmatrix}6.6335209\\\ 13.914984\\\
21.801714\end{pmatrix}\times 10^{-26}\textnormal{
kg},\qquad\begin{pmatrix}d_{\textnormal{Ar}}\\\ d_{\textnormal{Kr}}\\\
d_{\textnormal{Xe}}\end{pmatrix}=\begin{pmatrix}3.659\\\ 4.199\\\
4.939\end{pmatrix}\times 10^{-10}\textnormal{ m}.$ (105)

The setup for the test is similar to that of Section 7.1: a mixture of gases
at rest is initialized with a high temperature on the left, and a low
temperature on the right of the domain. Specifically, the mixture of Ar, Kr,
and Xe is initialized on the left half of the 1-dimensional spatial domain
with number densities
$(n_{\textnormal{Ar}}^{0},n_{\textnormal{Kr}}^{0},n_{\textnormal{Xe}}^{0})=(5,5,0.5)$
$\times 10^{25}$ m-3, and common temperature $T_{\textnormal{L}}=10$; on the
right half, the mixture is initialized with number densities
$(n_{\textnormal{Ar}}^{0},n_{\textnormal{Kr}}^{0},n_{\textnormal{Xe}}^{0})=(0.5,0.5,5)$
$\times 10^{25}$ m-3, and the common temperature $T_{\textnormal{R}}=1$. The
initial kinetic distributions are given by

$\displaystyle f_{\textnormal{Ar}}(x,\boldsymbol{v},0)$
$\displaystyle=\chi_{[-1,0]}(x)M_{5,\mathbf{0},10/m_{\textnormal{Ar}}}^{(3)}(\boldsymbol{v})+\chi_{(0,1]}(x)M_{0.5,\mathbf{0},1/m_{\textnormal{Ar}}}^{(3)}(\boldsymbol{v}),$
(106a) $\displaystyle f_{\textnormal{Kr}}(x,\boldsymbol{v},0)$
$\displaystyle=\chi_{[-1,0]}(x)M_{5,\mathbf{0},10/m_{\textnormal{Kr}}}^{(3)}(\boldsymbol{v})+\chi_{(0,1]}(x)M_{0.5,\mathbf{0},1/m_{\textnormal{Kr}}}^{(3)}(\boldsymbol{v}),$
(106b) $\displaystyle f_{\textnormal{Xe}}(x,\boldsymbol{v},0)$
$\displaystyle=\chi_{[-1,0]}(x)M_{0.5,\mathbf{0},10/m_{\textnormal{Xe}}}^{(3)}(\boldsymbol{v})+\chi_{(0,1]}(x)M_{5,\mathbf{0},1/m_{\textnormal{Xe}}}^{(3)}(\boldsymbol{v}),$
(106c)

where $\chi_{\mathcal{S}}$ is the indicator function on the set $\mathcal{S}$,
and $M_{n,\boldsymbol{u},\theta}^{(d)}$ denotes the $d$-dimensional
Maxwellian, defined in (3).

The numerical solution is computed on a discretization of the phase space
$X\times V=[-1,1]\times[-30,30]$, using reflective boundaries; the mesh
consists of $N_{x}=256\times N_{v}=192$ cells of uniform size. The IMEX method
(using the GST solver) uses a time step of size $\Delta t=0.9\Delta
t_{\textnormal{IMEX}}$, and the fully explicit method uses a time step of size
$\Delta t^{\mathfrak{n}}=0.9\Delta t_{\textnormal{EX}}^{\mathfrak{n}}$; the
simulation is run to a final time of $t_{\textnormal{F}}=0.1$ $\mu$s.

Moment plots for numerical solutions computed with the IMEX and fully explicit
method are compared in Figure 2 for $\varepsilon\approx 0.001$. The computed
solutions align closely. However, the explicit method requires $7089$ time
steps, while the IMEX method only requires $854$ steps. In practice, as
$\varepsilon$ decreases, the number of time steps for the explicit method
grows, slowing the simulation, while the IMEX method experiences no slow down
in computation time.

(a) Density of Ar

(b) Bulk velocity of Ar

(c) Temperature of Ar

(d) Density of Kr

(e) Bulk velocity of Kr

(f) Temperature of Kr

(g) Density of Xe

(h) Bulk velocity of Xe

(i) Temperature of Xe

Figure 2: Noble gas mixture of Ar, Kr, and Xe: Moment plots of species number
density, bulk velocity, and temperature, are given at $t_{\textnormal{F}}=0.1$
$\mu$s, for a mixture with $\varepsilon\approx 0.00103$. M-BGK uses the GST
algorithm for the moment solve. M-BGK${}_{\textnormal{EX}}$ uses a fully
explicit method.

##  8 Conclusions

In this paper, we have presented a simple iterative approach for solving a set
of implicit moment equations as part of a larger IMEX-based strategy for
computing solutions of a multi-species BGK model. The main challenge is that
the collision frequencies in the moment equations themselves depend on the
moments in a complicated way. The key contribution is a rigorous contraction
mapping proof, which shows that the iterative method converges to a unique
solution under reasonable time steps that do not degenerate to zero when the
collision parameter $\varepsilon\to 0$. This implicit solve is a key component
in the IMEX discretization of the multi-species BGK model, the implementation
of which is done in slab geometry but closed in three-dimensional velocity
space using a Chu reduction technique [11].

In future work, we will examine similar problems for multi-species ellipsoidal
statistical models (ES-BGK) [19, 2] that, unlike standard BGK models, recover
in highly collisional regimes, the correct Prandtl number for the compressible
Navier-Stokes equations. We will also extend the current code structure used
here to generate results for more general problems in multi-dimensional
domains. We will also explore a Fokker-Planck analog of the BGK model, often
referred to as the Lenard-Bernstein model [21], that is useful in plasma
physics simulations [13, 14, 29].

##  Appendix A Inverse Matrix Norms

The following lemma is used in the proofs of Lemma 4.6 and Lemma 5.11.

######  Lemma A.1.

The matrices $M^{\ell}$ and $\widehat{M}^{\ell}$, defined in Lemma 4.5 and
Lemma 5.10, respectively, satisfy the bounds

$\|M^{\ell}\|_{2}\leq\frac{\varepsilon}{\varepsilon+\Delta
tz_{1}^{\ell}}\leq\frac{\varepsilon}{\varepsilon+\Delta
tz_{\min}^{\mathfrak{n}}}\qquad\mbox{and}\qquad\|\widehat{M}^{\ell}\|_{2}\leq\frac{\varepsilon}{\varepsilon+\Delta
t\widehat{z}_{1}^{\ell}}\leq\frac{\varepsilon}{\varepsilon+\Delta
t\widehat{z}_{\min}^{\mathfrak{n}}},$ (107)

where $z_{1}^{\ell}\geq z_{\min}^{\mathfrak{n}}$ and
$\widehat{z}_{1}^{\ell}\geq z_{\min}^{\mathfrak{n}}$ are defined in
Proposition 3.4.

###### Proof.

We verify the bound for $M^{\ell}$; the proof of the bound on
$\widehat{M}^{\ell}$ is nearly identical. Let $M^{\ell}=(I+J^{\ell})^{-1}$,
where $J^{\ell}=\frac{\Delta t}{\varepsilon}V^{\top}Z^{\ell}V$ and
$z_{1}^{\ell}>0$ is the minimum positive eigenvalue of $Z^{\ell}$. Then
$J^{\ell}$ is symmetric, and for any $\boldsymbol{x}\in\mathbb{R}^{N-1}$,

$\boldsymbol{x}^{\top}J^{\ell}\boldsymbol{x}\geq\frac{\Delta
t}{\varepsilon}z_{1}^{\ell}\|\boldsymbol{x}\|_{2}^{2}\qquad\mbox{and}\qquad\|J^{\ell}\boldsymbol{x}\|^{2}_{2}=\boldsymbol{x}^{\top}(J^{\ell})^{2}\boldsymbol{x}\geq\left(\frac{{\Delta
t}}{\varepsilon}z_{1}^{\ell}\right)^{2}\|\boldsymbol{x}\|_{2}^{2}.$ (108)

Given $\boldsymbol{y}\in\mathbb{R}^{N-1}$, let
$\boldsymbol{x}=M^{\ell}\boldsymbol{y}=(I+J^{\ell})^{-1}\boldsymbol{y}$. Then
$\boldsymbol{y}=(I+J^{\ell})\boldsymbol{x}$ and, using (108),

$\|\boldsymbol{y}\|^{2}_{2}=\|\boldsymbol{x}\|^{2}_{2}+2\boldsymbol{x}^{\top}(J^{\ell})\boldsymbol{x}+\|J^{\ell}\boldsymbol{x}\|^{2}_{2}\geq\|\boldsymbol{x}\|^{2}_{2}+2\frac{\Delta
t}{\varepsilon}z_{1}^{\ell}\|\boldsymbol{x}\|^{2}_{2}+\left(\frac{{\Delta
t}}{\varepsilon}z_{1}^{\ell}\right)^{2}\|\boldsymbol{x}\|^{2}_{2}=\left(1+\frac{\Delta
t}{\varepsilon}z_{1}^{\ell}\right)^{2}\|\boldsymbol{x}\|^{2}_{2}.$ (109)

Since $\boldsymbol{y}$ is arbitrary and $z_{\min}^{\mathfrak{n}}\leq
z_{1}^{\ell}$, the bounds on $M^{\ell}$ follow. ∎

##  Appendix B Bounding Differences of $Z^{\ell}$ and $\widehat{Z}^{\ell}$
Matrices

The following lemma is a consequence of the difference of squares formula.

######  Lemma B.1.

Suppose that $x,y,\kappa\in\mathbb{R}$ such that $0<\kappa\leq x,y$. Then

$\left|\sqrt{x}-\sqrt{y}\right|\leq\frac{1}{2\sqrt{\kappa}}|x-y|.$ (110)

######  Lemma B.2 (Bounding $|\Delta^{\ell-1}A|$ and $|\Delta^{\ell-1}B|$ ).

Suppose $T_{i}^{\ell+1}\geq\min_{k}\\{T_{k}^{\mathfrak{n}}\\}\eqqcolon
T_{\min}^{\mathfrak{n}}$. Then

$\displaystyle|\Delta^{\ell-1}(A_{i,j})|$ $\displaystyle\leq
C^{\mathfrak{n}}_{A}\left(|\Delta^{\ell-1}(T_{i})|+|\Delta^{\ell-1}(T_{j})|\right),$
(111a) $\displaystyle|\Delta^{\ell-1}(B_{i,j})|$ $\displaystyle\leq
C_{B}^{\mathfrak{n}}\left(|\Delta^{\ell-1}(T_{i})|+|\Delta^{\ell-1}(T_{j})|\right),$
(111b)

where

$C_{A}^{\mathfrak{n}}=\frac{\max_{i,j}\\{c_{i,j}^{A}\\}\;m_{\max}^{\nicefrac{{1}}{{2}}}}{2\sqrt{2}m_{\min}\sqrt{T_{\min}^{\mathfrak{n}}}}\qquad\mbox{and}\qquad
C_{B}^{\mathfrak{n}}=\frac{\max_{i,j}\\{c_{i,j}^{B}\\}\;m_{\max}^{\nicefrac{{1}}{{2}}}}{2\sqrt{2}m_{\min}\sqrt{T_{\min}^{\mathfrak{n}}}}.$
(112)

###### Proof.

We prove only the bound in (111a). The proof of (111b) is nearly identical.
Let $\Psi^{\ell}_{i,j}=\Psi_{i,j}(T_{i}^{\ell},T_{j}^{\ell})$ where
$\Psi_{i,j}$ is given in (9). Then due to the lower bound on the temperature,

$x\coloneqq(\Psi^{\ell}_{i,j})^{2},\qquad
y\coloneqq(\Psi^{\ell-1}_{i,j})^{2},\quad\mbox{and}\quad\kappa\coloneqq\frac{2T_{\min}^{\mathfrak{n}}}{m_{\max}}$
(113)

satisfy the conditions of Lemma B.1. Therefore,

$\left|\Delta^{\ell-1}(\Psi_{i,j})\right|\leq\frac{1}{2}\sqrt{\frac{m_{\max}}{2T_{\min}^{\mathfrak{n}}}}\;\left|\frac{\Delta^{\ell-1}(T_{i})}{m_{i}}+\frac{\Delta^{\ell-1}(T_{j})}{m_{j}}\right|,$
(114)

and hence, using 3.3,

$|\Delta^{\ell-1}(A_{i,j})|=c_{i,j}^{A}|\Delta^{\ell-1}(\Psi_{i,j})|\leq
C^{\mathfrak{n}}_{A}\left(|\Delta^{\ell-1}(T_{i})|+|\Delta^{\ell-1}(T_{j})|\right).$
(115)

∎

######  Proposition B.3 (Bounding $\|\Delta^{\ell-1}Z\|_{2}$ and
$\|\Delta^{\ell-1}\widehat{Z}\|_{2}$).

Suppose that $T_{i}^{\ell+1}\geq\min_{k}\\{T_{k}^{\mathfrak{n}}\\}\eqqcolon
T_{\min}^{\mathfrak{n}}$. Then

$\|\Delta^{\ell-1}Z\|_{2}\leq
C_{Z}^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2}\qquad\mbox{and}\qquad\|\Delta^{\ell-1}\widehat{Z}\|_{2}\leq
C_{\widehat{Z}}^{\mathfrak{n}}\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2},$ (116)

where

$C_{Z}^{\mathfrak{n}}=\frac{2(N+N^{\nicefrac{{1}}{{2}}})C_{A}^{\mathfrak{n}}}{\rho_{\min}n_{\min}^{\nicefrac{{1}}{{2}}}}\qquad\mbox{and}\qquad
C_{\widehat{Z}}^{\mathfrak{n}}=\frac{2(N+N^{\nicefrac{{1}}{{2}}})C_{B}^{\mathfrak{n}}}{n_{\min}^{\nicefrac{{3}}{{2}}}}.$
(117)

###### Proof.

We prove the first inequality; the second is similar. Since
$Z=P^{-\nicefrac{{1}}{{2}}}(D-A)P^{-\nicefrac{{1}}{{2}}}$ is positive semi-
definite, then for $\boldsymbol{y}=P^{-\nicefrac{{1}}{{2}}}\boldsymbol{x}$,

$0\leq\boldsymbol{x}^{\top}\Delta^{\ell-1}Z\boldsymbol{x}=\boldsymbol{y}^{\top}\Delta^{\ell-1}(D-A)\boldsymbol{y}\qquad\Longrightarrow\qquad\frac{\boldsymbol{y}^{\top}\Delta^{\ell-1}A\boldsymbol{y}}{\boldsymbol{y}^{\top}\boldsymbol{y}}\leq\frac{\boldsymbol{y}^{\top}\Delta^{\ell-1}D\boldsymbol{y}}{\boldsymbol{y}^{\top}\boldsymbol{y}}.$
(118)

Thus $\|\Delta^{\ell-1}A\|_{2}\leq\|\Delta^{\ell-1}D\|_{2}$, and using Lemma
B.2 in the second line below gives,

$\begin{split}&\|\Delta^{\ell-1}Z\|_{2}\leq\frac{1}{\rho_{\min}}\|\Delta^{\ell-1}(D-A)\|_{2}\leq\frac{2}{\rho_{\min}}\|\Delta^{\ell-1}D\|_{2}\\\
&\quad=\frac{2}{\rho_{\min}}\max_{i}\big{|}\sum_{k}\Delta^{\ell-1}(A_{i,k})\big{|}\leq\frac{2C_{A}^{\mathfrak{n}}}{\rho_{\min}}\max_{i}\sum_{k}(|\Delta^{\ell-1}(T_{i})|+|\Delta^{\ell-1}(T_{k})|)\\\
&\quad\leq\frac{2C_{A}^{\mathfrak{n}}}{\rho_{\min}}N\|\Delta^{\ell-1}\boldsymbol{T}\|_{2}+\frac{2C_{A}^{\mathfrak{n}}}{\rho_{\min}}N^{\nicefrac{{1}}{{2}}}\|\Delta^{\ell-1}\boldsymbol{T}\|_{2}\leq\frac{2C_{A}^{\mathfrak{n}}}{\rho_{\min}}(N+N^{\nicefrac{{1}}{{2}}})\frac{\|\Delta^{\ell-1}\boldsymbol{\eta}\|_{2}}{n_{\min}^{\nicefrac{{1}}{{2}}}},\end{split}$
(119)

where the last line follow from the relation
$\boldsymbol{\eta}=Q^{\nicefrac{{1}}{{2}}}\boldsymbol{T}$ with
$Q_{j,k}=n_{k}\delta_{j,k}$. ∎

##  Appendix C Invariance of $\alpha_{i,j}$ and $\beta_{i,j}$

######  Lemma C.1.

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

$\alpha_{i,j}^{\ell}=\frac{\rho_{i}\lambda_{i,j}^{\ell}}{\rho_{i}\lambda_{i,j}^{\ell}+\rho_{j}\lambda_{j,i}^{\ell}}\qquad\mbox{and}\qquad\beta_{i,j}^{\ell}=\frac{n_{i}\lambda_{i,j}^{\ell}}{n_{i}\lambda_{i,j}^{\ell}+n_{j}\lambda_{j,i}^{\ell}},$
(120)

where (cf. (8)) $\lambda^{\ell}_{i,j}=c_{i,j}n_{j}\Psi^{\ell}_{i,j}$. Then

$\alpha_{i,j}^{\ell+1}-\alpha_{i,j}^{\ell}=0=\beta_{i,j}^{\ell+1}-\beta_{i,j}^{\ell}.$
(121)

###### Proof.

Only the first equality is shown. The other case is similar. Direct
calculation gives

$\alpha_{i,j}^{\ell+1}-\alpha_{i,j}^{\ell}=\frac{\rho_{i}\rho_{j}\left(\lambda_{i,j}^{\ell+1}\lambda_{j,i}^{\ell}-\lambda_{i,j}^{\ell}\lambda_{j,i}^{\ell+1}\right)}{(\rho_{i}\lambda_{i,j}^{\ell+1}+\rho_{j}\lambda_{j,i}^{\ell+1})(\rho_{i}\lambda_{i,j}^{\ell}+\rho_{j}\lambda_{j,i}^{\ell})},$
(122)

It is easy to verify that the product
$\lambda_{i,j}^{\ell+1}\lambda_{j,i}^{\ell}$ is symmetric in $i,j$. Hence the
difference in (122) is identically zero. ∎

##  Appendix D Temperature Source Terms Are Positive

This appendix contains three technical lemmas used in the proof of (33) in
Proposition 3.8.

######  Lemma D.1.

$s_{i}^{n}-s_{i}^{\ell+1}=\left|\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}_{i}^{\mathfrak{n}}\right|^{2}+\frac{\Delta
t}{\varepsilon}\sum_{j}\lambda_{i,j}^{\ell}\alpha_{j,i}\left[2s_{i}^{\ell+1}-2\left(\boldsymbol{u}_{i}^{\ell+1},\boldsymbol{u}_{j}^{\ell+1}\right)_{2}\right].$
(123)

###### Proof.

Rearrange the iterative step from (18) to obtain

$\displaystyle\boldsymbol{u}_{i}^{n}-\boldsymbol{u}_{i}^{\ell+1}$
$\displaystyle=\frac{\Delta
t}{\varepsilon}\sum_{j}\lambda_{i,j}^{\ell}\alpha_{j,i}(\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}_{j}^{\ell+1})$
(124a)
$\displaystyle\text{and}\qquad\boldsymbol{u}_{i}^{n}+\boldsymbol{u}_{i}^{\ell+1}$
$\displaystyle=2\boldsymbol{u}_{i}^{\ell+1}+\frac{\Delta
t}{\varepsilon}\sum_{j}\lambda_{i,j}^{\ell}\alpha_{j,i}(\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}_{j}^{\ell+1}).$
(124b)

Take the inner product of these to obtain the result. ∎

######  Lemma D.2.

$\displaystyle
m_{i}\beta_{i,j}(s_{i}^{\ell+1}-S_{i,j}^{\ell+1})+m_{j}\beta_{j,i}(s_{j}^{\ell+1}-S_{i,j}^{\ell+1})=m_{i}\alpha_{j,i}\beta_{i,j}\left|\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}_{j}^{\ell+1}\right|^{2},$
(125a) $\displaystyle\frac{\Delta t}{\varepsilon
d}\sum_{j}\lambda_{i,j}^{\ell}\beta_{j,i}\left[m_{j}\left(s_{j}^{\ell+1}-S_{i,j}^{\ell+1}\right)-m_{i}\left(s_{i}^{\ell+1}-S_{i,j}^{\ell+1}\right)\right]$
(125b) $\displaystyle\qquad=\frac{\Delta t}{\varepsilon
d}\sum_{j}m_{i}\lambda_{i,j}^{\ell}\alpha_{j,i}\bigg{[}\beta_{i,j}\left|\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}_{j}^{\ell+1}\right|^{2}-(1+\alpha_{i,j})s_{i}^{\ell+1}+\alpha_{j,i}s_{j}^{\ell+1}+2\alpha_{i,j}\left(\boldsymbol{u}_{i}^{\ell+1},\boldsymbol{u}_{j}^{\ell+1}\right)_{2}\bigg{]}.$

###### Proof.

The proof of (125a) is a calculation using the fact that
$\alpha_{i,j}+\alpha_{j,i}=1=\beta_{i,j}+\beta_{j,i}$, and the identities

$m_{j}\alpha_{i,j}\beta_{j,i}=m_{i}\alpha_{j,i}\beta_{i,j}\qquad\iff\qquad\alpha_{j,i}(m_{i}\beta_{i,j}+m_{j}\beta_{j,i})=m_{j}\beta_{j,i}.$
(126)

To prove (125b), use (125a). ∎

######  Proposition D.3.

$\frac{m_{i}}{d}(s_{i}^{\mathfrak{n}}-s_{i}^{\ell+1})+\frac{\Delta
t}{\varepsilon
d}\sum_{j}\lambda_{i,j}^{\ell}\beta_{j,i}\left[m_{j}(s_{j}^{\ell+1}-S_{i,j}^{\ell+1})-m_{i}(s_{i}^{\ell+1}-S_{i,j}^{\ell+1})\right]\geq
0.$ (127)

###### Proof.

Using Lemmas D.1 and D.2, and $\alpha_{i,j}+\alpha_{j,i}=1$,

$\begin{split}\frac{m_{i}}{d}&(s_{i}^{\mathfrak{n}}-s_{i}^{\ell+1})+\frac{\Delta
t}{\varepsilon
d}\sum_{j}\lambda_{i,j}^{\ell}\beta_{j,i}\left[m_{j}(s_{j}^{\ell+1}-S_{i,j}^{\ell+1})-m_{i}(s_{i}^{\ell+1}-S_{i,j}^{\ell+1})\right]\\\
&=\frac{m_{i}}{d}\left|\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}_{i}^{\mathfrak{n}}\right|^{2}+\frac{\Delta
t}{\varepsilon}\frac{m_{i}}{d}\sum_{j}\lambda_{i,j}^{\ell}\alpha_{j,i}(\alpha_{j,i}+\beta_{i,j})\left|\boldsymbol{u}_{i}^{\ell+1}-\boldsymbol{u}_{j}^{\ell+1}\right|^{2}\geq
0.\end{split}$ (128)

∎

## References

  * [1] Pierre Andries, Kazuo Aoki, and Benoit Perthame.  A consistent BGK-type model for gas mixtures.  Journal of Statistical Physics, 106:993–1018, 2002.
  * [2] Pierre Andries, Patrick Le Tallec, Jean-Philippe Perlat, and Benoıt Perthame.  The Gaussian-BGK model of Boltzmann equation with small Prandtl number.  European Journal of Mechanics-B/Fluids, 19(6):813–830, 2000.
  * [3] Uri M Ascher, Steven J Ruuth, and Raymond J Spiteri.  Implicit-explicit Runge–Kutta methods for time-dependent partial differential equations.  Applied Numerical Mathematics, 25(2-3):151–167, 1997.
  * [4] P.L. Bhatnagar, E.P. Gross, and M. Krook.  Model for collision processes in gases I: Small amplitude processes in charged and neutral one-component systems.  Phys. Rev., 94:511, 1954.
  * [5] G. Bird.  Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Volume 1.  Clarendon Press, 1994.
  * [6] 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.
  * [7] Sebastiano Boscarino, Lorenzo Pareschi, and Giovanni Russo.  A unified IMEX Runge–Kutta approach for hyperbolic systems with multiscale relaxation.  SIAM Journal on Numerical Analysis, 55(4):2085–2109, 2017.
  * [8] Stéphane Brull, Vincent Pavan, and Jacques Schneider.  Derivation of a BGK model for mixtures.  European Journal of Mechanics-B/Fluids, 33:74–86, 2012.
  * [9] Carlo Cercignani, Reinhard Illner, and Mario Pulvirenti.  The mathematical theory of dilute gases, volume 106.  Springer Science & Business Media, 2013.
  * [10] S. Chapman and T.G. Cowling.  The Mathematical Theory of Non-uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases.  Cambridge Mathematical Library. Cambridge University Press, 1990.
  * [11] C. K. Chu.  Kinetic‐theoretic description of the formation of a shock wave.  The Physics of Fluids, 8(1):12–22, 1965.
  * [12] F. Coron and B. Perthame.  Numerical passage from kinetic to fluid equations.  SIAM J. Numer. Anal., 28(1):26–42, 1991.
  * [13] 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.
  * [14] 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.
  * [15] Jeff Haack, Cory Hauck, Christian Klingenberg, Marlies Pirner, and Sandra Warnecke.  A consistent BGK model with velocity-dependent collision frequency for gas mixtures.  Journal of Statistical Physics, 184:1–17, 2021.
  * [16] J.R. Haack, C.D. Hauck, and M.S. Murillo.  A conservative, entropic multispecies BGK model.  J. Stat Phys, 168:826–856, 2017.
  * [17] 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.
  * [18] Evan Habbershaw and Steven M Wise.  A progress report on numerical methods for BGK-type kinetic equations.  Technical report, University of Tennessee, 2022.
  * [19] Lowell H Holway Jr.  Kinetic theory of shock structure using an ellipsoidal distribution function.  Rarefied Gas Dynamics, Volume 1, 1:193, 1965.
  * [20] 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.
  * [21] Andrew Lenard and Ira B Bernstein.  Plasma oscillations with diffusion in velocity space.  Physical Review, 112(5):1456, 1958.
  * [22] Richard L Liboff.  Kinetic theory: classical, quantum, and relativistic descriptions.  Springer Science & Business Media, 2003.
  * [23] Lorenzo Pareschi and Giovanni Russo.  Implicit–explicit Runge–Kutta schemes and applications to hyperbolic systems with relaxation.  Journal of Scientific computing, 25:129–155, 2005.
  * [24] B. Perthame and M. Pulvirenti.  Weighted L∞ bounds and uniqueness for the Boltzmann BGK model.  Archive Rat. Mech. Anal., 125(3):289–295, 1993.
  * [25] Benoît Perthame.  Global existence to the BGK model of Boltzmann equation.  J. Diff. Eq., 82:191–205, 1989.
  * [26] Gabriella Puppo.  Kinetic models of BGK type and their numerical integration.  arXiv preprint arXiv:1902.08311, 2019.
  * [27] Gary A Sod.  A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws.  Journal of computational physics, 27(1):1–31, 1978.
  * [28] Eleuterio F Toro.  Riemann solvers and numerical methods for fluid dynamics: a practical introduction.  Springer Science & Business Media, 2013.
  * [29] 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.