∎

11institutetext: Erica R. Johnson 22institutetext: Southwest Research
Institute, Space Science Department,
6220 Culebra Road, San Antonio, Texas 78238, USA
22email<EMAIL_ADDRESS>33institutetext: James A. Rossmanith
44institutetext: Iowa State University, Department of Mathematics,
411 Morrill Road, Ames, Iowa 50011, USA
44email<EMAIL_ADDRESS>55institutetext: Christine Vaughan
66institutetext: Physicians Health Plan
1400 East Michigan Ave, Lansing, Michigan 48912, USA
66email<EMAIL_ADDRESS>
# Positivity-Preserving Lax-Wendroff Discontinuous Galerkin Schemes for
Quadrature-Based Moment-Closure Approximations of Kinetic Models††thanks: This
research was partially funded by NSF Grants DMS–1620128 and DMS–2012699.

Erica R. Johnson  James A. Rossmanith111Corresponding author Christine Vaughan

###### Abstract

The quadrature-based method of moments (QMOM) offers a promising class of
approximation techniques for reducing kinetic equations to fluid equations
that are valid beyond thermodynamic equilibrium. In this work, we study a
particular five-moment variant of QMOM known as HyQMOM and establish that this
system is moment-invertible over a convex region in solution space. We then
develop a high-order discontinuous Galerkin (DG) scheme for solving the
resulting fluid system. The scheme is based on a predictor-corrector approach,
where the prediction is a localized space-time DG scheme. The nonlinear
algebraic system in this prediction is solved using a Picard iteration. The
correction is a straightforward explicit update based on the time-integral of
the evolution equation, where the space-time prediction replaces all instances
of the exact solution. In the absence of limiters, the high-order scheme does
not guarantee that solutions remain in the convex set over which HyQMOM is
moment-realizable. To overcome this, we introduce novel limiters that
rigorously guarantee that the computed solution does not leave the convex set
of realizable solutions, thus guaranteeing the hyperbolicity of the system. We
develop positivity-preserving limiters in both the prediction and correction
steps and an oscillation limiter that damps unphysical oscillations near
shocks. We also develop a novel extension of this scheme to include a BGK
collision operator; the proposed method is shown to be asymptotic-preserving
in the high-collision limit. The HyQMOM and the HyQMOM-BGK solvers are
verified on several test cases, demonstrating high-order accuracy on smooth
problems and shock-capturing capability on problems with shocks. The
asymptotic-preserving property of the HyQMOM-BGK solver is also numerically
verified.

###### Keywords:

discontinuous Galerkin hyperbolic conservation laws moment closure positivity-
preserving limiters

###### MSC:

65M60 82C80

††journal:

###### Contents

  1. 1 Introduction
    1. 1.1 Fluid models and the moment-closure problem
    2. 1.2 Moment-closure methods
    3. 1.3 Scope of this work
  2. 2 Brief review of moment-realizability
    1. 2.1 Existence: Moment-realizability
    2. 2.2 Example: $S=4$ case
    3. 2.3 Techniques for moment-closure
  3. 3 HyQMOM: Hyperbolic quadrature-based moment closure
    1. 3.1 Classical QMOM Approach
      1. 3.1.1 Example: $N=2$ case
      2. 3.1.2 Weak hyperbolicity and linear degeneracy of QMOM for all $N\geq 1$
    2. 3.2 Pressure Regularized QMOM
    3. 3.3 HyQMOM: Density Regularized QMOM
  4. 4 Locally-implicit Lax-Wendroff discontinuous Galerkin
    1. 4.1 Discontinuous Galerkin finite elements
    2. 4.2 Prediction step
    3. 4.3 Correction step
  5. 5 HyQMOM limiters
    1. 5.1 Positivity of the Rusanov scheme
    2. 5.2 Limiter I: Positivity-at-points in the prediction step
    3. 5.3 Limiter II: Positivity-in-the-mean in the correction step
    4. 5.4 Limiter III: Positivity-at-points in the correction step
    5. 5.5 Limiter IV: Unphysical oscillation limiter
  6. 6 Collisionless HyQMOM numerical examples
    1. 6.1 Smooth solution convergence test
    2. 6.2 Shock tube problem #1
    3. 6.3 Shock tube problem #2
    4. 6.4 Double rarefaction vacuum problem
  7. 7 Extension to HyQMOM-BGK
    1. 7.1 1D1V Boltzmann-BGK equation
    2. 7.2 HyQMOM-BGK and the asymptotic-preserving property
    3. 7.3 Prediction
    4. 7.4 Post-prediction BGK source evaluation
    5. 7.5 Correction
    6. 7.6 Asymptotic-preserving condition
  8. 8 HyQMOM-BGK numerical examples
    1. 8.1 Manufactured solution convergence test
    2. 8.2 BGK shock tube problem
  9. 9 Conclusions
  10. A Appendix

##  1 Introduction

Kinetic Boltzmann equations model the non-equilibrium dynamics of a wide
variety of fluids, including gases, multiphase flows, and plasma. These
equations have the following general form:

$f_{,t}+{\underline{v}}\cdot{\underline{\nabla}}_{{\underline{x}}}f+{\underline{\mathcal{F}}}\cdot{\underline{\nabla}}_{{\underline{v}}}f=\mathbb{C}(f),$
(1)

where $f\left(t,{\underline{x}},{\underline{v}}\right):\mathbb{R}_{\geq
0}\times\mathbb{R}^{D}\times\mathbb{R}^{V}\mapsto\mathbb{R}_{\geq 0}$ is the
distribution function that describes the state of the fluid,
$t\in\mathbb{R}_{\geq 0}$ is time, ${\underline{x}}\in\mathbb{R}^{D}$ is the
spatial coordinate, and ${\underline{v}}\in\mathbb{R}^{V}$ is the velocity
coordinate. Additionally,
${\underline{\mathcal{F}}}\left(t,{\underline{x}},{\underline{v}}\right)\in\mathbb{R}_{\geq
0}\times\mathbb{R}^{D}\times\mathbb{R}^{V}\mapsto\mathbb{R}^{V}$ is the
forcing term that could include lift, drag, gravity, and other forces acting
on the particles, and $\mathbb{C}(f):\mathbb{R}_{\geq 0}\mapsto\mathbb{R}$ is
the collision term that describes direct particle-particle interactions.

Kinetic models of the form (1) offer two desirable features: (1) the evolution
equations have a relatively simple form (i.e., advection in phase space); and
(2) the models are capable of accurately describing a large class of physical
phenomena that are important in application problems. However, the main
difficulty with kinetic models is that their solutions live in high-
dimensional phase space, which means that high fidelity numerical computations
are very expensive.

###  1.1 Fluid models and the moment-closure problem

One approach for reducing the computational complexity of kinetic models is to
replace them with so-called fluid models, which means that instead of evolving
the distribution function directly, one evolves a finite set of moments of the
distribution function. For example the $(\ell_{1},\ell_{2},\ell_{3})$ moment
of the distribution function, $f$, is defined as follows:

${\mathsf{M}}_{(\ell_{1},\ell_{2},\ell_{3})}:=\int_{\mathbb{R}^{3}}v_{1}^{\ell_{1}}v_{2}^{\ell_{2}}v_{3}^{\ell_{3}}f\left(t,{\underline{x}},{\underline{v}}\right)\,dv_{1}\,dv_{2}\,dv_{3},$
(2)

where ${\underline{v}}=(v_{1},v_{2},v_{3})$. If moments of the kinetic
equation (1) are computed, we arrive at three-dimensional evolution equations
of the following form:

${\mathsf{M}}_{(\ell_{1},\ell_{2},\ell_{3}),t}+{\mathsf{M}}_{(\ell_{1}+1,\ell_{2},\ell_{3}),x_{1}}+{\mathsf{M}}_{(\ell_{1},\ell_{2}+1,\ell_{3}),x_{2}}+{\mathsf{M}}_{(\ell_{1},\ell_{2},\ell_{3}+1),x_{3}}=0,$
(3)

where for simplicity, we have set the forcing and collision operator to zero:
${\underline{\mathcal{F}}}\equiv{\underline{0}}$ and $\mathbb{C}\equiv 0$.

The key benefit of considering a finite set of evolution equations of the form
(3) over the fully kinetic equation (1) is the reduction in the number of
independent variables from $1+D+V\leq 7$ to $1+D\leq 4$. However, we observe
from (3) the key challenge in fluid model approximations of the kinetic
equation: to evolve the moment ${\mathsf{M}}_{(\ell_{1},\ell_{2},\ell_{3})}$,
we need to know higher-order moments:
${\mathsf{M}}_{(\ell_{1}+1,\ell_{2},\ell_{3})}$,
${\mathsf{M}}_{(\ell_{1},\ell_{2}+1,\ell_{3})}$, and
${\mathsf{M}}_{(\ell_{1}+1,\ell_{2},\ell_{3}+1)}$. This issue is known simply
as the moment-closure problem. And in particular, to obtain a closed fluid
system, one needs to somehow approximate the highest moments in the system.
Furthermore, different choices lead to systems of differential (or integro-
differential) equations with vastly different mathematical properties.

###  1.2 Moment-closure methods

A standard approach for developing moment-closure approximation for (1) is to
assume a specific ansatz for the distribution function:

$f(t,{\underline{x}},{\underline{v}})\approx\sum_{\ell=1}^{M}\psi_{\ell}\left({\underline{v}},\beta_{\ell}\left(t,{\underline{x}}\right)\right)\quad\text{or}\quad
f(t,{\underline{x}},{\underline{v}})\approx\prod_{\ell=1}^{M}\psi_{\ell}\left({\underline{v}},\beta_{\ell}\left(t,{\underline{x}}\right)\right).$
(4)

The first systematic attempt at developing a moment-closure approach is due to
the seminal work of Grad article:grad49 , in which he proposed a moment-
closure that assumed the distribution was a Maxwell-Boltzmann distribution
multiplied by a polynomial in ${\underline{v}}$. Since Grad’s work, a vast
body of literature has developed on various moment closures, including modern
modifications of Grad’s closure (e.g., article:CaFaLi13 ; article:CaFaLi14 ;
article:Koellermeier2021 ; article:Koellermeier2014 ; article:KoTo17 ),
maximum entropy article:Dreyer87 ; article:Lev96 ; article:MuRu93  and its
numerous variants (e.g., article:AbdBru2016 ; article:BohmTorr2019 ), and
quadrature-based moment closures (e.g., article:desjardins2008 ; article:Fox08
; article:fox09 ; article:FoxLauVie2018 ; article:marchisio2005 ;
article:PaDeFo2019 ). A full review of all methods is well beyond the scope of
the current work, but can be found in the paper of Torrilhon
article:Torrilhon2009  and the references therein.

###  1.3 Scope of this work

In this work, we study a particular five-moment moment-closure known as the
hyperbolic quadrature-based moment closure (HyQMOM), originally due to Fox,
Laurent, and Vié article:FoxLauVie2018 . Our focus here is only on the one-
dimensional version of kinetic equation (1) (i.e., 1D1V). We begin in §2 with
a brief review of the moment-closure problem and a few strategies for
producing fluid approximations with desirable mathematical properties. In §3
we provide a brief review of the classical quadrature-based moment closure
(QMOM), show its shortcomings, and then establish that the HyQMOM system is
moment-invertible over a convex set in solution space. In §4 we introduce a
novel high-order Lax-Wendroff discontinuous Galerkin scheme for solving the
HyQMOM fluid system. The scheme is based on a predictor-corrector approach,
where the prediction step is based on a localized space-time discontinuous
Galerkin scheme. The nonlinear algebraic system that arises in this prediction
step is solved using a Picard iteration. The correction is a straightforward
explicit update based on the time-integral of the evolution equation, where
the space-time prediction replaces all instances of the exact solution. In the
absence of additional limiters, the proposed high-order scheme does not
guarantee that the numerical solution remains in the convex set over which
HyQMOM is moment-realizable. To overcome this challenge, we introduce in §5
novel limiters that rigorously guarantee that the computed solution does not
leave the convex set over which moment-invertible and hyperbolicity of the
fluid system is guaranteed. We develop positivity-preserving limiters in both
the prediction and correction steps and an oscillation limiter that damps
unphysical oscillations near shocks. In §6 we perform convergence tests to
verify the order of accuracy of the scheme and test the scheme on Riemann data
to demonstrate the shock-capturing and robustness of the method. In §7 we
develop an asymptotic-preserving article:Jin12  extension of the proposed
scheme that allows us to solve a five-moment fluid model with a Bhatnagar-
Gross-Krook (BGK) article:BGK54  collision operator. Finally, in §8 we perform
convergence tests to verify the order of accuracy of the scheme and verify the
method on Riemann data with different Knudsen numbers. Conclusions are
provided in §9.

##  2 Brief review of moment-realizability

For the remainder of the present work, we focus exclusively on the one-
dimensional version of (1). In particular, in this section, we focus on the
transport portion of the kinetic equation, and thus restrict ourselves to the
1D1V collisionless Boltzmann equation (aka Vlasov equation):

$f_{,t}+vf_{,x}=0,$ (5)

where $f(t,x,v):\mathbb{R}_{\geq
0}\times\mathbb{R}\times\mathbb{R}\mapsto\mathbb{R}_{\geq 0}$ is the
probability distribution function that describes the state of the fluid.

The moments of $f$ are defined as follows:

${\mathsf{M}}_{\ell}:=\int_{\mathbb{R}}v^{\ell}\,f\,dv,\quad\text{for}\quad\ell\in{\mathbb{N}}_{0}.$
(6)

A simple calculation reveals that for each $\ell\in{\mathbb{N}}_{0}$, the
moments satisfy the following equation:

${\mathsf{M}}_{\ell,t}+{\mathsf{M}}_{\ell+1,x}=0.$ (7)

The key difficulty is that the evolution of the $\ell^{\text{th}}$ moment
depends on the $(\ell+1)^{\text{st}}$ moment, meaning that the moment
expansion does not produce a closed system. Therefore, the key challenge for
developing fluid approximations of kinetic models is the following question:
How does one close the moment hierarchy?

######  Definition 1 (Univariate moment-closure problem)

Let $S\in{\mathbb{N}}_{0}$. Given only the first $S+1$ moments of a univariate
distribution function $f(v)$:

${\mathsf{M}}_{\ell}=\int_{-\infty}^{\infty}v^{\ell}f(v)\,dv,\quad\text{for}\quad\ell=0,1,\ldots,S,$
(8)

find an approximation of the next moment, ${\mathsf{M}}_{S+1}$, in terms of
the given moments.

The basic strategy in most moment-closure approaches is as follows: (1) start
with a finite set of moments (e.g., $\ell=0,1,\ldots,S$); (2) assume a form of
the distribution function with several free parameters (typically, the number
of free parameters is the same as the number of moments that will be tracked);
(3) determine the free parameters in the assumed distribution function so that
its moments match all of the known moments (this part is called moment-
inversion); and (4) compute the next moment of the assumed distribution
function, which is then used to provide the flux in the evolution equation for
${\mathsf{M}}_{S}$:

${\mathsf{M}}_{S,t}+\mkern
1.5mu\overline{\mkern-1.5mu{\mathsf{M}}\mkern-1.5mu}\mkern 1.5mu_{S+1,x}=0,$
(9)

where $\mkern 1.5mu\overline{\mkern-1.5mu{\mathsf{M}}\mkern-1.5mu}\mkern
1.5mu_{S+1}=\mkern 1.5mu\overline{\mkern-1.5mu{\mathsf{M}}\mkern-1.5mu}\mkern
1.5mu_{S+1}\left({\mathsf{M}}_{0},{\mathsf{M}}_{1},\ldots,{\mathsf{M}}_{S}\right)$.

###  2.1 Existence: Moment-realizability

Before considering specific strategies for approximating the missing moment,
${\mathsf{M}}_{S+1}$, it is worthwhile to discuss the general existence
problem first. We begin by defining some important quantities relevant
throughout this work, namely the mass density, macroscopic velocity, pressure,
heat flux, and modified kurtosis:

$\begin{gathered}\rho:=\int_{-\infty}^{\infty}f\,dv,\qquad
u:=\frac{1}{\rho}\int_{-\infty}^{\infty}vf\,dv,\qquad
p:=\int_{-\infty}^{\infty}(v-u)^{2}f\,dv,\\\
h:=\int_{-\infty}^{\infty}(v-u)^{3}f\,dv,\qquad\text{and}\qquad
k:=\int_{-\infty}^{\infty}(v-u)^{4}f\,dv-\left(\frac{p^{3}+\rho
q^{2}}{\rho\,p}\right),\end{gathered}$ (10)

where we assume that $\rho,p>0$. These primitive variables are directly linked
to the moments ${\mathsf{M}}_{\ell}$:

$\begin{gathered}\rho={\mathsf{M}}_{0},\qquad
u=\frac{{\mathsf{M}}_{1}}{{\mathsf{M}}_{0}},\qquad
p={\mathsf{M}}_{2}-\frac{{\mathsf{M}}_{1}^{2}}{{\mathsf{M}}_{0}},\qquad
h={\mathsf{M}}_{3}-\frac{3{\mathsf{M}}_{1}{\mathsf{M}}_{2}}{{\mathsf{M}}_{0}}+\frac{2{\mathsf{M}}_{1}^{3}}{{\mathsf{M}}_{0}^{2}},\\\
\text{and}\qquad
k=\frac{{\mathsf{M}}_{2}^{3}-2{\mathsf{M}}_{1}{\mathsf{M}}_{2}{\mathsf{M}}_{3}+{\mathsf{M}}_{0}{\mathsf{M}}_{3}^{2}+{\mathsf{M}}_{1}^{2}{\mathsf{M}}_{4}-{\mathsf{M}}_{0}{\mathsf{M}}_{2}{\mathsf{M}}_{4}}{{\mathsf{M}}_{1}^{2}-{\mathsf{M}}_{0}{\mathsf{M}}_{2}}.\end{gathered}$
(11)

Using these we define the normalized velocity variable, $s$, and the
normalized moments: $\widetilde{\mathsf{M}}_{j}$:

$s:=\frac{v-u}{\sqrt{T}}\qquad\text{and}\qquad\widetilde{\mathsf{M}}_{j}:=\frac{\sqrt{T}}{\rho}\int_{-\infty}^{\infty}s^{j}\,f\left(v(s)\right)\,ds=\int_{-\infty}^{\infty}s^{j}\,\widetilde{f}\left(s\right)\,ds,$
(12)

where $T=p/\rho$ is the temperature. The moments and the normalized moments
are related as follows:

$\displaystyle{{\mathsf{M}}}_{\ell}$
$\displaystyle=\rho\,\sum_{j=0}^{\ell}\begin{pmatrix}\ell\\\
j\end{pmatrix}T^{\frac{j}{2}}\,u^{\ell-j}\,\widetilde{\mathsf{M}}_{j}\qquad\text{and}\qquad\widetilde{\mathsf{M}}_{j}={\rho}^{-1}\,T^{-\frac{j}{2}}\sum_{\ell=0}^{j}\begin{pmatrix}j\\\
\ell\end{pmatrix}\left(-u\right)^{j-\ell}\,{{\mathsf{M}}}_{\ell}.$ (13)

By construction, the normalized moments have the following property:

$\widetilde{\mathsf{M}}_{0}=1,\qquad\widetilde{\mathsf{M}}_{1}=0,\qquad\text{and}\qquad\widetilde{\mathsf{M}}_{2}=1.$
(14)

######  Definition 2 (Realizable moments)

The following rescaled moments:

$\widetilde{\mathsf{M}}_{0}=1,\quad\widetilde{\mathsf{M}}_{1}=0,\quad\widetilde{\mathsf{M}}_{2}=1,\quad\widetilde{\mathsf{M}}_{3},\quad\ldots,\quad\widetilde{\mathsf{M}}_{S},$

