{'fraction_non_alphanumeric': 0.047284083943842856, 'fraction_numerical': 0.017629374507457458, 'mean_word_length': 4.447924948322468, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
thornado-transport: IMEX schemes for two-moment neutrino transport respecting Fermi-Dirac statistics

We develop implicit-explicit (IMEX) schemes for neutrino transport in a background material in the context of a two-moment model that evolves the angular moments of a neutrino phase-space distribution function. Considering the upper and lower bounds that are introduced by Pauli’s exclusion principle on the moments, an algebraic moment closure based on Fermi-Dirac statistics and a convex-invariant time integrator both are demanded. A finite-volume/first-order discontinuous Galerkin(DG) method is used to illustrate how an algebraic moment closure based on Fermi-Dirac statistics is needed to satisfy the bounds. Several algebraic closures are compared with these bounds in mind, and the Cernohorsky and Bludman closure, which satisfies the bounds, is chosen for our IMEX schemes. For the convex-invariant time integrator, two IMEX schemes named PD-ARS have been proposed. PD-ARS denotes a convex-invariant IMEX Runge-Kutta scheme that is high-order accurate in the streaming limit, and works well in the diffusion limit. Our two PD-ARS schemes use second-and third-order, explicit, strong-stability-preserving Runge-Kutta methods as their explicit part, respectively, and therefore are second-and third-order accurate in the streaming limit, respectively. The accuracy and convex-invariance of our PD-ARS schemes are demonstrated in the numerical tests with a third-order DG method for spatial discretization and a simple Lax-Friedrichs flux. The method has been implemented in our high-order neutrino-radiation hydrodynamics (thornado) toolkit. We show preliminary results employing tabulated neutrino opacities.


Introduction
Core-collapse supernovae (CCSNe) are the explosions of massive stars that end their lives. They are directly or indirectly responsible for the lion's share of elements heavier than oxygen and play important roles in many astrophysical phenomena, such as neutron star and black hole formation. Furthermore, these explosions occur at energies and densities relevant to address fundamental questions in nuclear, particle, and gravitational physics. A solid theoretical framework for the CCSN explosion mechanism may help answer important questions in fundamental physics [1].
One essential part of the explosion mechanism is neutrino transport. Neutrino energy deposition is believed to be the major driver of CCSN explosions, except in peculiar cases where rapid rotation is present and magnetohydrodynamic effects may dominate (for reviews, see [2,3,4,5]). Ideally, neutrino transport would be modeled by the Boltzmann transport equation, which is an integro-partial-differential equation evolving a phase-space distribution function f (e.g., see [6,7,8,9,10,11,12,13,14,15,16,17,18] for studies of CCSN with Boltzmann transport in various approximate settings). Simulating the neutrino transport implies finding a solution of the Boltzmann equation for a specific domain and period, with acceptable accuracy.
However, solving the Boltzmann transport equation with sufficient phase-space resolution and full weak interaction physics is at present too expensive. To balance physical fidelity and computational expediency, an approximate method called the two-moment method has been adopted (e.g., see [19,20,21,22]). Using the two-moment method, the evolved variables are the zeroth and first angular moments of the distribution function f -the spectral particle density J and flux H, respectively. However, the equation of H includes the second angular moment K. Knowledge of K is needed to close the two-moment system. Therefore, a closure that gives K consistent with J and H is needed. The better the closure predicts K, the more accurate the two-moment method will be. The two-moment method has been widely applied in the CCSN modeling community with different algebraic closures, such as the Minerbo [23] closure (e.g. O'Connor and Couch [24], Pan and et al. [25], Glas et al. [26], and Just et al. [21]) and the Levermore [27] closure (e.g. Vartanyan et al. [22], Cabezon et al. [28], and Kuroda et al. [19]).
Applying the two-moment method does simplify the problem, but doesn't guarantee an affordable solution. To be precise, how to discretize the continuous system of equatio ns given by the two-moment method and solve the discretized system efficiently remains a question. In fact, the time scales of neutrino interactions with the background (they can be ∼ O(10 −13 ) second) is short compared to the duration of the CCSN explosion (∼ O(1) second). This means that ∼ O(10 13 ) time steps could be needed for solving the system fully explicitly. On the other hand, solving the moment equations fully implicitly requires inverting global band-structured matrices whose sizes depend on the phase-space discretization. Such a global inversion is both expensive and unfriendly to parallelization. To circumvent these challenges, implicit-explicit (IMEX) methods are taken into consideration. By treating the transport terms in the two-ASTRONUM IOP Conf. Series: Journal of Physics: Conf. Series 1225 (2019) 012013 IOP Publishing doi: 10.1088/1742-6596/1225/1/012013 3 moment equations explicitly and the collision terms implicitly, IMEX methods are subject only to a time step governed by the explicit transport terms, and the matrices to be inverted are block diagonal. Therefore, IMEX methods require far fewer time steps compared with a fully explicit method, and the computation for each step is easily parallelizable. For the relativistic setting that we have, where the fluid and the neutrinos have comparable propagation speeds, IMEX methods can be efficient.
To model neutrino transport using a two-moment method, two (or at least two) things need to be chosen carefully: an algebraic closure based on Fermi-Dirac statistics for closing the twomoment equations and a convex-invariant, diffusion-accurate IMEX scheme to ensure a physical result. A convex-invariant scheme has the following property: if the solution u n ∈ W and W is a convex set, then u n+1 ∈ W . Since the neutrino distribution function is bounded (f ∈ [0, 1]) by the Pauli exclusion principle, its moments as weighted integrals of a bounded function over the domain ω ∈ S 2 are also bounded. We call the moments satisfying the constraints due to Fermi-Dirac statistics realizable moments and the realizable J and H define a convex set [29]. The algebraic closure should give a realizable K, and the well-posedness of the closure requires realizable J and H. This explains why an algebraic closure based on Fermi-Dirac statistics is needed. Realizability of J and H after each time step requires a convex-invariant IMEX scheme. Since the realizable moments form a convex set, it is possible to construct a realizability-preserving method using a convex-invariant IMEX scheme for two-moment neutrino transport. In addition, the physics of neutrino transport in CCSNe requires the IMEX scheme to be diffusion-accurate.
The study of moment realizability and realizability-preserving methods with diffusionaccurate IMEX schemes motivates this work. Gottlieb et al. [30] showed that standard strongstability-preserving IMEX schemes cannot have an order higher than first without a restricted time step. One way to obtain the second-order (or higher-order) accuracy is to add some correction steps after the standard step [31,32]. Unfortunately, the correction steps can deteriorate the accuracy of the IMEX scheme in the diffusion limit or restrict the time step. To keep things simple, we focus on IMEX schemes without correction steps and require them to be high-order (second or higher order) in the streaming limit and diffusion-accurate. We call these IMEX schemes PD-ARS.
thornado is our toolkit for high-order neutrino-radiation hydrodynamics based on high-order Runge-Kutta Discontinuous Galerkin (RKDG) methods. It is being developed at the University of Tennessee, Knoxville and Oak Ridge National Laboratory. It currently includes solvers for the Euler equations for fluid dynamics and the two-moment approximation of the radiative transfer equation [33]. In this paper, we focus on the transport methods in thornado with emphasis on IMEX.
This paper is organized as follows: Section 2 discusses the mathematical model, algebraic closures, and the constraints on the moments and algebraic closures imposed by Fermi-Dirac statistics; Section 3 gives a first-order finite-volume spatial discretization and shows how the spatial discretization preserves constraints in an IMEX step; Section 4 discusses how to use convex combination to construct two PD-ARS schemes, one with second-order accuracy in the streaming limit and the other with third-order accuracy in the same limit; Section 5 presents the results of the numerical tests, which demonstrate the properties of the PD-ARS schemes; Section 6 summarizes the achievements of this paper and discusses future work.

Transport Equations
We consider neutrino transport through a static background in a non-relativistic framework. For simplicity, we only include neutrino-matter interactions due to emission, absorption, and isoenergetic scattering. It is possible to include other interactions, such as neutrino-electron scattering and pair processes, in the two-moment model. However, these interactions result in dense coupling across neutrino energies and across neutrino species, and we have not yet considered the realizability-preserving property of our scheme in this case. We leave this for future work.
After scaling to dimensionless units, the Boltzmann equation can be written as where the distribution function f = f (ω, ε, x, t) gives the number of neutrinos propagating in the direction ω ∈ S 2 , with energy ε ∈ R + , at position x ∈ R 3 and time t ∈ R + . (ω) ∈ R 3 is the unit vector parallel to the neutrino three-momentum: p = ε . On the right-hand side, τ is a collision time scale. In opaque regions, where neutrinos have frequent interactions with the background, τ 1. In transparent regions, where neutrinos rarely interact and stream freely, τ 1. The collision term, C(f ), which models emission, absorption, and isoenergetic scattering is given by where ξ = σ A /(σ A + σ S ) is the ratio of the absorption opacity σ A to the total opacity. The scattering opacity is σ S . The limit ξ = 1, when σ S = 0, corresponds to pure absorption, while ξ = 0, when σ A = 0, corresponds to pure scattering. The equilibrium distribution function for neutrinos is given by the Fermi-Dirac distribution: where z := {ε, x}, T is the material temperature in energy units and µ is the neutrino chemical potential. Both T and µ depend on the spatial position x.

Two-Moment Model
Approximate solutions to the Boltzmann equation, Eq. (1), can be found by solving the twomoment model. To this end, define the angular moments of the distribution function as follows: The zeroth moment, J , is referred to as the particle density. The first moment, H, is the particle flux, and the second moment, K, is proportional to the stress tensor. By integrating Eq. (1) over the momentum-space angular dimension we obtain equations for the zeroth and the first moments: with M = (J , H) T , F = (H, K) T , and where η = (ξ f 0 , 0) T , D = diag(ξ, I), and I is the identity matrix. Hence, the process of solving the Boltzmann equation, Eq. (1), for the neutrino distribution function f (ω, z, t), is replaced by solving the two-moment equations for the neutrino number density, J (z, t), and flux, H(z, t).

Algebraic Closures
The moment equation for H involves the higher moment K and the two-moment model is open. To close the two-moment model, we consider algebraic closures. For the two-moment model, algebraic closures give an approximation to K using the lower moments: where k is the Eddington tensor. By assuming that the distribution function is symmetric about a preferred direction h = H/|H|, Levermore [27] proposed a simple form for the Eddington tensor: where χ = χ(J , |H|) is the Eddington factor. Thus, the two-moment model is closed by specifying the scalar χ in terms of J and |H|.

Constraints on the Moments
Neutrinos are fermions and obey the Pauli exclusion principle. Because of this, the neutrino distribution function is bounded; i.e. f ∈ [0, 1]. As a result, the angular moments J and H and the Eddington factor χ satisfy the following bounds [27,34,35,36]: where h = |H|/J is the flux factor. The inequalities in Eqs. (9) and (10) hold for moments as defined in Eq. (4) with f ∈ [0, 1]. In a relativistic model (e.g., [37,38]), the evolved quantities -functions of the moments in Eq. (4) -satisfy different bounds. We are currently exploring extensions of the realizability-preserving scheme presented here to the relativistic case. The constraints in Eq. (9) define realizable moments M. For fermions, realizable moments can only be constructed from a distribution satisfying the bounds f ∈ [0, 1]. Moreover, the set of realizable moments is convex: let R be the realizability set and M 1 , M 2 ∈ R, then λM 1 + (1 − λ)M 2 ∈ R for any λ ∈ [0, 1] [29]. As we will see later in Section 3, this convexity makes it possible to design a realizability-preserving discretization for solving the two-moment model numerically.
The inequalities in Eq. (10) deserve further attention. They are as important as those in Eq. (9), in maintaining consistency of the two-moment model with respect to Fermi-Dirac statistics. When designing a numerical scheme for the two-moment model that maintains realizable M, it is also necessary for the Eddington factor to satisfy the bounds in Eq. (10).
However, recently reported CCSN simulations using two-moment neutrino transport with algebraic closures have employed Eddington factors that can violate the bounds in Eq. (10). As examples, we consider the Eddington factors discussed in [39], where the suitability of several algebraic closures for two-moment neutrino transport was evaluated. In Fig. 1, we plot the Eddington factor versus the flux factor for two occupancies: J = 0.1 (low occupancy) and J = 0.9 (high occupancy). Of the algebraic closures plotted, few satisfy the bounds on the Eddington factor in Eq. (10). Kershaw [35], Wilson [40,41], Levermore [27], Minerbo [23], and Janka 2 [42] closures may work fine when the occupancy is low. When the occupancy is high, the Eddington factor due to these closures exceeds the upper bound for Fermi-Dirac statistics. The Eddington factor of Janka 1 [43] [45] use the Minerbo closure, while Vartanyan et al. [22], Cabezon et al. [28], Kuroda et al. [19], and Fujibayashi et al. [46] use the Levermore closure. When employing these closures in conditions of high occupancy, a numerical scheme may evolve the moments outside the realizable domain of Fermi-Dirac statistics given in Eq. (9). If this were to happen, the update step may give J > 1. Considering the fact that the collision term contains blocking factors -i.e. (1 − J )× something positive -J > 1 would change the blocking factors' sign, and it would be difficult to predict the impact of the subsequent induced errors on the simulation outcome. Besides, even with the simplified collision term Eq 6, J > 1 could result in information loss (such as loss to the neutrino distribution for the groups that exceed the bounds), closure failure, and other unpredictable consequences. If treatments are developed to map the unrealizable moments into the realizable domain, they should be developed to conserve lepton number, energy, and momentum.

Spatial Discretization
We discretize the two-moment model with a simple first-order finite volume method to illustrate how the closure affects the realizability-preserving property of the method. By assuming that the moments at time level t n (M n ) satisfy the bounds in Eq. (9), our goal is to identify sufficient conditions such that the moments at time level t n+1 (M n+1 ) also satisfy the bounds. To simplify, we limit the discussion to one spatial dimension and employ a uniform Cartesian mesh. (The extension to multiple spatial dimensions and high-order discretization using the discontinuous Galerkin method is given in [29]. ) We divide the spatial domain D into N uniform cells and denote the i-th cell by K i , with i = 1, . . . , N ; i.e., Integrating Eq. (5) over each cell K i gives where F (M a , M b ) is the numerical flux and 1 τ C(M i ) is the collision term evaluated with M i . In this paper we use the global Lax-Friedrichs flux (setting the largest absolute eigenvalue of the flux Jacobian to one): By treating the transport term explicitly with forward Euler and the collision term implicitly with backward Euler, we have where we have defined Considering Eqs. (14) and (15) and assuming that M n i is realizable for all i, it can be shown (Lemma 3 in [29]) that M n+1 i is realizable provided ∆t τ > 0. In Eq. (15), if β ∈ [0, 1], M n i is expressed as a convex combination of M n i and the expression in the square brackets on the right-hand side of Eq. (15). It follows that M n i is realizable if the expression inside the square brackets is realizable. It can be shown (Lemma 2 in [29]) that the expression in square brackets is realizable for a distribution satisfying f ∈ [0, 1]. For the two-moment model considered here, realizability depends on the algebraic closure. Specifically, if M n i is realizable and the Eddington factor satisfies the bounds in Eq.

Time Integration
Suppose that an algebraic closure based on Fermi-Dirac statistics is used (i.e., the Eddington factor satisfies Eq. (10)). Here we consider the construction of an Implicit-Explicit (IMEX) time integration scheme that maintains the bounds in Eq. (9). The semi-discretization of the two-moment model results in a system of ordinary differential equations of the forṁ where the solution vector u(t) = (M 1 (t), . . . , M N (t)) T is the collection of all cell-averaged moments, T is the transport operator, corresponding to the first term on the right-hand side of Eq. (12), and Q is the collision operator, corresponding to the second term on the right-hand side of Eq. (12).
Since the set of realizable moments is convex, convex-invariant schemes, which map the initial values into this set, can be used to design realizability-preserving methods for the twomoment model. Ideally, the scheme should also be high-order accurate and work well in the asymptotic diffusion limit (characterized by frequent collisions and long time scales). The following discussion considers the construction of such convex-invariant schemes.

Standard IMEX Schemes
Treating the transport operator explicitly and the collision operator implicitly, a standard s-stage IMEX scheme takes the following form [47]: where (ã ij ) and (a ij ), coefficients of the i-th stage, are elements of matricesÃ and A, respectively. The matricesÃ and A are lower triangular. (Ã is strictly lower triangular so that the transport part is explicit.) The vectorsw = (w 1 , . . . ,w s ) T and w = (w 1 , . . . , w s ) T are the weights in the assembly step in Eq. (19). These coefficients and weights must satisfy certain order conditions for consistency, accuracy, and other properties. For second-order temporal accuracy, the following conditions are required [48]: and wherec i = s j=1ã ij and c i = s j=1 a ij . The IMEX scheme is called globally stiffly accurate (GSA) if the coefficients satisfy [49]: Then, u n+1 = u (s) , which is simplifying because the assembly step in Eq. (19) is omitted. IMEX schemes are further classified by the structure of the implicit matrix A. If A is invertible, the IMEX scheme is of type A [47]. If a i1 = 0 for i = 1, . . . , s, w 1 = 0, and the submatrix consisting of the last s − 1 rows and columns is invertible, the IMEX scheme is of type ARS [50,47].

Convex-Invariant IMEX Schemes
To be convex-invariant, the coefficients and weights defining the IMEX scheme must satisfy additional constraints. Our goal is to find constraints on a ij ,ã ij ,w i , and w i that enable each u (i) in Eq. (18) to be expressed as a convex combination of realizable states. Following Hu et al. [32], the stage values in Eq. (18) can be rewritten as c ij u (j) +ĉ ij ∆t T (u (j) ) + a ii ∆t Q(u (i) ), i = 1, . . . , s, where c ij andĉ ij ≡c ij /c ij are defined in terms of a ij andã ij . For IMEX schemes of type ARS, c ij andc ij are given by [32] Note that c i1 =c i1 = 0 in Eq. (24), so that i−1 j=0 c ij = 1. If the IMEX scheme is GSA, u n+1 = u (s) . Moreover, if c ij ,c ij ≥ 0 and a ii > 0, each stage in Eq. (23) is a convex combination of explicit Euler steps (with time stepĉ ij ∆t), followed by an implicit Euler step. Each of the explicit Euler steps has a time step condition that ensures its realizability given byĉ ij ∆t ≤ ∆x; the CFL condition of the scheme. Using results proved in [29] and discussed in Section 3, the IMEX scheme is convex-invariant and realizability-preserving for the two-moment model in Section 3 provided (This CFL condition becomes more restrictive with high-order DG spatial discretization [29].)

Diffusion Accurate, Convex-Invariant IMEX Schemes
Accuracy in the diffusion limit is another important property to consider when an IMEX scheme is applied to the two-moment model. In the diffusion limit, the distribution function is nearly isotropic, so K ≈ 1 3 J I and H ≈ − 1 3 τ ∇J , and the two-moment model is approximately governed by (e.g., [51]) In the context of IMEX schemes, the above relationships imply that the following relations should hold [29]: where e i is the ith column of the s × s identity matrix, e is the vector of ones, and A andÃ are the matrices of the coefficients (ã ij ) and (a ij ). Eq. (28) implies: We have proved in [29] that only IMEX schemes of type ARS can be both diffusion accurate and convex-invariant. (Another short proof follows from the fact that IMEX schemes of type A havec 1 = 0 while c i = 0.)  10 However, the correction steps can impose time step constraints for realizability or accuracy in the diffusion limit that ruin the efficiency gains expected from the IMEX scheme. Because of this, we sacrifice overall high-order accuracy, and aim for IMEX schemes that are high-order accurate in the streaming limit, diffusion accurate, and convex-invariant. Combining these requirements we seek GSA IMEX schemes of type ARS with coefficients satisfying the following constraints [32,29]:

PD-ARS IMEX schemes
(i) Consistency of the implicit coefficients: (ii) High-order accuracy in the streaming limit. For second-order accuracy: For third-order accuracy: (iii) Diffusion accuracy: (iv) Convex-invariance: Fortunately, these IMEX schemes are easy to find. (The constraint in (v) is introduced from efficiency considerations to limit the number of implicit solves.) We call the IMEX schemes satisfying the above conditions PD-ARS (see also Definition 3 in [29]), and we provide two optimal PD-ARS schemes below: PD-ARS2 and PD-ARS3, each limiting to the optimal secondorder and third-order SSPRK schemes from [53], respectively. Its explicit tableau is SSPRK3. This scheme requires three implicit solves per time step, and c Sch = 1. Since PD-ARS3 is not more accurate than PD-ARS2 in collision-dominated regions (see our results in Section 5), it may not offer any practical advantage over PD-ARS2.

Numerical Tests
In this section we present numerical results obtained with the PD-ARS schemes in Section 4.
The tests in Section 5.1 are designed to compare the accuracy of the schemes in streaming, absorption, and scattering-dominated regimes in one spatial dimension. The test in Section 5.2 demonstrates the convex-invariance of PD-ARS schemes. All the tests in this section were computed with third-order accurate spatial discretization (polynomials of degree k = 2) and time step ∆t = 0.1 × ∆x, using the DG scheme from [29].

Accuracy Tests
To compare the accuracy of the IMEX schemes, we applied our PD-ARS schemes and the SSP2332 scheme from Pareschi & Russo [47] to problems with known smooth solutions in streaming, absorption (damping), and scattering-dominated (diffusing) regimes in one spatial dimension. All the tests in this subsection were computed with the maximum entropy closure in the low-occupancy limit (i.e., the Minerbo closure). In the streaming test, the second-and thirdorder accurate explicit strong-stability-preserving Runge-Kutta methods from [30] (SSPRK2 and SSPRK3, respectively) are also included. To compare the numerical results to analytic solutions, the averaged absolute error or the averaged relative error are computed in the L 1 -error norm. We compute the absolute error for the streaming and diffusion tests and the relative error for the damping test. They are averaged over the cell with an equal-weight quadrature for the cell integrals. To examine the convergence, we let the number of elements (N ) vary from 8 to 128.

Sine Wave Streaming
The sine wave streaming test is designed to test accuracy in the free-steaming regime; i.e. σ A = σ S = 0. A periodic domain of unit length is used and the initial condition is J 0 = H 0 = 0.5 + 0.49 × sin 2π x . We evolve the test until the sine wave has completed 10 crossings of the domain. Figure 2 plots the absolute error for the number density versus the number of elements N . We see the errors obtained with SSPRK3 and PD-ARS3 are smallest and decrease as N −3 , as expected for a scheme combining third-order accurate time stepping with third-order accurate spatial discretization. For all the other schemes, using second-order accurate explicit time stepping, the error decreases as N −2 . Among the secondorder schemes, SSP2332 has the smallest error. In the streaming limit, the PD-ARS schemes reduce to SSPRK schemes -PD-ARS2 to SSPRK2 and PD-ARS3 to SSPRK3, respectively. Therefore, the absolute errors of PD-ARS schemes and SSPRK schemes are indistinguishable.

Sine Wave Damping
The next test, adapted from [54], is designed for absorption- For σ A = 0.1, 1 and 10 we evolve the test until the initial condition has been damped by a factor e −10 . Figure 3 shows convergence results of the test in the relative error. Results for σ A = 0.1, 1, and 10 are plotted with red, green, and blue lines, respectively. SSP2332 is the most accurate among these schemes for σ A = 1 and 10. For σ A = 0.1, PD-ARS2 is the most accurate scheme for N = 8 and N = 16. We have seen the same behavior for the scheme proposed by McClarren et al. [55] (PC2 in [29]). Since these are special cases for N = 8 and N = 16, we do not recommend PD-ARS2 instead of SSP2332 for the damping problem. Only SSP2332, a second-order accurate scheme, displays a second-order convergence rate. The PD-ARS schemes are first-order accurate.

Sine Wave Diffusion
The final test with known smooth solutions, adopted from [56], is the sine wave diffusion test; i.e.
We evolve with σ S = 10 2 , 10 3 , and 10 4 , and adjust the end time so that t end /σ S = 1, at which time the amplitude of the sine wave has been reduced by a factor e −π 2 /27 ≈ 0.694 for all values of σ S . Figure 4 shows the absolute error, obtained using different values of σ S , for various IMEX schemes at t = t end , versus N . Results for σ S = 10 2 , 10 3 , and 10 4 are plotted with red, green, and blue lines, respectively. SSP2332 and PD-ARS schemes display third-order accuracy for the number density, J , and second-order accuracy for H x , and their errors are difficult to distinguish. For σ = 10 2 , the errors in the number density J do not drop below 10 −6 because of differences between the two-moment model and the diffusion equation used to obtain the analytic solution.

Neutrino Stationary State Test
In this section we consider a more "realistic" test: two-dimensional multigroup neutrino transport with emission, absorption, and isoenergetic scattering through a stationary background. This test is designed to test the realizability-preserving properties of the PD-ARS schemes. In the left panel in Figure 5 we plot the thermal state of the background, which mimics a collapsed stellar core: Mass Density: ρ = 4 × 10 14 × 7.5 (7.5 + (r/5 km) 4 Temperature: T = 1.5 × 10 11 × 1 1 + (r/50 km) 2 Electron Fraction: Y e = 0.25 × 1 + where the radius r = (x 1 ) 2 + (x 2 ) 2 is in kilometers. In the right panel in Figure 5 we plot the neutrino opacities, computed by interpolation in a pre-calculated opacity table based on [57]. This test is computed using Cartesian coordinates on a two-dimensional domain D = {x ∈ R 2 : x 1 ∈ [0, 200] km, x 2 ∈ [0, 200] km}, using a grid of 128 elements in each direction, 10 energy groups covering 0-300 MeV, reflecting inner boundaries, and outflow outer boundaries. Because we use Cartesian coordinates in two spatial dimensions, this problem has cylindrical geometry, and results in an artificial stratification of the radiation quantities. We initialize the neutrino number density to J = 10 −99 , the number flux density to H = 0, and evolve until an approximate steady state is reached (t = 5 ms). The background is kept fixed during the entire run. For this test, we employed both CB and Minerbo closures. We attempted to run this test with our PD-ARS schemes, SSP2332 from [47], IMEXRKCB2 proposed by Cavaglieri & Bewley [58], and the IMEX PC2 scheme proposed by McClarren et al. [55]. Only the PD-ARS schemes produce realizable moments and are able to evolve to a steady state with either CB or Minerbo closure. SSP2332 and IMEX PC2 failed after a few time steps with either CB or Minerbo closure because of the development of unrealizable moments. Even though IMEXRKCB2 with Minerbo closure can run and reach a steady state in this test, its results are different from that of PD-ARS schemes with CB closure, and there is no guarantee of stability.
Results obtained with the IMEX PD-ARS schemes are plotted in Figure 6 for various times: t = 0.01 ms (top panels), 0.35 ms (middle panels), and 5.0 ms (bottom panels). In the left column we plot the solution in the |H|J -plane. In the middle column we show scatter plots of the number density J versus radius for select neutrino energies: 5 MeV (red lines), 16 MeV (magenta lines), and 93 MeV (blue lines). In the right column we plot the flux factor |H|/J versus radius for the same neutrino energies as in the middle column. In the left panels, each solution point in the domain is marked as a red dot in the |H|J -plane, and the realizable domain is shown as the light blue region. The figures show that all the states in the simulation with the PD-ARS schemes are realizable. In the middle and right column, we can see how the neutrinos are generated near the core, stream out, and eventually reach an equilibrium distribution over the phase space. The oscillations in the flux factor seen in the right columns are associated with steep gradients in the radiation field as the initial transient propagates through the computational domain. Note that we do not apply any limiters to prevent oscillations in the numerical solution, and these will likely go away when we implement slope limiters; e.g., as described in [59]. The fact that we can still evolve the solution to a steady state speaks to the robustness of the scheme.

Summary and Outlook
We have developed IMEX schemes suitable for a two-moment model of neutrino transport that obey Fermi-Dirac statistics. The scheme employs algebraic closure based on Fermi-Dirac statistics, high-order discontinuous Galerkin methods for spatial discretization, and convexinvariant time integration to maintain realizability of the moments. Since the realizable domain is convex and its convexity can be inherited by a convex combination, a scheme having convex combinations as its stages can preserve the realizable domain. This encouraged us to construct realizability-preserving time integrators, realizability-preserving IMEX schemes, and a method with a realizability-preserving IMEX time integrator and high-order DG method.
In the applications that motivate this work, the neutrino distribution function can vary from 0 to 1. Hence, we have considered algebraic closures based on Fermi-Dirac statistics for both low and high occupancy. Among the seven algebraic closures we considered -Kershaw [35], Wilson [40,41], Levermore [27], Minerbo [23], Janka 1 [43], Janka 2 [42], and Cernohorsky & Bludman [44] -only the Cernohorsky & Bludman closure obeys Fermi-Dirac statistics for all occupancies. As a result, we employed the Cernohorsky & Bludman closure for the neutrino stationary state test in Section 5.2. We also ran our code with Minerbo closure, IMEX PC2, IMEX SSP2332 and IMEXRKCB2 schemes, and the results show that only PD-ARS schemes have stability. In addition, closures have impact on the simulation result. As we observed, both using PD-ARS scheme, there are ∼ 30% difference in the neutrino number densities (and relaxation time) between the results obtained with Minerbo closure and that with CB closure. Even though RKCB2 with Minerbo closure luckily survived our test, the results it gave were compromised: they were closer to the results of the PD-ARS scheme with Minerbo closure than to the results of the PD-ARS scheme with CB closure. In what way and to what degree the results are in fact compromised either by the closure or by a particular correction step for unrealizable moments are difficult to determine fully and is left for further study.
Two PD-ARS schemes are proposed. The one with SSPRK2 has second-order accuracy while the other with SSPRK3 has third-order accuracy, and both have the strong-stability preserving property in the streaming limit. Their accuracy was demonstrated on problems with known smooth solutions in streaming, absorption, and scattering-dominated regimes. The neutrino transport test with emission, absorption, and isoenergetic scattering through a stationary background, was designed to test the convex-invariance of our PD-ARS schemes. The neutrino stationary state test shows that a method combining an algebraic closure based on Fermi-Dirac statistics and convex-invariant time integration is promising for robust CCSN simulation.
In this work, we adopted Cartesian coordinates, a linear collision term, and a fixed material