# High order conservative LDG-IMEX methods for the degenerate nonlinear non-
equilibrium radiation diffusion problems***The first and the fourth authors
are supported by National Key R&D Program of China No. 2022YFA1004500, NSFC
No. 11971025 and No. 92270112, NSF of Fujian Province No. 2023J02003. The
second author is partially supported by NSFC No. 12031013, Shanghai pilot
innovation project No. 21JC1403500 and the Strategic Priority Research Program
of Chinese Academy of Sciences Grant No. XDA25010401. The third author is
supported by NSFC No. 12071214.

Shaoqin Zheng<EMAIL_ADDRESS>Min Tang<EMAIL_ADDRESS>Qiang
Zhang<EMAIL_ADDRESS>Tao Xiong<EMAIL_ADDRESS>
###### Abstract

In this paper, we develop a class of high-order conservative methods for
simulating non-equilibrium radiation diffusion problems. Numerically, this
system poses significant challenges due to strong nonlinearity within the
stiff source terms and the degeneracy of nonlinear diffusion terms. Explicit
methods require impractically small time steps, while implicit methods, which
offer stability, come with the challenge to guarantee the convergence of
nonlinear iterative solvers. To overcome these challenges, we propose a
predictor-corrector approach and design proper implicit-explicit time
discretizations. In the predictor step, the system is reformulated into a
nonconservative form and linear diffusion terms are introduced as a
penalization to mitigate strong nonlinearities. We then employ a Picard
iteration to secure convergence in handling the nonlinear aspects. The
corrector step guarantees the conservation of total energy, which is vital for
accurately simulating the speeds of propagating sharp fronts in this system.

For spatial approximations, we utilize local discontinuous Galerkin finite
element methods, coupled with positive-preserving and TVB limiters. We
validate the orders of accuracy, conservation properties, and suitability of
using large time steps for our proposed methods, through numerical experiments
conducted on one- and two-dimensional spatial problems. In both homogeneous
and heterogeneous non-equilibrium radiation diffusion problems, we attain a
time stability condition comparable to that of a fully implicit time
discretization. Such an approach is also applicable to many other reaction-
diffusion systems.

###### keywords:

non-equilibrium radiation diffusion , predictor-corrector procedure ,
conservative , high order , local discontinuous Galerkin method , IMEX

[1]organization=School of Mathematical Sciences, addressline=Xiamen
University, city=Xiamen, postcode=Fujian 361005, country=PR China
[2]organization=Institute of Natural Sciences, addressline=Shanghai Jiao Tong
University, city=Shanghai, postcode=200240, country=PR China
[3]organization=Department of Mathematics, addressline=Nanjing University,
city=Nanjing, postcode=210093, country=PR China [4]organization=School of
Mathematical Sciences, Fujian Provincial Key Laboratory of Mathematical
Modeling and High-Performance Scientific Computing, addressline=Xiamen
University, city=Xiamen, postcode=Fujian 361005, country=PR China

##  1 Introduction

In scenarios where radiation interacts dynamically with a material, but does
not fully reach thermodynamic equilibrium, a commonly employed description
involves a system of coupled time-dependent nonlinear diffusion equations.
This system is known for its strong nonlinearity and close coupling, and it
finds extensive applications across various fields, including inertial
confinement fusion [47], astrophysics [8], and Z-pinch experiments [38]. In
this work, we specifically explore a two-temperature (2T) model. This model
comprises a radiation diffusion equation within the framework of a gray
approximation and a material energy balance equation [8, 36, 50]. The 2T model
writes:

$\begin{cases}\dfrac{\partial E}{\partial t}-\nabla\cdot\left(D_{r}\nabla
E\right)=\sigma\left(T^{4}-E\right),\\\\[5.0pt] \dfrac{\partial T}{\partial
t}-\nabla\cdot\left(D_{t}\nabla
T\right)=\sigma\left(E-T^{4}\right),\end{cases}$ (1.1)

where $E({\bf x},t)$ is the radiation energy and $T({\bf x},t)$ is the
material temperature. The energy exchange between materials and photons is
controlled by the photon absorption cross-section $\sigma$, which is modeled
by

$\sigma(T)=\frac{z({\bf x})^{3}}{T^{3}},$ (1.2)

where $z({\bf x})$ is a spatial dependent material coefficient that represents
an atomic mass number. The following flux-limited energy radiation diffusion
coefficient $D_{r}$ is widely used [8, 23]:

$D_{r}=\frac{1}{3\sigma+\frac{|\nabla E|}{E}}.$ (1.3)

Moreover, the material conduction coefficient has the following form [43]

$D_{t}=\kappa T^{\frac{5}{2}},$ (1.4)

where $\kappa$ is a constant.

Adding the two equations of (1.1) together, one gets

$\dfrac{\partial(E+T)}{\partial t}=\nabla\cdot\left(D_{r}\nabla
E\right)+\nabla\cdot\left(D_{t}\nabla T\right),$ (1.5)

which indicates that the system (1.1) conserves the energy $E+T$ over the
space. Energy conservation is not only physically important, but also crucial
to capture a useful simplification of the 2T model. As $\sigma$ approaches
$\infty$, the radiation energy tends to approximate the thermal equilibrium,
namely $E\approx T^{4}$, (1.5) would lead to the gray radiative diffusion
equation [26]

$\dfrac{\partial(T^{4}+T)}{\partial t}=\nabla\cdot\left(D_{r}\nabla
T^{4}\right)+\nabla\cdot\left(D_{t}\nabla T\right).$ (1.6)

The gray radiation equation is an important approximation widely employed for
studying diverse radiative heat transfer phenomena, including radiative
transfer in stellar atmospheres.

Non-equilibrium radiation diffusion equations have attracted extensive
research efforts. Marshak, for instance, developed a time-dependent radiative
transfer model to investigate the impact of radiation on shock wave behavior
[30]. Analytical solutions have also been provided for specific Marshak wave
problems in previous studies [37, 6]. However, solving the system (1.1)
numerically presents substantial challenges, mainly due to the following
several reasons:

  * 1.

The source terms might be very stiff. Specifically, the coefficient $\sigma$,
as defined in (1.2), tends to be very large for large values of $z$ or small
values of $T$. This leads to pronounced energy exchanges and a close coupling
between $E$ and $T^{4}$. Dealing with such stiff source terms when $\sigma$ is
large requires implicit time discretizations, resulting in a highly nonlinear
system with multiscale coefficients.

  * 2.

The diffusion terms $\nabla\cdot\left(D_{r}\nabla E\right)$ and
$\nabla\cdot\left(D_{t}\nabla T\right)$ are nonlinear and degenerate. The
coefficients $D_{r}$ and $D_{t}$ defined in (1.3)-(1.4) nonlinearly depend on
$T$, and they may vary by several orders of magnitude when the temperature $T$
varies over a wide range. Moreover, when $T$ is small, both $D_{r}$ and
$D_{t}$ are close to $0$, which result in degenerate diffusion. Consequently,
solutions exhibit sharp fronts that propagate at finite speeds, akin to shock
wave solutions in hyperbolic systems [46].

  * 3.

A challenging condition number for a nonlinear iteration. When the temperature
$T$ is low, energy $E$ and temperature $T$ operate at significantly different
scales, with $E$ roughly proportionate to $T^{4}$. Consequently, using a fully
implicit time discretization for both $E$ and $T$, combined with Newton
iteration, can lead to a coefficient matrix with a challenging condition
number. It is essential to devise suitable iterative solvers and efficient
preconditioners to tackle this challenge effectively.

When solving (1.1) with explicit schemes, it is relatively straightforward to
maintain the conservation of total energy $E+T$, however a significant
drawback is the requirement of very small time steps. A parabolic time step
$\Delta t=\mathcal{O}(h^{2})$ with $\Delta t$ being the time step and $h$
being the mesh size is needed due to the nonlinear diffusion terms, and the
stiff source terms lead to a severe time step constraint as well [55]. Fully
implicit schemes offer the advantage of using large time steps. However, the
presence of strong nonlinearity and multiscale coefficients, when employing a
large time step, the solutions from the previous time step may not serve as a
suitable initial guess. As a result, iterative solvers used in fully implicit
methods may struggle to converge effectively [1, 20]. In past years,
significant research efforts have been dedicated to developing efficient
iterative solvers with effective preconditioners. One notable approach is the
Jacobian-free Newton-Krylov method [23, 24], which combines Newtonian external
iterations and conjugate gradient-like (Krylov) internal iterations, resulting
in superlinear convergence without the need to form Jacobian matrices.
Additionally, a physical-based preconditioning Newton-Krylov method was
explored in [33], and an operator-split preconditioner was investigated in
[31]. Various other efficient time discretization methods have also been
examined in [21, 29, 9, 32, 35, 22] and references therein. Two semi-implicit
schemes allowing for large time steps for the gray radiation diffusion
equation (1.6) have been developed in [46]. Furthermore, radiation diffusion
models are frequently coupled with hydrodynamic equations, and large
deformations of complex fluid flows lead to mesh distortions. Consequently,
for spatial discretization, there have been some works built upon distorted
meshes [19, 40, 56, 62], as well as moving meshes [54, 55].

Schemes that efficiently handle non-equilibrium radiation diffusion equations,
while being robust and capable of accommodating large time steps for accurate
sharp front capturing, remain relatively rare in the existing literature. To
design such a scheme, several requirements have to be met:

  * I

Preservation of the energy equation (1.5) at the discrete level. This
indicates that the source terms in the two equations in (1.1) are in balance
with each other. As far as (1.5) is preserved at the discrete level, when
$\sigma$ is large and the system reaches the thermodynamic equilibrium
$E=T^{4}$, the scheme’s accuracy can be controlled by the gray radiation
equation (1.6);

  * II

Conservation of energy. Because of the presence of the nonlinear degenerate
diffusion terms, the solution exhibits hyperbolic properties at the front. As
is well-established in the study of hyperbolic simulations, preserving the
conservation properties of the scheme’s discretization is of paramount
importance to accurately capture finite-speed front propagations.

  * III

Positivity Preserving. The system is only physically meaningful for positive
radiation energy $E$ and positive material temperature $T$. The radiation
temperature $T_{r}=E^{1/4}$ and all coefficients in (1.1) depend on $T$.
Hence, it would be crucial to maintain positivity of $E$ and $T$ in order to
ensure the robustness of the scheme and obtain physically meaningful
solutions.

  * IV

Efficient iterative solvers which can guarantee convergence and the providing
of a good initial guess. Fully implicit time discretizations coupled with
proper spatial discretizaitons may possibly meet the above mentioned
requirements, but due to strong nonlinearities and multiscale variations of
coefficients, one has to provide a good initial guess and a suitable
preconditioner in order to make an iterative solver converging properly.

We note that radiation diffusion equations belongs to reaction-diffusion
equations which appear in many other fields, such as in material sciences
[34], chemical reactions [14], ecology and biology systems [18], and many
numerical approaches have been developed for reaction-diffusion equations,
e.g. [63, 61, 27] and reference therein. However, for the radiation diffusion
equations (1.1), the main challenges are the degeneracy in the nonlinear
diffusion operators and stiffness in the source terms. Both are caused by the
low material temperature $T$, making degeneracy, stiffness and nonlinearity
strongly coupled, while standard reaction diffusion equations usually involve
only linear diffusion but various nonlinear reactions. The design of an
efficient iterative solver for radiation diffusion equations is very
challenging. In this paper, to address such difficulties, we present a novel
and efficient method for (1.1) utilizing implicit-explicit (IMEX) time
discretizations, designed to satisfy all four of the aforementioned criteria.
Some techniques we develop can also be applied to other reaction diffusion
systems. We introduce a predictor-corrector procedure as our proposed
approach. In the predictor step, we multiply both sides of the second equation
in (1.1) with $4T^{3}$, introduce $B=T^{4}$, and obtain a simpler system for
$B$ and $E$ [16]. This new system does not meet requirements I, II, and III.
However, we can design an efficient iterative solver capable of ensuring
convergence even for large time steps. The predictor step can provide a good
initial guess. Then, in the corrector step we solve the energy conservative
equation (1.5), and substitute $E+T$ into the first equation of (1.1) to
further update $E$. Requirements I and II are satisfied in the corrector step
and the bad condition number for solving $E$ and $T$ together is avoided.
Finally, requirement III can be attained through the utilization of suitable
spatial discretizations.

For requirement IV, several strategies are employed to alleviate the
challenges associated with solving a highly nonlinear system. First of all,
implicit treatments of nonlinear diffusion terms are avoided by using the idea
in [51, 49] from adding the same linear diffusion terms on both sides of those
equations, which however are discretized explicitly and implicitly
respectively. For both the predictor and corrector steps, only a mildly
nonlinear system needs to be solved, and a simple Picard iteration is adopted.
Secondly, a local discontinuous Galerkin (LDG) finite element method is
employed for spatial discretiztions, due to a shock wave-like behavior in the
solutions of radiation [30]. The LDG method was introduced by Cockburn and Shu
in [12, 13] for convection-diffusion problems, which is well-suited for h-p
adaptivity and is also very good at shock or sharp gradient capturing. For a
first order in space LDG method, the iterative matrix for $E$ and $B$ can be
shown to be an M-matrix [16]. This can guanrantee the convergence of the
Picard iteration in the predictor step. Finally, for second and third orders
of discretizations, appropriate spatial limiters are applied to control
numerical oscillations for sharp gradient solutions and preserve the
positivity of radiation energy $E$ and material temperature $T$. After spatial
discretizations, the nonlinearity iteration can be written in the form of a
linear system, with nonlinearities mainly appear in the diagonal part of the
resulting system, a fast convergence can be obtained.

The rest of this paper is organized as follows. In section 2, an efficient
IMEX time discretization is introduced via a predictor-corrector procedure,
followed by an LDG space discretization in section 3. In section 4, one and
two spatial dimensional numerical examples are performed to verify the high
order accuracy, conservation and large time step conditions, and good
performances for capturing sharp fronts of the radiation energy in both
homogeneous and heterogeneous media. Finally, a conclusion is drawn in section
5.

##  2 Time discretization

In this section, we present an efficient IMEX time discretization for (1.1),
employing a predictor-corrector procedure. The predictor step, though non-
conservative, ensures fast and robust nonlinear convergence. In contrast, the
corrector step maintains conservation, which is crucial for accurately
propagating sharp fronts. This novel approach significantly mitigates
nonlinearity as compared to a fully implicit time discretization, while
maintaining a generous time step stability condition and ensuring fast and
robust convergence. We begin with a first-order IMEX scheme, which serves as a
general framework that can be readily extended to higher orders by integrating
a multistage Runge-Kutta (RK) or multistep IMEX time-marching approach. Here,
we adopt an IMEX RK scheme. The IMEX time discretization can be combined with
any front-capturing spatial discretization. In this section, we keep space
continuous, with a detailed description of the spatial discretization provided
in the next section.

###  2.1 First order IMEX scheme

Let the discrete time steps be $t_{n}\,(n=0,1,\ldots)$ and the time step size
be $\Delta t=t_{n+1}-t_{n}$. Instead of solving (1.1) directly, we propose a
predictor-corrector procedure:

  * 1.

The predictor step:
To mitigate the stiffness of the source terms caused by $\sigma=z^{3}/T^{3}$,
as in [16], we multiply both sides of the second equation in (1.1) by
$4T^{3}$. (1.1) can be rewritten as:

$\begin{cases}\dfrac{\partial E}{\partial t}-\nabla\cdot\left(D_{r}\nabla
E\right)=\sigma\left(B-E\right),\\\\[10.0pt] \dfrac{\partial B}{\partial
t}-4\kappa\left[\nabla\cdot\left(T^{\frac{11}{2}}\nabla
T\right)-3T^{\frac{9}{2}}|\nabla
T|^{2}\right]=4z^{3}\left(E-B\right),\\\\[10.0pt] T=B^{1/4}.\end{cases}$ (2.1)

As discussed in the introduction, $T^{4}$ in the source term of (1.1) needs to
be treated implicitly to ensure thermal dynamic equilibrium ($E=T^{4}$) in the
stiff regime when $T$ is small. To achieve this without resorting to nonlinear
iterations for a given $\sigma$, we introduce a new variable, $B=T^{4}$, and
update a system for both $E$ and $B$ in this predictor step. Utilizing the
equation for $B$ not only avoids the stiffness associated with $\sigma$ in the
second equation but also maintains balance between $E$ and $B$, ensuring they
remain at the same scale when $T$ is small, thus preventing unfavorable
condition numbers. Furthermore, in this equation, we continue to use the
material temperature $T$ rather than $B$ for the diffusion terms, thereby
preventing negative powers of $B$ which can lead to stiffness when $B=T^{4}$
is small. Once we have determined $B$, we can calculate $T$ using the
relationship $B=T^{4}$. It is noteworthy that throughout this process and in
the following, we consistently use $E$ and $T$ as our input and output
variables.

To update (2.1) for $E$ and $B$, the diffusion terms remain nonlinear. An
implicit treatment would lead to a complex nonlinear system. Following the
approach presented in [51, 49], we introduce two linear diffusion terms
$\alpha_{0}\Delta E$ and $\beta_{0}\Delta B$ with constant coefficients
$\alpha_{0}$, $\beta_{0}$ on both sides of (2.1) for $E$ and $B$,
respectively. Starting from $E^{n}$ and $T^{n}$, a first-order IMEX scheme for
updating predicted solutions $\tilde{E}^{n+1}$ and $\tilde{T}^{n+1}$ is
defined as follows:

$\begin{cases}\frac{\tilde{E}^{n+1}-E^{n}}{\Delta
t}-\nabla\cdot\left(D^{n}_{r}\nabla E^{n}\right)+\alpha_{0}\,\Delta
E^{n}=\alpha_{0}\,\Delta\tilde{E}^{n+1}+\tilde{\sigma}^{n+1}\left(\tilde{B}^{n+1}-\tilde{E}^{n+1}\right),\\\\[10.0pt]
\frac{\tilde{B}^{n+1}-B^{n}}{\Delta t}-\mathcal{H}^{n}+\beta_{0}\,\Delta
B^{n}=\beta_{0}\,\Delta\tilde{B}^{n+1}+4z^{3}\left(\tilde{E}^{n+1}-\tilde{B}^{n+1}\right),\\\\[10.0pt]
\tilde{T}^{n+1}=(\tilde{B}^{n+1})^{1/4},\end{cases}$ (2.2)