where $S\in{\mathbb{Z}}_{\geq 3}$ and
$\bigl{|}\widetilde{\mathsf{M}}_{j}\bigr{|}<\infty$ $\forall
j\in{\mathbb{Z}}_{\geq 0}$, are called realizable if there exists a
probability density function,
$\widetilde{f}(s):\mathbb{R}\mapsto\mathbb{R}_{\geq 0}$, such that

$\widetilde{\mathsf{M}}_{j}=\int_{-\infty}^{\infty}s^{j}\,\widetilde{f}(s)\,ds\qquad\text{for}\quad
j=0,1,\ldots,S.$

This leads to an obvious question: for a given $S\in{\mathbb{Z}}_{\geq 3}$,
under what conditions is the set of moments,
$\\{\widetilde{\mathsf{M}}_{0}=1,\,\widetilde{\mathsf{M}}_{1}=0,\,\widetilde{\mathsf{M}}_{2}=1,\,\widetilde{\mathsf{M}}_{3},\ldots,\,\widetilde{\mathsf{M}}_{S}\\}$,
realizable? This question is the celebrated truncated Hamburger moment problem
(e.g., see Chapter 9 of book:schmudgen2017 ) for which we can state the
following result. Note that we state only the case where $S$ is even, although
a similar result also exists when $S$ is odd book:schmudgen2017 .

######  Theorem 2.1 (Truncated Hamburger moment problem (adapted from Theorem
9.27 of book:schmudgen2017 ))

Let $S\in{\mathbb{Z}}_{>3}$ be an even integer. The set of moments:

$\\{\widetilde{\mathsf{M}}_{0}=1,\,\widetilde{\mathsf{M}}_{1}=0,\,\widetilde{\mathsf{M}}_{2}=1,\,\widetilde{\mathsf{M}}_{3},\ldots,\,\widetilde{\mathsf{M}}_{S}\\}$

is realizable if all of the Hankel determinants for $m=0,1,\ldots,S/2$ are
positive:

$D_{m}:=\begin{vmatrix}1&0&1&\widetilde{\mathsf{M}}_{3}&\cdots&\widetilde{\mathsf{M}}_{m}\\\
0&1&\widetilde{\mathsf{M}}_{3}&\widetilde{\mathsf{M}}_{4}&\cdots&\widetilde{\mathsf{M}}_{m+1}\\\
1&\widetilde{\mathsf{M}}_{3}&\widetilde{\mathsf{M}}_{4}&\widetilde{\mathsf{M}}_{5}&\cdots&\widetilde{\mathsf{M}}_{m+2}\\\
\widetilde{\mathsf{M}}_{3}&\widetilde{\mathsf{M}}_{4}&\widetilde{\mathsf{M}}_{5}&\widetilde{\mathsf{M}}_{6}&\cdots&\widetilde{\mathsf{M}}_{m+3}\\\
\vdots&\vdots&\vdots&\vdots&&\vdots\\\
\widetilde{\mathsf{M}}_{m}&\widetilde{\mathsf{M}}_{m+1}&\widetilde{\mathsf{M}}_{m+2}&\widetilde{\mathsf{M}}_{m+3}&\cdots&\widetilde{\mathsf{M}}_{2m}\end{vmatrix}>0.$

###  2.2 Example: $S=4$ case

As an example, which will become relevant later in this work, consider the
case $S=4$, where the three relevant Hankel determinants are

$D_{0}=1,\quad D_{1}=1,\quad D_{2}=\begin{vmatrix}1&0&1\\\
0&1&\widetilde{\mathsf{M}}_{3}\\\
1&\widetilde{\mathsf{M}}_{3}&\widetilde{\mathsf{M}}_{4}\end{vmatrix}=\widetilde{\mathsf{M}}_{4}-\widetilde{\mathsf{M}}_{3}^{2}-1>0.$
(15)

Thus, the realizability condition in the univariate $S=4$ case is

$\widetilde{\mathsf{M}}_{4}>\widetilde{\mathsf{M}}_{3}^{2}+1,$ (16)

which is depicted in Figure 1. Also shown in this figure is the location of
the distribution in thermodynamic equilibrium (i.e., the Maxwellian
distribution):

$f(v)=\frac{\rho}{\sqrt{2\pi
T}}e^{-\frac{(v-u)^{2}}{2T}}\quad\Longrightarrow\quad\widetilde{f}(s)=\frac{1}{\sqrt{2\pi}}e^{-\frac{s^{2}}{2}}\quad\Longrightarrow\quad\begin{cases}\widetilde{\mathsf{M}}_{3}=0\\\
\widetilde{\mathsf{M}}_{4}=3\end{cases}.$ (17)

Finally, we note that the realizability condition (16) in terms of the
primitive variables can be written as follows:

$\frac{\rho\,k}{p^{2}}>0,$ (18)

which is satisfied if $\rho>0$, $p>0$, and $k>0$.

$-2$$-1$$0$$1$$2$$0$$1$$2$$3$$4$$-1$$0$$1$$1$$2$$3$Maxwellian$\widetilde{\mathbb{M}}_{3}$$\widetilde{\mathbb{M}}_{4}$Region
of Moment-Realizability for $S=4$ Figure 1: Region of the moment-realizability
when $S=4$. In other words, given the five rescaled moments:
$\widetilde{\mathbb{M}}_{0}=1$, $\widetilde{\mathbb{M}}_{1}=0$,
$\widetilde{\mathbb{M}}_{2}=1$, $\widetilde{\mathbb{M}}_{3}$, and
$\widetilde{\mathbb{M}}_{4}$, there exists a positive distribution function
matching the given moments provided that these moments lie within the shaded
pink region of the above graph. Thermodynamic equilibrium (i.e., the
Maxwellian distribution) occurs at $\widetilde{\mathbb{M}}_{3}=0$ and
$\widetilde{\mathbb{M}}_{4}=3$, which is in the interior of the shaded pink
region.

###  2.3 Techniques for moment-closure

There are many different moment closure methods for approximating the final
moment, $\mkern 1.5mu\overline{\mkern-1.5mu{\mathsf{M}}\mkern-1.5mu}\mkern
1.5mu_{S+1}$, and each method has its own merits and challenges. A full review
of all methods is well beyond the scope of the current work, but can be found
in the paper of Torrilhon article:Torrilhon2009  and the references therein.
In this work, we settle for a brief summary of three broad classes of the most
commonly used closures.

Grad closure:

This approach was originally developed by Grad article:grad49  in 1949, but
variants with improved hyperbolicity properties have been introduced more
recently article:CaFaLi13 ; article:CaFaLi14 ; article:Koellermeier2021 ;
article:Koellermeier2014 ; article:KoTo17 . The basic idea is that the
distribution function is approximated as a Maxwellian multiplied by a
polynomial in $v$:

$\widetilde{f}(s)\sim\frac{e^{-\frac{s^{2}}{2}}}{\sqrt{2\pi
T}}\left[\rho+\sum_{\ell=3}^{S}T^{-\frac{\ell}{2}}\beta_{\ell}\psi_{\ell}\left(s\right)\right],$
(19)

where $\psi_{\ell}\left(s\right)$ is the Hermite polynomial of degree $\ell$
and ${\underline{\beta}}$ are coefficients chosen so that $\widetilde{f}(s)$
matches the first $S+1$ input moments. A related approach widely used in
applications is the R13 model, which regularizes the 13-moment Grad closure
through additional terms from the Chapman-Enskog expansion; an excellent
review of the R13 model can be found in Torrilhon article:Torrilhon2016 .

Maximum entropy closure:

The maximum entropy closure article:Dreyer87 ; article:Lev96 ; article:MuRu93
and its numerous variants (e.g., see article:AbdBru2016 ; article:BohmTorr2019
) formulate the moment-inversion problem as an optimization problem to
maximize the entropy under some assumed form of the distribution function. In
the original formulation, the distribution function is approximated as an
exponential of a polynomial:

$\widetilde{f}(s)\sim
e^{\,{\underline{\beta}}\cdot{\underline{{\Phi}}}(s)}\qquad\text{where}\qquad{\underline{\Phi}}(s)=[1,s,s^{2},\dots]^{T},$
(20)

and the coefficients, ${\underline{\beta}}$, are chosen so that
$\widetilde{f}(s)$ matches the first $S+1$ input moments.

Quadrature-based moment closures:

In the quadrature-based method of moments (QMOM) (e.g., article:desjardins2008
; article:Fox08 ; article:fox09 ; article:FoxLauVie2018 ;
article:marchisio2005 ; article:PaDeFo2019 ), the distribution function is
represented as a sum of Dirac delta functions. This closure will be discussed
in detail in §3.

##  3 HyQMOM: Hyperbolic quadrature-based moment closure

In this section, we review the quadrature-based moment-closure (QMOM) approach
and describe in detail the five-moment hyperbolic regularization of QMOM
(HyQMOM), which is the main focus of the current work.

###  3.1 Classical QMOM Approach

The classical quadrature-based moment-closure (QMOM) approach is widely used
in modeling multiphase flows; key developments in this methodology have been
developed over the course of the past several years, e.g., see
article:ChKaMa12 ; article:desjardins2008 ; book:Fox2003 ; article:Fox08 ;
article:fox09 ; article:marchisio2005 ; article:Vikas11 .

The key idea is to assume that the distribution is a sum of Dirac delta
functions whose locations (abscissas) and strengths (weights) are free
parameters:

$\displaystyle{f}(t,x,v)\approx
f^{\star}(t,x,v):=\sum_{j=1}^{N}\omega_{j}\delta\left(v-\mu_{j}\right),$ (21)

where $\delta(v)$ is the Dirac delta and the quadrature weights, $\omega_{j}$,
and abscissas, $\mu_{j}$, are all functions of $t$ and $x$. This approach is
reminiscent of other discrete velocity models such as the Broadwell model
article:broad64a ; article:broad64b ; article:platkowski88 ; however, a key
difference is that the discrete velocities, $\mu_{\ell}$, change with the
solution.

The moment inversion problem requires us to find $(\omega_{j},\mu_{j})$ for
$j=1,\ldots,N$ by matching the first $2N$ moments of $f(t,x,v)$:

${\mathsf{M}}_{\ell}=\sum_{j=1}^{N}\omega_{j}\mu_{j}^{\ell}\quad\text{for}\quad\ell=0,1,\ldots,2N-1.$
(22)

The closure then comes from taking the next moment as follows:

${\mathsf{M}}^{\star}_{2N}=\sum_{j=1}^{N}\omega_{j}\mu_{j}^{2N}.$ (23)

####  3.1.1 Example: $N=2$ case

As a simple example, let us consider the $N=2$ case. The first four moments of
$f^{\star}$ in (21) with $N=2$ are

${\mathsf{M}}_{\ell}=\omega_{1}\mu^{\ell}_{1}+\omega_{2}\mu^{\ell}_{2}\quad\text{for}\quad\ell=0,1,2,3.$
(24)

The moment inversion problem is then this: given
$({\mathsf{M}}_{0},{\mathsf{M}}_{1},{\mathsf{M}}_{2},{\mathsf{M}}_{3})$, find
the parameters $(\mu_{1},\mu_{2},\omega_{1},\omega_{2})$ such that (24) is
satisfied.

This inversion problem is equivalent to finding the quadrature points and
weights for the following weighted Gaussian quadrature rule:

$\int_{-\infty}^{\infty}g(v)\,f(v)\,dv\,\approx\,{\omega}_{1}\,g\left({\mu}_{1}\right)+{\omega}_{2}\,g\left({\mu}_{2}\right),$
(25)

where $f(v)$ is a probability density function with moments
$({\mathsf{M}}_{0},{\mathsf{M}}_{1},{\mathsf{M}}_{2},{\mathsf{M}}_{3})$. If we
attempt to make this quadrature rule exact with $g(v)=1$, $v$, $v^{2}$, and
$v^{3}$, we again arrive at (24).

To find the correct Gaussian quadrature rule, we invoke results from classical
numerical analysis and look for polynomials of degree up to two that are
orthogonal in the weighted $L^{2}(-\infty,\infty)$ inner product:
$\langle\cdot,\cdot\rangle_{f}$. Such polynomials are easily obtained via
Gram-Schmidt, and indeed the relevant one here is the quadratic polynomial:

$\displaystyle\psi_{2}(v)$ $\displaystyle=p^{2}-p\rho\left(v-u\right)^{2}+\rho
h\left(v-u\right),$ (26)

where $\rho$, $u$, $p$, and $h$ are defined by (10)–(11). The abscissas are
the two distinct real roots of ${\psi}_{2}(v)$ and the weights can easily be
obtained by enforcing (24):

$\displaystyle{\mu}_{1},\,{\mu}_{2}=u+\frac{h}{2p}\mp\sqrt{\frac{p}{\rho}+\left(\frac{h}{2p}\right)^{2}},\quad{\omega}_{1},\,\omega_{2}=\frac{\rho}{2}\left[1\pm\frac{\left(\frac{h}{2p}\right)}{\sqrt{\frac{p}{\rho}+\left(\frac{h}{2p}\right)^{2}}}\right].$
(27)

####  3.1.2 Weak hyperbolicity and linear degeneracy of QMOM for all $N\geq 1$

While the above-described process can be used for any $N\geq 1$, it turns out
that the resulting fluid equations are always only weakly hyperbolic.
Furthermore, all of the waves in the Riemann problem solution are linearly
degenerate. We state these facts in the form of theorem A.1 in appendix A; for
completeness, we also provide a full proof of this theorem in appendix A. The
theorem results are well-known in the QMOM literature (e.g., see Chalons et
al. article:ChKaMa12 ), although previously, no theorem had been presented to
rigorously show both the weak hyperbolicity and linear degeneracy of all the
waves for all $N\geq 1$.

###  3.2 Pressure Regularized QMOM

To overcome the weak hyperbolicity present in classical QMOM, Chalons, Fox,
and Massot article:Chalons10  proposed to replace the delta function ansatz
(21) with a multi-Gaussian ansatz of the form:

$\displaystyle{f}(t,x,v)\approx
f^{\star}(t,x,v):=\frac{1}{\sqrt{2\pi\sigma}}\sum_{j=1}^{N}\omega_{j}\exp\left[-\frac{\left(v-\mu_{j}\right)^{2}}{2\sigma}\right],$
(28)

where the free parameters are now the quadrature weights, $\omega_{k}$, the
abscissas, $\mu_{k}$, and the additional parameter $\sigma$. A similar
approach using B-splines was also considered by Cheng and Rossmanith
article:ChRo13 . The additional parameter $\sigma$ allows this closure to
match an additional moment (i.e., a total of $2N+1$ moments can now be
matched), but more importantly, it provides a pressure regularization that
restores strong hyperbolicity. Unfortunately, this closure exhibits a
singularity in the limit of thermodynamic equilibrium, since in that limit,
all the quadrature points collapse to the macroscopic velocity:

$\sigma\rightarrow
T,\qquad\sum_{j=1}^{N}\omega_{j}=\rho,\qquad\text{and}\qquad\mu_{j}\rightarrow
u\quad\forall\,j\in\left[1,N\right].$ (29)

This type of singularity is also evident in other closures, most notably the
maximum entropy closure article:Junk98 .

###  3.3 HyQMOM: Density Regularized QMOM

As an alternative to the above-described pressure regularized QMOM, Fox et al.
article:FoxLauVie2018 ; article:PaDeFo2019  developed a density regularized
version, which they refer to as HyQMOM (hyperbolic quadrature-based method of
moments). This approach was also studied by Johnson thesis:Johnson2017  and
Wiersma thesis:Wiersma2019 . The five-moment HyQMOM system is the subject of
the current work, and we briefly review it in this section.

The idea of HyQMOM is as follows: approximate the distribution as a sum of
delta functions (as in classical QMOM), but place one (or more) of these delta
functions at known locations. This converts the quadrature rule from classical
Gaussian quadrature to something akin to Gauss-Radau quadrature. The version
of this idea relevant to the current work is the case of three delta
functions:

$f\approx
f^{\star}=\omega_{1}\delta(v-\mu_{1})+\omega_{2}\delta(v-u)+\omega_{3}\delta(v-\mu_{3}),$
(30)

where two delta distributions are at unknown locations $\mu_{1},\mu_{3}$ and
the last delta distribution is fixed at the velocity, $u$. Each of the
distributions is weighted by $\omega_{1},\omega_{2},\omega_{3}$. This results
in the following moment-inversion problem:

$\begin{split}\omega_{1}+\omega_{3}&=\widetilde{\rho},\\\
\omega_{1}\mu_{1}+\omega_{3}\mu_{3}&=\widetilde{\rho}u,\\\
\omega_{1}\mu_{1}^{2}+\omega_{3}\mu_{3}^{2}&=\widetilde{\rho}u^{2}+p,\\\
\omega_{1}\mu_{1}^{3}+\omega_{3}\mu_{3}^{3}&=\widetilde{\rho}u^{3}+3pu+h,\end{split}$
(31)

and

$\widetilde{\rho}:=\rho-\omega_{2}=\omega_{1}\left(\frac{\mu_{1}}{u}\right)^{4}+\omega_{3}\left(\frac{\mu_{3}}{u}\right)^{4}-\frac{6\rho
p^{2}u^{2}+4\rho phu+p^{3}+\rho h^{2}+\rho pk}{\rho pu^{4}}.$ (32)

System (31) can be solved in the same way that the $N=2$ classical QMOM system
was solved, namely by constructing a quadratic polynomial:

$\psi_{2}(v)=p^{2}-p\widetilde{\rho}\left(v-u\right)^{2}+h\widetilde{\rho}\left(v-u\right),$
(33)

which is just (26) with $\rho\rightarrow\widetilde{\rho}$. The roots of (33)
provide $\mu_{1}$ and $\mu_{3}$, and the corresponding weights, $\omega_{1}$
and $\omega_{3}$, can easily be computed from (31) once $\mu_{1}$ and
$\mu_{3}$ are known:

$\displaystyle{\mu}_{1},\,{\mu}_{3}=u+\frac{h}{2p}\mp\sqrt{\frac{p}{\widetilde{\rho}}+\left(\frac{h}{2p}\right)^{2}},\quad{\omega}_{1},\,\omega_{3}=\frac{\widetilde{\rho}}{2}\left[1\pm\frac{\left(\frac{h}{2p}\right)}{\sqrt{\frac{p}{\widetilde{\rho}}+\left(\frac{h}{2p}\right)^{2}}}\right].$
(34)

Finally, we can obtain $\omega_{2}$ and fully solve the moment inversion
problem by inserting the above expressions for $\omega_{1}$, $\omega_{3}$,
$\mu_{1}$, and $\mu_{3}$ into (32):

$\widetilde{\rho}=\rho-\omega_{2}=\frac{p^{3}}{p\left(k+\frac{p^{2}}{\rho}+\frac{h^{2}}{p}\right)-h^{2}}.$
(35)

Putting all of these results together yields:

$\begin{split}{\mathsf{M}}_{5}^{\star}=\,\rho
u^{5}+10pu^{3}+10hu^{2}+5ru+\frac{2hr}{p}-\frac{h^{3}}{p^{2}},\quad\text{where}\quad
r:=\frac{p^{2}}{\rho}+\frac{h^{2}}{p}+k.\end{split}$ (36)

From this, we can now assemble the full fluid approximation implied by the
5-moment HyQMOM approximation (30).

######  Definition 3 (5-moment HyQMOM fluid approximation)

The 5-moment HyQMOM approximation can be written either in conservative or
primitive form:

${\underline{q}}_{,t}+{\underline{f}}\left({\underline{q}}\right)_{,x}=\,{\underline{0}}\,\qquad\text{or}\qquad{\underline{\alpha}}_{,t}+{\underline{\underline{B}}}\left({\underline{\alpha}}\right)\,{\underline{\alpha}}_{,x}=\,{\underline{0}}\,,$
(37)

