# Numerical schemes for a multi-species quantum BGK model

Gi-Chan Bae , Marlies Pirner , Sandra Warnecke Seoul National
UniversityUniversity of MuensterUniversity of Wuerzburg

Abstract: We consider a kinetic model of an N-species gas mixture modeled with
quantum Bhatnagar-Gross-Krook (BGK) collision operators. The collision
operators consist of a relaxation to a Maxwell distribution in the classical
case, a Fermi distribution for fermions and a Bose-Einstein distribution for
bosons. In this paper we present a numerical method for simulating this model,
which uses an Implicit-Explicit (IMEX) scheme to minimize a certain potential
function. This is motivated by theoretical considerations coming from entropy
minimization. We show that theoretical properties such as conservation of
mass, total momentum and total energy as well as positivity of the
distribution functions are preserved by the numerical method presented in this
paper, and illustrate its usefulness and effectiveness with numerical examples

##  1 Introduction

In a kinetic description, the state of a dilute gas or plasma is given by a
distribution function that prescribes the density of particles at each point
in position-momentum phase space. In a time-dependent setting, the evolution
of this distribution function is due to a balance of particle advection and
binary collisions. Perhaps the most well-known model for collisions is the
Boltzmann collision operator, an integral operator that preserves collision
invariants (mass, momentum and energy) and dissipates the mathematical entropy
of the system. Unfortunately, the expense of evaluating this operator can be
prohibitive. Indeed, its evaluation requires the calculation of a five-
dimensional integral at every point in phase-space. Thus even with fast
spectral methods [35, 36, 16, 15], the collision operator is typically the
dominant part of a kinetic calculation. The quantum modification of the
celebrated Boltzmann equation was made in [39, 40] to incorporate the quantum
effect that cannot be neglected for light molecules (such as Helium) at low
temperature. Quantum Boltzmann equation is now fruitfully employed not just
for low temperature gases, but in various circumstances such as scattering
problem in solid [3, 13] and electrons on energy band structure in
semiconductor [27].

In the classical case, the Bhatnagar-Gross-Krook (BGK) operator is a widely
used surrogate for the Boltzmann operator that models collisions by a simple
relaxation mechanism. This simplification brings significant computational
advantages while also maintaining the conservation and entropy dissipation
properties of the Boltzmann operator. As in the classical case, the quantum
BGK models are widely used in place of the quantum Boltzmann equation [3, 13,
27, 31, 37, 32].

In the quantum case, an extension of the single-species BGK model to the
multi-species setting was recently developed in [7]. There a sufficient
condition is proven that guarantees the existence of equilibrium coefficients
so that the model shares the same conservation laws and H-theorem with the
quantum Boltzmann equation. Unlike the classical BGK model for gas mixtures
[20, 25, 19, 18, 38, 28, 21, 10, 1, 11], the equilibrium coefficients of the
local equilibrium for quantum multi-species gases are defined through highly
nonlinear relations that are not explicitly solvable. So it was necessary to
verify that such nonlinear relations uniquely determine the equilibrium
coefficients, leading to the well-definedness of the model.

In this paper we present a numerical implementation of the quantum multi-
species BGK model developed in [7]. The main obstacle of the quantum mixture
BGK model is that the equilibrium parameters can not be written by explicit
relation with mass, momentum, and energy. Here, we present a method that
enables an implicit treatment of the quantum BGK operator following the
recently developed method on the BGK model with velocity-dependent collision
frequency [22, 23]. The implementation is a discrete velocity method that
relies on standard spatial and temporal discretizations from the literature.
The crucial step involves the formulation of a convex entropy minimization
problem. In particular, the solver uses a numerical minimization procedure in
order to determine the coefficients of the attractors. The numerical method
was originally developed for a BGK model for gas mixtures with velocity
depedent collision frequencies [22, 23]. In this paper, we will see that this
method can be also adapted for quantum BGK models for gas mixtures since both
models have structural similarities. Both models have the difficulty that the
dependency of the attractors on the solution can not be given in an explicit
way. The implementation is a discrete velocity method that relies on standard
spatial and temporal discretizations from the literature. The key new
ingredient is a solver which enables an implicit treatment of the BGK
operator. The crucial step involves the formulation of a convex entropy
minimization problem. In particular, the solver uses a numerical minimization
procedure in order to determine the coefficients of the attractors. This
construction guarantees conservation and entropy properties at the discrete
level, up to numerical tolerances, even when using a discrete velocity mesh.

The remainder of this paper is organized as follows. In Section 2, 3 and 4, we
recall the multi-species quantum BGK model from [7] and its main important
properties. In Section 5.1, we present the first- and second-order implicit-
explicit time discretizations that are used in the paper. We also introduce
the optimization-based approach for the implicit evaluation of the BGK
operator. In Section 5.2, we describe the space discretization. In Section
5.3, we verify some structure preserving properties of the semi-discrete
scheme. In Section 5.4, we introduce the momentum discretization and summarize
the numerical implementation of the optimization algorithm introduced in
Section 5.1. In Section 6, we provide an array of numerical results that
illustrate the properties of our scheme.

##  2 A consistent multi-species quantum BGK model

We consider two distribution functions $f_{1}=f_{1}(x,p,t)\geq 0$ and
$f_{2}=f_{2}(x,p,t)\geq 0$ for species with masses $m_{1}>0$ and $m_{2}>0$,
respectively, with the phase space variables (position and momentum)
$x\in\Omega$ and $p\in\mathbb{R}^{3}$ and time $t\geq 0$. To be as general as
possible, we generalize the quantum mixture BGK model in [7] describing a
mixture of bosons and fermions to the mixture of gases including the
interaction of quantum-classical particles:

$\displaystyle\begin{split}\partial_{t}f_{1}+\frac{p}{m_{1}}\cdot\nabla_{x}f_{1}=\nu_{11}n_{1}(\mathcal{K}_{11}-f_{1})+\nu_{12}n_{2}(\mathcal{K}_{12}-f_{1}),\\\
\partial_{t}f_{2}+\frac{p}{m_{2}}\cdot\nabla_{x}f_{2}=\nu_{22}n_{2}(\mathcal{K}_{22}-f_{2})+\nu_{21}n_{1}(\mathcal{K}_{21}-f_{2}),\end{split}$
(1)

where $\mathcal{K}_{ij}$ is the local equilibrium describing the interactions
of $i$th and $j$th component and $\nu_{ij}n_{j}$ for $i,j=1,2$ are the
collision frequencies. More explicitly, there can be the following cases for
fermion $\tau=+1$, for boson $\tau=-1$, and for classical particle $\tau=0$:

$\displaystyle\mathcal{K}_{11}=\frac{1}{e^{m_{1}a_{1}\big{|}\frac{p}{m_{1}}-b_{1}\big{|}^{2}+c_{1}}+\tau},\quad\mathcal{K}_{12}=\frac{1}{e^{m_{1}a\big{|}\frac{p}{m_{1}}-b\big{|}^{2}+c_{12}}+\tau},$
(2)
$\displaystyle\mathcal{K}_{22}=\frac{1}{e^{m_{2}a_{2}\big{|}\frac{p}{m_{2}}-b_{2}\big{|}^{2}+c_{2}}+\tau^{\prime}},\quad\mathcal{K}_{21}=\frac{1}{e^{m_{2}a\big{|}\frac{p}{m_{2}}-b\big{|}^{2}+c_{21}}+\tau^{\prime}}.$
(3)

We denote $\mathcal{K}$ as Fermi-Dirac distribution, Bose-Einstein
distribution and Maxwellian for the case $\tau=+1,-1,0$, respectively. The
equilibrium parameters $(a_{i},b_{i},c_{i})$ and $(a,b,c_{12},c_{21})$ will be
determined later to satisfy the conservation laws and the entropy principle.
Note that the model includes the following cases depending on the types of the
particles

