Stellar core collapse in full general relativity with microphysics – Formulation and Spheical collapse test –

# Stellar core collapse in full general relativity with microphysics – Formulation and Spheical collapse test –

Yuichiro Sekiguchi

## 1 Introduction

### 1.1 Motivation

Gravitational collapse of massive stellar core to a neutron star or a black hole and the associated supernova explosion are one of the important and interesting events in the universe. From observational view point, they are among the most energetic events in astrophysics, producing a wide variety of observable signatures, namely, electromagnetic radiation, neutrinos, and gravitational radiation.

Most of the energy liberated in the collapse is eventually carried away by neutrinos from the system. The total energy of neutrinos emitted is several times ergs, where and are the mass and radius of the neutron star. Observations of gravitational collapse by neutrino detectors will provide important information of the deep inside of the core, because neutrinos can propagate from the central regions of the stellar core almost freely due to their small cross-sections with matters. Electromagnetic radiation, by contrast, interacts strongly with matters and thus gives information of collapse coming only from lower-density regions near the surface of the star. Bursts of neutrinos were first detected simultaneously by the Kamiokande and Irvine-Michigan-Brookhaven facilities in the supernova SN1987A, which occurred on February 23, 1987 in the Large Magellanic Cloud (for a review, see Ref. ?). Future detection of neutrinos will provide a direct clue to reveal the physical ingredient for the supernova explosion mechanism.

Gravitational wave astronomy will start in this decade. The first generation of ground-based interferometric detectors (LIGO, VIRGO, GEO600, TAMA300) are now in the scientific search for gravitational waves. Stellar core collapse is one of the important sources for these observatories. Observations of gravitational collapse by gravitational-wave detectors will provide unique information, complementary to that derived from electromagnetic and neutrino detectors, because gravitational waves can propagate from the innermost regions of a progenitor star to the detectors without attenuation by matters. Combination of the signatures of neutrinos and gravitational waves will provide much information about processes of the core collapse and ultimately, the physics that governs the stellar core collapse.

To obtain physically valuable information from these observations, it is necessary to connect the observed data and the physics behind it. For this purpose, a numerical simulation is the unique approach. However, simulating the stellar core collapse is one of the challenging problems because a rich diversity of physics has to be taken into account. All four known forces in nature are involved and play important roles during the collapse. General relativistic gravity plays an essential role in formation of a black hole and a neutron star. The weak interactions govern energy and lepton-number losses of the system. In particular, neutrinos transport most of the energy released during the collapse to the outside of the system. The electromagnetic and strong interactions determine the ingredient of the collapsing matter and the thermodynamical properties of the dense matter. Strong magnetic fields, if they are present, would modify the dynamics of the collapse, subsequent supernova explosion, and evolution of proto-neutron stars.

Due to these complexities, the explosion mechanism of core-collapse supernovae has not been fully understood in spite of the elaborate effort in the past about 40 years. Recent numerical studies have clarified that on the assumption of the spherical symmetry, the explosion does not succeed for the iron core collapse with the currently most elaborate input physics (neutrino interactions, neutrino transfer, and equation of states of the dense matter) on the basis of the standard “neutrino heating mechanism” (but see Ref. ? for successful explosion in O-Ne-Mg core collapse). To increase the neutrino-heating efficiency, a wide variety of multi-dimensional effects have been explored (for recent reviews, see e.g., Refs. ?), ?) and also Refs. ? and ? for simulations where successful explosions are obtained). However, it has not been completely clarified yet whether the increase of the heating efficiency due to such multi-dimensional effects suffices for yielding successful explosion, because the explosion energy resulting from these works is too low ergs.

Similarly, accurate predictions of gravitational waveforms are still hampered by the facts that reliable estimates of waveforms require a general relativistic treatment, and that appropriate treatments of microphysics such as nuclear equation of state (EOS), the electron capture, and neutrino emissions and transfers. Indeed, previous estimates of waveforms have relied either on Newtonian simulations with including microphysics to some extent, or general relativistic simulations with simplified microphysics. Recently, gravitational waveforms emitted in the rotating core collapse were derived by multidimensional simulations in general relativistic frameworks adopting a finite-temperature nuclear EOS and the electron capture. In their studies, however, the electron capture rate was not calculated in a self-consistent manner. Instead, they adopted a simplified prescription proposed in Ref. ? which is based on the result of a spherically symmetric simulation. However, it is not clear whether this treatment is justified for non-spherical collapse or not. Moreover, they did not take account of emission processes of neutrinos. More sophisticated simulations including microphysics are required to make accurate predictions of gravitational waveforms.