where

$\mathcal{H}=4\kappa\left[\nabla\cdot\left(T^{\frac{11}{2}}\nabla
T\right)-3T^{\frac{9}{2}}|\nabla T|^{2}\right].$ (2.3)

Here, the superscript $n$ or $n+1$ denotes the corresponding values at time
step $t_{n}$ or $t_{n+1}$ respectively, similarly in the following. As we can
see, in each equation, two added identical terms are discretized differently,
one explicitly and the other implicitly. This approach allows us to achieve
time stability close to that of an implicit scheme while only needing to solve
linearly implicit diffusion operators [49, 51]. From the above, we observe
that by introducing $B$ and adding the linear diffusion terms, when we solve
(2.2) with an IMEX method, the only nonlinearity arises from the source term
$\sigma(E-B)$ due to the coefficient $\sigma$ appearing in the first equation.
If $\sigma$ is constant or a predefined function, (2.2) becomes a simple
linear system. Additionally, if $\kappa=0$ for $\mathcal{H}$ in (2.3), (2.2)
is in a conservative form, and such a 2T model can be efficiently solved.

  * 2.

The corrector step:
However, if $\kappa$ is not zero, the term $\mathcal{H}$ in (2.3) within the
equation of $B$ is in a non-conservative form. This non-conservative form can
result in incorrect sharp front propagation, as will be demonstrated in our
numerical examples. To address this issue, a corrector step is required. To
ensure the conservation of total energy, as shown in (1.5), we begin by
replacing the second equation in the original equation (1.1) with (1.5),
resulting in the following:

$\begin{cases}\dfrac{\partial E}{\partial t}-\nabla\cdot\left(D_{r}\nabla
E\right)=\sigma\left(T^{4}-E\right),\\\\[10.0pt] \dfrac{\partial Q}{\partial
t}=\nabla\cdot\left(D_{r}\nabla E\right)+\nabla\cdot\left(D_{t}\nabla
T\right),\\\\[10.0pt] T=Q-E.\end{cases}$ (2.4)

The equation for $Q$ is in a conservative form, allowing for easy conservation
of total energy across space. Following (2.2), we also introduce two linear
diffusion terms, $\alpha_{0}\Delta E$ and $\gamma_{0}\Delta Q$, on both sides
of (2.4). Using a first-order IMEX time discretization, we obtain:

$\begin{cases}\dfrac{E^{n+1}-E^{n}}{\Delta t}-\nabla\cdot\left(D^{n}_{r}\nabla
E^{n}\right)+\alpha_{0}\,\Delta E^{n}=\alpha_{0}\,\Delta
E^{n+1}+\tilde{\sigma}^{n+1}\left((T^{n+1})^{4}-E^{n+1}\right),\\\\[10.0pt]
\dfrac{Q^{n+1}-Q^{n}}{\Delta t}-\nabla\cdot\left(D^{n}_{r}\nabla
E^{n}\right)-\nabla\cdot\left(D^{n}_{t}\nabla T^{n}\right)+\gamma_{0}\,\Delta
Q^{n}=\gamma_{0}\,\Delta Q^{n+1},\\\\[10.0pt]
T^{n+1}=Q^{n+1}-E^{n+1}.\end{cases}$ (2.5)

Solving for $E$ and $Q$ ensures good convergence and the conservation of $Q$.
As we can see, in (2.5), adding these linear diffusion terms is crucial to
simplify the nonlinear system. With this approach, we can first solve a linear
system for $Q$ from the second equation. Then, we can use $T=Q-E$ to replace
the $T^{4}$ term in the first equation, allowing us to solve this mildly
nonlinear equation to further update $E$. The nonlinearity is simply
$(Q-E)^{4}$ for $E$ and only appears in the diagonal part of the mass matrix
after spatial discretization.

######  Remark 1.

In the corrector step, we utilize the values obtained from the predictor step
to determine the coefficients $\sigma$ in the source term and to provide an
initial guess for the corrector step. One could contemplate simplifying the
scheme by eliminating the predictor step and directly employing a Picard
iteration in the corrector step to handle these coefficients. However, we
would mention that $\sigma=z^{3}/T^{3}$ and $T^{4}=(Q-E)^{4}$ appears as a
nonlinear term. In case of $T$ being small, with a stiff coefficient, the
nonlinear iteration for solving $E$ in (2.5) may not converge well, for
example, the stand Marshak wave problem in Example 4.4. Instead, in the
predictor step, except $\sigma$, others linearly depend on $E$ and $B$. Even
with a stiff $\sigma$, it appears in the diagonal part of the mass matrix, so
that a fast and robust convergence can be obtained. Starting from a good
initial guess provided by the predictor step, the corrector step can also
converge well.

###  2.2 High order IMEX scheme

The first-order IMEX scheme with a predictor-corrector procedure has offered a
highly efficient and versatile framework for solving the 2T model (1.1). To
attain high-order accuracy in time while preserving these desirable
properties, we employ a globally stiffly accurate IMEX RK time discretization
[7]. However, a multistep IMEX method can also be utilized [2].

We consider a system of additive ordinary differential equations:

$\frac{\mathrm{d}\bm{y}}{\mathrm{d}t}=L(t,\bm{y})+N(t,\bm{y}),\quad\bm{y}\left(t_{0}\right)=\bm{y}_{0},$
(2.6)

where $\bm{y}=(y_{1},y_{2},\ldots,y_{d})^{T}$. $L(t,\bm{y})$ and $N(t,\bm{y})$
are linear and nonlinear operators, respectively. $L(t,\bm{y})$ will be
discretized implicitly, while $N(t,\bm{y})$ will be discretized explicitly. An
$s$-stage IMEX RK time discretization can be represented by a double Butcher
tableau

$\begin{array}[]{c|c}\hat{c}&\hat{A}\\\ \hline\cr\vspace{-0.25cm}\hfil\\\
&\hat{b}^{T}\end{array},\ \ \ \qquad\begin{array}[]{c|c}{c}&{A}\\\
\hline\cr\vspace{-0.25cm}\hfil\\\ &{b^{T}}\end{array}.$

Here $A=\left(a_{ij}\right)$ and
$\hat{A}=\left(\hat{a}_{ij}\right)\in\mathbf{R}^{s\times s}$. $\hat{A}$ is a
strictly lower triangular matrix for explicit parts. For the implicit part,
$A$ can be taken as a lower triangular matrix with a nonzero diagonal to get
an efficient implementation, which is usually referred to as a diagonally
implicit RK (DIRK) scheme. The vectors are $b^{T}=(b_{1},b_{2},\ldots,b_{s})$,
$\hat{b}^{T}=(\hat{b}_{1},\hat{b}_{2},\ldots,\hat{b}_{s})$,
$c^{T}=(c_{1},c_{2},\ldots,c_{s})$, and
$\hat{c}^{T}=(\hat{c}_{1},\hat{c}_{2},\ldots,\hat{c}_{s})$, where
$c_{i}=\sum\limits_{j=1}^{i}a_{ij}$ and
$\hat{c}_{i}=\sum\limits_{j=1}^{i-1}\hat{a}_{ij}$. Denoting $t_{n}^{(j)}=$
$t_{n}+c_{j}\Delta t$, $\hat{t}_{n}^{(j)}=$ $t_{n}+\hat{c}_{j}\Delta t$, the
solution of (2.6) can be updated from time level $t^{n}$ to $t^{n+1}$ in the
following way:

$\begin{cases}\bm{Y}^{(1)}=\bm{y}_{n},\\\\[5.0pt]
\bm{Y}^{(i)}=\bm{y}_{n}+\Delta
t\sum\limits_{j=1}^{i-1}\hat{a}_{ij}N\left(\hat{t}_{n}^{(j)},\bm{Y}^{(j)}\right)+\Delta
t\sum\limits_{j=1}^{i}a_{ij}L\left(t_{n}^{(j)},\bm{Y}^{(j)}\right),\quad 2\leq
i\leq s,\\\\[5.0pt] \bm{y}_{n+1}=\bm{y}_{n}+\Delta
t\sum\limits_{i=1}^{s}\hat{b}_{i}N\left(\hat{t}_{n}^{(i)},\bm{Y}^{(i)}\right)+\Delta
t\sum\limits_{i=1}^{s}b_{i}L\left(t_{n}^{(i)},\bm{Y}^{(i)}\right).\end{cases}$
(2.7)

The IMEX RK scheme is called to be globally stiffly accurate, if the
coefficients satisfy [7] are required to satisfy:

$\hat{c}_{s}=c_{s}=1,\text{ and }a_{sj}=b_{j},\,\hat{a}_{sj}=\hat{b}_{j},\quad
j=1,2,\ldots,s.$

With such an IMEX scheme, the final updating of $\bm{y}_{n+1}$ in (2.7)
coincides with the last stage of updating $\bm{Y}^{(s)}$, so that we can take
$\bm{y}_{n+1}=\bm{Y}^{(s)}$ and avoid the last cumulative step. In A, Butcher
tableaux from first order to third order, which are adopted in this work, are
provided.

If we choose $L$ and $N$ in (2.6) based on the first order scheme (2.2) and
(2.5), with (2.7), the updating of the solutions at $t^{n+1}$ from $t^{n}$ can
be presented as follows:

$\displaystyle\frac{E^{n+1}-E^{n}}{\Delta t}$
$\displaystyle=\sum_{i=1}^{s-1}\hat{b}_{i}\left[\nabla\cdot(D^{(i)}_{r}\nabla
E^{(i)})-\alpha_{0}\Delta
E^{(i)}\right]+\sum_{i=1}^{s}b_{i}\left[\alpha_{0}\Delta
E^{(i)}+{\sigma}^{(i)}(B^{(i)}-E^{(i)})\right],$ (2.8a)
$\displaystyle\frac{Q^{n+1}-Q^{n}}{\Delta t}$
$\displaystyle=\sum_{i=1}^{s-1}\hat{b}_{i}\left[\nabla\cdot\left(D^{(i)}_{r}\nabla
E^{(i)}\right)+\nabla\cdot\left(D^{(i)}_{t}\nabla
T^{(i)}\right)-\gamma_{0}\Delta
Q^{(i)}\right]+\gamma_{0}\sum_{i=1}^{s}b_{i}\Delta Q^{(i)},$ (2.8b)
$\displaystyle T^{n+1}$ $\displaystyle=Q^{n+1}-E^{n+1}.$ (2.8c)

Here similarly the superscript $(i)$ denotes variables at the time stage
$t_{n}^{(i)}$ or $\hat{t}_{n}^{(i)}$, and the intermediate stage values for
$2\leq i\leq s$ are obtained from:

  * 1.

the predictor step:

$\displaystyle\frac{\tilde{E}^{(i)}-E^{n}}{\Delta t}$
$\displaystyle=\text{RHS}_{E}^{(i)}+a_{ii}\left[\alpha_{0}\Delta\tilde{E}^{(i)}+\tilde{\sigma}^{(i)}(\tilde{B}^{(i)}-\tilde{E}^{(i)})\right],$
(2.9a) $\displaystyle\frac{\tilde{B}^{(i)}-B^{n}}{\Delta t}$
$\displaystyle=\text{RHS}_{B}^{(i)}+a_{ii}\left[\beta_{0}\Delta\tilde{B}^{(i)}+z^{3}(\tilde{E}^{(i)}-\tilde{B}^{(i)})\right],$
(2.9b) $\displaystyle\tilde{T}^{(i)}$ $\displaystyle=(\tilde{B}^{(i)})^{1/4},$
(2.9c)

  * 2.

the corrector step:

$\displaystyle\frac{E^{(i)}-E^{n}}{\Delta t}$
$\displaystyle=\text{RHS}_{E}^{(i)}+a_{ii}\left[\alpha_{0}\Delta
E^{(i)}+\tilde{\sigma}^{(i)}(B^{(i)}-E^{(i)})\right],$ (2.10a)
$\displaystyle\frac{Q^{(i)}-Q^{n}}{\Delta t}$
$\displaystyle=\text{RHS}_{Q}^{(i)}+a_{ii}\gamma_{0}\Delta Q^{(i)},$ (2.10b)
$\displaystyle T^{(i)}$ $\displaystyle=Q^{(i)}-E^{(i)}.$ (2.10c)

The shorthand notations in (2.9)-(2.10) are defined as:

$\begin{cases}\text{RHS}_{E}^{(i)}=\sum\limits_{j=1}^{i-1}\left[\hat{a}_{ij}\left(\nabla\cdot(D^{(j)}_{r}\nabla
E^{(j)})-\alpha_{0}\Delta E^{(j)}\right)+a_{ij}\left(\alpha_{0}\Delta
E^{(j)}+\sigma^{(j)}(B^{(j)}-E^{(j)})\right)\right],\\\\[5.0pt]
\text{RHS}_{B}^{(i)}=\sum\limits_{j=1}^{i-1}\left[\hat{a}_{ij}\left(\mathcal{H}^{(j)}-\beta_{0}\Delta
B^{(j)}\right)+a_{ij}\left(\beta_{0}\Delta
B^{(j)}+z^{3}(E^{(j)}-B^{(j)})\right)\right],\\\\[5.0pt]
\text{RHS}_{Q}^{(i)}=\sum\limits_{j=1}^{i-1}\left[\hat{a}_{ij}\left(\nabla\cdot\left(D^{(j)}_{r}\nabla
E^{(j)}\right)+\nabla\cdot\left(D^{(j)}_{t}\nabla
T^{(j)}\right)-\gamma_{0}\Delta Q^{(j)}\right)+a_{ij}\gamma_{0}\Delta
Q^{(j)}\right].\end{cases}$

For the first stage, we take
$\tilde{E}^{(1)}=E^{n},\,\tilde{B}^{(1)}=(T^{n})^{4},\,E^{(1)}=E^{n},\,Q^{(1)}=E^{n}+T^{n}$.

###  2.3 Picard iteration

For the first-order IMEX scheme (2.2) and (2.5), or the high-order IMEX scheme
in the intermediate stages (2.9)-(2.10), each system is mildly nonlinear. Here
we will describe how to solve those mildly nonlinear systems with a simple
Picard iteration. The iteration does not rely on any specific spatial
discretization, so we keep space continuous first.

Taking the high-order IMEX scheme for the predictor step (2.9) as an example,
the updating $\tilde{E}$ and $\tilde{B}$ can be rewritten as:

$\displaystyle\left(\frac{1}{\Delta
t}-a_{ii}\alpha_{0}\Delta+a_{ii}\tilde{\sigma}^{(i)}\right)\tilde{E}^{(i)}-a_{ii}\tilde{\sigma}^{(i)}\tilde{B}^{(i)}=\text{RHS}_{E}^{(i)}+\frac{1}{\Delta
t}E^{n},$ (2.11a) $\displaystyle\left(\frac{1}{\Delta
t}-a_{ii}\beta_{0}\Delta+a_{ii}z^{3}\right)\tilde{B}^{(i)}-a_{ii}z^{3}\tilde{E}^{(i)}=\text{RHS}_{B}^{(i)}+\frac{1}{\Delta
t}B^{n}.$ (2.11b)

For a Picard iteration, starting from the iterative number $l=0$, we set
$\tilde{E}^{i,0}=E^{(i-1)}$ and $\tilde{B}^{i,0}=B^{(i-1)}$, and update
$\tilde{E}^{i,l+1}$ and $\tilde{B}^{i,l+1}$ from $\tilde{E}^{i,l}$ and
$\tilde{B}^{i,l}$ iteratively as:

$\displaystyle\left(\frac{1}{\Delta
t}-a_{ii}\alpha_{0}\Delta+a_{ii}\tilde{\sigma}^{i,l}\right)\tilde{E}^{i,l+1}-a_{ii}\tilde{\sigma}^{i,l}\tilde{B}^{i,l+1}=\text{RHS}_{E}^{(i)}+\frac{1}{\Delta
t}E^{n},$ (2.12a) $\displaystyle\left(\frac{1}{\Delta
t}-a_{ii}\beta_{0}\Delta+a_{ii}z^{3}\right)\tilde{B}^{i,l+1}-a_{ii}z^{3}\tilde{E}^{i,l+1}=\text{RHS}_{B}^{(i)}+\frac{1}{\Delta
t}B^{n}.$ (2.12b)

As observed, when $\tilde{\sigma}^{i,l}$ is set based on the previous
iterative step, (2.12) gives rise to a linear system for $\tilde{E}^{i,l+1}$
and $\tilde{B}^{i,l+1}$ that exhibits diagonal dominance. This property arises
from the positivity of $\sigma$ and $z({\bf x})$, resulting in rapid
convergence when employing iterative methods to solve it.

Similar to (2.11), the corrector step (2.10) for $E$ and $Q$ can be rewritten
as:

$\displaystyle\left(\frac{1}{\Delta
t}-a_{ii}\alpha_{0}\Delta+a_{ii}\tilde{\sigma}^{(i)}\right)E^{(i)}-a_{ii}\tilde{\sigma}^{(i)}B^{(i)}=\text{RHS}_{E}^{(i)}+\frac{1}{\Delta
t}E^{n},$ (2.13a) $\displaystyle\left(\frac{1}{\Delta
t}-a_{ii}\gamma_{0}\Delta\right)Q^{(i)}=\text{RHS}_{Q}^{(i)}+\frac{1}{\Delta
t}Q^{n}.$ (2.13b)

Here, $Q^{(i)}$ can be readily obtained by solving the linear system (2.13b),
which also possesses diagonal dominance. Subsequently, we solve (2.13a)
through a Picard iteration, with initial values
$B^{i,0}=\tilde{B}^{(i)}=(\tilde{T}^{(i)})^{4}$, as follows:

