Since $f\in\mathfrak{R}$ and $v^{i}\,\ell_{i}\in(-1,1)$, it follows that
$\mathsf{f}(\omega):=\big{(}\,1+v^{i}\,\ell_{i}(\omega)\,\big{)}\,f(\omega)\in\mathfrak{R}$
and thus $\boldsymbol{\mathcal{U}}\in\mathcal{R}$. ∎

######  Lemma 6.

Suppose Assumption 1 holds, $v<1$, and
$\boldsymbol{\mathcal{U}}\in\mathcal{R}$. Then there exists some
$\boldsymbol{\mathcal{M}}\in\mathcal{R}$ that satisfies Eq. (19).

###### Proof.

Let $\mathsf{f}\in\mathfrak{R}$ denote the underlying distribution for
$\boldsymbol{\mathcal{U}}\in\mathcal{R}$. Then the components of
$\boldsymbol{\mathcal{U}}$ can be written as

$\big{(}\mathcal{N},\mathcal{G}_{j}\big{)}^{\intercal}=\frac{1}{4\pi}\int_{\mathbb{S}^{2}}\mathsf{f}(\omega)\,\big{(}1,\ell_{j}(\omega)\big{)}^{\intercal}\,d\omega.$
(107)

Since $\mathsf{f}\in\mathfrak{R}$ and $v^{i}\,\ell_{i}\in(-1,1)$, it follows
that
$f(\omega):=\big{(}\,1+v^{i}\,\ell_{i}(\omega)\,\big{)}^{-1}\,\mathsf{f}(\omega)\in\mathfrak{R}$.
Taking the moments of $f$ leads to $\boldsymbol{\mathcal{M}}\in\mathcal{R}$.
Using the relation between $f$ and $\mathsf{f}$ it is then straightforward to
verify that $\boldsymbol{\mathcal{M}}$ satisfies Eq. (19). ∎

Lemma 6 shows the existence of realizable primitive moments corresponding to
given conserved moments. However, it does not provide guarantees on the
convergence of the iterative solver we use to find the primitive moments. In
the remainder of this subsection, we prove that the iterative solver in Eq.
(57) guarantees the convergence to a realizable moment
$\boldsymbol{\mathcal{M}}$. To start, in the following lemma we show that
realizability is guaranteed at each iteration of the solver in Eq. (57).

######  Lemma 7.

Let $\boldsymbol{\mathcal{U}}\in\mathcal{R}$ and $\lambda\leq\frac{1}{1+v}$ in
Eq. (56). Then, the solver in Eq. (57) guarantees that
$\boldsymbol{\mathcal{M}}^{[k+1]}=(\mathcal{D}^{[k+1]},\boldsymbol{\mathcal{I}}^{[k+1]})^{\intercal}\in\mathcal{R}$,
provided that
$\boldsymbol{\mathcal{M}}^{[k]}=(\mathcal{D}^{[k]},\boldsymbol{\mathcal{I}}^{[k]})^{\intercal}\in\mathcal{R}$.

###### Proof.

We write the iterative update in Eq. (57) as

$\displaystyle\boldsymbol{\mathcal{M}}^{[k+1]}=\left(\begin{array}[]{c}\mathcal{D}^{[k+1]}\\\
\mathcal{I}_{j}^{[k+1]}\end{array}\right)$
$\displaystyle=(1-\lambda)\left(\begin{array}[]{c}\mathcal{D}^{[k]}-\frac{\lambda}{1-\lambda}v^{i}\mathcal{I}_{i}^{[k]}\\\
\mathcal{I}_{j}^{[k]}-\frac{\lambda}{1-\lambda}v^{i}\mathsf{k}_{ij}^{[k]}\mathcal{D}^{[k]}\end{array}\right)+\lambda\left(\begin{array}[]{c}\mathcal{N}\\\
\mathcal{G}_{j}\end{array}\right)$ (108)
$\displaystyle=:(1-\lambda)\,\widetilde{\boldsymbol{\mathcal{M}}}^{[k]}+\lambda\,{\boldsymbol{\mathcal{U}}}.$

Since the realizable set $\mathcal{R}$ is convex and
${\boldsymbol{\mathcal{M}}}^{[k+1]}$ is a convex combination of
$\widetilde{\boldsymbol{\mathcal{M}}}^{[k]}$ and
$\boldsymbol{\mathcal{U}}\in\mathcal{R}$, it suffices to show that
$\widetilde{\boldsymbol{\mathcal{M}}}^{[k]}\in\mathcal{R}$. We observe that
the entries in $\widetilde{\boldsymbol{\mathcal{M}}}^{[k]}$ takes the exact
same form as the ones on the right-hand side of Eq. (19), except with
$\boldsymbol{v}$ replaced by $-\frac{\lambda}{1-\lambda}\boldsymbol{v}$. It
then follows from Lemma 5 that
$\widetilde{\boldsymbol{\mathcal{M}}}^{[k]}\in\mathcal{R}$ if
$\frac{\lambda}{1-\lambda}v\leq 1$, i.e., $\lambda\leq\frac{1}{1+v}$. ∎

It is well-known that, when solving a fixed-point problem defined by a
contraction operator, the Picard iteration converges to the unique fixed point
(see, e.g., [51]). We show below in Proposition 9 that the fixed-point
operator $\boldsymbol{\mathcal{H}}_{\boldsymbol{\mathcal{U}}}$ defined in Eq.
(56) is a contraction under mild assumptions on $v_{h}$, which thus guarantees
the convergence of the iterative solver in Eq. (57). The proof of Proposition
9 uses results from the following two technical lemmas.

######  Lemma 8.

For any $\boldsymbol{\mathcal{M}}\in\mathcal{R}$,
$\|\partial_{\mathcal{D}}(v^{i}\mathsf{k}_{ij}\mathcal{D})\|\leq v$.

###### Proof.

See A.2 for the proof. ∎

######  Lemma 9.

For any $\boldsymbol{\mathcal{M}}\in\mathcal{R}$,
$\|\nabla_{\boldsymbol{\mathcal{I}}}(v^{i}\mathsf{k}_{ij}\mathcal{D})\|\leq
2v$.

###### Proof.

See A.3 for the proof. ∎

We now state and prove Proposition 9.

######  Proposition 9.

Suppose $v<\sqrt{2}-1$ and $\lambda\in(0,1]$. Then,
$\boldsymbol{\mathcal{H}}_{\boldsymbol{\mathcal{U}}}$ defined in Eq. (56) is a
contraction operator, i.e., there exists some $L<1$ such that

$\|\boldsymbol{\mathcal{H}}_{\boldsymbol{\mathcal{U}}}(\boldsymbol{\mathcal{M}}^{(1)})-\boldsymbol{\mathcal{H}}_{\boldsymbol{\mathcal{U}}}(\boldsymbol{\mathcal{M}}^{(2)})\|\leq
L\|\boldsymbol{\mathcal{M}}^{(1)}-\boldsymbol{\mathcal{M}}^{(2)}\|\>,\quad\forall\boldsymbol{\mathcal{M}}^{(1)},\boldsymbol{\mathcal{M}}^{(2)}\in\mathcal{R}\>.$
(109)

###### Proof.

First, for convenience, we denote
$\Delta\mathcal{D}=\mathcal{D}^{(1)}-\mathcal{D}^{(2)}$ and
$\Delta\mathcal{I}_{j}=\mathcal{I}^{(1)}_{j}-\mathcal{I}^{(2)}_{j}$. It then
follows from the definition of
$\boldsymbol{\mathcal{H}}_{\boldsymbol{\mathcal{U}}}$ and the triangle
inequality that

$\|\boldsymbol{\mathcal{H}}_{\boldsymbol{\mathcal{U}}}(\boldsymbol{\mathcal{M}}^{(1)})-\boldsymbol{\mathcal{H}}_{\boldsymbol{\mathcal{U}}}(\boldsymbol{\mathcal{M}}^{(2)})\|\leq(1-\lambda)\left\|\left(\begin{array}[]{l}\Delta\mathcal{D}\\\
\Delta\mathcal{I}_{j}\end{array}\right)\right\|+\lambda\left\|\left(\begin{array}[]{c}v^{i}\,\Delta\mathcal{I}_{i}\\\
v^{i}(\mathsf{k}^{(1)}_{ij}\mathcal{D}^{(1)}-\mathsf{k}^{(2)}_{ij}\mathcal{D}^{(2)})\end{array}\right)\right\|$
(110)

Thus, it suffices to show that, there exists some $\tilde{L}<1$ such that

$\left\|\left(\begin{array}[]{c}v^{i}\,\Delta\mathcal{I}_{i}\\\
v^{i}(\mathsf{k}^{(1)}_{ij}\mathcal{D}^{(1)}-\mathsf{k}^{(2)}_{ij}\mathcal{D}^{(2)})\end{array}\right)\right\|\leq\tilde{L}\left\|\left(\begin{array}[]{l}\Delta\mathcal{D}\\\
\Delta\mathcal{I}_{j}\end{array}\right)\right\|,\quad\forall\boldsymbol{\mathcal{M}}^{(1)},\boldsymbol{\mathcal{M}}^{(2)}\in\mathcal{R}.$
(111)

Lemmas 8 and 9 imply that the gradients of $v^{i}\mathsf{k}_{ij}\mathcal{D}$
in the $\mathcal{D}$ and $\boldsymbol{\mathcal{I}}$ directions are bounded.
Thus, we have

$\displaystyle\|v^{i}(\mathsf{k}^{(1)}_{ij}\mathcal{D}^{(1)}-\mathsf{k}^{(2)}_{ij}\mathcal{D}^{(2)})\|$
$\displaystyle\leq\|\partial_{\mathcal{D}}(v^{i}\mathsf{k}_{ij}\mathcal{D})\|\|\Delta\mathcal{D}\|+\|\nabla_{\boldsymbol{\mathcal{I}}}(v^{i}\mathsf{k}_{ij}\mathcal{D})\|\|\Delta\mathcal{I}_{j}\|$
(112) $\displaystyle\leq
v\|\Delta\mathcal{D}\|+2v\|\Delta\mathcal{I}_{j}\|\>,$

which leads to

$\displaystyle\|v^{i}(\mathsf{k}^{(1)}_{ij}\mathcal{D}^{(1)}-\mathsf{k}^{(2)}_{ij}\mathcal{D}^{(2)})\|^{2}$
$\displaystyle\leq
v^{2}\|\Delta\mathcal{D}\|^{2}+4v^{2}\|\Delta\mathcal{D}\|\|\Delta\mathcal{I}_{j}\|+4v^{2}\|\Delta\mathcal{I}_{j}\|^{2}$
(113)
$\displaystyle\leq(3+2\sqrt{2})v^{2}\|\Delta\mathcal{D}\|^{2}+(2+2\sqrt{2})v^{2}\|\Delta\mathcal{I}_{j}\|^{2}\>,$

where the second inequality follows from the inequality,
$2ab\leq(\sqrt{2}+1)a^{2}+(\sqrt{2}-1)b^{2}$, with
$a=\sqrt{2}v\|\Delta\mathcal{D}\|$ and $b=\sqrt{2}v\|\Delta\mathcal{I}_{j}\|$.
Taking the square of the left-hand side in Eq. (111) and applying the
inequality in Eq. (113) gives

$\displaystyle\left\|\left(\begin{array}[]{c}v^{i}\,\Delta\mathcal{I}_{i}\\\
v^{i}(\mathsf{k}^{(1)}_{ij}\mathcal{D}^{(1)}-\mathsf{k}^{(2)}_{ij}\mathcal{D}^{(2)})\end{array}\right)\right\|^{2}$
$\displaystyle=v^{2}\|\Delta\mathcal{I}_{j}\|^{2}+\|v^{i}(\mathsf{k}^{(1)}_{ij}\mathcal{D}^{(1)}-\mathsf{k}^{(2)}_{ij}\mathcal{D}^{(2)})\|^{2}$
(114)
$\displaystyle\leq(3+2\sqrt{2})v^{2}(\|\Delta\mathcal{D}\|^{2}+\|\Delta\mathcal{I}_{j}\|^{2})\>.$

Let $\tilde{L}:=\sqrt{3+2\sqrt{2}}\,v$, the claim then holds when
$v<\big{(}\sqrt{3+2\sqrt{2}}\,\big{)}^{-1}=\sqrt{2}-1$. ∎

######  Theorem 2.

Suppose $v<\sqrt{2}-1$ and $\lambda\leq\frac{1}{1+v}$ in Eq. (56). Then, for
any given conserved moment
$\boldsymbol{\mathcal{U}}=(\mathcal{N},\boldsymbol{\mathcal{G}})^{\intercal}\in\mathcal{R}$
and initial primitive moment $\boldsymbol{\mathcal{M}}^{[0]}\in\mathcal{R}$,
the iterative solver in Eq. (57) converges to the unique realizable primitive
moment $\boldsymbol{\mathcal{M}}$ that satisfies Eq. (19) as $k\to\infty$.

###### Proof.

This theorem is a direct consequence from the realizability-preserving
property of the solver shown in Lemma 7 and the contraction property of
$\boldsymbol{\mathcal{H}}_{\boldsymbol{\mathcal{U}}}$ proved in Proposition 9.
∎

The results in Theorem 2 lead to the following corollary on the uniqueness of
realizable primitive moments associated with realizable conserved moments.

######  Corollary 1.

Suppose $v<\sqrt{2}-1$. For any conserved moment
$\boldsymbol{\mathcal{U}}\in\mathcal{R}$, there exists a unique realizable
primitive moment $\boldsymbol{\mathcal{M}}$ that satisfies Eq. (19).

###  5.4 Implicit Collision Update

In this section, we prove Proposition 6, which states that, under Assumption
1, the implicit update in Eq. (62c) preserves realizability when the iterative
solver in Eq. (61) is used. We first show in the following lemma that
realizability is preserved in each iteration when starting from a realizable
moment.

######  Lemma 10.

Let
$\boldsymbol{\mathcal{U}}^{{(\ast)}}=(\mathcal{N}^{{(\ast)}},\boldsymbol{\mathcal{G}}^{{(\ast)}})^{\intercal}\in\mathcal{R}$,
$\lambda\leq\frac{1}{1+v}$, and $\kappa\geq\chi\geq 0$ in Eq. (60). Then, the
solver in Eq. (61) guarantees that
$\boldsymbol{\mathcal{M}}^{[k+1]}=(\mathcal{D}^{[k+1]},\boldsymbol{\mathcal{I}}^{[k+1]})^{\intercal}\in\mathcal{R}$,
provided that
$\boldsymbol{\mathcal{M}}^{[k]}=(\mathcal{D}^{[k]},\boldsymbol{\mathcal{I}}^{[k]})^{\intercal}\in\mathcal{R}$.

###### Proof.

We follow the approach in the proof of Lemma 7 and write Eq. (61) as

$\displaystyle\boldsymbol{\mathcal{M}}^{[k+1]}$
$\displaystyle=(1-\lambda)\,\Lambda\left(\begin{array}[]{c}\mathcal{D}^{[k]}-\frac{\lambda}{1-\lambda}v^{i}\mathcal{I}_{i}^{[k]}\\\
\mathcal{I}_{j}^{[k]}-\frac{\lambda}{1-\lambda}v^{i}\mathsf{k}_{ij}^{[k]}\mathcal{D}^{[k]}\end{array}\right)+\lambda\,\Lambda\left(\begin{array}[]{c}\mathcal{N}^{{(\ast)}}+\Delta
t\chi\,\mathcal{D}_{0}\\\ \mathcal{G}_{j}^{{(\ast)}}\end{array}\right)$ (115)
$\displaystyle=:(1-\lambda)\,\Lambda\,\tilde{\boldsymbol{\mathcal{M}}}^{[k]}+\lambda\,\Lambda\,\tilde{\boldsymbol{\mathcal{U}}}^{{(\ast)}}.$

In the proof of Lemma 7, we have shown that
$\tilde{\boldsymbol{\mathcal{M}}}^{[k]}\in\mathcal{R}$ when
$\lambda\leq\frac{1}{1+v}$. Also, it is clear that
$\tilde{\boldsymbol{\mathcal{U}}}^{{(\ast)}}\in\mathcal{R}$ because
${\boldsymbol{\mathcal{U}}}^{{(\ast)}}\in\mathcal{R}$, $\chi\geq 0$ and
$\mathcal{D}_{0}\geq 0$. Since $\kappa\geq\chi\geq 0$, we have
$\mu_{\chi}\geq\mu_{\kappa}\geq 0$, implying that
$\Lambda\,\boldsymbol{\mathcal{M}}\in\mathcal{R}$ for all
$\boldsymbol{\mathcal{M}}\in\mathcal{R}$ based on the definition of
$\mathcal{R}$ in Eq. (22). Therefore,
$\Lambda\,\tilde{\boldsymbol{\mathcal{M}}}^{[k]}$ and
$\Lambda\,\tilde{\boldsymbol{\mathcal{U}}}^{{(\ast)}}$ are both realizable,
which, together with the convexity of $\mathcal{R}$, completes the proof. ∎

Similar to the moment conversion problem considered in Section 5.3, we next
show in the following proposition that the operator $\boldsymbol{\mathcal{Q}}$
in Eq. (60) is a contraction, which implies convergence of the Picard
iteration method in Eq. (61).

######  Proposition 10.

Suppose $v<\sqrt{2}-1$, $\lambda\in(0,1]$, and $\kappa\geq\chi\geq 0$. Then,
$\boldsymbol{\mathcal{Q}}$ is a contraction operator, i.e., there exists some
$L<1$ such that

$\|\boldsymbol{\mathcal{Q}}(\boldsymbol{\mathcal{M}}^{(1)})-\boldsymbol{\mathcal{Q}}(\boldsymbol{\mathcal{M}}^{(2)})\|\leq
L\|\boldsymbol{\mathcal{M}}^{(1)}-\boldsymbol{\mathcal{M}}^{(2)}\|\>,\quad\forall\boldsymbol{\mathcal{M}}^{(1)},\boldsymbol{\mathcal{M}}^{(2)}\in\mathcal{R}\>.$
(116)