The gravitational collapse of massive star is also the primary mechanism of black hole formation. Understanding the process of black hole formation is one of the most important issues in the theory of the stellar core collapse. A wide variety of recent observations have shown that black holes actually exist in the universe (e.g., see Ref. ?), and so far, about 20 stellar-mass black holes for which the mass is determined within a fairly small error have been observed in binary systems of our Galaxy and the Large Magellanic Clouds. The formation of a black hole through the gravitational collapse is a highly nonlinear and dynamical phenomenon, and therefore, numerical simulation in full general relativity is the unique approach to this problem. In spherical symmetry, fully general relativistic simulations of stellar core collapse to a black hole have been performed in a state-of-the-art manner, i.e., employing realistic EOSs, implementing microphysical processes, and the Boltzmann transfer of neutrinos. In the multidimensional case, by contrast, simulations only with simplified microphysics have been performed. Because multidimensional effects such as rotation and convection are likely to play an important role, multidimensional simulations in full general relativity employing a realistic EOS and detailed microphysics are necessary for clarifying the formation process of black holes.

Furthermore, recent observations have found the spectroscopic connections between several SNe and long gamma-ray bursts (GRBs) and clarified that some of long GRBs are associated with the collapse of massive stars. Supported by these observations, the collapsar model is currently one of promising models for the central engine of GRBs. In this model, one assumes that the central engine of the long GRBs is composed of a rotating black hole and a hot, massive accretion disk. Such a system may be formed as a result of the collapse of rapidly rotating massive core. In this model, one of the promising processes of the energy-deposition to form a GRB fireball is the pair annihilation of neutrinos emitted from the hot, massive disk (). The collapsar model requires the progenitor core to be rotating rapidly enough that the massive accretion disk can be formed around the black hole. Recent general relativistic numerical analyses have shown that if a progenitor of the collapse is massive and the angular momentum is large enough, a black hole surrounded by a massive disk will be formed. However, the formation mechanism of such system has not been clarified in detail. These also enhance the importance of exploring the stellar core collapse to a black hole taking account of microphysical processes.

As reviewed above, multidimensional simulations of stellar collapse in full general relativity including microphysics is currently one of the most required subjects in theoretical astrophysics. However, there has been no multidimensional code in full general relativity that self-consistently includes microphysics such as realistic EOS, electron capture, and neutrino emission. There have only existed fully general relativistic codes in spherical symmetry or Newtonian codes in multidimension. We have developed a fully general relativistic multidimensional code including a finite-temperature nuclear EOS, self-consistent treatment of the electron capture, and a simplified treatment of neutrino emission for the first time. In this code, by contrast with the previous ones, the electron capture rate is treated in a self-consistent manner and the neutrino cooling is taken into account for the first time. Because it is not currently feasible to fully solve the neutrino transfer equations in the framework of general relativity in multidimension because of restrictions of computational resources, it will be reasonable to take some approximation for the transfer equations at the current status. In this paper, the so-called neutrino leakage scheme is adopted as an approximate treatment of neutrino cooling, and a general relativistic version of the leakage scheme is developed.

### 1.2 The leakage schemes

The neutrino leakage schemes as an approximate method for the neutrino cooling have a well-established history (e.g. Ref. ?)). The basic concept of the original neutrino leakage schemes is to treat the following two regions in the system separately: one is the region where the diffusion timescale of neutrinos is longer than the dynamical timescale, and hence, neutrinos are ’trapped’ (neutrino-trapped region); the other is the region where the diffusion timescale is shorter than the dynamical timescale, and hence, neutrinos stream out freely out of the system (free-streaming region). The idea of treating the diffusion region separately has been applied to more advanced methods for the neutrino transfer (see e.g., Ref. ? and references therein).

Then, electron neutrinos and electron anti-neutrinos in the neutrino-trapped region are assumed to be in the -equilibrium state. The net local rates of lepton-number and energy exchange with matters are set to be zero in the neutrino-trapped region. To treat diffusive emission of neutrinos leaking out of the neutrino-trapped region, simple phenomenological source terms based on the diffusion theory are introduced. In the free-streaming region, on the other hand, it is assumed that neutrinos escape from the system without interacting with matters. Therefore, neutrinos carry the lepton number and the energy according to the local weak-interaction rates. Note that the neutrino fractions are not solved in the original version of the leakage scheme: Only the total lepton fraction (from which the neutrino fractions are calculated under the -equilibrium condition) is necessary in the neutrino-trapped region, and the neutrino fractions are set to be zero in the free-streaming region. As a result, neutrino quantities and the electron fraction are discontinuous at the boundary the neutrino-trapped and free-streaming regions.

