The equilibrium-diffusion limit for radiation hydrodynamics

TL;DR

This paper derives and analyzes the equilibrium-diffusion approximation (EDA) for radiation hydrodynamics, demonstrating its first-order accuracy and how common approximations preserve this accuracy in different frames, with practical implications for radiative-shock solutions.

Contribution

The paper provides a rigorous asymptotic derivation of the EDA and shows that it and related approximations maintain first-order accuracy in both comoving and lab frames.

Findings

EDA is first-order accurate with second-order transport corrections.

Grey nonequilibrium-diffusion and Eddington approximations preserve first-order accuracy.

Radiative-shock solutions approach equilibrium-diffusion solutions as the asymptotic parameter diminishes.

Abstract

The equilibrium-diffusion approximation (EDA) is used to describe certain radiation-hydrodynamic (RH) environments. When this is done the RH equations reduce to a simplified set of equations. The EDA can be derived by asymptotically analyzing the full set of RH equations in the equilibrium-diffusion limit. We derive the EDA this way and show that it and the associated set of simplified equations are both first-order accurate with transport corrections occurring at second order. Having established the EDA's first-order accuracy we then analyze the grey nonequilibrium-diffusion approximation and the grey Eddington approximation and show that they both preserve this first-order accuracy. Further, these approximations preserve the EDA's first-order accuracy when made in either the comoving-frame (CMF) or the lab-frame (LF). While analyzing the Eddington approximation, we found that the CMF and LF radiation-source equations are equivalent when neglecting ${\cal O}(\beta^2)$ terms and compared in the LF. Of course, the radiation pressures are not equivalent. It is expected that simplified physical models and numerical discretizations of the RH equations that do not preserve this first-order accuracy will not retain the correct equilibrium-diffusion solutions. As a practical example, we show that nonequilibrium-diffusion radiative-shock solutions devolve to equilibrium-diffusion solutions when the asymptotic parameter is small.

Figures

Equations

$${\partial }_t\rho +{\partial }_i\left(\rho u_i\right)=0,$$

$${\partial }_t\left(\rho u_i\right)+{\partial }_j\left(\rho u_i{u}_j+p_{ij}\right)=-S_{\text{rp},i},$$

$${\partial }_{t}E+{\partial }_i\left[u_j\left(E_{ij}+p_{ij}\right)\right]=-S_{\text{re}},$$

$$\frac{1}{c}{\partial }_t{I}_{\nu }+{\Omega }_i{\partial }_i{I}_{\nu }=Q_{\nu },$$

$$\begin{gathered} Q_{\nu }=-\left(\frac{{\nu }_{\text{o}}}{\nu }\right){\sigma }_{\text{t},{\nu }_{\text{o}}}{I}_{\nu }+{\left(\frac{\nu }{{\nu }_{\text{o}}}\right)}^2{\sigma }_{\text{a},{\nu }_{\text{o}}}{B}_{{\nu }_{\text{o}}}(T) \\ +{\left(\frac{\nu }{{\nu }_{\text{o}}}\right)}^2\frac{{\sigma }_{\text{s}}}{4\pi }{\int }_{4\pi }\left(\frac{{\nu }_{\text{o}}}{{\nu }^{\prime }}\right)I_{{\nu }^{\prime }}\left({\Omega }^{\prime }\right)d{\Omega }^{\prime }, \end{gathered}$$

$$B_{{\nu }_{\text{o}}}(T)=\frac{2h{\nu }_{\text{o}}^3}{c^2}{\left[e^{h{\nu }_{\text{o}}/k_{\text{B}}T}-1\right]}^{-1},$$

$$S_{\text{re}}\equiv {\int }_{4\pi }{\int }_0^{\infty }{Q}_{\nu }d\Omega d\nu ={\partial }_t\mathcal{E}+{\partial }_i{\mathcal{F}}_i,$$

$$\begin{gathered} S_{\text{rp},i}\equiv \frac{1}{c}{\int }_{4\pi }{\int }_0^{\infty }{\Omega }_i{Q}_{\nu }d\Omega d\nu \\ =\frac{1}{c^2}{\partial }_t{\mathcal{F}}_i+{\partial }_j{\mathcal{P}}_{ij}, \end{gathered}$$

$$\mathcal{E}\equiv \frac{1}{c}{\int }_{4\pi }{\int }_0^{\infty }{I}_{\nu }d\Omega d\nu ,$$

$${\mathcal{F}}_i\equiv {\int }_{4\pi }{\int }_0^{\infty }{\Omega }_i{I}_{\nu }d\Omega d\nu ,$$

$${\mathcal{P}}_{ij}\equiv \frac{1}{c}{\int }_{4\pi }{\int }_0^{\infty }{\Omega }_i{\Omega }_j{I}_{\nu }d\Omega d\nu ,$$

$$\mathcal{E}=a_{\text{R}}{T}^4,$$

$${\mathcal{F}}_i=-\frac{a_{\text{R}}c}{3{\sigma }_{\text{t},\text{R}}}{\partial }_i{T}^4+\frac{4}{3}{u}_i{a}_{\text{R}}{T}^4,$$

$${\mathcal{P}}_{ij}=\frac{1}{3}{a}_{\text{R}}{T}^4{\delta }_{ij},$$

$${\partial }_t\rho +{\partial }_i\left(\rho u_i\right)=0,$$

$$\begin{gathered} {\partial }_t\left(\rho u_i\right) \\ +{\partial }_j\left(\rho u_i{u}_j+p_{ij}+\frac{1}{3}{a}_{\text{R}}{T}^4{\delta }_{ij}\right)=0, \end{gathered}$$

