A Neutrino free emission rates and opacity sources

# Implementation of a simplified approach to radiative transfer in general relativity

## Abstract

We describe in detail the implementation of a simplified approach to radiative transfer in general relativity by means of the well-known neutrino leakage scheme (NLS). In particular, we carry out an extensive investigation of the properties and limitations of the NLS for isolated relativistic stars to a level of detail that has not been discussed before in a general-relativistic context. Although the numerous tests considered here are rather idealized, they provide a well-controlled environment in which to understand the relationship between the matter dynamics and the neutrino emission, which is important in order to model the neutrino signals from more complicated scenarios, such as binary neutron-star mergers. When considering nonrotating hot neutron stars we confirm earlier results of one-dimensional simulations, but also present novel results about the equilibrium properties and on how the cooling affects the stability of these configurations. In our idealized but controlled setup, we can then show that deviations from the thermal and weak-interaction equilibrium affect the stability of these models to radial perturbations, leading models that are stable in the absence of radiative losses, to a gravitational collapse to a black hole when neutrinos are instead radiated.

gravitation – relativistic processes – methods: numerical

## I Introduction

Neutron stars (NS) are a unique astrophysical laboratory to study the microphysics of superdense matter. From their births in a supernova explosion until their deaths in a merger of compact objects or in a gravitational collapse to a black hole (BH), NSs undergo the most extreme physical conditions with respect to density, temperature and composition. For example, the temperature can vary from a surface temperature of a few , during the inspiral of two old NSs in a binary system (BNSs), up to at the time of merger Bauswein et al. (2012). Similarly, a NS can reach up to times nuclear saturation density  Baiotti et al. (2005a). The merger of two NSs is considered one of the most promising astrophysical events in which to test the models of nuclear matter in regimes that are unattainable in laboratory experiments. In addition to the bulk properties of the NSs, such as mass and radius, also the details of the NS equation of state (EOS), are expected to be imprinted in the electromagnetic emission (EM) and gravitational wave (GW) signal produced by these events.

In particular, the GW signal that accompanies the merger of BNSs is expected to be the prime target for the advanced generation of kilometer length ground-based detectors, such as LIGO Abbott et al. (2009), Virgo Acernese et al. (2009), and KAGRA Kuroda et al. (2010). Individual merger signals can reveal the NS radii with an accuracy of  km together with their masses Hinderer et al. (2010). These observations will strongly constrain the EOS that describes the NS matter at zero temperature. Moreover, GWs, together with their EM counterparts, will be of great relevance to determine the properties of the engine behind short gamma-ray bursts (GRBs) which are thought to arise from NS/NS and BH/NS mergers Eichler et al. (1989); Rezzolla et al. (2011). The incorporation of microphysical nuclear EOSs, general relativity, magnetic fields, neutrino radiation, viscosity and other transport coefficients in the numerical modeling of those events, is essential to increase the realism of the simulations and help us explore the possible outcomes from a large set of initial conditions and parameters.

A great deal of effort has been spent in the field of nuclear physics on the development of an EOS appropriate to describe the behavior of matter at both supranuclear and subnuclear densities (see Lattimer and Prakash (2001); Haensel et al. (2007) for reviews). Unfortunately, only few constraints on the properties of the EOS are available from relativistic heavy-ion collision experiments and astronomical observations (see Lattimer (2012) for a review on the topic). The recent discoveries of two massive NSs in a binary system with a mass of Demorest et al. (2010); Antoniadis et al. (2013) seem to rule out some of the NS matter models. Nevertheless, the lack of precise estimates of their radii did not allow to significantly reduce the parameter space for the NS EOSs. Combined measurements of mass and radius of NSs can only be obtained from x-ray bursts and from the NS surface thermal emission, something unlikely to happen in quiescent systems as the one described in Demorest et al. (2010) and Antoniadis et al. (2013). These measurements are subject to many uncertainties, as the distance of the NS and its magnetic field, that can spoil the accuracy of the estimates of the radius Özel et al. (2010). More importantly, the estimates of the physical parameters from the observations are sensitive to the models adopted for their analysis  Steiner et al. (2010).

As a result of the complexity involved in the theoretical models of NS matter, only few EOSs cover temperatures up to , densities that reach several times nuclear saturation density and a range of composition that goes from pure neutron matter to symmetric matter (equal neutron and proton number, ). Amongst the few EOSs that have gained popularity for astrophysical simulations, particularly well studied are the EOS of Lattimer and Swesty (1991) (hereafter, LS-EOS) and the EOS proposed by Shen et al. (1998a, b). More recently, a new EOS has been derived by Shen et al. (2010, 2011) (hereafter, SHT-EOS), in which the meson mean fields and the Dirac wave function describing the nucleons are solved self-consistently for each Wigner-Seitz cell. Despite the accuracy reached in the theoretical models, these EOSs not always satisfy the constraints imposed by both experiments (symmetry energy and incompressibility) and astronomical observations (mass and radius) Steiner et al. (2010); Hebeler et al. (2013).

