A DG-IMEX Method for Two-moment Neutrino Transport: Nonlinear Solvers for Neutrino–Matter Coupling^*

In nuclear astrophysics applications, neutrino transport can be modeled by the Boltzmann equation. In this paper, we consider a nonrelativistic Boltzmann equation,

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mathsf{c}}}{\partial }_{t}f+{ \mathcal T }(f)={ \mathcal C }(f),\end{eqnarray} \tag{ 1 }$$

which governs the particle distribution function f( x , p , t) that describes the density of particles at position ${\boldsymbol{x}}\,=({x}^{1},{x}^{2},{x}^{3})\in {{\mathbb{R}}}^{3}$ with momentum ${\boldsymbol{p}}=({p}^{1},{p}^{2},{p}^{3})\in {{\mathbb{R}}}^{3}$ at time $t\in {{\mathbb{R}}}^{+}$, where c is the speed of light. Adopting curvilinear phase-space coordinates, the advection operator ${ \mathcal T }$ takes the form (see, e.g., Endeve et al. 2015)

$$\begin{eqnarray}&&{ \mathcal T }(f):= \displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \sum _{i=1}^{3}{\partial }_{{x}^{i}}\left(\sqrt{\gamma }\,{H}_{{\boldsymbol{x}}}^{i}f\right)+\displaystyle \frac{1}{\sqrt{\lambda }}\displaystyle \sum _{i=1}^{3}{\partial }_{{p}^{i}}\left(\sqrt{\lambda }\,{H}_{{\boldsymbol{p}}}^{i}f\right),\end{eqnarray} \tag{ 2 }$$

where ${H}_{{\boldsymbol{x}}}^{i}f$ and ${H}_{{\boldsymbol{p}}}^{i}f$ denote the position space flux and momentum space flux in the corresponding ith direction, respectively, and γ, λ are the determinants of the position space and momentum space metric tensors, respectively. While the form of the advection operator in Equation (2) holds for more general (nonorthogonal) spacetime and momentum space bases (see, e.g., Cardall et al. 2013a), we restrict ourselves to orthogonal bases in this paper. The position space flux considered in this paper takes the form ${H}_{{\boldsymbol{x}}}^{i}f=({p}^{i}/| {\boldsymbol{p}}| )f$, which is proportional to the particle propagation direction. We will restrict ourselves to spherical polar momentum coordinates, and in this case, ${H}_{{\boldsymbol{p}}}^{i}f$ can be obtained from Equations (A15) and (A16) in Endeve et al. (2015). The collision operator ${ \mathcal C }$ models the interactions between particles and a material background and includes emission, absorption, elastic scattering on nucleons and nuclei, NES, and thermal pair processes from electron–positron creation and annihilation. This is the neutrino opacity set described in Bruenn (1985; see Table 1 therein).

In this work, we consider the transport of electron neutrinos (ν_e ) and antineutrinos (${\bar{\nu }}_{e}$), which results in the coupled equations

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mathsf{c}}}{\partial }_{t}f+{ \mathcal T }(f)={ \mathcal C }(f,\bar{f}),\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mathsf{c}}}{\partial }_{t}\bar{f}+{ \mathcal T }(\bar{f})=\bar{{ \mathcal C }}(\bar{f},f),\end{eqnarray} \tag{ 3b }$$

where the neutrino and antineutrino distribution functions are denoted with f and $\bar{f}$, respectively, and the collision terms ${ \mathcal C }$ and $\bar{{ \mathcal C }}$ both depend on f and $\bar{f}$, since they include the thermal pair production and annihilation processes of neutrino–antineutrino pairs. Before giving a detail formulation of the collision terms, we first change to spherical polar momentum coordinates (ε, ω ) and decompose the neutrino three-momentum as p = ε ℓ ( ω ), where $\varepsilon \in {{\mathbb{R}}}^{+}$ denotes neutrino energy and the unit vector ${\boldsymbol{\ell }}({\boldsymbol{\omega }})\in {{\mathbb{R}}}^{3}$ only depends on the angular direction ${\boldsymbol{\omega }}\in {{\mathbb{S}}}^{2}$ relative to a local orthonormal basis. The momentum space volume element is then decomposed into d p = dV_ε d ω , where dV_ε $:=$ ε² d ε is the spherical shell energy volume element and d ω is the momentum space angular element. With this notation, the particle distribution f( x , p , t) can be written as f( x , ω , ε, t). At each x and t, the neutrino collision operator then takes the form (Bruenn 1985)

$$\begin{eqnarray}\,\begin{array}{rcl}{ \mathcal C }(f,\bar{f}) & = & \left(1-f({\boldsymbol{\omega }},\varepsilon )\right)\hat{\eta }(\varepsilon )-\hat{\chi }(\varepsilon )f({\boldsymbol{\omega }},\varepsilon )\\ & & +\ \displaystyle \frac{{\varepsilon }^{2}}{{\mathsf{c}}{\left({\mathsf{hc}}\right)}^{3}}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{R}^{{\rm\small{IS}}}({\boldsymbol{\ell }}\cdot {\boldsymbol{\ell }}^{\prime} ,\varepsilon )f({\boldsymbol{\omega }}^{\prime} ,\varepsilon )d{\boldsymbol{\omega }}^{\prime} \\ & & -\ \displaystyle \frac{{\varepsilon }^{2}}{{\mathsf{c}}{\left({\mathsf{hc}}\right)}^{3}}f({\boldsymbol{\omega }},\varepsilon ){\displaystyle \int }_{{{\mathbb{S}}}^{2}}{R}^{{\rm\small{IS}}}({\boldsymbol{\ell }}\cdot {\boldsymbol{\ell }}^{\prime} ,\varepsilon )\,d{\boldsymbol{\omega }}^{\prime} \\ & & +\ \displaystyle \frac{1}{{\mathsf{c}}{\left({\mathsf{hc}}\right)}^{3}}\left(1-f({\boldsymbol{\omega }},\varepsilon )\right){\displaystyle \int }_{{{\mathbb{R}}}^{+}}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}\\ & & \times \ {R}^{{\rm{I}}{\rm\small{N}}}({\boldsymbol{\ell }}\cdot {\boldsymbol{\ell }}^{\prime} ,\varepsilon ,\varepsilon ^{\prime} )f({\boldsymbol{\omega }}^{\prime} ,\varepsilon ^{\prime} )d{\boldsymbol{\omega }}^{\prime} \,{{dV}}_{\varepsilon ^{\prime} }\\ & & -\ \displaystyle \frac{1}{{\mathsf{c}}{\left({\mathsf{hc}}\right)}^{3}}f({\boldsymbol{\omega }},\varepsilon ){\displaystyle \int }_{{{\mathbb{R}}}^{+}}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}{R}^{{\rm{O}}{\rm\small{UT}}}\\ & & \times \ ({\boldsymbol{\ell }}\cdot {\boldsymbol{\ell }}^{\prime} ,\varepsilon ,\varepsilon ^{\prime} )\left(1-f({\boldsymbol{\omega }}^{\prime} ,\varepsilon ^{\prime} )\right)d{\boldsymbol{\omega }}^{\prime} \,{{dV}}_{\varepsilon ^{\prime} }\\ & & +\ \displaystyle \frac{1}{{\mathsf{c}}{\left({\mathsf{hc}}\right)}^{3}}\left(1-f({\boldsymbol{\omega }},\varepsilon )\right){\displaystyle \int }_{{{\mathbb{R}}}^{+}}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}\\ & & \times \ {R}^{{\rm{P}}{\rm\small{R}}}({\boldsymbol{\ell }}\cdot {\boldsymbol{\ell }}^{\prime} ,\varepsilon ,\varepsilon ^{\prime} )\left(1-\bar{f}({\boldsymbol{\omega }}^{\prime} ,\varepsilon ^{\prime} )\right)d{\boldsymbol{\omega }}^{\prime} \,{{dV}}_{\varepsilon ^{\prime} }\\ & & -\ \displaystyle \frac{1}{{\mathsf{c}}{\left({\mathsf{hc}}\right)}^{3}}f({\boldsymbol{\omega }},\varepsilon ){\displaystyle \int }_{{{\mathbb{R}}}^{+}}{\displaystyle \int }_{{{\mathbb{S}}}^{2}}\\ & & \times \ {R}^{{\rm{A}}{\rm\small{N}}}({\boldsymbol{\ell }}\cdot {\boldsymbol{\ell }}^{\prime} ,\varepsilon ,\varepsilon ^{\prime} )\bar{f}({\boldsymbol{\omega }}^{\prime} ,\varepsilon ^{\prime} )d{\boldsymbol{\omega }}^{\prime} \,{{dV}}_{\varepsilon ^{\prime} },\end{array}\,\end{eqnarray} \tag{ 4 }$$

with h the Planck constant, $\hat{\eta }$ the emissivity, $\hat{\chi }$ the absorption opacity, R ^IS the elastic (isoenergetic) scattering kernel (due to scattering with nucleons and nuclei), R^IN and R^OUT the NES kernels, and R^PR and R^AN the thermal production and annihilation kernels due to pair processes. The antineutrino collision operator $\bar{{ \mathcal C }}$ is defined analogously with $\hat{\eta }$, $\hat{\chi }$, and the opacity kernels R^OP, OP = IS, IN, OUT, PR, AN, replaced by their antineutrino counterparts.

For thermal emission and absorption, we follow the approach in, e.g., Burrows et al. (2006) and rewrite

$$\begin{eqnarray}&&\left(1-f\right)\hat{\eta }-\hat{\chi }\,f=\chi \left({f}_{0}-f\right),\end{eqnarray} \tag{ 5 }$$

where the effective opacity and equilibrium distribution are defined as

$$\begin{eqnarray}&&\chi :=\hat{\eta }+\hat{\chi }\,{\rm{a}}{\rm{n}}{\rm{d}}\,{f}_{0}:=\displaystyle \frac{\hat{\eta }}{\hat{\eta }+\hat{\chi }},\end{eqnarray} \tag{ 6 }$$

respectively. Similarly, the antineutrino emission and absorption term is written as $\bar{\chi }\left({\bar{f}}_{0}-\bar{f}\right)$, with $\bar{\chi }$ and ${\bar{f}}_{0}$ defined analogously. The equilibrium distributions f₀ and ${\bar{f}}_{0}$ are given by the Fermi–Dirac distribution, which is isotropic in angle ω . Specifically,

$$\begin{eqnarray}&&{f}_{0}(\varepsilon )=\displaystyle \frac{1}{{e}^{(\varepsilon -\mu )/{{\mathsf{k}}}_{{\mathsf{B}}}T}+1}\,{\rm{a}}{\rm{n}}{\rm{d}}\,{\bar{f}}_{0}(\varepsilon )=\displaystyle \frac{1}{{e}^{(\varepsilon -\bar{\mu })/{{\mathsf{k}}}_{{\mathsf{B}}}T}+1},\end{eqnarray} \tag{ 7 }$$

where k _B is the Boltzmann constant, T is the matter temperature, $\mu := {\mu }_{{e}^{-}}+({\mu }_{p}-{\mu }_{n})$ is the neutrino chemical potential, and $\bar{\mu }:= {\mu }_{{e}^{+}}-({\mu }_{p}-{\mu }_{n})$ is the antineutrino chemical potential. Here ${\mu }_{{e}^{-(+)}}$ is the electron (positron) chemical potential, μ_n is the neutron chemical potential, and μ_p is the proton chemical potential, which are evaluated from an appropriate EoS. Since ${\mu }_{{e}^{+}}=-{\mu }_{{e}^{-}}$, we have $\bar{\mu }=-\mu$.

Following Bruenn (1985), the neutrino scattering and pair process kernels are approximated with L-term Legendre expansions in the cosine of the scattering angle $\alpha := {\boldsymbol{\ell }}\cdot {\boldsymbol{\ell }}^{\prime}$, with expansion coefficients Φ depending on the neutrino energy. Specifically,

$$\begin{eqnarray}\,\begin{array}{rcl}{R}^{{\rm{I}}{\rm{S}}}(\alpha ,\varepsilon ) & \approx & \displaystyle \sum _{l=0}^{L}{{\rm{\Phi }}}_{l}^{{\rm{I}}{\rm{S}}}(\varepsilon ){P}_{l}(\alpha )\,{\rm{a}}{\rm{n}}{\rm{d}}\,\\ {R}^{{\rm{O}}{\rm{P}}}(\alpha ,\varepsilon ,{\varepsilon }^{{\rm{{\prime} }}}) & \approx & \displaystyle \sum _{l=0}^{L}{{\rm{\Phi }}}_{l}^{{\rm{O}}{\rm{P}}}(\varepsilon ,{\varepsilon }^{{\rm{{\prime} }}}){P}_{l}(\alpha ),\end{array}\,\end{eqnarray} \tag{ 8 }$$

where OP = IN, OUT, PR, and AN; P_l denotes the Legendre polynomial of degree l; and the expansion coefficients of degree l are given by

$$\begin{eqnarray}\,\begin{array}{rcl}{{\rm{\Phi }}}_{l}^{{\rm{I}}{\rm{S}}}(\varepsilon ) & = & \displaystyle \frac{1}{{C}_{l}}{\displaystyle \int }_{-1}^{1}{R}^{{\rm{I}}{\rm{S}}}(\alpha ,\varepsilon ){P}_{l}(\alpha )d\alpha \,{\rm{a}}{\rm{n}}{\rm{d}}\,\\ {{\rm{\Phi }}}_{l}^{{\rm{O}}{\rm{P}}}(\varepsilon ,{\varepsilon }^{{\rm{{\prime} }}}) & = & \displaystyle \frac{1}{{C}_{l}}{\displaystyle \int }_{-1}^{1}{R}^{{\rm{O}}{\rm{P}}}(\alpha ,\varepsilon ,{\varepsilon }^{{\rm{{\prime} }}}){P}_{l}(\alpha )d\alpha .\end{array}\,\end{eqnarray} \tag{ 9 }$$

Here ${C}_{l}:= {\int }_{-1}^{1}{P}_{l}^{2}(\alpha )d\alpha$, l = 0, ..., L, are normalization constants. In this work, we use P₀(α) = 1 and P₁(α) = α; thus, the normalization constants are C₀ = 2 and ${C}_{1}=\tfrac{2}{3}$. The antineutrino opacity kernels are approximated analogously. Due to particle conservation and detailed balance (e.g., Cernohorsky 1994), the NES and pair process kernels satisfy

$$\begin{eqnarray}\,\begin{array}{rcl}{R}^{{\rm{I}}{\rm\small{N}}}(\alpha ,\varepsilon ,\varepsilon ^{\prime} ) & = & {R}^{{\rm{O}}{\rm\small{UT}}}(\alpha ,\varepsilon ^{\prime} ,\varepsilon ),\quad \\ {R}^{{\rm{O}}{\rm\small{UT}}}(\alpha ,\varepsilon ,\varepsilon ^{\prime} ) & = & {R}^{{\rm{O}}{\rm\small{UT}}}(\alpha ,\varepsilon ^{\prime} ,\varepsilon ){e}^{(\varepsilon -\varepsilon ^{\prime} )/{{\mathsf{k}}}_{{\mathsf{B}}}T},\\ {R}^{{\rm{A}}{\rm\small{N}}}(\alpha ,\varepsilon ,\varepsilon ^{\prime} ) & = & {\bar{R}}^{{\rm{A}}{\rm\small{N}}}(\alpha ,\varepsilon ^{\prime} ,\varepsilon ),\qquad \\ {R}^{{\rm{P}}{\rm\small{R}}}(\alpha ,\varepsilon ,\varepsilon ^{\prime} ) & = & {R}^{{\rm{A}}{\rm\small{N}}}(\alpha ,\varepsilon ,\varepsilon ^{\prime} ){e}^{-(\varepsilon +\varepsilon ^{\prime} )/{{\mathsf{k}}}_{{\mathsf{B}}}T}.\end{array}\,\end{eqnarray} \tag{ 10 }$$

Note that these symmetry properties are enforced in the energy space; thus, they are preserved in the angular approximation for the kernels in Equation (8).

In the following, as a first approximation, we will only include the isotropic part of the kernels (i.e., L = 0). More realistic treatments would include linear corrections (L = 1), while further corrections (L > 1) have been shown to result in only minor differences for NES (e.g., Smit & Cernohorsky 1996).

The need for high spatial resolution and unconstrained spatial dimensionality, e.g., to capture fluid dynamics in our target applications, renders direct solutions of the Boltzmann equation too expensive. However, neutrino heating rates are sensitive to the neutrino energy distribution, which demands retention of the energy dimension of momentum space. Therefore, to balance computational cost with physical fidelity, we settle for solving for a finite number of angular moments of the distribution function by adopting a two-moment model. Two-moment models are widely used to model neutrino transport in CCSNe (e.g., Just et al. 2015; O'Connor 2015; Kuroda et al. 2016; Roberts et al. 2016; Skinner et al. 2019). In the spectral two-moment model, we solve for the zeroth and first moments of the neutrino distribution function, while the second moments are obtained from a closure procedure. These moments are defined respectively as

$$\begin{eqnarray}&&\,\begin{array}{l}\left\{{ \mathcal J },{{ \mathcal H }}^{i},{{ \mathcal K }}^{{ij}}\right\}({\boldsymbol{x}},\varepsilon ,t)=\displaystyle \frac{1}{4\pi }{\displaystyle \int }_{{{\mathbb{S}}}^{2}}f({\boldsymbol{x}},{\boldsymbol{\omega }},\varepsilon ,t)\\ \quad \ \times \ \left\{1,{{\ell }}^{i}({\boldsymbol{\omega }}),{{\ell }}^{i}({\boldsymbol{\omega }}){{\ell }}^{j}({\boldsymbol{\omega }})\right\}\,d{\boldsymbol{\omega }},\end{array}\,\end{eqnarray} \tag{ 11 }$$

with moments for antineutrinos ($\bar{{ \mathcal J }}$, ${\bar{{ \mathcal H }}}^{i}$, and ${\bar{{ \mathcal K }}}^{{ij}}$) defined analogously. Then, the zeroth moment ${ \mathcal J }$ ($\bar{{ \mathcal J }}$) is the spectral number density of neutrinos (antineutrinos), the first moment ${{ \mathcal H }}^{i}$ (${\bar{{ \mathcal H }}}^{i}$) is the spectral number flux density of neutrinos (antineutrinos), and the second moment ${{ \mathcal K }}^{{ij}}$ (${\bar{{ \mathcal K }}}^{{ij}}$) is proportional to the spectral pressure tensor of neutrinos (antineutrinos). Since the neutrino distribution function is bounded between 0 and 1, it can be shown that the moments must satisfy the bounds (Larecki & Banach 2011)

$$\begin{eqnarray}&&0\leqslant { \mathcal J }\leqslant 1\quad {\rm{and}}\quad | {\boldsymbol{ \mathcal H }}| \leqslant \left(1-{ \mathcal J }\right){ \mathcal J },\end{eqnarray} \tag{ 12 }$$

where ${\boldsymbol{ \mathcal H }}:= ({{ \mathcal H }}_{1},{{ \mathcal H }}_{2},{{ \mathcal H }}_{3})$ and $| {\boldsymbol{ \mathcal H }}| =\sqrt{{{ \mathcal H }}_{i}{{ \mathcal H }}^{i}}$, with ${{ \mathcal H }}_{i}$ associated to ${{ \mathcal H }}^{k}$ via the spatial metric γ_ik , i.e., ${{ \mathcal H }}_{i}={\gamma }_{{ik}}{{ \mathcal H }}^{k}$. (For notational convenience in this section, we sometimes use Einstein's summation convention where repeated Latin indices imply summation from 1 to 3.) Bounds equivalent to those in Equation (12) also hold for the antineutrino moments.

Taking the zeroth and first moments of Equation 3(a) leads to

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mathsf{c}}}{\partial }_{t}{ \mathcal J }+\displaystyle \frac{1}{4\pi }{\int }_{{{\mathbb{S}}}^{2}}{ \mathcal T }(f)\,d{\boldsymbol{\omega }}=\displaystyle \frac{1}{4\pi }{\int }_{{{\mathbb{S}}}^{2}}{ \mathcal C }(f,\bar{f})\,d{\boldsymbol{\omega }},\end{eqnarray} \tag{ 13a }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mathsf{c}}}{\partial }_{t}{{ \mathcal H }}_{j}+\displaystyle \frac{1}{4\pi }{\int }_{{{\mathbb{S}}}^{2}}{ \mathcal T }(f){{\ell }}_{j}\,d{\boldsymbol{\omega }}=\displaystyle \frac{1}{4\pi }{\int }_{{{\mathbb{S}}}^{2}}{ \mathcal C }(f,\bar{f}){{\ell }}_{j}\,d{\boldsymbol{\omega }},\end{eqnarray} \tag{ 13b }$$