$\displaystyle(\tau,\tau^{\prime})=\left\\{\begin{array}[]{ll}(+1,+1)&\text{(fermion-
fermion)}\cr(-1,-1)&\mbox{(boson-boson)}\cr(+1,-1)&\mbox{(fermion-
boson)}\cr(+1,~{}0~{})&\mbox{(fermion-
classical)}\cr(~{}0,-1~{})&\mbox{(classical-
boson)}\cr(~{}0~{},~{}0~{})&\mbox{(classical-classical)}\end{array}\right.$

We define the number density of particles $n_{i}$, momentum $P_{i}$, energy
$E_{i}$ of each species as

$\displaystyle n_{i}=\int_{\mathbb{R}^{3}}f_{i}dp,\quad
P_{i}=\int_{\mathbb{R}^{3}}f_{i}p~{}dp,\quad
E_{i}=\int_{\mathbb{R}^{3}}f_{i}\frac{|p|^{2}}{2m_{i}}dp,$ (4)

and $m_{i}n_{i}=N_{i}$. The parameters $(a_{i},b_{i},c_{i})$ in
$\mathcal{K}_{ii}$ are chosen such that the local equilibrium satisfies the
following identities:

$\displaystyle\int_{\mathbb{R}^{3}}\mathcal{K}_{ii}dp=n_{i},\quad\int_{\mathbb{R}^{3}}\mathcal{K}_{ii}p~{}dp=P_{i},\quad\int_{\mathbb{R}^{3}}\mathcal{K}_{ii}\frac{|p|^{2}}{2m_{i}}dp=E_{i},\quad(i=1,2).$
(5)

This ensures conservation of the number of particles, momentum and energy in
interactions of a species with itself.

The equilibrium parameters $(a,b,c_{12},c_{21})$ of $\mathcal{K}_{ij}$ for
$(ij)=(12),(21)$ are determined to satisfy the following identities:

$\displaystyle\begin{split}&\int_{\mathbb{R}^{3}}\mathcal{K}_{12}dp=n_{1},\qquad\int_{\mathbb{R}^{3}}\mathcal{K}_{21}dp=n_{2},\cr&\nu_{12}n_{2}\left(\int_{\mathbb{R}^{3}}\mathcal{K}_{12}p~{}dp-
P_{1}\right)+\nu_{21}n_{1}\left(\int_{\mathbb{R}^{3}}\mathcal{K}_{21}p~{}dp-
P_{2}\right)=0,\cr&\nu_{12}n_{2}\left(\int_{\mathbb{R}^{3}}\mathcal{K}_{12}\frac{|p|^{2}}{2m_{1}}dp-
E_{1}\right)+\nu_{21}n_{1}\left(\int_{\mathbb{R}^{3}}\mathcal{K}_{21}\frac{|p|^{2}}{2m_{2}}dp-
E_{2}\right)=0.\end{split}$ (6)

This ensures conservation of the number of particles, momentum and energy in
interactions of a species with the other one.

The existence of these parameters $(a_{i},b_{i},c_{i})$ and
$(a,b,c_{12},c_{21})$ is proven in [7] for quantum-quantum mixture case with
unit collision frequencies. We will show the proof for any choice of mixture
between classical particle, fermion, and boson for more general collision
frequency the correspondence with the classical case.

######  Remark 1.

For the correspondence of classical case and quantum case, let us assume that
the velocity distribution function $\bar{f}(x,v,t)$ and the momentum
distribution function $f(x,p,t)$ satisfy the following relation for $(i=1,2)$:

$\displaystyle\bar{f}_{i}(x,v,t)=\bar{f}_{i}\Big{(}x,\frac{p}{m_{i}},t\Big{)}=m_{i}^{3}f_{i}(x,p,t).$

This relation connects the conservation laws of classical case and quantum
case as follows:

$\displaystyle\int_{\mathbb{R}^{3}}\bar{f}_{i}(x,v,t)\left(\begin{array}[]{c}1\cr
m_{i}v\cr\frac{m_{i}}{2}|v|^{2}\end{array}\right)dv$
$\displaystyle=\int_{\mathbb{R}^{3}}\frac{1}{m^{3}_{i}}\bar{f}_{i}\left(x,\frac{p}{m_{i}},t\right)\left(\begin{array}[]{c}1\cr
p\cr\frac{|p|^{2}}{2m_{i}}\end{array}\right)dp=\int_{\mathbb{R}^{3}}f_{i}(x,p,t)\left(\begin{array}[]{c}1\cr
p\cr\frac{|p|^{2}}{2m_{i}}\end{array}\right)dp.$

Thus the macroscopic fields of the classical one

$\displaystyle n_{i}(x,t)$
$\displaystyle=\int_{\mathbb{R}^{3}}\bar{f}_{i}(x,v,t)dv,\quad
m_{i}n_{i}=\rho_{i}$ $\displaystyle U_{i}(x,t)$
$\displaystyle=\frac{1}{n_{i}}\int_{\mathbb{R}^{3}}\bar{f}_{i}(x,v,t)vdv,$
$\displaystyle T_{i}(x,t)$
$\displaystyle=\frac{1}{3n_{i}}\int_{\mathbb{R}^{3}}\bar{f}_{i}(x,v,t)m_{i}|v-U_{i}|^{2}dv,$

and the quantum one

$\displaystyle n_{i}(x,t)$
$\displaystyle=\int_{\mathbb{R}^{3}}f_{i}(x,p,t)dp,\quad m_{i}n_{i}=N_{i}$
$\displaystyle P_{i}(x,t)$
$\displaystyle=\int_{\mathbb{R}^{3}}f_{i}(x,p,t)pdp,$ $\displaystyle
E_{i}(x,t)$
$\displaystyle=\int_{\mathbb{R}^{3}}f_{i}(x,p,t)\frac{|p|^{2}}{2m_{i}}dp,$

have the following relation:

$\displaystyle\begin{split}\rho_{i}=&m_{i}\int_{\mathbb{R}^{3}}\bar{f}dv=m_{i}\int_{\mathbb{R}^{3}}fdp=N_{i},\cr\rho_{i}U_{i}=&m_{i}\int_{\mathbb{R}^{3}}\bar{f}vdv=\int_{\mathbb{R}^{3}}fpdp=P_{i},\cr\frac{3}{2}n_{i}T_{i}+\frac{1}{2}\rho_{i}|U_{i}|^{2}&=\int_{\mathbb{R}^{3}}\bar{f}\frac{m_{i}}{2}|v|^{2}dv=\int_{\mathbb{R}^{3}}f\frac{|p|^{2}}{2m_{i}}dp=E_{i}.\end{split}$
(7)

Having introduced our notation, we can deep dive into the model [7].

##  3 Properties of the model

In this section, we present the main properties of the quantum multi-species
BGK model (1).

###  3.1 Conservation properties, H-Theorem and Boundedness

First of all, the model satisfies the conservation properties and an
$H$-Theorem if the macroscopic fields have some boundedness conditions. This
is already proved in [7] for the Fermion-Fermion, Fermion-Boson and Boson-
Boson case with collision frequency $\nu_{ij}n_{j}=1,i,j=1,2$. We generalize
the result to include the Fermion-classical, Boson-classical and classical-
classical case with $(x,t)$ depending general collision frequency
$\nu_{ij}n_{j}$.

One of the most characterizing features of quantum mechanics is that the
number of particles cannot exceed the specific ratio of the energy. If the
number of particles exceeds that critical point, then Bose-Einstein
condensation occurs for bosons or saturated state occurs for fermions.
Consequently, the local equilibrium of the one-species quantum BGK model is
only defined below the specific ratio between number of particles and energy
(see [12, 4, 29, 30]). Hence we introduce the function $j_{\tau}$ by

$\displaystyle
j_{\tau}(x)=\frac{\int\frac{1}{e^{|p|^{2}+x}+\tau}dp}{\left(\int\frac{|p|^{2}}{e^{|p|^{2}+x}+\tau}dp\right)^{3/5}}.$

In the following, we will describe how the equilibrium parameters
$(a_{i},b_{i},c_{i})$ and $(c_{12},c_{21},a,b)$ in (3) will be determined in
order to ensure the conservation properties (5) and (6). We start withe the
equilibrium parameters in the one-species term. The equilibrium parameter
$c_{i}$ is uniquely determined by the following implicit relation,

$\displaystyle
j_{\tau}(c_{i})=\frac{n_{i}}{\left({\color[rgb]{0,0,0}2m_{i}E_{i}}-{|P_{i}|^{2}}/{n_{i}}\right)^{\frac{3}{5}}},$

if the right-hand-side is bounded as follows:

$\displaystyle\frac{n_{i}}{\left({\color[rgb]{0,0,0}2m_{i}E_{i}}-{|P_{i}|^{2}}/{n_{i}}\right)^{\frac{3}{5}}}\leq\begin{cases}j_{+1}(-\infty)\qquad&\mbox{for
fermion},\cr j_{-1}(0)\qquad&\mbox{for boson}.\end{cases}$

(with notation
$j_{+1}(-\infty)=\displaystyle{\lim_{x\rightarrow-\infty}j_{+1}(x)}$. The
reason is that the functions $j_{+1}$ and $j_{-1}$ are strictly decreasing on
$(-\infty,\infty)$ and $[0,\infty)$, respectively).

In the classical case, we can explicitly compute
$j_{0}(x)=\left(\frac{2\pi}{3}\right)^{\frac{3}{5}}e^{-\frac{2}{5}x}$, which
is strictly decreasing. Thus the maximum value of $j_{0}(x)$ is equal to
$\lim_{x\rightarrow-\infty}j_{0}(x)=\infty$, and the macroscopic fields of the
classical case has no restriction from above:

$\displaystyle
0\leq\frac{n_{i}}{\left({\color[rgb]{0,0,0}2m_{i}E_{i}}-{|P_{i}|^{2}}/{n_{i}}\right)^{\frac{3}{5}}}<\infty.$

We define $l_{+1},l_{0}$ and $l_{-1}$ to be such points where the respective
functions $j_{+1},j_{0}$ and $j_{-1}$ are maximal and use the following
notation

$\displaystyle l:\\{+1,-1,0\\}\rightarrow[-\infty,\infty],\quad
l_{\tau}=\begin{cases}l_{+1}=-\infty,\cr l_{-1}=0,\cr
l_{0}=-\infty.\end{cases}$

Then, if the specific ratio of macroscopic fields has the following bound

$\displaystyle\frac{n_{i}}{\left({\color[rgb]{0,0,0}2m_{i}E_{i}}-{|P_{i}|^{2}}/{n_{i}}\right)^{\frac{3}{5}}}\leq
j_{\tau}(l_{\tau}),$

then the equilibrium parameter $c_{i}$ is uniquely determined by

$\displaystyle
c_{i}=j_{\tau}^{-1}\left(\frac{n_{i}}{\left({\color[rgb]{0,0,0}2m_{i}E_{i}}-{|P_{i}|^{2}}/{n_{i}}\right)^{\frac{3}{5}}}\right).$

Finally, we can define $a_{i}$, and $b_{i}$ by

$\displaystyle
a_{i}=\left({\int_{\mathbb{R}^{3}}\frac{1}{e^{|p|^{2}+c_{i}}+\tau}dp}\right)^{\frac{2}{3}}n_{i}^{-\frac{2}{3}},\qquad
b_{i}=\frac{P_{i}}{n_{i}},$ (8)

for $(i=1,2).$ This choice ensures the conservation properties (5) and is
proven in [4]. To summarize, we have the following theorem.

######  Theorem 3.1.1.

Let the macroscopic fields satisfy

$\displaystyle\frac{n_{1}}{\left({\color[rgb]{0,0,0}2m_{1}E_{1}}-{|P_{1}|^{2}}/{n_{1}}\right)^{\frac{3}{5}}}\leq
j_{\tau}(l(\tau)),\qquad\frac{n_{2}}{\left({\color[rgb]{0,0,0}2m_{2}E_{2}}-{|P_{2}|^{2}}/{n_{2}}\right)^{\frac{3}{5}}}\leq
j_{\tau^{\prime}}(l(\tau^{\prime})).$

Then the parameters of $\mathcal{K}_{11}$ and $\mathcal{K}_{22}$ are uniquely
determined as

$\displaystyle
c_{1}=j_{\tau}^{-1}\left(\frac{n_{1}}{\left({\color[rgb]{0,0,0}2m_{1}E_{1}}-{|P_{1}|^{2}}/{n_{1}}\right)^{\frac{3}{5}}}\right),\qquad
c_{2}=j_{\tau^{\prime}}^{-1}\left(\frac{n_{2}}{\left({\color[rgb]{0,0,0}2m_{2}E_{2}}-{|P_{2}|^{2}}/{n_{1}}\right)^{\frac{3}{5}}}\right).$

and

$\displaystyle
a_{1}=\left({\int_{\mathbb{R}^{3}}\frac{1}{e^{|p|^{2}+c_{1}}+\tau}dp}\right)^{\frac{2}{3}}n_{1}^{-\frac{2}{3}},\qquad
a_{2}=\left({\int_{\mathbb{R}^{3}}\frac{1}{e^{|p|^{2}+c_{2}}+\tau^{\prime}}dp}\right)^{\frac{2}{3}}n_{2}^{-\frac{2}{3}},$

and

$\displaystyle b_{1}=\frac{P_{1}}{n_{1}},\quad b_{2}=\frac{P_{2}}{n_{2}}.$

Then with this choice of the equilibrium parameters of the one species quantum
equilibrium $\mathcal{K}_{11}$ and $\mathcal{K}_{22}$ the conservation
properties (5) are satisfied.

We note that the local Maxwellian for the classical case corresponds to the
original definition in [9]:

$\displaystyle\mathcal{M}_{ii}=\frac{N_{i}}{\pi^{\frac{3}{2}}}\left(\frac{3}{2}\frac{1}{\frac{E_{i}}{N_{i}}-\big{|}\frac{P_{i}}{N_{i}}\big{|}^{2}}\right)^{\frac{3}{2}}e^{-\frac{3}{2}\frac{\big{|}p-\frac{P_{i}}{N_{i}}\big{|}^{2}}{\frac{E_{i}}{N_{i}}-\big{|}\frac{P_{i}}{N_{i}}\big{|}^{2}}}=\frac{n_{i}}{\sqrt{2\pi\frac{T_{i}}{m_{i}}}^{3}}\exp\left(-\frac{|v-U_{i}|^{2}}{2\frac{T_{i}}{m_{i}}}\right),$

where we applied the correspondence of the macroscopic fields (7). Then we can
see that with such choice of $a_{i}$, $b_{i}$ and $c_{i}$, the local
equilibria $\mathcal{K}_{11}=\mathcal{M}_{11}$ and
$\mathcal{K}_{22}=\mathcal{M}_{22}$ satisfy (5).
Before stating the theorem for the mixture interaction terms, we introduce
some notations. We define

$\displaystyle\eta_{\tau}(x)=\int_{\mathbb{R}^{3}}\frac{1}{e^{|p|^{2}+x}+\tau}dp.$

Note that $\eta_{\tau}^{-1}$ always exist since $\eta_{\tau}$ is strictly
decreasing. We define $k_{\tau,\tau^{\prime}}$ by

$\displaystyle
k_{\tau,\tau^{\prime}}(x,y)=\frac{m_{1}^{3/2}\int_{\mathbb{R}^{3}}\frac{1}{e^{|p|^{2}+x}+\tau}dp}{\left(\frac{\nu_{12}}{m_{2}}\left(\frac{m_{1}^{3/2}}{2N_{1}}\int_{\mathbb{R}^{3}}\frac{|p|^{2}}{e^{|p|^{2}+x}+\tau}dp\right)+\frac{\nu_{21}}{m_{1}}\left(\frac{m_{2}^{3/2}}{2N_{2}}\int_{\mathbb{R}^{3}}\frac{|p|^{2}}{e^{|p|^{2}+y}+\tau^{\prime}}dp\right)\right)^{\frac{3}{5}}}.$

Using $\eta_{\tau}$ and $k_{\tau}$, we define $g_{\tau,\tau^{\prime}}$, which
is defined as a composite function of $k_{\tau}$ and $\eta^{-1}_{\tau}$, as
follows:

$\displaystyle
g_{\tau,\tau^{\prime}}(x)=k_{\tau,\tau^{\prime}}\big{(}x,y(x)\big{)}=\frac{m_{1}^{3/2}\int_{\mathbb{R}^{3}}\frac{1}{e^{|p|^{2}+x}+\tau}dp}{\left(\frac{\nu_{12}}{m_{2}}\left(\frac{m_{1}^{3/2}}{2N_{1}}\int_{\mathbb{R}^{3}}\frac{|p|^{2}}{e^{|p|^{2}+x}+\tau}dp\right)+\frac{\nu_{21}}{m_{1}}\left(\frac{m_{2}^{3/2}}{2N_{2}}\int_{\mathbb{R}^{3}}\frac{|p|^{2}}{e^{|p|^{2}+y(x)}+\tau^{\prime}}dp\right)\right)^{\frac{3}{5}}},$
(9)

where $y(x)$ denotes

$\displaystyle
y(x)=\eta_{\tau^{\prime}}^{-1}\left(\frac{m_{1}^{3/2}n_{2}}{m_{2}^{3/2}n_{1}}\eta_{\tau}(x)\right).$

######  Theorem 3.1.2.

Assume $\nu_{12}n_{2}=\nu_{21}n_{1}$. Assume further that

$\displaystyle\frac{n_{1}}{\left({\color[rgb]{0,0,0}2E_{1}}+{\color[rgb]{0,0,0}2E_{2}}-\frac{|P_{1}+P_{2}|^{2}}{{\color[rgb]{0,0,0}N_{1}+N_{2}}}\right)^{\frac{3}{5}}}\leq{\color[rgb]{0,0,0}g_{\tau,\tau^{\prime}}}\left(\max\left\\{l(\tau),h_{\tau}^{-1}\left(\frac{{\color[rgb]{0,0,0}m_{2}^{\frac{3}{2}}}n_{1}}{{\color[rgb]{0,0,0}m_{1}^{\frac{3}{2}}}n_{2}}h_{\tau^{\prime}}(l(\tau^{\prime}))\right)\right\\}\right).$

Then $c_{12}$, $c_{21}$ are defined as a unique solution of the following
relations:

$\displaystyle\frac{{\color[rgb]{0,0,0}m_{1}^{\frac{3}{2}}}h_{\tau}(c_{12})}{{\color[rgb]{0,0,0}m_{2}^{\frac{3}{2}}}h_{\tau^{\prime}}(c_{21})}=\frac{n_{1}}{n_{2}},\quad
k_{\tau,\tau^{\prime}}(c_{12},c_{21})=\frac{n_{1}}{\left({\color[rgb]{0,0,0}2E_{1}}+{\color[rgb]{0,0,0}2E_{2}}-\frac{|P_{1}+P_{2}|^{2}}{{\color[rgb]{0,0,0}N_{1}+N_{2}}}\right)^{\frac{3}{5}}}.$
(10)

With such $c_{12}$ and $c_{21}$, we define $a$ and $b$ by

$\displaystyle
a=\left(\frac{m_{1}^{\frac{3}{2}}\int_{\mathbb{R}^{3}}\frac{|p|^{2}}{e^{|p|^{2}+c_{12}}+\tau}dp+m_{2}^{\frac{3}{2}}\int_{\mathbb{R}^{3}}\frac{|p|^{2}}{e^{|p|^{2}+c_{21}}+\tau^{\prime}}dp}{2E_{1}+2E_{2}-\frac{|P_{1}+P_{2}|^{2}}{N_{1}+N_{2}}}\right)^{\frac{2}{5}},\quad
b=\frac{P_{1}+P_{2}}{N_{1}+N_{2}},$ (11)

Then with this choice of the equilibrium parameters of the mixture species
quantum equilibrium $\mathcal{K}_{12}$ and $\mathcal{K}_{21}$ the conservation
properties (6) are satisfied.

Another property of the model is that the distribution function in the fermion
case remains bounded from above by $1$ for all times $t\geq 0$. The proof can
be found in [7].

######  Lemma 3.1.3.

Let $f_{i}$ be a distribution function for fermions and $f_{i}(x,p,0)<1$. Then
we have $f_{i}(x,p,t)<1$ for $t\geq 0$.

Now we consider the entropy principle. For the convenience of notation, we
denote

$\displaystyle H_{\tau}(f)$
$\displaystyle=\begin{cases}\int_{\Omega}\int_{\mathbb{R}^{3}}f\ln
fdpdx\quad\mbox{for}\quad\tau=0\\\ \int_{\Omega}\int_{\mathbb{R}^{3}}f\ln
f+\tau^{-1}(1-\tau f)\ln(1-\tau f)dpdx\quad\mbox{for}\quad\tau=\pm
1\end{cases}$

($\tau=+1$ for fermion, $\tau=-1$, for boson) By abusing the notation
$\lim_{\tau\rightarrow 0^{+}}\tau^{-1}(1-\tau f)\ln(1-\tau f)=0$, we also
denote the integrand by $h_{\tau}$, namely

$\displaystyle h_{\tau}(z)=z\ln z+\tau^{-1}(1-\tau z)\ln(1-\tau z)$ (12)

for $z>0$ if $\tau=0,-1$ and $0<z<1$ if $\tau=+1$. Then, we define the entropy
for the gas mixture as follows:

$\displaystyle H_{\tau,\tau^{\prime}}(f_{1},f_{2})$
$\displaystyle=H_{\tau}(f_{1})+H_{\tau^{\prime}}(f_{2})$

Th function $H_{\tau,\tau^{\prime}}$ can be proven to be a non-increasing
function for the model (1).

######  Theorem 3.1.4.

Let $(f_{1},f_{2})$ be a solution of the equation (1). If the $i$-th component
is a fermion, we give an additional assumption that $f_{i}<1$. Then we have

$\displaystyle\frac{d}{dt}H_{\tau,\tau^{\prime}}(f_{1},f_{2})$
$\displaystyle\leq 0.$

The equality is characterized by $f_{1}$ and $f_{2}$ being two Fermion
distributions in the Fermion-Fermion case, two Boson distributions in the
Boson-Boson case, two Maxwell distributions in the classical-classical case, a
Fermi distribution and a Bose distribution in the Fermion-Boson case, a Fermi
distribution and a Maxwell distribution in the Fermion-classical case and a
Bose distribution and a Maxwell distribution in the Bose- classical case. In
all cases, these equilibrium distributions have the same $a$ and $b$.

The proof is a straightforward extension of the proof of Theorem 2.1(3) in
[7].

###  3.2 Macroscopic equations

From the quantum BGK model (1), one can derive the macroscopic equations for a
mixture. We first denote the macroscopic fields of the inter-species local
equilibrium $\mathcal{K}_{12}$ and $\mathcal{K}_{21}$:

$\displaystyle\begin{split}&P_{12}=\int_{\mathbb{R}^{3}}\mathcal{K}_{12}~{}p~{}dp,\quad
P_{21}=\int_{\mathbb{R}^{3}}\mathcal{K}_{21}~{}p~{}dp,\cr&E_{12}=\int_{\mathbb{R}^{3}}\mathcal{K}_{12}\frac{|p|^{2}}{2m_{1}}dp,\quad
E_{21}=\int_{\mathbb{R}^{3}}\mathcal{K}_{12}\frac{|p|^{2}}{2m_{2}}dp.\end{split}$
(13)

######  Theorem 3.2.1.

Assume $\nu_{12}n_{1}=\nu_{21}n_{2}$. Let $(f_{1},f_{2})$ be a solution to
(1), then we obtain the following formal conservation laws

$\displaystyle\begin{split}&\partial_{t}n_{1}+\frac{1}{m_{1}}\nabla_{x}\cdot
P_{1}=0,\quad\partial_{t}n_{2}+\frac{1}{m_{2}}\nabla_{x}\cdot P_{2}=0,\\\
&\partial_{t}P_{1}+\frac{1}{m_{1}}\nabla_{x}\cdot\int p\otimes
pf_{1}dp=\nu_{12}n_{2}(P_{12}-P_{1})\\\
&\partial_{t}P_{2}+\frac{1}{m_{2}}\nabla_{x}\cdot\int p\otimes
pf_{1}dp=\nu_{21}n_{1}(P_{21}-P_{2})\\\
&\partial_{t}E_{1}+\frac{1}{2m_{1}^{2}}\nabla_{x}\cdot\int|p|^{2}pf_{1}dp=\nu_{12}n_{2}(E_{12}-E_{1}),\cr&\partial_{t}E_{2}+\frac{1}{2m_{2}^{2}}\nabla_{x}\cdot\int|p|^{2}pf_{2}dp=\nu_{21}n_{1}(E_{21}-E_{2}),\end{split}$
(14)

where the exchange terms of momentum can be computed as

$\displaystyle
P_{12}-P_{1}=-(P_{21}-P_{2})=\frac{N_{1}N_{2}}{N_{1}+N_{2}}\left(\frac{P_{2}}{N_{2}}-\frac{P_{1}}{N_{1}}\right)$
(15)

Furthermore, we define the function
$\eta^{E}_{\tau}(c)=\int\frac{|p|^{2}}{e^{|p|^{2}+c}+\tau}dp,$ and obtain for
the exchange of energy

$\displaystyle\begin{split}&E_{12}-E_{1}=-(E_{21}-E_{2})\\\
&=\frac{1}{2}\frac{N_{1}|P_{1}+P_{2}|^{2}}{(N_{1}+N_{2})^{2}}+\frac{(E_{1}+E_{2})-\frac{1}{2}\frac{|P_{1}+P_{2}|^{2}}{N_{1}+N_{2}}}{m_{1}^{3/2}\eta^{E}_{\tau}(c_{12})+m_{2}^{3/2}\eta^{E}_{\tau}(c_{21})}m_{1}^{3/2}\eta^{E}_{\tau}(c_{12})-E_{1}\end{split}$
(16)

###### Proof.

We multiply the first equation of (1) by $(1,p,\frac{|p|^{2}}{2m_{1}})$, and
the second one by $(1,p,\frac{|p|^{2}}{2m_{2}}).$ Then we integrate them with
respect to the momentum $p$ to obtain (LABEL:macrosequ1) after a straight-
forward computation on the left-hand side.

The exchange of momentum can be computed as follows. After computing the
integral in the definition of $P_{12}$ and $P_{21}$ in (13), we observe that
$P_{12}=bN_{1}$ and $P_{21}=bN_{2}$. Substituting the quantity of $b$ in (11),
we obtain

$\displaystyle\frac{P_{12}}{N_{1}}-\frac{P_{1}}{N_{1}}$
$\displaystyle=\frac{P_{1}+P_{2}}{N_{1}+N_{2}}-\frac{P_{1}}{N_{1}}=\frac{N_{2}}{N_{1}+N_{2}}\left(\frac{P_{2}}{N_{2}}-\frac{P_{1}}{N_{1}}\right),$

and

$\displaystyle\frac{P_{21}}{N_{2}}-\frac{P_{2}}{N_{2}}$
$\displaystyle=\frac{P_{1}+P_{2}}{N_{1}+N_{2}}=\frac{N_{1}}{N_{1}+N_{2}}\left(\frac{P_{1}}{N_{1}}-\frac{P_{2}}{N_{2}}\right).$

Then, we can compute (similar as in section 2 in [4])

$\displaystyle
E_{12}-\frac{1}{2}\frac{|P_{12}|^{2}}{N_{1}}=\frac{1}{2}a^{-5/2}m_{1}^{3/2}\eta^{E}_{\tau}(c_{12})$
(17)

We can replace $a^{-5/2}$ with the formula (11):

$\displaystyle
E_{12}-\frac{1}{2}\frac{|P_{12}|^{2}}{N_{1}}=\frac{(E_{1}+E_{2})-\frac{1}{2}\frac{|P_{1}+P_{2}|^{2}}{N_{1}+N_{2}}}{m_{1}^{3/2}\eta^{E}_{\tau}(c_{12})+m_{2}^{3/2}\eta^{E}_{\tau^{\prime}}(c_{21})}m_{1}^{3/2}\eta^{E}_{\tau}(c_{12})$

So we obtain

$\displaystyle\begin{split}&E_{12}-E_{1}=-(E_{21}-E_{2})\\\
&=\frac{1}{2}{\color[rgb]{0,0,0}\frac{N_{1}|P_{1}+P_{2}|^{2}}{(N_{1}+N_{2})^{2}}}+\frac{(E_{1}+E_{2})-\frac{1}{2}\frac{|P_{1}+P_{2}|^{2}}{N_{1}+N_{2}}}{m_{1}^{3/2}\eta^{E}_{\tau}(c_{12})+m_{2}^{3/2}\eta^{E}_{\tau^{\prime}}(c_{21})}m_{1}^{3/2}\eta^{E}_{\tau}(c_{12})-E_{1}\end{split}$
(18)

∎

######  Remark 2.

For later purposes, we remark that in the space-homogeneous case the system of
equations reduces to

$\displaystyle\begin{split}\partial_{t}n_{1}&=0,\quad\hskip
71.13188pt\partial_{t}n_{2}=0,\\\
\partial_{t}P_{1}&=\nu_{12}n_{2}(P_{12}-P_{1}),\quad\partial_{t}P_{2}=\nu_{21}n_{1}(P_{21}-P_{2}),\\\
\partial_{t}E_{1}&=\nu_{12}n_{2}(E_{12}-E_{1}),\quad\partial_{t}E_{2}=\nu_{21}n_{1}(E_{21}-E_{2}).\end{split}$
(19)

######  Remark 3.

In the classical case $(\tau=\tau^{\prime}=0)$ we get from the relationship
(10)

$\displaystyle\frac{n_{1}}{n_{2}}=\frac{m_{1}^{3/2}\eta_{0}(c_{12})}{m_{2}^{3/2}\eta_{0}(c_{21})}=\frac{m_{1}^{3/2}\int_{\mathbb{R}^{3}}\frac{1}{e^{|p|^{2}+c_{12}}}dp}{m_{2}^{3/2}\int_{\mathbb{R}^{3}}\frac{1}{e^{|p|^{2}+c_{21}}}dp}=\frac{m_{1}^{3/2}}{m_{2}^{3/2}}\frac{e^{-c_{12}}}{e^{-c_{21}}}$

by computing the integrals explicitly. Using this, we can calculate

$\displaystyle\frac{m_{1}^{3/2}\eta_{0}^{E}(c_{12})}{m_{1}^{3/2}\eta_{0}^{E}(c_{12})+m_{2}^{3/2}\eta^{E}_{0}(c_{21})}=\frac{m_{1}^{3/2}e^{-c_{12}}}{m_{1}^{3/2}e^{-c_{12}}+m_{2}^{3/2}e^{-c_{21}}}=\frac{n_{1}}{n_{1}+n_{2}}$

and obtain

$\displaystyle E_{12}-E_{1}$
$\displaystyle=\frac{n_{1}n_{2}}{n_{1}+n_{2}}\left(\frac{E_{2}}{n_{2}}-\frac{E_{1}}{n_{1}}+\frac{m_{1}-m_{2}}{(N_{1}+N_{2})^{2}}\frac{1}{2}|P_{1}+P_{2}|^{2}\right)$
$\displaystyle=\frac{n_{1}n_{2}}{n_{1}+n_{2}}\Bigg{(}\frac{E_{2}}{n_{2}}-\frac{1}{2}\frac{|P_{2}|^{2}}{n_{2}N_{2}}-\frac{E_{1}}{n_{1}}+\frac{|P_{1}|^{2}}{n_{1}N_{1}}+m_{1}m_{2}\frac{n_{1}N_{1}+2n_{1}N_{2}+n_{2}N_{2}}{(N_{1}+N_{2})^{2}}\frac{1}{2}\frac{|P_{2}|^{2}}{N_{2}^{2}}$
$\displaystyle-
m_{1}m_{2}\frac{n_{1}N_{1}+2N_{1}n_{2}+n_{2}N_{2}}{(N_{1}+N_{2})^{2}}\frac{1}{2}\frac{|P_{1}|^{2}}{N_{1}^{2}}+m_{1}m_{2}\frac{(m_{1}-m_{2})n_{1}n_{2}}{(N_{1}+N_{2})^{2}}\frac{P_{1}}{N_{1}}\cdot\frac{P_{2}}{N_{2}}\Bigg{)}$
$\displaystyle=\frac{n_{1}n_{2}}{n_{1}+n_{2}}\Bigg{(}\frac{E_{2}}{n_{2}}-\frac{1}{2}\frac{|P_{2}|^{2}}{n_{2}N_{2}}-\frac{E_{1}}{n_{1}}+\frac{|P_{1}|^{2}}{n_{1}N_{1}}+m_{1}m_{2}\frac{n_{1}N_{1}+n_{2}N_{2}}{(N_{1}+N_{2})^{2}}\frac{1}{2}\left(\frac{|P_{2}|^{2}}{N_{2}^{2}}-\frac{|P_{1}|^{2}}{N_{1}^{2}}\right)$
$\displaystyle+\left(\frac{P_{2}}{N_{2}}-\frac{P_{1}}{N_{1}}\right)\cdot\left(\frac{P_{1}}{n_{1}}+\frac{P_{2}}{n_{2}}\right)\Bigg{)}.$

###  3.3 Decay rate for the velocities and kinetic temperatures in the space-
homogeneous case

In equilibrium, both distribution functions $f_{1}$ and $f_{2}$ will finally
share the same velocity and the same temperature. We can make those decay
rates explicit in the space-homogeneous case.

######  Theorem 3.3.1.

Assume $\nu_{12}n_{1}=\nu_{21}n_{2}$ and $\nu_{12},\nu_{21}$ independent of
$t$. In the space-homogeneous case of (1), we have the following convergence
rate for the momentum

$\frac{P_{1}}{N_{1}}-\frac{P_{2}}{N_{2}}=e^{-\frac{{\color[rgb]{0,0,0}\nu_{12}n_{2}N_{2}+\nu_{21}n_{1}N_{1}}}{N_{1}+N_{2}}t}\left(\frac{P_{1}(0)}{N_{1}}-\frac{P_{2}(0)}{N_{2}}\right).$
(20)

###### Proof.

We start with calculating

$\displaystyle\partial_{t}\left(\frac{P_{1}}{N_{1}}\right)=\frac{1}{N_{1}}\partial_{t}P_{1}={\color[rgb]{0,0,0}\nu_{12}n_{2}}\frac{1}{N_{1}}(P_{12}-P_{1}).$

The first equality uses that $N_{1}=m_{1}n_{1}$ is constant in the space-
homogeneous case (19). Then, we used the macroscopic equation for the time
evolution of $P_{1}$ from (19) Using the expression (15) for the exchange of
momentum leads to

$\displaystyle\partial_{t}\left(\frac{P_{1}}{N_{1}}\right)={\color[rgb]{0,0,0}\nu_{12}n_{2}}\frac{N_{2}}{N_{1}+N_{2}}\left(\frac{P_{2}}{N_{2}}-\frac{P_{1}}{N_{1}}\right).$
(21)

In a similar way, we can compute

$\displaystyle\partial_{t}\left(\frac{P_{2}}{N_{2}}\right)={\color[rgb]{0,0,0}\nu_{21}n_{1}}\frac{N_{1}}{N_{1}+N_{2}}\left(\frac{P_{1}}{N_{1}}-\frac{P_{2}}{N_{2}}\right).$

If we subtract the two equations we obtain

$\displaystyle\partial_{t}\left(\frac{P_{1}}{N_{1}}-\frac{P_{2}}{N_{2}}\right)=-\frac{{\color[rgb]{0,0,0}\nu_{12}n_{2}N_{2}+\nu_{21}n_{1}N_{1}}}{N_{1}+N_{2}}\left(\frac{P_{1}}{N_{1}}-\frac{P_{2}}{N_{2}}\right)$

It follows the result

$\displaystyle\frac{P_{1}}{N_{1}}-\frac{P_{2}}{N_{2}}=e^{-\frac{{\color[rgb]{0,0,0}\nu_{12}n_{2}N_{2}+\nu_{21}n_{1}N_{1}}}{N_{1}+N_{2}}t}\left(\frac{P_{1}(0)}{N_{1}}-\frac{P_{2}(0)}{N_{2}}\right).$

∎

######  Remark 4.

Using the relationship $P_{j}=N_{j}b_{j}$ from (8), one can equivalently write

$\displaystyle
b_{1}-b_{2}=e^{-\frac{\nu_{12}n_{1}N_{2}+\nu_{21}n_{2}N_{1}}{N_{1}+N_{2}}t}(b_{1}(0)-b_{2}(0)).$

We continue with the convergence rates of the quantities
$\frac{E_{1}}{n_{1}}-\frac{1}{2}\frac{|P_{1}|^{2}}{n_{1}N_{1}}$ and
$\frac{E_{2}}{n_{2}}-\frac{1}{2}\frac{|P_{2}|^{2}}{n_{2}N_{2}}$. In the
classical case, this quantity corresponds to the temperature.

######  Theorem 3.3.2.

Let $\nu_{12}n_{1}=\nu_{21}n_{2}=:\tilde{\nu}$ and $\tilde{\nu}$ be
independent of $t$. In the space-homogeneous case of (1), it is

$\displaystyle\begin{split}&\left(\frac{E_{1}}{n_{1}}-\frac{1}{2}\frac{|P_{1}|^{2}}{n_{1}N_{1}}\right)-\left(\frac{E_{2}}{n_{2}}-\frac{1}{2}\frac{|P_{2}|^{2}}{{\color[rgb]{0,0,0}n_{2}N_{2}}}\right)\\\
&=e^{-\tilde{\nu}t}\left(\left(\frac{E_{1}(0)}{n_{1}}-\frac{1}{2}\frac{|P_{1}(0)|^{2}}{n_{1}N_{1}}\right)-\left(\frac{E_{2}(0)}{n_{2}}-\frac{1}{2}\frac{|P_{2}(0)|^{2}}{n_{2}N_{2}}\right)\right)+\frac{1}{2}m_{1}m_{2}\frac{n_{2}N_{2}-n_{1}N_{1}}{(N_{1}+N_{2})^{2}}e^{-\tilde{\nu}t}(1-e^{-\tilde{\nu}t})\Bigg{|}\frac{P_{2}(0)}{N_{2}}-\frac{P_{1}(0)}{N_{1}}\Bigg{|}^{2}\\\
&+\left(E_{1}(0)+E_{2}(0)-\frac{1}{2}\frac{|P_{1}(0)+P_{2}(0)|^{2}}{N_{1}+N_{2}}\right)e^{-\tilde{\nu}t}\int_{0}^{t}e^{\tilde{\nu}s}\left[\frac{\frac{m_{1}^{3/2}H^{E}_{\tau}(c_{12}(s))}{n_{1}}-\frac{m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21}(s))}{n_{2}}}{m_{1}^{3/2}H^{E}_{\tau}(c_{12})+m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21})}\right]ds.\end{split}$
(22)

