Multidimensional Boltzmann Neutrino Transport Code in Full General Relativity for Core-collapse Simulations

Massive stars (≳8 M_⊙) end their lives with gravitational collapse and eventually form compact objects such as neutron stars (NSs) and black holes (BHs; Janka 2012; Burrows 2013; Foglizzo et al. 2015; O'Connor 2017, for reviews). Such a core-collapse event is a complicated phenomenon governed by nonlinear equations, and its accurate modeling is obtained only by numerical simulations. Ever since the first detailed study of the core-collapse events by Colgate & White (1966), many attempts have been made to model their dynamics and observational signals. Although successfully exploding models are not rare these days, the results do not agree with one another among different groups; moreover, the explosion energy obtained by the numerical simulations is commonly small, only ∼10% of the typical observed values (O'Connor & Couch 2018; Vartanyan et al. 2018; Burrows et al. 2019a, 2019b; Kuroda et al. 2020; Müller & Varma 2020; but see also Bollig et al. 2020; Burrows & Vartanyan 2021). Such discrepancies may be attributed to various approximations employed in the simulations, and there is a need for first-principles calculations, particularly of neutrino transfer without artificial approximations.

The numerical handling of neutrino transport is indeed the important component for the core-collapse simulations. In fact, neutrinos carry most (∼99%) of the energy released by the gravitational collapse, and the interactions between neutrinos and matter are believed to be crucial for the supernova dynamics (Martinez-Pinedo et al. 2016; Halzen & Scholberg 2017; Mezzacappa et al. 2020). In modern simulations, the shock wave generated by the core bounce stalls on the way to the core surface because of the photodissociation of nuclei, and most researchers think that the neutrino heating is the key process for shock revival to obtain successful explosions.

Neutrinos are not in thermal equilibrium except deep inside the core, and the neutrino distribution in momentum space needs to be calculated in principle. The neutrino transfer is described with the Boltzmann equations in the 6D phase space (Lindquist 1966; Ehlers 1971). There have been some attempts to directly solve them. The discrete ordinate (S_N ) method (Mezzacappa & Bruenn 1993; Liebendörfer et al. 2001; Sumiyoshi et al. 2005; Sumiyoshi & Yamada 2012; Nagakura et al. 2014) discretizes the entire phase space and solves the equation with a finite difference. On the other hand, Chan & Müller (2020) proposed an "intuitive particle-like approach" for the calculation of the angular advection in momentum space, instead of the direct finite differencing. Codes based on some spectral methods are also constructed; Radice et al. (2013) employed the spherical harmonics (P_N ) to decompose the distribution function in momentum space, whereas the distribution functions in both the configuration and momentum spaces are expanded with some basis functions by Peres et al. (2014).

The Monte Carlo method (Fleck & Cummings 1971; Abdikamalov et al. 2012; Richers et al. 2017) adopts a probabilistic description of neutrino transport. Monte Carlo transport is known to be noisy if the number of sample particles N is not large enough, with the error scaling as $1/\sqrt{N}$. On the other hand, it can treat sharp angular distributions, unlike the finite-difference schemes. Richers et al. (2017) made detailed comparisons between the Monte Carlo and discrete ordinate methods and summarized the merits and demerits of the two methods.

