SOME RECENT ADVANCES ON THE METHOD OF MOMENTS IN KINETIC THEORY

. This review presents some recent works on the construction of closure relations for moment systems extracted from a kinetic equation. A rough construction of those closures and the main properties of the resulting systems are described. Especially, based on the underlying kinetic equation, the main properties desired for such moment systems are the realizability, i


Introduction
Kinetic theory is widely used to describe, at mesoscopic level, clouds of objects interacting with each other in a medium.Numerical methods to solve those equations can be classified into two types: the direct methods that solve the kinetic equation directly, typically the discrete Monte Carlo methods and the discrete velocity methods; and the moment methods that are focused on in this note, which commonly preserve some structure of the underlying kinetic model.
The method of moments in kinetic theory can be viewed either as a Petrov-Galerkin discretization with respect to the kinetic variable, or it can be interpreted as a model reduction leading toward a macroscopic description of the flow as it describes the evolution of macroscopic quantities.
The first issue when constructing such a discretization or model consists in defining the set of solutions in a Galerkin framework or a closure relation in a macroscopic description.In practice, the equations resulting from a moment method are under-determined and one needs additional relations before solving them.
For a moment model, the closure is not unique, and it needs to be chosen carefully.This choice is generally driven by the kinetic properties that one aims at retrieving after the moment's extraction.Among those properties, we mainly focus on three of them which are most taken into account in those constructions: 1the existence of a positive underlying kinetic distribution, also called realizability, 2-the entropy decay which describes the trend of the solution toward an equilibrium, and 3-the hyperbolicity of the moment system.
The closure problem has been widely studied over the last decades, and only a non-exhaustive list of the recent constructions of closures is presented below.In the next section, the kinetic and models are presented in more detail together with those three properties and their applications.The following four sections describe families of closures from the literature and their specificities, namely the method constructed after Harold Grad's technique in Section 3, the ones based on quadrature techniques in Section 4, the ones constructed out of the study of the entropy decay of the kinetic solution in Section 5 and those that aim at approximating positive distribution function in Section 6. Section 7 gathers conclusive comments.

Kinetic equations and moment extraction
The generic kinetic equation and its properties is first presented, before its moment version is deduced.

Kinetic equation and properties
Consider the kinetic equation where the unknown f (t, x, v) is a density of particles depending on time t ∈ R + , position x ∈ R d , and velocity v ∈ R d , where the space dimension is generally d = 1, 2 or 3.This equation is assumed to satisfy three properties that one commonly aims at retrieving under a certain form after the moment extraction: Well-posedness and non-negativity: Together with appropriate initial and boundary conditions, this equation is assumed to possess a unique solution that is non-negative f (t, x, v) ≥ 0 and depends continuously on the data.
Furthermore, the entropy dissipation term D(f ) = 0 if and only if f is of the form of a Maxwellian, i.e. if where the mean velocity u ∈ R d , density ρ ∈ R * + and temperature They can be interpreted from the moments of f of order 0 to 2 where the power v i is understood here in a tensorial sense, i.e. v 0 = 1, v 1 = v, v 2 = vv T .Hyperbolicity: The left-hand side of ( 1) is an hyperbolic transport operator.The solution of such systems commonly takes the form of waves propagating in a medium, and even if the collision term Q impacts the shape of the solution, this property also holds for the moment equation at equilibrium, i.e. when approximating f by M and extracting its first three moments.

