A New MultiEnergy Neutrino RadiationHydrodynamics Code in Full General Relativity and Its Application to Gravitational Collapse of Massive Stars
Abstract
We present a new multidimensional radiationhydrodynamics code for massive stellar corecollapse in full general relativity (GR). Employing an M1 analytical closure scheme, we solve spectral neutrino transport of the radiation energy and momentum based on a truncated moment formalism. Regarding neutrino opacities, we take into account a baseline set in stateoftheart simulations, in which inelastic neutrinoelectron scattering, thermal neutrino production via pair annihilation and nucleonnucleon bremsstrahlung are included. While the Einstein field equations and the spatial advection terms in the radiationhydrodynamics equations are evolved explicitly, the source terms due to neutrinomatter interactions and energy shift in the radiation moment equations are integrated implicitly by an iteration method. To verify our code, we first perform a series of standard radiation tests with analytical solutions that include the check of gravitational redshift and Doppler shift. A good agreement in these tests supports the reliability of the GR multienergy neutrino transport scheme. We then conduct several test simulations of corecollapse, bounce, and shockstall of a 15M star in the Cartesian coordinates and make a detailed comparison with published results. Our code performs quite well to reproduce the results of fullBoltzmann neutrino transport especially before bounce. In the postbounce phase, our code basically performs well, however, there are several differences that are most likely to come from the insufficient spatial resolution in our current 3DGR models. For clarifying the resolution dependence and extending the code comparison in the late postbounce phase, we discuss that nextgeneration Exaflopsclass supercomputers are at least needed.
1 Introduction
Neutrino transport is an essential ingredient for numerical simulations to clarify the theory of corecollapse of massive stars and the formation mechanisms of compact objects (see, e.g., Foglizzo et al. (2015); Mezzacappa et al. (2014); Burrows (2013); Janka (2012); Kotake et al. (2012b) for recent reviews). In the collapsing iron core, neutrinos play crucial roles in every phase; deleptonization in the core determines the time of bounce and the mass of the protoneutron star (PNS) (e.g., Langanke et al. (2003)); the gigantic internal energy tapped in the PNS is almost completely carried away by neutrinos, a tiny fraction of which contributes to the heating of the postshock material, leading to the onset of corecollapse supernovae (CCSNe) in the context of the neutrinoheating mechanism (Wilson, 1985; Bethe & Wilson, 1985). Since these SN neutrinos are generally not in both thermal and chemical equilibrium, the distribution of neutrinos in the phase space should be computationally determined. This is a seven dimensional (7D) problem; (3D+1D in spacetime and 3D in momentum space), the reason why CCSN simulations are considered as one of the most challenging subjects in computational astrophysics.
Primarily due to neutrinodriven convection (e.g., Bethe (1990); Herant et al. (1994); Burrows et al. (1995); Janka & Müller (1996); Müller & Janka (1997)) and the standingaccretionshock instability (SASI, e.g., Blondin et al. (2003); Foglizzo et al. (2006, 2007); Ohnishi et al. (2006); Blondin & Mezzacappa (2007); Iwakami et al. (2008, 2009); Fernández & Thompson (2009); Guilet & Foglizzo (2012); Hanke et al. (2012); Foglizzo et al. (2012); Couch (2013); Fernández et al. (2014)), the postbounce cores are essentially of multidimensional (multiD) nature.^{1}^{1}1Recently multidimensionalities in the precollapse core (e.g., Meakin et al. (2011)) are also attracting much attention (e.g., Couch & Ott (2013); Müller & Janka (2014)). Due to the high compactness of the PNS, these multiD fluid motions are all under the influence of the general relativistic (GR) gravity, the consideration of which used to be standard in the pioneering era of CCSN simulations (e.g., May & White (1966); Schwartz (1967)). In rapidly rotating supernova cores (e.g., Woosley & Bloom (2006)), magnetohydrodynamics (MHD) instabilities naturally make the core dynamics intrinsically nonspherical (e.g., Ardeljan et al. (2000); Kotake et al. (2004, 2006b); Obergaulinger et al. (2006, 2014); Burrows et al. (2007); Masada et al. (2012, 2014); Sawai et al. (2013); Nakamura et al. (2014a); Iwakami et al. (2014)). All in all, in order to have the final word on the CCSN mechanisms quantitatively, one needs to deal with the 7D Boltzmann neutrino transport simulations in full GR MHD. Unfortunately, however, this still remains to be computationally very demanding even using exascale computing platforms in the next decade(s) to come (see discussions in Kotake et al. (2012a))^{2}^{2}2We here mean the feasibility of 7DGR Boltzmann neutrino transport simulation following 1 s after bounce with sufficient numerical resolutions in both space and momentum space..
Since the late 1990s, the ultimate sphericallysymmetric (1D) simulations where the GR Boltzmann transport is coupled to 1DGR hydrodynamics have been made feasible by Liebendörfer et al. (2001) and Sumiyoshi et al. (2005) (see Mezzacappa & Matzner (1989); Mezzacappa & Bruenn (1993b, c, a); Yamada (1997); Yamada et al. (1999); Bruenn et al. (2001); Liebendörfer et al. (2001, 2004) for the code developments). Unfortunately, however, these fullfledged 1D simulations fail to produce explosions except for superAGB stars at the lowmass end (Kitaura et al., 2006). In the context of the full Boltzmann calculations, Livne et al. (2004) were the first to perform twodimensional (2D) multiangle (i.e., assuming axisymmetry in both space and momentum space) neutrino hydrodynamics simulations using the discrete ordinate () method. Then it was demonstrated by Ott et al. (2008) that the multiangle transport is really important especially when the neutrino radiation field deviates significantly from spherical symmetry such as in the rapidly rotating cores (see also Brandt et al. (2011)). Going beyond the assumption of axisymmetry in the multiangle transport, Sumiyoshi & Yamada (2012) were the first to develop the fully multiangle Boltzmann transport code and then apply it for static backgrounds (Sumiyoshi et al., 2014). More recently, Nagakura et al. (2014) extended the code to include special relativistic (SR) corrections and showed the ability of the code by performing 1D corecollapse simulation of a 15 model. Albeit not yet implemented in hydrodynamics simulations, several novel formulations and schemes of the full Boltzmann equation have recently been reported in Cardall et al. (2013b, a); Shibata et al. (2014), and Peres et al. (2014).
At present, direct discretization of the Boltzmann transport equation fully into the neutrino angle and energy (such as by the method mentioned above) is still computationally very expensive. An approximation often made in previous works is to remove the full angular dependence of the Boltzmann equation by expanding the neutrino distribution function as a series of moments. The simplest version, in which one closes the moment expansion after the zeroth angular moment, is multigroup flux limited diffusion (MGFLD) scheme (e.g., Bruenn (1985); Livne et al. (2004); Kotake et al. (2006a); Swesty & Myra (2009); Zhang et al. (2013); Bruenn et al. (2013)). In FLD schemes, the basic equation is a diffusion equation for the neutrino energy density. In solving it, an appropriate fluxlimiter should be employed (e.g., Minerbo (1978); Pomraning (1981); Levermore (1984); Janka (1992)) to ensure the causality of the radiation field in the transparent regions where the diffusion approximation breaks down. The isotropic diffusion source approximation (IDSA) scheme developed by Liebendörfer et al. (2009) is basically positioned at a similar approximation level compared to the MGFLD scheme. In the IDSA, the neutrino distribution function is divided into two components, the one which is trapped with matter and has isotropic distribution function and the other in the freestreaming limit, each of which is solved independently, while satisfying the 1D Boltzmann equation as a whole. Due to the high computational efficiency, the IDSA has been extensively employed in both 2D (Suwa et al., 2010, 2011, 2013, 2014; Nakamura et al., 2014b) and 3D simulations (Takiwaki et al., 2012, 2014). One can also truncate the angular moment at the second order and transport the zeroth and first order angular moment. In this case, higher or equal to the second order moment needs to be determined independently to close the set of the transport equations. In the M1 moment scheme (e.g., Pons et al. (2000); Shibata et al. (2011)), one applies an analytic formula for the closure relation (see examples applied in postNewtonian MHD simulations (Obergaulinger et al., 2014) and GR simulations in 1D (O’Connor & Ott, 2013; O’Connor, 2014) and in 3D (Kuroda et al., 2012, 2014)). In contrast, one can selfconsistently determine the closure relation by the variable Eddington factor (VEF) method (e.g., Müller et al. (2010) and see references therein). In these cases, a model Boltzmann equation is integrated to iteratively obtain the solution up to the higher moments (i.e., the Eddington tensor) until the system is converged. Currently, the state of the art in multiD simulations of CCSNe is defined by multigroup (spectral) neutrino hydrodynamics simulations. More severe approximations include gray transport (Fryer et al., 1999; Scheck et al., 2006) or the lightbulb and leakage schemes (e.g., Janka & Müller (1996); Ruffert et al. (1996); Rosswog & Liebendörfer (2003); Kotake et al. (2003); Murphy & Burrows (2008); O’Connor & Ott (2011); Perego et al. (2014)), which have been often employed in many recent studies of multiD instabilities and the MHD mechanism of CCSNe.
In addition to the multiD and multiangle/truncated neutrino transport, the accurate treatment of GR is highly ranked among the todo lists towards the ultimate CCSN simulations. In most of previous multiD models with multigroup neutrino transport, GR effects are attempted to be modelled by using a modified gravitational potential that takes into account a 1D GR effect by replacing the monopole term of Newtonian gravity with the TOV potential (Buras et al., 2006b, a; Marek & Janka, 2009; Bruenn et al., 2009; Hanke et al., 2013). While there are a number of GR corecollapse simulations in 2D (e.g., Dimmelmeier et al. (2002); Shibata & Sekiguchi (2004); Müller et al. (2012b)) and in 3D (e.g., Shibata & Sekiguchi (2005); Ott et al. (2012); Kuroda et al. (2012, 2014); Mösta et al. (2014)), many of them especially in 3D have been made, for the sake of computational cost, to employ a simplified microphysics such as by the parametrization scheme (Liebendörfer et al., 2005) or by the neutrino leakage scheme (Sekiguchi, 2010). In our previous study (Kuroda et al., 2012), we performed 3D GR/SR hydrodynamics simulations of a star with the gray M1 scheme. We demonstrated that due to deeper gravitational well of GR the neutrino luminosity and the average neutrino energy in the postbounce phase increase when switching from SR to GR hydrodynamics. Since the neutrino heating rates in the postshock regions are sensitively affected by the emergent neutrino spectra, the GR effect whether it will or will not help the onset of neutrinodriven explosions needs to be investigated by multiD GR simulations with more sophisticated neutrino transport scheme.^{3}^{3}3It is worth mentioning that 2D GR models with the VEF method tend to explode more easily than the corresponding 2D Newtonian models with and without the GR correction (e.g. Müller et al., 2010, 2012b; Müller & Janka, 2014). This may support the speculation that GR is helpful for the working of the neutrino mechanism in multiD simulations.
In this paper, we present a new 3DGR radiation hydrodynamics code that is meant to apply for stellar corecollapse simulations. The spacetime treatment is based on the ArnowittDeserMisner 3+1 formalism and we employ the BaumgarteShapiroShibataNakamura (BSSN) formalism (Shibata & Nakamura, 1995; Baumgarte & Shapiro, 1999) to evolve the metric variables. The full GR radiationhydrodynamics equations are evolved in a conservative form, in which we solve the energydependent set of radiation moments up to the first order with the M1 moment scheme. This part is based on the partial implementation of the Thorne’s moment formalism (Thorne, 1981), which is extended by Shibata et al. (2011) in a more suitable manner applicable to the neutrino transport problem. Regarding the neutrinomatter interaction terms, we employ a baseline set of weak interactions as given in Bruenn (1985) and Rampp & Janka (2002) where nucleonnucleon bremsstrahlung is additionally taken into account. Our newly developed code is designed to evolve the Einstein field equation together with the GR radiation hydrodynamic equations in a selfconsistent manner while satisfying the Hamiltonian and momentum constraints. A nested structure embedded in the 3D Cartesian computational domain enables us to follow the dynamics starting from the onset of gravitational collapse of a 15 star (Woosley & Weaver, 1995), through bounce, up to about 50 ms postbounce. Since, it is still computationally too expensive to follow longterm evolution in full 3D until the neutrinodriven explosion takes place (e.g., at the earliest ms after bounce (Bruenn et al., 2009; Marek & Janka, 2009) or ms in 2D GR calculation (Müller & Janka, 2014)), we mainly focus on detailed comparisons between our pseudo1D neutrino profiles computed in the 3D Cartesian coordinates and previous 1D results to check the validity of our new code in the early postbounce phase.
This paper is organized as follows: In section 2, after we shortly introduce the formulation of the GR transport scheme, we describe the governing equations of hydrodynamics and neutrino transport in detail. Some practical implementation schemes how to satisfy important conservative quantities such as lepton number, energy, and momentum are given in section 3. The main results and detailed comparisons with previous studies are presented in section 5. Note that geometrized unit system is used in sections 2 and 3, i.e. the speed of light, the gravitational constant and the Planck constant are set to unity: , and cgs unit is used in section 5. Greek indices run from 0 to 3 and Latin indices do from 1 to 3.
2 Formalism
This section starts with a brief summary of the basic equations and the numerical schemes of GR radiation hydrodynamics.
Following our previous work (Kuroda et al., 2012), our code consists of the following three parts, where the evolution equations of metric, hydrodynamics, and neutrino radiation are solved, respectively. Each of them is solved in an operatorsplitting manner, but the system evolves selfconsistently as a whole satisfying the Hamiltonian and momentum constraints. Regarding the metric evolution, the spatial metric (in the standard (3+1) form: with and being the lapse and shift, respectively) and its extrinsic curvature are evolved using the BSSN formulation (Shibata & Nakamura, 1995; Baumgarte & Shapiro, 1999) (see Kuroda et al. (2012, 2014) for more details).
2.1 Radiation Hydrodynamics
Major differences compared to our previous code (Kuroda et al., 2012) are, on the one hand we evolved energy integrated (“gray”) neutrino radiation field with an approximate description of neutrinomatter interaction based on the neutrino leakage scheme, on the other hand we now solve the spectral neutrino transport where the source terms are treated selfconsistently following a standard procedure of the M1 closure scheme (Shibata et al., 2011). For convenience, we briefly summarize the basic equations of our newly developed code in the following (see Shibata et al. (2011) and Cardall et al. (2013a) for the complete derivation).
The total stressenergy tensor is expressed as
(1) 
where and is the stressenergy tensor of fluid and energydependent neutrino radiation field, respectively. Note in the above equation, summation is taken for all species of neutrinos () with representing heavylepton neutrinos (i.e. and their antiparticles), and represents neutrino energy measured in the comoving frame with the fluid. For simplicity, the neutrino flavor index is omitted below.
With introducing radiation energy (), radiation flux () and radiation pressure (), measured by an Eulerian observer or (, and ) measured in a comoving frame, can be written in covariant form as
(2)  
(3) 
In the above equations, is a unit vector orthogonal to the spacelike hypersurface and is the four velocity of fluid. In the truncated moment formalism (Thorne, 1981; Shibata et al., 2011), one evolves radiation energy () and radiation flux () in a conservative form and is determined by an analytic closure relation (e.g., Eq.(6)). The evolution equations for and are given by
(4) 
and
(5) 
respectively. Here is the determinant of the three metric and is the source term for neutrino matter interactions (see appendix A for the currently implemented processes). is defined by where denotes the third rank moment of neutrino distribution function (Shibata et al. (2011) for the explicit expression).
By adopting the M1 closure scheme, the radiation pressure can be expressed as
(6) 
where represents the variable Eddington factor, and corresponds to the radiation pressure in the optically thin and thick limit, respectively. They are written in terms of and (Shibata et al., 2011). Following Minerbo (1978); Cernohorsky & Bludman (1994) and Obergaulinger & Janka (2011), we take the variable Eddington factor as
(7) 
where
(8) 
In Eq.(8), is the projection operator. As we will discuss later, by the definition of in Eq.(8), one can appropriately reproduce several important neutrino behaviours, for example, neutrino trapping in the rapidly collapsing opaque core. We iteratively solve the simultaneous equations (78) to find the converged solution of .
The hydrodynamic equations are written in a conservative form as,
(9)  
(10)  
(11)  
(12) 
where , , , , and . is the rest mass density, is the Lorentz factor, is the specific enthalpy, , , is the electron fraction ( is the number density of ), and are the specific internal energy and pressure of matter, respectively and is the atomic mass unit. and are given by an equation of state (EOS) with denoting the entropy per baryon. We employ an EOS by Lattimer & Swesty (1991) (LS220, see section section 5.1 for more detail). In the right hand side of Eq.(11), represents the covariant derivative with respect to the three metric .
2.2 Conservation of Energy and Lepton Number
As explained in Section 2.1, the formalism of our code that treats the radiationhydrodynamics equations in a conservative form is suitable to satisfy the energy conservation of the total system (neutrinos and matters, see also Kuroda & Umeda (2010); Kuroda et al. (2012) for more details). Let us first show how the energy conservation law is obtained in our code.
To focus only on the energy exchange between the matter and neutrino radiation field, we omit the gravitational source term in the following discussion. Then, the equations of energy conservation of matter and neutrinos (e.g., Eqns. (4) and (11)) become
(13)  
(14) 
From the above two equations, one can readily see that the total energy (sum of matter and neutrinos) contained in the computational domain, , is conserved in our basic equations as long as there is no net energy flux through the numerical and momentum space boundaries (i.e. ).
The lepton number conservation needs to be satisfied with good accuracy because it determines the PNS mass and the postbounce supernova dynamics. We here explain how we treat it in our code. As for the electron and neutrino number conservation, the basic equations are given by
(15)  
(16) 
where
(17)  
(18) 
with . The conservation equation for neutrino number (16) corresponds to Eq.(3.22) (divided by ) in Shibata et al. (2011). Since the neutrino number density measured by an Eulerian observer is expressed as
(19) 
the equation of total lepton number conservation becomes
(20) 
Here the total lepton fraction is defined by
(21)  
From Eq.(20), one can readily see that the total lepton number is conserved irrespective of the included neutrino matter interaction processes in case that there is no net flux through the numerical and energy space boundaries. The distribution of into the each component (e.g., ) is determined by the details of the implemented microphysics, which should be checked carefully and will be reported in Sec. 5.
2.2.1 Neutrino Number Transport in the Diffusion Limit
In the collapsing iron core, it is well known that the central core becomes opaque to neutrinos due mainly to scattering off heavy nuclei when the central density exceeds g cm (Sato, 1975) and neutrinos are trapped with matter afterward. In the diffusion limit at large neutrino opacities, the trapped neutrinos move with the matter velocity for an Eulerian observer. Thus their advection equation can be described with the same form of as
(22) 
Because the source terms in the above equation amount equal to the negative value of the source terms in the electron number conservation equation, the core lepton number is conserved in a good accuracy until it gradually decreases by diffusion in the PNS cooling phase. Since the central core mass depends on the core lepton number, the CCSN simulation should be able to capture this important phenomena appropriately.
In our formalism, however, this is not a trivial problem because we solve the energymomentum conservation equations (4)(5) and not the neutrino number conservation equation (16). In this section we check whether our basic equations can reproduce the neutrino diffusion limit adequately, i.e., they reach asymptotically to Eq.(22).
For simplicity, we assume in the following that the spacetime is flat and the typical velocity of the matter field () is much smaller than the speed of light (slow motion limit; neglecting terms higher than the second order with respect to ).
Let us first check whether Eq.(16) can satisfy the local neutrino number conservation in the trapped region. In this limit, neutrino number density at each energy bin and its flux approach
(23)  
(24) 
From these relations, it is obvious that the neutrino number density at each energy bin is transferred with the matter velocity plus the diffusion velocity and the equation of the total lepton number (Eq. (20)) in the slow motion limit becomes
(25)  
demonstrating that Eq.(20) satisfies the local lepton number conservation in the trapped region.
Next, we take the diffusion limit of Eq.(4). In this limit, should approach 0 (i.e., the radiation flux () in the comoving frame vanishes). From this, the following relation can be derived,
(26) 
where we take a simple closure relation in the diffusion limit. Inserting this into the left hand side of Eq.(4), one can get
(27)  
(28)  
(29) 
Moving from the right hand side of Eq.(27) to that of Eq.(28), we assumed that has almost no spatial gradient well below the neutrino spheres (at high opacities) in the prebounce core. This is satisfied quite well as shown in the literature (Bruenn (1985)), in which a nearly flat profile is shown in their standard models within the central core with mass 0.50.6M after the central density exceeds g cm. Note that, in Eq. (28), the third term which is originally a part of the advection terms in the energy space balances with a part of the spatial advection term and they lead to the second term in Eq. (29). From this, one can clearly see how the apparent advection speed of at each energy bin approaches the matter velocity in the diffusion limit. Then, assuming that the neutrino number density at each energy bin can be approximately expressed as
(30) 
in the slow motion limit, when we divide Eq. (29) by and integrate it in energy space, it can be summarized as
(31) 
From this, one can clearly see that the neutrino number density is transported with the same matter velocity in the diffusion limit.
3 Numerical Method
In this section we describe how to evolve the radiationhydrodynamics variables^{4}^{4}4We omit our numerical method to evolve the spacetime variables that is essentially the same as in Kuroda et al. (2012).. As we explained in the previous section, we solve Eqs. (4), (5) and (9)(12) as our basic equations which are collectively expressed as
(32) 
where denotes conservative variables
(33) 
In Eq.(32), , , and denote advection term in space, advection term in momentum space, gravitational source and neutrinomatter interaction term, respectively. Throughout this paper other than Appendix B, we divide this equation into the following two parts which are expressed in the finite difference expression as,
(34) 
for the explicit part and
(35) 
for the implicit part, respectively. In Eqs.(34)(35), is the time step size between th and th time steps and the upper indices “” represents th time step. Variables with “” denote the time updated variables during an operatorsplitting procedure.
In Eq.(34), and represent the terms with respect to advection in space and gravitational fields at th time step, both of which are added first in an explicit manner to obtain conservative variables at a middle time step . Next in Eq. (35), the rest of terms, advection in energy space () and neutrinomatter interaction terms () at th time step are added to in order to find the converged solution of by an iterative method. The reason why we separate source terms into the two parts, explicit and implicit ones, is that the typical time step size, s, which is determined by the speed of light , typical minimum grid width in our calculation m and the CourantFriedrichsLewy number CFL, e.g., 0.25, is sufficiently short for the advection term in space and the gravitational source term as well as for all the geometrical variables by an explicit update. However it is too long to follow, e.g., the weakinteraction term, which has significantly shorter time scale (s). We thus need to treat these terms in an implicit way through Eq.(35) to ensure a numerical convergence and stability.
Here, we should comment on the treatment of the advection terms in energy space . As we mentioned above, we solve them implicitly in time as Eq.(35) in our main results. It means that there is a time gap between moment of evaluation for advection terms in real space and that in energy space . Although we can generally obtain consistent results with previous studies regardless of the time gap, it is noted that the treatment of the energy advection either by the implicit or explicit scheme leads to visible changes in the postbounce features. We will discuss this point further in Appendix B.
In the following sections, we are going to describe how to evaluate the advection terms in space (Sec. 3.1) and in energy space (Sec. 3.2), and then move on to describe the implicit time update in Sec. 3.3.
3.1 Advection in Space
We employed a standard highresolutionshockcapturing scheme and utilize the HLL (HartenLaxvan Leer) scheme (Harten et al., 1983) to evaluate the numerical fluxes in space (Kuroda & Umeda, 2010). As for the fastest and slowest characteristic wave speeds of the radiation field system (Eqs.(4) and (5)), we again use the same definition as in Kuroda et al. (2012) (see also Shibata et al. (2011)) and connect and smoothly via the variable Eddington factor as
(36) 
where and are the wave speed in the optically thin and thick limit, respectively.
To enforce the numerical flux of the radiation field in the opaque region asymptotically approach to the diffusion limit, we evaluate the energy flux () and the momentum flux () as
(37) 
and
(38) 
respectively (Audit et al., 2002; O’Connor & Ott, 2013). Here, and are the conservative variables and their corresponding fluxes, respectively, with denoting the left/right states for the Riemann problem. All the radiationhydrodynamical variables are defined at the cell center. For those cell centered (primitive) variables, we use a monotonized central reconstruction, which has second order accuracy in space, and obtain left/right states at the cell surface (see Kuroda & Umeda, 2010, for more detailed explanation). After the reconstruction, all the characteristic wave speeds of the matter and radiation fields are evaluated.
3.2 Advection in Energy Space
Regarding the advection term in energy space in all conservation equations (Eqns.(4), (5) and (16)), we define the advection fluxes at the interface of the energy bin as the same as in Müller et al. (2010) so that all energy integrated advection term will vanish. This can be achieved since both terms , appearing in energy and momentum conservation equations, and in number conservation equation are expressed in terms of linear combinations of radiation momenta, , , , and we therefore can define their cell surfaced values with an appropriate weighting function to suppress violation of all conservations simultaneously.
For all orders of the radiation momenta , the following conditions need to be satisfied for number conservation,
(40) 
and for energy and momentum conservations,
(41) 
respectively. The RHS of the above equations represent the finite difference expressions with and denoting th energy bin and the interface between  and th energy bins, respectively. is energy grid width . It is straightforwad to show that Eq.(40) can be automatically satisfied for any cellsurfaced quantities as long as they vanish at the outer boundary in the energy space. By introducing a definition , the RHS of Eq.(41) can be expressed as,
(42) 
As in Müller et al. (2010), we can get the solution of Eq.(42) as
(43)  
(44) 
where is a weighting function and is expressed as
(45) 
In this study, we used a “Harmonic” interpolation () for the energy density () as
(46) 
with (see Müller et al. (2010) for more details). Like these, just only by evaluating an appropriate radiation momenta at the energy bin surface, we can simultaneously vanish Eqns.(40)(41) numerically and maintain energy, momentum and number conservations.
3.3 Implicit Time Update
After the explicit update, we solve the following simultaneous equation (e.g., Eq. (35));
(47) 
where , and are expressed in terms of the primitive variables
(48) 
at th time step. Obviously, the baryon number density does not change in this step, i.e., . To get the solutions of the above simultaneous equation, we employ the NewtonRaphson iteration scheme with the inversion of the following matrix until a sufficient convergence is achieved;
(49)  
(50) 
for with the initial condition . As for the convergence criteria for the NewtonRaphson iteration, we monitor
(51) 
where represents a tolerance and we typically set .
We note that, even if we evaluate the advection terms in energy space explicitly in time as , Eqs.(47)(51) do not change except the term in Eq.(47) is now moved to the explicit part.
Our method for time update from th to th time step is summarized in order as follows and schematically drawn in Figure 1.