Unfortunately, the Boltzmann transport schemes mentioned above are all very computationally demanding. Hence, most core-collapse simulations have been performed by employing some phenomenological approximations: the flux-limited diffusion approximation (Arnett 1977), the isotropic diffusion source approximation (Liebendörfer et al. 2009), the ray-by-ray-plus approximation (Buras et al. 2006), and the truncated moment approximation (Thorne 1981; Shibata et al. 2011), with some analytical closure relations such as the M1 closure (Levermore 1984) being imposed by hand. Roughly speaking, these schemes reduce the computational cost by neglecting some information on the distribution in momentum space. Several studies were conducted to compare these approximations (Skinner et al. 2016; Cabezón et al. 2018; Just et al. 2018; O'Connor et al. 2018; Pan et al. 2018) and demonstrated that the different approximations may lead to quantitatively different results. This indicates that we need core-collapse simulations without appealing to these approximations, even if they are costly.

We should also remind ourselves that general relativity is another important ingredient. Core-collapse simulations with general relativistic neutrino transport, however, have been performed mostly with the approximate transport schemes (Müller et al. 2010, 2012; Kuroda et al. 2012; Ott et al. 2013; Kuroda et al. 2016; Roberts et al. 2016; O'Connor & Couch 2018; Kuroda et al. 2020). In fact, Boltzmann simulations with the general relativistic discrete ordinate method are limited to spherical symmetry (Yamada et al. 1999; Liebendörfer et al. 2001; Sumiyoshi et al. 2005). Although Monte Carlo neutrino transport codes in full general relativity have been developed in the context of NS merger and collapsar simulations (Foucart et al. 2018; Miller et al. 2019a, 2019b; Foucart et al. 2020), they have yet to be applied to the core-collapse simulations.

The general relativistic radiative transport is a common problem for computational high-energy astrophysics. General relativistic photon transport schemes have been developed in the context of BH accretion, and Monte Carlo photon transport codes have been constructed (Dolence et al. 2009; Schnittman & Krolik 2013). The long (Takahashi & Umemura 2017 and references therein) or short (Mihalas et al. 1978; Kunasz & Auer 1988) characteristics methods are an alternative that utilizes the transfer equation along the geodesic curves. Although the former is known to be very accurate, it is also computationally demanding, and the latter is sometimes preferred. A hybrid method (Zhu et al. 2015) was also developed.

We have developed a multidimensional Boltzmann radiation hydrodynamics code based on the discrete ordinate method. Sumiyoshi & Yamada (2012) developed a Newtonian transport code for the static background. It was later extended to incorporate special relativity to all orders of v/c implementing the two-energy-grid scheme in Nagakura et al. (2014). Nagakura et al. (2017) further extended the code to track the proper motion of a proto-NS by installing a moving-coordinates technique that makes use of the gauge degree of freedom in general relativity, since the code employs the conservative form of the general relativistic Boltzmann equation written in the 3 + 1 decomposition formulation in Shibata et al. (2014). Note, however, that this was realized only in the flat spacetime by a time-dependent uniform shift vector. We have run a series of core-collapse supernova (CCSN) simulations with this code (Nagakura et al. 2018, 2019; Harada et al. 2019, 2020; Iwakami et al. 2020), although they are limited to the flat spacetime.

In this paper, we report the results of some tests conducted for the fully general relativistic version of our code. We assess its capabilities regarding both the advection in phase space and the handling of interactions. Although the tests are limited to 1D (spherical symmetry) or 2D (axisymmetry) in this paper, our code is applicable to 3D calculations. Although we employ the spacetime metrics given analytically in this paper, it is applicable to any dynamical spacetime, such as those encountered in BH formations.

This paper is organized as follows. In Section 2, we describe the 3 + 1 formulation of the Boltzmann equation. The numerical method employed in our code is explained in Section 3. The results of the advection tests in the Schwarzschild and Kerr spacetimes are presented in Sections 4 and 5, respectively. We then show the results of the tests on the collision terms in Section 6. Finally, we summarize our findings and discuss future prospects in Section 7. Throughout the paper, we use the metric signature − +++; Greek (α, β, μ, ν) and Latin (i, j, k) indices run over 0–3 and 1–3, respectively.

The neutrino distribution in the curved spacetime is described by the general relativistic Boltzmann equation (Lindquist 1966; Ehlers 1971 for references),

$$\begin{eqnarray}&&{p}^{\mu }\displaystyle \frac{\partial f}{\partial {x}^{\mu }}-{{\rm{\Gamma }}}_{\mu \nu }^{i}{p}^{\mu }{p}^{\nu }\displaystyle \frac{\partial f}{\partial {p}^{i}}=-{p}^{\mu }{u}_{\mu }{S}_{\mathrm{rad}},\end{eqnarray} \tag{ 1 }$$

where f is the distribution function in phase space, and x^α and p^μ are the spacetime coordinates and four-momentum, respectively. Note that the spatial components of the four-momentum serve as the coordinates in momentum space; ${{\rm{\Gamma }}}_{\mu \nu }^{\alpha }$ is the Christoffel symbol, u^μ is the four-vector of matter, and S_rad is the collision term. The first and second terms on the left-hand side of the equation describe the neutrino advection in space and momentum space, respectively.

From a numerical point of view, it is desirable to cast the Boltzmann equation in the conservative form. Shibata et al. (2014) gave such a formulation to the general relativistic Boltzmann equation,

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

where the energy is defined as $\epsilon =-{p}_{\mu }{e}_{(0)}^{\mu };$ θ_ν and ϕ_ν are the zenith and azimuth angles in momentum space, respectively, as depicted in Figure 1; g is the determinant of the metric tensor; and the tetrad basis is specified by ${e}_{(\mu )}^{\alpha }$. In this paper, neutrinos are assumed to be massless, with their minuscule masses being neglected. Directional cosines in momentum space l_(i) are expressed as

$$\begin{eqnarray*}\begin{array}{rcl}{l}_{(1)} & = & \cos {\theta }_{\nu },\\ {l}_{(2)} & = & \sin {\theta }_{\nu }\cos {\phi }_{\nu },\\ {l}_{(3)} & = & \sin {\theta }_{\nu }\sin {\phi }_{\nu }.\end{array}\end{eqnarray*}$$

The factors ω₀, ${\omega }_{({\theta }_{\nu })}$, and ${\omega }_{({\phi }_{\nu })}$ are defined as

$$\begin{eqnarray*}\begin{array}{rcl}{\omega }_{(0)} & \equiv & {\epsilon }^{-2}{p}^{\mu }{p}_{\nu }{{\rm{\nabla }}}_{\mu }{e}_{(0)}^{\mu },\\ {\omega }_{({\theta }_{\nu })} & \equiv & \displaystyle \sum _{i=1}^{3}{\omega }_{(i)}\displaystyle \frac{\partial {l}_{(i)}}{\partial {\theta }_{\nu }},\\ {\omega }_{({\phi }_{\nu })} & \equiv & \displaystyle \sum _{i=2}^{3}{\omega }_{(i)}\displaystyle \frac{\partial {l}_{(i)}}{\partial {\phi }_{\nu }},\\ {\omega }_{i} & \equiv & {\epsilon }^{-2}{p}^{\mu }{p}_{\nu }{{\rm{\nabla }}}_{\mu }{e}_{(i)}^{\nu }.\end{array}\end{eqnarray*}$$

A natural choice of the tetrad basis may be given as

$$\begin{eqnarray}\begin{array}{rcl}{e}_{(0)}^{\mu } & = & {n}^{\mu },\\ {e}_{(1)}^{\mu } & = & {\gamma }_{{rr}}^{-1/2}{\left(\displaystyle \frac{\partial }{\partial r}\right)}^{\mu },\\ {e}_{(2)}^{\mu } & = & -\displaystyle \frac{{\gamma }_{r\theta }}{\sqrt{{\gamma }_{{rr}}({\gamma }_{{rr}}{\gamma }_{\theta \theta }-{\gamma }_{r\theta }^{2})}}{\left(\displaystyle \frac{\partial }{\partial r}\right)}^{\mu }\\ & & +\sqrt{\displaystyle \frac{{\gamma }_{{rr}}}{{\gamma }_{{rr}}{\gamma }_{\theta \theta }-{\gamma }_{r\theta }^{2}}}{\left(\displaystyle \frac{\partial }{\partial \theta }\right)}^{\mu },\\ {e}_{(3)}^{\mu } & = & \displaystyle \frac{{\gamma }^{r\phi }}{\sqrt{{\gamma }^{\phi \phi }}}{\left(\displaystyle \frac{\partial }{\partial r}\right)}^{\mu }+\displaystyle \frac{{\gamma }^{\theta \phi }}{\sqrt{{\gamma }^{\phi \phi }}}{\left(\displaystyle \frac{\partial }{\partial \theta }\right)}^{\mu }\\ & & +\sqrt{{\gamma }^{\phi \phi }}{\left(\displaystyle \frac{\partial }{\partial \phi }\right)}^{\mu }.\end{array}\end{eqnarray} \tag{ 3 }$$

Here n^μ and γ_ij are the time-like unit vector normal to the time-constant hypersurface and the spatial metric on it, respectively. We choose this normal vector n^μ as the time-like base, since the zeroth components of all of the space-like bases are vanishing then. Two of the space-like bases are chosen so that they span the meridian plane, with one of them being aligned with the radial coordinate; the last one is perpendicular to the meridian plane. The components of the normal vector can be explicitly written as

$$\begin{eqnarray}&&{n}^{\mu }=(1/\alpha ,-{\beta }^{i}/\alpha ),\end{eqnarray} \tag{ 4 }$$

where α is the lapse function and βⁱ is the shift vector. In the 3 + 1 formulation, the line element is expressed as

$$\begin{eqnarray}&&{{ds}}^{2}=(-{\alpha }^{2}+{\beta }^{k}{\beta }_{k}){{dt}}^{2}+2{\beta }_{i}{{dtdx}}^{i}+{\gamma }_{{ij}}{{dx}}^{i}{{dx}}^{j}.\end{eqnarray} \tag{ 5 }$$

In numerical computations, the values of the metric variables need to be evaluated both on the cell interfaces and in the centers. For the metrics given analytically as in Sections 4 and 5, they are just evaluated from those analytic formulae either on the interfaces or at the centers. For the test in Section 6 in which the metric is given numerically only at the cell interfaces from another simulation, the cell center values are calculated by linear interpolation.

Zoom In Zoom Out Reset image size

3.1. Extensions

3.2. Matter Background

In this and the next sections, we test our code's capability to treat neutrino advection in both space and momentum space. We switch off all neutrino reactions throughout these sections; hence, the background gas plays no role. The β-parameter, which is defined in Equation A3 to control the amount of upwinding in the numerical derivatives, becomes β = 1 in this case. It means that the spatial advection terms are completely upwind-differenced (see Appendix A for details).

In this section, we employ a Schwarzschild exterior solution, i.e., a spherically symmetric stationary solution of the Einstein equation for a vacuum, as the background spacetime. It is given as

$$\begin{eqnarray}\begin{array}{rcl}{{ds}}^{2} & = & -\left(1-\displaystyle \frac{2{GM}}{{c}^{2}r}\right){c}^{2}{{dt}}^{2}+{\left(1-\displaystyle \frac{2{GM}}{{c}^{2}r}\right)}^{-1}{{dr}}^{2}\\ & & +{r}^{2}(d{\theta }^{2}+{\sin }^{2}\theta d{\phi }^{2}),\end{array}\end{eqnarray} \tag{ 6 }$$

where t, r, θ, and ϕ are the spacetime coordinates; M is the mass of the central object; and c and G are the light speed and gravitational constant, respectively. The geodesic curves along which noninteracting neutrinos move in the Schwarzschild solution are described by a simple equation as given in Appendix B. We numerically evaluate Equation (3) to obtain the tetrad components, which are analytically reduced in the present case to

$$\begin{eqnarray}&&{e}_{(0)}^{0}={\alpha }^{-1},\quad {e}_{(1)}^{1}={\gamma }_{{rr}}^{-1},\quad {e}_{(2)}^{2}=\displaystyle \frac{1}{r},\quad {e}_{(3)}^{3}=\displaystyle \frac{1}{r\sin \theta }.\end{eqnarray} \tag{ 7 }$$

We set the gravity source to be a sphere with a constant density of ρ = 1 × 10¹⁵ g cm⁻³ and a radius of R = 12 km; hence, the mass is M = 3.62 M_⊙, whose Schwarzschild radius is r = 10.8 km, and the photon sphere radius is 16 km. Inside 12 km, we use the Schwarzschild interior solution for the uniform sphere of the same mass. The modification is not important, since we consider the neutrino propagation only outside this region in this section.

4.1. Spherically Symmetric Tests

We present the results of advection tests under the assumption of spherical symmetry in the neutrino distribution function; i.e., it depends only on r, , and θ_ν . Although our code is multidimensional, we intentionally suppress the angular degrees of freedom in space and treat the radial advection alone in space (plus the advections in momentum space) for the tests in Section 4.1. We discuss the energy and angular advection in momentum space in Sections 4.1.1 and 4.1.2, respectively.

4.1.1. Energy Advection Tests

In the gravitational field, neutrinos change their energy as they move. In the computation of neutrino transport, it is important to take such effects into account because the neutrino interactions strongly depend on the neutrino energy.

In order to study the capability of our code to treat those energy changes, we conduct the following test. We fix the neutrino distribution function to f = 1 on a single energy bin at a certain radius, which serves as a monochromatic and spherical neutrino source. We set f = 0 elsewhere initially. Then, neutrinos will flow out of this source and fill the space, the evolution of which we will compute with our new code. As for the initial angular distribution in momentum space, we assume that the neutrinos move in a single direction with the zenith angle θ_ν = 0 (outward) or π (inward). We can test redshift in the former and blueshift in the latter. Note, however, that it is difficult in our code to strictly set the single-angle distributions given above. We hence set f = 1 on either the first or the last angular bin and f = 0 on other bins in the numerical test. In order to focus on the energy advection, we switch off angular advection in this test.

The results are compared with the analytical formula for the neutrino energy _ana as a function of radius,

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{ana}}(r)=\left(1-\displaystyle \frac{2{GM}}{{c}^{2}{r}_{\mathrm{source}}}\right){\left(1-\displaystyle \frac{2{GM}}{{c}^{2}r}\right)}^{-1}{\epsilon }_{\mathrm{source}},\end{eqnarray} \tag{ 8 }$$

where _source and r_source are the energy and radius of the source, respectively.

It should be stressed that for mesh-based codes like ours, this is a very challenging problem, with sharp edges existing both in the energy and angular distributions. As will be witnessed later, rather large numerical diffusions occur inevitably. We choose this test, though, since this enables us to see most clearly if the code can reproduce the redshift/blueshift of neutrino energy as it moves in the gravitational well; the resolution dependence also manifests itself.

Throughout this test, we deploy the radial mesh with N_r = 128 grid points that cover the range r ∈ [0, 100] km. It is finer in the region r ∈ [10, 50] km. The number of angular mesh points in momentum space is ${N}_{{\theta }_{\nu }}=20$. The energy mesh has logarithmically spaced grid points and covers the range ∈ [0, 50] MeV. We vary the number of energy mesh points as N = 20, 30, 40, and 60 to study the resolution dependence of the result.

Figure 2 shows the energy distributions as a function of radius at a certain time for N = 60. The white dashed curves depict the analytical solutions. The left panel presents the result of the gravitational redshift test at the coordinate time t = 2 × 10⁻⁴ s. In this calculation, the neutrino source is located at r = 15 km and emits monochromatic neutrinos with = 20 MeV outward. The energy distributions obtained numerically trace the analytical curve, although they are somewhat broadened due to numerical diffusions arising from the finite energy resolution in the simulation. We stress again that this is actually a very challenging problem for finite-difference methods like ours. In fact, the single energy bin has a finite width and cannot express the monochromatic energy distribution very well in the first place. The same trend is found for the blueshift test, as exhibited in the right panel. The neutrino source is located at r = 40 km and emits neutrinos with = 10 MeV inward in this test. The energy of the neutrinos increases indeed as they propagate radially inward. It is also observed that the numerical diffusion is weaker at large radii, as the advection is rather small there. In both panels, the analytic curves are drawn from the source positions to the points that the massless neutrino reaches at the given time. It is obvious that the terminal points are well reproduced by the numerical computations.

Zoom In Zoom Out Reset image size

In order to see the resolution dependence of the numerical diffusion quantitatively, we repeat the redshift tests with different numbers of energy bins. We quantify the numerical diffusion by defining the following error function:

$$\begin{eqnarray}&&{E}_{\epsilon }(r)\equiv \displaystyle \frac{{\displaystyle \sum }_{n=1}^{{N}_{\epsilon }}f(r,{\epsilon }_{n}){\left({\epsilon }_{n}-{\epsilon }_{\mathrm{ana}}(r)\right)}^{2}d{\epsilon }_{n}}{{\left({\epsilon }_{\mathrm{ana}}(r)\right)}^{2}{\displaystyle \sum }_{n=1}^{{N}_{\epsilon }}f(r,{\epsilon }_{n})d{\epsilon }_{n}},\end{eqnarray} \tag{ 9 }$$

where _n and d _n are the value of energy at the nth cell center and the width of the same cell, respectively.

In Figure 3, we show the radial profiles of the error function rescaled by the number of energy mesh points N . It is seen that the rescaled error functions for the different energy resolutions almost coincide with one another, except for N = 20. This indicates that the error function is inversely proportional to N , roughly implying the first-order convergence. This is expected, since the energy advection term is evaluated with a first-order finite-difference scheme as described in Appendix A.

Zoom In Zoom Out Reset image size

4.1.2. Angular Advection Tests

The greatest advantage of directly solving the Boltzmann equation is that we are able to obtain information not only on energy but also on the angular distribution in momentum space. The direction of the neutrino momentum is specified by the zenith and azimuth angles, (θ_ν , ϕ_ν ; see Figure 1). Note that the distribution function depends on θ_ν alone in the spherical symmetry assumed in this section. As a neutrino moves nonradially, the zenith angle θ_ν , which is measured from the local radial direction, changes even in the flat spacetime. This angular advection is shown schematically in Figure 4. The blue curve is one of the geodesic curves along which the free neutrino moves in the Schwarzschild spacetime. Note that it is no longer a straight line due to gravity. In this example, the neutrino moves outward and the zenith angle approaches θ_ν = 0, i.e., the outward radial direction, with the increasing radius r. Since the geodesic curve is bent inward by gravity, the approach is slower for the Schwarzschild spacetime than in the flat spacetime. In this subsection, we test the capability of our code to reproduce this angular advection.

Zoom In Zoom Out Reset image size

The numerical setting is essentially the same as in the previous test for the energy advection; we put the monochromatic neutrino source uniformly on a sphere with a certain radius by setting f = 1 on an single energy bin there and f = 0 otherwise. The difference is that we choose θ_ν = π/2, which corresponds to a direction of p^r = 0, at the source. Note that we actually set f = 1 on a single angular bin nearest to θ_ν = π/2 for numerical convenience. We vary the source radius to investigate the angular advections both inside and outside the photon sphere, i.e., the circular orbit. The neutrinos emitted outside the photon sphere with p^r = 0 propagate outward with θ_ν decreasing monotonically to zero, whereas those emitted inside the photon sphere go inward with θ_ν increasing as they propagate. In this test, we switch off the energy advection, which should occur simultaneously in reality. This is to avoid numerical diffusions in both angle and energy at the same time. The results are compared with those for the flat spacetime, as well as with the reference solution obtained in Appendix B.

We employ the same radial mesh with N_r = 128 as in the previous tests. As for the mesh in momentum space, we vary the number of grid points as ${N}_{{\theta }_{\nu }}=20$, 30, and 40 to see the resolution dependence. We set ${N}_{{\phi }_{\nu }}=2$ for numerical convenience, although the distribution does not depend on ϕ_ν in the present case.

Figure 5 shows the angular distributions of neutrinos in momentum space as a function of the radius r for the angular advection tests with ${N}_{{\theta }_{\nu }}=40$. The white dashed curves depict the reference geodesic curves, truncated at the radius that the massless particles reach at the time of the snapshot. The left panel shows the result for the flat spacetime at t = 1 × 10⁻⁴ s, with the source placed at r = 20 km. As mentioned earlier, there is a nontrivial angular advection even in the flat spacetime, since we deploy the polar coordinates in space. Neutrinos emitted with p_r = 0 always move outward in the flat spacetime, and the zenith angle monotonically converges to θ_ν = 0 as they go outward. The middle panel presents the result at t = 2.5 × 10⁻⁴ s for the Schwarzschild spacetime. Note that the neutrino source is located outside the photon sphere, the radius of which is 16 km in the present test. The neutrinos emitted from this source with p_r = 0 move outward. The outward radial propagation of neutrinos is slower in this case than in the flat case because of the geodesic deflection and gravitational time delay. As a result, the radius-angle curve for the Schwarzschild spacetime is less steep than that for the flat spacetime. The right panel is the result at t = 2 × 10⁻⁴ s with the neutrino source located at r = 15 km, i.e., inside the photon sphere in the Schwarzschild spacetime. This time, the trajectory is directed radially inward, and, as a result, θ_ν approaches θ_ν = π, instead of 0. The distributions are consistent with the reference curves, although there are some numerical diffusions. Figure 6 shows the radial slice of angular distribution at r = 40 km for the test in Schwarzschild spacetime. The numerical result obtained with ${N}_{{\theta }_{\nu }}=20$ is shown in blue, while the analytical value is shown in green. The height of the peak is lowered due to the numerical diffusion; however, the position of the peak is consistent with the analytical value.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Just as in the energy advection tests in Section 4.1.1, we quantify the numerical diffusion as follows. We define the error function as

$$\begin{eqnarray}&&{E}_{{\theta }_{\nu }}(r)\equiv \displaystyle \frac{{\displaystyle \sum }_{n=1}^{{N}_{{\theta }_{\nu }}}f(r,{\left({\theta }_{\nu }\right)}_{n}){\left({\left({\theta }_{\nu }\right)}_{n}-{\left({\theta }_{\nu }\right)}_{\mathrm{ref}}\right)}^{2}{\left(d{\theta }_{\nu }\right)}_{n}}{{\displaystyle \sum }_{n=1}^{{N}_{{\theta }_{\nu }}}f(r,{\left({\theta }_{\nu }\right)}_{n})d{\theta }_{\nu }},\end{eqnarray} \tag{ 10 }$$

where ${({\theta }_{\nu })}_{n}$ and ${(d{\theta }_{\nu })}_{n}$ are the value of θ_ν at the nth cell center and the width of the same cell, respectively; ${({\theta }_{\nu })}_{\mathrm{ref}}$ is the zenith angle for the reference geodesic curve. We evaluate this function for both the flat and Schwarzschild spacetimes.

Figure 7 shows the radial profiles of the rescaled error function, i.e., the error function multiplied by the number of angular mesh points ${N}_{{\theta }_{\nu }}$. The left and right panels show the results for the flat and Schwarzschild spacetimes, respectively. For both cases, the rescaled errors for the different numbers of mesh points are close to each other, implying that our code is of first order in the angular advection. This is just as expected, since we adopt the first-order finite-difference scheme for the angular advection terms.

Zoom In Zoom Out Reset image size

4.1.3. Tests of Energy and Angular Advection Combined

In this section, we conduct the advection tests, taking both the energy and angular advections into account, since they occur simultaneously in reality. First, we redo the same test as in Section 4.1.2 with the energy advection turned on. We deploy the same mesh with N_r = 128, N = 20, and ${N}_{{\theta }_{\nu }}=20$ and place the neutrino steady source at r = 20 km that emits neutrinos in a single angular bin closest to θ_ν = π/2 and a single energy bin at = 20 MeV. The angular advection error ${E}_{{\theta }_{\nu }}$ is calculated for the energy-integrated distribution function and compared with the previous result without energy advection (see Section 4.1.2). Figure 8 shows ${E}_{{\theta }_{\nu }}$ as a function of radius for both cases. The results are very close to each other, implying that the inclusion of the energy advection does not affect the angular advection.

Zoom In Zoom Out Reset image size

Next, we consider a neutrino source that has smooth, extended energy and angular distributions in order to demonstrate that the numerical diffusion behavior is reduced for smoother distributions. We set the following distribution at the source:

$$\begin{eqnarray}f(x)=\left\{\begin{array}{ll}\displaystyle \frac{1}{1+{e}^{(\epsilon -\mu )/{kT}}}\displaystyle \frac{1+\cos 2{\theta }_{\nu }}{2} & \left({\theta }_{\nu }\leqslant \displaystyle \frac{\pi }{2}\right)\\ 0 & \left({\theta }_{\nu }\gt \displaystyle \frac{\pi }{2}\right)\end{array}\right.,\end{eqnarray} \tag{ 11 }$$

where we choose the parameters as μ = 20 MeV and kT = 10 MeV. We place this steady source at r = 20 km. We employ the same radial mesh with N_r = 128 as in the previous tests, whereas we vary the number of grid points in the energy mesh that covers the range ∈ [0, 300] MeV as N = 10 and 20; the cell number in the angular mesh is chosen to be either ${N}_{{\theta }_{\nu }}=20$ or 40.

Figure 9 shows the angular-integrated energy distribution (top left panel) and angular distribution for the neutrino energy of = 5 MeV (top right panel) at r = 40 km. The dashed lines are the distributions at the source position, while the solid lines give the analytic solutions (green) and the numerical results (blue for the lower resolution and red for the higher resolution). The bottom panels present the absolute values of the errors. As can be seen, the energy spectrum is already well reproduced with 10 energy-grid points, which is actually smaller than the standard number employed in our recent CCSN simulations, and the numerical result gets even closer to the analytic solution with 20 energy-grid points. Such a converging feature is also seen for the angular distributions. By comparing the right panel of Figure 9 (result for smooth distribution) and Figure 6 (result for single angular bin test; Section 4.1.2), it is intuitively understandable that numerical diffusion is much smaller for the extended distribution assumed here.

Zoom In Zoom Out Reset image size

The results of the above test for the smooth distribution suggest that 20 energy bins employed in our CCSN simulations are large enough, whereas 10 angular bins in momentum space are not sufficient at large radii, where the angular distribution becomes forward-peaked as assumed in this test. This is actually a well-known problem and is consistent with the previous investigation by Richers et al. (2017). We note, however, that the number of these mesh points can be increased by a factor of 2 or more (depending on which number is increased) when the latest Japanese flagship supercomputer Fugaku is available soon, which is roughly ∼40 times faster than the K supercomputer, which we used for the SN simulations so far.

4.2. Axisymmetric Tests

Our code is multidimensional in space. Here we test our code's capability to deal with the angular advection in space by again calculating the nonradial streaming of neutrinos in the Schwarzschild spacetime with this θ advection explicitly taken into account.

We hence run the code in 2D under axisymmetry in this section. We compute the distribution function of neutrinos on the ϕ constant meridional plane. The initial condition is as follows: we set f = 1 for a single cell at r = 30 km and θ = π/2 with θ_ν = π/2 and ϕ_ν = π (i.e., moving in the positive z direction) in momentum space. Note that the neutrinos move along a geodesic curve in this case that is essentially the same as the one in the previous test. We again switch off the energy advection in this test. Also, the ϕ_ν advection is switched off, since neutrinos following the geodesic do not advect in that direction in the current setting.

The radial mesh is the same as in Section 4.1, with N_r = 128. In addition to this, we deploy the θ mesh that has N_θ = 128 bins, covering the range θ ∈ [0, π]. As for the angular mesh in momentum space, we adopt ${N}_{{\theta }_{\nu }}=10$, 20, and 40 to see the resolution dependence; the number of ϕ_ν mesh points fixed to ${N}_{{\phi }_{\nu }}=10$.

Figure 10 shows the neutrino number density, i.e., the distribution function integrated over the momentum space, on the meridional plane for ${N}_{{\theta }_{\nu }}=40$. The white dashed curve depicts the geodesic curve drawn from the source to the point that the massless particles reach at the time of the snapshot. The left panel is the result for the flat spacetime at the coordinate time t = 1 × 10⁻⁴ s, where the geodesic is the straight line parallel to the z-axis. The right panel is for the Schwarzschild spacetime at t = 1.3 × 10⁻⁴ s. It is apparent that the geodesic is deflected by gravity in this case. In both cases, the numerically obtained distributions are consistent with the analytical curves. On the other hand, the broadening of the beam is also evident. Just as in the previous tests, this is partly because the beam has a finite width from the beginning and partly because there are numerical diffusions.

Zoom In Zoom Out Reset image size

We now quantify the numerical diffusion. This time, we look at the neutrino propagation speed, which should be the speed of light, but the diffusion will affect it. For this purpose, we first evaluate the number densities on the geodesic as a function of θ by linearly interpolating the values on the neighboring radial cells as follows:

$$\begin{eqnarray}&&{N}_{\mathrm{int}}(\theta )=\displaystyle \frac{({r}_{n}-r(\theta )){N}_{n-1}(\theta )+(r(\theta )-{r}_{n-1}){N}_{n}(\theta )}{{r}_{n}-{r}_{n-1}},\end{eqnarray} \tag{ 12 }$$

where N_n (θ) is the neutrino number density on the nth radial mesh point for the given θ, and r(θ) is the radial coordinate of the geodesic curve at the same angle θ. We parameterize the geodesic not with θ but with the "light-traveling distance," defined as

$$\begin{eqnarray}&&\int \sqrt{{g}_{00}^{-1}({g}_{11}{{dr}}^{2}+{g}_{22}d{\theta }^{2})},\end{eqnarray} \tag{ 13 }$$

where g_{μ ν} is the spacetime metric and the integration runs from the source to a point on the geodesic curve. This quantity has a simple physical interpretation: the light speed times the coordinate time it takes a massless particle to reach the point.

Figure 11 shows the number density profile for different time steps along the geodesic obtained with Equation (12) and parameterized with the light-traveling distance in Equation (13). The upper and lower panels show the results for the flat and Schwarzschild spacetimes, respectively. The dashed lines indicate the exact results. Although the number density is not constant owing to the beam broadening, and there are some superluminal diffusions, the distribution declines rapidly ahead of the exact position, and we may hence say that the propagation velocity is roughly consistent with the speed of light.

Zoom In Zoom Out Reset image size

We finally study the resolution dependence, defining the following error function to quantitatively estimate the numerical diffusion:

$$\begin{eqnarray}&&{E}_{r\theta }(\theta )\equiv \displaystyle \frac{{\displaystyle \sum }_{m=1}^{{N}_{r}}{N}_{n}(\theta ){\left({r}_{n}-r(\theta )\right)}^{2}{{dr}}_{n}}{{\left(r(\theta )\right)}^{2}{\displaystyle \sum }_{n=1}^{{N}_{r}}N({r}_{n},\theta ){{dr}}_{n}},\end{eqnarray} \tag{ 14 }$$

where r_n and dr_n are the radial coordinate at the center and the width of the nth radial cell, respectively, and r(θ) is the radius of the point on the geodesic curve at θ.

Figure 12 shows the rescaled error function, which is defined as the error function in Equation (14) multiplied by ${N}_{{\theta }_{\nu }};$ although Equation (14) is a function of θ, we employ the light-traveling distance to parameterize the geodesic in this figure. The left and right panels again present the results for the flat and Schwarzschild spacetimes, respectively. It is recognized that both cases roughly show the first-order convergence.

Zoom In Zoom Out Reset image size

We now move on to the advection tests in the Kerr spacetime. They are intended to demonstrate the applicability of our code to rotating BH spacetimes. We again ignore all neutrino reactions in this section.

We employ the Kerr–Schild coordinates,

$$\begin{eqnarray}\begin{array}{rcl}{{ds}}^{2} & = & -\left(1-\displaystyle \frac{2{GMr}}{{c}^{2}{\rm{\Sigma }}}\right){c}^{2}{{dt}}^{2}\\ & & +\displaystyle \frac{4{GMr}}{{c}^{2}{\rm{\Sigma }}}\left[\displaystyle \frac{{\rm{\Sigma }}}{{r}^{2}+{a}^{2}}{dr}-a{\sin }^{2}\theta d\phi \right]{cdt}\\ & & +\left(\displaystyle \frac{{\rm{\Sigma }}}{{r}^{2}+{a}^{2}}+\displaystyle \frac{2{GMr}{\rm{\Sigma }}}{{c}^{2}{\left({r}^{2}+{a}^{2}\right)}^{2}}\right){{dr}}^{2}\\ & & +{\rm{\Sigma }}d{\theta }^{2}+\displaystyle \frac{{\rm{\Xi }}}{{\rm{\Sigma }}}{\sin }^{2}\theta d{\phi }^{2}-\displaystyle \frac{4{GMra}}{{c}^{2}({r}^{2}+{a}^{2})}{\sin }^{2}\theta {drd}\phi ,\end{array}\end{eqnarray} \tag{ 15 }$$

where a is the BH spin parameter, and

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Sigma }} & = & {r}^{2}+{a}^{2}{\cos }^{2}\theta ,\\ {\rm{\Xi }} & = & ({r}^{2}+{a}^{2}){\rm{\Sigma }}+2{{GMa}}^{2}r{\sin }^{2}\theta /{c}^{2}.\end{array}\end{eqnarray} \tag{ 16 }$$