The boundary was given by hand as a single ’neutrino-trapping’ density () without calculating the optical depths of neutrinos in the previous studies. However, the location at which the neutrino trapping occurs in fact depends strongly on the neutrino energies () as , and hence, there are different neutrino-trapping densities for different neutrino energies. In the previous leakage schemes, on the other hand, all neutrinos were emitted in one moment irrespective of their energy. Consequently in the case of the so-called neutrino burst emission (e.g., Ref. ?)), for example, the duration in which the neutrinos are emitted was shortened and the peak luminosity at the burst was overestimated. The dependence of the neutrino-trapping densities and the neutrino diffusion rates on the neutrino energies are approximately taken into account in the recent simulations of mergers of binary neutron star. However, the lepton-number conservation equations for neutrinos are not solved, which is important to estimate the phase space blocking due to neutrinos.

Transfer equations of neutrinos are not solved in the leakage schemes. Therefore, the leakage schemes cannot treat non-local interactions among the neutrinos and matters. For example, the so-called neutrino heating and the neutrino pair annihilation cannot be treated in the leakage scheme. Nevertheless, we believe a detailed general relativistic leakage scheme presented in this paper to be a valuable approach because even by this approximated approach it is possible to incorporate the effects of neutrinos semi-quantitatively as shown in this paper. Also, the neutrino leakage scheme is an appropriate method for studying a number of phenomena for which the neutrino heating and neutrino transfer are expected to be not very important, e.g., prompt formation of a black hole and compact binary mergers. Both of these phenomena are the targets of the present code.

A first attempt towards a general relativistic leakage scheme was done in the previous study. In that study, not the region of the system but the energy momentum tensor of neutrinos was decomposed into two parts; ’trapped-neutrino’ and ’streaming-neutrino’ parts. However the source terms of hydrodynamic and lepton-number-conservation equations were determined using the single neutrino-trapping density as in the case of the previous leakage schemes. In this paper, we develop a new code implementing the microphysical processes in the general relativistic framework based on the previous study. As an application of the code, we perform simulations of stellar core collapse.

A lot of improved ingredients are installed into the present code: (1) The dependence of the neutrino diffusion rates on the neutrino energies are approximately taken into account following the recent study with detailed cross sections, instead of adopting the single neutrino-trapping density (see Appendix C). (2) The lepton-number conservation equations for neutrinos are solved to calculate self-consistently the chemical potentials of neutrinos. Then, the blocking effects due to the presence of neutrinos and the -equilibrium condition can be taken into account more accurately (see §3). (3) A stable explicit method for solving the equations of hydrodynamics, the lepton-number conservations, and neutrinos are developed. Such a special procedure is necessary because the characteristic timescale of the weak-interaction processes (hereafter referred to as the WP timescale ) is much shorter than the dynamical timescale in hot, dense matter regions. Note that in the previous leakage schemes, the -equilibrium was assumed to be achieved in such regions (i.e. ) and no such special treatments are required. See §2 for further discussions and §3 for details of the method. (4) The electron capture rate are calculated in a detailed manner including effects of the so-called thermal unblocking (see Appendix A).

The paper is organized as follows. First, issues in implementation of weak interactions and neutrino cooling in full general relativistic simulation is briefly summarized in §2. Then, framework of the general relativistic leakage scheme is described in detail in § 3. In § 4, EOSs employed in this paper are described in some details. Basic equations and numerical methods of the simulations are described in § 5. Numerical results obtained in this paper are shown in § 6. We devote § 7 to a summary and discussions. In appendices, details of the microphysics adopted in the present paper are summarized for the purpose of convenience. Throughout the paper, the geometrical unit is used otherwise stated.

## 2 Issues in implementation of weak interactions and neutrino cooling in fully general relativistic simulation

Because the characteristic timescale of the weak-interaction processes (the WP timescale ) is much shorter than the dynamical timescale in hot dense matters, the explicit numerical treatments of the weak interactions are computationally expensive in simple methods, as noted in the previous pioneering work by Bruenn: A very short timestep ( ) will be required to solve the equations explicitly.

The net rates of lepton-number and energy exchanges between matters and neutrinos may not be large, and consequently, an effective timescale may not be as short as the dynamical timescale. However, this does not immediately imply that one can solve the equations explicitly without employing any prescription. For example, the achievement of -equilibrium, where is the consequence of the cancellation of two very large weak interaction processes (the electron and the electron-neutrino captures, see Eq. (3.2)) and of the action of the phase space blocking. Note that the weak interaction processes depend enormously both on the temperature and the lepton chemical potentials. Therefore, small error in the evaluation of the temperature and a small deviation from the -equilibrium due to small error in calculation of the lepton chemical potentials will result in huge error. Then, stiff source terms appear and explicit numerical evolution often becomes unstable. Indeed, we found that a straightforward, explicit solution of the equations did not work.

