Perturbed, Entropy-Based Closure for Radiative Transfer

We derive a hierarchy of closures based on perturbations of well-known entropy-based closures; we therefore refer to them as perturbed entropy-based models. Our derivation reveals final equations containing an additional convective and diffusive term which are added to the flux term of the standard closure. We present numerical simulations for the simplest member of the hierarchy, the perturbed M1 or PM1 model, in one spatial dimension. Simulations are performed using a Runge-Kutta discontinuous Galerkin method with special limiters that guarantee the realizability of the moment variables and the positivity of the material temperature. Improvements to the standard M1 model are observed in cases where unphysical shocks develop in the M1 model.


Introduction
In this paper, we derive a new hierarchy of kinetic moment models in the context of frequencyintegrated (grey) photon transport. These new models are perturbations of well known entropy-based models; we therefore refer to them as perturbed entropy-based or PEB models. We present numerical simulations for the simplest member of the hierarchy, the perturbed M 1 or PM 1 model, in one spatial dimension. In this setting, the PM 1 model approximates the evolution of the photon radiation energy E and radiation flux F through a material medium with slab geometry. The photons interact with the material through scattering and emission/absorption processes.
Entropy-based (EB) models have been studied extensively in areas such as extended thermodynamics [20,45], gas dynamics [26,29,32,33,35,38,55], semiconductors [2-5, 27, 31, 34, 36, 51], quantum fluids [19,21], radiation transport [10, 11, 13, 14, 22-24, 28, 30, 43, 44, 59, 63], and phonon transport in solids [21]. In the context of radiative transfer, entropy-based models are commonly referred to as M N models, where N is order of the expansion. The M 1 model dates back to [43], where it was first derived using Maxwell-Boltzmann statistics. For problems with Bose-Einstein statistics, formal theoretical properties such as hyperbolicity and entropy dissipation were first reported in [22] for arbitrary N . However, computational studies have focused primarily on properties of the M 1 model and its extensions, including multigroup equations [62] and partial moment models [23,24]. In related work, one may find simulations of M 1 models based on other statistics, including Maxwell-Boltzmann [9][10][11] and Fermi-Dirac [8,58]. This attachment to M 1 is due to the fact that the higher order members of the M N hierarchy require the repeated solution of expensive numerical optimization problems. However, simulations of the M 2 model [44,63] (the next member in the hierarchy) have been performed for Bose-Einstein statistics and for M N up to order N = 15 for special benchmark problems using Maxwell-Boltzmann statistics [1,28].
There are several reasons to consider perturbative modifications to EB models. First, it is more economical to improve the model with perturbative corrections than to increase the number of moments, since the latter increases the memory footprint and makes the defining optimization problem more difficult to solve. In this case of grey photon transport, the flux in the M 1 model can be expressed analytically, i.e., no direct solution of the optimization problem is required. Thus, in this case, the argument against increasing the number of moments is especially compelling. A second reason is that perturbations add (among other things) diffusive terms to the EB model. It is hoped that these terms will smooth out nonphysical shocks which are known to exist in EB models. These shocks are a generic artifact of the modeling procedure that result when approximating linear transport in phase space by a nonlinear hyperbolic balance law for a set of moments. A third reason is that the specification of boundary conditions for moment equations that are consistent with the underlying kinetic boundary conditions is an open problem. However, at least in the case of linear moment equations, recent efforts [37] have shown the potential for well-posed boundary conditions for models with perturbative corrections. A fourth and final reason is that entropy-based closure do not depend on the properties of the material. Perturbations on the other hand can couple material properties into the closure.
Our goal in this work is to assess the qualitative behavior of the PM 1 model relative to the original M 1 model. We consider several test cases and find that the PM 1 model gives mixed results. Roughly speaking, it does quite well for shock problems that the M 1 model cannot handle. However, for more regular solutions, the two models perform comparably; and in some cases, the M 1 model performs slightly better.
One of the fundamental questions associated with any moment model is the issue of realizability. In the context of the M 1 and P M 1 models, the two unknowns E and F are called realizable if and only if they are the first two moments of an underlying kinetic distribution. This requirement on E and F is mathematically equivalent to the condition which must be satisfied point-wise in space and time. Here, c is the speed of light. It is expected that the solutions of the M 1 model will satisfy (1) because it (like all EB models) is derived assuming an ansatz for the kinetic distribution which is positive. However, the underlying ansatz for the PM 1 model is a perturbation of the EB ansatz that is no necessarily positive. Therefore, a modification of the PEB ansatz is needed which controls the contribution of the perturbative term. Even for the M 1 model, the realizability condition (1) can be destroyed by a numerical method unless special care is taken to enforce it. To address this issue in the current setting, we build on previous work with the M 1 model [48], using a Runge-Kutta discontinuous Galerkin (RKDG) method that is equipped with a special slope limiter [66,67] in the spatial variable. For implementation of the PM 1 model, this special limiter must be applied in combination with a control parameter that limits the size of the perturbations in the underlying ansatz of the PEB closure. The RKDG method [6] is a natural discretization here because we deal with a hyperbolic system of equations that is augmented by a diffusive term.
The remainder of the paper is organized as follows. In Section 2, we introduce the radiative transfer equation and moment model framework. In Section 3, we derive perturbed entropybased closures and give explicit expressions for the perturbed M 1 model. In Section 5, we focus on the P M 1 model in slab geometry and give details of the discontinuous Galerkin method used for simulation. In Section 6, we present numerical results. Section 7 is for discussion and conclusions. Several calculations and proofs are relegated to the Appendix.