Note that there is no (apparent) singularity at the event horizon in these coordinates. Throughout this section, we choose M = 5 M_⊙ and a = 0.5GM/c². The tetrad in Equation (3) is explicitly given as follows:

$$\begin{eqnarray}\begin{array}{rcl}{e}_{(0)}^{\mu } & = & \displaystyle \frac{1}{\alpha }\left[{\left(\displaystyle \frac{\partial }{\partial t}\right)}^{\mu }-{\beta }^{r}{\left(\displaystyle \frac{\partial }{\partial r}\right)}^{\mu }-{\beta }^{\phi }{\left(\displaystyle \frac{\partial }{\partial \phi }\right)}^{\mu }\right],\\ {e}_{(1)}^{\mu } & = & \displaystyle \frac{1}{\sqrt{{\gamma }_{{rr}}}}{\left(\displaystyle \frac{\partial }{\partial r}\right)}^{\mu },\\ {e}_{(2)}^{\mu } & = & \displaystyle \frac{1}{\sqrt{{\rm{\Sigma }}}}{\left(\displaystyle \frac{\partial }{\partial \theta }\right)}^{\mu },\\ {e}_{(3)}^{\mu } & = & \displaystyle \frac{{\gamma }^{r\phi }}{\sqrt{{\gamma }^{\phi \phi }}}{\left(\displaystyle \frac{\partial }{\partial r}\right)}^{\mu }+\displaystyle \frac{1}{\sqrt{{\gamma }^{\phi \phi }}}{\left(\displaystyle \frac{\partial }{\partial \phi }\right)}^{\mu },\end{array}\end{eqnarray} \tag{ 17 }$$

