#  An asymptotic-preserving method for the three-temperature radiative
transfer model

Ruo Li111CAPT, LMAM & School of Mathematical Sciences, Peking University,
Beijing, China, email<EMAIL_ADDRESS>Weiming Li222Institute of Applied
Physics and Computational Mathematics, Beijing, China,
<EMAIL_ADDRESS>Shengtong Liang333School of Mathematical
Sciences, Peking University, Beijing, China, email<EMAIL_ADDRESS>Yuehan Shao444School of Mathematical Sciences, Peking University, Beijing,
China, email<EMAIL_ADDRESS>Min Tang555Institute of Natural Sciences,
School of Mathematics, Shanghai Jiaotong University, Shanghai, China, email:
<EMAIL_ADDRESS>Yanli Wang666Beijing Computational Science Research
Center, Beijing, China, email<EMAIL_ADDRESS>
###### Abstract

We present an asymptotic-preserving (AP) numerical method for solving the
three-temperature radiative transfer model, which holds significant importance
in inertial confinement fusion. A carefully designed splitting method is
developed that can provide a general framework of extending AP schemes for the
gray radiative transport equation to the more complex three-temperature
radiative transfer model. The proposed scheme captures two important limiting
models: the three-temperature radiation diffusion equation (3TRDE) when
opacity approaches infinity and the two-temperature limit when the ion-
electron coupling coefficient goes to infinity. We have rigorously
demonstrated the AP property and energy conservation characteristics of the
proposed scheme and its efficiency has been validated through a series of
benchmark tests in the numerical part.

###### keywords:

radiative transfer equation; three-temperature; explicit; implicit;
asymptotic-preserving

##  1 Introduction

The radiative transfer equation holds significant importance in inertial
confinement fusion, and accurately modeling thermal radiative transfer
presents a formidable challenge. In many radiative transfer scenarios, the
timescale required for electron and ion temperatures to approach thermal
equilibrium exceeds the timescale needed to achieve quasi-neutral states, as
noted by Enaux et al. in [6]. Consequently, it becomes necessary to employ a
three-temperature system, wherein radiation, ion, and electron energies are
considered separately. For example, in the study conducted by Shiva et al.
[25], a three-temperature radiation hydrodynamic model is employed to show the
influence of electron thermal radiation on laser ablative shock waves
propagating through ambient atmospheric air.

Numerical simulations of the three-temperature radiation transport model
(3TRTM) pose substantial challenges, as highlighted by Evans et al. in [8].
These challenges arise primarily due to the following reasons: 1) The
intrinsic nonlinearity of the problem: The radiation intensity exhibits a
nonlinear dependence on electron temperature, and the governing equations for
electrons and ions manifest nonlinearity in diffusion and source terms. 2)
High dimensionality of radiation intensity: The radiation intensity is a
function of seven independent variables. 3) Restrictive time step constraints:
The speed of light imposes stringent limitations on the time step length
during simulations, as discussed by Enaux et al. in [6]. 4) Multiscale nature
of physical parameters: The physical parameters, including opacity and the
ion-electron coupling coefficient, exhibit significant variability across
several orders of magnitude within the computational domain of certain
benchmark tests.

The 3TRTM has two important limits: the three-temperature radiation diffusion
equation (3TRDE) emerges as opacity approaches infinity, and the two-
temperature limit arises when the ion-electron coupling coefficient tends
towards infinity. The 3TRDE depends solely on space and time, resulting in a
reduction in the dimension of independent variables. Nevertheless, it inherits
the same challenges owing to nonlinearity and the multiscale ion-electron
coupling coefficient. In the two-temperature limit, where ion and electron
temperatures coincide, the model shares a similar form as the gray radiation
transport model (GRTM), which only considers radiation-electron interactions.

Both of these limiting models, while simpler than the 3TRTM, present
formidable challenges and have gathered significant attention from
researchers. The 3TRDE comprises three reaction-diffusion-type equations with
nonlinear diffusion and stiff nonlinear reactions. Efforts are made to develop
robust and efficient schemes capable of handling nonlinearity in arbitrary
geometries. For example, a positivity-preserving and conservative scheme is
proposed in [5, 22] and the discretization on unstructured grids can be found
in [26, 27]. The pursuit of maximum principle-preserving schemes is discussed
in [32], while semi-analytic solutions for the 3TRDE are explored in [20]. To
enhance computational efficiency, adaptation, and multigrid methods [14, 21],
and certain iterative nonlinear solvers [1, 2, 7, 29] have been proposed.

GRTM, being a high-dimensional system, faces similar challenges as the 3TRTM.
Numerical methods for solving GRTM have been extensively studied and usually
fall into two categories: stochastic and deterministic methods. The widely
adopted stochastic method is the implicit Monte Carlo method (IMC), originally
introduced for the two-temperature model by Fleck and Cummings (1971) and
subsequently improved in [23]. It has also been extended to the three-
temperature model in [8, 24]. However, the IMC method is limited by the
statistical noise, and its ability to handle multiscale parameters for GRTE
and 3TRTM has not been thoroughly explored. An AP scheme means that the
discrete numerical scheme of the microscopic model will converge to a
discretization of the macroscopic model as the scale parameter approaches the
limit. To handle multiscale parameters and the high speed of light, numerous
AP schemes for GRTE have been proposed in the literature. For instance, in
[18], McClarren et al. designed a semi-implicit time integration method that
explicitly treats the streaming term while handling the material coupling term
implicitly with a linearized emission source. The unified gas kinetic scheme,
as developed in [28], extends the unified gas kinetic scheme for the linear
radiative transport equation to the nonlinear GRTE by employing a linearized
iterative solver for the material thermal temperature. An AP scheme based on a
prediction-correction approach that allows for larger time steps and can
capture the front position of the degenerate case is proposed in [30]. A high-
order AP scheme for GRTE is detailed in [4]. However, the literature on
efficient AP solvers for 3TRTM is much less.

In this work, we introduce an AP splitting method to address the complexities
of the 3TRTM. A carefully designed splitting method is developed that can
provide a general framework for extending AP schemes for the GRTE to the more
complex 3TRTM. The main idea involves dividing the nonlinear interactions
between radiation, electrons, and ions temperatures in a specific way, such
that the resulting split systems can effectively capture the diffusion limit
when opacity tends towards infinity and the two-temperature limits when the
collision time coefficient approaches infinity. More precisely, the original
3TRTM is divided into two parts. The first part consists of a seven-
dimensional equation governing radiation density flux, featuring a modified
source term, along with two ordinary differential equations (ODEs) for
electron and ion temperatures. The second part is composed of three equations
governing radiation density, electron, and ion temperatures, with all three
functions being macroscopic, depending solely on space and time. To relax the
time step requirement, a fully implicit solver is employed for both parts,
solving nonlinear systems through carefully designed alternating iterations.
While the first system is characterized by high dimensionality and
nonlinearity, it is almost the same as the GRTE and can be effectively
addressed by extending the efficient AP solvers proposed in the existing
literature for GRTE, as discussed in the preceding paragraph. In this paper,
we extend the AP scheme introduced in [10] for GRTE to 3TRTM and apply it to
the alternating iteration method, eliminating the need to solve a large
nonlinear system. On the other hand, the structure of the second system shows
similarity to the 3TRDE, making it both easier to solve and the various robust
and efficient existing schemes [30] can be utilized to solve it. For
simplicity, we employ a fully implicit solver. The nonlinear system obtained
by implicit discretization is also solved using the alternating iteration
method. We subsequently establish the AP properties in two distinct limits:
the 3TRDE and GRTE, and validate the energy conservation of the numerical
scheme.

The rest of this paper is organized as follows. In Sec. 2, the details of
3TRTM and the two asymptotic limits are discussed. The AP temporal splitting
method is proposed in Sec. 3. The nonlinear iterative solvers for both systems
are discussed in Sec. 4 with particular emphasis on the first high-dimensional
system. The spatial discretization together with a discussion of the AP
properties is shown in Sec. 5. To validate the performance of the scheme,
numerical simulations are conducted in Sec. 6, including tests of the AP
properties and tests on the homogeneous model problem, the Marshak wave
problem with and without conduction terms, and a two-dimensional Riemann
problem. Finally, some conclusion remarks are made in Sec. 7.

##  2 Three-temperature radiative transfer model

In this section, we will present the detailed formulation of 3TRTM and focus
on the derivation of the two limiting models: 3TRDE and the two-temperature
limit.

###  2.1 The 3TRTM

In 3TRTM, the temperatures for ions ($T_{i}$), electrons ($T_{e}$), and the
radiation field ($\psi$) are coupled together to represent the evolution of
the radiation field, electron energy, and ion energy density in a hot plasma.
The detailed model after nondimensionalization is as follows [8].

$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial\psi}{\partial
t}+\epsilon\Omega\cdot\nabla\psi=\frac{1}{|\Omega|}\sigma
acT_{e}^{4}-\sigma\psi+\frac{Q_{r}}{|\Omega|},$ (2.1a)
$\displaystyle\frac{\epsilon^{2}\partial E_{e}}{\partial
t}=\epsilon^{2}\nabla\cdot D_{e}\nabla
T_{e}+c\kappa(T_{i}-T_{e})+\sigma\int_{\mathbb{S}^{2}}\psi\mathop{}\\!\mathrm{d}\Omega-\sigma
acT_{e}^{4}+Q_{e},$ (2.1b) $\displaystyle\frac{\epsilon^{2}\partial
E_{i}}{\partial t}=\epsilon^{2}\nabla\cdot D_{i}\nabla
T_{i}+c\kappa(T_{e}-T_{i})+Q_{i}.$ (2.1c)

Here $\psi(t,{\bm{x}},\Omega)$ represents the radiation intensity, where
$\Omega$ denotes the angular variable located on the unit sphere
$\mathbb{S}^{2}$, and $|\Omega|$ is the surface area of $\Omega$, equaling
$4\pi$. The dimensionless parameter $\epsilon$ is analogous to the Knudsen
number in the Boltzmann equation, representing the ratio of the mean free path
to the characteristic length. $\sigma({\bm{x}},T_{e})$ denotes the opacity, a
material-specific coefficient. $c$ stands for the speed of light. The
radiation constant $a$ is defined as

$a=\frac{8\pi k^{4}}{15h^{3}c^{3}},$ (2.2)

where $h$ represents Planck’s constant, and $k$ denotes the Boltzmann
constant. $E_{e}(t,{\bm{x}})$ and $T_{e}(t,{\bm{x}})$ refer to the energy and
temperature of the electrons, while $E_{i}(t,{\bm{x}})$ and
$T_{i}(t,{\bm{x}})$ represent the energy and temperature of the ions. The
relationship between material temperature $T_{s}(t,{\bm{x}})$ and material
energy density $E_{s}(t,{\bm{x}})$ is

$E_{s}=C_{v_{s}}T_{s}>0,\qquad s=e,i,$ (2.3)

where $C_{v_{s}}>0$ represents the heat capacity. $Q_{r},Q_{e}$, and $Q_{i}$
are the source terms for radiation, electrons, and ions, respectively. The
thermal diffusion coefficients $D_{s}({\bm{x}},T_{s})$ in the conduction terms
$\nabla\cdot D_{s}\nabla T_{s}$ (where $s=e,i$) exhibit nonlinear dependence
on temperature as

$D_{s}=K_{s}T_{s}^{5/2},\qquad s=e,i,$ (2.4)

where $K_{s}$ represents the plasma-dependent coefficient [6]. The timescale
of ion-electron coupling is characterized by the collision time coefficient
$\kappa$, which can take various forms, as documented in the literature [3, 8,
26]. In this paper, $\kappa$ is treated as a constant but can vary
significantly over several orders of magnitude. The radiation energy $E_{r}$
is defined as

$E_{r}=\frac{1}{c}\int_{\mathbb{S}^{2}}\psi\mathop{}\\!\mathrm{d}\Omega,$
(2.5)

and the radiation temperature $T_{r}$ follows the relation

$E_{r}=aT_{r}^{4},\qquad T_{r}\geqslant 0.$ (2.6)

In Eq. (2.1), (2.1a) describes the evolution of radiation intensity, while
(2.1b) and (2.1c) govern the evolution of electron and ion temperatures. These
equations couple conduction and energy exchange among radiation, electrons,
and ions. Integrating (2.1a) over $\Omega$ and combining it with (2.1b) and
(2.1c) yields the conservation of total energy, expressed as

$\frac{\partial E_{e}}{\partial t}+\frac{\partial E_{i}}{\partial
t}+\frac{\partial E_{r}}{\partial t}+\nabla\cdot\bm{F}_{r}=\nabla\cdot
D_{e}\nabla T_{e}+\nabla\cdot D_{i}\nabla
T_{i}+Q_{r}/\epsilon^{2}+Q_{e}/\epsilon^{2}+Q_{i}/\epsilon^{2}.$ (2.7)

Here, $\bm{F}_{r}(t,{\bm{x}})$ represents the radiation flux

$\bm{F}_{r}(t,{\bm{x}})=\int_{\mathbb{S}^{2}}\Omega\psi\mathop{}\\!\mathrm{d}\Omega.$
(2.8)

###  2.2 Asymptotic behavior

In the 3TRTM (2.1), two limiting models emerge as the dimensionless parameter
$\epsilon$ approaches zero or when the collision time coefficient $\kappa$
tends to infinity. These limits provide valuable approximations to the
original model and have been thoroughly investigated in [6]. For the sake of
generality, we assume the source term $Q_{s}$ (where $s=r,e,i$) to be zero. In
this section, we will present a detailed discussion of these two limits.

#### Diffusion limit

As $\epsilon$ approaches zero, the radiation intensity $\psi$ becomes
independent of the angular variable $\Omega$ and can be described by the
radiation temperature $T_{r}$. Let $T_{r}^{(0)},T_{e}^{(0)}$, and
$T_{i}^{(0)}$ represent the leading-order terms of $T_{r},T_{e}$, and $T_{i}$,
respectively. This leads to an equilibrium diffusion limit, where the
corresponding local temperature $T^{(0)}$ satisfies a nonlinear diffusion
equation, expressed as follows [9].

$\displaystyle T_{r}^{(0)}=T_{e}^{(0)}=T_{i}^{(0)}=T^{(0)},\qquad
E_{r}(T^{(0)})=a\left(T^{(0)}\right)^{4},$ (2.9a)
$\displaystyle\begin{aligned} &\frac{\partial E_{e}(T^{(0)})}{\partial
t}+\frac{\partial E_{i}(T^{(0)})}{\partial t}+\frac{\partial
E_{r}(T^{(0)})}{\partial t}\\\ &\hskip
56.9055pt+\nabla\cdot\bm{F}(t,{\bm{x}})=\nabla\cdot D_{e}\nabla
T^{(0)}+\nabla\cdot D_{i}\nabla T^{(0)}.\end{aligned}$ (2.9b)

In this approximation, the radiative flux $\bm{F}(t,{\bm{x}})$ is related to
the temperature through Fick’s law of diffusion, expressed as

$\bm{F}(t,{\bm{x}})=\int_{\mathbb{S}^{2}}\Omega\psi\mathop{}\\!\mathrm{d}\Omega=-\frac{ac}{3\sigma}\nabla(T^{(0)})^{4}.$
(2.10)

Eq. (2.9) is referred to as the equilibrium diffusion limit [6]. This limit
provides a simplified description of the behavior of the system when the ratio
of the mean free path to the characteristic length $\epsilon$ becomes very
small.

#### Two temperature limit

When the collision time coefficient $\kappa$ approaches infinity, the system
converges to a two-temperature limit, resulting in a simplified two-
temperature model. In this limit, the electron temperature $T_{e}$ and the ion
temperature $T_{i}$ become equal, denoted as the material temperature $T$, and
the radiation transfer equations (2.1) can be reduced to a 2T (two-
temperature) radiation transfer equation by summing up (2.1b) and (2.1c). The
resulting 2T radiation transfer equation is

$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial\psi}{\partial
t}+\epsilon\Omega\cdot\nabla\psi=\frac{1}{4\pi}\sigma acT^{4}-\sigma\psi,$
(2.11a) $\displaystyle\epsilon^{2}C_{v}\frac{\partial T}{\partial
t}=\epsilon^{2}\nabla\cdot D_{T}\nabla
T+\sigma\int_{\mathbb{S}^{2}}\psi\mathop{}\\!\mathrm{d}\Omega-\sigma acT^{4},$
(2.11b) $\displaystyle C_{v}=C_{v_{e}}+C_{v_{i}},\qquad D_{T}=D_{e}+D_{i},$
(2.11c)

where $C_{v}$ ($D_{T}$) represents the heat capacity (thermal diffusion
coefficient), considering both the heat capacity (diffusion coefficient) of
electrons $C_{v_{e}}$ ($D_{e}$) and that of ions $C_{v_{i}}$ ($D_{i}$).

##  3 Asymptotic preserving time splitting scheme