###### Proof.

We compute

$\displaystyle\begin{split}\partial_{t}\left(\frac{E_{1}}{n_{1}}-\frac{1}{2}\frac{|P_{1}|^{2}}{n_{1}N_{1}}\right)=\frac{1}{n_{1}}\partial_{t}\left(E_{1}\right)-\frac{P_{1}}{n_{1}}\partial_{t}\left(\frac{P_{1}}{N_{1}}\right)=\tilde{\nu}\left(\frac{E_{12}}{n_{1}}-\frac{E_{1}}{n_{1}}-\frac{P_{1}}{n_{1}}\frac{N_{2}}{N_{1}+N_{2}}\left(\frac{P_{2}}{N_{2}}-\frac{P_{1}}{N_{1}}\right)\right).\end{split}$
(23)

using that $n_{1}$ is constant in the space-homogeneous case, equation (19)
and inserting (21).

Applying (LABEL:ex_en), it follows

$\displaystyle\begin{split}&\partial_{t}\left(\frac{E_{1}}{n_{1}}-\frac{1}{2}\frac{|P_{1}|^{2}}{n_{1}N_{1}}\right)\\\
&=\tilde{\nu}\left(\frac{1}{2}\frac{m_{1}|P_{1}+P_{2}|^{2}}{(N_{1}+N_{2})^{2}}+\frac{(E_{1}+E_{2})-\frac{1}{2}\frac{|P_{1}+P_{2}|^{2}}{N_{1}+N_{2}}}{m_{1}^{3/2}H^{E}_{\tau}(c_{12})+m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21})}\frac{m_{1}^{3/2}H^{E}_{\tau}(c_{12})}{n_{1}}-\frac{E_{1}}{n_{1}}-\frac{P_{1}}{n_{1}}\frac{N_{2}}{N_{1}+N_{2}}\left(\frac{P_{2}}{N_{2}}-\frac{P_{1}}{N_{1}}\right)\right).\end{split}$
(24)

Analogously for species 2 we obtain

$\displaystyle\partial_{t}\left(\frac{E_{2}}{n_{2}}-\frac{1}{2}\frac{|P_{2}|^{2}}{n_{2}N_{2}}\right)$
$\displaystyle=\tilde{\nu}\left(\frac{1}{2}\frac{m_{2}|P_{1}+P_{2}|^{2}}{(N_{1}+N_{2})^{2}}+\frac{(E_{1}+E_{2})-\frac{1}{2}\frac{|P_{1}+P_{2}|^{2}}{N_{1}+N_{2}}}{m_{1}^{3/2}H^{E}_{\tau}(c_{12})+m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21})}\frac{m_{2}^{3/2}H^{E}_{\tau}(c_{21})}{n_{2}}-\frac{E_{2}}{n_{2}}-\frac{P_{2}}{n_{2}}\frac{N_{1}}{N_{1}+N_{2}}\left(\frac{P_{1}}{N_{1}}-\frac{P_{2}}{N_{2}}\right)\right).$

Subtracting both, we get

$\displaystyle\begin{split}&\partial_{t}\left(\left(\frac{E_{1}}{n_{1}}-\frac{1}{2}\frac{|P_{1}|^{2}}{n_{1}N_{1}}\right)-\left(\frac{E_{2}}{n_{2}}-\frac{|P_{2}|^{2}}{n_{2}N_{2}}\right)\right)\\\
&=\tilde{\nu}\large(\frac{E_{2}}{n_{2}}-\frac{E_{1}}{n_{1}}+\frac{\frac{1}{2}(m_{1}-m_{2})n_{1}N_{1}+N_{2}(N_{1}+N_{2})}{(N_{1}+N_{2})^{2}}\frac{|P_{1}|^{2}}{n_{1}N_{1}}+\frac{\frac{1}{2}(m_{1}-m_{2})n_{2}N_{2}-N_{1}(N_{1}+N_{2})}{(N_{1}+N_{2})^{2}}\frac{|P_{2}|^{2}}{n_{2}N_{2}}\\\
&+\frac{(n_{1}N_{1}-n_{2}N_{2})}{(N_{1}+N_{2})^{2}n_{1}n_{2}}P_{1}\cdot
P_{2}+\frac{(E_{1}+E_{2})-\frac{1}{2}\frac{|P_{1}+P_{2}|^{2}}{N_{1}+N_{2}}}{m_{1}^{3/2}H^{E}_{\tau}(c_{12})+m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21})}\left[\frac{m_{1}^{3/2}H^{E}_{\tau}(c_{12})}{n_{1}}-\frac{m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21})}{n_{2}}\right]\large).\end{split}$

This can be written as

$\displaystyle\begin{split}&\partial_{t}\left(\left(\frac{E_{1}}{n_{1}}-\frac{1}{2}\frac{|P_{1}|^{2}}{n_{1}N_{1}}\right)-\left(\frac{E_{2}}{n_{2}}-\frac{1}{2}\frac{|P_{2}|^{2}}{n_{2}N_{2}}\right)\right)\\\
&=\tilde{\nu}\large(-\left(\left(\frac{E_{1}}{n_{1}}-\frac{1}{2}\frac{|P_{1}|^{2}}{n_{1}N_{1}}\right)-\left(\frac{E_{2}}{n_{2}}-\frac{1}{2}\frac{|P_{2}|^{2}}{N_{2}}\right)\right)+\frac{1}{2}\frac{m_{2}(n_{2}N_{2}-n_{1}N_{1})}{(N_{1}+N_{2})^{2}}\frac{|P_{1}|^{2}}{n_{1}N_{1}}+\frac{1}{2}\frac{m_{1}(n_{2}N_{2}-n_{1}N_{1})}{(N_{1}+N_{2})^{2}}\frac{|P_{2}|^{2}}{n_{2}N_{2}}\\\
&+P_{1}\cdot
P_{2}\frac{n_{1}N_{1}-n_{2}N_{2}}{(N_{1}+N_{2})^{2}n_{1}n_{2}}+\frac{(E_{1}+E_{2})-\frac{1}{2}\frac{|P_{1}+P_{2}|^{2}}{N_{1}+N_{2}}}{m_{1}^{3/2}H^{E}_{\tau}(c_{12})+m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21})}\left[\frac{m_{1}^{3/2}H^{E}_{\tau}(c_{12})}{n_{1}}-\frac{m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21})}{n_{2}}\right]\large)\\\
&=\tilde{\nu}\large(-\left(\left(\frac{E_{1}}{n_{1}}-\frac{|P_{1}|^{2}}{n_{1}N_{1}}\right)-\left(\frac{E_{2}}{n_{2}}-\frac{|P_{2}|^{2}}{n_{2}N_{2}}\right)\right)+\frac{1}{2}m_{1}m_{2}\frac{n_{2}N_{2}-n_{1}N_{1}}{(N_{1}+N_{2})^{2}}\Bigg{|}\frac{P_{1}}{N_{1}}-\frac{P_{2}}{N_{2}}\Bigg{|}^{2}\\\
&+\frac{(E_{1}+E_{2})-\frac{1}{2}\frac{|P_{1}+P_{2}|^{2}}{N_{1}+N_{2}}}{m_{1}^{3/2}H^{E}_{\tau}(c_{12})+m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21})}\left[\frac{m_{1}^{3/2}H^{E}_{\tau}(c_{12})}{n_{1}}-\frac{m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21})}{n_{2}}\right]\large).\end{split}$

Now, Duhamels formula gives