which, after plugging in the definitions of ${ \mathcal T }$ and ${ \mathcal C }$ in Equations (2) and (4), results in the moment equations for the number density and number flux,

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mathsf{c}}}{\partial }_{t}{ \mathcal J }+\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \sum _{i=1}^{3}{\partial }_{{x}^{i}}\left(\sqrt{\gamma }\,{{ \mathcal H }}^{i}\right)={\eta }_{{\rm{T}}{\rm\small{OT}}}-{\chi }_{{\rm{T}}{\rm\small{OT}}}\,{ \mathcal J },\end{eqnarray} \tag{ 14a }$$

$$\begin{eqnarray}&&\,\begin{array}{l}\displaystyle \frac{1}{{\mathsf{c}}}{\partial }_{t}{{ \mathcal H }}_{j}+\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \sum _{i=1}^{3}{\partial }_{{x}^{i}}\left(\sqrt{\gamma }\,{{ \mathcal K }}_{\,j}^{i}\right)=\displaystyle \frac{1}{2}\,{{ \mathcal K }}^{{ik}}\,{\partial }_{{x}^{j}}{\gamma }_{{ik}}\\ \qquad \ -\ ({\chi }_{{\rm{T}}{\rm\small{OT}}}+{\sigma }_{{\rm\small{IS}}}){{ \mathcal H }}_{j},\end{array}\,\end{eqnarray} \tag{ 14b }$$

where we have used the facts that the position space flux p /∣ p ∣f = ℓ ( ω )f and that contributions from momentum space fluxes vanish due to boundary conditions. (An analogous set of moment equations is derived for antineutrinos.) Equation (14) holds for general, time-independent curvilinear spatial coordinates encoded in the spatial metric γ_ij , including the commonly used Cartesian, spherical polar, and cylindrical coordinates. The metric tensor is used to raise and lower indices on vectors and tensors; e.g., ${{ \mathcal H }}_{j}={\gamma }_{{jk}}{{ \mathcal H }}^{k}$, ${{ \mathcal K }}_{j}^{i}={\gamma }_{{jk}}{{ \mathcal K }}^{{ik}}$.

Remark 1. We note that Equation (14) can also be obtained from the nonrelativistic (i.e., zero fluid velocity and no gravitational fields) limit of the moment equations in Shibata et al. (2011) and Cardall et al. (2013b; see also the number conservative two-moment model discussed by Mezzacappa et al. 2020; their Equations (123) and (125)). Moreover, when including relativistic effects, one is, among other things, confronted with choosing appropriate momentum space coordinates. The most common (and perhaps most natural) choice is to use momentum space coordinates in the frame of reference of the inertial observer instantaneously comoving with the fluid (i.e., the comoving frame), as opposed to the so-called laboratory frame (see, e.g., Mihalas & Mihalas 1999). (In the absence of fluid motion and gravitational fields, which is assumed here, there is no distinction between the comoving and laboratory frames.) Specifically, the choice of comoving frame momentum coordinates provides the most straightforward framework for describing neutrino–matter interactions, and this is the choice we intend to make when including relativistic effects in the future. However, this choice complicates the advection operator associated with the moment equations, which then includes Doppler and/or gravitational frequency shift terms. The inclusion of these terms is beyond the scope of the present paper but will be considered in a future study.

On the right-hand side of Equation (14), the elastic scattering opacity is given by ${\sigma }_{{\rm\small{IS}}}=\tfrac{4\pi {\varepsilon }^{2}}{{\mathsf{c}}{({\mathsf{hc}})}^{3}}{{\rm{\Phi }}}_{0}^{{\rm\small{IS}}}(\varepsilon )$, and we define the total emissivity and total opacity as

$$\begin{eqnarray}\,\begin{array}{rcl}{\eta }_{{\rm{T}}{\rm{O}}{\rm{T}}}({ \mathcal J },\bar{{ \mathcal J }}) & = & \chi \,{{ \mathcal J }}_{0}+{\eta }_{{\rm{S}}{\rm{C}}}({ \mathcal J })+{\eta }_{{\rm{T}}{\rm{P}}}(\bar{{ \mathcal J }})\,{\rm{a}}{\rm{n}}{\rm{d}}\,\\ {\chi }_{{\rm{T}}{\rm{O}}{\rm{T}}}({ \mathcal J },\bar{{ \mathcal J }}) & = & \chi +{\chi }_{{\rm{S}}{\rm{C}}}({ \mathcal J })+{\chi }_{{\rm{T}}{\rm{P}}}(\bar{{ \mathcal J }}),\end{array}\,\end{eqnarray} \tag{ 15 }$$

where ${{ \mathcal J }}_{0}=\tfrac{1}{4\pi }{\int }_{{{\mathbb{S}}}^{2}}{f}_{0}\,d{\boldsymbol{\omega }}={f}_{0}$. Let $d{\tilde{V}}_{\varepsilon }:= 4\pi {\varepsilon }^{2}d\varepsilon$ denote the scaled spherical shell energy volume element; then the opacity terms in Equation (15) are defined, respectively, as the scattering emissivity

$$\begin{eqnarray}&&{\eta }_{{\rm{S}}{\rm\small{C}}}({ \mathcal J })=\displaystyle \frac{1}{{\mathsf{c}}{\left({\mathsf{hc}}\right)}^{3}}{\int }_{{{\mathbb{R}}}^{+}}{{\rm{\Phi }}}_{0}^{{\rm{I}}{\rm\small{N}}}(\varepsilon ,\varepsilon ^{\prime} ){ \mathcal J }(\varepsilon ^{\prime} )d{\tilde{V}}_{\varepsilon ^{\prime} },\end{eqnarray} \tag{ 16 }$$

the scattering opacity

$$\begin{eqnarray}\,\begin{array}{rcl}{\chi }_{{\rm{S}}{\rm\small{C}}}({ \mathcal J }) & = & \displaystyle \frac{1}{{\mathsf{c}}{\left({\mathsf{hc}}\right)}^{3}}{\displaystyle \int }_{{{\mathbb{R}}}^{+}}\left[\,{{\rm{\Phi }}}_{0}^{{\rm{I}}{\rm\small{N}}}(\varepsilon ,\varepsilon ^{\prime} ){ \mathcal J }(\varepsilon ^{\prime} )\right.\\ & & +\ \left.{{\rm{\Phi }}}_{0}^{{\rm\small{OUT}}}(\varepsilon ,\varepsilon ^{\prime} )\left(1-{ \mathcal J }(\varepsilon ^{\prime} )\right)\right]\,d{\tilde{V}}_{\varepsilon ^{\prime} },\end{array}\,\end{eqnarray} \tag{ 17 }$$

the emissivity due to thermal pair processes

$$\begin{eqnarray}&&{\eta }_{{\rm{T}}{\rm\small{P}}}(\bar{{ \mathcal J }})=\displaystyle \frac{1}{{\mathsf{c}}{\left({\mathsf{hc}}\right)}^{3}}{\int }_{{{\mathbb{R}}}^{+}}{{\rm{\Phi }}}_{0}^{{\rm{P}}{\rm\small{R}}}(\varepsilon ,\varepsilon ^{\prime} )\left(1-\bar{{ \mathcal J }}(\varepsilon ^{\prime} )\right)d{\tilde{V}}_{\varepsilon ^{\prime} },\end{eqnarray} \tag{ 18 }$$

and the opacity due to thermal pair processes

$$\begin{eqnarray}\,\begin{array}{rcl}{\chi }_{{\rm{T}}{\rm\small{P}}}(\bar{{ \mathcal J }}) & = & \displaystyle \frac{1}{{\mathsf{c}}{\left({\mathsf{hc}}\right)}^{3}}{\displaystyle \int }_{{{\mathbb{R}}}^{+}}\left[\,{{\rm{\Phi }}}_{0}^{{\rm{P}}{\rm\small{R}}}(\varepsilon ,\varepsilon ^{\prime} )\left(1-\bar{{ \mathcal J }}(\varepsilon ^{\prime} )\right)\right.\\ & & +\ \left.{{\rm{\Phi }}}_{0}^{{\rm{A}}{\rm\small{N}}}(\varepsilon ,\varepsilon ^{\prime} )\bar{{ \mathcal J }}(\varepsilon ^{\prime} )\right]\,d{\tilde{V}}_{\varepsilon ^{\prime} }.\end{array}\,\end{eqnarray} \tag{ 19 }$$

We make the following remarks on Equation (14): (i) the scattering and pair process opacities depend on ${ \mathcal J }$ and $\bar{{ \mathcal J }}$, respectively, due to the Fermi blocking factors; (ii) since ${{\rm{\Phi }}}_{0}^{{\rm\small{OP}}}\geqslant 0$ and ${ \mathcal J }$, $\bar{{ \mathcal J }}$ are between 0 and 1, we have σ _IS , η_SC, χ_SC, η_TP, and χ_TP ≥ 0; and (iii) the emissivities and opacities depend on the neutrino energy ε and local matter states (e.g., density ρ, temperature T, and electron fraction Y_e ).

To close the two-moment model in Equation (14), we adopt an algebraic closure of the form (Levermore 1984)

$$\begin{eqnarray}&&{{ \mathcal K }}_{\,j}^{i}=\displaystyle \frac{1}{2}\left[\left(1-\psi \right){\delta }_{\,j}^{i}+\left(3\,\psi -1\right){\widehat{h}}^{i}\,{\widehat{h}}_{j}\,\right]\,{ \mathcal J },\end{eqnarray} \tag{ 20 }$$

where ${\widehat{h}}^{i}={{ \mathcal H }}^{i}/| {\boldsymbol{ \mathcal H }}|$ are components of a unit vector parallel to ${\boldsymbol{ \mathcal H }}$, and the Eddington factor $\psi =\psi ({ \mathcal J },h)$ with $h=| {\boldsymbol{ \mathcal H }}| /{ \mathcal J }$. We use the maximum entropy Eddington factor of Cernohorsky & Bludman (1994),

$$\begin{eqnarray}&&\psi ({ \mathcal J },h)=\displaystyle \frac{1}{3}+\displaystyle \frac{2(1-{ \mathcal J })(1-2\,{ \mathcal J })}{3}\,{\rm{\Theta }}\left(\displaystyle \frac{h}{1-{ \mathcal J }}\right),\end{eqnarray} \tag{ 21 }$$

where the closure polynomial is given by

$$\begin{eqnarray}&&{\rm{\Theta }}(x)=\displaystyle \frac{1}{5}(3-x+3\,{x}^{2}){x}^{2}.\end{eqnarray} \tag{ 22 }$$

We point out that this closure is based on Fermi–Dirac statistics and suitable for designing numerical methods for the two-moment model satisfying the bounds in Equation (12) (e.g., Chu et al. 2019).

For the antineutrino transport (Equation 3(b)), a two-moment model on the antineutrino number density $\bar{{ \mathcal J }}$ and number flux $\bar{{\boldsymbol{ \mathcal H }}}$ can be derived following similar procedure as in Equations (13)–(19). We then close the resulting two-moment model using a closure analogous to Equation (20). To simplify the notation, we denote the neutrino and antineutrino moments as ${\boldsymbol{ \mathcal M }}:= ({ \mathcal J },{\boldsymbol{ \mathcal H }},\bar{{ \mathcal J }},\bar{{\boldsymbol{ \mathcal H }}})$ and write the coupled two-moment models in operator form as

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mathsf{c}}}\partial t{\boldsymbol{ \mathcal M }}+\displaystyle \frac{1}{\sqrt{\gamma }}\displaystyle \sum _{i=1}^{3}{\partial }_{{x}^{i}}\left(\sqrt{\gamma }\,{{\boldsymbol{ \mathcal F }}}^{i}({\boldsymbol{ \mathcal M }})\right)={\boldsymbol{ \mathcal G }}({\boldsymbol{ \mathcal M }})+{\boldsymbol{ \mathcal C }}({\boldsymbol{ \mathcal M }}),\end{eqnarray} \tag{ 23 }$$

where ${\boldsymbol{ \mathcal F }}$, ${\boldsymbol{ \mathcal G }}$, and ${\boldsymbol{ \mathcal C }}$ are the position flux operator, geometry source operator, and collision operator, respectively. In particular,

$$\begin{eqnarray}\,\begin{array}{rcl}{\boldsymbol{ \mathcal C }}({\boldsymbol{ \mathcal M }}) & := & \left({\eta }_{{\rm{T}}{\rm\small{OT}}}-{\chi }_{{\rm{T}}{\rm\small{OT}}}\,{ \mathcal J },({\chi }_{{\rm{T}}{\rm\small{OT}}}+{\sigma }_{{\rm\small{IS}}}){\boldsymbol{ \mathcal H }},{\bar{\eta }}_{{\rm{T}}{\rm\small{OT}}}\right.\\ & & -\ \left.{\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}\,\bar{{ \mathcal J }},({\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}+{\bar{\sigma }}_{{\rm\small{IS}}})\bar{{\boldsymbol{ \mathcal H }}}\right).\end{array}\,\end{eqnarray} \tag{ 24 }$$

Neutrino–matter interactions mediate the exchange of lepton number, momentum, and energy between matter and neutrinos. As a first approximation, we assume that the fluid remains static, and momentum exchange is ignored. Under these assumptions, the matter is described by the mass density ρ( x ) (fixed in time), temperature T( x , t), and electron fraction Y_e ( x , t), and neutrino–matter interactions result in changes to the electron fraction,

$$\begin{eqnarray}\,\begin{array}{rcl}\displaystyle \frac{1}{{\mathsf{c}}}\,\displaystyle \frac{{{dY}}_{e}}{{dt}} & = & -\left(\displaystyle \frac{{{\mathsf{m}}}_{{\mathsf{b}}}}{\rho }\right)\left(\displaystyle \frac{1}{{\left({\mathsf{hc}}\right)}^{3}}\right){\displaystyle \int }_{{{\mathbb{R}}}^{+}}\left({\eta }_{{\rm{T}}{\rm\small{OT}}}-{\chi }_{{\rm{T}}{\rm\small{OT}}}\,{ \mathcal J })\right.\\ & & -\ \left.({\bar{\eta }}_{{\rm{T}}{\rm\small{OT}}}-{\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}\,\bar{{ \mathcal J }}\right)\,d{\tilde{V}}_{\varepsilon },\end{array}\,\end{eqnarray} \tag{ 25 }$$

and specific internal energy,

$$\begin{eqnarray}\,\begin{array}{rcl}\displaystyle \frac{1}{{\mathsf{c}}}\,\displaystyle \frac{d\epsilon }{{dt}} & = & -\left(\displaystyle \frac{1}{\rho }\right)\left(\displaystyle \frac{1}{{\left({\mathsf{hc}}\right)}^{3}}\right){\displaystyle \int }_{{{\mathbb{R}}}^{+}}\left({\eta }_{{\rm{T}}{\rm\small{OT}}}-{\chi }_{{\rm{T}}{\rm\small{OT}}}\,{ \mathcal J })\right.\\ & & +\ \left.({\bar{\eta }}_{{\rm{T}}{\rm\small{OT}}}-{\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}\,\bar{{ \mathcal J }}\right)\varepsilon \,d{\tilde{V}}_{\varepsilon },\end{array}\,\end{eqnarray} \tag{ 26 }$$

where m _b is the average baryon mass. The specific internal energy ( x , t) is defined such that ρ is the internal energy density. We note that, given any ρ and Y_e , the specific internal energy is a one-to-one (injective) function of the temperature T; i.e., one can map a given T to a unique , and vice versa. Specifically, ${(\partial \epsilon /\partial T)}_{\rho ,{Y}_{e}}\gt 0$.

Together with the number density evolution equations (see Equation (14a)) for neutrinos and antineutrinos, Equations (25) and (26) lead to the conservation of lepton number

$$\begin{eqnarray}&&\,\begin{array}{l}{\displaystyle \int }_{D}\left[\displaystyle \frac{\rho }{{{\mathsf{m}}}_{{\mathsf{b}}}}\displaystyle \frac{{{dY}}_{e}}{{dt}}+\displaystyle \frac{1}{{\left({\mathsf{hc}}\right)}^{3}}\displaystyle \frac{d}{{dt}}{\displaystyle \int }_{{{\mathbb{R}}}^{+}}\left({ \mathcal J }(\varepsilon )-\bar{{ \mathcal J }}(\varepsilon )\right)\,d{\tilde{V}}_{\varepsilon }\right]d{\boldsymbol{x}}\\ \quad \ =\ -{\oint }_{\partial D}\left[\displaystyle \frac{{\mathsf{c}}}{{\left({\mathsf{hc}}\right)}^{3}}{\displaystyle \int }_{{{\mathbb{R}}}^{+}}({\boldsymbol{ \mathcal H }}-\bar{{\boldsymbol{ \mathcal H }}})\,d{\tilde{V}}_{\varepsilon }\right]d{\boldsymbol{x}}\end{array}\,\end{eqnarray} \tag{ 27 }$$

and energy

$$\begin{eqnarray}&&\,\begin{array}{l}{\displaystyle \int }_{D}\left[\rho \,\displaystyle \frac{d\epsilon }{{dt}}+\displaystyle \frac{1}{{\left({\mathsf{hc}}\right)}^{3}}\displaystyle \frac{d}{{dt}}{\displaystyle \int }_{{{\mathbb{R}}}^{+}}\left({ \mathcal J }(\varepsilon )+\bar{{ \mathcal J }}(\varepsilon )\right)\varepsilon \,d{\tilde{V}}_{\varepsilon }\right]d{\boldsymbol{x}}\\ \quad \ =\ -{\oint }_{\partial D}\left[\displaystyle \frac{{\mathsf{c}}}{{\left({\mathsf{hc}}\right)}^{3}}{\displaystyle \int }_{{{\mathbb{R}}}^{+}}({\boldsymbol{ \mathcal H }}+\bar{{\boldsymbol{ \mathcal H }}})\varepsilon d{\tilde{V}}_{\varepsilon }\right]d{\boldsymbol{x}},\end{array}\,\end{eqnarray} \tag{ 28 }$$