$$\begin{gathered} {\partial }_t\left(E+\mathcal{E}\right) \\ +{\partial }_i\left[u_j\left(E_{ij}+p_{ij}+\frac{4}{3}{a}_{\text{R}}{T}^4{\delta }_{ij}\right)\right]= \\ {\partial }_i\left(\frac{a_{\text{R}}c}{3{\sigma }_{\text{t},\text{R}}}{\partial }_i{T}^4\right), \end{gathered}$$

$$x=\hat{x}{l}_{\infty },t=\hat{t}\frac{l_{\infty }}{a_{\infty }},u=\hat{u}{u}_{\infty },$$

$$\rho =\hat{\rho }{\rho }_{\infty },p=\hat{p}{\rho }_{\infty }{a}_{\infty }^2,e=\hat{e}{a}_{\infty }^2,$$

$$T=\hat{T}{T}_{\infty },{\sigma }_{\text{t}}={\hat{\sigma }}_{\text{t}}{\sigma }_{\text{t},\infty },{\sigma }_{\text{s}}={\hat{\sigma }}_{\text{s}}{\sigma }_{\text{s},\infty },$$

$$I_{\nu }={\hat{I}}_{\hat{\nu }}\frac{a_{\text{R}}ch{T}_{\infty }^3}{k_{\text{B}}},\nu =\hat{\nu }\frac{k_{\text{B}}{T}_{\infty }}{h},$$

$$\mathcal{E}=\hat{\mathcal{E}}{a}_{\text{R}}{T}_{\infty }^4,$$

$${\mathcal{F}}_i={\hat{\mathcal{F}}}_i{a}_{\text{R}}c{T}_{\infty }^4,$$

$${\mathcal{P}}_{ij}={\hat{\mathcal{P}}}_{ij}{a}_{\text{R}}{T}_{\infty }^4.$$

$${\partial }_t\rho +\mathcal{M}{\partial }_i\left(\rho u_i\right)=0,$$

$$\begin{gathered} \mathcal{M}{\partial }_t\left(\rho u_i\right)+{\partial }_j\left({\mathcal{M}}^2\rho u_i{u}_j+p_{ij}\right)= \\ -P_0\left(\mathcal{U}{\partial }_t{\mathcal{F}}_i+{\partial }_j{\mathcal{P}}_{ij}\right), \end{gathered}$$

$$\begin{gathered} {\partial }_t\left(\frac{1}{2}{\mathcal{M}}^2\rho u^2+\rho e\right) \\ +{\partial }_i\left[\mathcal{M}{u}_j\left(\frac{1}{2}{\mathcal{M}}^2\rho u_i{u}_j+\rho e{\delta }_{ij}+p_{ij}\right)\right] \\ =-P_0\left({\partial }_t\mathcal{E}+\mathcal{C}{\partial }_i{\mathcal{F}}_i\right), \end{gathered}$$

$$\mathcal{M}\equiv \frac{u_{\infty }}{a_{\infty }},P_0\equiv \frac{a_{\text{R}}{T}_{\infty }^4}{{\rho }_{\infty }{a}_{\infty }^2},$$

$$\mathcal{U}\equiv \frac{a_{\infty }}{c},\mathcal{C}\equiv \frac{c}{a_{\infty }},$$

$$\mathcal{L}\equiv \frac{l_{\infty }}{{\lambda }_{\text{t},\infty }}=l_{\infty }{\sigma }_{\text{t},\infty },{\mathcal{L}}_{\text{s}}\equiv \frac{{\sigma }_{\text{s},\infty }}{{\sigma }_{\text{t},\infty }},$$

$$\mathcal{M}=\mathcal{O}(1),P_0=\mathcal{O}(1),$$

$$\mathcal{U}=\mathcal{O}(ϵ),\mathcal{C}=\mathcal{O}({ϵ}^{-1}),$$

$$\mathcal{L}=\mathcal{O}({ϵ}^{-1}),{\mathcal{L}}_{\text{s}}=\mathcal{O}(ϵ).$$

$${\partial }_t\rho +{\partial }_i\left(\rho u_i\right)=0,$$

$$\begin{gathered} {\partial }_t\left(\rho u_i\right)+{\partial }_j\left(\rho u_i{u}_j+p_{ij}+{\mathcal{P}}_{ij}\right)= \\ -ϵ\frac{1}{c^2}{\partial }_t{\mathcal{F}}_i, \end{gathered}$$

$${\partial }_t\left(E+\mathcal{E}\right)+{\partial }_i\left[u_j\left(E_{ij}+p_{ij}\right)\right]=-\frac{1}{ϵ}{\partial }_i{\mathcal{F}}_i,$$

$$\rho =\sum _{n=0}^{\infty }{\rho }^{(n)}{ϵ}^n.$$

$${\left[{\partial }_t\rho +{\partial }_i\left(\rho u_i\right)\right]}^{(0)}=0,$$

$${\left[{\partial }_t\left(\rho u_i\right)+{\partial }_j\left(\rho u_i{u}_j+p_{ij}+{\mathcal{P}}_{ij}\right)\right]}^{(0)}=0,$$

$$\begin{gathered} {\left\{{\partial }_t\left(E+\mathcal{E}\right)+{\partial }_i\left[u_j\left(E_{ij}+p_{ij}\right)\right]\right\}}^{(0)} \\ =-{\partial }_i{\mathcal{F}}_i^{(1)}, \end{gathered}$$

$${\left[{\partial }_t\rho +{\partial }_i\left(\rho u_i\right)\right]}^{(1)}=0,$$

