A New Multi-Energy Neutrino Radiation-Hydrodynamics Code in Full General Relativity and Its Application to Gravitational Collapse of Massive Stars

Takami Kuroda1 , Tomoya Takiwaki1 and Kei Kotake1 1 Department of Physics, University of Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland Astrophysical Big Bang Laboratory, RIKEN, Saitama, 351-0198, Japan Department of Applied Physics, Fukuoka University, 8-19-1, Jonan, Nanakuma, Fukuoka, 814-0180, Japan Division of Theoretical Astronomy, National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo, 181-8588, Japan

We present a new multi-dimensional radiation-hydrodynamics code for massive stellar core-collapse 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 state-of-the-art simulations, in which inelastic neutrino-electron scattering, thermal neutrino production via pair annihilation and nucleon-nucleon bremsstrahlung are included. While the Einstein field equations and the spatial advection terms in the radiation-hydrodynamics equations are evolved explicitly, the source terms due to neutrino-matter 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 multi-energy neutrino transport scheme. We then conduct several test simulations of core-collapse, bounce, and shock-stall 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 full-Boltzmann 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 3D-GR models. For clarifying the resolution dependence and extending the code comparison in the late postbounce phase, we discuss that next-generation Exaflops-class supercomputers are at least needed.

hydrodynamics — methods: numerical — neutrinos — radiation: dynamics — supernovae: general

1 Introduction