Construction and properties of the moment system
Such a system is obtained by integrating (1) against a vector of basis function b(v), commonly chosen to be a polynomial basis (e.g.b 3 ) T ).For the sake of clarity in the following, the maximum degree of these polynomials is denoted N while the number of basis functions is denoted r (e.g. in the examples of b above N, d, r = 2, 3, 5 and N, d, r = 2, 3, 10).This yields the system where the moment vector f , flux F and collision operator Q are The unknown in the moment system (6a) is chosen to be f and one needs to define F and Q as functions of f instead of f .This is the so-called closure problem.One common strategy to close this system consists in interpreting (6a) as one particular discretization of (1), and the flux and collision terms are defined by replacing f in (6b) by one chosen reconstruction f R satisfying The reconstruction f R is defined from integral constraints.The sense given to such integrals (L 1 functions, distributions or measures) may impact on the well-definition of the flux and the collision term and on the properties satisfied by the moment system.Especially, f R needs to be integrable against |v| N +1 for the flux F to be well-defined, and one needs to give a sense to Q(f R ) when f R is not a function but a distribution or a measure.
With such a kinetic closure of the moment system, the properties of the kinetic equation turn into the following ones at the moment level: Hyperbolicity: The left-hand side of (6a) is an hyperbolic operator, i.e. the flux of Jacobian has real eigenvalues and a complete set of eigenvectors, for all directions n ∈ S d−1 .Entropy: There exists a set of entropy-entropy flux pair (h, g) satisfying for all directions n ∈ S d−1 .Ideally, one would like this hyperbolic entropy-entropy flux to be defined out of the kinetic entropy η(f ).The entropy dissipation can be obtained from the symmetrizability of the system ( [26]) which consists in writing the ansatz f as a function of parameters α and requiring that J α f (α) is symmetric positive definite and J α (F(α)n) is symmetric.Well-posedness and realizability: The moment solution is assumed to be unique, depending continuously on the data, and belongs to the following convex cone, so-called realizability domain where the integral can be understood in the sense of functions ).These three properties are rarely all satisfied by one choice of reconstruction f R , and compromises need to be made.The following sections aim to describe the axes of research that led to recent developments of closures and the properties they provide.In addition to those, one also needs to include in the balance: Numerical cost: The reconstruction, or directly the closure itself, needs to be numerically tractable and the numerical cost to obtain it needs to be put in balance with the brute force kinetic solvers.The first model emerging from the moment extraction is the system of Euler obtained by choosing for reconstruction f R = M the Maxwellian that minimizes Boltzmann entropy defined in (3) under the constraint of satisfying (4).This choice of reconstruction also provides a pressure law (here the closure) in these equations, and it leads naturally to a hyperbolic, entropy dissipated moment system with a realizable solution.All the next sections describe alternative closures or extensions of such a model when considering moments of higher order.

Necessity of these properties?
The method of moments applies to various kinetic equations, which possess similar features as (1).However, depending on the considered physics, some of the aforementioned properties are more important to preserve than others.None of the closures satisfy all the properties presented above, choices were made adapted to the considered physics or application.Therefore, we discuss the specificity of those properties and when they can be relaxed.

Positivity
All the density functions are positive in all the considered applications.But preserving this property and even identifying realizable vectors is difficult (see e.g.[29,42]).So it is rarely an objective on its own, and positivity is often enforced in order to obtain other properties.Typically, realizable closures are able to capture singular solutions (Dirac measures), while non-realizable ones are commonly simpler to compute in such singular regimes.Such regimes often appear in spray models or in radiative transfer and those physics often exploit realizable closures, while only relatively small perturbations away from the equilibrium are sufficient in rarefied gases and for some applications in plasma physics.

Entropy decay
This property is closely related to 1-the considered collision operator and 2-some equilibrium regime the solution tends to.For the first, some physics are modeled without collisions, and those models generally satisfy some other form of stability, which can be more meaningful to preserve when constructing closure relations.This can be the case for collisionless plasma or for the population balance equation, where the interactions between particles emerge from other effects than collisions.For the second, equilibrium distributions are often exploited, when available, when deriving moment closures and such a distribution needs to be accurately approximated by the considered reconstruction.

Hyperbolicity
Two types of hyperbolicity loss can be observed when constructing moment closures: First, the eigenvalues of the Jacobian of the flux can turn complex.This generally leads to non-physical oscillations, which need to be avoided.However, such a loss of hyperbolicity may appear in regimes that are not reached during a simulation.Second, the flux Jacobian may lack eigenvectors.Such weakly hyperbolic systems may create singular solutions, so-called delta-shocks, propagated with the flow (see e.g.[20,30]), which may be acceptable depending on the considered physics.Typically, in population balance or in radiative transfer, such solutions may occur.In such cases, the moments are understood in the sense of measures or of distributions.

Applications and their specificities
Here, we discuss the importance of the properties of the last subsection depending on the considered physical application and other specific properties.

Rarefied gases
The field of rarefied gases is probably the most classical application of kinetic theory.This is typically modeled exactly by (1) where the collision operator needs to be specified.
The study of moment models in this field is closely related to the construction of BGK relaxation operators.Indeed, relaxation operators are designed to be as simple as possible, exploiting only moment data and preserving the main properties of Boltzmann's term.Mainly, the entropy decay which drives the trend toward thermodynamical equilibrium is carefully studied for this purpose.