In many physically meaningful benchmark tests, $\epsilon$ can vary from very
small values to $\mathcal{O}(1)$, and $\kappa$ can span several orders of
magnitude, ranging from $\mathcal{O}(1)$ to very large values. As a result,
the numerical scheme should maintain uniform accuracy for all these extreme
cases. Therefore, it is important to design an efficient scheme that is AP in
both the diffusion limit as $\epsilon$ approaches zero and the 2T limit as
$\kappa$ tends towards infinity. Moreover, the preservation of total energy
(2.7) is of most importance for some benchmark tests, especially when the
regularity of the solution is low [6]. In this section, we propose an AP
conservative time-splitting method.

In this method, the 3TRTM is split into two systems. The first part is close
to the GRTE and can be effectively addressed by extending the existing
efficient AP solvers. The second part is a macroscopic model, where the
radiation density, electron temperature, and ion temperature depend solely on
space and time.

###  3.1 The AP splitting scheme

The splitting scheme is proposed in the framework of micro-macro decomposition
[12, 13]. The radiation intensity is decomposed into a macroscopic density
part that is independent of the angular variable $\Omega$ and the remaining
part:

$\psi(t,{\bm{x}},\Omega)=\rho(t,{\bm{x}})+g(t,{\bm{x}},\Omega),$ (3.1)

with

$\rho(t,{\bm{x}})=\langle\psi(t,{\bm{x}},\Omega)\rangle,\qquad
g(t,{\bm{x}},\Omega)=\psi(t,{\bm{x}},\Omega)-\rho(t,{\bm{x}}),\qquad\langle\cdot\rangle=\frac{1}{|\Omega|}\int_{\mathbb{S}^{2}}\cdot\mathop{}\\!\mathrm{d}\Omega.$
(3.2)

By integrating (2.1a) over the angular space, one can derive the equation for
$\rho$ that depends on $g$. Subtracting the equation for $\rho$ from (2.1a)
leads to the following micro-macro decomposition system:

$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial\rho}{\partial
t}+\epsilon\nabla\cdot\langle\Omega
g\rangle=-\sigma\left(\rho-\frac{acT_{e}^{4}}{|\Omega|}\right),$ (3.3a)
$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial g}{\partial
t}+\epsilon\nabla\cdot[\Omega(\rho+g)-\langle\Omega g\rangle]=-\sigma g,$
(3.3b) $\displaystyle\frac{\epsilon^{2}\partial E_{e}}{\partial
t}=\epsilon^{2}\nabla\cdot D_{e}\nabla
T_{e}+c\kappa(T_{i}-T_{e})+\sigma(|\Omega|\rho-acT_{e}^{4}),$ (3.3c)
$\displaystyle\frac{\epsilon^{2}\partial E_{i}}{\partial
t}=\epsilon^{2}\nabla\cdot D_{i}\nabla T_{i}+c\kappa(T_{e}-T_{i}).$ (3.3d)

The computational time $[0,t_{\rm end}]$ is divided into time intervals
$[t^{n},t^{n+1}]$ for $n=0,1,2,\cdots$. Utilizing this micro-macro
decomposition, the AP time-splitting scheme is formulated as follows.

  1. 1.

With the approximated solution $(\rho^{n},g^{n},E_{e}^{n},E_{i}^{n})$ at time
$t^{n}$, we solve the first system

$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial\rho}{\partial
t}+\epsilon\nabla\cdot\langle\Omega
g\rangle=-\frac{\sigma}{2}\left(\rho-\frac{1}{|\Omega|}acT_{e}^{4}\right),$
(3.4a) $\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial g}{\partial
t}+\epsilon\nabla\cdot[\Omega(\rho+g)-\langle\Omega g\rangle]=-\sigma g,$
(3.4b) $\displaystyle\epsilon^{2}\frac{\partial E_{e}}{\partial
t}=\frac{1}{2}c\kappa(T_{i}-T_{e})+\frac{\sigma}{2}(|\Omega|\rho-
acT_{e}^{4}),$ (3.4c) $\displaystyle\epsilon^{2}\frac{\partial E_{i}}{\partial
t}=\frac{1}{2}c\kappa(T_{e}-T_{i}),$ (3.4d)

for one time step $\Delta t=t^{n+1}-t^{n}$ and get
$(\rho^{n+1,\ast},g^{n+1,\ast},E_{e}^{n+1,\ast},E_{i}^{n+1,\ast})$.

  2. 2.

Using $(\rho^{n+1,\ast},g^{n+1,\ast},E_{e}^{n+1,\ast},E_{i}^{n+1,\ast})$ as
the initial condition, we update the second system

$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial\rho}{\partial
t}=-\frac{\sigma}{2}\left(\rho-\frac{acT_{e}^{4}}{|\Omega|}\right),$ (3.5a)
$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial g}{\partial t}=0,$ (3.5b)
$\displaystyle\epsilon^{2}\frac{\partial E_{e}}{\partial
t}=\epsilon^{2}\nabla\cdot D_{e}\nabla
T_{e}+\frac{1}{2}c\kappa(T_{i}-T_{e})+\frac{\sigma}{2}(|\Omega|\rho-
acT_{e}^{4}),$ (3.5c) $\displaystyle\epsilon^{2}\frac{\partial E_{i}}{\partial
t}=\epsilon^{2}\nabla\cdot D_{i}\nabla T_{i}+\frac{1}{2}c\kappa(T_{e}-T_{i}).$
(3.5d)

for one time step $\Delta t=t^{n+1}-t^{n}$ and get
$(\rho^{n+1},g^{n+1},E_{e}^{n+1},E_{i}^{n+1})$.

In the following, we will demonstrate that this splitting approach effectively
captures both the diffusion limit and the two-temperature limit while
conserving the total energy.

#### Diffusion limit

It is straightforward to demonstrate that as $\epsilon$ tends to zero in
(3.4), the leading-order terms $T_{r}^{(0)},T_{e}^{(0)}$, and $T_{i}^{(0)}$,
satisfy

$T_{e}^{(0)}=T_{i}^{(0)}=T^{(0)},\quad
E_{r}(T_{r})=\frac{1}{c}\int_{\mathbb{S}^{2}}\psi\mathop{}\\!\mathrm{d}\Omega=\frac{|\Omega|}{c}\rho=aT_{e}^{4}+\mathcal{O}(\epsilon)=E_{r}(T^{(0)})+\mathcal{O}(\epsilon),\quad
g=\mathcal{O}(\epsilon).$ (3.6)

Thus, from the relation between $E_{r}$ and $T_{r}$ in (2.6), we have
$T^{(0)}_{r}=T^{(0)}$. By adding up (3.4a), (3.4c), and (3.4d) together, the
source terms on the right-hand sides of these three equations cancel each
other, resulting in the following conservative equation:

$\frac{E_{r}(T^{(0)})}{\partial t}+\frac{\partial E_{e}(T^{(0)})}{\partial
t}+\frac{\partial E_{i}(T^{(0)})}{\partial
t}=-\frac{1}{\epsilon}\nabla\cdot\langle\Omega
g\rangle+\mathcal{O}(\epsilon).$ (3.7)

From (3.4b), one gets

$g=-\frac{\epsilon}{\sigma}\nabla\cdot[\Omega(\rho+g)-\langle\Omega
g\rangle]-\frac{\epsilon^{2}}{c\sigma}\frac{\partial g}{\partial t}.$ (3.8)

Multiplying both sides of (3.8) by $\Omega$ and integrating with respect to
angular variable, one obtains

$-\frac{1}{\epsilon}\nabla\cdot\langle\Omega
g\rangle=\frac{1}{\sigma}\nabla\cdot[\langle\Omega\Omega\rangle\rho]+\mathcal{O}(\epsilon).$
(3.9)

Here, we have used the fact that $g(t,{\bm{x}},\Omega)$ is of
$\mathcal{O}(\epsilon)$. Therefore, (3.7) can be expressed as follows.

$\displaystyle\frac{E_{r}(T^{(0)})}{\partial t}+\frac{\partial
E_{e}(T^{(0)})}{\partial t}+\frac{\partial E_{i}(T^{(0)})}{\partial
t}=\nabla\cdot\left(\frac{ac}{3\sigma}\nabla(T^{(0)})^{4}\right).$ (3.10)

Similarly, as $\epsilon\to 0$, the leading-order terms of (3.5) yield

$T_{r}^{(0)}=T_{e}^{(0)}=T_{i}^{(0)}=T^{(0)},\qquad\frac{|\Omega|}{c}\rho=E_{r}(T^{(0)})+\mathcal{O}(\epsilon).$
(3.11)

By adding up (3.5a), (3.5c), and (3.5d), the interaction terms cancel out, and
the limiting model for (3.5) can be expressed as

$\displaystyle\frac{E_{r}(T^{(0)})}{\partial t}+\frac{\partial
E_{e}(T^{(0)})}{\partial t}+\frac{\partial E_{i}(T^{(0)})}{\partial
t}=\nabla\cdot D_{e}\nabla T^{(0)}+\nabla\cdot D_{i}\nabla T^{(0)}.$ (3.12)

It is evident that the combination of (3.10) and (3.12) constitutes a time-
splitting scheme for the diffusion limit equation (2.9).

#### Two temperature limit

As $\kappa$ tends toward infinity, let $T$ be the material temperature in the
limiting case. and it holds that

$T_{e}=T+\mathcal{O}(\kappa^{-1}),\quad T_{i}=T+\mathcal{O}(\kappa^{-1}).$
(3.13)

We can derive the limiting equations for (3.4) by summing up (3.4c) and
(3.4d). The resulting limiting equations are

$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial\rho}{\partial
t}+\epsilon\nabla\cdot\langle\Omega
g\rangle=-\frac{\sigma}{2}\left(\rho-\frac{acT^{4}}{|\Omega|}\right),$ (3.14)
$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial g}{\partial
t}+\epsilon\nabla\cdot[\Omega(\rho+g)-\langle\Omega g\rangle]=-\sigma g,$
$\displaystyle\epsilon^{2}\frac{\partial E_{T}}{\partial
t}=\frac{\sigma}{2}\big{(}|\Omega|\rho-acT^{4}\big{)},$

where $E_{T}=(C_{v_{e}}+C_{v_{i}})T$. For the second system (3.5), as
$\kappa\to\infty$, one obtains (3.13). By summing up (3.5c) and (3.5d), the
limiting system of (3.5) becomes

$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial\rho}{\partial
t}=-\frac{\sigma}{2}\left(\rho-\frac{acT^{4}}{|\Omega|}\right),$ (3.15)
$\displaystyle\epsilon^{2}\frac{\partial E_{T}}{\partial
t}=\epsilon^{2}\nabla\cdot D_{T}\nabla T+\frac{\sigma}{2}(|\Omega|\rho-
acT^{4}),$

where $E_{T}=\big{(}C_{v_{e}}+C_{v_{i}}\big{)}T$ and
$D_{T}=D_{e}(T)+D_{i}(T)$. It is evident that the combination of equations
(3.14) and (3.15) constitutes a time-splitting method for the two-temperature
limit model (2.11).

###  3.2 Comparison with another splitting method

Compared to the original system, in the first system (3.4), spatial
derivatives only exist in (3.4a) and (3.4b), while the equations for the
electron and ion temperatures reduce into ODEs as in (3.4c) and (3.4d). On the
other hand, the second part (3.5) is a macroscopic system. In some benchmark
problems, the conduction terms $\nabla\cdot D_{s}\nabla T_{s},s=e,i$ are set
to be zero [24] and the second system (3.5) reduces to an ODEs system that can
be easily solved. The specific choices of 1) dividing the term
$\sigma\big{(}\rho-\frac{acT_{e}^{4}}{|\Omega|}\big{)}$ in (3.3a) and (3.3c)
into two $\frac{\sigma}{2}\big{(}\rho-\frac{acT_{e}^{4}}{|\Omega|}\big{)}$,
and placing them in (3.4a), (3.4c), (3.5a), and (3.5c); 2) dividing the term
$c\kappa(T_{i}-T_{e})$ in (3.3c) and (3.3d) into two
$\frac{c\kappa}{2}(T_{i}-T_{e})$, and placing them in (3.4c), (3.4d), (3.5c),
and (3.5d); are crucial for achieving the AP properties.

A splitting scheme for the macroscopic diffusion model is proposed in [12,
13]. Extending this approach, one can design a splitting scheme for the 3T
model by a straightforward extension. We split the 3T model into a classical
two-temperature radiative transfer system and a macroscopic system for the
electron and ion temperatures, as shown below:

$\displaystyle\text{Microscopic part:}\left\\{\begin{aligned}
&\frac{\epsilon^{2}}{c}\frac{\partial\psi}{\partial
t}+\epsilon\Omega\cdot\nabla\psi=-\sigma(\psi-\frac{1}{|\Omega|}acT_{e}^{4}),\\\
&\epsilon^{2}\frac{\partial E_{e}}{\partial t}=\sigma(|\Omega|\rho-
acT_{e}^{4}),\\\ &\epsilon^{2}\frac{\partial E_{i}}{\partial
t}=0.\end{aligned}\right.$ (3.16) $\displaystyle\text{Macroscopic
part:}\left\\{\begin{aligned} &\frac{\partial\psi}{\partial t}=0,\\\
&\epsilon^{2}\frac{\partial E_{e}}{\partial t}=\epsilon^{2}\nabla\cdot
D_{e}\nabla T_{e}+c\kappa(T_{i}-T_{e}),\\\ &\epsilon^{2}\frac{\partial
E_{i}}{\partial t}=\epsilon^{2}\nabla\cdot D_{i}\nabla
T_{i}+c\kappa(T_{e}-T_{i}).\end{aligned}\right.$ (3.17)

Here (3.16) has the same form as the GRTE, while (3.17) is a macroscopic
system that depends solely on $T_{e}$ and $T_{i}$. However, as $\epsilon$
approaches zero, the diffusion limit of the splitting system (3.16) and (3.17)
becomes

$\displaystyle\left\\{\begin{aligned} &\frac{\partial E_{r}(T^{(0)})}{\partial
t}+\frac{\partial E_{e}(T^{(0)})}{\partial
t}=\nabla\cdot\frac{ac}{3\sigma}\nabla\left(T^{(0)}\right)^{4},\\\
&T_{e}=T_{r}=T^{(0)},\quad
E_{r}(T^{(0)})=a\left(T^{(0)}\right)^{4},\end{aligned}\right.$ (3.18)
$\displaystyle\left\\{\begin{aligned} &\frac{\partial E_{e}(T^{(0)})}{\partial
t}+\frac{\partial E_{i}(T^{(0)})}{\partial t}=\nabla\cdot D_{e}\nabla
T_{e}+\nabla\cdot D_{i}\nabla T_{i},\\\
&T_{e}=T_{i}=T^{(0)}.\end{aligned}\right.$ (3.19)

This method is no longer a time-splitting approach for the diffusion limit
equation (2.9).

##  4 Temporal and angular discretization for the subsystems

In this section, we introduce the temporal and angular discretization based on
the splitting method in Sec. 3. We consider the nonlinearity in the
coefficients $\sigma$ and $D_{s}$ and employ a fully implicit solver for both
systems to ensure accuracy through the AP properties discussed in the previous
section. However, the computational cost of nonlinear iterations for a high-
dimensional system is extremely high, and efficient iterative solvers must be
designed. We adopt two alternating iteration methods, so that, under a larger
time step, these iterative methods can still converge quickly.

###  4.1 Splitting of $P_{N}$ system for angular discretization

We discretize the angular variable by $P_{N}$ method, where the radiative
intensity is approximated by a series of basis functions [15, 18, 19]. Without
loss of generality, we consider the following 3TRTM in slab geometry:

$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial\psi}{\partial
t}+\epsilon\mu\frac{\partial\psi}{\partial x}=-\sigma\psi+\frac{1}{2}\sigma
acT_{e}^{4},$ (4.1a) $\displaystyle\frac{\epsilon^{2}\partial E_{e}}{\partial
t}=\epsilon^{2}\frac{\partial}{\partial x}\left(D_{e}\frac{\partial
T_{e}}{\partial
x}\right)+c\kappa(T_{i}-T_{e})+\sigma\int_{-1}^{1}\psi\mathop{}\\!\mathrm{d}\mu-\sigma
acT_{e}^{4},$ (4.1b) $\displaystyle\frac{\epsilon^{2}\partial E_{i}}{\partial
t}=\epsilon^{2}\frac{\partial}{\partial x}\left(D_{i}\frac{\partial
T_{i}}{\partial x}\right)+c\kappa(T_{e}-T_{i}),$ (4.1c)

where $\mu\in[-1,1]$. The details of how to project a spatially 3D RTE with
angular variables defined on the surface of a 3D sphere to slab geometry can
be found in [16], and we omit the details here.

In the framework of the $P_{N}$ method [19], the radiative intensity $\psi$ is
approximated as

$\psi(t,x,\mu)\approx\sum_{l=0}^{M}\frac{2l+1}{2}\psi_{l}(t,x)P_{l}(\mu),$
(4.2)

where $P_{l}(\mu)$ represents the Legendre polynomial of degree $l$ and
$M\in\mathbb{N}$ is the expansion order. Utilizing the following orthogonality
property of the Legendre polynomial:

$\int_{-1}^{1}P_{m}(\mu)P_{n}(\mu)\mathop{}\\!\mathrm{d}\mu=\frac{2}{2n+1}\delta_{m,n},\quad
m,n\in\mathbb{N},$ (4.3)

the expansion coefficients $\psi_{l}(t,x)$ are

$\psi_{l}(t,x)=\int_{-1}^{1}P_{l}(\mu)\psi(t,x,\mu)\mathop{}\\!\mathrm{d}\mu,\qquad
l=0,\cdots,M.$ (4.4)

Therefore, it holds that

$\psi_{0}=2\rho.$ (4.5)

By substituting (4.2) into the governing equations of the 1D model (4.1a) and
utilizing the orthogonality of the basis functions (4.3), as well as the
following recursive property:

$P_{k+1}(\mu)=\frac{2k+1}{k+1}\mu P_{k}(\mu)-\frac{k}{k+1}P_{k-1}(\mu),\quad
k\in\mathbb{N}^{+},$ (4.6)

the $P_{N}$ system can be written as follows.

$\displaystyle\frac{2\epsilon^{2}}{c}\frac{\partial\rho}{\partial
t}+\epsilon\frac{\partial\psi_{1}}{\partial x}+2\sigma\rho=\sigma
acT_{e}^{4},$ (4.7a)
$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial\psi_{1}}{\partial
t}+\epsilon\left(\frac{2}{3}\frac{\partial\rho}{\partial
x}+\frac{2}{3}\frac{\partial\psi_{2}}{\partial x}\right)+\sigma\psi_{1}=0,$
(4.7b) $\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial\psi_{l}}{\partial
t}+\epsilon\left(\frac{l}{2l+1}\frac{\partial\psi_{l-1}}{\partial
x}+\frac{l+1}{2l+1}\frac{\partial\psi_{l+1}}{\partial
x}\right)+\sigma\psi_{l}=0,\quad 2\leqslant l\leqslant M-1,$ (4.7c)
$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial\psi_{M}}{\partial
t}+\frac{\epsilon M}{2M+1}\frac{\partial\psi_{M-1}}{\partial
x}+\sigma\psi_{M}=0.$ (4.7d)