Radiative Transfer and Moment Equations
We consider a collection of photons which move at the speed of light c through a static material medium. In engineering and physics applications, the fundamental quantity of interest is the radiation intensity ψ = ψ(x, Ω, ν, t) which is a function of position x ∈ K ⊂ R 3 , direction Ω ∈ S 2 , frequency ν ∈ (0, ∞), and time t ∈ (0, ∞). Roughly speaking, ψ is the flux of energy through a surface normal to Ω. If f is the kinetic density of photons-that is, the number density with respect to the Lebesgue measure dxdΩdν-then ψ = hνcf , where h is Planck's constant.
The material is characterized by a temperature T = T (x), an equation of state for the energy e = e(T ), and by scattering, absorption, and total cross-sections: σ s , σ a , and σ t = σ a + σ s that depend on x directly and also indirectly through the material temperature.

The Radiative Transfer Equation
The radiative transfer equation, which approximates the evolution of ψ, is given by The collision operator C models interactions of photons with the medium. For the purposes of this paper, we assume C has the form where φ is the angular integral of ψ: and the Planckian models blackbody radiation from the material. The constant k is Boltzmann's constant. The first term in C accounts for the loss of photons at a given frequency and angle due to both out-scattering and absorption by the material. The second group of terms gives the gain of photons due to in-scattering from other angles, re-emission by the material, and a generic external source s. To make calculations, s is assumed to be isotropic. However, such an assumption is not necessary. The evolution of the material temperature is determined by a balance of emitting and absorbed photons: where angle brackets are used as a shorthand notation for integration over angle and frequency: and the T 4 term in the first equation comes from the Stefan-Boltzmann Law: The constant a is the radiation constant. Though the material equation (6) plays an important role, we will focus here on simulating the transport equation (2).

Moment Equations
The large phase space on which (2) is defined makes direct numerical simulation prohibitively expensive. Thus, approximate models are needed to reduce the size of the system. A common and well-known approach is the method of moments, for which the angular and/or frequency dependency of ψ is approximated using a finite number of weighted averages. Derivation of any moment system begins with the choice of a vector-valued function m : S 2 → R n , Ω → [m 0 (Ω), . . . , m n−1 (Ω)] T , whose n components are linearly independent functions of Ω. Evolution equations for the moments u(x, t) := mψ(x, ·, t) are found by multiplying the transport equation by m and integrating over all angles to give 1 c ∂ t u + ∇ x · Ωmψ = mC(ψ; T ) .
The system (9) is not closed; a recipe, or closure, must be prescribed to express unknown quantities in terms of the given moments. Often this is done via an approximation for ψ in (9) that depends on u, and satisfies the consistency relation The resulting moment system is 1 In general, a closure is required to evaluate both the flux terms and the collision terms in (9). However for the collision operator in (3), no closure is required. Indeed, it is straightforward to show that mC(E(u); T ) = mC(ψ; T ) for any reconstruction that satisfies the consistency relation. Thus, for the purposes of this paper, we will be focused on closure of the flux term. As one might expect, the behavior of a moment system-and in particular its ability to capture fundamental features of the kinetic description-depends heavily on the form of the reconstruction.