Geometrical Update. We first evolve all the BSSN and the gauge variables G={, , , , , , } from th to th time step.

Explicit Update. In the second step, all the advection in space and gravitational source terms are added to obtain the fractional timestep values and their consistent primitive variables . Advection of total lepton number (Eq.(20)) is also performed simultaneously to evaluate .

Implicit Update. Finally, quantities in the fractional timestep () are implicitly updated to those in the th timestep () by the NewtonRaphson iteration until a sufficient convergence is obtained with a constraint that the local lepton fraction does not change (i.e. since ) in the diffusion region.
Here we shortly summarize some technical details to achieve a sufficient lepton number conservation in practical corecollapse simulations. As we have already mentioned in Sec. 2.2, neutrino number conservation is formally satisfied by solving the energymomentum conservation equations (4)(5). Especially, in the neutrino trapping regime, solving the energy conservation equation (4) is practically identical to solving the neutrino number conservation equation (16) which was proven in Sec. 2.2.1. However, in the finite difference method, it does not guarantee a perfect match between equations (4)(5) and equation (16) due to discretization error. To minimize the difference, we add the following constraint to Eq.(51)
(52) 
as another criterion to exit the NewtonRaphson iteration. Because of this additional criterion, the lepton number conservation is also satisfied when the NewtonRaphson iteration converges. Note that the local baryon number do not change in Eq.(47), i.e., should be met if the lepton number is conserved.
We also adopt another prescription in which is used to achieve a better convergence for the NewtonRaphson method. In some cases, happens to go beyond the range of the employed EOS table during the iteration, especially when the Jacobian matrix is not well evaluated. In those cases, we use to reset as below
(53) 
With this reset, we found that hardly goes beyond the range of the EOS table and also does not take unreasonable value for the supernova core. In addition, the number of the NewtonRaphson iteration can sometime be reduced. This is because the reset value automatically satisfies the lepton number conservation which leads to the faster convergence.
Without these two prescriptions, we observed that the central lepton fraction at bounce deviates maximally from the value immediately after neutrino trapping sets in (). Note that these adhoc constraints are introduced just to find a more efficient path toward convergence in the implicit time update. The criterion for convergence is solely given by Eq. (51). Since we do not include Eq.(20) into Eq.(47), there is no inconsistency between the number of simultaneous equations and the unknown variables .
4 Radiation Tests
In this section, we show results of several simple test problems to check the validity of our M1 radiation transport code. Except for Sec.4.4, a flat spacetime is assumed throughout this section.
4.1 Diffusion Wave Test
We first perform a diffusion wave test, by which we check the validity of our flux implementation in the diffusion limit. Following Pons et al. (2000), a Dirac function type radiation source is initially located at . Then we follow the diffusion of the source into the optically thick medium with zero absorptivity () and high scattering opacity (, ). The source term in Eqs.(4)(5) thus becomes
(54) 
The 3D computational domain is covered by 100 equidistant Cartesian zones in each direction () for a model with . This corresponds to a Peclet number . While, in the other model with , we raise the nested level to 2 to achieve sufficient resolution at the center and to check that the nested structure does not interfere the radiation propagation in the diffusion limit. In this model, numerical cells cover the base domain with 2 nested level. The minimum grid width and Peclet number at the center are thus and , respectively. Since is much higher than unity, in Eq.(39) approaches zero which enables us to check whether the radiation transport in the diffusion limit is solved appropriately.
4.2 Shadow Casting Test
Next, we move on to a shadow casting test to check the ability of our code whether anisotropic radiation field can be appropriately evolved in the free streaming limit. The initial setup is essentially the same as in Kanno et al. (2013). We set a perfect absorbing region at and inside the numerical domain and . We impose a constant radiation from the left boundary. Numerical resolution is and we set . Within the perfect absorbing region, the absorbing opacity is set as , otherwise . The source term in Eqs.(4)(5) thus becomes
(57) 
with and .
In Fig. 3, we show three different time snapshots (=3.5, 7.5 and 15.5) of the radiation energy density in a logarithmic scale. The absorbing region is expressed by whitedashed line and the radiation shock front () is denoted by vertical green line. From this figure, we see that the absorbing condition works appropriately in our scheme and the radiation front propagates with the light speed . Furthermore, the region beyond the absorbing box is barely contaminated by the radiation. These features are in good agreement with Kanno et al. (2013).
4.3 Propagation in Free Streaming Regime
The next test is a spherical expansion of radiation field from a pointlike source into optically thin medium. The aim of this test is to check whether our code using the Cartesian coordinates can maintain the sphericity of the field during the expansion. By following Just et al. (2015), we define the radiation source with radius which centres at the origin of the 2D Cartesian domain with . We also define the purely absorbing region centered at with radius . Inside those two regions, the absorptivity and equilibrium energy density are set as