$\displaystyle\left(\frac{1}{\Delta
t}-a_{ii}\alpha_{0}\Delta+a_{ii}\tilde{\sigma}^{(i)}\right){E}^{i,l+1}=a_{ii}\tilde{\sigma}^{(i)}{B}^{i,l}+\text{RHS}_{E}^{(i)}+\frac{1}{\Delta
t}E^{n},$ (2.14a) $\displaystyle B^{i,l+1}=(Q^{(i)}-E^{i,l+1})^{4}.$ (2.14b)

We solve the linear system (2.14a) to get $E^{i,l+1}$, and then update
$B^{i,l+1}$ from (2.14b). Since the initial values are set as
$E^{i,0}=\tilde{E}^{(i)}$ and $B^{i,0}=(Q^{(i)}-E^{i,0})^{4}$ from the
predictor step, a fast convergence can be obtained for the corrector step.

The above procedures are similar for the first order IMEX scheme (2.2) and
(2.5), we omit them to save space.

##  3 LDG spatial discretization

For the first-order IMEX scheme (2.2)-(2.5), or the high-order IMEX scheme
(2.8)-(2.10), we can couple them with any front capturing spatial
discretizations [28, 5, 1, 49, 58]. In this work, we utilize an LDG finite
element method. The LDG method offers great flexibility for h-p adaptivity and
excels at capturing sharp gradient propagations.

###  3.1 Some notations

For a computational domain $\Omega\subseteq\mathbb{R}^{2}$, we consider a
partition $\mathcal{T}_{h}$ of $\Omega$ with a set of non-overlapping
rectangular elements $\\{I_{i,j}\\}$, which can cover the whole domain
$\Omega$. Here $I_{ij}=I_{i}\times I_{j}$,
$I_{i}=\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]$ and
$I_{j}=\left[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}\right]$ for
$i=1,2,\cdots,N_{x},j=1,2,\cdots,N_{y}$. We denote the element length and
width as $h_{i}^{x}=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}$,
$h_{j}^{y}=y_{j+\frac{1}{2}}-y_{j-\frac{1}{2}}$, respectively.
$h=\max\limits_{i,j}\\{h_{i}^{x},h_{j}^{y}\\}$ is the maximum edge size of
these elements. The center of the element $I_{i,j}$ is $(x_{i},y_{j})$. We
also assume that $\mathcal{T}_{h}$ is quasi-uniform, namely,
$\max\limits_{i}\\{h/h_{i}^{x}\\}$ and $\max\limits_{j}\\{h/h_{j}^{y}\\}$ are
upper bounded by a given positive constant.

With the above partition, we follow [51] to give some notations which will be
used in the following. Given any non-negative integer vector
$\mathbf{k}=\left(k_{1},k_{2}\right)$, we define a finite-dimensional discrete
piecewise polynomial space as follows

$W_{h}^{\mathbf{k}}=\\{u\in
L^{2}(\Omega):u|_{K}\in\mathcal{Q}^{\mathbf{k}}(K),\forall
K\in\mathcal{T}_{h}\\},$

where $\mathcal{Q}^{\mathbf{k}}(K)$ consists of tensor product polynomials of
degree not exceeding $k_{\ell}$ along the $\ell$-th direction on each element
$K$, for $\ell=1,2$. Besides, we denote
${\mathbf{W}}_{h}^{\mathbf{k}}=W_{h}^{\mathbf{k}}\times W_{h}^{\mathbf{k}}$ as
a vector space, where each component belongs to $W_{h}^{\mathbf{k}}$. We
define a unit normal vector $\mathbf{n}^{e}$ on each edge $e$ of
$\mathcal{T}_{h}$ as follows: if $e\in\partial\Omega$, $\mathbf{n}^{e}$ is
defined as the unit normal vector pointing outside of $\Omega$; for an
interior edge $e=\partial K^{+}\cap\partial K^{-}$, the outward unit normal
vectors of $e$ taken from the elements $K^{+}$ and $K^{-}$ are denoted by
$\mathbf{n}^{+}$ and $\mathbf{n}^{-}$, respectively. Here we fix
$\mathbf{n}^{e}$ as one of $\mathbf{n}^{\pm}$. If we denote $u^{+}$ and
$u^{-}$ as the values of a function $u$ on $e$, taken from $K^{+}$ and $K^{-}$
respectively, then the jump $[[u]]$ over an edge $e$ for a scalar-valued
function $u$ is defined as

$[[u]]|_{e}=-\left(u^{+}\mathbf{n}^{+}+u^{-}\mathbf{n}^{-}\right)\cdot\mathbf{n}^{e}.$

For a vector-valued function $\mathbf{v}$, the jump
$[[\mathbf{v}\cdot\mathbf{n}]]$ is defined as

$[[\mathbf{v}\cdot\mathbf{n}]]|_{e}=-\left(\mathbf{v}^{+}\cdot\mathbf{n}^{+}+\mathbf{v}^{-}\cdot\mathbf{n}^{-}\right)\mathbf{n}^{e}\cdot\mathbf{n}^{+}.$

Accordingly, we express the averages of $u$ and $\mathbf{v}\cdot\mathbf{n}$ as

$\\{\\{u\\}\\}|_{e}=-\frac{1}{2}\left(u^{+}\mathbf{n}^{+}-u^{-}\mathbf{n}^{-}\right)\cdot\mathbf{n}^{e},\quad\\{\\{\mathbf{v}\cdot\mathbf{n}\\}\\}|_{e}=-\frac{1}{2}\left(\mathbf{v}^{+}\cdot\mathbf{n}^{+}-\mathbf{v}^{-}\cdot\mathbf{n}^{-}\right)\mathbf{n}^{e}\cdot\mathbf{n}^{+}.$

In this work, we take $\mathbf{n}^{e}=\mathbf{n}^{-}$, then

$[[u]]|_{e}=u^{+}-u^{-},\quad\\{\\{u\\}\\}|_{e}=\frac{1}{2}\left(u^{+}+u^{-}\right),$

and

$[[\mathbf{v}\cdot\mathbf{n}]]|_{e}=-(\mathbf{v}^{+}-\mathbf{v}^{-})\cdot\mathbf{n}^{-},\quad\\{\\{\mathbf{v}\cdot\mathbf{n}\\}\\}|_{e}=-\frac{1}{2}\left(\mathbf{v}^{+}+\mathbf{v}^{-}\right)\cdot\mathbf{n}^{-}.$

In our implementation, we use an orthogonal basis of $W_{h}^{\bm{k}}$ with a
uniform rectangular partition $h_{x}=h_{x}^{i}$ for $1\leq i\leq N_{x}$ and
$h_{y}=h_{y}^{j}$ for $1\leq j\leq N_{y}$. In this case, the numerical
solution can be expressed as

$u_{h}(x,y)=\sum\limits_{i=1}^{N_{x}}\sum\limits_{j=1}^{N_{y}}\sum\limits_{m=1}^{k_{1}}\sum\limits_{n=1}^{k_{2}}u_{i,j}^{m,n}H_{i}^{m}(x)H_{j}^{n}(y),\quad(x,y)\in\Omega,$

$u_{h}$ is $E_{h}$, $T_{h}$, or $B_{h}$ respectively. The local basis of
$\mathcal{Q}^{\mathbf{k}}(K)$ on each element $K$ is denoted as
$H_{i}^{m}(x)H_{j}^{n}(y)$ for $K=I_{i,j}\in\mathcal{T}_{h}$, and
$\mathbf{k}=(k_{1},k_{2})$. For example, up to third order, the local bases
are chosen as follows

$\displaystyle
H_{i}^{1}(x)=1,\,H_{i}^{2}(x)=\frac{x-x_{i}}{h_{x}},\,H_{i}^{3}(x)=\left(\frac{x-x_{i}}{h_{x}}\right)^{2}-\frac{1}{12},\quad
x\in I_{i},$ $\displaystyle
H_{j}^{1}(y)=1,\,H_{j}^{2}(y)=\frac{y-y_{j}}{h_{y}},\,H_{j}^{3}(y)=\left(\frac{y-y_{j}}{h_{y}}\right)^{2}-\frac{1}{12},\quad
y\in I_{j},$

with zero extension outside the cell $I_{i}$ or $I_{j}$ respectively.

###  3.2 First order IMEX-LDG scheme

With the above notations, a fully-discrete LDG scheme utilizing a first-order
IMEX scheme (2.2) and (2.5) is defined as follows. First, for the predictor
step (2.2), the scheme reads: given $E_{h}^{n},T_{h}^{n},B_{h}^{n}\in
W_{h}^{\mathbf{k}}$ and $\mathbf{p}^{n}_{h}$, $\mathbf{q}^{n}_{h}$,
$\mathbf{r}^{n}_{h}\in{\mathbf{W}}_{h}^{\mathbf{k}}$, we find
$\tilde{E}_{h}^{n+1},\tilde{B}_{h}^{n+1}\in W_{h}^{\mathbf{k}}$ and
$\tilde{\mathbf{p}}^{n+1}_{h}$,
$\tilde{\mathbf{q}}^{n+1}_{h}\in{\mathbf{W}}_{h}^{\mathbf{k}}$, such that for
any $\mu$, $\nu$, $\zeta\in W_{h}^{\mathbf{k}}$, and $\boldsymbol{\upsilon}$,
$\boldsymbol{\xi}\in{\mathbf{W}}_{h}^{\mathbf{k}}$, such that

$\displaystyle\frac{1}{\Delta
t}\left(\tilde{E}_{h}^{n+1}-E_{h}^{n},\mu\right)$
$\displaystyle=\mathcal{G}_{h}^{n}(\mu)-\alpha_{0}\mathcal{L}^{n}_{h,\mathbf{p}}(\mu)+\alpha_{0}\mathcal{L}^{n+1}_{h,\tilde{\mathbf{p}}}(\mu)+\left(\tilde{\sigma}_{h}^{n+1}(\tilde{B}^{n+1}_{h}-\tilde{E}^{n+1}_{h}),\mu\right),$
(3.1a) $\displaystyle\frac{1}{\Delta
t}\left(\tilde{B}_{h}^{n+1}-B_{h}^{n},\nu\right)$
$\displaystyle=\mathcal{H}_{h}^{n}(\nu)-\beta_{0}\mathcal{L}^{n}_{h,\mathbf{q}}(\nu)+\beta_{0}\mathcal{L}^{n+1}_{h,\tilde{\mathbf{q}}}(\nu)+\left(z^{3}(\tilde{B}^{n+1}_{h}-\tilde{E}^{n+1}_{h}),\nu\right),$
(3.1b) $\displaystyle(\tilde{T}_{h}^{n+1},\zeta)$
$\displaystyle=((\tilde{B}_{h}^{n+1})^{1/4},\zeta),$ (3.1c)
$\displaystyle(\tilde{\mathbf{p}}_{h}^{n+1},\boldsymbol{\upsilon})$
$\displaystyle=\boldsymbol{\mathcal{K}}_{h}(\tilde{E}_{h}^{n+1},\boldsymbol{\upsilon}),$
(3.1d) $\displaystyle(\tilde{\mathbf{q}}_{h}^{n+1},\boldsymbol{\xi})$
$\displaystyle=\boldsymbol{\mathcal{K}}_{h}(\tilde{B}_{h}^{n+1},\boldsymbol{\xi}),$
(3.1e)

with

$\mathcal{G}_{h}^{n}(\mu)=\boldsymbol{\mathcal{C}}_{h}((D_{r}\mathbf{p}_{h})^{n},\mu),\quad\mathcal{H}_{h}^{n}(\nu)=\boldsymbol{\mathcal{D}}_{h}(T_{h}^{n},\mathbf{r}_{h}^{n},\nu),\quad$

where

$\mathcal{L}^{n+1}_{h,\tilde{\mathbf{p}}}(\mu)=\boldsymbol{\mathcal{L}}_{h}(\tilde{\mathbf{p}}_{h}^{n+1},\mu),\quad\mathcal{L}^{n+1}_{h,\tilde{\mathbf{q}}}(\nu)=\boldsymbol{\mathcal{L}}_{h}(\tilde{\mathbf{q}}_{h}^{n+1},\nu),\quad\mathcal{L}^{n}_{h,\mathbf{p}}(\mu)=\boldsymbol{\mathcal{L}}_{h}(\mathbf{p}_{h}^{n},\mu),\quad\mathcal{L}^{n}_{h,\mathbf{q}}(\nu)=\boldsymbol{\mathcal{L}}_{h}(\mathbf{q}_{h}^{n},\nu).$

Correspondingly, for the corrector step (2.5), the scheme is defined as: we
look for $E_{h}^{n+1},Q_{h}^{n+1}\in W_{h}^{\mathbf{k}}$ and
$\mathbf{p}^{n+1}_{h}$,
$\mathbf{w}^{n+1}_{h}\in{\mathbf{W}}_{h}^{\mathbf{k}}$, for any $\mu$,
$\phi\in W_{h}^{\mathbf{k}},$ and $\boldsymbol{\upsilon}$,
$\boldsymbol{\eta}\in{\mathbf{W}}_{h}^{\mathbf{k}}$, such that

$\displaystyle\frac{1}{\Delta t}\left(E_{h}^{n+1}-E_{h}^{n},\mu\right)$
$\displaystyle=\mathcal{G}_{h}^{n}(\mu)-\alpha_{0}\mathcal{L}^{n}_{h,\mathbf{p}}(\mu)+\alpha_{0}\mathcal{L}^{n+1}_{h,\mathbf{p}}(\mu)+\left(\tilde{\sigma}_{h}^{n+1}(B^{n+1}_{h}-E^{n+1}_{h}),\mu\right),$
(3.2a) $\displaystyle\frac{1}{\Delta
t}\left(Q_{h}^{n+1}-Q_{h}^{n},\phi\right)$
$\displaystyle=\mathcal{E}_{h}^{n}(\phi)-\gamma_{0}\mathcal{L}^{n}_{h,\mathbf{w}}(\phi)+\gamma_{0}\mathcal{L}^{n+1}_{h,\mathbf{w}}(\phi),$
(3.2b) $\displaystyle T^{n+1}_{h}$ $\displaystyle=Q^{n+1}_{h}-E^{n+1}_{h},$
(3.2c) $\displaystyle(\mathbf{p}_{h}^{n+1},\boldsymbol{\upsilon})$
$\displaystyle=\boldsymbol{\mathcal{K}}_{h}(E_{h}^{n+1},\boldsymbol{\upsilon}),$
(3.2d) $\displaystyle(\mathbf{w}_{h}^{n+1},\boldsymbol{\eta})$
$\displaystyle=\boldsymbol{\mathcal{K}}_{h}(Q_{h}^{n+1},\boldsymbol{\eta}),$
(3.2e)

with

$\mathcal{E}_{h}^{n}(\phi)=\boldsymbol{\mathcal{F}}_{h}((D_{t}\mathbf{r})_{h}^{n},\phi)+\boldsymbol{\mathcal{C}}_{h}((D_{r}\mathbf{p}_{h})^{n},\phi),$

where

$\mathcal{L}^{n}_{h,\mathbf{p}}(\mu)=\boldsymbol{\mathcal{L}}_{h}(\mathbf{p}_{h}^{n},\mu),\quad\mathcal{L}^{n}_{h,\mathbf{w}}(\phi)=\boldsymbol{\mathcal{L}}_{h}(\mathbf{w}_{h}^{n},\phi).$

Those operators in (3.1)-(3.2) are defined as:

$\displaystyle\boldsymbol{\mathcal{L}}_{h}\left(\mathbf{p}_{h},\mu\right)=$
$\displaystyle-\sum\limits_{K\in\mathcal{T}_{h}}\int_{K}\mathbf{p}_{h}\cdot\nabla\mu
d{\bf x}+\sum\limits_{K\in\mathcal{T}_{h}}\int_{\partial
K}\mathbf{n}^{e}\cdot\hat{\mathbf{p}}_{h}\mu ds,$
$\displaystyle\boldsymbol{\mathcal{K}}_{h}\left(E_{h},\boldsymbol{\upsilon}\right)=$
$\displaystyle-\sum\limits_{K\in\mathcal{T}_{h}}\int_{K}E_{h}\nabla\cdot\boldsymbol{\upsilon}d{\bf
x}+\sum\limits_{K\in\mathcal{T}_{h}}\int_{\partial
K}\mathbf{n}^{e}\cdot\boldsymbol{\upsilon}\hat{E}_{h}ds,$
$\displaystyle\boldsymbol{\mathcal{C}}_{h}(D_{r}\mathbf{p}_{h},\mu)=$
$\displaystyle-\sum\limits_{K\in\mathcal{T}_{h}}\int_{K}D_{r,h}\mathbf{p}_{h}\cdot\nabla\mu
d{\bf x}+\sum\limits_{K\in\mathcal{T}_{h}}\int_{\partial
K}\hat{D}_{r,h}\mathbf{n}^{e}\cdot\hat{\mathbf{p}}_{h}\mu ds,$
$\displaystyle\boldsymbol{\mathcal{F}}_{h}(D_{t}\mathbf{r}_{h},\phi)=$
$\displaystyle-\sum\limits_{K\in\mathcal{T}_{h}}\int_{K}D_{t,h}\mathbf{r}_{h}\cdot\nabla\phi
d{\bf x}+\sum\limits_{K\in\mathcal{T}_{h}}\int_{\partial
K}\hat{D}_{t,h}\mathbf{n}^{e}\cdot\hat{\mathbf{r}}_{h}\phi ds.$
$\displaystyle\boldsymbol{\mathcal{D}}_{h}(T_{h},\mathbf{r}_{h},\nu)=$
$\displaystyle-4\kappa\sum\limits_{K\in\mathcal{T}_{h}}\int_{K}\left(T^{\frac{11}{2}}_{h}\mathbf{r}_{h}\cdot\nabla\nu+3T^{\frac{9}{2}}_{h}|\mathbf{r}_{h}|^{2}\nu\right)d{\bf
x}$ $\displaystyle+4\kappa\sum\limits_{K\in\mathcal{T}_{h}}\int_{\partial
K}\hat{d}_{h}(T_{h};\mathbf{n}^{e})\mathbf{n}^{e}\cdot\hat{\mathbf{r}}_{h}\nu
ds.$

