A realizability-preserving high-order kinetic scheme using WENO reconstruction for entropy-based moment closures of linear kinetic equations in slab geometry

We develop a high-order kinetic scheme for entropy-based moment models of a one-dimensional linear kinetic equation in slab geometry. High-order spatial reconstructions are achieved using the weighted essentially non-oscillatory (WENO) method, and for time integration we use multi-step Runge-Kutta methods which are strong stability preserving and whose stages and steps can be written as convex combinations of forward Euler steps. We show that the moment vectors stay in the realizable set using these time integrators along with a maximum principle-based kinetic-level limiter, which simultaneously dampens spurious oscillations in the numerical solutions. We present numerical results both on a manufactured solution, where we perform convergence tests showing our scheme converges of the expected order up to the numerical noise from the numerical optimization, as well as on two standard benchmark problems, where we show some of the advantages of high-order solutions and the role of the key parameter in the limiter.

1. Introduction. In recent years many approaches have been considered for the solution of time-dependent linear kinetic transport equations, which arise for example in electron radiation therapy or radiative heat transfer problems. Many of the most popular methods are moment methods, also known as moment closures because they are distinguished by how they close the truncated system of exact moment equations. Moments are defined through angular averages against basis functions to produce spectral approximations in the angle variable. A typical family of moment models are the so-called P N -methods [16,25] which are pure spectral methods. However, many high-order moment methods, including P N , do not take into account that the original kinetic density to be approximated must be nonnegative. The moment vectors produced by such models are therefore often where σ a and σ s are the nonnegative absorption and scattering interaction coefficients, and S a source. The operator C is a collision operator, which in this paper we assume to be linear and have the form (2.2) We assume that the kernel T > 0 is normalized to 1 −1 T (µ , µ)dµ ≡ 1. A typical example is isotropic scattering, where T (µ, µ ) ≡ 1/2. Equation (2.1) is supplemented by initial and boundary conditions: where ψ L , ψ R , and ψ t=0 are given.

2.2.
Moment equations and entropy-based closures. Moment equations are an angular discretization for (2.1), where the moments themselves are defined by angular averages against a set of basis functions. We use the following notation for angular integrals: for any integrable function φ = φ(µ); and therefore if we collect the basis functions into a vector m = m(µ) = (m 0 (µ), m 1 (µ), . . . , m N (µ)) T , the moments of a kinetic density φ are given by u = mφ . In this paper we consider the monomial moments m i (µ) = µ i , though all results can be extended to other bases, including, for example, partial [11,12] or mixed moments [13,32].
The closed system of moment equations is a system of partial differential equations of the form ∂ t u + ∂ x f (u) + σ a u = σ s r(u) + mS , (2.4) where the moment vector u(t, x) approximates mψ for the kinetic density ψ satisfying (2.1). In an entropy-based closure (commonly referred to as the M N model or the Levermore closure after he exposed their general structure in [24]), the functions f and r have the form f (u) := µmψ u and r(u) := mC(ψ u ) .
Hereψ u is an ansatz density reconstructed from the moments u by solving the constrained optimization problem: where the kinetic entropy density η is strictly convex and the minimum is simply taken over functions φ = φ(µ) such that η(φ) and mφ are well defined. This problem is typically solved through its strictly convex, unconstrained, finitedimensional dual,α (u) := argmin where η * is the Legendre dual of η. The first-order necessary conditions forα (u) show that the solution to (2.5) has the form ψ u = η * m Tα (u) (2.7) where η * is the derivative of η * . The kinetic entropy density η can be chosen according to the physics being modeled.
While in a linear setting such as ours, indeed any convex entropy η is dissipated by (2.1). As in [18] we focus on the Maxwell-Boltzmann entropy, because not only is it physically relevant for a wide variety of problems, but in particular it gives a positive ansatzψ u , since η * (y) = η * (y) = exp(y), and thuŝ ψ u = exp m Tα (u) . (2.8) Another standard choice for the entropy is η(z) = 1 2 z 2 , which yields the wellknown P N equations [18,25]. The P N equations are linear, and since η * (y) = y, the optimization problem (2.6) can be solved by hand. However the resulting ansatz is simply a linear combination of the basis polynomials and thus is not necessarily nonnegative. 1 The incorporation of the boundary conditions (2.3) is neither obvious nor trivial and is still an open problem [22,23,29,35]. This is not a focus of our work here, and therefore we only use a simple approach from previous work on entropy-based moment closures [2,18], which we discuss below in Section 3.2.
2.3. Moment realizability. Since the underlying kinetic density we are trying to approximate is nonnegative, a moment vector only makes sense physically if it can be associated with a nonnegative density. In this case the moment vector is called realizable. Additionally, since the entropy ansatz has the form (2.8), the optimization problem (2.5) only has a solution if the moment vector lies in the ansatz moment space In our case, where the domain of angular integration is bounded, the ansatz moment space A is exactly equal to the set of realizable moment vectors [19]. Therefore we can focus simply on realizable moments: Any φ such that u = mφ is called a representing density.
The realizable set is a convex cone. In the monomial basis, a moment vector is realizable if and only if its corresponding Hankel matrices are positive definite [9,33].
In general, angular integrals cannot be computed analytically. We define a quadrature for functions φ : Then the numerically realizable set is [3] Indeed, when replacing the integrals in the optimization problem (2.5) with quadrature, a minimizer can only exist when u ∈ R Q m . Below we often abuse notation and write φ when in implementation we mean its approximation by quadrature. We also liberally use the term realizable either to mean realizability with respect to R m or with respect to R Q m , where the specific meaning depends on whether exact integrals or those approximated by quadrature are meant in the context.