where the lapse function, shift vector, and inverse of the spatial metric are written as

$$\begin{eqnarray}\begin{array}{rcl}\alpha & = & \sqrt{\displaystyle \frac{{c}^{2}{\rm{\Sigma }}}{{c}^{2}{\rm{\Sigma }}+2{GMr}}},\\ {\beta }^{r} & = & \displaystyle \frac{2{GMr}}{{c}^{2}{\rm{\Sigma }}+2{GMr}},\\ {\beta }^{\phi } & = & -\displaystyle \frac{2{GMar}}{({c}^{2}{\rm{\Sigma }}+2{GMr})({r}^{2}+{a}^{2})},\\ {\gamma }^{{rr}} & = & \displaystyle \frac{{c}^{2}{\rm{\Xi }}}{{\rm{\Sigma }}({c}^{2}{\rm{\Sigma }}+2{GMr})},\\ {\gamma }^{r\phi } & = & \displaystyle \frac{2{GMar}}{({c}^{2}{\rm{\Sigma }}+2{GMr})({r}^{2}+{a}^{2})},\\ {\gamma }^{\phi \phi } & = & \displaystyle \frac{{\rm{\Sigma }}({r}^{2}+2{GMr}/{c}^{2}+{a}^{2})}{{\left({r}^{2}+{a}^{2}\right)}^{2}({\rm{\Sigma }}+2{GMr}/{c}^{2}){\sin }^{2}\theta },\end{array}\end{eqnarray} \tag{ 18 }$$