In the above and below, $(\cdot,\cdot)$ denotes the usual $L^{2}$ inner
product on $\Omega$. We can find that $\mathbf{p}_{h}$, $\mathbf{q}_{h}$,
$\mathbf{r}_{h}$, and $\mathbf{w}_{h}$ approximate $\nabla E$, $\nabla B$,
$\nabla T$, and $\nabla Q$ respectively.

The variables with $\hat{\cdot}$ are numerical fluxes. For such a diffusion
system, we adopt an alternating left-right flux, namely

$\displaystyle\hat{u}_{h}=u_{h}^{-},\quad\hskip 39.83385pt\text{for a scalar-
valued function }\,u=E,B,T,Q,$
$\displaystyle\mathbf{n}^{e}\cdot\hat{\mathbf{v}}_{h}=\mathbf{n}^{-}\cdot\mathbf{v}_{h}^{+},\quad\text{for
a vector-valued function
}\,\mathbf{v}=\mathbf{p},\mathbf{q},\mathbf{r},\mathbf{w}.$

The alternating right-left flux or central fluxes as in [10, 17, 48, 51] can
also be used. As used in [12, 52, 53],
$\hat{D}_{r,h}=\\{\\{D_{r,h}\\}\\}|_{e}$ is defined as a central numerical
flux in approximating of $D_{r,h}$, and $\hat{d}_{h}(T_{h};\mathbf{n}^{e})$ in
$\boldsymbol{\mathcal{D}}_{h}(T_{h},\mathbf{r}_{h},\nu)$ is an approximation
to $T^{\frac{11}{2}}$,

$\hat{d}_{h}(T_{h};\mathbf{n}^{e})=\left\\{\begin{array}[]{ll}[[\frac{2}{13}T_{h}^{\frac{13}{2}}]]|_{e}/[[T_{h}]]|_{e},&\text{
if }[[T_{h}]]|_{e}\neq 0,\\\\[10.0pt] (T_{h}^{\frac{11}{2}})^{-},&\text{
otherwise. }\end{array}\right.$

$\hat{D}_{t,h}$ is chosen to approximate $D_{t}$ in the same manner as
$\hat{d}_{h}(T_{h};\mathbf{n}^{e})$.

###  3.3 High order IMEX-LDG scheme

By applying a high-order IMEX RK scheme (2.7) to (3.1)-(3.2), similar to the
first-order scheme, the updating of the solutions at $t^{n+1}$ from $t^{n}$ is
given as follows:

$\displaystyle\frac{1}{\Delta t}\left(E_{h}^{n+1}-E_{h}^{n},\mu\right)$
$\displaystyle=\sum\limits_{i=1}^{s}\left[\hat{b}_{i}\left(\mathcal{G}_{h}^{(i)}(\mu)-\alpha_{0}\mathcal{L}^{(i)}_{h,\mathbf{p}}(\mu)\right)+b_{i}\left(\alpha_{0}\mathcal{L}^{(i)}_{h,\mathbf{p}}(\mu)+\left({\sigma}_{h}^{(i)}(B^{(i)}_{h}-E^{(i)}_{h}),\mu\right)\right)\right],$
(3.3a) $\displaystyle\frac{1}{\Delta
t}\left(Q_{h}^{n+1}-Q_{h}^{n},\phi\right)$
$\displaystyle=\sum\limits_{i=1}^{s}\left[\hat{b}_{i}\left(\mathcal{E}_{h}^{(i)}(\phi)-\gamma_{0}\mathcal{L}^{(i)}_{h,\mathbf{w}}(\phi)\right)+b_{i}\gamma_{0}\mathcal{L}^{(i)}_{h,\mathbf{w}}(\phi)\right],$
(3.3b) $\displaystyle T^{n+1}_{h}$ $\displaystyle=Q^{n+1}_{h}-E^{n+1}_{h},$
(3.3c) $\displaystyle(\mathbf{p}_{h}^{n+1},\boldsymbol{\upsilon})$
$\displaystyle=\boldsymbol{\mathcal{K}}_{h}(E_{h}^{n+1},\boldsymbol{\upsilon}),$
(3.3d) $\displaystyle(\mathbf{w}_{h}^{n+1},\boldsymbol{\eta})$
$\displaystyle=\boldsymbol{\mathcal{K}}_{h}(Q_{h}^{n+1},\boldsymbol{\eta}).$
(3.3e)

The intermediate stage values for $2\leq i\leq s$ are obtained from:

  * 1.

the predictor step:

$\displaystyle\frac{1}{\Delta
t}\left(\tilde{E}_{h}^{(i)}-E_{h}^{n},\mu\right)$
$\displaystyle=\text{RHS}_{E}^{(i)}(\mu)+a_{ii}\left[\alpha_{0}\mathcal{L}^{(i)}_{h,\tilde{\mathbf{p}}}(\mu)+\left(\tilde{\sigma}_{h}^{(i)}(\tilde{B}^{(i)}_{h}-\tilde{E}^{(i)}_{h}),\mu\right)\right],$
(3.4a) $\displaystyle\frac{1}{\Delta
t}\left(\tilde{B}_{h}^{(i)}-B_{h}^{n},\nu\right)$
$\displaystyle=\text{RHS}_{B}^{(i)}(\nu)+a_{ii}\left[\beta_{0}\mathcal{L}^{(i)}_{h,\tilde{\mathbf{q}}}(\nu)+\left(z^{3}(\tilde{B}^{(i)}_{h}-\tilde{E}^{(i)}_{h}),\nu\right)\right],$
(3.4b) $\displaystyle(\tilde{T}_{h}^{(i)},\zeta)$
$\displaystyle=((\tilde{B}_{h}^{(i)})^{1/4},\zeta),$ (3.4c)
$\displaystyle(\tilde{\mathbf{p}}_{h}^{(i)},\boldsymbol{\upsilon})$
$\displaystyle=\boldsymbol{\mathcal{K}}_{h}(\tilde{E}_{h}^{(i)},\boldsymbol{\upsilon}),$
(3.4d) $\displaystyle(\tilde{\mathbf{q}}_{h}^{(i)},\boldsymbol{\xi})$
$\displaystyle=\boldsymbol{\mathcal{K}}_{h}(\tilde{B}_{h}^{(i)},\boldsymbol{\xi});$
(3.4e)

  * 2.

the corrector step:

$\displaystyle\frac{1}{\Delta t}\left(E_{h}^{(i)}-E_{h}^{n},\mu\right)$
$\displaystyle=\text{RHS}_{E}^{(i)}(\mu)+a_{ii}\left[\alpha_{0}\mathcal{L}^{(i)}_{h,\tilde{\mathbf{p}}}(\mu)+\left(\tilde{\sigma}_{h}^{(i)}(B^{(i)}_{h}-E^{(i)}_{h}),\mu\right)\right],$
(3.5a) $\displaystyle\frac{1}{\Delta
t}\left(Q_{h}^{(i)}-Q_{h}^{n},\phi\right)$
$\displaystyle=\text{RHS}_{Q}^{(i)}(\phi)+a_{ii}\gamma_{0}\mathcal{L}^{(i)}_{h,\mathbf{w}}(\phi),$
(3.5b) $\displaystyle T^{(i)}_{h}$ $\displaystyle=Q^{(i)}_{h}-E^{(i)}_{h},$
(3.5c) $\displaystyle(\mathbf{p}_{h}^{(i)},\boldsymbol{\upsilon})$
$\displaystyle=\boldsymbol{\mathcal{K}}_{h}(E_{h}^{(i)},\boldsymbol{\upsilon}),$
(3.5d) $\displaystyle(\mathbf{w}_{h}^{(i)},\boldsymbol{\eta})$
$\displaystyle=\boldsymbol{\mathcal{K}}_{h}(Q_{h}^{(i)},\boldsymbol{\eta}).$
(3.5e)

The shorthand notations in (3.4)-(3.5) are defined as:

$\begin{cases}\text{RHS}_{E}^{(i)}(\mu)=\sum\limits_{j=1}^{i-1}\left[\hat{a}_{ij}\left(\mathcal{G}_{h}^{(j)}(\mu)-\alpha_{0}\mathcal{L}^{(j)}_{h,\mathbf{p}}(\mu)\right)+a_{ij}\left(\alpha_{0}\mathcal{L}^{(j)}_{h,\mathbf{p}}(\mu)+\left(\tilde{\sigma}_{h}^{(j)}(B^{(j)}_{h}-E^{(j)}_{h}),\mu\right)\right)\right],\\\\[5.0pt]
\text{RHS}_{B}^{(i)}(\nu)=\sum\limits_{j=1}^{i-1}\left[\hat{a}_{ij}\left(\mathcal{H}_{h}^{(j)}(\nu)-\beta_{0}\mathcal{L}^{(j)}_{h,\mathbf{q}}(\nu)\right)+a_{ij}\left(\beta_{0}\mathcal{L}^{(j)}_{h,\mathbf{q}}(\nu)+\left(z^{3}(\tilde{B}^{(j)}_{h}-\tilde{E}^{(j)}_{h}),\nu\right)\right)\right],\\\\[5.0pt]
\text{RHS}_{Q}^{(i)}(\phi)=\sum\limits_{j=1}^{i-1}\left[\hat{a}_{ij}\left(\mathcal{E}_{h}^{(j)}(\phi)-\gamma_{0}\mathcal{L}^{(j)}_{h,\mathbf{w}}(\phi)\right)+a_{ij}\gamma_{0}\mathcal{L}^{(j)}_{h,\mathbf{w}}(\phi)\right].\end{cases}$

As we can observe, each of the intermediate steps (3.4)-(3.5) in a high-order
scheme shares the same structure as the first-order scheme (3.1)-(3.2). Thus,
the first-order scheme provides a general framework for coupling with a multi-
stage IMEX RK method, making it convenient to extend to higher orders.
Furthermore, the choice of operators and numerical fluxes remains consistent
with the first-order scheme.

###  3.4 Picard iteration for a full scheme

To represent our scheme as a mildly nonlinear system more effectively, we
introduce notations for matrices and vectors. We define a long vector
$\bm{u_{c}}$ to represent the coefficients of a two-dimensional numerical
solution $u_{h}$, as follows:

$\bm{u_{c}}=\left(u_{1,1}^{1,1},\cdots,u_{1,1}^{k_{1},k_{2}},u_{2,1}^{1,1},\cdots,u_{2,1}^{k_{1},k_{2}},\cdots,u_{N_{x},1}^{1,1},\cdots,u_{N_{x},1}^{k_{1},k_{2}},u_{1,2}^{1,1},\cdots,u_{N_{x},N_{y}}^{k_{1},k_{2}}\right)^{T}.$
(3.6)

Such a vector $\bm{u_{c}}$ is the value to be updated by a numerical method in
the LDG framework. For instance, when taking the nonlinear integrals
$\left(\sigma_{h}E_{h},\mu\right)$ and $\left(\sigma_{h}B_{h},\mu\right)$ in
(3.1) and (3.2), practical approximations are employed. These approximations
involve Gaussian quadrature integration along each direction, facilitated by
the use of an interpolation operator $\mathcal{I}_{h}$ within each cell
$K=I_{i}\times I_{j}$, i.e.

$\left(\sigma_{h}E_{h},\mu\right)_{K}\approx\int_{K}\mathcal{I}_{h}\left(\sigma_{h}E_{h}\mu\right)d{\bf
x}=\dfrac{h_{x}h_{y}}{4}\sum\limits_{G=1}^{k_{1}+1}\sum\limits_{G^{\prime}=1}^{k_{2}+1}\hat{\omega}_{G}\hat{\omega}_{G^{\prime}}(\sigma_{h}E_{h}\mu)\left(\frac{h_{x}}{2}\eta_{G}+x_{i},\frac{h_{y}}{2}\eta_{G^{\prime}}+y_{j}\right),$

where $\eta_{G}$ are Gaussian quadrature points on the interval $[-1,1]$ with
corresponding weights $\hat{\omega}_{G}$. We define
$k_{1,2}=(k_{1}+1)(k_{2}+1)$, therefore, by taking $\mu$ over the local bases
$H^{m}_{i}(x)H^{n}_{j}(y)$ in each cell, $(\sigma_{h}E_{h},\mu)$ can be
formatted as a matrix $\Lambda_{h,T}$ multiplied by a long vector
$\bm{E_{c}}$, that is $\Lambda_{h,T}\bm{E_{c}}$, where $\Lambda_{h,T}$ is a
$k_{1,2}$-block at most tri-diagonal matrix in the size of $N_{x}N_{y}$. Here
$\sigma$ is $T$-dependent, so that $\Lambda_{h,T}$ depends on $T$, and
$\Lambda_{h,T}\bm{E_{c}}$ is nonlinear.

With similar matrix-by-vector notations, the high order fully-discrete IMEX-
LDG scheme (3.4), corresponding to (2.11), can be written in the following
form

$\displaystyle\left(\frac{1}{\Delta
t}M_{h}-a_{ii}\alpha_{0}S_{h}+a_{ii}\Lambda_{h,\tilde{T}}^{(i)}\right)\tilde{\bm{E}}_{\bm{c}}^{(i)}-a_{ii}\Lambda_{h,\tilde{T}}^{(i)}\tilde{\bm{B}}_{\bm{c}}^{(i)}=F_{E}^{(i)},$
(3.7a) $\displaystyle\left(\frac{1}{\Delta
t}M_{h}+a_{ii}Z_{h}-a_{ii}\beta_{0}S_{h}\right)\tilde{\bm{B}}_{\bm{c}}^{(i)}-a_{ii}Z_{h}\tilde{\bm{E}}_{\bm{c}}^{(i)}=F_{B}^{(i)},$
(3.7b)

where $M_{h}\bm{E_{c}}$, $M_{h}\bm{B_{c}}$, $S_{h}\bm{E_{c}}$,
$S_{h}\bm{B_{c}}$, $Z_{h}\bm{E_{c}}$ and $Z_{h}\bm{B_{c}}$ are matrix-by-
vector forms of $(E_{h},\mu)$, $(B_{h},\nu)$,
$(\mathcal{L}_{h,\mathbf{p}},\mu)$, $(\mathcal{L}_{h,\mathbf{q}},\nu)$,
$(z^{3}E_{h},\nu)$ and $(z^{3}B_{h},\nu)$, respectively. The right-hand-side
long vectors $F_{E}^{(i)}$ and $F_{B}^{(i)}$ are corresponding to
$(\text{RHS}_{E}^{(i)},\mu)+(E_{h},\mu)/\Delta t$ and
$(\text{RHS}_{B}^{(i)},\nu)+(B_{h},\nu)/\Delta t$ respectively.

Similar to (3.7), the corrector step (3.5), corresponding to (2.13), can be
written as

$\displaystyle\left(\frac{1}{\Delta
t}M_{h}-a_{ii}\alpha_{0}S_{h}+a_{ii}\Lambda_{h,\tilde{T}}^{(i)}\right)\bm{E_{c}}^{(i)}-a_{ii}\Lambda_{h,\tilde{T}}^{(i)}\bm{B_{c}}^{(i)}=F_{E}^{(i)},$
(3.8a) $\displaystyle\left(\frac{1}{\Delta
t}M_{h}-a_{ii}\gamma_{0}S_{h}\right)\bm{Q_{c}}^{(i)}=F^{(i)}_{Q}+\frac{1}{\Delta
t}M_{h}\bm{Q_{c}}^{n}.$ (3.8b)

For the predictor step, (3.7) can be expressed as a nonlinear algebraic system

$A(\tilde{U}^{(i)})\tilde{U}^{(i)}=F(U^{n},U^{(1)},U^{(2)},\ldots,U^{(i-1)}),$
(3.9)

with

$U=\left[\begin{array}[]{c}\bm{E_{c}}\\\ \bm{B_{c}}\end{array}\right],\quad
F=\left[\begin{array}[]{c}F_{E}^{(i)}\\\ F_{B}^{(i)}\end{array}\right],\quad
A(U)=\left[\begin{array}[]{cc}M_{1}(U)&M_{2}(U)\\\
M_{3}&M_{4}\end{array}\right].$

The corresponding submatrices in $A(U)$ are given by

$\displaystyle M_{1}(U)=\frac{1}{\Delta
t}M_{h}-a_{ii}\alpha_{0}S_{h}+a_{ii}\Lambda_{h,T},\quad
M_{2}(U)=-a_{ii}\Lambda_{h,T},$ $\displaystyle M_{3}=-a_{ii}Z_{h},\quad
M_{4}=\frac{1}{\Delta t}M_{h}+a_{ii}Z_{h}-a_{ii}\beta_{0}S_{h}.$

As we can see, only $M_{1}(U)$ and $M_{2}(U)$ mildly depend on $U$ due to
$\Lambda_{h,T}$. Hence an iteration procedure, corresponding to (2.12), can be
written as

$A(U^{i,l})U^{i,l+1}=F(U^{n},U^{(1)},U^{(2)},\ldots,U^{(i-1)}),$ (3.10)

with the initial value taking to be $U^{i,0}=U^{(i-1)}$.

A detailed procedure of the Picard iteration for solving (3.10) is as follows:
starting from an iterative number $l=0$, we set $E_{h}^{i,0}=E_{h}^{(i-1)}$,
$B_{h}^{i,0}=B_{h}^{(i-1)}$, and update the unknowns $B_{h}^{i,l+1}$ and
$E_{h}^{i,l+1}$ iteratively through the following three steps:

  * 1.

Step 1: compute $A(U^{i,l})$

$A(U^{i,l})=\left[\begin{array}[]{cc}M_{1}(U^{i,l})&M_{2}(U^{i,l})\\\
M_{3}&M_{4}\end{array}\right];$

  * 2.

Step 2: if $z({\bf x})$ is piecewise constant, $M_{3}$ is an easy invertible
diagonal matrix. From (3.10), by Gauss elimination, we update $E_{h}^{i,l+1}$
and $B_{h}^{i,l+1}$ with their long vectors by

$\begin{cases}\left(M_{1}{(M_{3})}^{-1}M_{4}-M_{2}\right)\bm{B_{c}}^{i,l+1}=M_{1}{(M_{3})}^{-1}F_{B}^{(i)}-F_{E}^{(i)},\\\
M_{3}\bm{E_{c}}^{i,l+1}=F_{B}^{(i)}-M_{4}\bm{B_{c}}^{i,l+1}.\end{cases}$
(3.11)

Otherwise, we solve (3.10) with the GMRES method [39].

  * 3.

Step 3: we take $\tilde{E}_{h}^{(i)}=E_{h}^{i,l+1}$ and
$\tilde{B}_{h}^{(i)}=B_{h}^{i,l+1}$ when iteration stops. The stop criteria is

$\|\bm{E_{c}}^{i,l+1}-\bm{E_{c}}^{i,l}\|<\delta.$

In our numerical tests, we take an $L^{2}$ norm and a threshold
$\delta=10^{-8}$.

For the corrector step (3.8), $\bm{Q_{c}}^{(i)}$ is solved directly from
(3.8b), and then obtain $\bm{E_{c}}^{(i)}$ and $\bm{B_{c}}^{(i)}$ iteratively
by

$\displaystyle\left(\frac{1}{\Delta
t}M_{h}-a_{ii}\alpha_{0}S_{h}+a_{ii}\Lambda_{h,\tilde{T}}^{(i)}\right)\bm{E_{c}}^{i,l+1}=a_{ii}\Lambda_{h,\tilde{T}}^{(i)}\bm{B_{c}}^{i,l}+F_{E}^{(i)},$
(3.12a)
$\displaystyle\bm{B_{c}}^{i,l+1}=(\bm{Q_{c}}^{(i)}-\bm{E_{c}}^{i,l+1})^{4}.$
(3.12b)

with the initial values $B_{h}^{i,0}=\tilde{B}_{h}^{(i)}$. The stop criteria
is the same as Step 3 in the predictor step.

For the Picard iteration in the predictor step, we can show that, if $E_{h}$
and $B_{h}\in W_{h}^{(0,0)}$, and the implicit part of the IMEX RK time
discretizations has a nonnegative diagonal, the resulting matrix
$A(\tilde{U}^{(i)})$ in (3.9) is an $\mathbf{M}$ matrix if
$\tilde{U}^{(i)}\geq 0$:

######  Theorem 1.

The matrix $A(\tilde{U}^{(i)})$ in (3.9) is an $\mathbf{M}$ matrix for
piecewise constant $\mathcal{Q}^{0}$ finite elements if $\tilde{U}^{(i)}\geq
0$.

We can easily show that the mass matrix $A(\tilde{U}^{(i)})$ has a nonnegative
diagonal, while off-diagonal entries are all non-positive. Besides,
$A(U^{(i)})$ is diagonally dominant, so $A(U^{(i)})$ is an $\mathbf{M}$
matrix, namely its inverse has all nonnegative arguments [4]. Thanks to this
property, rapid convergence can be assured. A similar fast convergence is also
observed in the iteration for the correction step (3.12), with
$\Lambda_{h,\tilde{T}}$ obtained from the predictor step and initial values
taken from $\tilde{B}$.

###  3.5 Positivity preserving and TVB limiters

For the solutions of (1.1), both the radiation energy $E({\bf x},t)$ and the
material temperature $T({\bf x},t)$ remain positive at all times [57].
Unfortunately, the schemes described above cannot preserve the positivity of
the solutions when they are close to zero. For solutions in $W_{h}^{(0,0)}$,
if the right-hand side terms $F_{E}^{(i)}$ and $F_{B}^{(i)}$ in (3.7) are non-
negative, the solutions without correction preserve positivity. However, this
requires a time step size on the order of $\mathcal{O}(h^{2})$. Furthermore,
an M-matrix cannot be guaranteed for higher-order schemes. In this work, to
develop an efficient scheme with larger time steps, we rely on the following
procedure to preserve positivity and control numerical oscillations.

Here we employ a positivity preserving limiter denoted as $P\Pi_{h}$ from
[59]. The $P\Pi_{h}$ limiter has also been applied to porous medium equations
in [49]. Additionally, TVB limiters [11] are utilized to control numerical
oscillations in the case of non-smooth solutions for second and third-order
schemes. As mentioned earlier, even for a first-order scheme with solutions in
$W^{0,0}_{h}$, the cell averages may not necessarily remain positive when
using large time step sizes. For diffusion systems, a cut-off limiter is
commonly employed to ensure a lower bound on the solutions [41, 51, 55].

For a one-dimensional case, we denote $\bar{u}_{j}$ as the cell average of the
numerical solution $u_{h}$ in the cell $I_{j}$ and $u_{\text{cold}}$ as the
minimum value of $u({\bf x},0)$. First, for the limiter $P\Pi_{h}$, it is
applied as follows:

  1. 1.

Check the cell average $\bar{u}_{j}$ in each cell $I_{j}$, if it is less than
$u_{\text{cold}}$, we set

$P\Pi_{h}u_{h}=u_{\text{cold }};$

  2. 2.

Then if $\bar{u}_{j}\geq u_{\text{cold}}$, but at least one value is less than
$u_{\text{cold}}$ appearing at two endpoints of the cell $I_{j}$ or at any
Gaussian point inside the cell $I_{j}$, we set
$P\Pi_{h}u_{h}=P\Pi^{1}_{h}u_{h}^{1}$, where $u_{h}^{1}$ is a local
$L^{2}$-projection of the solution $u_{h}$ to a linear polynomial within this
cell, where

$P\Pi^{1}_{h}u_{h}^{1}=\begin{cases}\left[1-2h_{j}^{-1}\left(x-x_{j}\right)\right]\bar{u}_{j}+2u_{\text{cold
}}h_{j}^{-1}\left(x-x_{j}\right),&\text{ if
}u_{h,j+\frac{1}{2}}^{1,-}<u_{\text{cold}},\\\ \,\\\
\left[1+2h_{j}^{-1}\left(x-x_{j}\right)\right]\bar{u}_{j}-2u_{\text{cold
}}h_{j}^{-1}\left(x-x_{j}\right),&\text{ if
}u_{h,j-\frac{1}{2}}^{1,+}<u_{\text{cold}};\end{cases}$

  3. 3.

Otherwise, we set $P\Pi_{h}u_{h}|_{I_{j}}=u_{h}|_{I_{j}}$.

After applying the PP limiter $P\Pi_{h}$, a generalized TVB limiter is further
applied to control numerical oscillations for second and third order schemes:

$\displaystyle\begin{cases}u_{j-\frac{1}{2}}^{+,(\bmod)}=\text{TVB minmod
}\left(\bar{u}_{j}-u_{j-\frac{1}{2}}^{+},\frac{\theta}{2}\left(\bar{u}_{j}-\bar{u}_{j-1}\right),\frac{\theta}{2}\left(\bar{u}_{j+1}-\bar{u}_{j}\right),\frac{1}{4}\left(\bar{u}_{j+1}-\bar{u}_{j-1}\right)\right),\\\
u_{j+\frac{1}{2}}^{-,(\bmod)}=\text{TVB minmod
}\left(u_{j+\frac{1}{2}}^{-}-\bar{u}_{j},\frac{\theta}{2}\left(\bar{u}_{j}-\bar{u}_{j-1}\right),\frac{\theta}{2}\left(\bar{u}_{j+1}-\bar{u}_{j}\right),\frac{1}{4}\left(\bar{u}_{j+1}-\bar{u}_{j-1}\right)\right),\end{cases}$

where

$\displaystyle\text{TVB minmod }(a,b,c,d)=\begin{cases}a,&\text{ if }|a|\leq
Mh^{2},\\\ \operatorname{minmod}(a,b,c,d),&\text{ otherwise, }\end{cases}$

and

$\displaystyle\operatorname{minmod}(a,b,c,d)=\begin{cases}\operatorname{sign}(a)\min(|a|,|b|,|c|),&\text{
if }ab>0,ac>0,\\\ 0,&\text{ otherwise. }\end{cases}$

If $u_{j-\frac{1}{2}}^{+,(\bmod)}$ is different from $u_{j-\frac{1}{2}}^{+}$,
or $u_{j+\frac{1}{2}}^{-,(\bmod)}$ is different from $u_{j+\frac{1}{2}}^{-}$,
the cell $I_{j}$ is identified as a troubled cell, and the polynomial in this
cell is modified to be

$u_{h}(x)=\bar{u}_{j}+\left(u_{j+\frac{1}{2}}^{-,(\bmod)}+u_{j-\frac{1}{2}}^{+,(\bmod)}\right)h_{j}^{-1}\left(x-x_{j}\right).$

In the generalized TVB limiter, we take the parameters $M=0.1$ and
$\theta=1.5$ in our numerical tests.

For the two-dimensional case, we denote $\bar{u}_{i,j}$ as the cell average of
the numerical solution $u_{h}$ in the cell $I_{ij}$ and $u_{\text{cold}}$ as
the minimum value of $u_{h}({\bf x},0)$.
$u_{i\pm\frac{1}{2},j\pm\frac{1}{2}}=u_{h}(x_{i\pm\frac{1}{2}},y_{j\pm\frac{1}{2}})$
are denoted the values of the endpoints within the cell $I_{i,j}$. The PP
limiter $P\Pi_{h}$ is applied as follows:

  1. 1.

Check the cell average $\bar{u}_{i,j}$ in each cell $I_{i,j}$, if it is less
than $u_{\text{cold}}$, we set

$P\Pi_{h}u_{h}=u_{\text{cold }};$

  2. 2.

Then if $\bar{u}_{i,j}\geq u_{\text{cold}}$, but at least one value is less
than $u_{\text{cold}}$ appearing at those endpoints of the cell $I_{ij}$ or at
any Gaussian point inside the cell $I_{ij}$, we set
$P\Pi_{h}u_{h}=P\Pi_{h}^{1}u_{h}^{1}$, where

$u_{h}^{1}=\bar{u}_{i,j}+\alpha_{1}\psi_{x}^{i}+\alpha_{2}\psi_{y}^{j},\quad\psi_{x}^{i}=\frac{x-x_{i}}{h_{x}^{i}},\,\psi_{y}^{j}=\frac{y-y_{j}}{h_{y}^{j}}.$

$u_{h}^{1}$ is a local $L^{2}$-projection of the solution $u_{h}$ to
$\mathcal{P}^{1}$ within this cell $I_{i,j}$. For $u^{1}_{h}$, the extremum
occurs only at four endpoints of the cell $I_{i,j}$, and due to
$\bar{u}_{i,j}\geq u_{\text{cold}}$ at most two adjacent points are less than
$u_{\text{cold}}$. We set those end-point values to be $u_{\text{cold}}$ if
they are less than $u_{\text{cold}}$. Next we only need to determine the
parameters $\alpha_{1}$, $\alpha_{2}$ by using two end-point values, and we
take two minimum point values at those end points. They are determined as
follows:

     * (a)

if only one end-point value at $(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}})$ is
modified, we have