Entropy-Based Closures
A general strategy for prescribing a closure is to use the solution of a constrained optimization problem min g∈Dom(H) where H(g) := η(g) and η : R → R is a strictly convex function that is related to the entropy of the system. For photons, the physically relevant entropy comes from Bose-Einstein statistics and is given by [49,52] η(g) = 2kν 2 c 3 [n g log(n g ) − (n g + 1) log(n g + 1)] , where the n g is the occupation number associated with g: The solution of (13) is expressed in terms of the Legendre dual Let Then we have the following.
Theorem 1. The solution of (13) is given by B(α), whereα =α(u) solves the dual problem It is also the Legendre dual variable of u with respect to the strictly convex entropy h(u) := H(B(α(u))), i.e.,α The moment system derived by setting E(u) = B(α) in (12) is hyperbolic and symmetric when expressed in theα variables and its solution formally dissipates h. Moreover, E is an inherently positive quantity.
Proof. The form of the minimizer in (18) can be derived formally using standard Lagrange multiplier techniques. However, a rigorous proof requires more technical arguments which can be found, for example, in [33] for the Maxwell-Boltzmann entropy and applied directly to the current setting. Once the existence of a minimizer is found, the other properties can be verified, as is done in [22,35]

Perturbed Entropy-Based (PEB) Closures
Perturbations to standard P N closures 1 have been derived for N = 3 in [47] and for general N in [54] (see also [27] and [60] ). The idea behind the derivation in [54] is to write ψ = ψ pn +ψ, where ψ pn is the standard P N expansion. The perturbationψ satisfies its own kinetic equation, which can be then used to approximateψ in terms of ψ pn . The resulting "D N " models gain a diffusive term in the equations for the highest order moments. uch an approach need not be restricted to the P N equations. Indeed, following this exact strategy, we define 1. The moment map M : g → u := mg ; 2. The expansion map E : u → B(α(u)); 3. The reconstruction R = E • M; 4. The kinetic perturbationψ = ψ − R(ψ).

The kinetic equation forψ is
where and we have used the relation to compute the matrix ∂α ∂u in (22). By operating withP u := I − P u on (2) , where P u := E ′ (u)M, we can write (21) as It should be noted, for future use, that the projection Q u , given by is self-adjoint in L 2 with respect to the positive weight W(u). Equation (25) for the perturbation is exact. To derive a closure, we neglect the time derivative and perturbative component of the flux to arrive at the following approximate balance equationP whereP In the appendix, we show that for the grey equations with Bose Einstein entropy, Knowing that this component will be integrated out in the final closure, we therefore solve (27) forP u E(u) +ψ in terms of a convective componentψ c and a diffusive componentψ d : where r s and r a are the scattering and absorption ratios, respectively: Inserting (30) back into the flux term of the moment equation (12) gives At this point, it is not clear whether this flux dissipates an entropy or if the convective flux f C is always hyperbolic. In general, the hyperbolicity of moment models closed by an entropy minimization principle follows from the fact that (in terms the Lagrange multiplierŝ α) the model can be written as a symmetric Lax-Friedrichs form [35]. This structure is not present here. However, at least for slab geometries, the convective flux in the P M 1 model is hyperbolic. (See Proposition 2 in the following section.) Moreover, in general, the diffusive flux satisfies a local dissipation law. Proof. A dissipation law for h is found by multiplying the closed moment system (12) where ∇ x acts on the components of Ω and the Lagrange multiplierα T on m. We only need to work with the term that is not in divergence form. We use the fact that whereQ u := Id − Q u and Q u is given in (26).

Controlling the Perturbations
While the entropy-based ansatz in (18) is positive for all Ω, the addition of the perturbation in (30) may lead to an ansatz which is not. As a consequence, the moments of the perturbed ansatz may not satisfy the realizability condition (1). To correct for this defect, we introduce a modification and approximate ψ with where δ(x, t) is a scalar control parameter. Several different choices for δ are possible. For example, one could select it to ensure that E(u) is positive everywhere. However, this choice requires pointwise evaluations with respect to Ω-a task we would like to avoid. Instead, we select δ in such a way as to preserve (1) in the numerical computation. While the exact form of δ depends on the details of the numerical method, the general framework relies on the realizability conditions for the moments. We call an array (Ψ 0 , Ψ 1 , . . . , Ψ N ) realizable with respect to (1, Ω, . . . , Ω ⊗N ) if there exists a non-negative measure on dΩdν with density Ψ(Ω, ν) such that Ψ k = Ω ⊗k Ψ for k = 1, . . . , N . The set R N of all such vectors is called the realizable set. Roughly speaking, we select δ to ensure that (Ψ 0 , Ψ 1 , . . . , Ψ N ) ∈ R N . Note that such a δ always exists: When δ = 0, there is no perturbative term and since the minimum entropy ansatz is always positive, (Ψ 0 , Ψ 1 , . . . , Ψ N ) ∈ R N . Details for the P M 1 are given in Section 5.3.

