Dynamics, nucleosynthesis, and kilonova signature of black hole - neutron star merger ejecta
Rodrigo Fernández, Francois Foucart, Daniel Kasen, Jonas Lippuner, Dhruv Desai, Luke F. Roberts
DDynamics, nucleosynthesis, and kilonova signatureof black hole - neutron star merger ejecta
Rodrigo Fern´andez , , , Francois Foucart , ,Daniel Kasen , , , Jonas Lippuner , Dhruv Desai ,Luke F. Roberts Department of Physics, University of Alberta, Edmonton, AB T6G 2E1,Canada. Department of Physics, University of California, Berkeley, CA 94720, USA. Department of Astronomy & Theoretical Astrophysics Center, University ofCalifornia, Berkeley, CA 94720, USA. Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley,CA 94720, USA. NASA Einstein Fellow TAPIR, Walter Burke Institute for Theoretical Physics, California Institute ofTechnology, Pasadena, CA 99125, USA. National Superconducting Cyclotron Laboratory and Department of Physicsand Astronomy, Michigan State University, East Lansing, MI 48824, USAE-mail: [email protected]
Abstract.
We investigate the ejecta from black hole - neutron star mergersby modeling the formation and interaction of mass ejected in a tidal tail and adisk wind. The outflows are neutron-rich, giving rise to optical/infrared emissionpowered by the radioactive decay of r -process elements (a kilonova ). Here weperform an end-to-end study of this phenomenon, where we start from the outputof a fully-relativistic merger simulation, calculate the post-merger hydrodynamicalevolution of the ejecta and disk winds including neutrino physics, determinethe final nucleosynthetic yields using post-processing nuclear reaction networkcalculations, and compute the kilonova emission with a radiative transfer code.We study the effects of the tail-to-disk mass ratio by scaling the tail density. Alarger initial tail mass results in fallback matter becoming mixed into the diskand ejected in the subsequent disk wind. Relative to the case of a disk withoutdynamical ejecta, the combined outflow has lower mean electron fraction, fasterspeed, larger total mass, and larger absolute mass free of high-opacity Lanthanidesor Actinides. In most cases, the nucleosynthetic yield is dominated by the heavy r -process contribution from the unbound part of the dynamical ejecta. A Solar-like abundance distribution can however be obtained when the total mass of thedynamical ejecta is comparable to the mass of the disk outflows. The kilonovahas a characteristic duration of 1 week and a luminosity of ∼ erg s − ,with orientation effects leading to variations of a factor ∼ < Keywords : accretion, accretion disks – dense matter – gravitational waves –hydrodynamics – neutrinos – nuclear reactions, nucleosynthesis, abundances. a r X i v : . [ a s t r o - ph . H E ] J un H-NS merger ejecta
1. Introduction
The recent success of Advanced LIGO in detecting gravitational waves (GWs) frombinary black hole (BH-BH) mergers [1, 2] has marked the onset of GW Astronomy. Asthe sensitivity of detectors improves, detection of double neutron star (NS-NS) andNS-BH mergers in GWs is expected within the next few years [3]. Mergers that involveNSs are of interest to a broad community because they can help to test the equationof state (EOS) of dense matter (e.g., [4]), are a prime candidate astrophysical site for r -process elements (e.g., [5, 6, 7]), and are expected to produce electromagnetic (EM)counterparts over a wide range of timescales and wavelengths (e.g., [8, 9]).The most easily detectable EM counterpart of NS-NS/NS-BH mergers is asupernova-like transient powered by the radioactive decay of r -process elementsproduced in the expanding ejecta, commonly known as a kilonova or macronova [10, 11, 12, 13, 14]. These transients arise from the sub-relativistic ejecta of themerger and are hence not affected by beaming like short gamma-ray bursts (SGRBs).The emission is detectable at optical and infrared wavelengths [15].The intimate relation between a kilonova and the production of r -process elementsmakes the transient a powerful diagnostic of the physical conditions in the merger.This diagnostic power arises from the sensitivity of the optical opacity to the type of r -process composition of the ejecta: even a small fraction of Lanthanides or Actinides(mass number A (cid:38) . − . c (e.g., [23]). For NS-BH mergers, ejectionis driven mainly by tidal forces, and therefore the ejecta is launched primarily onthe equatorial plane. For NS-NS mergers additional ejection occurs from the contactinterface, leading to an outflow more widely distributed in solid angle. The tidallyejected material is mostly neutron-rich, generating a robust abundance of elementswith A >
130 that follow the Solar System r -process distribution [24, 25]. A shock-heated and/or neutrino-irradiated component of the interface ejecta can have highelectron fraction Y e , by virtue of weak interactions, thereby producing lighter r -process elements [26]. The prevalence of the high- Y e component is an open researchquestion, however, as it is sensitive to the treatment of neutrinos and the nuclear EOS[27, 28, 29, 6, 30, 31, 32].The second source of ejecta is the remnant accretion disk. Outflows can belaunched on the thermal time ( ∼ −
100 ms) if there is sufficient neutrino irradiationfrom a hypermassive NS [33, 34, 35], or on longer times (cid:38) ∼ . r -process elements. Whether heavier elements ( A >
H-NS merger ejecta . M (cid:12) ) relative to the bound dynamicalejecta (0 . M (cid:12) ).Here we improve upon the work of [41] by using initial conditions mapped froma fully relativistic NS-BH simulation that includes neutrino irradiation and thereforeprovides more realistic initial conditions. In addition, we explore the effect of varyingthe relative initial masses of dynamical ejecta and disk, by scaling the initial dynamicalejecta density. Given that the dynamical ejecta moves nearly ballistically, this is agood first approximation to test the impact the ratio of dynamical- to disk ejectawithout having to perform a large number of costly merger simulations. We conductthe post-merger simulations in axisymmetry, and include the dominant energy andlepton number source terms, in addition to passive tracer particles. The output is post-processed with a nuclear reaction network and a radiative transfer code to compute r -process nucleosynthesis yields and kilonova light curve and spectral predictions,respectively.The structure of this paper is the following. Section 2 describes our computationalmethod, Section 3 presents the results of our hydrodynamic simulations of the remnantevolution, Section 4 describes our nucleosynthesis results, Section 5 our kilonovapredictions, and Section 6 presents our summary and discussion.
2. Methods
The initial condition for our models is the final snapshot of a BH-NS post-mergerremnant simulation presented in [28]. This simulation evolves the accretion diskresulting from the merger of a 7 M (cid:12) BH with a 1 . M (cid:12) NS (case M14-7-S8 from [42]),up to 20 ms after the merger. The dimensionless spin of the BH before merger is χ = 0 .
8, aligned with the orbital angular momentum of the system. After merger, theBH spin is χ = 0 . H-NS merger ejecta K = 220 MeV. The computational domain in SpEC encompasses thecentral object and the accretion disk, while dynamical ejecta with a fallback timelonger than ∼
20 ms is allowed to leave the domain. Up to the disruption of theneutron star by the gravitational field of the black hole, the SpEC evolution explicitlyimposes symmetry of the fluid variables across the equatorial plane, while the diskformation and post-merger evolution is performed without that symmetry requirement(although SpEC preserves quite accurately the original symmetry of the system). Wenote that while for the system studied here the exact solution is symmetric acrossthe equatorial plane, explicitly imposing this condition in simulations can create somenumerical artifacts in the results. In particular, small errors in the vertical componentof the fluid velocity in the equatorial plane can push material away from that plane,leading to density maxima offset from the equator. While this effect decreases as thenumerical resolution increases, it can be easily observed in lower-resolution regions -e.g. the low density regions of the tidal tail, shown in Fig. 1.We generate initial conditions for our 2D models by axisymmetrizing the matterdistribution from the SpEC simulation, keeping only the data with density higherthan 10 g cm − , as lower density material is not reliably evolved in the generalrelativistic merger simulation. The axisymmetric baryon density ˜ ρ is obtained byintegrating the baryon density ρ over a curve L ( r, z ) of constant coordinate height z and constant cylindrical radius r , and then dividing by the proper length of L . Forother axisymmetric quantities ˜ X , we use density-weighted averages, i.e.˜ X ( r, z ) = (cid:72) L ( r,z ) XρW √ gdl (cid:72) L ( r,z ) ρW √ gdl (1)with W the Lorentz factor, and g the determinant of the spatial 3-metric.Theaxisymmetrized disk used in this work is constructed from (˜ ρ, ˜ P , ˜ v xT , ˜ v yT , ˜ v zT ), whereP is the pressure, and v iT = u i /u t is the 3-velocity of the fluid with respect to thenumerical grid (“transport” velocity). We also constructed an axisymmetrized diskfrom the average temperature ˜ T , and found no significant differences in the subsequentevolution of the system. We can assess the effect of this axisymmetrization procedurefrom Figs. 5, 7 and 9 of Foucart et al. 2015 [28], which show the density, velocity,temperature and electron fraction of the fluid in the equatorial plane of the system20 ms after merger. The less symmetric variable is the temperature, which varies byup to a factor of 2 at constant radius in the equatorial plane. Other fluid variablesare slightly more axisymmetric.The dynamical ejecta is reconstructed from the evolution of tracer particlesinitialized from an earlier snapshot of the evolution, at t = 4 ms after merger.This is early enough that all of the dynamical ejecta material is still within thecomputational domain, but late enough that a Newtonian ballistic and adiabaticapproximation for the subsequent evolution of the particles is acceptable (see e.g. [6],where smoothed particle hydrodynamics simulations of the tidal ejecta of a BH-NSmerger were performed with and without the inclusion of pressure terms, with littledifferences in the evolution of the ejecta). Using these assumptions we determine thelocation of the bound and unbound components of the dynamical ejecta at t = 20 msafter merger, and axisymmetrize that mass distribution. For this component we onlyconsider matter which is located at a radius r (cid:38) GM BH c − when the particles areintroduced. This includes all of the dynamical ejecta material which leaves the SpECcomputational grid during the post-merger evolution. To limit the error introducedby the transition from a general relativistic simulation to Newtonian trajectories, we H-NS merger ejecta v/c ) / − GM BH / ( rc ) = − u t , where u µ is the 4-velocity of thefluid in the general relativistic simulation.The disk data is mapped onto a pseudo-Newtonian potential for long-termevolution ( § r cut = 260 km is made to prevent spuriousjumps in the density distribution when combining the two components, resulting ina gap of ∼
50 km in radius between disk and dynamical ejecta. This gap is filledon a timescale of ∼ − s once evolution begins, and its introduction does not havea significant outcome in the evolution. For both components, the interpolation isperformed after shifting the radial coordinate by 0 . r g everywhere, which yields a diskmass that agrees well with that computed in the GR simulation. A similar shift wasused in [6] to provide a better agreement between GR and Newtonian potentials whenmapping data from SpEC simulations of a BH-disk system in the damped harmonicgauge to SPH simulations using a Paczynski-Witta potential.The baseline configuration obtained in this way has a central BH mass of 8 . M (cid:12) and spin χ = 0 . . M (cid:12) , gravitationally bound dynamical ejecta mass0 . M (cid:12) (hereafter called fallback ) and unbound dynamical ejecta mass 0 . M (cid:12) (hereafter called unbound tail ). The degree of gravitational binding is computedrelative to the pseudo-Newtonian potential. Passive scalars are assigned to each ofthese components in order to track their evolution and interaction as a function oftime (Figure 1). Additional models are evolved which scale the dynamical ejecta massor remove it altogether, as described in § We solve the time-dependent hydrodynamic equations using FLASH3 [50, 51] inaxisymmetric spherical polar coordinates ( r, θ ). The public version of the code hasbeen modified to include the physics required to model the long-term evolution ofcompact object merger remnants [38, 52, 53]. Angular momentum transport is treatedwith an α viscosity prescription [54], and neutrinos are included via a leakage schemewith lightbulb-type self-irradiation, including only charged-current weak interactions.The equation of state is that of Timmes et al. [55], and gravity is modeled with thepseudo-Newtonian potential of Artemova et al. [56].The computational domain extends from r in = 24 . r in , covering the full range of polar angles. Theradial grid is logarithmic with cell size ∆ r/r (cid:39) . θ , achieving ∆ θ = ∆ r/r (cid:39) ◦ on the equator. We employ standardoutflow boundary conditions in radius, and reflecting in polar angle at the axis.During the first 1 /
10 of an orbit ( ∼ × − s), the system is evolved without H-NS merger ejecta . . . . . . . . . x [cm] × − . − . − . . . . . × . . . . . . . x [cm] × − . − . − . . . . . z [ c m ] × . . disk [g cm − ] 10 . . fallback [g cm − ] 10 . . unbound tail [g cm − ]0 . . . . . . . . . x [cm] × − . − . − . . . . . × . . . . . . . x [cm] × − . − . − . . . . . z [ c m ] × .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 .
25 0 .
30 0 .
35 0 .
40 0 .
45 0 . Figure 1:
Top:
Partial densities (= ρX i , with X i the mass fraction) of disk, bounddynamical ejecta (“fallback”), and unbound dynamical ejecta (“unbound tail”) componentsat t = 0 for the baseline model F0. Note that the color scale is not the same for allcomponents; here it is chosen to maximize contrast. Bottom:
Electron fraction at t = 0in model F0. The dynamical ejecta density is reconstructed from tracer particles and re-mapped onto the spherical grid. Rough edges are smoothed out by evolving the entire systemfor ∼ − s without source terms. The gap between disk and dynamical ejecta – introducedto eliminate overlap between the disk distribution and dynamical ejecta distribution – is alsofilled within that timescale. H-NS merger ejecta ρ amb ∝ r − , normalized so that it is nearly 10 orders ofmagnitude lower than the torus density peak ( ∼
60 g cm − ) at the same location.The density floor is set to be 90% of the initial ambient density. As time elapses, thefloor is gradually decreased interior to 200 km until it becomes constant and equal to10 g cm − . Source terms are suppressed when the density is less than 10 times higherthan the floor.In order to generate homologously expanding ejecta for radiative transfercalculations, the material flowing out is sampled at r map = 10 cm throughout thesimulation. Following the approach of [40] and [41], the sampled material is injectedinto a new computational domain with an inner radius at r = r map and outer radius r = 10 r map . In these new hydrodynamic simulations, material is evolved withoutenergy source terms except for radioactive heating from r -process nucleosynthesis [24],and with a much lower ambient density, 10 − g cm − . By the time the forward shockreaches r = 3 × cm, the velocity profile is proportional to radius and the kineticenergy is dominant in most of the material. For each simulation, 10 tracer particles are inserted in the domain, with randompositions that follow the mass distribution (i.e., all particles represent the same amountof mass). Thermodynamic trajectories as a function of time are then generated fromeach particle by interpolating the corresponding variables from the grid at each time.In addition to density, temperature, and composition, we sample neutrino and viscoussource terms.Nucleosynthesis calculations on each trajectory are carried out with the nuclearreaction network code SkyNet [57]. The network includes 7843 isotopes up to Cn.Forward strong rates are taken from the JINA REACLIB database [58], with inverserates computed assuming detailed balance. Spontaneous and neutron-induced fissionrates are taken from [59], [60], [61], and [62]. Most weak rates come from [63], [64],and [65], or otherwise from REACLIB. Nuclear masses and partition functions aretaken from REACLIB, which contains experimental data where available or Finite-Range Droplet Macroscopic model (FRDM, e.g. [66]) values otherwise. The networkevolves the temperature taking into account source terms due to nuclear reactions ascomputed by SkyNet, as well as viscous heating and neutrino interactions (chargedcurrent emission/absorption on free nucleons) from the thermodynamic trajectories.Computations for each trajectory start from nuclear statistical equilibrium (NSE).Most particles that initially reside in the disk reach temperatures T ≥ K, andtheir nucleosynthesis is computed from the last time they reach this temperature. Mostdynamical ejecta particles have temperatures lower than 3 × K at all times, makingcomputation of NSE infeasible (and a bad approximation). We therefore ignore allparticles for which the maximum temperature is less than 5 × K, and instead settheir abundance to the average final yield (at a time 10 s) of model M14-7-S8 from [6].That calculation uses the same GR merger simulation to produce tracer particles onlyfrom dynamical ejecta material, and starts the evolution at an earlier time such thatparticles satisfy T ≥ × K. For the remaining particles, for which the maximumtemperature satisfies 5 ≤ ( T / K) ≤
10, nucleosynthesis starts from the time atwhich the maximum temperature is reached. The network calculations performed
H-NS merger ejecta × s, which is sufficiently long to obtain relatively stableabundances as a function of mass number A except for the heaviest nuclei ( A > A = 130 is sub-dominant, we consider this to be a reasonable approximation (we ignore elementalabundance changes as a function of Z that take place at times longer than 10 s). We calculate synthetic light curves and spectra of the models using the MonteCarlo radiative transfer code SEDONA [67]. The calculations use the axisymmetrichomologous ejecta profiles from the hydrodynamical calculations described in § r -process elementsis approximated statistically based on the atomic data of [16], and using the Y e from the hydrodynamical simulations as a proxy for the composition. This is a fairapproximation given the very sensitive dependence of the production of heavy r-processin Y e [40, 57]. We use Y e < .
25 as a threshold for the production of lanthanides.The radiative transfer calculations assume local thermodynamic equilibrium forthe atomic level populations and ionization state. We use the average time-dependentradioactive heating rate from [12] for all points of the ejecta except those with Y e > .
4, which we asssume have zero heating. We assume 100% thermalizationof the radioactive decay products – in reality, non-efficient thermalization will reducethe kilonova luminosity, especially for times after peak [21]. We compute wavelength-dependent light curves for 20 polar angles, equally-spaced in their cosine, over therange 0 to π . We average the results over three wavelength regions: blue optical(3500 − A ), red optical (5000 − A ), and infrared (1 − µ m) to construct thespecific luminosity in different bands. Table 1 summarizes the properties of the eight models computed in this study. Thefirst group includes a simulation of the disk (Fdisk) or dynamical ejecta alone (Fdyn),as well as a case of disk alone with Y e = 0 . Y e = 0 . α = 0 .
03 is used (c.f. [38]), and the initial blackhole mass and spin are used in the pseudo-Newtonian potential and kept constant intime.The second group of models evolves the combined dynamical ejecta plus disk.The fiducial case (F0) simulates the system with the physical parameters from theGR merger simulation, whereas the remaining models (Ft0.1, Ft0.3, and Ft3.0) scalethe dynamical ejecta mass by factors of 1 /
10, 1 /
3, and 3. While the range of scalingfactors used here is very wide compared to the numerical error in the determinationof the mass of the disk, bound tail, and unbound tail, they are fairly representativeof potential uncertainties in the relative mass of these components when varying theparameters of the binary (i.e., for black hole masses in the range M BH ∼ [5 − M (cid:12) ).In particular, models which reduce the ratio of the mass of the dynamical ejecta tothe mass of the disk are more representative of low-mass BH-NS mergers, for which H-NS merger ejecta Table 1:
Models evolved and summary of results. Columns from left to right show modelname, initial masses of disk, fallback, and unbound tail, fraction of the disk and fallbackejected in the wind, total ejected fraction f w , tot (eq. [4]), mass-weighted radial velocity andelectron fraction in the wind (including disk and fallback), mass ejected with Y e > .
25, andaccreted fractions of disk and fallback.
Model M d M f M ut f w , d f w , f f w , tot ¯ v r, w /c ¯ Y e, w M w , . f acc , d f acc , f (10 − M (cid:12) ) (%) (%) (%) (10 − ) (10 − M (cid:12) ) (%) (%)Fdisk 6 0 0 8 ... 8 3.9 0.35 5 92 ...Fdisk-Y0.1 8 ... 8 3.9 0.31 4 92 ...Fdisk-ns a a Model Fdisk-ns is evolved without momentum or energy source terms.the disk mass often is an order of magnitude larger than the mass of the dynamicalejecta (see e.g. Fig.11 of [68]). ‡ The scaling factor could be even smaller for lowermass BH-NS mergers ( M BH < M (cid:12) ) or NS-NS mergers.The evolution of the accretion disk and dynamical ejecta including all sourceterms is carried out for 3000 orbits at the initial density peak, or about ∼
10 seconds.This timescale corresponds to several viscous times, and is determined by demandingsaturation in the mass ejection as measured at the radius at which the disk outflow issampled ( r map = 10 cm). The second simulation step that uses the sampled outflowand evolves it into homomology is run until about 100 −
200 s. This time is chosen sothat most of the mass distribution achieves homology while satisfying the requirementthat the swept up mass (in ambient medium) remains lower than 1% of the disk ejectamass.
3. Dynamics of disk wind and dynamical ejecta
The initial condition for the baseline model F0 is shown in Figure 1. Passive scalarstrack the evolution of the disk material as well as the gravitationally bound andunbound parts of the dynamical ejecta ( fallback and unbound tail ). For clarity, westart by discussing the behavior of models that include only the disk or the dynamicalejecta, and then proceed to describe their combined evolution.At t = 0, the electron fraction in the disk and dynamical ejecta is very different.With the exception of very low density regions near its inner edge, most of thedynamical ejecta has Y e (cid:39) .
05. In contrast, the bulk of the disk has Y e (cid:39) . − . ‡ We note that this statement only applies to the ratio of the disk and dynamical ejecta masses. Inabsolute terms, and at constant black hole spin, low mass black holes are not necessarily unfavorableto the production of unbound material. In fact, for low spins, high mass black holes are unable todisrupt the neutron star and do not produce any ejecta at all [69].
H-NS merger ejecta − − − time [s]10 − − − − − acc r e ti on r a t ea t I S C O [ M (cid:12) s − ] (a) F0 totalF0 diskF0 fall F0 tailFdynFdisk 10 − time [s]10 − − − − ou t fl o w r a t ea t r = c m [ M (cid:12) s − ] (b) F0 totalF0 diskF0 fallF0 tail FdynFdiskFdisk-ns Figure 2:
Left:
Mass accretion rate at the ISCO for the different components of the baselinemodel F0 shown as solid lines: disk (red), fallback (blue), and unbound tail (green). Forcomparison, dotted lines show accretion of disk material for the model without dynamicalejecta (Fdisk, red) and of fallback material for the model without disk (Fdyn, blue).
Right:
Same as the left panel but now showing total mass ejection at r = 10 cm. The disk-onlymodel without source terms (Fdisk-ns) is shown with a thin grey line. In the absence of dynamical ejecta (model Fdisk), the disk follows the usualstages of neutrino-cooled disks (e.g., [70]). Two time scales govern the evolution ofthe system: the orbital time t orb (cid:39) (cid:18) M (cid:12) M (cid:19) / (cid:16) r
55 km (cid:17) / ms , (2)where M is the black hole mass, and the viscous time t visc (cid:39) (cid:18) . α (cid:19) (cid:18) . H/R (cid:19) (cid:18) M (cid:12) M (cid:19) / (cid:16) r
55 km (cid:17) / ms , (3)where H/R is the height-to-radius ratio of the disk. Initially, neutrino cooling isimportant, balancing viscous heating. As the disk evolves and spreads, neutrinoemission gradually decreases. At the same time, nuclear recombination adds energyin the outer parts of the disk. By t = t visc (cid:39)
200 ms, the neutrino luminosity hasdecreased by two orders of magnitude from its initial value, and the disk has reachedan advective state, leading to vigorous convection and mass ejection. As shown inTable 1 however, most of the initial disk mass is accreted onto the BH (Figure 2a),with only ∼
8% ejected as a wind. The bulk of mass ejection reaches 10 cm at times (cid:38) t (cid:46)
30 ms), this outflow is also present in acomparison model evolved without any source terms other than gravity (Fdisk-ns).In the absence of an accretion disk (model Fdyn), most of the fallback mattermoves quickly towards the BH, while the unbound tail moves outward. Bothcomponents also expand in polar angle, due to small pressure gradients neglected inthe ballistic evolution of the dynamical ejecta material. Initially, the fallback materialhas nearly constant specific angular momentum, in excess of the value at the innermoststable circular orbit (ISCO) by a factor of ∼ .
3. Therefore, the part of this materialthat falls towards the BH forms an accretion disk, and the resulting accretion rateat the ISCO evolves in a way similar to the original accretion disk after fallback material has reached the ISCO (Figure 2a). Meanwhile, the unbound tail and part
H-NS merger ejecta fallback matter move continuously outward. Over the course of a few orbits,the outward-moving dynamical ejecta establishes a power-law density profile in radiusthat transitions smoothly between bound and unbound components. A very smallamount of unbound tail material is mixed within the fallback at small radius, andmoves towards the BH due to lack of angular momentum. As shown in Table 1, 12%of the initial fallback mass reaches a radius of 10 cm (9% gravitationally unbound),while the rest is accreted. A fraction 99.6% of the initial unbound tail mass leaves thesystem.The early evolution of the combined disk and dynamical ejecta (model F0) isillustrated in Figure 3. The gap that initially separates the outer edge of the disk andthe inner edge of the bound tail is filled on a time t < t orb by expansion of the diskand inward motion of the bound tail. Fallback matter continues to penetrate the diskon the equatorial plane at later times, with mixing of the two components occurringon a timescale of ∼
10 orbits at the disk density peak. For comparison, model Ft0.1(dynamical ejecta density scaled down by a factor of 10 relative to the baseline modelF0) does not experience this violent mixing. Instead, the fallback matter mixes onlyaround the periphery of the disk and does not penetrate regions of higher density.This behavior confirms that the initial gap between disk and dynamical ejecta doesnot influence the dynamics significantly. Instead, the relevant parameter is the ratioof dynamical ejecta mass to disk mass (which is related to their densities).A snapshot of the longer-term expansion of the system – compared to the casesin which the disk and dynamical ejecta evolve independently – is shown in Figure 5a.The bulk of the unbound tail material moves outward identically regardless of whetherthe disk is present. By t = t visc (cid:39)
200 ms, the fallback material has fully mixed intothe original accretion disk, as shown by the similar shape of the density profiles ofthe two components inside r ∼
200 km. The original disk material is compressedby the infalling matter, maintaining a sharp density drop at its initial outer edge(200 km), and a higher density peak relative to the case without dynamical ejecta atthe same evolution time. Nevertheless, the very outer edge of the disk has still madeits way into the dynamical ejecta, reaching a similar radius ( ∼ fallback in higher abundance than disk material outside r ∼
200 km. These results indicate that (a) with the exception of the outermostregions of the dynamical ejecta, the evolution of which is unaffected by the disk, the fallback and disk components mix efficiently (for model F0), becoming effectively asingle disk, and that (b) the wind ejection mechanism is the same as in the case withoutdynamical ejecta. Changes in the wind properties must therefore be a consequenceof the different initial binding energies and masses of the disk and bound dynamicalejecta components.
The overall trends in mass ejection as a function of initial fallback and disk massesare summarized in Table 1 and Figure 6a. Ejection is quantified by the ratio oftotal ejected mass ( M w,i ) at r = 10 cm to the initial mass ( M i ) in disk or fallbackcomponents , f w,i = M w,i /M i , with i = { d, f } for disk or fallback, respectively. Wealso define a total ejected fraction of disk plus fallback material, f w,tot = f w,d M d + f w,f M f M d + M f . (4) H-NS merger ejecta -1012 z [ c m ] (a) t = . . . x [10 cm]-2-101 z [ c m ] (e) 11 ms 1 2 3 x [10 cm](f) 15 ms 1 2 3 x [10 cm](g) 22 ms 1 2 3 4 x [10 cm](h) 43 ms10 disk [g cm − ] 10 fallback [g cm − ] 10 unbound tail [g cm − ] Figure 3:
Snapshots of the early evolution of the baseline model F0, showing how the fallback material mixes into the disk. The color maps show partial densities ρX i , with X i themass fraction of disk, fallback , or unbound tail material, with an opacity that varies linearlyfrom 0 to 1 alongside the color table. Note that the density range for the color scale differsfrom Figure 1; here it is the same for all components. The orbital time at the location of theinitial disk density peak is 2 . -1012 z [ c m ] (a) t = . . . x [10 cm]-2-101 z [ c m ] (e) 11 ms 1 2 3 x [10 cm](f) 15 ms 1 2 3 x [10 cm](g) 22 ms 1 2 3 4 x [10 cm](h) 43 ms10 disk [g cm − ] 10 fallback [g cm − ] 10 unbound tail [g cm − ] Figure 4:
Same as Figure 3, but for model Ft0.1. The fallback matter has a lower densitythan in model F0, and only mixes in the periphery of the disk.
H-NS merger ejecta radius [cm]10 d e n s it y [ g c m − ] F0 diskF0 fallF0 tail FdiskFdyn fallFdyn tail radius [cm] d e n s it y [ g c m − ] Ft0.3 diskFt0.3 fallFt0.3 tail FdiskFdyn fall ( ÷ d e n s it y [ g c m − ] Ft3.0 diskFt3.0 fallFt3.0 tail FdiskFdyn fall ( × Figure 5:
Left:
Angle-averaged density profiles for the different components of the disk-dynamical ejecta remnant system at t = 216 ms: disk (red), fallback (blue), and unbound tail (green). For comparison, dotted lines show profiles for the models with disk alone (Fdisk,red) and dynamical ejecta alone (Fdyn, blue and green) at the same time. Right:
Same asleft panel but for models with scaled dynamical ejecta at t = 216 ms: Ft3.0 (top) and Ft0.3(bottom). Models with disk and dynamical ejecta alone (Fdisk and Fdyn, respectively) areshown with dotted lines. The curves from model Fdyn are scaled by the same factor asmodels Ft0.3 and Ft3.0 (as labeled). Note that in computing the ejected mass, we include all material that reaches r = 10 cm, irrespective of whether it is gravitationally unbound or not. Equation 4is equivalent to adding up the mass of disk and fallback material ejected and dividingthe result by the sum of initial disk and fallback masses.The fractions of ejected disk and fallback material both decrease for increasing M f /M d , while the total fraction increases to a peak of f w,tot ∼
15% around M f ∼ M d .This behavior arises from the large fraction of fallback material that is turned aroundwhenever a disk is present ( > fallback densities for models with small and large M f /M d ratios (Ft0.3 and Ft3.0, respectively). Also shown for comparison in each panel is the fallback profile in the model with no disk (Fdyn), scaled to the appropriate dynamicalejecta normalization, and the disk profile in the model with no tail (Fdisk). Whenthe fallback mass is small, bound tail material causes only minor modifications to theradial disk profile, with most of the change occurring in the outer layers. In contrast,for a large initial M f /M d (model Ft3.0), fallback material mixes very effectively intothe disk, compressing it significantly relative to the case without dynamical ejecta.The presence of the disk stops fallback from immediately reaching the BH anddelays the onset of accretion of this material (as in e.g. Figure 2b). This partiallyaccounts for the fact that the ratio of ejected to accreted matter in fallback materialincreases as M f /M d decreases. Additionally, the (negative) net specific energy of thefluid is initially an order of magnitude larger at the disk density peak than within thebound dynamical ejecta material. The fallback component is thus less gravitationallybound than the disk component, and is relatively easier to eject. H-NS merger ejecta ∞ M f / M d . . . . e j ec t e d fr ac ti on diskfallback ( × . ∞ M f / M d − d i s k a nd f a ll b ac k m a ss e j ec t e d [ M (cid:12) ] total Y e > . Figure 6:
Left:
Fraction of the initial mass in disk (red), fallback (blue), and total disk plus fallback material (black, eq. [4]) that is ejected as a wind at r = 10 cm, as a function ofthe ratio of initial fallback to disk masses, for constant initial disk mass (c.f. Table 1). Notethat the fraction of fallback matter ejected is scaled down by a factor of 10 to fit in the plot. Right:
Total mass (black) and mass with Y e > .
25 (green) ejected in the disk wind as afunction of the ratio of initial fallback to disk masses. The mass with Y e > .
25 is a proxyfor the amount of Lanthanide- and Actinide-free mass. See Table 1 for numerical values.
Table 1 shows the mass-flux-weighted radial velocity and electron fraction of the diskand fallback material ejected at r = 10 cm for all models evolved. For increasing M f /M d , the outflow is faster and more neutron rich. This trend derives from the factthat the fallback is initially very neutron rich and less gravitationally bound (per unitmass) than the disk material. When both components mix, a larger mass-weightedoutflow velocity is obtained for fixed energy input per unit mass. Also, materialnear the transition between bound and unbound dynamical ejecta moves outward andcontributes to increase the total outflow momentum.For a more detailed breakdown of the outflow properties, Figure 7 shows masshistograms of electron fraction and radial velocity for the baseline model F0, and themodels without dynamical ejecta (Fdisk) and without disk (Fdyn), for comparison.The bulk of the unbound tail satisfies Y e (cid:46) . fallback component has asimilar electron fraction distribution as the unbound tail , and its velocity distributionis a continuation of the unbound tail distribution to lower velocities. When the disk ispresent, the fallback material acquires a broad Y e distribution that stretches fromthe initial unbound tail region to the disk Y e distribution. Material that mixesdeep into the disk is subject to weak interactions, acquiring a similar compositionas disk material, while matter that remains outermost in radius preserves the originaldynamical ejecta composition, with a smooth transition arising from the gradualdecrease in strength of weak interactions with increasing radius (temperature). Thedisk Y e distribution acquires a low- Y e tail stretching to Y e ∼ .
06, and the high- Y e cutoff is decreased from 0 .
45 to 0 .
4. The velocity distribution of the fallback materialin the baseline model F0 is essentially a scaled-up version of the distribution of thedisk component in the same model, again pointing to mixing. The high velocity end isdominated by turned-around fallback matter , which replaces matter at the outer edgeof the disk as the fastest and earliest ejecta in the polar direction (relative to model
H-NS merger ejecta . . . . . . Y e − − − − − m a ss i nb i n [ M (cid:12) ] (a) F0 totalF0 diskF0 fall F0 tailFdiskFdyn fall − − v r / c − − − m a ss i nb i n [ M (cid:12) ] (b) Figure 7:
Mass histograms of material ejected at r = 10 cm as a function of electronfraction (left) and radial velocity (right). Solid lines and dashed lines correspond to thebaseline model F0, while dotted lines show the disk-only case (Fdisk, red) and dynamicalejecta-only model (Fdyn, blue), for comparison. Fdisk). The disk component of model F0 loses the low-end of its velocity distributionrelative to model Fdisk. The low-velocity end is now dominated by outward-moving fallback matter far from the bound-unbound transition.Figure 6b shows ejected mass with Y e > .
25, which is indicative of material thatcannot form Lanthanides and Actinides and therefore has an optical opacity similar toiron-group elements, thus contributing to an optical kilonova signal (e.g. [40]). Eventhough the outflow becomes on average less neutron rich with increasing M f /M d , thereis more Lanthanide-free mass ejected until M f ∼ M d , due to the increasing total massof the outflows.Finally, we explore the effect of the initial conditions on the composition of thedisk outflow. Figure 8 compares the Y e distribution of two models that evolve thedisk alone: one using the electron fraction from the merger simulation (Fdisk), andanother one in which Y e = 0 . .
35 and 0 .
31 for Fdisk and Fdisk-Y0.1, respectively.Aside from the slightly lower mean value for model Fdisk-Y0.1, its distribution hasan extended low- Y e tail and a sharper high- Y e cutoff relative to model Fdisk.The insensitivity to initial value of Y e can be contrasted with the results ofsimulations starting from less compact accretion disks, as is obtained in NS-NSmergers. For example, [38] simulated disks with an initial density peak at r peak ∼ (10 − GM BH /c , as opposed to r peak ∼ GM BH /c here, with the same codeand physical model, finding an average electron fraction Y e ∼ . − .
2. Thisis significantly different than the results obtained here. The compactness of thedisk is determined primarily by the mass of the black hole. General relativisticsimulations yield density peaks of the post-merger disk at r peak ∼
50 km (once amostly axisymmetric disk forms at times 10 −
20 ms after merger), without muchdependence on the mass of the black hole [72, 73]. Accordingly, more massive blackholes have more compact disks.We conclude that the bulk of the disk outflow acquires its composition from the
H-NS merger ejecta . . . . . . Y e − − − m a ss i nb i n [ M (cid:12) ] FdiskFdisk-Y0.1 A − − − − − − a bund a n ce X ( A ) / A solarFdiskFdisk-Y0.1 Figure 8:
Left:
Mass histogram as a function of electron fraction for outflow materialmeasured at r = 10 cm in disk-only models that use an initial Y e from the merger simulation(Fdisk) or that impose Y e = 0 . Right:
Nucleosynthesis abundancesobtained with thermodynamic trajectories from models Fdisk and Fdisk-Y0.1, normalized sothat the mass fractions X ( A ) integrate to unity. The Solar System abundances from [74] arescaled to model Fdisk at A = 130. compactness of the disk, which regulates the temperature and therefore the weakinteraction timescale and beta equilibrium Y e . Memory of the initial conditionsimparts modifications at the ∼
10% level in the electron fraction distribution.
4. Nucleosynthesis
Figure 9a shows final abundances for model F0 as a whole and separated by ejectacomponent. The unbound tail produces heavy r -process elements, with a deficit below A ∼
90 (c.f. [6]), while the disk outflow generates mostly elements with
A <
A >
140 thanthe disk material: although part of it undergoes reprocessing by weak interactionsupon mixing, it is initially more neutron-rich. Given the relative masses of the ejectacomponents (Table 1), yields are dominated by the unbound tail for A (cid:38)
90. Thecombined abundance distribution under-produces elements with A (cid:46) A = 130 when normalizing to the rare-Earth peak. In particular,model Ft0.1 (smallest dynamical ejecta mass) is the only case that reaches relativeabundance values comparable to Solar at the first r -process peak ( A = 80), with goodoverall agreement elsewhere except for some overproduction at A = 100. On the otherhand, model Ft3.0 (highest dynamical ejecta mass) is hardly different than a casewithout disk (Fdyn) for A (cid:38) Y e = 0 . t = 0. Abundances below A = 130 are very similar in both cases, with a ratioof first to second peak abundances that compares well with the Solar System. Thisrobustness below A = 130 is consistent with previous nucleosynthesis calculations fromthe disk outflow alone [5, 7]. The low- Y e component in model Fdisk-Y0.1 (Fig. 8a) H-NS merger ejecta A − − − − − − a bund a n ce X ( A ) / A solartotaldisk falltail 0 50 100 150 200 250mass number A − − − − − − a bund a n ce X ( A ) / A solarF0Fdyn Ft0.1Ft0.3Ft3.0 Figure 9:
Left:
Abundances for model F0 broken down by components and weighted by theirrelative masses.
Right:
Combined abundance distributions for models with scaled dynamicalejecta. Model abundance curves and Solar System values are normalized to the rare-Earthpeak ( A = 164) of model F0. results in a higher production of elements with A >
130 than model Fdisk, with anexcellent ratio of third to second peak abundances and underproduction of rare-Earthelements. Both models display an abundance excess at A = 100 relative to the SolarSystem distribution. Also, both cases contain an abundance spike at A = 132 whichis related to the convective character of the disk outflow and an incomplete treatmentof nuclear heating in the simulations [7].
5. Kilonova Signature
Figure 10 shows the synthetic kilonova light curves for model F0. The bolometricluminosity peaks ≈ t (cid:46) Y e > .
25 (Lanthanide-free) wind ejecta has only a minor impact on the light curves. Figure 10 shows thatthis high- Y e material is surrounded by lanthanide-rich material in all directions (thehigh- Y e ‘jet’ along the polar axis has very low density).The spectral evolution of the F0 model (shown in Figure 11) helps clarify thebroadband light curve behavior. At early times ( t = 0 . (cid:38) t (cid:38) ≈ ≈ µ m), showing little color evolutionthereafter. This evolution of the SEDs is similar to that shown for 1D parameterizedmodels by [16], who also considered how the early optical emission is subject touncertainties in the lanthanide opacities. In addition, because the photosphere atearly times forms in the outermost layers of ejecta, the early optical emission may besensitive to how well the model resolves the low density surface layers.The broadband light curves vary by a factor of ≈ H-NS merger ejecta time (days) ν L ν ( e r g s − ) blueredIR -3 -2.5 -2 -1.5 -1 -0.5 0 x [10 cm] -3-2-10123 z [ c m ] v r / c x [10 cm] Y e . . . . . Figure 10:
Left:
Broadband light curves from model F0 in the wave bands 3500 − A (blue), 5000 − A (red), and 1 − µ m (green). For each color, 10 viewing angles equallyspaced in their cosine are shown with different shades, spanning the range θ = 0 (light,rotation axis) to θ = 90 ◦ (dark, equatorial plane). Right:
Snapshot of the radial velocity(left panel) and electron fraction (right panel) in model F0 at time t = 150 s, when most ofthe matter distribution has reached homology. This snapshot is used as input for radiativetransfer calculations. The white contour corresponds to a density 10 − g cm − . wavelength ( ˚ A ) . . . . . . . . s p ec i fi c l u m i no s it y ( e r g ss − ˚ A − ) θ = ◦ θ = ◦ wavelength ( ˚ A ) . . . . . . . . θ = ◦ θ = ◦ Figure 11:
Spectra of model F0 at 0.5 days (left panel) and 5.0 days (right panel) aftermerger, as observed from both polar ( θ = 0 ◦ ) and equatorial ( θ = 90 ◦ ) viewing angles. Thehigher equatorial expansion velocity of unbound tail ejecta leads to a greater blueshift of thespectrum for the θ = 90 ◦ viewing angles. The high frequency variations are Monte Carlonoise. H-NS merger ejecta time (days) ν L ν ( e r g s − ) blue Ft3.0F0Ft0.3Ft0.1 time (days) infrared Figure 12:
Angle-averaged, broadband light curves from models that scale the dynamicalejecta mass, as labeled. The left panel shows the blue optical band (3500 − A ) and theright panel the infrared band (1 − µ m). directions closer the polar axis ( θ ≈ ◦ ), due to the fact that the projected surfacearea of the unbound tail is larger by a factor of ≈ ≈ . c , for polar viewing anglesand of ≈ . c for equatorial viewing angles (Figure 10). Because a greater blueshiftmoves more of the blue edge of SED into the optical bands, the optical light curvesare always brightest from equatorial viewing angles (see Figure 11).Figure 12 shows the blue and infrared light curves of the entire set of models thatscale the dynamical ejecta mass, averaged over all viewing angles. Peak luminositiesvary from 6 × erg s − for the dimmest blue optical case to 7 × erg s − for thebrightest infrared luminosity. The peak time and luminosity of the infrared emissionare monotonic functions of the dynamical ejecta mass, consistent with the largeramount of radioactive material and longer diffusion time. The SEDs and orientationeffects of these models are similar to those described for model F0.Kasen et al. [40] studied the light curves of kilonovae from disk winds, and showedthe degree of blue optical emission depended on the amount of lanthanide free ( Y e (cid:38) .
25) material ejected. They further showed that blue wind emission may be obscuredfrom certain viewing angles if a small mass of lanthanide-containing dynamical ejectaoverlaid the wind. In the models of this paper, the lanthanide-free disk ejecta isgenerally obscured from all viewing angles due to the near complete covering factorof low- Y e unbound tail material (Figure 10). The light curves of these models aretherefore mostly insensitive to the amount of lanthanide-free ejecta. As discussedabove, the optical luminosity arises from the lanthanide-containing dynamical ejecta,which is hot enough at early times to include some optical component. H-NS merger ejecta
6. Summary
We have investigated the long-term evolution of the remnant system left behind aftera BH-NS merger, and its electromagnetic and nucleosynthetic signatures. Using initialconditions from a GR merger simulation with neutrino transport, the remnant compo-nents are evolved in axisymmetry, accounting for the dominant nuclear and neutrinosource terms. Nucleosynthesis is computed from passive tracer particles injected inthe simulation, and optical/infrared kilonova light curves are generated using a MonteCarlo radiative transfer code once the ejecta has reached homologous expansion. Wehave characterized the dynamics of the system and its observable signatures when therelative initial masses of disk and dynamical ejecta are varied. Our main results arethe following:1. – The presence of the disk prevents fallback material from directly reaching the BH.Instead, fallback can mix efficiently into the disk, increasing the effective disk mass anddecreasing the net specific binding energy (Fig. 5). The resulting disk outflow is fasterand more neutron-rich as the dynamical ejecta mass increases, at fixed disk mass (Ta-ble 1). The unbound dynamical ejecta evolves independently of the other components.2. – The total amount of mass ejected in the disk wind outflow increases as the dy-namical ejecta mass increases – at fixed disk mass – up to the point where the bounddynamical ejecta mass and disk mass are approximately equal at t = 0 (Fig. 6). Themass in Lanthanide- and Actinide-free composition also increases up to this point.3. – The composition of the disk outflow is determined primarily by the initial com-pactness of the disk (mass of the BH divided by the radius of the disk density peak).The compactness sets the strength of the gravitational field, and therefore the disktemperature, the weak interaction timescale, and the beta equilibrium Y e of the disk.Memory of the initial composition contributes at the ∼
10% level (Fig. 8).4. – The nucleosynthesis output in the model mapped from the relativistic mergersimulation (F0) is dominated by the unbound tail ( A (cid:38) A (cid:46) r -process distribution is obtained for small dynamicalejecta masses (model Ft0.1), more suitable for small BH masses (Fig. 9). Excess pro-duction of elements around A = 100 and at A = 132 is obtained.5. – The kilonova signature produces optical emission for the first day after merger,then evolves to the infrared. The peak luminosities are of order 10 erg s − , with avariation of a factor ∼ detailed properties ofthe early optical emission are sensitive to the outer edges of the disk and dynamicalejecta, and thus will depend on how well resolved are low-density regions in simula-tions. Nevertheless we consider the general properties of our results are robust.Our results on the long-term ejecta dynamics of the model with physical H-NS merger ejecta fallback and disk material in the outerparts of the disk outflow (see e.g. Figure 3 in that paper and the insensitivity of theelectron fraction of the disk outflow to the presence or absence of dynamical ejectain their Table 2). This difference can be attributed to the different relative masses ofdisk and bound dynamical ejecta, and on the way the initial conditions are obtainedhere. In our model F0, the ratio of fallback-to-disk masses is 0.038/0.060, whereasin the baseline model C2d of Ref. [41] this ratio is 0.02/0.20. The bound dynamicalejecta reaches a peak density a factor of ∼
16 lower in C2d than in our model F0,hence it is not able to penetrate the disk and mix in efficiently (in contrast to ourFigure 3). A similar behavior is observed in model Ft0.1, for which the bulk of thedisk remains unaffected by fallback material, which only mixes in on its periphery(Figure 4). Still, model Ft0.1 has a mean Y e of the disk outflow (0.30) that is lowerby 14% relative to the model without dynamical ejecta (Fdisk, 0.35), whereas thecorresponding fractional decrease in the models of Ref. [41] is about half (0.29 to0.27, respectively). This lower mixing efficiency could result from a number of factors,including a higher concentration of dynamical ejecta in the equatorial plane in thegeneral relativistic simulation, and on the fact that the disk and dynamical ejectaused by Ref. [41] are not computed separately, hence some mixing is already builtinto it (the mean Y e of the resulting outflow compares favorably to the average betweenour models Ft0.1 and Ft0.3, which bracket the initial fallback-to-disk mass ratio inmodel C2d of [41]).An important observational implication of our results is that a substantial amountof blue optical emission can still be generated by the Lanthanide-rich ejecta at earlytimes when the temperatures are high. The duration of this signal is (cid:46) erg s corresponds toan absolute magnitude M B (cid:39) −
14 (apparent magnitude m B (cid:39) . ≈
2. A similar magnitudeof viewing angle effects was found by [75]. Kyutoku et al. [76] emphasized theasymmetries intrinsic in the BH-NS merger dynamical ejecta due to the dominanceof tidal forces on the ejection. They also estimated a factor of a few variation in thekilonova light curve properties at peak as a function of viewing angle. Tanaka et al.[77] carried out multi-dimensional Monte Carlo radiative transfer, and showed thatin addition to variations in the line of sight, the asymmetries in the dynamical ejectacause the emission to be systematically bluer than in the case of NS-NS mergers. Incontrast to their results, however, we find that asymmetries due to the line of sightpersist beyond 10 days, particularly in the optical band. Finally, Kawaguchi et al.[78] conducted a parameter space study of kilonova properties from BH-NS mergersby developing a fitting formula for the dynamical ejecta mass and for the kilonovalight curves as a function of binary mass ratio, BH spin, and NS size (EOS). Ourscaling of dynamical ejecta mass at fixed disk mass is a reasonable approximation to
H-NS merger ejecta ∼ t − / obtained whenassuming a flat energy distribution of material in Keplerian orbits [79, 80, 81]. Areliable prediction of this time-dependence requires carrying out disk simulations inmagnetohydrodynamics (MHD). For instance, the normalization and time-dependenceof the disk accretion that we obtain in our models (e.g., Figure 11 of [41]) wouldimply insufficient accretion energy as required to power the X-ray emission fromGRB 130603B at a time of ∼ r -process elements, which when combined with coupling to the actual composition ofthe outflow (instead of just Y e ), would yield spectral predictions alongside broadbandlight curves. The ejecta and disk properties can be sensitive to parameter dependenciesother than just the mass ratio, such as the rotation of the neutron star (e.g., [84])or effects due to non-circularity of the binary (such as e.g., formation in dynamicalcaptures, [85]). The development of open-source radiative transfer codes would alsoallow end-to-end studies to be carried out by a wider community; to our knowledge,the only open-source tool available to date is SNEC [86]. Further improvements inthe calculations will reveal whether the trends found here are robust. Acknowledgments
We thank Kenta Hotokezaka, Masaomi Tanaka, Kunihito Ioka, and Shinya Wanajofor useful discussions. We also thank the anonymous referee for helpful commentsthat improved the presentation of this paper. RF acknowledges support from theUniversity of California Office of the President, from NSF grant AST-1206097, andfrom the Faculty of Science at the University of Alberta. Support for this workwas provided by NASA through Einstein Postdoctoral Fellowship grant numberedPF4-150122 (FF) awarded by the Chandra X-ray Center, which is operated by theSmithsonian Astrophysical Observatory for NASA under contract NAS8-03060. DK issupported in part by a Department of Energy Office of Nuclear Physics Early CareerAward, and by the Director, Office of Energy Research, Office of High Energy andNuclear Physics, Divisions of Nuclear Physics, of the U.S. Department of Energyunder Contract No. DE-AC02-05CH11231. The software used in this work was inpart developed by the DOE NNSA-ASC OASCR Flash Center at the University ofChicago. This research used resources of the National Energy Research ScientificComputing Center (NERSC), which is supported by the Office of Science of the U.S.
H-NS merger ejecta
Edison (repositories m1186 and m2058). This research also used the
Savio computational cluster resource provided by the Berkeley Research Computingprogram at the University of California, Berkeley (supported by the UC BerkeleyChancellor, Vice Chancellor of Research, and Office of the CIO). We also thank thehospitality of the Yukawa Institute for Theoretical Physics through the workshopNuclear Physics, Compact Stars, and Compact Star Mergers 2016, during which partof this work was carried out.
References [1] Abbott B P et al.
PRL
Preprint )[2] Abbott B P et al.
PRL et al.
Classical and Quantum Gravity Preprint )[4] Bauswein A, Baumgarte T W and Janka H T 2013
PRL
Preprint )[5] Just O, Bauswein A, Pulpillo R A, Goriely S and Janka H T 2015
MNRAS
Preprint )[6] Roberts L F, Lippuner J, Duez M D, Faber J A, Foucart F, Lombardi Jr J C, Ning S, Ott C Dand Ponce M 2017
MNRAS
Preprint )[7] Wu M R, Fern´andez R, Mart´ınez-Pinedo G and Metzger B D 2016
MNRAS
Preprint )[8] Rosswog S 2015
Int. J. of Mod. Phys. D Preprint )[9] Fern´andez R and Metzger B D 2016
Annu. Rev. Nuc. Part. Sci.
23 (
Preprint )[10] Li L X and Paczy´nski B 1998
ApJ
L59–L62 (
Preprint astro-ph/9807272 )[11] Metzger B D, Mart´ınez-Pinedo G, Darbha S, Quataert E, Arcones A, Kasen D, Thomas R,Nugent P, Panov I V and Zinner N T 2010
MNRAS
Preprint )[12] Roberts L F, Kasen D, Lee W H and Ramirez-Ruiz E 2011
ApJ
L21 (
Preprint )[13] Tanaka M 2016
Advances in Astronomy
Preprint )[14] Metzger B D 2016 preprint, arXiv:1610.09381 ( Preprint )[15] Metzger B D and Berger E 2012
ApJ
48 (
Preprint )[16] Kasen D, Badnell N R and Barnes J 2013
ApJ
25 (
Preprint )[17] Fontes C J, Fryer C L, Hungerford A L, Hakel P, Colgan J, Kilcrease D P and Sherrill M E 2015
High Energy Density Physics ApJ
18 (
Preprint )[19] Tanaka M and Hotokezaka K 2013
ApJ
113 (
Preprint )[20] Hotokezaka K, Wanajo S, Tanaka M, Bamba A, Terada Y and Piran T 2016
MNRAS
Preprint )[21] Barnes J, Kasen D, Wu M R and Mart´ınez-Pinedo G 2016
ApJ
110 (
Preprint )[22] Rosswog S, Feindt U, Korobkin O, Wu M R, Sollerman J, Goobar A and Martinez-Pinedo G2016 preprint, arXiv:1611.09822 ( Preprint )[23] Hotokezaka K, Kiuchi K, Kyutoku K, Okawa H, Sekiguchi Y i, Shibata M and Taniguchi K 2013
Phys. Rev. D Preprint )[24] Korobkin O, Rosswog S, Arcones A and Winteler C 2012
MNRAS
Preprint )[25] Bauswein A, Goriely S and Janka H T 2013
ApJ
ApJ
L39(
Preprint )[27] Sekiguchi Y, Kiuchi K, Kyutoku K and Shibata M 2015
PRD Preprint )[28] Foucart F, O’Connor E, Roberts L, Duez M D, Haas R, Kidder L E, Ott C D, Pfeiffer H P,Scheel M A and Szilagyi B 2015
PRD Preprint )[29] Palenzuela C, Liebling S L, Neilsen D, Lehner L, Caballero O L, O’Connor E and Anderson M2015
PRD Preprint )[30] Lehner L, Liebling S L, Palenzuela C, Caballero O L, O’Connor E, Anderson M and Neilsen D2016
Classical and Quantum Gravity Preprint )[31] Radice D, Galeazzi F, Lippuner J, Roberts L F, Ott C D and Rezzolla L 2016
MNRAS
Preprint )[32] Foucart F, O’Connor E, Roberts L, Kidder L E, Pfeiffer H P and Scheel M A 2016
Phys. Rev. D,in press ( Preprint ) H-NS merger ejecta [33] Ruffert M, Janka H T, Takahashi K and Schaefer G 1997 A&A
Preprint arXiv:astro-ph/9606181 )[34] Dessart L, Ott C D, Burrows A, Rosswog S and Livne E 2009
ApJ
Preprint )[35] Perego A, Rosswog S, Cabez´on R M, Korobkin O, K¨appeli R, Arcones A and Liebend¨orfer M2014
MNRAS
Preprint )[36] Metzger B D, Piro A L and Quataert E 2009
MNRAS
Preprint )[37] Lee W H, Ramirez-Ruiz E and L´opez-C´amara D 2009
ApJ
L93–L96 (
Preprint )[38] Fern´andez R and Metzger B D 2013
MNRAS
Preprint )[39] Kiuchi K, Sekiguchi Y, Kyutoku K, Shibata M, Taniguchi K and Wada T 2015
PRD Preprint )[40] Kasen D, Fern´andez R and Metzger B D 2015
MNRAS
Preprint )[41] Fern´andez R, Quataert E, Schwab J, Kasen D and Rosswog S 2015
MNRAS
Preprint )[42] Foucart F, Deaton M B, Duez M D, O’Connor E, Ott C D, Haas R, Kidder L E, Pfeiffer H P,Scheel M A and Szilagyi B 2014
PRD Preprint )[43] [44] Lindblom L, Scheel M A, Kidder L E, Owen R and Rinne O 2006
Class. and Quantum Grav.
447 (
Preprint gr-qc/0512093 )[45] Duez M D, Foucart F, Kidder L E, Pfeiffer H P, Scheel M A and Teukolsky S A 2008
Phys.Rev.D Preprint )[46] Foucart F, Deaton M B, Duez M D, Kidder L E, MacDonald I, Ott C D, Pfeiffer H P, ScheelM A, Szilagyi B and Teukolsky S A 2013
Phys.Rev.D Preprint )[47] Thorne K S 1980
Rev. Mod. Phys. Progress of Theoretical Physics
Preprint )[49] Lattimer J M and Swesty F D 1991
Nucl. Phys.
A535
ApJS
J. Par. Comp.
512 – 522[52] Fern´andez R and Metzger B D 2013
ApJ
MNRAS
A&A ApJS
ApJ
ApJ
82 (
Preprint )[58] Cyburt R H, Amthor A M, Ferguson R, Meisel Z, Smith K, Warren S, Heger A, Hoffman R D,Rauscher T, Sakharuk A, Schatz H, Thielemann F K and Wiescher M 2010
ApJS
Physical Review A&A
A61 (
Preprint )[61] Mamdouh A, Pearson J M, Rayet M and Tondeur F 2001
Nuclear Physics A
Preprint nucl-th/0010093 )[62] Wahl A C 2002
Systematics of Fission-Product Yields, Technical Report LA-13928 (Los Alamos,NM: Los Alamos National Laboratory)[63] Fuller G M, Fowler W A and Newman M J 1982
ApJS Atomic Data and Nuclear Data Tables Nuclear Physics A
Preprint nucl-th/0001018 )[66] M¨oller P, Nix J R, Myers W D and Swiatecki W J 1995
Atomic Data and Nuclear Data Tables
185 (
Preprint nucl-th/9308022 )[67] Kasen D, Thomas R C and Nugent P 2006
ApJ
Preprint astro-ph/0606111 )[68] Kyutoku K, Ioka K, Okawa H, Shibata M and Taniguchi K 2015
PRD Preprint )[69] Foucart F 2012
Phys.Rev.D Preprint )[70] Metzger B D, Piro A L and Quataert E 2008
MNRAS
Preprint )[71] Fern´andez R, Kasen D, Metzger B D and Quataert E 2015
MNRAS
Phys.Rev.D Preprint ) H-NS merger ejecta [73] Foucart F, Duez M D, Kidder L E, Scheel M A, Szilagyi B and Teukolsky S A 2012 Phys.Rev.D Preprint )[74] Goriely S 1999
A&A
MNRAS
Preprint )[76] Kyutoku K, Ioka K and Shibata M 2013
PRD Preprint )[77] Tanaka M, Hotokezaka K, Kyutoku K, Wanajo S, Kiuchi K, Sekiguchi Y and Shibata M 2014
ApJ
31 (
Preprint )[78] Kawaguchi K, Kyutoku K, Shibata M and Tanaka M 2016
ApJ
52 (
Preprint )[79] Rees M J 1988
Nature
MNRAS
L48–L51 (
Preprint arXiv:astro-ph/0611440 )[81] Metzger B D, Arcones A, Quataert E and Mart´ınez-Pinedo G 2010
MNRAS
Preprint )[82] Fong W, Berger E, Metzger B D, Margutti R, Chornock R, Migliori G, Foley R J, Zauderer B A,Lunnan R, Laskar T, Desch S J, Meech K J, Sonnett S, Dickey C, Hedlund A and Harding P2014
ApJ
118 (
Preprint )[83] Kisaka S, Ioka K and Nakar E 2016
ApJ
104 (
Preprint )[84] Dietrich T, Bernuzzi S, Ujevic M and Tichy W 2017
PRD Preprint )[85] Paschalidis V, East W E, Pretorius F and Shapiro S L 2015
PRD Preprint )[86] Morozova V, Piro A L, Renzo M, Ott C D, Clausen D, Couch S M, Ellis J and Roberts L F2015
ApJ
63 (
Preprint1505.06746