$\displaystyle\left(\frac{E_{1}}{N_{1}}-\frac{1}{2}\frac{|P_{1}|^{2}}{m_{1}N_{1}^{2}}\right)-\left(\frac{E_{2}}{N_{2}}-\frac{1}{2}\frac{|P_{2}|^{2}}{n_{2}N_{2}}\right)$
$\displaystyle=e^{-\tilde{\nu}t}\left(\left(\frac{E_{1}(0)}{n_{1}}-\frac{1}{2}\frac{|P_{1}(0)|^{2}}{n_{1}N_{1}}\right)-\left(\frac{E_{2}(0)}{n_{2}}-\frac{1}{2}\frac{|P_{2}(0)|^{2}}{n_{2}N_{2}}\right)\right)+\frac{1}{2}m_{1}m_{2}\frac{n_{2}N_{2}-n_{1}N_{1}}{(N_{1}+N_{2})^{2}}e^{-\tilde{\nu}t}\int_{0}^{t}e^{\tilde{\nu}s}\Bigg{|}\frac{P_{2}(s)}{N_{2}}-\frac{P_{1}(s)}{N_{1}}\Bigg{|}^{2}ds$
$\displaystyle+\left(E_{1}(0)+E_{2}(0)-\frac{1}{2}\frac{|P_{1}(0)+P_{2}(0)|^{2}}{N_{1}+N_{2}}\right)e^{-\tilde{\nu}t}\int_{0}^{t}e^{\tilde{\nu}s}\left[\frac{\frac{m_{1}^{3/2}H^{E}_{\tau}(c_{12}(s))}{n_{1}}-\frac{m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21}(s))}{n_{2}}}{m_{1}^{3/2}H^{E}_{\tau)}(c_{12})+m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21})}\right]ds$
$\displaystyle=e^{-\tilde{\nu}t}\left(\left(\frac{E_{1}(0)}{n_{1}}-\frac{1}{2}\frac{|P_{1}(0)|^{2}}{n_{1}N_{1}}\right)-\left(\frac{E_{2}(0)}{n_{2}}-\frac{1}{2}\frac{|P_{2}(0)|^{2}}{n_{2}N_{2}}\right)\right)+\frac{1}{2}m_{1}m_{2}\frac{n_{2}N_{2}-n_{1}N_{1}}{(N_{1}+N_{2})^{2}}e^{-\tilde{\nu}t}\int_{0}^{t}e^{\tilde{\nu}s}e^{-2\tilde{\nu}s}\Bigg{|}\frac{P_{2}(0)}{N_{2}}-\frac{P_{1}(0)}{N_{1}}\Bigg{|}^{2}ds$
$\displaystyle+\left(E_{1}(0)+E_{2}(0)-\frac{1}{2}\frac{|P_{1}(0)+P_{2}(0)|^{2}}{N_{1}+N_{2}}\right)e^{-\tilde{\nu}t}\int_{0}^{t}e^{\tilde{\nu}s}\left[\frac{\frac{m_{1}^{3/2}H^{E}_{\tau}(c_{12}(s))}{n_{1}}-\frac{m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21}(s))}{n_{2}}}{m_{1}^{3/2}H^{E}_{\tau}(c_{12}(s))+m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21}(s))}\right]ds$
$\displaystyle=e^{-\tilde{\nu}t}\left(\left(\frac{E_{1}(0)}{n_{1}}-\frac{1}{2}\frac{|P_{1}(0)|^{2}}{n_{1}N_{1}}\right)-\left(\frac{E_{2}(0)}{n_{2}}-\frac{1}{2}\frac{|P_{2}(0)|^{2}}{n_{2}N_{2}}\right)\right)+\frac{1}{2}m_{1}m_{2}\frac{n_{2}N_{2}-n_{1}N_{1}}{(N_{1}+N_{2})^{2}}e^{-\tilde{\nu}t}(1-e^{-\tilde{\nu}t})\Bigg{|}\frac{P_{2}(0)}{N_{2}}-\frac{P_{1}(0)}{N_{1}}\Bigg{|}^{2}$
$\displaystyle+\left(E_{1}(0)+E_{2}(0)-\frac{1}{2}\frac{|P_{1}(0)+P_{2}(0)|^{2}}{N_{1}+N_{2}}\right)e^{-\tilde{\nu}t}\int_{0}^{t}e^{\tilde{\nu}s}\left[\frac{\frac{m_{1}^{3/2}H^{E}_{\tau}(c_{12}(s))}{n_{1}}-\frac{m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21}(s))}{n_{2}}}{m_{1}^{3/2}H^{E}_{\tau}(c_{12}(s))+m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21}(s))}\right]ds.$

We get the last but one equality by using theorem 3.3.1 and the last equality
by computing the integral.

∎

######  Remark 5.

In the classical case we get, as in remark 3, the relationship

$\displaystyle\frac{n_{1}}{n_{2}}=\frac{m_{1}^{3/2}}{m_{2}^{3/2}}\frac{e^{-c_{12}}}{e^{-c_{21}}}.$

If we then compute the bracket

$\displaystyle\left[\frac{\frac{m_{1}^{3/2}H^{E}_{\tau}(c_{12}(s))}{n_{1}}-\frac{m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21}(s))}{n_{2}}}{m_{1}^{3/2}H^{E}_{\tau}(c_{12})+m_{2}^{3/2}H^{E}_{\tau^{\prime}}(c_{21})}\right]$

for $\tau=\tau^{\prime}=0$, we obtain by an explicit computation of
$H^{E}_{\tau}$

$\displaystyle\frac{\frac{m_{1}^{3/2}}{n_{1}}\frac{3}{2}e^{-c_{12}}\pi^{3/2}-\frac{m_{2}^{3/2}}{n_{2}}\frac{3}{2}e^{-c_{21}}\pi^{3/2}}{m_{1}^{3/2}\frac{3}{2}e^{-c_{12}}\pi^{3/2}+m_{2}^{3/2}\frac{3}{2}e^{-c_{21}}\pi^{3/2}}=\frac{\frac{m_{1}^{3/2}}{n_{1}}e^{-c_{12}}-\frac{m_{2}^{3/2}}{n_{2}}e^{-c_{21}}}{m_{1}^{3/2}e^{-c_{12}}+m_{2}^{3/2}e^{-c_{21}}}=0$

So in the classical-classical case, the last term of (22) vanishes, and we
obtain the temperature convergence rate

$\displaystyle\left(\frac{E_{1}}{n_{1}}-\frac{1}{2}\frac{|P_{1}|^{2}}{n_{1}N_{1}}\right)-\left(\frac{E_{2}}{n_{2}}-\frac{1}{2}\frac{|P_{2}|^{2}}{n_{2}N_{2}}\right)$
$\displaystyle=e^{-\tilde{\nu}t}\left(\left(\frac{E_{1}(0)}{n_{1}}-\frac{1}{2}\frac{|P_{1}(0)|^{2}}{n_{1}N_{1}}\right)-\left(\frac{E_{2}(0)}{n_{2}}-\frac{1}{2}\frac{|P_{2}(0)|^{2}}{n_{2}N_{2}}\right)\right)+\frac{1}{2}m_{1}m_{2}\frac{n_{2}N_{2}-n_{1}N_{1}}{(N_{1}+N_{2})^{2}}e^{-\tilde{\nu}t}(1-e^{-\tilde{\nu}t})\Bigg{|}\frac{P_{2}(0)}{N_{2}}-\frac{P_{1}(0)}{N_{1}}\Bigg{|}^{2}.$

##  4 Motivation of the structure of the local equilibria

In this section, we motivate the form of the local equilibria in (1) This will
motivate later the choice of our numerical scheme. For this, it will be
convenient to define the following notations. We denote

$\displaystyle\begin{split}\mathcal{K}_{11}=(e^{-\alpha_{11}\cdot\mathbf{p}_{1}}+\tau)^{-1},\quad\mathcal{K}_{22}=(e^{-\alpha_{22}\cdot\mathbf{p}_{2}}+\tau^{\prime})^{-1},\\\
\mathcal{K}_{12}=(e^{-\alpha_{12}\cdot\mathbf{p}_{1}}+\tau)^{-1},\quad\mathcal{K}_{12}=(e^{-\alpha_{21}\cdot\mathbf{p}_{2}}+\tau^{\prime})^{-1}\end{split}$
(25)

where

$\displaystyle\mathbf{p}_{k}(p):=(1,p,\frac{|p|^{2}}{2m_{k}})^{\top},\quad
k=1,2$ (26)

and the parameters $(a_{1},a_{2},a,b_{1},b_{2},b,c_{1},c_{2},c_{12},c_{21})$
can be mapped one-to-one to
$\alpha_{kj}=(\alpha^{0}_{k}j,\alpha^{1}_{kj},\alpha^{2}_{kj}),~{}k,j=1,2$.
Further, we recall the function  $h_{\tau}$ defined by (12) is convex,
therefore it follows

$h_{\tau}(z)\geq h_{\tau}(y)+\ln(y)(z-y),$ (27)

for all $y,z>0$ if $\tau=0,-1$ and $0<y,z<1$ if $\tau=+1$.

###  4.1 The one species local equilibria

We seek a solution of the entropy minimization problem

$\displaystyle\min_{g\in\chi_{k}}\int{\color[rgb]{0,0,0}h_{\tau}(g)}dv,\quad
k\in\\{1,2\\},$ (28)

where

$\displaystyle\chi_{k}=\left\\{g~{}\Big{|}~{}g\geq 0,\,(1+|p|^{2})g\in
L^{1}(\mathbb{R}^{3}),\,\int\mathbf{p}_{k}(p)(g-f_{k})dp=0\right\\}.$ (29)

The choice of the set $\chi_{k}$ ensures the conservation properties (5) for
intra-species collisions. Indeed, by standard optimization theory, any
critical point $(\mathcal{K}_{kk},\lambda^{kk})$ of the Lagrange functional
$L_{k}\colon\chi_{k}\times\mathbb{R}^{5}\to\mathbb{R}$, given by

$\displaystyle
L_{k}(g,\alpha)=\int{\color[rgb]{0,0,0}h(g)}dp-\alpha\cdot\int\mathbf{p}_{k}(p)(g-f_{k})dp,$
(30)

satisfies the first-order optimality condition

$\displaystyle\frac{\delta L_{k}}{\delta
g}(\mathcal{K}_{kk},\lambda^{kk})=(\ln\frac{\mathcal{K}_{kk}}{1-\tau_{k}\mathcal{K}_{kk}}-\lambda^{kk}\cdot\mathbf{p}_{k}(p))=0,$
(31)

with $\tau_{1}=\tau$ and $\tau_{2}=\tau^{\prime}$. This implies then that

$\displaystyle\mathcal{K}_{kk}=(\exp{(\lambda^{kk}\cdot\mathbf{p}_{k}(p))}+\tau_{k})^{-1}.$
(32)

Theorem 3.1.1 shows in a rigorous way that there exists a unique function of
the form (32) that satisfies these constraints. Therefore, we can prove the
following theorem.

######  Theorem 4.1.1.

The local equilibrium $\mathcal{K}_{kk}$ is the unique minimizer of (28).

###### Proof.

According to (27)

${\color[rgb]{0,0,0}h_{\tau_{k}}(g)}\geq{\color[rgb]{0,0,0}h_{\tau_{k}}}(\mathcal{K}_{kk})+\lambda^{kk}\cdot\mathbf{p}_{k}(g-\mathcal{K}_{kk}),$
(33)

point-wise in $p$. Thus it follows that for all $g\in\chi_{k}$,

$\displaystyle\int{\color[rgb]{0,0,0}h_{\tau_{k}}(g)}dp\geq\int{\color[rgb]{0,0,0}h_{\tau_{k}}}(\mathcal{K}_{kk})dp+\int\lambda^{kk}\cdot\mathbf{p}_{k}(g-\mathcal{K}_{kk})dp=\int{\color[rgb]{0,0,0}h_{\tau_{k}}}(\mathcal{K}_{kk})dp$
(34)

Hence $\mathcal{K}_{kk}$ is a minimizer of (28), and uniqueness follows
directly from the strict convexity of  $h_{\tau_{k}}$. ∎

###  4.2 The mixture local equilibria

For interactions between species, we seek a solution of the entropy
minimization problem

$\displaystyle\min_{g_{1},g_{2}\in\chi_{12}}\int h_{\tau}(g_{1})dp+\int
h_{\tau^{\prime}}(g_{2})dp,$ (35)

where

$\displaystyle\begin{split}\chi_{12}=\Bigg{\\{}(g_{1},g_{2})~{}\Big{|}&~{}g_{1},g_{2}>0,\,(1+|p|^{2})g_{1},\,(1+|p|^{2})g_{2}\in
L^{1}(\mathbb{R}^{3}),\\\ &\int g_{1}dp=\int f_{1}dp,\quad\int g_{2}dp=\int
f_{2}dp,\\\ &\int\begin{pmatrix}p\\\
\frac{|p|^{2}}{2m_{1}}\end{pmatrix}(g_{1}-f_{1})dp+\int\begin{pmatrix}p\\\
\frac{|p|^{2}}{2m_{2}}\end{pmatrix}(g_{2}-f_{2})dp=0\Bigg{\\}}.\end{split}$
(36)

Here, $\chi_{12}$ is chosen such that the constraints (6) for inter-species
collisions are satisfied. Similar to the case of intra-species collisions, we
consider the Lagrange functional
$L\colon\chi\times\mathbb{R}\times\mathbb{R}\times\mathbb{R}^{3}\times\mathbb{R}\to\mathbb{R}$

$\displaystyle\begin{split}L(g_{1},g_{2},\alpha_{0}^{1},\alpha_{0}^{2},\alpha_{1},\alpha_{2})&=\int
h(g_{1})dp+\int h(g_{2})dp\\\
&-\alpha_{0}^{1}\int(g_{1}-f_{1})dp-\alpha_{0}^{2}\int(g_{2}-f_{2})dp\\\
&-\alpha_{1}\cdot\left(\int p(g_{1}-f_{1})dp+\int p(g_{2}-f_{2})dp\right)\\\
&-\alpha_{2}\left(\int\frac{|p|^{2}}{2m_{1}}(g_{1}-f_{1})dp+\int\frac{|p|^{2}}{2m_{2}}(g_{2}-f_{2})dp\right).\\\
\end{split}$ (37)

Any critical point
$(\mathcal{K}_{12},\mathcal{K}_{21},\lambda_{0}^{1},\lambda_{0}^{2},\lambda_{1},\lambda_{2})$
of $L$ satisfies the first-order optimality conditions

$\displaystyle\frac{\delta L}{\delta
g_{1}}(\mathcal{K}_{12},\mathcal{K}_{21},\lambda_{0}^{1},\lambda_{0}^{2},\lambda_{1},\lambda_{2})=\ln\frac{\mathcal{K}_{12}}{1-\tau\mathcal{K}_{12}}-\lambda^{12}\cdot\mathbf{p}_{1}(p)=0,$
(38) $\displaystyle\frac{\delta L}{\delta
g_{2}}(\mathcal{K}_{12},\mathcal{K}_{21},\lambda_{0}^{1},\lambda_{0}^{2},\lambda_{1},\lambda_{2})=\ln\frac{\mathcal{K}_{21}}{1-\tau\mathcal{K}_{21}}-\lambda^{21}\cdot\mathbf{p}_{2}(p)=0,$
(39)

where $\lambda^{12}=(\lambda_{0}^{1},\lambda_{1},\lambda_{2})$ and
$\lambda^{21}=(\lambda_{0}^{2},\lambda_{1},\lambda_{2})$. Therefore

$\displaystyle\mathcal{K}_{12}$
$\displaystyle=(\exp(\lambda^{12}\cdot\mathbf{p}_{1}(p))+\tau)^{-1},$ (40)
$\displaystyle\mathcal{K}_{21}$
$\displaystyle=(\exp(\lambda^{21}\cdot\mathbf{p}_{2}(p))+\tau)^{-1}.$ (41)

Since we only require conservation of the combined momentum and kinetic
energy, there is only one Lagrange multiplier for the momentum constraint and
one Lagrange multiplier for the energy constraint. Hence,
$\lambda_{1}^{12}=\lambda_{1}^{21}$ and $\lambda_{2}^{12}=\lambda_{2}^{21}$ in
(41). When we are in the classical case ($\tau=0$), this restriction is the
same as the one used in [21], but more restrictive than the model in [28].

Theorem 3.1.2 shows the existence of functions of the form (41) which satisfy
the constraints in (5) and (6). As in the single species case, it follows that
these functions are unique minimizer of the corresponding minimization
problem.

######  Theorem 4.2.1.

$(\mathcal{K}_{12},\mathcal{K}_{21})$ as defined in (41) is the unique
minimizer of (35).

###### Proof.

According to (27)

$h_{\tau}(g)\geq
h_{\tau_{k}}(\mathcal{K}_{kj})+\lambda^{kj}\cdot\mathbf{p}_{k}(g-\mathcal{K}_{kj}),$
(42)

point-wise in $p$, for any measurable function $g$ and $k,j\in\\{1,2\\}$.
Therefore it follows that for any measureable functions $g_{1}$ and $g_{2}$,

$\displaystyle\int h_{\tau}(g_{1})dp+\int$ $\displaystyle
h_{\tau^{\prime}}(g_{2})dp\geq\int h_{\tau}(\mathcal{K}_{12})dp+\int
h_{\tau^{\prime}}(\mathcal{K}_{21})dp$
$\displaystyle\quad+\lambda^{12}\cdot\int\mathbf{p}_{1}(g_{1}-\mathcal{K}_{12})dp+\lambda^{21}\cdot\int\mathbf{p}_{2}(g_{2}-\mathcal{K}_{21})dp.$
(43)

Since $\lambda_{1}^{12}=\lambda_{1}^{21}$ and
$\lambda_{2}^{12}=\lambda_{2}^{21}$,

$\displaystyle\begin{split}\lambda^{12}\cdot\int\mathbf{p}_{1}&(g_{1}-\mathcal{K}_{12})dp+\lambda^{21}\cdot\int\mathbf{p}_{2}(g_{2}-\mathcal{K}_{21})dp=\lambda_{0}^{12}\int(g_{1}-\mathcal{K}_{12})dp+\lambda_{0}^{21}\int(g_{2}-\mathcal{K}_{21})dp\\\
&\quad+\lambda_{1}^{12}\cdot\left(\int p(g_{1}-\mathcal{K}_{12})dp+\int
p(g_{2}-\mathcal{K}_{21})dp\right)+\lambda_{2}^{12}\left(\int\frac{|p|^{2}}{2m_{1}}(g_{1}-\mathcal{K}_{12})dp+\int\frac{|p|^{2}}{2m_{2}}(g_{2}-\mathcal{K}_{21})dp\right).\end{split}$
(44)

If $(g_{1},g_{2})$ and $(\mathcal{K}_{12},\mathcal{K}_{21})$ are elements of
$\chi_{12}$, then the constraints in (36) imply that each of the terms above
is zero. In such cases, (4.2) reduces to

