An implementation of the microphysics in full general relativity : General relativistic neutrino leakage scheme
aa r X i v : . [ a s t r o - ph . H E ] S e p An implementation of the microphysics in fullgeneral relativity : General relativistic neutrinoleakage scheme
Yuichiro Sekiguchi
Division of Theoretical Astronomy, National Astronomical Observatory of Japan,Mitaka, Tokyo 181-8588, JapanE-mail: [email protected]
Abstract.
Performing fully general relativistic simulations taking account of microphysicalprocesses (e.g., weak interactions and neutrino cooling) is one of long standing problemsin numerical relativity. One of main difficulties in implementation of weak interactionsin the general relativistic framework lies on the fact that the characteristic timescaleof weak interaction processes (the WP timescale, t wp ∼ | Y e / ˙ Y e | ) in hot dense mattersis much shorter than the dynamical timescale ( t dyn ). Numerically this means that stiff source terms appears in the equations so that an implicit scheme is in generalnecessary to stably solve the relevant equations. Otherwise a very short timestep(∆ t < t wp ≪ t dyn ) will be required to solve them explicitly, which is unrealistic inthe present computational resources. Furthermore, in the relativistic framework, theLorentz factor is coupled with the rest mass density and the energy density. Thespecific enthalpy is also coupled with the momentum. Due to these couplings, it is verycomplicated to recover the primitive variables and the Lorentz factor from conservedquantities. Consequently, it is very difficult to solve the equations implicitly in the fullygeneral relativistic framework. At the current status, no implicit procedure have beenproposed except for the case of the spherical symmetry. Therefore, an approximate,explicit procedure is developed in the fully general relativistic framework in this paperas an first implementation of the microphysics toward a more realistic sophisticatedmodel. The procedure is based on the so-called neutrino leakage schemes which isbased on the property that the characteristic timescale in which neutrinos leak outof the system (the leakage timescale, t leak ) is much longer than the WP timescale.In the previous leakage schemes, however, the problems of the stiff source terms areavoided in an artificial manner. In this paper, I present a detailed neutrino leakagescheme and a simple and stable method for solving the equations explicitly in the fullygeneral relativistic framework. The drawback of the artificial treatment of the stiffsource terms is improved. I also perform a test simulation to check the validity of thepresent method, showing that it works fairly well. eneral relativistic leakage scheme
1. Introduction
Numerical relativity is the unique and powerful tool to explore dynamical phenomenain which strong gravity plays important roles. Stellar core collapse and mergers ofcompact star binaries are among the most important and interesting events in thefield. In theoretical view points, performing simulations of these phenomena is oneof challenging problems because a rich diversity of physics has to be taken into account.All four known forces of nature are involved and play important roles during the collapse.General relativistic gravity plays essential roles in formation of a black hole. Note alsothat general relativity may play important roles in supernova explosion as previouspioneering works [1, 2] showed. The weak interactions and emission of neutrinos governenergy and lepton-number losses, and hence driving the thermal and chemical evolutionsof the system. The strong interactions determine ingredients and properties of densematters. Strong magnetic fields, if they are present, may modify the dynamics.There is a long list of studies which explore these phenomena in the framework ofnumerical relativity (see [3] and references therein for resent simulations of stellar corecollapse, and see [4, 5] and references therein for those of compact binary merger). Inmost of the previous studies, however, treatments of microphysics are very simplifiedand more sophisticated studies are necessary.Furthermore, recent observations [6, 7, 8, 9, 10, 11, 12] have discovered thespectroscopic connections between several supernovae and long gamma-ray bursts(GRBs), clarifying that at least some of long GRBs are associated with the collapseof massive stars. Also, there are theoretical models that a short GRB occurs as a resultof a binary neutron star merger [13, 14]. The relevant process of the energy depositionto form a GRB fireball may be pair annihilation of neutrinos emitted from a hot massivedisk around a black hole formed after the collapse or the merger. These also enhancethe importance of exploring stellar core collapse and coalescence of compact star binaryin full general relativity taking account of microphysical processes.Gravitational wave astronomy will start in this decade. The first generation ofground-based interferometric detectors (LIGO [15], VIRGO [16], GEO600 [17]) arenow in the scientific search for gravitational waves. To obtain physically valuableinformation from these observations, it is necessary to connect the observed data andthe physics behind it. For this purpose, performing numerical simulation is the uniqueapproach. However, accurate predictions of gravitational waveforms are still hamperedby the facts that reliable estimates of waveforms require a general relativistic treatment[18, 19], and that appropriate treatments of microphysics such as a nuclear equationof state (EOS), the electron capture, and neutrino emissions and transfers. Generalrelativistic simulations including microphysics are required to make accurate predictionsof gravitational waveforms.As described above, to perform multidimensional simulations in the frame work ofnumerical relativity implementing microphysics is currently one of the most important eneral relativistic leakage scheme
The leakage schemes [28, 29, 30, 31, 32] as an approximate method for the neutrinocooling has a well-established history (e.g. [31]). The basic concept of the originalneutrino leakage schemes [28, 29] is to treat the following two regions in the systemseparately: one region is where the diffusion timescale of neutrinos is longer than thedynamical timescale, and hence neutrinos are ’trapped’ (neutrino-trapped region); theother region is where the diffusion timescale is shorter than the dynamical timescale, andhence neutrinos stream out freely out of the system (free-streaming region). The ideaof treating the diffusion region separately has been applied to more advanced methodsfor the neutrino transfer (see e.g., [33] and references therein).Then, electron neutrinos and anti-neutrinos in the neutrino-trapped region areassumed to be in β -equilibrium state. The net local rates of lepton-number and energyexchange with matters are set to be zero in the neutrino-trapped region. To treatdiffusive ’leakage’ of neutrinos out of the neutrino-trapped region, phenomenologicalsource terms based on the diffusion theory are introduced [28, 29]. In the free-streamingregion, on the other hand, it is assumed that neutrinos escape from the system withoutinteracting with matter. Therefore, neutrinos carry the lepton number and the energyaccording to the local weak interaction rates. The neutrino fractions are not solved inthe original version of the leakage scheme: Only the total lepton fraction is necessary inthe free-streaming region and the neutrino fractions are set to zero in the free-streamingregion. Note that there is a sharp discontinuity between the two regions. Consequently,thermodynamical quantities, in particular those of neutrinos and the electron fraction, eneral relativistic leakage scheme transfer of neutrinos are not solved in the leakage schemes. Therefore, theycannot treat non-local interactions among the neutrinos and matters; for example, theso-called neutrino heating [34] and the neutrino pair annihilation cannot be treated in theleakage scheme. Nevertheless, I consider a detailed general relativistic leakage schemepresented in this paper to be an important step towards more reliable and sophisticatedmodels, since the simulated physical timescales in the case of compact binary mergerswill be order of 10 ms and neutrino transfer is expected to be unimportant [35], andsince the neutrino heating would be not very important in the case of prompt black holeformation.Usually, the boundary between the neutrino-trapped and free-streaming regions isgiven by hand as a single ’neutrino-trapping’ density ( ρ trap ) in the previous simulationsof stellar core collapse [28, 29, 32]. In fact, however, the location at which the neutrinotrapping occurs depends strongly on the neutrino energies ( ǫ ν ), and hence, there aredifferent neutrino-trapping densities for different neutrino energies. The neutrino-trapping densities depend strongly on the neutrino energies as ρ trap ∝ ǫ − ν [36]. Thisimplies that neutrinos with lowest energy leave their corresponding neutrino-trappingregion first, and neutrinos with higher energy are emitted later.In the previous leakage schemes [28, 29, 32], on the other hand, all neutrinos areemitted in one moment irrespective of their energy. Consequently in the case of theso-called neutrino burst emission (e.g., [36]), for example, the duration in which theneutrinos are emitted is shortened and the peak luminosity at the burst is overestimatedin the previous leakage schemes [29, 27]. The dependence of the neutrino-trappingdensities and neutrino diffusion rates on the neutrino energies are approximately takeninto account in the recent simulations of binary neutron star mergers [37, 35]. Howeverthe lepton-number conservation equations for neutrinos are not solved [37].Recently, a numerical code based on a relativistic extension of the leakage schemeswas developed in [27], where not the region of the system but the energy momentumtensor of neutrinos are decomposed into two parts; ’trapped-neutrino’ and ’streaming-neutrino’ parts. However the source terms of hydrodynamic and the lepton-number-conservation equations are determined using the single neutrino-trapping density asin the case of the previous leakage schemes. More recently, Liebend¨orfer et al. [38]proposed a scheme, which they call the isotropic diffusion source approximation, wherethe neutrino distribution function is decomposed into an isotropic distribution functionof trapped neutrinos and a distribution function of streaming neutrinos.The present work is based on the previous studies described above. The frameworkof general relativistic extension of leakage scheme is based on my previous study in[27]. The treatment of neutrino diffusion rates is based on the recent work by Rosswogand Liebend¨orfer [35] where the neutrino-energy dependences are taken into account.Thus the remaining main problem to implement the relevant microphysics is thatstraightforward explicit scheme cannot be adopted to solve the equations [40] sincethe characteristic timescale of weak interaction processes ( t wp ∼ | Y e / ˙ Y e | , hereafter the eneral relativistic leakage scheme t dyn ) in hot dense regions,as described in Sec. 2. Note that the WP timescale is different from the so-called weaktimescale. In this paper I present a simple and stable method in which the equations aresolved explicitly in the dynamical timescale in the fully general relativistic framework.The paper is organized as follows. First, main difficulties of implementation of weakinteractions and neutrino cooling in full general relativity compared to implementationof them in Newtonian framework are briefly summarized in Sec. 2. Then, frameworkof the implementation of the microphysics is described in detail in Sec. 3. Some detailsof microphysics and numerics are described in Sec. 4, although GR leakage frameworkis independent of specific implementations of microphysics. In Sec. 5, results of a testsimulation is briefly described illustrating good ability of the present implementation.Section 6 is devoted to the summary and discussions. Throughout the paper, thegeometrical unit c = G = 1 is used otherwise stated.
2. Difficulties of implementation of weak interactions and neutrino coolingin full general relativity
Since the characteristic timescale of weak interaction precesses (the WP timescale t wp ∼ | Y e / ˙ Y e | ) is much shorter than the dynamical timescale t dyn in hot dense matters[40, 35], the numerical treatment of the weak interactions cannot be explicit, as notedin the previous pioneering work by Bruenn [40]. Otherwise a very short timestep (∆ t 3. General relativistic neutrino leakage scheme In the following, I describe in some detail a method for solving all of the equation in anexplicit manner. As described in the previous section, since t wp ≪ t dyn in the hot densematter regions, the source terms in the equations become too stiff to be solved explicitly.The characteristic timescale of leakage of neutrinos from the system t leak , however, ismuch longer than t wp in the hot dense matter region. Note that t leak ∼ L/c ∼ t dyn where L is the characteristic length scale of the system. On the other hand, t leak is comparableto t wp in the free-streaming regions where the WP timescale is longer than or comparablewith the dynamical timescale. Utilizing these facts, I approximate the original matterequations and reformulate them so that the source terms are to be characterized by theleakage timescale t leak . eneral relativistic leakage scheme Now, the problem is that the source term Q α in Eqs. (14) and (15) becomes too stiff tosolve explicitly in hot dense matter regions where t wp ≪ t dyn . To overcome the situation,the following procedures are adopted.First, it is assumed that the energy-momentum tensor of neutrinos are bedecomposed into ’trapped-neutrino’ (( T ν, T ) αβ ) and ’streaming-neutrino’ (( T ν, S ) αβ ) partsas [27], ( T ν ) αβ = ( T ν, T ) αβ + ( T ν, S ) αβ . (16)Here, the trapped-neutrinos phenomenologically represent neutrinos which interactsufficiently frequently with matter and are thermalized, while the streaming-neutrinopart describes a phenomenological flow of neutrinos streaming out of the system [27](see also [38] where a more sophisticate method based on the distribution function isadopted 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 α . (17)Then, the equation of the trapped-neutrino part becomes ∇ β ( T ν, T ) βα = Q α − Q leak α . (18)Third, the trapped-neutrino part is combined with the fluid part to give T αβ ≡ ( T F ) αβ + ( T ν, T ) αβ , (19)and Eqs. (14) and (18) are combined to give ∇ β T βα = − Q leak α . (20)Thus the equations to be solved is changed from Eqs. (14) and (15) to Eqs. (20)and (17). Note that the new equations only include the source terms Q leak α which ischaracterized by the leakage timescale t leak . Definition of Q leak α will be found in Sec. 3.3.The energy-momentum tensor of the fluid and trapped-neutrino parts ( T αβ ) istreated as that of the perfect fluid T αβ = ( ρ + ρε + P ) u α u β + P g αβ , (21)where ρ and u α are the rest mass density and the 4-velocity. The specific internal energydensity ( ε ) and the pressure ( P ) are the sum of the contributions from the baryons (freeprotons, free neutrons, α -particles, and heavy nuclei), leptons (electrons, positrons, and trapped-neutrinos ), and the radiation as, P = P B + P e + P ν + P r , (22) ε = ε B + ε e + ε ν + ε r , (23)where subscripts ’ B ’, ’ e ’, ’ r ’, and ’ ν ’ denote the components of the baryons, electronsand positrons, radiation, and trapped-neutrinos, respectively. eneral relativistic leakage scheme T ν, S ) αβ = En α n β + F α n β + F β n α + P αβ , (24)where F α n α = P αβ n α = 0. In order to close the system, we need an explicit expressionof P αβ . In this paper, I adopt a rather simple form P αβ = χEγ αβ with χ = 1 / 3. Thisapproximation may work well in high density regions but will violate in low densityregions. However, the violation will not affect the dynamics since the total amount ofstreaming-neutrinos emitted in low density regions will be small. Of course, a moresophisticated treatment will be necessary in a future study. The conservation equations of the lepton fractions can be written schematically as dY e dt = − γ e , (25) dY νe dt = γ νe , (26) dY ¯ νe dt = γ ¯ νe , (27) dY νx dt = γ νx , (28)where Y e , Y νe , Y ¯ νe , and Y νx denote the electron fraction, the electron neutrino fraction,the electron anti-neutrino fraction, and µ and τ neutrino and anti-neutrino fractions,respectively. It should be addressed that, in the previous simulations based on theleakage schemes [28, 29, 32, 37], the neutrino fractions are not solved.The source terms of neutrino fractions can be written, on the basis of the presentleakage scheme, as γ νe = γ local νe − γ leak νe , (29) γ ¯ νe = γ local¯ νe − γ leak¯ νe , (30) γ νx = γ local νx − γ leak νx , (31)where γ local ν and γ leak ν are the local production and the leakage rates of neutrinos,respectively (see Sec. 3.3). Note that only the trapped-neutrinos are responsible forthe neutrino fractions. The thermodynamical quantities (e.g., the pressure and thechemical potentials) of neutrinos can be calculated from the neutrino fractions on theassumption of thermalization of the trapped neutrinos.The source term for the electron fraction conservation can be written γ e = γ local νe − γ local¯ νe . (32)Since γ local ν s are characterized by by the WP timescale t wp , some procedures arenecessary to solve the lepton conservation equations explicitly. The following simpleprocedures are proposed to solve the equation stably. eneral relativistic leakage scheme n , the conservation equation of the total lepton fraction( Y l = Y e − Y νe + Y ¯ νe ), dY l dt = − γ l , (33)is solved together with the conservation equation of Y νx , Eq. (28), in advance of solvingwhole of the lepton conservation equations (Eqs. (25) – (28)). Note that the sourceterm γ l = γ leak νe − γ leak¯ νe is characterized by the leakage timescale t leak so that this equationmay 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 evolved Y l .Second, regarding 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 in thetimestep n + 1 do not exceed Y βν ’s. Then, the whole of the lepton conservation equations(Eqs. (25) – (28)) are solved explicitly utilizing the limiters.Third, the following conditions are checked µ p + µ e < µ n + µ νe , (34) µ n − µ e < µ p + µ ¯ νe , (35)where µ p , µ n , µ e , µ νe and µ ¯ νe are the chemical potentials of protons, neutrons, electrons,electron neutrinos, and electron anti-neutrinos, respectively. If both conditions aresatisfied, the values of the lepton fractions in the timestep n + 1 is set to be those in the β -equilibrium value; Y βe , Y βνe , and Y β ¯ νe . On the other hand, if either or both conditionsare not satisfied, the lepton fractions in the timestep n + 1 is set to be those obtainedby solving whole of the lepton-number conservation equations.A limiter for the evolution of Y νx may be also necessary in some case where the pairprocesses are dominant, for example, in simulations of collapse of population III stellarcore. In this case, the value of Y νx at the pair equilibrium (i.e. at µ νx = 0), Y pair νx maybe used to limit the source term.In the present implementation it is not necessary to somewhat artificially dividethe system into neutrino-trapped and free-streaming regions. Therefore there is nodiscontinuous boundary which existed in the previous leakage schemes [28, 29, 32].I found that simulations of the collapse of population III stellar core and theformation of a black hole, in which very high temperatures ( T > 100 MeV) are achieved,can be stably performed using the simple procedure presented in this paper. In this subsection the definitions of the leakage rates Q leak α and γ leak ν are presented.Because Q leak ν may be regarded as the emissivity of neutrinos measured in the fluid restframe , Q leak α is defined as [50] Q leak α = Q leak ν u α . (36) eneral relativistic leakage scheme H α which satisfies H α u α = 0, Eq. (3.3) may be the best choice in the present framework.The leakage rates Q leak ν and γ leak ν are assumed to satisfy the following properties.(i) The leakage rates approach the local rates Q local ν and γ local ν in the low density,transparent region.(ii) The leakage rates approach the diffusion rates Q diff ν and γ diff ν in the high density,opaque region.(iii) The above two limits should be connected smoothly.Here, the local rates can be calculated based on the theory of weak interactions (see Sec.4.3 for the local rates adopted in this paper) and the diffusion rates can be determinedbased on the diffusion theory (see Sec. 4.4 for the definition of the diffusion rate adoptedin this paper). There will be several prescriptions to satisfy the requirement (iii) [37, 35].In this paper, the leakage rates are defined as Q leak ν = (1 − e − bτ ν ) Q diff ν + e − bτ ν Q local ν , (37) γ leak ν = (1 − e − bτ ν ) γ diff ν + e − bτ ν γ local ν , (38)where τ ν is the optical depth of neutrinos and b is a parameter which is typically set as b − = 2 / 3. The optical depth can be computed from the cross sections in a standardmanner [37, 35]. The basic equations for the general relativistic hydrodynamics are the continuityequation, the lepton-number conservation equations, and the local conservation equationof the energy-momentum. The explicit forms of the equations are presented in thissubsection for the purpose of convenience. The continuityequation is ∇ α ( ρu α ) = 0 . (39)As fundamental variables for numerical simulations, the following quantities areintroduced: ρ ∗ ≡ ρwe φ and v i ≡ u i u t where w ≡ αu t . Then, the continuity equation iswritten as ∂ t ( ρ ∗ √ η ) + ∂ k ( ρ ∗ v k √ η ) = 0 , (40)where √ η ≡ q det η ij is the volume element of the flat space in the curvilinearcoordinates.The lepton-number conservation equations (25) – (28) can be abbreviated as dY L dt = γ L , (41)Using the continuity equation, they become ∂ t ( ρ ∗ Y L √ η ) + ∂ k ( ρ ∗ Y L v k √ η ) = ρ ∗ γ L . (42) eneral relativistic leakage scheme As discussed in Sec. 3.1, I solve the followingequations. ∇ β T βα = − Q leak α , (43) ∇ β ( T ν, S ) βα = Q leak α , (44)where T αβ and ( T ν, S ) αβ are given by Eqs. (21) and (24), respectively. The source term Q leak α is defined by Eq. (37).As fundamental variables for numerical simulations, I define the quantities ˆ u i ≡ hu i and ˆ e ≡ hw − P ( ρw ) − . Then, the Euler equation ( γ αi ∇ β T βα = γ αi Q α ), and the energyequation ( n α ∇ β T αβ = n α Q α ) can be written as ∂ t ( ρ ∗ ˆ u A √ η ) + ∂ k hn ρ ∗ ˆ u A v k + P αe φ δ kA o √ η i = − ρ ∗ " wh∂ A α − ˆ u i ∂ A β i + αe − φ wh ˆ u k ˆ u l ∂ A ˜ γ kl − αh ( w − w ∂ A φ + P ∂ A ( αe φ ) + P αe φ δ ̟A ̟ + αe φ Q A , (45) ∂ t ( ρ ∗ ˆ u ϕ √ η ) + ∂ k (cid:16) ρ ∗ ˆ u ϕ v k √ η (cid:17) = αe φ Q ϕ , (46) ∂ t ( ρ ∗ ˆ e √ η ) + + ∂ k hn ρ ∗ v k ˆ e + P e φ √ η ( v k + β k ) o √ η i = αe φ √ ηP K + ρ ∗ u t h ˆ u k ˆ u l K kl − ρ ∗ ˆ u i γ ij D j α + αe φ Q α n α , (47)where the subscript A denotes ̟ or z component.The evolution equations of streaming-neutrinos ( E and F i ) are written as ∂ t ( √ γE ) + ∂ k h √ γ ( αF k − β k E ) i = √ γ (cid:16) αP kl K kl − F k ∂ k α + αQ leak a n a (cid:17) , (48) ∂ t ( √ γF i ) + ∂ k h √ γ ( αP ki − β k F i ) i = √ γ (cid:18) − E∂ i α + F k ∂ i β k + α P kl ∂ i γ kl + αQ leak i (cid:19) . (49) In each numerical timestep, the so-called primitive variables ( ρ , Y L , T , and v i ) and theLorentz factor w = αu t = q γ ij u i u j must be calculated from the conserved quantities( ρ ∗ , ρ ∗ Y L , ˆ e , and ˆ u i ), where Y L is the representative of the lepton fractions. Since theequation of state (EOS) of the nuclear matter are usually tabularized in terms of theargument quantities ( ρ , Y p (= Y e ), T ), I am devoted to the cases of the tabularized EOSin the following.In the case where the whole of the lepton-number conservation equations are solved(see Sec. 3.2), the argument quantities ( ρ , Y e , T ) are calculated from the conservedquantities as follows.(i) Give a trial value, ˜ w , of the Lorentz factor. Then, one obtains a trial value, ˜ ρ , ofthe rest mass density: ˜ ρ = ρ ∗ / ( ˜ we φ ).(ii) A trial value, ˜ T , of the temperature can be obtained by solvingˆ e = ˆ e EOS ( ˜ ρ, Y e , ˜ T , Y νe , Y ¯ νe , Y νx ) , (50) eneral relativistic leakage scheme e EOS is constructed from EOS table. Note that ˆ e and ˆ e EOS in general containcontributions from trapped-neutrinos. One dimensional search over the EOS tableis required to obtain ˜ T .(iii) The next trial value of the Lorentz factor is given by solving ˜ w = q e − φ ˜ γ ij ˆ u i ˆ u j ˜ h − , where the specific enthalpy ˜ h is calculated from EOS tableas ˜ h = ˜ h ( ˜ ρ, Y e , ˜ T ).(iv) Repeat the procedures (i)–(iii) until a required degree of convergence is achieved.Convergent solutions of the temperature and w are obtained typically within 10iterations.On the other hand, in the case where the total lepton fraction is evolved, theargument quantities ( ρ, Y e , T ) must be recovered from the conserved quantities and Y l under the assumption of the β -equilibrium. In this case, two-dimensional reconstruction( Y l , ˆ e ) = ⇒ ( Y e , T ) (51)would be required for a given ˜ w . In this case, there may be in general more than onecouple of ( Y e , T ) which gives the same Y l and ˆ e . Therefore, I adopt a different methodto recover ( ρ, Y e , T ) [27].Under the assumption of the β -equilibrium, the electron fraction is related to thetotal lepton fraction as Y e = Y e ( ρ, Y l , T ). Using this relation, the original EOS tablecan be reconstructed in terms of the argument quantities of ( ρ , Y l , T ). Then, the samestrategy as in the above can be adopted. Namely,(i) Give a trial value, ˜ w of w . Then, one obtains a trial value, ˜ ρ , of the rest massdensity.(ii) A trial value, ˜ T , of the temperature can be obtained by solvingˆ e = ˆ e EOS ( ˜ ρ, Y l , ˜ T , Y νx ) . (52)One dimensional search over the EOS table is required to obtain ˜ T .(iii) The next trial value of w is given by solving ˜ w = q e − φ ˜ γ ij ˆ u i ˆ u j ˜ h − .(iv) Repeat the procedures (i)–(iii) until a required degree of convergence is achieved.The electron and electron neutrino fractions are given as Y e = Y e ( ρ, Y l , T ), Y νe = Y νe ( ρ, Y l , T ), Y ¯ νe = Y ¯ νe ( ρ, Y l , T ) in the new EOS table.The construction of EOS table in terms of the argument variables of ( ρ , Y l , T ) isimportant in the present implementation.In the case of a simplified or analytic EOS, the Newton-Raphson method may beapplied to recover the primitive variables. In the case of tabulated EOS, by contrast,the Newton-Raphson method may not good approach since it requires derivativesof thermodynamical quantities which in general cannot be calculated precisely fromtabulated EOS by the finite differentiating method (see also Sec. 4.2). eneral relativistic leakage scheme 4. Specific details of microphysics . While any EOS table can be used in the present code, an EOS [51, 52] basedon the relativistic mean field theory is adopted for baryon EOS (hereafter denoted byShen EOS) in the present version of our code. Note that the causality is guaranteed tobe satisfied in the relativistic EOS, whereas the sound velocity sometimes exceeds thevelocity of the light in the non-relativistic framework, e.g., in the EOS by Lattimer andSwesty [53]. Electrons and Positrons . If a EOS table for baryons does not include thecontributions of the leptons (electrons, positrons, and neutrinos if necessary) andphotons, one has to consistently include these contributions to the table. Electronsand positrons are described as ideal Fermi gases.To consistently calculate the contribution of the electrons, the charge neutralitycondition Y p = Y e should be solved in terms of the electron chemical potential µ e , foreach value of the baryon rest-mass density ρ and the temperature T in the EOS table: n e ( µ e , T ) ≡ n − − n + = ρY p m u (53)in terms of µ e for given ρ, T , and Y p . Here, m u = 931 . n − and n + are the total number densities (i.e., including electron-positron pairs) ofelectrons and positrons, respectively. Then all other quantities can be calculated from T and µ e . Radiations . The contribution of radiations is included in a standard manner: Theradiation pressure and the specific internal energy density are given by ε r = a r T ρ , P r = a r T , (54)where a r is the radiation constant a r = (8 π k B )(15 c h P ) − and k B and h P are theBoltzmann’s and the Plank’s constants respectively. Trapped-Neutrinos . In this paper, the trapped-neutrinos are assumed to interactsufficiently frequently with matter that be thermalized. Therefore they are described asideal Fermi gases with the matter temperature. Then, from the neutrino fractions Y ν ,the chemical potentials of neutrinos are calculated by solving Y ν = Y ν ( µ ν , T ) . (55)Using the chemical potentials and the matter temperature, the pressure and the internalenergy of the trapped-neutrinos are calculated. In high-resolution shock-capturing schemes, it is in general necessary to evaluate thesound velocity c s , c s = 1 h ∂P∂ρ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ǫ + Pρ ∂P∂ǫ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ . (56) eneral relativistic leakage scheme ∂P∂ρ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ǫ = X i = B,e,r,ν ∂P i ∂ρ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T − ∂P i ∂T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ X j = B,e,r,ν ∂ǫ j ∂ρ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T X k = B,e,r,ν ∂ǫ k ∂T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ − , (57) ∂P∂ǫ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ = X i = B,e,r,ν ∂P i ∂T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ X j = B,e,r,ν ∂ǫ j ∂T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ − , (58)where ’ B ’, ’ e ’, ’ r ’, and ν in the sum denote contributions of the baryon, the electrons,radiations, and neutrinos, respectively.Since there are in general the phase transition regions in a EOS table for baryons andthe EOS moreover may contain some non-smooth spiky structures, careful treatmentsare necessary when evaluating the derivatives of thermodynamical quantities. In thepresent EOS table, the derivatives are carefully evaluated so that there are no spikybehaviors in the resulting sound velocities. In this paper, the electron and positron captures ( γ ec ν and γ pc ν ) [54], the electron-positronpair annihilation ( γ pair ν ) [55], the plasmon decays ( γ plas ν ) [37], and the Bremsstrahlungprocesses ( γ Brems ν ) [56] are considered as the local production reactions of neutrinos.Then, the local rates for lepton fractions are γ local e = γ ec ν − γ pc ν , (59) γ local νe = γ ec ν + γ pair ν + γ plas ν + γ Brems ν , (60) γ local¯ νe = γ pc ν + γ pair ν + γ plas ν + γ Brems ν , (61) γ local νx = γ pair ν + γ plas ν + γ Brems ν . (62)Similarly, the local energy emission rate Q local ν is the sum of the contributions of theelectron and positron captures ( Q ec ν and Q pc ν ) , the electron-positron pair annihilation( Q pair ν ), the plasmon decays ( Q plas ν ), and Bremsstrahlung processes ( Q Brems ν ). I follow [35] for the neutrino diffusion rates γ diff ν and Q diff ν . I present the forms ofthe diffusion rates in the following for convenience. An alternative definition of thediffusion rates will be found in [37]. The cross sections adopted in this paper are thoseof neutrino-nucleus and neutrino-nucleon scattering, and neutrino absorptions on freenucleons. Explicit forms of these cross sections will be found in [56]Ignoring the higher order corrections, the neutrino cross sections can be written ingeneral as σ ( E ν ) = E ν ˜ σ, (63)where E ν is the neutrino energy and ˜ σ is ’cross section’ in which E ν dependence isfactored out. Similarly, the opacity and the optical depth are written as κ ( E ν ) = X κ i ( E ν ) = E ν X ˜ κ i = E ν ˜ κ, (64) eneral relativistic leakage scheme τ ( E ν ) = Z κ ( E ν ) ds = E ν Z ˜ κds = E ν ˜ τ . (65)Now the neutrino diffusion time may be defined by [37, 35] T diff ν ( E ν ) ≡ ∆ x ( E ν ) c τ ( E ν ) = E ν a diff ˜ τ c ˜ κ = E ν ˜ T diff ν , (66)where the distance parameter ∆ x ( E ν ) is set to be∆ x ( E ν ) = a diff τ ( E ν ) κ ( E ν ) . (67)Here a diff is a parameter which controls the diffusion rates. In this paper, I adopt a diff = 3 as suggested in [37].Finally, the neutrino diffusion rates are defined as N diff ν ≡ Z n ν ( E ν ) T diff ν ( E ν ) dE ν = 1 a diff πcg ν ( h P c ) ˜ κ ˜ τ T F ( η ν ) , (68) Q diff ν ≡ Z E ν n ν ( E ν ) T diff ν ( E ν ) dE ν = 1 a diff πcg ν ( h P c ) ˜ κ ˜ τ T F ( η ν ) , (69)from which the diffusion rates Q diff ν and γ diff ν are easily calculated. Here, g ν is thestatistical weight factor for neutrinos and n ( E ν ) dE ν is the number density of neutrinoin the range from E ν to E ν + dE ν under the Fermi-Dirac distribution. 5. Validity of general relativistic leakage scheme The numerical schemes for solving the Einstein’s equations are essentially same as thosein [4]; We adopt so-called BSSN formulation [58, 59] and use 4th-order finite differencescheme in the spatial direction and the 3rd-order Runge-Kutta scheme in the timeintegration. The advection terms such as β i ∂ i φ are evaluated by a 4th-order upwindscheme. The hydrodynamic equations, the lepton-number conservation equations, andModel Φ c ≤ . ≤ Φ c ≤ . ≤ Φ c ≤ . ≤ Φ c ≤ . c ≥ . S15 ∆ x η N 444 668 924 1212 1532 L (km) x low η resolution N 316 444 636 828 1020 L (km) Table 1. Summary of the regridding procedure. The values of the minimum gridspacing ∆ x (in units of km), the non-uniform-grid factor η , and the grid number N for each range of Φ c = 1 − α c are listed. eneral relativistic leakage scheme 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 )] Figure 1. The radial profiles of the infall velocity, the density, the entropy per baryonand the total lepton fraction at bounce, 2 ms and 6 ms after bounce. The results forthe finer grid resolution (solid curve) and for the coarser grid resolution (the dottedcurves) are shown together while they are almost identical. equations of streaming-neutrinos are solved using the high-resolution centered scheme[62]. A nonuniform grid is adopted in the numerical simulation, in which the grid spacingincreases as dx j +1 = ηdx j , dz l +1 = ηdz l (70)where dx j ≡ x j +1 − x j , dz l ≡ z l +1 − z l and η is a constant. The regridding procedure[60, 61] is furthermore used to compute the collapse accurately and to save theCPU time efficiently. For the regridding, I define an effective gravitational potentialΦ c ≡ − α c (Φ c > 0) where α c is the central value of the lapse function. In Table 1,parameters of the regridding procedure are summarized. More detailed set up of thesimulation will be found elsewhere [63].As a test problem, I performed a collapse simulation of spherical presupernovacore. A presupernova model (S15) of 15 M ⊙ with solar metallicity computed in [57] isadopted as the initial condition. I follow the dynamical evolution of central part whichis composed of the Fe core and some part of the Si-shell. The density, the electronfraction, and the temperature are used to calculate other thermodynamical quantitiesusing the EOS table.To check the validity of the code, the results are compared with those in the state-of- eneral relativistic leakage scheme Figure 2. Snapshots of the contours of the density (top left panels), the electronfraction Y e (top right panels), the entropy per baryon (bottom left panels), and thelocal neutrino energy emission rate (bottom right panels) in the x - z plane at selectedtime slices. the-art one-dimensional simulations (hereafter, the reference simulations) in full generalrelativity [64, 21, 65, 22], where one dimensional general relativistic Boltzmann equationis solved for neutrino transfer with relevant weak interaction processes. Since neutrinoheating processes ( ν e + n → p + e − and ¯ ν e + p → n + e + ) are not include in the presentimplementation, and multidimensional effects such as convection cannot be followedin the one-dimensional reference simulations, I pay particular attention in comparingresults during the collapse and the early phase ( ∼ 10 ms) after the bounce. After that,direct comparison cannot be done since in the present multidimensional code convectiveactivities set in. As shown below, results in the present simulation and in the referencesimulations agree well. The collapse proceeds until the nuclear density is reached in the central part of theiron core. Then, the inner core experiences the bounce due to the nuclear repulsiveforces, forming a strong shock wave at the edge of the inner core. The shock wavepropagates outward and when it crosses the neutrino-sphere, spiky burst emissions ofneutrinos occur (neutrino bursts): Neutrinos in hot post-shock region are copiously eneral relativistic leakage scheme -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 Figure 3. Time evolution of the neutrino luminosities. The results in the finer gridresolution (solid curves) and in the coarser grid resolution (dashed curves) are showntogether. The two results are almost identical until the convective phase sets in, whilethey are not in the convective phase. emitted without interacting matters. Eventually, negative gradients of the total leptonfraction are formed behind the shock since neutrinos carry away the lepton number. InFig. 1, we show the radial profiles of the infall velocity, the density, the entropy perbaryon and the total lepton fraction at selected time slices.The results agree at least semi-quantitatively with those in [64, 21, 65, 22]. Inparticular, the radial profiles of the infall velocity, the density, and the entropy perbaryon show good agreements. No such good agreements was reported in the previousNewtonian simulations in which leakage schemes are adopted [28, 29, 31, 32, 30]. Thenegative gradients quantitatively are little bit steeper in the present simulation. Thereason may be partly because the transfers of lepton-number and energy are not solvedin the present leakage scheme. Except for this quantitative difference, the two resultsagree well. It is found that the difference can be reduced by adjusting the parameter a diff introduced in Sec. 4.4.Recall that regions of negative Y l gradient are known to be convectively unstable[36]. Convective activities indeed set in in the present simulation as shown in Fig. 2. Comparisons of the neutrino luminosities are particularly important since they dependon both implementations of weak interactions (especially electron capture in the present eneral relativistic leakage scheme L ν = Z αe φ u t ˙ Q leak ν d x, (71)as functions of t − t bounce where t bounce is time at the bounce. The result also agreesapproximately with that in the reference simulations. The neutrino bursts occur whenthe shock wave crosses the neutrino-sphere soon after the bounce. The peak luminosityat the neutrino burst is L ν e ≈ . × ergs/s in the present simulation, which agreeswell with that in the reference simulations. The peak luminosity and the duration(width) of the neutrino burst emission can be improved by adjusting the parameter a diff . The modulation found in the later phase t − t bounce > 10 ms is due to convectiveactivities driven by negative gradients of electron fraction and entropy per baryon.Thus, Fig. 3 illustrates that the present detailed leakage scheme works fairly welland may be applied to simulations of rotating core collapse to a black hole and mergersof binary neutron stars.In the previous simulations based on the leakage scheme [32, 27] where the single’neutrino-trapping’ density is adopted, the luminosities do not agree with that in thereference simulations. In particular, the luminosities at the neutrino bursts are quitedifferent. In Figs. 1 and 3, I show results in the higher resolution (solid curves) and the lowerresolution (dashed curves). The radial profiles of the two resolutions are almost identical,showing that convergent results are obtained in the present simulation (see Fig. 1). Inthe time evolution of neutrino luminosities (see Fig. 3), the two results are almostidentical before the convective activities set in. In the later phase, on the other hand,the two results shows slight difference. Since the convection and the turbulence canoccurs in a infinitesimal scale length, the smaller-scale convection and turbulence arecaptured in the finer grid resolution. Further discussions associated with the convergenceand numerical accuracy will be found in [63]. 6. Summary and Discussions In this paper, I presented an implementation of the weak interactions and the neutrinocooling in the framework of full general relativity. Since the characteristic timescale ofweak interaction processes t wp ∼ | Y e / ˙ Y e | is much shorter than the dynamical timescale t dyn in hot dense matters, stiff source terms appears in the equations. In general, an eneral relativistic leakage scheme t leak is much longer than t wp and iscomparable to t dyn .By decomposing the energy tensor of neutrino into the trapped-neutrino andthe streaming-neutrino parts, the equations for the energy momentum tensor can berewritten so that the source terms are characterized by the leakage timescale t leak (seeEqs. (20) and (17)). The lepton-number conservation equations, on the other hand,include the source terms characterized by the WP timescale. Therefore the limiters forthe stiff source terms are introduced to solve the lepton-number conservation equationsexplicitly (see Sec. 3.2).In the numerical relativistic hydrodynamics, it is required to calculate the primitivevariables and the Lorentz factor from the conserved quantities. In this paper, I developa robust and stable procedure for it (Sec. 3.5).Finally, to check the validity of the present implementation, I performed acollapse simulation of spherical presupernova core and compared the results with thoseobtained in the state-of-the-art one-dimensional simulations in full general relativity[64, 21, 65, 22]. As shown in this paper, results in this paper agree well with thosein the state-of-the-art simulations. Thus the present implementation will be applied tosimulations of rotating core collapse to a black hole and mergers of binary neutron stars. Since the present implementation of the microphysics is simple and explicit, it hasadvantage that the individual microphysical processes can be easily improved andsophisticated.For example, the neutrino emission via the electron capture can be easilysophisticated as follows. To precisely calculate the electron capture rate, the completeinformation of the parent and daughter nuclei are required. In the nuclear equations ofstate currently available, however, a representative single-nucleus average for the trueensemble of heavy nuclei is adopted. The representative is usually the most abundantnuclei. The problem in evaluating the capture rate is that the nuclei which cause thelargest changes in Y e are neither the most abundant nuclei nor the nuclei with the largestrates, but the combination of the two. In fact, the most abundant nuclei tend to havesmall rates since they are more stable than others, and the fraction of the most reactivenuclei tend to be small [66, 67]. Assuming that the nuclear statistical equilibrium (NSE)is achieved, the electron capture rates under the NSE ensemble of heavy nuclei may becalculated for given ( ρ , Y e , T ). Such a numerical rate table can be easily employed in eneral relativistic leakage scheme τ ν ∼ σT ∼ 1) changes as T ∼ σ − / . Although thecorrection terms are in general very complicated, it is straightforward to include thecorrections in the present implementation. Note that the corrections become moreimportant for higher neutrino energies. Therefore, the correction terms might play rolesin the collapse of population III stellar core and the formation of a black hole in whichvery high temperatures ( T > 100 MeV) are achieved. I already started studies to explorethe importance of these corrections in the case of black hole formation.As briefly described in the introduction, one of main drawbacks of the presentimplementation of the neutrino cooling is that the transfer of neutrinos are not solved.Although to fully solve the transfer equations of neutrinos is far beyond the scope ofthis paper, there are a lot of rooms for improvements in the treatment of the neutrinocooling. For example, the relativistic moment formalism [69, 70], in particular the so-called M1 closure formalism, may be adopted. For this purpose, a more sophisticatedtreatment of the closure relation for P αβ is required. For example, one may adopt theEddington tensor of the form [71] P αβ = " − χ γ αβ + 3 χ − F α F β F γ F γ E, (72)where χ is the (variable) Eddington factor. In the diffusion limit where the neutrinopressure is isotropic χ = 1 / 3, while in the free streaming limit χ = 1. We plan toimplement a relativistic M1 closure formalism for the neutrino transfer in the nearfuture.To conclude, the present implementation of microphysics in full general relativityworks fairly well. We are now in the standpoint where simulations of stellar core collapseto a black hole and merger of compact stellar binaries can be performed includingmicrophysical processes. Fruitful scientific results will be reported in the near future. acknowledgments I thank M. Shibata and L. Rezzolla for valuable discussions, and the referees for valuablecomments. I also thank T. Shiromizu and T. Fukushige for their grateful aids. Numericalcomputations were performed on the NEC SX-9 at the data analysis center of NAOJand on the NEC SX-8 at YITP in Kyoto University. This work is partly supported bythe Grant-in-Aid of the Japanese Ministry of Education, Science, Culture, and Sport(21018008,21105511). eneral relativistic leakage scheme References [1] van Riper K A 1988 Astrophys. J. Prog. Theor. Phys. Class. Quant. Grav. Phys. Rev. D Class. Quant. Grav. Nature Nature Nature Astrophys. J. L17[10] Kawabata K S et al. 2003 Astrophys. J. L19[11] Modjaz M et al. 2006 Astrophys. J. L21[12] Sollerman J et al. 2006 Astron. Astrophys. Phys. Rep. New J. Phys. Sience Class. Quant. Grav. S385[17] Willke B et al. 2002 Class. Quant. Grav. Astron. Astrophys. Phys. Rev. D Astron. Astrophys. Astrophys. J. Suppl. Astrophys. J. Astrophys. J. Phys. Rev. Lett. Phys. Rev. Lett. Astrophys. J. PhD thesis (University of Tokyo)[28] Epstein R I and Pethick C J 1981 Astrophys. J. Astrophys. J. Nucl. Phys. A Phys. Rep. Astrophys. J. Astrophys. J. Astrophys. J. Mon. Not. R. Astron. Soc. Rev. Mod. Phys. Astron. Astrophys. Astrophys. J. Astrophys. J. Suppl. Astrophys. J. Suppl. Astrophys. J. Astrophys. J. Astron. Astrophys. Astrophys. J. Astron. Astrophys. Astrophys. J. Astrophys. J. J. Comp. Appl. Math. Phys. Rev. B571 eneral relativistic leakage scheme [50] Shibata M Sekiguchi Y and Takahashi R 2007 Prog. Theor. Phys. Nucl. Phys. A Prog. Theor. Phys. Nucl. Phys. A Astrophys. J. Astrophys. J. Nucl. Phys. A Rev. Mod. Phys. Phys. Rev. D Phys. Rev. D Astrophys. J. L39[61] Sekiguchi Y and Shibata M 2005 Phys. Rev. D J. Comput. Phys. Prog. Theor. Phys. submitted[64] Liebend¨orfer M et al. 2001 Phys. Rev. D Astrophys. J. Astrophys. J. Suppl. Phys. Rept. Phys. Rev. D Astrophys. J. Mon. Not. R. Astron. Soc. Astrophys. J.248