ThreeDimensional BoltzmannHydro Code for corecollapse in massive stars I. Special Relativistic Treatments
Abstract
We propose a novel numerical method for solving multidimensional, special relativistic Boltzmann equations for neutrinos coupled to hydrodynamics equations. It is meant to be applied to simulations of corecollapse supernovae. We handle special relativity in a nonconventional way, taking account of all orders of . Consistent treatment of advection and collision terms in the Boltzmann equations is the source of difficulties, which we overcome by employing two different energy grids: Lagrangian remapped and laboratory fixed grids. We conduct a series of basic tests and perform a onedimensional simulation of corecollapse, bounce and shockstall for a progenitor model with a minimum but essential set of microphysics. We demonstrate in the latter simulation that our new code is capable of handling all phases in corecollapse supernova. For comparison, a nonrelativistic simulation is also conducted with the same code, and we show that they produce qualitatively wrong results in neutrino transfer. Finally, we discuss a possible incorporation of general relativistic effects in our method.
Subject headings:
supernovae: general—neutrinos—hydrodynamicsYukawa Institute for Theoretical Physics, Kyoto University, Oiwakecho, Kitashirakawa, Sakyoku, Kyoto, 6068502, Japan
Numazu College of Technology, Ooka 3600, Numazu, Shizuoka 4108501, Japan
Advanced Research Institute for Science & Engineering, Waseda University, 341 Okubo, Shinjuku, Tokyo 1698555, Japan \addressDepartment of Science and Engineering, Waseda University, 341 Okubo, Shinjuku, Tokyo 1698555, Japan
1. Introduction
Quantitative studies on the mechanism of corecollapse supernovae (CCSNe) require detailed numerical simulations. Except for lowmass () progenitors, elaborate onedimensional (1D) simulations under spherical symmetry have not reproduced the supernova explosion (Sumiyoshi et al., 2005; Liebendörfer et al., 2005; Kitaura et al., 2006; Burrows et al., 2007). Last decade, most of supernova modelers have focused on the multidimensional (MultiD) aspects of dynamics (see e.g., Kotake et al. (2012a); Janka (2012); Burrows (2013) for recent review). In the postbounce phase, instabilitities drive postshock accretion flows into turbulence, making dynamics intrinsically multiD. This may be crucial for the supernova explosion, since the nonspherical turbulent motions increase the dwell time of material in the gain region, enhancing its absorption of hot neutrinos, boosting the post shock pressure, and eventually pushing the shock wave outwards (Takiwaki et al., 2012; Dolence et al., 2013).
As a matter of fact, we have recently witnessed shock revival in some of the currently most advanced simulations (Burrows et al., 2006; Marek & Janka, 2009; Suwa et al., 2010; Lentz et al., 2012; Müller et al., 2012a, b; Takiwaki et al., 2013), which has raised our hope that we will finally unveil the mechanism of CCSNe. Unfortunately, however, success or failure of the supernova explosion is a delicate problem. In fact, the latest results of MultiD simulations by different groups are still at odds with one another and no consensus has yet emerged concerning which ingredient(s) is (are) essential for explosion. Although various approaches, both phenomenological and ab initio, are being undertaken at present, only better simulations possibly with a Boltzmannequation solver that incorporate detailed microphysics and general relativity (GR) may give the conclusive answer.
Towards this goal, we are developing a numerical code for neutrino transfer, which solves the Boltzmann equations (Sumiyoshi & Yamada, 2012). Our code is based on the discreteordinate method, which finitedifferences the Boltzmann equations, deploying multiangle and multienergy bins in momentum space. Using some snapshots from threedimensional (3D) supernova simulations, Sumiyoshi & Yamada (2012) demonstrated the capabilities of this new code, which implements the minimum set of neutrino reactions (see also Sumiyoshi et al. (2014)). These simulations concerned neutrino transfer in static backgrounds, however, and no backreactions to matter were taken into account.
The next step should be a coupling of this code with a hydrodynamical code. This may not be so simple, though. Spherically symmetric 1D computations may be easier, since they can adopt Lagrangian formulations both for neutrino transfer and hydrodynamics (Mezzacappa & Bruenn, 1993; Mezzacappa et al., 2001; Liebendörfer et al., 2005; Sumiyoshi et al., 2005, 2007). Such formalisms as they are could not be applied in MultiD, however, and different formulations should be developed for the MultiD BoltzmannHydro simulations, i.e. the simulations that solve the Boltzmann equations and hydrodynamical equations simultaneously in multidimensions.
Unlike the previous 1D codes, we adopt an Eulerian picture in this paper. There are several reasons for this choice. Among other things, we have in mind that the Boltzmann solver will be coupled with a MultiD Eulerian hydrodynamics and gravity solvers, which have been well established and widely used in the highenergy astrophysical community. In addition, the Eulerian picture has a benefit to easily handle the left hand side of Boltzmann equation, i.e., advection terms. In general, Lagrangian formulations need to treat derivatives with respect to neutrino energy, which correspond to the Doppler’s effect caused by spatial and/or temporal variations in fluid velocities. This may cause problems particularly at a shock wave, where fluid velocities are discontinuous. For these reasons, we have opted the Eulerian approach.
It should be noted, however, that the Eulerian approach has its own demerits. For example, we normally need to handle transformations between the laboratory frame and the fluidrest frame defined locally, which are nothing but Lorentz transformations for flat spacetime, since neutrinomatter interactions are best described in the fluidrest frame. Physically consistent treatments of both advection and collision terms in the MultiD Boltzmann solver are technically difficult particularly in the method, since special relativistic (SR) effects such as Doppler shifts and aberrations should be handled on a rather coarse grid in the momentum space (see Section 3 for more details). Previous studies have attempted to alleviate this difficulty by employing an expansion of basic equations up to (see e.g., Hubeny & Burrows (2007)). In CCSNe, the maximum fluid velocity is around of the speed of light and such a firstorder approximation may be justified. The resultant equations are fairly complex and not easy to treat numerically, however, and the formulation is certainly not applicable to highly relativistic phenomena.
It should be also mentioned that several groups (Cardall et al., 2005, 2013; Peres et al., 2013) are developing different formulations for MultiD BoltzmannHydro simulations, which are yet to be implemented. Ott et al. (2008) performed detailed 2D BoltzmannHydro supernova simulations in the post bounce phase but they ignored SR effects. As shown later, consistent treatments of SR effects are indispensable to obtain correct behaviors in neutrino transfer. They are also the first step towards fully GR Boltzmann simulations, which will be needed to study more extreme phenomena such as black hole formations.
In this paper, we propose a novel formulation for the numerical computations of MultiD SR Boltzmann transfer based on the method, which treats SR effects to all orders of , where and denote the speed of light and fluid velocity, respectively. The accuracy of our new method is checked by the standard tests as well as by a realistic simulation of spherical collapse of a progenitor. As explained in the next section, the appropriate treatment of SR effects is crucial for numerically capturing the neutrinotrapping and the subsequent evolution up to bounce. In this paper, we particularly focus on this issue, and more detailed quantitative analyses of realistic supernova simulations by our BoltzmannHydro code will be reported later separately.
This paper is organized as follows. To facilitate readers’ understanding, we first give intuitive arguments on the importance of SR effect from the perspective of phase space (in Section 2), which will make clear why nonrelativistic BoltzmannHydro simulations fail to capture neutrino advections with matter and yield qualitatively wrong distributions of neutrinos. In Section 3, it is also emphasized that the treatment of SR effects is not so easy practically, and it is explained what is the main obstacle. Then the basic equations and formulations are presented in Section 4. After introducing two independent energygrids (which are essential for our SR treatment) in Section 5, the numerical algorithms are given in Section 6. We examine the accuracy of our new method by a series of SR Boltzmann and BoltzmannHydro simulations in Section 7. Finally we conclude the paper with a summary and discussions on further extensions of our code to a GR version in Section 8. Throughout this paper, Greek and Latin subscripts denote spacetime and space components, respectively. We use the metric signature of . Unless otherwise stated, we work in units with , where is the gravitational constant.
2. SR effects and neutrino trapping
Before going to details of our SR Boltzmann formulation and its numerical algorithm, we first give intuitive arguments on the importance of SR effects. As will be observed below, the ignorance of SR effects yields qualitatively wrong behaviors in neutrino distributions. The key ingredient is the angular distribution in phase space: isotropic distributions in the fluidrest frame become anisotropic in the laboratory frame after Lorentz transformations, the fact that ensures advection of neutrinos with matter and eventually the neutrino trapping.
In the following, we will explain this in a simplified and idealized setup. We consider neutrino transfer in moving matter that has uniform velocity and thermodynamic quantities. We assume in addition that neutrinos and matter are strongly coupled with each other via scattering and, as a result, the neutrino distribution function in phase space, , is isotropic in the fluidrest frame. This is a situation similar to the ones we see locally in the neutrino trapping phase. Then the neutrino flux at each point vanishes in the fluidrest frame and there is no net flux traversing fluid elements (see the left panel in Figure 1). The neutrino number in each fluid element is conserved as the fluid element moves at the finite velocity. This is the advection of neutrinos with matter and it should be evident why the Lagrangian approach is advantageous in dealing with it.
For comparison, the right panel in Figure 1 describes the same situation except in the laboratory frame. Here, we assume that the fluid is advected inwards (or leftwards in the figure). Since the neutrinos should be advected in the same direction as the fluid in the laboratory frame, the incoming neutrino flux is larger than the outgoing one, which means that the angular distribution of neutrinos is anisotropic in this frame. From the SR point of view, such anisotropies arise from the Doppler shift and relativistic beaming by Lorentz transformations. The mathematical expression of SR Boltzmann equations will be given in Section 4.
If we neglected all SR effects, not distinguishing the laboratory and fluidrest frames, we would not obtain the neutrino advection with matter, which is crucial for the neutrino trapping in the collapsing phase. In fact, neutrinos would be left behind as fluids are advected. The supernova core is not homogeneous in reality and both matter and neutrino densities are highest at the center. In the absence of advection, neutrinos would always flow outwards when they actually should move inwards, keeping pace with matter, and be effectively trapped in the core. As we will show later in Section 7.5, the number density of electrontype neutrinos becomes significantly smaller near the center for nonrelativistic simulations. It affects in turn the evolution of electron fraction and the size of inner core and eventually all the supernova dynamics thereafter.
3. Difficulties in handling SR effects
We give intuitive explanations more in detail on why SR treatments are not easy in the method, which we employ in this paper. The main source of the difficulties is scatterings, particularly those between neutrinos and nucleons (and nuclei). Other reactions such as neutrino absorptions and emissions have no technical difficulties
As mentioned in the previous sections, our BoltzmannHydro code is based on the Eulerian picture, and we discretize 6D phase space in the laboratory frame, as shown in the the left panel in Figure 2. In this picture, spherical coordinates in momentum space are adopted with the azimuthal dimension being collapsed. The radial direction corresponds to neutrino energy. Although the picture is drawn that way, gridding in each dimension is not necessarily uniform.
We first consider the isoenergetic scattering under the condition that fluid is at rest and, as a consequence, the laboratory fame coincides with the fluidrest frame. When a neutrino undergoes the isoenergetic scattering, it changes its flight direction specified by two angles, preserving energy. In the discretized momentum space, the neutrino moves from a bin to another with the same radialgrid number. The important thing is that only the angulargrid number is changed. In this case, there is no difficulty and, indeed, it has been implemented in Sumiyoshi & Yamada (2012); Sumiyoshi et al. (2014).
In the presence of nonvanishing fluid velocities, the problem becomes qualitatively different. In this case, the laboratory frame is different from the fluidrest frame and they are related with each other via a Lorentz transformation. The point is that the Lorentz transformation induces changes in both energy and angles. These energy shift and aberration are determined by the Doppler factor, which depends on the fluid velocity and neutrino angles (see Section 4). This is most clearly demonstrated in the right panel of Figure 2, in which the spherical coordinates given in the laboratory frame are Lorentztransformed to the fluidrest frame. It is evident that they are no longer spherically symmetric and distorted in the latter frame. This picture summarizes the difficulties in the treatment of scatterings even if they are isoenergetic. As is well known, the neutrino distribution function is a Lorentz invariant and its values at corresponding points in different frames are identical. The important point, however, is the fact that grid points are shifted by Lorentz transformations and concentric (equivalently isoenergetic) spheres in the laboratory frame are no longer spheres in the fluidrest frame. As a consequence, some interpolations are inevitable in evaluating the collision terms for scatterings in the fluidrest frame if one were to avoid the expansion. There are several difficulties to carry out this interpolation particularly in neutrino energy, though. The reasons are described shortly.
The rather low energy resolution we can afford in the Boltzmann code is one of the reasons. We can deploy at most energy bins (see Kotake et al. (2012b)). The distribution function depends strongly on the neutrino energy in general. In particular, it decreases almost exponentially at high energies. On the numerical mesh, may change several orders of magnitude between adjacent energygrid points. Highly accurate interpolations of are hence required on the coarse mesh. Note that since the isoenergetic scatterings between neutrinos and nucleons and/or nuclei dominate other reactions in CCSNe, the time step () of simulations is mostly determined by these processes. If the interpolations of are not accurate at high energies, we might find that becomes unreasonably small because of a large number of artificial scatterings. The fact that high energy neutrinos have larger cross sections makes matter worse. Not to mention, in the interpolation we further have to care about the conservation of neutrino numbers in scatterings.
After giving the SR Boltzmann equations in the next section, we present our idea to overcome these difficulties. We then demonstrate our successful handling of the isoenergetic scatterings in the realistic supernova simulations (see Section 7).
4. SR Boltzmann equations for neutrinos
We start with the covariant form of Boltzmann equation:
(1) 
which is valid even in curved spacetime. In the above expression, () denotes the neutrino distribution function in phase space; and are spacetime coordinates and fourmomentum of neutrino, respectively; since the latter satisfies the onshell condition: , in which is a neutrino mass, only three of four components are independent and this is why only spatial components appear in the second term on the left hand side; stands for the affine parameter of neutrino trajectory. The left hand side of Eq. (1) expresses a geodesic motion in the phase space, while the right hand side denotes symbolically the socalled collision terms, i.e., the terms that give the rate of changes in due to neutrinomatter interactions.
On the spherical coordinates in flat spacetime, which are the coordinates we employ for the laboratory frame in our Eulerian approach, Eq. (1) is cast into the following conservation form:
(2) 
where , , denote the spatial variables; as three independent components of neutrino fourmomentum, we do not use its spacial components but adopt energy and two angles, and (see Figure 3); is defined as . In Eq. (2) and the rest of this paper, we assume that neutrinos are massless, which is well justified as long as neutrino oscillations are ignored.
The collision term in Eq. (2), which is expressed with the laboratory time , is related with the original collision term in equation (1) as
(3) 
where denotes the neutrino energy measured in the laboratory frame. Similarly, the collision term in the fluidrest frame can be expressed with the proper time of each fluid element () as
(4) 
where denotes the neutrino energy in the fluidrest frame. Here is the fourvelocity of matter.
The Lorentz transformation of fourmomentum gives the relation of neutrino energies in the fluidrest and laboratory frames as
(5) 
where , denote the threevelocity and corresponding Lorentz factor of matter and is the unit vector that indicates the flight direction of neutrino in the laboratory frame. The factor in Eq. (5) expresses the Doppler shift of neutrino energy. From Eqs. (3)  (5), we can obtain the relation between the collision terms in the two frames as
(6) 
The Lorentz transformation also gives the relation between the flight directions in the fluidrest and laboratory frames as
(7) 
Here denotes the unit vector that specifies the flight direction of neutrino in the fluidrest frame. Using the Doppler factor , we obtain
(8) 
Note that this relation no longer contains neutrino energy and the angletransformations are decoupled from the energy transformations. This is a great simplification, which we make full use of in the following, and is a consequence of the assumption that neutrinos are massless. The solidangle element is then transformed as
(9) 
In the Boltzmann equation, neutrinomatter interactions are described in the collision terms. As is well known, they are obtained most easily in the fluidrest frame. We hence evaluate the collision term in this frame and use Eq. (6) to obtain the expression in the laboratory frame. The interactions that we take into account in this paper are the same as those in Sumiyoshi & Yamada (2012), the minimum set for supernova simulations. Since Sumiyoshi & Yamada (2012) worked in the Newtonian approximation, we need the following replacements to employ their collision terms:
(10) 
where denotes reaction kernels.
Here we take the collision terms for the isoenergetic scatterings in the laboratory frame and see how the neutrinonumber conservation is ensured, which will be useful in the next section. Following Sumiyoshi & Yamada (2012) and implementing the above replacements, we can write them as
(11)  
where and denote the isoenergetic scatteringkernel and neutrino distribution function in the fluidrest frame, respectively. The integration of Eq. (11) over the solid angle vanishes due to symmetric properties of scattering kernel: . This represents the conservation of neutrino number for the isoenergetic scatterings at each energy in the fluidrest frame.
5. Two energygrids
The origin of difficulties in the SR treatments is the fact that the neutrino momentum space is distorted by Lorentz transformations, i.e., the isoenergy surfaces in the laboratory frame do not coincide with the counterparts in the fluidrest frame. We then need highly accurate interpolations in energy of , taking care of neutrinonumber conservation, whose difficulties in the method were elucidated in Section 3.
We overcome these difficulties by employing not the grid shown in the left panel of Fig. 2 but the socalled Lagrangian remapped grid (hereafter LRG) in the laboratory frame, which is Lorentztransformed from the fluidrest frame. It is emphasized that LRG is the one we mainly use in our Eulerian approach. In Figure 4, we display the schematic picture of our LRG (see also Figure 2 for comparison). In this method, the energy grid is constructed so that it should be isotropic in the fluidrest frame.
As a consequence, it becomes anisotropic in the laboratory frame as observed in the left panel. The energy grids obtained in the laboratory frame that way are different from point to point at each time and change also in time because of inhomogeneous fluid motions. Thanks to the isotropic energy grids in the fluidrest frame, no special care is needed in the treatment of the isoenergetic scatterings on this grid. Note that the angular mapping is independent of energy. The angular grid is constructed, on the other hand, so that it should be uniform in the laboratory frame. It implies that the angular mesh is nonuniform in the fluidrest frame as shown in the right panel. In contrast to the energy grid, the nonuniform angular grid in the fluidrest frame causes no practical problems (see Eq. (32) for the correction by angular aberration).
One may say that the Lagrangian remapping method is nothing but the canonical Lagrangian approach, but there are several differences between the two. One of the important differences lies in the treatment of advection terms on the left hand side of the Boltzmann equation. As we have already mentioned in Section 1, the advection terms are fairly complicated in the Lagrangian approach due to the spatial inhomogeneities and temporal changes of matter velocity. We demonstrate this in a simplified situation in Figure 5. Here, we consider neutrinos propagating outwards (or rightwards in the figure) in optically thin matter. We further assume that the matter is moving inwards (or leftwards in the figure) at velocities that are piecewise constant with (see the bottom panel in this figure). The discontinuity may be regarded as a standing shock wave.
In this situation, the neutrinoenergy spectrum in the laboratory frame is uniform in space, since neutrinos are not interacting with matter at all (see the upper picture)
It is easily understood that the use of LRG, which is anisotropic and spatially nonuniform, complicates the calculation of spatial and neutrinoangular advection, a problem similar to the one in the ordinary Lagrangian method. This is mitigated in our method, however, by the introduction of yet another energy grid, which is isotropic and identical at all spatial grid points in the laboratory frame (referred hereafter to as the laboratoryfixed grid or LFG; see also Section 6.4 for more details). LFG is employed only to calculate the neutrino advection. Note that as long as we work in the laboratory frame, energyderivative terms do not appear explicitly on the left hand side of the Boltzmann equation and the advection on the LFG is particularly simple. It should be repeated that LFG is a grid only for temporary use to treat the neutrino advection. Accordingly on the LFG, which is obtained by interpolation in our method, is also a temporal variable. Instead, on LRG is the quantity to be solved and stored in our code.
6. Numerical Implementations
In this section, we explain the detailed numerical algorithm to implement various elements described above to our BoltzmannHydro solver, paying particular attention to the usage of different energy grids. Figure 6 summarizes multiple steps needed to update a numerical solution from to , where the superscripts represent the time steps. In the following, we describe each step in order in detail.
6.1. Step 1: Hydrodynamical evolutions
In our BoltzmannHydro solver, the operator splitting is employed: we first compute hydrodynamics, neglecting neutrino interactions, i.e., in the adiabatic manner; then from Step 2 through Step 4 we perform neutrino transfer for the matter distribution given in the first step as described below; feedbacks from neutrino interactions to the internal energy, velocity and electron fraction of matter are taken into account in Step 5.
The numerical code for the hydrodynamical evolution is essentially the same as that in Nagakura et al. (2013). It is based on the socalled central scheme with an explicit time evolution (Kurganov & Tadmor, 2000; Nagakura & Yamada, 2008; Nagakura et al., 2011). The code was successfully applied to the simulations of Standing Accretion Shock Instability (SASI) in the postbounce phase in our previous study (Nagakura et al., 2013). It is also noted that a series of standard tests for hydrodynamical schemes (e.g., shock tube problems) were carried out in Nagakura et al. (2011).
Although our Boltzmann solver is fully SR, the hydrodynamics solver is Newtonian. As a matter of fact, it is fully general relativistic (Nagakura & Yamada, 2008) except for its gravity solver, which is Newtonian and based on the MICCG technique (Nagakura et al., 2011). The implementation of an Einstein equation solver is currently underway, the perspective of which will be mentioned in Section 8.
The basic equations of Newtonian hydrodynamics in spherical coordinates are written in the following form:
(12) 
where each term is given as
(13)  
(14)  
(15)  
(16) 
Note that corresponds to the interactions between neutrinos and matter (the explicit expressions will be presented in Step 5) and denotes the volume factor in the spherical coordinates. Other variables, , , , , , and , are the mass density, pressure, internal energy density, electron fraction, fluid velocity, and Newtonian gravitational potential, respectively. The Newtonian selfgravity is solved with the Poisson’s equation,
(17) 
In our central scheme, the above system of equations is finitedifferenced in space with the piecewise parabolic (PPM) interpolation and the totalvariation diminishing (TVD) RungeKutta method is employed for time integration, which achieves secondorder accuracy in both space and time. We adopt the procedure proposed by Müller et al. (2010) in solving the energy equation (the 5th component in Eq. (12)), which reduces secular errors in the energy conservation.
Throughout this paper, we use Shen’s equation of state (EOS) (Shen et al., 2011) with lepton and photon contributions being added. (see e.g., Nagakura et al. (2013)). The original EOS table is rather coarse for the simulation of CCSNe. Indeed, we have found that trilinear interpolations in the original table reduce the accuracy of simulations particularly at the transition from inhomogeneous to homogeneous nuclear matter. We have hence reconstructed a new EOS table by interpolating all quantities with the tricubic Hermite functions. It is several times finer in and than the original table.
6.2. Step 2: Reconstruction of subgrid energy spectrum
In our Boltzmann solver, transformations between different energy grids are frequently performed. As mentioned in Section 3, we will be able to deploy at most energy bins, a rather coarse resolution. We hence need a subgrid modeling of neutrinoenergy spectrum. It is also important for computing fluxes at grid boundaries. As a matter of fact, if we did not take into account such subgrid distributions and assumed instead that neutrinos are populated uniformly in each grid, then a large number of neutrinos could artificially leak to neighboring grids either by inaccurate numerical fluxes or by imprecise interpolations (see also Step 4 on this issue).
In the reconstruction one should pay an adequate attention to the following two conditions:

monotonicity

Conservation of neutrino numbers
The first condition is familiar in the numerical treatment of hyperbolic systems and necessary to avoid artificial generation of extrema in spectra, which may cause numerical instabilities. The importance of the second condition is rather obvious. In fact, if it were violated, neutrinos would appear or disappear just by changing energy grids. As shown later, this condition is particularly important in the evaluation of on LFG. Note that the value of on each grid point actually represents the average in the energy bin in our formulation.
The reconstruction procedures are schematically shown in Figure 7, in which we are going to construct the subgrid energy spectrum for energy bin A in LRG. In so doing, not only grid point A but also neighboring grid points B and C are utilized. We distinguish two cases: (1) takes an extremal value locally on grid point A, i.e., both of ’s on grid points B and C are either larger or smaller than on grid point A; (2) otherwise.
The left panels in Figure 7 correspond to the first case. As shown there, we assume in this case a flat spectrum in the energy bin. This is not a bad approximation, since the actual spectrum is indeed nearly flat in the vicinity of a local extremum. In the second case, in which changes monotonically over the neighboring three grid points, we reconstruct a subgrid spectrum as follows, which is shown in the right panels in Figure 7.
We first determine the value of on the left and right interfaces of energy bin A as the averages of adjacent gridpoint values in the logarithmic scale. They are referred to as (), respectively. We use the grid widths for the weight in the average. We also define and as the largest (smallest) of and . Denoting the gridpoint value of by , we first construct a trial spectrum as follows:
(18)  
(19) 
where , , and are the energies at the left and right interfaces and grid point, respectively; , and are the corresponding values of given as
(20) 
It is clear that this expression becomes exactly correct if neutrinos are in thermal equilibrium and take a FermiDirac distribution.
It is obvious, however, that does not ensure the conservation of neutrino number. We hence need to modify . We first integrate in the energy bin to obtain the neutrino number, in it. This should have been equal to , the true value. Using the ratio,
(21) 
we scale the temporary spectrum as to obtain a new subgrid spectrum, which by definition satisfies the conservation of neutrino number exactly.
The new spectrum so obtained do not satisfy the monotonicity condition in general, which requires that should always lie between and . We hence apply a limiter if exceeds and/or : are modified so that they would lie in between. Owing to this limiter, obtained by integration of the new subgrid spectrum, again deviates from . We repeat the above procedure until the following condition is satisfied:
(22) 
where is a measure of convergence and is set to . It is important that the convergence of this iteration is guaranteed mathematically and that no artificial extremum emerges in the reconstructed spectrum.
6.3. Step 3: Lagrangian Remapping
Here we carry out the Lagrangian remapping of neutrinoenergy grids and compute the change in on LRG. The subgrid energy spectrum, which is obtained in the previous step, play an important in this process.
The procedure is summarized in Figure 8. Suppose that time integrations have been finished and all quantities associated with matter and neutrinos have been obtained at . The upper panel shows the gridpoint values of as well as subgrid spectra in three consecutive energy binges on LRG at . Note that the angular dimensions are suppressed in the figure. In Step 1, matter velocities are changed. As explained in Section 5, LRG is determined so that the neutrinoenergy grid be identical in each instantaneous fluidrest frame. We hence need to update LRG accordingly as shown with green lines in the middle panel of Figure 8. It is then evident that should be also changed on account of the shifts of the grid boundaries. As shown in the figure, neutrinos in the shaded areas determine the change of due to this remapping.
To evaluate the numbers of neutrinos in these regions, we use , the interface value of on the old LRG at . It is obtained as the smaller one of the two interface values derived from the subgrid spectra in the adjacent grids in order to prevent the moving of too many neutrinos. With this and the energies at the grid interfaces on the old () and new () LRG’s, the number of neutrinos to be moved to the adjacent grid is given by
(23) 
where is the extent of solid angle of the bin
It is emphasized that no interactions of neutrinos with matter have been taken into account yet up to this point. The change in considered above is induced by the acceleration or deceleration of matter (and hence of the local fluidrest frame). As explained in Section 5 in detail, if such a change in matter velocity occurred in optically thin matter, the neutrino energy spectrum should not change in the laboratory frame. The red or blue shifts of the energy spectrum in the fluidrest frame are compensated for by the Lagrangian remapping in our method. It should be also noted that the energy shift is proportional to the time step and is small as a consequence, the fact which is true even in the shock wave. This not only justifies the above estimation of but also is a huge advantage in the use of LRG compared with other Lagrangian formulations, in which large energy shifts may occur.
6.4. Step 4: Evaluations of Advection and Collision Terms
Now that the energy shift induced by matter motions has been treated, the remaining task is to consider the spatial advections and collisions of neutrinos. The latter is easy on LRG, which is essentially the comoving grid, and will be briefly explained at the end of this section. The former, on the other hand, is fairly complicated, since LRG is not uniform in space. In contrast to the ordinary Lagrangian formulation, in which the spatial advection is expressed as complicated partial derivatives, our method utilizes the fact that the advection is very simple in the laboratory frame. LFG defined in Section 5 is the main tool here.
As explained in Section 5, LFG is defined so that the energy grids are identical everywhere in space and it does not depend on the flight direction of neutrino. In addition LFG should have the following properties:

LFG covers the union of all energy ranges in LRG
^{5} . 
Each energy bin in LRG is covered by more than one energy bins in LFG.
These conditions are important to ensure the accuracy in the evaluation of the advection terms. Figure 9 displays the example of the relation between the two grids in the laboratory frame.
Given LFG, we evaluate the advection term as follows. We suppose that the nth time step has been completed and is given on each LRG and the subgrid spectrum has also been constructed according to Step 2 (see panel a in Figure 10). Using this subgrid spectrum, which is denoted by in this section, we assign to each grid point in LFG, which is denoted as :
(24) 
where corresponds to the neutrino energy at the grid point in LFG (see panel b in the same figure, in which LFG is presented in green while LRG is shown with black dots.). Here we would like to emphasize again that the spatial and angular advection
Once the numerical fluxes for the spatial and angular advections are obtained on LFG, we then calculate the corresponding numerical fluxes for LRG as follows. We take as an example panel b in Figure 10. As mentioned earlier, LFG is finer than LRG. Energy bin A in LRG, for instance, is covered by three bins C’, D’ and E’ in LFG. Let us look at bin E’ in LFG, which overlaps with bins A and C in LRG. The numerical flux in LFG should be hence shared with the latter two bins in LRG. For that purpose, we introduce a factor, , as defined shortly and divide the numerical flux into and . is defined as
(25) 
with
(26) 
where is the neutrino energy at the interface of bins A and C in LRG and is the energy at the left (right) boundary of bin E’ in LFG; is the gridpoint value of for bin A(B) in LRG. The numerical flux for bin A of LRG is a sum of the contributions from bins C’, D’ and E’ in LFG, each obtained in this manner.
So far we have explained our treatment of the spatial and angular advection terms as if they were finitedifferenced explicitly in time. As a matter of fact, we treat them implicitly in our method. This is important from a point of view of numerical stability and computational times. In the following, the detailed procedure is described.
We first rewrite the Boltzmann equation (Eq. (2)) as
(27) 
in which the second term on the left hand side is the sum of spatial and angular advection terms. Since the following treatment is common to each term in the sum, we drop the subscript hereafter.
In the implicit approach, both the advection and collision terms are evaluated at and the finitedifferenced equation is written as
(28) 
It is noted, however, that in SR is evaluated via the complex interpolation of between LFG and LRG and is nonlinear and highly complicated. This prevents us even from linearizing the equation
It is noteworthy, however, that in the Newtonian approximation, in which no distinction is made between the fluidrest and laboratory frames, the advection term is linear and can be treated completely implicitly. In fact the resultant nonrelativistic (NR) equation can be cast into the following form (see also Sumiyoshi & Yamada (2012)):
(29) 
where is the NR advection term and the subscript indicates the spatial and angular grid points symbolically. Coefficients , and are written as a function of space and flight directions. This fact suggests the following scheme for the advection term in SR:
(30) 
The first term on the right hand side is the relativistic advection term evaluated with and the second term in the parenthesis is a correction term. Eq. (30) is only semiimplicit as it is. We hence replace in Eq. (30) by a trial value, and repeatedly solve the Boltzmann equation (Eq. (27)), updating with obtained for until they coincide with each other within a certain error. This ensures fullimplicitness of our method as explained shortly.
The idea behind this method should be clear. If matter motion is not very fast compared with the speed of light, which is indeed the case in CCSNe, is almost equal to , and will be dominated by . Then this scheme is close to the NR implicit scheme. In addition, there are other important properties in the prescription: firstly if coincides with , as is the case at the end of iterations, the correction term vanishes and only remains. This property guarantees the full implicitness of our scheme. It is also mentioned that is actually limited most of times by the requirement that should not change by more than a few percent in a single time step, and the correction term is a small correction to the first term in most situations (but see below for the exceptional case.). In spite of this, we find that this correction term significantly improves the numerical stability, enabling us to take larger .
Although the above method works fairly well most of times, it fails sometimes when the matter velocity reaches several tens percent of the speed of light and the correction term becomes significantly larger than . In such cases the iteration does not converge unless we reduce . To avoid too small a value of , however, we introduce a limiter
(31) 
in which is determined so that the correction term should not exceed the first term. It is stressed that such a prescription is just a technique to improve the convergence and does not affect the final result, since the second term vanishes in the end anyway.
We turn to the collision terms before closing this section. There are no new difficulties in their treatment, since LRG is essentially the same as the fluidrest frame employed in the ordinary Lagrangian methods. There is one feature, however, which is original in our method.
The angular dimensions in momentum space are discretized in the same manner everywhere on LRG (see the left panel of Figure 4). This means that the angular gridding is not uniform in the fluidrest frame due to aberration by Lorentz transformation. The angular integration as shown in Eq. (11) for the isoenergetic scattering is normally performed in the fluidrest frame. In our approach, however, that is done on LRG in the laboratory frame. In so doing, the aberration effect is taken account of as the Jacobian as follows
(32) 
where is an arbitrary function of and is defined as . The Jacobian () has already been derived in Eq. (9).
As mentioned earlier, the collision terms are treated fully implicitly. The point is that the matrix structure originated from the collision terms is exactly the same as the one for the NR case, which implies that the numerical tools developed for our Newtonian code (Sumiyoshi & Yamada, 2012) can be utilized for the present code as they are. As a matter of fact, thanks to this implicit treatment of collision terms the time steps in the 1D test simulation of CCSNe (see Section 7.5) are comparable to those in our previous simulations with a 1D implicit GR Lagrangian BoltzmannHydro code (Sumiyoshi et al., 2005).
6.5. Step 5: Feedbacks to matter
Solving the Boltzmann equations in the previous step, we now treat feedbacks from the neutrinomatter interactions to hydrodynamics. The hydrodynamics equations and the conservation equation for electron number are written as
(33)  
(34) 
where the right hand sides correspond to the feedbacks and are written as
(35)  
(36)  
(37)  
(38) 
In these expressions, the invariant volume in the momentum space is denoted by and the subscript ”” indicates the neutrino species.
At the very end of all steps, we again perform Steps 2 and 3, since matter velocities are changed due to the momentum exchange between matter and neutrinos. This closes the update from to . We iterate these steps as many times as needed.
7. Validation
In order to validate our new formulation of SR Boltzmann RadiationHydrodynamics, we carry out a series of code tests. We first focus on the Boltzmann solver, i.e., the feedbacks to hydrodynamics are ignored. We test the advections and collisions separately in idealized setups in order to see the code performance in each sector clearly. In these tests, only electrontype neutrinos are taken into account, since the treatments of SR effects are common to other species.
We then perform SR BoltzmannHydro simulations of 1D spherical core collapse for the progenitor. In these test runs, we consider 3 species of neutrinos (, , and ) and implement minimal but essential microphysics. For comparison, a NR simulation is also performed. Based on the two results, we discuss the importance of SR effects.
7.1. Collision term: isoenergetic scattering
As discussed in Section 3, the isoenergetic scattering is the primary source of difficulties in the method for the SR Boltzmann equation. This test is meant to see whether our code can properly handle this process. This is a single zone calculation, in which we deploy only one spatial grid and the advection term is neglected. We are concerned only with the collision term. Hydrodynamical quantities are assumed to be constant in time and set as , MeV, and , where , , and denote the density, temperature and electron fraction, respectively. Under this thermodynamical condition, free nucleons and nuclei are both existent. We hence consider the following isoenergetic scatterings:
(39)  
(40) 
Although we drop the advection term in this test, we set a nonvanishing velocity as follows:
(41)  
(42)  
(43) 
where , , and denote the radial, –, and – components, respectively. They are assumed to be constant in time and are controlled by three parameters, , and . In this test, we set , , and , respectively. Note that this velocity is considerably large by the CCSNe standard.
We assume that neutrinos are distributed isotropically in the laboratory frame initially, and they have FermiDirac distributions in energy. The neutrino chemical potential can be obtained by the assumption that neutrinos are chemical equilibrium with matter. Since matter has a nonvanishing velocity, neutrinos are initially anistropic in the fluidrest frame. Then, should evolve towards an isotropic distribution in the latter frame by the isoenergetic scattering.
For this test, momentum space is covered with a grid of points in energy and points in angles. The gridding of LRG has been explained in detail in Section 5 and Figure 4 (see also Sumiyoshi & Yamada (2012) for the construction of angular grid).
We show the numerical results in Figures 11 to 15. Figure 11 displays the evolutions of for different angles but with the same energy ( MeV) in the fluidrest frame. As is expected, initially different values of are changed by the isoenergetic scatterings and converge to a certain value by the time s. Note that we work on LRG in the laboratory frame and these results are obtained by the Lorentz transformation to the fluidrest frame. The final isotropic distribution, , can be obtained analytically from the initial condition, , since isoenergetic scatterings do not change the number of neutrinos:
(44) 
where the subscript specifies the angular grid points. It is evident from the figure that the correct results are obtained numerically.
Figure 12 shows the angular evolution of on an isoenergy surface with MeV in the fluidrest frame. In the figure, wireframed pictures are drawn as follows: for each angular grid point, a node is placed at a distance proportional to the value of in the corresponding direction; the neighboring nodes are then connected by lines; we use the normalization, in which the maximum distance should be unity. As a consequence, an isotropic distribution corresponds to the unit sphere in this figure.
At the beginning (top left panel), the wire frame is elongated in one direction, indicating that the angular distribution is highly anisotropic. As time passes, however, it changes shapes and eventually (s) becomes isotropic although it may not appear so. This is due to the rather low resolution in this computation. Indeed, Figure 13, which presents the result for a higher resolution (, ), more clearly isotropy of the final distribution. It is also reminded that the angular grid is not uniform in the fluidrest frame due to aberration (it is uniform in the laboratory frame. See Section 5).
As an alternative presentation of isotropization, we show in Figure 14 the initial and final energy spectra for two different angles in the fluidrest frame. As demonstrated evidently in this figure, the initially different spectra converge at the end, implying that neutrino distributions become isotropic at all energies. In Figure 15 we show the same evolution in the laboratory frame, which is actually the frame we use for simulations. Contrary to the previous case, the initially identical spectrum for different angles is separated as time passes in the laboratory frame, indicating that the final distribution is anisotropic in this frame as it should.
7.2. Collision term: emission, absorption and isoenergetic scattering combined
To the isoenergetic scatterings, we add emissions and absorptions on nucleons and nuclei:
(45)  
(46)  
(47) 
The initial condition and computational setup are the same as those in the previous test.
Figure 16 shows the evolution of for different angles but with the same energy ( MeV) in the fluid rest frame, which corresponds to Figure 11. At first the isoenergetic scatterings isotropize the distribution in the fluidrest frame just as in the previous case. By the time s, the neutrino distribution is almost isotropic. Note that it is not a FermiDirac distribution at this time, i.e., neutrinos have not yet achieved chemical equilibrium with matter via emissions and absorptions. It is eventually established at s for this energy of neutrinos. Neutrinos with other energies undergo similar evolutions and reach FermiDirac distributions at different times. It is stressed again that this computation is done on LRG and the distribution in the fluidrest frame is obtained via a Lorentz transformation.
7.3. Advection term: 1D advection through a discontinuity in matter
We now turn to the advection term. Note that this is the main source of difficulties in the ordinary Lagrangian methods. Our formulation treats the spatial and angular advections in the laboratory frame but employs interpolations between two grids (LRG and LFG) as detailed in Section 6.4. It is hence important to confirm that the scheme indeed works properly.
Here, we consider the advection in matter that has a discontinuous velocity distribution. This is certainly the most difficult situation for our method. In contrast to the previous tests, we cut off all neutrinomatter reactions, assuming that the matter is optically thin, and consider the advection term alone. Note that in this case, the energy spectrum of neutrinos is unchanged as they propagate through the discontinuity in the laboratory frame whereas it undergoes a discontinuous change there in the fluidrest frame (see also Figure 5). We further assume spherical symmetry in this test, i.e., we omit derivative terms with respect to , , and in Eq. (2).
The computational setup is as follows. To avoid geometrical complications, we consider the advection in a waferthin spherical shell: the computational domain covers the range of cm by a uniform radial grid of 6 bins. The matter velocity is piecewise constant with a discontinuity between the 3rd and 4th grid points: v=0 for the first 3 grid points and for the rest of grid points. These velocities are again fixed during the computation. We inject outgoing neutrinos from the radial inner boundary with the FermiDirac distribution that is the same as in the previous tests and follow the subsequent evolution until a steady state is obtained. We deploy an LRG of and .
Figure 17 shows that the energy spectra for outgoing neutrinos () in the vicinity of the velocity discontinuity in the laboratory frame (left panel) and in the fluidrest frame (right panel). As is expected, neutrinos advect without any change of their spectrum when they pass through the discontinuity in the laboratory frame. We can also see in this figure that the energy bins for the outer two radial grid points are shifted from those for the inner two. The right panel shows the same spectrum but observed in the fluidrest frame. Due to the negative radial velocity in the outer region, neutrinos are blue shifted there (see also Figure 5). As demonstrated clearly in this test, our new formulation can reproduce the results just as expected without any numerical problems.
7.4. Advection term: 3D advection
This test is meant to check the MultiD advection in the optically thin matter with an inhomogeneous nonradial velocity distribution. We assume that the neutrino distribution is spherically symmetric in space. This is no problem, since the matter is optically thin and there is no interaction between the matter and neutrinos. This poses a challenge in our method, however, since LRG is not spherically symmetric in space and, as a consequence, there is no guarantee that the neutrino distribution remains spherically symmetric in our formulation. This test is hence good diagnostics on our handling of the spatial advection.
The 3D velocity distribution is set in a similar way to that in the previous test, Eq. (43), but with an additional spatial dependence:
(48)  
(49)  
(50) 
We again set a nonvanishing nonradial velocity by choosing , and . is given as follows:
(51) 
where and denote, respectively, the maximum and minimum radii of the computational region, which is the spherical shell with cm, and . We deploy to this computational domain an LRG with . In the following we demonstrate that the neutrino distribution remains spherically symmetric with this small number of spatial and angular grid points. We inject from the inner boundary outgoing neutrinos with the FermiDirac distribution employed in the previous tests. The simulation is continued until the neutrino distribution becomes steady.
We summarily display the results of this test in Figure 18. The upper left panel shows the energy spectra for different ’s (with being fixed) at and in the laboratory frame. Note that if the neutrino distribution is exactly spherically symmetric, these spectra should coincide with each other. As seen in this figure, they agree quite well despite they are computed on the LRG, which is not spherically symmetric. The upper right panel is the same as the upper left, but for . Note that these neutrinos propagate in a nonradial direction. Again their spectra depend on very weakly. Finally, we display in the bottom panel the energy spectra at a different radial location . This time, and are fixed and is varied. We can confirm also in this case that all energy spectra are in good agreement. It is emphasized again that these results are not trivial and, in fact, the test is very severe, since we assume here very fast matter motions ( of the speed of light) with large inhomogeneities. We hence think that our new method works satisfactorily.
7.5. SR BoltzmannHydro Simulations: the spherical collapse of progenitor
So far we have tested the advection and collision separately in simplified setups. In reality, however, they are nonlinearly coupled with each other and dictate the neutrino transfer and, as a consequence, the dynamics of CCSNe. In order to confirm that our new method is indeed applicable to realistic simulations of CCSNe, we conduct here a 1D spherically symmetric BoltzmannHydro simulation for the collapse of progenitor (a nonrotating star with the solar metallicity referred to as s15.0 in Woosley et al. (2002)). We employ an LRG with covering the computational domain of cm. For comparison, we also perform a NR simulation for the same setup. Although the simulation is continued after bounce until the shock wave is stalled, we focus here on the collapsing phase, since the infall velocity is largest and SR effects are most clearly discernible.
Figure 19 shows that the evolution of the number density of at the center for both the SR and NR simulations. Initially these two simulations follow almost the same evolutionary path. After the central density reaches , however, they start to deviate and become different by more than a factor of at . During the latter period, neutrinos undergo isoenergetic scatterings on nuclei called coherent scatterings and, as shown shortly, this is the source of discrepancy in fact.
In order to clearly see the SR effects by the matter motion, the left panels of Figure 20 show as a function of radius the radial component of the number flux, i.e., the energyintegrated firstangular moment of in the laboratory frame:
(52) 
where denotes the volume element of energy space in the laboratory frame. The upper panel corresponds to the time when the central density reaches whereas the bottom one shows the result at the time of , respectively. On the right panels, matter velocities are displayed as a function of time for the same times.
As is evident in the left panels, the number flux behaves qualitatively differently in the SR and NR cases: in the SR simulation is negative in the inner region (km), whereas it is positive everywhere in the NR. Simply put, neutrinos are moving in the opposite direction if SR is ignored. This is understood as follows (see also Section 2): matter is optically thick to neutrinos in the inner region and neutrinos tend to diffuse outwards as observed in the NR simulation; the matter is infalling, on the other hand, and tends to drag neutrinos inwards; this is made possible by frequent interactions between the matter and neutrinos; in fact, as demonstrated in Sections 7.1 and 7.2, scatterings and emissions/absorptions render the neutrino distribution isotropic in the fluidrest frame and, as a consequence, produce a flux in the direction of velocity in the laboratory frame after Lorentz transformation; if SR is neglected, neutrinos are isotropically distributed even in the laboratory frame and no dragging occurs; this is the cause for the discrepancy. Note that this dragging (and hence SR) is crucially important for neutrino trapping as shown next.
In Figure 21, we display the radial distribution of lepton fraction at two different times when the central density reaches and . The left panel presents the results of the SR simulation, while the right panel gives the NR counterpart. We can immediately recognize a remarkable difference. In the SR simulation, two lines are almost the same, in particular for , where denotes the mass coordinate. This means that the lepton number is conserved in each fluid element as it should after neutrino trapping. For the NR case, on the other hand, the lepton fraction is decreased even in the central region, while it is increased in the outer region. This means that neutrinos are diffusing outwards in the Lagrangian frame even after neutrino trapping, which is consistent with what we observed in the number flux above.
In the Lagrangian method, the leptonnumber conservation after neutrino trapping is handled almost trivially. Our BoltzmannHydro solver is based on the Eulerian picture, in which does evolve as a function of radius even after neutrino trapping. Only when SR is taken into account appropriately, can we reproduce the correct evolutions. The results of this test simulation are hence the clearest evidence that our new method properly handles the neutrino advection. Incidentally, we have also made a comparison with the result of 1D Lagrangian GR simulations (Sumiyoshi et al., 2005)
Finally, we present the mass shell trajectories during the postbounce phase (until 150 ms after the bounce) in Figure 22. After bounce, the shock wave propagates outwards initially through optically thick matter and generates a neutronization burst of when it breaks out of the neutrino sphere and is eventually stagnated at a certain radius in optically thin matter. These phases are, in general, difficult numerically for our method, since neutrino distributions evolve quite rapidly, and both optically thick and thin regions are involved, and a shock wave, i.e. a discontinuity in matter velocities, exists. In spite of these difficulties, our SR BoltzmannHydro code has run stably without problems. Although we have to wait for more detailed quantitative analyses of this model and others in multiD, which will be reported elsewhere as a sequel, the results shown so far indicate that our new code will be applicable to realistic CCSNe simulations.
8. Summary and possible extensions of the formulation
In this paper we have presented a novel method to solve numerically the SR Boltzmann equation in the laboratory frame based on the method, which overcomes technical difficulties inherent to the conventional approaches irrespective of the Lagrangian or Eulerian pictures. Our method is hybrid, deploying the Lagrangian remapped grids in the Laboratory frame. The employment of LRG simply solves the difficulties in the treatment of scatterings, which plague the conventional Eulerian approaches. As a tradeoff, the numerical treatment of the advection term becomes complicated as in the ordinary Lagrangian approaches. This problem is mitigated by the use of LFG, which is nothing but the ordinary grid fixed to the laboratory frame and adopted in the conventional Eulerian approaches. The advection becomes simplest on LFG. We have developed a scheme for the interpolation between LRG and LFG, which ensures the neutrino number conservation.
By carrying out a series of code tests, we have demonstrated that our new method works as expected, correctly handling both the collision and advection terms. With the same code, we have also conducted a 1D CCSNe simulations from core collapse through bounce till shock stall for a realistic progenitor model of with the minimal set of microphysics. We have paid particular attention to the collapsing phase, in which matter velocities reach maximum and our code faces the greatest challenge. We have found that the neutrinodragging due to matter motions, which is crucially important in neutrino trapping, is correctly captured in the SR simulation but not in the NR one. We have also observed only in the SR computation that the lepton fraction as a function of the Lagrangian mass coordinates does not change in time in the optically thick region. These results clearly indicate that the adequate treatment of SR effects is critically important to obtain the lepton fraction correctly.
The simulation was continued until the shock wave generated at bounce is stalled in the core. We have found no problem in the later phase, either, and we are now confident that our new method is applicable to the realistic simulation of CCSNe. In fact, we have already started such simulations in 2D and their results will be reported together with further tests in multiD in our forthcoming paper. It is finally stressed that our method could be applied to other more relativistic phenomena such as photon transfer in AGN or GRBs, since SR effects are taken into account to all orders of in our Boltzmann code. These possibilities will be studied in the future.
At the very end of the paper we comment on the extension of our formulation to GR BoltzmannHydro simulations. We have recently published a paper on the conservative form of GR Boltzmann equation (Shibata et al., 2014), which, in flat space time, would reduce to the one used in the current study. It turns out that our Lagrangian remapping method can be extended to this form of GR Boltzmann equation with some modifications. As shown in Eq. (21) of that paper, GR modifies only the advection terms with the collision terms being essentially unchanged from the SR case. In the GR case, the choice of LFG is nontrivial. We may be able to use a local tetrad with a timelike unit vector, , orthogonal to the spatial hypersurface of const. Then one important difference from the SR case is that depends on space and time, which implies that the GR Boltzmann equation has energyderivative terms on the left hand side even in the laboratory frame, which is nothing but gravitational redshift. It is noted, however, these terms may not pose problems, since gravitational field change only gradually both in time and space. Such extension is currently underway and will be published elsewhere.
Footnotes
 Of course, nonisoenergetic scatterings on electrons and neutrinos and pair processes are another complication, which will be addressed in the future.
 It is assumed here that the boundary condition is fixed and a steady state has been established.
 In numerical simulations, such a discontinuity is somewhat smeared and the flux becomes always finite.
 Note again that we suppress the angular dimension in Figure 8.
 Note that the energy range covered by LRG depends on the spatial position and flight direction.
 Hereafter the angular advections mean the advection of neutrinos in the twodimensional momentum subspace spanning all the flight directions.
 In this method the upwind and central finitedifferences are interpolated according to the optical depth. In so doing, we introduce the weighting factor, , which is linearly interpolated from LRG to LFG in the present formulation. Since takes a value in the range of 0.5 to 1 and does not strongly depend on the neutrino energy unlike , the simple linear interpolation is justified.
 But the collision terms can be easily treated implicitly. See the discussion of the end of this section.
 See Eqs. (12)–(16)
 In the comparison, we turn off the electron scattering in the Lagrangian GR simulation. Note that GR effect is negligible for before bounce.
