On the HolwayWeiss Debate: Convergence of the GradMomentExpansion in Kinetic Gas Theory
Abstract
Moment expansions are used as model reduction technique in kinetic gas theory to approximate the Boltzmann equation. Rarefied gas models based on socalled moment equations became increasingly popular recently. However, in a seminal paper by Holway [Phys. Fluids 7/6, (1965)] a fundamental restriction on the existence of the expansion was used to explain subshock behavior of shock profile solutions obtained by moment equations. Later, Weiss [Phys. Fluids 8/6, (1996)] argued that this restriction does not exist. We will revisit and discuss their findings and explain that both arguments have a correct and incorrect part. While a general convergence restriction for moment expansions does exist, it cannot be attributed to subshock solutions. We will also discuss the implications of the restriction and give some numerical evidence for our considerations.
0cm
1 Introduction
Rarefied gas dynamics is encountered in a number of modern technological fields such as highaltitude spacecrafts and microelectromechanical systems, and the modeling of rarefied gases has been of interest for more than one century. It is generally agreed that the Boltzmann equation, the fundamental equation in gas kinetic theory [8], provides an accurate description for rarefied gases in most applications. However, simulations using the Boltzmann equation require to solve a sixdimensional distribution function, which is computationally expensive and often unaffordable. Therefore, researchers have been trying to derive cheaper models from the Boltzmann equation by model reduction. A classical approach is the moment method, which was first introduced to gas kinetic theory by H. Grad in [11]. Grad suggested to expand the distribution function in velocity space using orthogonal polynomials with the local Maxwellian as the weight function, and derived a 13moment model by a truncation of such expansion.
In the past three decades many successful results in the modeling of rarefied gases could be obtained based on Grad’s approximation. Starting with the close relation to phenomenological extended thermodynamics (ET), Grad’s moment equations were used to compute sound waves, light scattering and shock profiles in the context of ET, see [19]. The regularized theory of [22] allows to formulate consistent boundary conditions [13, 27] based on Grad’s distribution. This leads to the successful application of Grad’s moment equations to a variety of boundary value problems [25]. Recently, Grad’s distribution has also been used with very high number of moments such that moment equations become a numerical technique to solve the Boltzmann equation efficiently and accurately [5, 26].
In the early days of Grad’s theory it was far from obvious that the theory was actually useful. In [11] H. Grad reports about artefacts in the shock profile computation obtained with his 13momentequations and concludes that the new equations do not show significant improvement over classical theories. The seminal paper of Holway [15] links the shock artefacts to a mathematical statement about the general convergence of moment expansions and concludes that moment equations cannot be used beyond a certain shock strength given by a critical Mach number. However, W. Weiss computed smooth shock profiles beyond the critical Mach number in [28] and subsequently published a paper [29] in which he disproves the statement of Holway about the usefulness of moment equations in shock waves.
In recent years computations of the authors of this paper made it clear that the HolwayWeissdebate deserves to be revisited, because the current representation in the literature is highly misleading. After introducing the foundations of Grad’s expansion in Sec. 2 we will carefully reframe the argument of Holway in Sec. 3 and discuss in detail in Sec. 3.2 what is right or wrong in both Holway’s statement and Weiss’ rebuttal. In Sec. 4 we discuss the consequences of the argument for nowadays moment models and give some numerical examples.
2 Basic Functional Analysis of the Grad Expansion
Given the density , velocity vector and temperature (in energy density units) of the gas the simplest approximation for the underlying distribution function is to assume equilibrium. The Maxwell distribution for the particle velocities
(1) 
then gives a distribution that fits the density, velocity and temperature. Here denotes the mass of a single gas molecule. For more general situations H. Grad proposed in 1949 to use an expansion which in the simplest form reads
(2) 
and became known as the Grad distribution. The indices represent multiindices such that are multivariate monomials up to degree and the coefficients parametrize the distribution. They are computed from the consistency requirement
(3) 
with given moments . Using this technique it is possible to reconstruct an approximation to the velocity distribution whenever a set of moments is known. For instance this approximative distribution can be used to close the transfer equations of the obtained from the Boltzmann equation which gives Grad’s moment equations.
2.1 Best Approximation
When using orthogonal polynomials instead of monomials the mathematical analysis of the Grad expansion becomes easier. In fact, the Grad distribution represents a best approximation in an appropriate subspace. We recall the following basic theorem, a variant of which could be found, e.g., in [9].
Theorem (Best Approximation).
Consider the weighted space with scalar product
(4) 
for functions on , and a set of orthonormal functions . For any function with we find the unique best approximation
(5) 
in the subspace spanned by .
In the Grad expansion the coefficients of are moments of which requires the choice with polynomial functions in order to find
(6) 
for the expansion. The natural choices for the inverse is a local Maxwellian with density , velocity , and temperature obtained from the respective local distribution . For more insight into how Grad’s expansion is used in modern gas modeling we refer to the textbook [21]. The polynomial functions satisfy the orthogonality condition
(7) 
such that they are nothing but Hermite polynomials with the Gaussian of the local equilibrium distribution as weight. The weighted functions are sometimes called Hermite functions.
The following theorem is the reason for a fundamental condition for the existence of a Grad expansion.
Theorem (Completeness).
The infinite set of Hermite functions form a complete orthogonal basis of the weighted space with .
This statement means any function with can be represented through a Grad expansion and vice versa. Hence, if is to be represented as an infinite Grad expansion, necessarily must hold.
2.2 Restrictions on the Distribution
The condition leads to
(8) 
hence the decay rate of must be faster than that of . In other words, if we let be the tail temperature of such that
(9) 
with some constants and , then for being the actual temperature of implied by ,
(10) 
must hold. In physical terms this means that the fast particles in the tail of the distribution should not be hotter than twice the average of all particles in the distribution.
2.3 Convergence Failure
However, (10) does not always hold. One typical example of distributions that may violate (10) is the superposition of two Gaussians:
(11) 
The velocity and temperature of such a distribution function can be calculated directly:
(12) 
Suppose . Then the tail of this distribution function is governed by Maxwellian with temperature . However, by (12), we can see that by decreasing , the value of can be set to be arbitrarily close to . This means that when , for any given and , we can always choose such that (10) does not hold, leading to the divergence of the Grad expansion.
Such an example is of interest since it comes from the MottSmith bimodal theory for the steady shock structure problem [18]. In such a theory, one considers a plane shock wave with Mach number , and it is assumed that the distribution function takes the form (11) everywhere, with the following dimensionless parameters:
(13) 
where varies with spatial location. To illustrate the divergence, we present in Fig. 1 the magnitude of the expansion coefficients for different and , from which one can observe that the coefficient does increase when is large and is small. In Fig. 2, we provide the comparison between the distribution function and the truncated Grad’s expansion. For and , Grad’s expansion converges, and the figure shows that the truncation at gives better approximation than . For and , Grad’s expansion diverges. Although the truncation at still seems to be approximating the original distribution function, the result given by clearly shows the failure of convergence.
3 A General Variant of Holway’s Argument
We will generalize the original argument of Holway to the case of a multidimensional domain in which we consider two arbitrary points connected by a straight line of length and direction unit vector as shown in Fig. 3. While Holway constructed his statement for the full Boltzmann collision operator, assuming some estimates for the gain part, we will directly use the BGK approximation to simplify the presentation. The extension to the Boltzmann operator is of purely technical nature and is not relevant for our discussion.
3.1 Derivation
The steady BGK equation is written with constant collision frequency for a particle velocity pointing from to
(14) 
which can be transformed into
(15) 
after multiplying with
(16)  
Replacing and we can compute the integral on the left hand side explicitly and find
(17) 
after multiplication with . Note that all terms in this equation are positive.
Holway then integrates from to with and writes
(18) 
to be used as lower bound of at . Using the mean value theorem with an , and we find
(19)  
(20) 
hence obtain
(21) 
with some positive constant independent of and for fixed . For close to either or the value of is very small and for at fixed collision frequency and length. If we assume an asymptotic decay of in the form
(22) 
the precise condition of Holway reads
(23) 
which implies the bound with the temperature at the position . This is because no matter how small becomes for large velocities, the exponentials are dominating and the inequality is only satisfied for all , if the exponential decays are behaving accordingly. Combining this with the decay restriction (10) of the Grad expansion at gives
(24) 
a condition involving the actual temperatures of the gas at two separated positions.
3.2 Holway’s Original Conclusions and Weiss’ Objection
In the original presentation of the argument in [15], Holway considered the onedimensional situation of a normal shock wave. Translated to our Fig. 3 he assumed a shock transition between and . He is using the coordinates and and considers particles moving in negative direction such that . His coordinate maybe identified with our position . Holway assumes that the shock is so thin that asymptotic shock conditions can be found already at and . Hence, the distributions at these points are given by the Maxwellians
(25)  
(26) 
where the parameters , and are defined in (13). Considering our setup in Fig. 3 this means that we find the equilibrium essentially in all positions , and , and at position . In particular Holway identified for the temperature and concluded from (24) that must hold in the boundary conditions for the applicability of Grad’s moment expansion.
The temperatures before and after the shock are connected by RankineHugoniot conditions depending on the shock’s Mach number. Assuming that belongs to the downstream or hot part after the shock and to the upstream or cold part before the shock, Holway computed a critical Mach number
(27) 
beyond which the temperatures before and after the shock would fail to satisfy condition (24). Actually, Holway [15] solved the equation incorrectly and got the limiting Mach number , as has been pointed out by Weiss in [29].
With this result Holway conjectured that when Grad’s moment expansion is truncated, the range of Mach numbers in which shock solutions exist is the interval , with being a number less than . Holway’s original phrasing was “When the expansion is truncated after a few terms, the region of convergence may be expected to be smaller than given by [the condition]”, where the condition refers to “” (as mentioned, that number is the result of a minor miscalculation). Although Holway did not explicitly explain what he meant by “the region of convergence”, he did provide a table showing “the ranges of convergence”, and the caption of the table is “range of Mach numbers for which continuous shock solutions exist”. Hence, he connected the convergence limit to the subshock artefacts of shock waves as reported by Grad [12].
In [28], Weiss found that the 21moment theory of extended thermodynamics [19] predicts smooth shock structures for any Mach number less than , which is in agreement with the theory of hyperbolic partial differential equations presented in [20]. However, this is beyond Holway’s proposed limit. This inspired Weiss to revisit Holway’s proof. Besides finding the calculational error, he also proposed another objection on Holway’s derivation. He claimed that Holway’s argument does not fix the direction of the shock wave and continues to derive the statement from the relation (21). Weiss concludes that this shows that one should set the boundary condition in the opposite way: should point to the hot fluid behind the shock, while should point to the cold fluid before the shock. In that case condition (24) is naturally satisfied, because in such a shock and the restriction on the Mach number is removed. In particular, Weiss concluded that Holway’s argument does not affect the existence of subshock artefacts in shock solutions of moment equations.
4 Further Discussion
4.1 What Holway and Weiss got right and what wrong
According to our derivation in Sec. 3.1, there is no restriction how to place a normal shock wave along the line from to . Weiss was right to point out that the shock direction could have been chosen differently from Holway. However, our multidimensional derivation also shows that the temperature condition (24) actually holds between the temperatures of any two position in any process.
Furthermore, in our view, one cannot conclude from (21) that . Note that the “constant” appearing in that statement actually depends on , that is, the scaled distance between and , with corresponding to the furthest away from . In fact, in the shock scenario of Holway we find
(28) 
hence, no contradiction to (21). In this sense, we support Holway’s argument on the convergence of Grad’s method, and in view of Sec. 3.1 we also conclude that Grad’s expansion might not converge in a process in which temperature ratios of more than a factor of 2 are present between any two points.
However, it remains unclear how this relates to the existence of smooth shock structure. Despite Holway’s conjecture, there is no clear evidence in his argument showing that a smooth shock structure does not exist for a Mach number larger than , especially for a lowmoment theory. Holway’s convergence argument can only explain the limiting diverging behavior of Grad’s theory, while it can be seen from Fig. 2 that sometimes an early truncation of Grad’s series can generate modest approximations of the distribution function. Although Weiss’ work [28] still focused only on Mach numbers less than , he has extended the same result to the 35moment theory, for which smooth shock structures exists for Mach number less than . The results are reported in [19, Chapter 12, Section 5]. In this sense, we support Weiss’ argument that the subshock problem cannot be related to the convergence restriction. Indeed, we tend to believe that lowmoment theories may still provide decent predictions for moderately rarefied gases, despite the possible divergence of Grad’s expansion in the limit of infinitely many moments.
Note, that the occurence of subshocks in shock solutions of moment equations is extensively explained, for example, in theoretical terms in the book [19] and on the basis of a model problem in [23]. The reason lies in the characteristics of the hyperbolic waves generated by the equation and their interaction with the relaxational part of the moment system. The highest characteristic velocity of the system gives the maximal speed with which infinitesimal disturbances can propagate. An inflow velocity exceeding this speed can only be sustained in a steady shock solution by introducing a discontinous subshock, which due to its nonlinearity may move faster than the characteristic limit.
4.2 Consequences for Moment Equations
Nowadays, a lot more numerical experiments have been done for Grad’s method with a large number of moments. To obtain results in strong nonequilibrium, another issue of Grad’s equations is the loss of hyperbolicity [19], that has to be fixed. Such an issue has been addressed systematically in [2, 3], which provides hyperbolic versions of Grad’s equations for any number of moments. However, such hyperbolicity fix does not change Grad’s ansatz, and therefore the convergence issue remains.
Based on the hyperbolic moment equations, several numerical experiments have been carried out in [6, 4, 7]. All the three works used a large number of moments in the simulation. However, after examining all the numerical tests in these three works, we find that for most cases, the maximum temperature ratio in the numerical solution is less than . There are only three exceptions, all found in [6]:

Shock tube problem with Knudsen number computed using 20/84 moments. See [6, Fig. 2].

Fourier flow with Knudsen number computed using 56 moments. See [6, Fig. 10(a)].

Fourier flow with Knudsen number computed using 56 moments. See [6, Fig. 10(b)].
In all the three cases, the Knudsen number is relatively small, so that the distribution function is close to the local equilibrium. Meanwhile, note that the numbers of moments presented in the above list correspond to the threedimensional case, meaning that an early truncation of Grad’s series is used equivalent to , and by our observation in Sec. 2.3, it can still be expected that one can get reasonable approximations of the distribution functions. Additionally, in all the above cases, the maximum temperature ratio does not exceed , so that the problem may still be manageable for a small number of moments. As a summary, existing numerical experiments support our conclusion that moderate nonequilibrium flow can still be well captured by Grad’s method with a moderate number of moments.
It is worth mentioning that the situation may change once Grad’s equations are linearized [24]. In the Grad expansion, the nonlinearity comes completely from the involvement of velocity and temperature in the distribution function. Therefore, the linearization of Grad’s equations changes the ansatz of the distribution function by replacing the local equilibrium by a global equilibrium:
(29) 
where
(30) 
with and being constants. Similar to (10), the convergence of the above expansion as requires that . For linearized Grad’s equations, this parameter is manually chosen by setting the global equilibrium about which the linearization is performed. Therefore, if is chosen sufficiently large such that throughout the computational domain, then the convergence can again be achieved regardless of whether the original Grad’s method converges. Such numerical results have appeared in a recent work [16] and [16, Figure 5.7] shows the result of the Fourier flow with Knudsen number and up to moments, where the temperature ratio is larger than . No divergence is observed since the value of is chosen to be even larger than the highest temperature in the numerical result. Nevertheless, such a method discards the consideration that the distribution functions are close to local equilibrium for dense gases, and therefore may lose efficiency of Grad’s method in the nearcontinuum regime. In fact, any equilibrium distribution with temperature different from will require a nontrivial expansion (29).
In his dissertation [14] Holway suggests a nonlinear alternative to Grad’s expansion specifically for the computation of normal shock profiles. He modifies the MottSmithansatz [18] of a bimodal distribution by replacing the Maxwellian with the lower temperature with a Grad expansion. Due to the superposition with the hotter Maxwellian a hot tail can be approximated at any point in the shock profile and convergence is recovered. Unfortunately, this approach hardly generalizes to other multidimensional processes in which the highest temperature is apriori not known.
4.3 Numerical Evidence
Since existing results do not explicitly show the failure of Grad’s method, we are going to support our argument by showing some results for a onedimensional problem () with boundary conditions. We assume that the fluid is located between two parallel diffusive walls with temperature and , and we want to solve for steady state. By approximating the collisions using the BGK operator, the governing equation can be written as
(31) 
with boundary conditions
(32)  
(33) 
where is the Knudsen number, and in the boundary conditions, and are chosen such that
(34) 
which ensures the mass conservation. Furthermore, we assume that the total mass equals :
(35) 
For this problem, we can find the exact solution in the limits and . When , the BGK equation shows that is independent of . Thus by the boundary condition, we see that
(36) 
The values of and can be solved from (35) and (34), and the result is
(37) 
The corresponding temperature is . Consider the distribution function on the boundaries. We see that a necessary condition for the convergence of the Grad expansion is and . This shows that if or , the convergence of the Grad expansion fails if is sufficiently large.
To study the limit of the solution as , we integrate the BGK equation with respect to , which gives us . By the boundary condition, we see that . Similarly, by multiplying the equation by and integrating with respect to , we know that is a constant. Now we use ChapmanEnskog expansion and assume
(38) 
By matching the terms, we get . We integrate the BGK equation against and approximate by , and the resulting equation, which is actually the Fourier equation, is
(39) 
Therefore is approximately a linear function. When , the distribution function tends to the local Maxwellian. Therefore by the boundary condition,
(40) 
Consequently,
(41) 
This shows that if or , we can find a sufficiently small such that the Grad expansion fails to converge.
To validate our argument that the Grad expansion may still work for small values of , in the numerical test, we choose the Knudsen number and , and we set the wall temperatures to be and , so that the Grad expansion will diverge for both small and large Knudsen numbers. In order to ensure that the number of boundary conditions matches the number of characteristics pointing insides the domain (which is necessary for the existence of the solution), we adopt the variation of Grad’s equations with global hyperbolicity [1], and we always choose an odd since the equations for even may have multiple exact solutions. For this onedimensional problem, the moment equations are solved with the shooting method with a quasiNewton method used for the iteration. The solutions for the temperature field are shown in Fig. 4 together with the result of the discrete velocity model (DVM) as reference solution.
From the DVM results, we see that for , the temperature on the right boundary exceeds twice the temperature on the left boundary. Therefore the Grad expansion is expected to be divergent. When , the temperature on the right boundary is less than , so that the Grad expansion also diverges. However, when we solve the moment equations, the quasiNewton iteration does not converge for some . For , the iteration converges only for and , while for , the iteration diverges for , but we can find solution for . For large , we cannot find solution for both Knudsen numbers.
From Fig. 4, one can also see that all the results are qualitatively correct. For smaller Knudsen number , in general, the result of gives better approximation. However, near the left boundary, where the Grad expansion is expected to be divergent, gives even larger error than . This can be observed more clearly from the approximation of the distribution functions, which are plotted in Fig. 5. It is shown that on the left boundary, the result of is closer to the DVM result, whereas on the right boundary, the result of is slightly better.
Similar plots for are provided in Fig. 6. In this example, due to the large discontinuity in the exact solution, it is harder for the moment method to get a good approximation. Nevertheless, when is not too large, the moment method still describes the general profile of the distribution functions.
5 Conclusion
We revisited a debate between L. H. Holway and W. Weiss on the convergence of moment approximations which keeps generating confusion in the community. In our view Holway was correct to show that there is a limit on the applicability of Grad’s moment expansion and in fact we generalized his argument to general multidimensional steady processes. However, Holway’s attempt to attribute the subshock behavior of moment equations when computing shock profiles to this convergence restriction is wrong. This misconception has been correctly pointed out by Weiss, whose smooth shock profile solutions based on moments remain valid. However, Weiss’ more substantial criticism of Holway’s argument turned out to be unfounded in our study.
Roughly speaking, Holway’s argument means that whenever there is a hot spot in a process, fast particles originating from that spot can be found anywhere in the domain, generating hot distribution tails that make Grad’s expansion diverge. While this bahavior is real and can be observed in specific computations, its implications for gas models based on a relatively small number of moments is probably negligible. We also discussed that even in simulations with many moments, convergence issues often remain undetected due to stabilizing effects like dissipation at moderate Knudsen numbers.
References
 [1] Z. Cai, Y. Fan, and R. Li, Globally hyperbolic regularization of Grad’s moment system in onedimensional space, Comm. Math. Sci., 11 (2013), pp. 547–571.
 [2] , Globally hyperbolic regularization of Grad’s moment system, Comm. Pure Appl. Math., 67 (2014), pp. 464–518.
 [3] , A framework on moment model reduction for kinetic equation, SIAM J. Appl. Math., 75 (2015), pp. 2001–2023.
 [4] Z. Cai, Y. Fan, R. Li, and Z. Qiao, Dimensionreduced hyperbolic moment method for the Boltzmann equation with BGKtype collision, Comm. Comput. Phys., 15 (2014), pp. 1368–1406.
 [5] Z. Cai and R. Li, Numerical regularized moment method of arbitrary order for BoltzmannBGK equation, SIAM J. Sci. Comput., 32 (2010), pp. 2875–2907.
 [6] Z. Cai, R. Li, and Z. Qiao, Globally hyperbolic regularized moment method with applications to microflow simulation, Comput. Fluids, 81 (2013), pp. 95–109.
 [7] Z. Cai and M. Torrilhon, Numerical simulation of microflows using moment methods with linearized collision operator, J. Sci. Comput., 74 (2018), pp. 336–374.
 [8] C. Cercignani, The Boltzmann Equation and its Applications, Applied Mathematical Sciences, Springer, New York, 1988.
 [9] E. W. Cheney, Introduction to Approximation Theory, Chelsea, New York, 1982.
 [10] W. Dreyer, Maximisation of the entropy in nonequilibrium, J. Phys. A, 20 (1987), pp. 6505–6517.
 [11] H. Grad, On the kinetic theory of rarefied gases, Comm. Pure Appl. Math., 2 (1949), pp. 331–407.
 [12] , The profile of a steady plane shock wave, Comm. Pure Appl. Math., 5 (1952), pp. 257–300.
 [13] X.J. Gu and D. Emerson, A computational strategy for the regularized 13 moment equations with enhanced wallboundary conditions, J. Comput. Phys., 225 (2007), pp. 263–283.
 [14] L. H. Holway, Approximation Procedures for Kinetic Theory, PhD Thesis, Harvard University, 1963.
 [15] , Existence of kinetic theory solutions to the shock structure problem, Phys. Fluids, 7 (1965), pp. 911–913.
 [16] Z. Hu, Z. Cai, and Y. Wang, Numerical simulation of microflows using Hermite spectral methods, 2019. arXiv:1807.06236.
 [17] C. Levermore, Moment closure hierarchies for kinetic theories, J. Stat. Phys., 83 (1996), pp. 1021–1065.
 [18] H. M. MottSmith, The solution of the boltzmann equation for a shock wave, Phys. Rev., 82 (1951), pp. 885–892.
 [19] I. Müller and T. Ruggeri, Rational Extended Thermodynamics, Second Edition, Springer, 1998.
 [20] T. Ruggeri, Breakdown of shock wave structure solutions, Phys. Rev. E, 47 (1993).
 [21] H. Struchtrup, Macroscopic Transport Equations for Rarefied Gas Flows, Interaction of Mechanics and Mathematics, Springer, New York, 2005.
 [22] H. Struchtrup and M. Torrilhon, Regularization of Grad’s 13 moment equations: Derivation and linear analysis, Phys. Fluids, 15 (2003), pp. 2668–2680.
 [23] M. Torrilhon, Characteristic waves and dissipation in the 13momentcase, Cont. Mech. Thermodyn., 12 (2000), p. 289.
 [24] , Convergence study of moment approximations for boundary value problems of the BoltzmannBGK equation, Commun. Comput. Phys., 18 (2015), pp. 529–557.
 [25] M. Torrilhon, Modeling nonequilibrium gas flow based on moment equations, Ann. Rev. Fluid Mech., 48 (2016), pp. 429–458.
 [26] M. Torrilhon and N. Sarna, Hierarchical boltzmann simulations and model error estimation, J. Comput. Phys., 342 (2017), pp. 66–84.
 [27] M. Torrilhon and H. Struchtrup, Boundary conditions for regularized 13momentequations for microchannelflows, J. Comput. Phys., 227 (2008), pp. 1982–2011.
 [28] W. Weiss, Continuous shock structure in extended thermodynamics, Phys. Rev. E, 52 (1995), pp. R5760–R5763.
 [29] , Comments on “Existence of kinetic theory solutions to the shock structure problem” [Phys. Fluids 7, 911 (1964)], Phys. Fluids, 8 (1996), pp. 1689–1690.