Over the last few years, a lot of progress has been made to understand the dynamics of BNS systems in general relativity with respect to the possible detection of GWs from such systems (for a review see Faber and Rasio (2012)). The level of sophistication reached by present-day general-relativistic simulations includes important microphysical aspects Sekiguchi et al. (2011a, b); Bauswein et al. (2012) that in the past were only accessible to Newtonian simulations Ruffert et al. (1996, 1997); Ruffert and Janka (2001); Rosswog and Liebendörfer (2003). Together with nuclear EOSs to describe the NS matter, current approaches also include weak-interaction processes, albeit modeled by means of approximate neutrino treatments. Under the extreme thermodynamic conditions reached during the merger, powerful neutrino bursts are produced from the hot and shock heated NS matter Sekiguchi et al. (2011a). When the nuclear matter is strongly compressed and its temperature reach values of about , weak-interaction processes become increasingly important, leading the material out of the original chemical equilibrium with respect to processes while copious amounts of neutrinos are emitted. Neutrinos are responsible for the luminosities of the order attained in BNS mergers and they may play a role in powering the relativistic jet needed for the beamed emission of a GRB. The energy required for the outflow can be efficiently deposited along the baryon-free axis of rotation of the BH by neutrino pair annihilation Ruffert et al. (1996).

The aim of this paper is to present the details of the implementation in the new Whisky code (WhiskyThermal) Baiotti et al. (2003, 2005b) of a ”neutrino leakage” scheme (NLS) along with state-of-the-art nuclear EOSs with density, temperature and composition dependence. The NLS is a pragmatic approach to neutrino radiation transfer that estimates the local changes in the lepton number and the associated energy losses via neutrino emission. The relative simplicity of the implementation along with the small computational costs turn the NLS into an interesting approach to upgrade numerical relativity codes that evolve the equations of general-relativistic hydrodynamics (GRHD) and magneto-hydrodynamics together with the Einstein equations for a dynamical spacetime.

Neutrino interactions are influenced by density, temperature and composition of the hadronic matter and by the energy of the emitted neutrinos. At rest-mass densities of few times and temperatures around , the scattering of neutrinos onto baryons is so efficient that neutrinos quickly reach the thermal equilibrium with the NS matter and can be considered as trapped (i.e., with mean free path ). When the density drops below , neutrinos with energy below interact rarely with the nuclear matter and can therefore be considered as free-streaming (i.e., with mean free path ). The NLS scheme is particularly suited to study the evolution of a NS because of the sharp density gradient at the surface where the density drops by several orders of magnitude. More accurate schemes such as the direct Monte Carlo method for solving the Boltzmann transport equation Abdikamalov et al. (2012), are able to capture properly the dynamic in the semitransparent regime. Unfortunately, the computational cost involved limits the use of such schemes to one-dimensional problems. Three-dimensional simulations are performed instead using approximate ”ray-by-ray” Scheck et al. (2006) multi-energy-group neutrino schemes [e.g., multigroup flux-limited diffusion Mezzacappa and Bruenn (1993) and isotropic diffusion source approximation Liebendörfer et al. (2009)]. Even more advanced, but still rudimentary, approaches to the solution of the radiative-transfer problem have been presented recently in Ref. Radice et al. (2013).

Besides carrying out an extensive investigation of the properties and limitations of the NLS to a level of detail that, to the best of our knowledge, has not discussed before in a general-relativistic context, we also present novel results regarding the equilibrium properties of nonrotating hot NSs and how the cooling affects the stability of these configurations. In particular, we use several temperature profiles close to the one reached by protoneutron stars (PNSs) during their formation in a supernova core collapse or by the hypermassive neutron star (HMNS) produced in a BNS merger and investigate how the energy losses can influence the stability of these objects. In our idealized but well-controlled setup we can then show that deviations from the thermal and weak-interaction equilibrium of these models will affect their stability to radial perturbations and thus lead to a neutrino-induced collapse of stable nonrotating NS models.

The organization of the paper is as follows: In Sec. II we summarize our mathematical framework regarding the gravitational field equations and the GRHD equations in the presence of neutrino emission. Sec. III describes the neutrino leakage scheme while Sec. IV discusses the nuclear equations of state for NS matter used in this work. Sec. V describes in detail the numerical framework adopted for our simulations. Sec. VI, presents a series of tests and results obtained with our code. In particular, we show simulations of nonrotating stable NS with nuclear EOS at finite temperature and in -equilibrium, comparing with predictions from linear-perturbation theory. Finally, Sec. VII reports the main result of this paper, namely, the neutrino-induced collapse of a stable nonrotating NS and the associated neutrino bursts when using different thermal EOSs. Appendix A reports the details of the free emission rates and opacity sources used in our code and that could be of direct use for anyone wanting to reproduce our implementation, while a detailed description of the convergence tests of our scheme are presented in Appendix B. A detailed description of a novel and robust technique for the transformation from conserved to primitive variables is finally presented in Appendix C. Unless stated otherwise, we use geometrized units in which . Latin indices run from 1 to 3 while greek indices run from 0 to 3.

## Ii Mathematical Framework

Although much of our mathematical and numerical framework for the solution of the Einstein-Euler equations has been presented in a scattered manner in a number of papers, i.e., Baiotti et al. (2005a, 2008), we include a short review here both for completeness and as an aid to anyone wanting to implement our approach to the handling of the radiative transfer.

### ii.1 Einstein equations