Plasma physics
The same kinetic equation is often used in plasma physics to model one or several populations of particles.These equations are generally coupled with an equation or a system for the electromagnetic fields.
Convergence to equilibrium is not only driven by the collisions (some models are collisionless), but also by electromagnetic effects.There are numerous plasma regimes and instabilities to capture described in the literature and as many models and collision operators.

Radiative transfer
Commonly, radiations are assumed to propagate exactly at the speed of light.Therefore, ( 1) is reduced into a model of the form where I is the radiation intensity, c is the speed of light, the normalized velocity Ω ∈ S 2 evolves on the unit sphere.Furthermore, the radiations interact only with the surrounding matter and not with each others.This leads to considering the collision operator L to be linear.These modifications have much impact on the construction of the moment model.On the one hand, the set of integration in the definition of the moments (5) becomes compact (it turns from R 3 in (1) into S 2 in ( 10)) which ensures that the realizability domain is open (see e.g.[37]).On the other hand, the collision operator LI (and the full kinetic equation) is linear, and any convex function η is an entropy dissipated through the equation, e.g.η(I) = I 2 /2 satisfies S 2 η ′ (I)LI = S 2 ILI ≤ 0 (see e.g.[2,24]).

Population balance equation
The population balance equation generically refers to the modeling of other clouds of objects interacting with each others.This includes e.g.liquid droplets or dust transported in another fluid or the motion of individuals in biological models.
The kinetic variable may include a velocity variable, as in (1), but also other state variables characterizing the object.These models also lead to various interaction operators (though commonly not collisions).Therefore, the impacts at the moment level are twofold: first, the moment extraction is performed with respect to all the kinetic variables and the set of integration in (2) is larger; second, the model (1) itself is modified, together with its interaction operator, and especially the properties that need to preserved at the moment level are various.

Grad's methods
Harold Grad was one of the first to propose an enriched moment closure to model out of equilibrium regimes.It is essentially applied for the construction of fluid models, typically rarefied gases.Though some recent works aim at extending the techniques and applying them in plasma physics ( [7]) or radiative transfer ( [16]).
Grad's work originally applied to rarefied gases.It aimed at capturing well the equilibrium represented by a Maxwellian and the small perturbations away from it, but it is not valid close to the boundary of the realizability domain.In this part, the moments are understood in the sense of L 1 functions.

Grad's closure (1949 ; [22])
The construction of Grad is a polynomial perturbation of the Maxwellian (3) where p is a polynomial of degree N .Its coefficients are typically computed using the change of unknown w : v → (v − u)/ √ T (with inverse v : w → √ T w + u), and exploiting expansions in some orthogonal polynomial basis, either Hermite (orthogonal with respect to the weight exp(−v 2 /2) on R) or Laguerre (with respect to exp(−v) on R + ).The reader is referred to [9,19,28] and references therein for computational details.This reconstruction takes the form of a (relatively) simple formula with analytic coefficients, so the numerical cost is very small.However, the other properties are not satisfied.Especially, f R is not positive on all the space v ∈ R 3 .A more problematic drawback is the appearance of non-hyperbolicity regions in the set of moments R that can be reached by this model.The appearance of such values during a simulation can be the source of non-physical artifacts.Grad also exhibited some perturbed entropy equality for this model.
In practice, this model is an extension of the Euler system that requires a low numerical cost to compute, and it is hyperbolic and entropy dissipating as long as the solution does not reach these pathologic regions.

Grad's regularizations
Several works aimed at extending the region of applicability of Grad's closure, mostly by compensating for the lack of hyperbolicity at the moment level.

Based on a Chapman-Enskog expansion
A first idea of regularization consisted in computing deviations from Grad's moment closure, by exploiting a Chapman-Enskog expansion near the equilibrium (see [44] and references therein).This expansion of the solution f in order of the Knudsen number commonly leads to diffusive (Navier-Stokes) or dispersive (Burnett or super-Burnett) systems.It is used here to correct the higher order moments away from the equilibrium and results in adding diffusion in the moment system.Such a regularization turns the model "valid" in a larger range of the attainable moment vectors, but the set of moment vectors satisfying positivity, hyperbolicity (for the flux part) and entropy decay remains a subset of the realizability domain.