$$\begin{gathered} {\left[{\partial }_t\left(\rho u_i\right)+{\partial }_j\left(\rho u_i{u}_j+p_{ij}+{\mathcal{P}}_{ij}\right)\right]}^{(1)} \\ =-\frac{1}{c^2}{\partial }_t{\mathcal{F}}_i^{(0)}, \end{gathered}$$

$$\begin{gathered} {\left\{{\partial }_t\left(E+\mathcal{E}\right)+{\partial }_i\left[u_j\left(E_{ij}+p_{ij}\right)\right]\right\}}^{(1)} \\ =-{\partial }_i{\mathcal{F}}_i^{(2)}. \end{gathered}$$

$${\mathcal{E}}^{(0)}={\left[a_{\text{R}}{T}^4\right]}^{(0)},$$

$${\mathcal{F}}_i^{(0)}=0,$$

Full text

††journal: Journal of Quantitative Spectroscopy and Radiative Transfer

1 Introduction

Radiation hydrodynamics (RH) describes how interactions between radiation and matter affect the thermodynamic states, and potentially, the dynamic flow characteristics of the matter-radiation system. Unfortunately, the full set of RH equations are computationally expensive and numerically difficult to solve, and various model approximations have been developed to aid their solution [1, 2, 3]. The Euler equations are typically assumed to provide a sufficient model of the material’s hydrodynamic response since the photon mean-free-path is generally much longer than the mean-free-path between material interactions, so that material viscosity and heat conduction may be neglected. For the radiation it is common to make the following assumptions, either independently or together: the system is absorption dominated, the system’s size is large compared to the photon mean-free-path, the radiation is in thermal equilibrium with the material, the radiation flux is diffusive, and the radiation pressure is isotropic. Taken together these assumptions are called the equilibrium diffusion approximation (EDA) [2]. When the EDA applies to a physical system it can be modeled by a simpler set of equations than the full set of RH equations. This simplified set of RH equations provides a reasonable description of stellar structure [4], high-temperature environments [1], fusion dominated energy sources [5, 6], a variety of astrophysical settings [2, 7, 8], and high-energy-density-physics [9].

As shown by Lowrie, Morel and Hittinger [10], the EDA can be derived via an asymptotic expansion of the RH equations, which is described as follows. First, an ansatz for the RH solution is made and the solution is expressed as an infinite power-series expansion in a small parameter, $\epsilon$. Then the RH equations are nondimensionalized. Next, nondimensional parameters that arise in the RH equations are chosen to scale by particular powers of $\epsilon$, which reflects the physical importance of those parameters, and effectively defines physical and mathematical limits. Finally, an infinite hierarchy of equations is obtained by equating all terms multiplied by the same power of $\epsilon$. This hierarchy is used to define the equilibrium-diffusion limit (EDL) since the scalings are chosen to ensure that the zeroth order RH solution satisfies the EDA. Thus, in the limit as $\epsilon$ approaches zero the RH solution transitions to the EDA solution. In this sense, the EDA is said to be accurate or exact to zeroth order in the EDL. However, it is possible for an approximation to be exact to higher than zeroth order. This has been shown to be true for the EDA by Lowrie, Morel and Hittinger [10], but they only analyzed the simplified RH equations at zeroth order and they did not determine at what order transport corrections modified the EDA. The main contribution of this paper is to show that the EDA and the simplified RH equations are both exact to first order in the EDL, and that transport corrections occur at second order. It is important to be clear that the EDA is a physical approximation to the full set of RH equations, whereas the asymptotic expansion of the RH equations is exact at each asymptotic level. If we summed over all asymptotic orders then the original transport content would be regained with no loss of information.

Asymptotic expansions have two very practical applications that are not necessarily obvious. The first is to demonstrate that analytic approximations to the full RH equations give proper accuracy in the EDL. For instance, one would expect the grey nonequilibrium-diffusion approximation of the RH equations to be “properly” accurate in the EDL. To determine if this is so, one performs an asymptotic expansion of these approximate RH equations. If this expansion yields the EDA to first order, then the solution of the approximate RH equations will properly transition to the EDA solution as $\epsilon\rightarrow 0$. Otherwise, the solution will approach the EDA solution at the wrong rate, or in the worst case, not approach it at all. The second application relates to numerical discretizations of the RH equations and their approximations. Numerical schemes are said to be “asymptotic preserving” when they give proper accuracy in an asymptotic limit at an appropriate computational cost. Numerical discretizations that are consistent, i.e., converge to the proper analytic solution as the mesh is refined, are not necessarily asymptotic preserving. To determine if a numerical discretization is asymptotic preserving, one performs the same asymptotic expansion for the discrete equations that was performed for the analytic equations. The discrete asymptotic expansion must approximate the analytic expansion in certain specific ways to be asymptotic preserving. If a consistent scheme is not asymptotic preserving, one can obtain accurate solutions in problems closely approaching the asymptotic limit, but the computational cost will generally be prohibitive. In this paper we analyze two approximations of the RH equations, but we do not analyze any discretization schemes. For future work on numerical discretizations, this paper defines the sense in which discretization schemes should, ideally, preserve the EDA. Specifically, to be asymptotic preserving, a numerical discretization must produce the EDA and a discretized set of its simplified RH equations through first order in the EDL.