respectively, where

${\underline{q}}=\begin{bmatrix}\rho\\\ \rho u\\\ \rho u^{2}+p\\\ \rho
u^{3}+3pu+h\\\ \rho
u^{4}+6pu^{2}+4hu+r\end{bmatrix},\quad{\underline{f}}\left({\underline{q}}\right)=\begin{bmatrix}\rho
u\\\ \rho u^{2}+p\\\ \rho u^{3}+3pu+h\\\ \rho u^{4}+6pu^{2}+4hu+r\\\
{\mathsf{M}}^{\star}_{5}\end{bmatrix},$ (38)

where $r$ and ${\mathsf{M}}^{\star}_{5}$ are defined by (36), and

${\underline{\alpha}}=\begin{bmatrix}[1.2]\rho\\\ u\\\ p\\\ h\\\
k\end{bmatrix},\quad{\underline{\underline{B}}}\left({\underline{\alpha}}\right)=\begin{bmatrix}[1.2]u&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&\rho&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0\\\
0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&u&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&\frac{1}{\rho}&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0\\\
0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&3p&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&u&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&1&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0\\\
-\frac{p^{2}}{\rho^{2}}&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&4h&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&-\frac{h^{2}}{p^{2}}-\frac{p}{\rho}&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&u+\frac{2h}{p}&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&1\\\
0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&5k&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&-\frac{2kh}{p^{2}}&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&\frac{2k}{p}&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>(39)

Furthermore, we note that that conservative flux Jacobian has the following
form:

${\underline{\underline{A}}}\left({\underline{q}}\right):={\underline{f}}\left({\underline{q}}\right)_{,{\underline{q}}}=\begin{bmatrix}[1.2]0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&1&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0\\\
0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&1&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0\\\
0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&1&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0\\\
0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&0&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&1\\\
\frac{\partial{\mathsf{M}}^{\star}_{5}}{\partial{\mathsf{M}}_{0}}&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&\frac{\partial{\mathsf{M}}^{\star}_{5}}{\partial{\mathsf{M}}_{1}}&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&\frac{\partial{\mathsf{M}}^{\star}_{5}}{\partial{\mathsf{M}}_{2}}&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&\frac{\partial{\mathsf{M}}^{\star}_{5}}{\partial{\mathsf{M}}_{3}}&\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}{-}&\frac{\partial{\mathsf{M}}^{\star}_{5}}{\partial{\mathsf{M}}_{4}}\end{bmatrix},$
(40)

where the details of the last row have been omitted for brevity. The matrices
${\underline{\underline{A}}}$ and ${\underline{\underline{B}}}$ as defined in
(40) and (39), respectively, are similar matrices, meaning that they share the
same eigenvalues.

######  Proposition 1 (Hyperbolicity of HyQMOM)

The HyQMOM fluid model as expressed either in conservative or primitive form
(37)–(39) is strictly hyperbolic for all $\rho$, $p$, $k$ $>0$. Note that
$\rho$, $p$, $k$ $>0$ defines a convex set in the solution space:
${\underline{q}}\in\mathbb{R}^{5}$. The resulting wave structure includes one
linearly degenerate wave and four nonlinear waves.

Note that the essential contents of this theorem were already known in Fox,
Laurent, and Vié article:FoxLauVie2018 , but never presented as a formal
proposition and proved. Therefore, we include the proof here.

###### Proof

The eigenvalues of the Jacobian matrix,
${\underline{\underline{B}}}\left({\underline{\alpha}}\right)$, in (39) can be
computed explicitly:

$\begin{gathered}\lambda_{1}=u+\frac{h}{2p}-\sqrt{a+b},\quad\lambda_{2}=u+\frac{h}{2p}-\sqrt{a-b},\quad\lambda_{3}=u,\\\
\lambda_{4}=u+\frac{h}{2p}+\sqrt{a-b},\quad\lambda_{5}=u+\frac{h}{2p}+\sqrt{a+b},\end{gathered}$
(41)

where

$a=\frac{p}{\rho}+\frac{k}{p}+\left(\frac{h}{2p}\right)^{2}\qquad\text{and}\qquad
b=\sqrt{\frac{k^{2}}{p^{2}}+\frac{k}{\rho}}.$ (42)

One can show via simple calculations that for all $\rho$, $p$, $k$ $>0$:

$\displaystyle
a>b>0,\quad\sqrt{a+b}>\frac{|h|}{2p},\quad\text{and}\quad\sqrt{a-b}>\frac{|h|}{2p}.$
(43)

Therefore, all five eigenvalues shown in (41)–(42) are real and distinct for
all $\rho$, $p$, $k$ $>0$, which is sufficient to show that system (37)–(39)
is strictly hyperbolic.

The eigenvectors of the flux Jacobian (40) can be written as follows:

${\underline{{\mathcal{R}}_{\ell}}}=\left[1,\lambda_{\ell},\lambda_{\ell}^{2},\lambda_{\ell}^{3},\lambda_{\ell}^{4}\right]^{T}\qquad\text{for}\qquad\ell=1,2,3,4,5.$
(44)

To determine whether the corresponding waves are linearly degenerate or not,
we need to compute the quantities:

$\frac{\partial\lambda_{\ell}}{\partial{\underline{q}}}\cdot{\underline{{\mathcal{R}}_{\ell}}}=\frac{\partial\lambda_{\ell}}{\partial{\underline{\alpha}}}\cdot\left(\frac{\partial{\underline{q}}}{\partial{\underline{\alpha}}}\right)^{-1}\cdot{\underline{{\mathcal{R}}_{\ell}}},\quad\forall\ell=1,2,3,4,5,$
(45)

where

$\frac{\partial\lambda_{\ell}}{\partial{\underline{q}}}:=\left[\frac{\partial\lambda_{\ell}}{\partial
q_{1}},\,\frac{\partial\lambda_{\ell}}{\partial
q_{2}},\,\frac{\partial\lambda_{\ell}}{\partial
q_{3}},\,\frac{\partial\lambda_{\ell}}{\partial
q_{4}},\,\frac{\partial\lambda_{\ell}}{\partial q_{5}}\right],$ (46)

and

$\frac{\partial{\underline{q}}}{\partial{\underline{\alpha}}}=\begin{bmatrix}[1.5]1&&0&&0&&0&&0\\\
u&&\rho&&0&&0&&0\\\ u^{2}&&2\rho u&&1&&0&&0\\\ u^{3}&&3\left(p+\rho
u^{2}\right)&&3u&&1&&0\\\
\frac{\rho^{2}u^{4}-p^{2}}{\rho^{2}}&&4\left(h+3pu+\rho
u^{3}\right)&&\frac{6\rho p^{2}u^{2}+2p^{3}-\rho h^{2}}{\rho
p^{2}}&&\frac{2h+4pu}{p}&&1\end{bmatrix}.$ (47)

We note that one of the waves is linearly degenerate while the remaining are
nonlinear:

$\frac{\partial\lambda_{3}}{\partial{\underline{q}}}\cdot{\underline{{\mathcal{R}}_{3}}}\equiv
0\qquad\text{and}\qquad\frac{\partial\lambda_{\ell}}{\partial{\underline{q}}}\cdot{\underline{{\mathcal{R}}_{\ell}}}\not\equiv
0\quad\text{for}\quad\ell=1,2,4,5.$ (48)

$\blacksquare$

##  4 Locally-implicit Lax-Wendroff discontinuous Galerkin

We consider generic one-dimensional conservation laws of the form:

${\underline{q}}_{,t}+{\underline{f}}\left({\underline{q}}\right)_{,x}={\underline{0}},$
(49)

where $t$ is time, $x$ is space,
$q(t,x):\in\mathbb{R}^{+}\times\mathbb{R}\mapsto\mathbb{R}^{M_{\text{eqn}}}$
is the vector of conserved variables, ${M_{\text{eqn}}}$ is the number of
equations, and
${\underline{f}}\left({\underline{q}}\right):\mathbb{R}^{M_{\text{eqn}}}\mapsto\mathbb{R}^{M_{\text{eqn}}}$
is the flux function. We assume that this system is hyperbolic, meaning that
the flux Jacobian:

${\underline{f}}\left({\underline{q}}\right)_{,{\underline{q}}}:\mathbb{R}^{M_{\text{eqn}}}\mapsto\mathbb{R}^{{M_{\text{eqn}}}\times{M_{\text{eqn}}}},$
(50)

has real eigenvalues and a complete set of eigenvectors over some convex
region ${\mathcal{D}}\subset\mathbb{R}^{{M_{\text{eqn}}}}$ in solution space
inside of which we are interested in solving the equation.

The Lax-Wendroff method article:LxW60  is a time discretization for hyperbolic
conservation laws based on the the Cauchy-Kovalevskaya article:Kovaleskaya1875
procedure to convert temporal derivatives into spatial derivatives. We begin
with a Taylor series in time:

${\underline{q}}(t+\Delta t,x)={\underline{q}}(t,x)+\Delta
t{\underline{q}}_{,t}(t,x)+\frac{1}{2}\Delta
t^{2}{\underline{q}}_{,t,t}(t,x)+\ldots,$ (51)

and then replace all time derivatives by spatial derivatives:

${\underline{q}}_{,t}=-{\underline{f}}\left({\underline{q}}\right)_{,x},\quad{\underline{q}}_{,t,t}=-{\underline{f}}\left({\underline{q}}\right)_{,t,x}=-\left[{\underline{\underline{f^{\prime}}}}\left({\underline{q}}\right)\,{\underline{q}}_{,t}\right]_{,x}=\left[{\underline{\underline{f^{\prime}}}}\left({\underline{q}}\right){\underline{f}}\left({\underline{q}}\right)_{,x}\right]_{,x},\quad\ldots,$
(52)

which results in the following:

${\underline{q}}(t+\Delta t,x)={\underline{q}}-\Delta
t{\underline{f}}\left({\underline{q}}\right)_{,x}+\frac{1}{2}\Delta
t^{2}\left[{\underline{\underline{f^{\prime}}}}\left({\underline{q}}\right){\underline{f}}\left({\underline{q}}\right)_{,x}\right]_{,x}+\ldots,$
(53)

where on the right-hand side, we have suppressed the evaluation at $(t,x)$.
The final step is to truncate the Taylor series at some finite number of
terms, and then replace all spatial derivatives by some discrete spatial
derivative operators. The above Lax-Wendroff formalism article:LxW60  has been
used in conjunction with with a variety of spatial discretizations, including
finite volume article:Le97 , weighted essentially non-oscillatory (WENO)
article:TitaToro02 , and discontinuous Galerkin article:Qiu05  operators.

In this work, we are concerned with the discontinuous Galerkin version of Lax-
Wendroff article:Qiu05 ; and in particular, we make use of the reformulation
of Gassner et al. article:GasDumHinMun2011  of the Lax-Wendroff discontinuous
Galerkin (LxW-DG) scheme in terms of a locally-implicit prediction step,
followed by an explicit correction step. The key advantage of this formulation
is that we do not need to explicitly compute the partial derivatives as shown
in (52); and instead, the locally-implicit solver automatically produces
discrete versions of these derivatives. The next challenge is to efficiently
solve the nonlinear algebraic equations arising from the locally-implicit
prediction step; we solve these equations by following Gassner et al.
article:GasDumHinMun2011  and making use of a Picard fixed point iteration. We
will follow the notational conventions of Guthrey and Rossmanith
article:GuRo2017  developed for locally-implicit and regionally-implicit LxW-
DG schemes. Note that the predictor-correct method is equivalent to the Lax-
Wendroff DG method for linear constant-coefficient hyperbolic systems. For
nonlinear systems, the predictor-correct method differs slightly from Lax-
Wendroff DG in that the Picard iteration inside the prediction step
approximates the direct computation of the nonlinear Taylor series expansion.

###  4.1 Discontinuous Galerkin finite elements

To discretize equation (49) in space, we use the discontinuous Galerkin (DG)
finite element method, which was first introduced by Reed and Hill
article:ReedHill73 . It was fully developed for time-dependent hyperbolic
conservation laws in a series of papers by Bernardo Cockburn, Chi-Wang Shu,
and collaborators (see article:CoShu98  and references therein for details).

We define the broken finite element space

$\mathcal{W}^{\Delta x}:=\left\\{w^{\Delta
x}\in[L^{\infty}(\Omega)]^{M_{\text{eqn}}}:w^{\Delta
x}|_{\mathcal{T}_{i}}\in[\mathbb{P}(M_{\text{deg}})]^{M_{\text{eqn}}}\quad\forall\mathcal{T}_{i}\right\\},$
(54)

where $\Delta x=(x_{\text{high}}-x_{\text{low}})/M_{\text{elem}}$ is a uniform
grid spacing with $M_{\text{elem}}$ being the number of elements.
Additionally, $M_{\text{eqn}}$ is the number of conserved variables,
$\mathbb{P}(M_{\text{deg}})$ is the set of all polynomials of degree at most
$M_{\text{deg}}$, and the computational mesh is described by non-overlapping
elements of width $\Delta x$ centered at the points $x_{i}$:

$\mathcal{T}_{i}=\left[x_{i}-\frac{\Delta x}{2},\,x_{i}+\frac{\Delta
x}{2}\right]\hskip 20.00003pt\text{for }i=1,\dots,M_{\text{elem}}.$ (55)

On each element, we define the local spatial variable, $\xi$:

$x=x_{i}+\left(\frac{\Delta x}{2}\right)\xi\quad\text{for}\quad\xi\in[-1,1].$
(56)

On each element, we approximate the solution by a finite expansion in terms of
the following orthonormal Legendre polynomial basis functions:

${\underline{\Phi}}=\left(1,\,\sqrt{3}\xi,\frac{\sqrt{5}}{2}\left(3\xi^{2}-1\right),\,\frac{\sqrt{7}}{2}\left(5\xi^{3}-3\xi\right),\,\frac{\sqrt{9}}{8}\left(35\xi^{4}-30\xi^{2}+3\right),\,\dots\right),$
(57)

with the orthonormality property:

$\frac{1}{2}\int_{-1}^{1}{\underline{\Phi}}\,{\underline{\Phi}}^{T}\,d\xi={\underline{\underline{\,\,{\mathbb{I}}\,\,}}},$
(58)

where ${\underline{\underline{\,\,{\mathbb{I}}\,\,}}}$ is the identity matrix.
On each element at time $t=t^{n}$, we approximate the solution as follows:

${\underline{q}}^{h}\left(t^{n},x_{i}+\frac{\Delta
x}{2}\xi\right):={\underline{\Phi}}(\xi)^{T}{\underline{\underline{Q_{i}^{n}}}}\quad\text{for}\quad\xi\in[-1,1],$
(59)

where

$\displaystyle{\underline{\Phi}}(\xi):[-1,1]\mapsto\mathbb{R}^{{M_{\text{C}}}}\quad\text{and}\quad{\underline{\underline{Q_{i}^{n}}}}\in\mathbb{R}^{{M_{\text{C}}}\times{M_{\text{eqn}}}}.$
(60)

The number of basis functions in 1D is ${M_{\text{C}}}={M_{\text{deg}}}+1$,
and the order of accuracy is ${M_{\text{O}}}={M_{\text{deg}}}+1$.

###  4.2 Prediction step

The prediction step is completely local to each element, and therefore, the
prediction step is inconsistent with the underlying conservation law. This
inconsistency allows us to freely choose updated variables; we use the
primitive variables for this step: ${\underline{\alpha}}=(\rho,u,p,h,k)$.

The prediction step is local on each space-element
$\left[t^{n},t^{n+1}\right]\times{\mathcal{T}}_{i}$, where
$t^{n+1}=t^{n}+\Delta t$. Let $t=t^{n}+\frac{\Delta t}{2}(1+\tau)$, for
$\tau\in[-1,1]$ and

${\underline{\alpha}}_{,\tau}={\underline{\Theta}}({\underline{\alpha}}):=-\frac{\Delta
t}{\Delta
x}\,{\underline{\underline{B}}}\left({\underline{\alpha}}\right){\underline{\alpha}}_{,\xi},$
(61)

where ${\underline{\underline{B}}}({\underline{\alpha}})$ is the (primitive
variable) flux Jacobian matrix defined by (39). We introduce a space-time
Legendre basis on each element:

$\Psi_{\ell}(\tau,\xi)=\Phi_{\ell_{1}}(\tau)\Phi_{\ell_{2}}(\xi),\quad\text{for}\quad\ell_{1}=1,\ldots,{M_{\text{O}}},\quad\ell_{2}=1,\ldots,{M_{\text{O}}}+1-\ell_{1},$
(62)

where

$\ell={M_{\text{O}}}\left(\ell_{1}-1\right)-\frac{\left(\ell_{1}-1\right)\left(\ell_{1}-2\right)}{2}+\ell_{2}\quad\text{such
that}\quad\ell=1,\ldots,{M_{\text{P}}},$ (63)

and ${M_{\text{P}}}={M_{\text{O}}}({M_{\text{O}}}+1)/2$ is the number of
space-time basis functions. These space-time basis functions are orthonormal
on $[-1,1]^{2}$:

$\frac{1}{4}\iint_{-1}^{1}{\underline{\Psi}}\,{\underline{\Psi}}^{T}\,d\tau\,d\xi={\underline{\underline{\,\,\mathbb{I}\,\,}}}.$
(64)

We write the predicted solution as follows:

${\underline{\alpha}}^{\text{ST}}\left(t^{n}+\frac{\Delta
t}{2}(1+\tau),\,x_{i}+\frac{\Delta
x}{2}\,\xi\right):={\underline{\Psi}}\left(\tau,\xi\right)^{T}{\underline{\underline{W^{n+\frac{1}{2}}_{i}}}},\qquad{\underline{\underline{W_{i}^{n+\frac{1}{2}}}}}\in\mathbb{R}^{{M_{\text{P}}}\times{M_{\text{eqn}}}},$
(65)

for $(\tau,\xi)\in[-1,1]^{2}$, where ${\underline{\underline{W}}}$ represents
the matrix of unknown coefficients.

We proceed by multiplying (61) by the test function ${\underline{\Psi}}$, then
integrating over $(\tau,\xi)\in[-1,1]^{2}$:

$\frac{1}{4}\iint_{-1}^{1}{\alpha}_{m,\tau}\,{\underline{\Psi}}\,d\tau\,d\xi=\frac{1}{4}\,\iint_{-1}^{1}{\Theta}_{m}({\underline{\alpha}})\,{\underline{\Psi}}\,d\tau\,d\xi,$
(66)

and applying integration-by-parts in $\tau$ only:

$\begin{gathered}\frac{1}{4}\int_{-1}^{1}{\alpha}^{\star}_{m}\left(\tau=1,\xi\right)\,{\underline{\Psi}}\left(\tau=1,\xi\right)\,d\xi-\frac{1}{4}\int_{-1}^{1}{\alpha}^{\star}_{m}\left(\tau=-1,\xi\right)\,{\underline{\Psi}}\left(\tau=-1,\xi\right)\,d\xi\\\
-\frac{1}{4}\iint_{-1}^{1}{\alpha}_{m}\,{\underline{\Psi}}_{,\tau}\,d\tau\,d\xi=\frac{1}{4}\,\iint_{-1}^{1}{\Theta}_{m}({\underline{\alpha}})\,{\underline{\Psi}}\,d\tau\,d\xi,\end{gathered}$
(67)

where the choice of values for ${\alpha}^{\star}_{m}$ at $\tau=1$ and
$\tau=-1$ still remains to be made. Before making this choice, however, let us
reverse integrate-by-parts, such that the newly introduced boundary terms are
always the internal values of the current space-time element:

$\begin{gathered}\frac{1}{4}\int_{-1}^{1}\Bigl{[}{\alpha}^{\star}_{m}\left(\tau=1,\xi\right)-{\alpha}_{m}\left(\tau=1,\xi\right)\Bigr{]}\,{\underline{\Psi}}\left(\tau=1,\xi\right)\,d\xi\\\
-\frac{1}{4}\int_{-1}^{1}\Bigl{[}{\alpha}^{\star}_{m}\left(\tau=-1,\xi\right)-{\alpha}_{m}\left(\tau=-1,\xi\right)\Bigr{]}\,{\underline{\Psi}}\left(\tau=-1,\xi\right)\,d\xi\\\
+\frac{1}{4}\iint_{-1}^{1}{\underline{\Psi}}\,{\alpha}_{m,\tau}\,d\tau\,d\xi=\frac{1}{4}\,\iint_{-1}^{1}{\Theta}_{m}({\underline{\alpha}})\,{\underline{\Psi}}\,d\tau\,d\xi.\end{gathered}$

Now we plug-in the following values for ${\alpha}^{\star}_{m}$ at $\tau=1$ and
$\tau=-1$:

${\alpha}^{\star}_{m}\left(\tau=1,\xi\right)={\alpha}_{m}\left(\tau=1,\xi\right)\quad\text{and}\quad{\alpha}^{\star}_{m}\left(\tau=-1,\xi\right)={\underline{\Phi}}\left(\xi\right)^{T}\,{\underline{A^{n}_{i\,(:,m)}}},$
(68)

as well as the following values for ${\alpha}_{m}$ on the interior of the
space-time element:

${\alpha}_{m}\left(\tau,\xi\right)={\underline{\Psi}}\left(\tau,\xi\right)^{T}{\underline{W^{n+\frac{1}{2}}_{i(:,m)}}}.$
(69)

The result is

$\begin{split}{\underline{\underline{L}}}\,\,{\underline{W^{n+\frac{1}{2}}_{i\,(:,m)}}}&=\frac{1}{4}\iint_{-1}^{1}{\Theta}_{m}\left({\underline{\Psi}}^{T}{\underline{\underline{W^{n+\frac{1}{2}}_{i}}}}\right){\underline{\Psi}}\,d\tau\,d\xi\\\
&+\left[\frac{1}{4}\int_{-1}^{1}\,{\underline{{\Psi}}}\left(-1,\xi\right)\,{\underline{\Phi}}\left(\xi\right)^{T}\,d\xi\right]\,{\underline{A^{n}_{i\,(:,m)}}},\end{split}$
(70)

for each equation $m=1,\ldots,{M_{\text{eqn}}}$, where

$\displaystyle{\underline{\underline{L}}}$
$\displaystyle=\frac{1}{4}\iint_{-1}^{1}{\underline{\Psi}}\,{\underline{\Psi}}^{T}_{,\tau}\,d\tau\,d\xi+\frac{1}{4}\int_{-1}^{1}{\underline{\Psi}}_{|_{\tau=-1}}{\underline{\Psi}}_{|_{\tau=-1}}^{T}\,d\xi,$
(71) $\displaystyle{\underline{\underline{A_{i}^{n}}}}$
$\displaystyle=\frac{1}{2}\sum_{a=1}^{{M_{\text{O}}}}\omega_{a}\,{\underline{\Phi}}\left(\mu_{a}\right)\left[{\underline{\alpha}}\left({\underline{\Phi}}\left(\mu_{a}\right)^{T}{\underline{\underline{Q_{i}^{n}}}}\right)\right]^{T},$
(72)

and ${\underline{\alpha}}({\underline{q}})$ gives the relationship between
conservative and primitive variables. Equation (70) is a nonlinear algebraic
equation that must be solved on each space-time element for the matrix of
unknown coefficients: ${\underline{\underline{W}}}$.

Following Gassner et al. article:GasDumHinMun2011 , instead of using Newton’s
method to solve the resulting non-linear equation, which involves inverting a
Jacobian matrix at every step, we used the much simpler Picard iteration.
There are two key advantages of the Picard iteration over Newton’s method.
First, since ${\underline{\underline{L}}}$ is independent of the solution and
the same on each space-time element, we only invert this relatively small
matrix once. Second, the Picard iteration converges to sufficiently high order
accuracy after ${M_{\text{O}}}-1$ iterations, so the need to compute residuals
is eliminated (see justification in article:GasDumHinMun2011 ). We can write
the Picard iteration as

$\begin{split}{\underline{W^{n+\frac{1}{2}(j)}_{i\,(:,m)}}}=&\frac{1}{4}\sum_{a=1}^{{M_{\text{O}}}}\sum_{b=1}^{{M_{\text{O}}}}\,\omega_{a}\,\omega_{b}\,{\underline{\hat{\Psi}}}\left(\mu_{b},\,\mu_{a}\right)\,{\Theta}_{m}\left({\underline{\Psi}}\left(\mu_{b},\mu_{a}\right)^{T}{\underline{\underline{W^{n+\frac{1}{2}(j-1)}_{i}}}}\right)\\\
+&\frac{1}{4}\sum_{b=1}^{{M_{\text{O}}}}\omega_{b}\,{\underline{\hat{\Psi}}}\left(-1,\xi_{b}\right)\,{\underline{\Phi}}\left(\xi_{b}\right)^{T}{\underline{A^{n}_{i\,(:,m)}}},\end{split}$
(73)

where $j=1,\ldots,{M_{\text{O}}}-1$ is the iteration count,
$m=1,\ldots,{M_{\text{eqn}}}$ is the equation index,
${\underline{\hat{\Psi}}}={\underline{\underline{L}}}^{-1}{\underline{\Psi}}$,
and $\omega_{a}$ and $\mu_{a}$ for $a=1,\ldots,{M_{\text{O}}}$ are the weights
and abscissas of the ${M_{\text{O}}}$-point Gauss-Legendre quadrature rule.
This gives a solution for the prediction step, which we know is inconsistent
with the conservation law. To make the final solution consistent (and high-
order) with the original conservation law, we next need to add a correction
step.

###  4.3 Correction step

The correction step is designed to work like a single forward Euler-like step
that uses the predicted solution. To perform this step, we begin with the
hyperbolic conservation law (49) and multiply by the spatial basis functions
defined in (57). Next, we integrate over the space-time element
$\left(\tau,\xi\right)\in\left[-1,1\right]^{2}$:

$\frac{1}{2}\iint_{-1}^{1}\left[{\underline{\Phi}}(\xi)\,{\underline{q}}_{,\tau}+\frac{\Delta
t}{\Delta
x}\,{\underline{\Phi}}(\xi)\,{\underline{f}}\left({\underline{q}}\right)_{,\xi}\right]\,d\xi\,d\tau={\underline{0}},$
(74)

which can be written as

$\begin{split}\frac{1}{2}\int_{-1}^{1}{\underline{\Phi}}(\xi)\,{\underline{q}}\left(t^{n+1},x_{i}+\frac{\Delta
x}{2}\xi\right)\,\,d\xi&=\frac{1}{2}\int_{-1}^{1}{\underline{\Phi}}(\xi)\,{\underline{q}}\left(t^{n},x_{i}+\frac{\Delta
x}{2}\xi\right)\,\,d\xi\\\ &-\frac{\Delta t}{2\Delta
x}\iint_{-1}^{1}{\underline{\Phi}}(\xi)\,{\underline{f}}\left({\underline{q}}\right)_{,\xi}\,d\xi\,d\tau.\end{split}$
(75)

We approximate ${\underline{q}}\left(t^{n+1},\cdot\right)$ and
${\underline{q}}\left(t^{n},\cdot\right)$ in (75) via appropriate versions of
ansatz (59). For the remaining term, we first apply integration-by-parts in
space, then replace the true solution $q$ by its space-time predicted
solution: (65), and replace exact integration by numerical quadrature. This
results in the following expression:

$\begin{split}{\underline{\underline{Q_{i}^{n+1}}}}&={\underline{\underline{Q^{n}_{i}}}}+\frac{\Delta
t}{2\Delta
x}\sum_{a=1}^{{M_{\text{O}}}}\sum_{b=1}^{{M_{\text{O}}}}\,\omega_{a}\,\omega_{b}\,{\underline{\Phi}}_{,\xi}\left(\mu_{a}\right)\left[{\underline{f}}\left({\underline{\Psi}}\left(\mu_{b},\mu_{a}\right)^{T}{\underline{\underline{W^{n+\frac{1}{2}}_{i}}}}\right)\right]^{T}\\\
&-{\frac{\Delta t}{\Delta
x}}\left(\,{\underline{\Phi}}(1)\left[{\underline{{\mathcal{F}}^{n+\frac{1}{2}}_{i+\frac{1}{2}}}}\right]^{T}-{\underline{\Phi}}(-1)\left[{\underline{{\mathcal{F}}^{n+\frac{1}{2}}_{i-\frac{1}{2}}}}\right]^{T}\,\right).\end{split}$
(76)

In the expressions after the approximation symbol, $\approx$, we replaced all
exact integration by a Gauss-Legendre quadrature, where $\omega_{a}$ and
$\mu_{a}$ for $a=1,\ldots,{M_{\text{O}}}$ are the weights and abscissas of the
${M_{\text{O}}}$-point Gauss-Legendre quadrature rule.

The time-integrated numerical fluxes are defined using the predicted solution
and the Rusanov article:Ru61  time-averaged flux:

$\displaystyle{\underline{{\mathcal{F}}^{n+\frac{1}{2}}_{i-\frac{1}{2}}}}:=\frac{1}{2}\sum_{a=1}^{{M_{\text{O}}}}\omega_{a}\,{\underline{\mathcal{F}}}\left(\mu_{a}\right),$
(77)

where the numerical flux at each temporal quadrature point is given by

${\underline{\mathcal{F}}}\left(\tau\right):=\frac{1}{2}\left\\{{\underline{f}}\left({\underline{W_{\text{R}}}}(\tau)\right)+{\underline{f}}\left({\underline{W_{\text{L}}}}(\tau)\right)\right\\}-\frac{1}{2}{\lambda_{\max}}(\tau)\left\\{{\underline{q}}\left({\underline{W_{\text{R}}}}(\tau)\right)-{\underline{q}}\left({\underline{W_{\text{L}}}}(\tau)\right)\right\\},$
(78)

where

$\displaystyle{\underline{W_{\text{L}}}}(\tau):={\underline{\Psi}}\left(\tau,1\right)^{T}{\underline{\underline{W^{n+\frac{1}{2}}_{i-1}}}},\qquad{\underline{W_{\text{R}}}}(\tau):={\underline{\Psi}}\left(\tau,-1\right)^{T}{\underline{\underline{W^{n+\frac{1}{2}}_{i}}}},$
(79)

and ${\lambda_{\max}}(\tau)$ is a local bound on the spectral radius of the
flux Jacobian,
${\underline{f}}({\underline{q}})_{,{\underline{q}}}={\underline{\underline{A}}}\left({\underline{q}}\right)$,
in the neighborhood of interface $x=x_{i-\frac{1}{2}}$ and at time
$t=t^{n}+(\tau+1)\Delta t/2$.

order | ${M_{\text{O}}}=1$ | ${M_{\text{O}}}=2$ | ${M_{\text{O}}}=3$ | ${M_{\text{O}}}=4$
---|---|---|---|---
CFL # used in practice | $0.90$ | $0.30$ | $0.14$ | $0.09$
Table 1: The CFL numbers used in all simulations as a function of the method
order.

These steps are all it takes to regain the coupling neglected in the
prediction step. We now have a solution that is not only consistent with the
conservation law but is also high order. The Courant-Friedrichs-Lewy (CFL)
number we use for each order is given in table 1. However, we still have some
work to do to ensure that the solution is physical. We must be careful to
maintain the positivity of the primitive variables $\rho$, $p$, $k$, as was
necessary for the system’s hyperbolicity and moment-realizability. We address
the limiters utilized to accomplish this in the next section.

##  5 HyQMOM limiters

The high-order numerical method as described in §4 does not guarantee that
density, pressure, and modified kurtosis remain positive throughout a time-
step:

$\rho^{n}>0,\quad p^{n}>0,\quad
k^{n}>0\quad\centernot\Longrightarrow\quad\rho^{n+1}>0,\quad p^{n+1}>0,\quad
k^{n+1}>0.$ (80)

Recall that positivity of these quantities is needed to guarantee the moment-
realizability of the moment-closure and strict hyperbolicity of the resulting
evolution equations. If we want positivity over a time step, we will need to
introduce positivity-preserving limiters. Additionally, if we want to control
unphysical oscillations near large gradients, shocks, and rarefactions, we
will also need non-oscillatory limiters. In this section, we derive all of
these limiters. In particular, we first need to establish that a simple first-
order scheme is positivity-preserving under some appropriate time-step
restriction; this is done in §5.1. Using this result we derive a suite of
limiters that ensure positivity: §5.2 (positivity at select points in the
prediction step), §5.3 (positivity of the average density, pressure, and
modified kurtosis in each element in the correction step), and §5.4
(positivity at select points in the correction step). We then develop an
unphysical oscillation limiter in §5.5.

###  5.1 Positivity of the Rusanov scheme

Before considering positivity limiters for the high-order method, we must
first establish that simple first-order schemes, in this case, we consider the
Rusanov (aka local Lax-Friedrichs) scheme article:Ru61 , are positivity-
preserving under some appropriate time-step restriction. This is established
in the theorem below, which is an extension of the result of Zhang and Shu
article:ZhangShu2011  for the compressible Euler equations.

###### Theorem 5.1

Let ${\underline{Q_{i}^{n}}}$ be some approximation of the element averages of
the conserved variables on element ${\mathcal{T}}_{i}$ at time $t^{n}$, and
let ${\underline{Q_{i}^{n+1}}}$ be the element averages produced by the
Rusanov scheme article:Ru61  at time $t^{n+1}=t^{n}+\Delta t$. Then

$\quad\rho_{i}^{n}>0,\,\,\,\,p_{i}^{n}>0,\,\,\,\,k_{i}^{n}>0\,\,\,\,\,\forall
i\quad\Longrightarrow\quad\rho^{n+1}_{i}>0,\,\,\,\,p^{n+1}_{i}>0,\,\,\,\,k^{n+1}_{i}>0\,\,\,\,\,\forall
i,$ (81)

under the CFL condition:

$\frac{\Delta t}{\Delta x}\max_{i}\left({\lambda_{i+\frac{1}{2}}}\right)<1,$
(82)

where

${\lambda_{i\pm\frac{1}{2}}}=\max\left\\{{\lambda_{\max}}\left({\underline{Q_{i}^{n}}}\right),\,{\lambda_{\max}}\left({\underline{Q^{n}_{i\pm
1}}}\right),\,{\lambda_{\max}}\left(\frac{1}{2}\left({\underline{Q_{i}^{n}}}+{\underline{Q^{n}_{i\pm
1}}}\right)\right)\right\\},$ (83)

and ${\lambda_{\max}}\left({\underline{Q_{i}^{n}}}\right)$ is a bound on the
spectral radius of the flux Jacobian (40) at state ${\underline{Q_{i}^{n}}}$.

###### Proof

Recall that the Rusanov scheme can be written as

${\underline{Q_{i}^{n+1}}}={\underline{Q_{i}^{n}}}-\frac{\Delta t}{\Delta
x}\left({\underline{{\mathcal{F}}_{i+\frac{1}{2}}}}-{\underline{{\mathcal{F}}_{i-\frac{1}{2}}}}\right),$
(84)

where the numerical fluxes are given by

${\underline{{\mathcal{F}}_{i\pm\frac{1}{2}}}}=\frac{1}{2}\left[{\underline{f}}({\underline{Q_{i}^{n}}})+{\underline{f}}({\underline{Q_{i\pm
1}^{n}}})\right]-\frac{1}{2}{\lambda_{i\pm\frac{1}{2}}}\left({\underline{Q_{i}^{n}}}-{\underline{Q^{n}_{i\pm
1}}}\right),$ (85)

where the flux function, ${\underline{f}}\left({\underline{q}}\right)$, is
defined by (38) and the local wave speed, ${\lambda_{i\pm\frac{1}{2}}}$, is
defined by (83). We can rewrite the above expression into the following
numerical update:

$\begin{split}{\underline{Q_{i}^{n+1}}}&=\left[1-\frac{\Delta t}{2\Delta
x}\left({\lambda_{i+\frac{1}{2}}}+{\lambda_{i-\frac{1}{2}}}\right)\right]{\underline{Q_{i}^{n}}}+\left[\frac{\Delta
t{\lambda_{i+\frac{1}{2}}}}{2\Delta
x}\right]{\underline{M^{+}_{i}}}+\left[\frac{\Delta
t{\lambda_{i-\frac{1}{2}}}}{2\Delta
x}\right]{\underline{M^{-}_{i}}},\end{split}$ (86)

where

$\displaystyle{\underline{M^{+}_{i}}}={\underline{Q_{i+1}^{n}}}-\left({\lambda_{i+\frac{1}{2}}}\right)^{-1}{{\underline{f}}\left({\underline{Q_{i+1}^{n}}}\right)}\quad\text{and}\quad{\underline{M^{-}_{i}}}={\underline{Q_{i-1}^{n}}}+\left({\lambda_{i-\frac{1}{2}}}\right)^{-1}{{\underline{f}}\left({\underline{Q_{i-1}^{n}}}\right)}.$
(87)

Under the CFL condition (82), we note that

$1-\frac{\Delta t}{2\Delta
x}\left({\lambda_{i+\frac{1}{2}}}+{\lambda_{i-\frac{1}{2}}}\right)\geq
0,\quad\frac{\Delta t{\lambda_{i+\frac{1}{2}}}}{2\Delta x}\geq
0,\quad\text{and}\quad\frac{\Delta t{\lambda_{i-\frac{1}{2}}}}{2\Delta x}\geq
0.$ (88)

Additionally, note that the coefficients in (86) sum to unity:

$\left[1-\frac{\Delta t}{2\Delta
x}\left({\lambda_{i+\frac{1}{2}}}+{\lambda_{i-\frac{1}{2}}}\right)\right]+\left[\frac{\Delta
t{\lambda_{i+\frac{1}{2}}}}{2\Delta x}\right]+\left[\frac{\Delta
t{\lambda_{i-\frac{1}{2}}}}{2\Delta x}\right]=1.$ (89)

Now let ${\mathcal{C}}_{\alpha}$ be any convex function of the conserved
variables:
${\underline{q}}=\left({\mathsf{M}}_{0},{\mathsf{M}}_{1},{\mathsf{M}}_{2},{\mathsf{M}}_{3},{\mathsf{M}}_{4}\right)$.
Then, applying the convex function to both sides of (86) we see that

$\begin{split}{\mathcal{C}}_{\alpha}\left({\underline{Q_{i}^{n+1}}}\right)\leq&\left[1-\frac{\Delta
t}{2\Delta
x}\left({\lambda_{i+\frac{1}{2}}}+{\lambda_{i-\frac{1}{2}}}\right)\right]{\mathcal{C}}_{\alpha}\left({\underline{Q_{i}^{n}}}\right)\\\
+&\left[\frac{\Delta t{\lambda_{i+\frac{1}{2}}}}{2\Delta
x}\right]{\mathcal{C}}_{\alpha}\left({\underline{M^{+}_{i}}}\right)+\left[\frac{\Delta
t{\lambda_{i-\frac{1}{2}}}}{2\Delta
x}\right]{\mathcal{C}}_{\alpha}\left({\underline{M^{-}_{i}}}\right),\end{split}$
(90)

which follows from conditions (88) and (89), as well as the property of convex
functions shown in lemma 5.