$\displaystyle P\Pi_{h}^{1}u^{1}_{h}=$
$\displaystyle\bar{u}_{i,j}+\left(\bar{u}_{i,j}-u_{\text{cold}}+\frac{u^{1}_{i+\frac{1}{2},j-\frac{1}{2}}-u^{1}_{i-\frac{1}{2},j+\frac{1}{2}}}{2}\right)\psi_{x}^{i}$
$\displaystyle+\left(\bar{u}_{i,j}-u_{\text{cold}}-\frac{u^{1}_{i+\frac{1}{2},j-\frac{1}{2}}-u^{1}_{i-\frac{1}{2},j+\frac{1}{2}}}{2}\right)\psi_{y}^{j},$

     * (b)

otherwise if two end-point values at $(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}})$
and $(x_{i-\frac{1}{2}},y_{j+\frac{1}{2}})$ are changed, we take

$P\Pi_{h}^{1}u^{1}_{h}=\bar{u}_{i,j}+2(\bar{u}_{i,j}-u_{\text{cold}})\psi_{x}^{i}.$

Other cases can be determined similarly, we omit them to save space;

  3. 3.

Otherwise, we set $P\Pi_{h}u_{h}|_{I_{ij}}=u_{h}|_{I_{ij}}$.

For the two-dimensional generalized TVB limiter, we define
$w^{1}(y)=u_{h}(x_{i},y)$ in the cell $I_{i,j}$ and then modify $w^{1}(y)$
along $x=x_{i}$ as in the one-dimensional case, that is,

$\displaystyle\begin{cases}u_{j-\frac{1}{2}}^{+,(\bmod)}=\text{TVB
minmod}\left(\bar{u}_{i,j}-w^{1}(y_{j-\frac{1}{2}}^{+}),\frac{\theta}{2}\left(\bar{u}_{i,j}-\bar{u}_{i,j-1}\right),\frac{\theta}{2}\left(\bar{u}_{i,j+1}-\bar{u}_{i,j}\right),\frac{1}{4}\left(\bar{u}_{i,j+1}-\bar{u}_{i,j-1}\right)\right),\\\
u_{j+\frac{1}{2}}^{-,(\bmod)}=\text{TVB
minmod}\left(w^{1}(y_{j+\frac{1}{2}}^{-})-\bar{u}_{i,j},\frac{\theta}{2}\left(\bar{u}_{i,j}-\bar{u}_{i,j-1}\right),\frac{\theta}{2}\left(\bar{u}_{i,j+1}-\bar{u}_{i,j}\right),\frac{1}{4}\left(\bar{u}_{i,j+1}-\bar{u}_{i,j-1}\right)\right).\end{cases}$

If $u_{j-\frac{1}{2}}^{+,(\bmod)}$ is different from
$w^{1}(y_{j-\frac{1}{2}}^{+})$, or $u_{j+\frac{1}{2}}^{-,(\bmod)}$ is
different from $w^{1}(y_{j+\frac{1}{2}}^{-})$, the cell $I_{i,j}$ is
identified as a troubled cell along $x=x_{i}$ and $w^{1}(y)$ is modified to be

$w^{1}(y)=\bar{u}_{i,j}+\left(u_{j+\frac{1}{2}}^{-,(\bmod)}+u_{j-\frac{1}{2}}^{+,(\bmod)}\right)h_{j}^{-1}\left(y-y_{j}\right).$

Similarly, along $y=y_{j}$ we can also define $w^{2}(x)=u_{h}(x,y_{j})$, and

$\displaystyle\begin{cases}u_{i-\frac{1}{2}}^{+,(\bmod)}=\text{TVB
minmod}\left(\bar{u}_{i,j}-w^{2}(x_{i-\frac{1}{2}}^{+}),\frac{\theta}{2}\left(\bar{u}_{i,j}-\bar{u}_{i-1,j}\right),\frac{\theta}{2}\left(\bar{u}_{i+1,j}-\bar{u}_{i,j}\right),\frac{1}{4}\left(\bar{u}_{i+1,j}-\bar{u}_{i-1,j}\right)\right),\\\
u_{i+\frac{1}{2}}^{-,(\bmod)}=\text{TVB
minmod}\left(w^{2}(x_{i+\frac{1}{2}}^{-})-\bar{u}_{i,j},\frac{\theta}{2}\left(\bar{u}_{i,j}-\bar{u}_{i-1,j}\right),\frac{\theta}{2}\left(\bar{u}_{i+1,j}-\bar{u}_{i,j}\right),\frac{1}{4}\left(\bar{u}_{i+1,j}-\bar{u}_{i-1,j}\right)\right).\end{cases}$

If $u_{i-\frac{1}{2}}^{+,(\bmod)}\neq w^{2}(x_{i-\frac{1}{2}}^{+})$, or
$u_{i+\frac{1}{2}}^{-,(\bmod)}\neq w^{2}(x_{i+\frac{1}{2}}^{-})$, the cell
$I_{i,j}$ is identified as a troubled cell along $y=y_{j}$ and $w^{2}(x)$ is
modified to be

$w^{2}(x)=\bar{u}_{i,j}+\left(u_{i+\frac{1}{2}}^{-,(\bmod)}+u_{i-\frac{1}{2}}^{+,(\bmod)}\right)h_{i}^{-1}\left(x-x_{i}\right).$

If $I_{i,j}$ is a troubled cell along either the $x$ or the $y$ direction, we
set $u_{h}(x,y)$ in the cell $I_{i,j}$ as
$u_{h}(x,y)=w^{1}(y)+w^{2}(x)-\bar{u}_{i,j}$.

###  3.6 Algorithm flowchart

Finally we present our high order IMEX-LDG scheme updating procedure from time
level $t_{n}$ to $t_{n+1}$ in Algorithm 1.

Input: $E^{n}$, $T^{n}$

Output: $E^{n+1}$, $T^{n+1}$

1 $E^{(1)}=E^{n}$, $B^{(1)}=(T^{n})^{4}$;

2 for  _$i=2$ to $s$_ do

3 $E^{i,0}=E^{(i-1)}$, $B^{i,0}=B^{(i-1)}$, $l=0$;

4 $M_{3}=-a_{ii}Z_{h},\quad M_{4}=\frac{1}{\Delta
t}M_{h}+a_{ii}Z_{h}-a_{ii}\beta_{0}S_{h}$;

5 while  _$\|E^{i,l+1}-E^{i,l}\|\geq\delta$ or $l=0$_ do

6 $M_{1}=\frac{1}{\Delta
t}M_{h}-a_{ii}\alpha_{0}S_{h}+a_{ii}\Lambda_{h,T},\quad
M_{2}=-a_{ii}\Lambda_{h,T}$;

7 update $E^{i,l+1}$ and $B^{i,l+1}$ by (3.11);

8 apply the limiters to $E^{i,l+1}$ and $B^{i,l+1}$;

9 $l=l+1$;

10 end while

11 $\tilde{E}^{(i)}=E^{i,l+1}$, $\tilde{B}^{(i)}=B^{i,l+1}$,
$\tilde{T}^{(i)}=(\tilde{B}^{(i)})^{1/4}$ ;

12 obtain $Q^{(i)}$ by (3.8b);

13 $E^{i,0}=\tilde{E}^{(i)}$, $l=0$;

14 while  _$\|E^{i,l+1}-E^{i,l}\|\geq\delta$ or $l=0$_ do

