A Detailed Comparison of Multidimensional Boltzmann Neutrino Transport Methods in Core-collapse Supernovae

Most massive stars ($M\gtrsim 10\,{M}_{\odot }$) end their lives in a cataclysmic core-collapse supernova (CCSN) explosion that releases around ${10}^{51}\,\mathrm{erg}$ of kinetic energy and around ${10}^{53}\,\mathrm{erg}$ of neutrino energy. The iron core begins to collapse when it exceeds its effective Chandrasekhar mass as degenerate electrons are captured onto nuclei and photodissociation breaks apart nuclei (e.g., Bethe 1990). Within a few tenths of a second after the onset of collapse, the inner core becomes very neutron-rich (${Y}_{e}\sim 0.3$) and exceeds nuclear densities ($\sim 2.7\,\times {10}^{14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$). At this point, the strong nuclear force kicks in, dramatically stiffening the equation of state (EOS) and abruptly stopping the collapse of the inner core within a few milliseconds. The inner core then rebounds, sending a shock wave through the supersonically infalling outer core. Neutrino cooling removes energy from the matter under the shock, and photodissociation of heavy nuclei weakens the shock. The shock subsequently stalls at around $150\,\mathrm{km}$ as it lacks pressure support from below to overcome the ram pressure of the accreting outer stellar core.

Understanding the mechanism that revives the shock's outward progress and results in a CCSN is presently the main target of CCSN theory. The canonical theory is the neutrino mechanism (Bethe & Wilson 1985), whereby neutrinos emitted from the dense inner core pass through the matter below the shock, depositing enough thermal energy to revive the shock via thermal support and by driving turbulence (e.g., Janka 2001; Burrows 2013; Müller 2016). However, the strongly nonlinear dynamics in this stage is inexorably coupled to a variety of microphysical processes. In particular, it is very sensitive to the properties of the neutrino field passing through the matter. During the stalled-shock phase, the star delicately straddles the line between explosion and total collapse, so small differences in how the neutrinos interact with the matter can be the difference between an explosion and a dud (e.g., Janka 2001; Murphy & Burrows 2008; O'Connor & Ott 2011; Couch & Ott 2015; Melson et al. 2015a; Burrows et al. 2016).

Computation has become the primary tool for studying these nonlinear processes, as it allows us to see detailed dynamics and make observable predictions (electromagnetic radiation, neutrinos, gravitational waves) under the assumptions imposed by the model. However, computational techniques and resources are still too primitive to allow for simulations complete with all of the required fidelity and involved physics. In general, simulations of CCSNe require a three-dimensional general relativistic (GR) treatment of magnetohydrodynamics, neutrino radiation transport, and a microphysical EOS (e.g., Janka 2012; Kotake et al. 2012; Ott 2016). The simulations also require sufficient resolution or subgrid modeling to capture everything from global dynamics to 100-meter-scale or smaller turbulence (Ott 2016; Radice et al. 2016).

Deep in the inner core, neutrinos are trapped and form an isotropic thermal distribution that slowly diffuses out. Outside the shock, the neutrinos are free streaming and move only radially outward. Though radiation transport methods are constructed to simulate these limits well, the intermediate semitransparent region is challenging to accurately simulate. This region is responsible for most of the dynamics that support the shock's progress due to neutrino heating. In addition, the neutrino opacity scales approximately as the square of the neutrino energy, causing the energy deposition rate and the location of the transition from trapped to free streaming to depend sensitively on neutrino energy. Hence, we require a means of simulating neutrinos of many energies in all regions of a CCSN.

A full treatment of classical neutrino radiation requires evolving the neutrino distribution function of each neutrino species according to the seven-dimensional Boltzmann equation (e.g., Lindquist 1966; Ehlers 1993; Mihalas & Weibel-Mihalas 1999) (three spatial dimensions, three momentum dimensions, time), which presents a significant computational challenge. A wide variety of methods have been used to capture the most important aspects of neutrino transport through the supernova that can be broadly categorized as either phenomenological, deterministic, or probabilistic methods. Though some methods are definitively more accurate than others, there is always a trade-off between efficiency and accuracy.