We evolve a conformal-traceless ”” formulation of the Einstein equations Nakamura et al. (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro (1998), BSSNOK. The spacetime is foliated into three-dimensional spacelike hypersurfaces, with a three-metric , extrinsic curvature , and the gauge functions (lapse) and (shift) to specify the coordinate frame. The three-metric is conformally transformed as

 ϕ=112ln(detγij),~γij=e−4ϕγij (1)

and the conformal factor is evolved as an independent variable, whereas is subject to the constraint . The extrinsic curvature is subject to the same conformal transformation and its trace is evolved as an independent variable, while in place of we evolve its conformal-traceless equivalent, i.e.,

 ~Aij=e−4ϕ(Kij−13γijK), (2)

with . In addition, new evolution variables defined in terms of the Christoffel symbols of the conformal three-metric are introduced, i.e.,

 ~Γi=~γjk~Γijk, (3)

and solved as independent evolution variables. As a result, the complete set of the evolution equations for the Einstein equations is given by

 (∂t−Lβ)~γij=−2α~Aij, (4) (∂t−Lβ)ϕ=−16αK, (5) (∂t−Lβ)~Aij=e−4ϕ[−DiDjα+α(Rij−8πSij)]TF +α(K~Aij−2~Aik~Akj), (6) (∂t−Lβ)K=−DiDiα +α[~Aij~Aij+13K2+4π(E+S)], (7) ∂t~Γi=~γjk∂j∂kβi+13~γij∂j∂kβk +23~Γi∂jβj−2~Aij∂jα+2α(~Γijk~Ajk+6~Aij∂jϕ −23~γij∂jK−8π~γijSj), (8)

where is the Lie derivative with respect to the shift vector, is the three-dimensional Ricci tensor, the covariant derivative associated with the three-metric , ”TF” indicates the trace-free part of tensor objects and , , and are the matter source terms defined as

 E ≡nαnβTαβ, (9) Si ≡−γiαnβTαβ, (10) Sij ≡γiαγjβTαβ, (11)

where is the future-pointing four-vector orthonormal to the spacelike hypersurface, is the energy-momentum tensor for a perfect fluid [cf., Eq. (22)], and . The Einstein equations also comprise a set of physical constraint equations that are satisfied within each spacelike slice,

 H ≡R+K2−KijKij−16πE=0, (12) Mi ≡Dj(Kij−γijK)−8πSi=0, (13)

which are usually referred to as Hamiltonian and momentum constraints. Here is the Ricci scalar on a three-dimensional time slice. Our specific choice of evolution variables introduces two additional constraints,

 det~γij =1, (14) ~Aii =0, (15)

which are enforced on each spacelike hypersurface. The remaining constraints, and are not actively enforced and can be used to monitor the accuracy of our numerical solution.

Our gauge choices reflect the experience matured over the last decade and are already discussed in detail in Baiotti et al. (2008). In particular, we evolve the lapse according to the ”” slicing condition Bona et al. (1995):

 ∂tα−βi∂iα=−2α(K−K0), (16)

where is the initial value of the trace of the extrinsic curvature and equals zero for the maximally sliced initial data we consider here. The shift is evolved using the hyperbolic -driver condition Alcubierre et al. (2003),

 ∂tβi−βj∂jβi = 34αBi, (17) ∂tBi−βj∂jBi = ∂t~Γi−βj∂j~Γi−ηBi, (18)

where is a parameter which acts as a damping coefficient. Following  Alcubierre et al. (2003) and because we are not considering binaries here, we set to be , where is the total gravitational mass of the system; a more sophisticated prescription in case of unequal-mass binaries can be found in Alic et al. (2010).

All the equations discussed above are solved using the McLachlan code, a three-dimensional finite-differencing code based on the Cactus Computational Toolkit Löffler et al. (2012).

### ii.2 Relativistic-hydrodynamic equations

In the following, we describe the matter evolution equations implemented in the new Whisky code. These are given by the general-relativistic-hydrodynamic equations describing a perfect compressible fluid governed by a generic EOS with temperature and composition dependence. In particular, the changes in composition, energy, and momentum due to neutrino radiation are modeled by introducing additional source terms.

When written in covariant form, the equations for the conservation of energy, momentum, baryon and lepton numbers are expressed as

 ∇αTαβ =Ψβ, (19) ∇α(nbuα) =0, (20) ∇α(neuα) =N, (21)

where and are the baryon and electron number densities, respectively. Note that the energy-momentum tensor accounts for the ordinary matter as well as for the trapped neutrinos and photons, but it does not include the nontrapped neutrinos. The corresponding energy-momentum tensor is taken into account only in the form of the source term modeling the radiative losses of energy and momentum. Furthermore, because we assume these nontrapped neutrinos to behave essentially as a test fluid, the associated energy momentum tensor is neglected when building the right-hand side of the Einstein equations.

As customary in a NLS, also in our implementation we do not evolve directly the number density and the energy distribution of the neutrinos. Instead, since the neutrinos are assumed to be in local thermal equilibrium with the baryonic matter, we can obtain direct estimates for the source terms and by simply considering the matter properties. Also, because we do not consider lepton species other than electrons in our fluid, the only degree of freedom in the composition is represented by the electron fraction defined as . The electron fraction is changed only by the production of electron neutrinos and antineutrinos as described by the source term .

The energy-momentum tensor used in Eq. (19) is that of a perfect fluid, given by

 Tαβ=ρ(1+ϵ+pρ)uαuβ+pgαβ, (22)

where is the baryon rest-mass density defined as , where denotes the nucleon mass1. The total isotropic pressure and the total specific internal energy contain contributions of baryons, electrons, photons, and trapped neutrinos

 p =pe+pb+pγ+pνe,¯νe, (23) ϵ =ϵe+ϵb+ϵγ+ϵνe,¯νe, (24)

where the contribution from trapped electron-neutrinos and antineutrinos in the dense baryonic component can be evaluated from the thermodynamic state of the fluid and assuming that the neutrinos are following a Fermi-Dirac distribution, i.e.,

 pνe,¯νe ≡ pνe+p¯νe=4π3T4[F3(ηνe)+F3(η¯νe)], (25) ϵνe,¯νe ≡ ϵνe+ϵ¯νe=13pνe,¯νeρ. (26)

Here, and are the degeneracy parameters for the electron-neutrinos and antineutrinos, the corresponding chemical potentials, is the temperature, which is expressed in throughout the paper, and is a Fermi integral [see Eq. (79) and Appendix A for a more extended discussion]. Note also that in Eq. (23) the electron contributions to the total pressure and internal energy are computed using the EOS for a semidegenerate gas of electrons as described in Ref. Timmes and Swesty (2000).

At rest-mass densities and temperatures , the contributions can be significant and around of the baryonic one. However, close to nuclear saturation density, it becomes less than of the total pressure and thus smaller than the typical uncertainties of the nuclear EOSs at such densities. Furthermore, the addition of a neutrino isotropic pressure to the fluid pressure is meaningful only in the trapped regime and cannot be made when the neutrinos are free streaming. For these reasons and to maintain a smooth transition between the regimes of trapped and free-streaming neutrinos, we neglect contributions of trapped neutrinos to pressure and internal energy, setting in Eqs. (23)–(24).

To compute the neutrino source term in Eq. (21), which describes the changes in electron fraction, we introduce the neutrino emission rates per baryon, and , for electron-neutrinos and -antineutrinos, respectively. As a result, in the fluid rest frame we can express the change in electron fraction as

 uα∇α(Ye) =R≡R¯νe−Rνe. (27)

Similarly, in order to determine the source terms describing the radiative losses of energy and momentum due to neutrinos in Eq. (19), we assume that the emission is isotropic in the fluid rest frame, so that no net momentum change can take place in the rest frame. We then obtain the following covariant equation that can be most easily verified in the comoving frame

 Ψβ =−ρm−1bQuβ≡−ρm−1b∑IQIuβ =−ρm−1b(Qνe+Q¯νe+Qντ,μ+Q¯ντ,μ)uβ, (28)

where we have introduced the neutrino emissivity as the sum of the emissivities of the various neutrino species with , which denote the radiated energy per unit time and baryon (i.e., it is expressed as in cgs units). Note that we collect the emissivity due to the and neutrinos into a single contribution . An important remark should be made about Eq. (21). Since the neutrino interaction time scale depends largely on the microphysics of the process described, it can differ significantly from the time scale over which the baryonic matter evolves, resulting in a problem with a stiff source term. When this happens, special numerical algorithms are needed to handle the evolution (see, e.g., Palenzuela et al. (2009)).

For numerical evolution, the equations are cast into a flux conservative formulation, based on the Valencia formulation originally developed by Martí et al. (1991); Banyuls et al. (1997). In particular, we write Eqs. (19)–(21) as balance laws

 ∂t(√γq)+∂i(√γf(i)(q)) =s(q), (29)

where is the determinant of the three-metric, while and are the flux vectors and source terms, respectively (see Font (2008) for details). The right-hand side contains the geometrical terms and the neutrino reaction rates that depend on the state of the fluid. The evolved conserved variables are given in terms of the primitive variables via the relations

 D ≡ρW, (30) τ ≡ρhW2−p−D, (31) Si ≡ρhW2vi, (32) ^Ye ≡DYe. (33)

In the equations above, is the three-velocity measured by an observer orthogonal to the hypersurface of constant coordinate time (i.e., a Eulerian observer), is the corresponding Lorentz factor, and is the specific enthalpy. When written out explicitly, Eq. (29) reads

 ∂t¯D =−∂i(wi¯D), (34) ∂t¯τ=−∂i(wi¯τ+√γαpvi)−αQm−1b¯D+α¯SklKkl−¯Si∂iα, (35) ∂t¯Sj=−∂i(wi¯Sj+√γαpδij)−αQm−1b¯Dvi+α2¯Skl∂jγkl+¯Sk∂jβk−(¯τ+¯D)∂jα, (36) ∂t¯Ye =−∂i(wi¯Ye)+¯DαWR, (37)

where , , , , , and denotes the advection speed with respect to the coordinates.

In order to close the system of equations, we need an EOS to relate the pressure to the rest of the primitive quantities (see Sec. IV). The new Whisky code implements several analytic EOSs, such as the polytropic EOS, the ideal-fluid (-law) EOS, as well as ”hybrid” EOSs Janka et al. (1993). Finally, the code can use microphysical finite-temperature EOSs in tabulated form to describe nuclear matter and the neutrino interactions, as will be discussed in Sec. IV.

The following sections are devoted to a description of our treatment of the radiative transport within the NLS and of the interactions between neutrinos and ordinary matter when the latter is modeled with the nuclear EOS described in Sec. IV.

## Iii Approximate radiative transfer: the neutrino leakage scheme

Neutrinos from hot, dense and neutron-rich nuclear material are produced via nonequilibrium weak-interaction processes during the first minutes (or months) of the life of a PNS and in the final stages of the merger of BNSs, when a HMNS is produced. The cooling produced by these neutrinos and the consequent influence on the stellar structure and equilibrium is of course of great importance to model these processes accurately. In order to study the possible influence of thermal effects on the onset of dynamical instabilities in PNSs or HMNSs, it is necessary to evolve the system for tens of dynamical time scales, which can be estimated to be ms, where and are the mass and size of the object. This time scale is short when compared with the typical diffusion time scale of neutrinos in dense matter, which is s. As a result, it is reasonable to use as a first approximation to the radiative transfer a simple leakage scheme that estimates the instantaneous energy loss via neutrino emission and by the total change in electron fraction of the nuclear matter. Clearly, this represents only a first step towards more sophisticated and accurate transport schemes such as those employed in stellar-core collapse simulations Scheck et al. (2006); Marek and Janka (2009); Liebendörfer et al. (2009); Abdikamalov et al. (2012).

The NLS was originally developed by van Riper et al. van Riper and Lattimer (1981), to describe the neutrino cooling through weak interactions during the core-collapse phase of a supernova. In Newtonian gravity, it was first applied in BNS simulations by Ruffert et al. Ruffert et al. (1996, 1997); Ruffert and Janka (1999, 2001) and more recently by Rosswog et al. (Rosswog and Davies, 2003; Rosswog and Liebendörfer, 2003; Rosswog et al., 2003). It has undergone a revival in general-relativistic three-dimensional supernova simulations with the work of Sekiguchi (2010) and O’Connor and Ott (2010); Ott et al. (2012, 2013) (see also Peres et al. (2013)). Finally, the NLS was also used recently to describe the neutrino cooling in general-relativistic BNS simulations (Sekiguchi et al., 2011a, b) and BH-NS mergers Deaton et al. (2013). The details of the implementation of the NLS differ among these authors and are sometimes fragmentary. For these reasons, and to help anyone wishing to reproduce our implementation, we give below a detailed account of our treatment.

In our scheme we account for three different neutrino species: electron neutrinos , electron-antineutrinos and heavy lepton neutrinos, that are considered as a single component with a statistical weight of 4. The creation of electron-neutrinos and antineutrinos by processes involving neutrons () and protons () leads to changes in the electron number, i.e.,

 e++n→p+¯νe, (38) e−+p→n+νe. (39)

that need to be taken into account. On the other hand, all three species of neutrinos and antineutrinos are created by annihilation of electron-positron pairs, i.e.,

 e++e−→¯νe+νe, (40) e++e−→¯ντ,μ+ντ,μ, (41)

and, especially at high temperatures, by plasmon decay, i.e.,

 γ→νe+¯νe, (42) γ→ντ,μ+¯ντ,μ. (43)

In order to compute the emission coefficients, we assume that the neutrinos are massless and in thermal equilibrium with the surrounding matter. As a result, the energy spectrum of the neutrinos follows a Fermi-Dirac distribution for ultrarelativistic particles at the temperature of the matter. A further assumption concerns the chemical potentials of the electron-neutrinos and of the electron-antineutrinos, which are assumed to be at -equilibrium Rosswog and Liebendörfer (2003), i.e.,

 μeqνe=μe−μn−μp=−μeq¯νe, (44)

where , , are the relativistic chemical potentials including the rest-mass of the particle. As a result, we can define , where now and we assume since the heavy neutrinos are only weakly interacting with the matter. The detailed derivation of the number and energy rates for the different emission processes is described in detail in Appendix A.

As mentioned in the Introduction, low-energy neutrinos (i.e., with ) interact very rarely at low densities, i.e., for . In this regime, if is the mean free path of a neutrino, the time needed by a neutrino of species to cross a slab of nuclear matter of depth , is simply given by

 tI∼L, (45)

where, we recall, . In this limit, the amount of energy loss by neutrinos can be directly computed using the free emission rates (see below and Appendix A for details). On the other hand, at high densities, i.e., for , the scattering and absorption of neutrinos onto baryons is extremely efficient already for neutrinos with and a purely diffusive regime can be reached, where . Under these conditions, neutrinos cannot freely escape but can only diffuse via a random walk over a time scale

 tI Extra open brace or missing close brace (46)

where is a coefficient depending on the optical depth evaluated along the direction of propagation Mihalas and Mihalas (1984). In order to compute the total mean free path, in our NLS we take into account scattering and absorption processes between neutrinos and the surrounding matter. In particular, for the calculation of the opacities we include the following processes:

• coherent neutrino scattering on heavy nuclei (with atomic mass number )

 νe+A →νe+A, ¯νe+A →¯νe+A, (47) ντ,μ+A →ντ,μ+A, ¯ντ,μ+A →¯ντ,μ+A, (48)
• neutrino scattering on free nucleons

 νe+[n,p] →νe+[n,p], (49) ¯νe+[n,p] →¯νe+[n,p], (50) ντ,μ+[n,p] →ντ,μ+[n,p], (51) ¯ντ,μ+[n,p] →¯ντ,μ+[n,p], (52)
• electron-flavor neutrinos absorption on free nucleons

 νe+n →p+e−,¯νe+p→n+e+ (53)

Summing all the various contributions to the mean free paths , we define the optical depth simply as the integral of the inverse of the mean free path along a line of sight which is specifically chosen for the geometry of the problem, i.e.,

 τI(EI)=∫x2x1dsλI(EI), (54)

where is the neutrino energy (see Appendix A for details), and is the proper length. Note we have to neglect the gravitational redshift of the neutrinos, as well as special-relativistic corrections due to the velocity of the absorbing matter. The knowledge of the optical depth has two important applications. First, it can be used to determine the extent of the optically thick regions, which are clearly different for different neutrino species. Second, it can be used to estimate the time scale a neutrino needs to diffuse out of the dense matter region. A simple estimate can be obtained from Eq. (46) by replacing with the optical depth , yielding

 tI(EI) =DλI(EI)τ2I(EI). (55)

This expression is also derived in Ref. Rosswog and Liebendörfer (2003) from a one-dimensional stationary diffusion model, and a similar formula is used by Ref. Ruffert et al. (1997). Results obtained from those prescriptions have been compared with a neutrino flux-limited diffusion transport scheme discussed in Ref. Baumgarte et al. (1996), providing also a calibration of the geometry parameter, which takes a value . Ref. Rosswog and Liebendörfer (2003) also proposes an alternative recipe given by

 tI(EI) =λI(EI)(τI(EI)+12)(2τI(EI)+1). (56)

In practice, it is hard to assess which of the two prescriptions is the most realistic one. The first one has the advantage that the energy dependency can be factored out, which will be crucial for the NLS described below. Considering the approximate nature of the NLS and the fact that the coefficient is itself determined empirically, we use hereafter only the prescription (55).

It should also be noted that the time scale over which our simulations are performed is typically at least 1 order of magnitude smaller than the typical diffusion time scale on which the neutrinos escape from the NS dense core, i.e.,  s. As a result, the contribution of the diffusive leakage of neutrinos is expected to be much smaller than the abundant emission of neutrinos at densities below the ones of the neutrinosphere. For this reason, the NLS can be considered as a reasonable first approximation to treat neutrino effects without resorting to more sophisticated transport schemes.

As shown in Appendix A, we can separate algebraically the neutrino energy dependence from the mean free path, introducing an inverse mean free path which does not depend on the energy, i.e.,

 ζI≡1E2IλI. (57)

Similarly, we can define an energy-independent optical depth as

 τI(EI)=E2IχI≡E2I∫x2x1ζI(x) ds. (58)

Using Eq. (55), we can now also factor out the energy dependence from the diffusion time scale. With these definitions, we can finally compute the diffusive number and energy emission rates per baryon as

 RDI(ηeqI) =ζIDχ2I∫∞0nI(E)E2dE, (59) QDI(ηeqI) =ζIDχ2I∫∞0nI(E)EdE, (60)

where, we recall, we assume the neutrino number densities to be described by the Fermi-Dirac distribution for ultrarelativistic particles. The integrals over energy can then be carried out analytically. The diffusive emission rates (59) and (60) approximate the real emission rates only in the diffusive regime. We will use them below in a prescription interpolating between optically thin and optically thick regimes.

The assumption of chemical equilibrium also impacts the opacities and the radii of the neutrinospheres, which are in general slightly overestimated in our approach. More specifically, as the region transparent to neutrinos is approached, the chemical potentials tend to vanish as a result of the decreased interaction of the nuclear matter with the neutrino flow. Clearly, this behavior cannot be reproduced when adopting the assumption of chemical equilibrium and in Ref. Ruffert et al. (1997) this problem was overcome by using an iterative procedure that adjusted the chemical potential of the neutrinos between the two regions. Here we follow instead the approach suggested more recently in Ref. Rosswog and Liebendörfer (2003) and use a simple interpolation of the emission rates at chemical equilibrium between the purely diffusive regime and the free-streaming one. This approach simplifies considerably the numerical scheme, allowing us to tabulate all the emission rates and opacities, leaving the computation of the optical depth as the only operation to be performed during the evolution (see Appendix A for details).

As a result, using the free neutrino emission rate per baryon, , and neutrino luminosity per baryon, , and the corresponding quantities for the diffusive emission and , we can now compute the effective rates, and via the following interpolation formulas

 RI(ηI) = RFI(ηeqI)(1+RFI(ηeqI)RDI(ηeqI))−1, (61) QI(ηI) = QFI(ηeqI)(1+QFI(ηeqI)QDI(ηeqI))−1. (62)

Alternatively, we have also implemented a more elaborated interpolation technique as suggested in Sekiguchi (2010)

 RI(ηI) = RFI(0)e−τI+RDI(ηeqI) (1−e−τI), (63) QI(ηI) = QFI(0)e−τI+QDI(ηeqI) (1−e−τI). (64)

This new prescription does not seem to lead to significant changes but we reserve a more careful comparison of the two prescriptions to future work.

After computing the effective emission rates, we can define a quantity that will be very useful as a diagnostic tool to estimate where the neutrino cooling is most effective. In particular, for each neutrino species, we define as the neutrinosphere the surface at which

 QIQFI=23. (65)

Clearly, inside this region the neutrino cooling is strongly suppressed and the time scale of the emission is dominated by the diffusion processes. Note that other definitions are possible and sometimes the neutrinosphere is defined as the surface at which the spectrally averaged optical depth is equal to (see, for instance, Rosswog and Liebendörfer (2003)).

Finally, the interpolated quantities can then be used to compute the net emission rate per baryon appearing in Eq. (37)

 R≡R¯νe(η¯νe)−Rνe(ηνe)=Rpc(η¯νe)−Rec(ηνe), (66)

and the luminosity per baryon appearing in Eqs. (35)–(36)

 Q≡∑IQI(ηI)=Qpc(η¯νe)+Qec(ηνe)+∑I[Qe+e−(ηI)+Qγ(ηI)], (67)

where are the emission rates from proton capture (), the emission rates from electron capture (), while and are the emission rates per unit mass from pair annihilation and plasmon decay, respectively. For the sum over these last two rates, here our notation is compact but not perfect; indeed the term in Eq. (67) should be read as

 ∑IQe+e−(ηI)=Qe+e−(ηνe,η¯νe)+Qe+e−(ηντ,μ,η¯ντ,μ), (68)

and the same for the plasmon-decay rates.

While the assumption of chemical equilibrium can be considered natural in the view of the many approximations already taken in a NLS, we should also caution that this assumption could overestimate considerably the emission rates, especially in the transition between the optically thin and thick region.

The effective local neutrino source terms are also used to estimate the total luminosity as the neutrino energy that is emitted between the coordinate times and , where the neutrino energy emitted at each point is the one measured by a local coordinate observer with worldline tangent to the four-vector . This observer is somewhat arbitrary, unless the spacetime is stationary with Killing vector . In this case, and assuming no further interactions, the energy of a neutrino measured by such observer is related to the energy measured when the neutrino reaches infinity by the simple relation , regardless of the direction of the emission. To compute the local luminosity, we use the effective emission rates, which are such that the contributions from the optically thick regions are strongly suppressed. The emission rates are given in the fluid rest frame, where a 4-momentum is emitted isotropically per baryon and unit proper time . Such a 4-momentum corresponds to an energy which is measured by our coordinate observers. The relation between the fluid proper time and the coordinate time is , and the number of baryons per coordinate volume is . Using these expressions, we define for each neutrino species the ”source luminosity” , which does not take into account the gravitational redshift, as

 LI =∫~fQI¯Dm−1bαd3x, (69)

where the integral is taken over the whole computational domain, and . In addition, we can compute the luminosity observed at infinity, , by replacing in Eq. (69) with , so as to account for the gravitational redshift. The expression for is strictly valid only for a stationary spacetime with Killing vector . For our nonrotating stationary stellar models, and thus . We will adopt this expression also for oscillating NSs, in which case the dominant error is caused by the fact that a significant fraction of the emitted neutrinos, even those emitted in the free-streaming regime, will hit the NS core and never reach infinity. This could easily reduce the luminosity by a factor around two.

## Iv Nuclear equations of state for neutron-star matter

The importance that the EOS has in our understanding of the properties of NSs is so large that its study requires no justification. Besides determining the composition and structure of NSs, the EOS also plays a fundamental role in defining the properties of the GW signal that is expected to be produced by NSs, either when isolated or in a binary system. A number of studies about the inspiral and merger of BNS have made this point very clearly, both when the EOS was idealized and analytic, e.g., Baiotti et al. (2007, 2008); Rezzolla et al. (2010), or when more realistic ones were considered, e.g., Freiburghaus et al. (1999); Shibata et al. (2005); Bauswein et al. (2010). Unfortunately, however, our knowledge of matter at nuclear densities is still plagued by too many uncertainties, making it impossible to derive, from robust first principles, what is the most realistic EOS for NS matter. These uncertainties leave ample room for a large variety of possible models for the nuclear matter at zero temperature, which is only mildly constrained by astronomical observations.

In practice, however, only a few nuclear EOSs are able to encompass the wide range of temperatures, compositions and rest-mass densities necessary to describe matter in supernova explosions and BNS mergers. Here, we will make use of two nuclear EOSs for hot dense matter: the Lattimer-Swesty (LS) EOS Lattimer and Swesty (1991), that adopts a compressible liquid drop model with Skyrme interaction Tondeur et al. (1984), and the Shen-Horowitz-Teige (SHT) EOS, that adopts a relativistic mean-field model for uniform matter with a modified NL3 set of interaction parameters (see Refs. Shen et al. (2010, 2011)). These two EOSs do not include contributions from exotic matter like quarks, hyperons, and kaon or pion condensates. Furthermore, they assume charge neutrality.

In these publicly available EOSs, the pressure and the specific internal energy are given as a function of three independent variables: the baryon number density , the electron/proton fraction , and the temperature . The SHT-EOS is given in a tabulated form, while the LS-EOS is provided in form of computer routines, which provide several choices for the incompressibility modulus, namely . For our simulations we concentrate on a value and cast also the LS-EOS into a tabular form using the routines described in Ref. O’Connor and Ott (2010). Note that the original table for the SHT-EOS has a fairly low resolution and we compensate this by computing a larger table via a thermodynamically consistent interpolation of the Helmholtz free energy Swesty (1996); Timmes and Swesty (2000). In a separate table the SHT-EOS provides also the slice. Since we store the logarithm of the temperature, we cannot use directly this slice, but we use it to extend the original temperature range via linear interpolation. In Table 1, we summarize the ranges covered by the two EOS tables used in our simulations and their resolution.

When it comes to the range in composition, the SHT-EOS can describe the properties of pure neutron matter, that is , but this case is provided as a separate table. To compute a regularly spaced table covering the full range, we fill the gap into the original table using a third-order polynomial interpolation. For the LS-EOS, on the other hand, the smallest available value for the electron fraction is and is taken from the table described in Ref. O’Connor and Ott (2010), and as produced by the routines provided by the LS-EOS. The regime of low electron fraction becomes relevant for NSs at low temperature and -equilibrium, which show a pronounced dip in the electron fraction at the rest-mass densities typical of the inner crust i.e., .

We note that the original table for the SHT-EOS contains some problematic regions where derived quantities, such as the speed of sound , are extremely noisy or even unphysical. Since our numerical evolution scheme requires a well-behaved sound speed, an additional smoothing step was necessary to remove these irregularities. In particular, we compute the sound speed at constant composition as

 c2s=1h⎡⎣∂p∂ρ∣∣∣T,Ye+pρ2∂p∂T∣∣∣ρ,Ye(∂ϵ∂T∣∣∣ρ,Ye)−1−∂p∂T∣∣∣ρ,Ye∂ϵ∂ρ∣∣∣T,Ye(∂ϵ∂T∣∣∣ρ,Ye)−1⎤⎦, (70)

where the required derivatives are computed using the limiting procedure described in Steffen (1990). On the other hand, the sound speed for the LS-EOS can be computed analytically from the derivatives of the Helmholtz free energy, yielding smoother results.

Both the LS-EOS and the SHT-EOS tables provide the necessary microphysical quantities to compute the neutrino interactions such as the nucleon and electron chemical potentials, the abundances of heavy nuclei and alpha particles, the proton to nucleon ratio , and the atomic mass number for the ”representative” heavy nuclei.

Fig. 1 illustrates the dependence of the pressure and specific energy as a function of the three independent variables for the SHT-EOS (the LS-EOS shows a similar behavior). The left panel shows the pressure as a function of the rest-mass density for three different temperatures (distinguished by different colors), while the right panel shows the specific internal energy as a function of the temperature for three different densities. Both panels report the influence of the electron fraction as shaded regions that represent the variation over the whole range of . As one can see, the influence of the composition on the pressure is quite large at low temperatures and densities below the crust-core transition, i.e., at . In the core region, on the other hand, the pressure is relatively unaffected by changes in composition and temperature. Note also that the specific internal energy becomes almost temperature independent at densities larger than nuclear saturation density, in stark contrast with simplified analytic EOSs such as the -law EOS.

We recall that the numerical scheme we use to solve the relativistic-hydrodynamic equations provides us at each time step with the rest-mass density , the electron fraction , and internal energy , but not the temperature . The EOS tables, on the other hand, use and as independent variables. As a result, a conversion is needed between and for any given couple of . This is done during the evolution by a discrete bisection algorithm that finds the nearest tabulated temperatures and then performs a linear interpolation.

At regions with high density but low temperature, the conversion becomes very inaccurate because of the very weak dependence of on (cf., right panel of Fig. 1), with the consequence that small errors in can be strongly amplified. Fortunately, this does not affect the hydrodynamic evolution since the pressure shows the same weak dependence on temperature. On the other hand, the local neutrino emission rates do depend strongly on the temperature, i.e., as , so that the free and the effective emission rates can be computed with a large relative error. Luckily, also in this case the errors do not have a serious impact because under these high-density conditions the neutrinos are trapped and the resulting effective luminosity is small. As a result, the influence of inaccurate emission rates on the total luminosity is small as long as the latter is dominated by the emission from hot matter outside the optically thick region.

## V Numerical Framework

Much of our numerical infrastructure has been discussed in other papers, e.g., Baiotti et al. (2008); Rezzolla et al. (2010), to which we refer the interested reader for details. Below, we only give a brief overview on the numerical methods used. For the time integration of the coupled set of the hydrodynamic and Einstein equations we use the method of lines (MOL) in conjunction with an explicit fourth-order Runge-Kutta method. In all our simulations we prescribe a Courant-Friedrichs-Lewy (CFL) factor of to compute the time step. The Einstein equations are spatially discretized using a fourth-order finite-difference operator, implemented by the publicly available McLachlan code. The hydrodynamic equations, on the other hand, are discretized in space using a finite-volume scheme based on a piecewise parabolic reconstruction (PPM) Colella and Woodward (1984) and the Harten-Lax-van Leer-Einfeldt (HLLE) Harten et al. (1983) approximate Riemann solver. Differently from the original PPM implementation reported in Ref. Baiotti et al. (2007), we reconstruct here the quantities instead of the three-velocities . This guarantees that the velocities reconstructed at the cell boundaries remain subluminal even under extreme conditions like those encountered near the central region of a BH. As is customary when solving the relativistic-hydrodynamic equations, the vacuum regions are treated by enforcing an artificial atmosphere with a rest-mass density and temperature close to the lowest values covered by the EOS considered (see Table 1) and with a constant electron fraction of .

A discussion on the global convergence order of the simulations reported here is presented in Appendix B and is generally of order , in agreement with what was found also in previous calculations using the same code with ideal-fluid EOS Baiotti et al. (2009).

### v.1 Adaptive-mesh refinement and singularity handling

The use of mesh-refinement techniques is of fundamental importance in our simulations. For this reason, we use the Carpet driver that implements a vertex-centered adaptive-mesh-refinement scheme adopting nested grids Schnetter et al. (2004). For the evolution of isolated NSs presented in this work, it is sufficient to use a fixed hierarchy of nested boxes centered around the origin, with a 2:1 refinement factor between successive grid levels.

In some of the simulations presented here, the final state of the evolution is a BH. We then use the isolated horizon framework Thornburg (2004) to measure its properties every few time steps. We do not make use of the excision technique Baiotti et al. (2005b) in our simulations, neither for the spacetime variables nor for the fluid. As described in Refs. Baiotti and Rezzolla (2006); Baiotti et al. (2007), in order to extend the simulations well past the BH formation, we add a small amount of dissipation to the evolution equations for the metric and gauge variables and rely on the singularity-avoiding gauge (16) (note that no dissipation is added to the evolution of matter variables). More specifically, we use an artificial dissipation of the Kreiss-Oliger-type Kreiss and Oliger (1973), i.e., we augment the right-hand side of the evolution equation of a given quantity with a term , where is the grid spacing. The dissipation is active on the whole computational domain with constant strength . The use of a numerical viscosity is necessary because all the field variables develop very steep gradients in the region near the BH center. Under these conditions, small high-frequency oscillations can lead to instabilities if not dissipated.

### v.2 optical depth grid

Different prescriptions are possible in order to estimate the optical depth that is needed in the NLS. A simple choice would be to perform the integrals (58) along the Cartesian coordinate lines and then take the minimum value over all directions. A better choice is to calculate the integrals along directions that are better adapted to the underlying geometry of the matter distribution. In this respect, since we are mainly interested in objects with an overall spherical geometry, a spherical polar coordinate system is better suited than a Cartesian one. Hence, we introduce a spherical grid and compute along a number of a radial coordinate directions the integrals of the energy-independent optical depths , i.e.,

 χI(r,θ,ϕ) =∫rmaxrζ