First-order quarter- and mixed-moment realizability theory and Kershaw closures for a Fokker-Planck equation in two space dimensions

Mixed-moment models, introduced before for one space dimension, are a modification of the method of moments applied to a (linear) kinetic equation, by choosing mixtures of different partial moments. They are well-suited to handle such equations where collisions of particles are modelled with a Laplace-Beltrami operator. We generalize the concept of mixed moments to two dimension. The resulting hyperbolic system of equations has desirable properties, removing some drawbacks of the well-known $\MN[1]$ model. We furthermore provide a realizability theory for a first-order system of mixed moments by linking it to the corresponding quarter-moment theory. Additionally, we derive a type of Kershaw closures for mixed- and quarter-moment models, giving an efficient closure (compared to minimum-entropy models). The derived closures are investigated for different benchmark problems.


Introduction
The full discretization of kinetic transport equations like the Fokker-Planck equation is in general very expensive since the discretized variable resides in X × S 2 × [0, t f ] where X ⊂ R 3 and S 2 denotes the unit sphere in R 3 . Thus, the solution of the Fokker-Planck equation is very high-dimensional.
A common approach to reduce the dimensionality is given by the method of moments [11,22]. One chooses a set of angular basis functions, tests the Fokker-Planck equation with it and integrates over the angular variable, removing the angular dependence while getting a system of equations in x and t. Assuming now a specific form of the underlying distribution function, different approximate systems arise. Typical examples are the well-known spherical harmonics or P N models [6,11,18] and their simplifications, the SP N [16] methods. These models are computationally inexpensive since they form an analytically closed system of hyperbolic differential equations. However, they suffer from severe drawbacks: The P N methods are generated by closing the balance equations with a distribution function which is a polynomial in the angular variable. This implies that this distribution function might be negative resulting in non-physical values like a negative particle density. Additionally, in many cases a very high number of moments is needed for a reasonable approximation of the transport solution. This is in particular true in beam cases, where the exact transport solution forms a Dirac delta. The entropy minimization M N -models [1,5,9,23,25] are expected to overcome this problem since their closure functions are always positive. In many situations these models perform very well. Still they produce unphysical steady-state shocks due to a zero netflux problem.
To improve this situation, half or partial moment models have been introduced in [10,13]. These models work especially well in one space dimension, because they capture the potential discontinuity of the probability density in the angular variable which in 1D is well-located. If however, a Laplace-Beltrami operator is used instead of the standard integral scattering operator, i.e. scattering is extremely forward-peaked [26], these half moment approximations fail [14,29].
To improve this situation a new model with mixed moments was proposed in [14,29] which is able to avoid this problem. Instead of choosing full or half moments, a mixture of both is used. Contrary to a typical half moment approximation, the lowest order moment (density) is kept as a full moment while all higher moments are averaged over half-spaces. This ensures the continuity of the underlying distribution function.
Since in one spatial dimension mixed moments perform very well, it seems reasonable to extend them to multi-D. It turns out, that here a mixture of full, half and quarter moments is necessary to derive the correct mixed-moment ansatz.
Realizability is the fact that a vector of moments is physically relevant, i.e. that it is the moment of a nonnegative distribution function. While in 1D realizability theory for full moments [8] and mixed moments [29] is completely solved, it remains an open problem in higher dimensions. We will give a realizability theory for quarter and mixed moments of order 1 and derive, as in 1D, a corresponding Kershaw closure which provides an analytically closed system of equations, in contrast to minimum-entropy models which requires the solution of a nonlinear system of equations.
The paper is organized as follows: First, we will give a short introduction to the method of moments and its consequences for the Fokker-Planck equation. After setting up the basic notations in Section 2, we provide necessary and sufficient conditions of order 1 for realizability in the case of quarter moments and mixed moments in Section 3. Section 4 deals with the construction of Kershaw closures. Finally, we present some numerical tests for the derived models in Section 5 as well as conclusions and outlook in Section 6.