References
 Bruenn, S. W. 1985, ApJS, 58, 771
 Burrows, A., Livne, E., Dessart, L., Ott, C. D., & Murphy, J. 2006, ApJ, 640, 878
 Burrows, A., Dessart, L., & Livne, E. 2007, Supernova 1987A: 20 Years After: Supernovae and GammaRay Bursters, 937, 370
 Burrows, A. 2013, Reviews of Modern Physics, 85, 245
 Cardall, C. Y., Lentz, E. J., & Mezzacappa, A. 2005, Phys. Rev. D, 72, 043007
 Cardall, C. Y., Endeve, E., & Mezzacappa, A. 2013, Phys. Rev. D, 88, 023011
 Dolence, J. C., Burrows, A., Murphy, J. W., & Nordhaus, J. 2013, ApJ, 765, 110
 Hubeny, I., & Burrows, A. 2007, ApJ, 659, 1458
 Janka, H.T. 2012, Annual Review of Nuclear and Particle Science, 62, 407
 Kitaura, F. S., Janka, H.T., & Hillebrandt, W. 2006, A&A, 450, 345
 Kotake, K., Takiwaki, T., Suwa, Y., et al. 2012, Advances in Astronomy, 2012,
 Kotake, K., Sumiyoshi, K., Yamada, S., et al. 2012, arXiv:1205.6284
 Kurganov, A., & Tadmor, E. 2000, Journal of Computational Physics, 160, 241
 Lentz, E., Bruenn, S. W., Harris, J. A., et al. 2012, Nuclei in the Cosmos (NIC XII),
 Liebendörfer, M., Rampp, M., Janka, H.T., & Mezzacappa, A. 2005, ApJ, 620, 840
 Marek, A., & Janka, H.T. 2009, ApJ, 694, 664
 Mezzacappa, A., & Bruenn, S. W. 1993, ApJ, 405, 669
 Mezzacappa, A., Liebendörfer, M., Messer, O. E., et al. 2001, Physical Review Letters, 86, 1935
 Misner, C. W., & Sharp, D. H. 1964, Physical Review, 136, 571
 Müller, B., Janka, H.T., & Dimmelmeier, H. 2010, ApJS, 189, 104
 Müller, B., Janka, H.T., & Heger, A. 2012, ApJ, 761, 72
 Müller, B., Janka, H.T., & Marek, A. 2012, ApJ, 756, 84
 Nagakura, H., & Yamada, S. 2008, ApJ, 689, 391
 Nagakura, H., Ito, H., Kiuchi, K., & Yamada, S. 2011, ApJ, 731, 80
 Nagakura, H., Yamamoto, Y., & Yamada, S. 2013, ApJ, 765, 123
 Ott, C. D., Burrows, A., Dessart, L., & Livne, E. 2008, ApJ, 685, 1069
 Peres, B., Penner, A. J., Novak, J., & Bonazzola, S. 2013, arXiv:1307.1666
 Shen, H., Toki, H., Oyamatsu, K., & Sumiyoshi, K. 2011, ApJS, 197, 20
 Shibata, M., Nagakura, H., Sekiguchi, Y., & Yamada, S. 2014, Phys. Rev. D, 89, 084073
 Sumiyoshi, K., Yamada, S., Suzuki, H., et al. 2005, ApJ, 629, 922
 Sumiyoshi, K., Yamada, S., & Suzuki, H. 2007, ApJ, 667, 382
 Sumiyoshi, K., & Yamada, S. 2012, ApJS, 199, 17
 Sumiyoshi, K., Takiwaki, T., Matsufuru, H., & Yamada, S. 2014, arXiv:1403.4476
 Suwa, Y., Kotake, K., Takiwaki, T., et al. 2010, PASJ, 62, L49
 Takiwaki, T., Kotake, K., & Suwa, Y. 2012, ApJ, 749, 98
 Takiwaki, T., Kotake, K., & Suwa, Y. 2013, arXiv:1308.5755
 Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Reviews of Modern Physics, 74, 1015