Second-order mixed-moment model with differentiable ansatz function in slab geometry

We study differentiable mixed-moment models (full zeroth and first moment, half higher moments) for a Fokker-Planck equation in one space dimension. Mixed-moment minimum-entropy models are known to overcome the zero net-flux problem of full-moment minimum entropy $M_N$ models. Realizability theory for these modification of mixed moments is derived for second order. Numerical tests are performed with a kinetic first-order finite volume scheme and compared with $M_N$, classical $MM_N$ and a $P_N$ reference scheme.


Introduction
We investigate time-dependent kinetic transport equations like the Fokker-Planck equation, arising from the Boltzmann equation [3,6] under the assumption of extremely forward-peaked scattering [30]. They describe the propagation of "radiation particles" like photons or electrons which travel at time t from their current position in a specific direction and how they interact with the surrounding matter. Without any assumptions or dimensional reductions this typically leads to a six-or seven-dimensional state space. Applications reach from electron transport in solids and plasmas, neutron transport in nuclear reactors, photon transport in superfluids and radiative transfer to the context of biological modelling, e.g. for studying cell movement (chemotaxis/haptotaxis) or wolf migration [7,18,21].
A common approach to reduce the dimensionality is given by the method of moments [12,26], which is a class of Galerkin methods for the approximation of such time-dependent kinetic transport equations. One chooses a set of angular basis functions, tests the kinetic equation against it and integrates over the angular variable, removing the angular dependence while getting a (potentially huge) system of differential equations in space and time. Well-known examples are the classical P N methods [5,12,22], their simplifications, the SP N [17] methods and entropy minimization M N models [1,4,9,28,29]. Especially the latter is favourable since the moment equations are always closed with a positive ansatz function, respecting the positivity of the kinetic distribution to be approximated. In many situations these models perform very well, but since they result from averaging over the complete velocity space, they can produce physically wrong steady-state shocks. It has been shown by Hauck [19] that these shocks exist for every odd order.
To improve this situation, half-or partial-moment models were introduced in [11,14]. These models work especially well in one space dimension since they capture the potential discontinuity of the probability density in the angular variable which in 1D is well-located. Unfortunately, in a Fokker-Planck operator is used instead of the standard integral-scattering operator (BGK type), these half-moment approximations fail significantly. A reason for this is that the domain of definition of the Laplace-Beltrami operator requires continuous functions in one dimension. [35].
An intermediate model respecting the continuity of full-moment models while allowing the flexibility of partial moments is the mixed-moment model, which was proposed in [16,35,37]. Contrary to a typical half-moment approximation, the lowest order moment (density) is kept as a full moment while all higher moments are half moments.
Although these MM N models satisfy the above-mentioned property of having a continuous ansatz function, the numerical discretization of it is highly non-trivial due to the appearance of microscopic terms (i.e. the moments of the Laplace-Beltrami operator depend on the values of the ansatz itself). Especially in multiple dimensions, naive implementations fail at discretizing the (semi-)microscopic quantities (line integrals over quadrant/octant boundaries) [34,37]. To overcome this (numerical) problem, we investigate a modification of the mixed-moment model. This new DMM N model has more regularity, i.e. its ansatz is differentiable, resulting in a more robust numerical implementation while maintaining most of the benefits of the classical MM N model.
The first part of the paper shortly reviews the method of moments and the minimum-entropy ansatz. Afterwards, the concept of realizability (the fact, that a moment vector is associated with a non-negative distribution function) is introduced and a concrete characterization of the realizable set for the DMM 2 model is derived. Furthermore, the eigenstructure of this model is explored. Then, the performance of the new model is investigated in two benchmark tests, showing that the DMM 2 is competitive compared to M N and MM N models with the same number of degrees of freedom. The paper is concluded by a summary and an outlook on future work.