Macroscopic Models
We consider the Fokker-Planck equation which describes the densities of particles with speed Ω ∈ S 2 at position x ∈ X ⊆ R 3 and time t under the events of scattering (proportional to σ s (t, x)), absorption (proportional to σ a (t, x)) and emission (proportional to Q (t, x, Ω)). The equation is supplemented with initial condition and Dirichlet boundary conditions: where n is the outward unit normal vector in x ∈ ∂X.
Similar to [31] we assume that geometry, initial and boundary conditions are independent of the z-direction, resulting in a solution ψ which is also z-independent. Therefore (2.1) can be reduced to x ∈ X ⊆ R 2 .
Parametrizing Ω in spherical coordinates and taking the symmetry reduction into account we obtain where ϕ ∈ [0, 2π] is the azimuthal and µ ∈ [−1, 1] the cosine of the polar angle. Then the Laplace-Beltrami operator on the unit sphere can be written as Definition 2.1. The vector of functions b : S 2 → R n consisting of n basis functions b i , i = 0, . . . n − 1 of maximal order N (in Ω) is called an angular basis.
The so-called moments u = (u 0 , . . . , u n−1 ) T of a given distribution function ψ are then defined by where the integration is performed component-wise.
Equations for u can then be obtained by multiplying (2.1) with b and integration over S 2 : b∂ t ψ + b∇ x · Ωψ + bσ a ψ = σ s bC (ψ) + bQ Collecting known terms, and interchanging integration and differentiation where possible, the moment system has the form Depending on the choice of b the terms Ω x bψ , Ω y bψ and in some cases even bC (ψ) cannot be given explicitly in terms of u. Therefore an ansatzψ has to be made for ψ closing the unknown terms. This is called the moment-closure problem.
In this paper the ansatz densityψ is reconstructed from the moments u by minimizing the entropy-functional under the moment constraints bψ = u. (2.7) 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, α(u) := argmiñ 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 [21].
The kinetic entropy density η can be chosen according to the physics being modelled. As in [17,21], Maxwell-Boltzmann entropy is used, thus η * (p) = η * (p) = exp(p). This entropy is used for non-interacting particles as in an ideal gas or an ensemble of photons.
A closed system of equations for u remains after substituting ψ in (2.5) withψ u : Also note that using the entropy η(ψ) = 1 2 ψ 2 the linear ansatẑ remains. If the angular basis is chosen as spherical harmonics of order N , (2.12) turns into the classical P N model [4,6,31].

Angular bases
This moment approach strongly depends on the choice of the ansatzψ and the angular basis b. In the following sections we will shortly derive the different angular bases for the models presented here. These bases will be generally collected in the basis-vector b(Ω). If we need to further distinguish between the models, the corresponding symbols defined in the following sections will be used. If a result is independent on the choice of the basis we will just use b as symbol.

Full moments
The full-moment basis f N of order N consists of the tensorial powers of Ω, i.e. (by abusing notation) f N = 1, Ω, Ω ⊗ Ω, Ω ⊗3 , . . . , Ω ⊗N T . (2.14) The corresponding tensorial moment of order i is given by which consists of the scalar moments Note that an equivalent system can be obtained by using the corresponding real-valued spherical harmonics of order N as angular basis [4,6,31].

Quarter moments
We follow the approach in [12] where general partial as well as the special case of quarter moments in two space-dimensions are treated. The main idea is not to integrate over the whole sphere S 2 but over subsets of it. As for full moments the corresponding components of the tensorial moments are given by In this paper, D will be one of the following quarterspaces or halfspaces Note that for D = S 2 we have that u |i| D = u |i| .
For pure quarter moments, we will have D = S ij for i, j ∈ {+, −}. The corresponding basis for quadrant ij is then given by where · should be understood as multiplication with every component. Consequently the complete set of basis functions is p N = p S++ , p S−+ , p S−− , p S+− .