Furthermore, we note that density, pressure, and modified kurtosis ($\rho$,
$p$, and $k$) are all convex functions of
${\underline{q}}=\left({\mathsf{M}}_{0},{\mathsf{M}}_{1},{\mathsf{M}}_{2},{\mathsf{M}}_{3},{\mathsf{M}}_{4}\right)$;
the density is trivially convex, while the pressure is convex if $\rho>0$, and
the modified kurtosis is convex if $\rho>0$ and $p>0$ (see definitions (11)).
Therefore, to prove result (81), all we need to show is that

${\mathcal{C}}_{\alpha}\left({\underline{M^{+}_{i}}}\right)>0\qquad\text{and}\qquad{\mathcal{C}}_{\alpha}\left({\underline{M^{-}_{i}}}\right)>0,$
(91)

with ${\mathcal{C}}_{\alpha}$ chosen as the density, pressure, and modified
kurtosis:

$\begin{gathered}{\mathcal{C}}_{\rho}\left({\underline{q}}\right):=q_{1},\quad{\mathcal{C}}_{p}\left({\underline{q}}\right):=q_{3}-\frac{q_{2}^{2}}{q_{1}},\\\
{\mathcal{C}}_{k}\left({\underline{q}}\right):=\frac{q_{3}^{3}-2q_{2}q_{3}q_{4}+q_{1}q_{4}^{2}+q_{2}^{2}q_{5}-q_{1}q_{3}q_{5}}{q_{2}^{2}-q_{1}q_{3}}.\end{gathered}$
(92)

  1. 1.

Density: We take $\alpha\equiv\rho$ and note that

${\mathcal{C}}_{\rho}\left({\underline{Q}}\pm\
\lambda^{-1}{\underline{f}}({\underline{Q}})\right)=\frac{\left(\lambda\pm
u\right)\rho}{\lambda}.$ (93)

Positivity of (93) follows from the fact that the wave speed, $\lambda$,
always exceeds the local fluid speed, $|u|$.

  2. 2.

Pressure: We take $\alpha\equiv p$ and note that

${\mathcal{C}}_{p}\left({\underline{Q}}\pm{\lambda^{-1}}{\underline{f}}({\underline{Q}})\right)=\frac{p\rho\left(\lambda\pm
u\right)^{2}+h\rho\left(u\pm\lambda\right)-p^{2}}{\lambda\left(\lambda\pm
u\right)\rho},$ (94)

for which we note that the numerator is a quadratic polynomial in $\lambda$.
The roots of this quadratic polynomial can easily be computed:

$\begin{split}z_{1}=\mp\left(u+\frac{h}{2p}\right)-\sqrt{\frac{p}{\rho}+\left(\frac{h}{2p}\right)^{2}},\quad
z_{2}=\mp\left(u+\frac{h}{2p}\right)+\sqrt{\frac{p}{\rho}+\left(\frac{h}{2p}\right)^{2}}.\end{split}$

Positivity of (94) follows from the fact that
$\lambda>\max\left\\{z_{1},\,z_{2}\right\\}$ (i.e., $\lambda$ is always to the
right of the roots) and $p\rho>0$ (i.e., the quadratic is concave up).

  3. 3.

Modified kurtosis: We take $\alpha\equiv k$ and note that

${\mathcal{C}}_{k}\left({\underline{Q}}\pm{\lambda^{-1}}{\underline{f}}({\underline{Q}})\right)=\left(\frac{k(\lambda\pm
u)}{\lambda}\right)\left(\frac{p\rho\left(\lambda\pm
u\right)^{2}+h\rho\left(u\pm\lambda\right)-p^{2}-k\rho}{p\rho\left(\lambda\pm
u\right)^{2}+h\rho\left(u\pm\lambda\right)-p^{2}}\right),$ (95)

where the numerator of the second fraction is again a quadratic polynomial in
$\lambda$. Showing that this quadratic is positive is sufficient to show that
the whole expression is positive since the remaining pieces are already
positive due to the previously established positivity of (93) and (94). The
roots of the quadratic are:

$\begin{split}z_{1}=\mp\left(u+\frac{h}{2p}\right)-\sqrt{\frac{k}{p}+\frac{p}{\rho}+\left(\frac{h}{2p}\right)^{2}},\quad
z_{2}=\mp\left(u+\frac{h}{2p}\right)+\sqrt{\frac{k}{p}+\frac{p}{\rho}+\left(\frac{h}{2p}\right)^{2}}.\end{split}$

Positivity of (95) follows from the fact that
$\lambda>\max\left\\{z_{1},\,z_{2}\right\\}$ (i.e., $\lambda$ is always to the
right of the roots) and $p\rho>0$ (i.e., the quadratic is concave up).

$\blacksquare$

We have now shown that the first-order method will maintain positivity from
one time step to the next under an appropriate CFL condition. However, higher-
order methods will not automatically guarantee positivity; we address this
issue in the subsequent subsections.

###  5.2 Limiter I: Positivity-at-points in the prediction step

The prediction step, as described in (73) requires numerical quadrature in
space-time in each Picard iteration. Furthermore, once the predicted solution
has been computed, it will again be integrated in space-time in the correction
step (i.e., see (76)–(79)). To guarantee that all of the numerical quadratures
in both the prediction and correction steps only use positive values of
density, pressure, and modified kurtosis, we introduce a prediction-step
positivity limiter.

Let the 1D Gauss-Legendre points internal to each element, augmented by the
element end-points (i.e., the location of the element faces), be defined as
follows:

${\mathbb{X}}_{{M_{\text{O}}}}:=\Bigl{\\{}-1,1\Bigr{\\}}\cup\Bigl{\\{}\text{roots
of the ${M_{\text{O}}^{\text{th}}}$ degree Legendre polynomial}\Bigr{\\}},$
(96)

where ${M_{\text{O}}}$ is the desired order of accuracy. Note that
${\mathbb{X}}_{{M_{\text{O}}}}$ contains a total of ${M_{\text{O}}}+2$ points.
We note that all of the quadratures in the prediction update (73) and
correction update (76)–(79) only depends on the predicted solution at a small
number of quadrature points, which are fully contained in the Cartesian
product of ${\mathbb{X}}_{{M_{\text{O}}}}$ with itself:

${\mathbb{X}}^{2}_{{M_{\text{O}}}}:={\mathbb{X}}_{{M_{\text{O}}}}\otimes{\mathbb{X}}_{{M_{\text{O}}}}.$
(97)

Therefore, ${\mathbb{X}}^{2}_{{M_{\text{O}}}}$ contains a total of
$({M_{\text{O}}}+2)^{2}$ points. Our goal is thus to enforce positivity at all
the space-time points $(\tau,\xi)\in{\mathbb{X}}^{2}_{{M_{\text{O}}}}$.

Following the strategy developed by Zhang and Shu article:ZhangShu11  for the
Runge-Kutta discontinuous Galerkin scheme, we apply the following procedure,
which is applied element-by-element.

Step 1.

On the current space-time element defined over

$(t,x)\in\left[t^{n},\,t^{n}+\Delta t\right]\times\left[x_{i}-\frac{\Delta
x}{2},\,x_{i}+\frac{\Delta x}{2}\right],$

the solution is given by (65). Find the minimum density, pressure, and
modified kurtosis of this solution over the points
$(\tau,\xi)\in{\mathbb{X}}^{2}_{{M_{\text{O}}}}$:

$\alpha_{i}^{(m)}:=\min_{(\tau,\xi)\in{\mathbb{X}}^{2}_{{M_{\text{O}}}}}\left\\{{\underline{\Psi}}\left(\tau,\xi\right)^{T}{\underline{W^{n+\frac{1}{2}}_{i\,(:,m)}}}\right\\},$
(98)

for $m=1$ (density), $m=3$ (pressure), and $m=5$ (modified kurtosis).

Step 2.

Rewrite the solution as

$\begin{split}{\underline{\alpha}}^{\text{ST}}\left(t^{n}+\frac{\Delta
t}{2}(1+\tau),\,x_{i}+\frac{\Delta
x}{2}\,\xi;\,\theta\right):=(1-\theta)\,{\underline{W^{n+\frac{1}{2}}_{i\,(1,:)}}}+\theta\,{\underline{\Psi}}\left(\tau,\xi\right)^{T}{\underline{\underline{W^{n+\frac{1}{2}}_{i}}}},\end{split}$
(99)

where $\theta\in[0,1]$, such that $\theta=1$ recovers the original solution
(65) and $\theta=0$ results in reducing the entire solution on to its space-
time average. We now choose the largest possible $\theta\in[0,1]$ so that (99)
is positive for components
$m=1\,(\text{density}),\,3\,(\text{pressure}),\,5\,(\text{modified kurtosis})$
at all the space-time points in ${\mathbb{X}}^{2}_{{M_{\text{O}}}}$. This is
achieved by taking

$\theta=\min_{m\in\left\\{1,3,5\right\\}}\min\left\\{1,\frac{{W^{n+\frac{1}{2}}_{i\,(1,m)}}-\varepsilon}{{W^{n+\frac{1}{2}}_{i\,(1,m)}}-\alpha_{i}^{(m)}}\right\\},$
(100)

where $\varepsilon>0$ is a preselected small constant (e.g., in this work we
select $\varepsilon=10^{-14}$).

###  5.3 Limiter II: Positivity-in-the-mean in the correction step

One of the key challenges in the correction step, as described by (76)–(79) is
to make sure that the element averages of the density, pressure, and modified
kurtosis remain positive at the end of the time-step: $\mkern
1.5mu\overline{\mkern-1.5mu\rho\mkern-1.5mu}\mkern 1.5mu_{i}^{n+1}>0$, $\mkern
1.5mu\overline{\mkern-1.5mup\mkern-1.5mu}\mkern 1.5mu_{i}^{n+1}>0$, and
$\mkern 1.5mu\overline{\mkern-1.5muk\mkern-1.5mu}\mkern 1.5mu_{i}^{n+1}>0$,
where the bar over each variable refers to the element average. The prediction
step limiter described in the previous subsection, §5.2, helps with this
positivity-in-the-mean but cannot guarantee it. Furthermore, without
positivity-in-the-mean, we cannot achieve positivity of the higher-order
polynomial inside the element (i.e., if the polynomial average is negative, a
significant portion of the polynomial must be negative inside the element). To
overcome this challenge, we extend the approach developed by Moe et al.
article:MoRoSe17 , which, for the element averages, blends the high-order
update described by (76)–(79) with a first-order Rusanov scheme. We have
already proved in Theorem 5.1 that the Rusanov scheme is guaranteed to
preserve positivity.

We begin by defining the Rusanov article:Ru61  (aka local Lax-Friedrichs)
update based on the element averages at $t=t^{n}$:

${\underline{Q_{i}^{\text{Rus}}}}\,:=\,{\underline{Q^{n}_{i(1,:)}}}-\frac{\Delta
t}{\Delta
x}\left(\,{\underline{{\mathcal{F}}_{i+\frac{1}{2}}^{\text{Rus}}}}-{\underline{{\mathcal{F}}_{i-\frac{1}{2}}^{\text{Rus}}}}\,\right),$
(101)

where the numerical flux is given by

${\underline{{\mathcal{F}}_{i-\frac{1}{2}}^{\text{Rus}}}}:=\frac{1}{2}\left[{\underline{f}}\left({\underline{Q^{n}_{i(1,:)}}}\right)+{\underline{f}}\left({\underline{Q^{n}_{i-1(1,:)}}}\right)\right]-\frac{1}{2}\Bigl{|}{\lambda_{i-\frac{1}{2}}}\Bigr{|}\left({\underline{Q^{n}_{i(1,:)}}}-{\underline{Q^{n}_{i-1(1,:)}}}\right),$
(102)

and $\Bigl{|}{\lambda_{i-\frac{1}{2}}}\Bigr{|}$ is a local bound of the flux
Jacobian spectral radius. Recall that ${\underline{Q_{i}^{\text{Rus}}}}$ is
guaranteed to have positive density, pressure, and modified kurtosis under a
time-step restriction (see Theorem 5.1).

Next, we write the update for the element averages of the high-order method in
terms of the low-order update (101):

${\underline{Q^{n+1}_{i\,(1,:)}}}={\underline{Q^{\text{Rus}}_{i\,(1,:)}}}-\frac{\Delta
t}{\Delta
x}\left(\,\theta_{i+\frac{1}{2}}\,{\underline{\Delta{\mathcal{F}}_{i+\frac{1}{2}}}}-\theta_{i-\frac{1}{2}}\,{\underline{\Delta{\mathcal{F}}_{i-\frac{1}{2}}}}\,\right),$
(103)

where the difference between the high and low-order fluxes is given by

${\underline{\Delta{\mathcal{F}}_{i-\frac{1}{2}}}}:={\underline{{\mathcal{F}}^{n+\frac{1}{2}}_{i-\frac{1}{2}}}}-{\underline{{\mathcal{F}}^{\text{Rus}}_{i-\frac{1}{2}}}},$
(104)

and $\theta_{i+\frac{1}{2}}\in[0,1]$ measures the amount of flux limiting,
where $\theta_{i+\frac{1}{2}}=0$ represents maximal limiting (i.e., no high-
order flux contributions) and $\theta_{i+\frac{1}{2}}=1$ represents no
limiting (i.e., no low-order flux contributions).

The strategy for the positivity-in-the-mean limiter is then to find the
maximum $\theta_{i+\frac{1}{2}}\in[0,1]$ such that $\forall i$

${\mathcal{C}}_{\alpha}\left({\underline{Q^{n+1}_{i\,(1,:)}}}\right)>0,$ (105)

for $\alpha\equiv\rho$ (density), $\alpha\equiv p$ (pressure), and
$\alpha\equiv k$ (modified kurtosis), as defined in (92). The strategy for
achieving this is outlined below and is applied element-by-element. The
process begins by initializing $\theta_{i+\frac{1}{2}}=1$ $\forall i$.

Step 1: (density)

Define

$\Gamma:=\frac{\Delta x}{\Delta
t}\left({Q^{\text{Rus}}_{i(1)}-\varepsilon}\right).$ (106)

Set $\Lambda_{\text{left}}=\Lambda_{\text{right}}=1$ (i.e., full high-order
flux), but modify these if there is any potential for the density to decrease
below zero.

Case 1.

If $\Delta{\mathcal{F}}_{i-\frac{1}{2}(1)}<0$ and
$\Delta{\mathcal{F}}_{i+\frac{1}{2}(1)}<0$, then

$\Lambda_{\text{left}}=\Lambda_{\text{right}}=\min\left\\{1,\frac{\Gamma}{\left|\Delta{\mathcal{F}}_{i-\frac{1}{2}(1)}\right|+\left|\Delta{\mathcal{F}}_{i+\frac{1}{2}(1)}\right|}\right\\}.$
(107)

Case 2.

If $\Delta{\mathcal{F}}_{i-\frac{1}{2}(1)}<0$ and
$\Delta{\mathcal{F}}_{i+\frac{1}{2}(1)}>0$ then

$\Lambda_{\text{left}}=\min\left\\{1,\frac{\Gamma}{\left|\Delta{\mathcal{F}}_{i-\frac{1}{2}(1)}\right|}\right\\}.$
(108)

Case 3.

If $\Delta{\mathcal{F}}_{i+\frac{1}{2}(1)}<0$ and
$\Delta{\mathcal{F}}_{i-\frac{1}{2}(1)}>0$ then

$\Lambda_{\text{right}}=\min\left\\{1,\frac{\Gamma}{\left|\Delta{\mathcal{F}}_{i+\frac{1}{2}(1)}\right|}\right\\}.$
(109)

Step 2: (pressure)

Compute the Rusanov pressure (which is guaranteed to be positive):

$p^{\text{Rus}}:={\mathcal{C}}_{p}\left({\underline{Q^{\text{Rus}}_{i}}}\right).$
(110)

Set $\mu_{11}=\mu_{10}=\mu_{01}=1$, but modify these if there is any potential
for the pressure to decrease below zero.

Part 2A.

Set

${\underline{Q^{\star}}}={\underline{Q^{n}_{i\,(1,:)}}}-\frac{\Delta t}{\Delta
x}\left(\Lambda_{\text{right}}\,{\underline{{\mathcal{F}}_{i+\frac{1}{2}}^{n+\frac{1}{2}}}}-\Lambda_{\text{left}}\,{\underline{{\mathcal{F}}_{i-\frac{1}{2}}^{n+\frac{1}{2}}}}\right),\quad
p^{\star}={\mathcal{C}}_{p}\left({\underline{Q^{\star}}}\right).$ (111)

If $p^{\star}<\varepsilon$, then we set
$\mu_{11}=\left({p^{\text{Rus}}-\varepsilon}\right)\bigl{/}\left(p^{\text{Rus}}-p^{\star}\right)$.

Part 2B.

Set

${\underline{Q^{\star}}}={\underline{Q^{n}_{i\,(1,:)}}}+\frac{\Delta t}{\Delta
x}\Lambda_{\text{left}}\,{\underline{{\mathcal{F}}_{i-\frac{1}{2}}^{n+\frac{1}{2}}}},\quad
p^{\star}={\mathcal{C}}_{p}\left({\underline{Q^{\star}}}\right).$ (112)

If $p^{\star}<\varepsilon$, then
$\mu_{10}=\left({p^{\text{Rus}}-\varepsilon}\right)/\left(p^{\text{Rus}}-p^{\star}\right)$.

Part 2C.

Set

${\underline{Q^{\star}}}={\underline{Q^{n}_{i\,(1,:)}}}-\frac{\Delta t}{\Delta
x}\Lambda_{\text{right}}\,{\underline{{\mathcal{F}}_{i+\frac{1}{2}}^{n+\frac{1}{2}}}},\quad
p^{\star}={\mathcal{C}}_{p}\left({\underline{Q^{\star}}}\right).$ (113)

If $p^{\star}<\varepsilon$, then
$\mu_{01}=\left(p^{\text{Rus}}-\varepsilon\right)/\left(p^{\text{Rus}}-p^{\star}\right)$.

Part 2D.

Set

$\mu=\min\bigl{\\{}\mu_{11},\mu_{10},\mu_{01}\bigr{\\}},\quad\Lambda_{\text{left}}\leftarrow\mu\,\Lambda_{\text{left}},\quad\Lambda_{\text{right}}\leftarrow\mu\,\Lambda_{\text{right}}.$
(114)

Step 3: (modified kurtosis)

Compute the Rusanov modified kurtosis (which is guaranteed to be positive):

$k^{\text{Rus}}:={\mathcal{C}}_{k}\left({\underline{Q^{\text{Rus}}_{i}}}\right).$
(115)

Set $\mu_{11}=\mu_{10}=\mu_{01}=1$, but modify these if there is any potential
for the pressure to decrease below zero.

Part 3A.

Set

${\underline{Q^{\star}}}={\underline{Q^{n}_{i\,(1,:)}}}-\frac{\Delta t}{\Delta
x}\left(\Lambda_{\text{right}}\,{\underline{{\mathcal{F}}_{i+\frac{1}{2}}^{n+\frac{1}{2}}}}-\Lambda_{\text{left}}\,{\underline{{\mathcal{F}}_{i-\frac{1}{2}}^{n+\frac{1}{2}}}}\right),\quad
k^{\star}={\mathcal{C}}_{k}\left({\underline{Q^{\star}}}\right).$ (116)

If $k^{\star}<\varepsilon$, then we set
$\mu_{11}=\left({k^{\text{Rus}}-\varepsilon}\right)\bigl{/}\left(k^{\text{Rus}}-k^{\star}\right)$.

Part 3B.

Set

${\underline{Q^{\star}}}={\underline{Q^{n}_{i\,(1,:)}}}+\frac{\Delta t}{\Delta
x}\Lambda_{\text{left}}\,{\underline{{\mathcal{F}}_{i-\frac{1}{2}}^{n+\frac{1}{2}}}},\quad
k^{\star}={\mathcal{C}}_{k}\left({\underline{Q^{\star}}}\right).$ (117)