Neutrino transport is an essential ingredient for numerical simulations to clarify the theory of core-collapse 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 proto-neutron 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 core-collapse supernovae (CCSNe) in the context of the neutrino-heating 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 space-time 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 neutrino-driven convection (e.g., Bethe (1990); Herant et al. (1994); Burrows et al. (1995); Janka & Müller (1996); Müller & Janka (1997)) and the standing-accretion-shock 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 multi-dimensional (multi-D) nature.111Recently multi-dimensionalities 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 multi-D 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 non-spherical (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 exa-scale computing platforms in the next decade(s) to come (see discussions in Kotake et al. (2012a))222We here mean the feasibility of 7D-GR 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 spherically-symmetric (1D) simulations where the GR Boltzmann transport is coupled to 1D-GR 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 full-fledged 1D simulations fail to produce explosions except for super-AGB stars at the low-mass end (Kitaura et al., 2006). In the context of the full Boltzmann calculations, Livne et al. (2004) were the first to perform two-dimensional (2D) multi-angle (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 multi-angle 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 multi-angle transport, Sumiyoshi & Yamada (2012) were the first to develop the fully multi-angle 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 core-collapse 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 multi-group 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 flux-limiter 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 free-streaming 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 post-Newtonian 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 self-consistently 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 multi-D simulations of CCSNe is defined by multi-group (spectral) neutrino hydrodynamics simulations. More severe approximations include gray transport (Fryer et al., 1999; Scheck et al., 2006) or the light-bulb 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 multi-D instabilities and the MHD mechanism of CCSNe.

In addition to the multi-D and multi-angle/truncated neutrino transport, the accurate treatment of GR is highly ranked among the to-do lists towards the ultimate CCSN simulations. In most of previous multi-D models with multi-group 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 core-collapse 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 neutrino-driven explosions needs to be investigated by multi-D GR simulations with more sophisticated neutrino transport scheme.333It 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 multi-D simulations.

In this paper, we present a new 3D-GR radiation hydrodynamics code that is meant to apply for stellar core-collapse simulations. The spacetime treatment is based on the Arnowitt-Deser-Misner 3+1 formalism and we employ the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism (Shibata & Nakamura, 1995; Baumgarte & Shapiro, 1999) to evolve the metric variables. The full GR radiation-hydrodynamics equations are evolved in a conservative form, in which we solve the energy-dependent 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 neutrino-matter interaction terms, we employ a baseline set of weak interactions as given in Bruenn (1985) and Rampp & Janka (2002) where nucleon-nucleon 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 self-consistent 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 long-term evolution in full 3D until the neutrino-driven 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 pseudo-1D 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 operator-splitting manner, but the system evolves self-consistently 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 neutrino-matter interaction based on the neutrino leakage scheme, on the other hand we now solve the spectral neutrino transport where the source terms are treated self-consistently 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 stress-energy tensor is expressed as


where and is the stress-energy tensor of fluid and energy-dependent neutrino radiation field, respectively. Note in the above equation, summation is taken for all species of neutrinos () with representing heavy-lepton neutrinos (i.e. and their anti-particles), 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


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




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


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




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 (7-8) to find the converged solution of .

The hydrodynamic equations are written in a conservative form as,


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 radiation-hydrodynamics 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


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




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


the equation of total lepton number conservation becomes


Here the total lepton fraction is defined by


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


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 energy-momentum 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


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


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,


where we take a simple closure relation in the diffusion limit. Inserting this into the left hand side of Eq.(4), one can get


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.5-0.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


in the slow motion limit, when we divide Eq. (29) by and integrate it in energy space, it can be summarized as


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 radiation-hydrodynamics variables444We omit our numerical method to evolve the space-time 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


where denotes conservative variables


In Eq.(32), , , and denote advection term in space, advection term in momentum space, gravitational source and neutrino-matter 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,


for the explicit part and


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 operator-splitting 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 neutrino-matter 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 Courant-Friedrichs-Lewy 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 weak-interaction 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 high-resolution-shock-capturing scheme and utilize the HLL (Harten-Lax-van 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


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




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 radiation-hydrodynamical 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.

is a modification parameter to fit the numerical flux to the diffusion limit, which we take as


where is the total opacity and is the grid width (Audit et al., 2002; O’Connor & Ott, 2013).

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,


and for energy and momentum conservations,


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 cell-surfaced 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,


As in Müller et al. (2010), we can get the solution of Eq.(42) as


where is a weighting function and is expressed as


In this study, we used a “Harmonic” interpolation () for the energy density () as


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));


where , and are expressed in terms of the primitive variables


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 Newton-Raphson iteration scheme with the inversion of the following matrix until a sufficient convergence is achieved;


for with the initial condition . As for the convergence criteria for the Newton-Raphson iteration, we monitor


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.

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

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

  3. Implicit Update. Finally, quantities in the fractional time-step () are implicitly updated to those in the -th time-step () by the Newton-Raphson 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.

Figure 1: A flowchart to visualize how to update variables (G, Q, and P) from -th to -th time step (see text).

Here we shortly summarize some technical details to achieve a sufficient lepton number conservation in practical core-collapse simulations. As we have already mentioned in Sec. 2.2, neutrino number conservation is formally satisfied by solving the energy-momentum 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)


as another criterion to exit the Newton-Raphson iteration. Because of this additional criterion, the lepton number conservation is also satisfied when the Newton-Raphson 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 Newton-Raphson 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


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 Newton-Raphson 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 ad-hoc 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


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.

The analytical solution of the zeroth and first order radiation momenta profiles at time (Pons et al., 2000) is given as


respectively. From Figure 2, it can be seen that our code can reproduce the analytical results quite well.

Figure 2: Diffusion tests with () in upper two panels and with () in bottom twos. Our results (crosses) for the evolutions of the energy density (left panels) and radial energy flux (right panels) are compared with the analytical ones (solid lines). Note that the simulations start at and 200 for models with and , respectively. The model with is performed with two nested level structure and we plot the spatial profile of in the bottom-left panel for reference. Note that we reduce number of plots for our results (crosses) to avoid messy overplots.

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


with and .

Figure 3: Shadow test at different time slices 3.5, 7.5 and 15.5 from top to bottom. Color contour represents radiation energy density in logarithmic scale. The absorbing region is expressed by white-dashed line and the expected radiation shock front (i.e., ) is denoted by a vertical green line.

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 white-dashed 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 point-like 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