###### Proof.

From the definitions of $\boldsymbol{\mathcal{H}}_{\boldsymbol{\mathcal{U}}}$
and $\boldsymbol{\mathcal{Q}}$ in Eqs. (56) and (60), we observe that

$\displaystyle\|\boldsymbol{\mathcal{Q}}(\boldsymbol{\mathcal{M}}^{(1)})-\boldsymbol{\mathcal{Q}}(\boldsymbol{\mathcal{M}}^{(2)})\|$
$\displaystyle=\|\Lambda(\boldsymbol{\mathcal{H}}_{\boldsymbol{\mathcal{U}}}(\boldsymbol{\mathcal{M}}^{(1)})-\boldsymbol{\mathcal{H}}_{\boldsymbol{\mathcal{U}}}(\boldsymbol{\mathcal{M}}^{(2)}))\|$
(117)
$\displaystyle\leq\|\Lambda\|\|\boldsymbol{\mathcal{H}}_{\boldsymbol{\mathcal{U}}}(\boldsymbol{\mathcal{M}}^{(1)})-\boldsymbol{\mathcal{H}}_{\boldsymbol{\mathcal{U}}}(\boldsymbol{\mathcal{M}}^{(2)})\|$

for all
$\boldsymbol{\mathcal{M}}^{(1)},\boldsymbol{\mathcal{M}}^{(2)}\in\mathcal{R}$.
Since $\kappa\geq\chi\geq 0$ and $\lambda>0$, we have
$0\leq\mu_{\kappa}\leq\mu_{\chi}\leq 1$, i.e., $\|\Lambda\|\leq 1$. The claim
is thus a direct consequence of Proposition 9. ∎

######  Theorem 3.

Suppose $v<\sqrt{2}-1$, $\lambda\leq\frac{1}{1+v}$, and $\kappa\geq\chi\geq 0$
in Eq. (60). Then, for any given conserved moment
$\boldsymbol{\mathcal{U}}^{{(\ast)}}=(\mathcal{N}^{{(\ast)}},\boldsymbol{\mathcal{G}}^{{(\ast)}})^{\intercal}\in\mathcal{R}$
and initial primitive moment $\boldsymbol{\mathcal{M}}^{[0]}\in\mathcal{R}$,
the iterative solver in Eq. (61) converges to a unique realizable primitive
moment $\boldsymbol{\mathcal{M}}$ as $k\to\infty$. Further, the conserved
moment $\boldsymbol{\mathcal{U}}$ associated to $\boldsymbol{\mathcal{M}}$ via
Eq. (19) is also realizable and solves the implicit system in Eq. (58).

###### Proof.

The convergence to a unique $\boldsymbol{\mathcal{M}}\in\mathcal{R}$ is given
by the realizability-preserving property in Lemma 10 and the contraction
property in Proposition 10. Realizability of $\boldsymbol{\mathcal{U}}$
follows from Lemma 5, and the formulation of the fixed-point problem in Eq.
(60) guarantees that $\boldsymbol{\mathcal{U}}$ solves the implicit system in
Eq. (58). ∎

###  5.5 Extension to Approximate Moment Closures

In the earlier sections, we have shown the realizability-preserving property
of the numerical scheme (62a)–(62c) under Assumption 1, in which the use of
exact moment closures is assumed. As discussed in Sections 2 and 3, the
approximate Minerbo closure is often used in practice to reduce the
computational cost, where the approximate Eddington factor $\psi_{\mathsf{a}}$
and heat-flux factor $\zeta_{\mathsf{a}}$, defined respectively in Eqs. (27)
and (28), are considered. In this section, we show that the realizability-
preserving and convergence analyses for the conserved-to-primitive moment
conversion (Eq. (19)) and the implicit update (Eq. (62c)) given in Sections
5.3 and 5.4 can be extended to the case when the approximate Minerbo closure
is used.

When the approximate Eddington factor $\psi_{\mathsf{a}}$ in Eq. (27) is used,
we replace Lemma 5 with the following lemma.

######  Lemma 11.

Suppose $v<1$. Let $\boldsymbol{\mathcal{U}}$ be given as in Eq. (19) with
$\boldsymbol{\mathcal{M}}\in\mathcal{R}$, then
$\boldsymbol{\mathcal{U}}\in\mathcal{R}$.

###### Proof.

Since
$\boldsymbol{\mathcal{M}}=\big{(}\mathcal{D},\boldsymbol{\mathcal{I}}\big{)}^{\intercal}\in\mathcal{R}$,
we know that $\mathcal{D}>0$ and $\mathcal{D}-\mathcal{I}\geq 0$. To show
$\boldsymbol{\mathcal{U}}=\big{(}\mathcal{N},\boldsymbol{\mathcal{G}}\big{)}^{\intercal}\in\mathcal{R}$,
we first prove $\mathcal{N}>0$. By definition,

$\mathcal{N}=\mathcal{D}+v^{i}\mathcal{I}_{i}\geq\mathcal{D}-v{\mathcal{I}}>\mathcal{D}-\mathcal{I}\geq
0\>,$ (118)

where the Cauchy-Schwartz inequality and the assumption that $v<1$ are used.
We next prove that $\mathcal{N}^{2}-\mathcal{G}^{2}\geq 0$ with
$\mathcal{G}:=|\boldsymbol{\mathcal{G}}|$, which implies
$\mathcal{N}-\mathcal{G}\geq 0$. Writing $\mathcal{N}^{2}$ and
$\mathcal{G}^{2}$ in terms of the primitive moments leads to

$\displaystyle\mathcal{N}^{2}$
$\displaystyle=\mathcal{D}^{2}+2(v^{i}\mathcal{I}_{i})\mathcal{D}+(v^{i}\mathcal{I}_{i})^{2}\>,$
(119)
$\displaystyle=\mathcal{D}^{2}\big{(}1+2\,v^{i}\,\hat{n}_{i}h+(v^{i}\,\hat{n}_{i})^{2}h^{2}\big{)}\>,$
$\displaystyle\mathcal{G}^{2}$
$\displaystyle=\mathcal{I}^{2}+2\mathcal{I}^{j}(v^{i}\,\mathsf{k}_{ij})\mathcal{D}+(v_{\ell}\,\mathsf{k}^{\ell
j})(v^{i}\,\mathsf{k}_{ij})\mathcal{D}^{2}\>.$
$\displaystyle=\mathcal{D}^{2}\big{(}h^{2}+2\,\hat{n}^{j}\,v^{i}\,\mathsf{k}_{ij}h+v_{\ell}\,\mathsf{k}^{\ell
j}\,v^{i}\,\mathsf{k}_{ij}\,\big{)}\>.$

Using the definition of $\mathsf{k}_{ij}$ in Eq. (10) we obtain

$\displaystyle\hat{n}^{j}\,v^{i}\,\mathsf{k}_{ij}$
$\displaystyle=\hat{n}^{j}\,v^{i}\,\frac{1}{2}\big{(}(1-\psi_{\mathsf{a}})\delta_{ij}+(3\psi_{\mathsf{a}}-1)\hat{n}_{i}\hat{n}_{j}\big{)}=\psi_{\mathsf{a}}\,v^{i}\,\hat{n}_{i}\>,$
(120) $\displaystyle v_{\ell}\,\mathsf{k}^{\ell j}\,v^{i}\,\mathsf{k}_{ij}$
$\displaystyle=\frac{1}{4}\big{(}(1-\psi_{\mathsf{a}})^{2}v^{2}+(1+\psi_{\mathsf{a}})(3\psi_{\mathsf{a}}-1)(v_{i}\,\hat{n}^{i})^{2}\big{)}\>.$
(121)

Plugging these terms into Eq. (119), denoting $s:=v^{i}\hat{n}_{i}$, and using
the assumption that $v<1$ leads to a sufficient condition for
$\mathcal{N}^{2}-\mathcal{G}^{2}\geq 0$: $\forall s\in[-1,1]$ and $\forall
h\in[0,1]$,

$\big{(}1-h^{2}-\frac{1}{4}(1-\psi_{\mathsf{a}})^{2}\big{)}+2(1-\psi_{\mathsf{a}})sh+\big{(}h^{2}-\frac{1}{4}(1+\psi_{\mathsf{a}})(3\psi_{\mathsf{a}}-1)\big{)}s^{2}\geq
0\>.$ (122)

From Lemma 12 (e), we have $1-\psi_{\mathsf{a}}\geq 0$. Thus, by applying the
inequality $2sh\geq-1-s^{2}h^{2}$ to the second term above, it suffices to
show that

$\big{[}\psi_{\mathsf{a}}-h^{2}-\frac{1}{4}(1-\psi_{\mathsf{a}})^{2}\big{]}+\big{[}\psi_{\mathsf{a}}h^{2}-\frac{1}{4}(1+\psi_{\mathsf{a}})(3\psi_{\mathsf{a}}-1)\big{]}s^{2}\geq
0\>.$ (123)

Since $\psi_{\mathsf{a}}-h^{2}-\frac{1}{4}(1-\psi_{\mathsf{a}})^{2}\geq 0$
(Lemma 12 (f)) and $s^{2}\in[0,1]$,

$\big{[}\psi_{\mathsf{a}}-h^{2}-\frac{1}{4}(1-\psi_{\mathsf{a}})^{2}+\psi_{\mathsf{a}}h^{2}-\frac{1}{4}(1+\psi_{\mathsf{a}})(3\psi_{\mathsf{a}}-1)\big{]}s^{2}\geq
0\>,$ (124)

which then becomes

$(\psi_{\mathsf{a}}-h^{2})(1-\psi_{\mathsf{a}})s^{2}\geq 0\>.$ (125)

With $h^{2}\leq\psi_{\mathsf{a}}\leq 1$ from Lemma 12 (e), the proof is
complete. ∎

Lemma 11 extends Lemma 5 by showing that the mapping from primitive moments to
conserved moments via Eq. (19) preserves realizability even when the
approximate Minerbo closure is used. However, there is not an analogous
extension of Lemma 6 to the case of approximate closures. To show that the map
from conserved to primitive moments is realizability-preserving with the
approximate closure, we verify that the analysis from Lemma 7 to Corollary 1
is still valid when the approximate closure is considered. Specifically, when
$\psi$ is replaced by $\psi_{\mathsf{a}}$, the result of Lemma 7 can be
obtained by invoking Lemma 11 rather than Lemma 5 in the proof, the results of
Lemmas 8 and 9 hold since it is shown in A.1 that $\psi_{\mathsf{a}}$ also
satisfies the required properties of $\psi$, and the remainder of the analysis
stays identical to the exact closure case considered in Section 5.3.
Therefore, we have shown that, when $\psi$ is replaced by $\psi_{\mathsf{a}}$,
the iterative solver in Eq. (57) converges to the unique, realizable primitive
moment that satisfies Eq. (19) for the given conserved moment, which implies
that the conserved to primitive moment map from Eq. (19) still preserves
realizability when $v<\sqrt{2}-1$. Further, we also verified that the
convergence and realizability-preserving properties of the iterative solver in
Eq. (61) for the implicit system in Eq. (60) given in Theorem 3 also hold in
the approximate closure case by applying the same arguments to the analysis in
Section 5.4.

##  6 Conservation Property

###  6.1 Simultaneous Number and Energy Conservation of the DG Scheme

It has been shown in Proposition 1 that the two-moment model in Eqs. (2)–(3)
conserves the Eulerian-frame energy up to $\mathcal{O}(v)$. In this section,
we discuss the simultaneous Eulerian-frame number and energy conservation
properties of the two-moment model with the discontinuous Galerkin phase-space
discretization presented in Section 4.1. We are primarily concerned with
consistency with the Eulerian-frame energy equation for the phase-space
advection problem. For this reason, we consider the collisionless case.

Eulerian-frame number conservation follows from the first component of the
semi-discrete DG scheme in Eq. (36) (treating the general case with
$d_{\boldsymbol{x}}=3$),

$\displaystyle\big{(}\,\partial_{t}\mathcal{N}_{h},\varphi_{h}\,\big{)}_{\boldsymbol{K}}=-\sum_{i=1}^{3}\int_{\tilde{\boldsymbol{K}}^{i}}\Big{[}\,\widehat{\mathcal{F}_{\mathcal{N}}^{i}}\big{(}\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}\big{)}\,\varphi_{h}|_{x_{\mbox{\tiny\sc
H}}^{i}}-\widehat{\mathcal{F}_{\mathcal{N}}^{i}}\big{(}\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}\big{)}\,\varphi_{h}|_{x_{\mbox{\tiny\sc
L}}^{i}}\,\Big{]}\,\tau\,d\tilde{\boldsymbol{z}}^{i}+\sum_{i=1}^{3}\big{(}\,\mathcal{F}_{\mathcal{N}}^{i}(\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}),\partial_{i}{\varphi_{h}}\,\big{)}_{\boldsymbol{K}}$
$\displaystyle\hskip
36.0pt-\int_{\boldsymbol{K}_{\boldsymbol{x}}}\Big{[}\,\varepsilon^{3}\,\widehat{\mathcal{F}_{\mathcal{N}}^{\varepsilon}}\big{(}\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}\big{)}\,\varphi_{h}|_{\varepsilon_{\mbox{\tiny\sc
H}}}-\varepsilon^{3}\,\widehat{\mathcal{F}_{\mathcal{N}}^{\varepsilon}}\big{(}\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}\big{)}\,\varphi_{h}|_{\varepsilon_{\mbox{\tiny\sc
L}}}\,\Big{]}\,d\boldsymbol{x}+\big{(}\,\varepsilon\,\mathcal{F}_{\mathcal{N}}^{\varepsilon}(\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}),\partial_{\varepsilon}{\varphi_{h}}\,\big{)}_{\boldsymbol{K}},$
(126)

where $\mathcal{F}_{\mathcal{N}}^{i}$ and
$\mathcal{F}_{\mathcal{N}}^{\varepsilon}$, respectively, are the first
component of the position and energy space fluxes, defined in Eq. (16), and
$\widehat{\mathcal{F}_{\mathcal{N}}^{i}}$ and
$\widehat{\mathcal{F}_{\mathcal{N}}^{\varepsilon}}$ are the corresponding
numerical fluxes, defined in Eqs. (41) and (45). Setting $\varphi_{h}=1$ as
the test function in Eq. (126) results in the equation for the element-
integrated Eulerian-frame number density. Note that the volume terms (the
second and fourth terms) on the right-hand side of Eq. (126) vanish when
$\varphi_{h}=1$. Then, because the numerical fluxes
$\widehat{\mathcal{F}_{\mathcal{N}}^{i}}$ and
$\widehat{\mathcal{F}_{\mathcal{N}}^{\varepsilon}}$ are continuous on element
interfaces, summation over all phase-space elements $\boldsymbol{K}\in D$
results in cancellation of all interior fluxes, and the resulting rate of
change in the total Eulerian-frame particle number is only due to the flow of
particles though the boundary of the domain $D$. That is, the DG scheme for
the Eulerian-frame particle number is conservative by construction.

As for Eulerian-frame energy conservation, similar to Eq. (7) in Proposition
1, the element-integrated Eulerian-frame energy equation can be derived by
adding the Eulerian-frame number equation in Eq. (126), with
$\varphi_{h}=\varepsilon$, and the sum of the three number flux equations in
Eq. (36) with test functions $\varphi_{h}=\varepsilon v_{h}^{j}$. To
accommodate this choice of test functions, the approximation space
$\mathbb{V}_{h}^{k}$ must include the piecewise linear function in the energy
dimension, i.e., $k\geq 1$. Let
$\mathcal{E}_{h}:=\varepsilon\,(\mathcal{N}_{h}+v^{j}_{h}\,{\mathcal{G}}_{h,j})$
denote the discretized Eulerian-frame energy density. Then, the resulting
equation for the element-integrated Eulerian-frame energy takes the form

$\displaystyle\big{(}\,\partial_{t}\mathcal{E}_{h}\,\big{)}_{\boldsymbol{K}}:=\big{(}\,\partial_{t}\mathcal{N}_{h},\varepsilon\,\big{)}_{\boldsymbol{K}}+\big{(}\,\partial_{t}\mathcal{G}_{j,h},\varepsilon
v_{h}^{j}\,\big{)}_{\boldsymbol{K}}$
$\displaystyle=-\sum_{i=1}^{3}\int_{\tilde{\boldsymbol{K}}^{i}}\Big{[}\,\widehat{\mathcal{F}_{\mathcal{E}}^{i}}\big{(}\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}\big{)}|_{x_{\mbox{\tiny\sc
H}}^{i}}-\widehat{\mathcal{F}_{\mathcal{E}}^{i}}\big{(}\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}\big{)}|_{x_{\mbox{\tiny\sc
L}}^{i}}\,\Big{]}\,\tau\,d\tilde{\boldsymbol{z}}^{i}-\int_{\boldsymbol{K}_{\boldsymbol{x}}}\Big{[}\,\varepsilon^{3}\,\widehat{\mathcal{F}_{\mathcal{E}}^{\varepsilon}}\big{(}\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}\big{)}|_{\varepsilon_{\mbox{\tiny\sc
H}}}-\varepsilon^{3}\,\widehat{\mathcal{F}_{\mathcal{E}}^{\varepsilon}}\big{(}\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}\big{)}|_{\varepsilon_{\mbox{\tiny\sc
L}}}\,\Big{]}\,d\boldsymbol{x}$ $\displaystyle\hskip
36.0pt+\big{(}\,\varepsilon\,\mathcal{F}_{\mathcal{N}}^{\varepsilon}(\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h})\,\big{)}_{\boldsymbol{K}}+\sum_{i=1}^{3}\big{(}\,\mathcal{F}_{\mathcal{G}_{j}}^{i}(\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}),\varepsilon\partial_{i}{v_{h}^{j}}\,\big{)}_{\boldsymbol{K}}+\mathcal{O}(v_{h}^{2}),$
(127)

where we have defined the position space numerical fluxes,