15 update $E^{i,l+1}$ by (3.12a);

16 apply the limiters to $E^{i,l+1}$;

17 $B^{i,l+1}=(Q^{(i)}-E^{i,l+1})^{4}$;

18 apply the limiters to $B^{i,l+1}$;

19 $l=l+1$;

20

21 end while

22 $E^{(i)}=E^{i,l+1}$, $T^{(i)}=Q^{(i)}-E^{(i)}$;

23

24 end for

25$E^{n+1}=E^{(s)}$, $T^{n+1}=T^{(s)}$;

26

Algorithm 1 The updating from $t^{n}$ to $t^{n+1}$ with an $s$-stage
$(k+1)$-th order IMEX method

######  Remark 2.

(Numerical boundary treatment) There are mainly three types of boundary
conditions in our numerical experiments: periodic, Neumann, and mixed boundary
conditions. The numerical fluxes at the boundary under periodic or Neumann
boundary conditions are chosen as in [12].

Taking 1D with $\Omega=[x_{L},x_{R}]$ as an example, for a periodic boundary
condition, we take $\phi_{L}^{-}=\phi_{R}^{-}$, $\phi_{R}^{+}=\phi_{L}^{+}$
where $\phi=u$ or $u_{x}$. For a Neumann boundary condition, we set
$\hat{u}_{L}=u_{L}^{+}$, $\hat{u}_{R}=u_{R}^{-}$, while
$\hat{v}_{L}=\hat{v}_{R}=0$ where $v=u_{x}$. Here $\phi^{\pm}_{L}$,
$\phi^{\pm}_{R}$, are the left and right limits at $x_{L}$ and $x_{R}$
respectively. Similarly for $u^{+}_{L}$ and $u^{-}_{R}$.

For a mixed boundary condition, e.g. in the Marshak wave problem in section 4,
the physical boundary and initial conditions on the domain
$x\in\Omega=[x_{L},x_{R}]=[0,1]$ are given as

$\begin{cases}\dfrac{E}{4}-\dfrac{1}{6\sigma}\dfrac{\partial E}{\partial
x}=1,\quad\dfrac{\partial T}{\partial x}=0,\quad x=0,\\\\[5.0pt]
\dfrac{E}{4}+\dfrac{1}{6\sigma}\dfrac{\partial E}{\partial
x}=0,\quad\dfrac{\partial T}{\partial x}=0,\quad x=1,\\\\[5.0pt]
E|_{t=0}=1.0\times 10^{-5},\quad T|_{t=0}=E^{1/4}|_{t=0}.\end{cases}$ (3.13)

Taking the left boundary $x_{L}=x_{\frac{1}{2}}=0$ as an example, we omit the
artificial term $\frac{|\nabla E|}{E}$ in $D_{r}$ and the boundary condition
is simplified to [60, 19, 16]

$\dfrac{E}{4}-\dfrac{1}{2}D_{r}\frac{\partial E}{\partial x}=1.$

Hence, in both the predictor and corrector steps, we set the numerical fluxes
as

$\hat{E}_{h}=E_{h}^{+},\quad\hat{D}_{r,h}\hat{p}_{h}=\dfrac{\hat{E}_{h}}{2}-2.$

Similarly at the right boundary $x_{N+\frac{1}{2}}=1$, we take

$\hat{E}_{h}=E_{h}^{-},\quad\hat{D}_{r,h}\hat{p}_{h}=-\dfrac{\hat{E}_{h}}{2}.$

In the corrector step, we set

$\hat{Q}_{h}=\hat{E}_{h}+\hat{T}_{h},\quad\hat{w}_{h}=3\sigma(\hat{T}_{h})\hat{D}_{r,h}\hat{p}_{h}.$

##  4 Numerical Examples

In this section, we conduct numerical experiments to validate the high-order
accuracy, conservation properties, suitability for large time steps, and
effectiveness in capturing sharp fronts in both homogeneous and heterogeneous
media using our proposed schemes. For the added diffusion terms, we take the
coefficients as $\alpha_{0}=\tau\max\\{\frac{1}{3\sigma}\\}$,
$\beta_{0}=\tau\max\\{D_{t}\\}$ and $\gamma_{0}=\alpha_{0}+\beta_{0}$ with
$\tau=0.6$ [49]. In space, the $\mathcal{Q}^{\bm{k}}$ basis with $k$-th
piecewise polynomial in each direction for $k=0,1,2$ is taken.
Correspondingly, an $s$-stage $(k+1)$-th order globally stiffly accurate IMEX
RK time discretization is employed in time with the double Butcher tableau
given in A. In the following, our schemes are denoted as $(k+1)$-th order
methods for $k=0,1,2$ with $s=2,3,5$, respectively.

######  Example 4.1.

(Accuracy test in 1D) First we consider a 1D example with smooth initial
values and periodic boundary conditions at the equilibrium, which are given by

$\begin{cases}T(x,0)=0.8+0.1\sin(x),\\\\[5.0pt] E(x,0)=T(x,0)^{4},\end{cases}$
(4.1)

on the computational domain $[-\pi,\pi]$. We take $\kappa=0.1$ and a
homogeneous medium $z(x)=1$. The problem is run to time $t=5$ using the 1st,
2nd, and 3rd order methods, respectively. Since the exact solution is not
available, we compute the numerical errors by comparing numerical solutions
with a reference solution. In this case, the source term is not stiff as $T$
is away from $0$, so we compute the reference solution by a 3rd order LDG
method with a 3rd order explicit strong-stability-preserving RK time
discretization [15], on a much refined mesh $N=1024$. In Tables 4.1-4.3, we
show the numerical $L^{2}$ errors and orders of accuracy for the 1st, 2nd, and
3rd order schemes with different time steps, respectively. From these tables,
we can see that our methods achieve the corresponding orders of accuracy when
the time step is $\Delta t=\mathcal{O}(h).$ From the numerical results, we
find that larger $\Delta t$ lead to larger errors. In Table 4.3, for the 3rd
order method, an order reduction can be observed for a large ratio of $\Delta
t/h$. The order increases with further mesh refinements. Such a phenomeno may
be due to a high order IMEX time discretization.

Table 4.1: The numerical $L^{2}$ errors and orders of accuracy for the 1st order scheme with different time steps for Example 4.1. $t=5$.  N | $\Delta t$ |  $L^{2}$ error | order |  $L^{2}$ error | order |  $L^{2}$ error | order
---|---|---|---|---|---|---|---
E | T | B
4 | $\frac{h}{2}$ | 3.46e-2 | - | 1.74e-2 | - | 3.51e-2 | -
8 | 1.73e-2 | 1.00 | 8.65e-3 | 1.01 | 1.76e-2 | 0.99
16 | 8.68e-3 | 1.00 | 4.34e-3 | 0.99 | 8.84e-3 | 1.00
32 | 4.34e-3 | 1.00 | 2.17e-3 | 1.00 | 4.42e-3 | 1.00
64 | 2.17e-3 | 1.00 | 1.09e-3 | 1.00 | 2.21e-3 | 1.00
128 | 1.08e-3 | 1.00 | 5.43e-4 | 1.00 | 1.10e-3 | 1.00
4 | $h$ | 3.47e-2 | - | 1.75e-2 | - | 3.52e-2 | -
8 | 1.74e-2 | 1.00 | 8.68e-3 | 1.01 | 1.77e-2 | 0.99
16 | 8.70e-3 | 1.00 | 4.36e-3 | 1.00 | 8.86e-3 | 1.00
32 | 4.35e-3 | 1.00 | 2.18e-3 | 1.00 | 4.43e-3 | 1.00
64 | 2.17e-3 | 1.00 | 1.09e-3 | 1.00 | 2.21e-3 | 1.00
128 | 1.09e-3 | 1.00 | 5.44e-4 | 1.00 | 1.11e-3 | 1.00
4 | $3h$ | 3.52e-2 | - | 1.77e-2 | - | 3.57e-2 | -
8 | 1.77e-2 | 0.99 | 8.85e-3 | 1.00 | 1.80e-2 | 0.99
16 | 8.86e-3 | 1.00 | 4.44e-3 | 1.00 | 9.02e-3 | 1.00
32 | 4.43e-3 | 1.00 | 2.22e-3 | 1.00 | 4.51e-3 | 1.00
64 | 2.21e-3 | 1.00 | 1.11e-3 | 1.00 | 2.25e-3 | 1.00
128 | 1.11e-3 | 1.00 | 5.55e-4 | 1.00 | 1.13e-3 | 1.00
4 | $5h$ | 3.53e-2 | - | 1.78e-2 | - | 3.58e-2 | -
8 | 1.79e-2 | 0.98 | 9.00e-3 | 0.98 | 1.82e-2 | 0.97
16 | 9.04e-3 | 0.99 | 4.55e-3 | 0.98 | 9.21e-3 | 0.99
32 | 4.56e-3 | 0.99 | 2.30e-3 | 0.99 | 4.64e-3 | 0.99
64 | 2.28e-3 | 1.00 | 1.15e-3 | 1.00 | 2.32e-3 | 1.00
128 | 1.14e-3 | 1.00 | 5.74e-4 | 1.00 | 1.16e-3 | 1.00
Table 4.2: The numerical $L^{2}$ errors and orders of accuracy for the 2nd order scheme with different time steps for Example 4.1. $t=5.$ N | $\Delta t$ |  $L^{2}$ error | order |  $L^{2}$ error | order |  $L^{2}$ error | order
---|---|---|---|---|---|---|---
E | T | B
4 | $\frac{h}{2}$ | 1.25e-2 | - | 6.18e-3 | - | 1.27e-2 | -
8 | 3.00e-3 | 2.06 | 1.57e-3 | 1.98 | 3.04e-3 | 2.06
16 | 7.36e-4 | 2.02 | 3.86e-4 | 2.02 | 7.48e-4 | 2.02
32 | 1.83e-4 | 2.01 | 9.60e-5 | 2.01 | 1.86e-4 | 2.01
64 | 4.58e-5 | 2.00 | 2.40e-5 | 2.00 | 4.66e-5 | 2.00
128 | 1.14e-5 | 2.00 | 5.99e-6 | 2.00 | 1.16e-5 | 2.00
4 | $h$ | 1.25e-2 | - | 6.19e-3 | - | 1.27e-2 | -
8 | 3.00e-3 | 2.06 | 1.57e-3 | 1.98 | 3.04e-3 | 2.06
16 | 7.36e-4 | 2.02 | 3.86e-4 | 2.02 | 7.48e-4 | 2.02
32 | 1.83e-4 | 2.01 | 9.60e-5 | 2.01 | 1.86e-4 | 2.01
64 | 4.58e-5 | 2.00 | 2.40e-5 | 2.00 | 4.66e-5 | 2.00
128 | 1.14e-5 | 2.00 | 5.99e-6 | 2.00 | 1.16e-5 | 2.00
4 | $3h$ | 1.25e-2 | - | 6.29e-3 | - | 1.27e-2 | -
8 | 3.03e-3 | 2.05 | 1.60e-3 | 1.97 | 3.10e-3 | 2.04
16 | 7.41e-4 | 2.03 | 3.90e-4 | 2.04 | 7.58e-4 | 2.03
32 | 1.85e-4 | 2.01 | 9.63e-5 | 2.02 | 1.87e-4 | 2.02
64 | 4.61e-5 | 2.00 | 2.41e-5 | 2.00 | 4.68e-5 | 2.00
128 | 1.15e-5 | 2.00 | 6.03e-6 | 2.00 | 1.17e-5 | 2.00
4 | $5h$ | 1.30e-2 | - | 6.29e-3 | - | 1.27e-2 | -
8 | 3.18e-3 | 2.03 | 1.93e-3 | 1.71 | 3.81e-3 | 1.74
16 | 9.11e-4 | 1.80 | 4.37e-4 | 2.14 | 8.77e-4 | 2.12
32 | 2.00e-4 | 2.19 | 9.99e-5 | 2.13 | 1.93e-4 | 2.18
64 | 4.81e-5 | 2.05 | 2.52e-5 | 1.99 | 4.92e-5 | 1.97
128 | 1.20e-5 | 2.01 | 6.28e-6 | 2.00 | 1.22e-5 | 2.00
Table 4.3: The numerical $L^{2}$ errors and orders of accuracy for the 3rd order scheme with different time steps for Example 4.1. $t=5.$ N | $\Delta t$ |  $L^{2}$ error | order |  $L^{2}$ error | order |  $L^{2}$ error | order
---|---|---|---|---|---|---|---
E | T | B
4 | $\frac{h}{2}$ | 1.71e-3 | - | 1.03e-3 | - | 1.78e-3 | -
8 | 1.97e-4 | 3.12 | 1.14e-4 | 3.18 | 2.04e-4 | 3.12
16 | 2.49e-5 | 2.99 | 1.47e-5 | 2.96 | 2.61e-5 | 2.97
32 | 3.11e-6 | 3.00 | 1.86e-6 | 2.98 | 3.31e-6 | 2.98
64 | 3.89e-7 | 3.00 | 2.34e-7 | 2.99 | 4.15e-7 | 2.99
128 | 4.86e-8 | 3.00 | 2.93e-8 | 3.00 | 5.20e-8 | 3.00
4 | $h$ | 1.72e-3 | - | 1.04e-3 | - | 1.78e-3 | -
8 | 1.97e-4 | 3.12 | 1.14e-4 | 3.18 | 2.04e-4 | 3.12
16 | 2.49e-5 | 2.99 | 1.47e-5 | 2.96 | 2.62e-5 | 2.96
32 | 3.12e-6 | 3.00 | 1.87e-6 | 2.97 | 3.33e-6 | 2.98
64 | 3.91e-7 | 3.00 | 2.36e-7 | 2.98 | 4.21e-7 | 2.98
128 | 4.90e-8 | 3.00 | 2.97e-8 | 2.99 | 5.29e-8 | 2.99
4 | 3$h$ | 1.98e-3 | - | 1.15e-3 | - | 2.06e-3 | -
8 | 2.21e-4 | 3.17 | 1.28e-4 | 3.17 | 2.36e-4 | 3.12
16 | 2.99e-5 | 2.89 | 1.78e-5 | 2.85 | 3.30e-5 | 2.84
32 | 4.28e-6 | 2.80 | 2.81e-6 | 2.66 | 5.30e-6 | 2.64
64 | 7.37e-7 | 2.54 | 6.50e-7 | 2.11 | 1.31e-6 | 2.01
128 | 1.25e-7 | 2.56 | 1.07e-7 | 2.61 | 2.16e-7 | 2.60
256 | 1.97e-8 | 2.66 | 1.60e-8 | 2.74 | 3.23e-8 | 2.74
512 | 2.85e-9 | 2.79 | 2.21e-9 | 2.85 | 4.48e-9 | 2.85
4 | 5$h$ | 2.13e-3 | - | 1.19e-3 | - | 2.20e-3 | -
8 | 4.29e-4 | 2.31 | 2.37e-4 | 2.32 | 4.82e-4 | 2.19
16 | 6.43e-5 | 2.74 | 3.65e-5 | 2.70 | 7.39e-5 | 2.71
32 | 1.06e-5 | 2.61 | 8.27e-6 | 2.14 | 1.68e-5 | 2.13
64 | 1.94e-6 | 2.44 | 1.64e-6 | 2.33 | 3.31e-6 | 2.35
128 | 3.70e-7 | 2.39 | 3.09e-7 | 2.41 | 6.18e-7 | 2.42
256 | 6.82e-8 | 2.44 | 5.69e-8 | 2.44 | 1.15e-7 | 2.43
512 | 1.11e-8 | 2.61 | 8.99e-9 | 2.66 | 1.82e-8 | 2.66

######  Example 4.2.

Next we consider the following initial boundary data [29], with a sharp
transition in the initial values

$\begin{cases}E(x,0)=E_{L}+(E_{R}-E_{L})\dfrac{1+\tanh[50(x-0.25)]}{2},\\\\[5.0pt]
T(x,0)=E(x,0)^{1/4},\\\\[5.0pt] \dfrac{E}{4}-\dfrac{1}{6\sigma}\dfrac{\partial
E}{\partial x}=1,\,\,\,\,\dfrac{\partial T}{\partial x}=0,x=0,\\\\[5.0pt]
\dfrac{E}{4}+\dfrac{1}{6\sigma}\dfrac{\partial E}{\partial
x}=V_{R},\dfrac{\partial T}{\partial x}=0,x=1,\end{cases}$ (4.2)

where $V_{R}=1\times 10^{-3}$ , $E_{L}=4$ and $E_{R}=4\times 10^{-3}.$ We
consider two cases $\kappa=0$ and $\kappa=0.1$ in $D_{t}$. For the case of
$\kappa=0$, the system of $E$ and $B$ itself is in a conservative form. We
take $z(x)=1$ and a time step size $\Delta t=\frac{1}{5}h$. The results at
$t=1$ for $\kappa=0$ and $t=0.5$ for $\kappa=0.1$ are shown in Fig. 4.1 and
4.2, respectively. Reference solutions are obtained by a 1st order explicit RK
LDG method on $1024$ elements with a small enough time step, which is denoted
as “ref”. The left column is the radiation temperature $T_{r}$ and the right
column is the material temperature $T$. Two different mesh sizes are
considered, $N=64$ and $N=128$. “WL” refers to numerical solutions obtained
with limiters, and correspondingly “NL” refers to numerical solutions without
limiters. It can be observed that the higher the order, the closer the
numerical solutions are, as compared to the reference solutions. In addition,
the 1st order solutions perform well without limiters, but the 2nd order and
3rd order solutions have small oscillations in front of the sharp gradient.
After applying limiters, these oscillations can be well controlled. As we can
see from Fig. 4.1 and Fig. 4.2, refining the mesh from $N=64$ to $N=128$, the
numerical solutions match the reference solutions better, no matter with or
without limiters. Especially, deviations due to the application of the
limiters are also reduced with the mesh refinement.

Figure 4.1: The numerical results for Example 4.2 with $\kappa=0$. 1st, 2nd,
and 3rd order schemes are used. “WL”: with limiters, “NL”: without limiters.
“ref” represents the reference solutions. Top row: $N=64$; Bottom row:
$N=128$. Left column: radiation temperature $T_{r}$; Right column: material
temperature $T$. $t=1$, $\Delta t=\frac{h}{5}.$