respectively. As mentioned in Shibata et al. (2014), there is no coordinate singularity if we use this tetrad in the Kerr–Schild coordinates.

As in Section 4.1.2, we perform the advection test by placing a point source in the Kerr spacetime. We treat the neutrino propagation only on the equatorial plane (θ = π/2), as our main concern here is the dragging of the inertial frame. We fix f = 1 for a single angular bin at the position of the point source and initially set f = 0 otherwise. The energy advection is turned off again in this case. We set the initial direction of neutrino momentum to ϕ_ν = 3π/2, i.e., the retrograde direction with respect to the BH spin, to maximize the frame-dragging effect. We further assume p^r = 0 to clearly distinguish the geodesic curves outside the photon sphere from those inside. There appears to be a bit of a complication then, because this does not correspond to θ_ν = π/2 for our choice of tetrad in the Kerr spacetime. That happens because the coordinate basis is not orthogonal. We need to find the value of θ_ν by solving the following equation:

$$\begin{eqnarray}&&{p}^{r}=\epsilon ({e}_{(0)}^{1}+{e}_{(1)}^{1}\cos {\theta }_{\nu }-{e}_{(3)}^{1}\sin {\theta }_{\nu })=0.\end{eqnarray} \tag{ 19 }$$

This yields, for example, θ_ν = 0.9668 for r = 28 km and θ_ν = 0.7695 for r = 22 km.

Throughout this test, the radial mesh is the same as in the previous tests, deploying N_r = 128 grid points. We vary the number of angular mesh points as ${N}_{{\theta }_{\nu }}=30$, 40, and 60 to see the resolution dependence.

Figure 13 shows the angular distribution in momentum space as a function of the radius r for ${N}_{{\theta }_{\nu }}=60$. The white dashed curves are the reference geodesic curves calculated in Appendix C, drawn from the source to the radius that the massless particles reach at the time of the snapshot. The left panel shows the result at the coordinate time t = 4 × 10⁻⁴ s with the source located at r = 28 km. In this case, the source is sitting outside the photon sphere (r = 26 km), and the neutrinos propagate radially outward. The right panel shows the result at t = 1.5 × 10⁻⁴ s with the source located at r = 22 km, i.e., inside the photon sphere, and the geodesic goes radially inward. The broadening of the beam is again apparent, which is inevitable, as the mesh size is finite, and in the latter case in particular, there are a fraction of neutrinos going radially outward due to the numerical diffusion. Nevertheless, most of the neutrinos propagate consistently with the geodesic curves in both cases. We may hence claim that our code can also handle the advection in the rotating spacetime.

Zoom In Zoom Out Reset image size

We look at the neutrino propagation speed. Since our computation does not take the ϕ advection into account directly thanks to the axisymmetry, we inspect the radial propagation. The geodesic is again parameterized by the light-traveling distance defined in the current case as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \int {\left({g}_{00}\right)}^{-1}\left[-({g}_{01}{dr}+{g}_{03}d\phi )\right.\\ \quad \left.-\sqrt{{\left({g}_{01}{dr}+{g}_{03}d\phi \right)}^{2}-{g}_{00}({g}_{11}{{dr}}^{2}+{g}_{33}d{\phi }^{2}+2{g}_{13}{drd}\phi )}\,\right].\end{array}\,\end{eqnarray} \tag{ 20 }$$