The Perturbed M 1 (PM 1 ) model
The perturbed M 1 model is based on the moments where E is the photon energy density and F is the energy flux density. The model approximates the evolution of E and F with the following system: where S(x, t) := ∞ 0 s(x, ν, t) dν and the closure for the pressure term is Here Π M 1 (u) is the term that comes from the entropy ansatz (the entropy-based term). The term Π C (u) is the convective correction and Π D (u) is the diffusive correction. These corrections can be expressed in terms of Π M 1 and which, in turn, can be expressed in terms of the unit vector n := F/|F | and the scalars Lemma 1. The correction terms Π C and Π D are given by Proof. see appendix.
Remark 1. The formula for the convective correction is independent of the specific form of E. In particular for the P 1 model, the pressure term is . This is consistent with the fact that the "D N " models in [54] contain only diffusive corrections.
Lemma 2. The entropy-based terms Π M 1 and Q M 1 are given by where the scalars χ 1 , χ 2 , and χ 3 are defined in (38).
Proof. See the appendix.
In slab geometry, we end up with the following expressions for the components of the pressure term Π = Π M 1 + Π C + Π D which can be computed from Lemma 1 and Lemma 2: where the convection and diffusion coefficients are given by and These coefficients are displayed in Figure 1. Note that χ, η, and D F are all even functions of the ratio F/(cE), while D E is odd.

Numerical Simulation using Discontinuous Galerkin
In slab geometries, the diffusion-corrected M 1 model and the material energy (6) reduce to ∂T is the specific heat at constant volume. In this setting, the convective flux of the P M 1 model is hyperbolic. Proposition 2. The perturbed M 1 system in slab geometry is hyperbolic if Π D = 0 and |F | < cE.
Proof. The proof is a direct calculation, given in the appendix. We simulate the system (47) using a Runge-Kutta discontinuous Galerkin (RKDG) method. The RKDG method is a method of lines: the DG discretization is only applied to spatial variables while time discretization is achieved by explicit Runge-Kutta time integrators. The presentation here is rather brief and relies on details found in [48], where the method was applied to the M 1 model. A general description of the RKDG method can be found, for example, in [16,17].

Spatial Discretization
We divide the computational domain (x L , x R ) into J cells with edges and let x j denote the center of each cell I j = (x j−1/2 , x j+1/2 ). We let h j := x j+1/2 − x j−1/2 be the length of the interval I j and set h := max j h j . We denote the finite-dimensional approximation space by where P k (I j ) is the space of polynomials of degree at most k on the interval I j .
The semidiscrete DG scheme is derived from a weak formulation of (47). However, following [18] we first reduce the convection-diffusion equations (47) to a system of first-order equations by introducing the auxiliary variable v: The exact solutions u(·, t), v(·, t) and T (·, t) are then replaced by approximations

and the resulting set of equations is required to hold for all test functions
Here we use the bracket notation: where and are the right and left limits of b h at the cell interfaces x j±1/2 . The term ub h (x) j is defined in an analogous fashion.
Since the components of u h (., t) and v h (., t) are piecewise polynomials, the edge values of u and v in (51) are not strictly defined. Thus, the nonlinear flux function f is replaced by a numerical fluxf which depends on the pointwise limits of u h , v h on either side of the edge at x j±1/2 :f The notations forû carry over analogously.
It remains to choose suitable numerical fluxesf andû. Since (47) has both a convective flux and a diffusive flux the choice is not obvious. Several approaches have been presented in literature [6,40,46,57]. In [6], the prescription for the diffusive term is given bŷ Combining this term with the Lax-Friedrichs flux for f C (u) gives the following total numerical flux:f where λ is the largest magnitude of any eigenvalue of the Jacobian associated with f C . These eigenvalues, in general, depend on material properties, the temperature T and the source term S. In contrast to the M1 model, they are not bounded by the speed of light c. For example, neglecting the temperature and source the maximum value is approximately 9.12 c.
We instead use the smaller value of λ = c, which is the particle speed in the transport equation and is consistent with the application of the control parameter to enforce realizability (see Section 3.3). Figure 2 shows the comparison of eigenvalues for the M1 and perturbed M1 modeling when c = 1, σ s = 1, σ t = 3, T = 0 = S. for P k (I j ) in each cell I j : The standard choice of basis for P(I j ) is generated by Legendre polynomials P l that are defined on the reference cell [−1, 1]: and normalized so that With ξ j (y) := x j + yh j /2, this gives a formulation defined on the reference cell: The remaining integrals are calculated by a quadrature rule. Note that (61b) can be solved locally for v j m (t) in each cell I j , which can then be substituted back into (61a).