Mixed moments
As will be shown below the quarter-moment basis exhibits undesired properties when applied to the Laplace-Beltrami operator. This has inspired the works in [14,29] where so-called mixed moments were developed. The main problem of the quarter-moment basis is that the ansatzψ is not continuous in Ω which is necessary for the solution of (2.1) in one space-dimension. There, mixed moments where constructed by starting from a half-moment ansatz and demanding continuity of the ansatz (2.10) with respect to this basis. There, it suffices to choose a full zeroth-order moment and half-moments for all higher order moments [29].
The construction of mixed moments in two dimensions works in the same spirit. We start with the general quarter-moment basis p N and demand continuity of the ansatzψ u .
Having (2.10) in mind, we obtain multipliers α ij for every quadrant S ij . For example, for N = 1 we havê At the poles of the sphere (µ = ±1) we have Ω = (0, 0) T , which implies thatψ u is continuous only if Similarly it holds that along the quarter-space boundaries (i.e. ϕ = k π 2 , k = 0, . . . , 3) some of the multipliers have to be the same (exactly those whose component of Ω does not vanish on the boundary, e.g. α (1,0) S+− since at ϕ = 0 we have Ω x = 1 = 0). Accordingly we obtain the moments In contrast to the one-dimensional setting one full moment, half moments for the basis functions contributing to either x or y direction, and quarter moments for the basis functions which contribute to both directions occur.
Note that in the fully three-dimensional setting the decomposition in z-direction has to be taken into account, leading to octants instead of quadrants.
To embed this in the framework we choose our basis function m N of order N as Formulas to efficiently calculate the appearing integrals in case of a linear ansatz can be found in Appendix Appendix A.

Moments of the Laplace-Beltrami operator
All that remains to obtain a closed set of equations in (2.12) is to correctly evaluate b∆ Ωψu . This can be done using the formal self-adjointness of the Laplace-Beltrami operator, using Due to this, the calculation of these integrals does not a priori depend on the choice of the ansatz but is true for every ψ.

Full moments f N -basis
Consequently, the corresponding moments are given by

Quarter moments p N -basis
We observe the following relations on the quadrant S ++ for i x , i y > 0, i x + i y = i: Similarly, the corresponding quantities in the other quarter-spaces can be obtained.
The moments of ∆ Ω ψ include the evaluation of the microscopic values ∂ ϕ ψ and ψ at the quarter-sphere boundaries. Note that, similar to the one-dimensional case, the mass-conservation property of the Laplace-Beltrami operator (i.e. ∆ Ω ψ = 0) is usually violated. This can be easily seen by summing up (2.24) over all quadrants, observing that ∂ ϕ ψ at the angles kπ 2 , k = 0, . . . , 3 remains. Also note that this quantity is not rigorously defined for minimum-entropy closures due to the discontinuity in the ansatzψ.

Mixed moments m N -basis
We observe for i x , i y > 0 the following: The calculations for ∆ Ω 1 Sij Ω ix x Ω iy y are equivalent to those of the quarter-moment basis.
Note that all these calculations are closure-independent. We therefore need to calculate u

Realizability
In this section we will define precisely the concept of realizability. Furthermore we give necessary and sufficient conditions for first order quarter and mixed moment models.
Definition 3.1 (Realizability). A moment vector u is said to be realizable with respect to basis b if there exists a non-negative distribution ψ(Ω) ≥ 0 such that u = bψ . The realizable set is defined as ψ is then called a realizing distribution.
Note that (2.9) is solvable if and only if u ∈ R b .
A standard example for realizability are the realizability conditions of first order in the full-moment setting.

Quarter moments
In this section first order necessary and sufficient conditions for realizability of a quarter-moment vector will be given. This will be used later to derive the corresponding conditions for mixed moments.
it is necessary and sufficient for the existence of a non-negative measure ψ which realizes u with respect to p Sij that and the normalized first moment Proof. We only prove the statement for S ++ . The proof for the other quadrants works similarly. Assume that ψ ≥ 0 in S ++ . Since ||Ω|| 2 ≤ 1 we obtain showing the necessity of (3.3).
Since Ω ∈ S ++ we obtain Ω x , Ω y ≥ 0 implying that u To show the sufficiency of (3.3) we give a realizing distribution function. A possible (but not necessarily unique) candidate is given by Therefore, ψ S++ is a realizing distribution for u under the given assumptions. 1 We assume for notational simplicity that δ has mass 1 even on the boundary of integration.