In the following of this section, we describe issues of implementation of weak interactions and neutrino cooling into the hydrodynamic equations in the conservative schemes in fully general relativistic simulations. Fiest, we illustrate that in the Newtonian framework, the equations may be solved implicitly in a relatively simple manner (see also Refs. ? and ?) and references therein). The equations of hydrodynamics, lepton-number conservations, and neutrino processes are schematically written as,

 ˙ρ = 0, (2.1) ˙vi = Svi(ρ,Ye,T,Qν), (2.2) ˙Ye = SYe(ρ,Ye,T,Qν), (2.3) ˙e = Se(ρ,Ye,T,Qν), (2.4) ˙Qν = SQν(ρ,Ye,T,Qν), (2.5)

where is the rest-mass density, is the velocity, is the electron fraction, is the (internal) energy of matter, is the temperature, and stands for the relevant neutrino quantities. We here omit the transport terms. ’s in the right-hand side stand for the relevant source terms. Comparing the quantities in the left-hand-side and the argument quantities in the source terms, only the relation between and is nontrivial. Usually, EOSs employed in the simulation are tabularized, and one dimensional search over the EOS table is required to solve them. Due to the relatively simple relations between the quantities to be evolved and the argument quantities, the above equations may be solved implicitly in a straightforward (although complicated) manner.

In the relativistic framework, the situation becomes much more complicated in conservative schemes, because the Lorentz factor () is coupled with rest-mass density and the energy density (see Eqs. (5.12) and (5.16) where is used instead of ), and because the specific enthalpy () is coupled with the momentum (see Eq. (5.14)).

It should be addressed that the previous fully general relativistic works in the spherical symmetry are based on the so-called Misner-Sharp coordinates. There are no such complicated couplings in these Lagrangian coordinates. Accordingly, the equations may be solved essentially in the same manner as in the Newtonian framework. Because no such simple Lagrangian coordinates are known in the multidimensional case, the complicated couplings inevitably appear in the relativistic framework.

Omitting the factors associated with the geometric variables (which are usually updated before solving the hydrodynamics equations) and the transport terms, the equations to be solved in the general relativistic framework are schematically written as,

 ˙ρ∗(ρ,\mathchar28928\relax) = 0, (2.6) ˙^ui(ui,h)=˙^ui(ui,ρ,Ye,T) = S^ui(ρ,Ye,T,Qν,\mathchar28928\relax), (2.7) ˙Ye = SYe(ρ,Ye,T,Qν,\mathchar28928\relax), (2.8) ˙^e(ρ,Ye,T,\mathchar28928\relax) = S^e(ρ,Ye,T,Qν,\mathchar28928\relax), (2.9) ˙Qν = SQν(ρ,Ye,T,Qν,\mathchar28928\relax), (2.10)

where is a weighted density, is a weighted four velocity, is a weighted energy density (see § 5.2 for the definition of these variables). The Lorentz factor is calculated by solving the normalization condition , which is rather complicated nonlinear equation schematically written as

 fnormalization(^ui,\mathchar28928\relax)=fnormalization(ui,ρ,Ye,T,\mathchar28928\relax)=0. (2.11)

The accurate calculation of the Lorentz factor and the accurate solution of the normalization condition are very important in the numerical relativistic hydrodynamics.

Now, it is obvious that the argument quantities in the source terms are not simply related with the evolved quantities in the left-hand-side of Eqs. (2)–(2). Solving the equations implicitly is not as straightforward as in the Newtonian case and no successful formulations have been developed. Moreover it might be not clear whether a convergent solution can be stably obtained numerically or not, because they are simultaneous nonlinear equations. Therefore, it may be not a poor choice to adopt an alternative approach in which the equations are solved explicitly with some approximations as described in the next section*)*)*)It should be stated that the implicit schemes are also approximated ones because a short WP timescale associated with the weak interaction is not fully resolved..

## 3 General relativistic neutrino leakage scheme

In this section, we describe a method for approximately solving hydrodynamic equations coupled with neutrino radiation in an explicit manner. As described in the previous section, because of the relation in the hot dense matter regions, the source terms in the equations become too stiff for the equations to be solved explicitly in the straightforward manner. The characteristic timescale of leakage of neutrinos from the system , by contrast, is much longer than in the hot dense matter region. Rather, where is the characteristic length scale of the system. On the other hand, is comparable to in the free-streaming regions but is longer than or comparable with there. All these facts imply that the WP timescale does not directly determine the evolution of the system but the leakage timescale does. Using this fact, we approximate some of original equations and reformulate them so that the source terms are to be characterized by the leakage timescale .

### 3.1 Decomposition of neutrino energy-momentum tensor

