Stellar core collapse in full general relativity with microphysics - Formulation and Spherical collapse test -
aa r X i v : . [ a s t r o - ph . H E ] S e p Stellar core collapse in full general relativity with microphysics– Formulation and Spheical collapse test –
Yuichiro
Sekiguchi
Division of Theoretical Astronomy, National Astronomical Observatory of Japan,Mitaka, Tokyo 181-8588, Japan
One of the longstanding issues in numerical relativity is to enable a simulation takingaccount of microphysical processes (e.g., weak interactions and neutrino cooling). We de-velop an approximate and explicit scheme in the fully general relativistic framework as afirst implementation of the microphysics toward a more realistic and sophisticated model-ing. In this paper, we describe in detail a method for implementation of a realistic equationof state, the electron capture and the neutrino cooling in a multidimensional, fully generalrelativistic code. The procedure is based on the so-called neutrino leakage scheme. To checkthe validity of the code, we perform a two dimensional (2D) simulation of spherical stel-lar core collapse. Until the convective activities set in, our results approximately agree, orat least are consistent, with those in the previous so-called state-of-the-art simulations. Inparticular, the radial profiles of thermodynamical quantities and the time evolution of theneutrino luminosities agree quantitatively. The convection is driven by negative gradients ofthe entropy per baryon and the electron fraction as in the previous 2D Newtonian simula-tions. We clarify which gradient is more responsible for the convection. Gravitational wavesfrom the convection are also calculated. We find that the characteristic frequencies of thegravitational-wave spectra are distributed for higher frequencies than those in Newtoniansimulations due to the general relativistic effects. §
1. Introduction
Motivation
Gravitational collapse of massive stellar core to a neutron star or a black holeand the associated supernova explosion are one of the important and interestingevents in the universe. From observational view point, they are among the mostenergetic 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 neu-trinos from the system. The total energy of neutrinos emitted is ≈ GM /R NS ∼ . M NS c ∼ several times 10 ergs, where M NS and R NS are the mass and radiusof the neutron star. Observations of gravitational collapse by neutrino detectorswill provide important information of the deep inside of the core, because neutrinoscan propagate from the central regions of the stellar core almost freely due to theirsmall cross-sections with matters. Electromagnetic radiation, by contrast, interactsstrongly 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 detectedsimultaneously by the Kamiokande and Irvine-Michigan-Brookhaven facilities inthe supernova SN1987A, which occurred on February 23, 1987 in the Large Magel-lanic Cloud (for a review, see Ref. 3)). Future detection of neutrinos will provide a typeset using PTP
TEX.cls h Ver.0.9 i Y. Sekiguchi direct clue to reveal the physical ingredient for the supernova explosion mechanism.Gravitational wave astronomy will start in this decade. The first generation ofground-based interferometric detectors (LIGO, VIRGO, GEO600, TAMA300 )are now in the scientific search for gravitational waves. Stellar core collapse is oneof the important sources for these observatories. Observations of gravitational col-lapse by gravitational-wave detectors will provide unique information, complemen-tary to that derived from electromagnetic and neutrino detectors, because gravita-tional waves can propagate from the innermost regions of a progenitor star to thedetectors without attenuation by matters. Combination of the signatures of neutri-nos and gravitational waves will provide much information about processes of thecore collapse and ultimately, the physics that governs the stellar core collapse.To obtain physically valuable information from these observations, it is neces-sary to connect the observed data and the physics behind it. For this purpose, anumerical simulation is the unique approach. However, simulating the stellar corecollapse is one of the challenging problems because a rich diversity of physics hasto be taken into account. All four known forces in nature are involved and playimportant roles during the collapse. General relativistic gravity plays an essentialrole in formation of a black hole and a neutron star. The weak interactions governenergy and lepton-number losses of the system. In particular, neutrinos transportmost of the energy released during the collapse to the outside of the system. Theelectromagnetic and strong interactions determine the ingredient of the collapsingmatter and the thermodynamical properties of the dense matter. Strong magneticfields, if they are present, would modify the dynamics of the collapse, subsequentsupernova explosion, and evolution of proto-neutron stars.Due to these complexities, the explosion mechanism of core-collapse supernovaehas not been fully understood in spite of the elaborate effort in the past about 40years. Recent numerical studies have clarified that on theassumption of the spherical symmetry, the explosion does not succeed for the ironcore collapse with the currently most elaborate input physics (neutrino interactions,neutrino transfer, and equation of states of the dense matter) on the basis of thestandard “neutrino heating mechanism” (but see Ref. 17) for successful explosionin O-Ne-Mg core collapse). To increase the neutrino-heating efficiency, a wide va-riety of multi-dimensional effects have been explored (for recent reviews, see e.g.,Refs. 9), 8) and also Refs. 18) and 19) for simulations where successful explosionsare obtained). However, it has not been completely clarified yet whether the increaseof the heating efficiency due to such multi-dimensional effects suffices for yieldingsuccessful explosion, because the explosion energy resulting from these works is toolow ∼ ergs.Similarly, accurate predictions of gravitational waveforms are still hampered bythe facts that reliable estimates of waveforms require a general relativistic treat-ment, and that appropriate treatments of microphysics such as nuclear equationof state (EOS), the electron capture, and neutrino emissions and transfers. Indeed,previous estimates of waveforms have relied either on Newtonian simulations withincluding microphysics to some extent, or general relativis-tic simulations with simplified microphysics. Recently, gravitational ull GR simulation with microphysics adopting a finite-temperature nu-clear EOS and the electron capture. In their studies, however, the electron capturerate was not calculated in a self-consistent manner. Instead, they adopted a simpli-fied prescription proposed in Ref. 36) which is based on the result of a sphericallysymmetric simulation. However, it is not clear whether this treatment is justifiedfor non-spherical collapse or not. Moreover, they did not take account of emissionprocesses of neutrinos. More sophisticated simulations including microphysics arerequired to make accurate predictions of gravitational waveforms.The gravitational collapse of massive star is also the primary mechanism of blackhole formation. Understanding the process of black hole formation is one of the mostimportant issues in the theory of the stellar core collapse. A wide variety of recentobservations have shown that black holes actually exist in the universe (e.g., seeRef. 37)), and so far, about 20 stellar-mass black holes for which the mass is de-termined within a fairly small error have been observed in binary systems of ourGalaxy and the Large Magellanic Clouds.
The formation of a black hole throughthe gravitational collapse is a highly nonlinear and dynamical phenomenon, andtherefore, numerical simulation in full general relativity is the unique approach tothis problem. In spherical symmetry, fully general relativistic simulations of stellarcore 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 Boltzmanntransfer of neutrinos.
In the multidimensional case, by contrast, simulationsonly with simplified microphysics have been performed.
Because multi-dimensional effects such as rotation and convection are likely to play an importantrole, multidimensional simulations in full general relativity employing a realistic EOSand detailed microphysics are necessary for clarifying the formation process of blackholes.Furthermore, recent observations have found the spectroscopic connec-tions between several SNe and long gamma-ray bursts (GRBs) and clarified thatsome of long GRBs are associated with the collapse of massive stars. Supported bythese observations, the collapsar model is currently one of promising models forthe central engine of GRBs. In this model, one assumes that the central engine ofthe long GRBs is composed of a rotating black hole and a hot, massive accretiondisk. Such a system may be formed as a result of the collapse of rapidly rotatingmassive core. In this model, one of the promising processes of the energy-depositionto form a GRB fireball is the pair annihilation of neutrinos emitted from the hot,massive disk ( ν e + ¯ ν e → e − + e + ). The collapsar model requires the progenitor coreto be rotating rapidly enough that the massive accretion disk can be formed aroundthe black hole. Recent general relativistic numerical analyses have shown that if aprogenitor of the collapse is massive and the angular momentum is large enough,a black hole surrounded by a massive disk will be formed. However, theformation mechanism of such system has not been clarified in detail. These alsoenhance the importance of exploring the stellar core collapse to a black hole takingaccount of microphysical processes.As reviewed above, multidimensional simulations of stellar collapse in full gen-
Y. Sekiguchi eral relativity including microphysics is currently one of the most required subjectsin theoretical astrophysics. However, there has been no multidimensional code infull general relativity that self-consistently includes microphysics such as realisticEOS, electron capture, and neutrino emission. There have only existed fully generalrelativistic codes in spherical symmetry or Newtonian codes in multidimen-sion.
We have developed a fully general relativistic multidimensional codeincluding a finite-temperature nuclear EOS, self-consistent treatment of the electroncapture, and a simplified treatment of neutrino emission for the first time. In thiscode, by contrast with the previous ones, the electron capture rate is treated ina self-consistent manner and the neutrino cooling is taken into account for the firsttime. Because it is not currently feasible to fully solve the neutrino transfer equa-tions in the framework of general relativity in multidimension because of restrictionsof computational resources, it will be reasonable to take some approximation for thetransfer equations at the current status. In this paper, the so-called neutrino leakagescheme is adopted as an approximate treatment of neutrino cooling, and a generalrelativistic version of the leakage scheme is developed.1.2.
The leakage schemes
The neutrino leakage schemes as an approximate methodfor the neutrino cooling have a well-established history (e.g. Ref. 57)). The basicconcept of the original neutrino leakage schemes is to treat the following tworegions in the system separately: one is the region where the diffusion timescale ofneutrinos is longer than the dynamical timescale, and hence, neutrinos are ’trapped’(neutrino-trapped region); the other is the region where the diffusion timescale isshorter than the dynamical timescale, and hence, neutrinos stream out freely outof the system (free-streaming region). The idea of treating the diffusion region sep-arately has been applied to more advanced methods for the neutrino transfer (seee.g., Ref. 59) and references therein).Then, electron neutrinos and electron anti-neutrinos in the neutrino-trappedregion 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-trappedregion. To treat diffusive emission of neutrinos leaking out of the neutrino-trappedregion, simple phenomenological source terms based on the diffusion theory are in-troduced. In the free-streaming region, on the other hand, it is assumed thatneutrinos escape from the system without interacting with matters. Therefore, neu-trinos carry the lepton number and the energy according to the local weak-interactionrates. Note that the neutrino fractions are not solved in the original version of theleakage scheme: Only the total lepton fraction (from which the neutrino fractions arecalculated under the β -equilibrium condition) is necessary in the neutrino-trappedregion, 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 theboundary the neutrino-trapped and free-streaming regions.The boundary was given by hand as a single ’neutrino-trapping’ density ( ρ trap )without calculating the optical depths of neutrinos in the previous studies. However, the location at which the neutrino trapping occurs in fact depends strongly ull GR simulation with microphysics E ν ) as ρ trap ∝ E − ν , and hence, there are differentneutrino-trapping densities for different neutrino energies. In the previous leakageschemes, on the other hand, all neutrinos were emitted in one momentirrespective of their energy. Consequently in the case of the so-called neutrino burstemission (e.g., Ref. 60)), for example, the duration in which the neutrinos are emit-ted was shortened and the peak luminosity at the burst was overestimated. The dependence of the neutrino-trapping densities and the neutrino diffusion rates onthe neutrino energies are approximately taken into account in the recent simulationsof mergers of binary neutron star.
However, the lepton-number conservationequations for neutrinos are not solved, which is important to estimate the phasespace blocking due to neutrinos.
Transfer equations of neutrinos are not solved in the leakage schemes. There-fore, the leakage schemes cannot treat non-local interactions among the neutrinosand matters. For example, the so-called neutrino heating and the neutrino pairannihilation cannot be treated in the leakage scheme. Nevertheless, we believe adetailed general relativistic leakage scheme presented in this paper to be a valuableapproach because even by this approximated approach it is possible to incorporatethe effects of neutrinos semi-quantitatively as shown in this paper. Also, the neu-trino leakage scheme is an appropriate method for studying a number of phenomenafor which the neutrino heating and neutrino transfer are expected to be not very im-portant, e.g., prompt formation of a black hole and compact binary mergers. Bothof these phenomena are the targets of the present code.A first attempt towards a general relativistic leakage scheme was done in theprevious study.
In that study, not the region of the system but the energy mo-mentum 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-trappingdensity as in the case of the previous leakage schemes. In this paper, we develop anew code implementing the microphysical processes in the general relativistic frame-work based on the previous study.
As an application of the code, we performsimulations of stellar core collapse.A lot of improved ingredients are installed into the present code: (1) The de-pendence of the neutrino diffusion rates on the neutrino energies are approximatelytaken into account following the recent study with detailed cross sections, insteadof adopting the single neutrino-trapping density (see Appendix C). (2) The lepton-number conservation equations for neutrinos are solved to calculate self-consistentlythe chemical potentials of neutrinos. Then, the blocking effects due to the pres-ence of neutrinos and the β -equilibrium condition can be taken into account moreaccurately (see § t wp ∼ | Y e / ˙ Y e | ) ismuch shorter than the dynamical timescale t dyn in hot, dense matter regions. Note that in the previous leakage schemes, the β -equilibrium was as-sumed to be achieved in such regions (i.e. ˙ Y e = 0) and no such special treatments Y. Sekiguchi are required. See § § including effects of theso-called thermal unblocking (see Appendix A).The paper is organized as follows. First, issues in implementation of weak inter-actions and neutrino cooling in full general relativistic simulation is briefly summa-rized in §
2. Then, framework of the general relativistic leakage scheme is describedin 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. Nu-merical results obtained in this paper are shown in §
6. We devote § c = G = 1 is used otherwise stated. §
2. Issues in implementation of weak interactions and neutrino coolingin fully general relativistic simulation
Because the characteristic timescale of the weak-interaction processes (the WPtimescale t wp ∼ | Y e / ˙ Y e | ) is much shorter than the dynamical timescale t dyn in hotdense matters, the explicit numerical treatments of the weak interactions arecomputationally expensive in simple methods, as noted in the previous pioneeringwork by Bruenn: A very short timestep ( ∆t < t wp ≪ t dyn ) will be required tosolve the equations explicitly.The net rates of lepton-number and energy exchanges between matters and neu-trinos may not be large, and consequently, an effective timescale may not be as shortas the dynamical timescale. However, this does not immediately imply that one cansolve the equations explicitly without employing any prescription. For example, theachievement of β -equilibrium, where ˙ Y e = 0 is the consequence of the cancellationof two very large weak interaction processes (the electron and the electron-neutrinocaptures, see Eq. (3.20)) and of the action of the phase space blocking. Note thatthe weak interaction processes depend enormously both on the temperature and thelepton chemical potentials. Therefore, small error in the evaluation of the tempera-ture and a small deviation from the β -equilibrium due to small error in calculationof the lepton chemical potentials will result in huge error. Then, stiff source termsappear and explicit numerical evolution often becomes unstable. Indeed, we foundthat a straightforward, explicit solution of the equations did not work.In the following of this section, we describe issues of implementation of weakinteractions and neutrino cooling into the hydrodynamic equations in the conserva-tive schemes in fully general relativistic simulations. Fiest, we illustrate that in theNewtonian framework, the equations may be solved implicitly in a relatively sim-ple manner (see also Refs. 74) and 59) and referencestherein). The equations of hydrodynamics, lepton-number conservations, and neu-trino processes are schematically written as,˙ ρ = 0 , (2.1)˙ v i = S v i ( ρ, Y e , T, Q ν ) , (2.2) ull GR simulation with microphysics Y e = S Y e ( ρ, Y e , T, Q ν ) , (2.3)˙ e = S e ( ρ, Y e , T, Q ν ) , (2.4)˙ Q ν = S Q ν ( ρ, Y e , T, Q ν ) , (2.5)where ρ is the rest-mass density, v i is the velocity, Y e is the electron fraction, e isthe (internal) energy of matter, T is the temperature, and Q ν stands for the relevantneutrino quantities. We here omit the transport terms. S ’s in the right-hand sidestand for the relevant source terms. Comparing the quantities in the left-hand-sideand the argument quantities in the source terms, only the relation between e and T is nontrivial. Usually, EOSs employed in the simulation are tabularized, and one di-mensional search over the EOS table is required to solve them. Due to the relativelysimple relations between the quantities to be evolved and the argument quantities,the above equations may be solved implicitly in a straightforward (although compli-cated) manner.In the relativistic framework, the situation becomes much more complicatedin conservative schemes, because the Lorentz factor ( Γ ) is coupled with rest-massdensity and the energy density (see Eqs. (5.12) and (5.16) where w ≡ αu t is usedinstead of Γ ), and because the specific enthalpy ( h = h ( ρ, Y e , T )) is coupled with themomentum (see Eq. (5.14)).It should be addressed that the previous fully general relativistic works in thespherical symmetry are based on the so-called Misner-Sharp coordinates. There are no such complicated couplings in these
Lagrangian coordinates. Accord-ingly, the equations may be solved essentially in the same manner as in the Newtonianframework. Because no such simple Lagrangian coordinates are known in the mul-tidimensional case, the complicated couplings inevitably appear in the relativisticframework.Omitting the factors associated with the geometric variables (which are usuallyupdated before solving the hydrodynamics equations) and the transport terms, theequations to be solved in the general relativistic framework are schematically writtenas, ˙ ρ ∗ ( ρ, Γ ) = 0 , (2.6)˙ˆ u i ( u i , h ) = ˙ˆ u i ( u i , ρ, Y e , T ) = S ˆ u i ( ρ, Y e , T, Q ν , Γ ) , (2.7)˙ Y e = S Y e ( ρ, Y e , T, Q ν , Γ ) , (2.8)˙ˆ e ( ρ, Y e , T, Γ ) = S ˆ e ( ρ, Y e , T, Q ν , Γ ) , (2.9)˙ Q ν = S Q ν ( ρ, Y e , T, Q ν , Γ ) , (2.10)where ρ ∗ is a weighted density, ˆ u α is a weighted four velocity, ˆ e is a weighted energydensity (see § u α u α = −
1, which is rather complicatednonlinear equation schematically written as f normalization ( ˆ u i , Γ ) = f normalization ( u i , ρ, Y e , T, Γ ) = 0 . (2.11)The accurate calculation of the Lorentz factor and the accurate solution of the nor-malization condition are very important in the numerical relativistic hydrodynamics. Y. Sekiguchi
Now, it is obvious that the argument quantities in the source terms are notsimply related with the evolved quantities in the left-hand-side of Eqs. (2.6)–(2.11).Solving the equations implicitly is not as straightforward as in the Newtonian caseand no successful formulations have been developed. Moreover it might be not clearwhether a convergent solution can be stably obtained numerically or not, becausethey are simultaneous nonlinear equations. Therefore, it may be not a poor choiceto adopt an alternative approach in which the equations are solved explicitly withsome approximations as described in the next section ∗ ) . §
3. General relativistic neutrino leakage scheme
In this section, we describe a method for approximately solving hydrodynamicequations coupled with neutrino radiation in an explicit manner. As described in theprevious section, because of the relation t wp ≪ t dyn in the hot dense matter regions,the source terms in the equations become too stiff for the equations to be solvedexplicitly in the straightforward manner. The characteristic timescale of leakage ofneutrinos from the system t leak , by contrast, is much longer than t wp in the hot densematter region. Rather, t leak ∼ L/c ∼ t dyn where L is the characteristic length scaleof the system. On the other hand, t leak is comparable to t wp in the free-streamingregions but t wp is longer than or comparable with t dyn there. All these facts implythat the WP timescale does not directly determine the evolution of the system but theleakage timescale does. Using this fact, we approximate some of original equationsand reformulate them so that the source terms are to be characterized by the leakagetimescale t leak .3.1. Decomposition of neutrino energy-momentum tensor
The basic equations of general relativistic hydrodynamics with neutrinos are ∇ α ( T Total ) αβ = ∇ α (cid:2) ( T F ) αβ + ( T ν ) αβ (cid:3) = 0 , (3.1)where ( T Total ) αβ is the total energy-momentum tensor, and ( T F ) αβ and ( T ν ) αβ arethe energy-momentum tensor of fluids and neutrinos, respectively. Equation (3.1)can be written as ∇ α ( T F ) αβ = Q β , (3.2) ∇ α ( T ν ) αβ = − Q β , (3.3)where the source term Q α is regarded as the local production rate of neutrinosthrough the weak processes.Now, the problem is that the source term Q α becomes too stiff to solve explicitlyin hot dense matter regions where t wp ≪ t dyn . To overcome this situation, thefollowing procedures are adopted.First, it is assumed that the energy-momentum tensor of neutrinos are be decom-posed into ’trapped-neutrino’ (( T ν, T ) αβ ) and ’streaming-neutrino’ (( T ν, S ) αβ ) parts ∗ ) It should be stated that the implicit schemes are also approximated ones because a short WPtimescale associated with the weak interaction is not fully resolved. ull GR simulation with microphysics ( T ν ) αβ = ( T ν, T ) αβ + ( T ν, S ) αβ . (3.4)Here, the trapped-neutrinos phenomenologically represent neutrinos which interactsufficiently frequently with matter and are thermalized. On the other hand, thestreaming-neutrino part describes a phenomenological flow of neutrinos streamingout of the system (see also Ref. 76) where a more sophisticate method in terms ofthe 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 Q leak α : ∇ β ( T ν, S ) βα = Q leak α . (3.5)Then, the equation of the trapped-neutrino part is ∇ β ( T ν, T ) βα = Q α − Q leak α . (3.6)Third, the trapped-neutrino part is combined with the fluid part as T αβ ≡ ( T F ) αβ + ( T ν, T ) αβ , (3.7)and Eqs. (3.2) and (3.6) are combined to give ∇ β T βα = − Q leak α . (3.8)Thus, the equations to be solved are changed from Eqs. (3.2) and (3.3) to Eqs. (3.8)and (3.5). Note that the new equations only include the source term Q leak α which ischaracterized by the leakage timescale t leak . Definition of Q leak α will be given in § T αβ ) istreated as that of the perfect fluid, T αβ = ( ρ + ρε + P ) u α u β + P g αβ , (3.9)where ρ and u α are the rest mass density and the 4-velocity. The specific inter-nal energy density ( ε ) and the pressure ( P ) are the sum of contributions from thebaryons (free protons, free neutrons, α -particles, and heavy nuclei), leptons (elec-trons, positrons, and trapped-neutrinos ), and the photons as, P = P B + P e + P ν + P ph , (3.10) ε = ε B + ε e + ε ν + ε ph , (3.11)where subscripts ’ B ’, ’ e ’, ’ ph ’, and ’ ν ’ denote the components of the baryons, elec-trons 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 F α n α = P αβ n α = 0. In order to close the system, we need an explicit expres-sion of P αβ . In this paper, we adopt a simple form P αβ = χEγ αβ with χ = 1 / Y. Sekiguchi
The lepton-number conservation equations
The conservation equations of the lepton fractions are written schematically as dY e dt = − γ e , (3.13) dY νe dt = γ νe , (3.14) dY ¯ νe dt = γ ¯ νe , (3.15) dY νx dt = γ νx , (3.16)where Y e , Y νe , Y ¯ νe , and Y νx denote the electron fraction, the electron neutrino frac-tion, the electron anti-neutrino fraction, and µ and τ neutrino and anti-neutrinofractions, respectively. We note that in the previous simulations based on the leak-age schemes, the neutrino fractions were not solved.The source terms of neutrino fractions can be written, on the basis of the presentleakage scheme, as γ νe = γ local νe − γ leak νe , (3.17) γ ¯ νe = γ local¯ νe − γ leak¯ νe , (3.18) γ νx = γ local νx − γ leak νx , (3.19)where γ local ’s and γ leak ’s are the local production and the leakage rates of each neu-trino, respectively (see § γ e = γ local νe − γ local¯ νe . (3.20)Because γ local ν ’s are characterized by the WP timescale t wp , some procedures arenecessary to solve the lepton conservation equations explicitly. The following simpleprocedures are proposed to solve the equations stably.First, in each timestep n , the conservation equation of the total lepton fraction( Y l = Y e − Y νe + Y ¯ νe ), dY l dt = − γ l , (3.21)is solved together with the conservation equation of Y νx , Eq. (3.16), in advance ofsolving whole of the lepton conservation equations (Eqs. (3.13) – (3.16)). Note thatthe source term γ l = γ leak νe − γ leak¯ νe is characterized by the leakage timescale t leak sothat 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 ( Y βe , Y βνe , and Y β ¯ νe ) are calculated from the evolved Y l . ull GR simulation with microphysics Y βνe and Y β ¯ νe as the maximum allowed values of the neutrinofractions in the next timestep n + 1, the source terms are limited so that Y ν ’s inthe timestep n + 1 cannot exceed Y βν ’s. Then, the whole of the lepton conservationequations (Eqs. (3.13) – (3.16)) 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 µ p , µ n , µ e , µ νe and µ ¯ νe are the chemical potentials of protons, neutrons,electrons, electron neutrinos, and electron anti-neutrinos, respectively. If both con-ditions are satisfied, the values of the lepton fractions in the timestep n + 1 are set tobe those in the β -equilibrium value; Y βe , Y βνe , and Y β ¯ νe . On the other hand, if eitheror both conditions are not satisfied, the lepton fractions in the timestep n + 1 is setto be those obtained by solving whole of the lepton-number conservation equations.A limiter for the evolution of Y νx may be also necessary in the case where the pairprocesses are dominant, for example, in the simulations for collapse of population IIIstellar core. In this case, the value of Y νx at the pair equilibrium (i.e. at µ νx = 0), Y pair νx may be used to limit the source term.3.3. Definition of leakage rates
In this subsection the definitions of the leakage rates Q leak α and γ leak ν are pre-sented. Because Q leak ν may be regarded as the emissivity of neutrinos measured inthe fluid rest frame , Q leak α is defined as Q leak α = Q leak ν u α . (3.24)The leakage rates Q leak ν and γ leak ν are assumed to satisfy the following properties.1. The leakage rates approach the local rates Q local ν and γ local ν in the low density,transparent region.2. The leakage rates approach the diffusion rates Q diff ν and γ diff ν 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 andthe diffusion rates can be determined based on the diffusion theory (see appendicesfor the definition of the local and diffusion rates adopted in this paper). There will beseveral prescriptions to satisfy the requirement (iii). In this paper, the leakagerates are defined by Q leak ν = (1 − e − bτ ν ) Q diff ν + e − bτ ν Q local ν , (3.25) γ leak ν = (1 − e − bτ ν ) γ diff ν + e − bτ ν γ local ν , (3.26)where τ ν is the optical depth of neutrinos and b is a parameter which is typicallyset as b − = 2 /
3. The optical depth can be computed from the cross sections in astandard manner. Y. Sekiguchi
In the present implementation, it is not necessary to artificially divide the systeminto neutrino-trapped and free-streaming regions by the single neutrino-trappingdensity. Therefore there is no discontinuous boundary which existed in the previousleakage schemes.
As the local production reactions of neutrinos, the electron and positron cap-tures ( γ ec νe and γ pc¯ νe ), the electron-positron pair annihilation ( γ pair ν e ¯ ν e for electron-type neutrinos and γ pair ν x ¯ ν x for the other type), the plasmon decays ( γ plas ν e ¯ ν e and γ plas ν x ¯ ν x ),and the Bremsstrahlung processes ( γ Brems ν e ¯ ν e and γ Brems ν x ¯ ν x ) are considered in this pa-per. Then, the local rates for the neutrino fractions are γ local νe = γ ec νe + γ pair ν e ¯ ν e + γ plas ν e ¯ ν e + γ Brems ν e ¯ ν e , (3.27) γ local¯ νe = γ pc¯ νe + γ pair ν e ¯ ν e + γ plas ν e ¯ ν e + γ Brems ν e ¯ ν e , (3.28) γ local νx = γ pair ν x ¯ ν x + γ plas ν x ¯ ν x + γ Brems ν x ¯ ν x . (3.29)Similarly, the local neutrino energy emission rate Q local ν is given by Q local ν = Q ec νe + Q pc¯ νe + 2 ( Q pair ν e ¯ ν e + Q plas ν e ¯ ν e + Q Brems ν e ¯ ν e )+ 4 ( Q pair ν x ¯ ν x + Q plas ν x ¯ ν x + Q Brems ν x ¯ ν x ) . (3.30)The explicit forms of the local rates in Eqs. (3.27)–(3.30) will be found in AppendicesA and B.We follow the recent work by Rosswog and Liebend¨orfer for the diffusiveneutrino emission rates γ diff ν and Q diff ν in Eqs (3.25) and (3.26). The explicit formsof γ diff ν and Q diff ν 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., whichis derived by the relativistic mean field theory based on the relativistic Br¨uckner-Hartree-Fock theory.
The so-called parameter set TM1 is adopted to reproducecharacteristic properties of heavy nuclei. The maximum mass of a cold spherical neu-tron star in this EOS is much larger than the canonical neutron star mass ≈ . M ⊙ as ≈ . M ⊙ . The framework of the relativistic mean field theory is extendedwith the Thomas-Fermi spherical cell model approximation to describe not only thehomogeneous matter but also an inhomogeneous one.The thermodynamical quantities of dense matter at various sets of ( ρ, Y p , T ) arecalculated to construct the numerical data table for simulation. The table covers awide range of density 10 . –10 . g/cm , electron fraction 0 . .
56, and temperature0–100 MeV, which are required for supernova simulation. It should be noted that thecausality is guaranteed to be satisfied in this framework, whereas the sound velocitysometimes 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 relativisticEOS. ull GR simulation with microphysics
Because the table of the original EOS by Shen et al. does not include thethermodynamical quantities of the leptons (electrons, positrons, and neutrinos ifnecessary) 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 andpositron, the charge neutrality condition Y p = Y e should be solved to determine theelectron chemical potential µ e for each value of the baryon rest-mass density ρ andthe temperature T in the EOS table. Namely, it is required to solve the equation n e ( µ e , T ) ≡ n − − n + = ρY e m u (4.1)in terms of µ e for given values of ρ , T , and Y e (= Y p ). Here, m u = 931 . n − and n + are the total number densities (i.e., includingthe electron-positron pair) of the electrons and positrons, respectively.Assuming that the electrons obey the Fermi-Dirac distribution (which is derivedunder the assumption of thermodynamic equilibrium), the number density ( n − ), thepressure ( P − ), and the internal energy density ( u − ) of the electron are written as n − = 1 π ~ Z ∞ p dp exp [ − η e + ˜ ǫ/k B T ] + 1 , (4.2) P − = 1 π ~ Z ∞ p ( ∂ ˜ ǫ/∂p ) dp exp [ − η e + ˜ ǫ/k B T ] + 1 , (4.3) u − = 1 π ~ Z ∞ p ˜ ǫdp exp [ − η e + ˜ ǫ/k B T ] + 1 . (4.4)Here ~ , k B , and η e ≡ µ e /k B T are the Planck’s constant, the Boltzmann’s constantand the so-called degeneracy parameter. ˜ ǫ ( p ) = p m e c + p − m e c is the kineticenergy of a free electron. If we choose the zero point of our energy scale for electronsat ˜ ǫ = 0, we have to assign a total energy of ˜ ǫ + = p m e c + p + m e c to a freepositron. Thus the number density ( n + ), the pressure ( P + ), and the internalenergy density ( u + ) of positrons are given by n + = 1 π ~ Z ∞ p dp exp [ − η + + ˜ ǫ + /k B T ] + 1 , (4.5) P + = 1 π ~ Z ∞ p ( ∂ ˜ ǫ + /∂p ) dp exp [ − η + + ˜ ǫ + /k B T ] + 1 , (4.6) u + = 1 π ~ Z ∞ p (˜ ǫ + 2 m e c ) dp exp [ − η + + ˜ ǫ + /k B T ] + 1 , (4.7)where η + = − η e is the degeneracy parameter of the positrons.4 Y. Sekiguchi
Photons
The pressure and the specific internal energy density of photons are given by P r = a r T , ε r = a r T ρ , (4.8)where a r is the radiation constant a r = ( π k B ) / (15 c ~ ) and c is the velocity of light.4.4. Trapped neutrinos
In this paper, the trapped-neutrinos are assumed to interact sufficiently fre-quently with matter that be thermalized. Therefore they are described as idealFermi gases with the matter temperature. Then, from the neutrino fractions Y ν , thechemical potentials of neutrinos are calculated by solving Y ν = Y ν ( µ ν , T ) . (4.9)Using the chemical potentials, µ ν , and the matter temperature, the pressure and theinternal energy of the trapped-neutrinos are calculated in the same manner as forelectrons.4.5. The sound velocity
In the high-resolution shock-capturing scheme for hydrodynamics, we in generalneed to evaluate the sound velocity c s , c s = 1 h " ∂P∂ρ (cid:12)(cid:12)(cid:12)(cid:12) ǫ + Pρ ∂P∂ǫ (cid:12)(cid:12)(cid:12)(cid:12) ρ . (4.10)The derivatives of the pressure are calculated by ∂P∂ρ (cid:12)(cid:12)(cid:12)(cid:12) ǫ = X i = B,e,r,ν ∂P i ∂ρ (cid:12)(cid:12)(cid:12)(cid:12) T − ∂P i ∂T (cid:12)(cid:12)(cid:12)(cid:12) ρ X j = B,e,r,ν ∂ǫ j ∂ρ (cid:12)(cid:12)(cid:12)(cid:12) T X k = B,e,r,ν ∂ǫ k ∂T (cid:12)(cid:12)(cid:12)(cid:12) ρ − , (4.11) ∂P∂ǫ (cid:12)(cid:12)(cid:12)(cid:12) ρ = X i = B,e,r,ν ∂P i ∂T (cid:12)(cid:12)(cid:12)(cid:12) ρ X j = B,e,r,ν ∂ǫ j ∂T (cid:12)(cid:12)(cid:12)(cid:12) ρ − , (4.12)where ’ B ’, ’ e ’, ’ ph ’ and ’ ν ’ in the sum denote the baryon, electron, photons, andneutrino quantities.The derivatives for the baryon parts are evaluated by taking a finite differenceof the table data. On the other hand, the derivatives for the electron parts can beevaluated semi-analytically. Because there are in general the phase transition regionsin an EOS table for baryons and moreover the EOS may contain some non-smoothspiky structures, careful treatments are necessary when evaluating the derivatives ofthermodynamical quantities. In the present EOS table, the derivatives are carefullyevaluated so that there are no spiky behaviors in the resulting sound velocities. ull GR simulation with microphysics §
5. Basic equations and Numerical methods
Einstein’s equation and gauge conditions
The standard variables in the 3+1 decomposition are the three-dimensional met-ric γ ij and the extrinsic curvature K ij on the three-dimensional hypersurface de-fined by γ µν ≡ g µν + n µ n ν , (5.1) K µν ≡ − L – n γ µν , (5.2)where g µν is the spacetime metric, n µ is the unit normal to a three-dimensionalhypersurface, and L – n is the Lie derivative with respect to the unit normal n µ . Thenwe can write the line element in the form ds = − α dt + γ ij ( dx i + β i dt )( dx j + β j dt ) , (5.3)where α and β i are the lapse function and the shift vector which describe the gaugedegree of freedom.In the BSSN reformulation, the spatial metric γ ij is conformally decom-posed as γ ij = e φ ˜ γ ij where the condition det(˜ γ ij ) = 1 is imposed for the conformalmetric ˜ γ ij . From this condition, the conformal factor is written as φ = ln γ and γ ≡ det( γ ij ). The extrinsic curvature K ij is decomposed into the trace part K and the traceless part A ij as K ij = A ij + (1 / γ ij K . The traceless part A ij isconformally decomposed as A ij = e φ ˜ A ij . Thus the fundamental quantities for theevolution equation are now split into φ, ˜ γ ij , K , and ˜ A ij . Furthermore, the auxiliaryvariable F i ≡ δ jk ∂ k ˜ γ ij is introduced in the BSSN reformulation. The basic equations to be solved are (cid:16) ∂ t − β k ∂ k (cid:17) φ = 16 (cid:16) − αK + ∂ k β k (cid:17) , (5.4) (cid:16) ∂ t − β k ∂ k (cid:17) ˜ γ ij = − α ˜ A ij + ˜ γ ik ∂ j β k + ˜ γ jk ∂ i β k −
23 ˜ γ ij ∂ k β k , (5.5) (cid:16) ∂ t − β k ∂ k (cid:17) K = − D k D k α + α (cid:20) ˜ A ij ˜ A ij + 13 K (cid:21) + 4 πα ( ρ h + S ) , (5.6) (cid:16) ∂ t − β k ∂ k (cid:17) ˜ A ij = e − φ (cid:20) α (cid:18) R ij − e φ ˜ γ ij R (cid:19) − (cid:18) D i D j α − e φ ˜ γ ij D k D k α (cid:19)(cid:21) + α (cid:16) K ˜ A ij − A ik ˜ A kj (cid:17) + ˜ A ik ∂ j β k + ˜ A jk ∂ i β k −
23 ˜ A ij ∂ k β k − πα (cid:18) e − φ S ij −
13 ˜ γ ij S (cid:19) , (5.7) (cid:16) ∂ t − β k ∂ k (cid:17) F i = − παj i +2 α (cid:26) f kj ∂ j ˜ A ik + ˜ A ik ∂ j f kj −
12 ˜ A jl ∂ i h jl + 6 ˜ A ki ∂ k φ − ∂ i K (cid:27) + δ jk n − A ij ∂ k α + (cid:16) ∂ k β l (cid:17) ∂ l h ij Y. Sekiguchi + ∂ k (cid:18) ˜ γ il ∂ j β l + ˜ γ jl ∂ i β l −
23 ˜ γ ij ∂ l β l (cid:19)(cid:27) , (5.8)where (3) R , (3) R ij , and D i are the Ricci scalar, the Ricci tensor, and the covariantderivative associated with three-dimensional metric γ ij , respectively. The mattersource terms, ρ h ≡ ( T Total ) αβ n α n β , j i ≡ − ( T Total ) αβ γ iα n β , and S ij ≡ ( T Total ) αβ γ iα γ jβ ,are the projections of the stress-energy tensor with respect to n µ and γ µν , and S ≡ γ ij S ij .We assume the axial symmetry of the spacetime and the so-called Cartoonmethod is adopted to avoid problems around the coordinate singularities ofthe cylindrical coordinates. Except for this, the numerical schemes for solving theEinstein’s equation are essentially the same as those in Ref. 90). We use 4th-order fi-nite difference scheme in the spatial direction and the 3rd-order Runge-Kutta schemein the time integration. The advection terms such as β i ∂ i φ are evaluated by a 4th-order upwind scheme.As the gauge conditions for the lapse, we use the so-called 1 + log slicing: ( ∂ t − L – β ) α = − Kα. (5.9)It is known that the 1 + log slicing enables to perform a long term evolution ofneutron stars as well as has strong singularity avoidance properties in the black holespacetime.The shift vector is determined by solving the following dynamical equation ∂ t β k = ˜ γ kl ( F l + ∆t∂ t F l ) . (5.10)Here the second term in the right-hand side is necessary for numerical stability. The hydrodynamic equations in leakage scheme
The basic equations for general relativistic hydrodynamics in our leakage schemeare the continuity equation, the lepton-number conservation equations, and the localconservation equation of the energy-momentum. We assume the axial symmetry ofthe spacetime and the hydrodynamics equations are solved in the cylindrical coordi-nates ( ̟, ϕ, z ) where ̟ = p x + y . In the axisymmetric case, the hydrodynamicsequations should be written in the cylindrical coordinate. On the other hand, in theCartoon method, Einstein’s equation are solved in the y = 0 plane for which x = ̟ , u ̟ = u x , u ϕ = xu y , and the other similar relations hold for vector and ten-sor quantities. Taking into these facts, the hydrodynamic equations may be writtenusing the Cartesian coordinates replacing ( ̟, ϕ ) by ( x, y ). In the following, we writedown explicit forms of the equations for the purpose of convenience. Numerical testsfor basic parts of the code of solving the hydrodynamics equations are extensivelyperformed in Ref. 89). The equations are solved using the third-order high-resolutioncentral scheme of Kurganov and Tadmor. ∇ α ( ρu α ) = 0 . (5.11) ull GR simulation with microphysics ρ ∗ ≡ ρwe φ and v i ≡ u i /u t where w ≡ αu t . Then, the continuity equationis written as ∂ t ( ρ ∗ ) + 1 x ∂ x ( ρ ∗ v x ) + ∂ z ( ρ ∗ v z ) = 0 . (5.12)Using the continuity equation, the lepton-number conservation equations (3.13)– (3.16) are written as ∂ t ( ρ ∗ Y L ) + 1 x ∂ x ( ρ ∗ Y L v x ) + ∂ z ( ρ ∗ Y L v z ) = ρ ∗ γ L , (5.13)where Y L and γ L are abbreviated expressions of the lepton fractions and the sourceterms.5.2.2. Energy-momentum conservationAs fundamental variables for numerical simulations, we define the quantitiesˆ u i ≡ hu i and ˆ e ≡ hw − P ( ρw ) − . Then, the Euler equation γ αi ∇ β T βα = − γ αi Q leak α ,and the energy equation n α ∇ β T αβ = − n α Q leak α can be written as ∂ t ( ρ ∗ ˆ u A ) + 1 x ∂ x h x n ρ ∗ ˆ u A v x + P αe φ δ xA oi + ∂ z h ρ ∗ ˆ u A v z + P αe φ δ zA i = − ρ ∗ (cid:20) wh∂ A α − ˆ u i ∂ A β i + αe − φ wh ˆ u k ˆ u l ∂ A ˜ γ kl − αh ( w − w ∂ A φ (cid:21) + P ∂ A ( αe φ ) + ( ρ ∗ u y v y + P αe φ ) δ xA x − αe φ Q leak A , (5.14) ∂ t ( ρ ∗ ˆ u y ) + 1 x ∂ x (cid:0) x ρ ∗ ˆ u y v y (cid:1) + ∂ z ( ρ ∗ ˆ u y v z ) = − αe φ Q leak y , (5.15) ∂ t ( ρ ∗ ˆ e ) + 1 x ∂ x h x n ρ ∗ v x ˆ e + P e φ ( v x + β x ) oi + ∂ z h ρ ∗ v z ˆ e + P e φ ( v z + β z ) i = αe φ P K + ρ ∗ u t h ˆ u k ˆ u l K kl − ρ ∗ ˆ u i γ ij D j α − αe φ Q leak α n α , (5.16)where the subscript A denotes x or z component.The evolution equation of streaming-neutrinos ∇ β ( T ν, S ) βα = Q leak α gives ∂ t ( ˆ E ) + 1 x ∂ x h x ( α ˆ F x − β x ˆ E ) i + ∂ z h ( α ˆ F z − β z ˆ E ) i = α ˆ EK − ˆ F k ∂ k α + αe φ Q leak a n a , (5.17) ∂ t ( ˆ F A ) + 1 x ∂ x (cid:20) x (cid:18) α ˆ Eδ xA − β x ˆ F A (cid:19)(cid:21) + ∂ z (cid:20)(cid:18) α ˆ Eδ zA − β z ˆ F A (cid:19)(cid:21) = − ˆ E∂ A α + ˆ F k ∂ A β k + 2 α ˆ E∂ A φ + ( ˆ E/ − ˆ F y β y ) δ xA x + αe φ Q leak A , (5.18) ∂ t ( ˆ F y ) − x ∂ x h x β x ˆ F y i − ∂ z h β z ˆ F y i = αe φ Q leak y , (5.19)where ˆ E = e φ E and ˆ F i = e φ F i , and the subscript A again denotes x or z compo-nent. The closure relation P αβ = Eγ αβ / Y. Sekiguchi
Recover of ( ρ , Y e / Y l , T ) The quantities numerically evolved in the relativistic hydrodynamics are theconserved quantities ( ρ ∗ , ˆ u i , ˆ e ) and the lepton fraction Y e or Y l . The argumentvariables, ( ρ , ( Y e or Y l ), T ), of the EOS table, together with the weight factor w = p γ ij u i u j , should be calculated from the conserved quantities at each timeslice. Note that the electron ( Y e ) or lepton fraction ( Y l ) is readily given by numericalevolution at each time slice whereas ρ , u i , and T are not. This fact requires us tofind an efficient method for determining w .5.3.1. Non- β -equilibrium caseIn the case that the β -equilibrium condition is not satisfied, the argument quan-tities ( ρ , Y e , T ) can be reconstructed from the conserved quantities in the followingstraightforward manner.1. Give a trial value of w , referred to as ˜ w . Then, one obtains a trial value of therest mass density from ˜ ρ = ρ ∗ / ( ˜ we φ ).2. A trial value of the temperature, ˜ T , can be obtained by solving the followingequation: ˆ e = ε + ˜ P ˜ ρ ! ˜ w − ˜ P ˜ ρ ˜ w ≡ ˜ e (˜ ρ, Y e , ˜ T ) . (5.20)Here, one dimensional search over the EOS table is required to obtain ˜ T .3. The next trial value of w is given by ˜ w = q e − φ ˜ γ ij ˆ u i ˆ u j ˜ h − , where thespecific enthalpy was calculated as ˜ h = ˜ h (˜ ρ, Y e , ˜ T ) in the step 2.4. Repeat the procedures (1)–(3) until a required degree of convergence is achieved.Convergent solutions of the temperature and w are obtained typically in 10iterations.5.3.2. The β -equilibrium caseOn the other hand, in the case that the β -equilibrium condition is satisfied, onehas to reconstruct the argument quantities ( ρ, Y e , T ) from the conserved quantitiesand Y l , under the assumption of the β -equilibrium. In this case, two-dimensionalrecover ( Y l , ˆ e ) = ⇒ ( Y e , T ) would be required for a given value of ˜ w . A seriousproblem is that in this case, there may be more than one combination of ( Y e , T )which gives the same values of Y l and ˆ e . Therefore, we have to adopt a differentmethod to recover ( ρ, Y e , T ). Under the assumption of the β -equilibrium, the electronfraction is related to the total lepton fraction: Y e = Y e ( ρ, Y l , T ). Using this relation,the EOS table can be rewritten in terms of the argument variables of ( ρ , Y l , T ).Then, the same strategy as in the non- β -equilibrium case can be adopted. Namely,1. Give a trial value ˜ w . Then one obtains a trial value of the rest mass density.2. A trial value of the temperature can be obtained by solving ˆ e = ˜ e (˜ ρ, Y l , ˜ T ), withone dimensional search over the EOS table.3. The next trial value of w is given by ˜ w = q e − φ ˜ γ ij ˆ u i ˆ u j ˜ h − .4. Repeat the procedures (1)–(3) until a required degree of convergence is achieved.The electron fraction is given as Y e = Y e ( ρ, Y l , T ) in the (new) EOS table. ull GR simulation with microphysics Model Φ c ≤ . ≤ Φ c ≤ . ≤ Φ c ≤ . ≤ Φ c ≤ . Φ c ≥ . ∆x δ N
444 668 924 1212 1532 L (km) 2330 2239 2188 2124 2103S15 ∆x δ N
316 444 636 828 1020 L (km) 2244 2151 2073 2043 2000Table I. Summary of the regridding procedure. The values of the minimum grid spacing ∆x (inunits of km), the non-uniform-grid factor δ , and the grid number N for each range of Φ c = 1 − α c are listed. In the case of a simplified or analytic EOS, the Newton-Raphson method maybe applied to recover the primitive variables. In the case of a tabulated EOS, by con-trast, the Newton-Raphson method may not be a good approach because it requiresderivatives of thermodynamical quantities which in general cannot be calculatedprecisely from a tabulated EOS by the finite differentiating method.5.4.
Grid Setting
In numerical simulations, we adopt a nonuniform grid, in which the grid spacingis increased as dx j +1 = (1 + δ ) dx j , dz l +1 = (1 + δ ) dz l (5.21)where dx j ≡ x j +1 − x j , dz l ≡ z l +1 − z l , and δ is a constant. In addition, a regriddingtechnique is adopted to assign a sufficiently large number of grid points insidethe collapsing core, saving the CPU time efficiently. The regridding is carried outwhenever the characteristic radius of the collapsing star decreases by a factor of a2–3. At each regridding, the minimum grid spacing is decreased by a factor of ∼ δ is unchanged (see Table I).All the quantities on the new grid are calculated using the fifth-order Lagrangeinterpolation. To avoid discarding the matter in the outer region, we also increasethe grid number at each regridding. For the regridding, we define a relativisticgravitational potential Φ c ≡ − α c ( Φ c >
0) where α c is the central value of thelapse function. Because Φ c is approximately proportional to M/R where M and R are characteristic mass and radius of the core, Φ − c can be used as a measure ofthe characteristic length scale of the stellar core for the regridding. In Table I, wesummarize the regridding parameters of each level of the grid. §
6. Results
As a test problem, we perform a collapse simulation of spherical presupernovacore and compare the results with those in the previous studies, to see the validity ofthe present code. Most of the following results are not novel astrophysically, but arenovel in the sense that stellar core collapse can be followed by a multidimensionalfully general relativistic simulation taking account of microphysical processes . In0
Y. Sekiguchi α Time after bounce t b [ms] D e n s it y ρ [ g / c m ] High resolutionLow resolution
Fig. 1. Evolution of the central density ρ c (upper panel) and the central value of the lapse function α c (lower panel). The solid curves are results for the finer grid resolution and the dotted curvesare results of the coarser grid resolution. § § § § Initial condition
In this paper, we adopt a recent presupernova model of massive star by Woosley,Heger, and Weaver: M ⊙ model with solar metallicity (hereafter S15 model).We follow the dynamical evolution of the central part which constitutes the Fe coreand the inner part of the Si-shell. We read in the density, the electron fraction,the temperature and the velocity ( v i ) of the original initial data and derive otherthermodynamical quantities using the EOS table.Note that the procedure of remapping the original initial data into the gridadopted in the numerical simulations is coordinate-dependent in general relativity.In this paper, we read in the original data as a function of the coordinate radius . Inthis case, the baryon rest-mass of the core is slightly larger than the original one,because it is defined by M ∗ = Z ρ ∗ dx = Z ρ ( we φ ) d x, (6.1)where we φ > ull GR simulation with microphysics Core bounce and shock formation
Figure 1 displays the time evolution of the central rest-mass density ρ and thecentral value of the lapse function. This figure shows that the stellar core collapse toa neutron star can be divided into three phases; the infall phase, the bounce phase,and the quasi-static evolution phase (see Refs. 21) and 24) for the case of rotationalcollapse). The general feature of the collapse is as follows.The infall phase sets in due to gravitational instability of the iron core trig-gered by the sudden softening of the EOS, which is associated primarily with theelectron capture and partially with the photo-dissociation of the heavy nuclei. Thecollapse in an early phase proceeds almost homologously. However, the collapse inthe central region is accelerated with time because the electron capture reduces thedegenerate pressure of electrons which provides the main part of the total pressure.Furthermore, the neutrino emission associated with the electron capture reduces thethermal pressure of the core. Here the inner part of the core, which collapses nearlyhomologously with a subsonic infall velocity, constitutes the inner core. On the otherhand, the outer region in which the infall velocity is supersonic constitutes the outercore.The collapse proceeds until the central part of the iron core reaches the nucleardensity ( ∼ × g/cm ), and then, the inner core experiences the bounce. Becauseof its large inertia and large kinetic energy induced by the infall, the inner coreovershoots its hypothetical equilibrium state. The stored internal energy of theinner core at the maximum compression is released through strong pressure wavesgenerated inside the inner core. The pressure waves propagate from the center to theouter region until they reach the sonic point located at the edge of the inner core.Because the sound cones tilt inward beyond the sonic point, the pressure disturbancecannot propagate further and forms a shock just inside the sonic point. A shock waveis formed at the edge of the inner core and propagates outward.After this phase, the proto-neutron star experiences the quasi-static evolutionphase. In this phase, the central value of density (the lapse function) increases(decreases) gradually, because the matter in the outer region falls into the proto-neutron star and because neutrinos are emitted carrying away the energy and lepton-number from the proto-neutron star.Figure 2 shows the radial profiles in the equator of the lepton fractions at thebounce. The central values of the electron, the electron-neutrinos, and the total-lepton fractions are ≈ .
32, 0.05, and 0.37, respectively. The electron-anti-neutrinofraction is almost zero through out the core because only very small amount ofpositrons exist due to the high degree of electron degeneracy.6.3.
Neutrino bursts and stall of shock
As the shock wave propagates outward, the kinetic energy of the infall matteris converted into the thermal energy behind the shock. The conversion rate of infallkinetic energy may be estimated approximately as L heat ∼ πR s ( ρ infall v / Y. Sekiguchi L e p t on F r ac ti on s R [km] electron fraction ν e fraction anti ν e fraction total lepton fraction Fig. 2. The radial profiles of the electron, ν e -, ¯ ν e -, and the total lepton fraction at the bounce. Theresults for the finer grid resolution (solid curve) and for the coarser grid resolution (the dottedcurves) are shown together. The two results are almost identical. ∼ . × ergs / s (cid:18) R s
100 km (cid:19) (cid:18) ρ infall g / cm (cid:19) (cid:16) v infall . c (cid:17) , (6.2)where R s and ρ infall are radius of the shock wave and the density of infall matter,and we here recover the velocity of the light ( c ). Here, we assume that all the kineticenergy is converted to the thermal energy.At the same time, the shock wave suffers from the energy loss by the photo-dissociation of the iron to α -particles and free nucleons. The fraction of this energyloss is ǫ diss ∼ . × ergs per 0 . M ⊙ . (6.3)Thus, the energy loss rate due to the photo-dissociation is L diss ∼ ˙ M shock ǫ diss ∼ . × ergs / s (cid:18) R s
100 km (cid:19) (cid:18) ρ infall g / cm (cid:19) (cid:16) v infall . c (cid:17) , (6.4)where ˙ M shock ∼ πR s ρ infall v infall is mass-accretion rate to the shock front.The ratio of L heat to L diss is L heat L diss ≈ . (cid:16) v infall . c (cid:17) . (6.5)Therefore the energy loss rate by the photo-dissociation will eventually overcome thehydrodynamic power, because the infall velocity, which is ≈ ( GM ic /R s ) / , decreasesas the shock wave propagates outward. ull GR simulation with microphysics -20 -10 0 10 20 30 40 50 60 L u m i no s it y [ e r g / s ] Time [ms] e neutrinoe antineutrino µ , τ neutrinos Fig. 3. Time evolution of the neutrino luminosities. The results in the finer grid resolution (solidcurves) and in the coarser grid resolution (dashed curves) are shown together. The two resultsare approximately identical until the convective phase sets in, whereas small disagreement isfound in the convective phase.
Furthermore, when the shock wave crosses the neutrino-sphere, spiky burst emis-sions of neutrinos, the so-called neutrino bursts, occur: Neutrinos in the hot post-shock region are copiously emitted without interacting core matter. Figure 3 showsthe neutrino luminosity as a function of time calculated by L ν = Z αe φ u t ˙ Q ν d x. (6.6)The peak luminosity is L ν e ≈ . × ergs/s. This neutrino burst significantlyreduces the thermal energy of the shock. Consequently, the shock wave stalls at ≈
80 km soon after the neutrino burst. The peak luminosity and the shock-stall ra-dius agree approximately with the previous one-dimensional fully general relativisticstudy.
When the shock wave stalls, negative gradients of the entropy per baryon and thetotal-lepton (electron) fraction appear because neutrinos carry away both the energyand the lepton number. Figure 4 shows the radial profiles of the infall velocity, thedensity, the entropy per baryon, and the total lepton fraction in the equator. Thisfigure clearly shows that negative gradients of the entropy per baryon and the totallepton fraction are formed above the neutrino sphere. As we shall see in § Y. Sekiguchi E n t r opy p e r B a r yon [ k B ] R [km] bounce2ms after bounce6ms after bounce
5 10 20 50 100 200 500 1000 0 0.1 0.2 0.3 0.4 0.5 T o t a l L e p t on F r ac ti on R [km] -8-6-4-2 0 2 V e l o c it y v R [ k m / s ] D e n s it y [ l og ( ρ g / c m )] Fig. 4. The radial profiles of the infall velocity, the density, the entropy per baryon, and the totallepton fraction at bounce, at 2 ms and 6ms after bounce. The results for the finer grid resolution(solid curves) and for the coarser grid resolution (the dotted curves) are shown together andthey are shown to be approximately identical.
Convective activities
Let us investigate the stability of the envelope of the proto-neutron star followingLattimer and Mazurek.
We consider the following parameter N ≡ g eff ρ (cid:20)(cid:18) dρdr (cid:19) amb − (cid:18) dρdr (cid:19) blob (cid:21) , (6.7)where g eff is the effective gravitational acceleration defined to be positive in thenegative radial direction, the subscript ’amb’ refers to the ambient core structure,and ’blob’ denotes the blob element which is under an isolated displacement. Thecondition N < (cid:18) dρdr (cid:19) blob = (cid:18) dρdP (cid:19) blob (cid:18) dPdr (cid:19) amb . (6.8) ull GR simulation with microphysics Fig. 5. Snapshots of the contours of the density (top left panels), the electron fraction Y e (top rightpanels), the entropy per baryon (bottom left panels), and the local neutrino energy emissionrate (bottom right panels) in the x - z plane at selected time slices. Using this relation, Eq. (6.7) is written as N = g eff ρ (cid:18) dρdP (cid:19) blob (cid:20)(cid:18) dPdρ (cid:19) blob (cid:18) dρdr (cid:19) amb − (cid:18) dPdr (cid:19) amb (cid:21) . (6.9)Because the pressure is a function of the entropy per baryon, the density, and thelepton fraction, ( dP/dr ) amb is rewritten to give N = g eff ρ (cid:18) ∂ρ∂P (cid:19) s,Y l "(cid:18) ∂P∂s (cid:19) ρ,Y l (cid:18) dsdr (cid:19) amb + (cid:18) ∂P∂Y l (cid:19) ρ,s (cid:18) dY l dr (cid:19) amb . (6.10)Here, we also assume that the blob elements do not interact the ambient mattersboth thermally and chemically, i.e. ds = dY l = 0 for the blob. Then, we have (cid:18) dPdρ (cid:19) blob = (cid:18) ∂P∂ρ (cid:19) s,Y l . (6.11)Equation (6.10) shows that when the pressure derivatives of given EOS (( ∂P/∂s ) ρY e and ( ∂P/∂Y l ) ρs ) are positive, configurations with negative gradients of entropy and Y l ( N <
0) are unstable. (Note that in the above treatment, we have ignored thedissociative effects caused by energy and lepton transports due to neutrinos.) Thus,6
Y. Sekiguchi
Fig. 6. Snapshots of the contours of gradients associated with the entropy per baryon( ∂P/∂s ) Y l ,ρ ( ds/dr ) (right panels) and associated with the lepton fraction ( ∂P/∂Y l ) s,ρ ( dY l /dr )(left panels) in the x - z plane at selected time slices. the negative gradients of the entropy per baryon and the total lepton fraction formedabove the neutrino sphere lead to the convective overturn (the proto-neutron starconvection). Indeed, convection occurs in our simulation.Figure 5 shows contours of the density, the electron fraction, the entropy perbaryon, and the neutrino energy-emission rate. Convective motions are activated atabout 8 ms after the bounce in the region located above the neutrino-sphere wherethe gradients of the entropy per baryon and Y l are imprinted (see Fig. 4). At about10 ms after the bounce, the lepton rich, hot blobs rise to form ’fingers’ (see in topleft panel in Fig. 5). Note that the neutrino energy emission rate in this finger isrelatively higher than that in other region. This is responsible for the small humpseen in the time-evolution of neutrino luminosity (see Fig. 3). Subsequently, the hotfingers expand to form ’mushroom structures’, and push the surface of the stalledshock (see top right panel in Fig. 5). At the same time, the lepton poor, coldermatters sink down to the proto-neutron star ( r .
20 km). The entropy per baryonjust behind the shock increases to be s & k B and the stalled shock gradually movesoutward to reach r ≈
200 km. As the hot, lepton rich matters are dug out from theregion below the neutrino-sphere, the neutrino luminosity is enhanced (see Fig. 3). ull GR simulation with microphysics
A more detailed comparison with the previous simulations is given in § ∂P/∂s ) Y l ,ρ ( ds/dr ) (right panelsin Fig. 6) and associated with the lepton fraction ( ∂P/∂Y l ) s,ρ ( dY l /dr ) (left panels inFig. 6). This figure clearly shows that negative gradient of the entropy per baryonis more important for the convection activated promptly after the bounce.6.5. Comparison with the previous studies
To check the validity of the code, the results presented in § § § in which 1D general relativisticBoltzmann equation is solved for neutrino transfer with relevant weak interactionprocesses. Because neutrino heating processes ( ν e + n → p + e − and ¯ ν e + p → n + e + ) are not included in the present implementation, and on the other hand,multidimensional effects such as convection cannot be followed in the one-dimensionalreference simulations, we pay particular attention to comparing results during thecollapse and the early phase ( ∼
10 ms) after the bounce (see results in § § where simple leakage schemes based on the singleneutrino-trapping density were adopted. Quantitatively, the negative gradients ofthe entropy per baryon and the lepton fraction are little bit steeper in the presentsimulation than those in 1D full Boltzmann simulations. The reason may be partlybecause the transfers of lepton-number and energy are not fully solved in the presentleakage scheme. Except for this small quantitative difference, the two results agreewell.For validating a scheme for the neutrino cooling, agreement of the neutrino lumi-nosities with those by 1D full Boltzmann simulation should be particularly checkedbecause they depend on both implementations of weak interactions (especially elec-tron capture in the present case) and treatments of neutrino cooling (the detailedleakage scheme). Also, accurate computation of the neutrino luminosities is requiredfor astrophysical applications, because neutrinos carry away the most of energy lib-erated during the collapse as the main cooling source and can be primary observable.8 Y. Sekiguchi
Our results, in particular the duration and the peak luminosity of the neutrino bursts,agree approximately with those in the previous simulations. Again, no such goodagreement was reported in the previous simulation.
The shock stall-radius is ≈
80 km. This value is consistent with (althoughslightly smaller than) that in Liebend¨orfer et al. ( R stall ≈
85 km) and smallerthan that in Sumiyoshi et al. ( R stall ≈
100 km). This is likely because in ourleakage scheme, neutrino heating is not taken into account.To summarize, the results in the present simulation agree well with those in theprevious 1D Boltzmann simulations qualitatively. Quantitatively, the present resultsagree approximately with those in the previous 1D Boltzmann simulations. We canobtain approximately correct results with a not computationally expensive schemewithout solving the Boltzmann equation. Thus, the present code may be adopted,as a first step, to other multidimensional simulations such as the rotating stellarcollapse to a black hole and mergers of compact binaries.6.5.2. Comparison of the results after the convection sets inIn this section, we compare our results in the convective phase with those in thetwo-dimensional (2D) Newtonian simulations in whicha wide variety of approximations were adopted for the treatment of neutrinos.In the present simulation, we have found both the vigorous convective activities(the proto-neutron star convection) and the enhancement of neutrino luminositiesdue to the convection. These features agree approximately with those in the previous2D simulations with a fluid-like treatment of neutrinos and with radial ray-by-ray, gray flux-limited diffusion approximation of neutrino transfers.
Ina spherically symmetric, gray flux-limited diffusion scheme, by contrast, onlymildly active convection was found and no enhancement in the neutrino luminositieswas observed.Note that the transport of energy and lepton number by neutrinos can flattenthe negative gradients of entropy and lepton fraction, and as a result, the convec-tion will be suppressed. In purely hydrodynamic simulations without neutrino pro-cesses (using a postbounce core obtained in 1D simulations with neutrinos),the proto-neutron star convection is strongly activated. In the radial ray-by-ray sim-ulations, the transfer of neutrinos in the angular direction is not taken intoaccount and the stabilizing effect is underestimated, resulting in the proto-neutronstar convection with the enhancement of neutrino luminosities. In the sphericallysymmetric simulation, the transfer of neutrino in the angular direction is as-sumed to occur fast enough to make the neutrino distribution function sphericallysymmetric, and consequently, the stabilizing effect is overestimated.Recently, Buras et al. performed simulations with a modified ray-by-ray,multi-group scheme in which some part of the lateral components are included, andfound that the proto-neutron star convection indeed sets in but has minor effectson the enhancement of the neutrino luminosities. Dessart et al. performed simu-lations employing a 2D multi-group flux-limited diffusion scheme and found similarresults as in Buras et al. Thus, although the proto-neutron star convection indeedoccurs, its influence on enhancing the neutrino luminosities may be minor. The ull GR simulation with microphysics -100-50 0 50 100-10 0 10 20 30 40 50 60 70 A [ c m ] Time after bounce t b [ms] high resolutionlow resolution Fig. 7. Gravitational wave quadrupole amplitude A due to the prompt convection as a functionof post bounce time t b . The results for the finer grid resolution (solid curve) and for the coarsergrid resolution (the dotted curves) are shown together. strong convective activities and the enhancement of neutrino luminosities found inthe present simulation should be considered as the maximum ones.Note that it is in intermediate regions ( τ ν ∼
1) that the stabilizing effect dueto the neutrino transfer works efficiently: At higher density region with τ ν ≫ τ ν ≪
1, on the other hand, neutrinoscarry away the energy and the lepton number without interacting with the mat-ter. Therefore a careful and detailed treatment of the neutrino transfer is requiredto clarify the degree of the stabilizing effect and the convection, although such acomputationally expensive sumulation is beyond the scope of this paper.The present result that the proto-neutron star convection occurs qualitativelyagrees with the recent simulations with detailed neutrino transfer.
If simula-tions are perfomed keeping in mind that the stabilizing effect due to the neutrinotransfer is not taken into account in the present scheme, the present code will beacceptable to explore the the rotating stellar collapse to a black hole and mergers ofcompact binaries.6.6.
Gravitational radiation
Associated with the convective motions, gravitational waves are emitted. Thegravitational waveforms are computed using a quadrupole formula.
In quadrupoleformulae, only the +-mode of gravitational waves with l = 2 and m = 0 is nonzero0 Y. Sekiguchi -23 -22 -21 -20
100 1000 h c h a r kp c frequency [Hz] high resolutionlow resolution Fig. 8. The frequency spectra of the characteristic gravitational-wave strain due to the promptconvection. The results for the finer grid resolution (solid curve) and for the coarser gridresolution (the dotted curves) are shown together. in axisymmetric spacetime and is written as h quad+ = ¨ I zz ( t ret ) − ¨ I xx ( t ret ) r sin θ ≡ A ( t ) r sin θ, (6.12)where I ij denotes a quadrupole moment, ¨ I ij its second time derivative, and t ret aretarded time. In fully general relativistic and dynamical spacetime, there is nounique definition for the quadrupole moment and nor is for ¨ I ij . Following Shibataand Sekiguchi, we choose the simplest definition of the form I ij = Z ρ ∗ x i x j d x. (6.13)Then, using the continuity equation, the first time derivative can be written as˙ I ij = Z ρ ∗ ( v i x j + x i v j ) d x, (6.14)and ¨ I ij is computed by the finite differencing of the numerical result for ˙ I ij . In thefollowing, we present A , which provides the amplitude of a given mode measuredby an observer located in the most optimistic direction (in the equatorial plane). Wealso calculate the characteristic gravitational-wave strain, h char ( f ) ≡ s π Gc D dEdf , (6.15) ull GR simulation with microphysics Fig. 9. Snapshots of the contours of √− N / π , in the x - z plane at selected time slices. where dEdf = 8 π c G f (cid:12)(cid:12)(cid:12) ˜ A ( f ) (cid:12)(cid:12)(cid:12) (6.16)is the energy power spectra of the gravitational radiation and˜ A ( f ) = Z A ( t ) e πift dt. (6.17)Figure 7 shows A ( t ). Because the system is initially spherically symmetric,no gravitational radiation is emitted before the onset of the convection. When theproto-neutron star convection sets in at ≈
10 ms after the bounce, gravitationalwaves start to be emitted. The peak amplitudes are A ∼
100 cm. After the peak isreached, gravitational waves generated by the smaller-scale convective motions areemitted with A ≈
50 cm.Figure 8 shows the spectra of h char due to the convective motions. In contrastto the spectra due to the core bounce (e.g. Refs. 20) and 30)), there is no dominantpeak frequency in the power spectra. Instead, several maxima for the frequencyrange 100–1000 Hz are present. Note that for gravitational waves due to the corebounce, the characteristic peak frequency is associated with the bounce timescaleof the core. The effective amplitude of gravitational waves observed in the most2 Y. Sekiguchi optimistic direction is h char ≈ × − for an event at a distance of 10 kpc, whichis as large as that emitted at the bounce of rotating core collapse. To check that gravitational waves are indeed originated by the convective mo-tions, we calculate the frequency √− N / π (see Eq. (6.7)) as shown in Fig. 9. Thisfrequency is in good agreement with the gravitational-wave frequency, implying thatgravitational waves are indeed due to the convective activities.M¨uller and Janka investigated gravitational waves due to the convective mo-tion inside the proto-neutron star. It is interesting to compare our results with theirs.They adopted a post-bounce model of Hillebrandt. They put an inner boundaryat radius r in = 15 km and assumed the hydrostatic equilibrium there. They do notinclude neutrino transfer while a sophisticated EOS is adopted. They found quali-tatively similar results to ours. According to their results, the maximum amplitudeof the quadrupole mode is A ≈
100 cm, which agrees well with our results. Thespectrum of the gravitational-wave strain has several maxima for f = 50–500 Hzwith the maximum value of h char ≈ × − . The peaks in h char are distributed forhigher frequency side in our results probably due to the general relativistic effects.We note that a similar general relativistic effect is observed for gravitational wavesat the bounce phase. These facts show that for deriving quantitatively correctspectra of gravitational waves, fully general relativistic simulations are necessary.6.7.
Numerical accuracy
In Figs. 1–3 we show the results both in the higher resolution (solid curves) andin the lower resolution (dashed curves). The radial profiles of the two resolutionsare almost identical, showing that convergent results are obtained in the presentsimulation (see Fig. 4). In the time evolution of neutrino luminosities (see Fig. 3),the two results are almost identical before the convective activities set in. In thelater phase, on the other hand, the two results show slight disagreement. Becausethe convection and the turbulence can occur in an infinitesimal scale length, thesmaller-scale convection and turbulence are captured in the finer grid resolution.However, the influence of the grid resolution on the neutrino luminosities is minorbecause the convection and turbulence are strongly activated in the region above theneutrino sphere (see the contours of the electron fraction and the entropy in Fig.5). On the other hand, most of the neutrinos are emitted from the region inside theneutrino sphere (see the contour of the local neutrino energy emission rate in Fig.5). The effect of the grid resolution can be seen in gravitational waves. In Figs. 7and 8 we show the quadrupole mode A ( t ) and the characteristic strain h char ( f ) bothin the higher resolution (solid curves) and in the lower resolution (dashed curves).After the formation of the lepton-rich, hot finger at ≈
10 ms after the bounce (see § A ( t ) between the finerand the coarser grid resolutions becomes noticeable (see Fig. 7). The characteristicpeaks of h char ( f ) in higher frequencies ( f ∼ ull GR simulation with microphysics -50 -40 -30 -20 -10 0 10 20 30 40 50 B a r yon m a ss [ M s o l a r ] Time after bounce t b [ms] -5 -4 -3 -2 -1 H - c on s t r a i n t High resolutionLow resolution
Fig. 10. Evolution of the averaged violation of the Hamiltonian constraint (upper panel) and baryonmass conservation (lower panel). constraint are calculated, which is written as H = − ψ − (cid:20) ˜ ∆ψ − ψ R + 2 πρ h ψ + ψ A ij ˜ A ij − ψ K (cid:21) , (6.18)where ψ ≡ e φ , and ˜ ∆ denotes the Laplacian with respect to ˜ γ ij . In this paper, theaveraged violation is defined according to ERROR = 1 M ∗ Z ρ ∗ | V | d x, (6.19)where M ∗ is the rest-mass density of the core (see Eq. 6.1) V = ˜ ∆ψ − ψ R + 2 πEψ + ψ A ij ˜ A ij − ψ K | ˜ ∆ψ | + (cid:12)(cid:12)(cid:12) ψ R (cid:12)(cid:12)(cid:12) + 2 πρ h ψ + ψ A ij ˜ A ij + ψ K . (6.20)Namely, we use ρ ∗ as a weight factor for the average. This weight factor isintroduced to monitor whether the main bodies of the system (proto-neutron starsand inner cores), in which we are interested, are accurately computed or not.We display the time evolution of the Hamiltonian-constraint violation and theconservation of the baryon mass of the system in Fig. 10. Several discontinuous4 Y. Sekiguchi changes in the Hamiltonian-constraint violation and the conservation of the baryonmass originate from the regridding procedures in which some matters of the outerregion are discarded.Before the bounce, the baryon mass is well conserved and the Hamiltonian-constraint violation is very small as ∼ − . After the bounce, the violation ofthe baryon-mass-conservation and the Hamiltonian constraint is enhanced due tothe existence of shock waves where the hydrodynamic scheme becomes essentiallya first-order scheme. The convergence of the baryon-mass-conservation and theHamiltonian-constraint violation also becomes worse in the convective phase. How-ever, the degree of violation of the Hamiltonian constraint and the baryon-mass-conservation is small and we may believe that the numerical results obtained in thepaper are reliable. §
7. Summary and Discussion
Summary
In this paper, we present a fully general relativistic hydrodynamic code in whicha finite-temperature EOS and neutrino cooling are implemented for the first time.Because the characteristic timescale of weak interaction processes t wp ∼ | Y e / ˙ Y e | (WPtimescale) is much shorter than the dynamical timescale t dyn in hot dense matters,stiff source terms appear in the equations. In general, an implicit scheme may berequired to solve them. However, it is not clear whether implicit schemes do workor not in the relativistic framework. The Lorentz factor is coupled with the rest-mass density and the energy density. The specific enthalpy is also coupled with themomentum. Due to these couplings, it is not straightforward to recover the primitivevariables and the Lorentz factor from conserved quantities. Taking account of thesefacts, we proposed an explicit method to solve all the equations noting that thecharacteristic timescale of neutrino leakage from the system t leak is much longerthan t wp and is comparable to t dyn .By decomposing the energy-momentum tensor of neutrinos into the trapped-neutrino and the streaming-neutrino parts, the hydrodynamic equations can berewritten so that the source terms are characterized by the leakage timescale t leak (see Eqs. (3.8) and (3.5)). The lepton-number conservation equations, on the otherhand, include the source terms characterized by the WP timescale. Taking accountof these facts, limiters for the stiff source terms are introduced to solve the lepton-number conservation equations explicitly (see § § M ⊙ spherical modelwith the solar metallicity computed by Woosley et al. After the shock formationand propagation, the shock wave suffers from the severe reduction of its energy dueto neutrino burst emission when the shock wave passes the neutrino-sphere. Even- ull GR simulation with microphysics Y l above the neutrinosphere. Because such configuration is convectively unstable, vigorous convective mo-tions are induced. All these properties agree qualitatively with those by the resent2D Newtonian simulations. We also compared our results with those in the previous simulations. Beforethe convection sets in, we compare our result with those in the state-of-the-art 1DBoltzmann simulations in full general relativity.
As shown in this paper,the radial structure of the core and the neutrino luminosities agree qualitatively wellwith those in their simulations. Quantitatively, they also agree approximately withthe previous results.After the convection sets in, we compare our result with those in 2D Newtoniansimulations.
Our result that the proto-neutron star con-vection occurs agree qualitatively with that in the previous simulations.
However, quantitative properties show disagreement because the transfer of neutri-nos are not fully solved in the present scheme. Note that the transport of energy andlepton-number by neutrinos can flatten the negative gradients of entropy and leptonfraction, stabilizing the convection. Therefore the convective activities obtained inthe present simulation should be considered as the maximum ones.If we keep in mind the above facts and note the good agreements of the radialstructure and neutrino luminosities, the present implementation will be applied tosimulations of rotating core collapse to a black hole and mergers of binary neutronstars as a first step towards more sophisticated models. A detailed treatment of theneutrino transfer is required to determine the degree of stabilizing effect, but this isfar beyond the scope of this paper.Gravitational waves emitted by the convective motions are also calculated. Thegravitational-wave amplitude is ≈ × − for an event of the distance 10 kpc. Re-flecting the contributions of multi-scale eddies with characteristic overturn timescale1–10 ms, the energy power spectrum shows several maxima distributed in f ≈ in which a similarcalculation (but in Newtonian gravity) is performed. The maximum amplitude ofgravitational waves in our results agrees well with that in M¨uller and Janka. Theseveral maxima in the energy power spectrum are distributed at higher-frequencyside in our results due to the general relativistic effects, showing that fully general rel-ativistic simulations are necessary for the accurate calculation of gravitational-wavespectra.7.2. Discussions
Because the present implementation of the microphysics is simple and explicit,it has advantage that the individual microphysical processes can be easily improvedand sophisticated. For example, the neutrino emission via the electron capturecan be easily sophisticated as follows. To precisely calculate the electron capturerate, the complete information of the parent and daughter nuclei is required. InEOSs currently available, however, a representative single-nucleus average for thetrue ensemble of heavy nuclei is adopted. The representative is usually the most6
Y. Sekiguchi abundant nucleus. The problem in evaluating the capture rate is that the nucleiwhich cause the largest changes in Y e are neither the most abundant nuclei nor thenuclei with the largest rates, but the combination of the two. In fact, the mostabundant nuclei tend to have small rates because they are more stable than others,and the fraction of the most reactive nuclei tend to be small. Assuming thatthe nuclear statistical equilibrium (NSE) is achieved, the electron capture rates underthe NSE ensemble of heavy nuclei may be calculated for a given set of ( ρ , Y e , T ).Such a numerical rate table can be easily employed in the present implementation.Also, the neutrino cross sections can be improved. As summarized in Ref. 108),there are a lot of corrections to the neutrino opacities. Note that small changesin the opacities may result in much larger changes in the neutrino luminosities,because the neutrino energy emission rates depend strongly on the temperature, andthe temperature at the last scattering surface ( τ ν ∼ σT ∼
1) changes as T ∼ σ − / .Although the correction terms are in general very complicated, it is straightforwardto include the corrections into our code. Note that the corrections become moreimportant for higher neutrino energies. Therefore, the correction terms might play acrucial role in the collapse of population III stellar core and the formation of a blackhole, in which very high temperatures ( T >
100 MeV) will be achieved. A studyto explore the importance of these corrections in the case of black hole formation isongoing.As briefly described in the introduction, one of the main drawbacks in the presentimplementation of the neutrino cooling is that the transfer of neutrinos are notsolved. Although fully solving the transfer equations of neutrinos is far beyond thescope of this paper, there are a lot of rooms for improvements in the treatment ofthe neutrino cooling. For example, the relativistic moment formalism, inparticular the so-called M1 closure formalism, may be adopted. For this purpose, amore sophisticated treatment of the closure relation for P αβ is required. We plan toimplement a relativistic M1 closure formalism for the neutrino transfer in the nearfuture.To conclude, the present implementation of microphysics in fully general rela-tivistic, multidimensional code works well and has a wide variety of applications. Weare now in the standpoint where simulations of stellar core collapse to a black holeand merger of compact stellar binaries can be performed including microphysicalprocesses. Fruitful scientific results will be reported in the near future. acknowledgments The author thanks M. Shibata for valuable discussions and careful reading ofthe manuscript, and L. Rezzolla, and K. Sumiyoshi for valuable discussions. Healso thanks T. Shiromizu and T. Fukushige for their grateful aids. He thanks S. E.Woosley, A. Heger and T. A. Weaver for providing the presupernova core used in thepresent simulation as an initial condition. Numerical computations were performedon the NEC SX-9 at the data analysis center of NAOJ and on the NEC SX-8 atYITP in Kyoto University. This work is partly supported by the Grant-in-Aid of theJapanese Ministry of Education, Science, Culture, and Sport (21018008,21105511). ull GR simulation with microphysics Appendix A
Electron and positron captures
In this section, we briefly summarize our treatments of electron and positroncaptures which are based on Ref. 66) and give the explicit forms of γ ec νe , γ pc¯ νe , Q ec νe ,and Q pc¯ νe appeared in Eqs. (3.27), (3.28), and (3.30), for the purpose of convenience.A.1. The electron and positron capture rates γ ec νe and γ pc¯ νe The ’net’ electron fraction is written as Y e = Y − − Y + where Y − ( Y + ) denotesthe number of electrons (positrons) per baryon including pair electrons. Then theelectron-neutrino number emission rate by the electron capture and the electron-anti-neutrino number emission rate by the positron capture are given by γ local νe = − ˙ Y − = − ( ˙ Y f − + ˙ Y h − ) , (A.1) γ local¯ νe = − ˙ Y + = − ( ˙ Y f + + ˙ Y h + ) , (A.2)where the electron and positron capture rates are decomposed into two parts, captureon by free nucleons (with superscript f ) and on heavy nuclei (with superscript h ).In the following, we present the explicit forms of ˙ Y f − , ˙ Y f + , ˙ Y h − , and ˙ Y h + .A.2. Capture on free nucleons ˙ Y f The electron capture rate (including the contribution of the inverse reaction ofthe neutrino capture) on free nucleons ( ˙ Y f − ) is given by˙ Y f − = X n λ ν e c ,f − X p λ ec ,f , (A.3)where λ ec ,f is the specific electron capture rate on free protons, λ ν e c ,f is the specificelectron-neutrino capture rate on free neutrons, and X p and X n are the mass fractionof free proton and neutron, respectively. Based on a balance argument, one canshow that λ ν e c ,f is related to λ ec ,f by λ ν e c ,f = exp (cid:18) η νe − η e − δmk B T (cid:19) λ ec ,f , (A.4)where η νe and η e are the chemical potentials of electron neutrinos and electrons inunits of k B T and δm = ( m n − m p ) c . Furthermore, we use the following relation fornon-degenerate free nucleons, X n ≈ X p exp (cid:18) η n − η p + δmk B T (cid:19) , (A.5)where η n and η p are the chemical potentials of free neutrons and free protons in unitsof k B T . Then we obtain˙ Y f − = [exp ( η νe − η e + η n − η p ) − X p λ ec ,f . (A.6)The positron capture rate (including the contribution of the inverse reaction) onfree nucleons is similarly given by˙ Y f + = X p λ ¯ ν e c ,f − X n λ pc ,f = [exp ( η ¯ νe + η e + η p − η n ) − X n λ pc ,f , (A.7)8 Y. Sekiguchi where η ¯ νe is the chemical potential of electron-anti-neutrinos in units of k B T , λ pc isthe specific positron capture rate on free neutrons, and λ ¯ ν e c ,f is the specific electron-anti-neutrino capture rate on free protons.A.3. Capture on heavy nuclei ˙ Y h The electron capture rate (including the contribution of the inverse reaction ofthe neutrino capture) on a heavy nucleus of mass number A ( ˙ Y h − ) is given by ˙ Y h − = X D A λ ν e c ,h − X P A λ ec ,h , (A.8)where λ ec ,h is the specific electron capture rate on the parent nucleus (mass fraction X P ), λ ν e c ,h is the specific electron neutrino capture rate on the daughter nucleus(mass fraction X D ), and A is the atomic mass of the parent and daughter nuclei.In the present simulations, we set X D = X P = X A . Then, under the assumptionof nuclear statistical equilibrium, one may approximate the capture rate on heavynuclei as, ˙ Y h − ≈ [exp ( η νe − η e + η n − η p ) − X A A λ ec ,h . (A.9)Similarly, the positron capture rate (including the contribution of the inversereaction) on heavy nuclei ( ˙ Y h + ) is given by˙ Y h + = X D A λ ¯ ν e c ,h − X P A λ pc ,h ≈ [exp ( η ¯ νe + η e + η p − η n ) − X A A λ pc ,h . (A.10)A.4. The specific capture rate λ The specific electron and positron capture rates on free nucleons and on heavynuclei and are written in the same form as λ ec ,f = ln 2 h f t i ec ,f eff I ec ,f , λ pc ,f = ln 2 h f t i pc ,f eff I pc ,f , (A.11) λ ec ,h = ln 2 h f t i ec ,h eff I ec ,h , λ pc ,h = ln 2 h f t i pc ,h eff I pc ,h , (A.12)where I ec ,f and I pc ,f are the phase space factors for the electron and positron cap-tures on free electrons, and I ec ,h and I pc ,h are those on heavy nuclei. h f t i eff ’s arethe effective f t -values introduced by Fuller et al., which is essentially the same asthe square of the nuclear transition matrix.The phase space factors are given by I ec ,f = (cid:18) k B Tm e c (cid:19) Z ∞ η η ( η + ζ ec ,f )
11 + e η − η e (cid:20) −
11 + e η − η νe + ζ ec ,f (cid:21) dη, (A.13) I pc ,f = (cid:18) k B Tm e c (cid:19) Z ∞ η η ( η + ζ pc ,f )
11 + e η + η e (cid:20) −
11 + e η − η ¯ νe + ζ pc ,f (cid:21) dη, (A.14) I ec ,h = (cid:18) k B Tm e c (cid:19) Z ∞ η η ( η + ζ ec ,h )
11 + e η − η e (cid:20) −
11 + e η − η νe + ζ ec ,h (cid:21) dη, (A.15) ull GR simulation with microphysics I pc ,h = (cid:18) k B Tm e c (cid:19) Z ∞ η η ( η + ζ pc ,h )
11 + e η + η e (cid:20) −
11 + e η − η ¯ νe + ζ pc ,h (cid:21) dη, (A.16)where ζ ec ,f , ζ pc ,f , ζ ec ,h , and ζ pc ,h are the nuclear mass-energy differences for electroncapture and positron capture in units of k B T . The nuclear mass-energy differencesfor capture on free nuclei are given by ζ ec ,f = − ζ pc ,fn ≈ η p − η n . (A.17)We follow Fuller et al. for the nuclear mass-energy differences for capture on heavynuclei: In the case of N <
40 or
Z >
20 (referred to as ’unblocked’ case), we set ζ ec ,h = − ζ pc ,h ≈ η p − η n . (A.18)In the case of N ≥
40 or Z ≤
20 (referred to as ’blocked’ case), on the other hand,we set ζ ec ,h ≈ η p − η n − k B T , (A.19) ζ pc ,h ≈ − η p + η n + 5(MeV) k B T . (A.20)Then, the threshold value of the electron and positron captures is given by η = m e c / ( k B T ) for ζ > − m e c / ( k B T ) and η = | ζ | for ζ < − m e c / ( k B T ) where wehave dropped the superscripts ’ec’, ’pc’, ’ f ’, and ’ h ’ in ζ for simplicity.The effective f t -value of electron or positron capture on free nuclei is given by(e.g. Ref. 66) log h f t i ec ,f eff = log h f t i pc ,f eff ≈ . . (A.21)We follow Fuller et al. for the effective f t -value of capture on heavy nuclei, whoproposed to uselog h f t i ec ,h eff ≈ . η e < | ζ ec ,h | . η e > | ζ ec ,h | . . T blocked , (A.22)log h f t i pc ,h eff ≈ . η e < | ζ pc ,h | . η e > | ζ pc ,h | . . T blocked , (A.23)where T = T / (10 K ). In this expression, the thermal unblocking effect is readilytaken into account. In the thermal unblocking, it costs ≈ .
13 MeV to pull a neutronout of a filled orbital 1 f / and place it in the gd -shell. A.5.
Energy emission rates Q ec νe and Q pc¯ νe The neutrino energy emission rates associated with electron and positron cap-tures in units of m e c s − are given by π ec = ln 2 J ec h f t i eceff , π pc = ln 2 J pc h f t i pceff , (A.24)0 Y. Sekiguchi where the phase space factors are given by J ec = (cid:18) k B Tm e c (cid:19) Z ∞ η η ( η + ζ ec )
11 + e η − η e (cid:20) − − e η − η νe + ζ ec (cid:21) dη, (A.25) J pc = (cid:18) k B Tm e c (cid:19) Z ∞ η η ( η + ζ pc )
11 + e η + η e (cid:20) − − e η − η ¯ νe + ζ pc (cid:21) dη. (A.26)In Eqs. (A.24)–(A.26), we have dropped the superscripts ’ f ’ and ’ h ’ in π ec , π pc , J ec , J pc , h f t i eceff , h f t i pceff , ζ ec , and ζ pc for simplicity.The average energy of electron neutrinos produced by electron and positroncaptures is defined, in units of m e c , as h ǫ νe i ec = J ec I ec , h ǫ ¯ νe i pc = J ec I pc . (A.27)Then, the local neutrino energy emission rates by the electron and positron capturesper unit volume is given by Q ec νe = ρm u (cid:20) X p h ǫ νe i ec ,f λ ec ,f + X A A h ǫ νe i ec ,h λ ec ,h (cid:21) , (A.28) Q pc¯ νe = ρm u (cid:20) X n h ǫ ¯ νe i pc ,f λ pc ,f + X A A h ǫ ¯ νe i pc ,h λ pc ,h (cid:21) . (A.29) Appendix B
Neutrino pair processes
In this section, we briefly summarize our treatment of pair processes of neutrinoemission and give the explicit forms of γ pair ν e ¯ ν e , γ plas ν e ¯ ν e , γ Brems ν e ¯ ν e , γ pair ν x ¯ ν x , γ plas ν x ¯ ν x , γ Brems ν x ¯ ν x , Q pair ν e ¯ ν e , Q plas ν e ¯ ν e , Q Brems ν e ¯ ν e , Q pair ν x ¯ ν x , Q plas ν x ¯ ν x , and Q Brems ν x ¯ ν x appeared in Eqs. (3.27), (3.28), (3.29) and(3.30), for the purpose of convenience.B.1. Electron-positron pair annihilation
We follow Cooperstein et al. for the rate of pair creation of neutrinos by theelectron-positron pair annihilation. The number emission rate of ν e or ¯ ν e by theelectron-positron pair annihilation can be written as γ pair ν e ¯ ν e = m u ρ C pair ν e ¯ ν e π σ cm e c ( k B T ) ( ~ c ) F ( η e ) F ( − η e ) h block i pair ν e ¯ ν e , (B.1)where σ ≈ . × − cm − and C pair ν e ¯ ν e = ( C V − C A ) + ( C V + C A ) with C V = + 2 sin θ W and C A = . The Weinberg angle is given by sin θ W ≈ .
23. Usingthe average energy of neutrinos produced by the pair annihilation, h ǫ ν e ¯ ν e i pair = k B T (cid:18) F ( η e ) F ( η e ) + F ( − η e ) F ( − η e ) (cid:19) , (B.2) ull GR simulation with microphysics h block i pair ν e ¯ ν e is evaluated as h block i pair νe ≈ (cid:20) (cid:18) η νe − h ǫ ν e ¯ ν e i pair k B T (cid:19)(cid:21) − (cid:20) (cid:18) η ¯ νe − h ǫ ν e ¯ ν e i pair k B T (cid:19)(cid:21) − . (B.3)The associated neutrino energy emission rate by the pair annihilation is given by Q pair ν e ¯ ν e = ρm u γ pair ν e ¯ ν e h ǫ ν e ¯ ν e i pair . (B.4)Similarly, the number emission rate of ν x or ¯ ν x by the electron-positron pairannihilation and the associated energy emission rate are given by γ pair ν x ¯ ν x = m u ρ C pair ν x ¯ ν x π σ cm e c ( k B T ) ( ~ c ) F ( η e ) F ( − η e ) h block i pair ν x ¯ ν x , (B.5) Q pair ν x ¯ ν x = ρm u γ pair ν x ¯ ν x h ǫ ν x ¯ ν x i pair , (B.6)where C ν x ¯ ν x = ( C V − C A ) + ( C V + C A − . The average neutrino energy and theblocking factor are given by h ǫ ν x ¯ ν x i pair = h ǫ ν e ¯ ν e i pair (B.7)and h block i pair νe ≈ (cid:20) (cid:18) η νx − h ǫ ν x ¯ ν x i pair k B T (cid:19)(cid:21) − (cid:20) (cid:18) η ¯ νx − h ǫ ν x ¯ ν x i pair k B T (cid:19)(cid:21) − , (B.8)where note that η ¯ νx = η νx .B.2. Plasmon decay
We follow Ruffert et al. for the rate of pair creation of neutrinos by the decayof transversal plasmons. The number emission rate of ν e or ¯ ν e can be written as γ plas ν e ¯ ν e = m u ρ C V π α fine σ cm e c ( k B T ) ( ~ c ) γ p e − γ p (1 + γ p ) h block i plas ν e ¯ ν e , (B.9)where α fine ≈ /
137 is the fine-structure constant and γ p ≈ p ( α fine / π )( π + 3 η e ).The blocking factor is approximately given by h block i plas ν e ¯ ν e ≈ (cid:20) (cid:18) η νe − h ǫ ν e ¯ ν e i plas k B T (cid:19)(cid:21) − (cid:20) (cid:18) η ¯ νe − h ǫ ν e ¯ ν e i plas k B T (cid:19)(cid:21) − , (B.10)where h ǫ ν e ¯ ν e i plas = k B T γ p γ p ! (B.11)is the average energy of neutrinos produced by the plasmon decay. The associatedneutrino energy emission rate is given by Q plas ν e ¯ ν e = ρm u γ plas ν e ¯ ν e h ǫ ν e ¯ ν e i plas . (B.12)2 Y. Sekiguchi
Similarly, the number emission rate of ν x or ¯ ν x by the plasmon decay and theassociated energy emission rate are given by γ plas ν x ¯ ν x = m u ρ ( C V − π α fine σ cm e c ( k B T ) ( ~ c ) γ p e − γ p (1 + γ p ) h block i plas ν x ¯ ν x , (B.13) Q plas ν x ¯ ν x = ρm u γ plas ν x ¯ ν x h ǫ ν x ¯ ν x i plas , (B.14)where the average neutrino energy is h ǫ ν x ¯ ν x i plas = h ǫ ν e ¯ ν e i plas and the blocking factoris given by h block i pair νe ≈ (cid:20) (cid:18) η νx − h ǫ ν x ¯ ν x i pair k B T (cid:19)(cid:21) − (cid:20) (cid:18) η ¯ νx − h ǫ ν x ¯ ν x i pair k B T (cid:19)(cid:21) − . (B.15)B.3. Nucleon-nucleon bremsstrahlung
We follow Burrows et al. for the rate of pair creation of neutrinos by thenucleon-nucleon bremsstrahlung radiation. They derived the neutrino energy emis-sion rate associated with the pair creation of ν x or ¯ ν x by the nucleon-nucleonbremsstrahlung radiation without the blocking factor: Q Brems , ν x ¯ ν x = 3 . × ζ Brems (cid:18) X n + X p + 283 X n X p (cid:19) ρ (cid:18) k B Tm e c (cid:19) . h ǫ ν x ¯ ν x i Brems , (B.16)where ζ Brems ∼ . h ǫ ν x ¯ ν x i Brems ≈ . k B T. (B.17)We multiply Q Brems , ν x ¯ ν x by the blocking factor, h block i Brems ν x ¯ ν x ≈ (cid:20) (cid:18) η νx − h ǫ ν x ¯ ν x i Brems k B T (cid:19)(cid:21) − (cid:20) (cid:18) η ¯ νx − h ǫ ν x ¯ ν x i Brems k B T (cid:19)(cid:21) − , (B.18)to obtain the ’blocked’ neutrino energy emission rate Q Brems ν x ¯ ν x = Q Brems , ν x ¯ ν x h block i Brems ν x ¯ ν x . (B.19)The number emission rate of ν x or ¯ ν x is readily given by γ Brems ν x ¯ ν x = 3 . × ζ Brems (cid:18) X n + X p + 283 X n X p (cid:19) m u ρ (cid:18) k B Tm e c (cid:19) . h block i Brems ν x ¯ ν x . (B.20)Noting that the weak interaction coefficients of the bremsstrahlung radiationare (1 − C V ) + (1 − C A ) for the pair creation of ν x ¯ ν x and C V + C A for the paircreation of ν e ¯ ν e , the number emission rate and the associated energy emission ratefor ν e or ¯ ν e may be written as γ Brems ν e ¯ ν e = C V + C A (1 − C V ) + (1 − C A ) γ Brems ν x ¯ ν x , (B.21) Q Brems ν e ¯ ν e = C V + C A (1 − C V ) + (1 − C A ) Q Brems ν x ¯ ν x . (B.22) ull GR simulation with microphysics Appendix C
Neutrino diffusion rates
We follow Ref. 62) for the diffusive neutrino-number emission rate γ diff ν and theassociated energy emission rate Q diff ν appeared in Eqs. (3.25) and (3.26). and presentthe explicit forms of them in the following for convenience. An alternative definitionof the diffusion rates are found in Ref. 63).C.1. Neutrino diffusion rates
To calculate the neutrino diffusion rates γ diff ν and Q diff ν , we first define neutrinodiffusion time. In this paper, we consider cross sections for scattering on nuclei( σ sc νA ), on free protons ( σ sc νp ), and on free neutrons ( σ sc νn ), and for absorption on freenucleons σ ab νn for electron neutrinos and σ ab νp for electron anti-neutrinos.Ignoring the higher order correction terms in neutrino energy E ν , these neutrinocross sections can be written in general as σ ( E ν ) = E ν ˜ σ, (C.1)where ˜ σ is a ’cross section’ in which E ν dependence is factored out. In practice, thecross sections contain the correction terms which cannot be expressed in the formof Eq. (C.1). We take account of these correction terms, approximating neutrino-energy dependence by temperature dependence according to E ν ≈ k B T F ( η ν ) F ( η ν ) . (C.2)The opacity is written as κ ( E ν ) = X κ i ( E ν ) = E ν X ˜ κ i = E ν ˜ κ, (C.3)and optical depth is calculated by τ ( E ν ) = Z κ ( E ν ) ds = E ν Z ˜ κds = E ν ˜ τ . (C.4)Then, we define neutrino diffusion time by T diff ν ( E ν ) ≡ a diff ∆x ( E ν ) c τ ( E ν ) = E ν a diff ˜ τ c ˜ κ = E ν ˜ T diff ν , (C.5)where the distance parameter ∆x ( E ν ) is given by ∆x ( E ν ) = τ ( E ν ) κ ( E ν ) . (C.6)Note that ˜ T diff ν can be calculated only using matter quantities. a diff is a parameterwhich controls how much neutrinos diffuse outward and we chose it to be 3 followingRef. 63). For a larger value of a diff , corresponding neutrino emission rate due todiffusion becomes smaller.4 Y. Sekiguchi
Finally, we define the neutrino diffusion rates by γ diff ν ≡ m u ρ Z n ν ( E ν ) T diff ν ( E ν ) dE ν = 1 a diff m u ρ πcg ν ( hc ) ˜ κ ˜ τ T F ( η ν ) , (C.7) Q diff ν ≡ Z E ν n ν ( E ν ) T diff ν ( E ν ) dE ν = 1 a diff πcg ν ( hc ) ˜ κ ˜ τ T F ( η ν ) . (C.8)C.2. Summary of cross sections
In this subsection, we briefly summarize the cross sections adopted in the presentneutrino leakage scheme.C.2.1. Neutrino nucleon scatteringThe total ν - p scattering cross section σ p for all neutrino species is given by σ sc νp = σ (cid:18) E ν m e c (cid:19) (cid:2) ( C V − + 3 g A ( C A − (cid:3) , (C.9)where g A is the axial-vector coupling constant g A ≈ − .
26. On the other hand, thetotal ν − n scattering cross section σ n is σ sc νn = σ (cid:18) E ν m e c (cid:19) (cid:2) g A (cid:3) . (C.10)C.2.2. Coherent scattering of neutrinos on nucleiThe differential cross section for the ν - A neutral current scattering is writtenas dσ sc A dΩ = σ π (cid:18) E ν m e c (cid:19) A [ WC FF + C LOS ] hS ion i (1 + cos θ ) , (C.11)where θ is the azimuthal angle of the scattering and W = 1 − ZA (1 − θ W ) . (C.12) hS ion i , C LOS , and C F F are correction factors due to the Coulomb interaction betweenthe nuclei, due to the electron polarization, and due to the finite size of heavynuclei.
Because it is known that the correction factor C LOS is important only forlow neutrino energies, we consider only hS ion i and C F F .The correction factor due to the Coulomb interaction between the nuclei is givenby hS ion i = 34 Z − d cos θ (1 + cos θ )(1 − cos θ ) S ion . (C.13)Itoh et al. presented a detailed fitting formula for the correction factor. However,the fitting formula is very complicated, we use a simple approximation based onRef. 117). S ion can be written approximately as S ion ≈ ( qa I ) Γ + f ( Γ )( qa I ) (C.14) ull GR simulation with microphysics -5 -4 -3 -2 -1 < S i on > log x simple approximationItoh et al. (2004) Fig. 11. Comparison of the ion-ion correction term, as a function of x ≡ E ν a I / ( ~ c ), between oursimplified estimation and a detailed fitting formula by Itoh et al. From the top to the bottom,the curves are for Γ = 1, 2, 5, 10, 20, 40, 80, 125, 160. for ( qa I ≪ q = (2 E ν / ~ c ) sin( θ/ a I = (4 πn A / − / is the ion-sphereradius, n A is the number density of a nucleus, and Γ = ( Ze ) / ( a I k B T ) is the con-ventional parameter that characterizes the strongness of the Coulomb interaction. f ( Γ ) is given by f ( Γ ) ≈ . − . Γ + 0 . Γ / + 0 . Γ − / . (C.15)The integration approximately gives for x ≡ E ν a I / ( ~ c ) < hS ion i ≈
16 1
Γ x − f ( Γ ) Γ x + 1135 ( f ( Γ )) Γ x − f ( Γ )) Γ x + 12268 ( f ( Γ )) Γ x . (C.16)To use this expression for the case of x ≥
1, we set the maximum value as hS ion i =max(1 , hS ion i ) where hS ion i = 1 corresponds to the case without the correction. Notethat x can be larger than unity ∗ ) . In this case, the simple approximation on averageunderestimates the effect of the Coulomb interaction between the nuclei (see Fig.11). Figure 11 compares the correction term as a function of the parameter x calcu-lated by our simple approximation and by a detailed fitting formula by Itoh et al. For smaller values of Γ , the disagreements become prominent. Note that the typicalvalue of Γ inside the collapsing core with T ∼ ρ ∼ g/cm , A = 56 and Z = 26 ( Fe) is Γ ∼ ∗ ) x ∼ ( E ν / ρ /A ) − / where ρ is the rest-mass density in units of 10 g/cm Y. Sekiguchi
C.2.3. Absorption on free neutronsThe total cross section of the absorption of electron neutrinos on free neutronsis given by σ ab n = σ (cid:18) g A (cid:19) (cid:18) E ν + ∆ np m e c (cid:19) (cid:20) − (cid:18) m e c E ν + ∆ np (cid:19)(cid:21) W M , (C.17)where ∆ np = m n c − m p c , and W M is the correction for weak magnetism and recoilwhich is approximately given by W M ≈ . E ν m n c . (C.18)Similarly, the total cross section of the absorption of electron anti-neutrinos onfree protons is given by σ ab p = σ (cid:18) g A (cid:19) (cid:18) E ¯ ν − ∆ np m e c (cid:19) (cid:20) − (cid:18) m e c E ¯ ν − ∆ np (cid:19)(cid:21) W ¯ M , (C.19)where W ¯ M is approximately given by W ¯ M ≈ − . E ¯ ν m p c . (C.20) References
1) K. Hirata, et al., Phys. Rev. Lett. (1987), 1490.2) R. M. Bionta et al., Phys. Rev. Lett. (1987), 1494.3) W. D. Arnett, J. N. Bahcall, R. P. Kirshner, and S. E. Woosley, Ann. Rev. Astron.Astrophys. (1989), 629.4) A. Abramovici, et al., Sience (1992), 325.5) F. Acernese, et al., Class. Quantum Grav. (2004), S385.6) B. Willke, et al., Class. Quantum Grav. (2002), 1377.7) M. Ando and the TAMA Collaboration, Class. Quantum Grav. (2002), 1409.8) K. Kotake, K. Sato, and K. Takahashi, Rep. Prog. Phys. (2006), 971.9) H.-Th. Janka, K. Langanke, A. Marek, G. Mart´ınez-Pinedo, and B. M¨uller, Phys. Rep. (2007), 38.10) H. A. Bethe, Rev. Mod. Phys. (1990), 801.11) M. Rampp and H.-Th. Janka, Astrophys. J. (2000), L33.12) A. Mezzacappa, M. Liebend¨orfer, O. E. B. Messer, W. R. Hix, F.-K. Thielemann, andS. W. Bruenn, Phys. Rev. Lett. (2001), 1935.13) T. A. Thompson, A. Burrows, and P. A. Pinto, Astrophys. J. (2003), 434.14) M. Rampp and H.-Th. Janka, Astron. Astrophys. (2002), 361.15) M. Liebend¨orfer, A. Mezzacappa, F.-K. Thielemann, O. E. B. Messer, W. R. Hix, andS. W. Bruenn, Phys. Rev. D (2001), 103004.16) K. Sumiyoshi, S. Yamada, H. Suzuki, H. Shen, S. Chiba, and H. Toki, Astrophys. J. (2005), 922.17) F. S. Kitaura, H.-Th. Janka, and W. Hillebrandt, Astron. Astrophys. (2006), 345.18) A. Marek and T. Th. Janka, Astrophys. J. (2009), 664.19) A. Burrows, E. Livne, L. Dessart, C. D. Ott, and J. Murphy, Astrophys. J. (2006),878.20) H. Dimmelmeier, J. A. Font, and E. M¨uller, Astron. Astrophys. (2002), 917; Astron.Astrophys. (2002), 523. ull GR simulation with microphysics
21) R. M¨onchmeyer, G. Sh¨afer, E. M¨uller, and R. E. Kates, Astron. Astrophys. (1991),417.22) A. Burrows and J. Hayes, Phys. Rev. Lett. (1996), 352.23) E. M¨uller and H.-Th. Janka, Astron. Astrophys. (1997), 140.24) T. Zwerger and E. M¨uller, Astron. Astrophys. (1997), 209.25) K. Kotake, S. Yamada, and K. Sato, Phys. Rev. D (2003), 044023; K. Kotake, S. Ya-mada, K. Sato, K. Sumiyoshi, H. Ono, and H. Suzuki, Phys. Rev. D (2004), 124004.26) C. D. Ott, A. Burrows, E. Livne, and R. Walder, Astrophys. J. (2004), 834.27) E. M¨uller, M. Rampp, R. Buras, H.-Th. Janka, and D. H. Shoemaker, Astrophys. J. (2004), 221.28) C. L. Fryer, D. E. Holz, and S. A. Hughes, Astrophys. J. (2004), 288.29) M. Shibata and Y. Sekiguchi, Phys. Rev. D (2004), 084024; Phys. Rev. D (2005),024014.30) Y. Sekiguchi and M. Shibata, Phys. Rev. D (2005), 084013.31) P. Cerd´a-Dur´an, G. Faye, H. Dimmelmeier, J. A. Font, J. Ma. Ib´a˜nez, E. M¨uller, andG. Sch¨afer, Astron. Astrophys. (2005), 1033.32) M. Shibata, Y. T. Liu, S. L. Shapiro, B. C. Stephens, Phys. Rev. D (2006), 104026.33) C. D. Ott, H. Dimmelmeier, A. Marek, H.-T. Janka, I. Hawke, B. Zink, and E. Schnetter,Phys. Rev. Lett. (2007), 261101.34) H. Dimmelmeier, C. D. Ott, H.-Th. Janka, A. Marek, and E. M¨uller, Phys. Rev. Lett. (2007), 251101.35) H. Shen, H. Toki, K. Oyamatsu, and K. Sumiyoshi, Nucl. Phys. A (1998), 435; Prog.Theor. Phys. (1998), 1013.36) M. Liebend¨orfer, Astrophys. J. (2005), 1042.37) M. J. Rees, in The future of theoretical physics and cosmology: celebrating Stephen Hawk-ing’s 60th birthday , edited by G. W. Gibbons et al. (Cambridge University Press, Cam-bridge, 2003), p.217, (astro-ph/0401365).38) J. E. McClintock and R. A. Remillard, in
Compact Stellar X-ray Sources , edited byW. H. G. Lewin and van der Klis (Cambridge University Press, Cambridge, 2006),(astro-ph/0306213).39) K. Sumiyoshi, S. Yamada, H. Suzuki, and S. Chiba, Phys. Rev. Lett. (2006), 091101;K. Sumiyoshi, S. Yamada, and H. Suzuki Astrophys. J. (2007), 382.40) K. Nakazato, K. Sumiyoshi, and S. Yamada, Astrophys. J. (2007), 1140.41) Y. Sekiguchi and M. Shibata, Prog. Theor. Phys. (2007), 1029.42) Y. T. Liu, S. L. Shapiro, and B. C. Stephens, Phys. Rev. D (2007), 084017.43) S. R. Kulkarni et al., Nature (1998), 663; T. J. Galama et al, Nature (1998), 670.44) J. Hjorth et al., Nature (2003), 847; K. Z. Stanek et al., Astrophys. J. (2003),L17; K. S. Kawabata et al., Astrophys. J. (2003), L19.45) M. Modjaz et al, Astrophys. J. (2006), L21; J. Sollerman et al, Astron. Astrophys. (2006), 503.46) S. E. Woosley, Astrophys. J. (1993), 273; B. Paczy´nski, Astrophys. J. (1998),L45; A. MacFadyen and S. E. Woosley, Astrophys. J. (1999), 262.47) M. Shibata and S. L. Shapiro, Astrophys. J. (2002), L39.48) Y. Sekiguchi and M. Shibata, Phys. Rev. D (2004), 084005.49) S. Yamada, H. Th. Janka, and H. Suzuki, Astron. Astrophys. (1999), 533.50) M. Liebend¨orfer et al., Astrophys. J. Suppl. (2004), 263.51) R. I. Epstein and C. J. Pethick, Astrophys. J. (1981), 1003.52) K. A. van Riper and J. M. Lattimer, Astrophys. J. (1981), 270.53) K. A. van Riper, Astrophys. J. (1982), 793.54) S. A. Bludman, I. Lichtenstadt, and G. Hayden, Astrophys. J. (1982), 661.55) A. Ray, S. M. Chitre, and K. Kar, Astrophys. J. (1984), 766.56) E. A. Baron, J. Cooperstein, and S. Kahana, Nucl. Phys. A (1985), 744.57) J. Cooperstein, Phys. Rep. (1988), 95.58) K. Kotake, S. Yamada, and K. Sato, Astrophys. J. (2003), 304.59) C. D. Ott et al., Astrophys. J. (2008), 1069.60) H. A. Bethe, Rev. Mod. Phys. (1990), 801.61) Y. Sekiguchi 2007 PhD thesis (University of Tokyo).62) S. Rosswog and M. Liebend¨orfer, Mon. Not. R. Astron. Soc. (2003), 673. Y. Sekiguchi
63) M. Ruffert, H.-Th. Janka, and G. Sch¨afer, Astron. Astrophys. (1996), 532.64) H. A. Bethe and J. R. Wilson, Astrophys. J. (1985), 14.65) S. W. Bruenn, Astrophys. J. Suppl. (1985), 771.66) G. M. Fuller, W. A. Fowler, and M. J. Newman, Astrophys. J. (1985), 1.67) J. Cooperstein and J. Wambach, Nucl. Phys. A (1984), 591.68) Bowers R L and Wilson J R 1982 Astrophys. J. Suppl. (1987), 744.70) A. Mezzacappa and S. W. Bruenn, Astrophys. J. (1993), 405.71) E. Livne et al., Astrophys. J. (2004), 277.72) R. Buras et al., Astron. Astrophys. (2006), 1049.73) A. Burrows et al., Astrophys. J. (2007), 416.74) A. Mezzacappa and . E. B. Messer, J. Comp. Appl. Math. (1999), 281.75) C. W. Misner and D. H. Sharp, Phys. Rev. (1964), B571.76) M. Liebend¨orfer, S. C. Whitehouse and T. Fischer, Astrophys. J. (2009), 1174.77) M. Shibata, Y. Sekiguchi, and R. Takahashi Prog. Theor. Phys. (2007), 257.78) J. Cooperstein et al., Astrophys. J. (1986), 653.79) A. Burrows et al., Nucl. Phys. A (2006), 356.80) Y. Sugahara and H. Toki, Nucl. Phys. A (1994), 557.81) e.g., R. Brockmann and R. Machliedt, Phys. Rev. C (1990), 1965.82) J. M. Lattimer and F. D. Swesty, Nucl. Phys. A (1991), 331.83) C. Ishizuka, A. Ohnishi, K. Tsunakihara, K. Sumiyoshi, and S, Yamada, J. Phys. G (2008), 085201.84) J. P. Cox and R. T. Giuli, Principles of Stellar Structure volume 2. , (Gordon & Breach,New York, 1968).85) J. W. York, Kinematics and Dynamics of General Relativity, in
Sources of gravitationalradiation , edited by Smarr, L., (Cambridge University Press, Cambridge, 1979).86) M. Shibata and T. Nakamura, Phys. Rev. D (1995), 5428.87) T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D (1999), 024007.88) M. Alcubierre et al, Int. J. Mod. Phys. D (2001), 273.89) M. Shibata, Phys. Rev. D (2003), 024033.90) K. Kiuchi et al., Phys. Rev. D (2009), 064037.91) M. Alcubierre and B. Br¨ugmann, Phys. Rev. D (2001), 104006.92) M. Shibata, Astrophys. J. (2003), 992.93) A. Kurganov and E. Tadmor, J. Comp. Phys. (2000), 214.94) S. E. Woosley, A. Heger, and T. A. Weaver, Rev. Mod. Phys. (2002), 1015.95) S. L. Shapiro and S. A. Teukolsky, Black Holes, White Dwarfs, and Neutron Stars (WileyInterscience, New York, 1983).96) M. Liebend¨orfer et al., Astrophys. J. (2005), 840.97) J. M. Lattimer and T. J. Mazurek, Astrophys. J. (1981), 955.98) A. Burrows and B. A. Fryxell, Astrophys. J. (1993), L33.99) M. Herant et al., Astrophys. J. (1994), 339.100) A. Burrwos, J. Hayes, and B. .A Fryxell, Astrophys. J. (1995), 830.101) W. Keil, H.-Th. Janka, and E. M¨uller, Astrophys. J. (1996), 111.102) A. Mezzacappa et al., Astrophys. J. (1998), 848.103) L. Dessart et al., Astrophys. J. (2006), 534.104) M. Shibata and Y. Sekiguchi, Phys. Rev. D (2003), 104020.105) E. E. Flanagan and S. A. Hughes, Phys. Rev. D (1998), 4535.106) W. Hillebrandt, in High Energy Phenomena around Collapsed Stars , ed. F. Pacini, (Rei-del, Dodrecht, 1987), p.73.107) M. B. Aufderheide et al., Astrophys. J. (1994), 389.108) C. J. Horowitz, Phys. Rev. D (2002), 043001.109) J. L. Anderson and E. A. Spiegel, Astrophys. J. (1972), 127.110) K. S. Thorne, Mon. Not. R. Astron. Soc. (1981), 439.111) N. Itoh et al. Astrophys. J. Suppl. (1996), 411.112) C. J. Horowitz, Phys. Rev. D (1997), 4577.113) L. B. Leinson, V. N. Oraevsky, and V. B. Semikoz, Phys. Lett. B (1988), 80.114) D. L. Tubbs and D. N. Schramm, Astrophys. J. (1975), 467; A. Burrows,T. J. Mazurek, and J. M. Lattimer, Astrophys. J. (1981), 325 ull GR simulation with microphysics (1997), 7529.116) N. Itoh et al. Astrophys. J. (2004), 1041.117) N. Itoh, Prog. Theor. Phys. (1975), 1580.118) A. Burrows and T. A. Thompson, in Stellar Collapse , ed. C. L. Fryer, (Kluwer Academic,Dordrecht, 2004); A. Burrows, S. Reddy, and T. A. Thompson, Nucl. Phys. A777