If $k^{\star}<\varepsilon$, then
$\mu_{10}=\left({k^{\text{Rus}}-\varepsilon}\right)/\left(k^{\text{Rus}}-k^{\star}\right)$.

Part 3C.

Set

${\underline{Q^{\star}}}={\underline{Q^{n}_{i\,(1,:)}}}-\frac{\Delta t}{\Delta
x}\Lambda_{\text{right}}\,{\underline{{\mathcal{F}}_{i+\frac{1}{2}}^{n+\frac{1}{2}}}},\quad
k^{\star}={\mathcal{C}}_{k}\left({\underline{Q^{\star}}}\right).$ (118)

If $k^{\star}<\varepsilon$, then
$\mu_{01}=\left(p^{\text{Rus}}-\varepsilon\right)/\left(k^{\text{Rus}}-k^{\star}\right)$.

Part 3D.

Set

$\mu=\min\bigl{\\{}\mu_{11},\mu_{10},\mu_{01}\bigr{\\}},\quad\Lambda_{\text{left}}\leftarrow\mu\,\Lambda_{\text{left}},\quad\Lambda_{\text{right}}\leftarrow\mu\,\Lambda_{\text{right}}.$
(119)

Step 4:

Set

$\theta_{i-\frac{1}{2}}\leftarrow\min\left\\{\theta_{i-\frac{1}{2}},\,\Lambda_{\text{left}}\right\\}\quad\text{and}\quad\theta_{i+\frac{1}{2}}\leftarrow\min\left\\{\theta_{i+\frac{1}{2}},\,\Lambda_{\text{right}}\right\\}.$
(120)

In all of the above formulas, we select in this work: $\varepsilon=10^{-14}$.

###  5.4 Limiter III: Positivity-at-points in the correction step

Once we have ensured that the element averages are positive, we then look to
enforce positivity of the corrected solution at spatial quadrature points:
$\xi\in{\mathbb{X}}_{{M_{\text{O}}}}$ as defined by (96). Following the ideas
developed by Zhang and Shu article:ZhangShu11  for the Runge-Kutta
discontinuous Galerkin scheme, we aim to find the maximum $\theta\in[0,1]$
such that

${\underline{q}}^{h}\left(t^{n+1},\,x_{i}+\frac{\Delta
x}{2}\,\xi;\theta\right):=\left(1-\theta\right)\,{\underline{Q^{n+1}_{i\,(1,:)}}}+\theta\,{\underline{\Phi}}(\xi)^{T}\,{\underline{\underline{Q^{n+1}_{i}}}}$
(121)

is positive at all points $\xi\in{\mathbb{X}}_{{M_{\text{O}}}}$ for every
space element ${\mathcal{T}}_{i}$. As in the prediction step limiter from
§5.2, $\theta=0$ means that the solution is limited fully down to its element
average, while $\theta=1$ means that no limiting is needed and the full high-
order approximation can be used. We apply the following procedure element-by-
element.

Step 1.

On the current element defined over

$x\in\left[x_{i}-\frac{\Delta x}{2},\,x_{i}+\frac{\Delta x}{2}\right],$

the solution is given by (59). Find the minimum density
$\forall\xi\in{\mathbb{X}}_{{M_{\text{O}}}}$ (see (96)) and compute the
corresponding damping parameter ($\theta$):

$\rho_{i}^{\text{min}}:=\min_{\xi\in{\mathbb{X}}_{{M_{\text{O}}}}}\left\\{{\underline{\Phi}}\left(\xi\right)^{T}{\underline{Q^{n+1}_{i\,(:,1)}}}\right\\},\quad\theta=\min\left\\{1,\frac{{Q^{n+1}_{i\,(1,1)}}-\varepsilon}{{Q^{n+1}_{i\,(1,1)}}-\rho_{i}^{{\text{min}}}}\right\\}.$
(122)

Finally, rescale the higher-order coefficients using the above calculated
damping parameter ($\theta$):

${\underline{Q^{n+1}_{i\,(:,\ell)}}}\leftarrow\theta\,{\underline{Q^{n+1}_{i\,(:,\ell)}}}\quad\forall\ell=2,\ldots,{M_{\text{C}}}.$
(123)

Step 2.

Now that the density is positive $\forall\xi\in{\mathbb{X}}_{{M_{\text{O}}}}$,
we repeat Step 1 for the pressure. That is, we find the average pressure and
the minimum pressure $\forall\xi\in{\mathbb{X}}_{{M_{\text{O}}}}$:

$\mkern 1.5mu\overline{\mkern-1.5muP\mkern-1.5mu}\mkern
1.5mu_{i}:={\mathcal{C}}_{p}\left({\underline{Q^{n+1}_{i\,(1,:)}}}\right),\quad
p_{i}^{\text{min}}:=\min_{\xi\in{\mathbb{X}}_{{M_{\text{O}}}}}\left\\{{\mathcal{C}}_{p}\left({\underline{\Phi}}\left(\xi\right)^{T}{\underline{\underline{Q^{n+1}_{i}}}}\right)\right\\},$
(124)

where ${\mathcal{C}}_{p}\left({\underline{q}}\right)$ is defined in (92). From
here, we compute the corresponding damping parameter and rescale the higher-
order coefficients:

$\theta=\min\left\\{1,\frac{\mkern
1.5mu\overline{\mkern-1.5muP\mkern-1.5mu}\mkern 1.5mu_{i}-\varepsilon}{\mkern
1.5mu\overline{\mkern-1.5muP\mkern-1.5mu}\mkern
1.5mu_{i}-p_{i}^{{\text{min}}}}\right\\},\quad{\underline{Q^{n+1}_{i\,(:,\ell)}}}\leftarrow\theta\,{\underline{Q^{n+1}_{i\,(:,\ell)}}}\quad\forall\ell=2,\ldots,{M_{\text{C}}}.$
(125)

Step 3.

Now that both density and pressure are positive
$\forall\xi\in{\mathbb{X}}_{{M_{\text{O}}}}$, we repeat Step 2 for the
modified kurtosis. That is, we find the average modified kurtosis and the
minimum modified kurtosis $\forall\xi\in{\mathbb{X}}_{{M_{\text{O}}}}$:

$\mkern 1.5mu\overline{\mkern-1.5muK\mkern-1.5mu}\mkern
1.5mu_{i}:={\mathcal{C}}_{k}\left({\underline{Q^{n+1}_{i\,(1,:)}}}\right),\quad
k_{i}^{\text{min}}:=\min_{\xi\in{\mathbb{X}}_{{M_{\text{O}}}}}\left\\{{\mathcal{C}}_{k}\left({\underline{\Phi}}\left(\xi\right)^{T}{\underline{\underline{Q^{n+1}_{i}}}}\right)\right\\},$
(126)

where ${\mathcal{C}}_{k}\left({\underline{q}}\right)$ is defined in (92). From
here, we compute the corresponding damping parameter and rescale the higher-
order coefficients:

$\theta=\min\left\\{1,\frac{\mkern
1.5mu\overline{\mkern-1.5muK\mkern-1.5mu}\mkern 1.5mu_{i}-\varepsilon}{\mkern
1.5mu\overline{\mkern-1.5muK\mkern-1.5mu}\mkern
1.5mu_{i}-k_{i}^{{\text{min}}}}\right\\},\quad{\underline{Q^{n+1}_{i\,(:,\ell)}}}\leftarrow\theta\,{\underline{Q^{n+1}_{i\,(:,\ell)}}}\quad\forall\ell=2,\ldots,{M_{\text{C}}}.$
(127)

In all the formulas presented above we use $\varepsilon=10^{-14}$.

###  5.5 Limiter IV: Unphysical oscillation limiter

The previously described limiters guarantee positivity for $\rho$, $p$, and
$k$, but there still may be unphysical oscillations near shocks, rarefactions,
or large gradients. We augment the method with one more limiter to eliminate
these oscillations: a variant of the strategy developed in Moe et al.
article:Moe15 . This limiter is applied once per time step and can remove
unphysical oscillations without overly diffusing the numerical solution. We
apply the following procedure.

Step 1.

Loop over each element ${\mathcal{T}}_{i}$ and compute the minimum and maximum
values of all of the following variables:
$w^{\ell}\in\left\\{\rho,u,p,h,r\right\\}$:

$w^{\ell}_{M_{i}}:=\max_{\xi\in{\mathbb{X}}_{{M_{\text{O}}}}}\left\\{w^{\ell}\left({\underline{q^{\Delta
x}}}(\xi)\right)\Bigl{|}_{{\mathcal{T}}_{i}}\right\\},\quad
w^{\ell}_{m_{i}}:=\min_{\xi\in{\mathbb{X}}_{{M_{\text{O}}}}}\left\\{w^{\ell}\left({\underline{q^{\Delta
x}}}(\xi)\right)\Bigl{|}_{{\mathcal{T}}_{i}}\right\\},$ (128)

for all $\ell=1,2,3,4,5$. Here ${\mathbb{X}}_{{M_{\text{O}}}}$ is taken to be
the ${M_{\text{O}}}$ roots of the ${M_{\text{O}}}^{\text{th}}$ Legendre
polynomial (i.e., Gauss-Legendre points) plus the element ends points (see eq.
96).

Step 2.

Compute upper and lower bounds over all neighborhoods,
$N_{{\mathcal{T}}_{i}}:=\left\\{{\mathcal{T}}_{i-1},{\mathcal{T}}_{i},\,{\mathcal{T}}_{i+1}\right\\}$:

$\begin{split}{M^{\ell}_{i}}&={\max\left\\{\mkern
1.5mu\overline{\mkern-1.5muw\mkern-1.5mu}\mkern
1.5mu^{\ell}_{i}+{\mathcal{A}}_{0}h^{1.5},\,\displaystyle\max_{j\in
N_{{\mathcal{T}}_{i}}}\left\\{w^{\ell}_{M_{j}}\right\\}\right\\}},\\\
{m^{\ell}_{i}}&={\min\left\\{\mkern
1.5mu\overline{\mkern-1.5muw\mkern-1.5mu}\mkern
1.5mu^{\ell}_{i}-{\mathcal{A}}_{0}h^{1.5},\,\displaystyle\min_{j\in
N_{{\mathcal{T}}_{i}}}\left\\{w^{\ell}_{m_{j}}\right\\}\right\\}},\end{split}$
(129)

where $\mkern 1.5mu\overline{\mkern-1.5muw\mkern-1.5mu}\mkern
1.5mu^{\ell}_{i}$ are the element averages of each variable, and
${\mathcal{A}}_{0}h^{1.5}$ is used to offset these averages to recover high-
order accuracy for smooth solutions in the limit $h\rightarrow 0$ (see Moe et
al. article:Moe15  for more details).

Step 3.

On each element ${\mathcal{T}}_{i}$, compute the largest damping parameters
between $[0,1]$ that guarantee that the high-order solution in
${\mathcal{T}}_{i}$ does not violate the maximum and minimum bounds defined by
(129):

$\theta=\min\left\\{1,\,\mu\cdot\min_{\ell}\left\\{\frac{M^{\ell}_{i}-\bar{w}^{\ell}_{i}}{w^{\ell}_{M_{i}}-\bar{w}^{\ell}_{i}}\right\\},\,\mu\cdot\min_{\ell}\left\\{\frac{m^{\ell}_{i}-\bar{w}^{\ell}_{i}}{w^{\ell}_{m_{i}}-\bar{w}^{\ell}_{i}}\right\\}\right\\},$
(130)

where the factor $\mu=10/11$ is introduced to slightly increase the
aggressiveness of the limiter (again, see Moe et al. article:Moe15  for more
details).

Step 4.

On each element ${\mathcal{T}}_{i}$, limit the conserved variables:

${\underline{Q^{n+1}_{i\,(\ell,:)}}}\leftarrow\theta\,{\underline{Q^{n+1}_{i\,(\ell,:)}}}\quad\forall\ell=2,\ldots,{M_{\text{C}}}.$
(131)

##  6 Collisionless HyQMOM numerical examples

In this section, we apply the proposed scheme and the corresponding limiters
to several test cases. In §6.1, we verify the claimed orders of accuracy on a
smooth exact solution. In §6.2 and §6.3 we apply the scheme to shock tube
initial data. These results demonstrate the ability of the non-oscillatory
limiter to control unphysical oscillations. Finally, in §6.4, we fully
validate the positivity limiters by applying the scheme to piecewise constant
initial data that lead to the formation of a vacuum. This example demonstrates
the ability of the positivity limiters to prevent negative states in density,
pressure, and modified kurtosis, both on the element average and on the
solution values internal to the element. In all the cases presented in §6.2,
§6.3, and §6.4, we compare the high-order scheme against a highly-resolved
first-order Rusanov scheme that is guaranteed to be positivity-preserving
without the need for any limiters.

###  6.1 Smooth solution convergence test

Consider the following exact solution to the 1D HyQMOM system (37)–(38) with
periodic boundary conditions on $x\in\bigl{[}-1,1\bigr{]}$:

$\begin{split}\rho(t,x)&=2+\sin(2\pi(x-t)),\\\ u(t,x)=1,\quad p(t,x)&=2,\quad
h(t,x)=4,\quad k(t,x)=8-4\left[{\rho(t,x)}\right]^{-1}.\end{split}$ (132)

The numerical solution is computed with grid resolutions of

${M_{\text{elem}}}=10\times 2^{\ell},\quad\text{for}\quad\ell=0,1,2,3,4,5,$
(133)

up to a final time of $t=1$. We verify the order of accuracy for the schemes
with orders of accuracy ${M_{\text{O}}}=2,3,4$.

The errors we report are based on the following error measure:

$\sum_{\ell=1}^{{M_{\text{eqn}}}}\frac{\Bigl{\|}f_{\ell}-g_{\ell}\Bigr{\|}_{\scriptstyle
L^{2}[-1,1]}}{\Bigl{\|}g_{\ell}\Bigr{\|}_{\scriptstyle
L^{2}[-1,1]}}=\sum_{\ell=1}^{{M_{\text{eqn}}}}\sqrt{\frac{\displaystyle\int_{-1}^{1}\Bigl{|}f_{\ell}(x)-g_{\ell}(x)\Bigr{|}^{2}\,dx}{\displaystyle\int_{-1}^{1}\Bigl{|}g_{\ell}(x)\Bigr{|}^{2}\,dx}},$
(134)

where ${\underline{f}}(x):[-1,1]\mapsto\mathbb{R}^{{M_{\text{eqn}}}}$ is the
approximate solution and
${\underline{g}}(x):[-1,1]\mapsto\mathbb{R}^{{M_{\text{eqn}}}}$ is the exact
solution. In practice, however, we replace the exact solution with a piecewise
Legendre polynomial approximation of degree ${M_{\text{O}}}+1$ on the
computational mesh. Repeated use of the orthonormality of the Legendre basis
functions yields the following (approximate) relative error on a mesh with $N$
elements and a numerical method of order ${M_{\text{O}}}$:

$e_{N}:=\sum_{\ell=1}^{{M_{\text{eqn}}}}\sqrt{\frac{\displaystyle\sum_{i=1}^{N}\left\\{\sum_{j=1}^{M_{C}}\Bigl{(}Q_{i(j,\ell)}-Q^{\star}_{i(j,\ell)}\Bigr{)}^{2}+\Bigl{(}Q^{\star}_{i(M_{C}+1,\ell)}\Bigr{)}^{2}\right\\}}{\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{M_{C}+1}\Bigl{(}Q^{\star}_{i(j,\ell)}\Bigr{)}^{2}}},$
(135)

where $Q$ and $Q^{\star}$ are the Legendre coefficients of the numerical and
exact solutions at the final time, respectively. The exact solution
coefficients are computed using Gaussian quadrature with $20$ quadrature
points per element:

${\underline{Q^{\star}_{i(k,:)}}}:=\frac{1}{2}\sum_{a=1}^{20}\omega_{a}^{\star}\,\phi_{k}\left(\mu_{a}^{\star}\right){\underline{q}}^{\star}\left(t=1,\,x_{i}+\frac{\Delta
x}{2}\mu_{a}^{\star}\right),$ (136)

where $\omega_{a}^{\star}$ and $\mu_{a}^{\star}$ for $a=1,\dots,20$ are the
weights and abscissas of the $20$ point quadrature rule, and
${\underline{q}}^{\star}$ is the exact solution. Gaussian quadrature rules
have been tabulated in many books and websites; we obtained our data from
web:gauss_quad .

The errors as defined by (135), as well as the base-2 logarithms of the ratio
of consecutive errors,

$\log_{2}\left(\frac{e_{N/2}}{e_{N}}\right)\approx\log_{2}\left(\frac{\left(N/2\right)^{-{M_{\text{O}}}}}{N^{-{M_{\text{O}}}}}\right)=\log_{2}\left(2^{{M_{\text{O}}}}\right)={M_{\text{O}}},$
(137)

are shown in Tables 2 (all limiters are turned off) and 3 (all limiters are
turned on). For the simulations that result in Table 3, none of the three
positivity limiters (i.e., Limiters I, II, and III) are active because the
solution is far away from positivity violations. In Table 3, the values
affected by Limiter IV are highlighted in red. Note that at low resolutions,
Limiter IV is active, and the results in Tables 2 and 3 differ slightly, but
that at higher resolutions, the effect of the limiter disappears. Note that
for all the simulations with limiters turned on, we used the value of
${\mathcal{A}}_{0}=5$ in formula eq. 129.

$N$ | $M_{\text{O}}=2$ | Eq. (137) | $M_{\text{O}}=3$ | Eq. (137) | $M_{\text{O}}=4$ | Eq. (137)
---|---|---|---|---|---|---
10 | 1.143e-01 | – | 1.171e-02 | – | 4.924e-03 | –
20 | 2.005e-02 | $2.511$ | 2.260e-03 | $2.374$ | 4.617e-04 | $3.415$
40 | 3.759e-03 | $2.415$ | 4.032e-04 | $2.487$ | 5.337e-06 | $6.435$
80 | 8.802e-04 | $2.095$ | 6.077e-05 | $2.730$ | 1.962e-07 | $4.765$
160 | 2.192e-04 | $2.006$ | 8.127e-06 | $2.903$ | 1.203e-08 | $4.028$
320 | 5.485e-05 | $1.999$ | 1.040e-06 | $2.966$ | 7.500e-10 | $4.004$
Table 2: (§6.1: smooth solution convergence test with limiters turned off) Relative $L^{2}$ errors for the HyQMOM equations with variable density, constant fluid velocity, pressure, heat flux, and fourth primitive moment, and periodic boundary conditions. In these simulations all four limiters were turned off.  $N$ | $M_{\text{O}}=2$ | Eq. (137) | $M_{\text{O}}=3$ | Eq. (137) | $M_{\text{O}}=4$ | Eq. (137)
---|---|---|---|---|---|---
10 | 3.154e-01 | – | 5.360e-02 | – | 4.924e-03 | –
20 | 4.887e-02 | $2.690$ | 2.260e-03 | $4.568$ | 4.617e-04 | $3.415$
40 | 3.759e-03 | $3.700$ | 4.032e-04 | $2.487$ | 5.337e-06 | $6.435$
80 | 8.802e-04 | $2.095$ | 6.077e-05 | $2.730$ | 1.962e-07 | $4.765$
160 | 2.192e-04 | $2.006$ | 8.127e-06 | $2.903$ | 1.203e-08 | $4.028$
320 | 5.485e-05 | $1.999$ | 1.040e-06 | $2.966$ | 7.500e-10 | $4.004$
Table 3: (§6.1: smooth solution convergence test with limiters turned on)
Relative $L^{2}$ errors for the one-dimensional HyQMOM equations with variable
density, constant fluid velocity, pressure, heat flux, and fourth primitive
moment, and periodic boundary conditions. In these simulations, all four
limiters were turned on.