The basic equations of general relativistic hydrodynamics with neutrinos are

 ∇α(TTotal)αβ=∇α[(TF)αβ+(Tν)αβ]=0, (3.1)

where is the total energy-momentum tensor, and and are the energy-momentum tensor of fluids and neutrinos, respectively. Equation (3.1) can be written as

 ∇α(TF)αβ = Qβ, (3.2) ∇α(Tν)αβ = −Qβ, (3.3)

where the source term is regarded as the local production rate of neutrinos through the weak processes.

Now, the problem is that the source term becomes too stiff to solve explicitly in hot dense matter regions where . To overcome this situation, the following procedures are adopted.

First, it is assumed that the energy-momentum tensor of neutrinos are be decomposed into ’trapped-neutrino’ () and ’streaming-neutrino’ () parts as,

 (Tν)αβ=(Tν,T)αβ+(Tν,S)αβ. (3.4)

Here, the trapped-neutrinos phenomenologically represent neutrinos which interact sufficiently frequently with matter and are thermalized. On the other hand, the streaming-neutrino part describes a phenomenological flow of neutrinos streaming out of the system (see also Ref. ? where a more sophisticate method in terms of the distribution function is adopted in the Newtonian framework).

Second, the locally produced neutrinos are assumed to leak out to be the streaming-neutrinos with a leakage rate :

 ∇β(Tν,S)βα=Qleakα. (3.5)

Then, the equation of the trapped-neutrino part is

 ∇β(Tν,T)βα=Qα−Qleakα. (3.6)

Third, the trapped-neutrino part is combined with the fluid part as

 Tαβ≡(TF)αβ+(Tν,T)αβ, (3.7)

and Eqs. (3.1) and (3.1) are combined to give

 ∇βTβα=−Qleakα. (3.8)

Thus, the equations to be solved are changed from Eqs. (3.1) and (3.1) to Eqs. (3.1) and (3.1). Note that the new equations only include the source term which is characterized by the leakage timescale . Definition of will be given in § 3.3.

The energy-momentum tensor of the fluid and trapped-neutrino parts () is treated as that of the perfect fluid,

 Tαβ=(ρ+ρε+P)uαuβ+Pgαβ, (3.9)

where and are the rest mass density and the 4-velocity. The specific internal energy density () and the pressure () are the sum of contributions from the baryons (free protons, free neutrons, -particles, and heavy nuclei), leptons (electrons, positrons, and trapped-neutrinos), and the photons as,

 P = PB+Pe+Pν+Pph, (3.10) ε = εB+εe+εν+εph, (3.11)

where subscripts ’’, ’’, ’’, and ’’ denote the components of the baryons, electrons and positrons, photons, and trapped-neutrinos, respectively.

The streaming-neutrino part, on the other hand, is set to be a general form of

 (Tν,S)αβ=Enαnβ+Fαnβ+Fβnα+Pαβ, (3.12)

where . In order to close the system, we need an explicit expression of . In this paper, we adopt a simple form with . This approximation may work well in high density regions but will violate in low density regions. However, the violation will not affect the dynamics because the total amount of streaming-neutrinos emitted in low density regions will be small. Of course, a more sophisticated treatment will be necessary in a future study.

### 3.2 The lepton-number conservation equations

The conservation equations of the lepton fractions are written schematically as

 dYedt=−γe, (3.13) dYνedt=γνe, (3.14) dY¯νedt=γ¯νe, (3.15) dYνxdt=γνx, (3.16)

where , , , and denote the electron fraction, the electron neutrino fraction, the electron anti-neutrino fraction, and and neutrino and anti-neutrino fractions, respectively. We note that in the previous simulations based on the leakage schemes, the neutrino fractions were not solved.

The source terms of neutrino fractions can be written, on the basis of the present leakage scheme, as

 γνe=γlocalνe−γleakνe, (3.17) γ¯νe=γlocal¯νe−γleak¯νe, (3.18) γνx=γlocalνx−γleakνx, (3.19)

where ’s and ’s are the local production and the leakage rates of each neutrino, respectively (see § 3.3). Note that only the trapped-neutrinos are responsible for the neutrino fractions. Assuming that the trapped neutrinos are thermalized and the distribution function is the equilibrium Fermi-Dirac one, the chemical potentials of neutrinos can be calculated from the neutrino fractions. Then the thermodynamical quantities of neutrinos can be also calculated.

The source term of the electron fraction conservation is

 γe=γlocalνe−γlocal¯νe. (3.20)

Because  ’s are characterized by the WP timescale , some procedures are necessary to solve the lepton conservation equations explicitly. The following simple procedures are proposed to solve the equations stably.

First, in each timestep , the conservation equation of the total lepton fraction (),

 dYldt=−γl, (3.21)