Time Discretization: Explicit SSP Runge-Kutta Schemes
The purpose of high-order, strong stability preserving (SSP) Runge-Kutta time integration methods is to achieve high-order accuracy in time while preserving desirable properties of the forward Euler method (for a review, see [25]). In this paper, we only use explicit schemes, which compute values of the unknowns at several intermediate stages. Each stage is a convex combination of forward Euler operators and this usually leads to modified CFL restrictions. The equations in (61) form a system of ODEs for the coefficients u j m (t) and T j m (t). For all j ∈ {1, . . . , J} and m ∈ {0, . . . , k}, we write this system in the abstract form: Here, L j u,m and L j T,m are the respective right-hand sides of the ODEs. Let {t n } N n=0 be an equidistant partition of [0, t final ] and set ∆t := t final /N . Let Λ denote the application of a generic slope limiter, and let π V m h be the orthogonal projection onto the finite dimensional space V m h . Note that Λ is applied at every Runge-Kutta stage. The algorithm for the optimal third-order SSP Runge-Kutta (SSPRK(3,3)) method [56] reads as follows: • m .
For the sake of completeness, we also state the optimal fourth order scheme SSPRK(5,4) [25]: Note that SSPRK(3,3) permits a timestep of the same size as forward Euler, while the SSPRK(5,4) method is less restrictive, allowing for a time step that is 1.508 times larger the forward Euler scheme.

Limiters
As in [48], two types of limiters are used. The first is standard; it is used to suppress spurious oscillations and maintain stability. There are many such limiters available. In this paper, we apply the moment limiter from [12], which is a modification of the original limiter in [7]. This limiter is applied to the variables u, but not the auxiliary variables v or the temperature T . Additional details can be found in [48].

Realizability-Preserving Limiter
The second limiter is a realizability-preserving limiter which is needed to ensure that the cell averages of E and F satisfy the realizability condition (1) at each stage of the numerical computation. The limiter is based on the work from [67] and [66] and is very similar to what was done in [48] for the M 1 model. The major difference here is the addition of the control parameter δ. An essential ingredient of the realizability limiter is the Gauss-Lobatto quadrature set where, for a spatial reconstruction of order k, M is the smallest integer such that 2M − 3 ≥ 2k + 1. This condition on M ensures accuracy of the scheme [15]. The weaker condition 2M − 3 ≥ k ensures that the quadrature integrates elements of the approximation space V k h exactly.
The realizability limiter is defined in order to ensure that u h (x ℓ j ) ∈ R 2 at each pointx ℓ j in the quadrature set. However, we enforce the convexity condition indirectly by requiring the positivity of the intermediate quantities 2 The inverse transformation that maps (Q, R) → (cE, F ) is given by An additional limiter is also used to enforce the positivity of the temperature reconstruction at each quadrature point. We now proceed to define the limiters. Let u j,n h = (cE j,n h , F j,n h ) and T j,n h be the approximations of u and T in cell I j at time t n , and letû j,n h andT j,n h denote the modifications of u j,n h and T j,n h that are generated by the limiting. We assume that the cell average of u j,n h , which we denote by u j,n h , is realizable, i.e., u j,n h ∈ R 2 . We also assume that the cell average of T j,n h , which we denote by T j,n h , is positive. Let Q j,n h (x) and R j,n h (x) be the approximations of Q and R, respectively, and define limited variables bŷ where θ j,n Q := min The parameter ε > 0 is chosen to maintain numerical stability with finite precision arithmetic; its value should be small relative to the magnitude of the variables in a given problem. The components ofû j,n h are then defined using (69). They satisfy the following property which is a key ingredient for maintaining realizability in the RKDG scheme.