$\displaystyle\int h_{\tau}(g_{1})dp+\int h_{\tau^{\prime}}(g_{2})dp\geq\int
h_{\tau}(\mathcal{K}_{12})dp+\int h_{\tau^{\prime}}(\mathcal{K}_{21})dp,$ (45)

which shows that $(\mathcal{K}_{12},\mathcal{K}_{21})$ solves (35). Since
$h_{\tau_{k}}$ is strictly convex, it follows that this solution is unique.

∎

##  5 Numerical scheme

###  5.1 Time discretization

Let $k,j=1,2$ and $k\neq j$. We write (1) as

$\displaystyle\partial_{t}f_{k}+\mathcal{T}_{k}(f_{k})=\mathcal{R}_{k}(f_{k},f_{j})$
(46)

with the combined relaxation operator

$\displaystyle\mathcal{R}_{k}(f_{k},f_{j})=\mathcal{R}_{kk}+\mathcal{R}_{kj}=\nu_{kk}n_{k}\left(\mathcal{K}_{kk}-f_{k}\right)+\nu_{kj}n_{j}\left(\mathcal{K}_{kj}-f_{k}\right)$
(47)

and the transport operator

$\displaystyle\mathcal{T}_{k}(f_{k})=\frac{p}{m_{k}}\cdot\nabla_{x}f_{k}.$
(48)

In the following, for simplicity, we assume that the collision frequencies
$\tilde{\nu}_{kk}:=\nu_{kk}n_{k}$ and $\tilde{\nu}_{kj}:=\nu_{kj}n_{k}$ are
constant in $x$ and $t$. But an extension to an $x$ and $t$ dependence of the
collision frequency would also be possible. Large collision frequencies result
in a stiff relaxation operator such that an implicit time discretization for
the relaxation part is a convenient choice. So we pursue implicit-explicit
(IMEX) schemes where $\mathcal{R}_{k}$ is treated implicitly and
$\mathcal{T}_{i}$ is treated explicitly.

Given $t_{\ell}=\ell\Delta t$ for $\ell\in\mathbb{N}_{0}$, a simple update of
$f_{k}^{\ell}\approx f_{k}(x,p,t_{\ell})$ from $t_{\ell}$ to $t_{\ell+1}$ uses
the approximation

$\displaystyle\mathcal{R}_{k}(f_{k}^{\ell+1},f_{j}^{\ell+1})\approx\tilde{\nu}_{kk}\left(\mathcal{K}_{kk}^{\ell+1}-f_{k}^{\ell+1}\right)+\tilde{\nu}_{kj}\left(\mathcal{K}_{kj}^{\ell+1}-f_{k}^{\ell+1}\right),$
(49)

where $\mathcal{K}_{kk}^{\ell+1}$ and $\mathcal{K}_{kj}^{\ell+1}$ are discrete
target functions that, as described in Section 5.1.3, depend on
$f_{k}^{\ell+1}$ and $f_{j}^{\ell+1}$ via the solution of a convex
minimization problem that is inspired by the work in [23]. By this procedure,
$\mathcal{K}_{kk}$ and $\mathcal{K}_{kj}$ are evaluated exactly at the next
time step (up to numerical tolerances) which results in the preservation of
conservation properties, and the first-order version inherits additional
properties from the continuum model.

####  5.1.1 First-order splitting

We split the relaxation and transport operators in (1).

##### Relaxation.

We perform the relaxation step in each spatial cell by a backward Euler method

$\displaystyle\frac{f_{k}^{\ast}-f_{k}^{\ell}}{\Delta
t}=\mathcal{R}_{k}(f_{k}^{\ast},f_{j}^{\ast}),$ (50)

which can be rewritten into the convex combination

$\displaystyle f_{k}^{\ast}=d_{k}f_{k}^{\ell}+d_{k}\Delta
t(\tilde{\nu}_{kk}\mathcal{K}_{kk}^{\ast}+\tilde{\nu}_{kj}\mathcal{K}_{kj}^{\ast})$
(51)

with

$\displaystyle d_{k}=\frac{1}{1+\Delta t(\tilde{\nu}_{kk}+\tilde{\nu}_{kj})}.$
(52)

The equation (51) represents an explicit update formula for $f_{k}^{\ast}$
provided that $\mathcal{K}_{kk}^{\ast}$ and $\mathcal{K}_{kj}^{\ast}$ can be
expressed as functions of $f_{k}^{\ell}$. In Section 5.1.3 we show how to
determine $\mathcal{K}_{kk}^{\ast}$ and $\mathcal{K}_{kj}^{\ast}$ in a
structure-preserving way.

##### Transport.

We compute the transport in $x$ for $f_{k}^{\ell+1}$ by a forward Euler method
with initial data $f_{k}^{\ast}$:

$\displaystyle\frac{f_{k}^{\ell+1}-f_{k}^{\ast}}{\Delta
t}+\mathcal{T}_{k}(f_{k}^{\ast})=0.$ (53)

Details on the numerical approximation of $\mathcal{T}_{k}$ are presented in
section 5.2.

####  5.1.2 Second-order IMEX Runge-Kutta

We use the following Butcher tableaux [2] for a second-order approach

$\displaystyle\begin{array}[]{c|ccc}0\\\ \gamma&0&\gamma\\\
1&0&1-\gamma&\gamma\\\ \hline\cr&0&1-\gamma&\gamma\end{array}\hskip
56.9055pt\begin{array}[]{c|ccc}0\\\ \gamma&\gamma\\\ 1&\delta&1-\delta&0\\\
\hline\cr&\delta&1-\delta&0\end{array}$ (62)

with

$\gamma=1-\frac{\sqrt{2}}{2}\quad\text{and}\quad\delta=1-\frac{1}{2\gamma}.$
(63)

The left table applies to the relaxation part, and the right table applies to
the transport terms. This IMEX Runge-Kutta scheme is L-stable and globally
stiffly accurate.
Applying this method to (1) and using the constants

$\displaystyle d_{k}=\frac{1}{1+\gamma\Delta
t(\tilde{\nu}_{kk}+\tilde{\nu}_{kj})},$ (64)

we can write the stages in the scheme as convex combination of three terms

$\displaystyle f_{k}^{(1)}$ $\displaystyle=d_{k}G_{k}^{(1)}$
$\displaystyle+d_{k}\gamma\Delta t\,\tilde{\nu}_{kk}\mathcal{K}_{kk}^{(1)}$
$\displaystyle+d_{k}\gamma\Delta t\,\tilde{\nu}_{kj}\mathcal{K}_{kj}^{(1)}$
(65a) $\displaystyle f_{k}^{(2)}$ $\displaystyle=d_{k}G_{k}^{(2)}$
$\displaystyle+d_{k}\gamma\Delta t\,\tilde{\nu}_{kk}\mathcal{K}_{kk}^{(2)}$
$\displaystyle+d_{k}\gamma\Delta t\,\tilde{\nu}_{kj}\mathcal{K}_{kj}^{(2)},$
(65b) $\displaystyle f_{k}^{\ell+1}$ $\displaystyle=f_{k}^{(2)}$ (65c)

where

$\displaystyle G_{k}^{(1)}$ $\displaystyle=f_{k}^{\ell}-\Delta
t\,\gamma\,\mathcal{T}_{k}(f_{k}^{\ell})$ (66a) $\displaystyle G_{k}^{(2)}$
$\displaystyle=f_{k}^{\ell}-\Delta
t\,\delta\,\mathcal{T}_{k}(f_{k}^{\ell})-\Delta
t\,(1-\delta)\mathcal{T}_{k}(f_{k}^{(1)})+\Delta
t\,(1-\gamma)\mathcal{R}_{k}(f_{k}^{(1)},f_{j}^{(1)})$ (66b)

depend on known data. For each stage, we have to determine the corresponding
values of the target functions in order to update the distribution functions.
In the following section, we explain how this can be done.

####  5.1.3 General implicit solver

We write the implicit updates in (51) and (65) in a generic steady state form

$\displaystyle\psi_{k}=d_{k}G_{k}+d_{k}\gamma\Delta
t(\tilde{\nu}_{kk}\mathcal{K}_{kk}+\tilde{\nu}_{kj}\mathcal{K}_{kj}).$ (67)

The functions $\mathcal{K}_{kk}$ and $\mathcal{K}_{kj}$ are the unique target
functions associated to $\psi_{k}$,

$d_{k}=\frac{1}{1+\gamma\Delta t(\tilde{\nu}_{kk}+\tilde{\nu}_{kj})},$ (68)

and $G_{k}$ is a known function. We want to express $\mathcal{K}_{kk}$ and
$\mathcal{K}_{kj}$ as functions of $G_{k}$ and $G_{j}$ so that (67) is an
explicit update formula for $\psi_{k}$. In Section 3, the existence and
uniqueness of $\mathcal{K}_{kk}$ and $\mathcal{K}_{kj}$ are proven by
algebraic considerations. In order to determine their values we follow a
different approach inspired by section 4. Applying the conservation properties
(5) and (6) to (67) leads to

$\displaystyle\begin{split}\int\tilde{\nu}_{11}\mathcal{K}_{11}\,\mathbf{p}_{1}dp+\int\tilde{\nu}_{22}\mathcal{K}_{22}\,\mathbf{p}_{2}dp+\int\tilde{\nu}_{12}\mathcal{K}_{12}\,\mathbf{p}_{1}dp+\int\tilde{\nu}_{21}\mathcal{K}_{21}\,\mathbf{p}_{2}dp\end{split}$
$\displaystyle\begin{split}\overset{\eqref{conserv 1},\eqref{conserv
2}}{=}\int\tilde{\nu}_{11}\psi_{1}\,\mathbf{p}_{1}dp+\int\tilde{\nu}_{22}\psi_{2}\,\mathbf{p}_{2}dp+\int\tilde{\nu}_{12}\psi_{1}\,\mathbf{p}_{1}dp+\int\tilde{\nu}_{21}\psi_{2}\,\mathbf{p}_{2}dp\end{split}$
$\displaystyle\begin{split}&\hskip
7.11317pt\overset{(\ref{eq:update_general})}{=}\int\tilde{\nu}_{11}d_{1}\left[G_{1}+\Delta
t\,\gamma\tilde{\nu}_{11}\mathcal{K}_{11}+\Delta
t\,\gamma\tilde{\nu}_{12}\mathcal{K}_{12}\right]\,\mathbf{p}_{1}dp+\int\tilde{\nu}_{22}d_{2}\left[G_{2}+\Delta
t\,\gamma\tilde{\nu}_{22}\mathcal{K}_{22}+\Delta
t\,\gamma\tilde{\nu}_{21}\mathcal{K}_{21}\right]\,\mathbf{p}_{2}dp\\\ &\hskip
42.67912pt+\int\tilde{\nu}_{12}d_{1}\left[G_{1}+\Delta
t\,\gamma\tilde{\nu}_{11}\mathcal{K}_{11}+\Delta
t\,\gamma\tilde{\nu}_{12}\mathcal{K}_{12}\right]\,\mathbf{p}_{1}dp+\int\tilde{\nu}_{21}d_{2}\left[G_{2}+\Delta
t\,\gamma\tilde{\nu}_{22}\mathcal{K}_{22}+\Delta
t\,\gamma\tilde{\nu}_{21}\mathcal{K}_{21}\right]\,\mathbf{p}_{2}dp\end{split}$
(69)

Sorting terms yields the following moment equations

$\displaystyle\begin{split}\int
d_{1}\left(\tilde{\nu}_{11}\mathcal{K}_{11}+\tilde{\nu}_{12}\mathcal{K}_{12}\right)\mathbf{p}_{1}dp+\int
d_{2}\left(\tilde{\nu}_{21}\mathcal{K}_{21}+\tilde{\nu}_{22}\mathcal{K}_{22}\right)\mathbf{p}_{2}dp\\\
=\int d_{1}(\tilde{\nu}_{11}+\tilde{\nu}_{12})G_{1}\mathbf{p}_{1}+\int
d_{2}(\tilde{\nu}_{21}+\tilde{\nu}_{22})G_{2}\mathbf{p}_{2}dp\end{split}$ (70)

which provide a set of constraints to determine $\mathcal{K}_{kk}$ and
$\mathcal{K}_{kj}$ from the given data $G_{k}$ and $G_{j}$. These represent
first-order optimality conditions associated to the minimization of the convex
potential function

$\varphi_{\rm
tot}(\alpha_{1},\alpha_{2},\alpha)=\int\left[d_{1}\tilde{\nu}_{11}w(\mathcal{K}_{11})+d_{2}\tilde{\nu}_{22}w(\mathcal{K}_{22})+d_{1}\tilde{\nu}_{12}w(\mathcal{K}_{12})+d_{2}\tilde{\nu}_{21}w(\mathcal{K}_{21})\right]dp\\\
-\mu_{1}\cdot\alpha_{1}-\mu_{2}\cdot\alpha_{2}-\mu\cdot\alpha$ (71)

where $\alpha_{k}=(\alpha^{0}_{k},\alpha^{1}_{k},\alpha^{2}_{k})$ and
$\alpha=(\alpha_{12}^{0},\alpha_{21}^{0},\alpha^{1},\alpha^{2})^{\top}$; the
auxiliary function reads

$\displaystyle
w(\mathcal{K}_{kj,\tau})=\begin{cases}\mathcal{K}_{kj,0}\quad&{\rm
for}\quad\tau_{k}=0,\\\ \log(1-\tau_{k}\mathcal{K}_{kj,\tau_{k}})\quad&{\rm
else};\end{cases}$ (72)

the given moments are

$\displaystyle\mu_{k}=\begin{pmatrix}\mu_{k}^{0}\\\ \mu_{k}^{1}\\\
\mu_{k}^{2}\end{pmatrix}=\int d_{k}\tilde{\nu}_{kk}G_{k}\mathbf{p}_{k}dp$ (73)

for $k=1,2$; and

$\displaystyle\mu=\begin{pmatrix}\mu_{12}^{0}\\\ \mu_{21}^{0}\\\ \mu^{1}\\\
\mu^{2}\end{pmatrix}$ $\displaystyle=\int\left[\begin{pmatrix}1\\\ 0\\\ p\\\
\frac{|p|^{2}}{2m_{1}}\end{pmatrix}d_{1}\tilde{\nu}_{12}G_{1}+\begin{pmatrix}0\\\
1\\\ p\\\
\frac{|p|^{2}}{2m_{2}}\end{pmatrix}d_{2}\tilde{\nu}_{21}G_{2}\right]dp.$ (74)

The minimization problem can be decoupled as follows:

######  Proposition 5.1.

The components of the minimizer of (71) can be found by minimizing the
following three convex potential functions independently:

$\displaystyle\varphi_{k}(\alpha_{k})$ $\displaystyle=\int
d_{k}\tilde{\nu}_{kk}\,w(\mathcal{K}_{kk})dp+\mu_{k}\cdot\alpha_{k}\quad\text{for}\quad
k=1,2\quad\text{and}$ (75) $\displaystyle\varphi(\alpha)$
$\displaystyle=\int\left[d_{1}\tilde{\nu}_{12}w(\mathcal{K}_{12})+d_{2}\tilde{\nu}_{21}w(\mathcal{K}_{21})\right]dp+\mu\cdot\alpha$
(76)

and the minimum of (71) is the sum of their minima.

###### Proof.

The statement is trivial because $\varphi_{\rm
tot}(\alpha_{1},\alpha_{2},\alpha)=\varphi_{1}(\alpha_{1})+\varphi_{2}(\alpha_{2})+\varphi_{1}(\alpha)$.
∎

The minimum of each potential function in (75) and (76) is found using
Newton’s method for convex optimization. More details are given in Section
5.4.

Actually, we can link these potential functions to dual problems when we
reformulate the modelling problem by using Lagrange functionals. For intra-
species interactions, the Lagrange functional reads

$\displaystyle L_{k}(g,\lambda)=\int
h_{\tau_{k}}(g)dp-\lambda\cdot\int\mathbf{p}_{k}(g-f_{k})dp$ (77)

using $h_{\tau_{k}}(g)$ given by (12). The first integral in (77) is the
entropy functional; the other integrals describe the conservation properties
as constraints. The Lagrange multipliers $\lambda$ solve the dual problem

$\displaystyle\alpha_{k}=\operatorname*{argmin}_{\lambda\in\Lambda_{k}}\int
w(\mathcal{K}_{kk}(\lambda))dp+\lambda\cdot\int\mathbf{p}_{k}f_{k}dp$ (78)

where
$\Lambda_{k}=\\{\lambda\in\mathbb{R}^{5}\,|\,\int\mathcal{K}_{kk}(\lambda)(1+|p|^{2})dp<\infty\\}$.
Analogously, we can formulate the dual problem for inter-species interactions:

$\displaystyle\begin{split}(\alpha_{12},\alpha_{21})=\operatorname*{argmin}_{(\lambda_{12},\lambda_{21})\in\Lambda}\Bigg{\\{}\int&w(\mathcal{K}_{12}(\lambda))+w(\mathcal{K}_{21}(\lambda))dp+\lambda_{12}^{0}\int
f_{1}dp+\lambda_{21}^{0}\int f_{2}dp\\\ &+\lambda^{1}\cdot\int
p(f_{1}+f_{2})dp+\lambda^{2}\int|p|^{2}\left(\frac{1}{2m_{1}}f_{1}+\frac{1}{2m_{2}}f_{2}\right)dp\Bigg{\\}}\end{split}$
(79)

for $\alpha_{kj}=(\alpha_{kj}^{0},\alpha^{1},\alpha^{2})$ and where
$\Lambda=\\{(\lambda_{12},\lambda_{21})\in\mathbb{R}^{5}\,|\,\int\mathcal{K}_{kj}(\lambda_{kj})(1+|p|^{2})dp<\infty\text{
for }k,j=1,2;k\neq j\\}$. We recognize the close relationship of (75) with
(78), respective of (76) with (79). The dual problems have unique solutions
according to Section 4. This is inherited to the potential functions because
$d_{k}\tilde{\nu}_{kj}$ is independent of $p$.

###  5.2 Space discretization