is solved together with the conservation equation of , Eq. (3.2), in advance of solving whole of the lepton conservation equations (Eqs. (3.2) – (3.2)). Note that the source term is characterized by the leakage timescale so that this equation can be solved explicitly in the hydrodynamic timescale. Then, assuming that the -equilibrium is achieved, values of the lepton fractions in the -equilibrium (, , and ) are calculated from the evolved .

Second, regarding and as the maximum allowed values of the neutrino fractions in the next timestep , the source terms are limited so that ’s in the timestep cannot exceed ’s. Then, the whole of the lepton conservation equations (Eqs. (3.2) – (3.2)) are solved explicitly using these limiters.

Third, the following conditions are checked

 μp+μe<μn+μνe, (3.22) μn−μe<μp+μ¯νe, (3.23)

where , , , and are the chemical potentials of protons, neutrons, electrons, electron neutrinos, and electron anti-neutrinos, respectively. If both conditions are satisfied, the values of the lepton fractions in the timestep are set to be those in the -equilibrium value; , , and . On the other hand, if either or both conditions are not satisfied, the lepton fractions in the timestep is set to be those obtained by solving whole of the lepton-number conservation equations.

A limiter for the evolution of may be also necessary in the case where the pair processes are dominant, for example, in the simulations for collapse of population III stellar core. In this case, the value of at the pair equilibrium (i.e. at ), may be used to limit the source term.

### 3.3 Definition of leakage rates

In this subsection the definitions of the leakage rates and are presented. Because may be regarded as the emissivity of neutrinos measured in the fluid rest frame, is defined as

 Qleakα=Qleakνuα. (3.24)

The leakage rates and are assumed to satisfy the following properties.

1. The leakage rates approach the local rates and in the low density, transparent region.

2. The leakage rates approach the diffusion rates and in the high density, opaque region.

3. The above two limits should be connected smoothly.

Here, the local rates can be calculated based on the theory of weak interactions and the diffusion rates can be determined based on the diffusion theory (see appendices for the definition of the local and diffusion rates adopted in this paper). There will be several prescriptions to satisfy the requirement (iii). In this paper, the leakage rates are defined by

 Qleakν=(1−e−bτν)Qdiffν+e−bτνQlocalν, (3.25) γleakν=(1−e−bτν)γdiffν+e−bτνγlocalν, (3.26)

where is the optical depth of neutrinos and is a parameter which is typically set as . The optical depth can be computed from the cross sections in a standard manner.

In the present implementation, it is not necessary to artificially divide the system into neutrino-trapped and free-streaming regions by the single neutrino-trapping density. Therefore there is no discontinuous boundary which existed in the previous leakage schemes.

As the local production reactions of neutrinos, the electron and positron captures ( and ), the electron-positron pair annihilation ( for electron-type neutrinos and for the other type), the plasmon decays ( and ), and the Bremsstrahlung processes ( and ) are considered in this paper. Then, the local rates for the neutrino fractions are

 γlocalνe=γecνe+γpairνe¯νe+γplasνe¯νe+γBremsνe¯νe, (3.27) (3.28) γlocalνx=γpairνx¯νx+γplasνx¯νx+γBremsνx¯νx. (3.29)

Similarly, the local neutrino energy emission rate is given by

 Qlocalν=Qecνe+Qpc¯νe + 2(Qpairνe¯νe+Qplasνe¯νe+QBremsνe¯νe) (3.30) + 4(Qpairνx¯νx+Qplasνx¯νx+QBremsνx¯νx) .

The explicit forms of the local rates in Eqs. (3.3)–(3.30) will be found in Appendices A and B.

We follow the recent work by Rosswog and Liebendörfer for the diffusive neutrino emission rates and in Eqs (3.3) and (3.3). The explicit forms of and are presented in Appendix C.

## 4 Equation of state

In this section we summarize details of EOSs adopted in our current code.

### 4.1 Baryons

In the present version of our code, we employ an EOS by Shen et al., which is derived by the relativistic mean field theory based on the relativistic Brückner-Hartree-Fock theory. The so-called parameter set TM1 is adopted to reproduce characteristic properties of heavy nuclei. The maximum mass of a cold spherical neutron star in this EOS is much larger than the canonical neutron star mass as . The framework of the relativistic mean field theory is extended with the Thomas-Fermi spherical cell model approximation to describe not only the homogeneous matter but also an inhomogeneous one.

The thermodynamical quantities of dense matter at various sets of are calculated to construct the numerical data table for simulation. The table covers a wide range of density g/cm, electron fraction , and temperature MeV, which are required for supernova simulation. It should be noted that the causality is guaranteed to be satisfied in this framework, whereas the sound velocity sometimes exceeds the velocity of the light in the non-relativistic framework, e.g., in the EOS by Lattimer and Swesty. This is one of the benefits of the relativistic EOS.

