Low mass binary neutron star mergers : gravitational waves and neutrino emission
Francois Foucart, Roland Haas, Matthew D. Duez, Evan O'Connor, Christian D. Ott, Luke Roberts, Lawrence E. Kidder, Jonas Lippuner, Harald P. Pfeiffer, Mark A. Scheel
LLow mass binary neutron star mergers : gravitational waves and neutrino emission
Francois Foucart, Roland Haas, Matthew D. Duez, Evan O’Connor, Christian D. Ott, LukeRoberts, Lawrence E. Kidder, Jonas Lippuner, Harald P. Pfeiffer,
7, 8 and Mark A. Scheel Lawrence Berkeley National Laboratory, 1 Cyclotron Rd, Berkeley, CA 94720, USA Max-Planck-Institut fur Gravitationsphysik, Albert-Einstein-Institut, D-14476 Golm, Germany Department of Physics & Astronomy, Washington State University, Pullman, Washington 99164, USA Department of Physics, North Carolina State University, Raleigh, North Carolina 27695, USA TAPIR, Walter Burke Institute for Theoretical Physics, MC 350-17, California Institute of Technology, Pasadena, California 91125, USA Cornell Center for Astrophysics and Planetary Science, Cornell University, Ithaca, New York, 14853, USA Canadian Institute for Theoretical Astrophysics, University of Toronto, Toronto, Ontario M5S 3H8, Canada Canadian Institute for Advanced Research, 180 Dundas St. West, Toronto, ON M5G 1Z8, Canada
Neutron star mergers are among the most promising sources of gravitational waves for advanced ground-based detectors. These mergers are also expected to power bright electromagnetic signals, in the form of shortgamma-ray bursts, infrared/optical transients powered by r-process nucleosynthesis in neutron-rich materialejected by the merger, and radio emission from the interaction of that ejecta with the interstellar medium.Simulations of these mergers with fully general relativistic codes are critical to understand the merger andpost-merger gravitational wave signals and their neutrinos and electromagnetic counterparts. In this paper, weemploy the Spectral Einstein Code (SpEC) to simulate the merger of low-mass neutron star binaries (two . M (cid:12) neutron stars) for a set of three nuclear-theory based, finite temperature equations of state. We show that thefrequency peaks of the post-merger gravitational wave signal are in good agreement with predictions obtainedfrom recent simulations using a simpler treatment of gravity. We find, however, that only the fundamental modeof the remnant is excited for long periods of time: emission at the secondary peaks is damped on a millisecondtimescale in the simulated binaries. For such low-mass systems, the remnant is a massive neutron star which,depending on the equation of state, is either permanently stable or long-lived (i.e. rapid uniform rotation issufficient to prevent its collapse). We observe strong excitations of l = 2 , m = 2 modes, both in the massiveneutron star and in the form of hot, shocked tidal arms in the surrounding accretion torus. We estimate theneutrino emission of the remnant using a neutrino leakage scheme and, in one case, compare these results witha gray two-moment neutrino transport scheme. We confirm the complex geometry of the neutrino emission,also observed in previous simulations with neutrino leakage, and show explicitly the presence of importantdifferences in the neutrino luminosity, disk composition, and outflow properties between the neutrino leakageand transport schemes. PACS numbers: 04.25.dg, 04.40.Dg, 26.30.Hj, 98.70.-f
I. INTRODUCTION
Compact binary mergers, including binary neutron stars(BNS), binary black holes (BBH), and neutron star-black hole(NS-BH) mergers are the primary targets of ground basedgravitational wave detectors such as Advanced LIGO [1], Ad-vanced VIRGO [2], and KAGRA [3]. In the presence of atleast one neutron star, the merger may also be accompanied bybright electromagnetic emission, whose detection can providevaluable complementary information on the properties of thesource, help characterize the merger environment, and pro-vide better localization than available from the gravitationalwave signal alone (see e.g. [4] for a review). Two of themost promising such counterparts are short gamma-ray bursts(SGRBs) [5–8] and radioactively powered transients originat-ing from r-process nucleosynthesis in the neutron-rich matterejected by the merger [9–13]. The latter would most likelypeak in the infrared about a week after merger [13–15]. Itcould also result in the production of many of the heavy ele-ments (
A > ∼ ) whose origin remains poorly understood to-day [11, 16, 17]. Finally, the ejecta could power long durationradio emission as it interacts with the interstellar medium.Numerical simulations of these mergers play an importantrole in the efforts to model the gravitational wave signal, pre- dict the properties of its electromagnetic counterparts, and es-timate the production of various elements. In this work, wefocus on the outcome of BNS mergers. The simulation ofBNS mergers with general relativistic hydrodynamics codeshas now been possible for about 15 years [18]. Nevertheless,and despite continuous improvements, current codes have notyet reached the accuracy required to model the gravitationalwave signal at the level required to extract as much informa-tion as possible from upcoming experiments (see e.g. [19]).Additionally, most codes do not take into account all of thephysics relevant to the evolution of the post-merger remnant,including at least a hot nuclear-theory based equation of state,a neutrino transport scheme accounting for neutrino-matterand neutrino-neutrino interactions [20, 21], and the evolutionof the magnetic fields with a high enough resolution to re-solve the growth of magnetohydrodynamics (MHD) instabili-ties [22, 23].Many recent simulations have taken steps toward an im-proved treatment of all relevant physics, and a better coverageof the available parameter space of BNS configurations. Fullygeneral relativistic simulations of the merger of two . M (cid:12) neutron stars with an approximate neutrino transport scheme,but no magnetic fields [20], have shown that neutrino-matterinteractions can lead to the ejection of material with a broad a r X i v : . [ a s t r o - ph . H E ] J a n distribution of temperature and composition, leading to theproduction of elements with abundances compatible with so-lar system observations [17]. The study of a similar merger atan unprecedented resolution with ideal MHD but using a sim-pler equation of state and no neutrinos [22] demonstrated mul-tiple mechanisms for the growth of magnetic field in globalsimulations of mergers – although, even at the highest resolu-tion currently achieved, the amplification of the magnetic fieldin the shear region between the two merging neutron stars isnot resolved [23], and the unresolved growth of magnetic fieldin that region could be very significant [23, 24]. Lower res-olution simulations including resistive MHD have shown thatdeviations from ideal MHD have a relatively modest effect onthe evolution of the post-merger remnant [25]. With or with-out resistivity, Dionysopoulou et al. [25] also find low-density,magnetically dominated outflow regions which may lead tothe formation of a relativistic jet. The combined effect of idealMHD with a subgrid model for the growth of the magneticfield, and a simple neutrino cooling scheme has also been in-vestigated [26, 27]. These results suggest that the growth ofstrong magnetic field could influence the mass and propertiesof the outflows, and the dynamics of the post-merger remnant.Large parameter space studies without microphysics [28–33] have also greatly improved our understanding of the gen-eral properties of the post-merger remnant. If the post-mergerremnant collapses to a black hole, the relevant timescale forthat collapse can be estimated from the total mass of the bi-nary. Order of magnitude estimates are available for the massejection from the merger, and we are beginning to understandthe main features of the post-merger gravitational wave sig-nal. More exotic merger scenarios, leading to more extrememass ejection, merger results, and orbital evolution have alsobeen investigated, including eccentric binaries [34, 35], andbinaries with rapidly rotating neutron stars [36–40]. Finally, alarge number of simulations have also been performed usingapproximate treatments of gravity, producing useful predic-tions for the properties of the post-merger gravitational wavesignal [41, 42], as well as the long term evolution of the post-merger remnant [43–47]. All of these, and many more earlierresults, provide us with a growing understanding of the be-havior of neutron star mergers (see also [48–50] for reviewsof earlier results), but certainly do not yet draw a completepicture.In this paper, we focus on an often neglected portion ofthe parameter space in general relativistic simulations: neu-tron star binaries at the very low end of the expected massdistribution [51]. We consider two . M (cid:12) neutron stars, andevolve them using the SpEC code [52], a set of three differenthot, nuclear theory based equations of state, and either a neu-trino cooling scheme or a gray two-moment neutrino trans-port scheme capable of properly taking into account neutrino-matter interactions [21]. A similar system has been evolvedusing ideal MHD and a simpler equation of state (leading,in particular, to an unrealistically large neutron star radius)in [24], with the aim to study the growth of magnetic field,and with similar equations of state but approximate gravityand no neutrinos in [53]. Nuclear theory-based equations ofstate and a neutrino cooling scheme have also recently been used for higher mass systems in [26, 27], providing us witha useful point of comparison for many of the qualitative fea-tures observed in our simulations. The post-merger gravita-tional wave signal can also be compared to numerical resultsusing nuclear theory-based equations of state but an approxi-mate treatment of gravity [41, 42].We describe our numerical methods in Sec. II, and initialconditions in Sec. III. Our simulations allow us to address anumber of important questions. First, in Sec. IV, we study thepost-merger gravitational wave signal of low-mass BNS sys-tems with hot nuclear-theory based equations of state. This al-lows us to compare our results with predictions coming fromsimulations using a more approximate treatment of gravity.Those simulations found clearly marked peaks in the gravi-tational wave signal at frequencies easily connected with theneutron star equation of state [41, 42]. We can also compareour results to differing predictions coming from general rela-tivistic simulations using simpler equations of state and moremassive neutron stars [31, 32]. The post-merger signal is ex-pected to depend on the equation of state of the neutron star,and its qualitative properties vary with the total mass of thesystem. Thus, our results using general relativistic simula-tions, a nuclear-theory based equations of state, and low-massneutron stars usefully complement existing studies in generalrelativity [29, 31, 32] and approximate gravity [42].Second, in Sec. V, we study the qualitative features of theremnant: its density, temperature, composition, rotation pro-file, and the excitation of l = 2 , m = 2 mode in the post-merger neutron star and the accretion torus. We compare ourresults with simulations of higher mass systems using a simi-lar level of physical detail [27].Third, in Sec. VI, we analyze neutrino emission, and for thesimulation in which a transport scheme is used, directly as-sess the impact of using the simpler neutrino cooling method,determine the changes in composition due to neutrino-matterinteractions, and estimate how these changes affect the prop-erties of the outflows (Sec. VII). This is the first time that sucha comparison is possible in binary neutron star mergers (we al-ready performed a similar study on the accretion disk formedin a NS-BH merger [21]). This provides us with better in-sights into the limits of the simple cooling scheme. We canalso qualitatively compare our results with the small set of ex-isting neutron star merger simulations using a similar neutrinotransport scheme [20]. II. NUMERICAL METHODSA. Evolution equations
We evolve Einstein’s equations and the general relativis-tic equations of ideal hydrodynamics using the Spectral Ein-stein Code (SpEC) [52]. SpEC evolves those equations ontwo separate grids: a pseudospectral grid for Einstein’s equa-tions, written in the generalized harmonic formulation [54],and a finite volume grid for the general relativistic equationsof hydrodynamics, written in conservative form. The lattermakes use of an approximate Riemann solver (HLL [55])and high-order shock capturing methods (fifth order WENOscheme [56, 57]), resulting in a second-order accurate evo-lution scheme. For the time evolution, we use a third-orderRunge-Kutta algorithm. Finally, after each time step, the twogrids communicate the required source terms, using a third-order accurate spatial interpolation scheme. Those sourceterms are the metric and its derivatives (from the pseudospec-tral grid to the finite volume grid) and the stress-energy tensorof the fluid (from the finite volume grid to the pseudospec-tral grid). In our current scheme, the radiation stress-energytensor does not self-consistently feed back onto the evolutionof Einstein’s equations. Direct measurements of the energyin the neutrino sector shows that the radiation energy densityis everywhere negligible as far as gravitational interactionsare concerned. The neutrino stress-energy tensor is, however,fully coupled to the general relativistic equations of hydrody-namics. More details about these numerical methods can befound in [58, 59].While the GR-hydrodynamics in SpEC was largely devel-oped for the evolution of NS-BH mergers, most of the algo-rithm carries over to the evolution of binary neutron star merg-ers. The main differences are in the choice of the evolutiongauge, and in the grid structure (Sec. II B). In the generalizedharmonic formalism, the evolution of the coordinates followsthe wave equation g ab ∇ c ∇ c x b = H a ( x ) , (1)where g ab is the spacetime metric, and H a ( x ) an arbitraryfunction. Before merger, we make the choice H a ( x ic , t ) = H a ( x ic ,
0) exp (cid:18) − t τ (cid:19) , (2)with τ = (cid:112) d /M , d the initial separation, and M the to-tal mass of the binary at infinite separation. Here, x ic are thecomoving spatial coordinates (which follow the rotation andinspiral of the binary) and H a ( x ic , are the value of H a atthe initial time, chosen so that ∂ t α ( x ic , t ) = ∂ t β i ( x ic , t ) = 0 .The lapse α and shift β i are obtained from the standard 3+1decomposition of the metric, g ab = − α dt + γ ij ( dx i + β i dt )( dx j + β j dt ) , (3)where γ ij is the 3-metric on a slice of constant t . After a tran-sient on the timescale τ , the gauge settles into the ”harmonicgauge” H a = 0 (see [60] for a discussion of the uses of thatgauge). During merger, we instead transition to the ”dampedharmonic” gauge condition [61], H DH a = (cid:20) log (cid:18) √ γα (cid:19)(cid:21) (cid:20) √ γα t a − γ ai β i α (cid:21) , (4)where γ is the determinant of the 3-metric γ ij . We do thisaccording to H a ( t ) = H DH a (cid:20) − exp (cid:18) − ( t − t DH ) τ m (cid:19)(cid:21) , (5)with t DH the time at which we turn on the damped harmonicgauge, discussed in Sec. II B, and τ m = 100 M . FIG. 1. Equatorial cut of the pseudospectral grid just before the twoneutron stars get into contact. The color scale shows the baryon den-sity (logarithmic scale).
We choose H a based on what we found most efficient inNS-BH mergers, in terms of the resolution required to reacha given accuracy. There, we find that the harmonic gauge is abetter choice for the evolution of neutron stars during inspiral.The damped harmonic gauge, on the other hand, is favoredfor the evolution of black holes and for very dynamical space-times. For the collapse of a merger remnant into a black hole,different gauge choices are necessary (R. Haas et al., in prep.).However, we do not have to take this into consideration here,since the low-mass systems that we study do not collapse to ablack hole within the duration of the simulation. B. Grid setup
Before the two neutron stars enter into contact, the pseu-dospectral grid on which we evolve Einstein’s equationstakes advantage of the approximate spherical symmetry inthe neighborhood of each star, and in the far-field region.The evolved spatial slice is decomposed into two small ballsaround the center of each neutron star, sets of spherical shellsaround each star and in the far-field region, and distortedcubes to connect the three spherically symmetric regions. Thegrid decomposition used in our simulations is pictured inFig. 1. The inner ball is expanded into Zernike polynomi-als [62, 63], the shells into Chebyshev polynomials (in radius)and spherical harmonics (in angle), and the distorted cubes inChebyshev polynomials. The grid follows the centers of theneutron stars, defined as the center of mass of the matter in the x < and x > half planes, through a simple rotation andscaling of the grid coordinates. We maintain this grid decomposition for the evolution ofEinstein’s equations up to the point at which the maximum Here and in the rest of the paper, the center of mass is defined in the coor-dinates of the simulation, and is clearly a gauge-dependent quantity. It is,however, a convenient quantity to use to define the map between the gridand physical coordinates and keep the neutron stars approximately fixed onthe grid.
FIG. 2. Equatorial cut of the pseudospectral grid just after the changeto a grid centered on the post-merger neutron star remnant. The colorscale shows the baryon density (logarithmic scale). density on the grid increases beyond the low-level oscillationsobserved during the inspiral. This rise in the density signifiesthe transition from two well-separated neutron star cores to asingle, more massive object. At that point, we switch to a gridwhich is fully centered on the coordinate center of mass of thesystem. This grid, pictured in Fig. 2, is made of a ball at theorigin of the coordinate system, surrounded by 59 sphericalshells extending to the outer edge of the computational do-main. Both before and after merger, the outer boundary is lo-cated at d , with d the initial separation of the binary. Forthe configurations considered here, the outer boundary is thusat a coordinate radius of M − M ∼ (1400 − .Both before and after merger, the resolution on the pseu-dospectral grid is chosen adaptively, in order to reach a targetaccuracy estimated from the convergence of the spectral ex-pansion of the solution [59, 64]. As in NS-BH mergers, ourstandard choice is to request a relative error of − in most ofthe computational domain, with a smooth transition to − close to the center of each of the neutron stars before the neu-tron stars get into contact (see [59, 64] for more detail on howthis is computed in practice). Larger and lower truncation er-rors can be imposed for convergence tests, chosen in order toconverge faster than the expected second-order convergenceof the finite volume code.The finite volume grid on which we evolve the general rel-ativistic equations of hydrodynamics is very simple. Beforethe two neutron stars get into contact, it is composed of twocubes, each with × cells each and centered on oneneutron star. The lower number of cells in the vertical dimen-sion is due to the assumption of equatorial symmetry beforemerger. The initial extent of the boxes is listed for each sim-ulation in Table I – it is chosen to cover a little more than theoriginal size of the neutron star. In the coordinate system co-moving with the neutron star centers, the neutron stars expandas the binary inspirals. To avoid losing matter to the outerboundary of the finite volume grid, we expand the grid by . every time the flux of matter across the outer boundaryexceeds . M (cid:12) s − . As the inspiral lasts less than
10 ms ,this implies a mass loss well below − M (cid:12) before merger.As the two neutron stars approach each other, the two finite volume boxes will eventually get into contact. During merger,we would like to follow the forming massive neutron star rem-nant, the tidal tails, the accretion disk, and any ejected mate-rial. We switch to a finite volume grid centered on the formingremnant and with − levels of refinement. Each level has, atour standard resolution, × cells, with the finest gridspacing listed in Table I and each coarser level increasing thegrid spacing by a factor of 2. The lower number of cells inthe vertical direction now reflects the fact that the remnant isless extended in that direction, and thus that we do not needthe finest grid to extend as far vertically as horizontally. Dur-ing merger, we no longer assume equatorial symmetry. Thisnumerical resolution, although insufficient to capture the evo-lution of magnetic fields [22, 23], has been shown to be suf-ficient to obtain reasonable accuracy in purely hydrodynamicsimulations [27]. We estimate the errors by performing a sim-ulation with × cells at each refinement level duringmerger ( × during the inspiral). These error estimateswill be discussed along the relevant results in the remainderof this paper.For these first simulations of neutron star mergers with mi-crophysics and a neutrino leakage scheme in SpEC, we onlyuse two levels of refinement during merger. This is far fromideal, as it allows us to follow material only up to slightlymore than
100 km from the central remnant (
50 km in the ver-tical direction). The reason for this choice is simply to keepthe simulations computationally affordable. With two levelsof refinement, the mergers use about cores for − days, or about , CPU-hours. Although not particu-larly large compared to some of the most advanced generalrelativistic magnetohydrodynamics simulations performed todate [22, 23, 27], this reaches the limits of our computationalresources. An important part of this work is to demonstratethe ability of the SpEC code to perform BNS mergers with mi-crophysics, and to run efficiently on a larger number of coresthan used in previous studies (we typically use ∼ − cores for NS-BH mergers). With our relatively small finitevolume grid, we can follow the formation of the post-mergerneutron star remnant. We can also extract the characteristicsof the remnant and of the surrounding accretion disk, observethe gravitational wave signal , and study neutrino effects. Themass of the outflows, on the other hand, is very approximate.For the simulation using neutrino transport, we are how-ever more interested in the properties of the outflows. Andthe composition of the ejected matter varies as it moves awayfrom the disk, due to neutrino absorption. Accordingly, forthat simulation, we choose to use a third level of refinement(i.e. we make the grid twice as large). C. Neutrino treatment
Most simulations presented here are performed using aleakage scheme to capture the cooling and composition evo- The pseudospectral grid on which we evolve Einstein’s equation extends tomuch larger distances. lution of the post-merger remnant due to neutrino emission.Our leakage scheme is described in more detail in [65, 66],and is based on previous work in Newtonian theory [67, 68]and core-collapse supernovae [69]. Effectively, it is a pre-scription for the energy and number of neutrinos leaving agiven point of the remnant as a function of the local proper-ties of the fluid and the optical depth between that point andthe outer boundary of the computational domain (computed asin [26]). The leakage scheme computes two rates of emission:a free-emission rate, which only depends on the local proper-ties of the fluid and is valid in the optically thin regime, andan approximate, steady state diffusion rate which depends onthe optical depth and is accurate within an order of magnitudein the optically thick regime. In regions of moderate opti-cal depth, an interpolation formula between those two emis-sion rates is used by the leakage scheme. For the emissionand absorption rates, we use the values derived by Rosswog& Liebendorfer [67] for the charged-current reactions (whichonly affect electron neutrinos and antineutrinos), p + e − ↔ n + ν e , (6) n + e + ↔ p + ¯ ν e , (7)electron/positron pair creation and annihilation, e + + e − ↔ ν + ¯ ν, (8)and plasmon decay, γ ↔ ν + ¯ ν. (9)We also use the rates of Burrows et al. [70] for nucleon-nucleon Bremsstrahlung N + N ↔ N + N + ν + ¯ ν. (10)Finally, we use the rates of Rosswog & Liebendorfer [67] forthe scattering of neutrinos on protons, neutrons and heavy nu-clei.To determine the errors in the leakage scheme and obtainmore accurate information about neutrino effects, we alsoperform a single simulation using a more advanced neutrinotransport method. We use the moment formalism, in whichthe energy-integrated (gray) energy density and momentumdensity of neutrinos is evolved on the grid [71 ? , 72]. Ananalytical closure (M1) is then used to compute the neutrinopressure tensor and close the system of equations. The detailsof our M1 algorithm can be found in [21]. Beyond providingmore accurate results and information about the spatial distri-bution of neutrinos, the transport scheme also allows us to takeinto account the absorption of neutrinos and antineutrinos inlow-density regions, which strongly affects the compositionof the fluid, and particularly of the disk and ejecta. These ef-fects were shown to be significant in black hole-neutron starmergers [21]. By explicitly comparing leakage and transportresults, we will show here that using a transport scheme is alsocritical for extracting the composition of the ejecta of binaryneutron star mergers. So far, only a few simulations have useda similar transport scheme for neutron star mergers [17, 20].The transport scheme used in [17, 20] could, however, dif-fer significantly from our algorithm in optically thick regions. Indeed, it relies partially on rates from the leakage schemeto compute the evolution of the component of the neutrinostress-energy tensor which is not trapped by the fluid. We in-stead evolve both the trapped and non-trapped components ofthe neutrino stress-energy tensor together within the momentformalism (but rely on the leakage scheme to estimate the av-erage energy of neutrinos in optically thin regions [21]).For both neutrino schemes, we consider three separate neu-trino species: the electron neutrinos ν e , electron antineutri-nos ¯ ν e , and a third species regrouping all heavy-lepton neutri-nos and antineutrinos ( ν µ , ¯ ν µ , ν τ , ¯ ν τ ), which we denote ν x .Grouping all heavy-lepton neutrinos and antineutrinos in onespecies is justified in our simulations because the temperaturesreached in BNS mergers are low enough to suppress the for-mation of the corresponding heavy leptons. Thus, we do nothave to take into account the charged-current reactions withmuon and tau leptons which would require us to differentiatebetween the four species regrouped in ν x .Finally, we note that the computation of the average en-ergy of the neutrinos differ from what was presented in ourprevious simulations [21, 65, 66]. This is due to the discov-ery of an error in the computation of the average energy ofneutrinos coming from optically thick regions in our leakagecode, and in the GR1D code upon which it is based. This er-ror leads to an overestimate of the neutrino energies. Interest-ingly, because the post-merger remnants in black hole-neutronstar and binary neutron star systems are in a quasi-equilibriumstate in which the neutrino luminosities and electron fractionof the fluid are adapting to the evolving condition in the rem-nant (density, temperature) on a timescale shorter than the dy-namical timescale of the remnant, we find that correcting thiserror does not significantly affect the composition evolution,the neutrino luminosities, or the dynamics of the remnant. Itdoes not even noticeably affect the absorption rate of neutrinosin optically thin regions when coupling our neutrino leakageand transport schemes, because the computation of the neu-trino temperature used to determine the absorption rate in thetransport scheme was previously corrected in optically thickregions to match known test results in spherical symmetry (asdescribed in [21]). This error in the neutrino energies wasdiscovered after the completion of the simulations presentedhere. We have performed significant portions of the simula-tions anew ( ∼ ), to verify that, as stated, our results arenot significantly modified by the use of an improved estimateof the neutrino energies. The neutrino energies quoted in thiswork use the improved estimate of the neutrino spectrum inoptically thick regions, and are typically ∼ lower thanbefore corrections. Updated estimates of the neutrino ener-gies in the systems considered in previous publications usingthe SpEC code [21, 65, 66] will be provided separately. III. INITIAL CONDITIONS AND EQUATIONS OF STATE
We obtain quasi-equilibrium initial conditions for neutronstar binaries by solving the constraints in Einstein’s equationsusing the spectral elliptic solver Spells [73]. The algorithmused to generate initial data for BNS is similar to that previ-
FIG. 3. Mass radius relationships for the three equations of state con-sidered here, for cold neutron stars in neutrino-less β -equilibrium.Horizontal dashed lines are provided at the initial mass of the indi-vidual neutron stars ( . M (cid:12) ) and total mass of the binaries ( . M (cid:12) )considered here. ously developed for NS-BH binaries [74], and spinning neu-tron star binaries [40], and will be discussed in more detail inan upcoming paper (Haas et al., in prep). We consider threedifferent physical configurations. Their main properties arelisted in Table I. In each case, the two neutron stars have equalgravitational masses in isolation, M NS = 1 . M (cid:12) M NS > ∼ M (cid:12) , as required by astronomicalobservations [80, 81]. The mass-radius relationships for theseequations of state are provided on Fig. 3 for cold, neutrino-less β -equilibrium neutron stars. The radii seen on Fig. 3for . M (cid:12) neutron stars are within the range deemed mostlikely by studies incorporating information about both nuclearphysics and astrophysical observations of neutron stars in qui-escent X-ray binaries [82, 83]. A different analysis of thesame astrophysical data however predicts significantly morecompact neutron stars [84], and the physical compactness TABLE I. Initial conditions and grid settings for the binary neu-tron star mergers studied here. “EoS” is the equation of state of theneutron star matter, R . M (cid:12) NS the radius of a neutron star of gravi-tational mass M NS = 1 . M (cid:12) , M . M (cid:12) b the baryonic mass of thesame neutron star, κ . M (cid:12) is the tidal coupling constant (see textand e.g. [85]) a , d the initial coordinate separation, ∆ x t =0FD the initialcoordinate grid spacing on the finite volume grid, and ∆ x mergerFD thegrid spacing in the finite volume grid used after the formation of asupermassive neutron star.EoS R . M (cid:12) NS M . M (cid:12) b κ . M (cid:12) d ∆ x t =0FD ∆ x mergerFD LS220 12.8 km 1.309 271 35.2 km 252 m 306 mLS220 12.8 km 1.309 271 35.2 km 198 m 234 mDD2 13.2 km 1.295 307 36.6 km 277 m 322 mSFHo 12.0 km 1.303 164 38.3 km 249 m 276 m a The authors thank Jan Steinhoff for providing the values of κ for theseequations of state. of neutron stars remains an important open question today.Before merger, the binaries considered here go through 2.5(LS220), 3 (DD2) and 4.5 (SFHo) orbits.The stiffest equation of state, DD2, has a maximum massabove . M (cid:12) even for cold, non-rotating neutron stars. Themerger of two . M (cid:12) neutron stars will thus result in a sta-ble remnant. For the other two equations of state, the mergerremnant will eventually collapse into a black hole. We willsee that a massive neutron star remnant does, however, sur-vive for at least
10 ms after the merger. Note that the maxi-mum baryon mass of a cold, uniformly rotating neutron star(as determined from mass-shedding sequences generated bythe code of Cook, Shapiro, and Teukolsky [86, 87]) is . M (cid:12) for the LS220 equation of state, . M (cid:12) for the SFHo equa-tion of state, and . M (cid:12) for the DD2 equation of state.The total baryon masses of the systems under considerationare . M (cid:12) (LS220), . M (cid:12) (SFHo), and . M (cid:12) (DD2).Hence, for the LS220 and SFHo equations of state, the rele-vant timescale for the remnant to collapse to a black hole isthe pulsar spin-down timescale, which is much longer thanthe timescales relevant to the study of BNS mergers or of theimmediate post-merger remnant evolution (see also [88] for amore detailed discussion of the fate of the post-merger rem-nant in BNS systems). IV. GRAVITATIONAL WAVES
We extract the gravitational wave signal from the sim-ulations following the method presented by Boyle &Mroue [89]: the gravitational waves are extracted on coordi-nate spheres at 24 radii equally spaced in r − within the range [140 , . We then extrapolate the signal to infinity byfitting the finite-radius data with a second-order polynomialin r − . For the extrapolation, we use finite radius values atthe same retarded time (as defined in [89]), to account for thefinite propagation speed of gravitational waves. The extrap-olation error can be estimated by comparing fits of differentpolynomial order to the finite radius data. FIG. 4. Dominant (2,2) mode of the gravitational wave strain for allthree configurations, at
100 Mpc from the source. The simulationwith the stiffest equation of state (DD2) shows the fastest inspiral.As opposed to other figures in this paper, we show the waveformsextrapolated to infinity, which explains the short length of the post-merger waveform.
A. Inspiral and tidal effects
The purpose of this work is to study the merger and post-merger dynamics of neutron star binaries. The simulationsperformed here are generally too short to perform detailedstudies of the gravitational waveform during inspiral. Wewill focus more on the post-merger signal in the next section.Qualitatively, however, we see that the late inspiral proceedsas expected. Figure 4 shows the dominant mode of the gravita-tional wave signal, approximately matched in time and phaseat the beginning of the LS220 simulation. We see that neu-tron stars with larger radii (LS220, DD2) inspiral faster thanneutron stars with small radii (SDHo), and merge at a lowerfrequency. The first effect is due to stronger tidal dissipationfor neutron stars of larger radii. Tides cause a phase differ-ence in the gravitational wave signal proportional to the tidalcoupling constant κ (see Table II), which is strongly corre-lated with the radius of the neutron star. The tidal couplingconstant is, for an equal mass binary, κ = k / (8 C ) with C NS = GM NS / ( R NS c ) the compactness of an isolated neu-tron star, and k the dimensionless Love number. The secondeffect is due to the fact that the merger frequency depends onthe size of the individual neutron stars, with more compactneutron stars merging at higher frequencies. B. Post-Merger gravitational wave signal
The post-merger signal provides more useful information.Previous studies have shown that the gravitational wave spec-trum of the post-merger remnant of a neutron star binaryshows clear peaks at frequencies dependent on the equation
FIG. 5. Power spectrum of the merger and post-merger gravitationalwave signal for optimally oriented mergers at a distance of
100 Mpc ,multiplied by f / . For reference, we also plot the design sensitivityof advanced LIGO (Zero-Detuned High Power detector strain noisespectrum, dashed blue curve) [96]. Note the clear dominant peaks at (2 . −
3) kHz . All spectra are Fourier transforms of the dominant(2,2) mode of the merger waveform in the time interval [ t peak − , t peak + 6 ms] , with a tapering window of width . used atthe beginning and end of the interval. of state of the neutron star [29, 31, 32, 36, 41, 53, 90–95].The strongest of those peaks ( f peak ) occurs at the frequency ofthe fundamental quadrupole mode of the remnant. Recently,Bauswein & Stergioulas [42] have associated the two largestpotential secondary peaks with rotating spiral structures in theremnant ( f spiral ), and a coupling between the fundamentalquadrupole mode and quasi-radial oscillations in the remnant( f − ). They also provide a fitting formula for the location ofthose peaks, based on simulations using an approximate treat-ment of gravity. These simulations assume that the metricis conformally flat, an assumption which can accommodateexactly non-spinning, isolated black holes and neutron stars,but not spinning compact objects or binary systems. In recentsimulation of BNS mergers for . M (cid:12) neutron stars, with ageneral relativistic code and nuclear-theory based equations ofstate, Palenzuela et al. [27] find deviations of ∼ . − . from the fits for f peak and f spiral provided in [42]. A minorpeak is sometimes also found close to the predicted location of f − , but other features of similar amplitudes are also foundat other frequencies.The spectra of the (2 , mode of the gravitational wave sig-nals observed in our simulations are shown in Fig. 5. For allequations of state, emission at the fundamental mode is clearlyvisible. A number of secondary peaks are also observed. Thepost-merger signals are weak when compared, for example, tothe design sensitivity of the Advanced LIGO detector: a veryclose event is required in order to detect the post-merger sig-nal [97]. In Table II, the frequencies of the fundamental andstrongest secondary peaks are compared with theoretical pre-dictions. For the fundamental mode we use the fitting formularecommended by Bauswein et al. 2012 [53], f fitpeak = ( − . R . + 6 . (cid:114) . . f peak < .
64 kHz]= ( − . R . + 8 . (cid:114) . . f peak > .
64 kHz] , with R . the radius of a neutron star of mass M NS = 1 . M (cid:12) [in kilometers], f peak given in kHz, and we assume a totalmass of . M (cid:12) . Here, we used the proportionality relation f peak ∝ (cid:112) M tot /R [53], with M tot the total mass of thebinary at infinite separation, and R max the radius of a neu-tron star at the maximum mass M max allowed for this equa-tion of state. We find disagreements of only ∼ (30 −
70) Hz between the numerical results and the fitting formula, whichwould translate to systematic errors of ∼ (100 − in theradius of a . M (cid:12) neutron star. We should also note that themerger of two . M (cid:12) neutron stars with the LS220 equationof state was studied with an approximate treatment of grav-ity in [53]. The dominant frequency of the post-merger signalwas .
55 kHz in that study. We find an extremely close valuefor that dominant frequency, .
56 kHz .For the secondary peak, we compare our results to the lin-ear fits to f spiral and f − provided in Fig. 4 of Bauswein &Stergoulias [42]. Bauswein & Stergoulias predict that the spi-ral mode should be the dominant secondary mode for mergersin which a stable or long-lived hypermassive neutron star isformed. This is in good agreement with our results as the firstsubdominant peak observed in our spectra agrees well withthe frequency of f spiral . The SFHo equation of state wave-form also shows a peak close to f − and the LS220 equationof state waveform has some extra power at f − . For the DD2equation of state waveform, the predicted location of the f − peak is too close to the merger frequency for any clear featureto be observed. We should note, however, that the predic-tions provided in [42] for f spiral and f − are not as powerfulas the unique fitting formula provided for f peak . This is be-cause a different linear relation between the frequency of themode and the compactness of the star has to be determinedfor each choice of neutron star masses. A universal relationbetween the secondary peak of the post-merger waveform andthe neutron star compactness has been proposed by Takami etal. [31, 32]. This universal relation would predict that the sec-ondary peak is at ∼ . for the LS220 and DD2 equationsof state and at ∼ . for the SFHo equation of state. TheSFHo has its third strongest peak at that frequency, and thatpeak is not much weaker than the secondary peak. The othertwo equations of state show a difference of (100 − between the theoretical predictions and the location of the sec-ondary peak, which would translate into (0 . − .
0) km errorsin the determination of the neutron star radius. Using the uni-versal formula from [31, 32] to infer the compactness of theneutron stars considered in this paper would thus lead to sig-nificant errors in the determination of neutron star radii. Simi-lar differences were observed by Bauswein & Stergoulias [42]for low-mass binaries. Like Takami et al. [31, 32] (but as op-posed to Bauswein & Stergoulias [42]), our code is fully gen-eral relativistic. The observed differences are thus not dueto the treatment of gravity. A more likely explanation is the use of a lower mass system combined with the use of nuclear-theory based equations of state. The choice of equation ofstate may also explain why Takami et al. [31] find significantdifferences between the frequency of the fundamental modeand f fitpeak , while we do not. The largest differences in [31]were observed for the simple Γ = 2 polytrope, while morerealistic equations of state performed better. The general rel-ativistic simulations of Palenzuela et al. [27], performed forhigher mass systems, found ∼ disagreement betweenthe numerical results and f fitpeak – a larger difference than inour simulations, but one that would still allow the recovery ofthe neutron star radius with systematic errors (cid:28) . . Itis quite likely that the fitted frequency f fitpeak is not universal,but nonetheless practically applicable to realistic neutron starequations of state.More insight can be gained in those post-merger featuresby considering a spectrogram of the gravitational wave signal,shown in Fig. 6. The quantity plotted there is ˜ h ( t , f ) = | (cid:82) h , ( t ) W ( t − t ) e ift dt | (cid:82) W ( t − t ) dt , (11)where we choose for the window function W the exact Black-man window, W ( t ) = 793818608 − (cid:18) π (cid:18) t ∆ t (cid:19)(cid:19) + 143018608 cos (cid:18) π (cid:18) t ∆ t (cid:19)(cid:19) (12)for | t | < ∆ t , and W ( t ) = 0 otherwise. The choice of thewindow size ∆ t is a trade-off between high time resolution(small ∆ t ) and high spectral resolution (large ∆ t ). For Fig. 6,we use ∆ t = 6 ms , which causes a noticeable smoothing ofthe spectrogram in time but is necessary to start resolving thesecondary peaks of emission.In the spectrograms, the fundamental peak is clearly visible,and strongest for the softest equation of state (SFHo). It isonly mildly damped over the short duration of the simulations,but varies in frequency as the structure of the remnant evolves.This shift can be as large as
200 Hz in the case of the DD2equation of state, and mostly occurs in the first ∼ aftermerger. Fig. 6 also clearly shows that the gravitational waveemission at the secondary peaks is extremely short-lived, witha decay timescale of ∼ − . The emission at thesecondary peak is both weaker and shorter-lived for the softestequation of state (SFHo). For low-mass systems, we thus seethat only the fundamental mode remains significantly excitedin the post-merger remnant. Gravitational wave emission atlower frequencies is largely coincident with the merger itself,and the secondary peaks are naturally broad in spectral spaceas the signal decays over only a few oscillation periods. Forthe low-mass systems studied here, there is thus a significantdifference between the strong, long-lived peak correspondingto the fundamental mode, and the weak, broad and short-livedpeaks at lower frequencies, which would naturally make thelatter difficult to observe or disentangle from detector noise.The interpretation of the emission of gravitational wavesat frequency f spiral as the result of rotating spiral structures FIG. 6. Spectrograms of the gravitational wave signal for the simulations using the SFHo (left), LS220 (middle) and DD2 (right) equations ofstate, for an optimally oriented binary at
100 Mpc . The dashed horizontal black curves show the location of the peaks predicted in [42]. Thetime is calculated as the difference between the center of the window function used to compute the spectrogram (see text) and the peak of thegravitational wave amplitude. The dominant peak is clearly visible at (2 . −
3) kHz . Secondary peaks in the (1 . − .
1) kHz range are poorlyresolved, and quickly damped. The strong emission at frequencies f < ∼ . is the merger signal itself.TABLE II. Frequency of the two strongest peaks in the spectrum ofthe post-merger gravitational wave signal, f , and f . The predictionfrom [42] for the dominant peaks, obtained from simulations usingan approximate treatment of gravity, are listed as f fitpeak , f fitspiral and f fit2 − . The prediction from [85] for f is f Bpeak .EoS f f f fitpeak f fitspiral f fit2 − f Bpeak
SFHo 2.96kHz 2.1kHz 3.03kHz 2.1kHz 1.7kHz 2.66kHzLS220 2.56kHz 1.8kHz 2.59kHz 1.9kHz 1.4kHz 2.27kHzDD2 2.35kHz 1.7kHz 2.39kHz 1.8kHz 1.2kHz 2.18kHz within the core of the remnant is partially supported by visu-alizations of the rest mass density in the equatorial plane (seeFig. 8). Right after merger, when the emission at f spiral is thestrongest, the LS220 and DD2 simulations show clear spiralarms, including in high-density regions. In the SFHo simula-tion, the spiral arms are significantly weaker. Later on, after merger, lower density spiral structures remain visible inthe LS220 and DD2 simulations, and are stronger in the lat-ter simulation. Again, this is in agreement with the observeddifference in the damping timescale of the emission at f spiral observed in the spectrograms (Fig. 6). Finally,
10 ms aftermerger, extended spiral structures remain visible, but they areconfined to low-density regions and are unlikely to signifi-cantly contribute to gravitational wave emission. Our resultsthus appear consistent with the interpretation of f spiral pro-posed by Bauswein & Stergoulias [42].An alternative universal relation, this time between the strongest peak of the post-merger gravitational wave signaland the tidal coupling constant κ , has been proposed byBernuzzi et al. [85]. Bernuzzi et al. [85] predict f Bpeak = 4 .
341 1 + 0 . κ . κ kHz . (13)We compare f Bpeak to our results in Table II. f Bpeak is system-atically (0 . − .
3) kHz below the peak frequency measuredin our simulations. The magnitude of the error is compara-ble to the scatter observed within the simulations presentedin Bernuzzi et al. [85]. It is unclear whether the system-atic underestimate of the peak frequency comes from our useof a hot nuclear theory based equation of state, the applica-tion of (13) to a low mass system, or the intrinsic scatteraround (13). Bernuzzi et al. use piecewise polytropic equa-tions of state with a Γ -law thermal component for the pres-sure, and f Bpeak is only fitted to mergers with total mass withinthe interval [2 . M (cid:12) , . M (cid:12) ] , so that comparison with ourresults requires extrapolation of their fitting formula to lowermass systems.Comparing the standard and high resolution simulations forthe LS220 equation of state does not reveal any significantresolution dependence in the spectrum of the gravitationalwave signal. Fig. 5 shows mild differences in the shape ofthe secondary peaks, which are however also observed whenchanging the exact interval over which we perform the Fouriertransform, and excellent agreement in the location and am-plitude of the primary peak. Even in the time domain, thelow-resolution and high-resolution waveforms appear to only0 FIG. 7. Dominant (2,2) mode of the gravitational wave strain forthe LS220 waveform at the standard and high resolution, at
100 Mpc from the source. We find excellent agreement even in the post-mergersignal. differ by a small time shift (see Fig. 7). We note that forboth Fig. 6 and Fig. 7, we use waveforms extracted at finiteradius ( r ∼ M ) to have access to a longer post-mergersignal. However, Fig. 6 shows that only a few millisecondsof post-merger signal are actually necessary to study the sec-ondary peaks in the post-merger spectrum, since the gravita-tional wave emission at those frequencies decays on a veryshort timescale. V. MERGER AND POST-MERGER REMNANT
We now consider the qualitative properties of the mergerand post-merger remnants. Figure 8 shows snapshot of thedensity profile in the equatorial plane of the binary . , and
10 ms after the peak of the gravitational wave am-plitude for all 3 equations of state, which we define as the timeof merger. Initially, the properties of the system are largely de-termined by the compactness of the pre-merger neutron stars.With the SFHo equation of state, i.e. for the most compactneutron stars, a compact core forms rapidly. For the otherequations of state (LS220, DD2), more strongly developedtidal features appear at merger. . after merger, we stillobserve two well-defined cores whose size and tidal distor-tion directly correlates to the pre-merger size of the individualneutron stars.After , the cores start to merge. The less compact neu-tron stars (LS220, DD2) produce clearly defined, high den-sity tidal tails. Finally,
10 ms after merger, the propertiesof the remnant are dominated by the high-density behaviorof the equation of state. For high neutron star masses, theLS220 equation of state is nearly as compact as the SFHoequation of state - and neither of them is able to support a sta-ble, non-rotating, cold . M (cid:12) neutron star. The LS220 and SFHo equations of state now have similarly compact cores,with central density rising slowly as the neutron star evolves.Since both post-merger remnants have baryon masses well be-low the maximum mass of a uniformly rotating cold neutronstar (at the mass-shedding limit), we expect the post-mergerneutron star to survive for much longer than the duration ofthe simulation. Collapse to a black hole should occur on thetimescale necessary for the post-merger neutron star to looseangular momentum through magnetically-driven spin-down.On the other hand, the DD2 equation of state forms a massiveneutron star which will remain stable even in the non-rotatinglimit. We see that the post-merger neutron star evolves onlyslowly at late times (see Fig. 9). All three configurations re-main far from axisymmetry. Strong spiral waves are drivenin the forming accretion disk, and a rapidly rotating bar moderemains present in the post-merger neutron star (see Fig. 10).Figure 9 also shows that for all three configurations, the peakof the gravitational wave emission is coincident with the corebounce of the merging neutron stars, i.e. the time at which thetwo cores touch. The core bounce causes a sharp increase inthe density and pressure of the cores, as well as strong shocks,shears, and heating at the interface between the merging neu-tron stars.The excitation of the post-merger remnant, shown inFig. 10, is consistent with the fundamental quadrupole mode,both in terms of the density and velocity patterns, and in termsof the measured gravitational wave frequency. Because ofthis excitation, the spatial variation of the instantaneous or-bital angular velocity, defined in the coordinates of the sim-ulation and with respect to the center of mass of the orig-inal binary, is fairly complex. We show the instantaneousorbital frequency of the fluid in Fig. 11. Although on aver-age the rotation frequency is higher at smaller radii, we ob-serve regions in which d Ω /dr > . There are thus regionsof the post-merger neutron star which may not be subject tothe axisymmetric magnetorotational instability (MRI) [98], atleast shortly after merger. However, all regions with densi-ties ρ < ∼ g cm − follow a typical disk profile for theorbital frequency, and should see rapid growth of the mag-netic field. The typical timescale for the growth of the ax-isymmetric MRI is τ MRI ∼ / (3Ω) ∼ [98]. Magneticfield should thus grow significantly in the disk over even theshort timescale of our simulation. So far, simulations haveonly been able to resolve the growth of the MRI in lower-density regions ( < ∼ g / cm ) [22], and indeed find rapidgrowth of the magnetic field at early times due to a combi-nation of magnetic instabilities and winding of the magneticfield. We should also note that the angular velocity in the rem-nant varies significantly when comparing a vertical slice alongthe direction of the bar mode in the remnant (as in Fig. 11),with a slice orthogonal to that direction (see, e.g., the az-imuthal dependence of the velocity on Fig. 10). Figure 11may give the impression that the bar mode is rotating at fre-quencies much lower than the expected f peak / . In fact, wecan check that the pattern speed of the bar mode is consistentwith f peak / . In the LS220 simulation, we measure a rota-tion of the bar of ∼ ◦ over the last . of evolution, or afrequency f bar ∼ .
25 kHz ∼ f peak / . We also find that the1 FIG. 8. Density in the equatorial plane of the post-merger remnant at three characteristic times: . after the peak of the gravitational waveemission (top), after the peak (middle), and
10 ms after the peak (bottom). We show results for the SFHo (left), LS220 (center) and DD2(right) equations of state. In each of those plots, we show density contours at (1 , , , , × g cm − (solid blue lines). Here and insubsequent figures, the coordinates are those of the simulation, as evolved in the generalized harmonic formulation of Einstein’s equations. orbital frequency in the core of the post-merger neutron star ishigher than the mass-shedding limit for a uniformly rotatingneutron star of the same baryon mass. Considering that the to-tal baryon mass of the system is below the maximum mass ofa cold, uniformly rotating neutron star (at the mass-sheddinglimit), this suggests that rotational support is initially suffi-cient to prevent the collapse of the remnant to a black hole.This should be true even if the differential rotation is rapidlyerased by angular momentum transport within the remnant.The properties of the post-merger remnant are strongly af-fected by the presence of a non-linear l = 2 perturbation. Vis-ible shock fronts remain in the disk at the end of the simu-lation, while a clear bar mode is observed in the post-merger neutron star. In Fig. 12 and Fig. 13, we plot the density, tem-perature, composition and entropy of the fluid in the equa-torial plane and in a meridional slice along the major axisof the bar observed in the post-merger neutron star. For allthree simulations, the neutron star remains neutron rich and,although hotter than the disk in terms of absolute temperature,with T ∼ (15 −
25) MeV , remains entirely supported by nu-clear forces and rotation [88]. The entropy per baryon in theneutron star is only S ∼ k B . The neutron star is thus colderthan the −
50 MeV temperatures observed in higher masssystems [29], presumably due to the fact that lower mass starsare less compact, and thus that shock heating of the neutronstar during the merger is not as important of an effect in the2
FIG. 9. Maximum value of the baryon density on the grid as a func-tion of time for all three equations of state. t peak is the time at whichthe amplitude of the gravitational wave signal is maximalFIG. 10. Density and velocity in the equatorial plane for the sim-ulation with the LS220 equation of state,
10 ms after merger. Theexcited quadrupole mode is still clearly visible at that time. low-mass systems considered here. This interpretation is par-tially confirmed by the fact that the softer equation of state(SFHo) has a core temperature significantly above the stifferequations of state (DD2, LS220). All three simulations how-ever show significant heating at the interface between the barand the disk, with T ∼ (30 −
40) MeV . At the lower densitiesobserved in that region, the thermal pressure now becomesimportant.The maximum density in the accretion disk remains high, ρ disk > ∼ g cm − . The shock front associated with the l = 2 perturbation of the disk heats the disk, causing tempera-tures of (8 −
10) MeV right at the shock. Accordingly, the diskis thick, with a scale height H ∼ R , and thermally supported,with P ( ρ , T, Y e ) (cid:29) P ( ρ , , Y e ) : for ρ = 10 g cm − , P ( ρ , T, Y e ) ≈ P ( ρ , , Y e ) for T = 2 MeV , and most ofthe disk is much hotter than that.We test the post-merger stars and disks for axisymmetricconvective instability using the relativistic Solberg-Hoiland- Ledoux condition ([99] generalized to include Y e gradientterms). This condition assumes an axisymmetric, stationarybackground, which is only approximately present, so we usedensity-weighted azimuthal average profiles. Consider theLS220 case, which is evolved both in leakage and M1. Con-centrating on the region near the equator, where the radial con-dition is most important, we find that the star and disk arestable everywhere except for a small region around cylindri-cal radius of 55 km, where the angular momentum gradientis Rayleigh unstable, an indication that matter at this distancehas not yet acheived equilibrium. The star and bulk of the diskare convectively stable even though the Y e gradient is unsta-ble through much of the star. This effect is counteracted bythe much larger stabilizing shear and entropy gradient terms.We can compare our results with recent simulations usingnuclear-theory based equations of state at similar compactnessand higher mass [26], similar radii and higher mass [20, 27],and simulations approximating the thermal dependence of theequation of state by a Γ -law for a wide range of systems ofhigher mass [29], and a single system with a large neutron starof similar mass [100]. The temperature of the remnant andtidal features are closest to those found by Neilsen et al. [26],confirming the correlation of those features with the compact-ness of the neutron star. More compact neutron stars generallyshow stronger shock heating in the post-merger neutron star,stronger excitation of the bar mode in the core, and weaker l = 2 features in the disk. The first two are due to smallerneutron stars merging at closer separation and thus higher ve-locities, causing a more violent merger event. The third ispresumably due to stronger tidal effects in less compact stars.The main features of the merger in our simulations are in factremarkably similar to those observed in [26, 27], especiallywhen one takes into account the expected dependence of theresults on the compactness of the star. The agreement betweenour results and [26, 27] is reassuring considering the use ofsimilar equations of state and a mostly identical neutrino leak-age scheme.This agreement does not, however, preclude systematic er-rors due to the limited physics included in the simulations.One such limitation is the use of a leakage scheme to treatneutrino cooling, and the absence of any treatment of neu-trino absorption. We can test the effect of these assumptionsby using a neutrino transport scheme. To do this, we evolvedthe merger with the LS220 equation of state using our M1transport scheme, which evolves the energy density and mo-mentum density of the neutrinos [21]. Figures 12-13 showthat the qualitative evolution of the neutron star remnant islargely unaffected by the treatment of the neutrinos. This isnot at all surprising, given that at the large densities exist-ing within the neutron star, the neutrinos are trapped and inequilibrium with the fluid (see also Sec. VI). As already ob-served in the disks resulting from NSBH mergers [21], theinclusion of neutrino transport causes a smoothing of the tem-perature profile. This is particularly visible in the equatorialplane at the bar-disk interface. Neutrino absorption also natu-rally causes cold low-density regions in the disk to be heated:most of the disk is ∼ hotter when using neutrino trans-port. The mild smoothing of the temperature profile, heating3 FIG. 11. Instantaneous orbital frequency of the fluid, in the coordinates of the simulation and
10 ms after merrger, in a slice orthogonalto the orbital plane and along the major axis of the bar in the post-merger neutron star. The black lines are density contours at ρ =(10 , , ) g / cm .FIG. 12. Hydrodynamic variables (in the equatorial plane)
10 ms after the peak of the gravitational wave signal. From top to bottom, we showthe density, temperature, and electron fraction for the four simulations considered in this paper – from left to right, SFHo, LS220 with neutrinoleakage, LS220 with neutrino transport, and DD2. of the outer disk, and shocked tidal arms can also be observedin the one-dimensional density and temperature profiles pre-sented on Fig. 14.These are fairly minor effects, with some impact onthe neutrino luminosity and neutrino-matter interactions (see Sec. VI), but not on the hydrodynamic properties of the post-merger remnant. The treatment of the neutrinos really beginsto matter when one considers the composition of the disk, asmeasured by its electron fraction Y e = n p / ( n p + n n ) , with n p and n n the number density of protons and neutrons respec-4 FIG. 13. Hydrodynamic variables
10 ms after the peak of the gravitational wave signal in a slice along the major axis of the remnant andorthogonal to the equatorial plane. From top to bottom, we show the density, entropy per baryon, temperature, and electron fraction for thefour simulations considered in this paper – from left to right, SFHo, LS220 with neutrino leakage, LS220 with neutrino transport, and DD2.FIG. 14. Density ( left ) and temperature ( right ) profiles in the equa-torial plane along the major axis of the merger remnant. We showresults for the merger using the LS220 equation of state and the neu-trino transport (solid lines) or neutrino leakage (dashed line) scheme.The two simulations are very similar, with a slightly smoother tem-perature profile and hotter outer disk when using neutrino transport.The hot shocked tidal arms are clearly visible around
40 km from thecenter in both plots. tively. In leakage simulations the neutrino emission, and thusthe composition of the disk, is set only by the local propertiesof the fluid and an estimated neutrino optical depth. The trans- port scheme, on the other hand, takes into account the irradi-ation of cold regions of the disk by the hot neutrino-emittingregions. Most of the neutrino emission occurs either at thecore-disk interface or at the tidal shocks in the disk, and thecomposition of the disk is largely set by the relative positionof a fluid element with respect to these two defining features.We will see in Sec. VI that the disk is mostly irradiated byelectron antineutrinos. Accordingly, the regions immediatelyadjacent to the neutrino-emitting regions, i.e. just next to thehot shock front in the disk or the surface of the post-mergerneutron star, are forced towards a very low electron fraction Y e ∼ . − . . Farther away from those emitting regions,we get Y e ∼ . − . . The leakage scheme, which doesnot take into account the absorption of electron antineutrinosin the disk, predicts an electron fraction Y e ∼ . − . inmost of the high-density regions of the disk, with smaller Y e in the corona and at low radii. There is thus a very significantdifference in the composition of the disk between leakage andtransport simulations, due to the fact that the composition ofmost of the disk is strongly affected by neutrino absorption.This can be contrasted with the remnant of NSBH mergers,where only the low-density corona is significantly affected byneutrino absorption [21]5 FIG. 15. Total neutrino luminosity for the three simulations using theleakage scheme and the LS220, SFHo, and DD2 equations of state(solid curves) as well as the simulation using the LS220 equation ofstate and neutrino transport (dashed red curve).
We also note that neutrino absorption in the transport sim-ulation causes the shocked regions in the tidal arms of thedisk to be about . hotter than predicted by the leak-age scheme (see Fig. 14). The shock front is energized byneutrinos emitted in the inner disk or in the core.Despite these effects, our simulations indicate that neutrinoleakage is a reasonable approximation as long as the exactcomposition of the fluid is not required, i.e. except when at-tempting to make predictions for r-process nucleosynthesis inthe dynamical ejecta and disk winds, and for the associatedradioactively powered electromagnetic signal.Finally, we note that for all the variables discussed in thissection, the low and high resolution simulations with theLS220 equation of state give results well below the differ-ences between simulations using different equations of state.We thus expect the results to be only mildly affected by thenumerical error in the simulations. VI. NEUTRINO EMISSION
When studying the neutrino emission in simulations usingthe leakage scheme, we can consider only the total energy andnumber of emitted neutrinos, as well as the predicted loca-tion of emission. We list the luminosity and average energyof neutrinos in all three simulations in Table IV, and providethe time dependence of the total luminosity in Fig. 15. For theaverage energy, we compute the energy-weighted root-mean-square energy, which is the important quantity when estimat-ing the energy deposited in optically thin regions through neu-trino absorption. For a blackbody spectrum, as assumed here,it is a factor of ∼ . larger than the number-weighted aver-age energy of the neutrinos. We find results similar to sim-ulations of higher mass neutron star mergers using leakage FIG. 16. Density (colour scale) and location of the surfaces of opticaldepth τ ∼ for ν e (solid black line), ¯ ν e (dashed red line), and ν x (dot-dashed black line) at the end of the simulation using the LS220equation of state and a neutrino transport scheme. schemes [26, 27] and leakage-transport hybrid [20] schemes.The electron antineutrinos have the highest luminosity, with L ¯ ν e ∼ (2 − × erg s − about
10 ms after merger. Theydominate the electron neutrinos by a factor of . − .To understand the geometry of the neutrino radiation, it isuseful to determine which regions of the remnant are opticallythick to neutrinos, and which are not. To do so, we show onFig. 16 the location of the surface of optical depth τ ∼ for allspecies of neutrinos, as computed by the leakage scheme. Theremnant neutron star is naturally optically thick to all neutri-nos. For the electron neutrinos and antineutrinos, this is alsothe case for most of the inner disk, up to or even further thanthe location of the shocked tidal arms.The electron neutrinos and antineutrinos are emitted fromdifferent regions of the remnant (see also [26]): as will bemade clearer when considering the results of the simulationusing neutrino transport, most electron neutrinos come fromthe polar regions of the core of the remnant, but most elec-tron antineutrinos are emitted in the hot shocked regions ofthe disk. This is largely due to the fact that these shocked re-gions are surrounded by material which has low optical depthto ¯ ν e , but is optically thick to ν e (see Fig. 16). The heavy lep-ton neutrinos, whose emissivity is more sensitive to the tem-perature of the fluid, are mostly emitted in the hottest regionsof the core. Accordingly, the luminosity in ¯ ν e scales with thetemperature of the shocked regions of the disk while the lumi-nosity in ν e and heavy-lepton neutrinos scales with the tem-perature of their respective neutrinosphere in the core (i.e. theregion of optical depth of order unity, which is closer to thesurface for ν e than for the heavy-lepton neutrinos). For allspecies, the highest luminosity is naturally reached in the caseof the most compact neutron stars (SFHo equation of state),given the stronger heating of the remnant occurring in thatcase. The LS220 equation of state, which forms a relativelycold core (see e.g. Fig. 13), shows a particularly large differ-ence between the emission of electron antineutrinos and elec-tron neutrinos.The average energy of the electron antineutrinos is set bythe equilibrium spectrum in the shocked regions of the disk,6 TABLE III. Neutrino luminosities and rms energy of the neutrinosaccording to the leakage scheme, after merger. The ν x valuesare for all heavy-lepton neutrinos and antineutrinos combined.EoS L ν e L ¯ ν e L ν x (cid:112) (cid:104) (cid:15) ν e (cid:105) (cid:112) (cid:104) (cid:15) ν e (cid:105) (cid:112) (cid:104) (cid:15) ν x (cid:105) Units erg s − MeVSFHo 1.9 3.0 2.2 14 21 29LS220 1.2 2.1 1.2 13 20 26DD2 1.6 2.2 0.9 13 20 24 since they are emitted there and then largely free-stream awayfrom the disk. The average energy of the other species is setby the temperature of the core at the point at which the neu-trinos decouple from the fluid, which occurs at higher density(and temperature) for the heavy-lepton neutrinos than for theelectron neutrinos.At the high luminosities observed here, we expect neutrinoabsorption in low-density regions to be important to the evo-lution of the composition of the disk and of the outflows. Thegeometrical properties of the neutrino distribution and the im-pact of the neutrinos on the evolution of the system may notbe suitably captured by the leakage scheme. We estimate theimpact of the use of the approximate leakage scheme by con-sidering the results of the simulation using the neutrino trans-port scheme. In that simulation, we evolve the neutrino energyand momentum density on our numerical grid, and take intoaccount neutrino absorption and scattering in low-density re-gions.Figures 17 and 18 show the energy density E ν of electronneutrinos and antineutrinos on the surface of matter density ρ = 10 g / cm , together with the neutrino ”transport” ve-locity v iT,ν = αF iν /E ν − β i , i.e. the velocity such that in theabsence of source terms ∂ t ( √ gE ν )+ ∂ i ( √ gE ν v iT,ν ) = 0 , andthe energy density E ν is simply transported through the gridat velocity v iT,ν . Here, F ν is the neutrino momentum den-sity. We confirm that most of the emission of ¯ ν e comes fromthe tidal arms in the disk (dark blue regions in Fig. 18), whilemost ν e escaping the system are emitted in the polar regions ofthe massive neutron star remnant (center of Fig. 17). We note,however, that even in the polar regions the electron neutrinosare outnumbered by the electron antineutrinos. Emission fromthe tidal arms is significantly beamed, due to fluid velocities v ∼ . c . The emission coming from the polar regions has awider opening angle.Figures 19 to 21 offer a different view of the same variables,through slices orthogonal to the equatorial plane and parallelto the major axis of the bar in the remnant. The differencesbetween the electron neutrinos, the electron antineutrinos, andthe heavy-lepton neutrinos are clearly visible. From the smallfluxes seen in Fig. 19, we can see that the electron neutrinosare trapped not only in the post-merger neutron star, but alsoin most of the disk – up to the second, weaker tidal shock at r ∼
55 km (also visible on Fig. 13). Electron neutrinos mostlyescape through the polar regions. Most of the electron neutri-nos emitted by the hot central neutron star toward the diskwill be absorbed in the disk, due to disk optical depth τ (cid:29) (see Fig. 16). The optical depth of the fluid to electron an- FIG. 17. Energy density (color scale) and normalized flux αF iν /E ν − β i (black arrows) of the electron neutrinos on the density isosurfacewith ρ = 10 g / cm . The region of brightest ν e emission, at thecenter of the plot, coincides with the polar region of the neutron starremnant.FIG. 18. Energy density (color scale) and normalized flux αF iν /E ν − β i (black arrows) of the electron antineutrinos on the density isosur-face with ρ = 10 g / cm . The regions of brightest ¯ ν e emission(dark blue) are in the hot, shocked tidal arms. tineutrinos is smaller. Only the regions inside of the first tidalshock in the disk are optically thick, as visible from the largeradial fluxes observed at radii r > ∼
40 km in Fig. 20 (but notin Fig. 19). Electron antineutrinos emitted from the shockedregions at r ∼
40 km are nearly free-streaming through theouter disk. Nearly all heavy-lepton neutrinos are emitted inthe hot regions of the core, and are free-streaming in most ofthe disk. Finally, we note that the core of the remnant is moreneutron rich than the equilibrium value at the temperature ofthe core ( Y eq e ∼ . − . at the center of the remnant,where we observe Y e ∼ . − . ). This results on a negativeequilibrium electron neutrino chemical potential, and leads toa suppression of the electron neutrino density compared to theelectron antineutrino density in the core. This is opposite to7 FIG. 19. Energy density and normalized flux αF iν /E ν − β i of theelectron neutrinos in the same vertical slice as in Fig. 13. The verticalstripes with higher neutrino energy in the polar regions are an artifactof the M1 closure.FIG. 20. Energy density and normalized flux αF iν /E ν − β i of theelectron antineutrinos in the same vertical slice as in Fig. 13. Thevertical stripes with higher neutrino energy in the polar regions arean artifact of the M1 closure. the situation found in the cores of protoneutron stars in a core-collapse supernova (see e.g. [101]).Figures 19 and 20 also show a known limitation of the M1transport scheme: its inability to handle crossing beams (seee.g. the code tests presented in [21]). Artificial shock frontsdevelop in the neutrino evolution in the regions in which neu-trinos from the disk and from the core converge. There, neutri-nos coming from different directions start to propagate in theiraverage direction of motion instead of crossing each others.This reduces the opening angle of electron neutrinos and an-tineutrinos leaving through the polar regions, and could plau-sibly affect the composition of matter outflows in those re-gions.From these observations, we can understand better the dif-ferences between the leakage and transport results discussedin Sec. V. In the post-merger neutron star and in the inner partsof the disk, neutrinos are either in equilibrium with the densematter or slowly diffusing through it without much of an effecton its composition over the short timescales considered here.Accordingly, transport and leakage results are very similar. Inthe outer disk and in the polar regions, low-density materialis strongly irradiated by neutrinos. If neutrino irradiation wasthe only driver of composition evolution (i.e. excluding neu-trino emissions and the advection of trapped neutrinos), wewould expect the fluid to have composition (assuming absorp- FIG. 21. Energy density and normalized flux αF iν /E ν − β i of theheavy-lepton neutrinos in the same vertical slice as in Fig. 13.FIG. 22. Expected equilibrium electron fraction of the fluid Y irr e ifthe composition was set solely by neutrino absorption. The whiteregion has Y irr e > . . tion cross-sections proportional to (cid:15) ν ) Y irr e ∼ (cid:104) (cid:15) ν e (cid:105) E ν e (cid:104) (cid:15) ν e (cid:105) E ν e + (cid:104) (cid:15) ¯ ν e (cid:105) E ¯ ν e , (14)a quantity plotted in Fig. 22 for our simulation. Otherwise, theeffect of absorption is at least to push the composition closerto Y irr e than what one would expect when neglecting absorp-tion.The impact of neutrino absorption can be observed by com-paring Fig. 22, which shows the electron fraction towardswhich the fluid is driven by neutrino absorption, and Fig. 13,which shows the actual electron fraction in the simulations(with both the leakage and the transport scheme). We seethat in the disk and right outside of the innermost shockedregion, neutrino absorption drives Y e down (electron antineu-trinos strongly dominate electron neutrinos). The closer to theshocked region a fluid element is, the more neutron-rich it be-comes. In the polar regions, on the other hand, there remainenough electron neutrinos to drive the composition to a higher Y e , at least close to the surface where the irradiation is thestrongest. Finally, outside of the weaker shock at r ∼
55 km (seen on Fig. 13), where the electron neutrinos decouple fromthe disk, the effect of neutrino absorption is also to raise Y e .Finally, we note that there are significant differences be-tween the leakage and transport simulations even in terms ofthe total neutrino luminosity in each species. Figure 23 showsthe total luminosity as a function of time for each species forthe LS220 simulation, either with neutrino leakage or with8 FIG. 23. Neutrino luminosities for the simulations using the LS220equation of state and the neutrino leakage scheme (solid curves), orthe neutrino transport scheme (dashed curves). We show the lumi-nosity in electron neutrinos (black), electron antineutrinos (red), andall heavy-lepton neutrinos and antineutrinos combined (green). neutrino transport. We observe differences of as much as afactor of two. The electron neutrino luminosity, which wasalready particularly low for the LS220 equation of state whenusing the leakage scheme (see Table IV), is a factor of foursmaller than the electron antineutrino luminosity when usingthe transport scheme. This can be explained by the fact that alarge fraction of the neutrinos, which the leakage scheme pre-dicts are escaping from the core, are emitted in the directionof the accretion disk. The disk is optically thick to electronneutrinos, and when using the transport scheme those neutri-nos are reabsorbed in the disk. Our leakage scheme does notaccount for this effect, because it allows neutrinos to movealong the path of smallest optical depth (in this case, alongthe polar axis), instead of following null geodesics when free-streaming.We also note that, even though the electron antineutrino lu-minosity appears to agree well in the leakage and transportsimulations, this is a mere coincidence. There are indeedtwo potential sources of disagreement between the leakageand transport results. The first is the error in the estimates ofthe leakage scheme for given physical conditions in the post-merger remnant. The second is that the physical properties ofthe post-merger remnant themselves evolve differently in theleakage and transport simulations. One important such differ-ence is that the disk is hotter in the simulation using neutrinotransport, due to neutrino absorption. So the predicted emis-sion rate of ¯ ν e is larger in the transport simulation. But theleakage scheme neglects neutrino absorption, which causesit to overpredict the total luminosity of ¯ ν e (albeit not by asmuch as for ν e ). In this specific case, the two effects appear tolargely cancel each other. But one can check that they are notnegligible by using the leakage scheme to predict the neutrinoluminosity at a given time in the simulation using a transport TABLE IV. Properties of the outflows measured within the first
10 ms following the merger. The three simulations using the LS220equation of state are the leakage simulation at our standard resolution(LS220-L0), the leakage simulation at high resolution (LS220-L1),and the simulation using the two-moment transport scheme (LS220-M1). M ejecta is the amount of mass which is flagged as unbound( hu t < − ) when leaving the computational grid, (cid:104) S (cid:105) and (cid:104) Y e (cid:105) aredensity-weighted averages of the entropy and electron fraction, and M polar is the mass of unbound material escaping through the upperand lower boundaries of the computational domain.Name M ejecta (cid:104) S (cid:105) (cid:104) Y e (cid:105) M polar /M ejecta Units − M (cid:12) k B baryon − SFHo 5 10 0.11 0.13LS220-L0 2 12 0.11 0.01LS220-L1 13 10 0.10 0.07LS220-M1 5 21 0.20 0.56DD2 13 11 0.11 0.20 scheme. The resulting ¯ ν e luminosity is a factor of two largerthan when using the transport scheme.In conclusion, it appears that using a neutrino transportscheme is necessary for reliable predictions regarding thecomposition of the outflows, the composition of the disk, orthe neutrino luminosity. The latter is particularly importantin the post-merger evolution of BNS mergers because of thecomplex geometry of neutrino emission, and the impact ofneutrino absorption on the evolution of the temperature andcomposition of the disk. VII. OUTFLOWS
As discussed in Sec. II, the simulations performed here usetoo small a numerical grid to accurately measure the mass ofunbound material leaving the system. For reference, however,we provide in Table IV the main properties of the materialwhich satisfies the approximate unbound condition hu t < − as it crosses the outer boundary of our grid, with u µ the 4-velocity of the fluid. The quantity hu t is conserved along afluid line if ∂ t is a Killing vector (see e.g. [102]), which weassume to be approximately true far from the remnant and af-ter merger. In Table IV and the rest of this section, all averagesrefer to density weighted averages, i.e. (cid:104) X (cid:105) = (cid:82) ρ √− gu t XdV (cid:82) ρ √− gu t dV , (15)with g the determinant of g µν and dV the flat-space volume el-ement. We note that for the total mass ejected the low and highresolution simulations with the LS220 equation of state givevery different results. The difference is due to variations in thesmall amount of mass which is ejected at the time of merger.This is the only quantity in our simulations for which the tworesolutions disagree, and a clear indicator that the mass of thedynamical ejecta is unreliable.The other notable result is that in all leakage simulationsthe ejected material has very similar properties: average spe-9 FIG. 24. Rate at which material flagged as unbound leaves the com-putational grid in the simulation using the LS220 equation of stateand the neutrino leakage scheme. We also show the flow of unboundmaterial at the same location in the simulation using the transportscheme. Note that due to neutrino absorptions at larger distances(the transport simulation uses a larger numerical grid), the true massoutflows in the simulation using a transport scheme are higher at latetimes, with ˙ M ∼ . M (cid:12) s − . The oscillations in the outflow ratesare largely due to the use of a rectangular outer boundary and theasymmetry of the outflows. cific entropy (cid:104) s (cid:105) ∼ k B baryon − , average electron frac-tion (cid:104) Y e (cid:105) ∼ . , matter ejection occurring mostly close to theequatorial plane, and no outflows at late times in the simu-lations. By comparison, the transport simulation ejects ma-terial with a higher electron fraction, twice the average en-tropy, and shows the ejection of material in the polar regions atlate times in the evolution (see Fig. 24). This neutrino-drivenwind has ˙ M ∼ . M (cid:12) s − . If maintained for −
100 ms ,it could be the dominant source of ejecta in the merger (seealso [103]). The wind is generally more proton rich andhotter than the earlier dynamical ejecta, with (cid:104) Y e (cid:105) ∼ . and (cid:104) s (cid:105) ∼ k B baryon − . In the M1 simulation, we find (cid:104) Y e (cid:105) ∼ . and (cid:104) s (cid:105) ∼ k B baryon − for the dynami-cal ejecta (material ejected in the first ∼ following themerger).The composition and entropy of the ejecta is important fortwo main reasons. The first is that neutron-rich material un-dergoes strong r-process nucleosynthesis, producing mostlyheavy elements with number of nucleons A > , whileless neutron-rich material will produce a larger fraction oflower mass elements with < A < [17]. The late-time wind observed here has an electron fraction which is atthe limit between the two regimes. Lippuner & Roberts [15]predict that, for the entropy seen in our simulations, the lim-iting composition is Y e ∼ . . This suggests the productionof a broad range of elements in the neutrino-driven wind, inagreement with what has been observed in higher mass neu-tron star mergers [17], and in the evolution of post-merger ac-cretion disks [45]. The second consequence is that the heavy- elements produced in neutron-rich material, particularly thelanthanides, have very high photon opacities. Their presencecauses the electromagnetic transients powered by r-processnucleosynthesis to peak on a timescale of a week instead ofa few days, and in the infrared instead of the optical [14]. Ifthere existed a lanthanide free region in the outflows, emissionfrom that region could power an earlier, optical signal [104].However, the conditions observed in our simulations do notfavor such a scenario.Finally, the strong wind generated in our simulation withneutrino transport also causes significant baryon loading ofthe polar regions, with a measured wind of ˙ M ∼ . M (cid:12) s − during the last of the simulation. This could hamperthe production of short-gamma ray bursts, at least until thecollapse of the post-merger neutron star to a black hole, whichwill not happen for the DD2 equation of state since the post-merger neutron star is stable in that case. VIII. CONCLUSIONS
We present a first set of simulations of neutron star merg-ers using the SpEC code, nuclear-theory based equations ofstate, and either a simple neutrino cooling prescription (leak-age) or a two-moment grey neutrino transport scheme. Forthis first study, we focus on equal mass neutron stars with M NS = 1 . M (cid:12) , at the low end of the expected range of neu-tron star masses. So far, only a few studies with fully gen-eral relativistic codes have included the effects of neutrinoswith leakage [26, 27] or transport schemes [20]. All have fo-cused solely on equal mass systems with M NS = 1 . M (cid:12) .These simulations allow us to study the properties of the post-merger gravitational waveform for BNS mergers with realisticequations of state, the qualitative properties of the post-mergerremnant, the effects of neutrinos on the post-merger evolution,and the impact of the choice of neutrino treatment on the re-sults of the simulations.The post-merger gravitational wave signal is known to bedominated by strong peaks at frequencies dependent on theneutron star equation of state, although there is some disagree-ment on the interpretation of the peaks and predicted emissionfrequency. The strongest post-merger peak, which is associ-ated with the fundamental l = 2 excitation of the post-mergerneutron star, largely dominates the post-merger signal in oursimulations. We find that the location of that peak is in goodagreement (to better than ∼
100 Hz ) with the predictions ofBauswein et al. [53], which were based on a large number ofsimulations using an approximate treatment of gravity. Othergeneral relativistic simulations have recently found large de-viations from those predictions [31], but for choices of equa-tion of state incompatible with the expected properties of neu-tron stars. General relativistic simulations with nuclear-theorybased equations of state of higher mass neutron star merg-ers [27] found agreement at the level ( [200 − differences) with the fitting formula provided by Bauswein etal. [53]. Although the accuracy of that fitting formula remainsan open question, our simulations tend to confirm that theyperform well for nuclear-theory based equations of state.0We are less optimistic when it comes to the determinationof the secondary peaks of the post-merger signal in low-masssystems. Although we do observe peaks at frequencies closeto those recently predicted by Bauswein & Stergoulias [42],and corroborated by the simulations of Palenzuela et al. [27],we find that their emission is limited to a very short time pe-riod of ∼ (1 −
3) ms , which causes the peaks to be very broadin frequency space. Their detectability when mixed with real-istic detector noise may be debatable.The properties of the post-merger remnant appear to fol-low the trends observed at higher mass when considering low-compactness neutron stars. The post-merger neutron star isnot as strongly heated as in higher mass, more compact sys-tems, but strong l = 2 modes survive for the entire durationof the simulation, up to
10 ms past merger. When treating theneutrinos through the simple leakage scheme, we find goodagreement in the basic properties of the remnant (density, tem-perature, composition) and the observed neutrino luminositywith recent results for higher mass systems with a similar levelof physical realism [26, 27].We also note that the neutron star remnants observed in oursimulations are differentially rotating, and that the surround-ing disks have rotation profiles close to equilibrium. The post-merger disk is unstable to the axisymmetric magnetorotationalinstability, while the neutron star remnant has a complex rota-tion profile in which parts of the neutron star have an angularfrequency profile satisfying d Ω /dr > . These regions arestable to the magnetorotational instability, at least in the first
10 ms following the merger.The inclusion of a more realistic neutrino transport methodsignificantly modifies the composition of the disk and out-flows. Neutrino energy deposition also powers sustained out-flows with moderate electron fraction (on average (cid:104) Y e (cid:105) ∼ . ) up to the end of the simulations. These outflows are ex-pected to produce a wide range of elements through r-processnucleosynthesis, as their properties are right at the limit be-tween neutron-rich material producing mostly heavy elements( A > ), and less neutron-rich material producing lowermass elements. Finally, absorption of electron neutrinos andantineutrinos in the optically thick disk causes heating in thedisk. It also leads to a disagreement in the predicted neu-trino luminosity between the leakage scheme, which does notaccount for absorption, and the transport scheme. None ofthese differences significantly affect the hydrodynamic of theremnant, but they show the importance of neutrino transportwhen assessing the mass, composition, and geometry of theoutflows, the resulting properties of radioactively poweredelectromagnetic transients, and the result of r-process nucle-osynthesis in the outflows. The fact that neutron star mergerscan produce a wide range of elements when neutrino trans-port is taken into account has been discussed in more detailby post-processing the results of a higher mass neutron starmerger simulation [17, 20], and of simulations of post-mergeraccretion disk models [45]. Our present simulations confirmthese results in low-mass neutron star mergers. Additionally,they explicitly show the difference between simulations us-ing a simple leakage scheme, and simulations using neutrinotransport. The results presented here are limited by a number of im-portant assumptions in our simulations. The first is the studyof equal mass binaries. Indeed, we know that unequal andequal mass BNS mergers show significant differences. In par-ticular, unequal mass mergers eject a larger amount of mate-rial [28]. The second is the fact that we ignore magnetic fields.Over the short time scales considered here, magnetohydrody-namics effects are not expected to affect the evolution of neu-tron star remnant, but could drive additional outflows from thedisk [22, 27] . Over longer time scales, magnetic fields wouldbe critical to the spin evolution of the remnant neutron star, an-gular momentum transport, heating in the disk, and possiblythe formation of relativistic jets and magnetically-driven out-flows. Finally, the relatively small numerical grid on whichwe evolve the equations of general relativistic hydrodynam-ics limits our ability to measure the mass and properties ofthe outflows accurately. We expect to address these issues inthe future. However, we do not expect these assumptions tosignificantly affect the main results of this work, i.e. the prop-erties of the post-merger gravitational wave signal and the im-portance of the neutrino-matter interactions in the first
10 ms following the merger.
ACKNOWLEDGMENTS
The authors wish to thank Jan Steinhoff for communicat-ing values for the Love number or the neutron stars used inour simulations, Rodrigo Fernandez and Dan Kasen for use-ful discussions over the course of this work, and the mem-bers of the SxS collaboration for their advice and support.Support for this work was provided by NASA through Ein-stein Postdoctoral Fellowship grants numbered PF4-150122(F.F.) and PF3-140114 (L.R.) awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophys-ical Observatory for NASA under contract NAS8-03060; andthrough Hubble Fellowship grant number 51344.001 awardedby the Space Telescope Science Institute, which is operatedby the Association of Universities for Research in Astronomy,Inc., for NASA, under contract NAS 5-26555. The authorsat CITA gratefully acknowledge support from the NSERCCanada. M.D.D. acknowledges support through NSF GrantPHY-1402916. L.K. acknowledges support from NSF grantsPHY-1306125 and AST-1333129 at Cornell, while the au-thors at Caltech acknowledge support from NSF Grants PHY-1404569, AST-1333520, NSF-1440083, and NSF CAREERAward PHY-1151197. Authors at both Cornell and Caltechalso thank the Sherman Fairchild Foundation for their sup-port. Computations were performed on the supercomputerBriar´ee from the Universit´e de Montr´eal, and Guillimin fromMcGill University, both managed by Calcul Qu´ebec and Com-pute Canada. The operation of these supercomputers is fundedby the Canada Foundation for Innovation (CFI), NanoQu´ebec,RMGA and the Fonds de recherche du Qu´ebec - Nature etTechnologie (FRQ-NT). Computations were also performedon the Zwicky cluster at Caltech, supported by the ShermanFairchild Foundation and by NSF award PHY-0960291. Thiswork also used the Extreme Science and Engineering Dis-1covery Environment (XSEDE) through allocation No. TG- PHY990007N, supported by NSF Grant No. ACI-1053575. [1] J. Aasi et al. (LIGO Scientific Collaboration), Class. QuantumGrav. , 074001 (2015), arXiv:1411.4547 [gr-qc].[2] F. Acernese et al. (Virgo Collaboration), Class. Quantum Grav. , 024001 (2015), arXiv:1408.3978 [gr-qc].[3] Y. Aso, Y. Michimura, K. Somiya, M. Ando, O. Miyakawa,T. Sekiguchi, D. Tatsumi, and H. Yamamoto (The KA-GRA Collaboration), Phys. Rev. D , 043007 (2013),arXiv:1306.6747 [gr-qc].[4] B. D. Metzger and E. Berger, Astrophys. J. , 48 (2012),arXiv:1108.6056 [astro-ph.HE].[5] R. Mochkovitch, M. Hernanz, J. Isern, and X. Martin, Nature , 236 (1993).[6] W. H. Lee and W. Klu´zniak, in Gamma-Ray Bursts, 4th Hun-stville Symposium , American Institute of Physics ConferenceSeries, Vol. 428, edited by C. A. Meegan, R. D. Preece, &T. M. Koshut (1998) p. 798.[7] H.-T. Janka, T. Eberl, M. Ruffert, and C. L. Fryer, Astrophys.J. , L39 (1999), astro-ph/9908290.[8] W. Fong and E. Berger, Astrophys. J. , 18 (2013),arXiv:1307.0819 [astro-ph.HE].[9] J. M. Lattimer and D. N. Schramm, Astrophys. J. , 549(1976).[10] L.-X. Li and B. Paczynski, Astrophys. J. , L59 (1998),arXiv:astro-ph/9807272 [astro-ph].[11] L. F. Roberts, D. Kasen, W. H. Lee, and E. Ramirez-Ruiz,Astrophys. J. Lett. , L21 (2011), arXiv:1104.5504 [astro-ph.HE].[12] D. Kasen, N. R. Badnell, and J. Barnes, Astrophys. J. , 25(2013), arXiv:1303.5788 [astro-ph.HE].[13] M. Tanaka and K. Hotokezaka, Astrophys. J. , 113 (2013),arXiv:1306.3742 [astro-ph.HE].[14] J. Barnes and D. Kasen, Astrophys. J. , 18 (2013),arXiv:1303.5787 [astro-ph.HE].[15] J. Lippuner and L. F. Roberts, Accepted for publication in theAstrophys. J.; arXiv:1508.03133 (2015), arXiv:1508.03133[astro-ph.HE].[16] O. Korobkin, S. Rosswog, A. Arcones, and C. Winteler,Mon. Not. Roy. Astr. Soc. , 1940 (2012), arXiv:1206.2379[astro-ph.SR].[17] S. Wanajo, Y. Sekiguchi, N. Nishimura, K. Kiuchi, K. Kyu-toku, and M. Shibata, Astrophys.J.Lett. , L39 (2014),arXiv:1402.7317 [astro-ph.SR].[18] M. Shibata and K. ¯o. Ury¯u, Phys. Rev. D , 064001 (2000),gr-qc/9911058.[19] K. Barkett, M. A. Scheel, R. Haas, C. D. Ott, S. Bernuzzi,D. A. Brown, B. Szil´agyi, J. D. Kaplan, J. Lippuner, C. D.Muhlberger, F. Foucart, and M. D. Duez, ArXiv e-prints(2015), arXiv:1509.05782 [gr-qc].[20] Y. Sekiguchi, K. Kiuchi, K. Kyutoku, and M. Shibata, Phys.Rev. D , 064059 (2015), arXiv:1502.06660 [astro-ph.HE].[21] F. Foucart, E. O’Connor, L. Roberts, M. D. Duez,R. Haas, L. E. Kidder, C. D. Ott, H. P. Pfeiffer, M. A.Scheel, and B. Szilagyi, Phys. Rev. D , 124021 (2015),arXiv:1502.04146 [astro-ph.HE].[22] K. Kiuchi, K. Kyutoku, Y. Sekiguchi, M. Shibata, andT. Wada, Phys. Rev. D , 041502 (2014), arXiv:1407.2660[astro-ph.HE]. [23] K. Kiuchi, P. Cerd´a-Dur´an, K. Kyutoku, Y. Sekiguchi, andM. Shibata, ArXiv e-prints (2015), arXiv:1509.09205 [astro-ph.HE].[24] B. Giacomazzo, J. Zrake, P. C. Duffell, A. I. MacFadyen,and R. Perna, Astrophys. J. , 39 (2015), arXiv:1410.0013[astro-ph.HE].[25] K. Dionysopoulou, D. Alic, and L. Rezzolla, ArXiv e-prints(2015), arXiv:1502.02021 [gr-qc].[26] D. Neilsen, S. L. Liebling, M. Anderson, L. Lehner,E. O’Connor, et al. , Phys. Rev. D , 104029 (2014),arXiv:1403.3680 [gr-qc].[27] C. Palenzuela, S. L. Liebling, D. Neilsen, L. Lehner, O. L.Caballero, E. O’Connor, and M. Anderson, Phys. Rev. D ,044045 (2015), arXiv:1505.01607 [gr-qc].[28] K. Hotokezaka, K. Kiuchi, K. Kyutoku, H. Okawa,Y. Sekiguchi, M. Shibata, and K. Taniguchi, Phys. Rev. D , 024001 (2013), arXiv:1212.0905 [astro-ph.HE].[29] K. Hotokezaka, K. Kiuchi, K. Kyutoku, T. Muranushi,Y. Sekiguchi, M. Shibata, and K. Taniguchi, Phys. Rev. D , 044026 (2013), arXiv:1307.5888 [astro-ph.HE].[30] S. Bernuzzi, A. Nagar, S. Balmelli, T. Dietrich, and M. Ujevic,Phys.Rev.Lett. , 201101 (2014), arXiv:1402.6244 [gr-qc].[31] K. Takami, L. Rezzolla, and L. Baiotti, Phys. Rev. Lett. ,091104 (2014), arXiv:1403.5672 [gr-qc].[32] K. Takami, L. Rezzolla, and L. Baiotti, Phys. Rev. D ,064001 (2015), arXiv:1412.3240 [gr-qc].[33] R. De Pietri, A. Feo, F. Maione, and F. L¨offler, ArXiv e-prints(2015), arXiv:1509.08804 [gr-qc].[34] R. Gold, S. Bernuzzi, M. Thierfelder, B. Brugmann, andF. Pretorius, Phys.Rev. D86 , 121501 (2012), arXiv:1109.5128[gr-qc].[35] W. E. East and F. Pretorius, Astrophys.J. , L4 (2012),arXiv:1208.5279 [astro-ph.HE].[36] S. Bernuzzi, T. Dietrich, W. Tichy, and B. Bruegmann, Phys.Rev. D , 104021 (2014), arXiv:1311.4443 [gr-qc].[37] W. Kastaun, F. Galeazzi, D. Alic, L. Rezzolla, and J. A. Font,Phys.Rev. D88 , 021501 (2013), arXiv:1301.7348 [gr-qc].[38] W. Kastaun and F. Galeazzi, Phys. Rev.
D91 , 064027 (2015),arXiv:1411.7975 [gr-qc].[39] T. Dietrich, N. Moldenhauer, N. K. Johnson-McDaniel,S. Bernuzzi, C. M. Markakis, B. Bruegmann, and W. Tichy,(2015), arXiv:1507.07100 [gr-qc].[40] N. Tacik, F. Foucart, H. P. Pfeiffer, R. Haas, S. Ossokine,J. Kaplan, C. Muhlberger, M. D. Duez, L. E. Kidder,M. A. Scheel, and B. Szil´agyi, ArXiv e-prints (2015),arXiv:1508.06986 [gr-qc].[41] A. Bauswein, N. Stergioulas, and H.-T. Janka, Phys. Rev. D , 023002 (2014), arXiv:1403.5301 [astro-ph.SR].[42] A. Bauswein and N. Stergioulas, Phys. Rev. D , 124056(2015), arXiv:1502.03176 [astro-ph.SR].[43] R. Fern´andez and B. D. Metzger, Mon. Not. Roy. Astr. Soc. , 502 (2013), arXiv:1304.6720 [astro-ph.HE].[44] A. Bauswein, S. Goriely, and H.-T. Janka, ApJ , 78(2013), arXiv:1302.6530 [astro-ph.SR].[45] O. Just, A. Bauswein, R. Ardevol Pulpillo, S. Goriely, andH.-T. Janka, Mon. Not. Roy. Astr. Soc. , 541 (2015),arXiv:1406.2687 [astro-ph.SR]. [46] A. Perego, S. Rosswog, R. M. Cabez´on, O. Korobkin,R. K¨appeli, A. Arcones, and M. Liebend¨orfer, Mon. Not. Roy.Astr. Soc. , 3134 (2014), arXiv:1405.6730 [astro-ph.HE].[47] R. Fern´andez, D. Kasen, B. D. Metzger, and E. Quataert,Mon. Not. Roy. Astr. Soc. , 750 (2015), arXiv:1409.4426[astro-ph.HE].[48] M. D. Duez, Class. Quantum Grav. , 114002 (2010),arXiv:0912.3529 [astro-ph.HE].[49] M. Shibata and K. Taniguchi, Living Reviews in Relativity ,6 (2011).[50] L. Lehner and F. Pretorius, Annu. Rev. Astron. Astr. , 661(2014).[51] F. ¨Ozel, D. Psaltis, R. Narayan, and A. Santos Villarreal, As-trophys. J. , 55 (2012), arXiv:1201.1006 [astro-ph.HE].[52] .[53] A. Bauswein, H.-T. Janka, K. Hebeler, and A. Schwenk, Phys.Rev. D , 063001 (2012), arXiv:1204.1888 [astro-ph.SR].[54] L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen, andO. Rinne, Class. Quantum Grav. , S447 (2006), arXiv:gr-qc/0512093v3 [gr-qc].[55] B. v. L. A. Harten, P. D. Lax, SIAM Rev. , 35 (1983).[56] X.-D. Liu, S. Osher, and T. Chan, Journal of ComputationalPhysics , 200 (1994).[57] G.-S. Jiang and C.-W. Shu, Journal of Computational Physics , 202 (1996).[58] M. D. Duez, F. Foucart, L. E. Kidder, H. P. Pfeiffer, M. A.Scheel, and S. A. Teukolsky, Phys. Rev. D , 104015 (2008),arXiv:0809.0002 [gr-qc].[59] F. Foucart, M. B. Deaton, M. D. Duez, L. E. Kidder, I. Mac-Donald, C. D. Ott, H. P. Pfeiffer, M. A. Scheel, B. Szi-lagyi, and S. A. Teukolsky, Phys. Rev. D , 084006 (2013),arXiv:1212.4810 [gr-qc].[60] L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen,and O. Rinne, Class. Quantum Grav. , 447 (2006), gr-qc/0512093.[61] B. Szil´agyi, L. Lindblom, and M. A. Scheel, Phys. Rev. D ,124010 (2009), arXiv:0909.3557 [gr-qc].[62] T. Matsushima and P. S. Marcus, Journal of ComputationalPhysics , 365 (1995).[63] W. T. M. Verkley, Journal of Computational Physics , 100(1997).[64] B. Szil´agyi, Int. J. Mod. Phys. D23 , 1430014 (2014),arXiv:1405.3693 [gr-qc].[65] M. B. Deaton, M. D. Duez, F. Foucart, E. O’Connor, C. D. Ott,L. E. Kidder, C. D. Muhlberger, M. A. Scheel, and B. Szi-lagyi, Astrophys. J. , 47 (2013), arXiv:1304.3384 [astro-ph.HE].[66] F. Foucart, M. B. Deaton, M. D. Duez, E. O’Connor, C. D.Ott, R. Haas, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, andB. Szilagyi, Phys. Rev. D , 024026 (2014), arXiv:1405.1121[astro-ph.HE].[67] S. Rosswog and M. Liebend¨orfer, Mon. Not. Roy. Astron. Soc. , 673 (2003), arXiv:astro-ph/0302301 [astro-ph].[68] M. Ruffert, H.-T. Janka, K. Takahashi, and G. Schaefer, As-tron. Astrophys. , 122 (1997), astro-ph/9606181.[69] E. O’Connor and C. D. Ott, Class. Quantum Grav. , 114103(2010), arXiv:0912.2393 [astro-ph.HE].[70] A. Burrows, S. Reddy, and T. A. Thompson, Nuc. Phys. A , 356 (2006), astro-ph/0404432.[71] K. S. Thorne, Mon. Not. Roy. Astr. Soc. , 439 (1981).[72] C. Y. Cardall, E. Endeve, and A. Mezzacappa, Phys. Rev. D , 103004 (2013), arXiv:1209.2151 [astro-ph.HE].[73] H. P. Pfeiffer, L. E. Kidder, M. A. Scheel, and S. A. Teukolsky,Comput. Phys. Commun. , 253 (2003), gr-qc/0202096. [74] F. Foucart, L. E. Kidder, H. P. Pfeiffer, and S. A. Teukolsky,Phys. Rev. D , 124051 (2008), arXiv:arXiv:0804.3787 [gr-qc].[75] H. P. Pfeiffer, D. A. Brown, L. E. Kidder, L. Lindblom,G. Lovelace, and M. A. Scheel, Class. Quantum Grav. ,S59 (2007), gr-qc/0702106.[76] A. W. Steiner, M. Hempel, and T. Fischer, Astrophys. J. ,17 (2013), arXiv:1207.2184 [astro-ph.SR].[77] M. Hempel, T. Fischer, J. Schaffner-Bielich, andM. Liebend¨orfer, Astrophys.J. , 70 (2012),arXiv:1108.0848 [astro-ph.HE].[78] J. M. Lattimer and F. D. Swesty, Nucl. Phys. A535 , 331(1991).[79] T. Fischer, M. Hempel, I. Sagert, Y. Suwa, and J. Schaffner-Bielich, European Physical Journal A , 46 (2014),arXiv:1307.6190 [astro-ph.HE].[80] P. Demorest, T. Pennucci, S. Ransom, M. Roberts, andJ. Hessels, Nature , 1081 (2010), arXiv:1010.5788 [astro-ph.HE].[81] J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S.Lynch, M. H. van Kerkwijk, M. Kramer, C. Bassa, V. S.Dhillon, T. Driebe, J. W. T. Hessels, V. M. Kaspi, V. I.Kondratiev, N. Langer, T. R. Marsh, M. A. McLaughlin,T. T. Pennucci, S. M. Ransom, I. H. Stairs, J. van Leeuwen,J. P. W. Verbiest, and D. G. Whelan, Science , 448 (2013),arXiv:1304.6875 [astro-ph.HE].[82] A. W. Steiner, J. M. Lattimer, and E. F. Brown, Astroph.J. , 33 (2010), arXiv:1005.0811 [astro-ph.HE].[83] A. W. Steiner, J. M. Lattimer, and E. F. Brown, ApJ Letters , L5 (2013), arXiv:1205.6871 [nucl-th].[84] S. Guillot, M. Servillat, N. A. Webb, and R. E. Rutledge,Astrophys. J. , 7 (2013), arXiv:1302.0023 [astro-ph.HE].[85] S. Bernuzzi, T. Dietrich, and A. Nagar, Physical Review Let-ters , 091101 (2015), arXiv:1504.01764 [gr-qc].[86] G. B. Cook, S. L. Shapiro, and S. A. Teukolsky, Astrophys. J. , 203 (1992).[87] G. B. Cook, S. L. Shapiro, and S. A. Teukolsky, Astrophys. J. , 227 (1994).[88] J. D. Kaplan, C. D. Ott, E. P. O’Connor, K. Kiuchi, L. Roberts,and M. Duez, Astrophys. J. , 19 (2014), arXiv:1306.4034[astro-ph.HE].[89] M. Boyle and A. H. Mrou´e, Phys. Rev. D , 124045 (2009),arXiv:0905.3177 [gr-qc].[90] X. Zhuge, J. M. Centrella, and S. L. W. McMillan, Phys. Rev.D , 6247 (1994), arXiv:gr-qc/9411029.[91] R. Oechslin, S. Rosswog, and F.-K. Thielemann, Phys. Rev.D , 103005 (2002), arXiv:gr-qc/0111005.[92] M. Shibata and K. Uryu, Prog. Theor. Phys. , 265 (2002).[93] L. Baiotti, B. Giacomazzo, and L. Rezzolla, Phys. Rev. D ,084033 (2008), arXiv:0804.0594 [gr-qc].[94] K. Kiuchi, Y. Sekiguchi, M. Shibata, and K. Taniguchi,Phys.Rev. D80 , 064037 (2009), arXiv:0904.4551 [gr-qc].[95] N. Stergioulas, A. Bauswein, K. Zagkouris, and H.-T. Janka,Mon. Not. Roy. Astr. Soc. , 427 (2011), arXiv:1105.0368[gr-qc].[96] D. Shoemaker (LIGO Collaboration), “Advanced LIGO an-ticipated sensitivity curves,” (2010), LIGO DocumentT0900288-v3.[97] J. Clark, A. Bauswein, L. Cadonati, H.-T. Janka, C. Pankow,and N. Stergioulas, Phys. Rev. D , 062004 (2014),arXiv:1406.5444 [astro-ph.HE].[98] S. A. Balbus and J. F. Hawley, Astrophys. J. , 214 (1991).[99] F. H. Seguin, apj , 745 (1975). [100] B. Giacomazzo and R. Perna, Astrophys.J. , L26 (2013),arXiv:1306.1608 [astro-ph.HE].[101] L. F. Roberts, G. Shen, V. Cirigliano, J. A. Pons, S. Reddy,and S. E. Woosley, Phys. Rev. Lett. , 061103 (2012),arXiv:1112.0335 [astro-ph.HE]. [102] E. Gourgoulhon, ArXiv e-prints (2010), arXiv:1003.5015 [gr-qc].[103] L. Dessart, C. D. Ott, A. Burrows, S. Rosswog, and E. Livne,Astrophys. J. , 1681 (2009), arXiv:0806.4380.[104] B. D. Metzger and R. Fern´andez, MNRAS441