Phenomenological approaches include the light bulb scheme and neutrino leakage and only very approximately account for neutrino effects. These schemes are very efficient, making them very conducive to parameter studies. In the light bulb scheme (e.g., Bethe & Wilson 1985; Janka & Müller 1996; Ohnishi et al. 2006; Murphy & Burrows 2008), the luminosity and temperature of each neutrino species are simply input parameters. All heating rates are based on this parameter, and cooling rates are estimated based on an approximate optical depth. The inner light bulb boundaries have also been combined with gray transport schemes in the semitransparent and transparent regimes (Scheck et al. 2006). In the leakage scheme, an approximate neutrino optical depth at each point is calculated, and this is used to set the cooling rate at each point (Ruffert et al. 1996; Rosswog & Liebendörfer 2003). Neutrino heating can be included approximately by assuming that the neutrino luminosity through a given point is determined by the energy leaking radially outward from below (O'Connor & Ott 2010; Fernández & Metzger 2013; Perego et al. 2016).

Approximate deterministic methods solve a simplified version of the Boltzmann equation in order to make a more tractable problem. The isotropic diffusion source approximation method evolves an isotropic trapped component and a free-streaming component of the distribution function (Liebendörfer et al. 2009). In truncated moment methods, the distribution function is discretized into an infinite list of angular moments, only the first few of which are directly evolved. Flux-limited diffusion (FLD; Levermore & Pomraning 1981; Mihalas & Klein 1982; Pomraning 1983; Castor 2004; Krumholz et al. 2007) is a one-moment method that evolves only the zeroth moment of the distribution function (energy density) and requires a closure relation to estimate the first moment (flux). However, FLD fails to capture much of the angular information about the distribution function and tends to smooth out angular variations (e.g., Janka 1992; Burrows et al. 2000; Liebendörfer et al. 2004; Ott et al. 2008; Zhang et al. 2013). In the two-moment method (Pomraning 1969; Anderson & Spiegel 1972; Thorne 1981; Dubroca & Feugeas 1999; Audit et al. 2002; Shibata et al. 2011; Vaytet et al. 2011; Cardall et al. 2013), the zeroth and first moments are evolved, and a closure relation is required to estimate the second moment (pressure tensor) and complete the system of equations. The closure can be provided by some ad hoc analytical function (e.g., Smit et al. 2000; Murchikova et al. 2017 and references therein). This is also known as the M1 method, though M1 confusingly refers to a specific closure as well (Levermore 1984; Dubroca & Feugeas 1999). The closure can also be more accurately determined using a direct solution of the Boltzmann equation, referred to as a variable Eddington tensor (VET) method (e.g., Stone et al. 1992; Hayes & Norman 2003).

In discrete ordinates (DO or S_n) methods, the distribution function is discretized into angular bins, each of which is directly evolved (e.g., Pomraning 1969; Mihalas & Weibel-Mihalas 1999 and references therein). In spherical harmonic (P_n) methods, the distribution function is decomposed in terms of spherical harmonics, and a small number of these are evolved (e.g., Pomraning 1973; Radice et al. 2013). Finally, fully spectral methods in all six space-momentum dimensions have been applied to stationary neutrino transport calculations (Peres et al. 2014).

Monte Carlo (MC) radiation transport (Fleck & Cummings 1971; Fleck & Canfield 1984; Densmore et al. 2007; Abdikamalov et al. 2012) is a probabilistic method that samples the trajectories of a finite number of individual neutrinos and assumes their behavior is representative of the rest of the bulk neutrino behavior. Tubbs (1978) applied MC methods to neutrino transport for the first time to study neutrino–matter equilibration in an infinite, uniform medium. MC transport has been long used in 1D steady-state transport calculations (Janka & Hillebrandt 1989; Janka 1991, 1992; Yamada et al. 1999; Keil et al. 2003; Abdikamalov et al. 2012), though the code of Abdikamalov et al. (2012) was also designed to perform time-dependent calculations. Richers et al. (2015) performed steady-state MC transport calculations on 2D snapshots of accretion disks from neutron star mergers, though the optical depths were much lower than in the CCSN context.

There is much more to a radiation transport code than the broad classes of methods mentioned above. Detailed Boltzmann transport using either DO or VET methods in dynamical simulations has been achieved in one (Mezzacappa & Bruenn 1993; Yamada et al. 1999; Burrows et al. 2000; Rampp & Janka 2002; Roberts 2012) and two (Livne et al. 2004; Ott et al. 2008; Nagakura et al. 2017a) spatial dimensions, but three-dimensional calculations are presently only possible when the fluid is assumed to be stationary (Sumiyoshi & Yamada 2012; Sumiyoshi et al. 2015). The local two-moment method is used in 1D (O'Connor 2015), 2D (Just et al. 2015a; O'Connor & Couch 2015), and 3D (Sekiguchi et al. 2015; Foucart et al. 2016a; Kuroda et al. 2016; Roberts et al. 2016, see also Müller & Janka 2015) core-collapse and neutron star merger simulations. FLD is also popular in 2D simulations (Dessart et al. 2006; Swesty & Myra 2009; Zhang et al. 2013). Various versions of the ray-by-ray (RbR) approximation can also be used to extend a 1D transport method to two (Burrows et al. 1995; Buras et al. 2006; Müller et al. 2010; Suwa et al. 2010; Pan et al. 2016) or three (Takiwaki et al. 2012; Hanke et al. 2013; Lentz et al. 2015; Melson et al. 2015b) spatial dimensions in an efficient manner by making transport along individual radial rays nearly independent of other rays or by solving a single spherically averaged 1D transport problem. Fully general relativistic neutrino radiation hydrodynamics simulations are now possible (Janka 1991; Yamada et al. 1999; Liebendörfer et al. 2001; Sumiyoshi et al. 2005; Müller et al. 2010; Kuroda et al. 2012, 2016; Shibata & Sekiguchi 2012; O'Connor 2015; Foucart et al. 2016a; Roberts et al. 2016), and many codes incorporate general relativistic effects with various levels of approximation (e.g., Müller et al. 2010; O'Connor & Couch 2015; Skinner et al. 2016). Special relativistic effects can be accounted for in full generality (Müller et al. 2010; Shibata & Sekiguchi 2012; Richers et al. 2015; Foucart et al. 2016a; Nagakura et al. 2017a) or by using only up to ${ \mathcal O }(v/c)$ terms (e.g., Rampp & Janka 2002; Dolence et al. 2015; Just et al. 2015a; Lentz et al. 2015; Skinner et al. 2016). Even if the transport method is equivalent, simulations differ in how neutrino–matter and neutrino–neutrino interactions are treated (e.g., Lentz et al. 2012a). In short, each piece of relevant physics can be simulated accurately, but 3D simulations containing all pieces remain a goal that has not yet been achieved.

Any computational method is an approximation of reality, and every method has strengths and weaknesses. It is therefore expected that computations performed by different codes should arrive at different solutions, though they should converge to the physical answer with increasing simulation fidelity. Understanding the weaknesses of a given method is a prerequisite to interpreting the physical meaning of simulation results. It is standard practice to test that codes produce known solutions to simple problems and to perform self-convergence tests to ensure that results are not mistakes or numerical artifacts. However, even with these practices in place, different codes produce different results, and independent verification is required to help determine which features of each are realistic (Calder et al. 2002).

Several works in the past have evaluated the accuracy of a low-order method like FLD or two-moment transport by comparing with a high-order DO, MC, or VET method (e.g., Janka 1992; Mezzacappa & Bruenn 1993; Messer et al. 1998; Burrows et al. 2000; Liebendörfer et al. 2004; Ott et al. 2008). However, in these comparisons, the low-order method is not expected to converge to the same result as the high-order method, which does not help to verify the high-order method. There have been few detailed comparisons between high-order methods that solve the same full Boltzmann equation (i.e., VET, DO, and MC methods), and none in more than one spatial dimension. Yamada et al. (1999) compared the results of a new DO implementation to the MC code of Janka (1991) in 1D GR snapshots of CCSN simulations, but they ignored fluid motion. Liebendörfer et al. (2005) performed a comparison of dynamical 1D CCSN simulations using the GR DO code Agile-BOLTZTRAN (Liebendörfer et al. 2004) with the Newtonian VET code VERTEX-PROMETHEUS (Rampp & Janka 2002). They found very good agreement once an effective potential was introduced to VERTEX-PROMETHEUS to account for GR effects (Marek et al. 2006). Several groups have since used the results of Liebendörfer et al. (2005) as a standard for comparison (e.g., Sumiyoshi et al. 2005; Müller et al. 2010; Suwa et al. 2011, 2016; Lentz et al. 2012a, 2012b; O'Connor & Ott 2013; Just et al. 2015b; O'Connor 2015). With the recent arrival of the many advanced multidimensional neutrino radiation hydrodynamics codes mentioned previously, continued independent verification is essential to interpreting simulation results.

In this paper, we perform the first detailed multidimensional comparison between fully special relativistic Boltzmann neutrino transport codes using a DO neutrino radiation hydrodynamics code (Nagakura et al. 2017a; hereafter NSY) and the MC radiation transport code Sedonu (Richers et al. 2015). We make the time-independent comparisons on spherically symmetric (1D) and in axisymmetric (2D) snapshots from CCSN simulations at around $100\,\mathrm{ms}$ after core bounce. Both codes are carefully configured to calculate the full steady-state neutrino distribution function from first principles in as similar a manner as possible. We find remarkably good agreement in all spectral, angular, and fluid interaction quantities, lending confidence to both methods. The MC method predicts sharper angular features due to the effectively infinite angular resolution, but struggles to drive down noise in quantities where subtractive cancellation is prevalent (e.g., net gain within the protoneutron star and off-diagonal components of the Eddington tensor). We test the importance of accounting for fluid velocities by setting all velocities to zero and find that the differences induced are much smaller than the errors due to finite momentum-space resolution. We compare directly computed second angular moments to those predicted by popular two-moment closures, and we find that the error from the approximate closure is comparable to the difference between the DO and MC methods.

This paper is organized as follows. In Section 2, we review the DO, MC, and two-moment methods. We present the results of the transport method comparisons in spherical symmetry in Section 4 and in axial symmetry in Section 5. We summarize our conclusions in Section 6. The MC transport code Sedonu is publicly available at https://bitbucket.org/srichers/sedonu, and the results obtained in this study from both transport codes are available at https://stellarcollapse.org/MCvsDO.

The transport of classical neutrinos is described in general by the Boltzmann equation (e.g., Lindquist 1966; Ehlers 1993; Mihalas & Weibel-Mihalas 1999):

$$\begin{eqnarray}&&\displaystyle \frac{{df}}{d\lambda }={\left[\displaystyle \frac{\partial f}{d\lambda }\right]}_{\mathrm{em} \mbox{-} \mathrm{abs}}+{\left[\displaystyle \frac{\partial f}{d\lambda }\right]}_{\mathrm{scat}}+{\left[\displaystyle \frac{\partial f}{d\lambda }\right]}_{\mathrm{pair}},\end{eqnarray} \tag{ 1 }$$

where λ is an affine parameter. In a coordinate basis, the geodesic equation gives

$$\begin{eqnarray}&&\displaystyle \frac{{df}}{d\lambda }={p}^{\alpha }\displaystyle \frac{\partial f}{\partial {x}^{\alpha }}-{{\rm{\Gamma }}}_{\mu \nu }^{\alpha }{p}^{\mu }{p}^{\nu }\displaystyle \frac{\partial f}{\partial {p}^{\alpha }},\end{eqnarray} \tag{ 2 }$$

where ${p}^{\alpha }$ are the neutrino four-momenta, ${x}^{\alpha }$ are the four spacetime coordinates, and ${{\rm{\Gamma }}}_{\mu \nu }^{\alpha }$ are Christoffel symbols of the spacetime metric. In an orthogonal coordinate system in flat spacetime, the distribution function f of each neutrino species is defined as

$$\begin{eqnarray}&&f({\boldsymbol{x}},{\boldsymbol{\Omega }},\epsilon ,t)=\displaystyle \frac{{dN}}{{dV}\,d{\boldsymbol{\Omega }}\,d\epsilon },\end{eqnarray} \tag{ 3 }$$

where N is the number of neutrinos, ${dV}({\boldsymbol{x}})$ is the volume element at position ${\boldsymbol{x}},t$ is the time, ${\boldsymbol{\Omega }}$ is the neutrino direction, and is the neutrino energy. The three source terms on the right-hand side of Equation (1) are interaction terms from emission and absorption, scattering, and pair processes, respectively. The neutrino propagation is generally calculated in reference to coordinates defined in the lab frame, but interactions between matter and neutrinos are formulated in the fluid rest frame (a.k.a. the comoving frame). It is thus important to very carefully keep track of the frame in which various quantities are defined. This is consistent with widely used conventions in the relativistic neutrino transport community (e.g., Shibata et al. 2011; Cardall et al. 2013).

The interaction terms are all local and formulated in a frame comoving with the underlying fluid. The emission and absorption terms are the simplest, as they depend linearly on the distribution function of a given neutrino species. The scattering term in general depends quadratically on the distribution function of a given species, since neutrinos are fermions and the reaction is inhibited by final-state neutrino blocking. However, under the assumption of isoenergetic scattering, it reduces to a linear dependence (see Appendix D). The pair term, which includes neutrino pair annihilation and creation and neutrino bremsstrahlung, depends on the product of the distribution functions of the species and anti-species. Our static transport calculations solve for the f that satisfies $\partial f/\partial t=0$.

There are six species of neutrinos in the standard model corresponding to the six leptonic species (${\nu }_{e},{\nu }_{\mu },{\nu }_{\tau }$ and their antiparticles). Electron neutrinos and antineutrinos interact with nucleons via both charged-current and neutral-current processes, while heavy lepton neutrinos interact only via neutral current processes. This makes each heavy lepton neutrino species individually less impactful than electron anti/neutrinos and makes all four species behave very similarly. In light of this, we simulate ${\nu }_{e}$ and ${\bar{\nu }}_{e}$ individually, but group all of the heavy lepton neutrinos into a single simulated species ${\nu }_{x}$ for computational efficiency.

Neutrino interaction rates depend on the properties of the fluid through which they traverse. In this study, we use the nonhyperonic EOS of Shen et al. (2011) to determine the abundances and chemical potentials of each constituent (i.e., leptons, nucleons, and nuclei) given the fluid density, temperature, and electron fraction. We consider the following minimum but essential sets of neutrino–matter interactions in the postbounce phase of CCSNe:

$$\begin{eqnarray*}{e}^{-}+p & \longleftrightarrow & {\nu }_{e}+n,\\ {e}^{+}+n & \longleftrightarrow & \bar{{\nu }_{e}}+p,\\ \nu (\bar{\nu })+N & \longleftrightarrow & \nu (\bar{\nu })+N,\\ {e}^{-}+{e}^{+} & \longleftrightarrow & \nu +\bar{\nu },\\ N+N & \longleftrightarrow & N+N+\nu +\bar{\nu },\end{eqnarray*}$$

where $N\in \{n,p\}$. From top to bottom, these processes are electron capture by free protons, positron capture by free neutrons, isoenergetic scattering with nucleons, electron–positron pair annihilation, and nucleon–nucleon bremsstrahlung, along with each of their inverse reactions.

Though a multitude of phenomenological, approximate, and exact transport methods exist in the literature, we will focus on three of them. The DO (Section 2.1) and MC (Section 2.2) methods described here both solve the Boltzmann equation directly in all three momentum dimensions and multiple spatial dimensions and so should converge to the same physical result. We also investigate how well approximate closure relations in the two-moment method (Section 3) compare to the solutions computed by DO and MC calculations.

2.1. Discrete Ordinates

The DO Boltzmann code of Nagakura, Sumiyoshi, and Yamada (hereafter NSY) is a grid-based multidimensional neutrino radiation hydrodynamics code that solves the conservative form of Equation (1) in the language of the 3 + 1 formulation of GR. The numerical method is essentially the same as described by Sumiyoshi & Yamada (2012), though it has since been extended to account for special relativistic effects and has been coupled with Newtonian hydrodynamics (Nagakura et al. 2014, 2017a). The newest version of this code was recently applied to axisymmetric CCSN simulations in Nagakura et al. (2017b).

The neutrino distribution function f is discretized onto a spherical-polar spatial grid described by radius r, polar angle θ, and azimuthal angle ϕ. The radial grid is constructed so as to provide good resolution where the density gradient is large. The radial mesh spacing is set to ${\rm{\Delta }}r=300\,{\rm{m}}$ at the center and decreases with increasing radius up to the location of the steepest density gradient at $r=10\,\mathrm{km}$, where ${\rm{\Delta }}r=104\,{\rm{m}}$. For $r\geqslant 10\,\mathrm{km}$, the spacing increases by $1.7 \%$ per zone up to $r=500\,\mathrm{km}$. For $r\geqslant 500\,\mathrm{km}$, the spacing increases by $3.8 \%$ per zone up to the outer boundary of $r=5000\,\mathrm{km}$. This results in 384 radial grid zones over the entire domain. The spatial angular grid is set to 128 Gaussian quadrature points in the range $0\leqslant \theta \leqslant \pi$. At each spatial location, f is discretized onto a spherical-polar momentum-space grid described by neutrino energy , neutrino polar angle $\bar{\theta }$ (where $\bar{\theta }=0$ is in the radial direction), and neutrino azimuthal angle $\bar{\phi }$. The first bin of the neutrino energy grid extends over 0–2 MeV in the 1D_1x and 2D calculations, 0–1 MeV in the 1D_2x calculations, and 0–0.5 MeV in the 1D_4x calculations. The rest of the energy bins are logarithmically spaced from 2 to $300\,\mathrm{MeV}$. The number of energy and direction bins used in each simulation is listed in Table 1.

Table 1.  List of Calculations

Problem Name	Special Relativity	Spatial Resolution	Angular Resolution	MC Particles	DO $\delta {E}_{\max }$
		$r\times \theta \times \epsilon$	$\bar{\theta }\times \bar{\phi }$	($\times {10}^{9}$)
1D Calculations
1D_1x	yes	$384\times 1\times 20$	$10\times 1$	1.82	$2.40\times {10}^{-4}$
1D_2x	yes	$384\times 1\times 40$	$20\times 1$	2.33	$5.50\times {10}^{-8}$
1D_4x	yes	$384\times 1\times 80$	$40\times 1$	2.67	$7.40\times {10}^{-8}$
1D_4x_nonrel	no	$384\times 1\times 80$	$40\times 1$	⋯	$9.80\times {10}^{-8}$
1D_4x_native	yes	$384\times 1\times 80$	⋯	2.96	⋯
1D_4x_native_nonrel	no	$384\times 1\times 80$	⋯	3.25	⋯
2D Calculations
2D_LR	yes	$270\times 128\times 20$	$10\times 6$	⋯	$7.00\times {10}^{-4}$
2D_LR_nonrel	no	$270\times 128\times 20$	$10\times 6$	⋯	$5.84\times {10}^{-4}$
2D_HR	yes	$270\times 128\times 28$	14 × 10	⋯	$4.84\times {10}^{-4}$
2D_HR_native	yes	$270\times 128\times 28$	⋯	63.4	⋯

Note. We list the numerical details of each calculation presented in this paper. The left-most column gives the name of each steady-state problem that is solved by one or both of Sedonu and the NSY code. The Special Relativity column indicates whether special relativistic effects are taken into account. The Resolution column describes the spatial resolution of the fluid data and the number of neutrino energy bins. The Angular Resolution column describes the number of discrete angular bins in the neutrino momentum space angular discretization of the distribution function. In the "native" MC calculations, neutrinos are accumulated directly into angular moments (${E}_{\epsilon },{F}_{\epsilon }^{r},{F}_{\epsilon }^{\theta }$, ${P}_{\epsilon }^{{rr}},{P}_{\epsilon }^{r\theta }$, ${P}_{\epsilon }^{\theta \theta },{P}_{\epsilon }^{\phi \phi }$) without making reference to any discrete angular representation of the distribution function. The final two columns are indicative of the fidelity of the simulation. The MC Particles column shows the number of MC neutrino packets that we simulate. The DO $\delta {E}_{\max }$ measures the time variability in the steady-state solution in DO calculations. Here, $\delta {E}_{\max }$ is defined as the maximum difference from the time-average energy density during the final $1\,\mathrm{ms}$ of the simulation in $r\lt 500$ km.

Download table as:  ASCIITypeset image

The NSY code treats the advection terms in the GR Boltzmann equation semi-implicitly. Both advection and collision terms are implemented self-consistently by using a mixed-frame approach with separate momentum-space grids in the lab and comoving frames. See Appendix A and Nagakura et al. (2017a) for implementation details.

Though the NSY code is capable of evolving coupled neutrino radiation hydrodynamics, we restrict the capability of the code in this study to evolve only the radiation field on top of a fixed fluid background with a flat spacetime metric until a nearly steady-state solution is reached. As we list in Table 1, the maximum time variability in the energy density at any spatial location relative to the value averaged over $1\,\mathrm{ms}$ is less than $0.1 \%$, which is significantly smaller than the difference between the DO and MC results.

2.2. Monte Carlo

The so-called local two-moment transport method (Pomraning 1969; Anderson & Spiegel 1972; Thorne 1981; Shibata et al. 2011; Cardall et al. 2013) is the current state-of-the-art method for time-dependent multidimensional simulations of neutrino radiation hydrodynamics (e.g., Just et al. 2015b; O'Connor & Couch 2015; Foucart et al. 2016b; Kuroda et al. 2016; Roberts et al. 2016; Radice et al. 2017).⁷ In the two-moment method, Equation (1) is again reformulated, this time in terms of specific moments of the distribution function. In an orthonormal coordinate system, these moments are defined by

$$\begin{eqnarray}{E}_{\epsilon }({\boldsymbol{x}},\epsilon ) & = & \displaystyle \int d{\boldsymbol{\Omega }}\epsilon f({\boldsymbol{\Omega }},\epsilon ),\\ {F}_{\epsilon }^{i}({\boldsymbol{x}},\epsilon ) & = & \displaystyle \int d{\boldsymbol{\Omega }}\epsilon f({\boldsymbol{\Omega }},\epsilon )[{\boldsymbol{\Omega }}\cdot {{\boldsymbol{e}}}_{(i)}],\\ {P}_{\epsilon }^{{ij}}({\boldsymbol{x}},\epsilon ) & = & \displaystyle \int d{\boldsymbol{\Omega }}\epsilon f({\boldsymbol{\Omega }},\epsilon )[{\boldsymbol{\Omega }}\cdot {{\boldsymbol{e}}}_{(i)}][{\boldsymbol{\Omega }}\cdot {{\boldsymbol{e}}}_{(j)}],\\ & \cdots & \end{eqnarray} \tag{ 4 }$$

where is the neutrino energy and ${{\boldsymbol{e}}}_{(i)}$ are the basis vectors. The ellipsis denotes that there is an infinite list of moments that can be used to reconstruct the two-dimensional angular dependence of the distribution function, much like an infinite list of terms in a Taylor series can be used to reconstruct a one-dimensional function. Each of these specific moments (${M}_{\epsilon }\in \{{E}_{\epsilon },{F}_{\epsilon }^{i},{P}_{\epsilon }^{{ij}},\ldots \}$) can be integrated in energy according to

$$\begin{eqnarray}&&M=\displaystyle \frac{1}{{({hc})}^{3}}\int d\left(\displaystyle \frac{{\epsilon }^{3}}{3}\right){M}_{\epsilon }\end{eqnarray} \tag{ 5 }$$

to get the total neutrino energy density, energy flux, and so on.

In the two-moment method, only the first two moments are evolved and are assumed to provide a good enough representation of the full distribution function. The evolution equations for each moment depend on higher-order moments. In a local two-moment scheme, the pressure tensor and higher-order moments are estimated based on the energy density and flux at the same spatial location. This estimate is referred to as a closure relation, many of which have been proposed based on various motivations (e.g., Smit et al. 2000; Murchikova et al. 2017 and references therein).

In this study, we do not perform two-moment radiation transport, but rather compare the radiation pressure tensor predicted by approximate closures to the actual radiation pressure tensor output from the MC and DO calculations. We consider three popular approximate closure relations based on different physical motivations: the maximum entropy closure of Minerbo (1978) in the classical limit (Minerbo), the isotropic rest-frame closure of Levermore (1984) (Levermore), and the closure of Janka (1991) in the form presented in Just et al. (2015b) that is empirically based on MC calculations of neutrino transport in protoneutron stars. In all three cases, the pressure tensor is expressed as an interpolation between optically thick and thin limits where the solution is known analytically:

$$\begin{eqnarray}&&{P}_{\epsilon }^{{ij}}=\displaystyle \frac{3{\chi }_{\epsilon }-1}{2}{P}_{\epsilon ,\mathrm{thin}}^{{ij}}+\displaystyle \frac{3(1-{\chi }_{\epsilon })}{2}{P}_{\epsilon ,\mathrm{thick}}^{{ij}},\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}{P}_{\epsilon ,\mathrm{thin}}^{{ij}} & = & {E}_{\epsilon }\displaystyle \frac{{F}_{\epsilon }^{i}{F}_{\epsilon }^{j}}{{{\boldsymbol{F}}}_{\epsilon }\cdot {{\boldsymbol{F}}}_{\epsilon }},\\ {P}_{\epsilon ,\mathrm{thick}}^{{ij}} & = & \displaystyle \frac{{E}_{\epsilon }}{3}{I}^{{ij}}+{ \mathcal O }\,\left(\displaystyle \frac{v}{c}\right),\end{eqnarray} \tag{ 7 }$$

and I is the 3 × 3 identity matrix. In our analysis, we ignore the ${ \mathcal O }(v/c)$ term for simplicity. The different closure relations are defined by

$$\begin{eqnarray}{\chi }_{\epsilon ,\mathrm{Minerbo}} & = & \displaystyle \frac{1}{3}+\displaystyle \frac{2}{15}{\zeta }_{\epsilon }^{2}(3-{\zeta }_{\epsilon }+3{\zeta }_{\epsilon }^{2}),\\ {\chi }_{\epsilon ,\mathrm{Levermore}} & = & \displaystyle \frac{3+4{\zeta }_{\epsilon }^{2}}{5+2\sqrt{4-3{\zeta }_{\epsilon }^{2}}},\\ {\chi }_{\epsilon ,\mathrm{Janka}} & = & \displaystyle \frac{1}{3}(1+{\zeta }_{\epsilon }^{1.31}+1.5{\zeta }_{\epsilon }^{3.56}).\end{eqnarray} \tag{ 8 }$$

Here, ${\zeta }_{\epsilon }=\sqrt{{{\boldsymbol{F}}}_{\epsilon }\cdot {{\boldsymbol{F}}}_{\epsilon }/{E}_{\epsilon }^{2}}$ is referred to as the the flux factor, ${P}^{{ij}}/E$ is referred to as the Eddington tensor, and in the context of spherically symmetric radiation transport, ${P}^{{rr}}/E$ is referred to as the Eddington factor.

We perform several spherically symmetric (1D) steady-state neutrino transport calculations using different momentum-space treatments listed in Table 1. In the 1D_1x, 1D_2x, and 1D_4x simulations, the neutrino radiation field is discretized into the number of angular and energy bins described in the "Spatial Resolution" and "Angular Resolution" columns. They differ only by the number of spatial and energy bins and are all run in each the NSY code and Sedonu. The 1D_4x_nonrel calculation (performed only by the NSY code) is identical to the 1D_4x calculation, except with all velocities set to zero. In the 1D_4x_native calculation performed only by Sedonu, MC particles are accumulated directly into angular moments rather into angular bins. The 1D_4x_native_nonrel calculation is identical to the 1D_native calculation, but with all velocities set to zero. Note that length contraction causes the simulations that include relativistic effects to have a slightly larger total rest mass, but only by $2\times {10}^{-5}\,{M}_{\odot }$. Throughout this section, we primarily compare the highest-resolution DO calculation (1D_4x) to the highest-resolution native-moment MC calculation (1D_4x_native). We use the lower-resolution calculations (1D_1x and 1D_2x) in resolution comparisons.

In Figure 1, we show the static fluid background that comes from a spherically symmetric neutrino radiation hydrodynamics simulation of the collapse of a $11.2\,{M}_{\odot }$ star (Woosley et al. 2002) at $100\,\mathrm{ms}$ after core bounce using the NSY code (Nagakura et al. 2017a). For the calculations in this paper, the NSY code is again used to calculate a steady-state solution of the neutrino radiation field on this background using the DO method. The opacities and emissivities are then exported to Sedonu, which computes a steady-state radiation field on the same background.

Zoom In Zoom Out Reset image size

In Figure 2 we show radial profiles of the total energy density of each neutrino species using both the 1D_4x DO calculation and the 1D_4x_native MC calculation. For each species, the energy density is approximately constant in the inner core ($r\lesssim 10\,\mathrm{km}$) as neutrinos are trapped and in equilibrium with the fluid. In this region, the neutrinos have a Fermi–Dirac distribution function, so the energy density is determined entirely by the fluid temperature and the electron and nucleon chemical potentials. In the outer (i.e., shock-processed) core at $r\gtrsim 10\,\mathrm{km}$, the temperature is very high ($T\approx 20\,\mathrm{MeV}$), and many more electron antineutrinos and heavy lepton neutrinos are emitted than in the cooler inner core. Beyond the energy-averaged neutrinospheres ($30\,\mathrm{km}\,\lesssim r\lesssim 60\,\mathrm{km}$, depending on the neutrino species), the neutrinos are only weakly coupled to the fluid, and the energy density decreases as $E\propto {r}^{-2}$.

Zoom In Zoom Out Reset image size

Both codes produce remarkably similar results, with differences in the energy density (Figure 2) smaller than 2% everywhere within the shock. The remaining differences near the energy-averaged neutrinospheres ($30\,\mathrm{km}\lesssim r\lesssim 60\,\mathrm{km}$, depending on species) are due to the MC random walk approximation and decrease when the critical optical depth is increased (see Appendix B). The differences outside the decoupling region ($r\gtrsim 40\mbox{--}80\,\mathrm{km}$, depending on species) decrease with increasing DO directional angular resolution. Outside of the shock, the difference in the energy density grows as the NSY code experiences a small departure from the ${r}^{-2}$ scaling of the energy density. This is an artifact of the finite spatial resolution. The size of the error visibly increases at $500\,\mathrm{km}$, where the radial resolution coarsens.

Figure 3 shows a quantity akin to the comoving-frame average neutrino energy,⁸ defined as

$$\begin{eqnarray}&&\bar{\epsilon }=\displaystyle \frac{\sum {E}_{i}}{\sum {E}_{i}/{\epsilon }_{i}},\end{eqnarray} \tag{ 9 }$$

where E_i is the lab-frame energy density in bin i and ${\epsilon }_{i}$ is the comoving-frame bin central energy. Because neutrinos are in equilibrium with the fluid below the decoupling region ($r=30\mbox{--}70\,\mathrm{km}$, depending on species), they have a Fermi–Dirac distribution function that depends only on the fluid temperature and electron and nucleon chemical potentials. The neutrino absorption opacity scales approximately as ${\kappa }_{a}\sim {\epsilon }^{2}$, so high-energy neutrinos are preferentially absorbed, causing the average neutrino energy to continuously decrease with radius. The average energy jumps at the shock front, since after passing the shock front, neutrinos are moving through matter falling with speeds of $| v| \sim 0.1c$. The comoving-frame neutrino energy is thus Doppler boosted even though the lab-frame energy density is constant across the shock.

Zoom In Zoom Out Reset image size

The differences in the average neutrino energy between the 1D_4x DO calculation and the 1D_4x_native MC calculation are smaller than about $0.1\,\mathrm{MeV}$ within the shock. Analyzing the various potential sources of errors, the differences at $r\lesssim 30\,\mathrm{km}$ are simply statistical noise that decreases with the square root of the number of MC neutrino packets simulated. The differences below the shock are primarily due to the MC random walk approximation error and decrease with an increased critical optical depth (see Appendix B) independent of the momentum-space resolution. The differences near and outside the shock are a result of finite energy resolution, which results in an interpolation error when the NSY code transforms energy and direction bins between the comoving and lab frames.

We plot the direction-integrated neutrino energy density spectra at $r=105\,\mathrm{km}$ for each neutrino species in Figure 4. This point is below the shock in the semitransparent region. The results of the 1D_4x DO calculation and the 1D_4x_native MC calculations agree in every bin with at most $\approx 1.5 \%$ of the peak value. In energy bins with little energy density, the relative errors become quite large, but bins with such small energy density have much less dynamical effect in CCSN simulations. Some statistical noise from the MC calculation can be seen in the range of $5\mbox{--}20\,\mathrm{MeV}$, but the small overall offsets are due to the finite neutrino energy and angular resolution.

Zoom In Zoom Out Reset image size

Figure 5 shows the lab-frame energy-integrated heavy lepton neutrino distribution function at three separate radii. The red curves are from the radius where $\rho =2\times {10}^{12}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and show that the neutrinos are nearly isotropic. The blue curves are from near the shock front and are nearly free streaming and very forward-peaked, as almost no neutrinos are moving inward. The green lines show the distribution function at an intermediate location of $r=68\,\mathrm{km}$ ($\rho =4.7\times {10}^{10}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$) between the trapped and free-streaming limits. In the plot, the distribution function is normalized by the largest value so the shapes can be easily compared. We assume a constant value for the distribution within the directional bin in the forward direction and linearly interpolate the distribution function for all other directions. This is done to ensure that in postprocessing the value of the distribution function never exceeds one. However, this gives rise to the artificially flattened nose of the distribution functions most apparent in the blue curves.

Zoom In Zoom Out Reset image size

The thickest lines in Figure 5 are from high-resolution 1D_4x MC and DO simulations, while thinner lines indicate lower resolution. The 1D_4x_native MC simulation does not collect data on a grid of discrete angular bins, so its results cannot be used to make such a plot. The importance of the angular resolution is very apparent for the blue curves at the shock front, since most of the neutrino energy is in a single angular bin. The MC results look remarkably similar to the DO results, though a lack of numerical diffusion in the MC calculations allows for slightly more sharply forward-peaked distribution functions for a given angular resolution. This angular dependence is reflected in all angular moments of the distribution function.

In Figure 6, we show the energy-integrated lab-frame radial flux factor (${F}^{r}/E$, see Equation (4)) of the distribution function of all three neutrino species using the 1D_4x DO calculation and the 1D_4x_native MC calculation. Below around $30\,\mathrm{km}$, the neutrinos are trapped and the distribution function is nearly isotropic, resulting in a minuscule flux relative to the energy density (corresponding to the red curves in Figure 5). In the transition region (corresponding to the green curves in Figure 5), an increasing fraction of the neutrino radiation energy is moving radially outward, causing the flux factor to approach one at large radii (corresponding to the blue curves in Figure 5). Here, ${F}^{\theta }$ and ${F}^{\phi }$ are identically zero due to spherical symmetry.

Zoom In Zoom Out Reset image size

The angular moments of the radiation field are naturally sensitive to the angular grid resolution. We see small differences of at most around 0.02 in the flux factor, but the size of this difference scales approximately linearly with the angular grid size for calculations with a coarser angular grid (not plotted). Sedonu consistently predicts a more rapid transition to free streaming than does the NSY code. Here the MC method shows a significant advantage in that by computing moments directly rather than postprocessing from an angular grid, we get angular moments with effectively infinite angular resolution. The NSY code comes very close to this solution but suffers from some angular diffusion. This causes the NSY code to predict distribution functions that are slightly, though artificially, more isotropic. The difference approaches a small but constant value at large radii, where almost all of the energy in the DO calculations is contained in the single outward-pointing angular bin. The Sedonu results are, however, visibly noisy in the difference plot, since subtractive cancellation tends to amplify statistical noise. There is a small hump visible in the heavy lepton neutrino difference plot just below $r=100\,\mathrm{km}$ that is a result of the MC random walk approximation. The size of this hump decreases when the critical optical depth is increased (see Appendix B), bringing it closer to the electron neutrino and electron antineutrino difference curves.

In Figure 7, we show the rr component of the energy-integrated lab-frame Eddington tensor (${P}^{{ij}}/E$, see Equation (4)) of the distribution function of all three neutrino species. Only the diagonal components of the Eddington tensor (${P}^{{rr}}/E,{P}^{\theta \theta }/E$, and ${P}^{\phi \phi }/E$) are nonzero in spherical symmetry. At $r\lesssim 30\,\mathrm{km}$, all diagonal components of the Eddington tensor are one-third because the radiation is nearly isotropic. After decoupling, the rr component approaches unity as all radiation is moving radially outward, while the $\theta \theta$ and $\phi \phi$ components (not plotted) approach zero.

Zoom In Zoom Out Reset image size

Once again, the differences between Sedonu and the NSY code scale approximately linearly with the neutrino direction angular zone sizes. However, the maximum difference of 0.03 is larger than the maximum flux factor difference of 0.02. Unlike with previously discussed radiation quantities, the random walk approximation does not add significant error to P^rr. Though we do not plot ${P}^{\theta \theta }$ or ${P}^{\phi \phi }$, the differences between the MC and DO results are similar to those of P^rr. Since the integral in Equation (4) contains a factor of ${({\boldsymbol{\Omega }}\cdot \hat{r})}^{2}$, the results do not suffer from subtractive cancellation, and the amount of statistical noise is significantly lower than that of the flux factor.

The dependence of the angular moments on angular resolution can be clearly seen in Figure 8, where we plot the rr component of the Eddington tensor for four MC calculations (solid lines) and three DO calculations (dashed lines). We first direct our attention to the blue lines, corresponding to the low-resolution 1D_1x DO and MC calculations. Even though both are postprocessed in the same way from the same angular grid, the MC results transition to large values of P^rr faster than do the DO results. At $r\approx 300\,\mathrm{km},{P}^{{rr}}/E$ saturates at the maximum value possible given the angular resolution, which the DO results approach at large radii. The same is true for the higher-resolution green and red curves, but the saturation occurs at a larger radius and is not so visibly obvious.

Zoom In Zoom Out Reset image size

Due to the effectively infinite angular resolution of the 1D_4x_native MC calculation, the corresponding black line in Figure 8 can be thought of as exact, modulo a small amount of MC noise. Going from coarsest to finest resolution, the maximum differences between the DO results and the black curve are 0.125, 0.057, and 0.028. This corresponds to a factor of 2.2 improvement when going from 1× to 2× resolution, and a factor of 2.0 when going from 2× to 4× resolution. Similarly, the maximum differences between the gridded MC and the native MC results are, in order of increasing resolution, 0.0896, 0.0243, and 0.0062. The accuracy improves by a factor of 3.7 when going from 1× to 2× resolution, and by a factor of 3.9 when going from 2× to 4× resolution. This trend, where DO results are near first-order convergence and gridded MC results are near second-order convergence, is apparent in the flux factor results as well. This is because the postprocessing angular integration scheme is second order (except in the forward-most bin, where it is first order), but the evolution scheme in the NSY code is only first order.

The use of an approximate analytic closure in two-moment radiation transport schemes is significantly faster than either the DO or MC methods. However, since there are many reasonable closures available, it is of great interest to evaluate how well these closures perform against our full Boltzmann results. We replot the electron neutrino P^rr curve (black line in Figure 8) from the 1D_4x_native MC calculation as a solid black line in Figure 9. We then use E and F^r from the same MC calculation to estimate P^rr using the three analytic closures given in Equation (8). The Janka and Minerbo closures perform similarly and have a maximum difference with MC of ∼0.03, which is comparable to the differences between the 1D_4x DO calculation and the 1D_4x_native MC calculation. The Levermore closure, however, performs better, with a maximum difference of ∼0.006. This is significantly better than the accuracy of any DO result and is comparable to the accuracy of the 1D_4x MC calculation. These results are also replicated in a similar analysis of ${P}^{\theta \theta }$ and ${P}^{\phi \phi }$ (not plotted). In short, analytic closures perform remarkably well in this particular steady-state spherically symmetric transport calculation.

Zoom In Zoom Out Reset image size

The primary role of neutrinos in the explosion mechanism of CCSNe is redistributing thermal energy from the protoneutron star region to the gain region that drives turbulence and pushes the shock outward. The relevance of these detailed transport calculations comes down to how the differences between the methods affect the heating and cooling of matter in the supernova. In Figure 10, we show the comoving-frame net gain, i.e., the heating rate less the cooling rate due to neutrinos. We show results from the 1D_4x MC and DO calculations (red line), the 1D_4x_nonrel MC and DO calculations (green line), and the 1D_1x MC and DO calculations (blue line). Below about $90\,\mathrm{km}$ the fluid is overall cooling, and the emitted neutrinos pass through the gain region from $90\,\mathrm{km}$ to $170\,\mathrm{km}$ and deposit energy. Below the shock, the net gains from the 1D_4x MC and DO calculations are very similar, with differences of $\lesssim 1 \%$ of the peaks in the gain curve.

Zoom In Zoom Out Reset image size

Just outside the shock, the fluid densities are low, and most nucleons are bound in nuclei. Because of this, the heating is due primarily to neutrino pair annihilation. The pair annihilation rates are underresolved in neutrino energy space even with 80 energy bins, resulting in significant differences between heating rates of different resolutions. However, the test results show only a 20% difference between the 1D_4x results and a test with an energy-space resolution of 160 bins. We can use the radial profiles of heating rate, density, and velocity to estimate the amount of energy per nucleon the fluid is heated before passing through the shock as

$$\begin{eqnarray}&&{\rm{\Delta }}E\approx \int \displaystyle \frac{{m}_{n}{ \mathcal H }(r)}{\rho (r)| v(r)| }{dr},\end{eqnarray} \tag{ 10 }$$

where ${ \mathcal H }$ is the heating rate. Using the heating rate from the highest-resolution simulations, this predicts a total heating of $\sim 0.1\,\mathrm{MeV}\,{\mathrm{nucleon}}^{-1}$. Compared to the postshock temperatures of $T\lesssim 20\,\mathrm{MeV}$ and the iron-56 binding energy of $8.8\,\mathrm{MeV}$, this preshock heating is insignificant.

The differences between MC and DO are amplified outside the shock, where we must divide a volumetric heating rate (in erg cm⁻³ s⁻¹) by the low density to get a specific heating rate. Also, recall that in our calculation of pair annihilation rates, we assume that the neutrino distribution functions are isotropic (see Appendix D for details). At large radii relative to the neutrinospheres and to leading order, however, the angular dependence actual reaction rates are (e.g., Bruenn 1985)

$$\begin{eqnarray}&&{R}_{\mathrm{pair},\mathrm{abs}}(\epsilon ,\bar{\epsilon },{\boldsymbol{\Omega }}\cdot \bar{{\boldsymbol{\Omega }}})\propto 1-{\boldsymbol{\Omega }}\cdot \bar{{\boldsymbol{\Omega }}}\sim {\left(\displaystyle \frac{r}{{r}_{\nu }}\right)}^{-2},\end{eqnarray} \tag{ 11 }$$

where r_ν is the neutrinosphere radius. The location of the neutrinosphere depends on the neutrino species and energy, but for a typical radius of ${r}_{\nu }=50\,\mathrm{km}$, this angular term scales the reaction rate by a factor of ∼0.1 at the shock. Thus, we expect the heating rates (and hence the heating rate differences) to be overestimated by a factor of ∼10 at the shock.

Including velocity dependence in neutrino transport algorithms is a complication that can significantly increase the complexity and cost of the transport calculation. It is natural to attempt to quantify the size of the error made in codes that neglect velocity dependence. We repeat the high-resolution calculations with the same rest-frame fluid profile shown in Figure 1, but set all velocities to zero. Velocity dependence changes the comoving frame neutrino energies and directions, modifying the rates at which neutrinos interact with the fluid. This has a very minor effect below the shock, but significantly changes the heating rate outside the shock where velocities are $\sim 0.1c$. However, the density drops by a factor of 10 across the shock, and the preshock fluid moves so quickly that the overall heating is negligible. These small errors outside the shock are unlikely to have a significant impact on simulation results. The volume-integrated net gain in the gain region (where there is net heating under the shock) is $2.16\times {10}^{51}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ in the 1D_4x DO calculation and $2.18\times {10}^{51}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ in the 1D_4x_native MC calculation, a difference of only $0.34 \%$. Compare this to the difference of the same quantity between the 1D_4x and 1D_1x DO calculations of $2.0 \%$ and between the 1D_4x and 1D_4x_nonrel DO calculations of $1.2 \%$. Though including velocity dependence impacts the heating rates more than the differences between the codes in the highest-resolution case, low resolution can cause significantly larger inaccuracies.

In this section, we describe results from the first multidimensional comparison of Boltzmann-level neutrino transport codes. We present a set of four axisymmetric time-independent neutrino transport calculations as listed in Table 1. Once again, the NSY code is used to calculate an approximately steady-state solution, and the opacities and emissivities are exported from the completed NSY calculations to Sedonu for the MC calculation. Due to computational cost, we only consider two resolutions in the DO code. The low-resolution (2D_LR) calculations have momentum-space resolution equivalent to the 1D_1x calculations, and the high-resolution (2D_HR) calculations have momentum-space resolution between that of the 1D_1x and 1D_2x calculations.

The rest-frame fluid profile used in all simulations shown in Figure 11 comes from a 2D simulation of the collapse of the same $11.2\,{M}_{\odot }$ star (Woosley et al. 2002) used in Section 4 (Nagakura et al. 2017a). In Figure 11 and in all other color-map plots in this section, data separated into quadrants show data from the northern hemisphere of the calculation only to ease visual comparison of data sets. Data on half-circles show the full simulation domain out to $r=210\,\mathrm{km}$. In Figure 11, a multidimensional structure in all fluid quantities is apparent and is due to neutrino-driven turbulent convection. The postshock velocity field in particular shows fluid speeds up to $0.037c$, compared to the maximum radial velocity of $0.015c$ in the 1D calculations. This multidimensional structure provides a challenge for any radiation transport method.

Zoom In Zoom Out Reset image size

We begin with a comparison of the spectral properties of the results from Sedonu and the NSY code. Figure 12 shows the comoving-frame average energy of each of the three simulated neutrino species for the 2D_HR DO calculation and the 2D_HR_native MC calculation. Just as in the 1D calculations, the electron neutrinos have the highest energy in the inner core due to the high electron chemical potential, and the lowest energy at large radii since they decouple at the largest radius and the lowest matter temperature. The DO and MC results are nearly identical, so differences can only be seen in the right half of each plot, where we subtract the DO results from the MC results. The average energies differ between the MC and DO results by at most 0.5 MeV, which is larger than in the 1D results in Figure 3 due to the lower energy resolution. When electron neutrinos and antineutrinos decouple from matter, the opacity is dominated by absorption. Because of this, MC packets that use the random walk approximation quickly lose energy, preventing errors from the random walk approximation from propagating outward through the rest of the domain. However, the heavy lepton neutrino opacity is dominated by scattering, so MC packets carrying errors from the random walk approximation retain their energy when traveling through the rest of the domain, causing errors to be slightly higher. Just as in the 1D results, the differences between the MC and DO average energy errors jump across the shock, since the NSY code suffers from numerical diffusion when transforming between grids in the two-grid approach (Appendix A). There are also a number of hot spots in the average energy differences within the shock, which correspond to regions of high velocity in Figure 11. These differences are also due to some numerical diffusion in the two-grid approach. Heavy lepton neutrinos are more strongly impacted by the protoneutron star convection, since they decouple at a smaller radius. The features visible in the heavy lepton neutrino average energy difference plot (the right panel of Figure 12) at small radii are diminished by reducing the MC random walk critical optical depth, independent of momentum-space resolution.

Zoom In Zoom Out Reset image size

Though the energy resolution is coarser than in the 1D calculations, we are able to compare the full spectra at a given location. For Figure 13, we choose the same radius of $105\,\mathrm{km}$ used for Figure 4 and an angle of 36° from the north pole. We plot the direction-integrated spectra of all three neutrino species using the 2D_HR DO calculation and the 2D_HR_native MC calculation. The neutrino energy density within each comoving-frame energy bin is measured in the lab frame and the individual neutrino energy in the comoving frame, resulting in a mixed-frame quantity. The results are remarkably similar and effectively reproduce the 1D results. The heights of the peaks differ by $\sim 5 \%$, which is comparable to the differences between MC and DO results in the energy density in the lower-resolution 1D results.

Zoom In Zoom Out Reset image size

In Figure 14, we plot the energy-integrated lab-frame electron neutrino flux factor and Eddington tensor using the 2D_HR DO calculation and the 2D_HR_native MC calculation. These plots effectively reproduce the 1D angular moment results in Figures 6 and 7. We again see that the MC results exhibit a quicker transition to a forward-peaked distribution, that the errors in the second moment (${P}^{{rr}}/E$) are larger than in the first (${F}^{r}/E$), and that the MC noise in the second moment is smaller than in the first. The magnitude of the differences in ${P}^{{rr}}/E$ of ∼0.1 at the shock can be compared to the 1D results in Figure 8. The 2D_HR differences are between the 1D_1x and 1D_2x differences, reflecting the fact that the 2D_HR angular resolution is between that of the 1D_1x and 1D_2x calculations. We also demonstrate this resolution dependence in the left (for ${F}^{r}/E$) and center (for ${P}^{{rr}}/E$) plots of Figure 15. The top left quadrant of each shows the difference between the moments calculated using the 2D_HR DO calculation and the 2D_HR_native MC calculation. The bottom left quadrant shows the difference between the 2D_LR DO calculation and the same MC calculation. The differences are significantly smaller for the higher-resolution DO calculation, indicating that the DO results are converging to the MC results with increasing resolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Unlike in spherical symmetry, ${P}^{r\theta }$ is not identically zero in axisymmetry and is thus a sensitive probe of multidimensional effects on the radiation field. Note that ${P}^{r\phi }$ and ${P}^{\theta \phi }$ are still identically zero, since we do not consider azimuthal fluid velocities. In Figure 16, we plot the energy-integrated lab-frame ${P}^{r\theta }/E$ for all three neutrino species using the 2D_HR DO calculation and the 2D_HR_native MC calculation. Since the dominant neutrino propagation direction is radial, the off-diagonal components of the pressure tensor are strongly correlated with the corresponding flux. In this particular snapshot, ${P}^{r\theta }$ happens to be overall mostly positive, and we find the morphology to be indeed very similar to ${F}^{\theta }$ (not shown). Generally, both positive and negative values are to be expected (see Nagakura et al. 2017a). It is interesting to note that the electron neutrino and electron antineutrino plots have complementary hot spots. Within the protoneutron star, nonradial neutrino fluxes are present due to turbulent fluid carrying trapped neutrinos. Outside of the convective zone of the protoneutron star (PNS) but still within the neutrinospheres, electron neutrinos tend to diffuse away from regions of high electron chemical potential while electron antineutrinos diffuse toward them. In tests where the inner $105\,\mathrm{km}$ is excluded from the calculation, the ${P}^{r\theta }$ distribution is much more uniform, suggesting that the hot spots are due to a combination of anisotropic neutrinos from the neutrinosphere interacting with multidimensional features in the fluid background.

Zoom In Zoom Out Reset image size

Once again, the MC and DO results for ${P}^{r\theta }$ look remarkably similar. Unlike for the diagonal moments, much subtractive cancellation occurs when computing ${P}^{r\theta }$, which in turn requires a large number of MC particles to drive down the noise. Similar to the other moments, the MC calculation tends to show larger values of ${P}^{r\theta }$, since its effectively infinite angular resolution is able to resolve finer angular structures. We demonstrate this resolution dependence in the right-most plot of Figure 15. The top left quadrant shows the difference between the electron neutrino ${P}^{r\theta }$ from the 2D_HR DO and 2D_HR_native MC calculations, while the bottom left quadrant compares the 2D_LR DO calculation to the same MC calculation. The differences are significantly larger for the lower-resolution calculation, indicating that the DO calculation is converging to the MC result with increasing angular resolution.

Figure 17 compares components of the electron neutrino Eddington tensor computed by the 2D_HR DO calculation to those predicted by the Levermore closure using E and Fⁱ from the same DO calculation. We demonstrated in Section 4 that in our spherically symmetric snapshot, the Levermore closure predicts ${P}^{{rr}}/E$ and ${P}^{\theta \theta }/E$ from only the flux factor with an accuracy within 0.01 of the actual Eddington tensor value. This result is reproduced in two dimensions for electron neutrinos in Figure 17. The left and center plots show ${P}^{{rr}}/E$ and ${P}^{\theta \theta }/E$, respectively. The top left quadrant of each shows the moment computed directly from the 2D_HR DO calculation (same data as depicted in Figure 14), and the bottom left quadrant shows the same moment predicted by the Levermore closure. They are visually identical, and the error plotted on the right side of each plot shows a maximum error of 0.014 in ${P}^{{rr}}/E$ and a maximum error of 0.0089 in ${P}^{\theta \theta }/E$. Though there is some multidimensional structure in how accurately the Levermore closure predicts the diagonal components, this effectively mirrors the results of the 1D calculations. The right-most plot shows ${P}^{r\theta }/E$ (same data as in Figure 16), multiplied by 10 for visibility on this color scale. The Levermore closure predicts this component of the moment within 0.0077. This is large compared to this component's maximum value of 0.012 and compared to a difference of ∼0.003 between DO and MC results. Thus, though this analytic closure has small relative errors for the diagonal components of the Eddington tensor, it has difficulty accurately predicting the small off-diagonal components in this CCSN snapshot. The Minerbo and Janka closures show errors at smaller radii, but the extrema of these errors are only slightly larger than those of the Levermore closure. Other neutrino species behave very similarly, except that the heavy lepton neutrino values for ${P}^{r\theta }$ (and hence errors) are significantly smaller.

Zoom In Zoom Out Reset image size

Ignoring special relativistic effects in radiation transport calculations greatly simplifies the problem. Fluid velocities are also generally only a few percent of the speed of light below the shock, but are larger in 2D ($\sim 0.037c$) than in 1D ($\leqslant 0.015c$). We test the effects of ignoring fluid velocities in Figure 15, where we plot the difference between the 2D_LR_nonrel and 2D_LR DO calculations. The error in ${F}^{r}/E$ and ${P}^{{rr}}/E$ from ignoring velocities is much smaller than the difference between MC and DO calculations or the difference between resolutions. The only exception is in the convective region of the protoneutron star, where the flux is determined entirely by the fluid velocity because the neutrinos are trapped. The magnitude of this error is at most comparable to the error induced by the coarse resolution, and is significantly smaller than the error in all components of the second moment predicted by the Levermore closure (Figure 17).

Finally, in Figure 18, we investigate how these different transport schemes affect the actual heating and cooling rates of the fluid. Once again, we show the results from the 2D_HR DO calculation in the top left quadrant, which outside the core appears visually identical to the 2D_HR_native MC calculation results in the bottom left quadrant. Just as in Figure 10, there is significant statistical noise within the core, where neutrinos are largely in equilibrium with the matter. We depict the difference between these results on the right half of the plot. The MC results show more rapid cooling in the cooling region and more rapid heating in the outer regions of the gain layer, but only by a few percent of the maximum net gain outside the protoneutron star. This is similar to the behavior of the lower-resolution 1D results in Figure 10, where MC calculations predict a smaller gain than do the DO calculations at $r\lesssim 125\,\mathrm{km}$ for the 1D_2x calculation and at $r\lesssim 140\,\mathrm{km}$ for the 1D_1x calculations. The 2D MC calculation also predicts larger heating than does the 2D_HR DO calculation in regions of high inward velocity. This is again an effect of the limited momentum-space resolution in the DO calculations that makes the two-grid approach somewhat diffusive in angle and energy. This is to be expected given that the average neutrino energies in these regions (Figure 12) are higher in the MC calculations. Overall, excluding the noisy region in the core, these errors are at most $\sim 2 \%$ of the amplitude of the net gain curve in Figure 10.

Zoom In Zoom Out Reset image size

The volume-integrated gain of the gain region (where there is net heating under the shock) is $9.00\times {10}^{51}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ in the 2D_HR DO calculation and $8.93\times {10}^{51}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ in the 2D_HR MC calculation, which is only a $0.35 \%$ difference. Compare this to the relative error of the same quantity between the 2D_LR and 2D_HR DO calculations of $1.3 \%$, which is smaller than in the 1D resolution comparison because our 2D resolutions are much more similar. Even so, the errors from low resolution are significantly more significant than the differences between the methods. The difference in integrated heating rate between the 2D_LR and 2D_LR_nonrel DO calculations is $2.0 \%$. This is larger than in the 1D relativity comparison because fluid velocities under the shock are larger in the 2D calculations than in the 1D calculations due to convective motion. Note that the integrated heating rate should not be compared with the 1D results because the fluid profiles are significantly different.

Neutrinos dominate energy transport in CCSNe and are a vital component of the CCSN explosion mechanism. It is therefore imperative to ensure neutrinos are simulated accurately in CCSN models. One means of checking the accuracy of a computational method is comparing against another accurate method. The grid-based DO method and the particle-based MC method both solve the full transport problem in very different ways. We perform the first detailed multidimensional comparison of special relativistic Boltzmann-level neutrino transport codes in the context of core-collapse supernovae using the grid-based DO code of Nagakura et al. (2017a) and the particle-based MC Sedonu code. We verify that both methods converge to the same result in the limit of large MC particle count and fine DO momentum-space resolution under the assumption of a static fluid background in spherical symmetry and in axisymmetry. This provides confidence in the accuracy of the results from these two completely different approaches.

We demonstrate an agreement of the average neutrino energy to within ∼0.1 MeV for 1D calculations and ∼0.5 MeV for coarser 2D calculations everywhere in the simulation domain for all three simulated neutrino species. We also demonstrate exquisite agreement in the local spectra of all three species. We find that numerical diffusion from a coarse momentum-space resolution dominates these small errors, though smaller errors result from finite spatial resolution and from the MC random walk approximation.

MC transport computes angular moments of the distribution function with great accuracy when the moments are computed natively during the calculation as opposed to postprocessed from an angular grid. The DO method requires a very high angular resolution of about 40 points in the polar direction to compute these moments with similar accuracy in 1D calculations, which is currently not possible in 2D calculations and certainly not possible in 2D time-dependent simulations. The MC method, however, requires a large number of particles to be simulated to reach low noise levels in moments that exhibit subtractive cancellation (i.e., Fⁱ in optically deep regions, ${P}^{r\theta }$ in 2D calculations). The present 2D calculation simulated 63 billion particles and still shows some noise in these quantities.

The approximate two-moment radiation transport scheme is significantly more efficient than either DO or MC transport by evolving on the energy density and flux. However, this method requires an ad hoc closure relation to complete the system of equations by making estimates of higher-order moments. We evaluate the performance of the Levermore, Janka, and Minerbo closures in predicting the second angular moments from the first. In the 1D calculations, the Levermore closure performs best, with an error comparable to the differences between the highest-resolution DO results and the MC results. In 2D calculations, this closure performs comparably well when predicting diagonal components of the second moment, but this accuracy is not sufficient for determining the very small off-diagonal components. Though careful tests would be required to assess the importance of these small off-diagonal components, these components reflect the multidimensional nature of CCSN dynamics.

Finally, we find that the difference in local heating and cooling rates between the DO and MC methods is at most 2% of the amplitude of the net gain curve in the cooling region of the CCSN in both 1D and 2D calculations. The volume-integrated gain in the gain region (where there is net heating under the shock), however, differs by only $0.3 \%$ in the highest-resolution 1D calculations and by $0.4 \%$ in the highest-resolution 2D calculations. In these cases, the DO and MC schemes share the same energy resolution, but the MC scheme has effectively infinite angular resolution. The differences in the same quantities due to changing the DO energy (and angular) resolution are larger than 1% in both 1D and 2D calculations, indicating that neutrino energy resolution is the dominant source of real error. Since both the MC and DO methods rely on opacities and emissivities discretized into energy bins, both suffer from this error. The errors in radiation quantities (energy density, angular moments, and average neutrino energies) below the shock are dominated by the finite momentum-space resolution of the DO calculations and the statistical noise and limited energy resolution of the MC calculations.

Though it is important to simulate all physics relevant to the CCSN mechanism, the numerical resolution can pose a significantly larger threat to simulation accuracy than a lack of physical elements. We test the effects of ignoring the special relativistic Lorentz transformation of neutrinos and find it to be severely subdominant to errors induced by low momentum-space resolution at the resolutions we use. The diagonal components of the Eddington tensor in the low-resolution DO calculations show resolution-induced errors of $\sim 15 \%$, and even the highest-resolution 1D calculations (which would be impossible in 2D) show errors of $\sim 3 \%$. This underscores the need for resolution tests in interpreting results of simulation results.

Though this study inspires much more confidence in both methods, we must mention several caveats. The largest is that opacities and emissivities are an extremely important component of radiation transport. In order to facilitate a detailed comparison of the methods themselves, we carefully configure both codes to use identical opacities at each spatial location and neutrino energy bin. However, we do not compare the effects of different approximations and physical processes present in the opacities that may overwhelm the small differences in the results between these codes.

Second, we must emphasize that our calculations are performed under the assumption of an unchanging fluid background and flat metric at one particular stage in the CCSN evolution. At different stages, especially during early postbounce prompt convection and the shock revival phase, the matter distribution and hence neutrino radiation fields are significantly different and would benefit from a similar analysis. We also use an approximate treatment of pair processes and neglect electron scattering. These simplifications are made to bring both codes to the same level, where we could be sure that they are simulating the same physics with the same level of approximation. This allows an isolated evaluation of the relative performance of both transport methods, but neglects many components of physics that should be included in realistic dynamical CCSN simulations.

The impact of the time-dependent features of the radiation field on the fluid dynamics will be the next necessary step in verifying neutrino radiation hydrodynamics codes. A similar careful verification of the choice and implementation of different neutrino interactions, the resolution and discretization scheme (including mesh geometry and refinement), the treatment of gravity, and the numerical hydrodynamics scheme would also greatly benefit the interpretation of simulation results. We leave this broader task of evaluating multidimensional time-dependent radiation hydrodynamics simulations of CCSNe to future work.

We release the data input to and output by both codes at http://www.stellarcollapse.org/MCvsDO. The Sedonu code is also open source and available at https://bitbucket.org/srichers/sedonu.git, along with a set of ready-to-run input data and parameter files to run the calculations in this paper. This Sedonu release contains many performance, usability, and flexibility changes implemented since previous releases. In addition, we incorporate a special relativistic, time-independent version of the MC random walk approximation that enables Sedonu to efficiently calculate neutrino transport through regions of arbitrarily large optical depth.

We would like to acknowledge Ryan Wollaeger, Kendra Keady, Adam Burrows, David Radice, Luke Roberts, and Yuki Nishino for many insightful discussions of radiation transport methods. S.R. was supported by the National Science Foundation (NSF) Blue Waters Graduate Fellowship. This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (Grants No. OCI-0725070 and No. ACI-1238993) and the State of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. This work benefited from access to the NSF XSEDE network under allocation TG-PHY100033 and to Blue Waters under PRAC award no. ACI-1440083. H.N. is supported in part by JSPS Postdoctoral Fellowships for Research Abroad No. 27-348. C.D.O. is supported in part by the International Research Unit of Advanced Future Studies at the Yukawa Institute for Theoretical Physics. This work is furthermore partially supported by the Sherman Fairchild Foundation and NSF under award nos. TCAN AST-1333520, CAREER PHY-1151197, and PHY-1404569. This work is supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan (15K05093, 16H03986, 26104006) and the HPCI Strategic Program of Japanese MEXT and K computer at the RIKEN and Post-K project (Project ID: hpci 170304, 170230, 170031). J.D. acknowledges support from the Laboratory Directed Research and Development program at Los Alamos National Laboratory (LANL). Work at LANL by S.R. and J.D. was done under the auspices of the National Nuclear Security Administration of the US Department of Energy. This paper has been assigned Yukawa Institute for Theoretical Physics report number YITP-17-61.

The NSY code solves the conservative form of the general relativistic Boltzmann equation, which can be written as (Nagakura et al. 2017a)

$$\begin{eqnarray}&&{\left.\displaystyle \frac{1}{\sqrt{-g}}\displaystyle \frac{\partial }{\partial {x}^{\alpha }}\right|}_{{\boldsymbol{q}}}\left[\sqrt{-g}f\left({e}_{(0)}^{\alpha }+\displaystyle \sum _{i=1}^{3}{{\ell }}_{(i)}{e}_{(i)}^{\alpha }\right)\right]\\ &&\quad -\,\displaystyle \frac{1}{{\epsilon }^{2}}\displaystyle \frac{\partial }{\partial \epsilon }({\epsilon }^{3}f{\omega }_{(0)})+\displaystyle \frac{1}{\sin \bar{\theta }}\displaystyle \frac{\partial }{\partial \bar{\theta }}(\sin \bar{\theta }\,f{\omega }_{(\bar{\theta })})\\ &&\quad +\,\displaystyle \frac{1}{{\sin }^{2}\bar{\theta }}\displaystyle \frac{\partial }{\partial \bar{\phi }}(f{\omega }_{(\bar{\phi })})={S}_{\mathrm{rad}},\end{eqnarray} \tag{ 12 }$$

where ${S}_{\mathrm{rad}}$ describes the collision term for neutrino–matter interactions, g is the determinant of the four-dimensional metric, and ${x}^{\alpha }$ are the spacetime coordinates; ${e}_{(\mu )}^{\alpha }$ ($\mu =0,1,2,3$) denote a set of tetrad bases for a local orthonormal frame on four-dimensional manifolds. In the present study, we assume that the spacetime is flat, so we simply use spherical-polar coordinates: ${x}^{0}=t,{x}^{1}=r,{x}^{2}=\theta$, ${x}^{3}=\phi ,{e}_{(0)}=\hat{t}$, ${e}_{(1)}=\hat{r},{e}_{(2)}=\hat{\theta },{e}_{(3)}=\hat{\phi }$, and $g=-{r}^{2}\sin \theta$. Here, ${x}^{1},{x}^{2}$, and x³ contain the same information as ${\boldsymbol{x}}$ in Section 2. The neutrino momentum space is also written in spherical-polar coordinates qⁱ ($i=1,2,3$). Here, ${q}^{1}=\epsilon \equiv -{p}_{\alpha }{e}_{(0)}^{\alpha }$ is the energy of a neutrino with four-momentum ${p}^{\alpha }$, ${q}^{2}=\bar{\theta }$ is the polar direction angle with respect to $\hat{r}$, and ${q}^{3}=\bar{\phi }$ is the azimuthal angle with respect to $\hat{\theta }$. Note that q² and q³ contain the same information as ${\boldsymbol{\Omega }}$ in Section 2. The derivative in the first term of Equation (12) is evaluated while holding the momentum coordinates constant. The direction cosines ℓ_i are equivalent to ${\boldsymbol{\Omega }}\cdot {{\boldsymbol{e}}}_{(i)}$ in Section 2 evaluated in the lab frame and can be written as

$$\begin{eqnarray}{{\ell }}_{(1)} & = & \cos \ \bar{\theta },\\ {{\ell }}_{(2)} & = & \sin \bar{\theta }\cos \ \bar{\phi },\\ {{\ell }}_{(3)} & = & \sin \bar{\theta }\sin \ \bar{\phi }.\end{eqnarray} \tag{ 13 }$$

The geometric coefficients ${\omega }_{(0)},{\omega }_{(\bar{\theta })}$, ${\omega }_{(\bar{\phi })}$ given in Nagakura et al. (2017a) reduce in flat spacetime to

$$\begin{eqnarray}{\omega }_{(0)} & = & 0,\\ {\omega }_{(\bar{\theta })} & = & -\displaystyle \frac{\sin \bar{\theta }}{r},\\ {\omega }_{(\bar{\phi })} & = & -\displaystyle \frac{\cot \theta }{r}{\sin }^{3}\bar{\theta }\sin \bar{\phi }.\end{eqnarray} \tag{ 14 }$$

The NSY code uses Lagrangian remapping grids (LRG) and laboratory fixed grids (LFG) in the fluid rest frame and the lab frame, respectively, to discretize the neutrino momentum space (Nagakura et al. 2014). The LFG are constructed so as to have an isotropic energy grid in the lab frame, while the LRG are constructed so as to have an isotropic energy grid in the fluid rest frame and a propagation angle-dependent energy grid in the lab frame. Here, isotropic means that each angular bin sees the same energy grid. This two-grid technique is essential to treating special relativistic effects in full generality in the DO method.

In regions where the scattering optical depth ${\tau }_{s}={\kappa }_{s}l$ is large, where l is the relevant length scale, direct MC radiation transport becomes very inefficient. The path length between scattering events is very small, so a great deal of computer time is spent performing these scattering events while there is little actual movement of energy and lepton number. In these regions, the neutrino transport is very well approximated as a diffusion process, a fact that we use to accelerate the computation.

In the past, MC neutrino transport schemes have excluded the inner regions of high optical depth in favor of an inner boundary condition (Janka 1991) or have employed the discrete diffusion MC approximation in these regions (Densmore et al. 2007; Abdikamalov et al. 2012). In order to keep the neutrino motion free of any specific grid geometry and to prevent a hard spatial boundary between two algorithms, we instead choose to implement the MC random walk approximation (Fleck & Canfield 1984). This treats neutrino motion over a specified distance as a diffusive process and relies on the assumption of isotropic, elastic scattering. In our implementation, we also assume the fluid is unchanging in space and in time during each diffusion step. Here we modify the method of Fleck & Canfield (1984) to treat static fluid backgrounds with relativistic fluid velocities.

The approximation accelerates MC transport in regions of high scattering optical depth using a solution to the diffusion equation:

$$\begin{eqnarray}&&{\partial }_{t}\psi ({\boldsymbol{r}},t)=D{{\rm{\nabla }}}^{2}\psi ({\boldsymbol{r}},t).\end{eqnarray} \tag{ 15 }$$

The diffusion constant can be shown to be $D=c/3{\kappa }_{s}$ by comparing the solution to the diffusion equation on an infinite, uniform background given initial conditions $\psi ({\boldsymbol{r}},0)={\delta }^{3}({\boldsymbol{r}})$ to the solution of a random walk process starting at ${\boldsymbol{r}}=0$ with step sizes determined by the scattering opacity ${\kappa }_{s}$ (Chandrasekhar 1943; Fleck & Canfield 1984). In the context of MC radiation transport, the solution $\psi ({\boldsymbol{r}},t)$ represents the probability density of the neutrino being at location ${\boldsymbol{r}}$ at time t.

Using the diffusion equation with this diffusion constant, we now specify a sphere of radius R in the comoving frame and derive the probability that a neutrino has escaped from the sphere after a certain time t. To do this, we again solve the diffusion equation, but this time with the boundary condition $\psi (R,t)=0$ to indicate that we are interested only in the first time a neutrino leaves the sphere, and we do not allow neutrinos to leave and then reenter the sphere. This can be solved via separation of variables and Sturm–Liouville orthogonality conditions to arrive at

$$\begin{eqnarray}&&\psi (r,t\gt 0)=\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{n}{2{R}^{2}}\displaystyle \frac{\sin (n\pi r/R)}{r}\exp \left[-{\left(\displaystyle \frac{n\pi }{R}\right)}^{2}{Dt}\right].\end{eqnarray} \tag{ 16 }$$

The probability that a neutrino has escaped the sphere after time t is represented by the volume integral of the diffusion solution (Figure 19):

$$\begin{eqnarray}&&{P}_{\mathrm{escape}}(R,t\gt 0)=1-{\displaystyle \int }_{0}^{R}\psi (r,t)4\pi {r}^{2}{dr}\\ &&\quad \quad =\,1-2\displaystyle \sum _{n=1}^{\infty }{(-1)}^{n-1}\exp \left[-{\left(\displaystyle \frac{n\pi }{R}\right)}^{2}{Dt}\right].\end{eqnarray} \tag{ 17 }$$

This solution is plotted in Figure 19.

Zoom In Zoom Out Reset image size

The diffusion equation is acausal in that there is a finite probability of a neutrino escaping at times less than the light travel time to the edge of the sphere. Because of this, we set ${P}_{\mathrm{escape}}(R,t\lt R/c)=0$. We can also use the escape probability at t = R/c as an estimate of the accuracy of the approximation. We only use the random walk approximation when

$$\begin{eqnarray}&&{P}_{\mathrm{escape}}(R,R/c)\lt \mathrm{tol}.\end{eqnarray} \tag{ 18 }$$

In this study, we use $\mathrm{tol}={10}^{-3}$, which corresponds to only using the random walk approximation when the scattering optical depth of the sphere is ${\kappa }_{s}R\geqslant 12$.

We tabulate ${P}_{\mathrm{escape}}(R,t)$, which can then be inverted via inverse transform sampling (e.g., Haghighat 2015) to randomly sample the escape time t_esc. The table extends over the range of $0\leqslant \chi \leqslant {\chi }_{\max }$ using 100 evenly spaced points in χ where $\chi ={Dt}/{R}^{2}$ and, in our calculations, ${\chi }_{\max }=2$ (corresponding to ${P}_{\mathrm{escape}}=0.997$). We evaluate the first 1000 terms in the sum in Equation (17), which is far more than is necessary for a converged solution, but tabulating P_escape is a very cheap one-time calculation.

We restrict the lab-frame radius of the sphere R_lab to the largest length scale between (1) the distance to the nearest grid cell boundary and (2) 1% of the grid cell's smallest dimension. However, since the sphere is defined in the fluid rest frame, its size must be further limited when the fluid is moving, since the sphere is effectively advected. The largest restriction occurs in the event that the displacement of the neutrino from its starting position to the surface of the sphere is parallel to the fluid velocity, so we will use this worst-case scenario to set the sphere size limiter. The four-vector ${d}_{\mathrm{com}}=\{t,{\boldsymbol{R}}\}$ connecting the neutrino's initial and final positions in the lab frame can be Lorentz transformed to give the displacement vector in the lab frame ${d}_{\mathrm{lab}}=\{\gamma (t+{vR}/{c}^{2}),\gamma ({\boldsymbol{R}}+{\boldsymbol{v}}t)\}$. The longest diffusion time the numerical scheme will allow is ${t}_{\max }={R}^{2}{\chi }_{\max }/D$, resulting in a maximum lab-frame displacement of ${R}_{\mathrm{lab},\max }=\gamma R(1+{Rv}{\chi }_{\max }/D)$. Inverting this, we set the comoving-frame radius to

$$\begin{eqnarray}&&R=\displaystyle \frac{2{R}_{\mathrm{lab}}}{\gamma }{\left(1+\sqrt{1+\displaystyle \frac{4{R}_{\mathrm{lab}}v{\chi }_{\max }}{\gamma D}}\right)}^{-1}.\end{eqnarray} \tag{ 19 }$$

The comoving-frame neutrino energy remains the same throughout the process since the scattering is assumed to be elastic. Absorption happens continuously throughout the diffusion process. The packet energy is decreased according to

$$\begin{eqnarray}&&{E}_{p}(t)={E}_{p}(0)\exp [-{\kappa }_{a}{ct}],\end{eqnarray} \tag{ 20 }$$

and ${E}_{p}(0)-{E}_{p}(t)$ is added to the fluid energy to account for neutrino absorption. The comoving frame packet energy averaged over the diffusion time is

$$\begin{eqnarray}{\bar{E}}_{p} & = & \displaystyle \frac{1}{t}{\displaystyle \int }_{0}^{t}{E}_{p,0}{e}^{-{\kappa }_{a}{ct}^{\prime} }{dt}^{\prime} \\ & = & \displaystyle \frac{{E}_{p,0}}{{ct}{\kappa }_{a}}(1-{e}^{-{\kappa }_{a}{ct}}).\end{eqnarray} \tag{ 21 }$$

If neutrino packets are created assuming the fluid emits for a time of $\delta {t}_{\mathrm{emit}}$, this means that the neutrino contributes ${\bar{E}}_{p}t/\delta {t}_{\mathrm{emit}}$ to the fluid cell's steady-state radiation energy content. Averaged over the diffusion process, most of the neutrino energy is distributed isotropically in terms of direction. However, there is a small asymmetry because the neutrino ends up at one point on the edge of the diffusion sphere. Averaged over the duration of the diffusion process, for a neutrino packet with energy E_p, there is a net energy flux of ${E}_{p}R/{ct}$ in the direction of the final displacement vector while ${E}_{p}({ct}-R)/{ct}$ is distributed isotropically in direction.

With the theoretical groundwork complete, we now describe the random walk algorithm itself. A comoving-frame diffusion sphere size R is first chosen according to Equation (19). If the scattering optical depth ${\kappa }_{s}R$ is sufficiently large (Equation (18)), the time the neutrino takes to reach the edge of the sphere t is sampled from Equation (17). A location at the edge of the comoving-frame sphere is randomly uniformly chosen, the displacement four-vector $\{t,{\boldsymbol{R}}\}$ is Lorentz transformed into the lab frame, and the neutrino is moved this distance. The new comoving neutrino direction is chosen uniformly in the $2\pi$ steradians moving strictly away from the diffusion sphere. The neutrino packet energy is decreased due to absorption according to Equation (20), and the absorbed energy is counted toward fluid heating. Comoving radiation energy in the amount of ${\bar{E}}_{p}R/c\delta {t}_{\mathrm{emit}}$ moving in the direction of the final displacement is Lorentz transformed into the lab frame and accumulated into the distribution function. The remaining ${\bar{E}}_{p}({ct}-R)/c\delta {t}_{\mathrm{emit}}$ of comoving radiation energy is divided evenly into N pieces, each is assigned an isotropically uniform random direction in the comoving frame, is Lorentz transformed into the lab frame, and is accumulated into the distribution function. This allows us to self-consistently treat both the isotropic and directional components of the radiation field without making reference to a particular grid structure. In this work, we found that N = 10 is a reasonable compromise between code performance and noise in the resulting radiation field.

The DO and MC methods are very different, so care is required to make meaningful comparisons between the two codes. The NSY code evolves the distribution function f, while Sedonu calculates the amount of neutrino energy in each spatial-energy-direction cell in nonnative calculations. The Sedonu distribution function value at the bin center $({r}_{a},{\theta }_{b},{\epsilon }_{k},{\bar{\theta }}_{m},{\bar{\phi }}_{n})$ is calculated using

$$\begin{eqnarray}&&{f}_{\mathrm{Sedonu},{abkmn}}=\displaystyle \frac{{\varepsilon }_{{akmn}}}{{V}_{{ab}}{\rm{\Delta }}{\left(\tfrac{{\epsilon }^{3}}{3}\right)}_{k}{\rm{\Delta }}\cos {(\bar{\theta })}_{m}{\rm{\Delta }}{\bar{\phi }}_{n}},\end{eqnarray} \tag{ 22 }$$

where ${\varepsilon }_{{akmn}}$ is the total neutrino energy content (units of erg) in spatial-direction-energy bin $\{{akmn}\},{V}_{a}$ is the spatial volume of the grid cell in the lab frame, $\bar{\theta }$ and $\bar{\phi }$ are the neutrino direction angles in the lab frame, and is the neutrino energy in the comoving frame.

In the 1D simulations where Sedonu collects radiation information on angular bins rather than native moments, we take care to ensure that the postprocessing for the two codes is equivalent. For both Sedonu and the NSY code, the distribution function is linearly interpolated to $\widetilde{f}$ on identical fine angular grids in $\{\cos (\bar{\theta }),\bar{\phi }\}$ of $\{80,40\}$ zones, respectively. Throughout this section, the subscript (a) refers to the spatial mesh index, the subscript (k) refers to the neutrino energy bin index, the subscripts (mn) refer to the direction indices on the coarse direction grid used in the transport calculation, the subscripts (pq) refer to the direction indices on the high-resolution postprocessing angular grid, and the superscripts (ij) refer to directions (i.e., $r,\theta$, or ϕ) in the lab frame.

The specific energy density (lab-frame energy density per comoving-frame neutrino energy) is computed as a sum over the coarse grid for Sedonu and over the fine grid in the NSY code. This takes advantage of the fact that Sedonu computes energy content directly and does not introduce interpolation error into the Sedonu results:

$$\begin{eqnarray}{E}_{\epsilon ,{ak},\mathrm{Sedonu}} & = & \displaystyle \frac{1}{{V}_{a}{\rm{\Delta }}{\left(\tfrac{{\epsilon }^{3}}{3}\right)}_{k}}\displaystyle \sum _{m}\displaystyle \sum _{n}{\varepsilon }_{{akmn}},\\ {E}_{\epsilon ,{ak},\mathrm{DO}} & = & {\epsilon }_{k}\displaystyle \sum _{p}\displaystyle \sum _{q}{\widetilde{f}}_{{akpq}}{\rm{\Delta }}\cos {({\bar{\theta }}_{\nu })}_{p}{\rm{\Delta }}{({\bar{\phi }}_{\nu })}_{q}.\end{eqnarray} \tag{ 23 }$$

For both Sedonu and the NSY code, the higher-order moments are evaluated as

$$\begin{eqnarray}{F}_{\epsilon ,{ak}}^{i} & = & {\epsilon }_{k}\displaystyle \sum _{p}\displaystyle \sum _{q}{\widetilde{f}}_{{akpq}}({{\boldsymbol{\Omega }}}_{{pq}}\cdot {{\boldsymbol{e}}}_{(i)}){\rm{\Delta }}\cos {(\bar{\theta })}_{p}{\rm{\Delta }}{(\bar{\phi })}_{q},\\ {P}_{\epsilon ,{ak}}^{{ij}} & = & {\epsilon }_{k}\displaystyle \sum _{p}\displaystyle \sum _{q}{\widetilde{f}}_{{akpq}}({{\boldsymbol{\Omega }}}_{{pq}}\cdot {{\boldsymbol{e}}}_{(i)})\\ & & \times ({{\boldsymbol{\Omega }}}_{{pq}}\cdot {{\boldsymbol{e}}}_{(j)}){\rm{\Delta }}\cos {(\bar{\theta })}_{p}{\rm{\Delta }}{(\bar{\phi })}_{q}.\end{eqnarray} \tag{ 24 }$$

The energy-integrated moments $M=\{E,{F}^{i},{P}^{{ij}}\}$ are computed using

$$\begin{eqnarray}&&{M}_{a}=\displaystyle \sum _{k}{\rm{\Delta }}{\left(\displaystyle \frac{{\epsilon }^{3}}{3}\right)}_{k}{M}_{\epsilon ,{ak}}.\end{eqnarray} \tag{ 25 }$$

Finally, the average neutrino energy is computed using

$$\begin{eqnarray}&&{\bar{\epsilon }}_{a}=\displaystyle \frac{{E}_{a}}{{N}_{a}}.\end{eqnarray} \tag{ 26 }$$

The three source terms on the right-hand side of Equation (1) each encapsulate multiple processes and are grouped into the mathematical nature of each term. In both Sedonu and the NSY code, all of these source terms are evaluated in the comoving frame. Details of how the NSY code computes reaction rates are explained by Bruenn (1985) and Sumiyoshi et al. (2005).

The emission and absorption term takes the form of

$$\begin{eqnarray}&&{\left[\displaystyle \frac{\partial f}{c\partial t}\right]}_{\mathrm{em} \mbox{-} \mathrm{abs}}={R}_{\mathrm{em}}(\epsilon )(1-f)-{R}_{\mathrm{abs}}(\epsilon )f,\end{eqnarray} \tag{ 27 }$$

where R_em and R_abs are the emission and absorption reaction rates, respectively. Sedonu takes advantage of the concept of stimulated absorption to account for final-state neutrino blocking (Burrows et al. 2006), in which the effective absorption reaction rate is ${\widetilde{R}}_{\mathrm{abs}}={R}_{\mathrm{abs}}+{R}_{\mathrm{emis}}$. This removes the need to treat final-state blocking explicitly in the neutrino emission process.

The scattering term accounts for neutrinos scattering into and out of a given direction according to

$$\begin{eqnarray}{\left[\displaystyle \frac{\partial f}{c\partial t}\right]}_{\mathrm{scat}} & = & \displaystyle \int d{\boldsymbol{\Omega }}^{\prime} \displaystyle \int d\left(\displaystyle \frac{\epsilon {{\prime} }^{3}}{3}\right)\\ & & \times [{R}_{\mathrm{scat}}(\epsilon ^{\prime} ,\epsilon ,{\boldsymbol{\Omega }}^{\prime} \cdot {\boldsymbol{\Omega }})f^{\prime} (1-f)\\ & & -{R}_{\mathrm{scat}}(\epsilon ,\epsilon ^{\prime} ,{\boldsymbol{\Omega }}\cdot {\boldsymbol{\Omega }}^{\prime} )f(1-f^{\prime} )].\end{eqnarray} \tag{ 28 }$$

The primed variables are the neutrino final-state quantities, and R_scat is the scattering reaction rate. Both Sedonu and the NSY code assume isoenergetic scattering, so the scattering reaction rate becomes ${R}_{\mathrm{scat}}(\epsilon ^{\prime} ,\epsilon ,{\boldsymbol{\Omega }}^{\prime} \cdot {\boldsymbol{\Omega }})=\delta (\epsilon ,\epsilon ^{\prime} ){\widetilde{R}}_{\mathrm{scat}}(\epsilon ,{\boldsymbol{\Omega }}\cdot {\boldsymbol{\Omega }}^{\prime} )$. Under this assumption, the scattering source term reduces to

$$\begin{eqnarray}&&{\left[\displaystyle \frac{\partial f}{c\partial t}\right]}_{\mathrm{scat}}=\int d{\boldsymbol{\Omega }}^{\prime} {\widetilde{R}}_{\mathrm{scat}}(\epsilon ,{\boldsymbol{\Omega }}\cdot {\boldsymbol{\Omega }}^{\prime} )(f^{\prime} -f).\end{eqnarray} \tag{ 29 }$$

Sedonu uses $\int d\mu {\widetilde{R}}_{\mathrm{scat}}(\epsilon ,\mu )$, where $\mu ={\boldsymbol{\Omega }}\cdot {\boldsymbol{\Omega }}^{\prime}$, as the scattering opacity directly.

Finally, pair annihilation and neutrino bremsstrahlung source terms take the form of

$$\begin{eqnarray}&&{\left[\displaystyle \frac{\partial f}{c\partial t}\right]}_{\mathrm{pair} \mbox{-} \mathrm{brem}}=\displaystyle \int d\bar{{\boldsymbol{\Omega }}}\displaystyle \int d\left(\displaystyle \frac{{\bar{\epsilon }}^{3}}{3}\right)\\ &&\qquad \times \,[{R}_{\mathrm{pair} \mbox{-} \mathrm{brem},\mathrm{emis}}(\epsilon ,\bar{\epsilon },{\boldsymbol{\Omega }}\cdot \bar{{\boldsymbol{\Omega }}})(1-f)(1-\bar{f})\\ &&\qquad -\,{R}_{\mathrm{pair} \mbox{-} \mathrm{brem},\mathrm{abs}}(\epsilon ,\bar{\epsilon },{\boldsymbol{\Omega }}\cdot \bar{{\boldsymbol{\Omega }}})f\bar{f}].\end{eqnarray} \tag{ 30 }$$

The barred variables are the neutrino antispecies quantities, and ${R}_{\mathrm{pair} \mbox{-} \mathrm{brem},\mathrm{emis}}$ is the reaction rate for pair and bremsstrahlung processes. In order to ensure the same assumptions go into both radiation transport schemes, the NSY code calculates these reactions assuming the antispecies is isotropic:

$$\begin{eqnarray}{\bar{f}}_{\mathrm{iso}} & = & \displaystyle \frac{1}{4\pi }\displaystyle \int d\bar{{\boldsymbol{\Omega }}}\bar{f},\\ {R}_{\mathrm{pair} \mbox{-} \mathrm{brem},\mathrm{iso}}(\epsilon ,\bar{\epsilon }) & = & \displaystyle \int d\bar{{\boldsymbol{\Omega }}}{R}_{\mathrm{pair} \mbox{-} \mathrm{brem}}(\epsilon ,\bar{\epsilon },{\boldsymbol{\Omega }}\cdot \bar{{\boldsymbol{\Omega }}}),\end{eqnarray} \tag{ 31 }$$

where ${\bar{f}}_{\mathrm{iso}}$ depends only on energy and not on direction. Under this assumption, the source term can be written as

$$\begin{eqnarray}\begin{array}{rcl}{\left[\displaystyle \frac{\partial f}{c\partial t}\right]}_{\mathrm{pair} \mbox{-} \mathrm{brem}} & = & \mathop{\overbrace{\displaystyle \int d\left(\displaystyle \frac{{\bar{\epsilon }}^{3}}{3}\right){R}_{\mathrm{pair} \mbox{-} \mathrm{brem},\mathrm{emis},\mathrm{iso}}(\epsilon ,\bar{\epsilon })(1-{\bar{f}}_{\mathrm{iso}})}}\limits^{{\widetilde{R}}_{\mathrm{pair} \mbox{-} \mathrm{brem},\mathrm{emis}}}\\ & & -f\mathop{\underbrace{\displaystyle \int d\left(\displaystyle \frac{{\bar{\epsilon }}^{3}}{3}\right)[{R}_{\mathrm{pair} \mbox{-} \mathrm{brem},\mathrm{emis},\mathrm{iso}}(\epsilon ,\bar{\epsilon })(1-{\bar{f}}_{\mathrm{iso}})+{R}_{\mathrm{pair} \mbox{-} \mathrm{brem},\mathrm{abs},\mathrm{iso}}(\epsilon ,\bar{\epsilon }){\bar{f}}_{\mathrm{iso}}]}}\limits_{{\widetilde{R}}_{\mathrm{pair} \mbox{-} \mathrm{brem},\mathrm{abs}}}.\end{array}\end{eqnarray} \tag{ 32 }$$

Sedonu uses ${\widetilde{R}}_{\mathrm{pair} \mbox{-} \mathrm{brem},\mathrm{abs}/\mathrm{emis}}$ in the same way as ${\widetilde{R}}_{\mathrm{abs}/\mathrm{emis}}$.

Since the NSY code evolves f, they use the reaction rates (units of cm⁻¹) directly, but Sedonu needs to convert the emission reaction rates to physical emissivities. For an emissivity η with units of (erg cm⁻³ s⁻¹),

$$\begin{eqnarray}&&\eta ={\widetilde{R}}_{\mathrm{emis}}\displaystyle \frac{{\epsilon }_{i}}{{c}^{2}{h}^{3}}{\rm{\Delta }}{\left(\displaystyle \frac{{\epsilon }^{3}}{3}\right)}_{i}.\end{eqnarray} \tag{ 33 }$$

Here, ${\rm{\Delta }}{({\epsilon }^{3}/3)}_{i}$ and ${\epsilon }_{i}$ are the momentum-space volume (normalized by $4\pi$) and center of energy bin i, respectively. The absorption and scattering reaction rates with tildes ($\widetilde{R}$) are already equivalent to the absorption and scattering opacities.