Although we employ the nuclear EOS by Shen et al. in this work, it is easy to replace the EOS. In the future we plan to implement other EOSs such as a hyperonic matter EOS.

Because the table of the original EOS by Shen et al. does not include the thermodynamical quantities of the leptons (electrons, positrons, and neutrinos if necessary) and photons, one has to consistently include them to the table.

### 4.2 Electrons and Positrons

To consistently calculate the pressure and the internal energy of the electron and positron, the charge neutrality condition should be solved to determine the electron chemical potential for each value of the baryon rest-mass density and the temperature in the EOS table. Namely, it is required to solve the equation

 ne(μe,T)≡n−−n+=ρYemu (4.1)

in terms of for given values of , , and . Here, MeV is the atomic mass unit, and and are the total number densities (i.e., including the electron-positron pair) of the electrons and positrons, respectively.

Assuming that the electrons obey the Fermi-Dirac distribution (which is derived under the assumption of thermodynamic equilibrium), the number density (), the pressure (), and the internal energy density () of the electron are written as

 n− = 1π2ℏ3∫∞0p2dpexp[−ηe+~ϵ/kBT]+1, (4.2) P− = 1π2ℏ3∫∞0p3(∂~ϵ/∂p)dpexp[−ηe+~ϵ/kBT]+1, (4.3) u− = 1π2ℏ3∫∞0p2~ϵdpexp[−ηe+~ϵ/kBT]+1. (4.4)

Here , , and are the Planck’s constant, the Boltzmann’s constant and the so-called degeneracy parameter. is the kinetic energy of a free electron. If we choose the zero point of our energy scale for electrons at , we have to assign a total energy of to a free positron. Thus the number density (), the pressure (), and the internal energy density () of positrons are given by

 n+ = 1π2ℏ3∫∞0p2dpexp[−η++~ϵ+/kBT]+1, (4.5) P+ = 1π2ℏ3∫∞0p3(∂~ϵ+/∂p)dpexp[−η++~ϵ+/kBT]+1, (4.6) u+ = 1π2ℏ3∫∞0p2(~ϵ+2mec2)dpexp[−η++~ϵ+/kBT]+1, (4.7)

where is the degeneracy parameter of the positrons.

### 4.3 Photons

The pressure and the specific internal energy density of photons are given by

 Pr=arT43,  εr=arT4ρ, (4.8)

where is the radiation constant and is the velocity of light.

### 4.4 Trapped neutrinos

In this paper, the trapped-neutrinos are assumed to interact sufficiently frequently with matter that be thermalized. Therefore they are described as ideal Fermi gases with the matter temperature. Then, from the neutrino fractions , the chemical potentials of neutrinos are calculated by solving

 Yν=Yν(μν,T). (4.9)

Using the chemical potentials, , and the matter temperature, the pressure and the internal energy of the trapped-neutrinos are calculated in the same manner as for electrons.

### 4.5 The sound velocity

In the high-resolution shock-capturing scheme for hydrodynamics, we in general need to evaluate the sound velocity ,

 c2s=1h[∂P∂ρ∣∣∣ϵ+Pρ∂P∂ϵ∣∣∣ρ]. (4.10)

The derivatives of the pressure are calculated by

 ∂P∂ρ∣∣∣ϵ = ∑i=B,e,r,ν⎡⎣∂Pi∂ρ∣∣∣T−∂Pi∂T∣∣∣ρ(∑j=B,e,r,ν∂ϵj∂ρ∣∣∣T)(∑k=B,e,r,ν∂ϵk∂T∣∣∣ρ)−1⎤⎦, (4.11) ∂P∂ϵ∣∣∣ρ = (∑i=B,e,r,ν∂Pi∂T∣∣∣ρ)(∑j=B,e,r,ν∂ϵj∂T∣∣∣ρ)−1, (4.12)

where ’’, ’’, ’’ and ’’ in the sum denote the baryon, electron, photons, and neutrino quantities.

The derivatives for the baryon parts are evaluated by taking a finite difference of the table data. On the other hand, the derivatives for the electron parts can be evaluated semi-analytically. Because there are in general the phase transition regions in an EOS table for baryons and moreover the EOS may contain some non-smooth spiky structures, careful treatments are necessary when evaluating the derivatives of thermodynamical quantities. In the present EOS table, the derivatives are carefully evaluated so that there are no spiky behaviors in the resulting sound velocities.

## 5 Basic equations and Numerical methods

### 5.1 Einstein’s equation and gauge conditions