Mixed moments
With the knowledge of Lemma 3.3 we are able to provide the realizability conditions for mixed moments of order 1.
Theorem 3.4 (First order necessary and sufficient conditions). For a vector of moments it is necessary and sufficient for the existence of a non-negative measure ψ which realizes u with respect to Proof. We start again with the necessity of (3.6): Note that the vector of component-wise absolute values |Ω| c : This can be reformulated to The left hand side will be extremal if ν is collinear to u, i.e. ν = u || u|| 2 . Thus (3.7) implies (3.6). The sign-constraints on the moments follow again from the signs of Ω x and Ω y in the corresponding halfspaces.
For sufficiency we give again a reproducing distribution. In the following we will make use that under (3.6) the reduced moment vector u is located in u (0,0) S ++ . Similarly the quantities will be in S ij . Reformulating (3.6) in terms of the normalized moments collected in φ |1| we see that . We now embed the moments φ |1| Sij ∈ R 2 into the normalized mixedmoment space (which is a subset of R 4 ) in the following way. Define Figure 1: Linear interpolation between two realizability boundaries in the projected three-space (φ and analogously for the other half-space combinations. This is visualized in Figure 1. It can be shown that φ |1| (and consequently u) can be written as a convex combination of the moments φ |1| S++ to φ |1| S+− . Indeed, defining we see that Therefore, due to the linearity of the problem, a non-negative realizing distribution for u is given by

Kershaw closures
A typical drawback of the minimum-entropy models defined by (2.10) is that the dual problem (2.9) cannot be solved analytically. The numerical solution, which has to be calculated at least once in every space-time cell, is challenging and expensive [3,17].
On the other hand, standard P N models may give physically irrelevant solutions since they do not ensure positivity of the underlying distribution function.
Due to this, Kershaw closures became recently a topic of increasing interest. They are constructed in such a way that they are automatically generated by a nonnegative distribution functionψ. Therefore, the moment vector including the unknown highest moment is also realizable with respect to the basis of one order higher. b is the isotropic moment. The last condition is also called the isotropic interpolation condition.
In one spatial dimension for a full-moment basis, the Kershaw closure is also exact on the realizability boundary, because there the realizing measure on the realizability boundary is unique. This property is no longer true in higher dimension (see e.g. [24]) or for other models. However, it gives an idea of how to construct such a closure.
The name of the closure is dedicated to David Kershaw who first proposed such an idea in [19]. For an introduction into Kershaw closures in one space dimension we refer to [29]. The construction of Kershaw closures of first order as done in [29] requires realizability information of second order. With this it is possible to linearly interpolate between different parts of the higher-order realizability boundaries (choosing the interpolation parameter in such a way that the isotropic moment is interpolated). The resulting model is then cheap to evaluate because it is analytically closed (in contrast to minimum-entropy models). Unfortunately, we were not able to provide a closed second-order realizability theory for mixed or quarter moments yet which implies that we can't identify the correct parts of this second-order realizability boundary for interpolation. However, it turns out that under some assumptions on φ |2| Sij the second-order realizability information from the half moments in one space dimension are sufficient to define the Kershaw closure for quarter moments in a similar fashion.
For mixed moments we choose a different approach to build the unknown second moment. Abusing the constructive procedure in Theorem 3.4 we are able to provide a closure for mixed moments by combining the quarter-moment closures accordingly.

Quarter moments
It was shown in [20] that, by assuming that the distribution function is symmetric around a preferred direction, the second moment of the M 1 closure can be decomposed into is the so-called Eddington factor. This implies that φ This is not necessarily true for all φ

|1|
Sij but has to be true on a subset of the quarterspace.
Then for any φ |2| Sij which is generated by a realizing distribution for φ |1| Sij we have that φ

|1|
Sij is an eigenvector Proof. We only prove the case S ij = S ++ . The other quadrants follow similarly. Assume that φ |1| S++ is generated by ψ ≥ 0.
Since Ω x , Ω y ≥ 0 we have the following. . On this support we have Ω x = 0, which implies that Similar to the one-dimensional case we observe that (keeping the support of ψ in mind) Therefore the eigenvalue λ := φ (0,2) S++ to the eigenvector φ |1| S++ has to satisfy φ |1| Sij it follows that ψ has support only on ϕ = π 4 and µ ∈ [−1, 1]. On this support we have Ω x = Ω y , which implies that In this case we have u (2,0) while the lower bound is u  (1,0) (d) This follows exactly in the same way as for full moments, see e.g. [24]. Note that here, φ |1| Sij Numerical tests suggest that for the QM 1 model the first normalized moment is indeed always close (but not equal) to an eigenvector of the second normalized moment, i.e. < ) u |1| Sij , v 1 remains small (where v 1 is the eigenvector of u |2| Sij with smaller enclosing angle between itself and u |1| Sij ), see Figure 2. This has been checked numerically using the QM 1 -tabulation strategy given in [12] Since the structure of eigenvectors of the QM 1 model is not immediately obvious to us, it seems reasonable to assume the simple form (4.1) for QK 1 as well: where the coefficients α 1 (u), α 2 (u) ∈ R ≥0 have to be chosen accordingly. The eigenvalue associated to φ |1| Sij is given by λ = α 1 + α 2 .
A short calculation shows that under this assumptions the bounds φ |1| Sij have to be fulfilled for all φ |1| Sij ∈ S ij , and not only on the subsets of S ij shown in Lemma 4.1. For every ν ∈ S ij it holds that where the last inequality is meant component-wise. This implies for ν = The lower bound on λ follows in the same way as for full moments observing that u The construction of the Kershaw closure in one spatial dimension linearly interpolates between the upper and lower boundary values for the second moment in such a way that the isotropic point is correctly hit. We define this interpolation not on the second moment but on its eigenvalues. The isotropic moment vector is given by φ . We therefore define 1] has to be chosen accordingly. Plugging in the isotropic point we get that Although other choices are possible we furthermore assume that α 1 + α 2 is rotationally symmetric, which implies that ζ = ζ φ |1| Sij 2 . Unfortunately, the simplest choice ζ ≡ ζ 1 2 , 1 2 leads to a moment system which is no longer hyperbolic. In fact, some of the eigenvalues of the flux Jacobian are imaginary, especially in those regions where the QM 1 eigenvectors deviate strongly from φ |1| Sij (compare Figure 2). We fixed this problem by choosing ζ to be linear in φ The proof of Lemma 4.1 reveals more information about the choice of α 1 . In situations (a), (b) and (d) (i.e. on the realizability boundary), α 1 has to be zero, while it has to be 1 3 − 2 3π at the isotropic point (recall that there α 2 = 2 3π due to the off-diagonals of the second moment). The simplest function fulfilling these two conditions is given by The corresponding closure relation is shown in Figure 3. Despite the fact that the assumption on the eigenvectors for QM 1 is not true for all φ |1| Sij the resulting second moment for QK 1 is very close to the minimum-entropy one. This is shown in Figure 4.
Since we are not able to provide a closed second-order realizability theory we cannot conclude immediately that this closure is in fact realizable with respect to p 2 . We therefore have to go a different way and, as before, provide a realizing distribution for this special case. We only prove the case S ij = S ++ and drop the lower index S ++ for readability issues. We make the ansatz 2 where ι ∈ {1, 2} and the parameters c ι± , ρ ι± ∈ R andv ι ∈ R 2 still have to be determined. Without loss of generality we assume u (0,0) = 1 for brevity. The moments of this ansatz are Hierarchically inserting the equations into the one of higher order we see that the coefficients c ι± , ρ ι± fulfil the system of equations Since ψ has to be supported in S ++ and should be nonnegative we also get the inequality-constraints Combining (4.6a) and (4.5b) shows that sign (ρ ι+ ) = sign (ρ ι− ), thus, without loss of generality, the quarter- sphere constraints (4.6b) can be rewritten as A visualization of these bounds can be found in Figure 5 for a specific example. The min function ensures that the first intersection with either the coordinate axes (ρ 1+ in the example) or the norm-bound (ρ 1− in the example) is taken. For ρ 2± there is no distinction necessary since we can only intersect with one of them.
We want to prove this exemplarily for ρ 2+ . After some symbolic simplifications we see that With this it immediately follows that 1 − φ |1| Note that the above choice of the coefficients ρ ι± is not unique. Since the corresponding values for c ι± are also non-negative, (4.4) is a non-negative distribution function with support in S ++ which realizes the desired moments. Therefore, the QK 1 closure is realizable.

