The equilibrium-diffusion limit for radiation hydrodynamics
TThe equilibrium-diffusion limit for radiation hydrodynamics
J.M. Ferguson a , J.E. Morel b , R.B. Lowrie c a XCP-Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA b Department of Nuclear Engineering, Texas A & M University, College Station, TX 77843, USA c CCS-Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
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 bederived by asymptotically analyzing the full set of RH equations in the equilibrium-diffusion limit. We derivethe EDA this way and show that it and the associated set of simplified equations are both first-order accuratewith transport corrections occurring at second order. Having established the EDA’s first-order accuracy wethen analyze the grey nonequilibrium-diffusion approximation and the grey Eddington approximation andshow that they both preserve this first-order accuracy. Further, these approximations preserve the EDA’sfirst-order accuracy when made in either the comoving-frame (CMF) or the lab-frame (LF). While analyzingthe Eddington approximation, we found that the CMF and LF radiation-source equations are equivalentwhen neglecting O ( β ) 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 notpreserve this first-order accuracy will not retain the correct equilibrium-diffusion solutions. As a practicalexample, we show that nonequilibrium-diffusion radiative-shock solutions devolve to equilibrium-diffusionsolutions when the asymptotic parameter is small. Keywords: asymptotics, equilibrium diffusion, radiation transport, radiation hydrodynamics, greynonequilibrium-diffusion approximation, grey Eddington approximation, radiative-shock solutions
1. Introduction
Radiation hydrodynamics (RH) describes howinteractions between radiation and matter affectthe thermodynamic states, and potentially, the dy-namic flow characteristics of the matter-radiationsystem. Unfortunately, the full set of RH equationsare computationally expensive and numerically dif-ficult to solve, and various model approximationshave been developed to aid their solution [1, 2, 3].The Euler equations are typically assumed to pro-vide a sufficient model of the material’s hydrody-namic response since the photon mean-free-path isgenerally much longer than the mean-free-path be-tween material interactions, so that material vis-cosity and heat conduction may be neglected. Forthe radiation it is common to make the followingassumptions, either independently or together: thesystem is absorption dominated, the system’s size islarge compared to the photon mean-free-path, theradiation is in thermal equilibrium with the ma- terial, the radiation flux is diffusive, and the ra-diation pressure is isotropic. Taken together theseassumptions are called the equilibrium diffusion ap-proximation (EDA) [2]. When the EDA applies toa physical system it can be modeled by a simplerset of equations than the full set of RH equations.This simplified set of RH equations provides a rea-sonable description of stellar structure [4], high-temperature environments [1], fusion dominated en-ergy sources [5, 6], a variety of astrophysical set-tings [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 expan-sion of the RH equations, which is described as fol-lows. First, an ansatz for the RH solution is madeand the solution is expressed as an infinite power-series expansion in a small parameter, (cid:15) . Thenthe RH equations are nondimensionalized. Next,nondimensional parameters that arise in the RHequations are chosen to scale by particular pow-
Preprint submitted to Journal of Quantitative Spectroscopy and Radiative Transfer October 15, 2018 a r X i v : . [ phy s i c s . c o m p - ph ] F e b rs of (cid:15) , which reflects the physical importance ofthose parameters, and effectively defines physicaland mathematical limits. Finally, an infinite hier-archy of equations is obtained by equating all termsmultiplied by the same power of (cid:15) . This hierarchy isused to define the equilibrium-diffusion limit (EDL)since the scalings are chosen to ensure that the ze-roth order RH solution satisfies the EDA. Thus,in the limit as (cid:15) approaches zero the RH solutiontransitions to the EDA solution. In this sense, theEDA is said to be accurate or exact to zeroth orderin the EDL. However, it is possible for an approx-imation to be exact to higher than zeroth order.This has been shown to be true for the EDA byLowrie, Morel and Hittinger [10], but they only an-alyzed the simplified RH equations at zeroth orderand they did not determine at what order transportcorrections modified the EDA. The main contribu-tion of this paper is to show that the EDA and thesimplified RH equations are both exact to first orderin the EDL, and that transport corrections occur atsecond order. It is important to be clear that theEDA is a physical approximation to the full set ofRH equations, whereas the asymptotic expansion ofthe RH equations is exact at each asymptotic level.If we summed over all asymptotic orders then theoriginal transport content would be regained withno loss of information.Asymptotic expansions have two very practicalapplications that are not necessarily obvious. Thefirst is to demonstrate that analytic approximationsto the full RH equations give proper accuracy inthe EDL. For instance, one would expect the greynonequilibrium-diffusion approximation of the RHequations to be properly accurate in the EDL. Todetermine if this is so, one performs an asymptoticexpansion of these approximate RH equations. Ifthis expansion yields the EDA to first order, thenthe solution of the approximate RH equations willproperly transition to the EDA solution as (cid:15) → O ( β ) terms. Of course, theradiation pressures are not equivalent, but theyonly differ by a symmetric-traceless term of O ( β ).The reason for analyzing these two approximationsis because the nonequilibrium-diffusion approxima-tion only applies to optically-thick systems, and isonly slightly less restrictive than the EDA. In con-tradistinction, the Eddington approximation maybe applied to highly rarefied systems where diffu-sion does not apply and which may be out of ther-mal equilibrium. Section 5 closes the main bodyof this paper with a summary and recommenda-tions for future work. In Appendix A we derivethe O ( β ) LF radiation-transport (RT) equation in-cluding scattering terms. We believe this is the firsttime that this equation has been presented in theliterature when retaining terms through O ( β ). InAppendix B we use results from Appendix A to de-rive the radiation intensity through O ( (cid:15) ). This isused to derive the zeroth- and first-order radiationvariables and the second-order radiation flux, andto show that the radiation energy and pressure con-tain transport corrections at second order.
2. The RH equations and the EDA
We present the LF RH equations in Subsection2.1, as well as the radiation-energy and radiation-momentum sources and the frequency-integratedangular moments of the radiation intensity. In Sub-section 2.2, we review the assumptions on which the EDA relies and then we present the LF EDA andits simplified RH equations.
The RH equations used here are the Euler equa-tions coupled to the radiation momentum andenergy sources, and the angle- and frequency-dependent RT equation: ∂ t ρ + ∂ i ( ρ u i ) = 0 , (1a) ∂ t ( ρ u i ) + ∂ j ( ρ u i u j + p ij ) = − S rp ,i , (1b) ∂ t E + ∂ i [ u j ( E ij + p ij )] = − S re , (1c)1 c ∂ t I ν + Ω i ∂ i I ν = Q ν , (1d)where Q ν = − (cid:16) ν o ν (cid:17) σ t ,ν o I ν + (cid:18) νν o (cid:19) σ a ,ν o B ν o ( T )+ (cid:18) νν o (cid:19) σ s π (cid:90) π (cid:16) ν o ν (cid:48) (cid:17) I ν (cid:48) (Ω (cid:48) ) d Ω (cid:48) , (1e)is the angle- and frequency-dependent LF radi-ation source. The angle-dependence resides inthe frequency ratios which are relativistically ex-act through all orders in β ≡ u/c , as well as theangle- and frequency-dependent LF radiation in-tensity, I ν = I ν (Ω), and Ω i is the LF photondirection of flight. The time and space deriva-tives, ∂ t ≡ ∂/∂t and ∂ i ≡ ∂x i , are with respectto LF time and space, where x i ∈ { x, y, z } , andwe are using the Einstein summation convention.Further, ρ is the mass density, u i is the materialvelocity, p ij is the material pressure, S rp ,i is theradiation-momentum source such that − S rp ,i is amaterial-momentum source, E = ρ u + ρ e and E ij ≡ ρ u i u j + ρ e δ ij are defined for notationalconvenience, S re is the radiation-energy source suchthat − S re is a material-energy source, c is the speedof light, and ν and ν o are the LF and CMF frequen-cies, respectively. We leave the Planck function inthe CMF, B ν o ( T ) = 2 h ν c (cid:104) e h ν o /k B T − (cid:105) − , (2)as a function of the CMF frequency, ν o , and thelocal material temperature, T , where h is Planck’sconstant and k B is Boltzmann’s constant. Finally, σ s , σ a ,ν o , and σ t ,ν o = σ a ,ν o + σ s , are the scatter-ing, absorption and total cross-sections in the CMF; σ t ,ν o and σ a ,ν o are functions of ν o while σ s is inde-pendent of frequency. We omit the subscript-o from3he Planck function and the cross sections, whichotherwise denotes a CMF variable.The LF radiation-energy and radiation-momentum sources, S re and S rp ,i , are the firsttwo frequency-integrated angular moments of Q ν ,respectively: S re ≡ (cid:90) π (cid:90) ∞ Q ν d Ω dν = ∂ t E + ∂ i F i , (3a) S rp ,i ≡ c (cid:90) π (cid:90) ∞ Ω i Q ν d Ω dν = 1 c ∂ t F i + ∂ j P ij , (3b)and E ≡ c (cid:90) π (cid:90) ∞ I ν d Ω dν , (4a) F i ≡ (cid:90) π (cid:90) ∞ Ω i I ν d Ω dν , (4b) P ij ≡ c (cid:90) π (cid:90) ∞ Ω i Ω j I ν d Ω dν , (4c)are the LF radiation energy density, radiation flux,and radiation-pressure, which are the first threefrequency-integrated angular moments of I ν , re-spectively. The angle- and frequency-integrals ofthe radiation sources (3) cannot be analytically per-formed since Q ν (1e) contains angle- and frequency-dependent variables whose functional form is notgenerally known. The EDA imposes four basic assumptions on aphysical system [2]: 1) the photon mean-free-pathis small compared to the size of the absorption-dominated system, 2) the matter-radiation systemis in thermal equilibrium, 3) the radiation flux isdiffusive, and 4) the radiation pressure is isotropic.These simplifications allow for a basic understand-ing of the underlying phenomena, the equations de-scribing them, and the form their solutions mighttake. In the LF the EDA and its simplified set ofRH equations are: E = a R T , (5a) F i = − a R c σ t , R ∂ i T + 43 u i a R T , (5b) P ij = 13 a R T δ ij , (5c) ∂ t ρ + ∂ i ( ρ u i ) = 0 , (6a) ∂ t ( ρ u i )+ ∂ j (cid:18) ρ u i u j + p ij + 13 a R T δ ij (cid:19) = 0 , (6b) ∂ t ( E + E )+ ∂ i (cid:20) u j (cid:18) E ij + p ij + 43 a R T δ ij (cid:19)(cid:21) = ∂ i (cid:18) a R c σ t , R ∂ i T (cid:19) , (6c)where a R is the radiation constant and σ t , R is theCMF Rosseland-averaged cross section [2]. Wepoint out that the time-derivative of the radiationflux in the total-momentum equation (6b) has beendropped in accordance with the diffusion approxi-mation.The strength of the EDA in solving RH problemsis that the radiation variables (5) in this simplifiedset of equations (6) are explicit functions of the hy-drodynamic variables. Thus, once an equation-of-state for the material is specified, these simplifiedRH equations represent a solvable system that areeasier to solve than the full set of RH equations (1)[30]. To be clear, the EDA neglects most of theradiation information that is contained in the RTequation (1d), and this is why we refer to any termsbeyond the EDA as transport corrections.
3. Asymptotic analysis of the RH equations
In this section, we present our main result whichis that the EDA (5) and its simplified set of RHequations (6) are first-order accurate in the EDLwith transport corrections beginning at second or-der. These results are derived in Appendix B. InSubsection 3.1, the RH variables and simplified RHequations are nondimensionalized and nondimen-sional parameters are defined. The EDL scalingswere originally presented by Lowrie, Morel and Hit-tinger [10], and we use these in Subsection 3.2 toscale the nondimensional RH equations. Finally,the scaled equations are asymptotically expandedand the main result of this paper is presented inSubsection 3.3. This analysis and its results onlyapply to the interior portion of a RH system, anddo not account for effects due to initial nor bound-ary conditions, nor boundary layers.
Each dimensional variable is decomposed into theproduct of a constant dimensional variable with4ubscript- ∞ , indicating a reference state, whichonly retains the dimension of the original variable,and a non-constant nondimensional variable with ahat, which retains the value of the original variable: x = ˆ x l ∞ , t = ˆ t l ∞ a ∞ , u = ˆ u u ∞ ,ρ = ˆ ρ ρ ∞ , p = ˆ p ρ ∞ a ∞ , e = ˆ e a ∞ ,T = ˆ T T ∞ , σ t = ˆ σ t σ t , ∞ , σ s = ˆ σ s σ s , ∞ ,I ν = ˆ I ˆ ν a R c h T ∞ k B , ν = ˆ ν k B T ∞ h , where l ∞ is the reference length of the system, a ∞ is a reference sound speed for the fluid, u ∞ is a ref-erence fluid velocity, ρ ∞ is a reference fluid massdensity, and T ∞ is a reference fluid temperature.The total and scattering reference cross sections are σ t , ∞ and σ s , ∞ , respectively. The radiation intensityand frequency, ˆ I ˆ ν and ˆ ν , are nondimensionalized insuch a way that the nondimensional radiation en-ergy density, radiation flux, and radiation pressureare consistent with their definitions (4) above: E = ˆ E a R T ∞ , (7a) F i = ˆ F i a R c T ∞ , (7b) P ij = ˆ P ij a R T ∞ . (7c)The nondimensionalized form of the Euler equa-tions coupled to the radiation energy and momen-tum sources are: ∂ t ρ + M ∂ i ( ρ u i ) = 0 , (8a) M ∂ t ( ρ u i ) + ∂ j (cid:0) M ρ u i u j + p ij (cid:1) = − P ( U ∂ t F i + ∂ j P ij ) , (8b) ∂ t (cid:18) M ρ u + ρ e (cid:19) + ∂ i (cid:20) M u j (cid:18) M ρ u i u j + ρ e δ ij + p ij (cid:19)(cid:21) = − P ( ∂ t E + C ∂ i F i ) , (8c)where the hats have been dropped for notationalconvenience. The nondimensional parameters are: M ≡ u ∞ a ∞ , P ≡ a R T ∞ ρ ∞ a ∞ , U ≡ a ∞ c , C ≡ ca ∞ , L ≡ l ∞ λ t , ∞ = l ∞ σ t , ∞ , L s ≡ σ s , ∞ σ t , ∞ , the first four of which are used in equations (8) andthe last two are from Appendix B; M is related tothe Mach number of the material flow, P is a mea-sure of the influence radiation has on the material, U is the ratio of the reference sound speed to thespeed of light, C is the inverse of U , L is a mea-sure of the system’s size compared to the radiationmean-free-path, and L s is a measure of whether thesystem is absorption or scattering dominated. The EDL scalings of the nondimensional param-eters were originally introduced by Lowrie, Moreland Hittinger [10]: M = O (1) , P = O (1) , U = O ( (cid:15) ) , C = O ( (cid:15) − ) , L = O ( (cid:15) − ) , L s = O ( (cid:15) ) . The first scaling implies that no assumption is madeabout whether the value of M should be small norlarge. As such, the EDL supports shock-wave so-lutions as shown by Lowrie and Rauenzahn [29],which is a result we will use in Subsection 4.3. Thesecond scaling similarly implies that no assump-tion is made about whether the system is radiationdominated or not. The third and fourth scalingsare consistent with nonrelativistic hydrodynamics.The fifth scaling implies that the reference lengthof the system is much larger than the photon mean-free-path, so that radiation may diffuse through thesystem. The sixth scaling implies that absorptiondominates scattering. These last two scalings areused in Appendix B. Reviewing the nondimension-alized equations (8) we see that only the radiationflux is affected by the scalings. The redimensional-ized and scaled version of equations (8) is: ∂ t ρ + ∂ i ( ρ u i ) = 0 , (11a) ∂ t ( ρ u i ) + ∂ j ( ρ u i u j + p ij + P ij ) = − (cid:15) c ∂ t F i , (11b) ∂ t ( E + E )+ ∂ i [ u j ( E ij + p ij )] = − (cid:15) ∂ i F i , (11c)where we have moved all terms without an (cid:15) to theleft-hand side. No approximations have gone intothese equations nor any of their variables yet.5 .3. Asymptotic expansion and results The Euler equations are first-order accurate us-ing scalings that agree with those above [31]. Thematerial and radiation variables in equations (11)are now expanded in powers of (cid:15) ; as an example,the expansion for the material mass density is: ρ = ∞ (cid:88) n =0 ρ ( n ) (cid:15) n . (12)Collecting equations at equal orders in (cid:15) forms aninfinite hierarchical set. The zeroth- and first-ordercontributions from equations (11) are:[ ∂ t ρ + ∂ i ( ρ u i )] (0) = 0 , (13a)[ ∂ t ( ρ u i ) + ∂ j ( ρ u i u j + p ij + P ij )] (0) = 0 , (13b) { ∂ t ( E + E ) + ∂ i [ u j ( E ij + p ij )] } (0) = − ∂ i F (1) i , (13c)[ ∂ t ρ + ∂ i ( ρ u i )] (1) = 0 , (13d)[ ∂ t ( ρ u i ) + ∂ j ( ρ u i u j + p ij + P ij )] (1) = − c ∂ t F (0) i , (13e) { ∂ t ( E + E ) + ∂ i [ u j ( E ij + p ij )] } (1) = − ∂ i F (2) i . (13f)The radiation variables E (0) , F (0) i , P (0) ij , E (1) , F (1) i , P (1) ij , and F (2) i are determined by analyzing theRT equation. This is done in Appendix B and wepresent the results here: E (0) = (cid:2) a R T (cid:3) (0) , (14a) F (0) i = 0 , (14b) P (0) ij = (cid:20) a R T δ ij (cid:21) (0) , (14c) E (1) = (cid:2) a R T (cid:3) (1) , (14d) F (1) i = (cid:20) − a R c σ t , R ∂ i T + 43 u i a R T (cid:21) (0) , (14e) P (1) ij = (cid:20) a R T δ ij (cid:21) (1) , (14f) F (2) i = (cid:20) − a R c σ t , R ∂ i T + 43 u i a R T (cid:21) (1) . (14g) Summing these variables through first order pro-duces the EDA as given in equations (5): E = a R T , (15a) F i = − a R c σ t , R ∂ i T + 43 u i a R T , (15b) P ij = 13 a R T δ ij , (15c)which are in agreement with the results presentedin [10]. The derivation in Appendix B shows thatthe radiation energy density and radiation pressurecontain transport corrections to the EDA at secondorder, and so the EDA is at most first-order ac-curate in the EDL. Summing the scaled equations(13) through first order, and using the results justpresented (14), produces the EDA’s simplified RHequations (6): ∂ t ρ + ∂ i ( ρ u i ) = 0 , (16a) ∂ t ( ρ u i )+ ∂ j (cid:18) ρ u i u j + p ij + 13 a R T δ ij (cid:19) = 0 , (16b) ∂ t (cid:0) E + a R T (cid:1) + ∂ i (cid:20) u j (cid:18) E ij + p ij + 43 a R T δ ij (cid:19)(cid:21) = ∂ i (cid:18) a R c σ t , R ∂ i T (cid:19) . (16c)Therefore, the EDA and its simplified RH equationsare first-order accurate with transport correctionsoccurring at second order. This is the main result ofthis paper. Both simplified models and numericaldiscretizations of the RH equations should preservethis asymptotic behavior. Models that have differ-ent asymptotic behavior may not attain the correctsolution in or near the EDL, while improper dis-cretizations may require an unreasonable computa-tional cost. These results only apply to the interiorsolution of an RH problem far from any boundaries,and sufficiently late in time, so that the initial andboundary conditions may be neglected.
4. Analysis of two grey RH approximations
We have established the first-order accuracy ofthe EDA and its simplified RH equations. As a re-minder, the assumptions comprising the EDA arelisted at the beginning of Subsection 2.2. In thissection, we analyze two RH approximations that6re commonly used to solve a wide span of RHproblems and we show that they both preservethe EDA’s first-order accuracy. The two approx-imations are the frequency-independent (“grey”)nonequilibrium-diffusion approximation and thegrey Eddington approximation. These approxima-tions are applied to the radiation energy and mo-mentum sources which are presented in Subsection4.1. The grey nonequilibrium-diffusion approxima-tion is analyzed in Subsection 4.2. Compared tothe EDA it relaxes the requirement that the radi-ation be in thermal equilibrium with the material,hence the moniker “nonequilibrium”. As a practicalexample, in Subsection 4.3, we present a radiative-shock solution from this approximation and showthat the EDA solution is obtained when (cid:15) is small.In contrast to the nonequilibrium-diffusion approxi-mation, the grey Eddington approximation only re-tains the EDA assumption that the radiation pres-sure is isotropic. As such, it can be used to de-scribe environments where diffusion does not apply,or that may be out of equilibrium, or which may behighly rarefied. This approximation is analyzed inSubsection 4.4.
The two approximations analyzed in this sectionare typically made in the CMF by modifying theradiation energy and momentum sources. Thesesources are equations (6.47) and (6.48) in [3]: ∂ t E o + 1 c ∂ t ( β i F o ,i ) + a i c F o ,i + ∂ i F o ,i + ∂ i ( u i E o ) + P o ,ij ∂ j u i = σ a c (cid:0) a R T − E o (cid:1) , (17a)1 c ∂ t F o ,i + a i c E o + 1 c ∂ t ( β j P o ,ij )+ ∂ j P o ,ij + 1 c F o ,j ∂ j β i + 1 c ∂ j ( β j F o ,i )= − σ t c F o ,i , (17b)where we have dropped the O ( β ) terms comingfrom the Lagrangian derivatives. In order to com-pare the CMF approximation to our results weLorentz transform the approximation to the LFand analyze it there. A discussion of how to per-form Lorentz transformations is outside the scopeof this paper and can be found in [2]. Althoughthe Lorentz transformation changes the coordinateframe being considered it does not alter the physi-cal content of the approximation. It should also be clear that it does not matter if the CMF radiationsources are scaled first or transformed first, so longas the scalings are also applied to the transforma-tion. Finally, in order to facilitate comparison ofthe radiation sources from the transformed approx-imation with the LF sources we record here the LFradiation energy and momentum sources: ∂ t E + ∂ i F i = σ a c (cid:0) a R T − E (cid:1) + β i ( σ a − σ s ) F i , (18a)1 c ∂ t F i + ∂ j P ij = − σ t c F i + β j (cid:0) σ t P ij + σ s E δ ij + σ a a R T δ ij (cid:1) , (18b)which are equations (2.29) and (2.30) in [32]. Thesimplest way to produce the scaled radiation-sourceequations is to make the following replacements, ∂ t → (cid:15) ∂ t , ∂ i → (cid:15) ∂ i , β i → (cid:15) β i , σ s → (cid:15) σ s , and then to divide the radiation-energy source by (cid:15) and the radiation-momentum source by (cid:15) . Thereason for dividing the radiation sources by (cid:15) and (cid:15) is to keep these sources at the corresponding orderof their associated hydrodynamic sources. The nonequilibrium-diffusion approximationmodifies the radiation-source equations by impos-ing the Eddington approximation, P o ,ij = E o δ ij or P o ,ii = E , dropping the time-derivative of theradiation flux and all acceleration terms, and drop-ping all other terms in the radiation-momentumsource until Fick’s First Law of Diffusion is ob-tained. Applied to the CMF radiation sources(17) it produces the following simplified sourceequations: ∂ t E o + ∂ i F o ,i + ∂ i ( u i E o ) + 13 E o ∂ i u i = σ a c (cid:0) a R T − E o (cid:1) , (19a)13 ∂ i E o = − σ t c F o ,i . (19b)Solving the radiation-momentum source for the ra-diation flux and plugging this into the radiation-energy source produces: F o,i = − c σ t ∂ i E o , (20a) ∂ t E o − ∂ i (cid:18) c σ t ∂ i E o (cid:19) + ∂ i ( u i E o ) + 13 E o ∂ i u i σ a c (cid:0) a R T − E o (cid:1) , (20b)which are equations (97.69) and (97.70) in [2].Lorentz transformation of the sources (19) alongwith the CMF Eddington approximation produces: ∂ t E + ∂ i F i − c ∂ t ( β i F i )= σ a c (cid:0) a R T − E (cid:1) + β i ( σ a − σ s ) F i , (21a)13 ∂ i E + β i c ∂ t E + β i c ∂ j F j = − σ t c F i + β i (cid:18) σ t E + σ s E + σ a a R T (cid:19) , (21b) P ij = 13 E δ ij + 1 c (cid:18) β j F i + β i F j − β k F k δ ij (cid:19) . (21c)The O ( β ) term on the radiation pressure (21c) issymmetric and traceless, so the Eddington approx-imation is retained after the Lorentz transforma-tion, P ii = E , since the trace of a traceless ob-ject is zero. We quickly compare the transformedsources, (21a) and (21b), with the LF sources (18).If the nonequilibrium-diffusion approximation wereapplied to the LF radiation sources (18) then theresult would differ from equations (21a) and (21b)by terms of O ( β ) and β i F i , and the radiation pres-sure would neglect the symmetric traceless portioncontained in equation (21c). However, we will findthat these differences between the CMF and LFsources, as well as the radiation pressures, are neg-ligible in the EDL. Scaling the radiation sources,(21a) and (21b), gives: ∂ t E + 1 (cid:15) ∂ i F i − (cid:15) c ∂ t ( β i F i )= 1 (cid:15) ( σ t − (cid:15) σ s ) c (cid:0) a R T − E (cid:1) + 1 (cid:15) β i ( σ t − (cid:15) σ s ) F i , (22a)13 ∂ i E + (cid:15) β i c ∂ t E + (cid:15) β i c ∂ j F j = − (cid:15) σ t c F i + β i (cid:16) σ t E + (cid:15) σ s E + ( σ t − (cid:15) σ s ) a R T (cid:17) . (22b)The radiation-energy source (22a) provides the so-lution for the scaled radiation energy density, E = a R T − (cid:15) (cid:20) σ s σ t (cid:0) a R T − E (cid:1) − β i c F i + 1 σ t c ∂ i F i (cid:21) , (23a)where terms of O ( (cid:15) ) have been discarded. Theradiation-momentum source (22b) provides the so-lution for the scaled radiation flux, F i = (cid:15) (cid:20) − c σ t ∂ i E + u i (cid:18) E + a R T (cid:19)(cid:21) + (cid:15) (cid:20) − σ t β i ∂ j F j + u i σ s σ t (cid:0) E − a R T (cid:1)(cid:21) , (23b)where terms of O ( (cid:15) ) have been discarded. Scalingthe radiation pressure (21c) gives: P ij = 13 E δ ij + 1 c (cid:15) (cid:18) β j F i + β i F j − β k F k δ ij (cid:19) . (23c)From these expressions (23) the zeroth- and first-order solutions are obtained along with the second-order radiation flux: E (0) = (cid:2) a R T (cid:3) (0) , (24a) F (0) i = 0 , (24b) P (0) ij = (cid:20) a R T δ ij (cid:21) (0) , (24c) E (1) = (cid:2) a R T (cid:3) (1) , (24d) F (1) i = (cid:15) (cid:20) − a R c σ t ∂ i T + 43 u i a R T (cid:21) (0) , (24e) P (1) ij = (cid:20) a R T δ ij (cid:21) (1) , (24f) F (2) i = (cid:15) (cid:20) − a R c σ t ∂ i T + 43 u i a R T (cid:21) (1) . (24g)The zeroth- and first-order contributions from thescaled radiation sources (22) are: S (0)re = ∂ t E (0) + ∂ i F (1) i , (25a) S (0)rp = 13 ∂ i E (0) , (25b) S (1)re = ∂ t E (1) + ∂ i F (2) i − ∂ t (cid:20) c β i F i (cid:21) (0) , (25c) S (1)rp = 13 ∂ i E (1) + (cid:20) β i c ∂ j F j (cid:21) (0) . (25d)8umming the zeroth- and first-order results for theradiation variables (24) reproduces the EDA (5): E = a R T , (26a) F i = − a R c σ t ∂ i T + 43 u i a R T , (26b) P ij = 13 a R T δ ij . (26c)As a reminder, the radiation sources are coupled tothe Euler equations; see equations (1). Summingthe zeroth- and first-order contributions from theradiation sources (25), coupled to the Euler equa-tions (1), and using the results in equations (24),reproduces the EDA’s simplified RH equations (6): ∂ t ρ + ∂ i ( ρu i ) = 0 , (27a) ∂ t ( ρ u i )+ ∂ j (cid:18) ρ u i u j + p ij + 13 a R T δ ij (cid:19) = 0 , (27b) ∂ t (cid:0) E + a R T (cid:1) + ∂ i (cid:20) u j (cid:18) E ij + p ij + 43 a R T δ ij (cid:19)(cid:21) = ∂ i (cid:18) a R c σ t ∂ i T (cid:19) . (27c)Thus, the nonequilibrium-diffusion approximationpreserves the EDA’s first-order accuracy, and thisresult holds if the approximation is made in theCMF or the LF. In this subsection we show that a particular so-lution of the nonequilibrium-diffusion approxima-tion transitions to the equilibrium-diffusion solu-tion when (cid:15) is small. Specifically, we briefly an-alyze the nonequilibrium-diffusion radiative-shocksolution method developed by Lowrie and Edwards[28] and we show that when (cid:15) = 0 .
001 that methodproduces the equilibrium-diffusion solution devel-oped by Lowrie and Rauenzahn [29].The 1D nondimensional equations solved byLowrie and Edwards, scaled in the EDL, are: ∂ x ( ρ u ) = 0 , (28a) ∂ x (cid:18) ρ u + p + 13 P E (cid:19) = 0 , (28b) ∂ x (cid:20) u (cid:18) ρ u + ρ e + p (cid:19)(cid:21) = − P (cid:15) ∂ x F , (28c) -0.01 0 0.00513.54.5 T , (cid:15) = 1 θ , (cid:15) = 1 T , (cid:15) = 0 . T , EDA Figure 1: Comparison of the M = 3 nonequilibrium-diffusion radiative-shock solutions for (cid:15) = 1 and (cid:15) = 0 . (cid:15) , as described in Subsection 4.3.The solutions associated with (cid:15) = 1 are the nonequilibrium-diffusion radiative-shock solutions, while the solution asso-ciated with (cid:15) = 0 .
001 is the equilibrium-diffusion solution.The equilibrium-diffusion solution produced by the methoddescribed by Lowrie and Rauenzahn is the dash-dotted linelabeled as the EDA solution. The EDA solution and thesolution using (cid:15) = 0 .
001 cannot be distinguished. M γ ( γ − ∂ x T + p ∂ x u = P σ t (cid:20) (cid:15) (cid:0) T − E (cid:1) + 2 β F (cid:21) . (28d)The fourth equation (28d) is the radiation inter-nal energy source, S rie = S re − β S rp . For amonatomic ideal-gas the adiabatic index is γ = 5 / M represents the initial Mach number of theunshocked gas, which we set to be M = 3 in thisexample. The nondimensional radiation flux for apurely absorbing system is F = − σ t ∂ x E + 13 β σ t (cid:0) E + T (cid:1) . (29)When (cid:15) = 1 we obtain the nonequilibrium-diffusionsolution and when (cid:15) = 0 .
001 we obtain theequilibrium-diffusion solution from the same solu-tion method; see Figure 1. For comparison withinthe figure, we also include the equilibrium-diffusionsolution produced by the method described byLowrie and Rauenzahn, which is the dash-dottedline labeled as the EDA solution. This EDA solu-tion and the nonequilibrium-diffusion solution us-ing (cid:15) = 0 .
001 cannot be distinguished. Thus, theseshock problems can serve as a test problem to en-sure that codes are asymptotic preserving. How-ever, other test problems are also needed that are9elevant to different portions of the RH problemspace.
In this subsection we analyze the grey Edding-ton approximation applied to the CMF radiationsources (17). The application of this approxima-tion produces the following CMF radiation sources: ∂ t E o + 1 c ∂ t ( β i F o ,i ) + a i c F o ,i + ∂ i F o ,i + ∂ i ( u i E o ) + 13 E o ∂ i u i = σ a c (cid:0) a R T − E o (cid:1) , (30a)1 c ∂ t F o ,i + a i c E o + 13 c ∂ t ( β i E o )+ 13 ∂ i E o + 1 c F o ,j ∂ j β i + 1 c ∂ j ( β j F o ,i )= − σ t c F o ,i . (30b)Lorentz transformation of these sources along withthe Eddington approximation to the LF produces: ∂ t E + ∂ i F i = σ a c (cid:0) a R T − E (cid:1) + β i ( σ a − σ s ) F i , (31a)1 c ∂ t F i + ∂ j P ij = − σ t c F i + β j (cid:0) σ t P ij + σ s E δ ij + σ a a R T δ ij (cid:1) , (31b) P ij = 13 E δ ij + 1 c (cid:18) β j F i + β i F j − β k F k δ ij (cid:19) . (31c)Again, the trace of the radiation pressure returnsthe Eddington approximation, P ii = E , since the O ( β ) term is traceless, as well as symmetric. Thesesource equations are identical to the LF radiationsources (18). Therefore, the analysis below appliesequally well when the Eddington approximation isapplied to the LF sources. Scaling these radiationsources, (31a) and (31b), gives: ∂ t E + 1 (cid:15) ∂ i F i = 1 (cid:15) ( σ t − (cid:15) σ s ) c (cid:0) a R T − E (cid:1) + 1 (cid:15) β i ( σ t − (cid:15) σ s ) F i , (32a) (cid:15) c ∂ t F i + ∂ j P ij = − (cid:15) σ t c F i + β j (cid:104) σ t P ij + (cid:15) σ s E δ ij + ( σ t − (cid:15) σ s ) a R T δ ij (cid:105) . (32b)The radiation-energy source (32a) provides the so-lution for the scaled radiation energy density, E = a R T − (cid:15) (cid:34) σ s σ t (cid:0) a R T − E (cid:1) − β i c F i + 1 σ t c ∂ i F i (cid:35) , (33a)where terms of O ( (cid:15) ) have been discarded. Theradiation-momentum source (32b) provides the so-lution for the scaled radiation flux, F i = (cid:15) (cid:20) − cσ t ∂ j P ij + u j (cid:0) P ij + a R T δ ij (cid:1)(cid:21) + (cid:15) (cid:20) − c σ t ∂ t F i + u i σ s σ t (cid:0) E − a R T (cid:1)(cid:21) . (33b)Scaling the radiation pressure (31c) gives: P ij = 13 E δ ij + (cid:15) c (cid:18) β j F i + β i F j − β k F k δ ij (cid:19) . (33c)From these expressions (33) the zeroth- and first-order solutions are obtained along with the second-order radiation flux: E (0) = (cid:2) a R T (cid:3) (0) , (34a) F (0) i = 0 , (34b) P (0) ij = (cid:20) a R T δ ij (cid:21) (0) , (34c) E (1) = (cid:2) a R T (cid:3) (1) , (34d) F (1) i = (cid:15) (cid:20) − a R c σ t ∂ i T + 43 u i a R T (cid:21) (0) , (34e) P (1) ij = (cid:20) a R T δ ij (cid:21) (1) , (34f) F (2) i = (cid:15) (cid:20) − a R c σ t ∂ i T + 43 u i a R T (cid:21) (1) . (34g)The zeroth- and first-order contributions from thescaled radiation sources (32) are: S (0)re = ∂ t E (0) + ∂ i F (1) i , (35a) S (0)rp = ∂ j P (0) ij , (35b) S (1)re = ∂ t E (1) + ∂ i F (2) i , (35c)10 (1)rp = 1 c ∂ t F (0) i + ∂ j P (1) ij . (35d)Summing the zeroth- and first-order results for theradiation variables (34) reproduces the EDA (5): E = a R T , (36a) F i = − a R c σ t ∂ i T + 43 u i a R T , (36b) P ij = 13 a R T δ ij . (36c)Summing the zeroth- and first-order contributionsfrom the radiation sources (35), coupled to the Eu-ler equations (1), and using the results (34), repro-duces the EDA’s simplified RH equations (6): ∂ t ρ + ∂ i ( ρ u i ) = 0 , (37a) ∂ t ( ρ u i )+ ∂ j (cid:18) ρ u i u j + p ij + 13 a R T δ ij (cid:19) = 0 , (37b) ∂ t (cid:0) E + a R T (cid:1) + ∂ i (cid:20) u j (cid:18) E ij + p ij + 43 a R T δ ij (cid:19)(cid:21) = ∂ i (cid:18) a R c σ t , R ∂ i T (cid:19) . (37c)Thus, the Eddington approximation preserves theEDA’s first-order accuracy. Further, it does notmatter whether the Eddington approximation is ap-plied to the CMF or to the LF radiation sources,the EDA’s first-order accuracy is preserved in bothsituations.
5. Summary
In this work we have derived the EDA fromthe RH equations via an asymptotic analysis.Our derivation showed that the EDA and itssimplified RH equations are first-order accurateand that transport corrections begin at secondorder. Since the EDA is a physical limit of the fullset of RH equations it is expected that simplifiedmodels of the RH equations should preserve theEDA’s first-order accuracy. We analyzed the greynonequilibrium-diffusion approximation and thegrey Eddington approximation and we showedthat they both preserved the EDA’s first-orderaccuracy. These approximations can be made inthe CMF or the LF, and we have shown that theEDA’s first-order accuracy is preserved in both cases. We also presented a test problem in whichan equilibrium-diffusion solution was capturedfrom a nonequilibrium-diffusion solver when (cid:15) was small. Other test problems that apply todifferent RH regimes are needed. Our results are inagreement with previous asymptotic analyses forneutron transport [11, 12] and radiative transfer[14]. Other analyses [17, 18, 13, 15] for neutrontransport and radiative transfer have discussedthe effects of initial and boundary conditions, aswell as boundary layers, on the asymptotic results.However, in this paper we have restricted ouranalysis to the interior solution sufficiently late intime and far away from any boundaries so thattheir effects on the analysis may be neglected.An analysis including the initial and boundaryconditions, and potentially boundary layers, shouldbe the subject of future work. Other work shouldanalyze other simplified models and numerical dis-cretizations, and present test problems confirmingthe analysis, when possible. It is expected thatnumerical discretizations which fail to preservethe EDA’s first-order accuracy will either fail toproduce accurate equilibrium-diffusion solutions orwill produce them at a prohibitive computationalcost. A different problem for future work isto investigate multigroup treatments of the RTequation to determine whether they preserve theEDA’s first-order accuracy.
Acknowledgements
One of us (JMF) would liketo thank Don Shirk and Bob Singleton for manyhelpful comments, as well as Scott Doebling for con-tinued support. This work was performed underthe auspices of the US Department of Energy un-der contract DE-AC52-06NA25396 as LA-UR-17-20878.
Appendix A. The O ( β ) LF RT equation This appendix is similar to Section 93 in [2],where the mixed-frame RT equation with certainCMF functions is presented. The purpose thereand here is to Taylor expand some of the CMFfunctions so that they depend on the LF frequencyinstead of the CMF frequency. However, we retain O ( β ) terms and scattering terms, both of whichare neglected there. We expand the frequency ra-tios, the cross sections, and the Planck functionthrough O ( β ). While the Planck function and thematerial cross sections are represented in the CMF,we drop the subscript-o for notational convenience.11e also expand the radiation intensity in the in-tegrand of equation (A.1) since it is a function of ν (cid:48) , which is a function of ν by equation (A.3). Ig-noring this expansion produces the wrong results.We refer to the resulting RT equation as the LFRT equation. The relativistically exact angle- andfrequency-dependent mixed-frame RT equation is:1 c ∂ t I ν + Ω i ∂ i I ν = − ν o ν σ t ,ν o I ν + (cid:18) νν o (cid:19) σ s π (cid:90) π ν o ν (cid:48) I ν (cid:48) (Ω (cid:48) ) d Ω (cid:48) + (cid:18) νν o (cid:19) σ a ,ν o B ν o . (A.1)The ratio of ν o to ν is a function of β and Ω i : ν o ν = γ u (1 − β i Ω i ) , (A.2)where the Einstein summation convention is used.The ratio of LF frequencies is: ν (cid:48) ν = 1 − β i Ω i − β i Ω (cid:48) i . (A.3)The Lorentz factor, γ u , expanded through O ( β ) is: γ u = 1 + 12 β , (A.4)so the O ( β ) expansion of the frequency ratios inequation (A.1) are: ν o ν = 1 − β i Ω i + 12 β , (A.5a) ν o ν (cid:48) = 1 − β i Ω (cid:48) i + 12 β , (A.5b) (cid:18) νν o (cid:19) = 1 + 2 β i Ω i + 3 ( β i Ω i ) − β . (A.5c)It is convenient to record here the O ( β ) expansionsof some identities that will be useful when Taylorexpanding our functions of interest: ν o − ν = ν (cid:18) − β i Ω i + 12 β (cid:19) , (A.6a)( ν o − ν ) = ν ( β i Ω i ) , (A.6b) ν = ν o (cid:20) β i Ω i + ( β i Ω i ) − β (cid:21) , (A.6c) ∂ ν o ν = 1 + β i Ω i + ( β i Ω i ) − β = νν o , (A.6d) ν ( ν o − ν ) ν o = ν (cid:20) − β i Ω i − ( β i Ω i ) + 12 β (cid:21) , (A.6e) ( ν o − ν ) f ( ν o ) = ( ν o − ν ) f ( ν ) . (A.6f)The Taylor-expansion of a general function of theCMF frequency, with respect to the LF frequency,through O (( ν o − ν ) ) ∼ O ( β ), is: f ( ν o ) = f + ( ν o − ν ) ∂ ν o f + 12 ( ν o − ν ) ∂ ν o f = f − β i Ω i ν∂ ν f + 12 β i β j (( δ ij − i Ω j ) ν∂ ν f + Ω i Ω j ν ∂ ν f (cid:1) . (A.7)The O ( β ) Taylor expansions of the total cross-section and the Planck function, and the product ofthe absorption cross-section with the Planck func-tion, are: σ t ,ν o = σ t ,ν − β i Ω i ν∂ ν σ t ,ν + 12 β i β j (( δ ij − i Ω j ) ν∂ ν σ t ,ν + Ω i Ω j ν ∂ ν σ t ,ν (cid:1) , (A.8a) B ν o = B ν − β i Ω i ν∂ ν B ν + 12 β i β j (( δ ij − i Ω j ) ν∂ ν B ν + Ω i Ω j ν ∂ ν B ν (cid:1) , (A.8b) σ a ,ν o B ν o = σ a ,ν B ν − β i Ω i ν∂ ν ( σ a ,ν B ν )+ 12 β i β j (( δ ij − i Ω j ) ν∂ ν ( σ a ,ν B ν )+ Ω i Ω j ν ∂ ν ( σ a ,ν B ν ) (cid:1) . (A.8c)The Taylor expansion of the LF radiation intensity,in the integrand of equation (A.1), proceeds alongthe same lines. However, the CMF frequency in theprevious expressions is now a LF frequency, ν (cid:48) , andthe necessary relations take a slightly different formthrough O ( β ): ν = ν (cid:48) (1 + β i (Ω (cid:48) i − Ω i )) , (A.9a)( ν (cid:48) − ν ) = νβ i (Ω (cid:48) i − Ω i ) , (A.9b)( ν (cid:48) − ν ) = ν β i β j (Ω (cid:48) i − Ω i ) (cid:0) Ω (cid:48) j − Ω j (cid:1) , (A.9c) ∂ ν (cid:48) ν = 1 + β i (Ω (cid:48) i − Ω i ) = νν (cid:48) , (A.9d)( ν (cid:48) − ν ) νν (cid:48) = νβ i (Ω (cid:48) i − Ω i ) . (A.9e)In arriving at the relations above, we have used thefact that β i β j Ω i (Ω j − Ω (cid:48) j ) is zero at O ( β ) since theangular variables are then the same. The O ( β )Taylor expansion is:12 ( ν (cid:48) ) = f + ( ν (cid:48) − ν ) (cid:16) νν (cid:48) (cid:17) ∂ ν f + 12 ( ν (cid:48) − ν ) (cid:18) ν ∂ ν f + ∂ ν f (cid:19) = f + β i (Ω (cid:48) i − Ω i ) ν∂ ν f , such that the Taylor expanded radiation intensity,through O ( β ), is I ν (Ω (cid:48) ) = I ν + β i (Ω (cid:48) i − Ω i ) ν∂ ν I ν . (A.10)We now combine the results in equations (A.8) and(A.10) to use in the three terms on the RHS ofequation (A.1). The first term is straight-forward: − ν o ν σ t ,ν o I ν = − σ t ,ν I ν + β i Ω i ( σ t ,ν I ν + I ν ν ∂ ν σ t ,ν ) − β i β j ( σ t ,ν I ν δ ij + I ν ν ∂ ν σ t ,ν δ ij + Ω i Ω j I ν ν ∂ ν σ t ,ν (cid:1) . (A.11)The second term is best analyzed by breaking itinto parts. The integrand is ν o ν (cid:48) I ν (cid:48) (Ω (cid:48) ) = I (cid:48) ν + β i ( − Ω (cid:48) i I (cid:48) ν + (Ω (cid:48) i − Ω i ) ν ∂ ν I (cid:48) ν )+ 12 β I (cid:48) ν , (A.12)where we have written I (cid:48) ν = I ν (Ω (cid:48) ) for notationalconvenience. The result of the angular integral is (cid:90) π ν o ν (cid:48) I ν (cid:48) (Ω (cid:48) ) d Ω (cid:48) = φ ν + β i ( − F ν,i + ν ∂ ν F ν,i − Ω i ν ∂ ν φ ν )+ 12 β φ ν . (A.13)The ratio ( ν/ν o ) multiplying the integral is givenin equation (A.5c), and the second term of equation(A.1) is now written as: (cid:18) νν o (cid:19) σ s π (cid:90) π ν o ν (cid:48) I ν (cid:48) (Ω (cid:48) ) d Ω (cid:48) = σ s π φ ν + σ s π β i (2 Ω i φ ν − Ω i ν ∂ ν φ ν − F ν,i + ν ∂ ν F ν,i )+ σ s π β i β j (cid:20)(cid:18) i Ω j − δ ij (cid:19) φ ν − i Ω j ν ∂ ν φ ν (cid:21) − σ s π β i β j Ω i ( F ν,j − ν ∂ ν F ν,j ) . (A.14)The third term is: (cid:18) νν o (cid:19) σ a ,ν o B ν o ( T ) = σ a ,ν B ν + β i Ω i [2 σ a ,ν B ν − ν ∂ ν ( σ a ,ν B ν )]+ β i β j (cid:34) (3 Ω i Ω j − δ ij ) σ a ,ν B ν + 12 Ω i Ω j ν ∂ ν ( σ a ,ν B ν )+ (cid:18) δ ij − i Ω j (cid:19) ν ∂ ν ( σ a ,ν B ν ) (cid:35) . (A.15)Collecting equations (A.11), (A.14) and (A.15), the O ( β ) LF RT equation is:1 c ∂ t I ν + Ω i ∂ i I ν = σ s π φ ν + σ a ,ν B ν − σ t ,ν I ν + β i Ω i (cid:104) σ t ,ν I ν + I ν ν∂ ν σ t ,ν + 2 σ s π φ ν + 2 σ a ,ν B ν − σ s π ν∂ ν φ ν − ν∂ ν ( σ a ,ν B ν ) (cid:105) − σ s π β i ( F ν,i − ν∂ ν F ν,i ) − β i β j (cid:0) σ t ,ν δ ij I ν + δ ij I ν ν∂ ν σ t ,ν + Ω i Ω j I ν ν ∂ ν σ t ,ν (cid:1) + σ s π β i β j (cid:34)(cid:18) i Ω j − δ ij (cid:19) φ ν − i Ω j ν∂ ν φ ν − i F ν,j + 2Ω i ν∂ ν F ν,j (cid:35) + β i β j (cid:34) (3Ω i Ω j − δ ij ) σ a ,ν B ν + 12 Ω i Ω j ν ∂ ν ( σ a ,ν B ν )+ (cid:18) δ ij − i Ω j (cid:19) ν∂ ν ( σ a ,ν B ν ) (cid:35) . (A.16)This result may be compared with equation 93.4 of[2], although there the O ( β ) terms are neglectedas are the scattering cross sections. We believe thisis the first time this equation has been presented inthe literature. Appendix B. The O ( (cid:15) ) analysis In this appendix we scale and analyze the O ( β )LF RT equation (A.16). First, we write σ a ,ν = σ t ,ν − σ s , and we reiterate that σ s is frequency-independent. The scaled LF RT equation is: (cid:15) c ∂ t I ν + (cid:15) Ω i ∂ i I ν = σ t ,ν ( B ν − I ν )+ (cid:15) (cid:26) σ s π ( φ ν − π B ν ) + β i Ω i (cid:104) σ t ,ν I ν + I ν ν ∂ ν σ t ,ν +2 σ t ,ν B ν − ν ∂ ν ( σ t ,ν B ν ) (cid:105)(cid:27) + (cid:15) (cid:40) σ s π β i (cid:104) − F ν,i + ν ∂ ν F ν,i
13 Ω i (2 φ ν − π B ν − ν ∂ ν φ ν + 4 π ν ∂ ν B ν ) (cid:105) + β i β j (cid:34) (3 Ω i Ω j − δ ij ) σ t ,ν B ν + 12 Ω i Ω j ν ∂ ν ( σ t ,ν B ν )+ (cid:18) δ ij − i Ω j (cid:19) ν ∂ ν ( σ t ,ν B ν ) − (cid:16) σ t ,ν I ν δ ij + δ ij I ν ν ∂ ν σ t ,ν + Ω i Ω j I ν ν ∂ ν σ t ,ν (cid:17)(cid:35)(cid:41) . (B.1a)This can be rearranged to produce the O ( (cid:15) ) solu-tion for the radiation intensity: I ν = B ν + (cid:15) (cid:26) − σ t ,ν Ω i ∂ i I ν + β i Ω i (cid:16) I ν + 2 B ν + νσ t ,ν (cid:2) I ν ∂ ν σ t ,ν − ∂ ν ( σ t ,ν B ν ) (cid:3)(cid:17) + σ s σ t ,ν (cid:18) φ ν π − B ν (cid:19)(cid:27) + (cid:15) (cid:26) − c σ t ,ν ∂ t I ν − σ s σ t ,ν β i (cid:20) π ( F ν,i − ν ∂ ν F ν,i ) − Ω i (cid:18) φ ν π − ν π ∂ ν φ ν − B ν + ν ∂ ν B ν (cid:19)(cid:21) + β i β j (cid:20) (3 Ω i Ω j − δ ij ) B ν − I ν (cid:18) νσ t ,ν ∂ ν σ t ,ν (cid:19) δ ij + 12 σ t ,ν ( δ ij − i Ω j ) ν ∂ ν ( σ t ,ν B ν )+ Ω i Ω j ν σ t ,ν (cid:2) ∂ ν ( σ t ,ν B ν ) − I ν ∂ ν σ t ,ν (cid:3)(cid:21)(cid:27) . (B.1b)The zeroth-, first- and second-order solutions ofthe radiation intensity can now be determined.Their frequency-integrated angular moments pro-duce the associated radiation variables. The zeroth-order solutions are: I (0) ν = B (0) ν , (B.2a) I (0) = (cid:90) ∞ I (0) ν dν = (cid:104) a R c π T (cid:105) (0) , (B.2b) E (0) = 1 c (cid:90) π I (0) d Ω = (cid:2) a R T (cid:3) (0) , (B.2c) F (0) i = (cid:90) π Ω i I (0) d Ω = 0 , (B.2d) P (0) ij = 1 c (cid:90) π Ω i Ω j I (0) d Ω = (cid:20) E δ ij (cid:21) (0) . (B.2e)The first-order solutions using these results are: I (1) ν = B (1) ν + Ω i (cid:104) − σ t ∂ i B ν + β i (3 B ν − ν ∂ ν B ν ) (cid:105) (0) , (B.3a) I (1) = (cid:90) ∞ I (1) ν dν = (cid:104) a R c π T (cid:105) (1) + Ω i π (cid:20) − a R cσ t , R ∂ i T + 4 u i a R T (cid:21) (0) , (B.3b) E (1) = 1 c (cid:90) π I (1) d Ω = (cid:2) a R T (cid:3) (1) , (B.3c) F (1) i = (cid:90) π Ω i I (1) d Ω= (cid:20) − a R c σ t , R ∂ i T + 43 a R u i T (cid:21) (0) , (B.3d) P (1) ij = 1 c (cid:90) π Ω i Ω j I (1) d Ω = (cid:20) E δ ij (cid:21) (1) . (B.3e)The Rosseland-averaged cross section [2], σ t , R , hasbeen used in passing from equation (B.3a) to(B.3b). The second-order radiation-intensity solu-tion, using the results in (B.2) and (B.3), is: I (2) ν = B (2) ν + Ω i (cid:26) − σ t ∂ i B ν + β i (3 B ν − ν ∂ ν B ν ) (cid:27) (1) + (cid:26) − c σ t ∂ t B ν − β (3 B ν − ν ∂ ν B ν )+ Ω i Ω j (cid:20) β j (cid:18) − σ t (cid:18) σ t ν ∂ ν σ t (cid:19) ∂ i B ν + β i (cid:18) B ν − ν ∂ ν B ν + 12 ν ∂ ν B ν (cid:19)(cid:19)(cid:21)(cid:27) (0) . (B.4a)This solution cannot be frequency integrated sincethe functional form of the total cross section is un-known. Previously, integration by parts was used tomove the frequency-derivatives, but that does notwork now. However, the frequency-dependent ra-diation flux can be constructed by taking the firstangular moment of (B.4a): F ν,i = (cid:90) π Ω i I (2) ν d Ω = 4 π (cid:20) − σ t ∂ i B ν + β i (3 B ν − ν ∂ ν B ν ) (cid:21) (1) . (B.5a)This expression can be frequency-integrated to pro-duce the second-order radiation flux: F (2) i = (cid:90) ∞ F (2) ν,i dν (cid:20) − a R c σ t , R ∂ i T + 43 u i a R T (cid:21) (1) . (B.5b)This completes the calculations that are used inSection 3. Explicit expressions for E (2) and P (2) ij cannot be produced because the functional form ofthe cross sections is not known. The zeroth andsecond angular moments of B (2) ν produce the EDA-like results for E (2) and P (2) ij . However, these angu-lar moments of the zeroth-order contribution to I (2) ν produce additional terms beyond the EDA expres-sions, which we have called transport corrections.Castor [3] provides some discussion of what theseterms might mean physically. References [1] Y. B. Zel‘dovich, Y. P. Raizer, Physics of ShockWaves and High-Temperature Hydrodynamic Phenom-ena, Dover Books on Physics Series, Dover Publications,Mineola, N.Y., 2002.[2] D. Mihalas, B. W. Mihalas, Foundations of RadiationHydrodynamics, Dover Books on Physics Series, Dover,Mineola, N.Y., 1999.[3] J. I. Castor, Radiation Hydrodynamics, CambridgeUniversity Press, New York, N.Y., 2007.[4] S. Chandrasekhar, Radiative Transfer, Dover Books onIntermediate and Advanced Mathematics, Dover Pub-lications, Mineola, N.Y., 1960.[5] S. Glasstone, R. Lovberg, Controlled ThermonuclearReactions: An Introduction to Theory and Experiment,R. E. Krieger Publishing Company, 1975.[6] S. Atzeni, J. Meyer-ter Vehn, The Physics of InertialFusion: BeamPlasma Interaction, Hydrodynamics, HotDense Matter, International Series of Monographs onPhysics, OUP Oxford, 2004.[7] F. H. Shu, The Physics of Astrophysics: Radiation, no.v. 1 in Series of books in astronomy, University ScienceBooks, 1991.[8] F. H. Shu, The Physics of Astrophysics: Gas dynam-ics, no. v. 2 in Series of books in astronomy, UniversityScience Books, 1991.[9] R. P. Drake, High-Energy-Density Physics: Fundamen-tals, Inertial Fusion, and Experimental Astrophysics,Shock Wave and High Pressure Phenomena, SpringerLondon, Limited, New York, N.Y., 2007.[10] R. B. Lowrie, J. E. Morel, J. A. Hittinger, The cou-pling of radiation and hydrodynamics, The Astrophys-ical Journal 521 (1) (1999) 432.[11] G. J. Habetler, B. J. Matkowsky, Uniform asymp-totic expansions in transport theory with small meanfree paths, and the diffusion approximation, Journal ofMathematical Physics 16 (4) (1975) 846–854.[12] E. W. Larsen, Neutron transport and diffusion in inho-mogeneous media. i, Journal of Mathematical Physics16 (7) (1975) 1421–1427.[13] E. W. Larsen, G. Pomraning, V. Badham, Asymp-totic analysis of radiative transfer problems, Journalof Quantitative Spectroscopy and Radiative Transfer29 (4) (1983) 285 – 310. [14] J. E. Morel, Diffusion-limit asymptotics of the trans-port equation, the P / equations, and two flux-limiteddiffusion theories, Journal of Quantitative Spectroscopyand Radiative Transfer 65 (5) (2000) 769 – 778.[15] F. Malvagie, G. Pomraning, Initial and boundary con-ditions for diffusive linear transport problems, Journalof Mathematical Physics 32 (1991) 805–820.[16] E. W. Larsen, J. B. Keller, Asymptotic solution ofneutron transport problems for small mean free paths,Journal of Mathematical Physics 15 (1) (1974) 75–81.[17] E. W. Larsen, J. E. Morel, W. F. Miller Jr., Asymp-totic solutions of numerical transport problems in op-tically thick, diffusive regimes I, Journal of Computa-tional Physics 69 (2) (1987) 283 – 324.[18] E. W. Larsen, J. E. Morel, Asymptotic solutions of nu-merical transport problems in optically thick, diffusiveregimes II, Journal of Computational Physics 83 (1)(1989) 212 – 236.[19] M. L. Adams, T. A. Wareing, W. F. Walters, Char-acteristic methods in thick diffusive problems, NuclearScience and Engineering 130 (1998) 18–46.[20] M. L. Adams, Discontinuous finite element transportsolutions in thick diffusive problems, Nuclear Scienceand Engineering 137 (2001) 298–333.[21] E. Larsen, G. Pomraning, Asymptotic analysis of non-linear marshak waves, SIAM Journal on Applied Math-ematics 39 (2) (1980) 201–212.[22] J. E. Morel, T. A. Wareing, K. Smith, A linear-discontinuous spatial differencing scheme for S n ra-diative transfer calculations, Journal of ComputationalPhysics 128 (1996) 445–462.[23] M. L. Adams, P. F. Nowak, Asymptotic analysisof a computational method for time- and frequency-dependent radiative transfer, Journal of ComputationalPhysics 146 (1) (1998) 366 – 403.[24] G. L. Olson, L. H. Auer, M. L. Hall, Diffusion, P , andother approximate forms of radiation transport, Jour-nal of Quantitative Spectroscopy and Radiative Trans-fer 64 (6) (2000) 619 – 634.[25] J. D. Densmore, E. W. Larsen, Asymptotic equilibriumdiffusion analysis of time-dependent monte carlo meth-ods for grey radiative transfer, Journal of Computa-tional Physics 199 (1) (2004) 175 – 204.[26] J. D. Densmore, Asymptotic analysis of the spatial dis-cretization of radiation absorption and re-emission inimplicit monte carlo, Journal of Computational Physics230 (4) (2011) 1116 – 1133.[27] S. Jin, Efficient asymptotic-preserving schemes for somemultiscale kinetic equations, SIAM Journal on ScientificComputing 21 (2) (1999) 441–454.[28] R. B. Lowrie, J. D. Edwards, Radiative shock solu-tions with grey nonequilibrium diffusion, Shock Waves18 (2008) 129–143.[29] R. B. Lowrie, R. M. Rauenzahn, Radiative shock solu-tions in the equilibrium diffusion limit, Shock Waves 16(2007) 445–453.[30] R. B. Lowrie, J. E. Morel, Issues with high-resolutiongodunov methods for radiation hydrodynamics, Jour-nal of Quantitative Spectroscopy and Radiative Trans-fer 69 (4) (2001) 475 – 489.[31] S. Chapman, T. Cowling, The Mathematical Theory ofNon-uniform Gases: An Account of the Kinetic The-ory of Viscosity, Thermal Conduction, and Diffusion inGases, Cambridge University Press, 1958.[32] D. Mihalas, R. I. Klein, On the solution of the time- ependent inertial-frame equation of radiative transferin moving media to O(v/c), Journal of ComputationalPhysics 46 (1982) 97–137.ependent inertial-frame equation of radiative transferin moving media to O(v/c), Journal of ComputationalPhysics 46 (1982) 97–137.