Approximating the $M_2$ method by the extended quadrature method of moments for radiative transfer in slab geometry

We consider the simplest member of the hierarchy of the extended quadrature 
method of moments (EQMOM), which gives equations for the zeroth-, first-, and 
second-order moments of the energy density of photons in the radiative 
transfer equations in slab geometry. 
First we show that the equations are well-defined for all moment vectors 
consistent with a nonnegative underlying distribution, and that the 
reconstruction is explicit and therefore computationally inexpensive. 
Second, we show that the resulting moment equations are hyperbolic. 
These two properties make this moment method quite similar to the attractive 
but far more expensive $M_2$ method. 
We confirm through numerical solutions to several benchmark problems that the 
methods give qualitatively similar results.

1. Introduction. The radiative transfer equations describe the interactions between photons and a background medium. The equations have been widely used in various applications such as atmospheric modeling, nuclear engineering and medical imaging. One challenge in simulating the radiative transfer equations is that the energy density of the photons depends on time, space, and an angle variable, and thus its state-space is high-dimensional. The equation for the energy density is a kinetic equation, and most numerical methods for kinetic equations fall in one of the following categories: direct simulation Monte Carlo (DSCMC), discrete-velocity models, and moment models. In this paper we will focus on moment models.
Moments, which are defined as angular averages of the energy density against polynomials, are natural quantities for the discretization of the angular dependence of the energy density because they are typically the quantities of interest from the solution. Moments are also the relevant quantities in a standard spectral discretization. The moment-closure problem is that of how to reconstruct the angular dependence of the energy density from a finite number of moments, thereby specifying a closed system of moment equations which approximate the exact moments of the radiative transfer equations. This reconstruction is called an ansatz. Ideally, the moment equations should have three properties: First, because the kinetic equation for the energy density is a continuum of coupled hyperbolic equations which satisfy certain conservation properties (for example conservation of mass), the moment equations should be a system of hyperbolic equations in conservative form. Notice that this also ensures local existence and uniqueness of solutions to the moment equations (at least in the one-dimensional case), and the plethora of numerical methods developed for simulation of hyperbolic conservation laws also become available. Second, because the energy density is inherently nonnegative, the moment equations should evolve in the realizable set, which is roughly defined as all those moment vectors which can be associated with a nonnegative density. Thus the reconstruction should yield a nonnegative ansatz for every realizable moment vector. Finally, the moment equations should satisfy an entropy dissipation law to be consistent with the second law of thermodynamics. 1 For the case of photons, the relevant entropy is the Bose-Einstein entropy.
The standard spectral method, colloquially known as the P n equations [13,10], where n indicates the degree of the highest-degree polynomials used to define the moments, indeed yields a set of hyperbolic equations in conservative form which dissipate the L 2 entropy, but their solution does not remain realizable [4,11]. Perhaps the only method known to possess all three properties (with dissipation of the Bose-Einstein entropy) is the M n method, where the energy density is reconstructed by solving a constrained strictly convex optimization problem (namely, that of minimizing the entropy subject to moment constraints) [9,7]. However, M n moment equations are expensive and difficult to implement because in general the defining optimization problem must be solved numerically and can be ill-posed [2,1], although they remain possibly interesting for large-scale parallel computation because this expense can be parallelized [8].
Another moment closure, the extended quadrature method of moments (EQ-MOM) was recently investigated in the context of radiative transfer in slab geometry [16]. The EQMOM uses an efficient procedure to a reconstruct convex combination of beta distributions to approximate the energy density from the moment vector, and good results were observed when the method was applied to standard test problems. However, it remained unclear whether the EQMOM is well-defined for all realizable moment vectors (that is, whether the reconstruction procedure always succeeded for realizable moment vectors) or whether the resulting set of moment equations is hyperbolic. In this work, we consider these questions for the simplest member of the EQMOM hierarchy, where a single beta distribution is used. The beta distribution has two degrees of freedom, so when combined with scaling (because the energy density need not be normalized), there are three degrees of freedom, and thus equations for the moments of zeroth-, first-, and second-order are derived. We show that the reconstruction procedure can be made explicit (that is, no iterative procedure is necessary), that the reconstruction is well-defined and nonnegative for all realizable moments, and finally that the resulting moment equations are hyperbolic. Thus this particular method shares two key properties with the M 2 method, and so we compare its solutions to those of the M 2 method.
The rest of this paper is arranged as follows: In Section 2 we introduce the basics of moment models, including the M n method and the EQMOM. Then in Section 3 we introduce the simplest member of the EQMOM hierarchy and analyze its properties. In Section 4 we compare this EQMOM method with M 2 using several benchmark problems. Finally in Section 5 we summarize and draw conclusions.
2. Preliminaries. We consider the one-dimensional case of slab geometry. Let I It, z, ν, µ be the specific intensity which depends on time t b R , the spatial coordinate z b z L , z R ¥ along the axis perpendicular to the slab, the angle variable µ b 1, 1¥, which gives the cosine of the angle between particle travel and the same axis, and frequency ν b R . Let T T t, z be the temperature of the background medium. The intensity and temperature are governed by the radiative transfer where c is the speed of light; σ a , σ s , and σ t ¢ σ a σ s are the absorption, scattering, and total interaction coefficients, respectively; is the Planckian; S is an external source; and For convenience we denote the sum of the terms on the right-hand side of (1a) Let m mµ m 0 µ, m 1 µ, . . . , m n µ T be a set of basis functions of the space H n of polynomials of µ up to degree n. The angular moments of the specific intensity I are mI e and they satisfy The term µmI e involves moments with respect to polynomials of order one larger than those in H n , and therefore (2) is not a closed system. A moment model is then defined by approximating the higher-order moments in terms of lower order moments to give a closed system of equations approximating (2). This closed system of equations has the form 1 c ∂E ∂t where 2 E¨mIe , f E¨µmI e , and rE¨mCI, T e .
The choice of f and r specify a closure.
The minimum entropy principle is one way of deriving a moment model which has several attractive theoretical properties. It is based on reconstructing an ansatẑ I M of I from the moments E by solving the following constrained variational minimization problem where HI is the Bose-Einstein entropy The solution of (5) has the form where α αE should be such that the polynomial α mµ d 0 for all µ b 1, 1¥, and α is the unique vector such that Î M me E. Then the M n method is defined by taking in (3).
The grey equations are obtained by integrating not just over angle but also over frequency, that is, replacing the e with t y in (4), (5), and (8) above. In this case, after integrating over frequency the entropy ansatzÎ M reduces tô where σ is the Stefan-Boltzmann constant. Properties of the M n model are discussed in [9,7], including a proof of its global hyperbolicity. Additionally, the M n ansatz is nonnegative, and therefore the moment vector E is always realizable. A moment vector E is said to be realizable if there exists a nonnegative distribution whose moments are equal to those in E. Finally, with the M n closure the moment equations (3) have an entropy dissipation law.
However, in general the M n closure is not given explicitly, that is, f E must be computed by solving the defining optimization problem (5) numerically. 3 This is typically solved through the unconstrained finite-dimensional dual to (5) which obtained by standard Lagrange-multiplier techniques. Unfortunately, this problem is expensive and difficult to solve numerically [2]. Due to these numerical difficulties, there has so far been no efficient general implementation of the M n model except for the simplest case, the M 1 method [3,12]. However, the M 1 method gives qualitatively incorrect solutions to some standard test problems, see for example [5].
The extended quadrature method of moments (EQMOM) [16] reconstructs an ansatzÎ B using a conical combination of beta distributions. The beta distribution on 1, 1¥ is given by where Bξ, η is the beta function. ThenÎ B is chosen of the form Notice that only a single value of the parameter δ is used. Therefore this ansatz has 2p 1 degrees of freedom, and so these parameters w i , γ i and δ can (theoretically) be chosen so that the moments ofÎ B up to order n 2p agree with those from a given vector E of moments. We will call this method the B n closure. An algorithm to determine the parameters w i , γ i , and δ from the moment vector E is given in [16], but it remains unclear whether this algorithm succeeds for all realizable moment vectors E. It is also unclear whether the resulting moment equations (3) with (or the corresponding definitions for the grey equations) are well-posed. We begin this analysis in the next section.
(a) The value of E 3 for normalized realizable moments using the B 2 closure (12).
(b) The value of E 3 for normalized realizable moments using the minimum entropy closure for the monochromatic case. This plot for the grey case is very similar.  1. If E 2 E 0 , then the specific intensity ansatz can only be I 1 Thus we must have E 3 E 1 which agrees with (12).
2. If E 2 E 0 E 2 1 , then the specific intensity can only be I E 0 δµE 1~E0 . This implies for E 3 E 3 1~E 2 0 , which also agrees with (12). Furthermore, this closure is correct in the isotropic case (that is, when E E 0 , 0, E 0~3 contains the moments of the constant density Iµ ¡ E 0~2 ).
In Figure 1(a) and Figure 1 We now show that the B 2 equations are globally hyperbolic in the interior of the realizable set M. We need the Jacobian matrix The characteristic polynomial of J is We then have the following theorem: Global strict hyperbolicity of B 2 ). The B 2 system (13) is globally strictly hyperbolic in the interior of the realizable set M, and its propagation speed is less than the speed of light.
thus pλ has one root in each intervals 1, Similar arguments work for E 1 d 0.
Denote λ 1 d λ 2 d λ 3 to be the three eigenvalues of the Jacobian matrix J , then it is clear that the eigenvector corresponding to λ j is R j ¡1 λ j λ 2 j ¦ T .
Theorem 3.2. The R 1 , R 3 -characteristic fields are genuinely non-linear, while the R 2 -characteristic field is neither genuinely non-linear nor linearly degenerate. 4 Proof. Let As ©λ j R j dp dλ λ λj 1 ∆λ, and dp dλ λ λj e 0, j 1, 3; dp dλ λ λ2 d 0, ©λ j R j 0 is equivalent to λ j being the common root of ∆λ and pλ. As the resultant of pλ and ∆λ is resp, ∆, λ 128 As Q and Z are both positive in the interior of the realizable region M, we have that resp, ∆, λ x 0 if E 1 x 0. Therefore when E 1 x 0, all characteristic fields satisfy ©λ j R j x 0.
In case that E 1 0, λ 2 0 is a common root of pλ and ∆λ, so Meanwhile we have thus ©λ j R j~ 0 for j 1, 3. Collecting the arguments above, one has that ©λ j R j~ 0, ¦E 0 , E 1 , E 2 b M, for j 1, 3; sign ©λ j R 2 signE 1 , under proper scaling of R 2 .
Let us summarize briefly features of the B 2 method: 1. The ansatz for the specific intensity preserves positivity and has explicit closure relationship. 2. The model is conservative and globally hyperbolic. 3. The signal speed of the B 2 model is less than the speed of light. 4. The closure is very close to that of M 2 . 5. Unlike M 2 , however, the closure can be stably computed for moment vectors up to and directly on ∂M, the boundary of the set of realizable moments.
Thus with the first two features, the B 2 system (13) has two of the important properties of M 2 . Unfortunately, it remains unclear whether the B 2 system also shares the last of these properties, that of entropy dissipation. 4. Numerical results. We solve equation (13) using the canonical finite volume scheme with the Lax-Friedrichs numerical flux, and the source term is treated implicitly. The number of cells are chosen to insure convergence of numerical solutions, and will be specified in the respective examples.
We derive boundary conditions on the kinetic level. On the left side of the domain, for example, where incoming kinetic data I L µ is provided, we define the components of the flux vectorf L whereÎ B is the ansatz associated with the moment vector from the first spatial cell on the left side of the domain. The flux vectorf R on the right side of the domain is analogously defined. We compute the integrals definingf L andf R using quadrature.
For the monochromatic cases below we simply set ν 1. We compare our model to numerical solutions of the M 2 method. We compute these solutions using the kinetic scheme and optimization techniques given in [1]. The entropy from that work is replaced with the Bose-Einstein entropy (6), and to avoid the singularity in the ansätze (7) or (9) when the polynomial α m passes through zero, we limit the step-size in the Armijo line search so that the polynomial α m remains negative at every angular quadrature point. These solutions are computed with 2000 cells except where stated otherwise. For solutions for the monochromatic case, where we do not consider the material temperature equation, we use the second-order kinetic scheme of [1]. For solutions to the grey equations, which we only include for the last test problem, we include the material temperature equation and update it implicitly as in [16]. Here we use the first-order version of the kinetic scheme.
All figures below present the numerical results for E 0 . This problem tests how a moment method handles crossing beams; it is the classical example where M 1 produces a nonphysical shock [5]. In Figure 3(a) we compare the steady-state solutions of the M 2 and B 2 methods. The B 2 solution is computed using 2000 cells. We see that the M 2 solution still appears to have the non-physical shock seen in M 1 solutions [5] while the the B 2 solution has no such shock.
Example 4.2 (Isotropic inflow into vacuum). In this example, we model a semiinfinite spatial domain z b ª, 1¥ with an isotropic inflow source imposed on the right boundary. We run the computation until t 0.8, so for numerical computations we use spatial domain 0, 1¥. We take σ a σ s 0. Initially, for all z, we take E 0 10 8 , E 1 0, E 2 E 0~3 . The isotropic inflow is specified at z 1. The specific intensity outside the right boundary is Iµ 0.5.
This problem tests how the moment method models simple kinetic transport from the boundary into the domain. The B 2 solution is computed using 16000  cells. The results are presented in Figure 3(b). Interestingly, the B 2 solution is smooth, while the M 2 solution has the expected structure of three waves. Example 4.3 (Plane source). In this test the spatial domain is unbounded (but 1, 1¥ for numerical simulations), the medium is purely scattering, σ s σ t 1, and the initial value is taken as E 0 z δz 10 8 , E 1 0 and E 2 E 0~3 .
This problem tests how the moment method handles strong spatial discontinuities. The B 2 solution is computed using 128000 cells. The numerical results are presented in Figure 3(c). While above the B 2 solutions have been smoother, here the B 2 solution has much sharper peaks at the fronts as the particles expand into the surrounding vacuum.
Because this is the only problem we are considering which neither has source terms nor is driven by boundary conditions, we computed the entropies of the numerical solutions and compare them in Figure 4. While we have not been able to prove that the B 2 solution dissipates entropy, we see in the figure that it does dissipate entropy for this problem, although not as quickly as the M 2 solution.
Example 4.4 (Su-Olson test case). In this problem, as in [16], we use pure absorption, σ a σ t 1, and eT aT 4 . The spatial domain is assumed to be unbounded, and in computing we take z b 0, 30¥, and imposing reflective boundary condition on z 0. The initial value is set as E 0 ac ! 10 10 , E 1 0, E 2 E 0~3 , and T 10 2.5 This problem is a benchmark commonly used due to the availability of a semianalytical solution [14]. The B 2 solution is computed using 16000 cells and the M 2 solution with 960 cells (as in [16]). The numerical results are presented in Figure  3(d). Here we see that the B 2 solution agrees with the benchmark solution values quantitatively somewhat better than the M 2 solution.

Conclusion.
We have considered the first member of the EQMOM hierarchy of moment methods for radiative transport proposed in [16], where ansätze for the specific intensity are constructed using the beta distribution. We showed that this method, which we call B 2 , is well-defined for all realizable moment vectors, has an explicit closure even on the boundary of realizability, and is hyperbolic. In these respects it is similar to the M n moment method, except that it is much less expensive to compute. We perform numerical simulations of standard benchmark problems to compare B 2 and M 2 (the corresponding member of the M n hierarchy) solutions and find good agreement between them.
For future work, the most important step is to extend this model beyond slab geometry to three-dimensional problems.