Similar analyses have been performed for neutron transport [11, 12] and radiative transfer [13, 14], which have shown that their analytic diffusion approximations are first-order accurate. Simplified physical models and numerical discretizations of these theories which fail to preserve the diffusion limit typically fail to obtain accurate diffusion solutions [15, 16, 12, 11, 17, 18, 19, 20, 13, 14, 21, 22, 23, 24, 25, 26]. Therefore, we expect that simplified physical models and numerical discretizations of the RH equations should preserve this asymptotic behavior in order to reasonably obtain correct equilibrium-diffusion solutions. When a physical model or numerical discretization preserves some diffusion limit they are generally referred to as being asymptotic preserving [27]. Previous analyses [17, 18, 13, 15] for neutron transport and radiative transfer have discussed the effects of initial and boundary conditions, as well as boundary layers, on the asymptotic results. However, in this paper we restrict our analysis to the interior solution sufficiently late in time and far from any boundaries so that their effects on the analysis may be neglected. An analysis including the initial and boundary conditions, and potentially boundary layers, should be the subject of future work.

The rest of this paper is organized as follows: In Section 2, the lab-frame (LF) RH equations and the EDA are presented. In Section 3, the main result of this paper is presented: the EDA and the simplified RH equations are first-order accurate in the EDL with transport corrections occurring at second-order. The derivation of these results is given in B. The EDA is a physical limit of the RH equations and we claim that its first-order accuracy should be preserved by simplified models of the RH equations. In Section 4, the grey nonequilibrium-diffusion approximation is analyzed and shown to preserve the EDA’s first-order accuracy. Further, we show that the nonequilibrium-diffusion approximation may be applied in either the comoving-frame (CMF) or the LF and the EDA’s first-order accuracy is preserved in both cases. As a practical example we show that a particular nonequilibrium-diffusion radiative-shock solution [28] devolves to the appropriate equilibrium-diffusion solution [29] when the asymptotic parameter is small. Then, the grey Eddington approximation is analyzed and shown to also preserve the EDA’s first-order accuracy. Again, the Eddington approximation may be applied in the CMF or the LF and the EDA’s first-order accuracy is preserved in both cases. Further, we show that when applying the Eddington approximation in the CMF or the LF the radiation-source equations are equivalent when neglecting ${\cal O}(\beta^{2})$ terms. Of course, the radiation pressures are not equivalent, but they only differ by a symmetric-traceless term of ${\cal O}(\beta)$. The reason for analyzing these two approximations is because the nonequilibrium-diffusion approximation only applies to optically-thick systems, and is only slightly less restrictive than the EDA. In contradistinction, the Eddington approximation may be applied to highly rarefied systems where diffusion does not apply and which may be out of thermal equilibrium. Section 5 closes the main body of this paper with a summary and recommendations for future work. In A we derive the ${\cal O}(\beta^{2})$ LF radiation-transport (RT) equation including scattering terms. We believe this is the first time that this equation has been presented in the literature when retaining terms through ${\cal O}(\beta^{2})$. In B we use results from A to derive the radiation intensity through ${\cal O}(\epsilon^{2})$. This is used to derive the zeroth- and first-order radiation variables and the second-order radiation flux, and to show that the radiation energy and pressure contain transport corrections at second order.

2 The RH equations and the EDA

We present the LF RH equations in Subsection 2.1, as well as the radiation-energy and radiation-momentum sources and the frequency-integrated angular moments of the radiation intensity. In Subsection 2.2, we review the assumptions on which the EDA relies and then we present the LF EDA and its simplified RH equations.

2.1 The RH equations and radiation variables

The RH equations used here are the Euler equations coupled to the radiation momentum and energy sources, and the angle- and frequency-dependent RT equation:

[TABLE]

is the angle- and frequency-dependent LF radiation source. The angle-dependence resides in the frequency ratios which are relativistically exact through all orders in $\beta\equiv u/c$, as well as the angle- and frequency-dependent LF radiation intensity, $I_{\nu}=I_{\nu}(\Omega)$, and $\Omega_{i}$ is the LF photon direction of flight. The time and space derivatives, $\partial_{t}\equiv\partial/\partial t$ and $\partial_{i}\equiv\partial x^{i}$, are with respect to LF time and space, where $x^{i}\in\{x,y,z\}$, and we are using the Einstein summation convention. Further, $\rho$ is the mass density, $u_{i}$ is the material velocity, $p_{ij}$ is the material pressure, $S_{\text{rp},i}$ is the radiation-momentum source such that $-S_{\text{rp},i}$ is a material-momentum source, $E=\tfrac{1}{2}\,\rho\,u^{2}+\rho\,e$ and $E_{ij}\equiv\tfrac{1}{2}\,\rho\,u_{i}\,u_{j}+\rho\,e\,\delta_{ij}$ are defined for notational convenience, $S_{\text{re}}$ is the radiation-energy source such that $-S_{\text{re}}$ is a material-energy source, $c$ is the speed of light, and $\nu$ and $\nu_{\text{o}}$ are the LF and CMF frequencies, respectively. We leave the Planck function in the CMF,

[TABLE]

as a function of the CMF frequency, $\nu_{\text{o}}$, and the local material temperature, $T$, where $h$ is Planck’s constant and $k_{\text{B}}$ is Boltzmann’s constant. Finally, $\sigma_{\text{s}}$, $\sigma_{\text{a},\nu_{\text{o}}}$, and $\sigma_{\text{t},\nu_{\text{o}}}=\sigma_{\text{a},\nu_{\text{o}}}+\sigma_{\text{s}}$, are the scattering, absorption and total cross-sections in the CMF; $\sigma_{\text{t},\nu_{\text{o}}}$ and $\sigma_{\text{a},\nu_{\text{o}}}$ are functions of $\nu_{\text{o}}$ while $\sigma_{\text{s}}$ is independent of frequency. We omit the subscript-o from the Planck function and the cross sections, which otherwise denotes a CMF variable.