$\widehat{\mathcal{F}_{\mathcal{E}}^{i}}\big{(}\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}\big{)}|_{x_{\mbox{\tiny\sc
H}/\mbox{\tiny\sc
L}}^{i}}=\varepsilon\,\big{[}\,\widehat{\mathcal{F}_{\mathcal{N}}^{i}}\big{(}\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}\big{)}+v_{h}^{j}\,\widehat{\mathcal{F}_{\mathcal{G}_{j}}^{i}}\big{(}\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}\big{)}\,\big{]}|_{x_{\mbox{\tiny\sc
H}/\mbox{\tiny\sc L}}^{i}},$ (128)

the energy space numerical fluxes,

$\widehat{\mathcal{F}_{\mathcal{E}}^{\varepsilon}}\big{(}\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}\big{)}|_{\varepsilon_{\mbox{\tiny\sc
H}/\mbox{\tiny\sc
L}}}=\varepsilon\,\widehat{\mathcal{F}_{\mathcal{N}}^{\varepsilon}}\big{(}\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}\big{)}|_{\varepsilon_{\mbox{\tiny\sc
H}/\mbox{\tiny\sc L}}},$ (129)

and $v_{h}^{2}:=|\boldsymbol{v}_{h}|^{2}$. Here,
$\mathcal{F}_{\mathcal{G}_{j}}^{i}$ and
$\widehat{\mathcal{F}_{\mathcal{G}_{j}}^{i}}$ represent the fluxes and the
corresponding numerical fluxes for the number flux equation, defined in Eqs.
(16) and (41), respectively. The third and fourth term on the right-hand side
of Eq. (127), which emanate from the energy derivative term of the number
equation and the spatial derivative of the number flux equations,
respectively, can be written as

$\big{(}\,\varepsilon\,\mathcal{F}_{\mathcal{N}}^{\varepsilon}(\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h})\,\big{)}_{\boldsymbol{K}}+\sum_{i=1}^{3}\big{(}\,\mathcal{F}_{\mathcal{G}_{j}}^{i}(\boldsymbol{\mathcal{U}}_{h},\boldsymbol{v}_{h}),\varepsilon\partial_{i}{v_{h}^{j}}\,\big{)}_{\boldsymbol{K}}=\int_{\boldsymbol{K}}\varepsilon\mathcal{K}^{i}_{\hskip
2.0ptj}\,\big{[}\,\partial_{i}v_{h}^{j}-(\partial_{i}v^{j})_{h}\,\big{]}\,\tau
d\boldsymbol{z}+\mathcal{O}(v_{h}^{2}),$ (130)

where $(\partial_{i}v^{j})_{h}$ is the discretized velocity derivative that
satisfies Eq. (44).

Provided (i) the numerical flux in Eq. (128) is uniquely defined on element
interfaces and (ii) the first term on the right-hand side of Eq. (130)
vanishes, Eq. (127) is, to $\mathcal{O}(v_{h}^{2})$, a phase-space
conservation law for the element-integrated Eulerian-frame energy, in
accordance with Proposition 1. These requirements — which generally require
the discrete velocity $\boldsymbol{v}_{h}$ to be continuous across the
elements — are not satisfied exactly by the DG scheme proposed here. Since the
components of $\boldsymbol{v}_{h}$ are represented by piecewise polynomials,
which are discontinuous on element boundaries, the violation in Eulerian-frame
energy conservation may be larger than what would be expected from
$\mathcal{O}(v_{h}^{2})$ contributions alone. We will investigate the
simultaneous conservation of Eulerian-frame number and energy numerically in
Section 8.

###  6.2 Energy Limiter

In addition to the potential violations of Eulerian-frame energy conservation,
beyond the $\mathcal{O}(v^{2})$ violations inherent to the model, from
discontinuous $\boldsymbol{v}_{h}$ as discussed in Section 6.1, the
realizability-enforcing limiter introduced in Section 5.2 can also affect the
conservation of energy. In fact, for small velocities (including $v=0$, when
the total energy should be preserved to machine precision) the realizability-
enforcing limiter is the dominant source of non-conservation of the Eulerian-
frame energy. To improve Eulerian-frame energy conservation, we propose an
“energy limiter” that corrects the change of energy induced by the
realizability-enforcing limiter via a redistribution of particles between
energy elements through a sweeping procedure. This approach maintains
Eulerian-frame number and energy conservation across all energy elements for a
given spatial element, at the expense of local number conservation in each
energy element. The energy limiter does not correct for Eulerian-frame energy
conservation violations inherent to the $\mathcal{O}(v)$ two-moment model or
due to discontinuous $\boldsymbol{v}_{h}$ (see Section 6.1).

To facilitate the discussion, we denote the element integrated Eulerian-frame
number and energy by N and E, respectively. Given evolved moments
$\boldsymbol{\mathcal{U}}_{h}=(\mathcal{N}_{h},\boldsymbol{\mathcal{G}}_{h})$
on element $\boldsymbol{K}$, the element-integrated number and energy can be
computed by

$\displaystyle\textsf{N}_{\boldsymbol{K}}$
$\displaystyle=\int_{\boldsymbol{K}}{\mathcal{N}}_{h}\varepsilon^{2}d\boldsymbol{z}:=|\boldsymbol{K}|\sum_{\boldsymbol{k}=1}^{|S_{\otimes}^{\boldsymbol{K}}|}w_{\boldsymbol{k}}^{(2)}{\mathcal{N}}_{\boldsymbol{k}},\quad\text{and}$
(131) $\displaystyle\textsf{E}_{\boldsymbol{K}}$
$\displaystyle=\int_{\boldsymbol{K}}({\mathcal{N}}_{h}+v^{j}\mathcal{G}_{h,j})\varepsilon^{3}d\boldsymbol{z}:=|\boldsymbol{K}|\sum_{\boldsymbol{k}=1}^{|S_{\otimes}^{\boldsymbol{K}}|}w_{\boldsymbol{k}}^{(3)}({\mathcal{N}}_{\boldsymbol{k}}+v^{j}\mathcal{G}_{\boldsymbol{k},j}),$
(132)

where $S_{\otimes}^{\boldsymbol{K}}$ denote the set of local DG nodes as
defined in Eq. (48), $\mathcal{N}_{\boldsymbol{k}}$ and
$\boldsymbol{\mathcal{G}}_{\boldsymbol{k}}$ denotes the nodal values at points
in $S_{\otimes}^{\boldsymbol{K}}$, and the weights $w_{\boldsymbol{k}}^{(2)}$
and $w_{\boldsymbol{k}}^{(3)}$ are given by the tensor product of the
$(k+1)$-point one-dimensional LG quadrature rules introduced in Section 4.1,
weighted by $\varepsilon^{2}$ and $\varepsilon^{3}$, respectively. Let
$\widetilde{\boldsymbol{\mathcal{U}}}_{h}:=\texttt{RealizabilityLimiter}(\widehat{\boldsymbol{\mathcal{U}}}_{h})$
be the output of the realizability-enforcing limiter given a potentially non-
realizable solution $\widehat{\boldsymbol{\mathcal{U}}}_{h}$, and let
$(\widetilde{\textsf{N}}_{\boldsymbol{K}},\widetilde{\textsf{E}}_{\boldsymbol{K}})$
and
$(\widehat{\textsf{N}}_{\boldsymbol{K}},\widehat{\textsf{E}}_{\boldsymbol{K}})$
denote the element-integrated number and energy, defined in Eqs. (131) and
(132), for $\widetilde{\boldsymbol{\mathcal{U}}}_{h}$ and
$\widehat{\boldsymbol{\mathcal{U}}}_{h}$, respectively. As discussed in
Section 5.2, the realizability-enforcing limiter gives a solution
$\widetilde{\boldsymbol{\mathcal{U}}}_{h}$ that is realizable on
$\widetilde{S}_{\otimes}^{\boldsymbol{K}}$ while maintaining number
conservation in each element; i.e.,
$\widetilde{\textsf{N}}_{\boldsymbol{K}}=\widehat{\textsf{N}}_{\boldsymbol{K}}$.
However, in part because of the additional factor of $\varepsilon$ in the
definition of the element-integrated energy in Eq. (132), the limiter results
in energy changes (i.e.,
$\widetilde{\textsf{E}}_{\boldsymbol{K}}\neq\widehat{\textsf{E}}_{\boldsymbol{K}}$),
which can lead to
$\sum_{{\boldsymbol{K}}\in\mathcal{T}}\widetilde{\textsf{E}}_{\boldsymbol{K}}\neq\sum_{{\boldsymbol{K}}\in\mathcal{T}}\widehat{\textsf{E}}_{\boldsymbol{K}}$;
i.e., a change in the global Eulerian-frame energy.

The proposed energy limiter corrects Eulerian-frame energy conservation
violations by redistributing particles via a sweeping procedure in the energy
dimension to produce
$\boldsymbol{\mathcal{U}}_{h}:=\texttt{EnergyLimiter}(\widetilde{\boldsymbol{\mathcal{U}}}_{h})$,
as detailed in Algorithm 2. Here we let $\mathcal{T}_{\boldsymbol{x}}$ denote
the collection of all spatial elements $\boldsymbol{K}_{\boldsymbol{x}}$ and
let
$\mathcal{T}_{\varepsilon}=\\{{K}_{\varepsilon,n}\\}_{n=1}^{N_{\varepsilon}}$
denote the collection of all energy elements ${K}_{\varepsilon}$ that cover
the energy domain $D_{\varepsilon}$. For a given spatial element
$\boldsymbol{K}_{\boldsymbol{x}}\in\mathcal{T}_{\boldsymbol{x}}$, the proposed
energy limiter sweeps through elements
$\boldsymbol{K}=K_{\varepsilon}\times\boldsymbol{K}_{\boldsymbol{x}}$ for all
$K_{\varepsilon}\in\mathcal{T}_{\varepsilon}$ in a user-prescribed order to
redistribute particles in a way that the number and energy are both conserved
for the given spatial element
$\boldsymbol{K}_{\boldsymbol{x}}\in\mathcal{T}_{\boldsymbol{x}}$, i.e.,

$\sum_{K_{\varepsilon}\in\mathcal{T}_{\varepsilon}}\textsf{N}_{K_{\varepsilon}\times\boldsymbol{K}_{\boldsymbol{x}}}=\sum_{K_{\varepsilon}\in\mathcal{T}_{\varepsilon}}\widehat{\textsf{N}}_{K_{\varepsilon}\times\boldsymbol{K}_{\boldsymbol{x}}}{\quad\text{and}\quad}\sum_{K_{\varepsilon}\in\mathcal{T}_{\varepsilon}}\widetilde{\textsf{E}}_{K_{\varepsilon}\times\boldsymbol{K}_{\boldsymbol{x}}}=\sum_{K_{\varepsilon}\in\mathcal{T}_{\varepsilon}}\widehat{\textsf{E}}_{K_{\varepsilon}\times\boldsymbol{K}_{\boldsymbol{x}}},$
(133)

which then leads to global number and energy conservation,
$\sum_{{\boldsymbol{K}}\in\mathcal{T}}\textsf{N}_{\boldsymbol{K}}=\sum_{{\boldsymbol{K}}\in\mathcal{T}}\widehat{\textsf{N}}_{\boldsymbol{K}}$
and
$\sum_{{\boldsymbol{K}}\in\mathcal{T}}\textsf{E}_{\boldsymbol{K}}=\sum_{{\boldsymbol{K}}\in\mathcal{T}}\widehat{\textsf{E}}_{\boldsymbol{K}}$,
by summing over all spatial elements.

The sweeping procedure and the particle redistribution strategy used in the
energy limiter are listed in Algorithm 2. Specifically, when the Eulerian-
frame energy conservation violation $\delta\textsf{E}$ is nonzero, the energy
limiter redistributes particles between elements in a pairwise manner to
correct $\delta\textsf{E}$. The pairwise redistribution strategy is detailed
in Algorithm 3, where a pair of scaling coefficients $(\theta_{1},\theta_{2})$
is computed by solving a linear system that requires the sum of the scaled
energies to correct $\delta\textsf{E}$ while preserving the sum of particles
(see Line 3). When at least one of the coefficients ($\theta_{1}$ and
$\theta_{2}$) is less than a prescribed threshold $\theta_{\min}>-1$, a
damping factor $\gamma$ is applied so that $\theta_{1}>\theta_{\min}$ and
$\theta_{2}>\theta_{\min}$, which preserves moment realizability. (The moment
realizability property is invariant to scaling by a positive scalar.) When the
linear system does not have a solution, or when
$\min(\theta_{1},\theta_{2})<\theta_{\min}$, the output of Algorithm 3 does
not fully correct $\delta\textsf{E}$, and the remainder is propagated to the
next pair of elements in the sweeping procedure. As shown in Algorithm 2,
beginning on Line 18, a backward sweep will be launched after the forward
sweep when $|\delta{\textsf{E}}|>\delta$; i.e., when the energy conservation
violation is not fully corrected. Here $\delta>0$ is a user-specified
tolerance on the energy conservation violation. In the implementation, we
choose to omit the condition in Line 19 and perform the full backward sweep in
order to improve the computational efficiency on GPUs. Moreover, to avoid
numerical issues, we restrict the damping factor $\gamma$ in Algorithm 3 such
that the resulting corrected moment ${\boldsymbol{\mathcal{U}}}_{h}$ remains a
strictly positive number density; i.e., $\mathcal{N}_{h}>0$. In the numerical
results reported in Section 8, we choose $\theta_{\min}=-0.5$ and permute the
energy elements in an ascending order based on the associated energy values.
We observe that the forward and backward sweeping procedure is sufficient for
correcting the energy conservation violations introduced by the realizability-
enforcing limiter — i.e., $\delta\textsf{E}\to 0$ during the sweeping
procedure — and that the additional scaling introduced in this energy
correction process has no noticeable adverse impact on the solution to the
two-moment system.

1

2Inputs: Discretized moments before and after the realizability-enforcing
limiter, i.e., $\widehat{\boldsymbol{\mathcal{U}}}_{h}$ and
$\widetilde{\boldsymbol{\mathcal{U}}}_{h}$; a permutation of the energy
elements, denoted as $K_{\varepsilon,n}$, $n=1,\dots,N_{\varepsilon}$.

3 Parameter: $\delta$ (Energy conservation violation tolerance)

${\boldsymbol{\mathcal{U}}}_{h}\leftarrow\widetilde{\boldsymbol{\mathcal{U}}}_{h}$
;

// Initialize the output moment

4 for  _each spatial element
$\boldsymbol{K}_{\boldsymbol{x}}\in\mathcal{T}_{\boldsymbol{x}}$_ do

5 for  _$n=1,\dots,N_{\varepsilon}$_ do

6 Compute
$({\textsf{N}}_{K_{\varepsilon,n}\times\boldsymbol{K}_{\boldsymbol{x}}},{\textsf{E}}_{K_{\varepsilon,n}\times\boldsymbol{K}_{\boldsymbol{x}}})$
from ${\boldsymbol{\mathcal{U}}}_{h}$ using Eqs. (131) and (132);

7 Compute
$\widehat{\textsf{E}}_{K_{\varepsilon,n}\times\boldsymbol{K}_{\boldsymbol{x}}}$
from $\widehat{\boldsymbol{\mathcal{U}}}_{h}$ using Eq. (132);

8

9
${\textsf{N}}_{n}\leftarrow{\textsf{N}}_{K_{\varepsilon,n}\times\boldsymbol{K}_{\boldsymbol{x}}}$,
${\textsf{E}}_{n}\leftarrow{\textsf{E}}_{K_{\varepsilon,n}\times\boldsymbol{K}_{\boldsymbol{x}}}$,
$\widehat{\textsf{E}}_{n}\leftarrow\widehat{\textsf{E}}_{K_{\varepsilon,n}\times\boldsymbol{K}_{\boldsymbol{x}}}$;

10

11 if
_$\sum_{n=1}^{N_{\varepsilon}}{\textsf{E}}_{n}\neq\sum_{n=1}^{N_{\varepsilon}}\widehat{\textsf{E}}_{n}$_
then

12 $\delta{\textsf{E}}\leftarrow{\textsf{E}}_{1}-\widehat{\textsf{E}}_{1}$;

/* Forward sweep */

13 for  _$n=1,\dots,N_{\varepsilon}-1$_ do

14
$\delta{\textsf{E}}\leftarrow\delta{\textsf{E}}+{\textsf{E}}_{n+1}-\widehat{\textsf{E}}_{n+1}$;

15 if  _$|\delta{\textsf{E}}| >\delta$_ then

16
$(\theta_{n},\theta_{n+1})=\texttt{ComputeCorrection}(\textsf{N}_{n},\textsf{E}_{n},\textsf{N}_{n+1},\textsf{E}_{n+1},\delta\textsf{E})$;

/* Update corrected moments, numbers, and energies */

17
${\boldsymbol{\mathcal{U}}}_{h}\leftarrow(1+\theta_{n})\,{\boldsymbol{\mathcal{U}}}_{h}$
on $K_{\varepsilon,n}\times\boldsymbol{K}_{\boldsymbol{x}}$, and
${\boldsymbol{\mathcal{U}}}_{h}\leftarrow(1+\theta_{n+1})\,{\boldsymbol{\mathcal{U}}}_{h}$
on $K_{\varepsilon,n+1}\times\boldsymbol{K}_{\boldsymbol{x}}$;

18
$({\textsf{N}}_{n},{\textsf{E}}_{n})\leftarrow(1+\theta_{n})\,({\textsf{N}}_{n},{\textsf{E}}_{n})$,
and
$({\textsf{N}}_{n+1},{\textsf{E}}_{n+1})\leftarrow(1+\theta_{n+1})\,({\textsf{N}}_{n+1},{\textsf{E}}_{n+1})$;

19
$\delta{\textsf{E}}\leftarrow{\textsf{E}}_{n}+{\textsf{E}}_{n+1}+\delta{\textsf{E}}$;

20

21

/* Backward sweep */

22 for  _$n=N_{\varepsilon}-1,\dots,2$_ do

23 if  _$|\delta{\textsf{E}}| >\delta$_ then