It has the same physical interpretation as in the Schwarzschild spacetime (Equation (13)). Figure 14 shows the profiles of the number density along the geodesic as a function of the light-traveling distance. Although the results are consistent with the propagation at the speed of light, the numerical diffusion is more remarkable than in the flat or Schwarzschild spacetime. This is due to the slower radial propagation of neutrinos in the current case, which is in turn caused by the frame-dragging in the rotating spacetime. Recall that we are considering the retrograde advection here.

Zoom In Zoom Out Reset image size

Using the error function given in Equation (10), we show in Figure 15 the rescaled error function for the advection test in the Kerr spacetime. It is apparent that the error scales linearly with the number of angular mesh points, just as in the previous tests in the flat or Schwarzschild spacetime.

Zoom In Zoom Out Reset image size

We have so far conducted various verification tests in the free-streaming limit alone. This is because the computation of advection is the most numerically challenging for mesh-based codes like ours. In the actual astrophysical simulations, the collision terms of the Boltzmann equation are also important. In this section, we perform a test that takes into account interactions of neutrinos with matter.

We fix the background matter distribution and compute with our new code the neutrino distribution that approaches a steady state. Since the collision terms are the focus of this test, we again assume spherical symmetry. As a reference, we employ the results obtained with another spherically symmetric, general relativistic Boltzmann solver (hereafter the 1D GR code) developed by Sumiyoshi et al. (2005). We run these two Boltzmann solvers for the same (fixed) matter distribution until the steady states are achieved. As for the matter distribution, we pick up a snapshot at 100 ms after bounce in the CCSN simulation by Sumiyoshi et al. (2005) with the 15 M_⊙ model of Woosley & Weaver (1995). It allows us to test the neutrino propagation in a wide range of mean free paths. The major neutrino–matter reactions are the same for both simulations, the rates of which are based on Bruenn (1985). The explicit expressions of the collision terms are identical to those described in Sumiyoshi & Yamada (2012). The chemical potential profile is also common because it affects the reaction rates; it is calculated from Shen's equation of state (Shen et al. 1998). In the 1D GR code, the hydrodynamics variables, as well as the spacetime metric, are the functions of the enclosed mass, as it is a Lagrangian code. However, the fluid velocity is fixed to zero. They are transformed into the functions of the radius when they are implemented in the new code.

The numbers of the radial mesh points are the same for the two codes: the radial mesh has N_r = 256 grid points, which covers the range r ∈ [0, 10⁴] km. As for the momentum space, we adopt ${N}_{{\theta }_{\nu }}=6$ and N = 14. The latter mesh covers ∈ [0, 300] MeV. The reader may be worried that the resolution employed here is relatively lower compared to the advection tests in Sections 4 and 5 or recent core-collapse simulations. This is no problem, since our purpose is the regression test, i.e., to confirm that the collision terms work properly.

Figure 16 shows the radial profiles of the electron-type neutrino number density for our new code and the 1D GR code (upper panel), as well as the relative difference with respect to the latter result (lower panel). The shock wave is located at r = 174 km, at which a small bump can be seen. It is apparent that they are in a reasonable agreement with the overall relative error of ∼10% (∼14% at the maximum). It is particularly small, ∼1%, near the center (r ≲ 10 km), where the neutrino distribution is close to that of thermal equilibrium.

Zoom In Zoom Out Reset image size

It should be noted that the above comparison does not tell which code is more accurate. To get some hints, we compare them with the equilibrium number densities, which are obtained by the energy integration of the Fermi–Dirac distribution. Figure 17 shows the relative deviation of the neutrino number density obtained with each of the two codes from the equilibrium density. As is clear, the difference is much smaller for the new code in the very optically thick regime. It is thought to come from the different way of evaluating the equilibrium distribution for the calculation of the reaction kernels. The 1D GR code simply calculates the equilibrium distribution by the value of energy at the center of the cell, whereas the new code derives the equilibrium value as the average over an energy bin by dividing it into subgrids. Incidentally, neither result is very accurate at large radii because of the rather coarse angular resolutions employed.

Zoom In Zoom Out Reset image size

Finally, we give some results of the comparison with the M1 closure, the currently most popular approximation in the literature. In the neutrino transport with the M1 closure, the second angular moment P^ij of the distribution function in momentum space is given as the interpolation of the optically thin and thick limits,

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

where the optically thin and thick limits are given as

$$\begin{eqnarray}&&{P}_{\mathrm{thin}}^{{ij}}(\epsilon )=E(\epsilon )\displaystyle \frac{{F}^{i}(\epsilon ){F}^{j}(\epsilon )}{{\gamma }_{{kl}}{F}^{k}(\epsilon ){F}^{l}(\epsilon )},\quad {P}_{\mathrm{thick}}^{{ij}}(\epsilon )=\displaystyle \frac{{\delta }^{{ij}}}{3},\end{eqnarray} \tag{ 22 }$$

respectively, where E() and Fⁱ () are the energy density and flux, respectively. Note that we only treat the case where fluid velocity is zero throughout this discussion. The Eddington factor χ is given, for example, as (Levermore 1984)

$$\begin{eqnarray}&&\chi =\displaystyle \frac{3+4{\bar{F}}^{2}}{5+2\sqrt{4-3{\bar{F}}^{2}}},\end{eqnarray} \tag{ 23 }$$

where the flux factor $\bar{F}$ may be calculated as (Shibata et al. 2011)

$$\begin{eqnarray}&&\bar{F}=\displaystyle \frac{\sqrt{{\gamma }_{{ij}}{F}^{i}(\epsilon ){F}^{j}(\epsilon )}}{E(\epsilon )}.\end{eqnarray} \tag{ 24 }$$

In the Boltzmann transport, the second moment can be directly calculated from the distribution function, and in spherical symmetry, the Eddington factor is χ = k^rr , where k^ij () = P^ij ()/E() is the Eddington tensor.

Figure 18 shows the radial profiles of the Eddington factors obtained directly by the Boltzmann solver (solid lines) and approximately with the M1 closure relation (dashed lines). It is observed that the M1 closure tends to give larger values than the Boltzmann solver; the latter, with lower angular resolutions, tends to underestimate the Eddington factor, particularly at large radii, where the neutrino angular distribution becomes forward-peaked. The M1 closure fails, on the other hand, in the semitransparent region ∼140 km. All of these results are quantitatively consistent with those found in our previous studies in the Minkowski spacetime (Harada et al. 2019). It is worth mentioning, however, that this is not a true comparison between the Boltzmann transport and the two-moment transport with the M1 closure. In this analysis of the M1 closure approximation, the Eddington factor is calculated from the energy density and the flux obtained by the Boltzmann solver. In actual two-moment transport, they should be computed on their own, and errors may be accumulated with time and could be larger than suggested in Figure 18. In the meantime, the 10 grid points in θ_ν are not sufficient in our Boltzmann solver and should be doubled or more, which is implied by this test.

Zoom In Zoom Out Reset image size

We have developed a numerical code to solve the multidimensional Boltzmann equation for neutrino transfer in full general relativity, employing the discrete ordinate method. We have performed a series of tests to confirm that our code can handle the general relativistic neutrino transport. We first conducted the tests of free streaming in the Schwarzschild and Kerr spacetimes. We computed the propagation of neutrinos from a source located either inside or outside the photon sphere in these spacetimes. We confirmed that the results are consistent with the reference geodesic curves. We also found that the deviation from the reference is inversely proportional to the number of mesh points either in energy or zenith angle, implying that the numerical accuracy of the energy or angular advection is of the first order, as expected from our finite-difference schemes adopted for these terms.

Next, we ran the code with the collision terms included for a fixed matter distribution extracted from our realistic 1D simulation of CCSNe. We demonstrated that the result is consistent with that obtained with another 1D GR code developed by Sumiyoshi et al. (2005). We also found that our new code can reproduce the thermal equilibrium distribution in the optically thick region more accurately than the 1D GR code.

In this paper, we ignore the evolution of matter configurations and spacetime entirely. Having completed the development of the Boltzmann solver, our next task is to combine it with the general relativistic hydrodynamics code and the solver of the Einstein equations, both of which have already been developed separately. The latter, in particular, is based on the Baumgarte–Shapiro–Shibata–Nakamura formalism (Shibata & Nakamura 1995; Baumgarte & Shapiro 1998) formulated for the polar coordinates, employing the partially implicit Runge–Kutta method (Baumgarte et al. 2013). Their performance tests will be reported elsewhere in the near future. Once completed, we will apply the integrated code to various general relativistic phenomena.