We assume a slab geometry, i.e.
$\partial_{x^{2}}f_{k}=\partial_{x^{3}}f_{k}=0$. So we reduce the physical
space dimension to one dimension and set $x:=x^{1}$ while the momentum domain
remains three dimensional $(p=(p^{1},p^{2},p^{3}))$. We divide the spatial
domain $[x_{\rm min},x_{\rm max}]$ into uniform cells
$I_{i}=[x_{i}-\frac{\Delta x}{2},x_{i}+\frac{\Delta x}{2}]$ for
$i\in\\{0,\dots,I\\}$.

We employ a second-order finite volume framework using approximate cell-
averaged quantities

$f_{k,i}^{\ell}\approx\frac{1}{\Delta x}\int_{I_{i}}f_{k}(x,p,t^{\ell})dx.$
(80)

The relaxation operators are approximated to second order by

$\mathcal{R}_{k,i}^{\ell}=\mathcal{R}_{k}(f_{k,i}^{\ell},f_{j,i}^{\ell})\approx\frac{1}{\Delta
x}\int_{I_{i}}\mathcal{R}\left(\,f_{k}(x,p,t^{\ell}),f_{j}(x,p,t^{\ell})\,\right)dx.$
(81)

Whereas the transport operator $\mathcal{T}_{k}$ is discretized with numerical
fluxes $\mathscr{F}_{i+\frac{1}{2}}$ by

$\displaystyle\mathcal{T}_{k}\left(g\right)\approx\mathcal{T}_{i;k}(g)=\frac{1}{\Delta
x}\left(\mathscr{F}_{k+\frac{1}{2}}(g)-\mathscr{F}_{i-\frac{1}{2}}(g)\right)$
(82)

for any grid function $g=\\{g_{i}\\}$. We follow [33] and use

$\displaystyle\mathscr{F}_{i+\frac{1}{2}}(g)=\frac{p^{1}}{2m}\left(g_{i+1}+g_{i}\right)-\frac{|p^{1}|}{2m}\left(g_{i+1}-g_{i}-\phi_{i+\frac{1}{2}}(g)\right)$
(83)

where $\phi_{i+\frac{1}{2}}$ is a flux limiter. The choice
$\phi_{i+\frac{1}{2}}=0$ leads to a first-order approximation, and a second-
order method is provided by

$\displaystyle\phi_{i+\frac{1}{2}}(g)=\operatorname{minmod}\left((g_{i}-g_{i-1}),(g_{i+1}-g_{i}),(g_{i+2}-g_{i+1})\right)$
(84)

where

$\displaystyle\operatorname{minmod}(a,b,c)=\begin{cases}s\min(|a|,|b|,|c|),\quad&\mathrm{sign}(a)=\mathrm{sign}(b)=\mathrm{sign}(c)=:s,\\\
0,\quad&\text{otherwise}.\end{cases}$ (85)

We guarantee positivity during a simple forward Euler update of (53) by
enforcing the CFL condition

$\displaystyle\Delta t<\alpha\frac{m\Delta x}{\max|p^{1}|}$ (86)

with $\alpha=1$ for the first-order flux and $\alpha=\frac{2}{3}$ for the
second-order flux. (See Proposition 5.2.)

###  5.3 Properties of the semi-discrete scheme

In this section, we review the positivity preservation, conservation
properties, and the entropy behavior of the semi-discrete scheme.

####  5.3.1 Positivity of distribution functions

The first-order time stepping scheme in Section 5.1.1 preserves positivity for
both first- and second-order numerical fluxes in space; see Proposition 5.2.
We discuss the positivity for the second-order scheme 5.1.2 in Proposition
5.3, and give a sufficient criterion for the space homogeneous case.
Additionally, we show that the upper bound for distribution functions of
fermions is preserved by our scheme; see Proposition 5.4.

######  Proposition 5.2.

The first-order time discretization in Section 5.1.1 together with the space
discretization described in Section 5.2 is positivity preserving, provided
that

$\displaystyle\Delta t\leq\beta\frac{m_{k}\Delta x}{\max|p^{1}|},$ (87)

with $\beta=1$ and $\beta=\frac{2}{3}$ for the first-order and second-order
fluxes, respectively.

###### Proof.

The proof can be performed analogously to the proof of proposition 5.1 in
[23]. ∎

Second-order time-stepping makes it more difficult to guarantee positivity.
Nevertheless, we derive some sufficient conditions on $\Delta t$ in order to
preserve positivity in the second-order scheme presented in Section 5.1.2.

######  Proposition 5.3.

For the space homogeneous case, the second-order IMEX scheme presented in
Section 5.1.2 is positivity preserving provided that

$\displaystyle\Delta
t\leq\frac{1}{(1-2\gamma)(\tilde{\nu}_{kk}+\tilde{\nu}_{kj})}$ (88)

for $k,j=1,2.$

###### Proof.

The proof can be performed analogously to the proof of proposition 5.2 in
[23]. ∎

For large collision frequencies $\nu_{kj}$, the time step condition (88) can
be restrictive. So one might be interested in enforcing the milder (but still
sufficient) local condition

$\displaystyle\Delta
t\leq\frac{f_{k}^{\ell}}{(1-\gamma)\left[(\tilde{\nu}_{kk}+\tilde{\nu}_{kj})f_{k}^{(1)}-(\tilde{\nu}_{kk}A_{kk}^{(1)}+\tilde{\nu}_{kj}A_{kj}^{(1)})\right]}.$
(89)

Large collision frequencies push the numerical kinetic distribution to the
corresponding target function. Hence, the denominator in (89) becomes large,
and the condition is not restrictive.

A distribution function of a fermion has the additional upper bound $f<1$. Our
scheme preserves this property which is shown in the following propositions.

######  Proposition 5.4.

If $f_{k}$ represents the distribution function of a fermion with
$f_{k}^{\ell}<1$, the time discretization in Section 5.1.1 together with the
space discretization described in Section 5.2 leads to $f_{k}^{\ell+1}<1$.

###### Proof.

Let $f_{k}^{\ell}<1$. The local equilibrium of a fermion is a Fermi-Dirac
distribution function $\mathcal{F}$ for which $0<\mathcal{F}<1$ by definition.
Hence, for the relaxation step it holds

$\displaystyle f_{k}^{*}=d_{k}f_{k}^{\ell}+d_{k}\Delta
t\,(\tilde{\nu}_{kk}\mathcal{F}_{kk}^{*}+\tilde{\nu}_{kj}\mathcal{F}_{kj}^{*})<d_{k}+d_{k}\Delta
t(\tilde{\nu}_{kk}+\tilde{\nu}_{kj})=1.$ (90)

Here, we used the definition of $d_{k}$ given by (64). For the transport step
(53), the first-order fluxes lead to

$\displaystyle f_{k,i}^{\ell+1}$ $\displaystyle=\left(1-\frac{\Delta
t}{m_{k}\Delta x}|p^{1}|\right)f_{k,i}^{\ast}+\frac{\Delta t}{m_{k}\Delta
x}|p^{1}|f_{k,i-\operatorname{sign}(p^{1})}^{\ast}\overset{\eqref{eq:relax-
smaller-1}}{\leq}\left(1-\frac{\Delta t}{m_{k}\Delta
x}|p^{1}|\right)+\frac{\Delta t}{m_{k}\Delta x}|p^{1}|=1.$

For the second-order fluxes, define
$\sigma:=\operatorname{sign}(f_{k,i}^{\ast}-f_{k,i-1}^{\ast})$. We conclude
that

$\displaystyle\phi_{i+\frac{1}{2}}(f_{k}^{\ast})\leq\begin{cases}0\quad&\text{if}\quad\sigma=-1\\\
f_{k,i+1}^{\ast}-f_{k,i}^{\ast}\quad&\text{if}\quad\sigma=+1\end{cases},$
$\displaystyle-\phi_{i-\frac{1}{2}}(f_{k}^{\ast})\leq\begin{cases}f_{k,i-1}^{\ast}-f_{k,i}^{\ast}\quad&\text{if}\quad\sigma=-1\\\
0\quad&\text{if}\quad\sigma=+1\end{cases}.$

It follows that

$\displaystyle f_{k,i}^{\ell+1}$
$\displaystyle\overset{\eqref{eq:transport_x}}{=}\left(1-\frac{\Delta
t}{m_{k}\Delta x}|p^{1}|\right)f_{k,i}^{\ast}+\frac{\Delta t}{m_{k}\Delta
x}|p^{1}|f_{k,i-\operatorname{sign}(p^{1})}^{\ast}+\frac{\Delta t}{m_{k}\Delta
x}\frac{|p^{1}|}{2}(\phi_{i+\frac{1}{2}}(f_{k}^{\ast})-\phi_{i-\frac{1}{2}}(f_{k}^{\ast}))$
$\displaystyle\leq\left(1-\frac{\Delta t}{m_{k}\Delta
x}|p^{1}|\right)f_{k,i}^{\ast}+\frac{\Delta t}{m_{k}\Delta
x}|p^{1}|f_{k,i-\operatorname{sign}(p^{1})}^{\ast}+\frac{\Delta t}{m_{k}\Delta
x}\frac{|p^{1}|}{2}\begin{cases}(f_{k,i-1}^{\ast}-f_{k,i}^{\ast})\;&\text{if}\quad\sigma=-1\\\
(f_{k,i+1}^{\ast}-f_{k,i}^{\ast})\;&\text{if}\quad\sigma=+1\end{cases}$
$\displaystyle=\left(1-\frac{3}{2}\frac{\Delta t}{m_{k}\Delta
x}|p^{1}|\right)f_{k,i}^{\ast}+\frac{\Delta t}{m_{k}\Delta
x}|p^{1}|f_{k,i-\operatorname{sign}(p^{1})}^{\ast}+\frac{\Delta t}{m_{k}\Delta
x}\frac{|p^{1}|}{2}\begin{cases}f_{k,i-1}^{\ast}\;&\text{if}\quad\sigma=-1\\\
f_{k,i+1}^{\ast}\;&\text{if}\quad\sigma=+1\end{cases}$
$\displaystyle\overset{\eqref{eq:relax-smaller-1}}{\leq}1.$

∎

######  Proposition 5.5.

If $f_{k}$ represents the distribution function of a fermion with
$f_{k}^{\ell}<1$, the time discretization in Section 5.1.2 leads to
$f_{k}^{\ell+1}<1$ for the space homogeneous case.

###### Proof.

Let $f_{k}^{\ell}<1$. The local equilibrium of a fermion is a Fermi-Dirac
distribution function $\mathcal{F}$ for which $0<\mathcal{F}<1$ by definition.
Hence,

$\displaystyle\begin{split}f_{k}^{\ell+1}&=d_{k}\left[f_{k}^{\ell}+\Delta
t(1-\gamma)(\tilde{\nu}_{kk}\mathcal{F}_{kk}^{(1)}+\tilde{\nu}_{kj}\mathcal{F}_{kj}^{(1)}-(\tilde{\nu}_{kk}+\tilde{\nu}_{kj})f_{k}^{(1)})\right]+\gamma\Delta
td_{k}(\tilde{\nu}_{kk}\mathcal{F}_{kk}^{(2)}+\tilde{\nu}_{kj}\mathcal{F}_{kj}^{(2)})\\\
&=d_{k}\left[f_{k}^{\ell}(1-2\Delta t(1-\gamma)d_{k})+\Delta
t(1-\gamma)d_{k}(\tilde{\nu}_{kk}\mathcal{F}_{kk}^{(1)}+\tilde{\nu}_{kj}\mathcal{F}_{kj}^{(1)})\right]+\gamma\Delta
td_{k}(\tilde{\nu}_{kk}\mathcal{F}_{kk}^{(2)}+\tilde{\nu}_{kj}\mathcal{F}_{kj}^{(2)})\\\
&<d_{k}\left[1-2\Delta t(1-\gamma)d_{k}+2\Delta t(1-\gamma)d_{k}+\gamma\Delta
t(\tilde{\nu}_{kk}+\tilde{\nu}_{kj})\right]=1.\end{split}$ (91)

∎

####  5.3.2 Conservation of mass, total momentum and total energy

In this section, we concern the conservation of mass, total momentum, and
total energy for the semi-discrete scheme. The proofs of the following
propositions work analogously as and can be found in the proofs of proposition
5.3 and 5.4 in [23].

######  Proposition 5.6.

The relaxation step in the first-order splitting scheme presented in Section
5.1.1 satisfies the conservation laws

$\displaystyle\int m_{1}f_{1}^{\ast}dp=\int m_{1}f_{1}^{\ell}dp,\quad\int
m_{2}f_{2}^{\ast}dp=\int m_{2}f_{2}^{\ell}dp,$ (92)
$\displaystyle\int\left(m_{1}pf_{1}^{\ast}+m_{2}pf_{2}^{\ast}\right)dp=\int\left(m_{1}pf_{1}^{\ell}+m_{2}pf_{2}^{\ell}\right)dp,$
(93)
$\displaystyle\int\left(\frac{|p|^{2}}{2m_{1}}f_{1}^{\ast}+\frac{|p|^{2}}{2m_{2}}f_{2}^{\ast}\right)dp=\int\left(\frac{|p|^{2}}{2m_{1}}f_{1}^{\ell}+\frac{|p|^{2}}{2m_{2}}f_{2}^{\ell}\right)dp.$
(94)

######  Proposition 5.7.

For each $i=1,2$, the transport step in the first-order splitting scheme in
Section 5.1.1, combined with the space discretization presented in Section 5.2
satisfies the conservation laws

$\displaystyle\sum_{i=0}^{I}\int\mathbf{p}_{k}f^{\ell+1}_{k,i}dp\Delta
x=\sum_{i=0}^{I}\int\mathbf{p}_{k}f^{\ast}_{k,i}dp\Delta x$ (95)

for periodic or zero boundary conditions.

Since the second-order time-stepping scheme in Section 5.1.2 can be broken
into relaxation and transport parts, each of which preserves the conservation
of mass, total momentum, and total energy, we can state the following:

######  Corollary 5.3.1.

For periodic or zero boundary conditions, any combination of temporal and
space discretization presented to Sections 5.1 and 5.2, respectively,
conserves mass, total momentum and total energy.

####  5.3.3 Entropy inequality

We study the entropy behavior for the first-order scheme in Section 5.1.1.
Both the relaxation and the transport step dissipate entropy; see Propositions
5.8 and 5.10. Moreover, the minimal entropy is reached for the relaxation step
if the distribution functions coincide with the corresponding target
functions; see Proposition 5.9.

######  Proposition 5.8.

Let $h_{\tau}$ given by (12). The relaxation step in the first-order splitting
scheme in Section 5.1.1 fulfills the discrete entropy inequality

$\displaystyle\int
h_{\tau}(f_{1}^{\ast})+h_{\tau^{\prime}}(f_{2}^{\ast})dp\leq\int
h_{\tau}(f_{1}^{\ell})+h_{\tau^{\prime}}(f_{2}^{\ell})dp.$ (96)

###### Proof.

By convexity

$h_{\tau_{k}}(f_{k}^{\ell})\geq
h_{\tau_{k}}(f_{k}^{\ast})+h_{\tau_{k}}^{\prime}(f_{k}^{\ast})(f_{k}^{\ell}-f_{k}^{\ast}).$
(97)

For $f\geq 0$ ($\tau\in\\{-1,0\\})$, respective $0\leq f<1$ ($\tau=+1$), the
derivative

$\displaystyle h_{\tau}^{\prime}(f)=\log\frac{1}{1-\tau f}$ (98)

is monotonically increasing such that

$\displaystyle(h_{\tau}^{\prime}(x)-h_{\tau}^{\prime}(y))(y-x)\leq 0$ (99)

for all $x,y\geq 0$ ($\tau\in\\{-1,0\\})$ and $0\leq x,y<1$ ($\tau=+1$),
respectively. Moreover, since

$\displaystyle
h_{\tau_{k}}^{\prime}(\mathcal{K}_{kk}^{\ast})=\alpha_{k}\cdot\mathbf{p}_{k},$
(100)

it holds

$\displaystyle\int
h_{\tau_{k}}^{\prime}(\mathcal{K}_{kk}^{\ast})\tilde{\nu}_{kk}(\mathcal{K}_{kk}^{\ast}-f_{k}^{\ast})dp=\int\alpha_{k}\cdot\mathbf{p}_{k}\,\tilde{\nu}_{kk}(\mathcal{K}_{kk}^{\ast}-f_{k}^{\ast})dp=0$
(101)

which vanishes as the conservation properties are satisfied at the semi-
discrete level as well by construction of the scheme. Analogously for the
inter-species terms,

$\displaystyle\int$ $\displaystyle
h_{\tau}^{\prime}(\mathcal{K}_{12}^{\ast})\tilde{\nu}_{12}(\mathcal{K}_{12}^{\ast}-f_{1}^{\ast})dp+\int
h_{\tau^{\prime}}^{\prime}(\mathcal{K}_{21}^{\ast})\tilde{\nu}_{21}(\mathcal{K}_{21}^{\ast}-f_{2}^{\ast})dp$
(102)
$\displaystyle=\alpha_{12}^{0}\int\tilde{\nu}_{12}(\mathcal{K}_{12}^{\ast}-f_{1}^{\ast})dp+\alpha_{21}^{0}\int\tilde{\nu}_{21}(\mathcal{K}_{21}^{\ast}-f_{2}^{\ast})dp$
$\displaystyle\quad+\begin{pmatrix}\alpha^{1}\\\
\alpha^{2}\end{pmatrix}\cdot\int(\tilde{\nu}_{12}(\mathcal{K}_{12}^{\ast}-f_{1}^{\ast})\begin{pmatrix}p\\\
\frac{|p|^{2}}{2m_{1}}\end{pmatrix}+\tilde{\nu}_{21}(\mathcal{K}_{21}^{\ast}-f_{2}^{\ast})\begin{pmatrix}p\\\
\frac{|p|^{2}}{2m_{2}}\end{pmatrix})dp$ $\displaystyle=0.$

The implicit step (51) is

$f_{k}^{\ast}-f_{k}^{\ell}=\Delta
t\,\tilde{\nu}_{kk}(\mathcal{K}_{kk}^{\ast}-f_{k}^{\ast})+\Delta
t\,\tilde{\nu}_{kj}(\mathcal{K}_{kj}^{\ast}-f_{k}^{\ast}).$ (103)