Setting the Control Parameter
We now define the control parameter δ, discussed in Section 3.3. Our definition is guided by the following result.
These relations are easily verified by applying the definition of δ and using the fact that (cE, F, cΠ M 1 ) ∈ R 3 .
With δ given by (71), one can show that cell averages of u h remain realizable and that the cell average of T h remains positive in a forward Euler step. Let Lemma 6. Assume that 2M − 3 ≥ k and for each ℓ = 1, . . . , M , that and Π j,n δ,ℓ satisfies (C1) and (C2). Assume further that ∆t satisfies the following conditions: where h := min j h j . Then after a forward Euler time step, Proof. We refer the reader to the proof of Theorem 3 in [48] for the M 1 model, which relies exactly on the conditions (A1)-(A3) and (C1)-(C2). The only difference is that (C1) and (C2) are assumed in Lemma 6, while in [48]  Proof. Application of the limiters in (70) ensures that the conditions of Lemma 6 hold at each stage in the SSP-RK scheme. Each successive stage is an application of the forward Euler operator to the current stage with an appropriately modified time step. Thus, the conclusions of Lemma 6 apply at the next stage, including the final stage, which gives u j,n+1 h .

Numerical Results
In this section, we present several numerical computations in slab geometry for a choice of test cases that are common for the M 1 model. The goal is to compare and contrast the perturbed M 1 model with the M 1 model and to point out benefits and drawbacks. Benchmark solutions are generated by the discrete ordinates method with an upwind scheme in space, high-order spherical harmonics or semi-analytic expressions. The RKDG implementation has been verified and benchmarked in [48]. The correct implementation of the additional perturbative terms has been checked by the method of manufactured solutions [53].
As in [48] our algorithm is implemented in MATLAB, and Gauss-Lobatto quadrature on [-1,1] is used. Additionally, the Runge-Kutta time integration methods as well as parameters for the admissibility limiter are applied in the same way. In order to satisfy the conditions of Lemma 6 and to guarantee stability, the time step is set to where k is the polynomial degree, h = min j h j , w min and w max are the minimum and maximum quadrature weights, respectively. The quantities σ t,max and τ max are the maximum values of σ t,ℓ and σ a,ℓ (T j,n ℓ ) 3 . The constant c 4 is not needed to preserve realizability of the moments but rather to enforce a parabolic CFL condition. Without this condition, unstable modes grow without bound until the control parameter δ turns on and damps them. Unfortunately, the parabolic CFL restriction leads to small time steps.
The stability parameter for the realizability limiter of E and F is set to ε = 10 −10 . The same value is also used to enforce conditions (C1) and (C2), i.e., the control parameter in (71) is chosen such that cΠ δ ≤ cE − ε and |F ± cΠ δ | ≤ cE ± F ± ε.
In Sections 6.1-6.3, we study simulations with c = 1 and neglect the energy equation which is included in the last two cases from Section 6.4. Unless otherwise stated, slope and realizability limiters are always turned on for all DG calculations. If transformation to characteristic variables for the slope limiter is used, it will be explicitly stated. For this problem, coupling to the material is ignored, and the material energy equation is not included in the simulation.

Two-Beam Instability
In Figure 3, one can observe the formation of a shock in the M 1 solution which persists at the steady-state. The perturbed M 1 model also develops an unphysical transient profile in which the particle number jumps in the center of the domain. While this artifact persists, the steady state solution (t = 3) appears continuous. The steady-state solution also has noticeable kinks in the at x ≈ ±0.3. For comparison, discrete ordinates solutions are plotted for which 256 discretization points in angle and 1000 points in space are used. The perturbed M 1 is throughout closer to the transport solution.
Remark 2. Precise explanations for the occurrence of shocks and kinks in the perturbed M 1 solution require an additional analysis. For example, an explanation for the formation of shocks in the M 1 model is given [10]. However, such an analysis goes beyond the purpose of this paper and will be postponed to future work.   3] hits an isotropic source S = 1/2 generating particles on the interval 1 ≤ x < 1.5. In order to avoid complications of spatial discontinuities in the fluxes [41,64,65], we smooth the source and material cross sections, which enter into the perturbative components of the flux. The source S is smoothed at the end points x = 1 and x = 1.5

Source-Beam Problem
Similarly, we design the material properties with the cross sections: The function p H is a Hermite polynomial of order 10 with p H (±1) = ±1 and p are set. Initially, there are (almost) no particles in the system, i.e., The value of c is again set to one. The perturbed M 1 results are compared to M 1 and transport solutions. Classic M 1 calculations are performed using the DG method from [48], with the same computational parameters as the P M 1 model and slope limiting performed in the characteristic variables. The transport solution is computed using the discrete ordinates method with 600 spatial cells and 256 discrete angles.
One can observe in Figure 4 that as time increases, particles entering from the left boundary encounter the source in the interior. As this happens the M 1 profile diverges from the transport solution. Even as steady state is achieved at t = 4, there is a large difference for x ≤ 1. The P M 1 profile agrees much better with the transport solution.