Mixed moments
Since the mixed-moment setting is even more complicated than the quarter-moment one, we want to follow a different approach to construct an approximate closure. As has been shown in Theorem 3.4 it suffices to have a non-negative quarter-moment distribution ψ Sij to construct a mixed-moment reproducing distribution. If these distributions fulfil the isotropic moment interpolation condition for the quarter-moment basis, by linearity, the resulting mixed-moment distribution (3.12) will fulfil it as well for the mixed-moment basis.
Such a distribution function for quarter moments has been found in (4.4). Thus setting φ |1| Sij as in (3.8) will provide the desired mixed-moment distribution function ψ m1 using the interpolation defined by (3.11).
A closure can now be generated by calculating the second moments with respect to those quarter-space distributions which have support in the corresponding set D. For example, Sij .

Remark 4.2.
Numerical experiments suggest that this closure is also (weakly) hyperbolic. No imaginary eigenvalues of the flux Jacobian have been found on a finely discretized grid of the realizable set. However, at the isotropic point both flux Jacobians have three times the eigenvalue 0. Still, the transformation matrix has full rank, i.e. the geometric multiplicity of the eigenvalue 0 is 3.

Remark 4.3.
Any other choice of ψ Sij which fulfils the desired properties with respect to the quarter-moment basis will be a feasible choice as well. Exemplarily, the quarter-moment minimum-entropy closure QM 1 shall be mentioned. An efficient implementation using tabulation is given in [12].
Note that this closure is not an approximation of the MM 1 model and might behave differently. This is shown exemplarily in Figure 6 where the second moment φ The corresponding shape of the Kershaw closure is very different to those of the MM 1 model.