Models
In slab geometry, the transport equation under consideration for the particle distribution ψ = ψ(t, x, µ) has the form The physical parameters are the absorption and scattering coefficient σ a , σ s : T × X → R ≥0 , respectively, and the emitting source Q : Collision of particles is modelled by the Laplace-Beltrami operator This operator appears, for example, as the result of an asymptotic analysis of the Boltzmann equation under the assumption of small energy loss and deflection, and forward-peaked scattering in the context of electron transport [16,20,30].
The transport equation (2.1) is supplemented by initial and boundary conditions: In general, solving equation (2.1) is very expensive in two and three dimensions due to the high dimensionality of the state space.
For this reason it is convenient to use some type of spectral or Galerkin method to transform the highdimensional equation into a system of lower-dimensional equations. Typically, one chooses to reduce the dimensionality by representing the angular dependence of ψ in terms of some basis b. The so-called moments of a given distribution function ψ with respect to b are then defined by where the integration · := 1 −1 · dµ is performed componentwise.
Collecting known terms, and interchanging integrals and differentiation where possible, the moment system has the form The solution of (2.5) is equivalent to the one of (2.1) if b is a basis of L 2 (S 2 , R).
Since it is impractical to work with an infinite-dimensional system, only a finite number of n < ∞ basis functions b of order N can be considered. Unfortunately, there always exists an index i ∈ {0, . . . , n−1} such that the components of b i · µ are not in the linear span of b. Therefore, the flux term cannot be expressed in terms of u without additional information. Furthermore, the same might be true for the projection of the scattering operator onto the moment-space given by bC (ψ) . This is the so-called closure problem. One usually prescribes some ansatz distributionψ u (t, x, µ) :=ψ(u(t, x), b(µ)) to calculate the unknown quantities in (2.5). Note that the dependence on the angular basis in the short-hand notationψ u is neglected for notational simplicity.
In this paper the ansatz densityψ is reconstructed from the moments u by minimizing the entropy-functional under the moment constraints The kinetic entropy density η : R → R is strictly convex and twice continuously differentiable and the minimum is simply taken over all functions ψ = ψ(µ) such that H(ψ) is well defined. The obtained ansatẑ ψ =ψ u , solving this constrained optimization problem, is given bŷ This problem, which must be solved over the space-time mesh, is typically solved through its strictly convex finite-dimensional dual, where η * is the Legendre dual of η. The first-order necessary conditions for the multipliers α(u) show that the solution to (2.8) has the formψ This approach is called the minimum-entropy closure [25]. The resulting model has many desirable properties: symmetric hyperbolicity, bounded eigenvalues of the directional flux Jacobian and the direct existence of an entropy-entropy flux pair (compare [25,34]).
The kinetic entropy density η can be chosen according to the physics being modelled. As in [19,25], Maxwell-Boltzmann entropy is used, thus η * (p) = η * (p) = exp(p). This entropy is used for non-interacting particles as in an ideal gas.
Substituting ψ in (2.5) withψ u yields a closed system of equations for u: For convenience, (2.12) can be written in the form of a usual first-order hyperbolic system of balance laws In this paper, a variant of the so-called mixed-moment basis [16,35] is used. This ansatz is a combination of the full-moment (b i = µ i ) and half-moment monomial basis [10,11].
The classical mixed-moment basis consists of a full zeroth moment and half moments for every higher moment. The resulting ansatz (2.10) is continuous but not continuously differentiable in µ = 0, leading to a microscopic term of the formψ u (0) in the scattering term b∆ µψu [16,35]. While this can be treated easily in one dimension, discretization problems arise in higher dimensions, where the microscopic quantity has to be replaced by an integration over a spherical arc of the unit sphere [34,37].
For this reason, we modify the mixed-moment basis in such a way that the ansatz is differentiable in µ = 0, removing the microscopic quantity. In one dimension, it suffices to choose a full zeroth and first moment to obtain the desired regularity. The corresponding moments have the form · dµ denote integration over the halfspaces. Accordingly, the angular basis Using this basis, it holds that Note that for l = 2 and l = 3 the quantities u 0± = ψ ± and u 1± = µψ ± , respectively, appear, which have to be determined using the closure relation (2.10).
Definition 2.2. The classical mixed-moment model will be referred to as the MM N model, while the differentiable mixed-moment model will be called the DMM N model.  Figure 1 shows typical ansatz functionsψ for the MM 2 and DMM 2 model. It can be seen that the MM 2 ansatz is only continuous, while the DMM 2 ansatz is also continuously differentiable in µ.