24
$(\theta_{n},\theta_{n-1})=\texttt{ComputeCorrection}(\textsf{N}_{n},\textsf{E}_{n},\textsf{N}_{n-1},\textsf{E}_{n-1},\delta\textsf{E})$;

/* Update corrected moments, numbers, and energies */

25
${\boldsymbol{\mathcal{U}}}_{h}\leftarrow(1+\theta_{n})\,{\boldsymbol{\mathcal{U}}}_{h}$
on $K_{\varepsilon,n}\times\boldsymbol{K}_{\boldsymbol{x}}$, and
${\boldsymbol{\mathcal{U}}}_{h}\leftarrow(1+\theta_{n-1})\,{\boldsymbol{\mathcal{U}}}_{h}$
on $K_{\varepsilon,n-1}\times\boldsymbol{K}_{\boldsymbol{x}}$;

26
$({\textsf{N}}_{n},{\textsf{E}}_{n})\leftarrow(1+\theta_{n})\,({\textsf{N}}_{n},{\textsf{E}}_{n})$,
and
$({\textsf{N}}_{n-1},{\textsf{E}}_{n-1})\leftarrow(1+\theta_{n-1})\,({\textsf{N}}_{n-1},{\textsf{E}}_{n-1})$;

27
$\delta{\textsf{E}}\leftarrow{\textsf{E}}_{n}+{\textsf{E}}_{n-1}+\delta{\textsf{E}}$;

28

29 else

30 break;

31

32

33

34

Algorithm 2
${\boldsymbol{\mathcal{U}}}_{h}=\texttt{EnergyLimiter}(\widetilde{\boldsymbol{\mathcal{U}}}_{h})$

1

2Inputs:
$\textsf{N}_{1},\textsf{E}_{1},\textsf{N}_{2},\textsf{E}_{2},\delta\textsf{E}$

3 Parameter: $\theta_{\min}>-1$

4