respectively. Note that (ρ/m _b )Y_e = n_e is the electron density (technically the electron minus positron density), and that neutrinos and electrons have lepton number 1, while antineutrinos and positrons have lepton number −1. Equation (27) implies that the total lepton number in domain D (left-hand side of Equation (27)) only changes due to fluxes through the domain boundary ∂D (right-hand side of Equation (27)). A similar conservation statement holds for the total energy in D as given in Equation (28).

We apply the nodal DG discretization (see, e.g., Hesthaven & Warburton 2008 for an overview) to Equation (23) in space ${\boldsymbol{x}}\in {{\mathbb{R}}}^{3}$ and energy $\varepsilon \in {{\mathbb{R}}}^{+}$, in which a logically Cartesian mesh with coordinate aligned elements is considered. To derive the discretized equation from the nodal DG method, we first divide the computational domain $D={D}_{{\boldsymbol{x}}}\times {D}_{\varepsilon }\subset {{\mathbb{R}}}^{3}\times {{\mathbb{R}}}^{+}$ into a disjoint union ${ \mathcal K }$ of open elements K , where each element takes the form

$$\begin{eqnarray}\,\begin{array}{rcl}{\boldsymbol{K}} & = & \{({\boldsymbol{x}},\varepsilon ):{x}^{i}\in {K}^{i}:= ({x}_{{\rm{L}}}^{i},{x}_{{\rm{H}}}^{i}),\\ i & = & 1,2,3,\varepsilon \in {K}^{\varepsilon }:= ({\varepsilon }_{{\rm{L}}},{\varepsilon }_{{\rm{H}}})\},\end{array}\,\end{eqnarray} \tag{ 29 }$$

with ${\rm{\Delta }}{x}^{i}={x}_{{\rm{H}}}^{i}-{x}_{{\rm{L}}}^{i}$ and Δε = ε_H − ε_L denoting the side lengths of K . For i = 1, 2, 3, the spatial surface elements in direction xⁱ are denoted as ${\tilde{{\boldsymbol{K}}}}^{i}={\times }_{j\ne i}{K}^{j}$, with space variables ${\tilde{{\boldsymbol{x}}}}^{i}:= \{{x}^{j}\in {\boldsymbol{x}}:j\ne i\}$ on ${\tilde{{\boldsymbol{K}}}}^{i}$, where × is the Cartesian product operator. We use V _K to denote the proper volume of the element

$$\begin{eqnarray}&&{V}_{{\boldsymbol{K}}}={\int }_{{\boldsymbol{K}}}{dV},\quad \mathrm{where}\quad {dV}=\sqrt{\gamma }\,d{\boldsymbol{x}}\,{{dV}}_{\varepsilon }=\sqrt{\gamma }\,d{\boldsymbol{x}}\,{\varepsilon }^{2}d\varepsilon .\end{eqnarray} \tag{ 30 }$$

We let the approximation space for the DG method, ${{\mathbb{V}}}^{{N}_{{\boldsymbol{x}}},{N}_{\varepsilon }}$, be constructed from the tensor product of one-dimensional polynomials of maximal degrees N _x and N_ε in space and energy, respectively. Note that functions in ${{\mathbb{V}}}^{{N}_{{\boldsymbol{x}}},{N}_{\varepsilon }}$ can be discontinuous across element interfaces. The semidiscrete DG problem is to find ${{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}\in {{\mathbb{V}}}^{{N}_{{\boldsymbol{x}}},{N}_{\varepsilon }}$ (which approximates ${\boldsymbol{ \mathcal M }}$ in Equation (23)) such that (see Cockburn & Shu 2001)

$$\begin{eqnarray}&&\,\begin{array}{l}\displaystyle \frac{1}{{\mathsf{c}}}{\partial }_{t}{\displaystyle \int }_{{\boldsymbol{K}}}{{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}\,v\,{dV}+\displaystyle \sum _{i=1}^{3}{\displaystyle \int }_{{K}^{\varepsilon }}{\displaystyle \int }_{{\tilde{{\boldsymbol{K}}}}^{i}}\left(\sqrt{\gamma }\,{\widehat{{\boldsymbol{ \mathcal F }}}}^{i}{\left.({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}})v\right|}_{{x}_{{\rm{H}}}^{i}}\right.\\ \quad \ -\ \left.\sqrt{\gamma }\,{\widehat{{\boldsymbol{ \mathcal F }}}}^{i}{\left.({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}})v\right|}_{{x}_{{\rm{L}}}^{i}}\,\right)d{\tilde{{\boldsymbol{x}}}}^{i}\,{{dV}}_{\varepsilon }\\ \quad \ -\ \displaystyle \sum _{i=1}^{3}{\displaystyle \int }_{{\boldsymbol{K}}}\left({{\boldsymbol{ \mathcal F }}}^{i}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}})\displaystyle \frac{\partial v}{\partial {x}^{i}}\right){dV}={\displaystyle \int }_{{\boldsymbol{K}}}{\boldsymbol{ \mathcal G }}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}})v\,{dV}\\ \quad \ +\ {\displaystyle \int }_{{\boldsymbol{K}}}{\boldsymbol{ \mathcal C }}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}})v\,{dV}\end{array}\,\end{eqnarray} \tag{ 31 }$$

for all $v\in {{\mathbb{V}}}^{{N}_{{\boldsymbol{x}}},{N}_{\varepsilon }}$ and all ${\boldsymbol{K}}\in { \mathcal K }$. Here the numerical flux approximating the flux on the spatial surface element ${\tilde{{\boldsymbol{K}}}}^{i}$ is denoted as ${\widehat{{\boldsymbol{ \mathcal F }}}}^{i}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}})$. In this work, we consider the Lax–Friedrichs (LF) flux,

$$\begin{eqnarray}\,\begin{array}{rcl}{\widehat{{\boldsymbol{ \mathcal F }}}}^{i}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}){| }_{{x}^{i}} & = & {{\boldsymbol{f}}}^{\mathrm{LF},i}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}{| }_{{x}^{i,-}},{{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}{| }_{{x}^{i,+}})\\ & := & \ \displaystyle \frac{1}{2}\left(\,{{\boldsymbol{ \mathcal F }}}^{i}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}{| }_{{x}^{i,-}})+{{\boldsymbol{ \mathcal F }}}^{i}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}{| }_{{x}^{i,+}})\right.\\ & & -\ \left.{\alpha }^{i}(\,{{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}{| }_{{x}^{i,+}}-{{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}{| }_{{x}^{i,-}}\,\right),\end{array}\,\end{eqnarray} \tag{ 32 }$$

where ${{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}{| }_{{x}^{i,\pm }}={\mathrm{lim}}_{\delta \to 0}{{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}{| }_{{x}^{i}\pm \delta }$ are the evaluations of ${{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}$ at the immediate right/left of xⁱ , which thus are functions of $({\tilde{{\boldsymbol{x}}}}^{i},\varepsilon ,t)$. The parameter ${\alpha }^{i}=| | \mathrm{eig}\left(\partial {{\boldsymbol{ \mathcal F }}}^{i}/\partial {\boldsymbol{ \mathcal M }}\right)| {| }_{\infty }$ is the largest eigenvalue of the flux Jacobian. For massless neutrinos, which propagate at the speed of light, we can take αⁱ = c (i.e., the global LF flux).

In each element K , we approximate the conserved variables ${\boldsymbol{ \mathcal M }}$ by

$$\begin{eqnarray}\,\begin{array}{rcl}{\boldsymbol{ \mathcal M }}({\boldsymbol{x}},\varepsilon ,t) & \approx & {{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}({\boldsymbol{x}},\varepsilon ,t)=\displaystyle \sum _{{n}_{\varepsilon }=1}^{{N}_{\varepsilon }}\displaystyle \sum _{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}={\bf{1}}}^{{{\boldsymbol{N}}}_{{\boldsymbol{x}}}}\\ & & \times \ {{\boldsymbol{ \mathcal M }}}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }}(t){{\ell }}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}({\boldsymbol{x}}){{\ell }}_{{n}_{\varepsilon }}(\varepsilon ),\end{array}\,\end{eqnarray} \tag{ 33 }$$

where ${\{{{\ell }}_{{n}_{\varepsilon }}\}}_{{n}_{\varepsilon }=1}^{{N}_{\varepsilon }}$ is chosen to be a collection of Lagrange polynomials in energy of degree N_ε , n _x $:=$ {n¹, n², n³} is a multi-index that goes from 1 = {1, 1, 1} to N _x = {N _x , N _x , N _x }, and ${\{{{\ell }}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}\}}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}={\bf{1}}}^{{{\boldsymbol{N}}}_{{\boldsymbol{x}}}}$ is the collection of multidimensional polynomials defined as ${{\ell }}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}({\boldsymbol{x}}):= {\prod }_{i=1}^{3}{{\ell }}_{{n}^{i}}({x}^{i})$, with $\{{{\ell }}_{{n}^{i}}\}{}_{{n}^{i}=1}^{{N}_{{\boldsymbol{x}}}}$ one-dimensional Lagrange polynomials in xⁱ of degree N _x . It then follows from these definitions that ${\{{{\ell }}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}({\boldsymbol{x}}){{\ell }}_{{n}_{\varepsilon }}(\varepsilon )\}}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}={\bf{1}},{n}_{\varepsilon }=1}^{{{\boldsymbol{N}}}_{{\boldsymbol{x}}},{N}_{\varepsilon }}$ forms a basis for ${{\mathbb{V}}}^{{N}_{{\boldsymbol{x}}},{N}_{\varepsilon }}$ on K . Motivated by the numerical experiments reported in Bassi et al. (2013), we consider here the Lagrange polynomials with Gauss–Legendre interpolation points (instead of Gauss–Legendre–Lobatto points); e.g., ${{\ell }}_{{n}_{\varepsilon }}$ on interval K^ε is defined as

$$\begin{eqnarray}\,\begin{array}{rcl}{{\ell }}_{{n}_{\varepsilon }}(\varepsilon ) & = & {{\ell }}_{{N}_{\varepsilon },{n}_{\varepsilon }}^{\mathrm{loc}}(\xi ):= \displaystyle \prod _{\displaystyle \genfrac{}{}{0em}{}{j=1}{j\ne {n}_{\varepsilon }}}^{{N}_{\varepsilon }}\displaystyle \frac{\xi -{\xi }_{j}^{\varepsilon }}{{\xi }_{{n}_{\varepsilon }}^{\varepsilon }-{\xi }_{j}^{\varepsilon }},\\ \forall \xi & \in & I:= [-1/2,1/2],\quad {n}_{\varepsilon }=1,\ \ldots ,\ {N}_{\varepsilon }.\end{array}\,\end{eqnarray} \tag{ 34 }$$

The local variable $\xi (\varepsilon ):= \displaystyle \frac{\varepsilon -{\varepsilon }_{{\rm{C}}}}{{\rm{\Delta }}\varepsilon }$ with the center ε_C $:=$ (ε_L + ε_H)/2, and interpolation points $\{{\xi }_{j}^{\varepsilon }\}{}_{j=1}^{{N}_{\varepsilon }}$ are given by the N_ε -point Gauss–Legendre quadrature abscissas on I. On interval Kⁱ , ${{\ell }}_{{n}^{i}}$ is defined analogously as ${{\ell }}_{{n}^{i}}({x}^{i})={{\ell }}_{{N}_{{\boldsymbol{x}}},{n}^{i}}^{\mathrm{loc}}(\xi )$, with Gauss–Legendre interpolation points $\{{\xi }_{j}^{i}\}{}_{j=1}^{{N}_{{\boldsymbol{x}}}}$ on I. With this choice of basis, it follows that the expansion in Equation (33) becomes a nodal representation of ${{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}$, i.e., ${{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}({{\boldsymbol{x}}}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}},{\varepsilon }_{{n}_{\varepsilon }},t)={{\boldsymbol{ \mathcal M }}}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }}(t)$, where ${{\boldsymbol{x}}}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}=\{{x}_{{n}^{1}}^{1},{x}_{{n}^{2}}^{2},{x}_{{n}^{3}}^{3}\}$ and ${\varepsilon }_{{n}_{\varepsilon }}$ are the global space and energy variables corresponding to the local interpolation points $\{{\xi }_{j}^{i}\}{}_{j=1}^{{N}_{{\boldsymbol{x}}}}$ and $\{{\xi }_{j}^{\varepsilon }\}{}_{j=1}^{{N}_{\varepsilon }}$ on element K , respectively. Figure 2 shows an example of nodal DG elements in a reduced space ${\boldsymbol{x}}\in {\mathbb{R}}$ and energy $\varepsilon \in {{\mathbb{R}}}^{+}$ with the interpolation points in both local and global coordinates.

Zoom In Zoom Out Reset image size

We then follow the standard practice and approximate Equation (31) by a semidiscrete system consisting of ${({N}_{{\boldsymbol{x}}})}^{3}{N}_{\varepsilon }$ equations, each on a nodal value ${{\boldsymbol{ \mathcal M }}}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }}$, n _x = 1, ..., N _x , n_ε = 1, ..., N_ε . To derive these equations, we approximate the integrals in Equation (31) using the N _x and N_ε point Gauss–Legendre quadrature rules in space and energy, respectively, with the associated weights ${\{{w}_{{n}^{i}}\}}_{{n}^{i}=1}^{{N}_{{\boldsymbol{x}}}}$ and ${\{{w}_{{n}_{\varepsilon }}\}}_{{n}_{\varepsilon }=1}^{{N}_{\varepsilon }}$ normalized such that ${\sum }_{{n}^{i}=1}^{{N}_{{\boldsymbol{x}}}}{w}_{{n}^{i}}=1$ and ${\sum }_{{n}_{\varepsilon }=1}^{{N}_{\varepsilon }}{w}_{{n}_{\varepsilon }}=1$. Specifically, for $v({\boldsymbol{x}},\varepsilon )={{\ell }}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}({\boldsymbol{x}}){{\ell }}_{{n}_{\varepsilon }}(\varepsilon )$, we have

$$\begin{eqnarray}&&\partial t{\int }_{{\boldsymbol{K}}}{{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}\,v\,{dV}\approx {w}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }}\,{\varepsilon }_{{n}_{\varepsilon }}^{2}\,{\sqrt{\gamma }}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}\,| {\boldsymbol{K}}| \,\partial t{{\boldsymbol{ \mathcal M }}}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }}\end{eqnarray} \tag{ 35 }$$

and

$$\begin{eqnarray}&&\,\begin{array}{l}{\displaystyle \int }_{{\boldsymbol{K}}}\left({\boldsymbol{ \mathcal G }}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}})+{\boldsymbol{ \mathcal C }}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}})\right)v\,{dV}\approx {w}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }}\\ \quad \ \times \ {\varepsilon }_{{n}_{\varepsilon }}^{2}\,{\sqrt{\gamma }}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}\,| {\boldsymbol{K}}| \left({\boldsymbol{ \mathcal G }}{\left({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}\right)}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }}+{\boldsymbol{ \mathcal C }}{\left({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}\right)}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }}\right),\end{array}\,\end{eqnarray} \tag{ 36 }$$

where $| {\boldsymbol{K}}| =({\prod }_{i=1}^{3}{\rm{\Delta }}{x}^{i}){\rm{\Delta }}\varepsilon$, ${w}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }}=({\prod }_{i=1}^{3}{w}_{{n}^{i}}){w}_{{n}_{\varepsilon }}$, and ${\gamma }_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}$ denotes the value of γ at ${{\boldsymbol{x}}}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}$. For the surface integrals and "volume terms," we denote the index ${\tilde{{\boldsymbol{n}}}}_{i}:= \{{n}^{j}\in {{\boldsymbol{n}}}_{{\boldsymbol{x}}}:j\ne i\}$ and obtain, e.g.,

$$\begin{eqnarray}&&\,\begin{array}{l}{\displaystyle \int }_{{K}^{\varepsilon }}{\displaystyle \int }_{{\tilde{{\boldsymbol{K}}}}^{i}}\sqrt{\gamma }\,{\widehat{{\boldsymbol{ \mathcal F }}}}^{i}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}})v{| }_{{x}_{{\rm{H}}}^{i}}\,d{\tilde{{\boldsymbol{x}}}}^{i}\,{{dV}}_{\varepsilon }\approx {w}_{{\tilde{{\boldsymbol{n}}}}_{i},{n}_{\varepsilon }}\\ \quad \ \times \ {\varepsilon }_{{n}_{\varepsilon }}^{2}{\left.\left({\sqrt{\gamma }}_{{\tilde{{\boldsymbol{n}}}}_{i}}{\widehat{{\boldsymbol{ \mathcal F }}}}_{{\tilde{{\boldsymbol{n}}}}_{i},{n}_{\varepsilon }}^{i}{{\ell }}_{{n}^{i}}\right)\right|}_{{x}_{{\rm{H}}}^{i}}\,| {\tilde{{\boldsymbol{K}}}}^{i}| \,| {K}^{\varepsilon }| \end{array}\,\end{eqnarray} \tag{ 37 }$$

and

$$\begin{eqnarray}&&\,\begin{array}{l}{\displaystyle \int }_{{\boldsymbol{K}}}{{\boldsymbol{ \mathcal F }}}^{i}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}})\displaystyle \frac{\partial v}{\partial {x}^{i}}\,{dV}\approx {w}_{{\tilde{{\boldsymbol{n}}}}_{i},{n}_{\varepsilon }}\,{\varepsilon }_{{n}_{\varepsilon }}^{2}\,\displaystyle \sum _{{k}^{i}=1}^{{N}_{{\boldsymbol{x}}}}\\ \quad \ \times \ {w}_{{k}^{i}}{\left.\left({\sqrt{\gamma }}_{{\tilde{{\boldsymbol{n}}}}_{i}}{{\boldsymbol{ \mathcal F }}}_{{\tilde{{\boldsymbol{n}}}}_{i},{n}_{\varepsilon }}^{i}\displaystyle \frac{\partial {{\ell }}_{{k}^{i}}}{\partial {x}^{i}}\right)\right|}_{{x}_{{k}^{i}}^{i}}\,| {\tilde{{\boldsymbol{K}}}}^{i}| \,| {K}^{\varepsilon }| ,\end{array}\,\end{eqnarray} \tag{ 38 }$$

