Direct Flux Gradient Approximation to Close Moment Model for Kinetic Equations
aa r X i v : . [ phy s i c s . c o m p - ph ] F e b Direct Flux Gradient Approximation to Close Moment Model forKinetic Equations
Ruo Li ∗ , Weiming Li † , Lingchao Zheng ‡ February 16, 2021
Abstract
To close the moment model deduced from kinetic equations, the canonical approach is to providean approximation to the flux function not able to be depicted by the moments in the reducedmodel. In this paper, we propose a brand new closure approach with remarkable advantages thanthe canonical approach. Instead of approximating the flux function, the new approach close themoment model by approximating the flux gradient. Precisely, we approximate the space derivativeof the distribution function by an ansatz which is a weighted polynomial, and the derivative of theclosing flux is computed by taking the moments of the ansatz. Consequently, the method providesus an improved framework to derive globally hyperbolic moment models, which preserve all thoseconservative variables in the low order moments. It is shown that the linearized system at the weightfunction, which is often the local equilibrium, of the moment model deduced by our new approach isautomatically coincided with the system deduced from the classical perturbation theory, which cannot be satisfied by previous hyperbolic regularization framework. Taking the Boltzmann equation asexample, the linearlization of the moment model gives the correct Navier-Stokes-Fourier law same asthat the Chapman-Enskog expansion gives. Most existing globally hyperbolic moment models arere-produced by our new approach, and several new models are proposed based on this framework.
Keywords:
Kinetic equation; moment model; global hyperbolicity; conservation law; Maxwellianiteration.
The kinetic equation is the evolution equation of particle distribution function in phase space, anddescribes the motion of particles and their interactions with each other or with the background medium.Common examples of the kinetic equation include the Boltzmann equation for rarefied gas, the radiativetransfer equation for photon transport, the vlasov equation for plasmas, etc. They have wide applicationsin fields like rarefied gases, microflow, semi-conductor device simulation, radiation astronomy, opticalimaging, and so on. Among the numerous methods developed to solve the kinetic equations [6, 4, 20, 23],moment methods are attractive due to their clear physical interpretation and high efficiency in thetransitional regimes [18, 31, 24, 30, 9, 22, 17].The moment method approximates the original kinetic equation by studying the evolution of a finitenumber of moments of the distribution function. In gas kinetic theory, it was first proposed in Grad’sseminal paper [18] in 1949, in which the notable Grad’s 13 moment system is also presented. For anymoment model, the evolution equations of the lower order moments rely on the higher order moments,hence a moment closure is needed, which is the central problem of the moment method. Different typesof moment closures have been developed, resulting in many different kinds of moment models. Thesemodels could roughly be divided into two types: the first type are called hyperbolic or first order PDEs,these include Grad’s 13 moment system [18], the maximum entropy moment system [24], approximatemaximum entropy moment systems [30], the P N and M N models for radiative transfer [35, 31], etc. Thesecond kind are parabolic type or second order PDEs, these include regularized moment methods ofvarious kinds [38, 36, 39] and the diffusion approximation for the radiative transfer equation. This paper ∗ CAPT, LMAM & School of Mathematical Sciences, Peking University, Beijing, China, email: [email protected] † Laboratory of Computational Physics, Institute of Applied Physics and Computational Mathematics, Beijing, China,email: [email protected] ‡ School of Mathematical Sciences, Peking University, Beijing, China, email: [email protected] M model, which is already globallyhyperbolic and based on which the MP N model is derived. On the other hand, in [13, 16], it was pointedout that the moment model derived by directly applying the framework in [11, 15] lead to incorrect NSFlaw and Eddington approximation, when a one-step Maxwellian iteration is applied. There has beensome recent progress in this direction. In [16], for the MP N model of the radiative transfer equation[17], a new hyperbolic regularization was proposed, in which the distribution function and its derivativewere approximated in different spaces. The HMP N model, based on this hyperbolic regularization, nolonger had the above disadvantages. This inspires us to propose a new hyperbolic regularization for themoment models of kinetic equation.In this paper, we propose a new approach to close hyperbolic moment models derived for kinetic equa-tions. We directly approximate the flux gradient to give the moment closure. Precisely, we approximatethe derivative of the distribution function by a weighted polynomial. In traditional moment methods,the moment closure is often given by approximating the distribution function, and some of them sufferfrom lack of hyperbolicity [10, 41, 17]. The new closure approach provide us an improved framework[11, 15] to carry out model reduction for generic kinetic equations. We first prove that using the improvedframework, the moment system is globally hyperbolic, as long as the weight function satisfies the basicrequirement of positivity. Moreover, for an N -th order moment system, moments with orders from 0 to N − P N model, the M N model, and the HMP N model for the radiative transferequation [35, 31, 16, 25] can be regarded as special cases derived by the new approach. Using the newclosure approach, we propose new globally hyperbolic moment systems for kinetic equations, such as HMP N model for monochromatic case and a moment model for the Boltzmann-Peierls equation. Theimproved framework is extended to 3-D case naturally, where the advantages discussed above are allpreserved. Using the improved framework, we derive a new hyperbolic moment model for the Grad’s13-moment system for quantum gas.The remaining parts of the paper is organized as follows. In Section 2, we introduce our new closureapproach in the 1-D case as the essential ingredients of the improved framework, and some propertiesare proved. In Section 3, we list some classical moment models and put them into the category of theimproved framework, while some new models are proposed. In Section 4, the framework is extended to3-D case, and the 3-D examples are shown in Section 5. The paper ends with a conclusion in Section 6.2 Model in one-dimensional case
In this section, we first introduce the moment method for kinetic equation and our new framework in1-D case, where the kinetic equation is often given by ∂f∂t + v ( ξ ) ∂f∂x = S ( f ) , t ∈ R + , x ∈ R , v ( ξ ) ∈ R , ξ ∈ G ⊂ R . (2.1) f ( t, x, ξ ) is the distribution function depending on time t , spatial variable x and velocity-related variable ξ . The k -th moment is defined as E k ( t, x ) := h τ ( ξ ) v k f i , where τ ( ξ ) >
0, and h·i := Z G · d ξ. An N -th order moment system can be written as ∂E k ∂t + ∂E k +1 ∂x = S k , ≤ k ≤ N, (2.2)where S k := h τ ( ξ ) v k Si . However, the moment system (2.2) is not closed, since there are N + 1 equations and N + 2 variables E , E , · · · , E N +1 . Thus a moment closure needs to be applied. A typical moment closure is to use afunction of E , E , · · · , E N instead of E N +1 in (2.2), which is given by E N +1 = E N +1 ( E , E , · · · , E N ) , (2.3)and a commonly used method to obtain (2.3) is to give an ansatz ˆ f , which satisfies h τ ( ξ ) v k ˆ f i = E k , ≤ k ≤ N. (2.4)Then E N +1 is taken as the ( N + 1)-th moment of ˆ f to finish the moment closure. Based on this idea,many practical moment models [18, 24, 41, 1, 17] for kinetic equations have been proposed.On the other hand, in order to make (2.2) closed, one only needs to give the moment closure on ∂E N +1 ∂x . This inspires us to take the moment closure as ∂E N +1 ∂x = ∂E N +1 ∂x (cid:18) E , E , · · · , E N ; ∂E ∂x , ∂E ∂x , · · · , ∂E N ∂x (cid:19) . (2.5)It is natural to expect a quasilinear system to be derived, that we require that ∂E N +1 ∂x relies linearlyon ∂E ∂x , · · · , ∂E N ∂x . Since we can use an ansatz ˆ f to approximate the distribution function f , we cansimilarly give an ansatz ˆ g on ∂f∂x , which satisfies h τ ( ξ ) v k ˆ g i = ∂E k ∂x , ≤ k ≤ N, (2.6)and the moment closure (2.5) is then given by ∂E N +1 ∂x = h τ ( ξ ) v N +1 ˆ g i . [18, 41, 17] suggest using a weighted polynomial to approximate the distribution function f to give themoment closure (2.3), which motivates us to use a weighted polynomial to approximate ∂f∂x to give themoment closure (2.5). Precisely, we take a positive weight function ω [ η ] , where η is a vector composedof some parameters depend on moments E , E , · · · , E N , and the ansatz is given byˆ g = ω [ η ] N X i =0 g i v i . (2.7)3enote E i = h τ ( ξ ) v i ω [ η ] i as the i -th moment of weight function ω [ η ] , and matrix D ∈ R ( N +1) × ( N +1) and K ∈ R ( N +1) × ( N +1) satisfy D i,j = E i + j , and K i,j = E i + j +1 , 0 ≤ i, j ≤ N , i.e. D = E E E · · · E N E E E · · · E N +1 E E E · · · E N +2 ... ... ... . . . ... E N E N +1 E N +2 · · · E N , K = E E E · · · E N +1 E E E · · · E N +2 E E E · · · E N +3 ... ... ... . . . ... E N +1 E N +2 E N +3 · · · E N +1 . (2.8)Since τ ( ξ ) > ω [ η ] ( ξ ) >
0, it is not difficult to prove that D is symmetric and positive definite, and K is symmetric. Furthermore, denote g = ( g , g , · · · , g N ) T ∈ R N +1 , then according to the constraints(2.6), we have D g = (cid:18) ∂E ∂x , ∂E ∂x , · · · , ∂E N ∂x (cid:19) T , (2.9)and the moment closure is given by (cid:18) ∂E ∂x , ∂E ∂x , · · · , ∂E N +1 ∂x (cid:19) T = K g = KD − (cid:18) ∂E ∂x , ∂E ∂x , · · · , ∂E N ∂x (cid:19) T . (2.10)Notice in (2.10), the gradient of the flux functions rely linearly on ∂E k ∂x , k = 0 , · · · , N , which satisfiesour previous requirement. Therefore, the moment system (2.2) can be written as ∂ E ∂t + KD − ∂ E ∂x = S , (2.11)where E = ( E , E , · · · , E N ) T ∈ R N +1 , and S = ( S , S , · · · , S N ) T ∈ R N +1 . The moment system (2.11) is globally hyperbolic, and the eigenvalues can be calculated. We first givethis theorem:
Theorem 1.
Moment system (2.11) is globally hyperbolic.Proof.
To obtain the hyperbolicity, one only needs to investigate whether KD − is real diagonalizable.Notice that D is symmetric and positive definite, and K is symmetric, thus KD − is real diagonalizable.Therefore, moment system (2.11) is globally hyperbolic. In order to investigate the eigenvalue of (2.11), we first introduce a series of monic orthogonal polynomials φ [ η ] i ( v ), which satisfy Z G τ ( ξ ) ω [ η ] ( ξ ) φ [ η ] i ( v ) φ [ η ] j ( v ) d ξ = 0 , when i = j. (2.12)The Gram-Schmidt procedure to obtain the orthogonal polynomial can be formulated as φ [ η ]0 = 1 , φ [ η ] i = v i − i − X j =0 K i,j K j,j φ [ η ] j , (2.13)where K i,j := h τ ( ξ ) ω [ η ] ( ξ ) v i φ [ η ] j i = Z G τ ( ξ ) ω [ η ] ( ξ ) v i φ [ η ] j d ξ. (2.14)According to the orthogonality (2.12), we have K i,i = h τ ( ξ ) ω [ η ] ( ξ )( φ [ η ] i ) i > , K k,i = 0 , when k < i. Furthermore, multiplying v k τ ( ξ ) ω [ η ] ( ξ ) by (2.13), and taking integration with respect to ξ on G , K k,i can be calculated by K k,i = E k + i − i − X j =0 K i,j K j,j K k,j . τ ( ξ ) > ω [ η ] ( ξ ) >
0, one can always define the orthogonal polynomials φ [ η ] i by (2.13).Then, on the eigenvalue of (2.11), we have the following theorem: Theorem 2.
The eigenvalues of the moment system (2.11) are the zeros of the ( N + 1) -th orthogonalpolynomial φ [ η ] N +1 ( v ) .Proof. Consider the characteristic polynomial of (2.11)Det( KD − − λ ) = 0 , which can be rewritten as Det( K − λ D ) = 0 . Since λ is a eigenvalue, there exists a vector r = ( r , r , . . . , r N ) T ∈ R N +1 , satisfying that ( K − λ D ) r = 0,i.e. N X i =0 ( K k,i − λD k,i ) r i = 0 , ≤ k ≤ N, which is N X i =0 ( E k + i +1 − λ E k + i ) r i = 0 = * v k ( v − λ ) N X i =0 r i v i τ ( ξ ) ω [ η ] ( ξ ) + , ≤ k ≤ N. Therefore, ( v − λ ) P Ni =0 r i v i is orthogonal to v k , with respect to τ ( ξ ) ω [ η ] ( ξ ), for any 0 ≤ k ≤ N . Notice( v − λ ) P Ni =0 r i v i is a polynomial of v whose degree is not greater than N + 1, thus we have r N = 0, and( v − λ ) N X i =0 r i v i = r N φ [ η ] N +1 . Therefore, λ is the root of the orthogonal polynomial φ [ η ] N +1 . The conventional moment system, derived by postulating an ansatz for the distribution function toobtain (2.3), can be formulated as a conservation law, without considering the right hand side. Due tothe conservation, it is a popular way to deduce a moment system. However, the conventional momentsystem often suffers from lack of hyperbolicity. As a globally hyperbolic moment system, in this section,we will show some properties of (2.11) in comparison with the conventional moment system.First, it is direct to show the conservation of the moments with orders from 0 to N −
1, i.e., we have
Theorem 3.
In the moment system (2.2) , the first N equations can be written into conservation form.Proof. The first N equations, corresponding to the momensts with orders from 0 to ( N − ∂E k ∂t + ∂E k +1 ∂x = S k , ≤ k ≤ N − , therefore, the governing equations of the moments with orders from 0 to ( N −
1) can be written intoconservation form.However, compared to the conventional moment system given by (2.3) and (2.4), the frameworkcannot ensure the conservation property of the entire system. Nevertheless, we have the followingtheorem:
Theorem 4.
For the one-dimensional kinetic equation (2.1) , if we use the ansatz ˆ f = h N X i =0 η i v i ! (2.15) to approximate the distribution function f and derive a moment system A . The resulting system isequivalent to the system (2.11) by taking the weight function ω [ η ] = h ′ N X i =0 η i v i ! , (2.16) which is denoted as A . Moreover, the system (2.11) in this situation is in conservation form. roof. We write the system derived by taking the ansatz (2.15) into quasilinear form, formulated as ∂ E ∂t + M ∂ E ∂x = S , (2.17)where the element in the i -th row and j -th column, which we hereafter will refer to as the ( i, j )-th elementof matrix M , is M i,j = ∂E i +1 ∂E j . Notice the moment closure E N +1 is given by E N +1 = h τ ( ξ ) v N +1 ˆ f i = * τ ( ξ ) v N +1 h N X i =0 η i v i !+ , (2.18)and the constraints provided by the first N moments E k = h τ ( ξ ) v k ˆ f i = * τ ( ξ ) v k h N X i =0 η i v i !+ , ≤ k ≤ N, (2.19)since τ ( ξ ) is independent with η i , thus when 0 ≤ i ≤ N, ≤ k ≤ N + 1, ∂E k ∂η i = ∂ D τ ( ξ ) v k h (cid:16)P Ni =0 η i v i (cid:17)E ∂η i = * τ ( ξ ) v k h ′ N X i =0 η i v i ! v i + = D τ ( ξ ) ω [ η ] v i + k E = E k + i . Therefore, ∂E k ∂η i is the ( k, i )-th element in D (2.8), 0 ≤ i, k ≤ N , and ∂η i ∂E k is the ( k, i )-th element of D − ,since D is symmetric and positive definite. Furthermore, ∂E i +1 ∂η j is ( i, j )-th element of K , 0 ≤ i, j ≤ N .Thus, ∂E i +1 ∂E j = N X k =0 ∂E i +1 ∂η k ∂η k ∂E j = N X k =0 K i,k D − k,j = ( KD − ) i,j . Therefore, M = KD − , and the system (2.17) is equivalent to (2.11).At last, since (2.17) can be written as * τ ( ξ ) v k ∂ ˆ f∂t + v ∂ ˆ f∂x !+ = (cid:10) τ ( ξ ) v k S (cid:11) , ≤ k ≤ N, (2.17) and (2.11) are in conservation form.Furthermore, when the number of parameters η is less than N , we have the following theorem Theorem 5.
For the one-dimensional kinetic equation (2.1) , if we use the weight function ˜ ω [ η ] = h n X i =0 η i v i ! , (2.20) and make a weighted polynomial ˆ f = P Ni =0 f i v i ˜ ω [ η ] to approximate the distribution function f and derivea moment system A , and denote the system (2.11) by taking the weight function ω [ η ] = h ′ n X i =0 η i v i ! , (2.21) as A . Then the results of one-step Maxwellian iteration of A and A are the same.Proof. We still use the notation in (2.8), and A can be written as ∂ E ∂t + KD − ∂ E ∂x = S . On the other hand, denote ˜ E n = h τ ( ξ ) v n ˜ ω [ η ] i , and ˜ D as ˜ D i,j = ˜ E i + j , and ˜ K as ˜ K i,j = ˜ E i + j +1 , 0 ≤ i, j ≤ N , then A can be written as ∂ E ∂t + ∂ ( ˜ K ˜ D − E ) ∂x = S . w = ( η , η ∗ ) ∈ R N +1 , where there is anone-to-one map between w and E . Additionally, we assume that when ˆ f = ˜ ω [ η ] , E k = ˜ E k , and η ∗ = 0.Without loss of generality, we assume that A := ∂ E ∂ w is invertible.Then A and A can be respectively written as A ∂ w ∂t + ∂ ( ˜ K ˜ D − E ) ∂ w ∂ w ∂x = S , A ∂ w ∂t + KD − A ∂ w ∂x = S , In order to prove that the one-step Maxwellian iteration of these two systems are equivalent, one onlyneeds to prove that when ˆ f = ˜ ω [ η ] , i.e. E k = ˜ E k for 0 ≤ k ≤ N ,( KD − A ) m,l (cid:12)(cid:12)(cid:12) η ∗ =0 = ∂ ( ˜ K ˜ D − E ) m ∂η l (cid:12)(cid:12)(cid:12) η ∗ =0 , ≤ m ≤ N, ≤ l ≤ n. (2.22)The right hand side can be calculated as ∂ ( P Ni,j =0 ˜ K m,i ˜ D − i,j E j ) ∂η l (cid:12)(cid:12)(cid:12) η ∗ =0 = ∂ [( P Ni,j =0 ˜ K m,i ˜ D − i,j E j ) (cid:12)(cid:12) η ∗ =0 ] ∂η l = ∂ ˜ K m, ∂η l = K m,l (2.23)where when ˆ f = ˜ ω [ η ] , E j = ˜ E j = ˜ D j, , thus N X j =0 ˜ D − i,j E j (cid:12)(cid:12)(cid:12) η ∗ =0 = ( , i = 0 , , i = 0 , and ∂ ˜ K i, ∂η l = K i,l , ∂ ˜ D i, ∂η l = D i,l , ≤ i ≤ N, ≤ l ≤ n. are applied. On the other hand, notice that A i,l | η ∗ =0 = ∂E i ∂η l (cid:12)(cid:12)(cid:12) η ∗ =0 = ∂ ˜ D i, ∂η l (cid:12)(cid:12)(cid:12) η ∗ =0 = D i,l | η ∗ =0 . Therefore, ( KD − A ) m,l (cid:12)(cid:12)(cid:12) η ∗ =0 = K m,l , ≤ m ≤ N, ≤ l ≤ n, (2.24)thus (2.22) holds. Remark 1.
The result of Maxwellian iteration for many kinetic equations and their moment modelsis important. For instance, the one-step Maxwellian iteration of Grad’s 13-moment equation leads toNavier-Stokes-Fourier law [37, 13]. However, as pointed out in [13, 16], a direct application of thehyperbolic regularization in [15, 11] may change the result of Maxwellian iteration. According to abovetheorem, the moment system (2.11) fixes this defect.
In this section, we first list some existing hyperbolic moment models, some of which can be regardedas special cases of our framework, i.e., they can be written as the form in (2.11) by taking a properweight function ω [ η ] . Then we will also propose some novel hyperbolic moment models based on our newframework.In order to put moment models into our framework, we need to prove they are equivalent to (2.11).Apart from Theorem 4, in order to judge whether two moment models are the same, we introduce thislemma. Lemma 6.
If two one-dimensional N -th order moment system with N + 1 equations, denoted as system A and A , satisfy these two conditions:1. For each system, the governing equations of the moments with orders from to ( N − can bewritten in conservation form. In other words, the first N equations are equivalent to ∂E k ∂t + ∂E k +1 ∂x = S k , ≤ k ≤ N − . . The characteristic polynomial of A is the same as the characteristic polynomial of A .Then these two moment models are equivalent.Proof. According to the first condition, the moment system can be written as these quasi-linear forms ∂ E ∂t + A ∂ E ∂x = S and ∂ E ∂t + A ∂ E ∂x = S , where E = ( E , E , · · · , E N ) T ∈ R N +1 , and S = ( S , S , · · · , S N ) T ∈ R N +1 . The coefficient matrix A and A can be written as A = · · ·
00 0 1 0 · · ·
00 0 0 1 · · · · · · c , c , c , c , · · · c ,N and A = · · ·
00 0 1 0 · · ·
00 0 0 1 · · · · · · c , c , c , c , · · · c ,N . Therefore, the characteristic polynomials are p ( λ ) = λ N +1 − N X i =0 c ,i λ i and p ( λ ) = λ N +1 − N X i =0 c ,i λ i . Noticing that p ( λ ) = p ( λ ), we have c ,i = c ,i , i = 0 , , · · · , N . Thus we have A = A , thus A and A are equivalent. This section considers the one-dimensional Boltzmann equation. The Boltzmann equation depicts themovement of the gas molecules from a statistical point of view [5, 37]. It reads ∂f∂t + ξ ∂f∂x = S ( f ) , (3.1)where f = f ( t, x, ξ ) is the distribution function, ξ ∈ R is the microscopic velocity, and the right-handside S is the collision term, depicting the interaction between gas particles.The thermodynamic equilibrium is f eq = ρ √ πθ exp (cid:18) − | ξ − u | θ (cid:19) , (3.2)where ρ , u and θ are macroscopic variables, satisfying that ρ = h f i , ρu = h ξf i and ρθ = h| ξ − u | f i .Furthermore, the commonly used definition of the moments are E k = h ξ k f i , k ≥ . Compared with (2.1), we have v ( ξ ) = ξ and τ ( ξ ) = 1. To derive a reduced model from (3.1), Grad [18] proposed the moment method, in which a weightedpolynomial is applied to approximate the distribution function, with the weight function set as f eq , i.e.˜ ω [ η ] = ρ √ πθ exp (cid:18) − θ ( ξ − u ) (cid:19) . (3.3)By taking η = (cid:18) ln (cid:18) ρ √ πθ (cid:19) − u θ , uθ , − θ (cid:19) T , The weight function (3.3) can be written as the form(2.20), with h ( ζ ) = exp( ζ ). 8owever, in [32, 10], it was pointed out the Grad’s moment method is not globally hyperbolic, whichlimits the applications of the Grad’s moment method. Therefore, the hyperbolic moment equation [8, 9]is proposed as the result of a globally hyperbolic regularization of Grad’s moment system.With respect to the weight function (3.3), the orthogonal polynomial is the generalized Hermitepolynomial, formulated asHe [ η ] k ( ξ ) = ( − k exp (cid:18) − ( ξ − u ) θ (cid:19) d k d ξ k exp (cid:18) − ( ξ − u ) θ (cid:19) , (3.4)it is not difficult to check that He [ η ] k ( ξ ) = θ − k/ He k (cid:18) ξ − u √ θ (cid:19) , where He k is the k -th order Hermite polynomial. Therefore, the zeros of He [ η ] k are u + c i √ θ , where c i are zeros of Hermite polynomial He k , 1 ≤ i ≤ k . In other words, when we take the weight function as(3.3) in our framework, the eigenvalues of the moment system (2.11) are the zeros of He [ η ] N +1 .On the other hand, according to the discussions in [9], the governing equation of the moments withorders from 0 to N − E k , k ≤ N −
1) in HME system can be written in conservation form, andthe eigenvalues of the HME system is the zeros of He [ η ] N +1 . According to lemma 6, HME system can beput into our framework, by taking the weight function as ω [ η ] = h ′ n X i =0 η i ξ i ! = exp n X i =0 η i ξ i ! = ρ √ πθ exp (cid:18) − θ ( ξ − u ) (cid:19) . The maximum entropy model [24] makes use of the following ansatz for the distribution function:ˆ f ( t, x, ξ ) = exp N X i =0 η i ( t, x ) ξ i ! , (3.5)which corresponds to h ( ζ ) = exp( ζ ) in Theorem 4. Therefore, according to Theorem 4, take the weightfunction as ω [ η ] = h ′ N X i =0 η i ξ i ! = exp N X i =0 η i ξ i ! , and we know the two moment models are the same. This section considers the radiative transfer equation (RTE) in slab geometry1 c ∂f∂t + ξ ∂f∂x = S ( f, T ) , (3.6)where c is the speed of light, ξ ∈ [ − ,
1] denotes the cosine of the angle between the photon velocitydirection and the x axis, and the right-hand side S [7, 28] depicts the interaction between photons andthe background medium which has material temperature T .Furthermore, the commonly used definition of the moments are E k = h ξ k f i , k ≥ . Compared with (2.1), we have v ( ξ ) = ξ and τ ( ξ ) = 1. P N model The P N model [21] suggests to approximate the distribution function by a polynomial, i.e.ˆ f ( t, x, ξ ) = N X i =0 η i ( t, x ) ξ i , (3.7)9hich corresponds to Theorem 4 with h ( ζ ) = ζ . Therefore, take the weight function as ω [ η ] = 1 , then according to Theorem 4, we know that the moment system derived by our framework is the sameas the P N model. M N model To derive the M N model [31, 24, 14] for gray RTE, one uses the following ansatz for the distributionfunction: ˆ f ( t, x, ξ ) = " N X i =0 η i ( t, x ) ξ i − , (3.8)which corresponds to h ( ζ ) = 1 ζ in Theorem 4.On the other hand, for the M N model for monochromatic radiative transfer, the ansatz to approximatethe distribution function is given byˆ f ( t, x, ξ ) = " exp N X i =0 η i ( t, x ) ξ i ! − − , (3.9)which corresponds to h ( ζ ) = 1exp( ζ ) − ω [ η ] = 1 (cid:16)P Ni =0 η i ξ i (cid:17) , and the weight function for the monochromatic RTE is ω [ η ] = exp (cid:16)P Ni =0 η i ξ i (cid:17)(cid:16) exp (cid:16)P Ni =0 η i ξ i (cid:17) − (cid:17) . Then the M N model for the grey RTE and the monochromatic RTE can also regarded as special casesof our new framework (2.11). HMP N model In [17], the researchers suggest to use a weighted polynomial to approximate the specific intensity f ,where the weight function is given by the ansatz of the M model in grey medium, i.e.ˆ f = 1( η + η ξ ) N X i =0 f i ξ i , where η , η are determined by the moments E and E , and η η ∈ ( − , MP N model is not globally hyperbolic for N ≥ N = 1, this system is different from the M model. Moreover, this system’s one-step Maxwellian changesthe result of the MP N model. Therefore, in [16] a new hyperbolic regularization was proposed to fixthese defects, and the resulting system is HMP N model. Now we claim that the HMP N model can alsobe regarded as an example of our new framework.Noticing that in MP N model, we use ˜ ω [ η ] = 1 / ( η + η ξ ) as the weight function to approximate thespecific intensity f . Thus it is natural to use its derivative to approximate ∂f∂x , i.e. the weight functionis taken as ω [ η ] = 1( η + η ξ ) . (3.10)10ctually, the weight function is the same as the weight function we used in the M N model in Subsec-tion 3.2.2, when N = 1.According to lemma 6, one only need to check1. According to [16], the governing equations of the moments with orders from 0 to N − HMP N model is the zeros of the ( N + 1)-thorthogonal polynomial with respect to the weight function (3.10).Therefore, we have the HMP N model [16] is also a special case of our framework, when the weightfunction is taken as (3.10). According to Theorem 4, we have that the moment system (2.11) is the sameas the M model. According to Theorem 5, the Maxwellian iteration of the HMP N model is the sameas the MP N model. HMP N model for monochromatic case Since the
HMP N model for grey medium leads to a globally hyperbolic models which also preserves nicephysical properties [16], it is natural to propose the HMP N model by taking the weight function as themonochromatic M model, i.e. ˜ ω [ η ] = 1exp( η + η ξ ) − , (3.11)which corresponds to (2.21) with h ( ζ ) = ζ ) − .Taking the weight function as ω [ η ] = exp( η + η ξ )(exp( η + η ξ ) − , (3.12)one can propose a globally hyperbolic moment system with the first ( N − N = 1, the resultingsystem is the M model. According to Theorem 5, the resulting system preserves the result of theMaxwellian iteration the MP N model for monochromatic case. The one-dimensional Boltzmann-Peierls equation [33], which characterizes phonon transport, reads ∂f∂t + v ( ξ ) ∂f∂x = S ( f ) , (3.13)where v ( ξ ) = d w d ξ , w ( ξ ) is the dispersion relationship. For instance, if we only consider harmonicinteraction between adjacent atoms [19], w ( ξ ) = 2 r Km (cid:12)(cid:12)(cid:12)(cid:12) sin ξa (cid:12)(cid:12)(cid:12)(cid:12) , (3.14)where ξ is the wave number of lattice vibration, while K , m and a are constants. Usually, ξ ∈ B =[ − π/a, π/a ) is from the first Brillouin zone. S ( f ) depicts the phonon scattering process.As a high order extension of the approach in [2], the moments are defined as E k = h wv k f i , i.e. τ ( ξ ) = w ( ξ ) in this case. E is proportional to the energy density, and E proportional to the heatflux. We consider only the single mode relaxation time approximation [34], then equilibrium distributionfunction is f eq = 1exp( ~ w ( ξ ) / ( k B T )) − , (3.15)where ~ is the reduced Planck constant, k B is the Boltzmann constant, and T is temperature defined byenergy conservation (cid:28) ~ w f eq − f ˜ τ (cid:29) = 0 , (3.16)11here ˜ τ is relaxation time.A natural idea is to use a weighted polynomial to approximate the distribution function f , with theweight function taken as the equilibrium f eq . Some recent works can be found in [2, 3].Using the conventional moment method, the weight function is taken as˜ ω [ η ] = 1exp ( η w ( ξ )) − , and the ansatz is given by ˆ f = ˜ ω [ η ] N X i =0 f i v i . The resulting system is then formulated as ∂ E ∂t + M ∂ E ∂x = S , (3.17)with M i,j = ∂E i +1 ∂E j , and E N +1 is given by E N +1 = h w ( ξ ) ˆ f v N +1 i . Under our new framework, the derivative ∂f∂x can be approximated by a weighted polynomial, with ω [ η ] = [exp( η w ( ξ )) − − exp( η w ( ξ )) w ( ξ ) . (3.18)Precisely, the ansatz is formulated as ∂f∂x ≈ ˆ g = ω [ η ] N X i =0 g i v i . (3.19)According to Theorem 1, the moment system is globally hyperbolic. On the Maxwellian iteration, noticethat ∂ ˜ E i ∂η = ∂ h w ( ξ )˜ ω [ η ] v i i ∂η = h w ( ξ ) ∂ ˜ ω [ η ] ∂η v i i = h w ( ξ ) ω [ η ] v i i = E i , and according to proof of Theorem 5, the moment system of the new framework has the same result ofMaxwellian iteration with the conventional moment system (3.17). This section considers the D-dimensional kinetic equation, which has the form ∂f∂t + v ( ξ ) · ∇ x f = S ( f ) , t ∈ R + , x ∈ R D , v ( ξ ) ∈ R D , ξ ∈ G ⊂ R D , (4.1)where the distribution function f = f ( t, x , ξ ), and v ( ξ ) ∈ R D is a function of ξ . The k -th moment isdefined by E k ( t, x ) := h τ ( ξ ) ψ k ( v ) f i , ≤ k ≤ M − , (4.2)where ψ k are chosen polynomials of v , M is the number of moments considered in the moment system,and τ ( ξ ) > ψ k could be different. For example, in Grad’s 13-moment model, D = 3, M = 13, τ ( ξ ) = 1, v ( ξ ) = ξ , and the sequence of ψ k is ψ = (1 , ξ , ξ , ξ , ξ , ξ ξ , ξ ξ , ξ , ξ ξ , ξ , k ξ k ξ , k ξ k ξ , k ξ k ξ ) T . (4.3)In the N -th order HME model [9], we have M = (cid:0) N + DD (cid:1) , τ ( ξ ) = 1, v ( ξ ) = ξ , and ψ N ( α ) = ξ α ,where ξ α = Q Dd =1 ξ α d d , and N ( · ) is a map from the multi-index group { α ∈ N D : | α | ≤ N } to the set { , , , · · · , M − } .As h τ ( ξ ) f i often corresponds to density, hereafter we assume that ψ = 1.12 moment system can be written as ∂E k ∂t + D X d =1 ∂ h τ ( ξ ) v d ψ k f i ∂x d = S k , ≤ k ≤ M − , (4.4)where S k := h τ ( ξ ) ψ k Si , It is obvious that the moment system (4.4) is not closed. Thus one has to seek a moment closure. Similaras the 1-D case, we define the moment closure ∂ h τ ( ξ ) v d ψ k f i ∂x d , d = 1 , , · · · , D , by imposing an ansatz oneach ∂f∂x d . Similar to 1-D case, we approximate ∂f∂x d by a weighted polynomial ˆ g d , d = 1 , , · · · , D . We take apositive weight function ω [ η ] , where η depends on E k , 0 ≤ k ≤ M −
1. Thus for different directions x d ,their weight function is the same. The ansatz is then given byˆ g d = ω [ η ] M − X i =0 g d,i ψ i , ≤ d ≤ D. (4.5)Denote E k = h τ ( ξ ) ψ k ω [ η ] i is the k -th moment of weight function ω [ η ] . Let matrix D ∈ R M × M and K d ∈ R M × M , satisfying that D i,j = h τ ( ξ ) ψ i ψ j ω [ η ] i , and K d,i,j = h τ ( ξ ) v d ψ i ψ j ω [ η ] i , 0 ≤ i, j ≤ M − τ ( ξ ) ω [ η ] >
0, we have D is symmetric and positive definite, and K d issymmetric for all d = 1 , · · · , D . Furthermore, for a given 1 ≤ d ≤ D , denote g d ∈ R M to satisfy that its k -th element is g d,k , 0 ≤ k ≤ M −
1, then according to the constraints h τ ( ξ ) ψ k ˆ g d i = ∂E k ∂x d , ≤ k ≤ M − , (4.6)we have D g d = (cid:18) ∂E ∂x d , ∂E ∂x d , · · · , ∂E M − ∂x d (cid:19) T , (4.7)and the moment closure is given by (cid:18) ∂ h τ ( ξ ) v d ψ f i ∂x d , · · · , ∂ h τ ( ξ ) v d ψ M − f i ∂x d (cid:19) T = K d g d = K d D − (cid:18) ∂E ∂x d , ∂E ∂x d , · · · , ∂E M − ∂x d (cid:19) T . (4.8)Therefore, the multi-dimensional moment system can be written as ∂ E ∂t + D X d =1 K d D − ∂ E ∂x d = S , (4.9)where E = ( E , E , · · · , E M − ) T ∈ R M , and S = ( S , S , · · · , S M − ) T ∈ R M . Theorem 7.
The moment system (4.9) is globally hyperbolic.Proof.
Notice that τ ( ξ ) ω [ η ] > D is symmetric and positive definite and K d is symmetricfor d = 1 , · · · , D , thus any linear combination D P d =1 n d K d D − is real diagonalizable. Therefore, momentsystem (4.9) is globally hyperbolic. 13 .3 Comparison with the conventional moment system It is easy to obtain this theorem
Theorem 8.
In the multi-dimensional moment system, the moment h τ ( ξ ) φf i satisfy an equation inconservation form if v d φ can be linearly expressed by ψ l , ≤ l ≤ M − , for any d = 1 , , · · · , D . Remark 2.
When the considered moments are the moments with orders from to N , i.e. ψ N ( α ) = v α ,then the moments with orders from to ( N − , h τ ( ξ ) ξ α f i , | α | ≤ N − are conservative. Remark 3.
In 13-moment system, i.e. ψ k is given by (4.3) , then h τ ( ξ ) f i , h τ ( ξ ) ξ f i , h τ ( ξ ) | ξ | f i , are conservative, which often correspond to density ρ , momentum ρ u , and second-order moment ρθ + ρ | u | . However, for the entire system, the conservation form is not always preserved. The following theoremdiscusses a special case.
Theorem 9.
For the multi-dimensional kinetic equation (4.1) , if we use the ansatz ˆ f = h M − X k =0 η k ψ k ! (4.10) to approximate the distribution function f and derive a moment system A . The resulting system isequivalent to the system (4.9) by taking the weight function ω [ η ] = h ′ M − X k =0 η k ψ k ! , (4.11) which is denoted as A . Moreover, the system (4.9) in this situation is in conservation form.Proof. We rewrite the moment system A into quasi-linear form, ∂ E ∂t + D X d =1 M d ∂ E ∂x d = S , (4.12)where the ( i, j )-th element of matrix M d is M d,i,j = ∂ h τ ( ξ ) v d ψ i ˆ f i ∂E j . Notice h τ ( ξ ) v d ψ i ˆ f i = * τ ( ξ ) v d ψ i h M − X k =0 η k ψ k !+ , ≤ i ≤ M − , (4.13)and the constraints provided by the moments satisfying α ∈ I are E i = h τ ( ξ ) ψ i ˆ f i = * τ ( ξ ) ψ i h M − X k =0 η k ψ k !+ , ≤ i ≤ M − , (4.14)we have that when d = 1 , , · · · , D and 0 ≤ i, j ≤ M − ∂E i ∂η j = ∂ D τ ( ξ ) ψ i h (cid:16)P M − k =0 η k ψ k (cid:17)E ∂η j = * τ ( ξ ) ψ i h ′ M − X k =0 η k ψ k ! ψ j + = D τ ( ξ ) ω [ η ] ψ i ψ j E = D i,j , and ∂ h τ ( ξ ) v d ψ i ˆ f i ∂η j = ∂ D τ ( ξ ) v d ψ i h (cid:16)P M − k =0 η k ψ k (cid:17)E ∂η j = * τ ( ξ ) v d ψ i h ′ M − X k =0 η k ψ k ! ψ j + = D τ ( ξ ) ω [ η ] v d ψ i ψ j E = K d,i,j . ∂η j ∂E i is the ( i, j )-th element of D − , since D is symmetric and positive definite. Furthermore,according to ∂ h τ ( ξ ) v d ψ i ˆ f i ∂E j = M − X k =0 ∂ h τ ( ξ ) v d ψ i ˆ f i ∂η k ∂η k ∂E j = M − X k =0 K d,i,k D − k,j = ( K d D − ) i,j . Therefore, M d = K d D − , for any d = 1 , , · · · , D , and the system (4.12) is equivalent to (4.9).At last, since (4.12) can be written as * τ ( ξ ) ψ k ∂ ˆ f∂t + v · ∇ x ˆ f !+ = h τ ( ξ ) ψ k Si , ≤ k ≤ M − , (4.12) and (4.9) are in conservation form.Next, we consider cases where the number of parameters η is less than the number of moments M .We calculate the result of one-step Maxwellian iteration for setting f (0) = ω [ η ] instead of f (0) = f eq .Similar as 1-D case, we have the following theorem Theorem 10.
For the multi-dimensional kinetic equation (4.1) , if we use the weight function ˜ ω [ η ] = h n X i =0 η i ψ i ! , (4.15) where n ≤ M − , and construct a weighted polynomial ˆ f = P M − i =0 f i ψ i ˜ ω [ η ] to approximate the distributionfunction f and derive a moment system A , and denote the system (4.9) by taking the weight function ω [ η ] = h ′ n X i =0 η i ψ i ! , (4.16) as A . Then the results of one-step Maxwellian iteration of A and A are the same.Proof. Denote ˜ D as ˜ D i,j = h τ ( ξ ) ψ i ψ j ˜ ω [ η ] i , and ˜ K d as ˜ K d,i,j = h τ ( ξ ) v d ψ i ψ j ˜ ω [ η ] i , 0 ≤ i, j ≤ M −
1, then A can be written as ∂ E ∂t + D X d =1 ∂ ( ˜ K d ˜ D − E ) ∂x d = S . Furthermore, we can rewrite these two systems by a variable w = ( η , η ∗ ) ∈ R M , where there is anone-to-one map between w and E . Additionally, we assume that when ˆ f = ˜ ω [ η ] , η ∗ = 0. Without lossof generality, we assume that A := ∂ E ∂ w is invertible.Then A and A can be respectively written as A ∂ w ∂t + D X d =1 ∂ ( ˜ K d ˜ D − E ) ∂ w ∂ w ∂x d = S , A ∂ w ∂t + D X d =1 K d D − A ∂ w ∂x d = S , In order to prove that the one-step Maxwellian iteration of these two systems are equivalent, one onlyneeds to prove that( K d D − A ) m,l (cid:12)(cid:12) η ∗ =0 = ∂ ( ˜ K d ˜ D − E ) m ∂η l (cid:12)(cid:12)(cid:12) η ∗ =0 , ≤ m ≤ M − , ≤ l ≤ n, ≤ d ≤ D. (4.17)The right hand side can be calculated as ∂ ( P M − i,j =0 ˜ K d,m,i ˜ D − i,j E j ) ∂η l (cid:12)(cid:12)(cid:12) η ∗ =0 = ∂ [( P M − i,j =0 ˜ K d,m,i ˜ D − i,j E j ) | η ∗ =0 ] ∂η l = ∂ ˜ K d,m, ∂η l (cid:12)(cid:12)(cid:12) η ∗ =0 = K d,m,l | η ∗ =0 , (4.18)where ψ = 1, thus ( M − X j =0 ˜ D − i,j E j ) | η ∗ =0 = ( , i = 0 , , i = 0 . ∂ ˜ K d,i, ∂η l = K d,i,l ∂ ˜ D i, ∂η l = D i,l , ≤ i ≤ M − , ≤ l ≤ n are applied. On the other hand, notice that A i,l | η ∗ =0 = ∂E i ∂η l (cid:12)(cid:12)(cid:12) η ∗ =0 = ∂ ˜ D i, ∂η l (cid:12)(cid:12)(cid:12) η ∗ =0 = D i,l | η ∗ =0 , ≤ i ≤ M − , ≤ l ≤ n. Therefore, ( K d D − A ) m,l (cid:12)(cid:12)(cid:12) η ∗ =0 = K d,m,l , ≤ m ≤ M − , ≤ l ≤ n, ≤ d ≤ D, (4.19)Thus (4.17) holds. In this section, we will list come existing hyperbolic moment models in multi-dimensional case. Theycan be regarded as special cases of our new framework. Furthermore, we will also propose some newhyperbolic models.
In Subsection 3.1, we consider Boltzmann equation in 1-D case and verify the equality of our frameworkand some existing hyperbolic moment models, such as one-dimensional HME and maximum entropymodel. In this section, we will prove that in multi-dimensional case, the HME and maximum entropymodel can still be included in our new framework (4.9).First, the multi-dimensional Boltzmann equation reads ∂f∂t + ξ · ∇ x f = S ( f ) , (5.1)where f = f ( t, x , ξ ) is the distribution function, with x ∈ R D , and ξ ∈ R D . Compared with (4.1), wehave v ( ξ ) = ξ .The collision term S ( f ) depicts the interactions between particles, whose form is given by [18, 37] S ( f ) = Z R D Z S D − ( f ′ f ′∗ − f f ∗ ) B ( | ξ − ξ ∗ | , σ ) d n d ξ ∗ , where f, f ′ , f ∗ , f ′∗ are the shorthand notations for f ( t, x , ξ ), f ( t, x , ξ ′ ), f ( t, x , ξ ∗ ) and f ( t, x , ξ ′∗ ), respec-tively. ( ξ , ξ ∗ ) and ( ξ ′ , ξ ′∗ ) are the velocities before and after collision, n represents the collision angle,and B ( | ξ − ξ ∗ | , σ ) is the collision kernel.The local equilibrium is given as Maxwellian [26], formulated as f eq = ρ √ πθ D exp (cid:18) − | ξ − u | θ (cid:19) , where ρ = h f i , ρ u = h ξ f i , ρ | u | + D ρθ = (cid:28) | ξ | f (cid:29) . (5.2) The Grad’s moment model [18] suggests to use a weighted polynomial to approximate the distributionfunction, where the weight function is taken as the local equilibrium,˜ ω [ η ] = ρ √ πθ D exp (cid:18) − | ξ − u | θ (cid:19) , By taking η = ln ρ √ πθ D ! − | u | θ , u θ , − θ ! T , one can rewrite the weight function into the form(4.15). 16he moments considered in N -th order HME model is E k , 0 ≤ k ≤ M −
1, where M = (cid:0) N + DD (cid:1) , and E N ( α ) = h ξ α f i , | α | ≤ N. Compared with (4.2), we have τ ( ξ ) = 1. Remark 4.
In this example, we require that N (cid:0)(cid:8) α ∈ N D : | α | ≤ N (cid:9)(cid:1) = { , , , · · · , M − } . For example, we can let [9] N ( α ) = D X i =1 (cid:18)P Dk = D − i +1 α k + i − i (cid:19) . In the following discussion of this paper, we will use ( · ) α to represent the N ( α ) of vector ( · ) , and ( · ) α,β to represent the ( N ( α ) , N ( β )) -th element of matrix ( · ) . The ansatz for f is ˆ f = X | α |≤ N f α H [ η ] α = ˜ ω [ η ] X | α |≤ N f α He [ η ] α , where H [ η ] α ( ξ ) = He [ η ] α ( ξ )˜ ω [ η ] ( ξ ) is orthogonal basis, and He [ η ] α ( ξ ) is generalized Hermite polynomial,which is orthogonal polynomials with respect to the weight function ˜ ω [ η ] . Actually, according to [12, 9],we have He [ η ] α ( ξ ) = D Y d =1 He [ η ] α d ( ξ d ) = D Y d =1 He α d (cid:18) ξ d − u d √ θ (cid:19) , where He [ η ] k and He k are defined in Subsection 3.1.1.With a little abuse of the notation, we assume He [ η ] α are unit orthogonal polynomials, and H [ η ] α isunit orthogonal basis. Property 11.
Multi-dimensional HME model is equivalent to the moment model (4.9) by taking theweight function as ω [ η ] = ρ √ πθ D exp (cid:18) − | ξ − u | θ (cid:19) . Proof.
According to the statements in [11], multi-dimensional HME system can be written as˜ D ∂ w ∂t + D X d =1 ˜ M d ˜ D − ∂ w ∂x d = ˜ S , (5.3)where ˜ S ∈ R M is a vector whose N ( α )-th element is˜ S α = h He [ η ] α Si = Z R D He [ η ] α S d ξ . Moreover, w ∈ R M is a vector corresponding to the ansatz of the distribution function ˆ f . Furthermore,we have [11] for | α | , | β | ≤ N and d = 1 , , · · · , D ,˜ D α,β = Z R D He [ η ] α ∂ ˆ f∂ w β d ξ , (5.4)and ˜ M d,α,β = Z R D ξ d He [ η ] α He [ η ] β ˜ ω [ η ] d ξ . (5.5)On the other hand, in our new framework, take the weight function as ω [ η ] = ρ √ πθ D exp (cid:18) − | ξ − u | θ (cid:19) , ψ N ( α ) = ξ α , | α | ≤ N , then the resulting system is (4.9). Now we prove that thesetwo system are the same. First, notice that He [ η ] α , | α | ≤ N and ψ N ( α ) = ξ α , | α | ≤ N are two bases of { ξ α : | α | ≤ N } , there exists a matrix T ∈ R M × M , which satisfies(He [ η ]0 , He [ η ] e , . . . , He [ η ] α ∗ ) T = T ( ξ , ξ e , . . . , ξ α ∗ ) T , where α ∗ is the last element in I = { α ∈ N D : | α | ≤ N } , i.e. N ( α ∗ ) = M − S and S , we have˜ S = T S . (5.6)Moreover, let matrix ˜ D ′ be defined as ˜ D ′ α,β = Z R D ξ α ∂ ˆ f∂ w β d ξ , (5.7)then we have ˜ D ′ ∂ w ∂s = ∂ E ∂s , s = t, x , x , · · · , x D , and ˜ D = T ˜ D ′ . Moreover,˜ M d = TK d T T , and D = T − T − T , where D and K d are as defined in (4.9). Therefore, we have˜ M d = TK d D − T − , and (5.3) can be rewritten as T ∂ E ∂t + D X d =1 TK d D − ∂ E ∂x d = T S , which is equivalent to (4.9). Therefore, multi-dimensional HME can be included by our framework. Grad’s 13-moment [18] is the most famous model in the Grad-type moment hierarchy, and the HMEmodel for 13-moment system is its hyperbolic regularized version [15]. Different from the Grad’s momentmethod mentioned in the previous section, the moments considered here do not correspond to h τ ( ξ ) ξ α f i ,with | α | ≤ N . In this section, we will introduce the HME model for 13-moment system and put it intoour new framework (4.9).In this case, the dimension D = 3, and the moments are defined as E = ( h f i , h ξ f i , h ξ f i , h ξ f i , h ξ f i , h ξ ξ f i , h ξ ξ f i , h ξ f i , h ξ ξ f i , h ξ f i , hk ξ k ξ f i , hk ξ k ξ f i , hk ξ k ξ f i ) T , which implies that τ ( ξ ) = 1, and ψ = (1 , ξ , ξ , ξ , ξ , ξ ξ , ξ ξ , ξ , ξ ξ , ξ , k ξ k ξ , k ξ k ξ , k ξ k ξ ) T .∂f∂x d is approximated by the ansatz ˆ g d , written asˆ g d = ω [ η ] M − X i =0 g i ψ [ η ] i . then the moment system can be written as ∂ E ∂t + X d =1 K d D − ∂ E ∂x d = S , (5.8)where the ( i, j )-th element of D and K d are D i,j = h ψ [ η ] i ψ [ η ] j ω [ η ] i , K d,i,j = h ξ d ψ [ η ] i ψ [ η ] j ω [ η ] i . Using the same technique as in the proof of Property 11, it is not difficult to see that (5.8) is equivalentto the HME13 system proposed in [15]. 18 .1.3 Maximum entropy model
The moments considered in the full moment maximum entropy model [24] for multi-dimensional case aredescribed by ψ N ( α ) = ξ α , | α | ≤ N .The maximum entropy model uses the following ansatz for the distribution function:ˆ f ( t, x , ξ ) = exp M − X k =0 η k ( t, x ) ψ k ! , (5.9)which corresponds to h ( ζ ) = exp( ζ ) in Theorem 9. Therefore, according to Theorem 9, take the weightfunction as ω [ η ] = h ′ M − X k =0 η k ψ k ! = exp M − X k =0 η k ψ k ! , and we know the two moment models are the same. Three-dimensional radiative transfer equation is written as1 c ∂I∂t + ξ · ∇ x I = S , (5.10)where the distribution function (specific intensity) I = I ( t, x , ξ ), and t ∈ R + , x ∈ R , ξ ∈ S . The righthand side S depicts the interactions with other photons and the background medium. In RTE, we have τ ( ξ ) = 1 and v ( ξ ) = ξ .In radiative transfer equation, notice that k ξ k = 1, we have that X d =1 h ξ α ξ d f i = h ξ α f i , α ∈ N . Thus the moments h ξ α f i , | α | ≤ N can be expressed by h ξ α f i , α ∈ I , with I = { α ∈ N : | α | ≤ N, α ≤ } , and M = I = ( N + 1) . A map from I to { , , , · · · , M − } can be defined as N ( α ) = ( α + α + α ) + ( α + α + α + 1) α + α , and the basis function is ψ N ( α ) = ξ α , α ∈ I . M N model To derive the M N model for gray RTE, one uses the following ansatz for the distribution function:ˆ f ( t, x , ξ ) = M − X k =0 η k ( t, x ) ψ k ! − , (5.11)which corresponds to h ( ζ ) = 1 ζ in Theorem 9.On the other hand, for the M N model for monochromatic radiative transfer, the ansatz to approximatethe distribution function is given byˆ f ( t, x , ξ ) = " exp M − X k =0 η k ( t, x ) ψ k ! − − , (5.12)which corresponds to h ( ζ ) = 1exp( ζ ) − ω [ η ] = 1 (cid:18) M − P k =0 η k ψ k (cid:19) , and the weight function for the monochromatic RTE is ω [ η ] = exp (cid:18) M − P k =0 η k ψ k (cid:19)(cid:18) exp (cid:18) M − P k =0 η k ψ k (cid:19) − (cid:19) . Then the M N model for the grey RTE and the monochromatic RTE can also regarded as examples ofour new framework. HMP N model In [25], the
HMP N model was extended to 3D case, and a globally hyperbolic moment model for radiativetransfer equation was carried out.By taking ω [ η ] = 1( η + η ξ + η ξ + η ξ ) , using the same technique in the proof of Property 11, we know (4.9) is equivalent to 3D HMP N modelin [25]. Furthermore, according to Theorem 8, the first ( N − HMP N model preserves the result of the one-step Maxwellian iteration of the 3D MP N model, whichwill be changed by a direct application of the hyperbolic regularization in [15, 11]. The quantum Boltzmann equation, as well as the well-known Uehling-Uhlenbeck equation [40], is writtenas ∂f∂t + ξ · ∇ x f = S ( f ) , where the collision term is defined as S ( f ) = Z R Z π Z π [(1 − θf )(1 − θf ∗ ) f ′ f ′∗ − (1 − θf ′ )(1 − θf ′∗ ) f f ∗ ] gσ sin χ d χ d ε d ξ ∗ . Here f, f ′ , f ∗ , f ′∗ are the shorthand notations for f ( t, x , ξ ), f ( t, x , ξ ′ ), f ( t, x , ξ ∗ ) and f ( t, x , ξ ′∗ ), respec-tively. ( ξ , ξ ∗ ) and ( ξ ′ , ξ ′∗ ) are the velocities before and after collision. ε is the scattering angle, χ isthe deflection angle, g = | ξ − ξ ∗ | , and σ is the differential cross section. θ = 1 , , − v ( ξ ) = ξ .The macroscopic variables can be defined as ρ = m ˆ h Z R f d ξ , ρ u = m ˆ h Z R ξ f d ξ , p = m h Z R | ξ − u | f d ξ , where m is the mass of the particle, ˆ h = h/m , and h is Planck’s constant.The thermodynamic equilibrium is f eq = 1 z − exp (cid:16) | ξ − u | RT (cid:17) + θ , (5.13)where z and RT is related to ρ and p as ρ = m ˆ h √ πRT Li , p = m ˆ h √ πRT RT Li , (5.14)20nd Li s := − θ Li s ( − θ z ) is the polylogarithm. For the special case θ = 0, let Li s = z .The considered moments in the quantum Grad’s 13-moment system [41] are E = ( h f i , h ξ f i , h ξ f i , h ξ f i , h ξ f i , h ξ ξ f i , h ξ ξ f i , h ξ f i , h ξ ξ f i , h ξ f i , hk ξ k ξ f i , hk ξ k ξ f i , hk ξ k ξ f i ) T , which implies that τ ( ξ ) = 1, and ψ = (1 , ξ , ξ , ξ , ξ , ξ ξ , ξ ξ , ξ , ξ ξ , ξ , k ξ k ξ , k ξ k ξ , k ξ k ξ ) T . On the other hand, in [41], the authors proposed an ansatz for the distribution f G = f eq X i =0 f i ψ i , where f i are coefficients. It was pointed out in [13] that quantum Grad’s 13-moment system is nothyperbolic and a direct application of the classical framework in [11, 15] will change the NSF law ofthe quantum Grad’s 13-moment system. Therefore, researchers in [13] proposed a new regularized 13-moment system, by splitting the expansion into the equilibrium and non-equilibrium part and thenapplying the framework in [11, 15]. However, this method can not be extended to higher order momentmodel, in which case the conservation of most moments can not be preserved. In this paper, based onthe method in Section 4, we propose a new hyperbolic regularization of the quantum Grad’s 13-momentsystem.The weight function is taken as ω [ η ] = z − exp (cid:16) | ξ − u | RT (cid:17)(cid:16) z − exp (cid:16) | ξ − u | RT (cid:17) + θ (cid:17) , a globally hyperbolic moment system is obtained. According to Theorem 8, we have the conservation of ρ , ρ u , and p . For higher order moment models, we can get that higher order moments satisfy equationsin conservation form, which can not be ensured by the moment system proposed in [13]. On the otherhand, according to Theorem 10, the NSF law of the new moment system is the same as the quantumGrad’s 13-moment system, since f eq and ω [ η ] satisfy the conditions in Theorem 10. We proposed an improved framework for globally hyperbolic model reduction of kinetic equations by anew closure approach. Different from previous works, most of which derive moment closure by specifyingan ansatz for the distribution function, our method directly gives the closing relationship of the fluxgradient by approximating the spatial derivative of the distribution function with weighted polynomials.Besides being globally hyperbolic, moment models derived in this way satisfy nice physical properties,such as except for the highest order moments, all the other moments satisfy equations in conservationform. Also, the results of the Maxwellian iteration after the first step remain unchanged, ensuring themodels satisfy important properties like the Navier-Stokes-Fourier law. As shown in previous sections,most of the globally hyperbolic moment models developed in the previous literature could be regarded asspecial cases of this new framework, and we also derived new moment models, such as for the monochro-matic radiative transfer equation, phonon Boltzmann equation and quantum gases, which are either thefirst moment models for such equations with global hyperbolicity, or improvements on previously devel-oped ones. This framework provides us with a new tool to derive in a routine way globally hyperbolicmoment models which preserve nice physical properties. Numerical schemes specially designed for thistype of moment model will be further investigated.
Acknowledgements
The work of R.L. and L.Z. is partially supported by Science Challenge Project, No. TZ2016002 and theNational Natural Science Foundation in China (Grant No. 11971041). The work of W.L. is partiallysupported by Science Challenge Project (NO. TZ2016002) and National Natural Science Foundation ofChina (12001051). 21 eferences [1] G. W. Alldredge, R. Li, and W. Li. Approximating the M method by the extended quadraturemethod of moments for radiative transfer in slab geometry. Kinetic & Related Models , 9(2), 2016.[2] Z. Banach and W. Larecki. Nine-moment phonon hydrodynamics based on the modified grad-typeapproach: formulation.
Journal of Physics A: Mathematical and General , 37(41):9805, 2004.[3] Z. Banach and W. Larecki. Nine-moment phonon hydrodynamics based on the modified grad-type approach: hyperbolicity of the one-dimensional flow.
Journal of Physics A: Mathematical andGeneral , 37(45):11053, 2004.[4] G. Bird.
Molecular Gas Dynamics and the Direct Simulation of Gas Flows . Oxford: ClarendonPress, 1994.[5] L. Boltzmann. Weitere studien ¨uber das w¨armegleichgewicht unter gas-molek¨ulen.
Wiener Berichte ,66:275–370, 1872.[6] J. E. Broadwell. Study of rarefied shear flow by the discrete velocity method.
Journal of FluidMechanics , 19(03):401–414, 1964.[7] T. A. Brunner. Forms of approximate radiation transport.
Tech. Rep SAND2002-1778 , 2002.[8] Z. Cai, Y. Fan, and R. Li. Globally hyperbolic regularization of Grad’s moment system in onedimensional space.
Comm. Math. Sci. , 11(2):547–571, 2013.[9] Z. Cai, Y. Fan, and R. Li. Globally hyperbolic regularization of Grad’s moment system.
Comm.Pure Appl. Math. , 67(3):464–518, 2014.[10] Z. Cai, Y. Fan, and R. Li. On hyperbolicity of 13-moment system.
Kinet. Relat. Mod. , 7(3):415–432,2014.[11] Z. Cai, Y. Fan, and R. Li. A framework on moment model reduction for kinetic equation.
SIAM J.Appl. Math. , 75(5):2001–2023, 2015.[12] Z. Cai and R. Li. Numerical regularized moment method of arbitrary order for Boltzmann-BGKequation.
SIAM J. Sci. Comput. , 32(5):2875–2907, 2010.[13] Y. Di, Y. Fan, and R. Li. 13-moment system with global hyperbolicity for quantum gas.
Journal ofStatistical Physics , 167(5):1280–1302, 2017.[14] B. Dubroca and J. Feugeas. Theoretical and numerical study on a moment closure hierarchy forthe radiative transfer equation.
Comptes Rendus de l’Academie des Sciences Series I Mathematics ,329(10):915–920, 1999.[15] Y. Fan, J. Koellermeier, J. Li, R. Li, and M. Torrilhon. Model reduction of kinetic equations byoperator projection.
Journal of Statistical Physics , 162(2):457–486, 2016.[16] Y. Fan, R. Li, and L. Zheng. A nonlinear hyperbolic model for radiative transfer equation in slabgeometry.
SIAM Journal on Applied Mathematics , 80(6):2388–2419, 2020.[17] Y. Fan, R. Li, and L. Zheng. A nonlinear moment model for radiative transfer equation in slabgeometry.
Journal of Computational Physics , 404:109128, 2020.[18] H. Grad. On the kinetic theory of rarefied gases.
Comm. Pure Appl. Math. , 2(4):331–407, 1949.[19] Y. Guo and M. Wang. Phonon hydrodynamics and its applications in nanoscale heat transport.
Physics Reports , 595:1–44, 2015.[20] C. K. Hayakawa, J. Spanier, and V. Venugopalan. Coupled forward-adjoint Monte Carlo simulationsof radiative transport for the study of optical probe design in heterogeneous tissues.
SIAM Journalon Applied Mathematics , 68(1):253–270, 2007.[21] J. H. Jeans. Stars, gaseous, radiative transfer of energy.
Monthly Notices of the Royal AstronomicalSociety , 78:28–36, 1917. 2222] J. Koellermeier, R. Schaerer, and M. Torrilhon. A framework for hyperbolic approximation of kineticequations using quadrature-based projection methods.
Kinet. Relat. Mod. , 7(3):531–549, 2014.[23] E. W. Larsen and J. E. Morel. Advances in discrete-ordinates methodology. In
Nuclear Computa-tional Science , pages 1–84. Springer, 2010.[24] C. D. Levermore. Moment closure hierarchies for kinetic theories.
Journal of Statistical Physics ,83(5-6):1021–1065, 1996.[25] R. Li, P. Song, and L. Zheng. A nonlinear moment model for radiative transfer equation. arXivpreprint arXiv:2005.13142 , 2020.[26] J. C. Maxwell. On stresses in rarefied gases arising from inequalities of temperature.
Proc. R. Soc.Lond. , 27(185–189):304–308, 1878.[27] J. Mc Donald and M. Torrilhon. Affordable robust moment closures for cfd based on the maximum-entropy hierarchy.
Journal of Computational Physics , 2013.[28] R. G. McClarren, T. M. Evans, R. B. Lowrie, and J. D. Densmore. Semi-implicit time integrationfor P N thermal radiative transfer. Journal of Computational Physics , 227(16):7561–7586, 2008.[29] J. McDonald, J. Sachdev, and C. Groth. Use of the gaussian moment closure for the modellingof continuum and micron-scale flows with moving boundaries. In
Computational Fluid Dynamics2006 , pages 783–788. Springer, 2009.[30] J. McDonald and M. Torrilhon. Affordable robust moment closures for CFD based on the maximum-entropy hierarchy.
J. Comput. Phys. , 251:500–523, 2013.[31] G. N. Minerbo. Maximum entropy eddington factors.
Journal of Quantitative Spectroscopy andRadiative Transfer , 20(6):541–545, 1978.[32] I. M¨uller and T. Ruggeri.
Rational Extended Thermodynamics, Second Edition , volume 37 of
Springer tracts in natural philosophy . Springer-Verlag, New York, 1998.[33] R. Peierls. On the kinetic theory of thermal conduction in crystals. In
Selected Scientific Papers OfSir Rudolf Peierls: (With Commentary) , pages 15–48. World Scientific, 1997.[34] J.-P. M. P´eraud and N. G. Hadjiconstantinou. Efficient simulation of multidimensional phonontransport using energy-based variance-reduced monte carlo formulations.
Physical Review B ,84(20):205331, 2011.[35] G. Pomraning.
The equations of radiation hydrodynamics . Pergamon Press, 1973.[36] H. Struchtrup. Stable transport equations for rarefied gases at high orders in the Knudsen number.
Phys. Fluids , 16(11):3921–3934, 2004.[37] H. Struchtrup.
Macroscopic Transport Equations for Rarefied Gas Flows: Approximation Methodsin Kinetic Theory . Springer, 2005.[38] H. Struchtrup and M. Torrilhon. Regularization of Grad’s 13 moment equations: Derivation andlinear analysis.
Phys. Fluids , 15(9):2668–2680, 2003.[39] M. Torrilhon and H. Struchtrup. Regularized 13-moment equations: shock structure calculationsand comparison to Burnett models.
J. Fluid Mech. , 513:171–198, 2004.[40] E. A. Uehling and G. Uhlenbeck. Transport phenomena in Einstein-Bose and Fermi-Dirac gases. i.
Physical Review , 43(7):552, 1933.[41] R. Yano. Semi-classical expansion of distribution function using modified Hermite polynomials forquantum gas.