5Compute $(\theta_{1},\theta_{2})$ by solving $\left\\{\begin{aligned}
\theta_{1}{\textsf{N}}_{1}+\theta_{2}{\textsf{N}}_{2}&=0\>\\\
\theta_{1}{\textsf{E}}_{1}+\theta_{2}{\textsf{E}}_{2}&=-\delta\textsf{E}\>\end{aligned}\right.$;

6 if  _no solution_ then

7 $(\theta_{1},\theta_{2})\leftarrow(0,0)$;

8

9if  _$\min(\theta_{1},\theta_{2}) <\theta_{\min}$_ then

10 $\gamma\leftarrow\frac{\theta_{\min}}{\min(\theta_{1},\theta_{2})}$;

$(\theta_{1},\theta_{2})\leftarrow(\gamma\theta_{1},\gamma\theta_{2})$ ;

// Limit for realizability

11

12

Algorithm 3
$(\theta_{1},\theta_{2})=\texttt{ComputeCorrection}(\textsf{N}_{1},\textsf{E}_{1},\textsf{N}_{2},\textsf{E}_{2},\delta\textsf{E})$

##  7 Implementation, Programming Models, and Portability

The DG-IMEX method proposed here has been implemented in the toolkit for high-
order neutrino radiation hydrodynamics (thornado). Here we briefly discuss
some considerations in this process.

Neutrino transport is only one component (along with, e.g., hydrodynamics,
nuclear reaction kinetics, and gravity) of a broader, multiphysics simulation
framework needed to model multiscale astrophysical systems, e.g., core-
collapse supernova explosions. However, the number of evolved degrees of
freedom is relatively high compared to other components. For example,
simulations incorporating a two-moment model (four moments), evolving three
independent neutrino flavors (six species), with 16 linear elements (k=1) to
discretize the energy dimension evolve $4\times 6\times 16\times(k+1)=768$
degrees of freedom per spatial point. As such, spectral neutrino radiation
transport represents the bulk of the computational load in such scientific
applications. With this in mind, node-level performance and portability for
heterogeneous computing systems are prioritized in thornado development as a
collection of modular physics components that can be incorporated into
distributed mutliphysics simulation codes (e.g., Flash-X [54]), which are
equipped with native infrastructure for distributed parallelism. In
particular, we target frameworks that utilize adaptive mesh refinement, where
simulation data is mapped to smaller grid blocks of relatively even size.

thornado uses a combination of compiler directives and optimized linear
algebra libraries to accelerate all components of the DG-IMEX method. All of
the solver components — e.g., the computation of numerical fluxes, evaluation
of phase-space divergences, and limiters — are reduced to small, discrete
kernels that can be executed either as collapsible, tightly-nested loops over
phase-space dimensions or basic linear algebra operations. In addition to
optimizing many key metrics for GPU performance (e.g., occupancy, register
pressure, and memory coalescence), this strategy naturally exposes vector-
level parallelism which also benefits performance on modern, multicore CPUs.
This is especially important when invoking iterative solvers, such as those
described in Sections 4.3.1 and 4.3.2, across many independent phase-space
points. Since iteration counts can vary, assigning an even number of phase-
space points to each thread can lead to severe load imbalance among GPU
threads. We address this problem by tracking the convergence of each point
independently, removing them from calculations in each kernel until all points
have converged.

Our portability strategy focuses on maintaining a single code-base that can
efficiently execute on different hardware architectures and software
environments. thornado contains three distinct implementations of compiler
directives that are managed with C preprocessor macros: traditional OpenMP
(CPU multi-core), OpenMP offload (GPU), and OpenACC (GPU). We refer to code
listings in [55] for specific examples.

Interfaces to optimized linear algebra routines are also written in a generic
way for portability across different libraries. Currently, thornado has linear
algebra interfaces supporting several LAPACK and BLAS [56] routines with GPU
implementations from NVIDIA, AMD, Intel, and MAGMA [57]. This approach hides
the complexities of managing different interfaces in a single thornado module
that can be easily used throughout the code. In addition to the individual
routine interfaces, each linear algebra package requires specific attention to
interoperability with the compiler directives to ensure correct
synchronization when using multiple execution streams per device. This is
managed during initialization with compiler directives and C preprocessor
macros.

We provide timing results and a breakdown of the computational cost associated
with key solver components for one of the numerical examples in Section 8.

##  8 Numerical Tests

In this section, we demonstrate the performance of our implementation of the
DG-IMEX method to solve the $\mathcal{O}(v)$ two-moment model. We consider
problems with and without collisions. For problems with collisions, we use the
IMEX scheme proposed in [34] (see also [35] for details). For problems without
collisions, we use the optimal second- and third-order accurate strong
stability-preserving Runge–Kutta methods from [58], referred to as SSPRK2 and
SSPRK3, respectively. For the tests in Sections 8.2 and 8.3, unless stated
otherwise, we set the time step to $\Delta
t=0.3\times|K_{\boldsymbol{x}}^{1}|/(k+1)$, where $k$ is the polynomial
degree. For the tests in Sections 8.4-8.6, we enforce the time step
restriction given in Theorem 1.

Collisions tend to drive the distribution towards isotropy in the angular
dimensions of momentum space (i.e., $|\boldsymbol{\mathcal{I}}|\to 0$), which
places the comoving-frame moments $\boldsymbol{\mathcal{M}}$ safely inside the
realizable domain. Therefore, to emphasize the improved robustness resulting
from our analysis, our main focus is on phase-space advection problems without
collisions, where the moments evolve close to the boundary of the realizable
domain.

###  8.1 Moment Conversion Solver

The solution of the conserved-to-primitive moment conversion problem in Eq.
(19) and the implicit system in Eq. (58) contribute the majority of the
computational cost of the realizability-preserving scheme. In this section, we
test the iterative solver for solving the moment conversion problem Eq. (19)
with various solver configurations, and the results reported provide guidance
for selecting iterative solver configurations for this critical part of the
algorithm.

As discussed in Section 4.3.1, we formulate the moment conversion problem in
Eq. (19) as a fixed-point problem on the primitive moments
$\boldsymbol{\mathcal{M}}=(\mathcal{D},\boldsymbol{\mathcal{I}})^{\intercal}$
of the form stated in Eq. (56). In Lemma 7, we have shown that the moment
realizability is preserved in the iterative procedure when Eq. (56) is solved
with the Picard iteration method in Eq. (57) and $\lambda\leq(1+v)^{-1}$ in
Eq. (56). The convergence of Picard iteration is guaranteed in Theorem 2 with
the additional assumption that $v<\sqrt{2}-1$.

We first compare the iteration counts required for convergence of the Picard
iteration solver and an Anderson acceleration (AA) solver, using two different
choices for $\lambda$. The AA technique was first proposed in [59] to
accelerate the convergence of fixed-point iterations by accounting for the
past iteration history to compute new iterates. Here we follow the formulation
and implementation in [60, 35] and apply the AA solver to the moment
conversion problem Eq. (56). In Figure 4, the iteration counts are reported
for the two iterative solvers applied to solve Eq. (56) at varying fluid speed
$v:=|\boldsymbol{v}|$ and flux factor
$h=|\boldsymbol{\mathcal{I}}|/\mathcal{D}$, with $\lambda$ chosen to be the
largest allowable value, i.e., $\lambda=(1+v)^{-1}$ and a more conservative
value $\lambda=0.5$. The AA solver uses the memory parameter $m=1$ (defined in
[35]), so that only information from the previous and current iterate is used.
The stopping criteria for both solvers are given as

$\|\boldsymbol{\mathcal{M}}^{[k]}-\boldsymbol{\mathcal{M}}^{[k-1]}\|\leq\texttt{tol}\,\|\boldsymbol{\mathcal{U}}\|,$
(134)

where we consider the norms in the $L^{2}$ sense and the tolerance
$\texttt{tol}=10^{-8}$. For each choice of $(v,h)$ in Figure 4, the fixed-
point problem is solved for 100 randomly generated
$\boldsymbol{\mathcal{U}}\in\mathcal{R}$ (varying the direction of
$\boldsymbol{v}$ and $\boldsymbol{\mathcal{I}}/\mathcal{D}$ randomly), and the
averaged iteration counts over these 100 problems are recorded. In each test,
the initial guess takes the form
$\boldsymbol{\mathcal{M}}^{[0]}=\boldsymbol{\mathcal{U}}$. The results in
Figure 4 illustrate that, for both the Picard iteration and the AA solvers,
choosing the parameter $\lambda$ to be the largest allowable value
$(1+v)^{-1}$ indeed reduces the iteration counts from the more conservative
choice $\lambda=0.5$, particularly in the low velocity regime. In addition, it
can be found from Figure 4 that AA solver consistently outperforms the Picard
iteration method, and the advantage of using AA grows as the velocity
increases. We note that the realizability-preserving and convergence
properties analyzed in Section 5.3 are only applicable to the Picard iteration
solver, and not to the AA solver.777The realizability-preserving and
convergence properties of the AA solver require additional conditions such as
boundedness of extrapolation coefficients, which we do not enforce in the
implementation. However, in the numerical results reported in Figure 4, we
have not observed convergence failure by any of the solvers, even when the
velocity is larger than the upper bound ($v=\sqrt{2}-1$; plotted as a red
vertical line in each panel in Figure 4) required in the convergence analysis
in Theorem 2.

In Figure 5, we show results from experimenting with two choices for the
initial guess,
$\boldsymbol{\mathcal{M}}^{[0]}=(\mathcal{N},\boldsymbol{0})^{\intercal}$ and
$\boldsymbol{\mathcal{M}}^{[0]}=\boldsymbol{\mathcal{U}}=(\mathcal{N},\boldsymbol{\mathcal{G}})^{\intercal}$,
for the AA solver with $\lambda=(1+v)^{-1}$, which is the best performing
configuration shown in Figure 4. As shown in Figure 5, initializing with the
conserved moment $\boldsymbol{\mathcal{U}}$ generally outperforms the
isotropic initial condition $(\mathcal{N},\boldsymbol{0})$, except for the
case when the flux factor $h=0$, for which the isotropic initial condition is
exactly the primitive moment. Since we expect moments with $h=0$ to be rarely
encountered in numerical simulations, adopting the AA solver with
$\lambda=(1+v)^{-1}$ and initial guess
$\boldsymbol{\mathcal{M}}^{[0]}=\boldsymbol{\mathcal{U}}$ appears to be the
best choice. This conjecture is confirmed in the performance comparison
reported in Section 8.7, where we observe a considerable improvement in terms
of computational time by using the initial guess
$\boldsymbol{\mathcal{M}}^{[0]}=\boldsymbol{\mathcal{U}}$.

|
---|---
|
Figure 4: Iteration counts for the Picard iteration (top panels) and AA (bottom panels) solvers with modified Richardson iteration parameter $\lambda=(1+v)^{-1}$ and $\lambda=0.5$ (right and left columns, respectively), applied to the moment conversion problems with various velocity $v$ and flux factor $h$. The reported iteration counts are the average over 100 randomly generated moment conversion problems at each $(v,h)$, where the randomness is applied to the directions of $\boldsymbol{v}$ and $\boldsymbol{\mathcal{I}}/\mathcal{D}$. In each panel, the red vertical line indicates the upper velocity bound for guaranteed convergence, $v=\sqrt{2}-1$. |
---|---
Figure 5: Iteration counts for the AA solver with modified Richardson
iteration parameters $\lambda=(1+v)^{-1}$ applied to the moment conversion
problems with various velocity $v$ and flux factor $h$ for two different
initial guesses. The reported iteration counts are the average over 100
randomly generated moment conversion problems at each $(v,h)$.

###  8.2 Sine Wave Streaming

The first test we consider that evolves the two-moment system models free-
streaming radiation through a background with a spatially (and temporally)
constant velocity field in the $x^{1}$-direction. That is, we set
$\chi=\sigma=0$, while $\boldsymbol{v}=[v,0,0]^{\intercal}$, with $v=0.1$. The
purpose of this test is to verify (i) the correct radiation propagation speed
in this idealized setting, and (ii) the expected order of accuracy of the
implemented method. We consider a periodic one-dimensional unit spatial domain
$D_{x^{1}}=[0,1]$. The initial number density and flux are set to
$\mathcal{D}(x^{1},0)=\mathcal{D}_{0}(x^{1})=0.5+0.49\times\sin(2\pi x^{1})$
and $\mathcal{I}^{1}(x^{1},0)=\mathcal{D}_{0}(x^{1})$, respectively. Then, the
flux factor is $h=1$, and the analytic solution is given by
$\mathcal{D}(x^{1},t)=\mathcal{I}^{1}(x^{1},t)=\mathcal{D}_{0}(x^{1}-t)$;
i.e., the initial profile propagates with unit speed, _independent_ of $v$.
(As noted by [23], dropping the velocity-dependent terms in the time
derivatives of Eqs. (2) and (3), as is done in [3, 4], the propagation speed
becomes $1+v$ for this test, which is unphysical.) Since the background
velocity is constant, there is no coupling in the energy dimension. Therefore,
this test is performed with a single energy. We run this test until $t=1$,
when the initial profile has crossed the grid once before returning to its
initial position.

Figure 6: Error in the $L^{2}$ norm versus number of spatial elements $N$ for
the Sine Wave Streaming test. Results obtained with second- and third-order
schemes (using $k=1$ polynomials and SSPRK2 time stepping and $k=2$
polynomials and SSPRK3 time stepping, respectively), along with dotted
reference lines proportional to $1/N^{k+1}$, are plotted using black and red
lines, respectively.

In Figure 6, we plot the error in the $L^{2}$ norm versus the number of
spatial elements for the second-order scheme, using linear polynomials ($k=1$)
and SSPRK2 time stepping, and the third-order scheme, using quadratic
polynomials ($k=2$) and SSPRK3 time stepping. The results confirm the expected
convergence rate to the exact solution.

###  8.3 Gaussian Diffusion

The next test we consider, adopted from [3], models the diffusion of particles
through a medium moving with constant velocity in the $x^{1}$-direction. We
consider a purely scattering medium and set $\sigma=3.2\times 10^{3}$ and
$\chi=0$, and let $\boldsymbol{v}=[v,0,0]^{\intercal}$, with $v=0.1$. We let
the spatial domain be periodic, $D_{x^{1}}=[0,3]$, with initial conditions
$\mathcal{D}_{0}(x^{1})=\exp[-(x^{1}-x_{0}^{1})^{2}/(4t_{0}\kappa_{\mathsf{D}})]$
and
$\mathcal{I}_{0}^{1}(x^{1})=-\kappa_{\mathsf{D}}\partial_{x^{1}}\mathcal{D}_{0}$,
where $\kappa_{\mathsf{D}}=(3\sigma)^{-1}$, and we set $x_{0}^{1}=1$ and
$t_{0}=5$. Then, the evolution of the number density is approximately governed
by the advection-diffusion equation

$\partial_{t}\mathcal{D}+\partial_{x^{1}}\big{(}\,\mathcal{D}v-\kappa_{\mathsf{D}}\partial_{x^{1}}\mathcal{D}\,\big{)}=0,$
(135)

whose analytical solution is given by [3]

$\mathcal{D}(x^{1},t)=\sqrt{\frac{t_{0}}{t_{0}+t}}\exp\Big{\\{}\,-\frac{((x^{1}-vt)-x_{0}^{1})^{2}}{4(t_{0}+t)\kappa_{\mathsf{D}}}\,\Big{\\}}.$
(136)

Since there is no coupling in the energy dimension ($v$ is constant), we
perform this test with a single energy. We use quadratic elements ($k=2$) and
the IMEX time stepping scheme from [34], integrating the collision term
implicitly. For this test, the time step is set to $\Delta t=C_{\rm
CFL}\times|K_{\boldsymbol{x}}^{1}|$, where $C_{\rm CFL}$ is specified below.
The purpose of this test is to investigate the performance of the DG-IMEX
scheme in a regime where both advection and diffusion contribute to the
evolution of the number density. For $t>0$, the Gaussian profile is advected
with the flow, while the amplitude decreases and the width increases due to
diffusion.

The left panel in Figure 7 shows the number density versus $x^{1}-vt$ for
various times as the Gaussian profile propagates once across the periodic
domain and returns to its initial position at $t=30$. At this time the
amplitude is reduced by a factor $\sqrt{5/35}\approx 0.378$. For this
simulation, the spatial domain is discretized using $96$ elements, so that the
radio of the element width to the mean-free path is
$|K_{\boldsymbol{x}}^{1}|\sigma=10^{2}$. The numerical solution (open circles)
agrees well with the expression in Eq. (136) (solid lines).

The right panel in Figure 7 shows the error in the $L^{2}$ norm at $t=5$
versus the number of spatial elements $N$. Since the expression given by Eq.
(136) is not an exact solution to the $\mathcal{O}(v)$ two-moment model in the
limit of high scattering opacity, we compare the numerical results to a high-
resolution reference solution, computed with $8192$ elements. (We have
confirmed that for fixed $N$ and $t$, and varying $v$, the difference between
the numerical solution and the expression in Eq. (136) increases as
$\mathcal{O}(v^{2})$.) The solid black curve with squares shows the error
obtained with the standard CFL number $C_{\rm CFL}:=C_{\rm
CFL}^{0}=0.3/(k+1)$. For smaller $N$, the error falls off as $N^{-3}$ (see
dashed black reference line), consistent with a third-order convergence rate,
while for larger $N$ the convergence rate transitions to first-order (see
dotted black reference line). The reason for this change in convergence rate
is because the IMEX scheme, taken from [34], is formally only first-order
accurate. Spatial discretization errors dominate for small $N$, but since
these errors decrease with the third-order rate, temporal errors become
dominant for large $N$. To verify this, we also plot convergence results
obtained after reducing the time step by a factor of $25$, $C_{\rm
CFL}:=C_{\rm CFL}^{0}/25$; see solid red curve with squares. For this case,
the error decreases with the third-order rate for all $N$. We also show the
error obtained with a second-order IMEX scheme (IMEXRKCB2 from [61]), using
the standard CFL number. Due to better temporal accuracy, the error decreases
with the third-order rate, but the scheme does not satisfy the convex-
invariant conditions delineated in [34], and can therefore not be guaranteed
to maintain moment realizability by our analysis.

|
---|---
Figure 7: Results for the Gaussian diffusion test. The left panel shows the
numerical solution (open circles) versus the shifted coordinate $x^{1}-vt$ for
various times, compared with the analytic solution (solid lines) to the
advection-diffusion equation in Eq. (135). The right panel shows the error in
the $L^{2}$ norm versus the number of elements $N$. The error is computed with
respect to a high-resolution reference run with $N=8192$. Convergence results
are shown for standard and reduced CFL number, using the first-order SSP-IMEX
scheme from [34] (solid black and red, respectively; see text for details),
and a second-order IMEX scheme from [61] with standard CFL number (dotted
blue). The dotted and dashed black reference lines are proportional to
$N^{-1}$ and $N^{-3}$, respectively.

###  8.4 Streaming Doppler Shift

This test, adopted from [24] (see also [3, 4]), models the propagation of
free-streaming radiation along the $x^{1}$-direction through a background with
a spatially varying velocity field. Because the two-moment model adopts
momentum-space coordinates associated with a comoving observer, the radiation
energy spectra will be Doppler shifted. We consider a one-dimensional spatial
domain $D_{x^{1}}=[0,10]$. Again, we set $\chi=\sigma=0$, while the velocity
field is set to $\boldsymbol{v}=(v,0,0)^{\intercal}$, where

$v(x^{1})=\left\\{\begin{array}[]{ll}0,&x^{1}\in[0,2)\\\
v_{\max}\times\sin^{2}[2\pi(x^{1}-2)/6],&x^{1}\in[2,3.5)\\\
v_{\max},&x^{1}\in[3.5,6.5)\\\
v_{\max}\times\sin^{2}[2\pi(x^{1}-2)/6],&x^{1}\in[6.5,8)\\\
0,&x^{1}\in[8,10]\end{array}\right.,$ (137)

and where we will vary $v_{\max}$. We set the energy domain to
$D_{\varepsilon}=[0,50]$. In this test, we discretize the spatial and energy
domains into $128$ and $32$ elements, respectively, and use quadratic elements
($k=2$) and SSPRK3 time stepping. In the computational domain, the moments are
initially set to $\mathcal{D}=1\times 10^{-40}$ and $\mathcal{I}^{1}=0$ for
all $(x^{1},\varepsilon)\in D_{x^{1}}\times D_{\varepsilon}$. At the inner
spatial boundary, we impose an incoming, forward-peaked radiation field with a
Fermi-Dirac spectrum; i.e., we set
$\mathcal{D}(\varepsilon,x^{1}=0)=1/[\exp(\varepsilon/3-3)+1]$ and
$\mathcal{I}^{1}(\varepsilon,x^{1}=0)=0.999\times\mathcal{D}(\varepsilon,x^{1}=0)$,
so that the flux factor $h\approx 1$. (We impose outflow boundary conditions
at $x^{1}=10$.) Then, for $t>0$, a radiation front propagates through the
computational domain, and a steady state is established for $t\gtrsim 10$,
where the spectrum is Doppler-shifted according to the velocity field. From
special relativistic considerations, similar to [3], the analytical spectral
number density in the steady state can be written as

$\mathcal{D}_{\rm A}=\frac{s^{2}}{\exp(s\varepsilon/3-3)+1},$ (138)

where $s=\sqrt{(1+v)/(1-v)}$. The purpose of this test is to (i) compare
steady state numerical solutions with the prediction from special relativity
given by Eq. (138), and (ii) investigate the simultaneous Eulerian-frame
number and energy conservation properties of the method as the initial
conditions are evolved to steady state. To reach an approximate steady state,
we run all models until $t=20$.

####  8.4.1 General Solution Characteristics

|
---|---
|
Figure 8: Steady state solutions ($t=20$) for the streaming Doppler shift
problem for various $v_{\max}\in\\{0.0,0.1,0.2,0.4\\}$. In the top panels, we
plot spectra at $x^{1}=5$ ($\mathcal{D}\varepsilon^{2}$ versus $\varepsilon$;
left) and the RMS energy, as defined in Eq. (139), versus position $x^{1}$
(right). In these panels, solid lines represent the analytic (A) solution from
special relativistic considerations, given by Eq. (138), while dotted lines
represent the numerical (N) results. In all panels, black, red, blue, and
magenta curves represent runs with $v_{\max}$ set to $0.0$, $0.1$, $0.2$, and
$0.4$, respectively. In the bottom left panel, the Eulerian-frame (solid
lines) and comoving-frame (dotted lines) number densities are plotted versus
position. In the bottom right panel, the Eulerian-frame (solid lines) and
comoving-frame (dotted lines) energy densities are plotted versus position.

Figure 8 displays steady state solution characteristics for models where we
have varied $v_{\max}\in\\{0.0,0.1,0.2,0.4\\}$. From the top left panel, we
see that the $\mathcal{O}(v)$ spectra (dotted) — which for $v_{\max}>0$ are
redshifted relative to the case with $v_{\max}=0$ — agree well with the
analytic, special relativistic results (solid) for the lower values of
$v_{\max}$, while the difference between the $\mathcal{O}(v)$ and the special
relativistic results are larger when $v_{\max}=0.4$, which is to be expected
since $\mathcal{O}(v^{2})$ terms are no longer negligible. The top right panel
shows the RMS energy, defined as

$\varepsilon_{\rm
RMS}=\sqrt{\int_{D_{\varepsilon}}\mathcal{D}\varepsilon^{5}d\varepsilon/\int_{D_{\varepsilon}}\mathcal{D}\varepsilon^{3}d\varepsilon},$
(139)

versus position, computed from the numerical, $\mathcal{O}(v)$ solution and
the analytic, special relativistic solution in Eq. (138). Indeed, at
$x^{1}=5$, the relative difference in $\varepsilon_{\rm RMS}$ between the
$\mathcal{O}(v)$ and the special relativistic result for $v_{\max}=0.1$ is
$4.8\times 10^{-3}$, while it is $2.0\times 10^{-2}$ and $9.1\times 10^{-2}$
for $v_{\max}=0.2$ and $v_{\max}=0.4$, respectively. That is, the relative
error in the RMS energy increases roughly as $\mathcal{O}(v^{2})$.

The bottom panels of Figure 8 show the Eulerian- and comoving-frame number
densities ($N$ and $D$, respectively; left), and the Eulerian- and comoving-
frame energy densities ($E$ and $J$, respectively; right) versus position.
Here,

$\big{\\{}\,N,\,D,\,E,\,J\,\big{\\}}=4\pi\int_{D_{\varepsilon}}\big{\\{}\,\mathcal{N},\,\mathcal{D},\,\mathcal{E},\,\mathcal{D}\varepsilon\,\big{\\}}\varepsilon^{2}d\varepsilon,$
(140)

where $\mathcal{N}$ and $\mathcal{E}$ are defined in Eqs. (4) and (5),
respectively. Relative to where $v=0$, both $D$ and $J$ are lower in the
region where $v>0$, which is expected from the redshifted spectra displayed in
the top left panel in Figure 8. The Eulerian-frame quantities, $N$ and $E$,
are practically unaffected by the velocity field, and remain relatively
constant throughout the spatial domain.

####  8.4.2 Simultaneous Number and Energy Conservation

Next, we investigate the simultaneous number and energy conservation
properties of the scheme as a function of time for $t\in[0,20]$. Here, a main
challenge stems from the fact that, since the flux factor $h\approx 1$, the
moments evolve close to the boundary of the realizable set $\mathcal{R}$. With
high-order polynomials ($k\geq 1$), the solution can become non-realizable in
one or more quadrature points in some elements, which then triggers the
realizability-enforcing limiter discussed in Section 5.2. The realizability-
enforcing limiter preserves the Eulerian-frame particle number, but not the
Eulerian-frame energy, which is the reason for introducing the ‘energy
limiter’ in Section 6.2. Here, we demonstrate the performance of the energy
limiter, and its effect on the simultaneous number and energy conservation
properties of the method. Recall from Proposition 1 that the $\mathcal{O}(v)$
two-moment model is conservative for the Eulerian-frame energy only to
$\mathcal{O}(v^{2})$.

In the context of the current test, the Eulerian-frame number density
satisfies a conservation law of the form

$\partial_{t}N+\partial_{1}F_{N}^{1}=0.$ (141)

Integration over space $D_{x^{1}}$ and time $[t_{0},t]$ gives

$\underbrace{\int_{D_{x^{1}}}[\,N(x^{1},t)-N(x^{1},0)]\,dx^{1}}_{\Delta N_{\rm
int}(t)}+\underbrace{\int_{0}^{t}[F_{N}^{1}|_{x^{1}=10}-F_{N}^{1}|_{x^{1}=0}]\,d\tau}_{\Delta
N_{\rm ext}(t)}=0,$ (142)

where $\Delta N_{\rm int}(t)$ and $\Delta N_{\rm ext}(t)$ represent the change
in the total number of particles, from $t_{0}$ to $t$, _interior_ and
_exterior_ to the domain $D_{x^{1}}$, respectively. Since there is no creation
or destruction of particles, the sum vanishes. We can obtain a similar
expression for the Eulerian-frame energy (with $E$ replacing $N$ in Eq.
(142)), but for the $\mathcal{O}(v)$ two-moment model considered here, by
Proposition 1, one would in general expect

$\Delta E_{\rm int}+\Delta E_{\rm ext}=\mathcal{O}(v^{2}),$ (143)

at the continuous level. At the discrete level, with consistent discretization
of the left-hand side of the two-moment model, Eulerian-frame energy
violations of $\mathcal{O}(v^{2})$ should be considered optimal. (Recall the
discussion on this issue specific to the DG scheme from Section 6.1.) For this
test, given our chosen spatial resolution and use of quadratic elements,
velocity jumps across elements are small, and we expect to observe near
optimal Eulerian-frame energy conservation properties. However, the acceptable
level of Eulerian-frame energy nonconservation is application dependent, and
should be considered on a case-by-case basis.

|
---|---
Figure 9: Plot of Eulerian-frame number (left) and energy (right) balance
versus time for the streaming Doppler shift problem for a run with
$v_{\max}=0.2$. In the left panel, $\Delta N_{\rm int}$ (solid), $\Delta
N_{\rm ext}$ (dash-dotted), and $\Delta N_{\rm int}+\Delta N_{\rm ext}$
(dotted), see Eq. (142), are plotted. Similarly, in the right panel, $\Delta
E_{\rm int}$ (solid), $\Delta E_{\rm ext}$ (dash-dotted), and $\Delta E_{\rm
int}+\Delta E_{\rm ext}$ (dotted), see Eq. (143), are plotted.

Figure 9 plots the number and energy balance versus time for a run with
$v_{\max}=0.2$. Initially, there are essentially no particles in the
computational domain $D_{x^{1}}$, and the flux at the outer boundary is zero.
For $t>0$, particles enter the domain through the inner boundary, and $\Delta
N_{\rm int}$ begins to increase linearly with time, while $\Delta N_{\rm ext}$
decreases at the same rate, and the sum $\Delta N_{\rm int}+\Delta N_{\rm
ext}$ remains zero to machine precision (see also Figure 10). Around $t=10$,
the particles that entered the domain at $t=0$ reach the outer boundary,
establishing a balance between particles entering and leaving the domain, and
the system reaches a steady state where both $\Delta N_{\rm int}$ and $\Delta
N_{\rm ext}$ remain unchanged. The evolution observed for the Eulerian-frame
energy quantities ($\Delta E_{\rm ext}$ and $\Delta E_{\rm ext}$) is similar
to that for the particle number, and, on the scale of the ordinate on the
right panel in Figure 9, the sum $\Delta E_{\rm ext}+\Delta E_{\rm ext}$
remains close to zero.

|
---|---
|
Figure 10: Results from the streaming Doppler shift problem for models with
various values of $v_{\max}$. In the upper left panel the relative change in
the Eulerian-frame particle number, defined as $(\Delta N_{\rm int}+\Delta
N_{\rm ext})/\Delta N_{\rm int}(t=20)$, is plotted versus time. Similarly, in
the upper right panel, the relative change in the Eulerian-frame energy,
defined as $(\Delta E_{\rm int}+\Delta E_{\rm ext})/\Delta E_{\rm int}(t=20)$,
is plotted. In these panels, solid lines represent results obtained with the
fiducial algorithm with the energy limiter on, while dotted lines represent
results obtained with the energy limiter turned off. In the lower left panel,
the minimum (solid) and mean (dashed) value of the limiter parameter
$\theta_{\boldsymbol{K}}^{\boldsymbol{\mathcal{U}}}$ (see Algorithm 1) are
plotted versus time for the fiducial models with the energy limiter on.
Results with $v_{\max}$ set to $0.1$, $0.2$, and $0.3$ are plotted with red,
blue, and magenta curves, respectively. In the lower right panel, the relative
change in the Eulerian-frame energy at $t=20$ is plotted versus $v_{\max}$ for
models where the energy limiter is on (solid lines) and off (dotted lines)
(the dashed black reference line is proportional to $v_{\max}^{2}$).

Figure 10 shows select results from runs with
$v_{\max}\in\\{0.05,0.1,0.2,0.3,0.4\\}$ and provides further details on the
simultaneous number and energy conservation properties of the scheme when
applied to the streaming Doppler shift problem. For comparison, we have also
run all the models with the energy limiter turned off. The relative change in
the total number of particles (top left panel) is at the level of machine
precision for all values of $v_{\max}$ (with and without the energy limiter),
which is expected since the two-moment model is formulated in number
conservative form.

As mentioned earlier, the moments evolve close to the boundary of the
realizable set in this test, and the realizability-enforcing limiter is
frequently invoked to enforce pointwise realizability in all elements. From
the lower left panel in Figure 10, we observe that the minimum value of the
limiter parameter $\theta_{\boldsymbol{K}}^{\boldsymbol{\mathcal{U}}}$ (solid
lines) varies between $0$ and $0.1$ for $t\lesssim 10$, while the average
value of $\theta_{\boldsymbol{K}}^{\boldsymbol{\mathcal{U}}}$ (dashed lines)
grows from about $0.8$ to $1$ during this time. For $t\leq 10$, a
discontinuity in the moments, driven by the inner boundary condition,
propagates through the domain and is mainly responsible for triggering the
realizability-enforcing limiter. Recall that
$\theta_{\boldsymbol{K}}^{\boldsymbol{\mathcal{U}}}=0$ implies full limiting
where all the moments within an element are set to their respective cell
average, while $\theta_{\boldsymbol{K}}^{\boldsymbol{\mathcal{U}}}=1$ implies
no limiting. After $t\approx 10$, the average
$\theta_{\boldsymbol{K}}^{\boldsymbol{\mathcal{U}}}$ value hovers around
unity, but a few elements still require significant limiting when $t>10$,
especially for the $v_{\max}=0.4$ model, with the minimum
$\theta_{\boldsymbol{K}}^{\boldsymbol{\mathcal{U}}}$ still dropping down close
to zero. Closer inspection reveals that a few elements around the location of
the velocity gradients ($x^{1}\in[2,3.5)$ and $x^{1}\in[6.5,8)$) require
limiting beyond $t=10$.

The relative change in the Eulerian-frame energy is plotted in the top right
panel of Figure 10, which reveals a significant improvement in conservation
for the fiducial models with the energy limiter turned on (solid lines), when
compared to models without the energy limiter (dotted lines). (We compute the
relative change in number and energy by normalizing by interior values at
$t=20$, when the system is in steady state, because the initial interior
values are close to zero.) For the models without the energy limiter, the
relative change in total energy immediately jumps to about $1\times 10^{-5}$,
and continues to grow for later times, while it remains relatively constant
for the fiducial models. This implies that the realizability-enforcing limiter
is the main driver of Eulerian-frame energy nonconservation for this test, and
not the inherent nonconservation properties of the continuum $\mathcal{O}(v)$
two-moment model. (Velocity jumps across elements are small, and we infer from
our results that their contribution to energy nonconservation is negligible.)
The relative change in the Eulerian-frame energy at $t=20$ is plotted versus
$v_{\max}$ in the lower right panel of Figure 10. For the models with the
energy limiter turned off, the relative change is essentially independent of
$v_{\max}$, while Eulerian-frame energy violations grow as $v_{\max}^{2}$ for
the fiducial models. Note that the energy limiter only recovers energy
conservation violations caused by the realizability-enforcing limiter. Energy
conservation violations caused by the use of the $\mathcal{O}(v)$ two-moment
model are unaffected by the energy limiter. Since we observe the
$v_{\max}^{2}$ scaling when using the energy limiter, we posit that the DG
discretization maintains consistency with the continuum model on this aspect.

###  8.5 Transparent Shock

In this test, we investigate the performance of the method when the background
velocity gradient is varied. We consider a one-dimensional spatial domain
$D_{x^{1}}=[0,2]$, set the opacities $\chi=\sigma=0$, and the velocity field
$\boldsymbol{v}=(v,0,0)^{\intercal}$, where

$v(x^{1})=\frac{1}{2}v_{\max}\times\big{[}\,1-\tanh\big{(}\,(x^{1}-1)/H\,\big{)}\,\big{]}.$
(144)

We will vary both the velocity magnitude $v_{\max}$ and gradient, parametrized
by the length scale $H$. The energy domain is again $D_{\varepsilon}=[0,50]$.
We discretize the spatial and energy domains using $80$ and $32$ elements,
respectively, and use quadratic elements ($k=2$) and SSPRK3 time stepping.
Then, $\Delta x^{1}=0.025$, and $\Delta x^{1}/H=5/6$, $2.5$, and $25$, for
$H=3\times 10^{-2}$, $10^{-2}$, and $10^{-3}$, respectively. We use the same
boundary conditions as in the Doppler shift test, and the moments are
initially set to $\mathcal{D}=1\times 10^{-8}$ and $\mathcal{I}^{1}=0$ for all
$(x^{1},\varepsilon)\in D_{x^{1}}\times D_{\varepsilon}$. With the given
initial and boundary conditions, the equations are integrated until $t=3$,
when an approximate steady state has been established.

Figure 11 shows velocity profiles and comoving-frame and Eulerian-frame number
densities versus position around the ‘shock’ ($x^{1}\in[0.9,1.1]$) for the
different values of $H$ for the case with $v_{\max}=-0.1$. (The markers
indicate locations of LG quadrature points in each spatial element.) For
$H=3\times 10^{-2}$, the shock is resolved by the spatial grid, for
$H=10^{-2}$ it is under-resolved, while for $H=10^{-3}$, the velocity profile
is discontinuous.

The comoving-frame number densities increase across the velocity gradient
because of the Doppler effect, increasing the particle energy measured by the
comoving observer, who is moving towards the inner boundary. Beyond the shock,
the values for the computed comoving-frame number densities (solid lines) are
about $0.5~{}\%$ higher than the analytic values obtained using Eq. (138)
(dashed lines), and this fact is independent of the value of $H$. The
Eulerian-frame number densities, which should remain unaffected by the
presence of the background velocity, are essentially constant across the
shock. These results indicate that the method is able to capture Doppler
shifts correctly, even when velocity gradients are large.

|
---|---
Figure 11: Results from the transparent shock problem for $v_{\max}=-0.1$ and
various values of the shock width parameter $H$, defined in Eq. (144). In the
left panel, velocity profiles are plotted versus position for $H=3\times
10^{-2}$, $10^{-2}$, and $10^{-3}$ (black, red, and blue curves,
respectively). In the right panel, the numerical and analytic (special
relativistic) comoving-frame number densities (solid lines with markers and
dashed lines, respectively), and the numerical Eulerian-frame number densities
(dotted lines with markers) are plotted versus position. For the results
displayed in the right panel, the line colors correspond to the velocity
profile with the matching color plotted in the left panel.

In Figure 12, the left panel displays the relative change in the Eulerian-
frame energy, defined by the left-hand side of Eq. (143), versus $|v_{\max}|$
at $t=3$, for various values of $H$. The right panel displays the relative
error in $\varepsilon_{\rm RMS}$. Both panels display results obtained with
and without the energy limiter. (The relative change in the Eulerian-frame
particle number, not shown, is at the level of machine precision for all
models.) Considering the results obtained with the energy limiter active, in
the left panel we observe that, for a given value of $H$, the relative change
in the total energy increases with increasing $|v_{\max}|$, roughly
proportional to $|v_{\max}|^{2}$. Also, for a given value of $|v_{\max}|$, the
relative change in total energy increases with decreasing shock width $H$. For
the models with the energy limiter turned off, the behavior is different in a
large region of the $(v_{\max},H)$-space. With $H=3\times 10^{-2}$ (dotted
black line), the relative change in the Eulerian-frame energy at $t=3$ is
around $2\times 10^{-5}$; independent of $v_{\max}$. This can be attributed to
the realizability-enforcing limiter. For the models with $H=10^{-2}$ (dotted
red line), the energy change is roughly constant until $|v_{\max}|=0.1$, when
the relative energy change begins to increase with $|v_{\max}|$ in a manner
similar to the models with the energy limiter activated (solid red line). The
models with the steepest velocity gradient ($H=10^{-3}$; dotted blue line)
follow the corresponding models with active energy limiter, and the relative
change in the Eulerian-frame energy increases as $|v_{\max}|^{2}$ for all
$|v_{\max}|$. From this we conclude that the energy limiter can help to
recover the $\mathcal{O}(v^{2})$ Eulerian-frame energy conservation property
of the $\mathcal{O}(v)$ two-moment model for small velocities and velocity
gradients. The relative error in $\varepsilon_{\rm RMS}$, which increases as
$|v_{\max}|^{2}$ for all models, is essentially unaffected by the energy
limiter.

|
---|---
Figure 12: Results for the transparent shock problem, plotted versus
$|v_{\max}|$. The left panel displays the relative change in the Eulerian-
frame energy at $t=3$ for models where the energy limiter is on (solid lines)
and off (dotted lines). The right panel displays the absolute relative
difference in $\varepsilon_{\rm RMS}$ with respect to the exact solution from
special relativity at $x^{1}=2$. In both panels, black, red, and blue curves
correspond to $H=3\times 10^{-2}$, $10^{-2}$, and $10^{-3}$, respectively, and
the dashed references line are proportional to $|v_{\max}|^{2}$.

###  8.6 Transparent Vortex

The final test, inspired by the test in Section 4.2.3 of [3], considers
evolution in a two-dimensional spatial domain $D_{x^{1}}\times
D_{x^{2}}=[-5,5]\times[-5,5]$. We set the opacities $\chi=\sigma=0$, and the
velocity field is given by $\boldsymbol{v}=[v^{1},v^{2},0]^{\intercal}$, where

$\displaystyle v^{1}(x^{1},x^{2})$
$\displaystyle=-v_{\max}\,x^{2}\,\exp\big{[}\,(1-r^{2})/2\,\big{]},$ (145a)
$\displaystyle v^{2}(x^{1},x^{2})$ $\displaystyle=\hskip
8.0ptv_{\max}\,x^{1}\,\exp\big{[}\,(1-r^{2})/2\,\big{]},$ (145b)

and $r=\sqrt{(x^{1})^{2}+(x^{2})^{2}}$. The energy domain is
$D_{\varepsilon}=[0,50]$, and we discretize the spatial and energy domains
using $48\times 48$ and $32$ elements, respectively. We use quadratic elements
($k=2$) and SSPRK3 time stepping. The upper left panel in Figure 13 shows the
velocity field for the case with $v_{\max}=0.1$. The main purpose of this test
is to investigate the robustness of the method in configurations where the
radiation field propagates through a spatially variable velocity field with
various relative angles between the radiation flux and velocity vectors. The
moments are initially set to $\mathcal{D}=1\times 10^{-8}$ and
$\mathcal{I}^{1}=\mathcal{I}^{2}=0$ for all $(x^{1},x^{2},\varepsilon)\in
D_{x^{1}}\times D_{x^{2}}\times D_{\varepsilon}$. At the inner $x^{1}$
boundary, we impose an incoming, radiation field with a Fermi-Dirac spectrum:
We set
$\mathcal{D}(\varepsilon,x^{1}=-5,x^{2})=0.05/[\exp(\varepsilon/3-3)+1]$,
$\mathcal{I}^{1}(\varepsilon,x^{1}=-5,x^{2})=0.95\times\mathcal{D}(\varepsilon,x^{1}=-5,x^{2})$,
and $\mathcal{I}^{2}(\varepsilon,x^{1}=-5,x^{2})=0$, so that the flux factor
is $h=0.95$. With these initial and boundary conditions, the moment equations
are integrated until a steady state is reached ($t=20$).

|
---|---
|
Figure 13: Results for the transparent vortex problem. In the upper left
panel, the magnitude of the velocity, for the case with $v_{\max}=0.1$, is
displayed in grayscale with velocity vectors overlaid. The black, magenta,
red, and blue markers indicate spatial positions for which we plot numerical
energy spectra in solid lines in the upper right panel, with line colors
corresponding to the marker colors in the upper left panel. The dotted lines
are analytic spectra obtained from special relativistic considerations using
the local three-velocity. In the lower left panel, the relative change in
Eulerian-frame energy is plotted versus time for models with $v_{\max}$ in Eq.
(145) set to $0.01$ (black lines), $0.03$ (red lines), and $0.1$ (blue lines).
Results obtained with and without the energy limiter are plotted using solid
and dotted lines, respectively. The lower right panel plots the relative
difference between the incoming and outgoing particle fluxes in the
$x^{1}$-direction.

In the upper right panel in Figure 13, the solid lines represent numerical
energy spectra at spatial locations indicated with solid markers of matching
color in upper left panel. At the location of the black marker the velocity is
close to zero, and thus the black line represents the spectrum of the incoming
radiation. The red and blue markers are located where
$\boldsymbol{v}=(v_{\max},0,0)^{\intercal}$ and
$\boldsymbol{v}=(-v_{\max},0,0)^{\intercal}$, respectively, and the spectra at
these locations are, respectively, red- and blue-shifted relative to the
spectrum sampled at the location of the black marker. Analytic spectra at the
locations of the black, red, and blue markers are plotted with dotted lines,
which indicate good agreement between numerical and analytical solutions
across all energies. At the locations of the black, red, and blue markers, we
find that $\varepsilon_{\rm RMS}$ is approximately $15.6$, $14.2$, and $17.2$,
respectively. At the location of the magenta marker, which is placed on the
opposite side of the vortex (relative to the black marker), the velocity is
again close to zero, and it is expected that the spectrum at this location
agrees with the spectrum at the location of the black marker. Comparing the
solid black and magenta lines in the upper right panel, we observe that the
spectral number density is consistently higher at the location of the magenta
marker (by a constant factor of about $1.07$). Comparing $\varepsilon_{\rm
RMS}$ at the two locations, we find that the relative difference is less than
$10^{-4}$.

The lower left panel in Figure 13 plots the relative change in the Eulerian-
frame energy versus time for models with $v_{\max}\in\\{0.01,0.03,0.1\\}$.
Results obtained with the energy limiter on are plotted with solid lines,
while dotted lines correspond to results with the energy limiter turned off.
For all models, the relative change in the Eulerian-frame energy is less than
$10^{-4}$. For the models with $v_{\max}=0.1$, the relative change reaches the
largest amplitudes for $t\in[4,7]$, when a radiation front, driven by the
boundary condition at $x^{1}=-5$, propagates through the vortex. The model
with the energy limiter makes a better recovery than the corresponding model
with the energy limiter turned off. For smaller $v_{\max}$, the relative
change in the Eulerian-frame energy is clearly much smaller when the energy
limiter is active. These results demonstrate the contribution to Eulerian-
frame energy nonconservation caused by the realizability-enforcing limiter.
For both suites of models (energy limiter on or off), the relative change in
the Eulerian-frame number (not shown) is at the level of machine precision for
all models.

The lower right panel in Figure 13, similar to Figure 6 (b) in [3], shows, for
$t=20$, the relative difference between the energy integrated
$x^{1}$-component of the number flux densities evaluated at the inner and
outer boundaries of $D_{x^{1}}$, defined as
$|I^{1}(5,x^{2})-I^{1}(-5,x^{2})|/I^{1}(-5,x^{2})$. As discussed by Just et
al. [3], this quantity should vanish for exact calculations, while errors of
$\mathcal{O}(v^{2})$ are to be expected for the $\mathcal{O}(v)$ two-moment
model. Comparing with their results, the curves plotted in our figure share
similar features. Moreover, for $v_{\max}=0.01$, the maximum relative
difference is $6.15\times 10^{-4}$, for $v_{\max}=0.03$ it is $4.87\times
10^{-3}$, while it is $5.68\times 10^{-2}$ for $v_{\max}=0.1$; i.e., the
maximum error grows as $v_{\max}^{2}$.

Despite the growing (with $v_{\max}$) relative difference between the number
fluxes at the inner and outer boundaries in the $x^{1}$-direction, we point
out that, due to number conservation, the integrated number fluxes through the
inner and outer boundaries balance each other. That is, in the steady state at
$t=20$,
$\int_{D_{x^{2}}}I^{1}(-5,x^{2})\,dx^{2}=\int_{D_{x^{2}}}I^{1}(5,x^{2})\,dx^{2}$.
However, the distribution of particles along the $x^{2}$-direction becomes
nonuniform in the wake of the vortex, while a uniform distribution is expected
as $|\boldsymbol{v}|\to 0$. We illustrate this further in Figure 14. The left
panel shows that, within the vortex ($r\lesssim 2$), the comoving-frame number
density is higher than the reference value $D_{0}$ for $x^{2}>0$, and lower
than $D_{0}$ for $x^{2}<0$, which is consistent with the Doppler shift of the
spectra in the respective regions. In the wake of the vortex, the comoving-
frame number density is relatively higher in the region centered around
$x^{2}=0$, while it is lower further away (compare red and blue regions for
$x^{1}\gtrsim 2$ in the left panel in Figure 14). The Eulerian-frame number
density is relatively unaffected by the vortex for $x^{1}<0$, but exhibits a
spatial distribution similar to the comoving-frame number density in the wake.
In contrast, the spatial distribution of the RMS energy is more consistent
with expectations: Within the vortex, $\varepsilon_{\rm RMS}>\varepsilon_{{\rm
RMS},0}$ for $x^{2}>0$, while $\varepsilon_{\rm RMS}<\varepsilon_{{\rm
RMS},0}$ for $x^{2}<0$. Moreover, the RMS energy returns to the reference
value in the wake of the vortex, with almost uniform distribution along the
$x^{2}$-direction. We do not have a complete theoretical explanation for the
spatial distribution of the number densities in the wake of the vortex, but
suspect that the two-moment approximation and the associated closure, which
assumes that the radiation field is axisymmetric about a preferred direction
in momentum space [46], is insufficient for capturing relativistic aberration
effects.

|  |
---|---|---
Figure 14: Results for the transparent vortex problem at $t=20$ for a model
with $v_{\max}=0.1$. The left panel shows the relative deviation in comoving-
frame number density from $D_{0}=D(x^{1}=-5,x^{2})$, $(D-D_{0})/D_{0}$, with
vectors of the comoving-frame number flux
$(I^{1}-I_{0}^{1},I^{2})^{\intercal}$ overlaid, where
$I_{0}^{1}=I^{1}(x^{1}=-5,x^{2})$ is the first component of the comoving-frame
number flux density at the inner boundary in the $x^{1}$-direction, which is
subtracted to better illustrate the flow, since $|I^{2}|\ll|I^{1}|$ generally
holds. Similarly, the middle panel shows the corresponding relative deviation
in the Eulerian-frame number density $(N-N_{0})/N_{0}$, with vectors of the
Eulerian-frame number flux $(F_{N}^{1}-F_{N,0}^{1},F_{N}^{2})^{\intercal}$.
The right panel shows the relative deviation in the RMS energy,
$(\varepsilon_{{\rm RMS},0}-\varepsilon_{\rm RMS})/\varepsilon_{{\rm RMS},0}$,
where $\varepsilon_{{\rm RMS},0}=\varepsilon_{\rm RMS}(x^{1}=-5,x^{2})$.

###  8.7 Performance Evaluation

To demonstrate the GPU functionality and performance characteristics of the
DG-IMEX method as implemented in thornado, we consider the Streaming Doppler
Shift test, described in Section 8.4, with $v_{\max}=0.1$. To more accurately
capture a production workload, the tests are performed in three spatial
dimensions, with the number of elements similar to what would be used for a
single process invoking thornado in a multiphysics simulation. The benchmark
is run in two configurations, using tensor product polynomials of degree $k=1$
and $k=2$, respectively. The SSPRK2 time stepper is used for both
configurations. For $k=1$, we use 16 energy elements and $96\times 3\times 3$
spatial elements, while 12 energy elements and $64\times 2\times 2$ spatial
elements are used for $k=2$ — thus keeping the total number of spatial degrees
of freedom the same. Our goal is to provide a high-level demonstration of
performance characteristics and the relative cost of main algorithmic
components, while we defer a rigorous performance analysis to future work.

The tests are performed on a single node of the Summit computer at the Oak
Ridge Leadership Computing Facility (OLCF). Each Summit node has 2 IBM POWER9
CPUs and 6 NVIDIA V100 GPUs, but here we limit our benchmarks to a single CPU
or GPU. For the CPU runs, we use seven cores with one thread per core as this
is the number of cores that would be available to one process if we divide the
resources equally with one GPU per process. All runs use version 22.5 of the
NVIDIA nvfortran compiler with standard -O2 optimizations. Optimized linear
algebra libraries are provided by IBM ESSL (v6.3.0) on the CPU and NVIDIA
cuBLAS (v11.0.3) on the GPU. For the GPU runs, all computations are done on
the GPU using OpenACC and libraries; the CPU process is only used to launch
kernels and manage data transfer. In both cases, the salient metric is wall-
time per time step (lower is better).

Figure 15 shows a breakdown of the relative cost associated with evaluating
the major components of the explicit phase-space advection operator. The
polynomial degree has little effect on the absolute wall-time, especially for
the GPU runs. For the CPU runs, the relative cost of linear algebra (MatMul)
is somewhat higher when $k=2$. As can be seen comparing the right and left
panels, the initial guess in the conserved-to-primitive calculation can have a
non-trivial impact on the total wall-time by reducing the total number of
solver iterations. We measure a total speedup factor of 8–10 for the V100
relative to the multi-core CPU runs on the POWER9. Notably, the relative cost
for linear algebra and limiters becomes negligible when using the GPU, and the
majority of the computational cost is shifted to the iterative conserved-to-
primitive calculations. We speculate that one approach to further improve the
performance would be to combine the calculation of all of the primitive
moments on the quadrature set $\widetilde{S}_{\otimes}^{\boldsymbol{K}}$,
defined in Eq. (52), into a single kernel, rather than to calculate them
separately for each evaluation of $\boldsymbol{\mathcal{F}}^{i}$ and
$\boldsymbol{\mathcal{F}}^{\varepsilon}$, defined in Eqs. (39)–(40), which
results in some duplicate evaluations. While these savings may be significant
for the phase-space advection problem considered here, refactoring will be
considered in the context of a more physics-complete implementation. With more
realistic collision terms included, the relative cost of the explicit phase-
space advection part is expected to be small (see, e.g., [55, 35]).

|
---|---
Figure 15:  Breakdown of normalized wall-time for components of the Streaming
Doppler Shift test problem as implemented in thornado. The left panel uses an
initial guess of
$\boldsymbol{\mathcal{M}}^{[0]}=(\mathcal{N},\boldsymbol{0})^{\intercal}$ in
the conserved-to-primitive calculation, and the right panel uses
$\boldsymbol{\mathcal{M}}^{[0]}=\boldsymbol{\mathcal{U}}=(\mathcal{N},\boldsymbol{\mathcal{G}})^{\intercal}$.
Absolute wall-clock times per time step are shown above each bar. Flux
captures the calculation of fluxes $\boldsymbol{\mathcal{F}}^{i}$ and
$\boldsymbol{\mathcal{F}}^{\varepsilon}$ in Eqs. (39)–(40). MatMul represents
the matrix-matrix multiplications used throughout the explicit operator, e.g.,
to evaluate polynomials $\boldsymbol{\mathcal{U}}_{h}$ in quadrature points on
all elements. Primitive captures the iterative conserved-to-primitive
calculation described in Section 4.3.1. Limiters includes the application of
the realizability-enforcing limiter described Section 5.2 and the energy
limiter described in Section 6.2. Other is used to capture all remaining wall-
time spent in the explicit step.

##  9 Summary and Conclusions

We have proposed and analyzed a realizability-preserving numerical method for
evolving a spectral two-moment model for neutral particles interacting with a
moving background fluid. This number-conservative moment model is based on
comoving-frame momentum coordinates, includes special relativistic corrections
to $\mathcal{O}(v)$, and, as a result, contains velocity-dependent terms
accounting for spatial advection, Doppler shift, and angular aberration. The
nonlinear two-moment model solves for comoving-frame angular moments,
representing number density and components of the number flux densitiy, and is
closed by expressing higher-order moments (rank-two and rank-three tensors) in
terms of the evolved moments using the maximum entropy closure (both exact and
approximate) due to Minerbo [19]. The two-moment model is closely related to
that promoted in [23], predicts wave speeds bounded by the speed of light
(Proposition 2), and is consistent, to $\mathcal{O}(v)$, with Eulerian-frame
energy and momentum conservation (Proposition 1).

The numerical method is based on the DG method to discretize phase-space, and
IMEX time stepping, where the phase-space advection part is integrated with
explicit methods, and the collision term is integrated with implicit methods.
The discretized spatial and energy derivative terms in the moment equations
have been equipped with tailored numerical fluxes, which in the case of exact
moment closure (Assumption 1) allow us to derive explicit time-step
restrictions that guarantee realizable cell-averaged moments due to these
terms, and $\mathcal{N}>0$ overall. Unfortunately, a corresponding time-step
restriction was not found for the source terms associated with phase-space
advection in the number flux equation to guarantee the second moment
realizability condition, in the sense of cell averages, for the evolved
moments (i.e., $|\boldsymbol{\mathcal{G}}|\leq\mathcal{N}$) in the general
multidimensional case. However, an analysis in the semi-discrete setting
revealed that the moments evolve tangentially to the boundary of the
realizable domain when $|\boldsymbol{\mathcal{G}}|=\mathcal{N}$, and we found
a sufficient time-step restriction to guarantee realizable cell averages in
the one-dimensional, planar geometry case. Given a positive cell-averaged
number density, a realizability-enforcing limiter is proposed to recover
pointwise moment realizability in each element. Specific properties of the
IMEX scheme (i.e., convex-invariance, as defined in [34]) extend the
applicability of our results beyond the forward-backward Euler sequence
analyzed in detail.

Retention of specific $\mathcal{O}(v)$ terms in the time derivative of the
two-moment system, motivated by the desire to maintain wave speeds bounded by
the speed of light and consistency with Eulerian-frame energy and momentum
conservation equations, results in increased computational complexity of the
numerical scheme in two (related) ways. First, since the evolved moments are
nonlinear functions of the primitive moments used to close the moment
equations, a nonlinear system must be solved to recover primitive moments from
evolved moments. Second, because the collision operators are formulated in
terms of primitive moments, the implicit collision update requires the
solution of a similar nonlinear system. For both cases, solution methods have
been formulated as fixed-point problems, and we have proposed tailored fixed-
point operators in Eqs. (56) and (60), for the primitive recovery and implicit
collision solve, respectively. The fixed-point operators are designed to
preserve moment realizability in each iteration (subject to mild conditions on
the step size), and we have proven convergence for cases with exact _and_
approximate moment closures, subject to the additional constraint
$|\boldsymbol{v}|\leq\sqrt{2}-1$, which is mild when considering the
applicability of the $\mathcal{O}(v)$ model. Numerically, we did not observe
convergence failures for the primitive recovery problem, even when violating
the condition on the velocity, or when combining the algorithm with Anderson
acceleration, which our analysis here did not consider.

The proposed algorithm has been implemented and tested against a series of
benchmark problems. Using two problems with a constant background velocity —
in the streaming and diffusion regimes, respectively — we demonstrate the
expected rate of error convergence in the $L^{2}$ norm. Additional tests with
spatially varying (smooth and discontinuous) background velocity fields — the
Streaming Doppler Shift, Transparent Shock, and Transparent Vortex tests —
were used to document the robustness of the proposed algorithm, and
qualitative accuracy with respect to special relativistic considerations
(e.g., correct Doppler shifts) for sufficiently small background velocities.
In these tests, the moments evolve close to the boundary of the realizable
domain, and the realizability-enforcing limiter is frequently triggered to
recover pointwise realizability from (guaranteed) realizable cell averages.
Without this recovery procedure, the algorithm fails invariably on these
challenging problems.

We have analyzed the simultaneous Eulerian-frame number and energy
conservation properties of the proposed method. While the DG method provides
flexibility in the approximation spaces to capture conservation properties
beyond those inherent to the model formulation (i.e., number conservation in
the present setting), the approximation of the background velocity by
piecewise polynomials from the DG approximation space, which accommodates
discontinuities, can result in Eulerian-frame energy conservation errors that
exceed the $\mathcal{O}(v^{2})$ scaling predicted by the continuum model.
However, we found that the realizability-enforcing limiter is the main
contributor to Eulerian-frame energy conservation violations when the
background velocity field is smooth and its magnitude is within the range of
applicability of an $\mathcal{O}(v)$ model. For this reason, an energy limiter
is proposed to recover conservation violations introduced by the
realizability-enforcing limiter. This limiter trades local number conservation
for number _and_ energy conservation after integration over the phase-space
energy dimension, and has no observed negative impact on solution accuracy,
while improving Eulerian-frame energy conservation properties of the method.
With the energy limiter active, we observe that energy conservation violations
scale as $\mathcal{O}(v^{2})$, in accordance with the continuum model. We
emphasize that the energy limiter introduces a rescaling of the moments, which
does not impact moment realizability. However, the proposed strategy to
promote Eulerian-frame energy conservation is not feasible without the
realizability-preserving property.

Our goal is to apply the proposed algorithm to model neutrino transport in
core-collapse supernova simulations. Several extensions are needed to achieve
this goal. First, the collision term must be extended to include a complete
set of neutrino weak interactions, and the model extended to include coupling
to dynamical equations for the background fluid. Second, because neutrinos are
Fermions, for which the Pauli exclusion principle implies an upper bound on
the phase-space density and associated bounds on the moments, the analysis
should be extended to apply to moment closures based on Fermi-Dirac
statistics. Third, because special _and_ general relativistic effects
contribute to the dynamics in nontrivial ways, further development and
analysis is required to design realizability-preserving methods for fully
relativistic moment models. We believe the methodologies developed in this
paper can be helpful in these endeavors, and hope to present progress on
addressing these challenges in future work.

##  Appendix A Technical Proofs

###  A.1 Various Bounds for the Exact and Approximate Eddington Factors

In the following lemma, we list several bounds on functions dependent on the
exact or approximate Eddington factors ($\psi$ or $\psi_{\mathsf{a}}$). These
bounds are used in the proofs of Lemmas 8 and 9 in A.2 and A.3, respectively,
as well of the proof of Lemma 11 in Section 5.5.

######  Lemma 12.

Let $\psi$ be the Eddington factor in the exact Minerbo closure as given in
Eq. (26) and let

$\phi_{1}:=3\psi-1-3\psi^{\prime}h\quad\text{and}\quad\phi_{2}:=(3\psi-1)h^{-1}\>.$
(146)

Then, the following bounds hold when $h\in[0,1]$.

  * (a)

$-4\leq\phi_{1}\leq 0$,

  * (b)

$\phi_{2}^{2}-\psi^{\prime}\phi_{2}\geq 0$,

  * (c)

$3(\psi^{\prime})^{2}-3\psi^{\prime}\phi_{2}+\phi_{2}^{2}\geq 0$,

  * (d)

$\partial_{h}(\phi_{2}^{2}-\psi^{\prime}\phi_{2}+(\psi^{\prime})^{2})>0$.

Moreover, Let $\psi_{\mathsf{a}}$ be the approximate Eddington factor defined
in Eq. (27) and let

$\phi_{\mathsf{a},1}:=3\psi_{\mathsf{a}}-1-3\psi_{\mathsf{a}}^{\prime}h\quad\text{and}\quad\phi_{\mathsf{a},2}:=(3\psi_{\mathsf{a}}-1)h^{-1}\>.$
(147)

Then the bounds (a)–(d) hold when $(\psi,\phi_{1},\phi_{2})$ are replaced by
$(\psi_{\mathsf{a}},\phi_{\mathsf{a},1},\phi_{\mathsf{a},2})$. In addition,
the following two bounds hold for the approximate Eddington factor when
$h\in[0,1)$.

  * (e)

$\psi_{\mathsf{a}}-h^{2}-\frac{1}{4}(1-\psi_{\mathsf{a}})^{2}\geq 0$,

  * (f)

$h^{2}\leq\psi_{\mathsf{a}}\leq 1$.

Since both $\psi$ and $\psi_{\mathsf{a}}$ are one-dimensional functions
defined between 0 and 1, the proof of the bounds are straightforward but are
rather tedious. Instead of giving rigorous proofs for these bounds, we plot
the functions of interest in Figure 16, from which the bounds can be visually
verified.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 16: Verification of the bounds (a)–(f) in Lemma 12 in the cases when
the exact Eddington factor $\psi$ and approximate $\psi_{\mathsf{a}}$ are
considered.

###  A.2 Proof of Lemma 8

###### Proof of Lemma 8..

Using the definition of the closure terms $\mathsf{k}_{ij}$ in Eq. (10), we
have from chain rule that

$\displaystyle v^{i}{\partial_{\mathcal{D}}(\mathsf{k}_{ij}\mathcal{D})}$
$\displaystyle=v^{i}\left(\frac{1}{2}\big{[}3\hat{n}_{i}\hat{n}_{j}-\delta_{ij}\big{]}\frac{\partial\psi}{\partial
h}\frac{\partial
h}{\partial\mathcal{D}}\mathcal{D}+\frac{1}{2}\big{[}(1-\psi)\delta_{ij}+(3\psi-1)\hat{n}_{i}\hat{n}_{j}\big{]}\right)$
$\displaystyle=v^{i}\left(-\frac{1}{2}\big{[}3\hat{n}_{i}\hat{n}_{j}-\delta_{ij}\big{]}\psi^{\prime}h+\frac{1}{2}\big{[}(1-\psi)\delta_{ij}+(3\psi-1)\hat{n}_{i}\hat{n}_{j}\big{]}\right)$
(148)
$\displaystyle=\frac{1}{2}(3\psi-1-3\psi^{\prime}h)\big{(}v^{i}\hat{n}_{i}\hat{n}_{j}-\frac{1}{3}v_{j}\big{)}+\frac{1}{3}v_{j}=\frac{1}{2}\phi_{1}\big{(}v^{i}\hat{n}_{i}\hat{n}_{j}-\frac{1}{3}v_{j}\big{)}+\frac{1}{3}v_{j}\>,$

where $\phi_{1}:=(3\psi-1-3\psi^{\prime}h)$ as defined in Eq. (146). Since
$\|\partial_{\mathcal{D}}(v^{i}\mathsf{k}_{ij}\mathcal{D})\|^{2}=\sum_{j}\left(v^{i}{\partial_{\mathcal{D}}(\mathsf{k}_{ij}\mathcal{D})}\right)^{2}$,
summing up the squares leads to

$\displaystyle\|\partial_{\mathcal{D}}(v^{i}\mathsf{k}_{ij}\mathcal{D})\|^{2}$
$\displaystyle=\frac{1}{4}\phi_{1}^{2}\sum_{j}\left(v^{i}\hat{n}_{i}\hat{n}_{j}-\frac{1}{3}v_{j}\right)^{2}+\frac{1}{3}\phi_{1}v^{j}\left(v^{i}\hat{n}_{i}\hat{n}_{j}-\frac{1}{3}v_{j}\right)+\frac{1}{9}v^{j}v_{j}$
$\displaystyle=\left(\frac{\phi_{1}^{2}}{12}+\frac{\phi_{1}}{3}\right)(v^{i}\hat{n}_{i})^{2}+\left(\frac{\phi_{1}^{2}}{36}-\frac{\phi_{1}}{9}+\frac{1}{9}\right)v^{i}v_{i}\>.$
(149)

It follows from Lemma 12 (a) that $\phi_{1}(h)\in[-4,0]$ for $h\in[0,1]$.
Therefore,
$\left(\frac{\phi_{1}^{2}}{12}+\frac{\phi_{1}}{3}\right)=\frac{\phi_{1}}{12}\left(\phi_{1}+4\right)\geq
0$ and we have

$\|\partial_{\mathcal{D}}(v^{i}\mathsf{k}_{ij}\mathcal{D})\|^{2}\leq\left(\big{(}\frac{\phi_{1}^{2}}{12}+\frac{\phi_{1}}{3}\big{)}+\big{(}\frac{\phi_{1}^{2}}{36}-\frac{\phi_{1}}{9}+\frac{1}{9}\big{)}\right)v^{i}v_{i}=\frac{1}{9}\left(\phi_{1}+1\right)^{2}v^{i}v_{i}\>.$
(150)

Since $\phi_{1}\in[-4,0]$, $\frac{1}{9}\left(\phi_{1}+1\right)^{2}\leq 1$,
which concludes the proof. ∎

###  A.3 Proof of Lemma 9

###### Proof of Lemma 9..

Using the definition of the closure terms $\mathsf{k}_{ij}$ in Eq. (10), we
have from chain rule that

$v^{i}{\partial_{\mathcal{I}_{k}}(\mathsf{k}_{ij}\mathcal{D})}=\frac{1}{2}\,{\psi^{\prime}}\,\Big{[}\,3\,s\,\hat{n}_{j}-v_{j}\,\Big{]}\,\hat{n}_{k}+\frac{(3\,\psi-1)}{2h}\,\Big{[}\,v_{k}\,\hat{n}_{j}+s\delta_{jk}\,-2\,s\,\hat{n}_{j}\,\hat{n}_{k}\,\Big{]}\>.$
(151)

To show
$\|\nabla_{\boldsymbol{\mathcal{I}}}(v^{i}\mathsf{k}_{ij}\mathcal{D})\|\leq
2v$, we prove in the following that

$\|\nabla_{\boldsymbol{\mathcal{I}}}(v^{i}\mathsf{k}_{ij}\mathcal{D})\boldsymbol{y}\|\leq
2vy\>,\quad\forall\boldsymbol{y}=\big{(}y^{1},y^{2},y^{3}\big{)}^{\intercal},$
(152)

where $y:=\sqrt{y^{i}y_{i}}\,$. Let $\phi_{2}:=(3\psi-1)h^{-1}$ as defined in
Eq. (146). Then,

$v^{i}{\partial_{\mathcal{I}_{k}}(\mathsf{k}_{ij}\mathcal{D})}y^{k}=\frac{1}{2}\,{\psi^{\prime}}\,\Big{[}\,3\,s\,\hat{n}_{j}-v_{j}\,\Big{]}\,(y^{k}\hat{n}_{k})+\frac{1}{2}\phi_{2}\,\Big{[}\,\hat{n}_{j}(y^{k}v_{k})\,+sy_{j}\,-2\,s\,\hat{n}_{j}\,(y^{k}\hat{n}_{k})\,\Big{]}\>.$
(153)

Summing up the squares leads to

$\displaystyle\|\nabla_{\boldsymbol{\mathcal{I}}}(v^{i}\mathsf{k}_{ij}\mathcal{D})\boldsymbol{y}\|^{2}=$
$\displaystyle\sum_{j}\left(v^{i}{\partial_{\mathcal{I}_{k}}(\mathsf{k}_{ij}\mathcal{D})}y^{k}\right)^{2}=\frac{1}{4}\phi_{2}^{2}s^{2}y^{2}+\frac{1}{4}\phi_{2}^{2}(y^{k}v_{k})^{2}+\frac{1}{4}(\psi^{\prime})^{2}v^{2}(y^{k}\hat{n}_{k})^{2}$
(154)
$\displaystyle\qquad+\frac{1}{4}\Big{[}3(\psi^{\prime})^{2}-2\psi^{\prime}\phi_{2}\Big{]}s^{2}(y^{k}\hat{n}_{k})^{2}+\frac{1}{2}\Big{[}\psi^{\prime}\phi_{2}-\phi_{2}^{2}\Big{]}s(y^{k}v_{k})(y^{k}\hat{n}_{k})\>.$

Since $\phi_{2}^{2}-\psi^{\prime}\phi_{2}\geq 0$ (Lemma 12 (b)), we apply the
inequality $-2ab\leq a^{2}+b^{2}$ and obtain

$\frac{1}{2}\Big{[}\psi^{\prime}\phi_{2}-\phi_{2}^{2}\Big{]}s(y^{k}v_{k})(y^{k}\hat{n}_{k})\leq\frac{1}{4}\Big{[}\phi_{2}^{2}-\psi^{\prime}\phi_{2}\Big{]}(y^{k}v_{k})^{2}+\frac{1}{4}\Big{[}\phi_{2}^{2}-\psi^{\prime}\phi_{2}\Big{]}s^{2}(y^{k}\hat{n}_{k})^{2}\>.$
(155)

Therefore,

$\displaystyle\|\nabla_{\boldsymbol{\mathcal{I}}}(v^{i}\mathsf{k}_{ij}\mathcal{D})\boldsymbol{y}\|^{2}\leq$
$\displaystyle\,\frac{1}{4}\phi_{2}^{2}s^{2}y^{2}+\frac{1}{4}\Big{[}2\phi_{2}^{2}-\psi^{\prime}\phi_{2}\Big{]}(y^{k}v_{k})^{2}+\frac{1}{4}(\psi^{\prime})^{2}v^{2}(y^{k}\hat{n}_{k})^{2}$
(156)
$\displaystyle+\frac{1}{4}\Big{[}3(\psi^{\prime})^{2}-3\psi^{\prime}\phi_{2}+\phi_{2}^{2}\Big{]}s^{2}(y^{k}\hat{n}_{k})^{2}\>.$

Further, since $\phi_{2}^{2}\geq 0$, $(\psi^{\prime})^{2}\geq 0$,
$2\phi_{2}^{2}-\psi^{\prime}\phi_{2}\geq 0$ (Lemma 12 (b)), and
$3(\psi^{\prime})^{2}-3\psi^{\prime}\phi_{2}+\phi_{2}^{2}\geq 0$ (Lemma 12
(c)), we can take the upper bounds $s^{2}\leq v^{2}$,
$(y^{k}v_{k})^{2}\leq(vy)^{2}$, and $(y^{k}\hat{n}_{k})^{2}\leq y^{2}$ to
obtain

$\|\nabla_{\boldsymbol{\mathcal{I}}}(v^{i}\mathsf{k}_{ij}\mathcal{D})\boldsymbol{y}\|^{2}\leq\,\Big{[}\phi_{2}^{2}-\psi^{\prime}\phi_{2}+(\psi^{\prime})^{2}\Big{]}{(vy)^{2}}\>.$
(157)

It follows from Lemma 12 (d) that
$\partial_{h}(\phi_{2}^{2}-\psi^{\prime}\phi_{2}+(\psi^{\prime})^{2})>0$,
which implies
$\max_{h\in[0,1]}\big{[}\phi_{2}^{2}-\psi^{\prime}\phi_{2}+(\psi^{\prime})^{2}\big{]}=\phi_{2}^{2}(1)-\psi^{\prime}(1)\phi_{2}(1)+(\psi^{\prime}(1))^{2}=4$.
Thus,

$\|\nabla_{\boldsymbol{\mathcal{I}}}(v^{i}\mathsf{k}_{ij}\mathcal{D})\boldsymbol{y}\|^{2}\leq\,4{(vy)^{2}}\>,\quad\forall\boldsymbol{y}=\big{(}y^{1},y^{2},y^{3}\big{)}^{\intercal},$
(158)

which proves the claim. ∎

## References

  * [1] D. Mihalas, B. W. Mihalas, Foundations of radiation hydrodynamics, Dover (New York), 1999.
  * [2] M. Rampp, H. T. Janka, Radiation hydrodynamics with neutrinos. Variable Eddington factor method for core-collapse supernova simulations, Astronomy & Astrophysics 396 (2002) 361–392.  doi:10.1051/0004-6361:20021398.
  * [3] O. Just, M. Obergaulinger, H.-T. Janka, A new multidimensional, energy-dependent two-moment transport code for neutrino-hydrodynamics, MNRAS 453 (2015) 3386–3413.
  * [4] M. A. Skinner, J. C. Dolence, A. Burrows, D. Radice, D. Vartanyan, FORNAX: A Flexible Code for Multiphysics Astrophysical Simulations, ApJS 241 (2019) 7.
  * [5] S. W. Bruenn, J. M. Blondin, W. R. Hix, E. J. Lentz, O. E. B. Messer, A. Mezzacappa, E. Endeve, J. A. Harris, P. Marronetti, R. D. Budiardja, M. A. Chertkow, C.-T. Lee, CHIMERA: A Massively Parallel Code for Core-collapse Supernova Simulations, ApJS 248 (1) (2020) 11.  doi:10.3847/1538-4365/ab7aff.
  * [6] A. Mezzacappa, E. Endeve, O. E. B. Messer, S. W. Bruenn, Physical, numerical, and computational challenges of modeling neutrino transport in core-collapse supernovae, Living Reviews in Computational Astrophysics 6 (1) (2020) 4.  doi:10.1007/s41115-020-00010-8.
  * [7] O. Just, A. Bauswein, R. Ardevol Pulpillo, S. Goriely, H. T. Janka, Comprehensive nucleosynthesis analysis for ejecta of compact binary mergers, MNRAS 448 (1) (2015) 541–567.  doi:10.1093/mnras/stv009.
  * [8] F. Foucart, Neutrino transport in general relativistic neutron star merger simulations, Living Reviews in Computational Astrophysics 9 (1) (2023) 1.  doi:10.1007/s41115-023-00016-y.
  * [9] S. Chapman, T. Cowling, The Mathematical Theory of Non-uniform Gases, Cambridge Mathematical Library, Cambridge University Press, 1970.
  * [10] J. I. Castor, Radiative Transfer in Spherically Symmetric Flows, Astrophysical Journal 178 (1972) 779–792.  doi:10.1086/151834.
  * [11] J. R. Buchler, Radiation hydrodynamics in the fluid frame., JQSRT 22 (1979) 293–300.  doi:10.1016/0022-4073(79)90119-5.
  * [12] J. R. Buchler, Radiation transfer in the fluid frame., JQSRT 30 (1983) 395–407.  doi:10.1016/0022-4073(83)90102-4.
  * [13] D. Mihalas, R. I. Klein, On the solution of the time-dependent inertial-frame equation of radiative transfer in moving media to O(v/c)., Journal of Computational Physics 46 (1982) 97–137.  doi:10.1016/0021-9991(82)90007-9.
  * [14] H. Nagakura, K. Sumiyoshi, S. Yamada, Three-dimensional Boltzmann Hydro Code for Core Collapse in Massive Stars. I. Special Relativistic Treatments, Astrophysical Journal Supplement Series 214 (2) (2014) 16.  doi:10.1088/0067-0049/214/2/16.
  * [15] A. Burrows, S. Reddy, T. A. Thompson, Neutrino opacities in nuclear matter, Nuclear Physics A 777 (2006) 356–394.  doi:10.1016/j.nuclphysa.2004.06.012.
  * [16] H.-T. Janka, K. Langanke, A. Marek, G. Martínez-Pinedo, B. Müller, Theory of core-collapse supernovae, Physics Reports 442 (1) (2007) 38–74.  doi:https://doi.org/10.1016/j.physrep.2007.02.002.
  * [17] H.-T. Janka, Explosion Mechanisms of Core-Collapse Supernovae, Annual Review of Nuclear and Particle Science 62 (1) (2012) 407–451.  doi:10.1146/annurev-nucl-102711-094901.
  * [18] T. Fischer, G. Guo, K. Langanke, G. Martinez-Pinedo, Y.-Z. Qian, M.-R. Wu, Neutrinos and nucleosynthesis of elements, arXiv e-prints (2023) arXiv:2308.03962arXiv:2308.03962, doi:10.48550/arXiv.2308.03962.
  * [19] G. N. Minerbo, Maximum entropy Eddington factors., JQSRT 20 (1978) 541–545.
  * [20] J. Cernohorsky, S. A. Bludman, Maximum entropy distribution and closure for bose-einstein and Fermi-Dirac radiation transport, Astrophysical Journal 433 (1) (1994) 450 – 455.
  * [21] M. Shibata, K. Kiuchi, Y. Sekiguchi, Y. Suwa, Truncated Moment Formalism for Radiation Hydrodynamics in Numerical Relativity, Progress of Theoretical Physics 125 (2011) 1255–1287.
  * [22] C. Y. Cardall, E. Endeve, A. Mezzacappa, Conservative 3+1 general relativistic Boltzmann equation, Physical Review D 88 (2013) 023011.
  * [23] R. B. Lowrie, D. Mihalas, J. Morel, Comoving-frame radiation transport for nonrelativistic fluid velocities, JQSRT 69 (2001) 291–304.
  * [24] N. Vaytet, E. Audit, B. Dubroca, F. Delahaye, A numerical model for multigroup radiation hydrodynamics, JQSRT 112 (2011) 1323–1335.
  * [25] B. Cockburn, C.-W. Shu, Runge-Kutta Discontinuous Galerkin Methods for Convection-Dominated Problems, Journal of Scientific Computing 16 (2001) 173–261.
  * [26] E. W. Larsen, J. E. Morel, Asymptotic Solutions of Numerical Transport Problems in Optically Thick, Diffusive Regimes II, Journal of Computational Physics 83 (1989) 212–236.
  * [27] M. L. Adams, Discontinuous finite element transport solutions in thick diffusive problems, Nuclear science and engineering 137 (3) (2001) 298–333.
  * [28] J.-L. Guermond, G. Kanschat, Asymptotic Analysis of Upwind Discontinuous Galerkin Approximation of the Radiative Transport Equation in the Diffusive Limit, SIAM J. Numer. Anal. 48 (2010) 53–78.
  * [29] E. Audit, P. Charrier, J. P. Chièze, B. Dubroca, A radiation-hydrodynamics scheme valid from the transport to the diffusion limit, arXiv e-prints (2002) astro–ph/0206281doi:10.48550/arXiv.astro-ph/0206281.
  * [30] B. Ayuso, J. A. Carrillo, S. C.-W., Discontinuous Galerkin Methods for the One-Dimensional Vlasov–Poisson System, Kinetic and Related Models 4 (4) (2011) 955–989.
  * [31] Y. Cheng, I. M. Gamba, P. J. Morrison, Study of conservation and recurrence of Runge–Kutta discontinuous Galerkin schemes for Vlasov–Poisson systems, Journal of Scientific Computing 56 (2013) 319–349.
  * [32] U. Ascher, S. Ruuth, R. Spiteri, Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations, Applied Numerical Mathematics 25 (1997) 151–167.
  * [33] L. Pareschi, G. Russo, Implicit-Explicit Runge-Kutta Schemes and Application to Hyperbolic Systems with Relaxation, Journal of Scientific Computing 25 (2005) 129–155.
  * [34] R. Chu, E. Endeve, C. Hauck, A. Mezzacappa, Realizability-preserving dg-imex method for the two-moment model of fermion transport, Journal of Computational Physics 389 (2019) 62–93.
  * [35] M. P. Laiu, E. Endeve, R. Chu, J. A. Harris, O. E. B. Messer, A DG-IMEX Method for Two-moment Neutrino Transport: Nonlinear Solvers for Neutrino-Matter Coupling, Astrophysical Journal Supplement Series 253 (2) (2021) 52.  doi:10.3847/1538-4365/abe2a8.
  * [36] C. D. Levermore, Moment closure hierarchies for kinetic theories, Journal of Statistical Physics 83 (1996) 1021–1065.
  * [37] G. W. Alldredge, M. Frank, C. D. Hauck, A regularized entropy-based moment method for kinetic equations, SIAM Journal on Applied Mathematics 79 (5) (2019) 1627–1653.
  * [38] E. Olbrant, C. D. Hauck, M. Frank, A realizability-preserving discontinuous galerkin method for the m1 model of radiative transfer, Journal of Computational Physics 231 (17) (2012) 5612–5639.
  * [39] C. D. Hauck, High-order entropy-based closures for linear transport in slab geometry, Communications in Mathematical Sciences 9 (1) (2011) 187–205.
  * [40] L. F. Richardson, Ix. the approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam, Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 210 (459-470) (1911) 307–357.
  * [41] Y. Saad, Iterative methods for sparse linear systems, SIAM, 2003.
  * [42] M. Liebendörfer, O. E. B. Messer, A. Mezzacappa, S. W. Bruenn, C. Y. Cardall, F. K. Thielemann, A Finite Difference Representation of Neutrino Radiation Hydrodynamics in Spherically Symmetric General Relativistic Spacetime, Astrophysical Journal Supplement Series 150 (1) (2004) 263–316.  doi:10.1086/380191.
  * [43] B. Müller, H.-T. Janka, H. Dimmelmeier, A New Multi-dimensional General Relativistic Neutrino Hydrodynamic Code for Core-collapse Supernovae. I. Method and Code Tests in Spherical Symmetry, Astrophysical Journal Supplement Series 189 (1) (2010) 104–133.  doi:10.1088/0067-0049/189/1/104.
  * [44] R. W. Lindquist, Relativistic transport theory, Annals of Physics 37 (1966) 487–518.
  * [45] E. Endeve, C. Y. Cardall, A. Mezzacappa, Conservative Moment Equations for Neutrino Radiation Transport with Limited Relativity, ArXiv e-prints (Dec. 2012).  arXiv:1212.4064.
  * [46] C. D. Levermore, Relating Eddington factors to flux limiters., JQSRT 31 (1984) 149–160.
  * [47] S. C. Noble, C. F. Gammie, J. C. McKinney, L. Del Zanna, Primitive variable solvers for conservative general relativistic magnetohydrodynamics, Astrophysical Journal 641 (2006) 626–637.
  * [48] R. E. Curto, L. A. Fialkow, Recursiveness, positivity, and truncated moment problems, Houston Journal of Mathematics 17 (4) (1991).
  * [49] J. S. Hesthaven, T. Warburton, Nodal discontinuous Galerkin methods: Algorithms, analysis and applications, Springer, 2008.
  * [50] M. Abramowitz, I. A. Stegun, R. H. Romer, Handbook of mathematical functions with formulas, graphs, and mathematical tables, American Journal of Physics 56 (10) (1988) 958–958.  arXiv:https://doi.org/10.1119/1.15378, doi:10.1119/1.15378.
URL https://doi.org/10.1119/1.15378

  * [51] E. Hairer, S. Nørsett, G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, Springer Series in Computational Mathematics, Springer Berlin Heidelberg, 1993.
  * [52] X. Zhang, C.-W. Shu, On maximum-principle-satisfying high order schemes for scalar conservation laws, Journal of Computational Physics 229 (2010) 3091–3120.
  * [53] X. Zhang, C.-W. Shu, On positivity preserving high order discontinuous galerkin schemes for compressible euler equations on rectangular meshes, Journal of Computational Physics 229 (2010) 8918–8934.
  * [54] A. Dubey, K. Weide, J. O’Neal, A. Dhruv, S. Couch, J. A. Harris, T. Klosterman, R. Jain, J. Rudi, B. Messer, M. Pajkos, J. Carlson, R. Chu, M. Wahib, S. Chawdhary, P. M. Ricker, D. Lee, K. Antypas, K. M. Riley, C. Daley, M. Ganapathy, F. X. Timmes, D. M. Townsley, M. Vanella, J. Bachan, P. M. Rich, S. Kumar, E. Endeve, W. R. Hix, A. Mezzacappa, T. Papatheodore, Flash-X: A multiphysics simulation software instrument, SoftwareX 19 (2022) 101168.  doi:10.1016/j.softx.2022.101168.
  * [55] M. P. Laiu, J. A. Harris, R. Chu, E. Endeve, thornado-transport: Anderson- and GPU-accelerated nonlinear solvers for neutrino-matter coupling, Journal of Physics: Conference Series 1623 (2020) 012013.  doi:10.1088/1742-6596/1623/1/012013.
  * [56] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, D. Sorensen, LAPACK Users’ Guide, 3rd Edition, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999.
  * [57] S. Tomov, J. Dongarra, M. Baboulin, Towards dense linear algebra for hybrid GPU accelerated manycore systems, Parallel Computing 36 (5-6) (2010) 232–240.  doi:10.1016/j.parco.2009.12.005.
  * [58] C.-W. Shu, O. S., Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes, Journal of Computational Physics 77 (1988) 439–471.
  * [59] D. G. Anderson, Iterative procedures for nonlinear integral equations, J. ACM 12 (4) (1965) 547–560.  doi:10.1145/321296.321305.
URL http://doi.acm.org/10.1145/321296.321305

  * [60] H. Walker, P. Ni, Anderson acceleration for fixed-point iterations, SIAM Journal on Numerical Analysis 49 (4) (2011) 1715–1735.  arXiv:https://doi.org/10.1137/10078356X, doi:10.1137/10078356X.
URL https://doi.org/10.1137/10078356X

  * [61] D. Cavaglieri, T. Bewley, Low-storage implicit/explicit Runge–Kutta schemes for the simulation of stiff high-dimensional ODE systems, Journal of Computational Physics 286 (2015) 172 – 193.