Based on the Jacobian of the flux
A second idea based on the study of the Jacobian of the flux F(f ) aimed at modifying the structure of the Jacobian.The system is written in a quasi-linear form and one aims at modifying the matrix A(f ).This idea has been fruitful in the literature, and only a few of these recent strategies are listed below: • A first modification (see [10]) consisted in adding a regularization term, an appropriate flux term, to correct the characteristic speed.• A second modification ( [27]) was justified by the use of an appropriate quadrature method together with the expansion of the form (11). • A third approach ( [11]) consisted in performing the expansion (11) separately to each of the terms in (6a).That corresponded to taking into account a contribution of higher degree into the flux.
Eventually, they provide globally hyperbolic moment systems based on a known approximation of the distribution function (though still not positive).However, the resulting moment system a priori remains in a non-conservative formalism.

Through projection techniques
Other interpretations of the third idea of the last section were afterward provided and extended.Grad's method is interpreted here as a polynomial expansion using orthogonality w.r.t. to the Gaussian measure.Then, [15] suggested exploiting projections onto a set of weighted polynomials.These projections P are applied to the kinetic solution Pf , on its derivatives P∂ t Pf and P∇ x Pf and onto the flux term Pv T ∇ x Pf .Using an appropriate combination of those leads to a closed reduced system from (1) which is globally hyperbolic and symmetrizable.The system is eventually written in a non-conservative form and it is a modification of (2) in the sense that it is not obtained directly by imposing the closure terms.
This work was also extended to radiative transfer (see e.g.[16] and references therein).

Quadrature methods
As the moments are weighted integrals of a distribution function, a natural approximation consists in using appropriate quadrature formulae.
4.1.Quadrature Method Of Moments (QMOM ; 1984) The QMOM closure was indeed based on this simple remark.It originates in the field of spray modeling ( [34]).The original method consists simply in approaching the integrals in the definition of the vector of moments f using quadrature.One therefore needs to define the positive weights and abscissa of this quadrature, then the same quadrature is used to approach the integral in the definition of the flux function F(f ).This is equivalent to choosing a reconstruction in the form of a sum of Dirac (measures or distributions) A wide variety of algorithms is available in the literature to compute such quadratures, and the computational cost is rather low to compute them.
In terms of modeling properties, the resulting moment system is only weakly hyperbolic as the Jacobian of the flux possesses n eigenvalues (these are the quadrature points v i ) with a double multiplicity and only one eigenvector per eigenvalue.Such weakly hyperbolic systems tend to enhance δ-shock solutions, i.e. each of the δ in this reconstruction is propagated at its speed v i .
Similarly, entropies are dissipated through this equation (any ), but this is not directly related to the one dissipated at the kinetic level.

Multi-Gaussian or extended QMOM
The extended QMOM ( [48]) or multi-Gaussian ( [12,45]) closure consists in replacing the Diracs by regular functions in the quadrature formula.This is typically done using Gaussians with the same dispersion.The purpose of this modification is mostly to extend this model to a wider range of physical phenomena.For instance, non-zero ansatz is required in some specific values when considering evaporation of droplets in a spray, while Dirac deltas are zero everywhere but at the quadrature points.Also, the dispersion of the Gaussian leads to a model with a non-zero temperature.
The definition of dispersion is possible using a clever change of variables.However, using the same dispersion for all the Gaussians is necessary to perform the computations.
Eventually, the moment system is based on a strictly positive ansatz.At this step, few things can be said about hyperbolicity and entropy dissipation.

Hyperbolic QMOM
A recent alternative aimed at turning the QMOM model strongly hyperbolic.This technique is based on the algorithm of Chebychev ( [13,47]).The characteristic polynomial of the Jacobian of the flux can be obtained at the last step of this sequence of polynomials.All the coefficients of the QMOM closure can also be determined from the coefficients of this polynomial sequence.The HyQMOM closure consists in following this sequence up to one higher degree, and fixing the coefficients at the last iteration such that the Jacobian J f F is diagonalizable.A reconstruction with N + 1 positive Dirac measures can be deduced out of this technique.

Entropy-based methods
In a similar manner as the Maxwellian minimizes the entropy D under the constraints of satisfying (6b), the entropy-based methods exploit the minima of an entropy function.