Realizability
Since the underlying kinetic density to be approximated is non-negative, a moment vector only makes sense physically if it can be associated with a non-negative distribution function. In this case the moment vector is called realizable.
If u ∈ R b , then u is called realizable. Any ψ such that u = bψ is called a representing density.

Remark 3.2.
(a) The realizable set is a convex cone, and (b) Representing densities are not necessarily unique.
Additionally, since the entropy ansatz has the form (2.10), in the Maxwell-Boltzmann case, the optimization problem (2.8) only has a solution if the moment vector lies in the ansatz space In the case of a bounded angular domain, the ansatz space A is equal to the set of realizable moment vectors [23]. Therefore, it is sufficient to focus on realizable moments only.
Unfortunately, the definition of the realizable set is not constructive, making it hard to check if a moment vector is realizable or not. Therefore, other characterizations of R b are necessary.
For example, in the classical mixed-moment problem of first order, the realizable set is characterized by the inequalities [16,35] u 1+ − u 1− ≤ u 0 and ± u 1± ≥ 0.
In this paper, we want to focus on the lowest-order non-trivial model of the differentiable mixed-moment hierarchy, i.e. N = 2.
Proof. At first, we want to show that (3.1) and (3.2) are necessary. Assume that ψ ≥ 0 is arbitrary but fixed and u = bψ . Note that u 0± ≥ ±u 1± ≥ u 2± ≥ 0 and u 0± u 2± ≥ u 2 1± due to the half-moment realizability conditions [8,35]. Then we have (using u 0− + u 0+ = u 0 ) that Since u 1− ≤ 0 it follows that The upper bound can be shown to be necessary in a similar way.
(3.2) follows from the positivity of 1 and µ 2 . We want to remark that the standard second-order full-moment realizability condition for u 2 = u 2+ + u 2− , namely u 0 (u 2+ + u 2− ) ≥ u 2 1 , is implied by (3.1) and (3.2). To show that the above inequalities are also sufficient, we provide a non-negative realizing distribution with support in [−1, 1]: In this case, ψ = u 0 δ(µ) (due to the quadratic term the second moment vanishes faster than the first moment so we have φ2± φ1± → 0 and This corresponds to the standard half-moment realizability conditions of second order. We first note that (3.1) and (3.2) imply that φ 2± ∈ [0, 1] since otherwise the bounds become complex. Second, we have that implying the classical full-moment realizability conditions φ 1 ∈ [−1, 1] and φ 2 ≤ 1.
We start the investigation at the different parts of the realizability boundary.
Plugging this into the definition of φ 1+ we get that, after some elementary transformations, Similarly, we obtain −φ 1− ≥ φ 2− and the same in the case Since the realizable set is always convex, the argumentation must also hold in the interior of the above set.
The normalized realizable set of the DMM 2 model, defined by (3.1) and (3.2), is shown in Figure 2.
gives a surprising insight into realizability of mixed-moment models. While the realizable set for full-moment and classical mixed-moment models can be characterized by inequalities with rational functions of the moments, the differential mixed-moment model requires non-linearities. This implies that it might be impossible to transfer the general mixed-moment structure (which uses a linearity argument) shown in [35] to the differentiable case.
The four eigenvalues λ 1 ≤ λ 2 ≤ λ 3 ≤ λ 4 of (4.2), which have been obtained numerically, are shown in Figures 3 to 6. In Figure 3, the eigenvalues are shown along the cut which is exactly the mean of the upper and lower bound on φ 1 .
It can be seen that the eigenvalues are discontinuous in the degenerate corners of the realizable set (e.g. λ 4 at φ 2− = 1 = −φ 1 , φ 2+ = 0). This property exists also for the classical mixed-moment MM 1 or the M 2 model [34].
We investigate the hyperbolicity of the moment system in Figure 4 by comparing the distances of adjacent eigenvalues. Figure 4a shows that all eigenvalues coincide only if φ 1 = φ 2+ = φ 2− = 0. Otherwise, at least two eigenvalues differ from each other.
The results in Figure 4b propose the strict hyperbolicity of the moment system in the interior of the realizable set (since all eigenvalues differ). However, at the realizability boundary (e.g. φ 2+ + φ 2− = 1) at least two eigenvalues coincide.
Since the optimization problem (2.9) is ill-conditioned close to the realizability boundary, the calculation of the multipliers α at this part of the realizable set is error-prone or impossible, resulting in meaningless pictures. We therefore investigate isotropically-regularized moments where an increase of the regularization parameter r ∈ [0, 1] moves the original moment vector φ towards the isotropic moment vector (in case of the DMM 2 model: φ iso = 0, 1 6 , 1 6 T ). Figure 5 shows the eigenvalues for r = 0.05 and φ ∈ ∂R b ρ=1 .
Similar to Figure 4b, Figure 6 shows the minimal distance of the 5% regularized boundary moments. It is visible that the minimal distance is attained at φ 2+ + φ 2− = 1, which indicates that on this part of the boundary the moment system is only weakly hyperbolic.
However, an analytical investigation of the eigenvalues in the limit cases

Numerical experiments
We use the first-order, realizability-preserving, implicit-explicit kinetic scheme derived in [31]. All results are computed on a grid with 1000 points. The reference solution is given by the P 99 model [27]. 0 0.5

Plane source
In this test case an isotropic distribution with all mass concentrated in the middle of an infinite domain x ∈ (−∞, ∞) is defined as initial condition, i.e.
where the small parameter ψ vac = 0.5 × 10 −8 is used to approximate a vacuum. In practice, a bounded domain must be used which is large enough that the boundary should have only negligible effects on the solution. For the final time t f = 1, the domain is set to X = [−1.2, 1.2] (recall that for all presented models the maximal speed of propagation is bounded in absolute value by one).
At the boundary the vacuum approximation is used again. Furthermore, the physical coefficients are set to σ s ≡ 1, σ a ≡ 0 and Q ≡ 0.
All solutions are computed with an even number of cells, so the initial Dirac delta lies on a cell boundary. Therefore it is approximated by splitting it into the cells immediately to the left and right. In

Source beam
We present a discontinuous version of the source-beam problem from [15], as in [2,36]. The spatial domain is X = [0, 3], and The final time is t f = 2.5. As above, the results for M N , MM N and DMM N models are shown in Figure 8.
As has been remarked in [34], the MM 1 and the M 2 model coincide well in this situation. Surprisingly, a similar statement is valid for the DMM 2 and M 3 model. As before, the DMM 2 model behaves qualitatively above the level of the MM 1 and M 2 model and similarly or slightly below those of the MM 2 model (which has the highest number of degrees of freedom).

Conclusions and outlook
We have derived the DMM N model and its associated realizability domain R b for N = 2. Numerical results suggest that, despite having one degree of freedom less, the DMM 2 model performs comparable to the MM 2 model. The key advantage of this class of moment models is that in the approximation of the Laplace-Beltrami operator only macroscopic quantities occur, whereas microscopic terms are present in the classical mixed-moment model. While this appears to have no significant impact in one dimension, where the position of the microscopic term is well-located, a more stable numerical approximation can be expected in two or three dimensions.
Future work should include the derivation of realizability theory for moment-orders N ≥ 3, to gain more insight into the arising non-linearities in this modified problem. Furthermore, the DMM N should be investigated in higher dimensions, especially in the context of the Fokker-Planck operator. The results in [34,37] indicate that mixed moments are hardly applicable in this framework due to the difficulty in the discretization of the Laplace-Beltrami operator. This should be avoidable using the differentiable basis functions. Finally, Kershaw closures [24,32,33,35] should be investigated to improve the efficiency of the DMM N model by avoiding the need to solve the moment system (2.7).