Figure 4.2: The numerical results for Example 4.2 with $\kappa=0.1,$ and
$N=128$. 1st, 2nd, and 3rd order schemes are used. “WL” : with limiters, “NL”:
without limiters. “ref” represents the reference solutions. Top row: $N=64$;
Bottom row: $N=128$. Left column: radiation temperature $T_{r}$; Right column:
material temperature $T$. $t=0.5$, $\Delta t=\frac{h}{5}.$

######  Example 4.3.

We now consider an example with a periodic boundary condition to verify the
conservation errors. The initial values are taken to be

$\begin{cases}E(x,0)=E_{L}+(E_{R}-E_{L})\dfrac{1+\tanh[200(x-0.5)^{2}]}{2},\\\\[5.0pt]
T(x,0)=E(x,0)^{1/4},\end{cases}$ (4.3)

with $E_{L}=1,E_{R}=0.0001$, which are shown in the left column of Fig. 4.3.
We run the solution with $\kappa=0.5$ up to time $t=1$. We take $z(x)=1$ and
mesh numbers $N=128$. 1st, 2nd, and 3rd order schemes are used. “c” denotes
numerical results obtained with the conservation corrector step, while “nc” is
without the corrector step. As we have observed numerically, without the
corrector step, a smaller time step is needed for the convergence of
iteration. We take $\Delta t=\frac{1}{30}h$ for the 2nd method, and $\Delta
t=\frac{1}{50}h$ for the 3rd order method, without a corrector step, and
$\Delta t=\frac{1}{5}h$ for all others. In Fig. 4.3, on the right column, we
show the time evolution of conservation errors for the total energy $Q=E+T$.
We can observe that with a corrector step, the errors are much smaller than
those without a corrector step. In Fig. 4.4, we show the numerical solutions
for different orders with or without a corrector step. As we can see, with a
corrector step, all results match each other well. Without a corrector step,
we can clearly observe a deviation, especially for first and second order
methods. However, we would note that with a corrector step, the method itself
is conservative for the total energy, although limiters for second and third
order methods would slightly destroy such a conservation.

Figure 4.3: The numerical results for Example 4.3 with $\kappa=0.5$ and
$N=128$. 1st, 2nd, and 3rd order schemes are used. “c”: with a conservation
corrector step; “nc”: without a corrector step. Left: initial values; Right:
time evolution of conservation errors. $\Delta t=\frac{1}{30}h$ for “2nd, nc”,
$\Delta t=\frac{1}{50}h$ for “3rd, nc” and $\Delta t=\frac{1}{5}h$ for all
others.

Figure 4.4: The numerical results for Example 4.3 with $\kappa=0.5$ and
$N=128$. 1st, 2nd, and 3rd order schemes are used. “c”: with a conservation
corrector step; “nc”: without a corrector step. Left: radiation energy $E$;
Right: material temperature $T$. $\Delta t=\frac{1}{30}h$ for “2nd, nc”,
$\Delta t=\frac{1}{50}h$ for “3rd, nc” and $\Delta t=\frac{1}{5}h$ for all
others.

######  Example 4.4.

Here we consider the standard homogeneous Marshak wave problem [30, 45, 59]
with the initial and boundary conditions (3.13), where the atomic mass number
$z(x)=1$. This is a benchmark problem for radiation diffusion problems. It is
very challenging since the left inflow boundary value is not consistent with
the initial datum. Similarly, we consider both $\kappa=0$ and $\kappa=0.1$ in
$D_{t}$, respectively. The reference solutions are computed by a 1st order
explicit RK LDG method on $1024$ elements with a small enough time step. In
Fig. 4.5 and Fig. 4.6, we show the results at different times $t=1,2$ and
$t=1,1.5$, respectively. The left column is for the radiation temperature
$T_{r}$ and the right column is for the material temperature $T$. A mesh
refinement from $N=64$ to $N=128$ is also considered. As we can see that the
1st order numerical results are very deviated from the reference solutions,
second and third order methods capture the sharp fronts more accurately. With
mesh refinement, all results are getting closer to the reference solutions,
indicating that our methods are convergent with mesh refinement. The results
are also consistent with those in [33, 22]. As compared to [60, 16], a
relatively larger time step size $\Delta t=\frac{1}{5}h$ can be used, where
$\Delta t=\frac{4}{125}h$ with $h=\frac{1}{80}$ in [60] and $\Delta
t=\frac{1}{8}h^{2}$ with $h=\frac{1}{16}$ in [16]. In Fig 4.7, we compare the
results with and without a corrector step. In the case of $\kappa=0$, the
system $E$ and $B$ in (2.1) itself is in a conservative form, although it only
conserves the total energy $Q=E+T$ up to an error of numerical precision, the
results with or without a corrector step are almost the same. For
$\kappa=0.1$, the results with a corrector step clearly match the reference
solutions better, especially for first and second order methods. This has
demonstrated that the corrector step is very necessary.

######  Remark 3.

As discussed in Remark 1, if we eliminate the predictor step and directly
employ a Picard iteration in the corrector step, it may not work well for some
challenge problems. For this standard Marshak wave problem, if we consider
$E(x,0)=10^{-7}$, $N=64$ and $\Delta t=\frac{h}{5}$, even a first order scheme
does not converge well. Such an approach with higher orders work even worse.
However, our predictor-corrector procedure works well for these test cases.
With the results in the previous Example 4.3, we have shown that our methods
can ensure both conservation and robustness.

Figure 4.5: The numerical results of 1st, 2nd, and 3rd order schemes for
Example 4.4 with $\kappa=0,$ $\Delta t=\frac{h}{5},$ time $t=1$ and $t=2.$
“ref” represents the reference solutions. Top: $N=64$; Bottom: $N=128$. Left
column: radiation temperature $T_{r}$; Right column: material temperature $T$.

Figure 4.6: The numerical results of 1st, 2nd, and 3rd order schemes for
Example 4.4 with $\kappa=0.1,$ $\Delta t=\frac{h}{5},$ time $t=1$ and $t=1.5.$
“ref” represents the reference solutions. Top: $N=64$; Bottom: $N=128$. Left
column: radiation temperature $T_{r}$; Right column: material temperature $T$.

Figure 4.7: The numerical results of 1st, 2nd and 3rd order schemes for
Example 4.4 with $t=1,$ $\Delta t=\frac{h}{5},$ $N=64$. “ref” means the
reference solutions. Left column: radiation temperature $T_{r}$; Right column:
material temperature $T$. Top: $\kappa=0$; Bottom: $\kappa=0.1$. “c” and “nc”
denote numerical results with or without a corrector step, respectively.

######  Example 4.5.

In this example, we consider a heterogeneous Marshak problem [30, 45, 59] with
the initial and boundary conditions (3.13) with $z=3$ inside the interval
$[\frac{1}{3},\frac{2}{3}]$ and $z=1$ elsewhere. We take $\kappa=0.1$ and two
different times $t=1,3$. In Fig. 4.8, we show the results of 2nd and 3rd order
methods, on a mesh $N=180$ and the time step is taken to be $\Delta
t=\frac{1}{15}h$. The 1st order results are not presented since they deviate a
lot as shown above. From the numerical results, we can observe that the
radiation temperature $T_{r}$ and the material temperature propagate much
slower with a large value of $z$ in the middle region. Besides, the radiation
temperature $T_{r}$ is approaching the material temperature $T$ in the middle
region, reaching a thermodynamic equilibrium. In Fig. 4.9, we compare the
numerical results between the 2nd and 3rd order methods. “ref” represents the
reference solution obtained by the 1st order explicit RK LDG method on $1024$
elements with a small enough time step. We can find that the solutions of the
3rd order method match the reference better than the 2nd order ones,
especially for the radiation temperature. From this example, we can see that
high order methods can capture fine structures than corresponding lower order
methods.

Figure 4.8: The numerical results of 2nd and 3rd order schemes for Example 4.5
at time $t=1$ and $t=3$. $N=180$ and $\Delta t=\frac{h}{15}$. Left: 2nd order;
Right: 3rd order.

Figure 4.9: The comparison between the results of 2nd and 3rd order schemes
for Example 4.5 at time $t=1$ and $t=3$. Left: radiation temperature $T_{r}$;
Right: material temperature $T$. $N=180$ and $\Delta t=\frac{h}{15}$.

######  Example 4.6.

(Accuracy test in 2D) In this example we test the errors and convergence
orders of accuracy at the equilibrium in the 2D case. We consider a 2D
radiation diffusion problem with two given source terms $f_{1}(x,y,t)$ and
$f_{2}(x,y,t)$, which are chosen properly so that exact solutions are
available for the following system:

$\begin{cases}\dfrac{\partial E}{\partial t}-\nabla\cdot\left(D_{r}\nabla
E\right)=\sigma\left(T^{4}-E\right)+f_{1}(x,y,t),\\\\[5.0pt] \dfrac{\partial
T}{\partial t}-\nabla\cdot\left(D_{t}\nabla
T\right)=\sigma\left(E-T^{4}\right)+f_{2}(x,y,t),\end{cases}$ (4.4)

where
$T(x,y,t)=\left(0.8+0.1\sin\left(2\pi\left(x-t\right)\right)\right)\left(0.8+0.1\sin\left(2\pi\left(y-t\right)\right)\right)$
and $E(x,y,t)=T(x,y,t)^{4}$, $\Omega=[0,1]^{2}$. For convenience we remove the
artificially added term in the energy radiation coefficient, that is
$D_{r}=\frac{1}{3\sigma}$.

We take $\kappa=0.01$ and $z(x)=1$. The problem is run to time $t=0.5$ using
the 1st, 2nd, and 3rd order methods, respectively. We list the numerical
$L^{2}$ errors and orders of accuracy with different time step sizes in Table
4.4. Similar results as the 1D case are obtained.

Table 4.4: The numerical $L^{2}$ errors and orders of accuracy for the 1st, 2nd and 3rd order schemes with different time steps for Example 4.6. $t=0.5.$ | N | $\Delta t$ |  $L^{2}$ error | order |  $L^{2}$ error | order |  $L^{2}$ error | order
---|---|---|---|---|---|---|---|---
| E | T | B
1st | $16\times 16$ | $\frac{h}{2}$ | 1.22e-2 | - | 1.02e-2 | - | 1.23e-2 | -
$32\times 32$ | 6.21e-3 | 0.98 | 5.17e-3 | 0.98 | 6.23e-3 | 0.98
$64\times 64$ | 3.13e-3 | 0.99 | 2.60e-3 | 0.99 | 3.14e-3 | 0.99
$16\times 16$ | $h$ | 1.85e-2 | - | 1.80e-2 | - | 2.10e-2 | -
$32\times 32$ | 9.54e-3 | 0.96 | 9.35e-3 | 0.95 | 1.08e-2 | 0.95
$64\times 64$ | 4.84e-3 | 0.98 | 4.77e-3 | 0.97 | 5.50e-3 | 0.98
$16\times 16$ | $2h$ | 3.67e-2 | - | 4.07e-2 | - | 4.48e-2 | -
$32\times 32$ | 1.89e-2 | 0.96 | 2.14e-2 | 0.93 | 2.33e-2 | 0.94
$64\times 64$ | 9.65e-3 | 0.97 | 1.10e-2 | 0.96 | 1.19e-2 | 0.97
2nd | $16\times 16$ | $\frac{h}{2}$ | 1.11e-3 | - | 1.04e-3 | - | 1.14e-3 | -
$32\times 32$ | 2.78e-4 | 2.00 | 2.63e-4 | 1.99 | 2.90e-4 | 1.98
$64\times 64$ | 6.93e-5 | 2.00 | 6.57e-5 | 2.00 | 7.27e-5 | 2.00
$16\times 16$ | $h$ | 2.56e-3 | - | 2.89e-3 | - | 2.51e-3 | -
$32\times 32$ | 6.33e-4 | 2.01 | 7.45e-4 | 1.95 | 6.40e-4 | 1.97
$64\times 64$ | 1.57e-4 | 2.01 | 1.88e-4 | 1.99 | 1.61e-4 | 1.99
$16\times 16$ | $2h$ | 9.68e-3 | - | 1.01e-2 | - | 8.91e-3 | -
$32\times 32$ | 2.34e-3 | 2.05 | 2.79e-3 | 1.86 | 2.33e-3 | 1.94
$64\times 64$ | 5.81e-4 | 2.01 | 7.22e-4 | 1.95 | 5.93e-4 | 1.97
3rd | $16\times 16$ | $\frac{h}{2}$ | 7.63e-5 | - | 2.40e-5 | - | 2.59e-5 | -
$32\times 32$ | 1.13e-5 | 2.75 | 3.26e-6 | 2.88 | 3.31e-6 | 2.97
$64\times 64$ | 1.59e-6 | 2.83 | 4.32e-7 | 2.91 | 4.22e-7 | 2.97
$16\times 16$ | $h$ | 4.56e-4 | - | 1.11e-4 | - | 1.23e-4 | -
$32\times 32$ | 7.17e-5 | 2.67 | 1.62e-5 | 2.78 | 1.50e-5 | 3.04
$64\times 64$ | 1.08e-5 | 2.73 | 2.37e-6 | 2.78 | 1.95e-6 | 2.94
$16\times 16$ | $2h$ | 2.80e-3 | - | 8.52e-4 | - | 1.11e-3 | -
$32\times 32$ | 4.55e-4 | 2.62 | 1.10e-4 | 2.95 | 1.21e-4 | 3.20
$64\times 64$ | 7.16e-5 | 2.67 | 1.61e-5 | 2.77 | 1.47e-5 | 3.04

######  Example 4.7.

Now we consider a 2D blast wave problem (1.1) with initial and boundary
conditions given by

$\begin{cases}E(x,y,0)=10^{-3}+100e^{-100((x-1)^{2}+(y-1)^{2})},\\\
T(x,y,0)=E(x,y,0)^{\frac{1}{4}},\\\ \frac{\partial
E}{\partial\bm{n}}|_{\partial\Omega}=0,\quad\frac{\partial
T}{\partial\bm{n}}|_{\partial\Omega}=0,\end{cases}$ (4.5)

where $\Omega=[0,1]\times[0,1]$, $\kappa=0.01$.

First we consider a homogeneous case with $z(x,y)=1$. In Fig. 4.10 and Fig.
4.11, we show the numerical solutions $T_{r}$ and $T$ at time $t=0.5,1.5,2$ on
a mesh of $64\times 64$, respectively. We take a time step $\Delta
t=\frac{1}{5}h$. From top to bottom, numerical results obtained by the 1st,
2nd, and 3rd order methods are presented, respectively. We can observe that
our methods can capture the sharp fronts of the blast wave well. Moreover, in
Fig. 4.12, we show the cutting plots for the numerical solutions along $y=x+1$
at time $t=2$ on a mesh of $N\times N$ for $N=64,128,512$. We can see that the
three methods converge and match each other. 2nd and 3rd order methods are
clearly better than the 1st order method. From zoom-in figures, we can see the
results of 3rd order method are slightly better than the 2nd order ones.

######  Example 4.8.

Finally we consider the problems (1.1)-(4.5) in a heterogeneous medium. The
atomic mass number $z(x,y)$ is $1$ everywhere except in two inner square
regions $\frac{3}{16}<x<\frac{7}{16}$, $\frac{9}{16}<y<\frac{13}{16}$ and
$\frac{9}{16}<x<\frac{13}{16}$, $\frac{3}{16}<y<\frac{7}{16}$, where the value
of $z(x,y)$ is $3$, as shown in Fig. 4.13. In Fig. 4.14 and Fig. 4.15, we
present the numerical solutions $T_{r}$ and $T$ at time $t=1,2,2.5$ on a mesh
of $64\times 64$, respectively. As compared to $\Delta t=\frac{1}{32}h$ in
[16, 62], we use a lager time step $\Delta t=\frac{1}{10}h$. From top to
bottom, numerical results obtained by the 1st, 2nd, and 3rd order methods are
presented, respectively. From the results, we find that the results of the 1st
order method are very smeared due to numerical viscosities. 2nd and 3rd order
methods perform better. In Fig. 4.16, we show the cutting plots for the
numerical solutions along $y=x+1$ at $t=2$ on a mesh of $N\times N$ for
$N=64,128,512$. We can observe that by a mesh refinement, the results of 2nd
and 3rd methods match each other. In addition, we observe that when the front
of the radiation temperature propagates to the interface of two different
materials, the process is hindered by a region with dense atoms. Higher energy
exchange $\sigma=z({\bf x})^{3}/T^{3}$ occurs in these regions, resulting in a
well balance between the material and radiation temperatures. The temperature
propagates faster in the lower $z({\bf x})$ region. Our results agree with
those in [57, 55].

######  Remark 4.

We are not able to prove the convergence of the Picard iteration for (3.9)
theoretically. For the 2nd and 3rd order methods in space, the matrix
$A(U^{(i)})$ may not be diagonally dominant, so that $A(U^{(i)})$ cannot be
guaranteed to be an $\mathbf{M}$ matrix. However, in our numerical examples,
all results converge quickly under a tolerance $\delta=10^{-8}$. In Table 4.5,
we show an average number of iterations per stage for the nonlinear algebraic
system (3.9) in the heterogeneous medium. We can observe that in such
challenging cases, our methods with large time steps converge around four
iterations, which show the robustness of our proposed methods.

Table 4.5: Average iteration numbers per stage for some examples in the heterogeneous medium. |  | mesh |  time step |  an average number of iteration per stage
---|---|---|---|---
1D heterogeneous Marshak wave(Example 4.5) | 2nd | 180 | $\frac{1}{15}h$ | 4.02
3rd | 5.05
2D heterogeneous Blast wave(Example 4.8) | 1st | $64\times 64$ | $\frac{1}{10}h$ | 4.12
2nd | 3.99
3rd | 4.00

Figure 4.10: The numerical results for Example 4.7. From top to bottom: the
radiation temperature $T_{r}$ for 1st, 2nd, and 3rd order schemes. From left
to right: time $t=0.5,1.5,2.$ $\Delta t=\frac{1}{5}h$.