Gaussian Source
The next test case simulates particles with an initial energy density that is a Gaussian distribution in space and a zero energy density flux: Periodic boundary conditions on [−L, L] are prescribed where L = t final + 1. The computational domain is always chosen large enough to ensure that a negligible number of particles reaches the boundaries. No internal source is assumed (so S = 0), and the medium is purely scattering with σ s = σ t = 1. The velocity c is also set to one and the material energy equation (6) is neglected. All DG results are computed with h = 0.01 and polynomial degree k = 2. For comparison, discrete ordinates solutions of the transport equations are obtained with h = 0.005 and 128 angular points. Figure 6 displays the solutions at t final = 1, 2, 3, 10. The M 1 model gives the expected wave effects that are washed out at larger times. These effects do not occur in the perturbed M 1 results. However, the P M 1 forms Gaussian bell that are higher and more narrow than the benchmark solution. At lower times, their maximum propagation speed is roughly half the correct velocity. Nevertheless, at t final = 10 the front of the P M 1 model catches up with the reference solution, at which point all three models agree reasonably well.
The perturbationψ from Section 3.2 is related to the difference between the M 1 and transport solution. Figure 6 indicates that this quantity is highly time-dependent. Additionally, the spatial gradient ofψ is large at shorter times. Hence, this numerical example violates the smallness assumptions made in the derivation of the perturbed M 1 model in Section 3.2. Thus the lack of accuracy is not surprising.

Including the Material Energy Equation
The next two examples involve (2) coupling to the energy equation (6). The linearized Marshak wave problem from [61] is analyzed first and then a Marshak wave with material parameters taken from [48].

Smoothed Su-Olson's Benchmark Problem
In [61], tabulated data is provided for analytic solutions to a linearized Marshak wave problem, which serves as a validation of numerical algorithms in the radiative transfer community. In particular, this semi-analytic benchmark is also compared to diffusion-corrected P N approximations in [54]. It is therefore of interest to study solutions of the perturbed M 1 model to this problem.
We compute approximations to (2) and (6) in slab geometry with the following physical data As in the source-beam problem, we seek to avoid the discontinuous material properties in [61] by introducing a smoothed version:  where p H is the Hermite polynomial described in Section 6.2. However, we still compare our numerical results to the semi-analytic solutions from [61] because the length of the smoothing window 2∆ = 0.02 is relatively small.
Initially, the medium is cold and there is no radiation: Only at t = 3.16228 solutions from both models are close to each other as well as to the semi-analytic points.

Thin Marshak Wave
In this problem, incoming radiation is prescribed on the left boundary by well-posed boundary conditions, and the material is assumed to be purely absorbing: Due to above incoming radiation on the left boundary, radiation propagates through the medium from left to the right. The material temperature T (Figure 8a) and the energy density E (Figure 8a) decay smoothly to zero. The M 1 and PM 1 models yield very similar solutions. For comparison, a reference solution is computed using a P 99 model that is calculated with the DG method from [42]. The simulation of the P 99 model uses linear elements and 800 spatial cells. The reference solution shows a much stronger decay in the energy and material temperature.

Discussion and Conclusions
In this paper, we have derived a hierarchy of closures based on perturbations of well-known entropy-based closures. The derivation has been done in the context of grey photon transport. Our derivations reveals final equations containing an additional convective and diffusive term which are added to the flux term of the standard closure. This is different to perturbations to standard P N closures [54] which only gain a diffusive component.
For the first member of the hierarchy, the PM 1 model, we compute explicit formulas for all terms. The resulting system of equations is a convection-diffusion system which, for slab geometries, is discretized by using a Runge-Kutta discontinuous Galerkin method. By introducing a special limiter and an additional control parameter to modify the pressure term, we ensure that cell averages of the moments satisfy the important realizability property (1).   We perform simulations to compare qualitative results with the M 1 model and with highlyresolved discretizations of the original transport equation. Improvements to the standard M 1 model are observed in cases where unphysical shocks develop in the classical M 1 model. However, for problems with continuous solutions, there is little or no improvement; and in some cases, the M 1 model perform marginally better.
Finally, we discuss some open problems in this framework which might be addressed in future: • Moment systems from entropy-based closures are known to be hyperbolic and satisfy a local dissipation law [22,35], but such results are not yet known for the perturbed models. The partial result in Proposition 1 confirms that the additional diffusive term dissipates the entropy. Moreover, neglecting the diffusion contribution the system indeed becomes hyperbolic for the special case of the M 1 model in one dimension (i.e., for slab geometry).
• Standard issues in analysis such as existence and uniqueness of solutions have yet to be investigated for the P M 1 model or for PEB models in general.
• In Section 5.3.2, the control parameter δ is chosen to guarantee conditions (C1)-(C2). However, this ansatz is a crude modification of the original perturbative model and could distort numerical solutions. A more subtle approach, possibly along the lines of flux-limited diffusion, e.g. [39], would be preferable.
• An undesirable aspect of the RKDG method is the time step restriction required by the explicit time integrator. For convection-diffusion equations, stiff sources, and/or long time scales, this time step restriction is very harsh. To lower the computational effort, (semi-)implicit time discretizations are therefore necessary. In addition, the method does not address the challenges of spatially discontinuous fluxes that arise from discontinuities in sources and material cross-sections.
• Another issue is the formation of unphysical shocks in the (perturbed) M 1 solution. Further analysis is needed to understand why and when they appear and how to further mitigate them.