Numerical computations were carried out in part on the Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan, and the supercomputer SX-Aurora of the Particle, Nuclear and Astrophysics Simulation Program (No. 2019-002, 2020-004) at the High Energy Accelerator Research Organization (KEK). This work is supported by Grants-in-Aid for Scientific Research (19K03837, 20H01905, 19H05811) and Research Activity Start-up (19K23435) from the Japan Society for the Promotion of Science (JSPS) and Grants-in-Aid for Scientific Research on Innovative areas "Gravitational wave physics and astronomy: Genesis" (17H06357, 17H06365) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. S.Y. is supported by the Institute for Advanced Theoretical and Experimental Physics, Waseda University, and the Waseda University Grant for Special Research Projects (project No. 2020-C273). For providing high-performance computing resources, the Computing Research Center, KEK, JLDG on SINET of NII, Research Center for Nuclear Physics, Osaka University, Yukawa Institute of Theoretical Physics, Kyoto University, Nagoya University, Hokkaido University, and Information Technology Center, University of Tokyo, are acknowledged. This work was supported by MEXT as "Program for Promoting Researches on the Supercomputer Fugaku" (Toward a unified view of the universe: from large scale structures to planets). This research used the high-performance computing resources of the K-computer/the supercomputer Fugaku provided by RIKEN, the FX10 provided by Tokyo University, the FX100 provided by Nagoya University, the Grand Chariot provided by Hokkaido University, and Oakforest-PACS provided by JCAHPC through the HPCI System Research Project (Project ID: hp130025, 140211, 150225, 150262, 160071, 160211, 170031, 170230, 170304, 180111, 180179, 180239, 190100, 190160, 200102, 200124).

In this appendix, we describe how the advection terms are evaluated in our finite-difference scheme. It is the general relativistic extension from the Newtonian counterpart given in Sumiyoshi & Yamada (2012). Throughout this section, variables with the lowercase subscript n (such as r_n ) are defined at the cell centers, whereas those with the uppercase subscript N (such as r_N ) are defined at the cell interfaces.

The evaluation of the advection terms in the discretized form is tricky when they are written in the conservative form. This is because there are some delicate cancellations among the terms that are a part of different advection terms. This may be understood if one considers a constant f in the phase space. Then all advections should be vanishing. This is apparent in the Boltzmann equation written as in Equation (1), since all derivatives of f are zero. This is not so obvious, however, if the Boltzmann equation is cast into the conservative form, as in Equation (2). In fact, the advection terms are not simply vanishing but rather are canceled among them in nontrivial ways. Analytically, these cancellations are no problem; it is difficult to enforce in the finite difference, however. If one treats them carelessly, neutrinos may appear from nowhere. The formulation described in Sumiyoshi & Yamada (2012) ensures the perfect cancellation of those terms for the flat spacetime. In the general relativistic case with an arbitrary metric, the formulation for the flat spacetime is not applicable, though. In this paper, we do not stick to the perfect cancellation, and the advection terms are evaluated in a rather straightforward way, as described below. As shown in the advection tests in Sections 4 and 5, however, we found no problem in the general relativistic cases even if the cancellation is not exactly enforced. We will investigate this issue further and, if necessary, revise the general relativistic formulation in future works.

The numerical grid employed in this paper is based on the original flat version of the code, and the details are described in Sumiyoshi & Yamada (2012). The original code employed $\mu \equiv \cos \theta$ and ${\mu }_{\nu }\equiv \cos {\theta }_{\nu }$ instead of θ and θ_ν in order to simplify the advection terms, and the six coordinate variables (r, μ, ϕ, , μ_ν , and ϕ_ν ) were discretized as described in the paper. In our general relativistic extension, on the other hand, the use of μ is not so useful, though. We hence use θ as an independent variable but give its values at the cell centers and interfaces according to ${\theta }_{i}={\cos }^{-1}({\mu }_{i})$ and ${\theta }_{I}={\cos }^{-1}({\mu }_{I})$ from the corresponding values of μ.

The spatial advection terms (r, θ, and ϕ) are finite-differenced as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{1}{\sqrt{-g}}\displaystyle \frac{\partial }{\partial {x}^{i}}\left[{L}^{i}\sqrt{-g}f\right]\\ \quad \longrightarrow \displaystyle \frac{1}{\sqrt{-{g}_{n}}{{dx}}_{n}^{i}}\left[\sqrt{-{g}_{N}}{L}_{N}^{i}{f}_{N}-\sqrt{-{g}_{N-1}}{L}_{N-1}^{i}{f}_{N-1}\right],\end{array}\end{eqnarray} \tag{ A1 }$$

where dxⁱ _n is the width of the nth cell, and we define ${L}^{\mu }\equiv {e}_{(0)}^{\mu }+{\sum }_{i=1}^{3}{l}_{i}{e}_{i}^{\mu }$ for notational simplicity. In the flat or Schwarzschild case, for example, ${L}^{1}=\cos {\theta }_{\nu }$, ${L}^{2}\,=\sin {\theta }_{\nu }\cos {\phi }_{\nu }$, ${L}^{3}=\sin {\theta }_{\nu }\sin {\phi }_{\nu }$. The value of L_N f_N on the Nth cell interface (${L}_{N}^{i}{f}_{N}$) is evaluated as

$$\begin{eqnarray}\begin{array}{rcl}{L}_{N}^{i}{f}_{N} & = & \displaystyle \frac{{L}_{n}^{i}-| {L}_{n}^{i}| }{2}[(1-{\beta }_{N}){f}_{n}+{\beta }_{N}{f}_{n+1}]\\ & & +\displaystyle \frac{{L}_{n+1}^{i}+| {L}_{n+1}^{i}| }{2}[{\beta }_{N}{f}_{n}+(1-{\beta }_{N}){f}_{n+1}],\end{array}\end{eqnarray} \tag{ A2 }$$

where β_N is introduced to smoothly switch from the upwind differencing in the free-streaming limit (β_N = 1) to the central differencing in the diffusion limit (β_N = 1/2). We use the following expression of β_N based on Mezzacappa & Bruenn (1993):

$$\begin{eqnarray}&&{\beta }_{N}=1-\displaystyle \frac{1}{2}\displaystyle \frac{{\alpha }_{p}{{dr}}_{N}{\lambda }_{N}}{1+{\alpha }_{p}{{dr}}_{N}{\lambda }_{N}},\end{eqnarray} \tag{ A3 }$$

where the mean free path on the cell interface is calculated as λ_N = (λ_n+1 + λ_n )/2. The adjustable parameter α_p is set to be 100, following Sumiyoshi & Yamada (2012). In the free-streaming case β_N = 1, the upwind direction is accounted for by the signature of Lⁱ .

The energy advection term is expressed as follows:

$$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{\partial }{\partial ({\epsilon }^{3}/3)}\left({\epsilon }^{3}f{\omega }_{(0)}\right)\\ \quad \longrightarrow -\displaystyle \frac{3}{{\epsilon }_{N}^{3}-{\epsilon }_{N-1}^{3}}\left[{\epsilon }_{N}^{3}{\omega }_{(0)N}{f}_{N}-{\epsilon }_{N-1}^{3}{\omega }_{(0)N-1}\,{f}_{N-1}\right],\end{array}\end{eqnarray} \tag{ A4 }$$

with

$$\begin{eqnarray}&&{\omega }_{(0)N}{f}_{N}=\displaystyle \frac{{\omega }_{(0)N}+| {\omega }_{(0)N}| }{2}{f}_{n+1}+\displaystyle \frac{{\omega }_{(0)N}-| {\omega }_{(0)N}| }{2}{f}_{n}.\end{eqnarray} \tag{ A5 }$$

Here the one-sided differencing is employed, and its direction is dictated by the signature of ω_(0)n ; the forward finite differencing is adopted for ω_(0)n > 0, i.e., in the case for redshift, whereas the backward differencing is employed for the blueshift case.

The angular advection terms are finite-differenced in a similar way, by switching the direction of the one-sided differencing according to the signatures of ω_(*); employing ${\mu }_{\nu }=\cos {\theta }_{\nu }$, we write the θ_ν advection term as