where ∣K^ε ∣ = Δε, $| {\tilde{{\boldsymbol{K}}}}^{i}| ={\prod }_{j\ne i}{\rm{\Delta }}{x}^{j}$, ${w}_{{\tilde{{\boldsymbol{n}}}}_{i},{n}_{\varepsilon }}=({\prod }_{j\ne i}{w}_{{n}^{j}}){w}_{{n}_{\varepsilon }}$, ${\gamma }_{{\tilde{{\boldsymbol{n}}}}_{i}}$ is the evaluation of γ at ${\tilde{{\boldsymbol{x}}}}_{{\tilde{{\boldsymbol{n}}}}_{i}}^{i}$, and ${\widehat{{\boldsymbol{ \mathcal F }}}}_{{\tilde{{\boldsymbol{n}}}}_{i},{n}_{\varepsilon }}^{i}$ and ${{\boldsymbol{ \mathcal F }}}_{{\tilde{{\boldsymbol{n}}}}_{i},{n}_{\varepsilon }}^{i}$ are the evaluations of ${\widehat{{\boldsymbol{ \mathcal F }}}}^{i}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}})$ and ${{\boldsymbol{ \mathcal F }}}^{i}({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}})$ at $({\tilde{{\boldsymbol{x}}}}_{{\tilde{{\boldsymbol{n}}}}_{i}}^{i},{\varepsilon }_{{n}_{\varepsilon }})$, respectively. Here ${\gamma }_{{\tilde{{\boldsymbol{n}}}}_{i}}$, ${\widehat{{\boldsymbol{ \mathcal F }}}}_{{\tilde{{\boldsymbol{n}}}}_{i},{n}_{\varepsilon }}^{i}$, and ${{\boldsymbol{ \mathcal F }}}_{{\tilde{{\boldsymbol{n}}}}_{i},{n}_{\varepsilon }}^{i}$ are functions of xⁱ .

Plugging the terms in Equations (35)–(38) into Equation (31) and dividing through by ${w}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }}\,{\sqrt{\gamma }}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}\,{\varepsilon }_{{n}_{\varepsilon }}^{2}\,| {\boldsymbol{K}}|$ leads to the semidiscrete form of the moment equations when the test function is $v={{\ell }}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}\,{{\ell }}_{{n}_{\varepsilon }}$,

$$\begin{eqnarray}&&\,\begin{array}{l}\displaystyle \frac{1}{{\mathsf{c}}}\partial t{{\boldsymbol{ \mathcal M }}}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }}=-\displaystyle \frac{1}{{\sqrt{\gamma }}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}}\displaystyle \sum _{i=1}^{3}\displaystyle \frac{1}{{w}_{{n}^{i}}{\rm{\Delta }}{x}^{i}}{\left.\left({\sqrt{\gamma }}_{{\tilde{{\boldsymbol{n}}}}_{i}}\,{\widehat{{\boldsymbol{ \mathcal F }}}}_{{\tilde{{\boldsymbol{n}}}}_{i},{n}_{\varepsilon }}^{i}\,{{\ell }}_{{n}^{i}}\right)\right|}_{{x}_{{\rm{L}}}^{i}}^{{x}_{{\rm{H}}}^{i}}\\ \quad \ +\ \displaystyle \frac{1}{{\sqrt{\gamma }}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}}}\displaystyle \sum _{i=1}^{3}\displaystyle \frac{1}{{w}_{{n}^{i}}{\rm{\Delta }}{x}^{i}}\displaystyle \sum _{{k}^{i}=1}^{{N}_{{\boldsymbol{x}}}}{w}_{{k}^{i}}{\left.\left({\sqrt{\gamma }}_{{\tilde{{\boldsymbol{n}}}}_{i}}{{\boldsymbol{ \mathcal F }}}_{{\tilde{{\boldsymbol{n}}}}_{i},{n}_{\varepsilon }}^{i}\displaystyle \frac{\partial {{\ell }}_{{k}^{i}}}{\partial {x}^{i}}\right)\right|}_{{x}_{{k}^{i}}^{i}}\\ \quad \ +\ {\boldsymbol{ \mathcal G }}{\left({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}\right)}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }}+{\boldsymbol{ \mathcal C }}{\left({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}\right)}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }},\end{array}\,\end{eqnarray} \tag{ 39 }$$

for n _x = 1, ..., N _x and n_ε = 1, ..., N_ε in all ${\boldsymbol{K}}\in { \mathcal K }$. Equation (39) defines the spatial and energy discretization of the moment equations and provides the basis for implementation in thornado. Before discussing the time discretization, we further simplify Equation (39) utilizing the structure of the collision operator ${\boldsymbol{ \mathcal C }}$. Here we introduce the collection of space and energy nodes from all elements, denoted as

$$\begin{eqnarray}&&S:= {\cup }_{{\boldsymbol{K}}\in { \mathcal K }}\{({{\boldsymbol{x}}}_{{{\boldsymbol{n}}}_{{\boldsymbol{x}}}},{\varepsilon }_{{n}_{\varepsilon }})\in {\boldsymbol{K}}:{{\boldsymbol{n}}}_{{\boldsymbol{x}}}={\bf{1}},\ \ldots ,\ {{\boldsymbol{N}}}_{{\boldsymbol{x}}},{n}_{\varepsilon }=1,\ \ldots ,\ {N}_{\varepsilon }\},\end{eqnarray} \tag{ 40 }$$

and we define the spatial component of S as S _x and the energy component of S as S_ε , i.e.,

$$\begin{eqnarray}&&{\boldsymbol{x}}\in {S}_{{\boldsymbol{x}}},\,\varepsilon \in {S}_{\varepsilon }\quad \ \mathrm{if}\ \mathrm{and}\ \mathrm{only}\ \mathrm{if}\ \quad ({\boldsymbol{x}},\varepsilon )\in S.\end{eqnarray} \tag{ 41 }$$