The LF radiation-energy and radiation-momentum sources, $S_{\text{re}}$ and $S_{\text{rp},i}$, are the first two frequency-integrated angular moments of $Q_{\nu}$, respectively:

[TABLE]

and

[TABLE]

are the LF radiation energy density, radiation flux, and radiation-pressure, which are the first three frequency-integrated angular moments of $I_{\nu}$, respectively. The angle- and frequency-integrals of the radiation sources (3) cannot be analytically performed since $Q_{\nu}$ (1e) contains angle- and frequency-dependent variables whose functional form is not generally known.

2.2 The Equilibrium-diffusion approximation

The EDA imposes four basic assumptions on a physical system [2]:

1. the photon mean-free-path is small compared to the size of the absorption-dominated system,
2. the matter-radiation system is in thermal equilibrium,
3. the radiation flux is diffusive, and
4. the radiation pressure is isotropic. These simplifications allow for a basic understanding of the underlying phenomena, the equations describing them, and the form their solutions might take. In the LF the EDA and its simplified set of RH equations are:

[TABLE]

[TABLE]

where $a_{\text{\tiny R}}$ is the radiation constant and $\sigma_{\text{t},\text{\tiny R}}$ is the CMF Rosseland-averaged cross section [2]. We point out that the time-derivative of the radiation flux in the total-momentum equation (6b) has been dropped in accordance with the diffusion approximation.

The strength of the EDA in solving RH problems is that the radiation variables (5) in this simplified set of equations (6) are explicit functions of the hydrodynamic variables. Thus, once an equation-of-state for the material is specified, these simplified RH equations represent a solvable system that are easier to solve than the full set of RH equations (1) [30]. To be clear, the EDA neglects most of the radiation information that is contained in the RT equation (1d), and this is why we refer to any terms beyond the EDA as transport corrections.

3 Asymptotic analysis of the RH equations

In this section, we present our main result which is that the EDA (5) and its simplified set of RH equations (6) are first-order accurate in the EDL with transport corrections beginning at second order. These results are derived in B. In Subsection 3.1, the RH variables and simplified RH equations are nondimensionalized and nondimensional parameters are defined. The EDL scalings were originally presented by Lowrie, Morel and Hittinger [10], and we use these in Subsection 3.2 to scale the nondimensional RH equations. Finally, the scaled equations are asymptotically expanded and the main result of this paper is presented in Subsection 3.3. This analysis and its results only apply to the interior portion of a RH system, and do not account for effects due to initial nor boundary conditions, nor boundary layers.

3.1 Nondimensionalization of the RH equations

Each dimensional variable is decomposed into the product of a constant dimensional variable with subscript-$\infty$, indicating a reference state, which only retains the dimension of the original variable, and a non-constant nondimensional variable with a hat, which retains the value of the original variable:

[TABLE]

where $l_{\infty}$ is the reference length of the system, $a_{\infty}$ is a reference sound speed for the fluid, $u_{\infty}$ is a reference fluid velocity, $\rho_{\infty}$ is a reference fluid mass density, and $T_{\infty}$ is a reference fluid temperature. The total and scattering reference cross sections are $\sigma_{\text{t},\infty}$ and $\sigma_{\text{s},\infty}$, respectively. The radiation intensity and frequency, $\hat{I}_{\hat{\nu}}$ and $\hat{\nu}$, are nondimensionalized in such a way that the nondimensional radiation energy density, radiation flux, and radiation pressure are consistent with their definitions (4) above:

[TABLE]

The nondimensionalized form of the Euler equations coupled to the radiation energy and momentum sources are:

[TABLE]

where the hats have been dropped for notational convenience. The nondimensional parameters are:

[TABLE]

the first four of which are used in equations (8) and the last two are from B; ${\cal M}$ is related to the Mach number of the material flow, $P_{0}$ is a measure of the influence radiation has on the material, ${\cal U}$ is the ratio of the reference sound speed to the speed of light, ${\cal C}$ is the inverse of ${\cal U}$, ${\cal L}$ is a measure of the system’s size compared to the radiation mean-free-path, and ${\cal L}_{\text{s}}$ is a measure of whether the system is absorption or scattering dominated.

3.2 Asymptotic scalings

The EDL scalings of the nondimensional parameters were originally introduced by Lowrie, Morel and Hittinger [10]:

[TABLE]

The first scaling implies that no assumption is made about whether the value of ${\cal M}$ should be small nor large. As such, the EDL supports shock-wave solutions as shown by Lowrie and Rauenzahn [29], which is a result we will use in Subsection 4.3. The second scaling similarly implies that no assumption is made about whether the system is radiation dominated or not. The third and fourth scalings are consistent with nonrelativistic hydrodynamics. The fifth scaling implies that the reference length of the system is much larger than the photon mean-free-path, so that radiation may diffuse through the system. The sixth scaling implies that absorption dominates scattering. These last two scalings are used in B. Reviewing the nondimensionalized equations (8) we see that only the radiation flux is affected by the scalings. The redimensionalized and scaled version of equations (8) is:

[TABLE]

where we have moved all terms without an $\epsilon$ to the left-hand side. No approximations have gone into these equations nor any of their variables yet.

3.3 Asymptotic expansion and results

The Euler equations are first-order accurate using scalings that agree with those above [31]. The material and radiation variables in equations (11) are now expanded in powers of $\epsilon$; as an example, the expansion for the material mass density is:

[TABLE]

Collecting equations at equal orders in $\epsilon$ forms an infinite hierarchical set. The zeroth- and first-order contributions from equations (11) are:

[TABLE]

The radiation variables ${\cal E}^{(0)}$, ${\cal F}_{i}^{(0)}$, ${\cal P}_{ij}^{(0)}$, ${\cal E}^{(1)}$, ${\cal F}_{i}^{(1)}$, ${\cal P}_{ij}^{(1)}$, and ${\cal F}_{i}^{(2)}$ are determined by analyzing the RT equation. This is done in B and we present the results here:

[TABLE]

Summing these variables through first order produces the EDA as given in equations (5):

[TABLE]

which are in agreement with the results presented in [10]. The derivation in B shows that the radiation energy density and radiation pressure contain transport corrections to the EDA at second order, and so the EDA is at most first-order accurate in the EDL. Summing the scaled equations (13) through first order, and using the results just presented (14), produces the EDA’s simplified RH equations (6):

[TABLE]

Therefore, the EDA and its simplified RH equations are first-order accurate with transport corrections occurring at second order. This is the main result of this paper. Both simplified models and numerical discretizations of the RH equations should preserve this asymptotic behavior. Models that have different asymptotic behavior may not attain the correct solution in or near the EDL, while improper discretizations may require an unreasonable computational cost. These results only apply to the interior solution of an RH problem far from any boundaries, and sufficiently late in time, so that the initial and boundary conditions may be neglected.

4 Analysis of two grey RH approximations

We have established the first-order accuracy of the EDA and its simplified RH equations. As a reminder, the assumptions comprising the EDA are listed at the beginning of Subsection 2.2. In this section, we analyze two RH approximations that are commonly used to solve a wide span of RH problems and we show that they both preserve the EDA’s first-order accuracy. The two approximations are the frequency-independent (“grey”) nonequilibrium-diffusion approximation and the grey Eddington approximation. These approximations are applied to the radiation energy and momentum sources which are presented in Subsection 4.1. The grey nonequilibrium-diffusion approximation is analyzed in Subsection 4.2. Compared to the EDA it relaxes the requirement that the radiation be in thermal equilibrium with the material, hence the moniker “nonequilibrium”. As a practical example, in Subsection 4.3, we present a radiative-shock solution from this approximation and show that the EDA solution is obtained when $\epsilon$ is small. In contrast to the nonequilibrium-diffusion approximation, the grey Eddington approximation only retains the EDA assumption that the radiation pressure is isotropic. As such, it can be used to describe environments where diffusion does not apply, or that may be out of equilibrium, or which may be highly rarefied. This approximation is analyzed in Subsection 4.4.

4.1 The grey radiation sources

The two approximations analyzed in this section are typically made in the CMF by modifying the radiation energy and momentum sources. These sources are equations (6.47) and (6.48) in [3]:

[TABLE]

where we have dropped the ${\cal O}(\beta^{2})$ terms coming from the Lagrangian derivatives. In order to compare the CMF approximation to our results we Lorentz transform the approximation to the LF and analyze it there. A discussion of how to perform Lorentz transformations is outside the scope of this paper and can be found in [2]. Although the Lorentz transformation changes the coordinate frame being considered it does not alter the physical content of the approximation. It should also be clear that it does not matter if the CMF radiation sources are scaled first or transformed first, so long as the scalings are also applied to the transformation. Finally, in order to facilitate comparison of the radiation sources from the transformed approximation with the LF sources we record here the LF radiation energy and momentum sources:

[TABLE]

which are equations (2.29) and (2.30) in [32]. The simplest way to produce the scaled radiation-source equations is to make the following replacements,

[TABLE]

and then to divide the radiation-energy source by $\epsilon^{2}$ and the radiation-momentum source by $\epsilon$. The reason for dividing the radiation sources by $\epsilon$ and $\epsilon^{2}$ is to keep these sources at the corresponding order of their associated hydrodynamic sources.

4.2 The nonequilibrium-diffusion approximation

The nonequilibrium-diffusion approximation modifies the radiation-source equations by imposing the Eddington approximation, ${\cal P}_{\text{o},ij}=\tfrac{1}{3}\,{\cal E}_{\text{o}}\,\delta_{ij}$ or ${\cal P}_{\text{o},ii}={\cal E}$, dropping the time-derivative of the radiation flux and all acceleration terms, and dropping all other terms in the radiation-momentum source until Fick’s First Law of Diffusion is obtained. Applied to the CMF radiation sources (17) it produces the following simplified source equations:

[TABLE]

Solving the radiation-momentum source for the radiation flux and plugging this into the radiation-energy source produces:

[TABLE]

which are equations (97.69) and (97.70) in [2]. Lorentz transformation of the sources (19) along with the CMF Eddington approximation produces:

[TABLE]

The ${\cal O}(\beta)$ term on the radiation pressure (21c) is symmetric and traceless, so the Eddington approximation is retained after the Lorentz transformation, ${\cal P}_{ii}={\cal E}$, since the trace of a traceless object is zero. We quickly compare the transformed sources, (21a) and (21b), with the LF sources (18). If the nonequilibrium-diffusion approximation were applied to the LF radiation sources (18) then the result would differ from equations (21a) and (21b) by terms of ${\cal O}(\beta^{2})$ and $\beta_{i}\,{\cal F}_{i}$, and the radiation pressure would neglect the symmetric traceless portion contained in equation (21c). However, we will find that these differences between the CMF and LF sources, as well as the radiation pressures, are negligible in the EDL. Scaling the radiation sources, (21a) and (21b), gives:

[TABLE]

The radiation-energy source (22a) provides the solution for the scaled radiation energy density,

[TABLE]