Using (103) and the convexity of $h_{\tau}$ leads to

$\displaystyle h_{\tau_{k}}(f_{k}^{\ast})-h(f_{k}^{\ell})$ $\displaystyle\leq
h^{\prime}(f_{k}^{\ast})(f_{k}^{\ast}-f_{k}^{\ell})$ (104)
$\displaystyle\overset{\eqref{eq:implicit_step_restated}}{=}\Delta
t\,h_{\tau_{k}}^{\prime}(f_{k}^{\ast})\tilde{\nu}_{kk}(\mathcal{K}_{kk}^{\ast}-f_{k}^{\ast})+\Delta
t\,h_{\tau}^{\prime}(f_{k}^{\ast})\tilde{\nu}_{kj}(\mathcal{K}_{kj}^{\ast}-f_{k}^{\ast}).$

Thus after integrating (104) with respect to $p$ and making use of (101) and
(102), we obtain

$\displaystyle\int h_{\tau}(f_{1}^{\ast})dp-\int h_{\tau}(f_{1}^{\ell})dp+\int
h_{\tau^{\prime}}(f_{2}^{\ast})dp-\int h_{\tau^{\prime}}(f_{2}^{\ell})dp$
(105) $\displaystyle\leq\Delta
t\int(h_{\tau}^{\prime}(f_{1}^{\ast})-h_{\tau}^{\prime}(\mathcal{K}_{11}^{\ast}))\tilde{\nu}_{11}(\mathcal{K}_{11}^{\ast}-f_{1}^{\ast})dp+\Delta
t\int(h_{\tau^{\prime}}^{\prime}(f_{2}^{\ast})-h_{\tau^{\prime}}^{\prime}(\mathcal{K}_{22}^{\ast}))\tilde{\nu}_{22}(\mathcal{K}_{22}^{\ast}-f_{2}^{\ast})dp$
$\displaystyle\,+\Delta
t\int(h_{\tau}^{\prime}(f_{1}^{\ast})-h_{\tau}^{\prime}(\mathcal{K}_{12}^{\ast}))\tilde{\nu}_{12}(\mathcal{K}_{12}^{\ast}-f_{1}^{\ast})dp+\Delta
t\int(h_{\tau^{\prime}}^{\prime}(f_{2}^{\ast})-h_{\tau^{\prime}}^{\prime}(\mathcal{K}_{21}^{\ast}))\tilde{\nu}_{21}(\mathcal{K}_{21}^{\ast}-f_{2}^{\ast})dp$
$\displaystyle\leq 0.$

The last inequality comes by (99). ∎

######  Proposition 5.9.

The inequality in Proposition 5.8 is an equality if and only if
$f_{1}^{\ell}=\mathcal{K}_{12}^{\ell}$ and
$f_{2}^{\ell}=\mathcal{K}_{21}^{\ell}$. In such cases
$f_{1}^{\ast}=\mathcal{K}_{12}^{\ast}$ and
$f_{2}^{\ast}=\mathcal{K}_{21}^{\ast}$.

###### Proof.

The proof works analogously as and can be found in [23]. ∎

######  Proposition 5.10.

Let $h_{\tau}$ by given by (12) .The transport step in the first-order
splitting scheme in Section 5.1.1 combined with the first-order spatial
discretization in Section 5.2 fulfills the discrete entropy inequality

$\displaystyle\sum_{i=0}^{I}\left\\{\int
h_{\tau}(f_{1,i}^{\ell+1})+h_{\tau^{\prime}}(f_{2,i}^{\ell+1})dp\right\\}\Delta
x\leq\sum_{i=0}^{I}\left\\{\int
h_{\tau_{1}}(f_{1,i}^{\ast})+h_{\tau^{\prime}}(f_{2,i}^{\ast})dp\right\\}\Delta
x$ (106)

for periodic or zero boundary conditions, provided that

$\displaystyle\Delta t\leq\frac{m_{k}\Delta x}{\max|p^{1}|}.$ (107)

###### Proof.

We can apply the same proof as in [23] because $h_{\tau}$ is convex. ∎

We combine the two propositions above and obtain the following:

######  Corollary 5.3.2.

For periodic or zero boundary conditions, the first-order splitting scheme
5.1.1 combined with the first-order numerical fluxes in Section 5.2 fulfills
the discrete entropy inequality

$\displaystyle\sum_{i=0}^{I}\left\\{\int
h_{\tau}(f_{1,i}^{\ell+1})+h_{\tau^{\prime}}(f_{2,i}^{\ell+1})dp\right\\}\Delta
x\leq\sum_{i=0}^{I}\left\\{\int
h(f_{1,i}^{\ast})+h(f_{2,i}^{\ast})dp\right\\}\Delta x$ (108)

provided that

$\displaystyle\Delta t\leq\frac{m_{i}\Delta x}{\max|p^{1}|}.$ (109)

###  5.4 Momentum discretization

Eventually, we discretize the momentum variable. We center the discrete
momenta $p_{q}=(p_{q_{1}}^{1},p_{q_{2}}^{2},p_{q_{3}}^{3})^{\top}$, with
$q=(q_{1},q_{2},q_{3})\in\mathbb{N}^{3}_{0}$, around $u_{\rm mix}$ with the
mixture mean velocity

$\displaystyle u_{\rm mix}=\frac{p_{1}+p_{2}}{N_{1}+N_{2}}$ (110)

and restrict them to a finite cube. This means, for each component
$r\in\\{1,2,3\\}$,

$\displaystyle p^{r}\in[m_{k}u_{{\rm
mix}}^{r}-6m_{k}v_{{\rm{th}},k},\,m_{k}u_{{\rm
mix}}^{r}+6m_{k}v_{{\rm{th}},k}]$ (111)

where $v_{{\rm{th}},k}=\sqrt{\frac{T_{\rm mix}}{m_{k}}}$ is the thermal
velocity of species $i$ and

$\displaystyle T_{\rm
mix}=\frac{n_{1}T_{1}+n_{2}T_{2}}{n_{1}+n_{2}}+\frac{1}{3}\frac{N_{1}N_{2}}{N_{1}+N_{2}}\frac{|\frac{P_{1}}{N_{1}}-\frac{P_{2}}{N_{2}}|^{2}}{n_{1}+n_{2}}$
(112)

is the mixture temperature. An adequate resolution is ensured by the momentum
mesh size $\Delta p_{k}=0.25m_{k}v_{{\rm{th}},k}$ in each direction, as in
[34].

We emphasize the advantage of the multi-species BGK model that it is possible
to use different grids for each species/equation. This feature becomes
beneficial when the species masses, and hence the thermal speeds, differ
significantly.

All momentum integrals are replaced by discrete sums using the trapezoidal
rule, i.e.

$\displaystyle\int(\cdot)dp\approx\sum_{q}\omega_{q}(\cdot)_{q}(\Delta
p_{k})^{3}$ (113)

where $\omega_{q}=\omega_{q_{1}}\omega_{q_{2}}\omega_{q_{3}}$ are the weights
and

$\displaystyle\omega_{q_{p}}=\begin{cases}1\quad\text{if
}\min(q_{p})<q_{p}<\max(q_{p}),\\\ \frac{1}{2}\quad\text{else}.\end{cases}$
(114)

We need to distinguish between discrete and continuous moments, especially
when determining the discrete local equilibria $\mathcal{K}_{kk,q}$ and
$\mathcal{K}_{kj,q}$. Since the minimization of (75) and (76) is solved using
a discrete momentum grid and discrete moments $\bar{\mu}_{k},\bar{\mu}$ as
input, the parameters $\alpha_{kk}$ and $\alpha_{kj}$ are determined such that
$\mathcal{K}_{kk,q}$ and $\mathcal{K}_{kj,q}$ have the desired discrete
moments. Thus, the conservation and entropy properties are fulfilled at the
discrete level. (A similar approach for the standard single-species BGK
equation is given in [34].)

######  Theorem 5.4.1.

Propositions 5.2, 5.3, and 5.8-5.10 all hold true after replacing continuous
integrals by their respective quadratures. Additionally, the scheme in Section
5.1.3 satisfies the following conservation properties for $\ell\geq 0$

$\displaystyle\sum_{i,q}\omega_{q}\left(f_{1,iq}^{\ell}\mathbf{p}_{1,q}(\Delta
p_{1})^{3}+f_{2,iq}^{\ell}\mathbf{p}_{2,q}(\Delta p_{2})^{3}\right)\Delta
x\sum_{i,q}\omega_{q}\left(f_{1,iq}^{0}\mathbf{p}_{1,q}(\Delta
p_{1})^{3}+f_{2,iq}^{0}\mathbf{p}_{2,q}(\Delta p_{2})^{3}\right)\Delta x$
(115)

with $\mathbf{p}_{k,q}=(1,p_{q},\frac{|p_{q}|^{2}}{2m_{k}})^{\top}$ and
$f_{k,iq}^{\ell}\approx f_{k,i}^{\ell}(p_{q})$.

##### Optimization algorithm

The minimization of (75) and (76) is solved by Newton’s method which requires
the evaluation of the gradients

$\displaystyle\nabla_{\alpha_{k}}\varphi_{k}$
$\displaystyle\approx-\sum_{q}\omega_{q}d_{k,q}\tilde{\nu}_{kk,q}\,\mathcal{K}_{kk,q}\,\mathbf{p}_{k,q}(\Delta
p_{k})^{3}+\bar{\mu}_{k}$ (116) $\displaystyle\nabla_{\alpha}\varphi$
$\displaystyle\approx-\sum_{q}\omega_{q}d_{1,q}\nu_{12,q}\,\mathcal{K}_{12,q}\,\mathbf{p}_{12,q}(\Delta
p_{1})^{3}-\sum_{q}\omega_{q}d_{2,q}\tilde{\nu}_{21,q}\,\mathcal{K}_{21,q}\,\mathbf{p}_{21,q}(\Delta
p_{2})^{3}+\bar{\mu},$ (117)

and the Hessians

$\displaystyle\nabla_{\alpha_{k}}^{2}\varphi_{k}$
$\displaystyle\approx\sum_{q}\omega_{q}d_{k,q}\tilde{\nu}_{kk,q}\,\zeta(\mathcal{K}_{kk,q},\tau_{k})\,\mathbf{p}_{k,q}\otimes\mathbf{p}_{k}(\Delta
p_{k})^{3}$ (118) $\displaystyle\nabla_{\alpha}^{2}\varphi$
$\displaystyle\approx\sum_{q}\omega_{q}d_{1,q}\tilde{\nu}_{12,q}\,\zeta(\mathcal{K}_{12,q},\tau_{1})\,\mathbf{p}_{12}\otimes\mathbf{p}_{12,q}(\Delta
p_{1})^{3}+\sum_{q}\omega_{q}d_{2,q}\tilde{\nu}_{21,q}\,\zeta(\mathcal{K}_{21,q},\tau_{2})\,\mathbf{p}_{21}\otimes\mathbf{p}_{21,q}(\Delta
p_{2})^{3}$ (119)

where $\mathbf{p}_{12,q}=(1,0,p_{1,q},\frac{|p_{1,q}|^{2}}{2m_{1}})^{\top}$,
$\mathbf{p}_{21,q}=(0,\mathbf{p}_{2,q})^{\top}$ and

$\displaystyle\zeta(g,\tau)=\begin{cases}g\quad&\text{for }\tau=0,\\\
g^{2}e^{-\alpha\cdot\mathbf{p}_{i}}\quad&\text{for }\tau=\pm 1.\end{cases}$
(120)

The input data in (73) is computed in a straight-forward way:

$\displaystyle\bar{\mu}_{k}\approx\sum_{q}\omega_{q}d_{k,q}\tilde{\nu}_{kk,q}\,G_{k,q}\,\mathbf{p}_{k,q}(\Delta
p_{k})^{3}.$ (121)

Analogously for the input data $\bar{\mu}$ in (74).

##  6 Numerical results

In this section, we present several numerical tests. We illustrate the
properties of our model and demonstrate the properties of our scheme.

###  6.1 Relaxation in a homogeneous setting

####  6.1.1 Decay rates and illustration of the schemes’ properties

We validate our numerical scheme for quantum particles and verify the decay
rates for the mean velocities and kinetic temperatures which are given
analytically in Section 3.2.

Initially, we set the distribution functions to Maxwellians

$\displaystyle f_{k}=\mathcal{M}[n_{k},U_{k},T_{k},m_{k}]=\frac{n_{k}}{(2\pi
T_{k}m_{k})^{3/2}}\exp\left(-\frac{|p-m_{k}U_{k}|^{2}}{2T_{k}m_{k}}\right)$
(122)

with

$\displaystyle m_{1}$ $\displaystyle=1.0,\quad n_{1}=1.0,\quad
U_{1}=(0.5,0,0)^{\top},\quad T_{1}=1.0,$ $\displaystyle m_{2}$
$\displaystyle=1.5,\quad n_{2}=1.2,\quad U_{2}=(0.1,0,0)^{\top},\quad
T_{2}=0.5.$

These initial data are chosen to only illustrate the basic properties of the
model and scheme, respectively, but we do not incorporate further physical
details (e.g. for a specific quantum regime). The collision frequencies are
set to $\tilde{\nu}_{kj}=1$.

For the simulation, we use a momentum grid with $48^{3}$ nodes and the first-
order splitting scheme from Section 5.1.1 with the time step $\Delta t=0.01$.

We study any combination of classical particles, fermions and bosons.
Exemplary for the interactions of fermions with fermions, we illustrate the
evolution of the entropy and the entropy dissipation in Figure 1. In Figure 2,
we demonstrate the conservation properties where the numerical oscillations in
mass, total momentum and total energy are only of the order $10^{-14}$.

In Figure 3, we verify the behavior of the mean velocities converging
exponentially fast to a common value. The numerical decay rate and the
analytical one (20) coincide very well. We only display the rate for the
interactions of fermions with fermions because the decay rate is independent
of the type of the species.

In Figure 4, we consider the behavior of the temperatures where we distinguish
between the kinetic temperatures

$\displaystyle
T_{k}=\frac{2}{3}\left(\frac{E_{k}}{N_{k}}-\frac{1}{2}\frac{|P_{k}|^{2}}{m_{k}N_{k}^{2}}\right)$
(123)

and the physical temperatures $\vartheta_{k}$ of the fluid. The latter ones
can be calculated via a nonlinear system of equations with the parameters
density, fugacity and kinetic temperature [26]. In the first column, we
observe that the kinetic temperatures do not converge to a common value
whenever a quantum particle is involved. This is also visible in the second
column. The numerical and analytical decay rates for the kinetic temperatures
coincide very well, and the difference converges to a constant value for
quantum particles. Such behavior of the kinetic temperatures for quantum
particles comes by an additional term for the decay rates (22) which vanishes
for classical-classical interactions, see Remark 5. Additionally, we compare
the results to the physical temperatures $\vartheta_{k}$. Even though the
kinetic temperatures behave differently for quantum particles, the physical
temperatures converge to a common value in all cases as predicted by the
theory.

(a) entropy

(b) entropy dissipation

Figure 1: Entropy and entropy dissipation for the test case in Section 6.1.1,
exemplary for fermion-fermion interactions. The entropy decays monotonically.

(a) species mass

(b) total momentum

(c) total energy

Figure 2: Illustration of the conservation properties for the test case in
Section 6.1.1, exemplary for fermion-fermion interactions. The mass densities
of each species ($\rho_{k}=m_{k}n_{k}$), the total momentum ($M$) and total
energy ($E$) have small oscillations of the order of $10^{-14}$.

Figure 3: Mean velocities for the test case in Section 6.1.1, exemplary for
fermion-fermion interactions. The mean velocities converge exponentially fast
to a common value, and the numerical decay rate coincides very well with the
analytical one.

fermion-fermion:

boson-fermion:

fermion-classic:

boson-boson:

boson-classic:

classic-classic:

Figure 4: Evolution of the temperatures for the test case in Section 6.1.1.
First column: kinetic temperatures $T_{k}$; whenever a quantum particle is
involved, the kinetic temperatures do not converge to a common value. Second
column: decay rates for kinetic temperatures in logarithmic scale — numerical
and analytical values coincide very well. Additionally, the difference between
the physical temperatures $\vartheta_{k}$ is displayed which decays
exponentially fast, whereas the kinetic temperatures $T_{k}$ behave
differently for quantum particles.

####  6.1.2 Sulfur-Flourine-electrons test case

We run a space homogeneous, 3-species test case inspired by [24]. In the
following, the index $S$ refers to sulfur ions, the index $F$ refers to
fluorine ions, and the index $e$ refers to electrons. For convenience, we
clarify in the Appendix how the model and the numerical scheme can be extended
straight-forwardly to more than two species.

We incorporate collision frequencies $\nu_{kj}=0.00753\frac{1}{{\rm fs}}$
which are approximately of the same order as those used in [24]. The masses of
the species are

$\displaystyle m_{S}$ $\displaystyle=32.07\text{u}-11m_{e},\quad
m_{F}=19\text{u}-7m_{e},\quad m_{e}=9.11\cdot 10^{-28}\text{g}$

with the atomic mass $\text{u}=1.6605\cdot 10^{-24}\,\text{g}$. The ions are
treated like classical particles and initialized by
$f_{k}=\mathcal{M}[n_{k},U_{k},T_{k},m_{k}]$ $(k=S,F)$ with

$\displaystyle n_{S}$ $\displaystyle=10^{19}\,\text{cm}^{-3},\quad
n_{F}=6\cdot 10^{19}\,\text{cm}^{-3},$ $\displaystyle U_{S}$
$\displaystyle=U_{F}=0\,\frac{\text{cm}}{\text{s}},$ $\displaystyle T_{S}$
$\displaystyle=T_{F}=15\,\text{eV},$

where $M$ is defined in (122). For the electrons, we compare the behavior when
they are treated like classical particles to the behavior when they are
treated like fermions. In the former case, we initialize
$f_{e}=\mathcal{M}[n_{e},U_{e},\vartheta_{e},m_{e}]$ with

$\displaystyle n_{e}=53\cdot 10^{19}\,\text{cm}^{-3},\quad
U_{e}=0\,\frac{\text{cm}}{\text{s}},\quad\vartheta_{e}=100\,\text{eV}.$