A simplified illustration of S _x and S_ε is given in Figure 2(b), in which ${\boldsymbol{x}}\in {\mathbb{R}}$. At each x _k ∈ S _x , the nodal value of ${{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}$ on all ε ∈ S_ε is then denoted as ${{\boldsymbol{ \mathcal M }}}_{{\boldsymbol{k}}}(t):= \{{{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}({{\boldsymbol{x}}}_{{\boldsymbol{k}}},{\varepsilon }_{q},t):{\varepsilon }_{q}\in {S}_{\varepsilon }\}.$ With these notations, we write Equation (39) in the operator form in the remainder of this paper as

$$\begin{eqnarray}&&\,\begin{array}{l}\displaystyle \frac{1}{{\mathsf{c}}}\partial t{{\boldsymbol{ \mathcal M }}}_{{\boldsymbol{k}}}={{\mathsf{F}}}_{{\mathsf{M}}}{\left({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}\right)}_{{\boldsymbol{k}}}+{{\mathsf{G}}}_{{\mathsf{M}}}{\left({{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}\right)}_{{\boldsymbol{k}}}+{{\mathsf{C}}}_{{\mathsf{M}}}({{\boldsymbol{ \mathcal M }}}_{{\boldsymbol{k}}}),\\ \quad \ \forall {\boldsymbol{k}}\ \mathrm{such}\ \mathrm{that}\ {{\boldsymbol{x}}}_{{\boldsymbol{k}}}\in {S}_{{\boldsymbol{x}}},\end{array}\,\end{eqnarray} \tag{ 42 }$$

where F _M , G _M , and C _M denote, respectively, the discrete position space flux operator, the discrete geometry source, and the discrete collision operator. Here the collision term ${{\mathsf{C}}}_{{\mathsf{M}}}({{\boldsymbol{ \mathcal M }}}_{{\boldsymbol{k}}})$ depends only on ${{\boldsymbol{ \mathcal M }}}_{{\boldsymbol{k}}}$ instead of the full discretized solution ${{\boldsymbol{ \mathcal M }}}_{{\rm{h}}}$, since the physical interactions modeled in the collision operator (see Equation (4)) are independent of the position x while coupled in the energy domain. Specifically, let $({{\mathsf{J}}}_{{\bf{k}}},{{\bf\mathsf{H}}}_{{\bf{k}}},{\bar{{\mathsf{J}}}}_{{\bf{k}}},{\bar{{\bf\mathsf{H}}}}_{{\bf{k}}})$ denote the discretized moments ${{\boldsymbol{ \mathcal M }}}_{{\boldsymbol{k}}}$; then it follows from Equation (24) that

$$\begin{eqnarray}&&\,\begin{array}{l}{{\mathsf{C}}}_{{\mathsf{M}}}({{\bf{ \mathcal M }}}_{{\bf{k}}}):= ({\eta }_{{\rm{T}}{\rm\small{OT}}}-{\chi }_{{\rm{T}}{\rm\small{OT}}}\,{{\mathsf{J}}}_{{\bf{k}}},({\chi }_{{\rm{T}}{\rm\small{OT}}}+{\sigma }_{{\rm\small{IS}}}){{\bf\mathsf{H}}}_{{\bf{k}}},\\ {\bar{\eta }}_{{\rm{T}}{\rm\small{OT}}}-{\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}\,{\bar{{\mathsf{J}}}}_{{\bf{k}}},({\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}+{\bar{\sigma }}_{{\rm\small{IS}}}){\bar{{\bf\mathsf{H}}}}_{{\bf{k}}}),\end{array}\,\end{eqnarray} \tag{ 43 }$$

where σ _IS and ${\bar{\sigma }}_{{\rm\small{IS}}}$ are the values of σ _IS and ${\bar{\sigma }}_{{\rm\small{IS}}}$ on the energy nodes S_ε , respectively, and η_TOT , χ_TOT , ${\bar{\eta }}_{{\rm{T}}{\rm\small{OT}}}$, and ${\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}$ are the discrete counterparts of η_TOT , χ_TOT , ${\bar{\eta }}_{{\rm{T}}{\rm\small{OT}}}$, and ${\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}$. These discrete opacities are given by replacing the energy integrals in Equations (15)–(19) with the numerical integrals using S_ε as quadrature points, for example,

$$\begin{eqnarray}\,\begin{array}{rcl}{\eta }_{{\rm{S}}{\rm\small{C}}}({ \mathcal J }) & = & \displaystyle \frac{1}{{\mathsf{c}}{\left({\mathsf{hc}}\right)}^{3}}{\displaystyle \int }_{{{\mathbb{R}}}^{+}}{{\rm{\Phi }}}_{0}^{{\rm{I}}{\rm\small{N}}}(\varepsilon ,\varepsilon ^{\prime} ){ \mathcal J }(\varepsilon ^{\prime} )d{\tilde{V}}_{\varepsilon ^{\prime} }\\ & \approx & \ \displaystyle \frac{1}{{\mathsf{c}}{\left({\mathsf{hc}}\right)}^{3}}\displaystyle \sum _{q=1}^{| {S}_{\varepsilon }| }{w}_{q}^{\varepsilon }{{\rm{\Phi }}}_{0}^{{\rm{I}}{\rm\small{N}}}(\varepsilon ,{\varepsilon }_{q}){{\mathsf{J}}}_{q}=: {\eta }_{{\rm{S}}{\rm\small{C}}}({\mathsf{J}}),\end{array}\,\end{eqnarray} \tag{ 44 }$$

where ${w}_{q}^{\varepsilon }$ denotes the weight associated with energy node ε_q , and ∣S_ε ∣ denotes the total number of energy nodes in S_ε . Specifically, for ε_q ∈ S_ε in some energy element K_ε ,

$$\begin{eqnarray}&&{w}_{q}^{\varepsilon }:= 4\pi {\varepsilon }_{q}^{2}{w}_{{n}_{\varepsilon }}| {K}_{\varepsilon }| ,\end{eqnarray} \tag{ 45 }$$

where ${w}_{{n}_{\varepsilon }}$ is the local Gauss–Legendre weight at ε_q in K_ε defined earlier in this section.

Finally, we apply this same nodal DG discretization on the electron fraction and specific internal energy evolution (Equations (25) and (26)). Augmenting Equation (42) to the resulting semidiscrete equations then gives

$$\begin{eqnarray}&&\,\begin{array}{l}\displaystyle \frac{1}{{\mathsf{c}}}\partial t{{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}={\mathsf{F}}{\left({{\boldsymbol{ \mathcal V }}}_{{\rm{h}}}\right)}_{{\boldsymbol{k}}}+{\mathsf{G}}{\left({{\boldsymbol{ \mathcal V }}}_{{\rm{h}}}\right)}_{{\boldsymbol{k}}}+{\mathsf{C}}({{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}),\\ \forall {\boldsymbol{k}}\ \mathrm{such}\ \mathrm{that}\ {{\boldsymbol{x}}}_{{\boldsymbol{k}}}\in {S}_{{\boldsymbol{x}}},\end{array}\,\end{eqnarray} \tag{ 46 }$$

where ${\boldsymbol{ \mathcal V }}:= ({Y}_{e},\epsilon ,{\boldsymbol{ \mathcal M }})$, ${{\boldsymbol{ \mathcal V }}}_{{\rm{h}}}:= ({Y}_{e,{\rm{h}}},{\epsilon }_{{\rm{h}}},{{\boldsymbol{ \mathcal M }}}_{{\rm{h}}})$ is the fully discretized version of ${\boldsymbol{ \mathcal V }}$, and ${{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}$ denotes the nodal value of ${{\boldsymbol{ \mathcal V }}}_{{\rm{h}}}$ at point x _k ∈ S _x . Here the operators F $:=$ (0, 0, F _M ), G $:=$ (0, 0, G _M ), and C is given by the discrete version of the right-hand sides in Equations (25) and (26) and C _M . Specifically, in C, the energy integrals in Equations (25) and (26) are evaluated using a quadrature with the energy nodes ε_q ∈ S_ε as abscissas and weights defined in Equation (45).

Remark 2. In multidimensional CCSN simulations, the use of curvilinear coordinates, e.g., spherical polar coordinates, suffers from an excessively stringent Courant–Friedrichs–Lewy (CFL) time-step restriction due to singularities at the origin and poles (see, e.g., Müller 2020, Section 2.1.1 and references therein). Although these singularities are not an issue in the spherically symmetric CCSN models considered here, they will become a problem when extending the nodal DG scheme to multidimensional CCSN models. Several techniques have been developed to address this issue, such as mesh coarsening (Skinner et al. 2019), element averaging/merging (Asaithambi & Mahesh 2017; Müller et al. 2019), and spectral filtering (Müller et al. 2019). For future multidimensional simulations, we will consider (i) adopting one of the aforementioned techniques or (ii) using Cartesian coordinates in combination with adaptive mesh refinement.

An IMEX time integration scheme (Ascher et al. 1997; Pareschi & Russo 2005) is considered here for solving the two-moment model in Equation (23). When applied to transport equations with collision terms, IMEX schemes usually handle the collision term with an implicit method while applying an explicit method on the advection term (Hu et al. 2018; see also Just et al. 2015; O'Connor 2015; Kuroda et al. 2016; Skinner et al. 2019 for applications of IMEX-type schemes to neutrino transport). This approach relaxes the excessive time-step restriction from an explicit and stiff collision term and avoids the spatially coupled nonlinear solves from an implicit advection term. While the class of IMEX schemes considered in this paper is detailed in Chu et al. (2019), we include it in the following paragraph for completeness. We also stress that even though the nonlinear solution strategies given in Section 4 are motivated from the implicit part of this class of IMEX schemes, they are general enough to be used with any IMEX scheme that treats the collision term implicitly.

To perform time integration, we discretize the time interval [t₀, t_f] into N time steps, t₀ = t⁰ < t¹ < ... < t^N = t_f, and denote ${{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}({t}^{n})$ as ${{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{n}$, n = 1, ..., N. At each spatial node x _k and time tⁿ , the IMEX scheme integrates the semidiscrete Equation (46) in time via

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{(0)}={{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{n},\end{eqnarray} \tag{ 47a }$$

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{(i)}=\displaystyle \sum _{j=0}^{i-1}{c}_{{ij}}\,{{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{({ij})}+{a}_{{ii}}\,{\mathsf{c}}\,{\rm{\Delta }}t\,{\mathsf{C}}({{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{(i)}\,),\quad i=1,\ \ldots ,\ s,\end{eqnarray} \tag{ 47b }$$

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{n+1}={{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{(s)},\end{eqnarray} \tag{ 47c }$$

where s is the number of stages, the parameters c_ij ≥ 0, ${\sum }_{j=0}^{i-1}{c}_{{ij}}=1$, a_ii > 0, and ${{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{({ij})}$ is given by an explicit update,

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{({ij})}={{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{(j)}+{\hat{c}}_{{ij}}\,{\mathsf{c}}\,{\rm{\Delta }}t\left(\,{\mathsf{F}}{\left({{\boldsymbol{ \mathcal V }}}_{{\rm{h}}}^{(j)}\right)}_{{\boldsymbol{k}}}+{\mathsf{G}}{\left({{\boldsymbol{ \mathcal V }}}_{{\rm{h}}}^{(j)}\right)}_{{\boldsymbol{k}}}\right),\end{eqnarray} \tag{ 48 }$$

with parameters ${\hat{c}}_{{ij}}\geqslant 0$. Thus, in each time step, the s-stage IMEX scheme requires s evaluations of the discrete flux and geometry operators F and G and s inversions of the discrete collision operator C. Here the number of evaluations of F and G is identical to the number of stages due to the fact that, by reusing values of F and G from earlier stages, each stage only requires one additional evaluation of F and G. In general, the inversion of C is the dominant cost in the IMEX scheme. In the remainder of this section, we present the details of the nonlinear system arising from inverting C.

At each stage i of the IMEX scheme, Equation 47(b) can be considered as the nonlinear system

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{(+)}={{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{(* )}+\tau \,{\mathsf{C}}({{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{(+)}),\end{eqnarray} \tag{ 49 }$$

where τ $:=$ a_ii c Δt > 0 denotes the effective time step, ${{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{(* )}$ denotes the weighted sum of explicit updates ${{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{({ij})}$, and ${{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}^{(+)}$ denotes the unknown nodal values to be solved.

Let (Y _{e, k} , _k ) denote the nodal value of (Y_e , ) at x _k , and let $({{\mathsf{J}}}_{{\bf{k}}},{{\bf\mathsf{H}}}_{{\bf{k}}},{\bar{{\mathsf{J}}}}_{{\bf{k}}},{\bar{{\bf\mathsf{H}}}}_{{\bf{k}}})$ denote the discrete moment ${{\boldsymbol{ \mathcal M }}}_{{\boldsymbol{k}}}$, which collects the values of moments $({ \mathcal J },{\boldsymbol{ \mathcal H }},\bar{{ \mathcal J }},\bar{{\boldsymbol{ \mathcal H }}})$ at x _k on all energy nodes in S_ε . Equation (49) can then be considered as a nonlinear system on ${{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}:= ({{\mathsf{Y}}}_{e,{\boldsymbol{k}}},{\epsilon }_{{\boldsymbol{k}}}$, ${{\mathsf{J}}}_{{\bf{k}}},{{\bf\mathsf{H}}}_{{\bf{k}}},{\bar{{\mathsf{J}}}}_{{\bf{k}}},{\bar{{\bf\mathsf{H}}}}_{{\bf{k}}})$. For notational simplicity, we suppress all k subscripts when denoting the nodal values of electron fraction, specific internal energy, and moments in the remainder of this paper. It follows from the definition of C that, at each x _k , the resulting nonlinear system is given by

$$\begin{eqnarray}\,\begin{array}{rcl}{{\mathsf{Y}}}_{e}^{(+)} & = & {{\mathsf{Y}}}_{e}^{(* )}-\tau \,\displaystyle \frac{{{\mathsf{m}}}_{{\mathsf{b}}}}{\rho }\,\displaystyle \sum _{q=1}^{| {S}_{\varepsilon }| }{w}_{q}^{(2)}\left({\eta }_{{\rm{T}}{\rm\small{OT}}}-{\chi }_{{\rm{T}}{\rm\small{OT}}}\,{{\mathsf{J}}}_{q}^{(+)})\right.\\ & & -\ \left.({\bar{\eta }}_{{\rm{T}}{\rm\small{OT}}}-{\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}\,{\bar{{\mathsf{J}}}}_{q}^{(+)}\right),\end{array}\,\end{eqnarray} \tag{ 50a }$$

$$\begin{eqnarray}\,\begin{array}{rcl}{\epsilon }^{(+)} & = & {\epsilon }^{(* )}-\tau \,\displaystyle \frac{1}{\rho }\,\displaystyle \sum _{q=1}^{| {S}_{\varepsilon }| }{w}_{q}^{(3)}\left({\eta }_{{\rm{T}}{\rm\small{OT}}}-{\chi }_{{\rm{T}}{\rm\small{OT}}}\,{{\mathsf{J}}}_{q}^{(+)})\right.\\ & & +\ \left.({\bar{\eta }}_{{\rm{T}}{\rm\small{OT}}}-{\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}\,{\bar{{\mathsf{J}}}}_{q}^{(+)}\right),\end{array}\,\end{eqnarray} \tag{ 50b }$$

$$\begin{eqnarray}&&{{\mathsf{J}}}^{(+)}={{\mathsf{J}}}^{(* )}+\tau \left({\eta }_{{\rm{T}}{\rm\small{OT}}}-{\chi }_{{\rm{T}}{\rm\small{OT}}}\,{{\mathsf{J}}}^{(+)}\right),\end{eqnarray} \tag{ 50c }$$

$$\begin{eqnarray}&&{\bar{{\mathsf{J}}}}^{(+)}={\bar{{\mathsf{J}}}}^{(* )}+\tau \left({\bar{\eta }}_{{\rm{T}}{\rm\small{OT}}}-{\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}\,{\bar{{\mathsf{J}}}}^{(+)}\right),\end{eqnarray} \tag{ 50d }$$

$$\begin{eqnarray}&&{{{\bf\mathsf{H}}}}^{(+)}={{{\bf\mathsf{H}}}}^{(* )}-\tau ({{\chi }}_{{\rm{T}}{\rm\small{OT}}}+{{\sigma }}_{{\rm\small{IS}}}){{{\bf\mathsf{H}}}}^{(+)},\end{eqnarray} \tag{ 50e }$$

$$\begin{eqnarray}&&{\bar{{\bf\mathsf{H}}}}^{(+)}={\bar{{\bf\mathsf{H}}}}^{(* )}-\tau ({\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}+{\bar{\sigma }}_{{\rm\small{IS}}}){\bar{{\bf\mathsf{H}}}}^{(+)}.\end{eqnarray} \tag{ 50f }$$

In Equations 50(a) and (b), a quadrature is used to evaluate the energy integrals in Equations (25) and (26), as described when defining C in Equation (46). Here $({{\mathsf{J}}}_{q},{\bar{{\mathsf{J}}}}_{q})$ denotes the value of $({\mathsf{J}},\bar{{\mathsf{J}}})$ at energy node ε_q ∈ S_ε . The physical constants are absorbed into the weights, i.e.,

$$\begin{eqnarray}&&{w}_{q}^{(2)}:=\displaystyle \frac{1}{{\left({\mathsf{h}}{\mathsf{c}}\right)}^{3}}{w}_{q}^{\varepsilon }\,{\rm{a}}{\rm{n}}{\rm{d}}\,{w}_{q}^{(3)}:=\displaystyle \frac{1}{{\left({\mathsf{h}}{\mathsf{c}}\right)}^{3}}{\varepsilon }_{q}{w}_{q}^{\varepsilon },\end{eqnarray} \tag{ 51 }$$

with ${w}_{q}^{\varepsilon }$ defined in Equation (45).

We take a two-step approach to solve Equation (50), which first solves the fully coupled nonlinear system in Equations 50(a)–(d); plug the solution $({{\mathsf{Y}}}_{e}^{(+)},{\epsilon }^{(+)},{{\mathsf{J}}}^{(+)},{\bar{{\mathsf{J}}}}^{(+)})$ into Equations 50(e)–(f) to compute $({\chi }_{{\rm{T}}{\rm\small{OT}}},{\sigma }_{{\rm\small{IS}}},{\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}},{\bar{\sigma }}_{{\rm\small{IS}}})$; and then update $({{\bf\mathsf{H}}}^{(+)},{\bar{{\bf\mathsf{H}}}}^{(+)})$. We note that, once $({\chi }_{{\rm{T}}{\rm\small{OT}}},{\sigma }_{{\rm\small{IS}}},{\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}},{\bar{\sigma }}_{{\rm\small{IS}}})$ are known, solving Equations 50(e)–(f) is straightforward. Thus, we focus on the solution procedure of the coupled system in Equations 50(a)–(d), where the opacities $({\eta }_{{\rm{T}}{\rm\small{OT}}},{\chi }_{{\rm{T}}{\rm\small{OT}}},{\bar{\eta }}_{{\rm{T}}{\rm\small{OT}}},{\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}})$ are functions of $({{\mathsf{Y}}}_{e},\epsilon ,{\mathsf{J}},\bar{{\mathsf{J}}})$. Specifically, while the opacities are written explicitly as functions of $({\mathsf{J}},\bar{{\mathsf{J}}})$ in Equations (15)–(19), they also depend on the matter state (Y _e , ) through the opacity kernels ${{\rm{\Phi }}}_{0}^{{\rm{I}}{\rm\small{N}}}$, ${{\rm{\Phi }}}_{0}^{{\rm{O}}{\rm\small{UT}}}$, ${{\rm{\Phi }}}_{0}^{{\rm{P}}{\rm\small{R}}}$, and ${{\rm{\Phi }}}_{0}^{{\rm{A}}{\rm\small{N}}}$ in Equations (16)–(19). In a fully implicit approach, these opacities need to be updated in the solution procedure of Equations 50(a)–(d) in order to remain consistent with $({{\mathsf{Y}}}_{e}^{(+)},{\epsilon }^{(+)},{{\mathsf{J}}}^{(+)},{\bar{{\mathsf{J}}}}^{(+)})$.

We proceed to investigate various approaches to solve the coupled nonlinear system in Equations 50(a)–(d).

To simplify the notation, we denote the matter states as ${\boldsymbol{u}}:= ({{\mathsf{Y}}}_{e}^{(+)}$, ⁽⁺⁾), the discretized neutrino and antineutrino distributions as ${\boldsymbol{ \mathcal U }}:= ({{\mathsf{J}}}^{(+)},{\bar{{\mathsf{J}}}}^{(+)})$, and the collection of all of the unknowns as ${\boldsymbol{U}}:= ({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})$ in the remainder of the paper. To formulate the system in Equation (52) as a fixed-point problem, we write it as

$$\begin{eqnarray}&&{\boldsymbol{U}}={\boldsymbol{G}}({\boldsymbol{U}}):= \left(\begin{array}{c}{\bf{g}}({\boldsymbol{ \mathcal U }})\\ {\boldsymbol{ \mathcal G }}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\end{array}\right),\end{eqnarray} \tag{ 54 }$$

where

$$\begin{eqnarray}&&{\bf{g}}({\boldsymbol{ \mathcal U }}):= \left(\begin{array}{c}-\displaystyle \frac{{{\mathsf{m}}}_{{\mathsf{b}}}}{\rho }\displaystyle \sum _{q=1}^{| {S}_{\varepsilon }| }{w}_{q}^{(2)}\left({{\mathsf{J}}}_{q}^{(+)}-{\bar{{\mathsf{J}}}}_{q}^{(+)}\right)+{C}_{{{\mathsf{Y}}}_{e}}\\ -\displaystyle \frac{1}{\rho }\displaystyle \sum _{q=1}^{| {S}_{\varepsilon }| }{w}_{q}^{(3)}\left({{\mathsf{J}}}_{q}^{(+)}+{\bar{{\mathsf{J}}}}_{q}^{(+)}\right)+{C}_{\epsilon }\end{array}\right)\end{eqnarray} \tag{ 55 }$$

with ${C}_{{{\mathsf{Y}}}_{e}}={{\mathsf{Y}}}_{e}^{(* )}$ + $\displaystyle \frac{{{\mathsf{m}}}_{{\mathsf{b}}}}{\rho }{\sum }_{q=1}^{| {S}_{\varepsilon }| }{w}_{q}^{(2)}\left({{\mathsf{J}}}_{q}^{(* )}\right.$ − $\left.{\bar{{\mathsf{J}}}}_{q}^{(* )}\right)$, ${C}_{\epsilon }\,={\epsilon }^{(* )}+\displaystyle \frac{1}{\rho }{\sum }_{q=1}^{| {S}_{\varepsilon }| }{w}_{q}^{(3)}\left({{\mathsf{J}}}_{q}^{(* )}+{\bar{{\mathsf{J}}}}_{q}^{(* )}\right)$, and

$$\begin{eqnarray}&&{\boldsymbol{ \mathcal G }}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }}):= \left(\begin{array}{c}\left({{\mathsf{J}}}^{(* )}+\tau \,{\eta }_{{\rm{T}}{\rm\small{OT}}}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\right)/\left(1+\tau \,{\chi }_{{\rm{T}}{\rm\small{OT}}}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\right)\\ \left({\bar{{\mathsf{J}}}}^{(* )}+\tau \,{\bar{\eta }}_{{\rm{T}}{\rm\small{OT}}}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\right)/\left(1+\tau \,{\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\right)\end{array}\right).\end{eqnarray} \tag{ 56 }$$

Here ${\boldsymbol{ \mathcal G }}$ is an equivalent form of Equations 52(c) and (d), which ensures that G is a contraction map, i.e., the Lipschitz constant of G is strictly less than 1.

The coupled fixed-point algorithm considers Equation (54) as a fixed-point problem with unknowns U ; e.g., applying the Picard iteration on Equation (54) leads to

$$\begin{eqnarray}&&{{\boldsymbol{U}}}^{[k+1]}={\boldsymbol{G}}({{\boldsymbol{U}}}^{[k]}),\end{eqnarray} \tag{ 57 }$$

where U ^[k] denotes the kth iterate of the unknowns U , starting from an initial guess U ^[0]. When G is a contraction mapping, the Picard iteration guarantees that, as k → ∞, the iterate U ^[k] converges to ${\boldsymbol{U}}=({{\mathsf{Y}}}_{e}^{(+)}$, ${\epsilon }^{(+)},{{\mathsf{J}}}^{(+)},{\bar{{\mathsf{J}}}}^{(+)})$, the solution to Equation (52). Here the opacities $(\eta ,\chi ,\bar{\eta },\bar{\chi })$ in G are updated at each iteration k using U ^[k] and thus are consistent with the solution. While the Picard iteration guarantees convergence when G is a contraction, the convergence could be slow. To achieve faster convergence, we implement Anderson acceleration (Anderson 1965; Walker & Ni 2011) to solve Equation (54). Anderson acceleration utilizes information from previous iterations to update the unknowns, which is expected to give faster convergence than the Picard iteration but at a cost of additional memory usage. Specifically, in iteration k, Anderson acceleration on the coupled problem first computes the residual

$$\begin{eqnarray}&&{{\boldsymbol{r}}}^{[k]}:= {\boldsymbol{G}}({{\boldsymbol{U}}}^{[k]})-{{\boldsymbol{U}}}^{[k]},\end{eqnarray} \tag{ 58 }$$

then solves a least-squares problem $\alpha * := {{argmin}}_{\alpha \in {{\mathbb{R}}}^{{m}_{k}+1}}$ $\left\{{\parallel {\sum }_{i=0}^{{m}_{k}}{\alpha }_{i}\,{{\boldsymbol{r}}}^{[k-i]}\parallel }_{2}^{2}\,:{\sum }_{i=0}^{{m}_{k}}{\alpha }_{i}=1\right\}$ with ${m}_{k}:= \min \{m,k\}$, and finally updates

$$\begin{eqnarray}&&{{\boldsymbol{U}}}^{[k+1]}={\sum }_{i=0}^{{m}_{k}}{\alpha }_{i}^{* }\,{\boldsymbol{G}}({{\boldsymbol{U}}}^{[k-i]}).\end{eqnarray} \tag{ 59 }$$

Here the truncation parameter m ≥ 0 is an integer that indicates the "memory" of Anderson acceleration, i.e., the maximum number of residuals kept in memory. When m = 0, the solver reduces to Picard iteration. For m > 0, Anderson acceleration updates U using a linear combination of the last m_k iterates that leads to the minimum residual. In the numerical tests in Section 5, we use m = 2, which we have found to significantly reduce the number of iterations when compared to Picard. For m = 2, the additional memory required for Anderson acceleration is small, since each implicit solve is local in space. In addition, the least-squares problem for ${\alpha }_{i}^{* }$ is small and can be written out explicitly or solved using LAPACK's DGELS.

The other solver we considered for the nonlinear coupled system in Equation (52) is Newton's method, which formulates Equation (52) as a root-finding problem,

$$\begin{eqnarray}&&\left(\begin{array}{c}{\boldsymbol{f}}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\\ {\boldsymbol{ \mathcal F }}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\end{array}\right)=: {\bf{F}}({\boldsymbol{U}})={\bf{0}},\end{eqnarray} \tag{ 60 }$$

where

$$\begin{eqnarray}\,\begin{array}{rcl}{\boldsymbol{f}}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }}) & := & \left(\begin{array}{c}{{\boldsymbol{f}}}_{{{\mathsf{Y}}}_{e}}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\\ {{\boldsymbol{f}}}_{\epsilon }({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\end{array}\right)\\ & := & \ \left(\begin{array}{c}{{\mathsf{Y}}}_{e}^{(+)}+\displaystyle \frac{{{\mathsf{m}}}_{{\mathsf{b}}}}{\rho }\displaystyle \sum _{q=1}^{| {S}_{\varepsilon }| }{w}_{q}^{(2)}\left({{\mathsf{J}}}_{q}^{(+)}-{\bar{{\mathsf{J}}}}_{q}^{(+)}\right)-{C}_{{{\mathsf{Y}}}_{e}}\\ {\epsilon }^{(+)}+\displaystyle \frac{1}{\rho }\displaystyle \sum _{q=1}^{| {S}_{\varepsilon }| }{w}_{q}^{(3)}\left({{\mathsf{J}}}_{q}^{(+)}+{\bar{{\mathsf{J}}}}_{q}^{(+)}\right)-{C}_{\epsilon }\end{array}\right)\end{array}\,\end{eqnarray} \tag{ 61 }$$

with constants ${C}_{{{\mathsf{Y}}}_{e}}$ and C defined as in Equation (55), and

$$\begin{eqnarray}&&\begin{array}{l}{\boldsymbol{ \mathcal F }}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }}):= \left(\begin{array}{c}{{\boldsymbol{ \mathcal F }}}_{{\mathsf{J}}}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\\ {{\boldsymbol{ \mathcal F }}}_{\bar{{\mathsf{J}}}}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\end{array}\right)\\ := \ \left(\begin{array}{c}\left(1+\tau \,{\chi }_{{\rm{T}}{\rm\small{OT}}}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\right){{\mathsf{J}}}^{(+)}-\left({{\mathsf{J}}}^{(* )}+\tau \,{\eta }_{{\rm{T}}{\rm\small{OT}}}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\right)\\ \left(1+\tau \,{\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\right){\bar{{\mathsf{J}}}}^{(+)}-\left({\bar{{\mathsf{J}}}}^{(* )}+\tau \,{\bar{\eta }}_{{\rm{T}}{\rm\small{OT}}}({\boldsymbol{u}},{\boldsymbol{ \mathcal U }})\right)\end{array}\right).\end{array}\end{eqnarray} \tag{ 62 }$$

Applying Newton's method to solve Equation (60) leads to the following update of U in iteration k,

$$\begin{eqnarray}&&{{\boldsymbol{U}}}^{[k+1]}={{\boldsymbol{U}}}^{[k]}+\delta {{\boldsymbol{U}}}^{[k]},\end{eqnarray} \tag{ 63 }$$

where the Newton step is given by

$$\begin{eqnarray}&&\delta {{\boldsymbol{U}}}^{[k]}=-{\boldsymbol{J}}{\left({{\boldsymbol{U}}}^{[k]}\right)}^{-1}\,{\bf{F}}({{\boldsymbol{U}}}^{[k]}),\quad \mathrm{with}\quad {\boldsymbol{J}}({{\boldsymbol{U}}}^{[k]}):= {\bf{F}}^{\prime} ({{\boldsymbol{U}}}^{[k]}).\end{eqnarray} \tag{ 64 }$$

Equation (64) shows that two key components are needed in Newton's method: evaluating the Jacobian J and solving the linear system for the Newton step. For the coupled problem in Equation (60), Jacobian evaluation at a given $\hat{{\boldsymbol{U}}}$ requires computing $\displaystyle \frac{\partial {\boldsymbol{f}}}{\partial {\boldsymbol{u}}}$, $\displaystyle \frac{\partial {\boldsymbol{f}}}{\partial {\boldsymbol{ \mathcal U }}}$, $\displaystyle \frac{\partial {\boldsymbol{ \mathcal F }}}{\partial {\boldsymbol{u}}}$, and $\displaystyle \frac{\partial {\boldsymbol{ \mathcal F }}}{\partial {\boldsymbol{ \mathcal U }}}$ at $(\hat{{\boldsymbol{u}}},\hat{{\boldsymbol{ \mathcal U }}})=: \hat{{\boldsymbol{U}}}$. From Equation (61), it is clear that evaluating $\displaystyle \frac{\partial {\boldsymbol{f}}}{\partial {\boldsymbol{u}}}$ and $\displaystyle \frac{\partial {\boldsymbol{f}}}{\partial {\boldsymbol{ \mathcal U }}}$ at $(\hat{{\boldsymbol{u}}},\hat{{\boldsymbol{ \mathcal U }}})$ is straightforward with minimal cost. However, it follows from Equation (62) that evaluating $\displaystyle \frac{\partial {\boldsymbol{ \mathcal F }}}{\partial {\boldsymbol{u}}}$ and $\displaystyle \frac{\partial {\boldsymbol{ \mathcal F }}}{\partial {\boldsymbol{ \mathcal U }}}$ at $(\hat{{\boldsymbol{u}}},\hat{{\boldsymbol{ \mathcal U }}})$ requires gradients of the opacities $({\eta }_{{\rm{T}}{\rm\small{OT}}},{\chi }_{{\rm{T}}{\rm\small{OT}}},{\bar{\eta }}_{{\rm{T}}{\rm\small{OT}}},{\bar{\chi }}_{{\rm{T}}{\rm\small{OT}}})$ with respect to the matter state u and the neutrino and antineutrino number densities ${\boldsymbol{ \mathcal U }}$, respectively. In particular, Equations (15)–(19) and (44) imply that $\displaystyle \frac{\partial {{\rm{\Phi }}}^{{\rm\small{OP}}}}{\partial {\boldsymbol{u}}}$, which is the gradient of opacity kernels with respect to the matter state, is involved in the computation of $\displaystyle \frac{\partial {\boldsymbol{ \mathcal F }}}{\partial {\boldsymbol{u}}}(\hat{{\boldsymbol{u}}},\hat{{\boldsymbol{ \mathcal U }}})$. As discussed earlier, tabulated opacity kernels are considered in this paper. Thus, we can only obtain an approximate $\displaystyle \frac{\partial {{\rm{\Phi }}}^{{\rm\small{OP}}}}{\partial {\boldsymbol{u}}}$ from the tabulated quantities. Further, when the opacity kernels are not tabulated in terms of u , the gradient $\displaystyle \frac{\partial {{\rm{\Phi }}}^{{\rm\small{OP}}}}{\partial {\boldsymbol{u}}}$ has to be approximated using the chain rule. The detailed calculations of these kernel derivatives are given in the Appendix. Once the approximate Jacobian is obtained, the linear solve in Equation (64) is performed via LAPACK's DGESV.

The next approach we consider is a nested algorithm, which formulates Equation (54) as a nested fixed-point problem with two layers,

$$\begin{eqnarray}&&{\boldsymbol{u}}={\bf{g}}(\hat{{\boldsymbol{ \mathcal U }}}({\boldsymbol{u}})),\end{eqnarray} \tag{ 65a }$$

$$\begin{eqnarray}&&\hat{{\boldsymbol{ \mathcal U }}}({\boldsymbol{u}})={\boldsymbol{ \mathcal G }}({\boldsymbol{u}},\hat{{\boldsymbol{ \mathcal U }}}({\boldsymbol{u}})),\end{eqnarray} \tag{ 65b }$$

where the outer layer, Equation 65(a), is a fixed-point problem on the matter states u , and the inner layer, Equation 65(b), is on the distributions $\hat{{\boldsymbol{ \mathcal U }}}$ for fixed matter states u . These two problems are nested in the sense that evaluating the right-hand side of Equation 65(a) at a given u requires solving Equation 65(b). For example, applying the Picard iteration to both Equations 65(a) and (b) gives the iterative scheme

$$\begin{eqnarray}&&{{\boldsymbol{u}}}^{[k+1]}={\bf{g}}(\hat{{\boldsymbol{ \mathcal U }}}({{\boldsymbol{u}}}^{[k]})),\end{eqnarray} \tag{ 66a }$$

where $\hat{{\boldsymbol{ \mathcal U }}}({{\boldsymbol{u}}}^{[k]})={{\boldsymbol{ \mathcal U }}}^{[k,* ]}$, the limit point of the inner Picard iteration

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal U }}}^{[k,{\ell }+1]}={\boldsymbol{ \mathcal G }}({{\boldsymbol{u}}}^{[k]},{{\boldsymbol{ \mathcal U }}}^{[k,{\ell }]}).\end{eqnarray} \tag{ 66b }$$

In practice, we use Anderson acceleration with m = 2, as described in Section 4.1, to accelerate both the outer and inner solves separately.

The nested approach was considered in Laiu et al. (2020) for relaxing the nonlinear coupling between the electric field and electron concentration when solving implicit systems for semiconductor models. Here the nested approach is motivated by the fact that in solving Equation (54), the most costly part is evaluating the opacity kernels Φ at a given matter state u , which is performed in ${\boldsymbol{ \mathcal G }}$ whenever u is updated. Therefore, while the coupled approach seems simple and straightforward, the nested structure in Equation (65) justifies the additional complexity by reducing the number of updates (on u ) in Equation 65(a) via a more accurate distribution update given by solving Equation 65(b) at the current matter state. Note that the matter state u is fixed in the solution procedure of the inner problem in Equation 65(b), which does not require opacity kernel evaluations and results in much cheaper inner iterations.

The nested Newton's method formulates the inner layer of the nested system in Equation (65) as a root-finding problem, resulting in

$$\begin{eqnarray}&&{\boldsymbol{u}}={\bf{g}}(\hat{{\boldsymbol{ \mathcal U }}}({\boldsymbol{u}})),\end{eqnarray} \tag{ 67a }$$

$$\begin{eqnarray}&&{\boldsymbol{ \mathcal F }}({\boldsymbol{u}},\hat{{\boldsymbol{ \mathcal U }}}({\boldsymbol{u}}))={\bf{0}},\end{eqnarray} \tag{ 67b }$$

where the outer layer in Equation 67(a) is still a fixed-point problem on u , and the inner layer in Equation 67(b) is a root-finding problem on $\hat{{\boldsymbol{ \mathcal U }}}$ for fixed u , with ${\boldsymbol{ \mathcal F }}$ defined in Equation (62). Here Equation 67(a) is solved using Anderson acceleration, and, whenever the right-hand side of Equation 67(a) is evaluated at some given u , the inner problem in Equation 67(b) is solved via Newton's method to obtain $\hat{{\boldsymbol{ \mathcal U }}}({\boldsymbol{u}})$ that is used to evaluate ${\bf{g}}(\hat{{\boldsymbol{ \mathcal U }}}({\boldsymbol{u}}))$. This nested solver is identical to the nested fixed-point algorithm in Section 4.3, except that the inner problem is solved via Newton's method instead of Anderson acceleration. Specifically, let u ^[k] be the kth iterate in the outer layer, then $\hat{{\boldsymbol{ \mathcal U }}}({{\boldsymbol{u}}}^{[k]})={{\boldsymbol{ \mathcal U }}}^{[k,* ]}$, which is the limit point of the Newton iterate

$$\begin{eqnarray}\,\begin{array}{rcl}{{\boldsymbol{ \mathcal U }}}^{[k,{\ell }+1]} & = & {{\boldsymbol{ \mathcal U }}}^{[k,{\ell }]}+\delta {{\boldsymbol{ \mathcal U }}}^{[k,{\ell }]},\quad \mathrm{with}\\ \delta {{\boldsymbol{ \mathcal U }}}^{[k,{\ell }]} & = & -{\left[\displaystyle \frac{\partial {\boldsymbol{ \mathcal F }}}{\partial {\boldsymbol{ \mathcal U }}}({{\boldsymbol{u}}}^{[k]},{{\boldsymbol{ \mathcal U }}}^{[k,{\ell }]})\right]}^{-1}\,{\boldsymbol{ \mathcal F }}({{\boldsymbol{u}}}^{[k]},{{\boldsymbol{ \mathcal U }}}^{[k,{\ell }]}).\end{array}\,\end{eqnarray} \tag{ 68 }$$

Since the matter state u is fixed in the inner iterations, Equations (15)–(19) and (44) imply that the Jacobian $\displaystyle \frac{\partial {\boldsymbol{ \mathcal F }}}{\partial {\boldsymbol{ \mathcal U }}}$ can be calculated with no additional opacity kernel evaluation, which justifies this nested approach. The reason we choose not to formulate the outer layer as a root-finding problem and solve it with Newton's method is to avoid the costly opacity kernel gradient approximation discussed in Section 4.2.

In the numerical tests discussed in this section, the DG-IMEX scheme and the nonlinear solvers are implemented following the specifics below, unless otherwise noted.

• 1.  
  DG scheme—We consider problems with one spatial dimension (imposing spherical symmetry). For the relaxation problems in Section 5.2, we solve the space-homogeneous problem in Equation (49) for a single spatial element. For the proto–neutron star deleptonization problems in Section 5.3, the spatial domain r ∈ [0, 300] km is divided into 128 geometrically progressing elements, with the first element of size Δx = 1 km and the last element of size Δx ≈ 4.54 km. In both tests, the energy domain covering ε ∈ [0, 300] MeV is divided into 16 geometrically progressing elements, where the first element has Δε = 4 MeV and the last element has Δε ≈ 50 MeV. The spatial and energy DG elements considered here are linear (k = 1).
• 2.  
  IMEX scheme—In the proto–neutron star deleptonization problem in Section 5.3, we use the IMEX scheme in Equations (47) and (48) with two stages (s = 2). When written in the so-called Shu–Osher form (as in Equations (47) and (48)), the coefficients are given by Chu et al. (2019):

  $$\begin{eqnarray}\,\begin{array}{rcl}\left[\begin{array}{cc}{c}_{10} & \\ {c}_{20} & {c}_{21}\end{array}\right] & = & \left[\begin{array}{cc}1 & \\ 1/2 & 1/2\end{array}\right],\quad \left[\begin{array}{cc}{\hat{c}}_{10} & \\ {\hat{c}}_{20} & {\hat{c}}_{21}\end{array}\right]=\left[\begin{array}{cc}1 & \\ 0 & 1\end{array}\right],\quad \mathrm{and}\\ \left[\begin{array}{c}{a}_{11}\\ {a}_{22}\end{array}\right] & = & \left[\begin{array}{c}1\\ 1/2\end{array}\right].\end{array}\,\end{eqnarray} \tag{ 69 }$$

  This scheme consists of two evaluations of the explicit part and two implicit solves. We use the realizability-enforcing limiter in Chu et al. (2019) to enforce realizable moments (see Equation (12)) after each stage.
• 3.  
  The EoS and neutrino opacity tables—In all tests, we use a tabulated version of the SFHo EoS (Steiner et al. 2013). Thermodynamic (dependent) variables are tabulated as a function of mass density, temperature, and electron fraction (ρ, T, and Y_e ). The EoS table covers the ranges ρ ∈ [1.66 × 10³, 3.16 × 10¹⁵] g cm⁻³ using N_ρ = 185 points (logarithmically spaced to achieve about 15 points per decade), T ∈ [1.16 × 10⁹, 1.84 × 10¹²] K using N_T = 81 points (logarithmically spaced to achieve about 25 points per decade), and Y_e ∈ [0.01, 0.6] using ${N}_{{Y}_{e}}=30$ points (linearly spaced). The neutrino opacities are taken from Bruenn (1985), with all input thermodynamic quantities computed with the SFHo EoS. The absorption and scattering opacities (χ and σ _IS ) are tabulated in terms of the neutrino energy ε in addition to ρ, T, and Y_e (using the same resolution as the EoS table). The neutrino energy range covers ε ∈ [0.1, 300] MeV using N_ε = 40 logarithmically spaced points. The NES and pair creation and annihilation kernels are tabulated in terms of neutrino energy pairs ε and $\varepsilon ^{\prime}$ (using the same points as for χ and σ _IS ), T (using the same points as in the EoS table), and the degeneracy parameter η = μ_e /(k _B T) ∈ [1 × 10⁻³, 2.5 × 10³] using N_η = 60 logarithmically spaced points. To evaluate dependent variables from the table following, e.g., Mezzacappa & Messer (1999), we use bilinear interpolation (or the higher-dimensional equivalent), while derivatives with respect to any of the independent variables are computed by taking the derivative of the interpolation formula. When interpolating the opacity kernels, we enforce the symmetries in Equation (10).
• 4.  
  Nonlinear solvers—The four nonlinear solvers are implemented following the description in Section 4, with one exception that, in the implementation, the effective emission and absorption opacity χ in Equation (15) is lagged. In other words, when solving Equation (52), χ is evaluated at the starting matter state $({{\mathsf{Y}}}_{e}^{(* )},{\epsilon }^{(* )})$ and is not being updated in the solution procedure. We choose to lag χ in the nonlinear solvers to simplify the Jacobian calculation in Newton's method. Since χ is usually varying slowly in time, lagging χ has minimal impact on the solution accuracy. For the fixed-point solvers, χ can be updated at each iteration at a minor additional cost. As mentioned in Section 4, we choose the truncation parameter in Anderson acceleration to be m = 2, unless otherwise specified. In this case, the least-squares problem for determining α* in Equation (59) becomes an inversion of a 3 × 3 matrix and is thus solved analytically. A numerical justification of this choice of m is given in Section 5.2. In Newton's method, the Jacobian matrix is constructed using the derivatives given in the Appendix.We also note that Jacobian-free Newton–Krylov (JFNK) methods, where the Newton step is computed by solving an approximate Newton system with a Krylov solver, are not well suited for this problem (see, e.g., Knoll & Keyes 2004 for a comprehensive survey of these methods). The reason is that, in the JFNK methods, one evaluation of the opacity kernels is needed in every Krylov iteration to approximate the Jacobian matrix. Thus, the JFNK methods require several opacity evaluations per Newton iteration, while the standard Newton's method only requires one, which makes the JFNK methods more expensive for solving these problems, where the opacity evaluation is a dominant computational cost.
• 5.  
  Nonlinear solver initial guess—When solving the coupled nonlinear system in Equation (52), a natural choice of initial guess for the unknowns $({{\mathsf{Y}}}_{e}^{(+)},{\epsilon }^{(+)},{{\mathsf{J}}}^{(+)},{\bar{{\mathsf{J}}}}^{(+)})$ is $({{\mathsf{Y}}}_{e}^{(* )},{\epsilon }^{(* )},{{\mathsf{J}}}^{(* )},{\bar{{\mathsf{J}}}}^{(* )})$, the weighted sum from explicit steps defined in Equation (49). In this work, we add a "presolve" step that aims to provide a better starting point for the iterative solvers and speed up the computation. Specifically, for a given $({{\mathsf{Y}}}_{e}^{(* )},{\epsilon }^{(* )},{{\mathsf{J}}}^{(* )},{\bar{{\mathsf{J}}}}^{(* )})$, the presolve step solves a subsystem of Equation (52), which is obtained by setting the opacities η_SC, η_TP, χ_SC, and χ_TP in Equation (53) to be zero, i.e., with only emission, absorption, and isoenergetic scattering. The solution of this simplified system then serves as the initial guess for $({{\mathsf{Y}}}_{e}^{(+)},{\epsilon }^{(+)},{{\mathsf{J}}}^{(+)},{\bar{{\mathsf{J}}}}^{(+)})$ in the nonlinear solvers for the system in Equation (52). This presolve step is computationally inexpensive (no opacity table interpolations are needed, since the emission and absorption opacity χ is lagged as discussed in the previous paragraph) while giving a reasonable initial guess for the full system. We observe that, without the presolve step, the coupled and nested Newton's methods in Sections 4.2 and 4.4 could diverge if the initial guess $({{\mathsf{Y}}}_{e}^{(* )},{\epsilon }^{(* )},{{\mathsf{J}}}^{(* )},{\bar{{\mathsf{J}}}}^{(* )})$ is too far away from the solution.
• 6.  
  Nonlinear solver convergence criteria—We set the convergence criteria for the iterative solvers based on the relative residual of the system in Equation (52) at the current iterate. Specifically, for the coupled solvers in Sections 4.1 and 4.2, the convergence criteria are

  $$\begin{eqnarray}&&\,\begin{array}{l}\left|{{\boldsymbol{f}}}_{{{\mathsf{Y}}}_{e}}({{\boldsymbol{u}}}^{[k]},{{\boldsymbol{ \mathcal U }}}^{[k]})\right|\leqslant {\mathtt{tol}}\,\left|{{\mathsf{Y}}}_{e}^{[0]}\right|,\\ \left|{{\boldsymbol{f}}}_{\epsilon }({{\boldsymbol{u}}}^{[k]},{{\boldsymbol{ \mathcal U }}}^{[k]})\right|\leqslant {\mathtt{tol}}\,\left|{\epsilon }^{[0]}\right|,\end{array}\,\end{eqnarray} \tag{ 70a }$$

  $$\begin{eqnarray}&&\,\begin{array}{l}\parallel {{\boldsymbol{ \mathcal F }}}_{{\mathsf{J}}}({{\boldsymbol{u}}}^{[k]},{{\boldsymbol{ \mathcal U }}}^{[k]})\parallel \leqslant {\mathtt{tol}}\,\parallel {{\mathsf{J}}}^{[0]}\parallel ,\\ \parallel {{\boldsymbol{ \mathcal F }}}_{\bar{{\mathsf{J}}}}({{\boldsymbol{u}}}^{[k]},{{\boldsymbol{ \mathcal U }}}^{[k]})\parallel \leqslant {\mathtt{tol}}\,\parallel {\bar{{\mathsf{J}}}}^{[0]}\parallel ,\end{array}\,\end{eqnarray} \tag{ 70b }$$

  where tol > 0 is a constant relative tolerance, and $({{\boldsymbol{f}}}_{{{\mathsf{Y}}}_{e}},{{\boldsymbol{f}}}_{\epsilon },{{\boldsymbol{ \mathcal F }}}_{{\mathsf{J}}},{{\boldsymbol{ \mathcal F }}}_{\bar{{\mathsf{J}}}})$ are the residuals as defined in Equations (61) and (62). As for the nested solvers in Sections 4.3 and 4.4, the outer layer (Equations 65(a) and 67(a)) uses the convergence criteria in Equation 70(a), while the convergence criteria in the inner layer (Equations 65(b) and 67(b)) are given by

  $$\begin{eqnarray}&&\,\begin{array}{l}\parallel {{\boldsymbol{ \mathcal F }}}_{{\mathsf{J}}}({{\boldsymbol{u}}}^{[k]},{{\boldsymbol{ \mathcal U }}}^{[k,{\ell }]})\parallel \leqslant {\mathtt{tol}}\,\parallel {{\mathsf{J}}}^{[k,0]}\parallel ,\\ \parallel {{\boldsymbol{ \mathcal F }}}_{\bar{{\mathsf{J}}}}({{\boldsymbol{u}}}^{[k]},{{\boldsymbol{ \mathcal U }}}^{[k,{\ell }]})\parallel \leqslant {\mathtt{tol}}\,\parallel {\bar{{\mathsf{J}}}}^{[k,0]}\parallel .\end{array}\,\end{eqnarray} \tag{ 71 }$$

  When all convergence criteria are satisfied, the solvers return the current iterate as the solution. In all numerical experiments, we choose the norm to be the discrete L² norm on the energy domain, i.e., $\parallel {\mathsf{J}}\parallel := {\left({\sum }_{q=1}^{| {S}_{\varepsilon }| }{w}_{q}^{(2)}{{\mathsf{J}}}_{q}^{2}\right)}^{1/2}$, and we set the relative tolerance to be tol = 10⁻⁸.

Remark 3. In the coupled fixed-point and the two nested solvers, the basic version of Anderson acceleration outlined in Section 4.1 is implemented. It is known (see, e.g., Walker & Ni 2011 and references therein) that the iterations in Anderson acceleration may suffer from "stagnation" when the least-squares problem is ill conditioned. Due to the choice of small truncation parameter (e.g., m = 2) and the contractive property of the collision operator, we do not encounter the stagnation issue in any of the numerical tests presented in this section. Stagnation may potentially become a practical concern when applying Anderson acceleration to solve more complicated systems, e.g., fully coupled neutrino radiation hydrodynamics. A common approach to mitigate the stagnation issue is to control the condition number of the least-squares problems. This can be achieved by (i) solving a proper reformulation of the least-squares problems using QR factorization (Ni & Walker 2010), (ii) modifying the truncation parameter m adaptively (Yang et al. 2009), and/or (iii) regularizing the least-squares problems (Scieur et al. 2020). Alternatively, one may consider the recently proposed globally convergent variant of Anderson acceleration (Zhang et al. 2020), in which ill conditioning is handled via regularization and global convergence is guaranteed using safeguarding steps.

The first class of test problems we consider here is the relaxation problem, where the neutrino and antineutrino transport Equations 3(a) and (b) are solved with only the collision terms considered, i.e., the space-homogeneous case where the advection terms ${ \mathcal T }(f)$ and ${ \mathcal T }(\bar{f})$ are zero. In these relaxation problems, the collision operator relaxes the distributions f and $\bar{f}$ to the equilibrium Fermi–Dirac distributions f₀ and ${\bar{f}}_{0}$ in Equation (7), respectively, as time evolves. Due to the lack of advection terms, there is no spatial coupling. Thus, these problems can be solved independently in space, which makes them ideal test cases for the nonlinear collision system solvers considered in this paper. In this setup, the semidiscrete moment and matter equations reduce from Equation (46) to

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\mathsf{c}}}\partial t{{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}={\mathsf{C}}({{\boldsymbol{ \mathcal V }}}_{{\boldsymbol{k}}}),\quad \forall {\boldsymbol{k}}\ \mathrm{such}\ \mathrm{that}\ {{\boldsymbol{x}}}_{{\boldsymbol{k}}}\in {S}_{{\boldsymbol{x}}}.\end{eqnarray} \tag{ 72 }$$

For the relaxation problems, we discretize Equation (72) in time with the backward Euler method. At each time step, this time discretization results in a coupled system that takes the form of Equation (50) with the effective time step τ = Δt. To solve this system, we apply the nonlinear solvers in Section 4 to the subsystem in Equation (52) (with τ = Δt) to obtain $({{\mathsf{Y}}}_{e}^{(+)},{\epsilon }^{(+)},{{\mathsf{J}}}^{(+)},{\bar{{\mathsf{J}}}}^{(+)})$, which are then used to compute $({{\bf\mathsf{H}}}^{(+)},{\bar{{\bf\mathsf{H}}}}^{(+)})$. We start with trivial initial states for H and $\bar{{\bf\mathsf{H}}}$, which remain unchanged with time.

We test the nonlinear solvers on the relaxation problem in Equation (72) with two initial matter states that present problems with different degrees of collisionality. The first state, which represents the high-density, strongly collisional region inside a proto–neutron star, is sampled at radius r⁽¹⁾ = 9.756 km from the center of a collapsed stellar core, and the second state, which represents the lower-density regions around the surface of a proto–neutron star with relatively weaker collisionality, is sampled at radius r⁽²⁾ = 39.52 km. We obtained the matter states at these two locations from a spherically symmetric CCSN simulation 50 ms after core bounce (Liebendörfer et al. 2005; results from the vertex code using a 15 M_⊙ progenitor). Specifically, the matter state at r⁽¹⁾ is given by ρ⁽¹⁾ = 1.084 × 10¹⁴ g cm⁻³, T⁽¹⁾ = 1.845 × 10¹¹ K, and ${Y}_{e}^{(1)}=0.2728;$ and the matter state at r⁽²⁾ is ρ⁽²⁾ = 1.032 × 10¹² g cm⁻³, T⁽²⁾ = 8.806 × 10¹⁰ K, and ${Y}_{e}^{(2)}=0.1347$. The inverse mean free paths associated with the neutrino opacities for these matter states are plotted versus neutrino energy in Figure 3. Considering the isoenergetic scattering opacity (which is equal for neutrinos and antineutrinos and increases with neutrino energy as ε²), the mean free path ${\lambda }_{{\rm\small{IS}}}={\sigma }_{{\rm\small{IS}}}^{-1}$ varies from about 10 to about 10⁻⁴ km in the high collisional state (left panel) in the energy range ε ∈ [1, 300] MeV. In the low collisional state, the scattering mean free path varies from λ _IS ≈ 10³ to 10⁻² km. Correspondingly, the collision time τ _IS = λ _IS /c varies from about 3 × 10⁻² to 3 × 10⁻⁷ ms in the high collisional state case and from about 3 to 3 × 10⁻⁵ ms in the low collisional state case.

Zoom In Zoom Out Reset image size

We run the simulations from initial time t₀ = 0 to final time t_f = 0.5 ms for the high-collision case and from t₀ = 0 to final time t_f = 30 ms for the low-collision case. This is to guarantee that the distributions have relaxed to the equilibrium distributions by the end of the simulations. Each simulation is solved on a single spatial element (the relaxation problem is space-homogeneous) with 16 geometrically progressing energy elements that divide the energy domain [0, 300] MeV, where the first element has Δε = 4 MeV and the last element has Δε ≈ 50 MeV. In the simulations, the initial matter states take the electron fraction ${Y}_{e,0}={Y}_{e}^{(i)}$ and specific internal energy ₀ calculated at $({\rho }^{(i)},{T}^{(i)},{Y}_{e}^{(i)})$ from the EoS with i = 1 and 2. The initial neutrino and antineutrino moments are given by

$$\begin{eqnarray}&&\,\begin{array}{r}\begin{array}{l}{{ \mathcal J }}_{0}({r}^{(i)},\varepsilon )={\bar{{ \mathcal J }}}_{0}({r}^{(i)},\varepsilon )=0.99\\ \,\,\times \,\exp \left(-\displaystyle \frac{{\left(\varepsilon -2{{\mathsf{k}}}_{{\mathsf{B}}}{T}^{(i)}\right)}^{2}}{2\times {10}^{2}}\right)\,{\rm{a}}{\rm{n}}{\rm{d}}\\ {{ \mathcal H }}_{0}({r}^{(i)},\varepsilon )={\bar{{ \mathcal H }}}_{0}({r}^{(i)},\varepsilon )=0\end{array}\end{array}\,\end{eqnarray} \tag{ 73 }$$

for i = 1 and 2, which were chosen such that the initial moments are away from the expected final equilibrium distributions while making sure the initial moments are realizable, i.e., satisfy Equation (12).

Figure 4 illustrates the initial and final equilibrium electron neutrino and antineutrino number densities for the two relaxation problems. Here the initial number densities are given in Equation (73), and the final number densities were generated by solving Equation (72) using the nodal DG and backward Euler discretization, together with the proposed nonlinear solvers. We validated these results by comparing them to number densities from simulations on finer energy-temporal meshes, in which no noticeable differences were observed. The number densities for the high-collision case are shown in Figure 4(a), and the ones for the low-collision case are shown in Figure 4(b). For the high-collision case, the temperature is 1.845 × 10¹¹ K at t₀ = 0 and increases by 7.3% to 1.979 × 10¹¹ K at t_f = 0.5 ms, while the electron fraction drops by 14.0% from 0.2728 to 0.2347. For the low-collision case, the temperature rises by 18.8% from 8.806 × 10¹⁰ K at t₀ = 0 to 1.046 × 10¹¹ K at t_f = 30 ms; meanwhile, the electron fraction increases by 2.2% from 0.1347 to 0.1376. From these plots, we observe that, in the high-collision case, the neutrino number density ${ \mathcal J }$ is at least 2 orders of magnitude higher than the antineutrino number density $\bar{{ \mathcal J }}$ at the final equilibrium; thus, they are further away from the initial densities than the ones in the low-collision case. This is one of the reasons that the high collision rate problems are more challenging than the low collision rate problems in the earlier stage, as discussed in the following paragraphs.

Zoom In Zoom Out Reset image size

The conservation of lepton number and energy (see Equations (27) and (28), respectively) in the relaxation tests is shown in Figure 5. Since the relaxation problem is space-homogeneous, the right-hand sides of Equations (27) and (28) are both zero. It then follows from the nonlinear system formulation in Equation (52) that the lepton number and energy are conserved if and only if Equations 52(a) and (b) are satisfied. Thus, the lepton number and energy are expected to be conserved up to the nonlinear solver tolerance at each point on the spacetime grid, which is confirmed in the results reported in Figure 5. In addition, we also observe from Figure 5 that, in both the high and low collision rate cases, tightening the nonlinear solver tolerance from 10⁻⁸ to 10⁻¹² indeed improves the conservation results.

Zoom In Zoom Out Reset image size

Figure 6 shows iteration counts versus time for each nonlinear solver on the two relaxation problems using various time-step sizes: Δt = 10⁻⁴, 10⁻³, and 10⁻² ms. These time-step sizes are motivated by the fact that in the context of the IMEX scheme in Equations (47) and (48), the maximum stable time step is ${\rm{\Delta }}t={{\mathsf{C}}}_{{\rm{CFL}}}\times 3\times {10}^{-3}\left({\rm{\Delta }}x/\text{1 km}\right)$ ms, where C _CFL ≲ 1 is the CFL number. For a spatial resolution ${\rm{\Delta }}x={ \mathcal O }(\text{1 km})$, our chosen time steps bracket what is typically used in CCSN simulations. The results for the high-collision problem are shown in Figure 6(a), while the ones for the low-collision problem are shown in Figure 6(b). In these figures, the top plot shows the "outer" iteration counts for each solver, while the bottom plot shows the averaged "inner" iteration counts for the two nested solvers. Here the outer iteration counts represent the number of iterations needed for solving Equations (54), (60), 65(a), and 67(a) in the coupled AA (Anderson acceleration), coupled Newton, nested AA, and Nested newton solvers, respectively; the inner iteration counts are the number of iterations needed for solving Equations 65(b) and 67(b) in the nested AA and nested Newton solvers, respectively. Since the inner Equations 65(b) and 67(b) are solved in every outer iteration, the reported inner iteration counts in Figures 6(a) and (b) are averaged over the number of times that Equations 65(b) and 67(b) were solved; i.e., the total number of inner iterations taken in the solver is the product of the outer iteration count and the averaged inner iteration count.

Zoom In Zoom Out Reset image size

Before comparing the solvers, we first observe from the results that the implicit system in Equation (52) in the high collision rate case indeed requires more iterations to reach convergence than the same system does in the low collision rate case. Also, the overall iteration count grows as Δt increases. These observations agree with our expectation, since increasing Δt effectively increases the collision rates, and higher collision rates result in stronger coupling of the number densities, which makes the system in Equation (52) harder to solve.

For the nonlinear solver performance, we observe that the coupled Newton solver requires more iterations than the coupled AA solver on harder problems, e.g., problems with higher collision rates, larger time steps, or at an earlier stage (further away from equilibrium); on easier problems, the coupled Newton solver converges in fewer iterations than the coupled AA solver. As expected, the nested solvers indeed reduce the number of outer iterations when compared to the coupled solvers, presumably by providing a more accurate update of the neutrino and antineutrino number densities from the inner iteration. The nested AA and nested Newton solvers share nearly identical outer iteration counts. This is due to the fact that both nested solvers use Anderson acceleration in the outer layer (Equations 65(a) and 67(a)), which takes number densities $\hat{{\boldsymbol{ \mathcal U }}}({\boldsymbol{u}})$ from the inner layers (Equations 65(b) and 67(b)) of the two solvers. Since the problems in the inner layers are equivalent, the solutions $\hat{{\boldsymbol{ \mathcal U }}}({\boldsymbol{u}})$ are identical up to the residual tolerance. The nested AA solver generally requires more (inner) iterations to converge than the nested Newton solver does, especially on harder problems. We also note that the comparison of iteration counts here does not fully reflect the performance of the solvers in terms of computational time. In particular, as discussed in Section 4.3, opacity kernel evaluations make outer iterations much more computationally expensive than inner iterations. In addition, the iterations in Anderson acceleration are computationally cheaper than the ones in Newton's method, since they do not require constructing and inverting the Jacobian matrix. We defer the comparison of computational times to Section 5.3, where a more realistic test problem is considered.

Next, we explore the effect of the truncation parameter m on the convergence of Anderson acceleration by comparing the iteration count for the coupled AA solver with different values of m on both the high and low collision rate relaxation problems considered in the previous test. In this comparison, the time-step size is fixed to Δt = 10⁻³ ms, and m varies from zero to 4, where m = 0 resembles the simple Picard iteration. From the results reported in Figure 7, we observe that, in these problems, Anderson acceleration does converge faster as the value of m increases, while the marginal benefit becomes insignificant for m ≥ 2. This result justifies our choice of m = 2 in the numerical tests throughout the paper, since higher values of m lead to larger memory footprints and more expensive least-squares solves in Anderson acceleration with minimal improvement in the iteration counts.

Zoom In Zoom Out Reset image size

We next investigate how early termination of the iterative solvers affects solution accuracy. We choose to test the nested AA solver with early termination of the outer loop and compare the resulting solution to a converged reference solution at the final time. Here the solution process is terminated either when the convergence criteria (Equation 70(a)) are satisfied or when the outer iteration count for solving Equation 65(a) reaches a preset maximum (MaxIter). We test the solver on the relaxation problem with a low collision rate and time step Δt = 10⁻¹ ms, which is much larger than the usual time step for a stable explicit scheme. The choice of a large time step makes the effect of early termination more pronounced. Figure 8 reports the iteration counts for MaxIter = 1 and 2 on the test problem, along with the electron neutrino and antineutrino (energy-integrated) number densities, temperatures, and electron fractions at the final time. These results are compared to a fully converged reference solution (MaxIter = 100), where the nested fixed-point solver converges well before the nominal maximal outer iteration is reached, as shown in Figure 8(a). Here the presolve step discussed in Section 5.1 is turned off. We note that when setting MaxIter = 1, the nested AA solver (without presolve) resembles an approach used in earlier works such as Just et al. (2015), where the radiation quantities are updated using opacities computed from the lagged matter states; i.e., matter states at time tⁿ are used to update the radiation quantities at tⁿ⁺¹. The results in Figure 8 show that when MaxIter = 1, while the relative differences in the earlier time steps are rather significant (from 30% to 0.5%), the solution still converges to an identical (up to the solver tolerance) equilibrium at the final time. However, when MaxIter = 2, it is clear that the solution converges to a different equilibrium. Indeed, it can be seen from Figure 8(e) that the lepton number and energy are not conserved in the solution with MaxIter = 2, while they are conserved up to the solver tolerance in the fully converged solution (MaxIter = 100). When MaxIter = 2, the changes in lepton number and energy lead to a different equilibrium. As for the solution with MaxIter = 1, despite the fact that the solution process of the nonlinear system (Equation (52)) is terminated early, the lepton number and energy are actually conserved in exact arithmetic. This is because, when MaxIter = 1, the early terminated nested iterates satisfy Equations 52(a) and (b) exactly, which enforces lepton number and energy convergences and leads to the correct equilibrium. However, the early terminated solutions generally do not satisfy Equations 52(c) and (d) and result in rather inaccurate solutions in the transient state. In Figure 9, we repeat the test but with the presolve step turned on. Here we also include the option MaxIter = 0 in the comparison, in which the radiation and matter quantities are only updated in the presolve step, i.e., with the NES and pair processes ignored. From Figure 9, it can be observed that even with MaxIter = 0, the presolve step gives a fairly accurate solution at the final time, as the lepton number and energy are conserved up to the solver tolerance in the presolve step. When the maximum iteration is allowed to be higher, the presolve step improves the solution accuracy, at least in the earlier stage. We observe that while the solution is still inaccurate when MaxIter = 2, it is much closer to the reference solution when the presolve is turned on.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We further investigate and compare the performance of the nonlinear solvers in a more realistic setting with a proto–neutron star deleptonization problem, using initial matter profiles from spherically symmetric CCSN simulations. In this test, we solve the full moment equations presented in Section 2.2 using the DG phase-space discretization in Section 3.1 and the IMEX time integration scheme in Section 3.2.

For this test, we adopt matter profiles from Liebendörfer et al. (2005). Specifically, we use profiles for mass density, temperature, and electron fraction obtained with the vertex code using a 15 M_⊙ progenitor from Woosley & Weaver (1995; model G15 in Liebendörfer et al. 2005). Other thermodynamics quantities (e.g., internal energy and electron, proton, and neutron chemical potentials) are obtained from the tabulated SFHo EoS (Steiner et al. 2013). To investigate the sensitivity to conditions encountered over an extended period covering the neutrino heating phase, we run the comparison on profiles taken at t = 50, 100, 150, and 250 ms after core bounce. The initial matter profiles are plotted in Figure 10.

Zoom In Zoom Out Reset image size

Since we do not have the radiation quantities from Liebendörfer et al. (2005), to initialize the radiation field, we adopt the analytical distribution function from the homogeneous sphere test, f_HS (e.g., Smit et al. 1997), which is a solution to the steady-state transport problem of radiation emanating from a sphere whose constant absorption opacity and emissivity inside a radius R₀ are χ₀ and χ₀ f₀, respectively,

$$\begin{eqnarray}&&{f}_{{\rm{H}}{\rm\small{S}}}[{f}_{0},{\chi }_{0},{R}_{0}](r,\mu )={f}_{0}\left(1-{e}^{-{\chi }_{0}s[{R}_{0}](r,\mu )}\right),\end{eqnarray} \tag{ 74 }$$

where f₀ is an isotropic equilibrium distribution,

$$\begin{eqnarray}&&\,\begin{array}{l}s[{R}_{0}](r,\mu )\\ \ =\ \left\{\begin{array}{ll}r\,\mu +{R}_{0}\,g[{R}_{0}](r,\mu ) & {\rm{if}}\quad r\lt {R}_{0},\mu \in [-1,+1],\\ 2\,{R}_{0}\,g[{R}_{0}](r,\mu ) & {\rm{if}}\quad r\geqslant {R}_{0},\mu \in [{\left(1-{\left({R}_{0}/r\right)}^{2}\right)}^{1/2},+1],\\ 0 & {\rm{otherwise}},\end{array}\right.\end{array}\,\end{eqnarray} \tag{ 75 }$$

and $g{[{R}_{0}]{(r,\mu )=[1-(r/{R}_{0})}^{2}(1-{\mu }^{2})]}^{1/2}$.

To adopt the homogeneous sphere distribution to the current setting, we first estimate the energy-dependent neutrinosphere radius R_ν (ε),

$$\begin{eqnarray}&&{\int }_{{R}_{\nu }(\varepsilon )}^{\infty }\chi (\varepsilon ,r){dr}=\displaystyle \frac{2}{3},\end{eqnarray} \tag{ 76 }$$

where χ is the absorption opacity. (Neutrinosphere radii for the various initial profiles used here are plotted versus neutrino energy in the lower right panel of Figure 10.) The neutrino distribution function is then set to

$$\begin{eqnarray}&&\,\begin{array}{l}f(r,\varepsilon ,\mu )\\ \ =\ \left\{\begin{array}{ll}{f}_{{\rm{H}}{\rm\small{S}}}[{f}_{0}(r,\varepsilon ),\chi (r,\varepsilon ),{R}_{\nu }(\varepsilon )](r,\mu ), & r\lt {R}_{\nu }(\varepsilon )\\ {f}_{{\rm{H}}{\rm\small{S}}}[{f}_{0}({R}_{\nu }(\varepsilon )),\chi ({R}_{\nu }(\varepsilon )),{R}_{\nu }(\varepsilon )](r,\mu ), & r\geqslant {R}_{\nu }(\varepsilon ),\end{array}\right.\end{array}\,\end{eqnarray} \tag{ 77 }$$

where f₀(r, ε) is taken to be the Fermi–Dirac distribution in Equation (7). Finally, the initial moments are computed as

$$\begin{eqnarray}&&\left\{{ \mathcal J },{ \mathcal H }\right\}(r,\varepsilon )=\displaystyle \frac{1}{2}{\int }_{-1}^{1}f(r,\varepsilon ,\mu ){\mu }^{\{\mathrm{0,1}\}}\,d\mu .\end{eqnarray} \tag{ 78 }$$

(The same procedure is adopted to initialize the antineutrinos.)

The evolution of the electron fraction in the deleptonization problem is illustrated in Figure 11, where we plot the electron fraction versus mass density over 10 ms of evolution starting from the 100 ms profile in Figure 10 (dashed lines). Here the evolution of the electron fraction was generated from the DG-IMEX scheme with the proposed nonlinear solvers as described in Section 5.1. We validated the result by comparing it against solutions from simulations on finer spatial energy–temporal meshes, in which no noticeable differences were observed. For higher densities (ρ ≳ 3 × 10¹² g cm⁻³), neutrinos are effectively trapped, and the electron fraction remains largely unchanged. For lower densities, neutrinos (created by electron capture on protons) are not trapped and escape the computational domain, which results in a lowering of the electron fraction (deleptonization). Figure 12 shows conservation of lepton number and energy over 5 ms of evolution in the deleptonization problem starting from the 100 ms postbounce matter profile. Here the lepton number and energy both comprise two parts: (i) the interior lepton number and energy in the computation domain, i.e., the time integrals of the left-hand sides in Equations (27) and (28), respectively, and (ii) the accumulated outflow lepton number and energy at the boundary, which are, respectively, the negations of the right-hand sides in Equations (27) and (28), integrated in time. The conservation results of the lepton number/energy, as well as the evolution of the interior lepton number/energy and the accumulated outflow lepton number/energy, are illustrated in Figure 12(a). The evolution of the individual components of the interior lepton number/energy, including matter lepton number (electron number), neutrino lepton number, internal energy, and neutrino energy, are shown in Figure 12(b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For each profile, the deleptonization problems are simulated from the profile time (50, 100, 150, and 250 ms after core bounce) to 5 ms after the profile time, which are referred to as the initial time t₀ = 0 and the final time t_f = 5 ms in the remainder of the paper. The IMEX time integration scheme discussed in Section 5.1 is used in the simulations, where the time step Δt is determined by the stability requirement for the explicit advection part. These problems are solved in the spatial domain r ∈ [0, 300] km and energy domain [0, 300] MeV, which are divided into 128 and 16 geometrically progressing elements, respectively. Here the first spatial element has Δx = 1 km, the last spatial element has Δx ≈ 4.54 km, the first energy element has Δε = 4 MeV, and the last energy element has Δε ≈ 50 MeV.

Figure 13 shows the iteration counts of the nonlinear solvers on the deleptonization problem for various profiles. The results for profiles taken at 50, 100, 150, and 250 ms after core bounce are illustrated in Figures 13(a)–(d), respectively. As in Figure 6, the top plot in these figures shows the "outer" iteration counts of each solver, while the bottom plot shows the averaged "inner" iteration counts of the nested solvers, where the inner iteration counts are averaged over the number of times that the inner Equations 65(b) or 67(b) were solved. In addition, here both the outer and inner iteration counts are averaged over time (from t₀ to t_f ), and the averaged iteration counts are plotted against the mass density, which corresponds to spatial locations for each profile. (We have found that the number of iterations at a given location varies little from t₀ to t_f .)

Zoom In Zoom Out Reset image size

From Figure 13, we observe that the simulations with all four profiles give consistent results; the nested solvers require fewer outer iterations than the coupled ones do, and the nested AA solver requires more inner iteration to converge than the nested Newton solver does, especially for harder problems (at higher mass density). This observation also agrees with the results reported in Section 5.2 on the relaxation problems. Computational times for these simulations are reported in Table 1, where tests 1, 5, 9, and 13 correspond to simulations shown in Figure 13(a), tests 2, 6, 10, and 14 correspond to simulations shown in Figure 13(b), and so on. In Tables 1 and 2, we report the total computation time t_Tot as well as the detailed timing measurements in each simulation, such as the computational time spent on (i) solving the nonlinear system in Equation (52) in the implicit step (t_Im); (ii) the opacity evaluations/interpolations when solving Equation (52) (t_Op); (iii) linear algebra operations, such as the least-squares solve in Anderson acceleration or the assembly and inversion of Jacobian matrices in Newton's method (t_LA); (iv) the presolve step (t_Ps); (v) the explicit update of the advection term (t_Ex); and (vi) the positivity limiter (t_PL). ⁷ The reported timing results in Tables 1 and 2 are all linearly scaled such that the highest reported computational time (16,085 s from test 2) is scaled to 100, which corresponds to the total computation time when the coupled Newton solver is used to solve the deleptonization problem with the 100 ms profile. For each simulation, we also record the solver configurations, such as the value of the truncation parameter m in Anderson acceleration, the maximum allowed iteration (MaxIter), and whether the presolve step is performed or not. We also note that these reported data are not mutually exclusive, e.g., ${t}_{\mathrm{Tot}}\approx {t}_{\mathrm{Im}}+{t}_{\mathrm{Ex}}+{t}_{\mathrm{PL}}$ and ${t}_{\mathrm{Im}}\approx {t}_{\mathrm{Op}}+{t}_{\mathrm{LA}}+{t}_{\mathrm{Ps}}$. Figure 14 provides a column chart that visualizes the results with the 100 ms profile reported in Table 1. There is no qualitative difference between the results from problems with different profiles. The column chart in Figure 14 confirms that the majority of the computational time in these simulations is spent on opacity evaluation/interpolation, which is proportional to the outer iteration counts. It also shows that the nested solvers indeed speed up the computations by taking inner iterations to reduce the number of outer iterations. Further, we observe that the lower computational cost on the linear algebra operations (t_LA) makes the nested AA solver outperform the nested Newton solver in terms of the total computation time by about 8%, despite the higher inner iteration counts. These results indicate that the nested AA solver leads to the least computation time for the deleptonization problems among all of the tested solvers.

Zoom In Zoom Out Reset image size

Table 1. Overview of Nonlinear Solver Timing Results for the Deleptonization Problem

No.	Solver	Profile	tTot	${t}_{\mathrm{Im}}$	tOp	tLA	tPs	tEx	tPL	m	Presolve	MaxIter
1	Coupled Newton	50 ms	98.8	96.8	78.6	11.7	2.7	0.8	1.0	⋯	Yes	100
2	Coupled Newton	100 ms	$\boxed{100}$	97.9	79.5	11.9	2.7	0.8	1.0	⋯	Yes	100
3	Coupled Newton	150 ms	97.5	95.4	77.5	11.5	2.7	0.8	1.0	⋯	Yes	100
4	Coupled Newton	250 ms	92.9	90.8	73.8	10.8	2.7	0.8	1.0	⋯	Yes	100
5	Coupled AA	50 ms	65.2	63.2	57.6	0.9	2.7	0.8	1.0	2	Yes	100
6	Coupled AA	100 ms	67.8	65.8	60.0	1.0	2.7	0.8	1.0	2	Yes	100
7	Coupled AA	150 ms	69.0	67.0	61.1	1.0	2.7	0.8	1.0	2	Yes	100
8	Coupled AA	250 ms	67.5	65.4	59.7	1.0	2.6	0.8	1.0	2	Yes	100
9	Nested Newton	50 ms	36.6	34.5	26.0	4.8	2.7	0.8	1.0	2	Yes	100
10	Nested Newton	100 ms	36.9	34.9	26.4	4.7	2.7	0.8	1.0	2	Yes	100
11	Nested Newton	150 ms	36.4	34.3	26.1	4.5	2.7	0.8	1.0	2	Yes	100
12	Nested Newton	250 ms	35.7	33.7	25.8	4.3	2.6	0.8	1.0	2	Yes	100
13	Nested AA	50 ms	33.7	31.6	26.0	0.9	2.7	0.8	1.0	2	Yes	100
14	Nested AA	100 ms	34.3	32.2	26.5	0.9	2.7	0.8	1.0	2	Yes	100
15	Nested AA	150 ms	33.8	31.7	26.0	0.9	2.7	0.8	1.0	2	Yes	100
16	Nested AA	250 ms	33.1	31.0	25.5	0.8	2.6	0.8	1.0	2	Yes	100

Notes. The detailed computational time is reported for the results shown in Figure 13. Here the four iterative solvers are compared on four initial matter profiles with identical solver parameters. The results are linearly scaled so that the highest measurement (boxed; 16,085 s) is scaled to 100.

$\begin{eqnarray*}&&\begin{array}{l}{t}_{\mathrm{Tot}}:\ {\rm{total}}\,{\rm{simulation}}\,{\rm{time;}}\ \quad {t}_{\mathrm{Im}}:\ {\rm{implicit}}\,{\rm{solution}}\,{\rm{time;}}\quad {t}_{\mathrm{Op}}:\ {\rm{opacity}}\,{\rm{interpolation}}\,{\rm{time;}}\\ {t}_{\mathrm{LA}}:\ {\rm{dense}}\,{\rm{linear}}\,{\rm{algebra}}\,{\rm{time;}}\quad {t}_{\mathrm{Ps}}:\ {\rm{initial}}\,{\rm{presolve}}\,{\rm{time;}}\quad {t}_{\mathrm{Ex}}:\ {\rm{explicit}}\,{\rm{update}}\,{\rm{time;}}\\ {t}_{\mathrm{PL}}:\ {\rm{positivity}}\,{\rm{limiter}}\,{\rm{time;}}\quad m:\ {\rm{Anderson}}\,{\rm{acceleration}}\,{\rm{truncation}}\,{\rm{parameter;}}\\ {\rm{Presolve:}}\,{\rm{whether}}\,{\rm{the}}\,{\rm{presolve}}\,{\rm{step}}\,{\rm{is}}\,{\rm{performed}}\,{\rm{in}}\,{\rm{solver}}\,{\rm{initialization;}}\\ {\rm{MaxIter:}}\,{\rm{maximum}}\,{\rm{allowed}}\,{\rm{(outer)}}\,{\rm{iteration}}\,{\rm{for}}\,{\rm{(nested)}}\,{\rm{iterative}}\,{\rm{solvers;}}\\ {\rm{Relations:}}\ \quad {t}_{\mathrm{Tot}}\approx {t}_{\mathrm{Im}}+{t}_{\mathrm{Ex}}+{t}_{\mathrm{PL}},{t}_{\mathrm{Im}}\approx {t}_{\mathrm{Op}}+{t}_{\mathrm{LA}}+{t}_{\mathrm{Ps}}.\end{array}\end{eqnarray*}$

Download table as:  ASCIITypeset image

Table 2. Overview of Nonlinear Solver Timing Results for the Deleptonization Problem

#	Solver	Profile	tTot	${t}_{\mathrm{Im}}$	tOp	tLA	tPs	tEx	tPL	m	Presolve	MaxIter
14	Nested AA	100 ms	$\boxed{100}$	94.0	77.1	2.6	7.8	2.4	2.9	2	Yes	100
17	Nested AA	100 ms	117.8	111.9	96.1	2.0	⋯	2.4	2.9	0	No	100
18	Nested AA	100 ms	99.9	94.0	81.5	1.5	⋯	2.4	2.9	1	No	100
19	Nested AA	100 ms	99.7	93.8	81.3	1.5	⋯	2.4	2.9	2	No	100
20	Nested AA	100 ms	118.2	112.3	91.4	3.3	7.7	2.4	2.9	0	Yes	100
21	Nested AA	100 ms	99.3	93.4	76.7	2.6	7.7	2.4	2.9	1	Yes	100
22	Nested AA	100 ms	47.2	41.2	34.0	0.8	⋯	2.4	2.9	2	No	1
23	Nested AA	100 ms	83.9	77.9	67.2	1.2	⋯	2.5	2.9	2	No	2
24	Nested AA	100 ms	100.9	94.9	82.3	1.5	⋯	2.4	2.9	2	No	100
25	Nested AA	100 ms	16.9	10.8	3.7	1.2	7.6	2.5	2.9	2	Yes	0
26	Nested AA	100 ms	54.2	48.2	36.2	1.9	7.7	2.4	2.9	2	Yes	1
27	Nested AA	100 ms	89.1	83.0	68.2	2.2	7.8	2.4	2.9	2	Yes	2

Notes. The detailed computational time is reported for the results shown in Figures 15–17. Here the nested AA solver is tested on the 100 ms profile with various solver parameters. The results are linearly scaled so that the total time in test 14 (boxed; 5518 s) is scaled to 100.

$\begin{eqnarray*}&&\begin{array}{l}{t}_{\mathrm{Tot}}:\ {\rm{total}}\,{\rm{simulation}}\,{\rm{time;}}\ \quad {t}_{\mathrm{Im}}:\ {\rm{implicit}}\,{\rm{solution}}\,{\rm{time;}}\quad {t}_{\mathrm{Op}}:\ {\rm{opacity}}\,{\rm{interpolation}}\,{\rm{time;}}\\ {t}_{\mathrm{LA}}:\ {\rm{dense}}\,{\rm{linear}}\,{\rm{algebra}}\,{\rm{time;}}\quad {t}_{\mathrm{Ps}}:\ {\rm{initial}}\,{\rm{presolve}}\,{\rm{time;}}\quad {t}_{\mathrm{Ex}}:\ {\rm{explicit}}\,{\rm{update}}\,{\rm{time;}}\\ {t}_{\mathrm{PL}}:\ {\rm{positivity}}\,{\rm{limiter}}\,{\rm{time;}}\quad m:\ {\rm{Anderson}}\,{\rm{acceleration}}\,{\rm{truncation}}\,{\rm{parameter;}}\\ {\rm{Presolve:}}\,{\rm{whether}}\,{\rm{the}}\,{\rm{presolve}}\,{\rm{step}}\,{\rm{is}}\,{\rm{performed}}\,{\rm{in}}\,{\rm{solver}}\,{\rm{initialization;}}\\ {\rm{MaxIter:}}\,{\rm{maximum}}\,{\rm{allowed}}\,{\rm{(outer)}}\,{\rm{iteration}}\,{\rm{for}}\,{\rm{(nested)}}\,{\rm{iterative}}\,{\rm{solvers;}}\\ {\rm{Relations:}}\ \quad {t}_{\mathrm{Tot}}\approx {t}_{\mathrm{Im}}+{t}_{\mathrm{Ex}}+{t}_{\mathrm{PL}},{t}_{\mathrm{Im}}\approx {t}_{\mathrm{Op}}+{t}_{\mathrm{LA}}+{t}_{\mathrm{Ps}}.\end{array}\end{eqnarray*}$

Download table as:  ASCIITypeset image

To further analyze the performance of the nested AA solver, we experiment with different solver parameters on the deleptonization problem using the 100 ms profile. Specifically, we tested the nested AA solver with the Anderson acceleration truncation parameter set to m = 0, 1, and 2 for solving the outer system in Equation 65(a). This experiment is to verify the benefit of Anderson acceleration on solving the smaller outer system in Equation 65(a), which has only two unknowns. In addition, for each choice of m, we initialize the solver with and without the presolve step introduced in Section 5.1, which helps us assess the effect of that step in a more realistic setting. The resulting iteration counts are shown in Figure 15, which are averaged over the time from t₀ to t_f as in Figure 13. The computation times are reported in Table 2 (tests 17–21). From these results, we observe that moving from Picard iteration (m = 0) to Anderson acceleration (m > 0) does reduce the outer iteration count and thus improves the computation time by around 18%, especially for harder problems (at higher mass density). However, unlike the results for the coupled AA solver reported in Figure 7, increasing the truncation parameter from m = 1 to 2 does not lead to any observable reduction in either the iteration counts or computation time. This result is not unexpected, since Anderson acceleration is applied here to solve the outer system in Equation 65(a) with only two unknowns, while the results reported in Figure 7 are from solving the fully coupled system in Equation (54). Another observation from Figure 15 is that there is no clear benefit in applying the presolve step on this problem. The presolve step does slightly reduce the iteration counts; however, the additional cost of the presolve step wipes out the gains from fewer iterations. We suspect that the diminished benefit of the presolve step is due to the fact that in the deleptonization problem, the explicit time-step restriction from the advection term limits the stiffness of the implicit problem and forces the implicit update to be small, which makes the effect of presolve insignificant. For deleptonization problems, larger implicit time steps could potentially be achieved by techniques such as subcycling the explicit steps or by adopting the general multirate framework proposed by Sandu (2019). We expect to see a greater impact of the presolve step on these problems with larger time steps, as is observed in the relaxation tests.

Zoom In Zoom Out Reset image size

Finally, we investigate the effect of early termination for the nested AA solver on the deleptonization problems. In these tests, the solver configurations are identical to the ones reported in Figures 8 and 9 for the relaxation problems; e.g., the outer loop of the nested AA solver is terminated early by restricting the maximum number of outer iterations (MaxIter). Here we test the solver on the deleptonization problem with the 100 ms initial matter profile for MaxIter = 1, 2, and 100 and the presolve step discussed in Section 5.1 turned off. We then repeat the test for MaxIter = 0, 1, 2, and 100 with the presolve step turned on. The results without and with the presolve step are shown in Figures 16 and 17, respectively, where the iteration counts are reported along with the electron neutrino and antineutrino (energy-integrated) number densities, temperatures, and electron fractions at the final time t = 5 ms. Here the fully converged solutions (MaxIter = 100) are considered as reference solutions, where the nested solver converges well before the nominal maximal outer iteration is reached, as shown in Figures 16(a) and 17(a). As mentioned in the earlier subsection, the nested AA solver with MaxIter = 1 and the presolve step turned off is effectively lagging the opacities by computing them from the matter states at the previous time step while updating the radiation quantities at the current time step (see, e.g., Just et al. 2015). From Figure 16, we observe that the early terminated solutions are in good agreement with the reference solution. However, allowing two outer iterations does not give a better solution than the one with only a single outer iteration, which is possibly due to the issue of lepton number and energy conservation for early terminated solutions, as discussed in Section 5.2. Another potential reason is that Anderson acceleration does not guarantee monotone decreasing of the residuals (see Toth & Kelley 2015; Kelley 2018; Pollock & Rebholz 2019; Evans et al. 2020 for convergence analysis of AA). The behavior of the residuals from Anderson acceleration alternately increasing and decreasing was observed in Pollock & Rebholz (2019), and a potential explanation is given from Pollock & Rebholz (2019), Theorem 4.5. Similar results can be observed in Figure 17, where the presolve step is in effect. Here the comparison includes the case that MaxIter = 0, where both the radiation and matter quantities are solely updated in the presolve step, with the NES and pair processes omitted. The relative difference in the solution with MaxIter = 0 is mostly around 10⁻²–10⁻³; however, the difference could go up to 10¹ for the energy-integrated antineutrino number density in the high mass density region. We also note that the presolve step reduces the difference in the solution with MaxIter = 2 by a few orders of magnitude in the low mass density region. The reason is that, in the low-density region, the nonlinear solve usually converges within two iterations with the presolve step. The results in Figures 16 and 17 suggest that, for problems requiring few iterations, limiting MaxIter to 1 can give sufficiently accurate solutions while reducing the computation time by roughly a factor of 2 from the fully converged case, as shown in the results for tests 22–27 in Table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size