Treatment of the Laplace-Beltrami operator
A clear drawback of the mixed-moment Kershaw closure is that the moments of the Laplace-Beltrami operator (see Section 2.2.3) require the evaluation of the ansatz function. When using the QK 1 distribution to realize the mixed-moment closure, the ansatz is a linear combination of Dirac deltas, which cannot be evaluated in a feasible way. One could follow the authors in [14] and replace the ansatz in the Laplace-Beltrami moments by the corresponding polynomial distribution (2.13). Then, we have (4.8) Unfortunately, this closure may lead to an undesired behaviour of the moment-system (2.12), since its solution may leave the realizable set (this can be easily shown even in one spatial dimension with the ansatz given in [14]). This follows from the fact that the tangent field spanned by Su may point outside of the realizable set at those parts of the boundary where the polynomial ansatz (2.13) is negative.
Consider for example u = (1, 1, 0, 0, 0) T . Then we have that We therefore modify this approach slightly and replace the Dirac ansatz with the tabulated QM 1 distribution. Since the fluxes for QM 1 and QK 1 are very similar, it can be expected that the error done with this replacement is negligible. Furthermore, since the ansatz is positive, the moment-system (2.12) will behave as expected. Indeed, calculating the approximation for the previous example we get approximately which points inside of the realizable set. Note the fourth and fifth component which converges to ±∞. This is consistent with the original problem since on the realizability boundary, the exponential ansatz converges to a combination of Dirac deltas.