###  6.2 Shock tube problem #1

Consider the Riemann problem for (37)–(39) with the following initial data at
$t=0$:

$\bigl{(}\rho,u,p,h,k\bigr{)}(t=0,x)=\begin{cases}\bigl{(}1.5,\,-0.5,\,1.5,\,1.0,\,2.3\bar{3}\bigr{)}\quad
x<0,\\\ \bigl{(}1.0,\,-0.5,\,1.0,\,0.5,\,1.75\bigr{)}\quad x>0,\end{cases}$
(138)

on $x\in[-1.2,1.2]$ with extrapolation boundary conditions.

Shown in Figure 2 are results from a simulation run with two distinct methods:
(1) the ${M_{\text{O}}}=4$ scheme with 200 elements and full limiters (shown
as blue dots), and (2) the first-order Rusanov scheme with 20,000 elements
(shown as a solid red line). For the ${M_{\text{O}}}=4$ scheme, we are
plotting four points per element in order to show the intra-element solution
structure. The panels show the primitive variables: (a) density: $\rho(t,x)$,
(b) macroscopic velocity: $u(t,x)$, (c) pressure: $p(t,x)$, (d) heat flux:
$h(t,x)$, (e) modified kurtosis: $k(t,x)$, and (f) primitive fourth-moment:
$r(t,x)$. Note that we used the value of ${\mathcal{A}}_{0}=5$ in formula eq.
129. These results clearly demonstrate the non-oscillatory limiters’ ability
to adequately control unphysical oscillations and produce accurate solutions.