It holds $T_{e}=\vartheta_{e}$ for classical particles. In the latter case —
electrons being treated as fermions — we initialize the distribution function
by a Fermi-Dirac function, but we keep the same macroscopic quantities, i.e.

$\displaystyle f_{e}=\left[\frac{(2\pi
m_{e}\vartheta_{e})^{3/2}}{\alpha\,n_{e}}e^{\frac{|p|^{2}}{2m_{e}\vartheta_{e}}}+1\right]^{-1}$
(124)

with the scaling factor $\alpha=1.061711634$ which leads to the desired $\int
f_{e}dp=n_{e}$.

We use momentum grids with $48^{3}$ nodes for each species, and we use the
second-order IMEX RK scheme from Section 5.1.2 with time step $\Delta t=0.1$
fs.

We illustrate the evolution of the temperatures in Figure 5. For the purely
classic test case, the physical and the kinetic temperatures coincide such
that the temperature in equilibrium $T_{\rm eq}$ can be precomputed from the
initial data [23]:

$\displaystyle T_{\rm eq}=T_{\rm
mix}(0)\overset{\eqref{eq:T_mix}}{=}\frac{n_{1}T_{1}(0)+n_{2}T_{2}(0)+n_{3}T_{3}(0)}{n_{1}+n_{2}+n_{3}}.$
(125)

In Figure 5, we observe that all species temperatures converge to that value
for the classical simulation. Additionally, we display the results when we
consider the electrons to be fermions instead. As predicted by the theory, the
physical temperatures converge to a common value. However, the physical
temperatures generally differ from the kinetic temperatures in the quantum
case. As a consequence, the physical temperature in equilibrium does not equal
$T_{\rm eq}$.

Figure 5: Evolution of the physical temperatures for the Sulfur-Fluorine-
electrons quantum test case in Section 6.1.2. When both the ions and the
electrons are treated classically (lines with dots), the physical temperatures
(which coincide with the kinetic temperatures (123)) converge to the mixture
temperature $T_{\rm eq}$ defined in (125). When the electrons are treated like
fermions instead, the physical temperatures do converge to a common value as
predicted by the theory. However, this value differs from $T_{\rm eq}$.

###  6.2 Sod problem

We run a quantum-kinetic version of the well-known Sod problem [17] in the
fluid regime for fermions. As carried out in [14], the limiting equations for
the kinetic equations in the fluid regime are the quantum Euler equations.

We implement a single-species test case with the multi-species model by
assuming $m_{1}=m_{2}=m$, $n_{1}=n_{2}=n$, $U_{1}=U_{2}=U$ and
$T_{1}=T_{2}=T$. We set $m=1$ and use $\tilde{\nu}_{kj}=2\cdot 10^{4}$ for
approaching the fluid regime. The initial data is given by
$f_{1}=f_{2}=\mathcal{M}[n,u,T,m]$ where $\mathcal{M}$ is defined in (122)
with

$\displaystyle n=1,\qquad U=0,\qquad T=1,$ (126)

for $x\leq 0$ and

$\displaystyle n=0.125,\qquad U=0,\qquad T=0.8$ (127)

for $x>0$.

The simulations are run using a velocity grid with $48^{3}$ points and 300
equally spaced cells in $x$. We use the second-order IMEX Runge-Kutta scheme
from Section 5.1.2 combined with the second-order finite volume scheme from
Section 5.2.

Numerical results of the macroscopic quantities are given in Figure 6. The
fluid limit is recovered fairly well by the density $n$, mean velocity $u$ and
kinetic temperature $T$. We see again that the physical temperature
$\vartheta$ deviates from the kinetic temperature.

(a) density

(b) mean velocity

(c) temperature

Figure 6: Numerical solution at $t=0.055$ of the Sod problem in Section 6.2.
We show results for a 2-species kinetic simulation for fermions. The solutions
for both species are identical; we show only the species 1 results. For
reference, the exact solution for the quantum Euler equations is also provided
(dotted gray line). The kinetic solution recovers the fluid limit fairly well.

##  7 Appendix

The two-species model can be extended to a system of $N$-species which undergo
binary interactions. For ease in notation, we illustrate here the 3-species
case. Each distribution function $f_{k}$, $k=1,\dots,3$, represents the
solution to

$\displaystyle\partial_{t}f_{k}+\frac{p}{m_{k}}\cdot\nabla_{x}f_{k}=\tilde{\nu}_{k1}(\mathcal{K}_{k1}-f_{k})+\tilde{\nu}_{k2}(\mathcal{K}_{k2}-f_{k})+\tilde{\nu}_{k3}(\mathcal{K}_{k3}-f_{k})$
(128)

with $\tilde{\nu}_{kj}=\nu_{kj}n_{j}$. Since we still consider only binary
interactions, the properties in section 2 and 3.1 are still satisfied.

The presented numerical scheme is based on the general implicit solver in
Section 5.1.3. Since the transport operators act only on the individual
species, we focus on and shortly illustrate the scheme of the relaxation
process.

As above, we write the implicit updates of the distribution functions in a
generic steady state form

$\displaystyle f_{k}=d_{k}G_{k}+d_{k}\gamma\Delta
t(\tilde{\nu}_{kk}\mathcal{K}_{kk,\tau_{k}}+\tilde{\nu}_{kj}\mathcal{K}_{kj,\tau_{j}}+\tilde{\nu}_{kl}\mathcal{K}_{kl,\tau_{l}})$
(129)

for $k,j,l\in\\{1,2,3\\}$, each of $k,j,l$ distinct, where
$\mathcal{K}_{kk,\tau_{k}}$, $\mathcal{K}_{kj,\tau_{j}}$ and
$\mathcal{K}_{kl,\tau_{l}}$ are the unique attractors associated to $f_{k}$,

$d_{k}=\frac{1}{1+\gamma\Delta
t(\tilde{\nu}_{kk}+\tilde{\nu}_{kj}+\tilde{\nu}_{kl})},$ (130)

and $G_{k}$ is a known function. When we can express $\mathcal{K}_{kk}$,
$\mathcal{K}_{kj}$ and $\mathcal{K}_{kl}$ as functions of $G_{k}$, $G_{j}$ and
$G_{l}$, (129) provides an explicit update formula for $f_{k}$.

We apply the conservation properties to (129). An analogous calculation as in
the 2-species case leads to a set of constraints to determine the attractors
from the given data:

$\displaystyle\begin{split}&\int
d_{1}\left(\tilde{\nu}_{11}\mathcal{K}_{11,\tau_{1}}+\tilde{\nu}_{12}\mathcal{K}_{12,\tau_{1}}+\tilde{\nu}_{13}\mathcal{K}_{13,\tau_{1}}\right)\mathbf{p}_{1}dp\int
d_{2}\left(\tilde{\nu}_{21}\mathcal{K}_{21,\tau_{2}}+\tilde{\nu}_{22}\mathcal{K}_{22,\tau_{2}}+\tilde{\nu}_{23}\mathcal{K}_{23,\tau_{2}}\right)\mathbf{p}_{2}dp\\\
+&\int
d_{3}\left(\tilde{\nu}_{31}\mathcal{K}_{31,\tau_{3}}+\tilde{\nu}_{32}\mathcal{K}_{32,\tau_{3}}+\tilde{\nu}_{33}\mathcal{K}_{33,\tau_{3}}\right)\mathbf{p}_{3}dp\\\
&=\int
d_{1}\left(\tilde{\nu}_{11}+\tilde{\nu}_{12}+\tilde{\nu}_{13}\right)G_{1}\mathbf{p}_{1}dp+\int
d_{2}\left(\tilde{\nu}_{21}+\tilde{\nu}_{22}+\tilde{\nu}_{23}\right)G_{2}\mathbf{p}_{2}dp+\int
d_{3}\left(\tilde{\nu}_{31}+\tilde{\nu}_{32}+\tilde{\nu}_{33}\right)G_{3}\mathbf{p}_{3}dp.\end{split}$
(131)

These constraints (131) represent first-order optimality conditions associated
to the minimization of the convex function

$\displaystyle\varphi_{\rm
tot}(\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{12},\alpha_{13},\alpha_{23})=\varphi_{1}(\alpha_{1})+\varphi_{2}(\alpha_{2})+\varphi_{3}(\alpha_{3})+\varphi(\alpha_{12})+\varphi(\alpha_{13})+\varphi(\alpha_{23})$

with

$\displaystyle\varphi_{k}(\alpha_{k})$ $\displaystyle=\int
d_{k}\tilde{\nu}_{kk}w[\mathcal{K}_{kk,\tau_{k}}]dp+\mu_{kk}\cdot\alpha_{k}$

and

$\displaystyle\varphi(\alpha_{kj})$
$\displaystyle=\int\left(d_{k}\tilde{\nu}_{kj}w[\mathcal{K}_{kj,\tau_{k}}]+d_{j}\tilde{\nu}_{jk}w[\mathcal{K}_{jk,\tau_{j}}]\right)dp+\mu_{kj}\cdot\alpha_{kj},$

where

$\displaystyle
w[\mathcal{K}_{kj,\tau_{k}}]=\frac{\log(1-\tau_{k}\mathcal{K}_{kj,\tau_{k}})}{\tau_{k}}=\begin{cases}-\mathcal{K}_{kj,\tau_{k}}\qquad&\text{for
}\tau_{k}=0,\\\ \log(1-\mathcal{K}_{kj,1})\quad&\text{for }\tau_{k}=+1,\\\
-\log(1+\mathcal{K}_{kj,-1})\quad&\text{for }\tau_{k}=-1.\end{cases}$

Moreover, $\alpha_{k}=(\alpha_{k}^{0},\alpha_{k}^{1},\alpha_{k}^{2})^{\top}$;

$\displaystyle\mu_{kk}=\begin{pmatrix}\mu_{kk}^{0}\\\ \mu_{kk}^{1}\\\
\mu_{kk}^{2}\end{pmatrix}=\int d_{k}\tilde{\nu}_{kk}G_{k}\mathbf{p}_{k}dp$

for $k=1,2,3$; for $k\neq j:$
$\alpha_{kj}=(\alpha_{kj}^{0},\alpha_{jk}^{0},\alpha_{kj}^{1},\alpha_{kj}^{2})^{\top}$;
and

$\displaystyle\mu_{kj}=\begin{pmatrix}\mu_{kj}^{0}\\\ \mu_{jk}^{0}\\\
\mu_{kj}^{1}\\\ \mu_{kj}^{2}\end{pmatrix}$
$\displaystyle=\int\left[\begin{pmatrix}1\\\ 0\\\ p\\\
\frac{|p|^{2}}{2m_{k}}\end{pmatrix}d_{k}\nu_{kj}G_{k}+\begin{pmatrix}0\\\ 1\\\
p\\\ \frac{|p|^{2}}{2m_{j}}\end{pmatrix}d_{j}\nu_{jk}G_{j}\right]dp.$

##  8 Acknowledgements

Marlies Pirner was funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy EXC 2044-390685587,
Mathematics Münster: Dynamics–Geometry–Structure, by the Alexander von
Humboldt foundation and the German Science Foundation DFG (grant no. PI
1501/2-1).

## References

  * [1] P. Andries, K. Aoki and B. Perthame,_A consistent BGK-type model for gas mixtures,_ Journal of Statistical Physics 106 (2002) 993-1018
  * [2] Ascher, U.M., Ruuth, S.J., Spiteri, R.J.: Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Applied Numerical Mathematics 25 (1997), no. 2, 151-167
  * [3] Ashcroft, N. W.: Solid State Physic. Thomson Press, 2003.
  * [4] Bae, G. C., Yun, S. B. (2020). Quantum BGK Model near a Global Fermi–Dirac Distribution. SIAM Journal on Mathematical Analysis, 52(3), 2313-2352.
  * [5] Bae, G. C., Klingenberg, C., Pirner, M., Yun, S. B. (2019). BGK model of the multi-species Uehling Uhlenbeck equation. arXiv preprint arXiv:1912.01677.
  * [6] Bae, G.-C., Yun, S.-B.: Stationary quantum BGK model for bosons and fermions in a bounded interval. J. Stat. Phys. 178 (2020), no. 4, 845-868.
  * [7] Bae, G. C., Klingenberg, C., Pirner, M., Yun, S. B.: BGK model of the multi-species Uehling-Uhlenbeck equation. Kinet. Relat. Models (2021), no.1, 25–44.
  * [8] Bae, G. C., Klingenberg, C., Pirner, M., Yun, S. B.: BGK model for two-component gases near a global Maxwellian.SIAM J. Math. Anal. 55 (2023), no.2, 1007–1047.
  * [9] Bae, G. C., Klingenberg, C., Pirner, M., Yun, S. B.: BGK model for two-component gases near a global Maxwellian. SIAM Journal of Mathematical Analysis 55 (2), 2023
  * [10] Bobylev, A. V., Bisi, M., Groppi, M., Spiga, G., Potapenko, I. F. (2018). A general consistent BGK model for gas mixtures. Kinetic and Related Models, 11(6).
  * [11] S. Brull, V. Pavan and J. Schneider, Derivation of a BGK model for mixtures, European Journal of Mechanics B/Fluids, 33 (2012) 74-86
  * [12] Escobedo, M., Mischler, S., Valle, M. A.: Entropy maximisation problem for quantum relativistic particles. Bull. Soc. Math. France. 133 (2005) no. 1, 87-120.
  * [13] Feng, D., Jin, G.: Introduction to Condensed Matter Physics, World Scientific Publishing Company, 2005.
  * [14] Francis Filbet, Jingwei Hu, Shi Jin (2012): A numerical scheme for the quantum Boltzmann equation with stiff collision terms. ESAIM: M2AN 46 (2) 443-463, DOI: 10.1051/m2an/2011051
  * [15] Gamba, Irene M and Haack, Jeffrey R and Hauck, Cory D and Hu, Jingwei, _A fast spectral method for the Boltzmann collision operator with general collision kernels_ , SIAM Journal on Scientific Computing 39 (4), pp. B658–B674 (2017)
  * [16] _Gamba, Irene M and Tharkabhushanam, Sri Harsha, Spectral-Lagrangian methods for collisional models of non-equilibrium statistical states_ , Journal of Computational Physics 228 (6), pp. 2012–2036 (2009)
  * [17] Gary A. Sod: A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. Journal of Computational Physics 27 (1978), no. 1, 1–31
  * [18] V. Garzó, A. Santos and J. J. Brey, A kinetic model for a multicomponent gas Physics of Fluids, 1 (1989) 380-383
  * [19] J. Greene, Improved Bhatnagar-Gross-Krook model of electron-ion collisions. Phys. Fluids 16, 2022– 2023 (1973)
  * [20] E. P. Gross and M. Krook, Model for collision processes in gases: small-amplitude oscillations of charged two-component systems, Physical Review 3 (1956) 593
  * [21] J. R. Haack, C.D. Haack, and M.S.Murillo . A conservative, entropic multispecies BGK model. Journal of Statistical Physics, 168 (2017), 826-856.
  * [22] Haack, J., Hauck, C., Klingenberg, C., Pirner, M., Warnecke, S.: A consistent BGK model with velocity-dependent collision frequency for gas mixtures. J. Stat. Phys. 184 (2021), no. 31.
  * [23] Haack, J., Hauck, C., Klingenberg, C., Pirner, M., Warnecke, S. (2022). Numerical schemes for multi-species BGK equations with velocity-dependent collision frequency. submitted. arXiv preprint arXiv:2202.05652.
  * [24] Haack, J., Hauck, C., Murillo, M.: A conservative, entropic multispecies BGK model. J. Stat. Phys. 168 (2017), no. 4, 826–856
  * [25] B. Hamel, Kinetic model for binary gas mixtures, Physics of Fluids 8 (1965) 418-425
  * [26] Hu, J., Jin, S. (2011). On kinetic flux vector splitting schemes for quantum Euler equations. Kinetic and Related Models.
  * [27] Jüngel, A.: Transport equations for semiconductors. Lecture Notes in Physics, 773. Springer-Verlag, Berlin, 2009.
  * [28] C. Klingenberg, M.Pirner, G.Puppo, _A consistent kinetic model for a two-component mixture with an application to plasma_ , Kinetic and related Models 10 (2017) 445-465
  * [29] Lu, X.: A Modified Boltzmann Equation for Bose-Einstein Particles: Isotropic Solutions and Long-Time Behavior. J. Statist. Phys. 98 (2000), no. 5-6, 1335–1394.
  * [30] Lu, X.: On Spatially Homogeneous Solutions of a Modified Boltzmann Equation for Fermi-Dirac Particles. J. Statist. Phys. 105 (2001), no. 1-2, 353-388.
  * [31] Madelung, O.: Introduction to solid-state theory. Vol. 2. Springer Science and Business Media, 2012.
  * [32] Markowich, P. A., Ringhofer, C. A., Schmeiser, C. : Semiconductor equations. Springer-Verlag, Vienna, 1990.
  * [33] Mieussens, Luc and Struchtrup, Henning. Numerical comparison of Bhatnagar–Gross–Krook models with proper Prandtl number. Physics of Fluids. 16. 10.1063/1.1758217 (2004)
  * [34] Mieussens, Luc. _Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamics_ ; Mathematical Models and Methods in Applied Sciences, Vol. 10, No. 08 (2000): 1121-1149
  * [35] Mouhot, Clément and Pareschi, Lorenzo, _Fast algorithms for computing the Boltzmann collision operator_ , Mathematics of computation 75, 256, pp. 1833–1852 (2006)
  * [36] Pareschi, Lorenzo and Russo, Giovanni, _Numerical solution of the Boltzmann equation I: Spectrally accurate approximation of the collision operator_ , SIAM journal on numerical analysis 37, 4, pp. 1217–1245 (2000)
  * [37] Reinhard, P. G., Suraud, E.: A quantum relaxation-time approximation for finite fermion systems. Annals of Physics 354 (2015), 183–202.
  * [38] V. Sofonea and R. Sekerka, BGK models for diffusion in isothermal binary fluid systems, Physica, 3 (2001), 494-520
  * [39] Uehling, E. A., Uhlenbeck, G. E.: Transport phenomena in einstein-bose and fermi-dirac gases. i. Physical Review, 43 (1933), no. 7, 552.
  * [40] Uehling, E. A.: Transport phenomena in einstein-bose and fermi-dirac gases. ii. Physical Review, 46 (1934), no. 10, 917.