A Appendix
Computation of equation (29). We first calculate two frequency integrals. Let κ := hc kα T m and θ := −κν. Then With these two integrals, it is easy to show that −α/4 is the unique solution to the linear system so that Using the definition of P u , and again the integrals above: (87) Proof of Lemma 1. The proof is a straight-forward calculation. It turns out to be more efficient to calculate Π D and Π C without directly using (22). Instead for any function g, we compute (Ω ∨ Ω)P u g = (Ω ∨ Ω)g − Using (88), we find for the diffusive correction, For the convection correction, we use the fact that φ, B and S are independent of Ω. This implies that m 1 φ ≡ Ωφ = 0 and similarly for B and S. We also use the Stefan Boltzmann-Law (8) and the identity u 0 = cE = ∞ 0 φ dν. This gives Proof of Lemma 2. Let {e 1 , e 2 , e 3 } be any orthogonal basis for R 3 . Then and Now set e 3 = n = F/|F | and note that, according to Lemma 7 below, E(u) depends on Ω only through Ω 3 . Thus, only the terms with even powers of Ω 1 and Ω 2 will survive. For k = 2, this means cΠ M 1 = Ω 2 1 E(u) e 1 ∨ e 1 + Ω 2 2 E(u) e 2 ∨ e 2 + Ω 2 3 E(u) n ∨ n, and for k = 3, Q M 1 = 3 Ω 2 1 Ω 3 E(u) e 1 ∨ e 1 ∨ n + 3 Ω 2 2 Ω 3 E(u) e 2 ∨ e 2 ∨ n + Ω 3 3 E(u) n ∨3 .
The goal then is to write these formulas in terms of Ω 3 only. Let us focus first on Π M 1 . Because E(u) depends only on Ω 3 , symmetry arguments can be used to conclude that first two terms in (93) are the same. Combined with the far right relation (91), this gives cΠ M 1 = Ω 2 1 E(u) (e 1 ∨ e 1 + e 2 ∨ e 2 ) + (Ω 2 3 E(u) n ∨ n = Ω 2 1 E(u) (e 1 ∨ e 1 + e 2 ∨ e 2 + n ∨ n) + (Ω 2 3 − Ω 2 1 )E(u) n ∨ n where we have used the fact that e 1 ∨ e 1 + e 2 ∨ e 2 + n ∨ n is the identity. From the definition of χ 2 , we conclude that Similarly for k = 3, Q M 1 = 3 Ω 2 1 Ω 3 E(u) (e 1 ∨ e 1 + e 2 ∨ e 2 + n ∨ n) ∨ n + (Ω 2 3 − 3Ω 2 1 )Ω 3 E(u) n ∨3 Lemma 7. For the M 1 model, the multiplierα 1 is co-linear with F , that iŝ Proof. If E(u) = η ′ * − hνc k (α 0 +α 1 m 1 ) solves the optimization problem (13), then by definition Let R be any orthogonal 3 × 3 matrix which preserves F . Then multiplying (99) by R gives where we have used the fact that the measure dΩ is invariant under the action of R. Because the solution of the optimization is unique, we conclude that Rα 1 =α 1 and therefore, since R is arbitrary,α 1 and F must be co-linear hjhkjhkj We show that the radical α + β 2 /4 in the formula for the eigenvalues is positive for all f = 1. Note that (39b) implies η = 1/3 + χ ′ f − χ and hence, The prime notation always refers to the derivative with respect to f . With this, we conclude