3.
A high-order kinetic scheme. A kinetic scheme for (2.4) can be thought of as first defining a spatial discretization for the underlying kinetic equation (2.1) and subsequently performing the angular discretization with the moment closure. We divide the spatial domain (x L , x R ) into a (for simplicity) uniform grid of J cells I j = (x j−1/2 , x j+1/2 ), where the cell edges are given by x j±1/2 = x j ± ∆x/2 for cell centers x j = x L + (j − 1/2)∆x, and ∆x = (x R − x L )/J. For (2.1) we define a finite-volume scheme for the cell means which with the Godunov (or 'upwind') numerical flux gives: where ψ − j±1/2 and ψ + j±1/2 in the flux terms denote the values of the approximate solution at the cell edges x j±1/2 from the left and right, respectively, and we generally use the bar with subsequent subscript j, i.e. · j , to indicate a cell average over the j-th cell as in (3.1).
To obtain a high-order scheme in space one only has to give a high-order reconstruction of the point-values of ψ, the distribution underlying the cell means, not only at the edge values for the flux terms, but also throughout the cells when σ a or σ s depends on x. We use the popular weighted essentially non-oscillatory (WENO) reconstruction method [34], which gives a polynomial reconstruction of ψ from the cell averages of the j-th cell and its neighbors. Now we perform the moment closure by replacing ψ with the entropy ansatzψ in (2.7), multiplying through by the angular basis functions, and integrating out the angle: where