From these expressions (23) the zeroth- and first-order solutions are obtained along with the second-order radiation flux:

[TABLE]

The zeroth- and first-order contributions from the scaled radiation sources (22) are:

[TABLE]

Summing the zeroth- and first-order results for the radiation variables (24) reproduces the EDA (5):

[TABLE]

As a reminder, the radiation sources are coupled to the Euler equations; see equations (1). Summing the zeroth- and first-order contributions from the radiation sources (25), coupled to the Euler equations (1), and using the results in equations (24), reproduces the EDA’s simplified RH equations (6):

[TABLE]

Thus, the nonequilibrium-diffusion approximation preserves the EDA’s first-order accuracy, and this result holds if the approximation is made in the CMF or the LF.

4.3 Confirmation of the asymptotic analysis by radiative-shock solutions

In this subsection we show that a particular solution of the nonequilibrium-diffusion approximation transitions to the equilibrium-diffusion solution when $\epsilon$ is small. Specifically, we briefly analyze the nonequilibrium-diffusion radiative-shock solution method developed by Lowrie and Edwards [28] and we show that when $\epsilon=0.001$ that method produces the equilibrium-diffusion solution developed by Lowrie and Rauenzahn [29].

The 1D nondimensional equations solved by Lowrie and Edwards, scaled in the EDL, are:

[TABLE]

The fourth equation (28d) is the radiation internal energy source, $S_{\text{rie}}=S_{\text{re}}-\beta\,S_{\text{rp}}$. For a monatomic ideal-gas the adiabatic index is $\gamma=5/3$, and ${\cal M}_{0}$ represents the initial Mach number of the unshocked gas, which we set to be ${\cal M}_{0}=3$ in this example. The nondimensional radiation flux for a purely absorbing system is

[TABLE]

When $\epsilon=1$ we obtain the nonequilibrium-diffusion solution and when $\epsilon=0.001$ we obtain the equilibrium-diffusion solution from the same solution method; see Figure 1. For comparison within the figure, we also include the equilibrium-diffusion solution produced by the method described by Lowrie and Rauenzahn, which is the dash-dotted line labeled as the EDA solution. This EDA solution and the nonequilibrium-diffusion solution using $\epsilon=0.001$ cannot be distinguished. Thus, these shock problems can serve as a test problem to ensure that codes are asymptotic preserving. However, other test problems are also needed that are relevant to different portions of the RH problem space.

4.4 The Eddington approximation

In this subsection we analyze the grey Eddington approximation applied to the CMF radiation sources (17). The application of this approximation produces the following CMF radiation sources:

[TABLE]

Lorentz transformation of these sources along with the Eddington approximation to the LF produces:

[TABLE]

Again, the trace of the radiation pressure returns the Eddington approximation, ${\cal P}_{ii}={\cal E}$, since the ${\cal O}(\beta)$ term is traceless, as well as symmetric. These source equations are identical to the LF radiation sources (18). Therefore, the analysis below applies equally well when the Eddington approximation is applied to the LF sources. Scaling these radiation sources, (31a) and (31b), gives:

[TABLE]

The radiation-energy source (32a) provides the solution for the scaled radiation energy density,

[TABLE]

From these expressions (33) the zeroth- and first-order solutions are obtained along with the second-order radiation flux:

[TABLE]

The zeroth- and first-order contributions from the scaled radiation sources (32) are:

[TABLE]

Summing the zeroth- and first-order results for the radiation variables (34) reproduces the EDA (5):

[TABLE]

Summing the zeroth- and first-order contributions from the radiation sources (35), coupled to the Euler equations (1), and using the results (34), reproduces the EDA’s simplified RH equations (6):

[TABLE]

Thus, the Eddington approximation preserves the EDA’s first-order accuracy. Further, it does not matter whether the Eddington approximation is applied to the CMF or to the LF radiation sources, the EDA’s first-order accuracy is preserved in both situations.

5 Summary

In this work we have derived the EDA from the RH equations via an asymptotic analysis. Our derivation showed that the EDA and its simplified RH equations are first-order accurate and that transport corrections begin at second order. Since the EDA is a physical limit of the full set of RH equations it is expected that simplified models of the RH equations should preserve the EDA’s first-order accuracy. We analyzed the grey nonequilibrium-diffusion approximation and the grey Eddington approximation and we showed that they both preserved the EDA’s first-order accuracy. These approximations can be made in the CMF or the LF, and we have shown that the EDA’s first-order accuracy is preserved in both cases. We also presented a test problem in which an equilibrium-diffusion solution was captured from a nonequilibrium-diffusion solver when $\epsilon$ was small. Other test problems that apply to different RH regimes are needed. Our results are in agreement with previous asymptotic analyses for neutron transport [11, 12] and radiative transfer [14]. Other analyses [17, 18, 13, 15] for neutron transport and radiative transfer have discussed the effects of initial and boundary conditions, as well as boundary layers, on the asymptotic results. However, in this paper we have restricted our analysis to the interior solution sufficiently late in time and far away from any boundaries so that their effects on the analysis may be neglected. An analysis including the initial and boundary conditions, and potentially boundary layers, should be the subject of future work. Other work should analyze other simplified models and numerical discretizations, and present test problems confirming the analysis, when possible. It is expected that numerical discretizations which fail to preserve the EDA’s first-order accuracy will either fail to produce accurate equilibrium-diffusion solutions or will produce them at a prohibitive computational cost. A different problem for future work is to investigate multigroup treatments of the RT equation to determine whether they preserve the EDA’s first-order accuracy.