Comparing (4.7) with (3.3a) and (3.3b), the equation for the microscopic part
$g$ becomes a system for
$(\psi_{1},\psi_{2},\cdots,\psi_{M})^{T}=(\psi_{1},\bar{{\bm{u}}}^{T})^{T}$.
Let coefficient matrices be ${\bf
A}=(A_{i,j})\in{\mathbb{R}}^{(M-1)\times(M+1)}$ and ${\bf
B}=(B_{i,j})\in{\mathbb{R}}^{(M-1)\times(M+1)}$, where

$A_{l-1,l}=\frac{l}{2l+1},\quad l=2,\cdots,M,\qquad
B_{l-1,l+2}=\frac{l+1}{2l+1},\quad l=2,\cdots,M-1,$ (4.8)

and all other entries are zero. Eq. (4.7) can then be rewritten as

$\displaystyle\frac{2\epsilon^{2}}{c}\frac{\partial\rho}{\partial
t}+\epsilon\frac{\partial\psi_{1}}{\partial x}=-\sigma(2\rho-acT_{e}^{4}),$
(4.9a) $\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial\psi_{1}}{\partial
t}+\frac{2\epsilon}{3}\frac{\partial\rho}{\partial
x}+\frac{2\epsilon}{3}\frac{\partial\psi_{2}}{\partial x}=-\sigma\psi_{1},$
(4.9b)
$\displaystyle\frac{\epsilon^{2}}{c}\frac{\partial\bar{\bm{u}}}{\partial
t}+\epsilon{\bf A}\frac{\partial{\bm{u}}}{\partial x}+\epsilon{\bf
B}\frac{\partial{\bm{u}}}{\partial x}=-\sigma\bar{\bm{u}},$ (4.9c)

where ${\bm{u}}=(2\rho,\psi_{1},\bar{\bm{u}}^{T})^{T}$. The governing equation
for the electron becomes

$\epsilon^{2}\frac{\partial E_{e}}{\partial
t}=\epsilon^{2}\frac{\partial}{\partial x}\left(D_{e}\frac{\partial
T_{e}}{\partial x}\right)+c\kappa(T_{i}-T_{e})+\sigma(2\rho-acT_{e}^{4}),$
(4.10)

and the equation for the conservation of the total energy (2.7) is written as

$\frac{\partial E_{e}}{\partial t}+\frac{\partial E_{i}}{\partial
t}+\frac{2}{c}\frac{\partial\rho}{\partial
t}=-\frac{1}{\epsilon}\frac{\partial\psi_{1}}{\partial
x}+\frac{\partial}{\partial x}\left(D_{e}\frac{\partial T_{e}}{\partial
x}\right)+\frac{\partial}{\partial x}\left(D_{i}\frac{\partial T_{i}}{\partial
x}\right).$ (4.11)

Eqs. (4.9), (4.10) and (4.1c) together form the $P_{N}$ system for the 1D
3TRTM.

Utilizing the splitting method introduced in the previous section, the two
parts of the splitting approach for this $P_{N}$ system are as follows.