FLORIAN SCHNEIDER, JOCHEN KALL AND GRAHAM ALLDREDGE
and s = mS . Hereψ + j±1/2 denote the evaluations at the cell edges of the WENO reconstructions made using the entropy ansätze of the neighboring cells evaluated from the right, and respectively forψ − j±1/2 evaluated from the left. The first step in the scheme, then, is to compute the ansatz for each cell.
3.1. Numerical optimization for angular reconstruction. In order to computeψ j at the cell means, we first compute the multipliersα(u j ) by solving the dual problem (2.6). The gradient and Hessian of the objective function are We use the numerical optimization techniques proposed in [3]. We assume that the moment vector u j has been scaled so that its zeroth component is one. The optimizer stops at the first iterate α which satisfies where · is the Euclidean and · 1 the 1−norm in R N +1 , τ and ε user-specified tolerances, and u 0 (α) = exp(m T α is the zero-th order moment associated with the multiplier vector α. This criterion is similar to the one in [3] but modified for the following reasons. Unlike the algorithm there, we modify the zero-th component of the final multipliers so that the zero-th moment (which gives the local density) is exactly matched. That is, while the optimizer stops at the first α which satisfies (3.4), the multiplier vector it returns (to define ψ j for the kinetic reconstructions in Section 3.2) is α := α 0 − log(u 0 (α)), α 1 , . . . , α N T . (3.5) Since m 0 ≡ 1, this gives m exp(m T α) = m exp(m T α) /u 0 (α), which ensures exp(m T α) = 1. The form of τ on the right-hand side of (3.4a) ensures that g(α) is bounded by τ : where we have used that (3.4a) implies |u 0 (α)−1| < τ and consequently 1/u 0 (α) < 1/(1 − τ ). The second condition (3.4b) in the stopping criterion enforces an approximate lower bound on the ratio 2ψ j /ψ j , whereψ j =ψ uj = exp(m Tα (u j )) is the entropy ansatz associated with an exact solution of (2.8) and ψ j = exp(m T α) is the ansatz associated with the Lagrange multiplier vector which satisfies (3.4). While of course we do not know the exact entropy ansatz, we can approximate the difference between the exact solutionα(u j ) and another multiplier vector α from the iterations of the optimizer by the Newton direction This leads to the approximate bound where m ∞ = max i,µ |m i (µ)| = 1. Finally, and exactly as in [3], we use a isotropic-regularization technique to return multipliers for nearby moments when the optimizer fails (for example, by reaching a maximum number of iterations or being unable to solve for the Newton direction). Isotropically regularized moments are defined by the convex combination v(u, r) : where u iso = 1 2 m is the moment vector of the normalized isotropic density φ(µ) ≡ 1/2. Then the optimizer moves through a sequence of values 0, r 1 , r 2 , . . . , r M , advancing in this sequence only if the optimizer fails to converge for v(u, r) after k reg iterations for the current value of r. It is assumed that r M is chosen large enough that the optimizer will always converge for v(u, r M ) for any realizable u.

3.2.
Spatial WENO reconstruction. We use the standard WENO reconstruction method given for example in [34,36]. For the unfamiliar reader, in this section we briefly introduce the method, while a more detailed documentation and demo implementations of the reconstruction procedures we implemented can be found on our webpage [1].
Here, as in the previous section, we use ψ j to indicate the entropy ansatz associated with the Lagrange multipliers returned by the optimizer to approximate the true multipliersα(u j ). Time dependence is again suppressed for clarity of exposition.
At x = x j−1/2 , the cell edge between the (j − 1)-th and j-th cells, for each µ we evaluate a weighted combination of polynomials of degree k − 1, p jm (·, µ) ∈ P k−1 (the space of polynomials up to degree k − 1), m = 0, 1, . . . , k, each solving the interpolation problem The WENO method then gives weights ω ± j−1/2,m to form the weighted averages and finally we approximate the values at the cell interfaces by The weights ω ± j−1/2,m are non-linear functions of the cell-averages and reflect the smoothness of each polynomial p jm . They are computed such that for smooth data the approximation order at the cell edge is maximized. This gives an order 2k − 1 approximation at the cell edge, while the overall order in the interior of the cell is k.
When at least one of the interaction coefficients σ a or σ s is spatially dependent, we must also specify the reconstruction inside each cell to compute, for example,

FLORIAN SCHNEIDER, JOCHEN KALL AND GRAHAM ALLDREDGE
the σ a u j term in (3.2). Here we must make a choice, because both p + j−1/2 and p − j+1/2 are order k reconstructions of the density ψ in the j-th cell. 3 We denote this polynomial ψ j (x, µ) and choose it to be This particular reconstruction allows us to derive a realizability-preserving time step in Theorem 4.1. Finally, it remains to incorporate boundary conditions. We define 'ghost cells' at the cell indices j ∈ {1 − k, . . . , 0, J + 1, . . . , J + k}, namely those indices which are used in (3.8) but have not yet been defined. We assume that we can smoothly extend ψ L and ψ R in µ to [−1, 1]. We then use the simplest possible approach 4 and set We note, however, that the validity of this approach is not entirely noncontroversial, but the question of appropriate boundary conditions for moment models is an open problem [22,23,29,35] which we do not explore here.
In case of periodic boundary conditions, we use the data from the physical cells j ∈ {J − k + 1, . . . , J} from the right side of the domain to fill the ghost cells j ∈ {1 − k, . . . , 0}. The ghost cells on the right side of the domain analogously take the values from the left side of the physical domain.
3.3. High-order time integration. If we collect the approximate cell means from each spatial cell into one long vector u h (t) := (u T 1 (t), u T 2 (t), . . . , u T J (t)) T , then equation (3.2) can be written as Since entropy-based moment closures are only defined on the realizable set, it is important to choose a time integrator for which we can prove that realizability is preserved. Therefore we follow [2,39] and use a strong stability-preserving (SSP) method whose stages and steps are convex combinations of forward Euler steps.
Since the realizable set is convex, the analysis of a forward Euler step then suffices to prove realizability preservation of the high-order method. When possible we use a SSP Runge-Kutta (SSP-RK) method, but such methods only exist up to order four [17,31]. For higher orders we use the so-called two-step Runge-Kutta (TSRK) SSP methods [20] as well as their generalizations, the multistep Runge-Kutta (MSRK) SSP methods [5]. They combine Runge-Kutta schemes with positive weights and high-order multistep-methods to achieve a total order which is higher than four while preserving the important SSP property. If we let Figure 1 -One possible startup procedure for SSP TSRK schemes. The first step from t0 to t1 is subdivided into substeps (here there are three substeps of sizes ∆t/4, ∆t/4, and ∆t/2). A one-step SSP Runge-Kutta scheme is used for the first substep, and subsequent substeps are taken with the TSRK scheme itself, doubling the step sizes until reaching t1. Illustration taken from [20].
the solution at the n-th time instant t n = n∆t, an s-stage TSRK method in the low-storage implementation has the following form [20]. For 0 ≤ ≤ s set and update using where the coefficients d , q m , ζ, η m , ρ define the scheme. In this work we only consider explicit schemes, where q m = 0 for ≥ m. The positive coefficient ρ is called the radius of absolute monotonicity and indicates how much we can scale our time step ∆t while fulfilling the CFL condition for forward Euler steps (see Section 4.1 below). Such schemes provide reasonably good effective CFL numbers, which are the ratios ρ/s for each method.
If the forward-Euler method is stable under the timestep restriction ∆t ≤ ∆t 0 , the high order scheme with radius of absolute monotonicity ρ is stable under the time-step restriction ∆t ≤ ρ∆t 0 . A larger effective CFL indicates that, in order to reach a given final time, the operator L h will need to be evaluated fewer times. The explicit Euler method, which serves as the reference for efficiency in this context, has an effective CFL condition of 1 (s = ρ = 1).
The effective CFL numbers of the integration schemes used here are given in Table 1.
Unfortunately, multistep methods are not self-starting, so they need a predictor for the first time-step. Since we can only use convex combinations of forward Euler steps for all time steps in order to prove that realizability is maintained we must use a lower-order method. We use the strategy given in [20]: First, we predict with a smaller step size ∆t = ∆t/2 q , for an integer q ≥ 1, with the ten-stage, fourth-order explicit SSPRK(1, 4, 10) method. Then we use the corresponding TSRK method and double the step-size after every iteration until we reach t = ∆t. This procedure is shown in Figure 1 for q = 2.
For the five-step method MSRK (5,7,12), which we use for seventh-order simulations, we have to initialize four steps. For these initialization steps we use the twostep method TSRK (2,7,12), whose initial step we predict using the same method given above. However the radius of absolute monotonicity ρ for TSRK (2,7,12 approximately 2.7659, while for MSRK (5,7,12), the radius of absolute monotonicity ρ is approximately 3.0886. This means the time steps we take after initialization will be longer than those we can take with the TSRK method during initializating without violating the realizability-preserving CFL condition. Therefore we stop increasing the step size in the initialization routine when we reach ∆t/2 (as opposed to ∆t), and then continue with the TSRK initialization steps of size ∆t/2 until we have computed u h (t 4 ) = u h (4∆t). At this point we have all the previous steps we need to compute u h (t 5 ) using the five-stage MSRK(5, 7, 12) method.

4.
Realizability preservation and limiting. The strategy to prove that the moments u j remain realizable at each time step and inner stage of the time integrator begins by proving that the kinetic density reconstructed for the cell means remains nonnegative after a forward Euler step with a certain time-step restriction. Since we use time integrators which are a convex combinations of Euler steps (see Section 3.3) this immediately gives nonnegativity for all time steps and internal stages of time integration. This proof, however, requires the assumption that the point-wise values of the current polynomial reconstruction ψ j are nonnegative at every spatial and angular quadrature point. One therefore introduces a limiter to enforce this nonnegativity.

4.1.
Realizability preservation of the cell means. We follow along the lines of the main proof in [2] and will provide weaker conditions which follow [37,39]. A spatial quadrature rule plays a crucial role here. We use Gauss-Lobatto rules which are exact for polynomials of degree 2Q − 3, where Q is the number of quadrature nodes. These rules are characterized by nodes y i and weightsŵ i for i ∈ {1, 2, . . . , Q} on the reference interval [− 1 2 , 1 2 ]: We let x ji := x j + ∆xy i denote the quadrature nodes shifted and scaled for cell I j , and note that since the Gauss-Lobatto rules include the endpoints, we have x j1 = x j−1/2 and x jQ = x j+1/2 . We also use the property that, since the collision kernel T in (2.2) is nonnegative, there exists a realizable moment vector u C such that (4.2) Theorem 4.1 (Main theorem). Assume that (i) for all cells j ∈ {1, 2, . . . , J} the coefficients 0 ≤ S(t n , x)| Ij , σ a (t n , x)| Ij , and σ s (t n , x)| Ij are in P k S −1 (I j ); (ii) the cell means u n j at time step t n are realizable; (iii) in each cell the ratio between the exact entropy ansatz and its approximation from the optimizer satisfiesψ for all µ ∈ [−1, 1] in the angular quadrature set; and (iv) the point-wise values of the polynomial reconstructionsψ j (x, µ) ∈ P k−1 at the quadrature nodes of the Q-point Gauss-Lobatto on each cell I j are nonnegative, for Q = (k + k S + 1)/2 . 5 Let σ t,max := max j∈{1,...,J} max i∈{1,...,Q} σ s (t n , x ji ) + σ a (t n , x ji ). (

4.4)
Then under the CFL condition the cell means u n+1 j after one forward Euler step are realizable.
Proof. For simplicity we will neglect the time-index for quantities from t = t n and use it only for time-level t = t n+1 . An Euler step is given by ∆t −σ a u j + σ s r(u) j + s j (4.6b) and consequently we have u n+1 j = mφ j for φ j given by where the total interaction coefficient is defined as σ t := σ a + σ s , and ψ C ≥ 0 is the entropy ansatz corresponding to the realizable part u C of the collision operator (see (4.2)). Note that φ j = φ j (µ) depends on µ. Let us first consider the case µ > 0. Stripping away positive terms and using µ ≤ 1 gives Next, we want to use our polynomial reconstruction (3.9) and the Q-point Gauss-Lobatto quadrature to relate the exact cell meanψ j to the reconstruction's values at the cell edge ψ − j+1/2 and the cell average σ t ψ j . But since only the approximations from the optimization are used to define the reconstruction, we must first multiply and divideψ j by ψ j . Note that this is possible since ψ j > 0. Then we use the assumed bound (4.3), note carefully that for µ > 0 indeed ψ − j+1/2 = ψ j (x jQ )), and apply L ∞ bounds on the total interaction coefficient. This gives where in the last line we have introduced the notationψ ji :=ψ j (x ji ). One can see that (4.5) ensures nonnegativity of both terms in the final expression. Recalling thatŵ 1 =ŵ Q , the cases µ < 0 and µ = 0 follow analogously, and together we have that φ j ≥ 0 at each µ, which shows that u n+1 j is realizable.

Remark 1. In practice, (4.3) is only approximately enforced by (3.4). However, we
have never had problems losing realizability in our numerical simulations, including those that approach the boundary of realizability.

Remark 2.
We have assumed that the reconstruction (3.9) is always used to compute σ t ψ j . However, when σ a and σ s are constant in each cell, one simply has σ t ψ j = σ tψj , and therefore the reconstruction is not needed. The proof can then be completed with an upper bound onψ j /ψ j , which is similarly easy to enforce approximately, and one ends up with a slightly less restrictive CFL condition. However, we did not use this in our implementation.
The time-step restriction (4.5) can then be scaled by ρ for the corresponding time-integration scheme to give realizability of every stage and step in the scheme.
In practice, to achieve an order k method for sources S or interaction coefficients σ a or σ s which are not piecewise degree k − 1 polynomials, one would approximate them using the same spatial reconstruction techniques that we use for the density to achieve an order k approximation of the corresponding terms. Thus in Theorem 4.1, one would not use a value of k S larger than k.

4.2.
Limiting. The first role of the limiter, then, is first to ensure that the pointwise values of the polynomial reconstructionsψ j are nonnegative at the spatial and angular quadrature points. However, as we will see in Section 5.3 below, the numerical solutions using a limiter which only ensures nonnegativity can still contain spurious oscillations. Therefore we extend the same limiter to enforce local maximum principles as well, thereby much more effectively dampening such oscillations.

Positivity-preserving limiter.
To preserve nonnegativity we can simply apply a linear scaling limiter. The limited spatial reconstruction is defined as ψ θ j (x) := (1−θ)ψ j +θψ j (x); notice that θ = 1 corresponds to no limiting. For each quadrature point (both in space and angle, though here we suppress the angular argument) we Then in each cell we set θ = θ j := min i=1,...,Q {θ ji } (where one should keep in mind that θ j still depends on µ). One immediately sees that this limiter ensures the positivity, preserves the cell means ψ j , and following arguments from [37,38], does not destroy accuracy of the scheme if ψ j > 0. However, it has been remarked in [41] that in some pathological situations this limiter may reduce the accuracy to second order. See also [4] where a similar observation has been made for a realizability limiter in a discontinuous-Galerkin scheme.

4.2.2.
Maximum principle-satisfying limiter. The limiter we introduce here is a slightly modified version of the maximum-principle limiter from [40].
Since we know a priori that ψ satisfies a strict maximum principle m ≤ ψ(x, µ) ≤ M for all x and µ, a natural strategy to dampen artificial oscillations in numerical solutions is to enforce a local maximum principle. Specifically, we would like the polynomial reconstruction ψ j (x) to be bounded by the data of those cells which influence it. The corresponding index set of influential nodes is N j,k = {j−k, . . . , j+ k} (cf. (3.8)), so the local maximum principle we would like to enforce is However, this tends to flatten smooth extrema, so, inspired by the modified minmod function in [8], we relax the strict maximum principle by setting the maximum principle bounds M j and m j locally as where c is a local bound on the relative derivative of ψ, i.e.
and the maximum and minimum in µ are taken over the angular quadrature nodes. 6 Therefore the maximum principle that we will actually enforce is for all spatial quadrature points x ji ∈ I j and all angular quadrature points. To enforce this maximum principle, for each spatial quadrature point we set 1 otherwise, (4.10) and finally for each cell we choose θ j := min i {θ ji }.

FLORIAN SCHNEIDER, JOCHEN KALL AND GRAHAM ALLDREDGE
As with the positivity-preserving limiter, it can be shown that this limiter does not destroy accuracy [38]. 5. Numerical results. In this section we present results to confirm that our scheme converges with the expected order, to show the effect of various parameters in the scheme, and to highlight some of the features of high-order solutions.
Except where otherwise noted, we used the following parameter values: τ = 10 −9 Optimization gradient tolerance, ε = 0.01 Optimization tolerance on 1 −ψ j /ψ j , For the value q determining the number of initialization steps for the multi-step time integrators we used two for fifth-order simulations and three for sixth-and seventh-order simulations.
The time step is chosen to fulfill (4.5) (replacing ∆t by ∆t/ρ for the appropriate time integrator) with equality.
In both benchmark problems we use isotropic scattering, C(ψ) = 1 2 ψ − ψ. For the first stage of the first time step, the initial multipliers for the optimizer are those associated with the normalized isotropic distribution. For the following stages and time steps, the initial multipliers are set to those from the previous stage or step at the same spatial cell. However, if in these later stages the optimizer cannot converge before initializing the regularization loop, we first switch back to the multipliers associated with the normalized isotropic distribution and restart the optimizer (this time allowing regularization if necessary). In our experience, the isotropic multipliers are the safest choice for the initial condition, and therefore this technique reduces the number of times the regularization must be used.
As in [2], if regularization has to be applied to u j , we replace it (in (4.6)) by the regularized moments v(u j , r) (3.7) for which the optimizer converged. Then Theorem 4.1 can be applied to ensure that the next iterate will be realizable as well.
In all figures below we plot the zeroth-order component of the reconstruction u j (x) = mψ j (x, µ) for x ∈ I j , see (3.9).

M N manufactured solution.
In general analytical solutions for minimumentropy models are not known. Therefore, to test the convergence and efficiency of our scheme, we use the method of manufactured solutions, and we follow the target solution given in [4] but add a spatially and temporally dependent absorption interaction coefficient. The solution is defined on the spatial domain X = (−π, π) with periodic boundary conditions. A kinetic density in the form of the entropy ansatz is given by A source term is defined by applying the transport operator to φ: where σ a (t, x) := 4 − 4 cos(x − t)). Thus by inserting this S into (2.1) (and taking σ s = 0) we have that φ is, by construction, a solution of (2.1).
A straightforward computation shows that S ≥ 0 (for any a or K), which means that Theorem 4.1 will apply to the resulting moment system. Furthermore we take so that the maximum value of φ for (t, x) ∈ [0, t f ] × X is one. As K is increased, φ converges to a Dirac delta at µ = 1.
The moment vector w = mφ is then a solution of (2.4) for M N models with N ≥ 1.
We used the final time t f = π/5 and chose K = 4, for which w 1 /w 0 ∈ [0.67, 0.8] (recall that |w 1 /w 0 | < 1 is necessary for realizability). We are using a fairly low value of K because, in order to show convergence for our scheme with the highestorder (k = 7), we need to be able to have a tighter control on the errors from the numerical optimization. When K is higher, these errors are larger and drown out the convergence in space and time. In the following, we used the M 3 model so that our results included the effects of the numerical optimization.
We compute errors in the zero-th moment of the solution, which we denote w 0 (t, x) = φ(t, x, ·) . Then L 1 and L ∞ errors for u 0,h (t, x) (that is, the zero-th component of a numerical solution u h ) are defined as respectively. We approximate u 0,h (t f , x) using the same reconstruction technique as for the scheme for the underlying kinetic density (3.9) and integrate with respect to µ. Then we approximate the integral in E 1 h using a 100-point Gauss-Lobatto quadrature rule over each spatial cell I j , and E ∞ h is approximated by taking the maximum over these quadrature nodes.
The observed convergence order ν is defined by where for i ∈ {1, 2}, E p hi is the error E p h for the numerical solution using cell size ∆x i , for p ∈ {1, ∞}.
Convergence tables are given in Table 2 for solutions with a tighter optimization tolerance of τ = 10 −11 . We observe that the scheme converges with at least its designed order until the errors are roughly O(τ ), where errors from the numerical optimization halt the convergence. For many of the solutions on the coarsest grids, the convergence is faster than designed, likely because the WENO reconstruction is order 2k − 1 at the cell interfaces. The effects are indeed more pronounced for higher orders.
In Figure 2 we plot the error of solutions of various orders against their computation time. Here we confirm the expectation that for smaller errors, higher-order solutions require less computation time.   We set σ s = 1 and σ a = 0. All solutions here are computed with an even number of cells, so the delta function in the initial condition lies on a cell boundary. Therefore we approximate the delta function by splitting it into the cells immediately to the left and right.
Since this is a highly non-smooth problem, the choice of c is not immediately obvious. Indeed, the solution contains such strong gradients, so choosing c too small can greatly reduce the accuracy of the simulation. After some numerical experimentation, we observed that in the 1000-cell seventh-order simulations the value c = 15 generally gave a good trade-off between smoothness in the solution without adding too much diffusivity. In the 4000-cell fourth-order simulations, where the non-smoothness is further resolved, c = 10 gave the best results.
In Figure 3 and Figure 4, we plot several solutions at the final time for the M 3 and M 7 models, respectively, as well as a reference solution 7 of the kinetic equation (2.1). We consider numerical solutions of orders k ∈ {1, 4, 7}. The first-order solutions are included to indicate our best guess of the true M N solutions, and they largely agree with those presented in [18] (although we have computed solutions on a much finer grid in order to better resolve sharp peaks in the solutions).
In Figure 3, we present seventh-order M 3 solutions. In Figure 3a, we see that while the seventh-order solution using only the positivity-preserving limiter closely matches the highly-resolved first-order solution, some spurious oscillations around x = ±0.25 remain. These oscillations are largely removed by using the maximumprinciple limiter with c = 15, as we see in Figure 3b, though in Figure 3c we see that this solution is also somewhat more diffusive.
The M 7 solution to the plane-source problem is more oscillatory. Figure 4a shows that our kinetic scheme with fourth-order reconstructions and the maximumprinciple limiter performs well. However, seventh-order solutions are notably more diffusive, as shown at the leading edge of the solution in Figure 4b.
That the high-order solutions are more diffusive than the true M N solutions is not necessarily a disadvantage: the more diffusive solutions are actually closer to the reference kinetic solution. Indeed, the M N equations are a spectral method for the original kinetic equation, and thus should not be applied to a non-smooth problem like this one without filtering to avoid the Gibbs phenomenon. Exactly how to apply such a filter for moment models of kinetic equations is a topic of ongoing research [26,30], but here it seems that the maximum-principle limiter already filters the solution somewhat. 5.3. Source-beam. Finally we present a discontinuous version of the source-beam problem from [14]. The spatial domain is X = [0, 3], and  Figure 5a using J = 150 cells and seventh-order reconstructions, using the maximum-principle limiter with c = 1, along with a reference solution. 8 We see that increasing the moment order to N = 3 qualitatively improves the solution significantly.   Figure 5b shows some strengths and weaknesses of our scheme for the M 1 model, where for this problem the shock is the strongest. When compared to a much more finely resolved low-order approximation, the seventh-order approximation on a coarse grid nicely fits the features of the solution despite the discontinuous physical parameters. However, the incoming beam can only be resolved up to first order (note the piecewise-constant approximations behind the beam), resulting in a more diffusive solution at the shock around x = 1.
The source-beam problem is particularly well-suited to test the oscillation-dampening effects of the limiters. In Figure 6 we compare several seventh-order solutions. First in Figure 6a we see in the M 1 solution that the positivity-preserving limiter does not dampen spurious oscillations while the maximum principle-preserving limiter does a better job even for a relatively large value of c. Second, in Figure 6b, we show with the M 2 model that for small values of c, the maximum principlepreserving limiter is too diffusive. 6. Conclusions and outlook. In this paper we describe how to implement a kinetic scheme of (in principle) arbitrarily high order for entropy-based moment closures of linear kinetic equations in one space dimension. For spatial reconstructions we use the well-known WENO method to reconstruct the underlying entropy ansätze using interpolating polynomials, and time integration is performed using multi-step SSP methods. These SSP time integrators play a key role in allowing us to give a time-step restriction which guarantees that the moments stay in the realizable set. The other key component is a limiter, which not only ensures positivity of the polynomial reconstructions on a spatial quadrature set, but also enforces a local maximum principle which dampens spurious oscillations in numerical solutions.
We performed convergence tests with a manufactured solution that included the effects of a space-and time-dependent absorption interaction coefficient, and these  results validated that the scheme is converging at least as fast as expected, and often faster at lower resolutions with higher orders. The convergence tests also showed that errors from the numerical optimization routine needed for the angular reconstructions limit the overall accuracy of the scheme. This indeed eliminates the benefit of going beyond a certain order (depending on the optimization tolerance τ ). However, our manufactured-solution tests showed that increases in efficiency can be obtained before the optimization errors dominate the solution.
Using the plane-source, a challenging highly non-smooth benchmark problem, we showed that with the maximum-principle limiter accurate solutions can be obtained which even limit some of the spurious oscillations due to Gibbs phenomena in the true M N solutions, thus pushing our solutions closer to the kinetic solution. However, the choice of the parameter c plays an important role in the approximation quality of the scheme, and a good value is not always available a priori. With another benchmark problem, the source beam, we also demonstrated the benefit of using the maximum principle-preserving limiter. Again, one must strike a balance between flattening smooth extrema and dampening spurious oscillations with the choice of its parameter.
Compared to the discontinuous-Galerkin implementation in [4], where realizability preservation is complex due to the structure of the set of realizable moments, realizability preservation on the kinetic level is much simpler. However, the wider stencils in the WENO reconstruction process will influence the overall parallelizability of the scheme.
Future work should continue to work towards practical implementations of entropy-based moment closures. Models in two and three spatial dimensions should be implemented, and here a notable challenge is the increasing number of angular quadrature points that will be needed. Indeed, our reconstructions are performed at every angular quadrature point, so more efficient WENO techniques will be necessary. Other collision models should also be considered. The Laplace-Beltrami operator in the Fokker-Planck equation does not fall under the types of collision operators considered here but is an important model for problems with forward-peaked scattering. This appears, for example, in important applications such as radiotherapy. Finally, since we have only considered explicit time-stepping schemes, our time-step restriction scales with the mean-free path. Implicit-explicit or asymptotic-preserving schemes should be developed to handle moment models near diffusive or fluid regimes without requiring extremely small time steps.