Entropy-minimizing closure
The idea of minimizing an entropy in order to construct a distribution was studied in various fields (typically in information theory with the work of Shannon), and is often associated to Dave Levermore in kinetic theory ( [31]).This reconstruction is defined as This choice offers most of the properties expected in Section 2.2.For most physically relevant η, the reconstruction f R > 0 is strictly positive, the resulting moment model is symmetric hyperbolic ( [18,21]) and dissipates the underlying kinetic entropy H(f R ) (see e.g.[26]).Euler's system can be interpreted as the first model in this hierarchy.However, two drawbacks of these models have been widely studied in the literature: The first one, closer to the application, is the relatively high computational cost to compute this closure.In order to compute this closure, one needs to solve the optimization problem (12) at every location at every time.Algorithms were developed for this specific minimization problem, to reduce the numerical costs ( [5,23]), especially close to the boundary of the realizability domain where the condition of the optimization problem deteriorates.
The second problem appears only when the set of integration is unbounded (e.g. in rarefied gases v ∈ R 3 , but not in radiative transfer Ω ∈ S 2 for instance).In such cases, Michael Junk ( [25]) showed that there exists realizable vectors f ∈ R that do not possess such a reconstruction.This may lead to ill-posed problems.

Reducing computational costs
Several works aimed at reducing the computational costs of the entropy-based closure.A first approach ( [6]) consisted in changing basis functions when using the minimization algorithm in order to solve better conditioned problems.However, the optimization problem remains singular along the boundary of the realizability domain, where the distribution function turns singular.Another approach ( [4]) based on a Tikhonov regularization of the optimization problem leads to a well-defined solution along this boundary.

Relaxed moment constraints
A first idea to circumvent the appearance of such a Junk line consisted in relaxing the last moment constraint into an inequality ( [43]).For all moments realized by a minimum-entropy reconstruction, this reconstruction is identical.For all moment vectors along a Junk line, this technique provides a reconstruction with a lower (finite) entropy, but that does not realize the higher order moment.

φ-divergence closure
The φ-divergence closure, developed through the thesis of M. Abdelmalik ( [1,3]), can be interpreted as a modification of the entropy to minimize.In this framework, the equilibrium represented by the Maxwellian M does not depend on the perturbation, and needs to be fixed a priori.Once the Maxwellian M is fixed, the minimizer of an approximation of the relative entropy is computed.This provides a symmetric hyperbolic model for the perturbation away from the equilibrium, then it dissipates some approximation of the relative entropy.With an appropriate choice of approximation of the relative entropy, this yields a positive reconstruction f R , and no Junk line can be created through this closure.

Realizability-based methods
Some closures were constructed based on the study of the realizability domain.In the monovariate case v ∈ [−1, +1], the reconstruction is known to be unique along the boundary of the realizability domain.With the knowledge of this unique representation, one can then decompose any realizable vector as a positive sum of vectors which all have a known reconstruction.

Interpolative closures
The interpolative closures were developed in the fields of CFD ( [33]) and of radiative transfer ( [32,39,41]).It aimed at reducing the computational costs of the entropy-based closure by using a high order polynomial approximation that interpolates the exact entropy-minimizing closure at specific locations in the realizability domain.From a theoretical point of view, none of the properties of realizability, entropy dissipation and hyperbolicity can be ensured globally in the realizability domain, but these closures approximate a closure that does satisfy all of them (except for [32] that lose these properties in certain regions of R) and numerical experiments showed no loss of these properties.

Projective closure
The projective closure ( [37,38]) consists in a reconstruction of the form f R = α 0 f eq + f QM OM , where f eq corresponds to a chosen equilibrium function, typically a Maxwellian M , and f QM OM is a quadraturebased reconstruction.This can be interpreted using projection techniques on the boundary of the realizability domain.
This closure ensures realizability and weak hyperbolicity.However, as for the φ-divergence technique, the numerical computation of this closure when f eq depends on f remains difficult and the entropic structure in that case is weaker than the symmetric hyperbolic structure of the entropy-minimizing closure.

Conclusive remarks
The construction of appropriate closure for moment equations remains an active field of research for accurately simulating many phenomena emerging in engineering and physics.Even though many solutions have been suggested, none possesses all the desired properties and choices still need to be made.
Other issues emerge when trying to use those moment models, and that are not discussed in this short note.First, most of the presented closures are appropriately defined for the flux term, but many of them still require to be extended to the collision operator (together with a proper entropy study).Second, if constructing appropriate moment fluxes is complicated, constructing boundary conditions for moment equations remains an open problem and well-posedness results including such boundaries are rare.