Acknowledgements One of us (JMF) would like to thank Don Shirk and Bob Singleton for many helpful comments, as well as Scott Doebling for continued support. This work was performed under the auspices of the US Department of Energy under contract DE-AC52-06NA25396 as LA-UR-17-20878.

Appendix A The ${\cal O}(\beta^{2})$ LF RT equation

This appendix is similar to Section 93 in [2], where the mixed-frame RT equation with certain CMF functions is presented. The purpose there and here is to Taylor expand some of the CMF functions so that they depend on the LF frequency instead of the CMF frequency. However, we retain ${\cal O}(\beta^{2})$ terms and scattering terms, both of which are neglected there. We expand the frequency ratios, the cross sections, and the Planck function through ${\cal O}(\beta^{2})$. While the Planck function and the material cross sections are represented in the CMF, we drop the subscript-o for notational convenience. We also expand the radiation intensity in the integrand of equation (38) since it is a function of $\nu^{\,\prime}$, which is a function of $\nu$ by equation (40). Ignoring this expansion produces the wrong results. We refer to the resulting RT equation as the LF RT equation. The relativistically exact angle- and frequency-dependent mixed-frame RT equation is:

[TABLE]

The ratio of $\nu_{\text{o}}$ to $\nu$ is a function of $\beta$ and $\Omega_{i}$:

[TABLE]

where the Einstein summation convention is used. The ratio of LF frequencies is:

[TABLE]

The Lorentz factor, $\gamma_{\text{u}}$, expanded through ${\cal O}(\beta^{2})$ is:

[TABLE]

so the ${\cal O}(\beta^{2})$ expansion of the frequency ratios in equation (38) are:

[TABLE]

It is convenient to record here the ${\cal O}(\beta^{2})$ expansions of some identities that will be useful when Taylor expanding our functions of interest:

[TABLE]

The Taylor-expansion of a general function of the CMF frequency, with respect to the LF frequency, through ${\cal O}((\nu_{\text{o}}-\nu)^{2})\sim{\cal O}(\beta^{2})$, is:

[TABLE]

The ${\cal O}(\beta^{2})$ Taylor expansions of the total cross-section and the Planck function, and the product of the absorption cross-section with the Planck function, are:

[TABLE]

The Taylor expansion of the LF radiation intensity, in the integrand of equation (38), proceeds along the same lines. However, the CMF frequency in the previous expressions is now a LF frequency, $\nu^{\,\prime}$, and the necessary relations take a slightly different form through ${\cal O}(\beta^{2})$:

[TABLE]

In arriving at the relations above, we have used the fact that $\beta_{i}\beta_{j}\Omega_{i}(\Omega_{j}-\Omega_{j}^{\,\prime})$ is zero at ${\cal O}(\beta^{2})$ since the angular variables are then the same. The ${\cal O}(\beta^{2})$ Taylor expansion is:

[TABLE]

such that the Taylor expanded radiation intensity, through ${\cal O}(\beta^{2})$, is

[TABLE]

We now combine the results in equations (45) and (48) to use in the three terms on the RHS of equation (38). The first term is straight-forward:

[TABLE]

The second term is best analyzed by breaking it into parts. The integrand is

[TABLE]

where we have written $I^{\,\prime}_{\nu}=I_{\nu}(\Omega^{\,\prime})$ for notational convenience. The result of the angular integral is

[TABLE]

The ratio $\left(\nu/\nu_{\text{o}}\right)^{2}$ multiplying the integral is given in equation (42c), and the second term of equation (38) is now written as:

[TABLE]

The third term is:

[TABLE]

Collecting equations (49), (52) and (53), the ${\cal O}(\beta^{2})$ LF RT equation is:

[TABLE]

This result may be compared with equation 93.4 of [2], although there the ${\cal O}(\beta^{2})$ terms are neglected as are the scattering cross sections. We believe this is the first time this equation has been presented in the literature.

Appendix B The ${\cal O}(\epsilon^{2})$ analysis

In this appendix we scale and analyze the ${\cal O}(\beta^{2})$ LF RT equation (54). First, we write $\sigma_{\text{a},\nu}=\sigma_{\text{t},\nu}-\sigma_{\text{s}}$, and we reiterate that $\sigma_{\text{s}}$ is frequency-independent. The scaled LF RT equation is:

[TABLE]

The zeroth-, first- and second-order solutions of the radiation intensity can now be determined. Their frequency-integrated angular moments produce the associated radiation variables. The zeroth-order solutions are:

[TABLE]

The first-order solutions using these results are:

[TABLE]

The Rosseland-averaged cross section [2], $\sigma_{\text{t},\text{\tiny R}}$, has been used in passing from equation (57a) to (57b). The second-order radiation-intensity solution, using the results in (56) and (57), is:

[TABLE]

This solution cannot be frequency integrated since the functional form of the total cross section is unknown. Previously, integration by parts was used to move the frequency-derivatives, but that does not work now. However, the frequency-dependent radiation flux can be constructed by taking the first angular moment of (58a):

[TABLE]

This completes the calculations that are used in Section 3. Explicit expressions for ${\cal E}^{(2)}$ and ${\cal P}_{ij}^{(2)}$ cannot be produced because the functional form of the cross sections is not known. The zeroth and second angular moments of $B_{\nu}^{(2)}$ produce the EDA-like results for ${\cal E}^{(2)}$ and ${\cal P}_{ij}^{(2)}$. However, these angular moments of the zeroth-order contribution to $I_{\nu}^{(2)}$ produce additional terms beyond the EDA expressions, which we have called transport corrections. Castor [3] provides some discussion of what these terms might mean physically.