(a) | (b)
---|---
(c) | (d)
(e) | (f)
Figure 2: (§6.2: shock tube problem #1) Numerical solution of shock tube
problem #1 on $x\in[-1.2,1.2]$ with initial conditions given by (138). Shown
are the results from a simulation run with two distinct methods: (1) the
${M_{\text{O}}}=4$ scheme with 200 elements and full limiters (shown as blue
dots), and (2) the first-order Rusanov scheme with 20,000 elements (shown as a
solid red line). For the ${M_{\text{O}}}=4$ scheme, we are plotting four
points per element in order to show the intra-element solution structure. The
panels show the primitive variables: (a) density: $\rho(t,x)$, (b) macroscopic
velocity: $u(t,x)$, (c) pressure: $p(t,x)$, (d) heat flux: $h(t,x)$, (e)
modified kurtosis: $k(t,x)$, and (f) primitive fourth-moment: $r(t,x)$.

###  6.3 Shock tube problem #2

Consider the Riemann problem for (37)–(39) with the following initial data at
$t=0$:

$\bigl{(}\rho,u,p,h,k\bigr{)}(t=0,x)=\begin{cases}\bigl{(}1.0,\,-0.7,\,1.5,\,1.5,\,1.75\bigr{)}&\quad
x<0,\\\ \bigl{(}0.5,\,-0.9,\,1.0,\,1.0,\,1.0\,\bigr{)}&\quad x>0,\end{cases}$
(139)

on $x\in[-1.2,1.2]$ with extrapolation boundary conditions.

Shown in Figure 3 are results from a simulation run with two distinct methods:
(1) the ${M_{\text{O}}}=4$ scheme with 200 elements and full limiters (shown
as blue dots), and (2) the first-order Rusanov scheme with 20,000 elements
(shown as a solid red line). For the ${M_{\text{O}}}=4$ scheme, we are
plotting four points-per-element in order to show the intra-element solution
structure. The panels show the primitive variables: (a) density: $\rho(t,x)$,
(b) macroscopic velocity: $u(t,x)$, (c) pressure: $p(t,x)$, (d) heat flux:
$h(t,x)$, (e) modified kurtosis: $k(t,x)$, and (f) primitive fourth-moment:
$r(t,x)$. Note that we used the value of ${\mathcal{A}}_{0}=5$ in formula eq.
129.

Again, just as in the previous example, these results demonstrate the ability
of the non-oscillatory limiters to adequately control unphysical oscillations
and produce accurate solutions.

(a) | (b)
---|---
(c) | (d)
(e) | (f)
Figure 3: (§6.3: shock tube problem #2) Numerical solution of shock tube
problem #2 on $x\in[-1.2,1.2]$ with initial conditions given by (139). Shown
are the results from a simulation run with two distinct methods: (1) the
${M_{\text{O}}}=4$ scheme with 200 elements and full limiters (shown as blue
dots), and (2) the first-order Rusanov scheme with 20,000 elements (shown as a
solid red line). For the ${M_{\text{O}}}=4$ scheme, we are plotting four
points per element in order to show the intra-element solution structure. The
panels show the primitive variables: (a) density: $\rho(t,x)$, (b) macroscopic
velocity: $u(t,x)$, (c) pressure: $p(t,x)$, (d) heat flux: $h(t,x)$, (e)
modified kurtosis: $k(t,x)$, and (f) primitive fourth-moment: $r(t,x)$.

###  6.4 Double rarefaction vacuum problem

In the final example, we solve a vacuum problem where the right and left
initial velocities are large and opposite, creating a vacuum state in the
center of the solution domain. The initial states are

$\bigl{(}\rho,u,p,h,k\bigr{)}(t=0,x)=\begin{cases}\bigl{(}1.0,\,-2.0,\,1.0,\,0.0,\,2.0\bigr{)}&\quad
x<0,\\\ \bigl{(}1.0,\,+2.0,\,1.0,\,0.0,\,2.0\bigr{)}&\quad x>0.\end{cases}$
(140)

The computational domain is $x\in[-1.2,1.2]$ with extrapolation boundary
conditions.

Shown in Figure 4 are results from a simulation run with two distinct methods:
(1) the ${M_{\text{O}}}=4$ scheme with 200 elements and full limiters (shown
as blue dots), and (2) the first-order Rusanov scheme with 20,000 elements
(shown as a solid red line). For the ${M_{\text{O}}}=4$ scheme, we are
plotting four points per element in order to show the intra-element solution
structure. The panels show the primitive variables: (a) density: $\rho(t,x)$,
(b) macroscopic velocity: $u(t,x)$, (c) pressure: $p(t,x)$, (d) heat flux:
$h(t,x)$, (e) modified kurtosis: $k(t,x)$, and (f) primitive fourth-moment:
$r(t,x)$. Note that we used the value of ${\mathcal{A}}_{0}=5$ in formula eq.
129.

We comment on two important findings from this simulation. First, this example
demonstrates the ability of the positivity limiters to prevent negative states
in density, pressure, and modified kurtosis, both on the element average and
the solution values internal to the element. In this simulation, all three
variables, $\rho$, $p$, and $k$, become very small, but all stay strictly
above zero. Because all three remain strictly positive, the moments remain
realizable, and the numerical simulation remains nonlinear stable. Second,
while the simulation results from the ${M_{\text{O}}}=4$ scheme with 200
elements do show some differences in the vacuum region with the highly
resolved Rusanov solution, especially in the density plot shown in Figure
4(a), the solution remains qualitatively correct. We can investigate this
further by increasing the grid resolution; in Figure 5 we show the density
plots at different grid resolutions: (a) $N=200$, (b) $N=400$, (c) $N=800$,
and (b) $N=1600$. These results verify that the differences between the
${M_{\text{O}}}=4$ scheme and the highly resolved Rusanov scheme disappear at
higher resolutions.

(a) | (b)
---|---
(c) | (d)
(e) | (f)
Figure 4: (§6.4: double rarefaction vacuum problem) Numerical solution of the double rarefaction vacuum problem on $x\in[-1.2,1.2]$ with initial conditions given by (140). Shown are the results from a simulation run with two distinct methods: (1) the ${M_{\text{O}}}=4$ scheme with 200 elements and full limiters (shown as blue dots), and (2) the first-order Rusanov scheme with 20,000 elements (shown as a solid red line). For the ${M_{\text{O}}}=4$ scheme, we are plotting four points per element in order to show the intra-element solution structure. The panels show the primitive variables: (a) density: $\rho(t,x)$, (b) macroscopic velocity: $u(t,x)$, (c) pressure: $p(t,x)$, (d) heat flux: $h(t,x)$, (e) modified kurtosis: $k(t,x)$, and (f) primitive fourth-moment: $r(t,x)$. (a) | (b)
---|---
(c) | (d)
Figure 5: (§6.4: double rarefaction vacuum problem) Numerical solution of the
double rarefaction vacuum problem on $x\in[-1.2,1.2]$ with initial conditions
given by (140). Shown are the densities at various grid resolutions: (a)
$N=200$, (b) $N=400$, (c) $N=800$, and (d) $N=1600$. In each panel, we compare
the ${M_{\text{O}}}=4$ scheme with full limiters (shown as blue dots) with the
first-order Rusanov scheme with 20,000 elements (shown as a solid red line).

##  7 Extension to HyQMOM-BGK

Up to this point, we have only considered the HyQMOM approximation applied to
the Vlasov model eq. 5; this allowed us to study the mathematical structure of
HyQMOM and to develop accurate high-order methods and limiters. On the other
hand, the practicality of HyQMOM is not for solving collisionless kinetic
models since, in this regime, it would be far better to directly solve the
Vlasov equation with Lagrangian or semi-Lagrangian approaches. Instead, the
true benefit of the HyQMOM approximation is in the approximation of kinetic
systems near thermodynamic equilibrium – a regime we study in this section.

In this section, we extend the previously developed numerical method to HyQMOM
with a BGK collision operator. Importantly, we develop this extension so that
the resulting HyQMOM solver adheres to the following two key design
parameters:

  1. 1.

The method should remain high-order accurate irrespective of the Knudsen
number: $\varepsilon>0$.

  2. 2.

For fixed mesh parameters (i.e., fixed $\Delta x$ and $\Delta t$), the method
should remain stable in the singular limit: $\varepsilon\rightarrow 0^{+}$.
This property is often referred to as the asymptotic-preserving (AP) property,
and a variety of schemes with this property can be found in the literature
(e.g., see article:Ben08 ; article:Caflish1997 ; article:Coron1991 ;
article:Gabetta1997 ; article:Jin1995 ; article:Jin99 ; article:Jin12 ;
article:Jin1996 ; article:Pieraccini2007 ; article:Xiong2017 ).

The specific approach we detail in this section is novel and directly relies
on the prediction and correction format of the method developed in section 4.

###  7.1 1D1V Boltzmann-BGK equation

Consider the 1D1V Boltzmann-BGK equation article:BGK54 :

$f_{,t}+vf_{,x}=\frac{1}{\varepsilon}\left({\mathcal{M}}-f\right),$ (141)

where $\varepsilon>0$ is the Knudsen number, which is a non-dimensional ratio
of the particle mean-free path to a characteristic length scale, and
${\mathcal{M}}(t,x,v):\mathbb{R}_{\geq
0}\times\mathbb{R}\times\mathbb{R}\mapsto\mathbb{R}_{\geq 0}$ is the Maxwell-
Boltzmann distribution:

${\mathcal{M}}(t,x,v):=\frac{\rho}{\sqrt{2\pi T}}e^{-\frac{(v-u)^{2}}{2T}}.$
(142)

In this expression, $\rho$ is density, $p$ is pressure, and $T=p/\rho$ is
temperature (e.g., see definitions (10)). For $\varepsilon\gg 1$, and for a
fixed $t$ and $x$, the collision operator is weak, and the solution behaves
similarly to Vlasov equation (5). For $\varepsilon\ll 1$, and for a fixed $t$
and $x$, the BGK collision operator forces $f$ towards the Maxwell-Boltzmann
distribution (i.e., thermodynamic equilibrium):

$f(t,x,v)\rightarrow{\mathcal{M}}(t,x,v)+{\mathcal{O}}(\varepsilon).$ (143)

###  7.2 HyQMOM-BGK and the asymptotic-preserving property

Relevant in this work are the first five moments of (141) with the HyQMOM
moment-closure (36):

$\displaystyle{\underline{q}}_{,t}+{\underline{f}}\left({\underline{q}}\right)_{,x}=\frac{1}{\varepsilon}\,{\underline{S^{\text{cons}}}}\left({\underline{q}}\right),\qquad{\underline{\alpha}}_{,t}+{\underline{\underline{B}}}\left({\underline{\alpha}}\right)\,{\underline{\alpha}}_{,x}=\frac{1}{\varepsilon}\,{\underline{S^{\text{prim}}}}\left({\underline{\alpha}}\right),$
(144)

where only the fourth and fifth components of the source terms are nonzero:

$S^{\text{cons}}_{4}=S^{\text{prim}}_{4}=-h,\quad
S^{\text{cons}}_{5}=-k+\frac{2p^{2}}{\rho}-\frac{4puh+h^{2}}{p},\quad
S^{\text{prim}}_{5}=-k+\frac{2p^{2}}{\rho}+\frac{h^{2}}{p}.$ (145)

In eq. 144 we are using definitions (36), (38), and (39). For $\varepsilon\ll
1$, and for a fixed $t$ and $x$, the BGK collision operator forces the heat
flux, $h$, and modified kurtosis, $k$, towards their Maxwell-Boltzmann values:

$h(t,x)=0+{\mathcal{O}}(\varepsilon)\qquad\text{and}\qquad
k(t,x)=\frac{2p(t,x)^{2}}{\rho(t,x)}+{\mathcal{O}}(\varepsilon).$ (146)

In particular, in the $\varepsilon\rightarrow 0^{+}$ limit, solutions of the
HyQMOM-BGK system converge to solutions of the 1D compressible Euler equations
at a convergence rate of ${\mathcal{O}}\left(\varepsilon\right)$:

$\begin{bmatrix}\rho\\\ \rho u\\\ \rho
u^{2}+p\end{bmatrix}_{,t}+\begin{bmatrix}\rho u\\\ \rho u^{2}+p\\\ \rho
u^{3}+3pu\end{bmatrix}_{,x}={\underline{0}},$ (147)

where $h\equiv 0$ and $k\equiv\frac{2p}{\rho}$. Furthermore, by including the
next order term in the Chapman-Enskog expansion, one can show that solutions
to HyQMOM-BGK converge to solutions of the 1D Navier-Stokes equations at a
convergence rate of ${\mathcal{O}}(\varepsilon^{2})$ article:Bardos1991 :

$\begin{bmatrix}\rho\\\ \rho u\\\ \rho
u^{2}+p\end{bmatrix}_{,t}+\begin{bmatrix}\rho u\\\ \rho u^{2}+p\\\ \rho
u^{3}+3pu\end{bmatrix}_{,x}=\begin{bmatrix}0\\\ 0\\\ \frac{3}{2}\varepsilon
pT_{,x}\end{bmatrix}_{,x},$ (148)

where $T=\frac{p}{\rho}$ and $q\equiv-\frac{3}{2}\varepsilon pT_{,x}$.

######  Definition 4 (Asymptotic-preserving (AP) property article:Jin12 )

Let ${\underline{q^{(\Delta t,\Delta x)}}}(t,x;\varepsilon)$ be an
approximation to the exact solution of eq. 144–eq. 145 as computed by a
numerical method with mesh parameters $\Delta t,\Delta x>0$. We assume that
for a fixed $\varepsilon>0$, this method is convergent to the exact solution
of eq. 144–eq. 145. This numerical method is said to be asymptotic-preserving
(AP) provided that the vanishing mesh parameter limit, $\Delta t,\Delta
x\rightarrow 0^{+}$, and the vanishing Knudsen number limit,
$\varepsilon\rightarrow 0^{+}$, commute:

$\lim_{\varepsilon\rightarrow 0^{+}}\left[\lim_{\Delta t,\Delta x\rightarrow
0^{+}}\,{\underline{q^{(\Delta t,\Delta
x)}}}(t,x;\varepsilon)\right]=\lim_{\Delta t,\Delta x\rightarrow
0^{+}}\left[\lim_{\varepsilon\rightarrow 0^{+}}\,{\underline{q^{(\Delta
t,\Delta x)}}}(t,x;\varepsilon)\right].$

Practically, this means that an AP scheme remains stable and accurate for a
fixed mesh, $\Delta t$, and $\Delta x$, for all $\varepsilon>0$, including in
the limit $\varepsilon\rightarrow 0^{+}$.

The goal of this section is to develop an extension of the Lax-Wendroff DG
scheme developed in section 4 and section 5 for the HyQMOM-BGK system eq.
144–eq. 145 that behaves, on the discrete level, as a consistent and stable
numerical method for eq. 147 and eq. 148 in the singular limit
$\varepsilon\rightarrow 0^{+}$. The key innovation in this work is that we
make use of the prediction-correction formulation of Lax-Wendroff DG to
incorporate the collision operator.

###  7.3 Prediction

To describe the HyQMOM-BGK prediction step, it is first useful to define the
following matrices that allow us to map Legendre coefficients to nodal space-
time Gauss-Legendre quadrature points and back again:

$\displaystyle C^{1}_{(a,b)}$
$\displaystyle:=\Psi_{b}\left(\tau_{a},\,\xi_{a}\right),\qquad$
$\displaystyle{\underline{\underline{C^{1}}}}\in\mathbb{R}^{M_{\text{O}}^{2}\times{M_{\text{P}}}},$
(149) $\displaystyle C^{2}_{(b,a)}$
$\displaystyle:=\frac{\omega_{a}}{4}\,\Psi_{b}\left(\tau_{a},\,\xi_{a}\right),\qquad$
$\displaystyle{\underline{\underline{C^{2}}}}\in\mathbb{R}^{{M_{\text{P}}}\times
M_{\text{O}}^{2}},$ (150)

where $\omega_{a}$ and $\left(\tau_{a},\,\xi_{a}\right)$ for
$a=1,\ldots,M_{\text{O}}^{2}$ are tensor product Gauss-Legendre weights and
abscissas. These two matrices satisfy

${\underline{\underline{C^{2}}}}\,{\underline{\underline{C^{1}}}}={\underline{\underline{\,\,{\mathbb{I}}\,\,}}}\in\mathbb{R}^{{M_{\text{P}}}\times{M_{\text{P}}}}.$
(151)

HyQMOM consists of five evolution equations, and the first three are
unaffected by the collision operator; therefore, the update inside the Picard
iteration for the three collision invariants (i.e., density, macroscopic
velocity, and pressure) remains the same as in the collisionless case: eq. 73.
On the other hand, the update for the heat flux, $h$, has a non-zero BGK
contribution; however, the BGK term is linear in the heat flux (see eqs. 144
and 145), which allows for simple treatment. The strategy we pursue here is to
include the BGK source term in the implicit portion of the Picard update to
remain uniformly stable in $\varepsilon>0$. After simple algebra, we arrive at
the following update for the heat flux, $h$:

$\displaystyle\begin{split}\left(\frac{\Delta
t}{2}\,{\underline{\underline{I}}}+\varepsilon\,{\underline{\underline{L}}}\,\right)\,{\underline{W^{n+\frac{1}{2}(j)}_{i\,(:,4)}}}=\,&\frac{\varepsilon}{4}\sum_{a=1}^{{M_{\text{O}}}}\sum_{b=1}^{{M_{\text{O}}}}\,\omega_{a}\,\omega_{b}\,{\underline{{\Psi}}}\left(\mu_{b},\,\mu_{a}\right)\,{\Theta}_{(a,b,4)}\\\
+&\frac{\varepsilon}{4}\sum_{b=1}^{{M_{\text{O}}}}\omega_{b}\,{\underline{{\Psi}}}\left(-1,\xi_{b}\right)\,{\underline{\Phi}}\left(\xi_{b}\right)^{T}{\underline{A^{n}_{i\,(:,4)}}},\end{split}$
(152)

where we are using the short-hand:

${\Theta}_{(a,b,m)}:={\Theta}_{m}\left({\underline{\Psi}}\left(\mu_{b},\mu_{a}\right)^{T}{\underline{\underline{W^{n+\frac{1}{2}(j-1)}_{i}}}}\right).$
(153)

The update for the final primitive variable, $k$ (modified kurtosis), requires
more work. The source term shown in eqs. 144 and 145 is linear in $k$, but it
also includes nonlinear terms from three previously updated quantities:
$\rho$, $p$, and $h$. To construct these nonlinear quantities, we first apply
the mapping from Legendre to nodal values via eq. 149, and then evaluate the
nonlinear portion of the source:

$\displaystyle m=1,3,4:$
$\displaystyle{\underline{\widehat{W}_{(:,m)}}}={\underline{\underline{C^{1}}}}\,{\underline{W_{i(:,m)}^{n+\frac{1}{2}(j)}}}\in\mathbb{R}^{M_{\text{O}}^{2}},$
(154) $\displaystyle a=1,\ldots,M_{\text{O}}^{2}:$
$\displaystyle\left\\{\rho_{a},p_{a},h_{a}\right\\}=\widehat{W}_{(a,\\{1,3,4\\})},\quad\widehat{\mathcal{S}}^{\,(j)}_{ia}=\frac{2p_{a}^{2}}{\rho_{a}}+\frac{h_{a}^{2}}{p_{a}}.$

From here, the update for the modified kurtosis inside the Picard iteration
looks very similar to the update for heat flux (see eq. 152), but with the
additional nonlinear terms computed from eq. 154, which now need to be mapped
back to Legendre coefficients via eq. 150. After some simple algebra, the
update takes the following form:

$\begin{split}\left(\,\frac{\Delta
t}{2}{\underline{\underline{I}}}+\varepsilon\,{\underline{\underline{L}}}\,\right)\,{\underline{W^{n+\frac{1}{2}(j)}_{i\,(:,5)}}}=\,&\frac{\varepsilon}{4}\sum_{a=1}^{{M_{\text{O}}}}\sum_{b=1}^{{M_{\text{O}}}}\,\omega_{a}\,\omega_{b}\,{\underline{{\Psi}}}\left(\mu_{b},\,\mu_{a}\right)\,{\Theta}_{(a,b,5)}\\\
+&\frac{\varepsilon}{4}\sum_{b=1}^{{M_{\text{O}}}}\omega_{b}\,{\underline{{\Psi}}}\left(-1,\xi_{b}\right)\,{\underline{\Phi}}\left(\xi_{b}\right)^{T}{\underline{A^{n}_{i\,(:,5)}}}+\frac{\Delta
t}{2}\,{\underline{\underline{C^{2}}}}\,\,{\underline{\widehat{\mathcal{S}}^{(j)}_{i}}}.\end{split}$
(155)

###  7.4 Post-prediction BGK source evaluation

Once the Picard iterations are complete and all five primitive variables have
been predicted, there is one final computation that must be completed to
prepare us for the correction step: we need to evaluate and project the BGK
source term, ${\underline{S^{\text{cons}}}}$, from eq. 145. This is done
similar to eq. 154 by first mapping the predicted solution from Legendre to
nodal values via eq. 149, then evaluating the source components at nodal
values, and finally mapping back to Legendre coefficients via eq. 150:

$\displaystyle m=1,\ldots,5:$
$\displaystyle{\underline{\widehat{W}_{(:,m)}}}={\underline{\underline{C^{1}}}}\,{\underline{W_{i(:,m)}^{n+\frac{1}{2}}}}\in\mathbb{R}^{M_{\text{O}}^{2}},$
(156) $\displaystyle a=1,\ldots,M_{\text{O}}^{2}:$
$\displaystyle\left\\{\rho_{a},\,u_{a},\,p_{a},\,h_{a},\,k_{a}\right\\}=\widehat{W}_{(a,1:5)},$
$\displaystyle\widehat{\Delta{\mathcal{M}}}_{(a,4)}=-h_{a},\,\,\,\,\,\widehat{\Delta{\mathcal{M}}}_{(a,5)}=-k_{a}+\frac{2p_{a}^{2}}{\rho_{a}}-\frac{4p_{a}u_{a}h_{a}+h_{a}^{2}}{p_{a}},$
$\displaystyle m=4,5:$
$\displaystyle{\underline{\Delta{\mathcal{M}}_{i(:,m)}}}={\underline{\underline{C^{2}}}}\,{\underline{\widehat{\Delta{\mathcal{M}}}_{(:,m)}}}\in\mathbb{R}^{{M_{\text{P}}}}.$

In the above expressions, we use the notation $\Delta{\mathcal{M}}$ to signify
that these BGK source terms are, in fact, measuring the deviations in the heat
flux, $h$, and the modified kurtosis, $k$, from their Maxwell-Boltzmann values
(e.g., see (146)).

We choose to do the above BGK source evaluation and projection, eq. 156, as a
separate step rather than just as part of the correction update since we need
to be extra careful in assuring that the final update is asymptotic-
preserving. Indeed, we show in the next section how to obtain a fully
asymptotic-preserving scheme.

###  7.5 Correction

Just as in the prediction step, we begin by defining matrices that allow us to
map Legendre coefficients to nodal space Gauss-Legendre quadrature points and
back again:

$\displaystyle C^{3}_{(a,b)}$
$\displaystyle:=\Phi_{b}\left(\xi_{a}\right),\qquad$
$\displaystyle{\underline{\underline{C^{3}}}}\in\mathbb{R}^{M_{\text{O}}\times{M_{\text{C}}}},$
(157) $\displaystyle C^{4}_{(b,a)}$
$\displaystyle:=\frac{\omega_{a}}{2}\,\Phi_{b}\left(\xi_{a}\right),\qquad$
$\displaystyle{\underline{\underline{C^{4}}}}\in\mathbb{R}^{{M_{\text{C}}}\times
M_{\text{O}}},$ (158)

where $\omega_{a}$ and $\xi_{a}$ for $a=1,\ldots,M_{\text{O}}$ are Gauss-
Legendre weights and abscissas. These two matrices satisfy

${\underline{\underline{C^{4}}}}\,{\underline{\underline{C^{3}}}}={\underline{\underline{\,\,{\mathbb{I}}\,\,}}}\in\mathbb{R}^{{M_{\text{C}}}\times{M_{\text{C}}}}.$
(159)

As far as the correction step is concerned, the only difference between the
collisionless update, as shown through eqs. 76, 77, 78 and 79 and the BGK
version is the additional BGK source integral needed in eq. 76:

$\begin{split}{\underline{\underline{Q_{i}^{n+1}}}}&={\underline{\underline{Q^{n}_{i}}}}\,+\,\ldots\,+\,\underset{\text{BGK
source}}{\underbrace{\frac{\frac{\Delta
t}{2}}{2\varepsilon}\iint_{-1}^{1}{\underline{\Phi}}\,\left[{\underline{S^{\text{cons}}}}\right]^{T}\,d\tau\,d\xi}}.\end{split}$
(160)

To eventually achieve the asymptotic-preserving (AP) property, we introduce
the following Legendre-in-space-Radau-in-time quadrature:

$\iint_{-1}^{1}g\left(\tau,\xi\right)\,d\tau\,d\xi\approx\sum_{k=1}^{M_{\text{O}}}\sum_{\ell=1}^{M_{\text{O}}}\omega_{k}\,\omega^{R}_{\ell}\,g\left(\tau^{R}_{\ell},\xi_{k}\right),$
(161)

where for $a=1,\ldots,{M_{\text{O}}}$, $(\omega_{a},\xi_{a})$ are again 1D
Gauss-Legendre weights/abscissas, while $(\omega^{R}_{a},\xi^{R}_{a})$ are 1D
Gauss-Radau weights/abscissas. In particular, what we aim to do here is to
handle the BGK source in eq. 160 using a strategy that replaces the actual
Legendre-in-space-Radau-in-time quadrature shown via eq. 161, by a version
where the function values at the $\tau=1$ quadrature points are replaced by
the unknown solution $Q^{n+1}$:

$\displaystyle
m=4,5:\quad\frac{1}{2\varepsilon}\iint_{-1}^{1}{\underline{\Phi}}\,S^{\text{cons}}_{m}\,d\tau\,d\xi\approx\underset{\text{explicit}}{\underbrace{\frac{1}{\varepsilon}\,{\underline{\underline{R}}}\,\,{\underline{\Delta{\mathcal{M}}_{i(:,m)}}}}}+\underset{\text{implicit}}{\underbrace{\frac{1}{r\varepsilon}\left({\underline{{S}_{(:,m)}}}-{\underline{Q^{n+1}_{i(:,m)}}}\right)}},$
(162)

where $\Delta{\mathcal{M}}$ is defined by eq. 156, ${S}$ are Maxwell-Boltzmann
moments (the precise definition is provided below in eq. 164), and

$\displaystyle{\underline{\underline{R}}}=\frac{1}{2}\sum_{k=1}^{M_{\text{O}}}\sum_{\ell=1}^{M_{\text{O}}-1}\omega_{k}\,\omega^{R}_{\ell}\,{\underline{\Phi}}\left(\xi_{k}\right){\underline{\Psi}}\left(\tau^{R}_{\ell},\xi_{k}\right)^{T}\in\mathbb{R}^{{M_{\text{C}}}\times{M_{\text{P}}}},\quad
r=1\bigl{/}\omega^{R}_{M_{\text{O}}}.$ (163)

This quadrature provides a strategy for implicitly handling the BGK collision
term, which is critical for achieving the asymptotic-preserving (AP) property.
We illustrate the modified Gauss-Radau quadrature strategy in fig. 6.

The full correction update is detailed below. The first three moments are
collision invariants and thus updated via eqs. 76, 77, 78 and 79. From these
updated moments, we compute the Maxwell-Boltzmann moments, $S$, that are
required in eq. 162:

$\displaystyle m=1,2,3:$
$\displaystyle{\underline{\widehat{M}_{(:,m)}}}={\underline{\underline{C^{3}}}}\,\,{\underline{Q_{i(:,m)}^{n+1}}}\in\mathbb{R}^{M_{\text{O}}},$
(164) $\displaystyle a=1,\ldots,M_{\text{O}}:$
$\displaystyle\left\\{\rho,\,u,\,p\right\\}=\left\\{\widehat{M}_{(a,1)},\,\frac{\widehat{M}_{(a,2)}}{\widehat{M}_{(a,1)}},\,{\widehat{M}_{(a,3)}}-\frac{\widehat{M}_{(a,2)}^{2}}{\widehat{M}_{(a,1)}}\right\\},$
$\displaystyle\left\\{\widehat{S}_{i(a,4)},\,\widehat{S}_{i(a,5)}\right\\}=\left\\{\rho
u^{3}+3pu,\,\,\rho u^{4}+6pu^{2}+\frac{3p^{2}}{\rho}\right\\}.$

We then update the final two moments in a two-step process, where the first
step is to apply a collisionless update:

$\displaystyle\begin{split}m=4,5:\quad{{\underline{\widetilde{Q}_{i(:,m)}^{n+1}}}}&={\underline{Q_{i(:,m)}^{n}}}-{\frac{\Delta
t}{\Delta
x}}\left(\,{\underline{\Phi}}(1)\,{{\mathcal{F}}^{n+\frac{1}{2}}_{i+\frac{1}{2}(m)}}-{\underline{\Phi}}(-1)\,{{\mathcal{F}}^{n+\frac{1}{2}}_{i-\frac{1}{2}\,(m)}}\,\right)\\\
&+\frac{\Delta t}{2\Delta
x}\sum_{a=1}^{{M_{\text{O}}}}\sum_{b=1}^{{M_{\text{O}}}}\,\omega_{a}\,\omega_{b}\,{\underline{\Phi}}_{,\xi}\left(\mu_{a}\right)f_{m}\left({\underline{\underline{W^{n+\frac{1}{2}}_{i}}}}\,{\underline{\Psi}}\left(\mu_{b},\mu_{a}\right)\right),\end{split}$
(165)

followed by a collision step:

${\underline{Q_{i(:,4:5)}^{n+1}}}=\left(\frac{r\varepsilon}{\frac{\Delta
t}{2}+r\varepsilon}\right){{\underline{\widetilde{Q}_{i(:,4:5)}^{n+1}}}}+\left(\frac{\frac{\Delta
t}{2}}{\frac{\Delta
t}{2}+r\varepsilon}\right)\Biggl{(}r{\underline{\underline{R}}}\,{\underline{\Delta{\mathcal{M}}_{i(:,4:5)}}}+{\underline{\underline{C^{4}}}}\,\,{\underline{\widehat{S}_{i(:,4:5)}}}\Biggr{)}.$
(166)

###### Note 1

All the limiters described in section 5 can still be applied to the HyQMOM-BGK
solver described in this section.

(a) $\widehat{W}^{n+\frac{1}{2}}_{1}$$\widehat{W}^{n+\frac{1}{2}}_{2}$$\widehat{W}^{n+\frac{1}{2}}_{3}$$\widehat{W}^{n+\frac{1}{2}}_{4}$$\widehat{W}^{n+\frac{1}{2}}_{5}$$\widehat{W}^{n+\frac{1}{2}}_{6}$$\widehat{W}^{n+\frac{1}{2}}_{7}$$\widehat{W}^{n+\frac{1}{2}}_{8}$$\widehat{W}^{n+\frac{1}{2}}_{9}$$\widehat{W}^{n+\frac{1}{2}}_{10}$$\widehat{W}^{n+\frac{1}{2}}_{11}$$\widehat{W}^{n+\frac{1}{2}}_{12}$$\widehat{W}^{n+\frac{1}{2}}_{13}$$\widehat{W}^{n+\frac{1}{2}}_{14}$$\widehat{W}^{n+\frac{1}{2}}_{15}$$\widehat{W}^{n+\frac{1}{2}}_{16}$ |  (b) $\widehat{W}^{n+\frac{1}{2}}_{1}$$\widehat{W}^{n+\frac{1}{2}}_{2}$$\widehat{W}^{n+\frac{1}{2}}_{3}$$\widehat{W}^{n+\frac{1}{2}}_{4}$$\widehat{W}^{n+\frac{1}{2}}_{5}$$\widehat{W}^{n+\frac{1}{2}}_{6}$$\widehat{W}^{n+\frac{1}{2}}_{7}$$\widehat{W}^{n+\frac{1}{2}}_{8}$$\widehat{W}^{n+\frac{1}{2}}_{9}$$\widehat{W}^{n+\frac{1}{2}}_{10}$$\widehat{W}^{n+\frac{1}{2}}_{11}$$\widehat{W}^{n+\frac{1}{2}}_{12}$$Q^{n+1}_{1}$$Q^{n+1}_{2}$$Q^{n+1}_{3}$$Q^{n+1}_{4}$
---|---
Figure 6: Numerical quadrature on the canonical space-time element
$(\tau,\xi)\in[-1,1]^{2}$ using a tensor product between 1D Gauss-Legendre
points in $\xi$ and Gauss-Radau points in $\tau$. Panel (a) shows the
Legendre-in-space-Radau-in-time points in the case ${M_{\text{O}}}=4$ and the
known solution values from the prediction step at those quadrature points.
Panel (b) shows the same thing, but the function values at the $\tau=1$
quadrature points are replaced by the unknown solution $Q^{n+1}$. This
strategy of replacing the $\tau=1$ quadrature point function values with the
unknown solution provides a strategy for implicitly handling the BGK collision
term.

###  7.6 Asymptotic-preserving condition

The advantage of the above-proposed scheme for the HyQMOM-BGK is that it
remains high-order accurate uniformly in $\varepsilon>0$ and is asymptotic-
preserving in the $\varepsilon\rightarrow 0^{+}$ limit. The first claim is
demonstrated via numerical examples in the next section; the second claim is
easily demonstrated in this section.

###### Lemma 1

The method LxW-DG method for HyQMOM-BGK described in section 7.3, section 7.4,
and section 7.5 is asymptotic-preserving in the $\varepsilon\rightarrow 0^{+}$
limit.

###### Proof

In the prediction step, the only updates directly affected by the Knudsen
number, $\varepsilon>0$, are the updates for the heat flux, $h$, and the
modified kurtosis, $k$. Taking the $\varepsilon\rightarrow 0^{+}$ limit of
both eq. 152 and eq. 155 yields:

${\underline{W^{n+\frac{1}{2}(j)}_{i\,(:,4)}}}\rightarrow{\underline{0}}\qquad\text{and}\qquad{\underline{W^{n+\frac{1}{2}(j)}_{i\,(:,5)}}}\rightarrow{\underline{\underline{C^{2}}}}\,\,{\underline{\widehat{\mathcal{S}}^{(j)}_{i}}}.$
(167)

This is precisely the desired effect: moments converge to their Maxwell-
Boltzmann values.

In the correction step, the only update directly affected by the Knudsen
number, $\varepsilon>0$, is eq. 166. Taking the $\varepsilon\rightarrow 0^{+}$
limit of this update yields:

${\underline{\Delta{\mathcal{M}}_{i(:,4:5)}}}\rightarrow{\underline{0}}\qquad\text{and}\qquad{\underline{Q_{i(:,4:5)}^{n+1}}}\rightarrow{\underline{\underline{C^{4}}}}\,\,{\underline{\widehat{S}_{i(:,4:5)}}}.$
(168)

Again, this is precisely the desired effect: moments converge to their
Maxwell-Boltzmann values. $\blacksquare$

##  8 HyQMOM-BGK numerical examples

In this section, we apply the proposed HyQMOM-BGK scheme to several test
cases. In §8.1 we verify the claimed orders of accuracy on a smooth
manufactured solution with different Knudsen numbers. These tests also show
the scheme’s uniform accuracy and order of accuracy as a function of the
Knudsen number. In §8.2 we apply the scheme to shock tube initial data with
different Knudsen numbers. These results demonstrate the ability of the non-
oscillatory limiter to control unphysical oscillations. Also shown by these
results is the asymptotic-preserving (AP) property of the scheme for small
$\varepsilon>0$; in particular, we include the exact Riemann solution for the
compressible Euler equations as a point of comparison.

###  8.1 Manufactured solution convergence test

We consider the following manufactured solution:

$\begin{split}{\underline{\alpha^{\text{ms}}}}(t,x)&:=\left(\rho^{\text{ms}},\,u^{\text{ms}},\,p^{\text{ms}},\,h^{\text{ms}},\,k^{\text{ms}}\right)(t,x),\\\
\left(\rho^{\text{ms}},\,p^{\text{ms}},\,h^{\text{ms}},\,k^{\text{ms}}\right)\left(t,x\right)&:=\left(\rho_{\varepsilon},\,p_{\varepsilon},\,h_{\varepsilon},\,k_{\varepsilon}\right)\sqrt{\pi}\left(2-\cos\left(2\pi\left(t-x\right)\right)\right),\\\
u^{\text{ms}}\left(t,x\right)&:=\frac{1-3\varepsilon}{4+8\varepsilon},\end{split}$
(169)

where

$\begin{split}\rho_{\varepsilon}=\frac{1+2\varepsilon}{2+2\varepsilon},\quad
p_{\varepsilon}=\frac{2+33\left(\varepsilon+\varepsilon^{2}\right)}{32(1+\varepsilon)(1+2\varepsilon)},\quad
h_{\varepsilon}=-\frac{125\varepsilon}{128\left(1+2\varepsilon\right)^{2}},\\\
r_{\varepsilon}=\frac{12+\left(\varepsilon+\varepsilon^{2}\right)\left(1021+2017\left(\varepsilon+\varepsilon^{2}\right)\right)}{512(1+\varepsilon)(1+2\varepsilon)^{3}},\quad
k_{\varepsilon}=r_{\varepsilon}-\frac{p^{2}_{\varepsilon}}{\rho_{\varepsilon}}-\frac{h_{\varepsilon}^{2}}{p_{\varepsilon}}.\end{split}$
(170)

Note that this solution is $\varepsilon$-dependent and well-defined for all