Numerical results
We present the derived models in several benchmark test-cases. The numerical discretizations are obtained with a two-dimensional generalization of the high-order realizability-preserving discontinuous-Galerkin scheme given in [2] for the Kershaw models and a generalization of the realizability-preserving kinetic scheme given in [30] for minimum-entropy models. Those generalizations are in principle straight-forward but might be topic of a follow-up paper. The spatial and temporal order is four where roughly 10000 rectangular (discontinuous-Galerkin scheme) and 20000 triangular elements (kinetic scheme) in space were used. The reference solution is discretized with the second-order scheme given in [7,27,28].

Line Source
The line-source test is a Green's function problem, where a pulse of particles is emitted in an infinite medium [6]. It has been widely investigated for isotropic scattering in [15]. We choose the following set of parameters for this problem, smoothing the initial Dirac delta in space: • Parameters: σ a = Q = 0, σ s = 1 • Initial condition: ψ t=0 (x, Ω) = max( exp − x 2 +y 2 2 σ 2 8 π σ 2 , 10 −4 ) with σ = 0.03 • Boundary conditions: Isotropic in Ω, consistent with initial conditions. We show several model solutions in Figures 7 and 8. We observe that convergence is achieved quickly with increasing moment order N for M N models which can be observed as well in one spatial dimension. This is a huge contrast to the integral-scattering operator used in [15], where only slow convergence can be observed. Note that, by construction, the MM 1 model has no rotational symmetry. The main directions of propagation follow the half-spaces, i.e. along the cartesian axes. Still, it performs significantly better than M 1 (especially along the diagonal cut shown in Figure 5.8b), but worse than M 2 which almost converged to the reference solution. We observed that the naive implementation of higher-order mixed-moment closures like in the MM 2 model are numerically instable due to errors in the non-linear closure of the Laplace-Beltrami operator (see Section 2.2.3). The polynomial closure (4.8) has a slightly lower peak speed as the tabulated closure but more interesting is that the latter is much closer to the MM 1 solution (Figure 5.7d). This shows, that in contrast to the one-dimensional situation in [14], the choice of the Laplace-Beltrami closure has a significant effect on the solution of the mixed-moment system.
Since σ s > 0, the QK 1 model cannot be applied here.

Two Beams
This is the typical situation where the M 1 model completely fails and is an often-used benchmark problem in 1D (see e.g. [14,29]). Instead of using opposing beams we investigate two orthogonal beams hitting each other in a void. We show two different variants of this. The first one is given by • Final time: t f = 1.2 • Parameters: σ s = σ a = Q = 0 • Initial condition: ψ t=0 (x, Ω) = 10 −4 4π • Boundary conditions: where σ 2 = 0.05. Figure 9 shows the solutions of full-and mixed-moment minimum-entropy models up to order N = 3, In the second version, the whole setup is rotated clockwise. The two beams sit in the upper-left and lowerleft corner, respectively, both pointing to the center of the domain. Due to the rotational symmetry of the full-moment models and the reference solution, only the MK 1 closure is shown in Figure 11. It is visible that in contrast to the first version of this problem (compare Figure 5.9h) no shock is produced when the beams hit.

Conclusions and future work
We developed a two-dimensional extension of mixed moments, leading to a generalization of the original minimum-entropy MM N model, proposed in [14,29]. Additionally, we provided a first-order realizability theory for mixed as well as quarter moments. Since the numerical inversion of the general minimum-entropy moment problem is very expensive we developed approximate, but still realizable, closures, the so-called quarter-and mixed-moment Kershaw closures. They have approximately the same cost as the corresponding polynomial closure while being realizable and non-linear, better adapting to the correct eigenvalues of the true Fokker-Planck solution. Although first-order mixed moments completely resolved the zero-netflux problem of the M 1 model, a first order approximation in two dimensions is in general not enough, as has been shown in the two-beams test-case. Under some assumptions on the alignment of the beams this drawback can be completely removed. On the other hand, going to higher-order approximations, like M N or MM N with N ≥ 2, a much better approximation can be done. Thus it seems reasonable to extend the concept of Kershaw closures for either K N or MK N for N ≥ 2 to obtain a cheap and accurate approximation of the Fokker-Planck equation. Furthermore, more research is necessary to obtain better-quality high-order numerical solutions of Kershaw closures.
where B (z, w) is the beta function.