$$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{\partial }{\partial {\mu }_{\nu }}(\sin {\theta }_{\nu }f{\omega }_{({\theta }_{\nu })})\\ \quad \longrightarrow -\displaystyle \frac{1}{{\left(d{\mu }_{\nu }\right)}_{n}}\left[\sin {\left({\theta }_{\nu }\right)}_{N}{\omega }_{({\theta }_{\nu })N}{f}_{N}-\right.\\ \quad \left.\sin {\left({\theta }_{\nu }\right)}_{N-1}{\omega }_{({\theta }_{\nu })N-1}{f}_{N-1}\right],\end{array}\end{eqnarray} \tag{ A6 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{\omega }_{({\theta }_{\nu })N}{f}_{N} & = & \displaystyle \frac{{\omega }_{({\theta }_{\nu })N}+| {\omega }_{({\theta }_{\nu })N}| }{2}{f}_{n+1}\\ & & +\displaystyle \frac{{\omega }_{({\theta }_{\nu })N}-| {\omega }_{({\theta }_{\nu })N}| }{2}{f}_{n}.\end{array}\end{eqnarray} \tag{ A7 }$$

The ϕ_ν advection term is finite-differenced as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{1}{{\sin }^{2}{\theta }_{\nu }}\displaystyle \frac{\partial }{\partial {\phi }_{\nu }}(f{\omega }_{({\phi }_{\nu })})\\ \quad \longrightarrow \displaystyle \frac{1}{{\sin }^{2}{\left({\theta }_{\nu }\right)}_{{n}_{{\theta }_{\nu }}}}\displaystyle \frac{1}{{\left(d{\phi }_{\nu }\right)}_{{n}_{{\phi }_{\nu }}}}\\ \quad \times \left[{\omega }_{({\phi }_{\nu }){N}_{{\phi }_{\nu }}}{f}_{{N}_{{\phi }_{\nu }}}-{\omega }_{({\phi }_{\nu }){N}_{{\phi }_{\nu }}-1}{f}_{{N}_{{\phi }_{\nu }}-1}\right],\end{array}\end{eqnarray} \tag{ A8 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{\omega }_{({\phi }_{\nu })N}{f}_{N} & = & \displaystyle \frac{{\omega }_{({\phi }_{\nu })N}-| {\omega }_{({\phi }_{\nu })N}| }{2}{f}_{n+1}\\ & & +\displaystyle \frac{{\omega }_{({\phi }_{\nu })N}+| {\omega }_{({\phi }_{\nu })N}| }{2}{f}_{n}.\end{array}\end{eqnarray} \tag{ A9 }$$

We summarize how we calculate the geodesic curves in the Schwarzschild spacetime, employed as the reference in Section 4. In the Schwarzschild (exterior) spacetime, the geodesic motion on the equatorial plane satisfies the following equation (see Section 25.6 in Misner 1973, for example):

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{dr}}{d\phi }\right)}^{2}={r}^{2}\left(\displaystyle \frac{{r}^{2}}{{b}^{2}}+\displaystyle \frac{2{GM}}{{c}^{2}r}-1\right),\end{eqnarray} \tag{ B1 }$$

where r and ϕ are the coordinate variables, and b is the impact parameter. When the trajectory is on the meridional plane, one may replace the azimuth angle ϕ with the zenith angle θ. The impact parameter b can be expressed in terms of the radius r₀ and zenith angle ${({\theta }_{\nu })}_{0}$ of a reference point on the geodesic curve as

$$\begin{eqnarray}&&b=\displaystyle \frac{{r}_{0}\sin {\left({\theta }_{\nu }\right)}_{0}}{1-2{GM}/({c}^{2}{r}_{0})}.\end{eqnarray} \tag{ B2 }$$

For numerical calculations, it is more useful to rewrite the above Equation (B1) as

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}u}{d{\phi }^{2}}=\displaystyle \frac{3{{GMu}}^{2}}{{c}^{2}}-u,\end{eqnarray} \tag{ B3 }$$

where the new variable is introduced as u = 1/r. The second derivative of u changes signature at the radius r = 3GM/c², which corresponds to the radius of the photon sphere. In this work, Equation (B3) is solved with the fourth-order explicit Runge–Kutta method by dividing it into the following two equations:

$$\begin{eqnarray}&&\displaystyle \frac{{du}}{d\phi }=v,\quad \displaystyle \frac{{dv}}{d\phi }=\displaystyle \frac{3{{GMu}}^{2}}{{c}^{2}}-u.\end{eqnarray} \tag{ B4 }$$

For the angular advection tests in Section 4.1.2, we need the zenith angle θ_ν in momentum space along the geodesic, which is given in Equation (13) in Shibata et al. (2014),

$$\begin{eqnarray}&&{\theta }_{\nu }={\tan }^{-1}\left(\displaystyle \frac{{p}_{(2)}}{{p}_{(1)}}\right)={\tan }^{-1}\left(r\sqrt{1-\displaystyle \frac{2{GM}}{{c}^{2}r}}\displaystyle \frac{{p}^{\theta }}{{p}^{r}}\right),\end{eqnarray} \tag{ B5 }$$

where p_(i) is the momentum component for the tetrad basis and given as ${p}_{(i)}={p}_{\mu }{e}_{(i)}^{\mu }$.

For the Kerr spacetime, there is no simple differential equation to describe the geodesic curve, unlike for the Schwarzschild spacetime. We hence solve the geodesic equation,

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{x}^{\mu }}{d{\lambda }^{2}}=-{{\rm{\Gamma }}}_{\alpha \beta }^{\mu }\displaystyle \frac{{{dx}}^{\alpha }}{d\lambda }\displaystyle \frac{{{dx}}^{\beta }}{d\lambda },\end{eqnarray} \tag{ C1 }$$

where λ is the affine parameter. Since we consider the geodesic curves only on the equatorial plane, we solve the equations only for t, r, and ϕ. We employ the fourth-order explicit Runge–Kutta for the following forms of the equations:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{dt}}{d\lambda } & = & {p}^{t},\qquad \displaystyle \frac{{{dp}}^{t}}{d\lambda }=-{{\rm{\Gamma }}}_{\alpha \beta }^{t}{p}^{\alpha }{p}^{\beta },\\ \displaystyle \frac{{dr}}{d\lambda } & = & {p}^{r},\qquad \displaystyle \frac{{{dp}}^{r}}{d\lambda }=-{{\rm{\Gamma }}}_{\alpha \beta }^{r}{p}^{\alpha }{p}^{\beta },\\ \displaystyle \frac{d\phi }{d\lambda } & = & {p}^{\phi },\qquad \displaystyle \frac{{{dp}}^{\phi }}{d\lambda }=-{{\rm{\Gamma }}}_{\alpha \beta }^{\phi }{p}^{\alpha }{p}^{\beta }.\end{array}\end{eqnarray} \tag{ C2 }$$

Here p^μ is the four-momentum.

In the Kerr spacetime, the photon sphere for the prograde orbits with respect to the BH spin is given as

$$\begin{eqnarray}&&{r}_{\mathrm{prog}}=\displaystyle \frac{2{GM}}{{c}^{2}}\left\{1+\cos \left[\displaystyle \frac{2}{3}{\cos }^{-1}\left(\displaystyle \frac{{c}^{2}| a| }{{GM}}\right)\right]\right\},\end{eqnarray} \tag{ C3 }$$

whereas for the retrograde orbits, it is given as

$$\begin{eqnarray}&&{r}_{\mathrm{ret}}=\displaystyle \frac{2{GM}}{{c}^{2}}\left\{1+\cos \left[\displaystyle \frac{2}{3}{\cos }^{-1}\left(-\displaystyle \frac{{c}^{2}| a| }{{GM}}\right)\right]\right\}.\end{eqnarray} \tag{ C4 }$$

For the metric parameters, we employ the same values as the advection tests in Section 5, with M = 5 M_⊙ and a = 0.5GM/c².

Figure 19 shows some geodesic curves on the equatorial plane. The left panel exhibits four prograde geodesic curves with p^r = 0 at r = 16, 17, 18, and 19 km. The first two radii are smaller than that of the photon sphere, which is r = 17.33 km in the present case, whereas the last two are larger. The photon sphere is indicated with the black circle in the figure. The right panel, on the other hand, presents some retrograde geodesic curves with p^r = 0 at r = 25, 26, 27, and 28 km. The photon sphere has a radius of r = 26.07 km in this case. It is found in both panels that the geodesic curves are confined either inside or outside the photon sphere.

Zoom In Zoom Out Reset image size

For the advection tests in the Kerr spacetime in Section 5, we need the zenith angle θ_ν of the four-momentum on the geodesic curve. Just as for the Schwarzschild spacetime, it is obtained from

$$\begin{eqnarray}&&{\theta }_{\nu }={\tan }^{-1}\left(\displaystyle \frac{{p}_{(3)}}{{p}_{(1)}}\right)={\tan }^{-1}\left(\displaystyle \frac{{p}_{r}{e}_{(3)}^{r}+{p}_{\phi }{e}_{(3)}^{\phi }}{{p}_{r}{e}_{(1)}^{r}}\right).\end{eqnarray} \tag{ C5 }$$