Figure 4.11: The numerical results for Example 4.7. From top to bottom: the
material temperature $T$ for 1st, 2nd, and 3rd order schemes. From left to
right: time $t=0.5,1.5,2.$ $\Delta t=\frac{1}{5}h$.

Figure 4.12: The numerical results for Example 4.7 at time $t=2$. Left: the
radiation temperature $T_{r}$; Right: the material temperature $T$. $\Delta
t=\frac{1}{5}h$.
$y$$x$01$\frac{3}{16}$$\frac{7}{16}$$\frac{9}{16}$$\frac{13}{16}$$\frac{3}{16}$$\frac{7}{16}$$\frac{9}{16}$$\frac{13}{16}$1313
Figure 4.13: The atomic mass number $z$ in the heterogeneous case for Example
4.8.

Figure 4.14: The numerical results for Example 4.8. From top to bottom: the
radiation temperature $T_{r}$ for 1st, 2nd, and 3rd order schemes. From left
to right: time $t=1.0,2.0,2.5.$ $\Delta t=\frac{1}{10}h$.

Figure 4.15: The numerical results for Example 4.8. From top to bottom: the
material temperature $T$ for 1st, 2nd, and 3rd order schemes. From left to
right: time $t=1.0,2.0,2.5.$ $\Delta t=\frac{1}{10}h$.

Figure 4.16: The numerical results for Example 4.8 at time $t=2$. Left: the
radiation temperature $T_{r}$; Right: the material temperature $T$. $\Delta
t=\frac{1}{10}h$.

##  5 Conclusion

In this paper, we have developed a class of high order conservative LDG-IMEX
methods for non-equilibrium radiation diffusion problems. The proposed scheme
is based on a predictor-corrector approach. During the predictor step, we
solve a reformulated system to improve convergence and approach thermodynamic
equilibrium. Subsequently, in the corrector step, we solve the original system
using nonlinear coefficients and initial estimates obtained from the predictor
step. This dual-step process ensures the conservation of total energy and
robust convergence. We incorporate linear diffusion terms to circumvent
implicit discretization for nonlinear diffusion terms. This approach permits
the use of larger time step sizes, in contrast to the restrictive parabolic
time step conditions $\Delta t=\mathcal{O}(h^{2})$ associated with purely
explicit discretizations. Local discontinuous Galerkin finite element
approximations in space are employed. Numerical examples in both 1D and 2D
illustrate the benefits of high-order conservative methods in accurately
capturing steep solution fronts within both homogeneous and heterogeneous
media. Our proposed methods exhibit robustness, as confirmed by numerical
experiments. However, theoretical analysis to ensure the convergence of such a
nonlinear iteration is a challenging task, which we plan to explore in our
future work. The extension to a 3D radiation diffusion model [25] or three-
temperature (3T) model [44, 56], and other reaction-diffusion equations [42]
will also be investigated.

##  Appendix A IMEX Butcher tableau

In this paper, the double Butcher tableaux we used from [3] are listed below:

######  Tableau A.1.

2 stage 1st order:

0 | 0 | 0 |  | 0 | 0 | 0
---|---|---|---|---|---|---
1 | 0 | 1 |  | 1 | 1 | 0
| 0 | 1 |  |  | 1 | 0

.

######  Tableau A.2.

3 stage 2nd order:

0 | 0 | 0 | 0 |  | 0 | 0 | 0 | 0
---|---|---|---|---|---|---|---|---
$\gamma$ | 0 | $\gamma$ | 0 |  | $\gamma$ | $\gamma$ | 0 | 0
1 | $\gamma$ | 0 | $\gamma$ |  | 1 | 0 | 1 | 0
| $\gamma$ | 0 | $\gamma$ |  |  | 0 | 1 | 0

,

where $\gamma=\frac{1}{2}$.

######  Tableau A.3.

5 stage 3rd order:

0 | 0 | 0 | 0 | 0 | 0 |  | 0 | 0 | 0 | 0 | 0 | 0
---|---|---|---|---|---|---|---|---|---|---|---|---
1/2 | 0 | 1/2 | 0 | 0 | 0 |  | 1/2 | 1/2 | 0 | 0 | 0 | 0
2/3 | 0 | 1/6 | 1/2 | 0 | 0 |  | 2/3 | 11/18 | 1/18 | 0 | 0 | 0
1/2 | 0 | -1/2 | 1/2 | 1/2 | 0 |  | 1/2 | 5/6 | -5/6 | 1/2 | 0 | 0
1 | 0 | 3/2 | -3/2 | 1/2 | 1/2 |  | 1 | 1/4 | 7/4 | 3/4 | -7/4 | 0
| 0 | 3/2 | -3/2 | 1/2 | 1/2 |  |  | 1/4 | 7/4 | 3/4 | -7/4 | 0

.

## References

  * [1] T. Arbogast, C.-S. Huang, and X. Zhao.  Finite volume WENO schemes for nonlinear parabolic problems with degenerate diffusion on non-uniform meshes.  Journal of Computational Physics, 399:108921, 2019.
  * [2] U. M. Ascher, S. J. Ruuth, and B. T. Wetton.  Implicit-explicit methods for time-dependent partial differential equations.  SIAM Journal on Numerical Analysis, 32:797–823, 1995.
  * [3] U.M. Ascher, S.J. Ruuth, and R.J. Spiteri.  Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations.  Applied Numerical Mathematics, 25(2-3):151–167, 1997.
  * [4] A. Berman and R.J. Plemmons.  Nonnegative Matrices in the Mathematical Sciences.  SIAM, 1994.
  * [5] Marianne Bessemoulin-Chatard and Francis Filbet.  A finite volume scheme for nonlinear degenerate parabolic equations.  SIAM Journal on Scientific Computing, 34(5):B559–B583, 2012.
  * [6] S. Bingjing and G.L. Olson.  Benchmark results for the non-equilibrium Marshak diffusion problem.  Journal of Quantitative Spectroscopy and Radiative Transfer, 56(3):337–351, 1996.
  * [7] S. Boscarino, L. Pareschi, and G. Russo.  Implicit-explicit Runge–Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit.  SIAM Journal on Scientific Computing, 35(1):A22–A51, 2013.
  * [8] R.L. Bowers and J.R. Wilson.  Numerical Modeling in Applied Physics and Astrophysics.  Jones and Bartlett Publishers, Inc., 1991.
  * [9] P.N. Brown, D.E. Shumaker, and C.S. Woodward.  Fully implicit solution of large-scale non-equilibrium radiation diffusion with high order time integration.  Journal of Computational Physics, 204(2):760–783, 2005.
  * [10] P. Castillo, B. Cockburn, D. Schötzau, and C. Schwab.  Optimal a priori error estimates for the $hp$-version of the local discontinuous Galerkin method for convection-diffusion problems.  Mathematics of Computation, 71(238):455–478, 2002.
  * [11] B. Cockburn and C.-W. Shu.  TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework.  Mathematics of Computation, 52(186):411–435, 1989.
  * [12] B. Cockburn and C.-W. Shu.  The local discontinuous Galerkin method for time-dependent convection-diffusion systems.  SIAM Journal on Numerical Analysis, 35(6):2440–2463, 1998.
  * [13] B. Cockburn and C.-W. Shu.  Runge-Kutta discontinuous Galerkin methods for convection-dominated problems.  Journal of Scientific Computing, 16(3):173–261, 2001.
  * [14] Irving R. Epstein and John A. Pojman.  An introduction to nonlinear chemical dynamics.  Oxford Academic, 1998.
  * [15] S. Gottlieb and C.-W. Shu.  Total variation diminishing Runge-Kutta schemes.  Mathematics of Computation, 67(221):73–85, 1998.
  * [16] Z. Huang and Y. Li.  Monotone finite point method for non-equilibrium radiation diffusion equations.  BIT Numerical Mathematics, 56(2):659–679, 2016.
  * [17] J. Jang, F. Li, J.-M. Qiu, and T. Xiong.  Analysis of asymptotic preserving DG-IMEX schemes for linear kinetic transport equations in a diffusive scaling.  SIAM Journal on Numerical Analysis, 52(4):2048–2072, 2014.
  * [18] A.M. Jokisaari, P.W. Voorhees, J.E. Guyer, J. Warren, and O.G. Heinonen.  Benchmark problems for numerical implementations of phase field models.  Computational Materials Science, 126:139–151, 2017.
  * [19] K.S. Kang.  ${P}_{1}$ Nonconforming Finite Element Multigrid Method for Radiation Transport.  SIAM Journal on Scientific Computing, 25(2):369–384, 2003.
  * [20] D.I. Ketcheson, C.B. Macdonald, and S. Gottlieb.  Optimal implicit strong stability preserving Runge–Kutta methods.  Applied Numerical Mathematics, 59(2):373–392, 2009.
  * [21] D.A. Knoll, L. Chacon, L.G. Margolin, and V.A. Mousseau.  On balanced approximations for time integration of multiple time scale systems.  Journal of Computational Physics, 185(2):583–611, 2003.
  * [22] D.A. Knoll, R.B. Lowrie, and J.E. Morel.  Numerical analysis of time integration errors for nonequilibrium radiation diffusion.  Journal of Computational Physics, 226(2):1332–1347, 2007.
  * [23] D.A. Knoll, W.J. Rider, and G.L. Olson.  An efficient nonlinear solution method for non-equilibrium radiation diffusion.  Journal of Quantitative Spectroscopy and Radiative Transfer, 63(1):15–29, 1999.
  * [24] D.A. Knoll, W.J. Rider, and G.L. Olson.  Nonlinear convergence, accuracy, and time step control in nonequilibrium radiation diffusion.  Journal of Quantitative Spectroscopy and Radiative Transfer, 70(1):25–36, 2001.
  * [25] X. Lai, Z. Sheng, and G. Yuan.  Monotone finite volume scheme for three dimensional diffusion equation on tetrahedral meshes.  Communications in Computational Physics, 21(1):162–181, 2017.
  * [26] E.W. Larsen, G.C. Pomraning, and V.C. Badham.  Asymptotic analysis of radiative transfer problems.  Journal of Quantitative Spectroscopy and Radiative Transfer, 29:285–310, 1983.
  * [27] Shu Liu, Siting Liu, Stanley Osher, and Wuchen Li.  A first-order computational algorithm for reaction-diffusion type equations via primal-dual hybrid gradient method.  Journal of Computational Physics, 500:112753, 2024.
  * [28] Yuanyuan Liu, Chi-Wang Shu, and Mengping Zhang.  High order finite difference WENO schemes for nonlinear degenerate parabolic equations.  SIAM Journal on Scientific Computing, 33(2):939–965, 2011.
  * [29] R.B. Lowrie.  A comparison of implicit time integration methods for nonlinear relaxation and diffusion.  Journal of Computational Physics, 196(2):566–590, 2004.
  * [30] R.E. Marshak.  Effect of radiation on shock wave behavior.  The Physics of Fluids, 1(1):24–29, 1958.
  * [31] V.A. Mousseau and D.A. Knoll.  New physics-based preconditioning of implicit methods for non-equilibrium radiation diffusion.  Journal of Computational Physics, 190(1):42–51, 2003.
  * [32] V.A. Mousseau and D.A. Knoll.  Temporal accuracy of the nonequilibrium radiation diffusion equations applied to two-dimensional multimaterial simulations.  Nuclear Science and Engineering, 154(2):174–189, 2006.
  * [33] V.A. Mousseau, D.A. Knoll, and W.J. Rider.  Physics-based preconditioning and the Newton–Krylov method for non-equilibrium radiation diffusion.  Journal of Computational Physics, 160(2):743–765, 2000.
  * [34] J. D. Murray.  Mathematical Biology II: Spatial Models and Biomedical Applications, volume 3rd Edition.  Springer New York, 2003.
  * [35] C.C. Ober and J.N. Shadid.  Studies on the accuracy of time-integration methods for the radiation–diffusion equations.  Journal of Computational Physics, 195(2):743–772, 2004.
  * [36] D.L. Peterson, R.L. Bowers, J.H. Brownell, A.E. Greene, K.D. McLenithan, T.A. Oliphant, N.F. Roderick, and A.J. Scannapieco.  Two-dimensional modeling of magnetically driven Rayleigh-Taylor instabilities in cylindrical Z pinches.  Physics of Plasmas, 3(1):368–381, 1996.
  * [37] G.C. Pomraning.  The non-equilibrium Marshak wave problem.  Journal of Quantitative Spectroscopy and Radiative Transfer, 21(3):249–261, 1979.
  * [38] A.C. Robinson and C.J. Garasi.  Three-dimensional z-pinch wire array modeling with ALEGRA-HEDP.  Computer Physics Communications, 164(1-3):408–413, 2004.
  * [39] Y. Saad and M.H. Schultz.  GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems.  SIAM Journal on scientific and Statistical Computing, 7(3):856–869, 1986.
  * [40] Z. Sheng, J. Yue, and G. Yuan.  Monotone finite volume schemes of nonequilibrium radiation diffusion equations on distorted meshes.  SIAM Journal on Scientific Computing, 31(4):2915–2934, 2009.
  * [41] Y. Shi, H. Yong, C. Zhai, J. Qi, and P. Song.  A functional expansion tally method for gray radiative transfer equations in implicit Monte Carlo.  Journal of Computational and Theoretical Transport, 47(7):581–598, 2018.
  * [42] J. Smoller.  Shock Waves and Reaction—Diffusion Equations, volume 258.  Springer Science & Business Media, 2012.
  * [43] L. Spitzer Jr and R. Härm.  Transport phenomena in a completely ionized gas.  Physical Review, 89(5):977, 1953.
  * [44] S. Su and J. Wu.  A vertex-centered and positivity-preserving finite volume scheme for two-dimensional three-temperature radiation diffusion equations on general polygonal meshes.  Numerical Mathematics: Theory, Methods and Applications, 13(1):220–252, 2020.
  * [45] M. Tang, L. Wang, and X. Zhang.  Accurate front capturing asymptotic preserving scheme for nonlinear gray radiative transfer equation.  SIAM Journal on Scientific Computing, 43(3):B759–B783, 2021.
  * [46] M. Tang and X. Zhang.  Semi-implicit front capturing schemes for the degenerate nonlinear radiative diffusion equation.  Journal of Computational Physics, 436:110290, 2021.
  * [47] N.J. Turner and J.M. Stone.  A module for radiation hydrodynamic calculations with ZEUS-2D using flux-limited diffusion.  The Astrophysical Journal Supplement Series, 135(1):95, 2001.
  * [48] H. Wang, C.-W. Shu, and Q. Zhang.  Stability and error estimates of local discontinuous Galerkin methods with implicit-explicit time-marching for advection-diffusion problems.  SIAM Journal on Numerical Analysis, 53(1):206–227, 2015.
  * [49] H. Wang, Q. Zhang, S. Wang, and C.-W. Shu.  Local discontinuous Galerkin methods with explicit-implicit-null time discretizations for solving nonlinear diffusion problems.  Science China Mathematics, 63(1):183–204, 2020.
  * [50] A.M. Winslow.  Multifrequency-gray method for radiation diffusion with Compton scattering.  Journal of Computational Physics, 117(2):262–273, 1995.
  * [51] T. Xiong, W. Sun, Y. Shi, and P. Song.  High order asymptotic preserving discontinuous Galerkin methods for gray radiative transfer equations.  Journal of Computational Physics, page 111308, 2022.
  * [52] Y. Xu and C.-W. Shu.  Local discontinuous Galerkin methods for high-order time-dependent partial differential equations.  Communications in Computational Physics, 7(1):1, 2010.
  * [53] J. Yan and C.-W. Shu.  Local discontinuous Galerkin methods for partial differential equations with higher order derivatives.  Journal of Scientific Computing, 17(1):27–47, 2002.
  * [54] X. Yang, W. Huang, and J. Qiu.  A moving mesh finite difference method for equilibrium radiation diffusion equations.  Journal of Computational Physics, 298:661–677, 2015.
  * [55] X. Yang, W. Huang, and J. Qiu.  Moving mesh finite difference solution of non-equilibrium radiation diffusion equations.  Numerical Algorithms, 82(4):1409–1440, 2019.
  * [56] Y. Yu, X. Chen, and G. Yuan.  A finite volume scheme preserving maximum principle for the system of radiation diffusion equations with three-temperature.  SIAM Journal on Scientific Computing, 41(1):B93–B113, 2019.
  * [57] G. Yuan, X. Hang, Z. Sheng, and J. Yue.  Progress in numerical methods for radiation diffusion equations.  Chinese Journal of Computational Physics, 26(4):475, 2009.
  * [58] Peng Zhang and Tao Xiong.  High order implicit finite difference schemes with a semi-implicit WENO reconstruction for nonlinear degenerate parabolic equations.  Journal of Computational Physics, 467:111442, 2022.
  * [59] Q. Zhang and Z.-L. Wu.  Numerical simulation for porous medium equation by local discontinuous Galerkin finite element method.  Journal of Scientific Computing, 38(2):127–148, 2009.
  * [60] R. Zhang, X. Yu, X. Cui, and T. Feng.  Discontinuous finite element method for 1D non-equilibrium radiation diffusion equations.  Chinese Journal of Computational Physics, 29(5):641, 2012.
  * [61] Rongpei Zhang, Xijun Yu, Jiang Zhu, and Abimael F.D. Loula.  Direct discontinuous galerkin method for nonlinear reaction–diffusion systems in pattern formation.  Applied Mathematical Modelling, 38:1612–1621, 2014.
  * [62] X. Zhao, Y. Chen, Y. Gao, C. Yu, and Y. Li.  Finite volume element methods for nonequilibrium radiation diffusion equations.  International Journal for Numerical Methods in Fluids, 73(12):1059–1080, 2013.
  * [63] Jianfeng Zhu, Yong-Tao Zhang, Stuart A. Newman, and Mark Alber.  Application of discontinuous galerkin methods for reaction-diffusion systems in developmental biology.  Journal of Scientific Computing, 40:391–418, 2009.