First part (microscopic part): $\displaystyle\left\\{\begin{aligned}
&\frac{2\epsilon^{2}}{c}\frac{\partial\rho}{\partial
t}+\epsilon\frac{\partial\psi_{1}}{\partial x}=-\frac{\sigma}{2}(2\rho-
acT_{e}^{4}),\\\ &\frac{\epsilon^{2}}{c}\frac{\partial\psi_{1}}{\partial
t}+\frac{2\epsilon}{3}\frac{\partial\rho}{\partial
x}+\frac{2\epsilon}{3}\frac{\partial\psi_{2}}{\partial x}=-\sigma\psi_{1},\\\
&\frac{\epsilon^{2}}{c}\frac{\partial\bar{\bm{u}}}{\partial t}+\epsilon{\bf
A}\frac{\partial{\bm{u}}}{\partial x}+\epsilon{\bf
B}\frac{\partial{\bm{u}}}{\partial x}=-\sigma\bar{\bm{u}},\\\
&\epsilon^{2}C_{v_{e}}\frac{\partial T_{e}}{\partial
t}=\frac{1}{2}c\kappa(T_{i}-T_{e})+\frac{\sigma}{2}(2\rho-acT_{e}^{4}),\\\
&\epsilon^{2}C_{v_{i}}\frac{\partial T_{i}}{\partial
t}=\frac{1}{2}c\kappa(T_{e}-T_{i}),\end{aligned}\right.$ (4.12a) Second part
(macroscopic part): $\displaystyle\left\\{\begin{aligned}
&\frac{2\epsilon^{2}}{c}\frac{\partial\rho}{\partial
t}=-\frac{\sigma}{2}(2\rho-acT_{e}^{4}),\\\
&\epsilon^{2}C_{v_{e}}\frac{\partial T_{e}}{\partial
t}=\epsilon^{2}\frac{\partial}{\partial x}\left(D_{e}\frac{\partial
T_{e}}{\partial
x}\right)+\frac{1}{2}c\kappa(T_{i}-T_{e})+\frac{\sigma}{2}(2\rho-
acT_{e}^{4}),\\\ &\epsilon^{2}C_{v_{i}}\frac{\partial T_{i}}{\partial
t}=\epsilon^{2}\frac{\partial}{\partial x}\left(D_{i}\frac{\partial
T_{i}}{\partial x}\right)+\frac{1}{2}c\kappa(T_{e}-T_{i}),\\\
&\frac{\partial\psi_{1}}{\partial
t}=0,\quad\frac{\partial\bar{\bm{u}}}{\partial
t}=\bm{0},\quad\bm{0}=(0,0,\dots,0)^{T}\in\mathbb{R}^{M-1}.\end{aligned}\right.$
(4.12b)

###  4.2 Fully temporal discretization and the iterative solver

In the split $P_{N}$ system, the nonlinear terms include the fourth power of
temperature and the nonlinear conduction term. Handling systems where both the
conduction term and the fourth power term are implicit presents challenges.
Explicitly treating the conduction term during the iteration process would
impose a time step requirement of $\Delta t=\mathcal{O}(\Delta x^{2})$. On the
other hand, explicitly treating the fourth power term during the iteration
process would result in a time step requirement dependent on $\epsilon$. While
the fully implicit scheme is utilized for both systems, designing an efficient
iterative solver is crucial to avoid solving a large nonlinear system,
especially for the first system (4.12a).

To iteratively solve these systems, we design an alternating iterative method.
It consists of two parts: one implicitly solves a linearized system, and the
nonlinear part is considered explicitly. The other part implicitly solves the
quartic equation at each space point, while the linearized system part is
considered explicitly. The iterations of these two parts are performed
alternately until certain convergence criteria are satisfied, and the
numerical solution of the implicit system can be obtained. An AP solver for
GRTE is developed in [10, 11], whose computational cost is comparable to an
explicit scheme, while its accuracy and stability are independent of the
multiscale parameters $\epsilon$. In the subsequent part, we propose an
efficient iteration method for (4.12a) based on the AP solver for GRTE in
[10].

According to the order analysis in [10], the orders of the expansion
coefficients are given by

$\psi_{l}=\mathcal{O}{(\epsilon^{l})},\qquad\frac{\partial\psi_{l}}{\partial
x}=\mathcal{O}{(\epsilon^{l})}.$ (4.13)

For the convenience of readers, we provide the details of order analysis in A.
In Eq. (4.12a), ${\bf A}\frac{\partial{\bm{u}}}{\partial x}$ has a low order
of $\epsilon$ and is treated implicitly, while the terms ${\bf
B}\frac{\partial{\bm{u}}}{\partial x}$ are treated explicitly since their
orders are higher. The terms on the right-hand side are treated implicitly,
while the coefficient $\sigma$, which may nonlinearly depend on $T_{e}$, is
treated explicitly to simplify the iteration process.

We now present the details of the iterative method for the microscopic part,
as displayed in Eqs. (4.14) and (4.15). Eq. (4.14) is a direct generalization
of the IMEX method presented in [10]. The approach outlined in [10] requires
solving a fourth-order polynomial equation for each spatial grid, offering
computational efficiency akin to an explicit scheme while benefiting from a
time step at the same order of the mesh size. The advection terms of Eq.
(4.15) are updated by the results obtained in (4.14), while only a local
nonlinear system has to be solved at each grid point.

$\displaystyle\frac{2\epsilon^{2}}{c}\frac{\rho^{n+1,\ast,2k+1}-\rho^{n}}{\Delta
t}+\epsilon\left(\frac{\partial\psi_{1}}{\partial
x}\right)^{n+1,\ast,2k}=-\frac{\sigma^{n+1,\ast,2k}}{2}(2\rho-
acT_{e}^{4})^{n+1,\ast,2k+1},$ (4.14a)
$\displaystyle\frac{\epsilon^{2}}{c}\frac{\psi_{1}^{n+1,\ast,2k+1}-\psi_{1}^{n}}{\Delta
t}+\frac{2\epsilon}{3}\left(\frac{\partial\rho}{\partial
x}\right)^{n+1,\ast,2k+1}+\frac{2\epsilon}{3}\left(\frac{\partial\psi_{2}}{\partial
x}\right)^{n+1,\ast,2k}=-\sigma^{n+1,\ast,2k+1}\psi_{1}^{n+1,\ast,2k+1},$
(4.14b)
$\displaystyle\frac{\epsilon^{2}}{c}\frac{\bar{\bm{u}}^{n+1,\ast,2k+1}-\bar{\bm{u}}^{n}}{\Delta
t}+\epsilon{\bf A}\left(\frac{\partial{\bm{u}}}{\partial
x}\right)^{n+1,\ast,2k+1}+\epsilon{\bf
B}\left(\frac{\partial{\bm{u}}}{\partial
x}\right)^{n+1,\ast,2k}=-\sigma^{n+1,\ast,2k+1}\bar{\bm{u}}^{n+1,\ast,2k+1}$
(4.14c)
$\displaystyle\epsilon^{2}C_{v_{e}}\frac{T_{e}^{n+1,\ast,2k+1}-T_{e}^{n}}{\Delta
t}=\frac{1}{2}c\kappa(T_{i}-T_{e})^{n+1,\ast,2k+1}+\frac{\sigma^{n+1,\ast,2k}}{2}(2\rho-
acT_{e}^{4})^{n+1,\ast,2k+1},$ (4.14d)
$\displaystyle\epsilon^{2}C_{v_{i}}\frac{T_{i}^{n+1,\ast,2k+1}-T_{i}^{n}}{\Delta
t}=\frac{1}{2}c\kappa(T_{e}-T_{i})^{n+1,\ast,2k+1},$ (4.14e)

$\displaystyle\frac{2\epsilon^{2}}{c}\frac{\rho^{n+1,\ast,2k+2}-\rho^{n}}{\Delta
t}+\epsilon\left(\frac{\partial\psi_{1}}{\partial
x}\right)^{n+1,\ast,2k+2}=-\frac{\sigma^{n+1,\ast,2k+1}}{2}\bigg{(}2\rho^{n+1,\ast,2k+2}-ac(T_{e}^{n+1,\ast,2k+1})^{3}T_{e}^{n+1,\ast,2k+2}\bigg{)},$
(4.15a)
$\displaystyle\frac{\epsilon^{2}}{c}\frac{\psi_{1}^{n+1,\ast,2k+2}-\psi_{1}^{n}}{\Delta
t}+\frac{2\epsilon}{3}\left(\frac{\partial\rho}{\partial
x}\right)^{n+1,\ast,2k+2}+\frac{2\epsilon}{3}\left(\frac{\partial\psi_{2}}{\partial
x}\right)^{n+1,\ast,2k+1}=-\sigma^{n+1,\ast,2k+1}\psi_{1}^{n+1,\ast,2k+2}$
(4.15b)
$\displaystyle\frac{\epsilon^{2}}{c}\frac{\bar{\bm{u}}^{n+1,\ast,2k+2}-\bar{\bm{u}}^{n}}{\Delta
t}+\epsilon{\bf A}\left(\frac{\partial{\bm{u}}}{\partial
x}\right)^{n+1,\ast,2k+2}+\epsilon{\bf
B}\left(\frac{\partial{\bm{u}}}{\partial
x}\right)^{n+1,\ast,2k+1}=-\sigma^{n+1,\ast,2k+2}\bar{\bm{u}}^{n+1,\ast,2k+2}$
(4.15c)
$\displaystyle\epsilon^{2}C_{v_{e}}\frac{T_{e}^{n+1,\ast,2k+2}-T_{e}^{n}}{\Delta
t}=\frac{1}{2}c\kappa(T_{i}-T_{e})^{n+1,\ast,2k+2}+\frac{\sigma^{n+1,\ast,2k+1}}{2}\bigg{(}2\rho^{n+1,\ast,2k+2}-ac(T_{e}^{n+1,\ast,2k+1})^{3}T_{e}^{n+1,\ast,2k+2}\bigg{)},$
(4.15d)
$\displaystyle\epsilon^{2}C_{v_{i}}\frac{T_{i}^{n+1,\ast,2k+2}-T_{i}^{n}}{\Delta
t}=\frac{1}{2}c\kappa(T_{e}-T_{i})^{n+1,\ast,2k+2},$ (4.15e)

where $()^{n+1,\ast,2k+1}$ is the numerical solution after the step (4.14) in
one iteration. When solving (4.14), the governing equation of $T_{e}$ is
obtained first. Specifically, by substituting the governing equation of $\rho$
(4.14a) and $T_{i}$ (4.14e) into that of $T_{e}$ (4.14d) in (4.14), a fourth-
order polynomial equation of $T_{e}^{n+1,\ast,2k+1}$ is derived as

$\mathscr{C}_{4}(T_{e}^{n+1,\ast,2k+1})^{4}+\mathscr{C}_{1}(T_{e}^{n+1,\ast,2k+1})+\mathscr{C}_{0}=0.$
(4.16)

Here, the coefficients are given by

$\begin{gathered}\mathscr{C}_{4}=\frac{\sigma^{n+1,\ast,2k}ac\Delta
t}{2\epsilon^{2}+\sigma^{n+1,\ast,2k}c\Delta
t},\qquad\mathscr{C}_{1}=C_{v_{e}}+\frac{C_{v_{i}}c\kappa\Delta
t}{2\epsilon^{2}C_{v_{i}}+c\kappa\Delta t},\\\
\mathscr{C}_{0}=-C_{v_{e}}T_{e}^{n}-\frac{C_{v_{i}}c\kappa\Delta
t}{2\epsilon^{2}C_{v_{i}}+c\kappa\Delta
t}T_{i}^{n}-\frac{2\sigma^{n+1,\ast,2k}\Delta
t}{2\epsilon^{2}+\sigma^{n+1,\ast,2k}\Delta
tc}\rho^{n}+\frac{\sigma^{n+1,\ast,2k}(\Delta
t)^{2}c}{\epsilon(2\epsilon^{2}+\sigma^{n+1,\ast,2k}\Delta
tc)}\left(\frac{\partial\psi_{1}}{\partial
x}\right)^{n+1,\ast,2k}.\end{gathered}$ (4.17)

This equation is a key component of the iterative method and allows for the
efficient update of the electron temperature $T_{e}$. The existence and
uniqueness of the solution to the fourth-order polynomial equation (4.16),
with the coefficients provided in (4.17), are ensured as long as
$\rho^{n+1,\ast,2k+1}$, $T_{e}^{n+1,\ast,2k+1}$, and $T_{i}^{n+1,\ast,2k+1}$
exist. In fact, let

$f(T)=\mathscr{C}_{4}T^{4}+\mathscr{C}_{1}T+\mathscr{C}_{0}.$ (4.18)

Then, it is easy to verify that $f^{\prime}(T)>0$ for all $T>0$ and $f(0)<0$.
Thus, the positive solution of (4.16) is unique.

Once $T_{e}^{n+1,\ast,2k+1}$ is obtained,
$\rho^{n+1,\ast,2k+1},\sigma^{n+1,\ast,2k+1}$ and $T_{i}^{n+1,\ast,2k+1}$ can
be directly calculated using (4.14a), (4.14d) and (4.14e). The high-order
terms, $\psi_{l}^{n+1,\ast,2k+1},l=1,\cdots,M$, can be solved successively.
Notably, this numerical scheme can be solved with the same computational cost
as an explicit scheme, even though the dominant terms are treated implicitly.
For example, when updating $\psi_{1}^{n+1,\ast,2k+1}$, the implicit term is
$\rho^{n+1,\ast,2k+1}$, which has already obtained by (4.14a). When solving
the system (4.15), a linear system of (4.15a), (4.15b), (4.15d) and (4.15e) is
first solved to obtain
$\rho^{n+1,\ast,2k+2},\psi_{1}^{n+1,\ast,2k+2},T_{e}^{n+1,\ast,2k+2}$ and
$T_{i}^{n+1,\ast,2k+2}$. Then, the high-order terms,
$\psi_{l}^{n+1,\ast,2k+2},l=2,\cdots,M$, are solved successively as in
(4.14c).

The alternative iteration method is also adopted for the temporal
discretization of the macroscopic part (4.12b). Assuming the numerical
solution after the iteration of the first part is $()^{n+1,\ast}$, the
iterative method for the macroscopic part is displayed in Eqs. (4.19) and
(4.20). In Eq. (4.19), the interaction terms on the right side are treated
implicitly to allow for larger time step length, while the conduction terms
including the coefficients such as the diffusion coefficients $D_{s},s=e,i$
and the opacity $\sigma$ are treated explicitly to simplify the iteration
process. Eq. (4.19) solves a local nonlinear system at each spatial grid where
the governing equation of $T_{e}$, which is a fourth-order polynomial
equation, is initially solved. Subsequently, in Eq. (4.20), the conduction
coefficients are updated by the results obtained in Eq. (4.19), and two linear
parabolic equations are solved.

$\displaystyle\frac{2\epsilon^{2}}{c}\frac{\rho^{n+1,2k+1}-\rho^{n+1,\ast}}{\Delta
t}=-\frac{\sigma^{n+1,2k}}{2}\left(2\rho-acT_{e}^{4}\right)^{n+1,2k+1},$
(4.19a) $\displaystyle\begin{aligned}
&\epsilon^{2}C_{v_{e}}\frac{T_{e}^{n+1,2k+1}-T_{e}^{n+1,\ast}}{\Delta
t}=\epsilon^{2}\frac{\partial}{\partial
x}\left(D_{e}^{n+1,2k}\left(\frac{\partial T_{e}}{\partial
x}\right)^{n+1,2k}\right)+\frac{1}{2}c\kappa(T_{i}-T_{e})^{n+1,2k+1}\\\
&\qquad\qquad\qquad\qquad+\frac{\sigma^{n+1,2k}}{2}(2\rho-
acT_{e}^{4})^{n+1,2k+1},\end{aligned}$ (4.19b)
$\displaystyle\epsilon^{2}C_{v_{i}}\frac{T_{i}^{n+1,2k+1}-T_{i}^{n+1,\ast}}{\Delta
t}=\epsilon^{2}\frac{\partial}{\partial
x}\left(D_{i}^{n+1,2k}\left(\frac{\partial T_{i}}{\partial
x}\right)^{n+1,2k}\right)+\frac{1}{2}c\kappa(T_{e}-T_{i})^{n+1,2k+1},$ (4.19c)

$\displaystyle\frac{2\epsilon^{2}}{c}\frac{\rho^{n+1,2k+2}-\rho^{n+1,\ast}}{\Delta
t}=-\frac{\sigma^{n+1,2k+1}}{2}\left(2\rho^{n+1,2k+2}-ac(T_{e}^{n+1,2k+1})^{3}T_{e}^{n+1,2k+2}\right),$
(4.20a) $\displaystyle\begin{aligned}
&\epsilon^{2}C_{v_{e}}\frac{T_{e}^{n+1,2k+2}-T_{e}^{n+1,\ast}}{\Delta
t}=\epsilon^{2}\frac{\partial}{\partial
x}\left(D_{e}^{n+1,2k+1}\left(\frac{\partial T_{e}}{\partial
x}\right)^{n+1,2k+2}\right)+\frac{1}{2}c\kappa(T_{i}-T_{e})^{n+1,2k+2}\\\
&\qquad\qquad\qquad\qquad+\frac{\sigma^{n+1,2k+1}}{2}(2\rho^{n+1,2k+2}-ac(T_{e}^{n+1,2k+1})^{3}T_{e}^{n+1,2k+2}),\end{aligned}$
(4.20b)
$\displaystyle\epsilon^{2}C_{v_{i}}\frac{T_{i}^{n+1,2k+2}-T_{i}^{n+1,\ast}}{\Delta
t}=\epsilon^{2}\frac{\partial}{\partial
x}\left(D_{i}^{n+1,2k+1}\left(\frac{\partial T_{i}}{\partial
x}\right)^{n+1,2k+2}\right)+\frac{1}{2}c\kappa(T_{e}-T_{i})^{n+1,2k+2}.$
(4.20c)

However, using solely (4.20) is not sufficient for the iteration process
involving the macroscopic part (4.12b). This limitation arises because the
linearization of the term $T^{4}$ in Eq. (4.20) will no longer match the
density $\rho$ as $\epsilon$ diminishes. Consequently, the alternating
iteration method is adopted, wherein the numerical solution $()^{n+1,2k+1}$
serves as a precondition. Additionally, numerical observations reveal that in
this alternating iteration method, the time step length remains independent of
$\epsilon$ and maintains the same order as the spatial step. Besides, when
$D_{e}=D_{i}=0$, (4.12b) reduces to an ODEs system, greatly simplifying the
numerical solution of the system.

##  5 Spatial discretization and the AP property

In this section, the spatial discretization of the system (4.12a) and (4.12b)
is introduced. The finite volume method is applied to (4.12a), while the
finite difference method is used for (4.12b). The AP property and the energy
conservation properties of the fully discretized scheme are discussed.

###  5.1 Spatial discretization

The spatial discretization for the model (4.14), (4.15), (4.19) and (4.20) in
the slab geometry is introduced, and the two-dimensional extension is
straightforward. Let the computational domain be $[0,1]$ and $\Delta
x=\frac{1}{N}$. $x_{j-1/2}=(j-1/2)\Delta x,j=1,\cdots,N$ represent the uniform
mesh in the Cartesian coordinates. Let $\Gamma_{j}$ denote the cell
$[x_{j-1/2},x_{j+1/2}]$. $T_{e,j}^{n}$, $T_{i,j}^{n}$, and ${\bm{u}}_{j}^{n}$
represent the average of $T_{e},T_{i}$ and ${\bm{u}}$ at time $t^{n}$ inside
cell $\Gamma_{j}$, respectively. The finite volume method is utilized to
approximate the convection term in (4.14) and (4.15). The revised Lax-
Friedrichs scheme described in [10] is adopted to determine the numerical
flux.

$\displaystyle\frac{l+1}{2l+1}\left(\frac{\partial\psi_{l+1}}{\partial
x}\right)_{j}^{\rm old}\approx\frac{1}{\Delta x}\left(\mathcal{F}_{l+1}^{\rm
old}\left({\bm{u}}_{j+1/2}^{L,\rm old},{\bm{u}}_{j+1/2}^{R,\rm
old}\right)-\mathcal{F}_{l+1}^{\rm old}\left({\bm{u}}_{j-1/2}^{L,\rm
old},{\bm{u}}_{j-1/2}^{R,\rm old}\right)\right),\ l=0,1,\cdots,M-1,$ (5.1a)
$\displaystyle\frac{l}{2l+1}\left(\frac{\partial\psi_{l-1}}{\partial
x}\right)_{j}^{\rm new}\approx\frac{1}{\Delta x}\left(\mathcal{F}_{l-1}^{\rm
new}\left({\bm{u}}_{j+1/2}^{L,\rm new},{\bm{u}}_{j+1/2}^{R,\rm
new}\right)-\mathcal{F}_{l-1}^{\rm new}\left({\bm{u}}_{j-1/2}^{L,\rm
new},{\bm{u}}_{j-1/2}^{R,\rm new}\right)\right),\quad l=1,2,\cdots,M.$ (5.1b)

where $()^{\rm new}$ is $()^{n+1,\ast,2k+1}$ and $()^{\rm old}$ is
$()^{n+1,\ast,2k}$ for the numerical flux in (4.14), while $()^{\rm new}$ is
$()^{n+1,\ast,2k+2}$ and $()^{\rm old}$ is $()^{n+1,\ast,2k+1}$ for the
numerical flux in (4.15). The exact forms of $\mathcal{F}_{l}^{\rm old}$ and
$\mathcal{F}_{l}^{\rm new}$ are given by

$\displaystyle\mathcal{F}_{l+1}^{\rm
old}({\bm{u}}^{L}_{j+\frac{1}{2}},{\bm{u}}^{R}_{j+\frac{1}{2}})=\frac{l+1}{2l+1}\frac{1}{2}(\psi_{l+1,j+\frac{1}{2}}^{L}+\psi_{l+1,j+\frac{1}{2}}^{R})-\frac{\alpha_{j+\frac{1}{2}}}{2}(\psi_{l,j+\frac{1}{2}}^{R}-\psi_{l,j+\frac{1}{2}}^{L}),\quad
l=0,1,\dots,M-1,$ (5.2a) $\displaystyle\mathcal{F}_{l-1}^{\rm
new}({\bm{u}}^{L}_{j+\frac{1}{2}},{\bm{u}}^{R}_{j+\frac{1}{2}})=\frac{l}{2l+1}\frac{1}{2}(\psi_{l-1,j+\frac{1}{2}}^{L}+\psi_{l-1,j+\frac{1}{2}}^{R}),\quad
l=1,2,\dots,M-1,$ (5.2b) $\displaystyle\begin{aligned} &\mathcal{F}_{M-1}^{\rm
new}({\bm{u}}^{L}_{j+\frac{1}{2}},{\bm{u}}^{R}_{j+\frac{1}{2}})=\frac{M}{2M+1}\frac{1}{2}(\psi_{M-1,j+\frac{1}{2}}^{L}+\psi_{M-1,j+\frac{1}{2}}^{R})-\frac{\alpha_{j+\frac{1}{2}}}{2}(\psi_{M,j+\frac{1}{2}}^{R}-\psi_{M,j+\frac{1}{2}}^{L}),\end{aligned}$
(5.2c)

with

$\alpha_{j+\frac{1}{2}}=\alpha(\sigma,\epsilon)=\frac{c}{\epsilon^{2}}\exp\left(-\frac{\sigma_{j}+\sigma_{j+1}}{2\epsilon^{2}}\right),$
(5.3)

and ${\bm{u}}^{L}$ and ${\bm{u}}^{R}$ are given by reconstruction. The
simplest way is

${\bm{u}}_{j+\frac{1}{2}}^{L}={\bm{u}}_{j},\quad{\bm{u}}_{j+\frac{1}{2}}^{R}={\bm{u}}_{j+1}.$
(5.4)

To improve the resolution of the solution, linear reconstruction or WENO
reconstruction can be used.

######  Remark 1.

The different forms of flux are due to the non-consistent governing equation
of $\psi_{M}$. The choice of $\alpha_{j}$ is necessary to ensuring the AP
property of the numerical scheme and we refer to [10] for more details.

For the diffusion terms in (4.19) and (4.20), the finite difference method is
employed and is expressed as

$\displaystyle\frac{\partial}{\partial x}\left(D_{s}^{p}\left(\frac{\partial
T_{s}}{\partial x}\right)_{j}^{q}\right)\approx\frac{1}{\Delta
x}\left(\frac{D_{s,j+1}^{p}+D_{s,j}^{p}}{2}\frac{T_{s,j+1}^{q}-T_{s,j}^{q}}{\Delta
x}-\frac{D_{s,j}^{p}+D_{s,j-1}^{p}}{2}\frac{T_{s,j}^{q}-T_{s,j-1}^{q}}{\Delta
x}\right),\quad s=e,i,$ (5.5)

where $p,q$ is the label for the time step.

###  5.2 AP property

To demonstrate the AP property, it is necessary to show that for a fixed mesh
size and time step, the numerical scheme should automatically converge to a
scheme for the diffusion limit (2.9) as the parameter $\epsilon$ approaches
zero. Simultaneously, it should also converge to a scheme for the two-
temperature limit (2.11) as the parameter $\kappa$ tends to infinity. In the
subsequent part, we illustrate the AP property of the numerical scheme
presented in Sec. 4.2 and Sec. 5.1.

#### AP property to the diffusion limit

As the parameter $\epsilon$ approaches zero, the numerical scheme described in
Sec. 5.1 converges to a consistent scheme for the nonlinear diffusion equation
(2.9).

The fully discretized scheme for the $\rho,T_{e}$ and $T_{i}$ are

$\displaystyle\frac{2\epsilon^{2}}{c}\frac{\rho^{n+1,\ast}_{j}-\rho^{n}_{j}}{\Delta
t}+\frac{\mathcal{F}_{1,j+1/2}^{n+1,\ast}-\mathcal{F}_{1,j-1/2}^{n+1,\ast}}{\Delta
x}=-\frac{\sigma^{n+1,\ast}_{j}}{2}(2\rho-acT_{e}^{4})_{j}^{n+1,\ast},$ (5.6)
$\displaystyle\epsilon^{2}C_{v_{e}}\frac{T_{e,j}^{n+1,\ast}-T_{e,j}^{n}}{\Delta
t}=\frac{1}{2}c\kappa(T_{i}-T_{e})^{n+1,\ast}_{j}+{\frac{\sigma_{j}^{n+1,\ast}}{2}}(2\rho-
acT_{e}^{4})^{n+1,\ast}_{j},$
$\displaystyle\epsilon^{2}C_{v_{i}}\frac{T_{i,j}^{n+1,\ast}-T_{i,j}^{n}}{\Delta
t}=\frac{1}{2}c\kappa(T_{e}-T_{i})^{n+1,\ast}_{j},$

where

$\mathcal{F}_{1,j+1/2}^{n+1,\ast}=\frac{\epsilon}{2}(\psi_{1,j}^{n+1,\ast}+\psi_{1,j+1}^{n+1,\ast})-\alpha_{j+1/2}^{n+1,\ast}\epsilon(\rho_{j+1}^{n+1,\ast}-\rho_{j}^{n+1,\ast}).$
(5.7)

Summing up the equations in (5.6), we have

$\displaystyle\frac{2\epsilon^{2}}{c}\frac{\rho_{j}^{n+1,\ast}-\rho_{j}^{n}}{\Delta
t}+\epsilon^{2}C_{v_{e}}\frac{T_{e,j}^{n+1,\ast}-T_{e,j}^{n}}{\Delta
t}+\epsilon^{2}C_{v_{i}}\frac{T_{i,j}^{n+1,\ast}-T_{i,j}^{n}}{\Delta
t}=-\frac{\mathcal{F}_{1,j+1/2}^{n+1,\ast}-\mathcal{F}_{1,j-1/2}^{n+1,\ast}}{\Delta
x}.$ (5.8)

The equation of $\psi_{1}$ in (4.12a) is

$\frac{\epsilon^{2}}{c}\frac{\psi_{1,j}^{n+1,\ast}-\psi_{1,j}^{n}}{\Delta
t}+\frac{\mathcal{F}_{0,j+1/2}^{n+1,\ast}-\mathcal{F}_{0,j-1/2}^{n+1,\ast}}{\Delta
x}+\frac{\mathcal{F}_{2,j+1/2}^{n+1,\ast}-\mathcal{F}_{2,j-1/2}^{n+1,\ast}}{\Delta
x}=-\sigma_{j}^{n+1,\ast}\psi_{1,j}^{n+1,\ast},$ (5.9)

where

$\mathcal{F}_{0,j+1/2}^{n+1,\ast}=\frac{2}{3}\frac{\epsilon}{2}(\rho_{j}^{n+1,\ast}+\rho_{j+1}^{n+1,\ast}),\quad\mathcal{F}_{2,j+1/2}^{n+1,\ast}=\frac{2}{3}\frac{\epsilon}{2}(\psi_{2,j}^{n+1,\ast}+\psi_{2,j+1}^{n+1,\ast})-\alpha_{j+1/2}^{n+1,\ast}\frac{\epsilon}{2}(\psi_{1,j+1}^{n+1,\ast}-\psi_{1,j}^{n+1,\ast}).$
(5.10)

In the limit as $\epsilon$ approaches zero, order analysis (as detailed in A)
yields

$\psi_{1,j}^{n+1,\ast}=\mathcal{O}(\epsilon),\qquad\psi_{2,j}^{n+1,\ast}=\mathcal{O}(\epsilon^{2}).$
(5.11)

With this, we can rewrite the source term in (5.9) as

$-\sigma^{n+1,\ast}_{j}\psi_{1,j}^{n+1,\ast}=\frac{2\epsilon}{3}\frac{\rho_{j+1}^{n+1,\ast}-\rho_{j-1}^{n+1,\ast}}{2\Delta
x}+\mathcal{O}(\epsilon^{2}).$ (5.12)

Substituting (5.12) into (5.7), one obtains

$\mathcal{F}^{n+1,\ast}_{1,j+1/2}=\frac{\epsilon^{2}}{2}\left(\frac{2}{3\sigma_{j}^{n+1,\ast}}\frac{\rho_{j+1}^{n+1,\ast}-\rho_{j-1}^{n+1,\ast}}{2\Delta
x}+\frac{2}{3\sigma_{j+1}^{n+1,\ast}}\frac{\rho_{j+2}^{n+1,\ast}-\rho_{j}^{n+1,\ast}}{2\Delta
x}\right)+\mathcal{O}(\epsilon^{3}).$ (5.13)

With equation (5.13), equation (5.8) is reduced to

$\displaystyle\frac{2\epsilon^{2}}{c}\frac{\rho_{j}^{n+1,\ast}-\rho_{j}^{n}}{\Delta
t}+\epsilon^{2}C_{v_{e}}\frac{T_{e,j}^{n+1,\ast}-T_{e,j}^{n}}{\Delta
t}+\epsilon^{2}C_{v_{i}}\frac{T_{i,j}^{n+1,\ast}-T_{i,j}^{n}}{\Delta t}$
(5.14) $\displaystyle=\frac{\epsilon^{2}}{2\Delta
x}\left(\frac{2}{3\sigma_{j+1}^{n+1,\ast}}\frac{\rho_{j+2}^{n+1,\ast}-\rho_{j}^{n+1,\ast}}{2\Delta
x}-\frac{2}{3\sigma_{j-1}^{n+1,\ast}}\frac{\rho_{j}^{n+1,\ast}-\rho_{j-2}^{n+1,\ast}}{2\Delta
x}\right)+\mathcal{O}(\epsilon^{3}).$

As $\epsilon$ tends to zero, the dominance of the sources in (5.6) leads to

$T_{e,j}^{n+1,\ast}=(T^{(0)})_{j}^{n+1,\ast}+\mathcal{O}(\epsilon),\qquad
T_{i,j}^{n+1,\ast}=(T^{(0)})_{j}^{n+1,\ast}+\mathcal{O}(\epsilon),\qquad
2\rho_{j}^{n+1,\ast}=ac((T^{(0)})^{4})_{j}^{n+1,\ast}+\mathcal{O}(\epsilon).$
(5.15)

By substituting (5.15) into (5.14), it holds that, as $\epsilon$ approaches
zero,

$\displaystyle\frac{a((T^{(0)})^{4})_{j}^{n+1,\ast}-a((T^{(0)})^{4})_{j}^{n}}{\Delta
t}+\frac{C_{v_{e}}(T^{(0)}_{j})^{n+1,\ast}-C_{v_{e}}(T^{(0)}_{j})^{n}}{\Delta
t}+\frac{C_{v_{i}}(T^{(0)}_{j})^{n+1,\ast}-C_{v_{i}}(T^{(0)}_{j})^{n}}{\Delta
t}$ (5.16) $\displaystyle=\frac{1}{2\Delta
x}\left(\frac{ac}{3\sigma((T^{(0)})_{j+1}^{n+1,\ast})}\frac{((T^{(0)})^{4})_{j+2}^{n+1,\ast}-((T^{(0)})^{4})_{j}^{n+1,\ast}}{2\Delta
x}-\frac{ac}{3\sigma((T^{(0)})_{j-1}^{n+1,\ast})}\frac{((T^{(0)})^{4})_{j}^{n+1,\ast}-((T^{(0)})^{4})_{j-2}^{n+1,\ast}}{2\Delta
x}\right),$

which is a numerical scheme for the nonlinear diffusion equation (2.9). We
have completed the proof of the first part.

In the second part, as $\epsilon$ tends to zero in (4.12b), it holds that

$T_{e,j}^{n+1}=(T^{(0)})_{j}^{n+1}+\mathcal{O}(\epsilon),\qquad
T_{i,j}^{n+1}=(T^{(0)})_{j}^{n+1}+\mathcal{O}(\epsilon),\qquad
2\rho_{j}^{n+1}=(ac(T^{(0)})^{4})_{j}^{n+1}+\mathcal{O}(\epsilon),$ (5.17)

Summing up the equations for $\rho$, $T_{e}$, and $T_{i}$ in (4.12b) and
employing (5.5), we can derive that

$\displaystyle\frac{a((T_{j}^{(0)})^{4})_{j}^{n+1}-a((T_{j}^{(0)})^{4})_{j}^{n+1,\ast}}{\Delta
t}+\frac{C_{v_{e}}(T^{(0)}_{j})^{n+1}-C_{v_{e}}(T^{(0)}_{j})^{n+1,\ast}}{\Delta
t}+\frac{C_{v_{i}}(T^{(0)}_{j})^{n+1}-C_{v_{i}}(T^{(0)}_{j})^{n+1,\ast}}{\Delta
t}$ (5.18) $\displaystyle\qquad=\frac{1}{\Delta
x}\left(\frac{D_{e}((T^{(0)})_{j+1}^{n+1})+D_{e}((T^{(0)})_{j}^{n+1})}{2}\frac{(T^{(0)})_{j+1}^{n+1}-(T^{(0)})_{j}^{n+1}}{\Delta
x}-\frac{D_{e}((T^{(0)})_{j}^{n+1})+D_{e}((T^{(0)})_{j-1}^{n+1})}{2}\frac{(T^{(0)})_{j}^{n+1}-(T^{(0)})_{j-1}^{n+1}}{\Delta
x}\right)$ $\displaystyle\qquad+\frac{1}{\Delta
x}\left(\frac{D_{i}((T^{(0)})_{j+1}^{n+1})+D_{i}((T^{(0)})_{j}^{n+1})}{2}\frac{(T^{(0)})_{j+1}^{n+1}-(T^{(0)})_{j}^{n+1}}{\Delta
x}-\frac{D_{i}((T^{(0)})_{j}^{n+1})+D_{i}((T^{(0)})_{j-1}^{n+1})}{2}\frac{(T^{(0)})_{j}^{n+1}-(T^{(0)})_{j-1}^{n+1}}{\Delta
x}\right).$

Considering (5.16), (5.17), and (5.18), we observe that the numerical scheme
(4.12), with the spatial discretization (5.1) and (5.5), reduces to an
implicit scheme of the diffusion limit (2.9). Furthermore, (4.12a) and (4.12b)
converge to the numerical schemes corresponding to the 1D form of the split
systems (3.10) and (3.12). This convergence completes the proof.

#### AP property to the two temperature limit

As the parameter $\kappa$ tends to infinity, the numerical scheme proposed in
Sec. 5.1 converges to a consistent numerical scheme for the two-temperature
model given by (2.11).

As $\kappa$ approaches infinity, it holds that

$T_{e,j}^{n+1,\ast}=T_{j}^{n+1,\ast}+\mathcal{O}(\kappa^{-1}),\quad
T_{i,j}^{n+1,\ast}=T_{j}^{n+1,\ast}+\mathcal{O}(\kappa^{-1}).$ (5.19)

With $C_{v}=C_{v_{e}}+C_{v_{i}}$, by summing up the governing equation for
$T_{e}$ and $T_{i}$ in (4.12a), we can derive

$\displaystyle\frac{2\epsilon^{2}}{c}\frac{\rho_{j}^{n+1,\ast}-\rho_{j}^{n}}{\Delta
t}+\frac{\mathcal{F}^{n+1,\ast}_{1,j+1/2}-\mathcal{F}^{n+1,\ast}_{1,j-1/2}}{\Delta
x}=-\frac{\sigma_{j}^{n+1,\ast}}{2}(2\rho-acT^{4})_{j}^{n+1,\ast},$ (5.20)
$\displaystyle\frac{\epsilon^{2}}{c}\frac{\psi_{1,j}^{n+1,\ast}-\psi_{1,j}^{n}}{\Delta
t}+\frac{\mathcal{F}_{0,j+1/2}^{n+1,\ast}-\mathcal{F}_{0,j-1/2}^{n+1,\ast}}{\Delta
x}+\frac{\mathcal{F}_{2,j+1/2}^{n+1,\ast}-\mathcal{F}_{2,j-1/2}^{n+1,\ast}}{\Delta
x}=-\sigma_{j}^{n+1,\ast}\psi_{1,j},$
$\displaystyle\frac{\epsilon^{2}}{c}\frac{\bar{\bm{u}}_{j}^{n+1,\ast}-\bar{\bm{u}}_{j}^{n}}{\Delta
t}+\frac{{\bf F}^{n+1,\ast}_{1,j+1/2}-{\bf F}^{n+1,\ast}_{1,j-1/2}}{\Delta
x}=-\sigma_{j}^{n+1,\ast}\bar{\bm{u}}_{j}^{n+1,\ast},$
$\displaystyle\epsilon^{2}C_{v}\frac{T_{j}^{n+1,\ast}-T_{j}^{n}}{\Delta
t}=\frac{\sigma_{j}^{n+1,\ast}}{2}(2\rho-acT^{4})_{j}^{n+1,\ast}.$

Here, $\bf F$ is a vector with $M-1$ entries, and the $l$-th entry represents
the corresponding numerical flux of $\psi_{l+1},l=1,\cdots,M-1$. When $\kappa$
approach infinity in (4.12b), one obtains

$T_{e,j}^{n+1}=T_{j}^{n+1}+\mathcal{O}(\kappa^{-1}),\quad
T_{i,j}^{n+1}=T_{j}^{n+1}+\mathcal{O}(\kappa^{-1}).$ (5.21)

By summing up the governing equations for $T_{e}$ and $T_{i}$ in (4.12b), we
have

$\displaystyle\frac{2\epsilon^{2}}{c}\frac{\rho_{j}^{n+1}-\rho_{j}^{n+1,\ast}}{\Delta
t}=-\frac{\sigma_{j}^{n+1}}{2}\left(2\rho-acT^{4}\right)_{j}^{n+1},$ (5.22)
$\displaystyle\epsilon^{2}C_{v}\frac{T_{j}^{n+1}-T_{j}^{n+1,\ast}}{\Delta
t}=\frac{\epsilon^{2}}{\Delta
x}\left(\frac{D_{T,j+1}^{n+1}+D_{T,j}^{n+1}}{2}\frac{T_{j+1}^{n+1}-T_{j}^{n+1}}{\Delta
x}\right.$
$\displaystyle\qquad\qquad\qquad\qquad\left.-\frac{D_{T,j}^{n+1}+D_{T,j-1}^{n+1}}{2}\frac{T_{j}^{n+1}-T_{j-1}^{n+1}}{\Delta
x}\right)+\frac{\sigma^{n+1}}{2}\left(2\rho-acT^{4}\right)_{j}^{n+1},$

with $C_{v}=C_{v_{e}}+C_{v_{i}}$, and $D_{T}=D_{e}(T)+D_{i}(T)$. It is evident
that (5.20) and (5.22) together constitute a splitting scheme of the two-
temperature model (2.11), thus completing the proof.

###  5.3 Conservation of total energy

In this section, we will discuss the energy-conserving property of this
numerical scheme. For the 3TRTE model (4.1), it conserves the total energy, as
demonstrated by the theorem below.

#### Energy-conserving property

For the 3TRTE system (4.1), with the periodic boundary conditions in the
spatial space and the absence of external source terms, the system conserves
total energy. This conservation is precisely expressed by

$\frac{\partial}{\partial t}E_{\rm total}(t)=0,$ (5.23)

where the total energy $E_{\rm total}(t)$ is defined as

$E_{\rm total}(t)=\int
E_{e}(t,x)+E_{i}(t,x)+E_{r}(t,x)\mathop{}\\!\mathrm{d}x.$ (5.24)

The proof is straightforward by integrating (4.11) in the spatial domain, and
we will omit it. For the numerical scheme (4.12), utilizing the spatial
discretization methods presented in (5.1) and (5.5), it effectively conserves
the discretized total energy. This conservation is expressed as follows.

#### Conservation of the discretized total energy

The numerical scheme (4.12), with the spatial discretization given by (5.1)
and (5.5), conserves the total energy

$\mathcal{E}_{\rm total}^{n+1}=\mathcal{E}_{\rm total}^{n}$ (5.25)

for the 3TRTE system (4.1) with periodic boundary condition. The discretized
total energy at time $t=t^{n}$ is defined as

$\mathcal{E}_{\rm
total}^{n}=\sum_{j=1}^{N}(E_{e,j}^{n}+E_{i,j}^{n}+E_{r,j}^{n}),$ (5.26)

with $N$ the total number of mesh points in the spatial domain.

To prove this, first, by summing up the governing equation for $\rho$,
$T_{e}$, and $T_{i}$ in (4.12), with the spatial discretization (5.1) and
(5.5), it holds that

$\displaystyle\frac{2\epsilon^{2}}{c}\frac{\rho_{j}^{n+1}-\rho_{j}^{n}}{\Delta
t}+\epsilon^{2}C_{v_{e}}\frac{T_{e,j}^{n+1}-T_{e,j}^{n}}{\Delta
t}+\epsilon^{2}C_{v_{i}}\frac{T_{i,j}^{n+1}-T_{i,j}^{n}}{\Delta t}$ (5.27)
$\displaystyle\qquad\qquad=-\frac{\mathcal{F}_{1,j+1/2}^{n+1,\ast}-\mathcal{F}_{1,j-1/2}^{n+1,\ast}}{\Delta
x}$ $\displaystyle\qquad\qquad+\frac{\epsilon^{2}}{\Delta
x}\left(D^{n+1}_{e,j+1/2}\frac{T^{n+1}_{e,j+1}-T^{e,n+1}_{e,j}}{\Delta
x}-D^{n+1}_{e,j-1/2}\frac{T^{n+1}_{e,j}-T^{e,n+1}_{e,j-1}}{\Delta x}\right)$
$\displaystyle\qquad\qquad+\frac{\epsilon^{2}}{\Delta
x}\left(D^{n+1}_{i,j+1/2}\frac{T^{n+1}_{i,j+1}-T^{i,n+1}_{i,j}}{\Delta
x}-D^{n+1}_{i,j+1/2}\frac{T^{n+1}_{i,j+1}-T^{i,n+1}_{i,j}}{\Delta x}\right),$

with

$D^{n+1}_{s,j+1/2}=\frac{1}{2}\left(D^{n+1}_{s,j+1}+D^{n+1}_{s,j}\right),\qquad
s=e,i.$ (5.28)

Summing up over the spatial space, we obtain

$\displaystyle\sum_{j=1}^{N}\left(\frac{2\epsilon^{2}}{c}\frac{\rho_{j}^{n+1}-\rho_{j}^{n}}{\Delta
t}+\epsilon^{2}C_{v_{e}}\frac{T_{e,j}^{n+1}-T_{e,j}^{n}}{\Delta
t}+\epsilon^{2}C_{v_{i}}\frac{T_{i,j}^{n+1}-T_{i,j}^{n}}{\Delta t}\right)$
(5.29)
$\displaystyle\qquad\qquad=\frac{\mathcal{F}_{1,1/2}^{n+1,\ast}-\mathcal{F}_{1,N+1/2}^{n+1,\ast}}{\Delta
x}$ $\displaystyle\qquad\qquad+\frac{\epsilon^{2}}{\Delta
x}\left(D^{n+1}_{e,N+1/2}\frac{T^{n+1}_{e,N+1}-T^{n+1}_{e,N}}{\Delta
x}-D^{n+1}_{e,1/2}\frac{T^{n+1}_{e,1}-T^{n+1}_{e,0}}{\Delta x}\right)$
$\displaystyle\qquad\qquad+\frac{\epsilon^{2}}{\Delta
x}\left(D^{n+1}_{i,N+1/2}\frac{T^{n+1}_{i,N+1}-T^{n+1}_{i,N}}{\Delta
x}-D^{n+1}_{i,1/2}\frac{T^{n+1}_{i,1}-T^{n+1}_{i,0}}{\Delta x}\right).$

With the periodic boundary condition, we have

$\mathcal{F}_{1,1/2}^{n+1,\ast}=\mathcal{F}_{1,N+1/2}^{n+1,\ast},\qquad
D^{n+1}_{s,1/2}=D^{n+1}_{s,N+1/2},\qquad T^{n+1}_{s,0}=T^{n+1}_{s,N},\qquad
T^{n+1}_{s,1}=T^{n+1}_{s,N+1},\qquad s=e,i.$ (5.30)

Substituting (5.30) into (5.29) yields

$\sum_{j=1}^{N}\left(\frac{2\epsilon^{2}}{c}\frac{\rho_{j}^{n+1}-\rho_{j}^{n}}{\Delta
t}+\epsilon^{2}C_{v_{e}}\frac{T_{e,j}^{n+1}-T_{e,j}^{n}}{\Delta
t}+\epsilon^{2}C_{v_{i}}\frac{T_{i,j}^{n+1}-T_{i,j}^{n}}{\Delta t}\right)=0.$
(5.31)

Thus, we have completed the proof.

###  5.4 The whole algorithm

In this section, for the completeness of the numerical scheme, we provide an
outline of the numerical algorithm as follows.

  1. 1.

Let $t=0$ and $n=0$, and calculate the initial value of $\rho_{j}^{0}$
,$\psi_{j}^{0}$, $\bar{\bm{u}}^{0}_{j}$, $T_{e,j}^{0}$, $T_{i,j}^{0}$;

  2. 2.

Set time step $\Delta t$ with the CFL condition (6.2);

  3. 3.

Solve (4.12) to obtain $\rho_{j}^{n+1},\psi_{j}^{n+1},\bar{\bm{u}}^{n+1}_{j}$,
$T_{e,j}^{n+1}$, $T_{i,j}^{n+1}$:

    1. (a)

Let $k=0$;

    2. (b)

Set
$\rho_{j}^{n+1,\ast,2k}=\rho_{j}^{n},\psi_{1,j}^{n+1,\ast,2k}=\psi_{1,j}^{n},\bar{\bm{u}}^{n+1,\ast,2k}_{j}=\bar{\bm{u}}_{j}^{n}$,
$T_{e,j}^{n+1,\ast,2k}=T_{e,j}^{n},T_{i,j}^{n+1,\ast,2k}=T_{i,j}^{n}$ and
$\sigma_{j}^{n+1,\ast,2k}=\sigma_{j}^{n}$;

    3. (c)

Obtain
$\rho_{j}^{n+1,\ast,2k+1},T_{e,j}^{n+1,\ast,2k+1},T_{i,j}^{n+1,\ast,2k+1}$
using (4.14) and (5.1);

    4. (d)

Update $\sigma_{j}^{n+1,\ast,2k+1}$ using $T_{e,j}^{n+1,\ast,2k+1}$;

    5. (e)

Obtain $\psi_{1,j}^{n+1,\ast,2k+1},\bar{\bm{u}}^{n+1,\ast,2k+1}_{j}$ using
(4.14) and (5.1);

    6. (f)

Obtain
$\rho_{j}^{n+1,\ast,2k+2},\psi_{1,j}^{n+1,\ast,2k+2},T_{e,j}^{n+1,\ast,2k+2},T_{i,j}^{n+1,\ast,2k+2}$
using (4.15) and (5.1);

    7. (g)

Update $\sigma_{j}^{n+1,\ast,2k+2}$ using $T_{e,j}^{n+1,\ast,2k+2}$;

    8. (h)

Obtain $\bar{\bm{u}}^{n+1,\ast,2k+2}_{j}$ using (4.15) and (5.1);

    9. (i)

Calculate the error $\mathcal{E}$ using (6.1a);

    10. (j)

If $\mathcal{E}^{n+1,\ast,2k+2}<\bar{\epsilon}$, set
$\rho_{j}^{n+1,\ast}=\rho_{j}^{n+1,\ast,2k+2},\psi_{1,j}^{n+1,\ast}=\psi_{1,j}^{n+1,\ast,2k+2},\bar{\bm{u}}^{n+1,\ast}_{j}=\bar{\bm{u}}_{j}^{n+1,\ast,2k+2}$,
$T_{e,j}^{n+1,\ast}=T_{e,j}^{n+1,\ast,2k+2}$ and
$T_{i,j}^{n+1,\ast}=T_{i,j}^{n+1,\ast,2k+2}$, else let $k\leftarrow k+1$, and
return to step (c);

    11. (k)

Set
$\psi_{1,j}^{n+1}=\psi_{1,j}^{n+1,*},\bar{\bm{u}}^{n+1}_{j}=\bar{\bm{u}}_{j}^{n+1,*}$.

    12. (l)

Let $k=0$;

    13. (m)

Set
$\rho_{j}^{n+1,2k}=\rho_{j}^{n+1,\ast},\bar{\bm{u}}^{n+1,2k}_{j}=\bar{\bm{u}}_{j}^{n+1,\ast}$,
$T_{e,j}^{n+1,2k}=T_{e,j}^{n+1,\ast},T_{i,j}^{n+1,2k}=T_{i,j}^{n+1,\ast}$ and
$\sigma_{j}^{n+1,2k}=\sigma_{j}^{n+1,\ast}$;

    14. (n)

Obtain
$\rho_{j}^{n+1,2k+1},\bar{\bm{u}}^{n+1,2k+1}_{j},T_{e,j}^{n+1,2k+1},T_{i,j}^{n+1,2k+1}$
using (4.19) and (5.5);

    15. (o)

Update $\sigma_{j}^{n+1,2k+1}$ and $D_{s}^{n+1,2k+1},s=e,i$ using
$T_{e,j}^{n+1,2k+1}$ and $T_{i,j}^{n+1,2k+1}$;

    16. (p)

Obtain
$\rho_{j}^{n+1,2k+2},\bar{\bm{u}}^{n+1,2k+2}_{j},T_{e,j}^{n+1,2k+2},T_{i,j}^{n+1,2k+2}$
using (4.20) and (5.5);

    17. (q)

Update $\sigma_{j}^{n+1,2k+2}$ and $D_{s}^{n+1,2k+2},s=e,i$ using
$T_{e,j}^{n+1,2k+2}$ and $T_{i,j}^{n+1,2k+2}$;

    18. (r)

Calculate the error $\mathcal{E}$ using (6.1b);

    19. (s)

If $\mathcal{E}^{n+1,2k+2}<\bar{\epsilon}$, set
$\rho_{j}^{n+1}=\rho_{j}^{n+1,2k+2}$, $T_{e,j}^{n+1}=T_{e,j}^{n+1,2k+2}$ and
$T_{i,j}^{n+1}=T_{i,j}^{n+1,2k+2}$, else let $k\leftarrow k+1$, and return to
step (n);

  4. 4.

Let $n\leftarrow n+1$, and return to step 2.

##  6 Numerical examples

In this section, the validation of the numerical method is conducted through a
series of tests. First, the AP property of this splitting scheme is verified.
Then, numerical examples including spatially homogeneous conditions, the 1D
Marshak wave problem, and the 2D Riemann problem are investigated. In the
first test, the periodic condition is applied, while the inflow boundary
condition in [10] is applied for the last three tests.

###  6.1 Stopping criteria

We discuss the stopping criteria for each iteration to solve (4.12). The key
stopping criterion is based on the error measure, denoted as $\mathcal{E}$,
which ensures the accuracy of the solution. Two forms of this error,
$\mathcal{E}^{n+1,\ast.2k+2}$ and $\mathcal{E}^{n+1,2k+2}$, are adopted to
monitor the convergence of the iterative process. The detailed form is defined
as

$\displaystyle\mathcal{E}^{n+1,\ast,2k+2}$ $\displaystyle=\frac{1}{\Delta
t}\left\\{\Delta
x\sum_{j=1}^{N}\left[\sum_{l=1}^{M}\frac{2}{2l+1}\left(\psi_{l,j}^{n+1,\ast,2k+2}-\psi_{l,j}^{n+1,\ast,2k}\right)^{2}+2\left(2\rho_{j}^{n+1,\ast,2k+2}-2\rho_{j}^{n+1,\ast,2k}\right)^{2}\right.\right.$
(6.1a)
$\displaystyle\qquad\left.\left.+\left(T_{e,j}^{n+1,\ast,2k+2}-T_{e,j}^{n+1,\ast,2k}\right)^{2}+\left(T_{i,j}^{n+1,\ast,2k+2}-T_{i,j}^{n+1,\ast,2k}\right)^{2}\right]\right\\}^{1/2};$
$\mathcal{E}^{n+1,2k+2}=\frac{1}{\Delta t}\left\\{\Delta
x\sum_{j=1}^{N}\left[\left(T_{e,j}^{n+1,2k+2}-T_{e,j}^{n+1,2k}\right)^{2}+\left(T_{i,j}^{n+1,2k+2}-T_{i,j}^{n+1,2k}\right)^{2}\right]\right\\}^{1/2}.$
(6.1b)

The parameter $\bar{\epsilon}$ is introduced as a stopping threshold for these
errors. It is a problem-specific value and set to $10^{-10}$ for both 1D and
2D problems in our simulation. The CFL (Courant-Friedrichs-Lewy) condition is
utilized to determine the time step length for this splitting method, which is
calculated as

$\Delta t=C\Delta x/c,$ (6.2)

where $\Delta x$ is the spatial mesh size, $c$ is the speed of light, and $C$
is the CFL number.

###  6.2 Tests of the AP property

We first investigate the AP property of this numerical scheme. The initial
conditions and parameters are similar to examples in [10], which is an
equilibrium as

$T_{r}=T_{e}=T_{i}=(3+\sin(\pi x))/4,\qquad x\in[0,2].$ (6.3)

Periodic condition is applied, with other parameters set as

$a=1,\quad c=1,\quad C_{v_{e}}=0.1,\quad C_{v_{i}}=0.2,\quad\sigma=10,\quad
K_{e}=0.01,\quad K_{i}=0.02.$ (6.4)

The AP property in the radiation diffusion limit with the parameter $\kappa$
fixed at $\kappa=1$ is first studied. The expansion order in the $P_{N}$
method is set to $M=7$. During the simulation, the mesh size $N$ is varied,
with values set as $50$, $100$, $200$, $400$, and $800$. The reference
solution is obtained using the same method with a mesh size of $N=1600$. The
CFL number is set to $C=0.1$, and the time step length is determined as (6.2).
The $l_{2}$ error between the numerical solution and the reference solution at
$t=0.1$ is depicted in Fig. 1, where $\epsilon$ is set to $1,0.1$, and
$0.001$. The results in Fig. 1 demonstrate that, for all three temperatures,
the $l_{2}$ error converges at the first order for all values of $\epsilon$,
confirming the AP property of this splitting scheme concerning $\epsilon$. In
Fig 2, the AP property related to $\kappa$ is investigated, with $\epsilon$
fixed at $\epsilon=1$, and the other parameters identical as (6.4). Fig. 2
displays the $l_{2}$ error between the numerical solution and the reference
solution at $t=0.1$, with $\kappa$ taking values of $\kappa=1,10$, and $100$.
The results in Fig. 2 demonstrate that, for all three temperatures, the
$l_{2}$ error converges at the first order for all values of $\kappa$,
providing evidence for the AP property of this splitting scheme concerning
$\kappa$.

(a) $T_{r}$ (b) $T_{e}$ (c) $T_{i}$

Figure 1: (Tests of the AP property in Sec. 6.2) The $l_{2}$ error between the
numerical solution with mesh size $N=50,100,200,400$ and $800$ and the
reference solution. Here, the reference solution is obtained by the same
numerical method with $N=1600$. The parameter $\epsilon$ is set as
$\epsilon=1,0.1$ and $0.001$ and $\kappa=1$.

(a) $T_{r}$ (b) $T_{e}$ (c) $T_{i}$

Figure 2: (Tests of the AP property in Sec. 6.2) The $l_{2}$ error between the
numerical solution with grid size $N=50,100,200,400$ and $800$ and the
reference solution. Here, the reference solution is obtained by the same
numerical method with $N=1600$. The parameter $\kappa$ is set as $\kappa=1,10$
and $100$ and $\epsilon=1$.

To investigate the impact of the CFL number on the efficiency of the splitting
method (4.12), we examine the iteration numbers for both steps in the
numerical scheme under various CFL numbers. Initially, with $\kappa=1$, we
test three different values of $\epsilon$ as $\epsilon=1,0.1$, and $0.001$.
Tab. 1 displays the iteration numbers for both steps of the numerical scheme
(4.12) when the CFL number $C=0.01,0.05$, and $0.1$, with a fixed mesh size of
$N=200$. The results reveal that, for the first part of the numerical scheme,
iteration numbers are consistently small across different CFL numbers and all
three values of $\epsilon$. In the second part of the scheme, the iteration
numbers are consistently small with the CFL number for all tested $\epsilon$
values too.

Subsequently, with $\epsilon=1$ held constant, we explore the effect of three
different $\kappa$ values as $\kappa=1,10$, and $100$. Tab. 2 illustrates the
iteration numbers for both steps of the numerical scheme (4.12) at the CFL
number $C=0.01,0.05$, and $0.1$, with a fixed mesh size of $N=200$. Similar to
the $\epsilon$ scenario, the iteration numbers for the first part of the
numerical scheme consistently stay small across various CFL numbers and all
three $\kappa$ values. In the second part of the scheme, iteration numbers
also consistently remain small.

In summary, across all tested regimes of $\epsilon$ and $\kappa$, the
iteration numbers for both steps in the numerical scheme (4.12) are
consistently small. Importantly, the CFL number shows minimal influence on the
iteration numbers, indicating the high efficiency of this asymptotic-
preserving splitting numerical method.

iteration number | first iteration | second iteration
---|---|---
$C=0.01$ | $C=0.05$ | $C=0.1$ |  $C=0.01$ | $C=0.05$ | $C=0.1$
$\epsilon=1$ | 6 | 6 | 6 | 6 | 8 | 8
$\epsilon=0.1$ | 6 | 8 | 8 | 8 | 8 | 8
$\epsilon=0.001$ | 8 | 8 | 8 | 8 | 8 | 8
Table 1: (Tests of the AP property in Sec. 6.2) Iteration numbers for the both steps in the numerical scheme (4.12) with CFL number $C=0.01,0.05$ and $0.1$ for $\epsilon=1,0.1$, and $0.001$ and $\kappa=1$.  iteration number | first iteration | second iteration
---|---|---
$C=0.01$ | $C=0.05$ | $C=0.1$ |  $C=0.01$ | $C=0.05$ | $C=0.1$
$\kappa=1$ | 6 | 6 | 6 | 6 | 8 | 8
$\kappa=10$ | 6 | 6 | 6 | 6 | 8 | 8
$\kappa=100$ | 6 | 6 | 6 | 6 | 8 | 8
Table 2: (Tests of the AP property in Sec. 6.2) Iteration numbers for the both
steps in the numerical scheme (4.12) with CFL number $C=0.01,0.05$ and $0.1$
for $\kappa=1,10$, and $100$ and $\epsilon=1$.

###  6.3 Homogeneous model problem

In this section, we examine two homogeneous model problems, each involving a
source affecting either the ions or the photons, respectively. In this case,
the system becomes an ODE system with three equations with strong stiffness.
We can verify that the splitting scheme can provide an accurate solution even
if the ODE is stiff. The parameters for these tests are consistent with those
in references like [5, 6, 8]. Detailed parameters for both tests are provided
in Tab. 3. The source term is defined as

$Q=\frac{\bar{\rho}\mathcal{C}}{\sqrt{2\pi}t_{w}}\exp\left(\frac{(t-t_{c})^{2}}{2t_{w}^{2}}\right)[{\rm
GJcm^{-3}ns^{-1}}],$ (6.5)

and the initial conditions for both problems are as follows.

$T_{i}=T_{e}=T_{r}=2.52487\times 10^{-5}{\rm keV}.$ (6.6)

| Problem I | Problem II
---|---|---
$Q_{i}$ |  (6.5) | $0$
$Q_{r}$ | $0$ |  (6.5)
$C_{v_{e}}$ | $0.1\bar{\rho}$ | $0.1T_{e}\bar{\rho}$
$C_{v_{i}}$ | $0.05\bar{\rho}$ | $0.05\bar{\rho}$
$\kappa$ | $C_{v_{e}}/(c\tau)$ | $0.01379T_{e}^{-0.5}$
$c~{}({\rm cmns}^{-1})$ | $29.979$ | $29.979$
$a~{}({\rm GJcm}^{-3}{\rm keV}^{4})$ | $0.01372$ | $0.01372$
$\tau~{}({\rm ns})$ | $0.1$ | $0.1$
$\epsilon$ | $1$ | $1$
$\mathcal{C}$ | 25.06628 | 25.06628
$t_{w}~{}({\rm ns})$ | 1.0 | 1.0
$t_{c}~{}({\rm ns})$ | 1.0 | 1.0
$\bar{\rho}~{}({\rm gcm}^{-3})$ | 3.0 | 3.0
$\sigma~{}({\rm cm}^{-1})$ | $\sigma_{0}/T_{e}^{2}$ | $\sigma_{0}/T_{e}^{2}$
$\sigma_{0}$ | 0.5 | 0.5
$M$ | 7 | 7
$\Delta t({\rm ns})$ | $0.0025$ | $0.0025$
Table 3: (Homogeneous model problem in Sec. 6.3) Parameters for the
homogeneous problem with source terms.

(a) Problem I

(b) Problem II

Figure 3: (Homogeneous model problem in Sec. 6.3) The time evolution of the
numerical solution, the reference solution by the same method with a time step
length of $\Delta t=0.01\times 2^{-10}$ns, and the reference solution from
[8]. The dot-dashed lines are the numerical solution by the AP $P_{N}$ method.
The solid lines are the reference solution by the same method with a time step
length of $\Delta t=0.01\times 2^{-10}$ns, and the symbols are the reference
solutions from [8]. (a) Solutions of Problem I, (b) Solutions of Problem II.

#### Problem I: Ion-source homogeneous problem

In Problem I, the source term of the ion temperature is introduced. The source
term is specified in (6.5), and the initial condition is presented in (6.6).
Additionally, the opacity depends on the temperature of electron $T_{e}$ and
follows the relation $\sigma=\frac{1}{2T_{e}^{2}}.$

During the simulation, the expansion number for the $P_{N}$ method is set to
$M=7$, and the time step length is defined as $\Delta t=0.0025$ns. The
simulation is run until the final time $t=20$ns. The time evolution of
$T_{r}$, $T_{e}$, and $T_{i}$ is illustrated in Fig. 3a, along with the
reference solution obtained using the same method with a time step length of
$\Delta t=0.01\times 2^{-10}$ ns and reference solutions from [8]. The
simulation results show close agreement with the reference solutions.
Initially, all three temperatures are identical, but they gradually diverge
over time. Due to the source term affecting the ions, the ion temperature
exhibits the highest value. The $l_{2}$ error between the numerical solution
and the reference solution is presented in Fig. 4a, demonstrating first-order
convergence in time. When the implicit midpoint scheme is employed for
temporal discretization, the convergence rate improves to the second order as
shown in Fig. 4b. Here, the implicit midpoint method with a time step length
of $\Delta t=0.01\times 2^{-10}$ ns is adopted as the reference solution.

(a) First order

(b) Second order

Figure 4: (Homogeneous model problem in Sec. 6.3) The $l_{2}$ error of the
numerical solution by the AP $P_{N}$ method and the reference solution by the
same method with a time step length of $\Delta t=0.01\times 2^{-10}$ns. (a)
The error for the first-order $P_{N}$ method, (b) The error for the second-
order $P_{N}$ method.

#### Problem II: Photons-source homogeneous problem

In Problem II, the source term in (6.5) is introduced for the photons, and the
identical initial condition (6.6) is employed. The specific parameters are
detailed in Tab. 3, where the specific heat of electrons $C_{v_{e}}$ depends
on the electron temperature and is defined as

$C_{v_{e}}=0.1\bar{\rho}T_{e},\qquad\bar{\rho}=3.0,$ (6.7)

and the coefficient $\kappa$ depends on $T_{e}$ as well

$\kappa=0.01379T_{e}^{-0.5}.$ (6.8)

In the simulation, the expansion number for the $P_{N}$ method is set to
$M=7$, and the time step length is configured as $\Delta t=0.0025$ns. The time
evolution of the radiative temperature $T_{r}$, electron temperature $T_{e}$,
and ion temperature $T_{i}$ is illustrated in Fig. 3b, alongside the reference
solutions obtained using the same method with a time step length of $\Delta
t=0.01\times 2^{-10}$ns and the reference solution from [5], respectively. The
simulation is conducted until the final time $t=20$ns. The results indicate
that the numerical solution matches the reference solution well. In this
problem, all three temperatures are nearly identical initially, with the
radiative temperature $T_{r}$ gradually becoming the highest over time due to
the influence of the source term.

###  6.4 1D Marshak wave problem without conduction terms

In this section, we investigate the Marshak problem without conduction terms.
This problem is a classical benchmark extensively studied in [6, 10, 28]. For
this specific example, the computational domain is $[0,0.5]$ and the
equilibrium initial condition is chosen as

$T_{e}=T_{i}=10^{-6},\quad\psi=\frac{1}{2}acT_{e}^{4}.$ (6.9)

The inflow boundary conditions are defined as

${\rm
left~{}boundary:}~{}T_{e}=T_{i}=1,\quad\psi=\frac{1}{2}acT_{e}^{4},\qquad{\rm
right~{}boundary:}~{}T_{e}=T_{i}=10^{-6},\quad\psi=\frac{1}{2}acT_{e}^{4},$
(6.10)

where the parameters are set as follows.

$\epsilon=1,\qquad a=0.01372,\qquad c=299.79,\qquad C_{v_{e}}=0.03,\qquad
C_{v_{i}}=0.27,\qquad\sigma=300T_{e}^{-3}.$ (6.11)

Detailed application of the inflow boundary condition within the framework of
the $P_{N}$ method is proposed in [10].

(a) $T_{r},\kappa=1$ (b) $T_{e},\kappa=1$ (c) $T_{i},\kappa=1$

Figure 5: (1D Marshak wave problem without conduction terms in Sec. 6.4) The
numerical solution of $T_{r},T_{e}$ and $T_{i}$ for $\kappa=1$ at different
times. The symbols are the numerical solution by the AP $P_{N}$ method, and
the solid lines are the reference solution obtained from the method in [6].

In the absence of the conduction terms, both $D_{e}$ and $D_{i}$ are set to
zero. A grid size of $N=200$ is employed, along with linear reconstruction. We
use $M=7$ in the $P_{N}$ method. The numerical solution for the radiative
temperature $T_{r}$, as well as the temperature of the electron and ion
$T_{e},T_{i}$ at time $t=0.074,0.37$, and $0.74$ for $\kappa=1$, is depicted
in Fig. 5. Additionally, the reference solution obtained using the method from
[6] is included. The results illustrate that all three temperatures propagate
forward and closely align with the reference solution.

(a) $T_{r},\kappa=0.001$ (b) $T_{e},\kappa=0.001$ (c) $T_{i},\kappa=0.001$

(d) $T_{r},\kappa=100$ (e) $T_{e},\kappa=100$ (f) $T_{i},\kappa=100$

Figure 6: (1D Marshak wave problem without conduction terms in Sec. 6.4) The
numerical solution of $T_{r},T_{e}$ and $T_{i}$ for $\kappa=0.001$ and $100$
at different times. The top column is for $\kappa=0.001$ and the bottom column
is for $\kappa=100$. The symbols are the numerical solution by the AP $P_{N}$
method, and the solid lines are the reference solution obtained from the
method in [6].

To demonstrate the AP property of the scheme proposed in Sec. 5.4, the
behavior of the numerical solution is examined for two additional values of
$\kappa$: $\kappa=0.001$ and $100$, as depicted in Fig. 6. For both small
$\kappa$ ($\kappa=0.001$) and large $\kappa$ ($\kappa=100$), Fig. 6 reveals
that the numerical solution closely aligns with the reference solution.
Importantly, it is clearly illustrated that the behavior of the ion
temperature $T_{i}$ differs from that of $T_{e}$ at $\kappa=0.001$, while the
two behave similarly for $\kappa=100$. For a more precise comparison of the
behavior of $T_{e}$ and $T_{i}$ for different $\kappa$, the numerical
solutions of $T_{e}$ and $T_{i}$ for various $\kappa$ at $t=0.74$ are
presented in Fig. 7. It is observed that when $\kappa=0.001$ in Fig. 7a, the
two temperatures exhibit significant differences. However, as $\kappa$
increases to $\kappa=1$, the numerical solutions of $T_{e}$ and $T_{i}$ become
nearly identical. Furthermore, the numerical solutions for these two
temperatures completely overlap when $\kappa$ increases to $100$.
Additionally, the same time step length, which is independent of $\kappa$, is
employed for all three cases, indicating the AP property of this new method
concerning $\kappa$.

(a) $\kappa=0.001$

(b) $\kappa=1$

(c) $\kappa=100$

Figure 7: (1D Marshak wave problem without conduction terms in Sec. 6.4) The
comparison of $T_{e}$ and $T_{i}$ for different $\kappa$ at different times.
The symbols are the numerical solution of $T_{e}$ and the solid lines are the
numerical solution of $T_{i}$. (a) $\kappa=0.001$, (b) $\kappa=1$, (c)
$\kappa=100$.

###  6.5 Marshak wave problem with conduction terms

In this section, the Marshak wave problem with conduction terms is examined.
The same initial and boundary conditions as (6.9) and (6.10) are utilized. The
parameters are set as

$\epsilon=1,\quad a=0.01372,\quad c=29.979,\quad C_{v_{e}}=0.3,\quad
C_{v_{i}}=0.27,\quad\sigma=300T_{e}^{-3}.$ (6.12)

Additionally, the conduction coefficients are set to $K_{e}=1.0$ and
$K_{i}=0.1$.

(a) $T_{r},\kappa=1$ (b) $T_{e},\kappa=1$ (c) $T_{i},\kappa=1$

Figure 8: (1D Marshak wave problem with conduction terms in Sec. 6.5) The
numerical solution of $T_{r},T_{e}$ and $T_{i}$ for $\kappa=1$ at different
times. The symbols are the numerical solution by the AP $P_{N}$ method, and
the solid lines are the reference solution obtained from the $S_{N}$ method.

In the simulation, the same grid size, expansion order, and time step length
as in Sec. 6.4 are adopted. The numerical solution of the radiative
temperature $T_{r}$, and the temperature of electron and ion $T_{e},T_{i}$ for
$\kappa=1$ at $t=0.01,0.05$ and $0.1$ is plotted in Fig. 8. It shows that the
numerical solution and the reference solution are almost on top of each other.
Here the reference solution is obtained by the $S_{N}$ method. For the $S_{N}$
method, the fully implicit scheme is applied with the mesh size in the spacial
space $N=1600$, and the CFL number $C=0.8$. The upwind scheme is applied to
solve the convection term, with $8$ Gaussian nodes utilized to discrete the
angular space. Comparing Fig. 5 and 8, the evolution of these three
temperatures is quite different with and without the conduction terms. With
the conduction terms, the evolution of these three temperatures is much
smoother, aligning with the analytical properties of the system (2.1).
Furthermore, significant differences between $T_{e}$ and $T_{i}$ are noticed
when $\kappa=1$.

(a) $T_{r},\kappa=0.001$ (b) $T_{e},\kappa=0.001$

(c) $T_{i},\kappa=0.001$

(d) $T_{r},\kappa=100$ (e) $T_{e},\kappa=100$ (f) $T_{i},\kappa=100$

Figure 9: (1D Marshak wave problem with conduction terms in Sec. 6.5) The
numerical solution of $T_{r},T_{e}$ and $T_{i}$ for $\kappa=0.001$ and $100$
at different times. The top column is for $\kappa=0.001$ and the bottom column
is for $\kappa=100$. The symbols are the numerical solution by the AP $P_{N}$
method, and the solid lines are the reference solution obtained from the
$S_{N}$ method.

To demonstrate the AP property of this scheme, we have conducted tests for two
additional values of $\kappa$ as $\kappa=0.001$ and $100$. The numerical
solutions for the three temperatures at $t=0.01,0.05$, and $0.1$ are presented
in Fig. 9, alongside the reference solution obtained by the $S_{N}$ method.
For both values of $\kappa$, the numerical solutions closely align with the
reference solutions, demonstrating the AP property of the numerical scheme
introduced in Sec. 5.4. Furthermore, Fig. 9b and 9c provide a clear depiction
of the different behaviors between electrons and ions for small $\kappa$.
However, as $\kappa$ increases, the behaviors of electrons and ions become
more similar. To facilitate a precise comparison of their behavior, Fig. 10
displays the temperature profiles of electrons and ions for different values
of $\kappa$. It is evident that when $\kappa=0.001$, the behaviors of the two
temperatures differ significantly. In contrast, when $\kappa=1$, the
differences diminish considerably, and as $\kappa$ increases to $100$, the
behaviors of these two temperatures nearly coincide. Notably, the same time
step length, independent of $\kappa$, is employed for all three cases,
indicating the AP property of this new method with respect to $\kappa$.

(a) $\kappa=0.001$

(b) $\kappa=1$

(c) $\kappa=100$

Figure 10: (1D Marshak wave problem with conduction terms in Sec. 6.5) The
comparison of $T_{e}$ and $T_{i}$ for different $\kappa$ at different times.
The symbols are the numerical solution of $T_{e}$ and the solid lines are the
numerical solution of $T_{i}$. (a) $\kappa=0.001$, (b) $\kappa=1$, (c)
$\kappa=100$.

###  6.6 2D Riemann problem

(a) Problem setting for the 2D Riemann problem

(b) Problem setting for the 1D Riemann problem

Figure 11: (2D Riemann problem in Sec. 6.6) Computational regime and problem
setting for the 2D Riemann problem. Here the blue and red line in Fig. 11b is
the computational regime in the corresponding 1D problem.

In this section, a 2D Riemann problem is examined using a similar setup to the
one described in [31]. This problem is defined on a square spatial domain
$[-1.5,1.5]\times[-1.5,1.5]$ as shown in Fig. 11a. In this example, the
constant parameters are set as follows.

$a=1,\quad c=1,\quad C_{v_{e}}=1,\quad C_{v_{i}}=0.1.$ (6.13)

As in Fig. 11, the opacity $\sigma$ is chosen based on the spatial region:

$\sigma(x,y)=\left\\{\begin{aligned} 10,\quad&x+y>0\text{ and }x<y,\quad\text{
(white region)},\\\ 1,\quad&x+y<0\text{ or }x>y,\quad\text{(black and gray
regions)}.\end{aligned}\right.$ (6.14)

Moreover, the initial conditions for the material temperatures are given as

$T_{e}(0,x,y)=T_{i}(0,x,y)=\left\\{\begin{aligned}
1,\quad&x+y<0,\quad\text{(black region)},\\\ 0.1,\quad&x+y>0,\quad\text{(gray
and white regions)}.\end{aligned}\right.$ (6.15)

The boundary conditions are set as

$T_{e}(t,x,y)=T_{i}(t,x,y)=\left\\{\begin{aligned} 1,\quad&x=-1.5\text{ or
}y=-1.5,\\\ 0.1,\quad&x=1.5\text{ or }y=1.5,\end{aligned}\qquad\forall
t>0.\right.$ (6.16)

The initial and boundary radiative intensity is in thermal equilibrium

$\begin{gathered}\psi(0,x,y,\Omega)=\frac{1}{4\pi}acT_{r}^{4}(0,x,y),\\\
\psi(t,x,y,\Omega)=\frac{1}{4\pi}acT_{r}^{4}(t,x,y),\quad(x,y)\in\partial[-1.5,1.5]^{2},\end{gathered}$
(6.17)

with $T_{r}(0,x,y)=T_{e}(0,x,y)$ given in (6.15), and
$T_{r}(t,x,y)=T_{e}(t,x,y)$ for boundary condition.

(a) $T_{r},P_{N}$ (b) $T_{e},P_{N}$ (c) $T_{i},P_{N}$

(d) $T_{r},S_{N}$ (e) $T_{e},S_{N}$ (f) $T_{i},S_{N}$

Figure 12: (2D Riemann problem in Sec. 6.6) The radiative temperature $T_{r}$
and the temperature of electron and ion $T_{e},T_{i}$ at $t=1$ for $\kappa=1$.
The first column is the numerical solution obtained by the AP $P_{N}$ method,
and the second column is the reference solution obtained by the $S_{N}$
method.

In the simulation, the $P_{N}$ method employs an expansion order of $M=30$,
and the grid size is set to $N_{x}=N_{y}=128$ with WENO reconstruction [17] in
B. First, we consider the case $\kappa=1$. The radiative temperature $T_{r}$,
electron temperature $T_{e}$, and ion temperature $T_{i}$ at $t=1$ are
presented in Fig. 12, where the reference solution is obtained using the
$S_{N}$ method with a mesh size of $256\times 256$. In the $S_{N}$ method, the
angular space is discreted with $1202$ Lebedev points (order $59$ polynomial),
and the upwind scheme is utilized for the convection term with the time step
length (6.2) with $\rm CFL=0.1$. It is evident that the numerical solution
closely matches the reference solution. Moreover, in the white region, the
radiative temperature is lower than in the gray region, consistent with the
opacity value $\sigma$. To provide a more detailed comparison between the
numerical and reference solutions, we plot the numerical solution along the
lines $y=x-0.5$ and $y=x+0.5$ in Fig. 13, demonstrating a close agreement.
Additionally, to verify the numerical solution more carefully, two one-
dimensional problems along the regime shown in Fig. 11b are tested with the
initial and boundary conditions the same as the 2D problem, where the same
$P_{N}$ method is utilized with the expansion number $M=30$, and the grid size
$N=1600$. The numerical results are shown in Fig. 14, where the behavior of
the numerical solution for the 1D problem is almost the same as the 2D
numerical solution along $y=x-0.5$ and $y=x+0.5$. This indicates that a
smaller opacity $\sigma$ accelerates the propagation of radiation, yet
decelerates the heating of the material by the radiation. Note that in the
problems considered in this section, the conduction term is absent, meaning
the change in material temperature is solely influenced by radiation.
Therefore, for a small opacity $\sigma$, the material temperature rises over a
large spatial range, but the increase is relatively small. This shows that the
numerical results are reasonable.

######  Remark 2.

In this example, the interface for the one-dimensional test is located at the
cell boundary, and the interface for the two-dimensional test passes through
the center of the cell.

To investigate the impact of $\kappa$ on the temperatures $T_{e}$ and $T_{i}$,
the numerical solutions for different $\kappa$ along the same lines $y=x-0.5$
and $y=x+0.5$ are presented in Fig. 15. When $\kappa=0$, the system reduces to
a two-temperature limit, where $T_{i}$ remains the same as the initial
condition and significantly differs from $T_{e}$ at $t=1$, as shown in Fig.
15a. As $\kappa$ increases, the variation in $T_{i}$ becomes more pronounced.
As illustrated in Fig. 15b and 15c, the numerical solution of $T_{i}$
gradually converges towards that of $T_{e}$. When $\kappa$ reaches $\kappa=1$,
the discrepancy between these two numerical solutions is small. Moreover, at
$\kappa=10$ and $100$, these two solutions nearly coincide. Fig. 15
demonstrates that the behavior of temperatures $T_{e}$ and $T_{i}$ aligns with
the AP property of the system (2.1) as $\kappa$ increases.

(a) $T_{r},\kappa=1$ (b) $T_{e},\kappa=1$ (c) $T_{i},\kappa=1$

Figure 13: (2D Riemann problem in Sec. 6.6) The radiative temperature $T_{r}$
and the temperature of electron and ion $T_{e},T_{i}$ at $t=1$ along $y=x-0.5$
and $y=x+0.5$ for $\kappa=1$. The solid line is the numerical solution
obtained by the AP $P_{N}$ method, and the symbol is the reference solution
obtained by the $S_{N}$ method.

(a) $T_{r},\kappa=1$ (b) $T_{e},\kappa=1$ (c) $T_{i},\kappa=1$

Figure 14: (2D Riemann problem in Sec. 6.6) Comparision of the numerical
solution with the 1D problem. The radiative temperature $T_{r}$ and the
temperature of electron and ion $T_{e},T_{i}$ at $t=1$ for $\kappa=1$. The
solid line is the numerical solution obtained by the AP $P_{N}$ method along
$y=x-0.5$ and $y=x+0.5$, and the symbol is the numerical solution obtained by
the 1D problem.

(a) $\kappa=0$

(b) $\kappa=0.01$

(c) $\kappa=0.1$

(d) $\kappa=1$

(e) $\kappa=10$

(f) $\kappa=100$

Figure 15: (2D Riemann problem in Sec. 6.6) The temperature $T_{e},T_{i}$ at
$t=1$ along $y=x-0.5$ and $y=x+0.5$ for $\kappa=0,0.01,0.1,1,10$ and $100$.
The solid line is the numerical solution of the temperature $T_{e}$, and the
symbol is the numerical solution of the temperature $T_{i}$.

To verify the AP property of this numerical scheme, the numerical solution is
examined for an extremely large $\kappa$ value, specifically $\kappa=1000$.
The same expansion order for the $P_{N}$ method, mesh size, and time step size
are utilized. The numerical solutions of $T_{r},T_{e}$, and $T_{i}$ at $t=1$
are presented in Fig. 16, along with the reference solution obtained using the
$S_{N}$ method. As depicted in Fig. 16, the numerical solution closely aligns
with the reference solution, even with $\kappa=1000$. Furthermore, Fig. 17
displays the numerical solution along $y=x-0.5$ and $y=x+0.5$. The comparison
with the reference solution further illustrates the AP property of the
numerical scheme introduced in Sec. 5.4.

(a) $T_{r},P_{N}$ (b) $T_{e},P_{N}$ (c) $T_{i},P_{N}$

(d) $T_{r},S_{N}$ (e) $T_{e},S_{N}$ (f) $T_{i},S_{N}$

Figure 16: (2D Riemann problem in Sec. 6.6) The radiative temperature $T_{r}$
and the temperature of electron and ion $T_{e},T_{i}$ at $t=1$ for
$\kappa=1000$. The first column is the numerical solution obtained by the AP
$P_{N}$ method, and the second column is the reference solution obtained by
the $S_{N}$ method.

(a) $T_{r},\kappa=1000$ (b) $T_{e},\kappa=1000$ (c) $T_{i},\kappa=1000$

Figure 17: (2D Riemann problem in Sec. 6.6) The radiative temperature $T_{r}$
and the temperature of electron and ion $T_{e},T_{i}$ at $t=1$ along $y=x-0.5$
and $y=x+0.5$ for $\kappa=1000$. The solid line is the numerical solution
obtained by the AP $P_{N}$ method, and the symbol is the reference solution
obtained by the $S_{N}$ method.

##  7 Conclusion

In this work, we have introduced an AP numerical method for solving the three-
temperature radiative transfer model. This method utilizes a splitting scheme
to ensure the accuracy and stability of the solution in both the diffusion
limit as $\epsilon$ approaches zero and the two-temperature limit as $\kappa$
tends to infinity. Additionally, the numerical scheme has been demonstrated to
possess the AP property and conserve energy. We have presented various
numerical examples, including homogeneous model problems, the 1D Marshak wave
problem, and the 2D Riemann problems, to validate the AP property and the
efficiency of the proposed numerical scheme. Future work may involve exploring
more complex numerical tests and applications of this method in different
physical contexts such as the frequency-dependent problems [33]. Furthermore,
this splitting method can be generalized to provide an asymptotic-preserving
coupling method of radiation and fluid, as well as radiation and plasma.

## Acknowledges

We thank Prof. Peng Song from the Institute of Applied Physics and
Computational Mathematics and Prof. Tao Xiong from Xiamen University for their
valuable suggestions. The work of Ruo Li is partially supported by the
National Natural Science Foundation of China (Grant No. 12288101). W. Li is
partially supported by the National Natural Science Foundation of China (Grant
No. 12001051). M. Tang is partially supported by the National Natural Science
Foundation of China (Grant No. 12031013), Shanghai pilot innovation project
21JC1403500, and the Strategic Priority Research Program of Chinese Academy of
Sciences Grant No. XDA25010401. This work of Y. Wang is partially supported by
the Foundation of the President of China Academy of Engineering Physics
(YZJJZQ2022017) and the National Natural Science Foundation of China (Grant
No. 12171026, U2230402, and 12031013).

##  Appendix A Order analysis

In this section, we present an order analysis of $\psi_{l}$, where
$l=0,1,\cdots$. A similar analysis is available in [10], and we included it
here for the completeness of the algorithm. As stated in [10], it is confirmed
that for $\psi_{l}$, the following relationship holds:

$\psi_{l}=\mathcal{O}(\epsilon^{l}),\qquad l=0,1,\cdots.$ (A.1)

To prove this, the Chapman-Enskog expansion is utilized, and the specific
intensity $\psi$ is expanded in power series of $\epsilon$ as

$\psi=\psi^{(0)}+\epsilon\psi^{(1)}+\epsilon^{2}\psi^{(2)}+\cdots.$ (A.2)

Define

$\psi_{l}^{(k)}=\int_{-1}^{1}P_{l}(\mu)\psi^{(k)}(t,x,\mu)\mathop{}\\!\mathrm{d}\mu,\qquad
l=0,1,\cdots.$ (A.3)

Then, it holds that

$\psi_{l}=\psi_{l}^{(0)}+\epsilon\psi_{l}^{(1)}+\epsilon^{2}\psi_{l}^{(2)}+\cdots.$
(A.4)

We will show that for any fixed $k$,

$\psi_{l}^{(k)}=0,\qquad\forall l>k.$ (A.5)

First, if $k=0$, matching the order of $\mathcal{O}(1)$ in (4.7a) we can
obtain

$\psi^{(0)}=acT_{e}^{4}.$ (A.6)

Then, by the orthogonality of the spherical harmonic functions, it holds that

$\psi_{l}^{(0)}=0,\qquad\forall l\geqslant 1.$ (A.7)

Hence, for the case of $k=0$, (A.5) remains valid. Subsequently, assuming that
(A.5) holds for all $k$ up to and including $n$, as stated

$\psi_{l}^{(k)}=0,\qquad\forall k\leqslant n,\qquad l>k.$ (A.8)

When considering the case where $k=n+1$, it is established that for $\psi_{l}$
with $l>0$ that

$\sigma\psi_{l}=\frac{\epsilon^{2}}{c}\frac{\partial\psi_{l}}{\partial
t}+\frac{\epsilon l}{2l+1}\frac{\partial\psi_{l-1}}{\partial
x}+\frac{\epsilon(l+1)}{2l+1}\frac{\partial\psi_{l+1}}{\partial x}.$ (A.9)

Matching the order of $\mathcal{O}(\epsilon^{n+1})$ yields

$\sigma\psi_{l}^{(n+1)}=\frac{\epsilon^{2}}{c}\frac{\partial\psi_{l}^{(n-1)}}{\partial
t}+\frac{\epsilon l}{2l+1}\frac{\partial\psi_{l-1}^{(n)}}{\partial
x}+\frac{\epsilon(l+1)}{2l+1}\frac{\partial\psi_{l+1}^{(n)}}{\partial x}.$
(A.10)

According to (A.8), it can be observed that $\psi_{l}^{(n+1)}=0$ for values of
$l$ greater than $n+1$. With this established, we have successfully
demonstrated (A.5), and consequently, (A.1) can be derived naturally.

##  Appendix B WENO reconstruction

The 3-order WENO reconstruction [17]:

$\displaystyle f^{L,1}=\frac{3}{2}f_{j}-\frac{1}{2}f_{j-1},\quad
f^{L,2}=\frac{1}{2}f_{j}+\frac{1}{2}f_{j+1},$ (B.1) $\displaystyle
f^{R,1}=\frac{3}{2}f_{j}-\frac{1}{2}f_{j+1},\quad
f^{R,2}=\frac{1}{2}f_{j}+\frac{1}{2}f_{j-1},$
$\displaystyle\omega_{L,1}=\frac{1/3}{\left[10^{-6}+(f_{j}-f_{j-1})^{2}\right]^{2}},\quad\omega_{L,2}=\frac{2/3}{\left[10^{-6}+(f_{j+1}-f_{j})^{2}\right]^{2}},$
$\displaystyle\omega_{R,1}=\frac{1/3}{\left[10^{-6}+(f_{j+1}-f_{j})^{2}\right]^{2}},\quad\omega_{R,2}=\frac{2/3}{\left[10^{-6}+(f_{j}-f_{j-1})^{2}\right]^{2}},$
$\displaystyle
f_{j+1/2}^{L}=\frac{\omega_{L,1}f^{L,1}+\omega_{L,2}f^{L,2}}{\omega_{L,1}+\omega_{L,2}},\quad
f_{j-1/2}^{R}=\frac{\omega_{R,1}f^{R,1}+\omega_{R,2}f^{R,2}}{\omega_{R,1}+\omega_{R,2}}.$

Then, the numerical flux can be obtained from $f_{j+1/2}^{L}$ and
$f_{j-1/2}^{R}$.

## References

  * [1] H. An, X. Jia, and H.F. Walker.  Anderson acceleration and application to the three-temperature energy equations.  Journal of Computational Physics, 347:1–19, 2017.
  * [2] H. An, Z. Mo, X. Xu, and X. Jia.  Operator-based preconditioning for the 2-D 3-T energy equations in radiation hydrodynamics simulations.  Journal of Computational Physics, 385:51–74, 2019.
  * [3] W. Bennett and R. McClarren.  Self-similar solutions for high-energy density radiative transfer with separate ion and electron temperatures.  Proceedings of the Royal Society A, 477:20210119, 2021.
  * [4] S. Boscarino, J.M. Qiu, G. Russo, and T. Xiong.  A high order semi-implicit IMEX WENO scheme for the all-Mach isentropic Euler system.  Journal of Computational Physics, 392:594–618, 2019.
  * [5] R. Chauvin, S. Guisset, B. Manach-Perennou, and L. Martaud.  A colocalized scheme for three-temperature grey diffusion radiation hydrodynamics.  Communications in Computational Physics, 31(04), 2022.
  * [6] C. Enaux, S. Guisset, C. Lasuen, and Q. Ragueneau.  Numerical resolution of a three temperature plasma model.  Journal of Scientific Computing, 82(3):1–29, 2020.
  * [7] C. Enaux, S. Guisset, C. Lasuen, and G. Samba.  Numerical methods for coupling multigroup radiation with ion and electron temperatures.  Communications in Applied Mathematics and Computational Science, 17:43–78, 2022.
  * [8] T.M. Evans and J.D. Densmore.  Methods for coupling radiation, ion, and electron energies in grey Implicit Monte Carlo.  Journal of Computational Physics, 225(2):1695–1720, 2007.
  * [9] J.A. Fleck Jr and J.D. Cummings Jr.  An implicit Monte Carlo scheme for calculating time and frequency dependent nonlinear radiation transport.  Journal of Computational Physics, 8(3):313–342, 1971.
  * [10] J. Fu, W. Li, P. Song, and Y. Wang.  An asymptotic-preserving IMEX method for nonlinear radiative transfer equation.  Journal of Scientific Computing, 92:1–35, 2022.
  * [11] J. Fu, W. Li, X. Tao, and Y. Wang.  Solving nonlinear radiative transfer equation using IMEX PN method.  Under Preparation, 2023.
  * [12] S. Ganesan and M. Singh.  An operator-splitting finite element method for the numerical solution of radiative transfer equation.  arXiv preprint arXiv:2122.07949, 2021.
  * [13] P.H. Hauschildt and E. Baron.  A 3D radiative transfer framework I. Non-local operator splitting and continuum scattering problems.  Astronomy & Astrophysics, 451:273–284, 2006.
  * [14] J. Jiang, Y. Huang, S. Shu, and S. Zeng.  Some new discretization and adaptation and multigrid methods for 2-D 3-T diffusion equations.  Journal of Computational Physics, 224:168–181, 2007.
  * [15] V. Laboure, R. McClarren, and C. Hauck.  Implicit filtered $P_{N}$ for high-energy density thermal radiation transport using discontinuous galerkin finite elements.  Journal of Computational Physics, 321:624–643, 2016.
  * [16] E.E. Lewis and W.F. Miller.  Computational methods of neutron transpot.  American Nuclear Society, La Grange Park, Illinois, 1 1984.
  * [17] X.D. Liu, S. Osher, and T. Chan.  Weighted essentially non-oscillatory schemes.  Journal of Computational Physics, 115:200–212, 1994.
  * [18] R. McClarren, T. Evans, R. Lowrie, and J. Densmore.  Semi-implicit time integration for PN thermal radiative transfer.  Journal of Computational Physics, 227(16):7561–7586, 2008.
  * [19] R. McClarren, J.P. Holloway, and T.A. Brunner.  On solutions to the Pn equations for thermal radiative transfer.  Journal of Computational Phyiscs, 227:2864–2885, 2008.
  * [20] R. McClarren and J. Wöhlbier.  Solutions for ion-electron-radiation coupling with radiation and electron diffusion.  Journal of Quantitative Spectroscopy and Radiative Transfer, 112:1119–130, 2011.
  * [21] Z. Mo, L. Shen, and G. Wittum.  Parallel adaptive multigrid algorithm for 2-D 3-T diffusion equations.  International Journal of Computer Mathematics, 81(3):361–374, 2004\.
  * [22] G. Peng, Z. Gao, W. Yan, and X. Feng.  A positivity-preserving finite volume scheme for three-temperature radiation diffusion equations.  Applied Numerical Mathematics, 152:125–140, 2020.
  * [23] Y. Shi, X. Han, W. Sun, and P. Song.  A continuous source tilting scheme for radiative transfer equations in implicit Monte Carlo.  Journal of Computational and Theoretical Transport, 50(1), 2021\.
  * [24] Y. Shi and H. Yong.  A maximum-principle preserving implicit Monte Carlo method for a three-temperature radiative transfer model.  Available at SSRN 4408635, 2023.
  * [25] S. Shiva, C. Leela, P. Kiran, C.D. Sijoy, and S. Chaturvedi.  The effects of electron thermal radiation on laser ablative shock waves from aluminum plasma into ambient air.  Physics of Plasmas, 23:053107, 2016.
  * [26] C.D. Sijoy and S. Chaturvedi.  TRHD: Three-temperature radiation-hydrodynamics code with an implicit non-equilibrium radiation transport using a cell-centered monotonic finite volume scheme on unstructured-grids.  Computer Physics Communications, 190:98–119, 2015.
  * [27] C.D. Sijoy and S. Chaturvedi.  Combining node-centered parallel radiation transport and higher-order multi-material cell-centered hydrodynamics methods in three-temperature radiation hydrodynamics code TRHD.  Computer Physics Communications, 203:94–109, 2016.
  * [28] W. Sun, S. Jiang, and K. Xu.  An asymptotic preserving unified gas kinetic scheme for gray radiative transfer equations.  Journal of Computational Physics, 285(15):265–279, 2015.
  * [29] 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.
  * [30] M. Tang and X. Zhang.  Semi-implicit front capturing schemes for the degenerate nonlinear radiative diffusion equation.  Journal of Computational Physics, 436:110290, 2021.
  * [31] 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, 463:111308, 2022.
  * [32] 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.
  * [33] X. Zhang, P. Song, Y. Shi, and M. Tang.  A fully asymptotic preserving decomposed multi-group method for the frequency-dependent radiative transfer equations.  Journal of Computational Physics, 491:112368, 2023.