The standard variables in the 3+1 decomposition are the three-dimensional metric and the extrinsic curvature on the three-dimensional hypersurface defined by

 γμν ≡ gμν+nμnν, (5.1) Kμν ≡ −12\@fontswitchLto0.0pt−−nγμν, (5.2)

where is the spacetime metric, is the unit normal to a three-dimensional hypersurface, and is the Lie derivative with respect to the unit normal . Then we can write the line element in the form

 ds2=−α2dt2+γij(dxi+βidt)(dxj+βjdt), (5.3)

where and are the lapse function and the shift vector which describe the gauge degree of freedom.

In the BSSN reformulation, the spatial metric is conformally decomposed as where the condition is imposed for the conformal metric . From this condition, the conformal factor is written as and . The extrinsic curvature is decomposed into the trace part and the traceless part as . The traceless part is conformally decomposed as . Thus the fundamental quantities for the evolution equation are now split into , , and . Furthermore, the auxiliary variable is introduced in the BSSN reformulation.

The basic equations to be solved are

 (∂t−βk∂k)ϕ=16(−αK+∂kβk), (5.4) (∂t−βk∂k)~γij=−2α~Aij+~γik∂jβk+~γjk∂iβk−23~γij∂kβk, (5.5) (∂t−βk∂k)K=−DkDkα+α[~Aij~Aij+13K2]+4πα(ρh+S), (5.6) +α(K~Aij−2~Aik~Ak j)+~Aik∂jβk+~Ajk∂iβk−23~Aij∂kβk −8πα(e−4ϕSij−13~γijS), (5.7) (∂t−βk∂k)Fi=−16παji +2α{fkj∂j~Aik+~Aik∂jfkj−12~Ajl∂ihjl+6~Ak i∂kϕ−23∂iK} +δjk{−2~Aij∂kα+(∂kβl)∂lhij +∂k(~γil∂jβl+~γjl∂iβl−23~γij∂lβl)}, (5.8)

where , , and are the Ricci scalar, the Ricci tensor, and the covariant derivative associated with three-dimensional metric , respectively. The matter source terms, , , and , are the projections of the stress-energy tensor with respect to and , and .

We assume the axial symmetry of the spacetime and the so-called Cartoon method is adopted to avoid problems around the coordinate singularities of the cylindrical coordinates. Except for this, the numerical schemes for solving the Einstein’s equation are essentially the same as those in Ref. ?. We use 4th-order finite difference scheme in the spatial direction and the 3rd-order Runge-Kutta scheme in the time integration. The advection terms such as are evaluated by a 4th-order upwind scheme.

As the gauge conditions for the lapse, we use the so-called slicing:

 (∂t−\@fontswitchLto0.0pt−−β)α=−2Kα. (5.9)

It is known that the slicing enables to perform a long term evolution of neutron stars as well as has strong singularity avoidance properties in the black hole spacetime.

The shift vector is determined by solving the following dynamical equation

 ∂tβk=~γkl(Fl+\mathchar28929\relaxt∂tFl). (5.10)

Here the second term in the right-hand side is necessary for numerical stability.

### 5.2 The hydrodynamic equations in leakage scheme

The basic equations for general relativistic hydrodynamics in our leakage scheme are the continuity equation, the lepton-number conservation equations, and the local conservation equation of the energy-momentum. We assume the axial symmetry of the spacetime and the hydrodynamics equations are solved in the cylindrical coordinates where . In the axisymmetric case, the hydrodynamics equations should be written in the cylindrical coordinate. On the other hand, in the Cartoon method, Einstein’s equation are solved in the plane for which , , , and the other similar relations hold for vector and tensor quantities. Taking into these facts, the hydrodynamic equations may be written using the Cartesian coordinates replacing by . In the following, we write down explicit forms of the equations for the purpose of convenience. Numerical tests for basic parts of the code of solving the hydrodynamics equations are extensively performed in Ref. ?). The equations are solved using the third-order high-resolution central scheme of Kurganov and Tadmor.

#### 5.2.1 The Continuity and lepton-number conservation equations

The continuity equation for the baryon rest mass is

 ∇α(ρuα)=0. (5.11)

As fundamental variables for numerical simulations, the following quantities are introduced: and where . Then, the continuity equation is written as

 ∂t(ρ∗)+1x∂x(ρ∗vx)+∂z(ρ∗vz)=0. (5.12)

Using the continuity equation, the lepton-number conservation equations (3.2) – (3.2) are written as

 ∂t(ρ∗YL)+1x∂x(ρ∗YLvx)+∂z(ρ∗YLvz)=ρ∗γL, (5.13)

where and are abbreviated expressions of the lepton fractions and the source terms.

#### 5.2.2 Energy-momentum conservation

As fundamental variables for numerical simulations, we define the quantities and . Then, the Euler equation , and the energy equation