Signatures of hypermassive neutron star lifetimes on r-process nucleosynthesis in the disc ejecta from neutron star mergers

Abstract

We investigate the nucleosynthesis of heavy elements in the winds ejected by accretion discs formed in neutron star mergers. We compute the element formation in disc outflows from hypermassive neutron star (HMNS) remnants of variable lifetime, including the effect of angular momentum transport in the disc evolution. We employ long-term axisymmetric hydrodynamic disc simulations to model the ejecta, and compute r-process nucleosynthesis with tracer particles using a nuclear reaction network containing ∼8000 species. We find that the previously known strong correlation between HMNS lifetime, ejected mass and average electron fraction in the outflow is directly related to the amount of neutrino irradiation on the disc, which dominates mass ejection at early times in the form of a neutrino-driven wind. Production of lanthanides and actinides saturates at short HMNS lifetimes (≲10 ms), with additional ejecta contributing to a blue optical kilonova component for longer-lived HMNSs. We find good agreement between the abundances from the disc outflow alone and the solar r-process distribution only for short HMNS lifetimes (≲10 ms). For longer lifetimes, the rare-earth and third r-process peaks are significantly underproduced compared to the solar pattern, requiring additional contributions from the dynamical ejecta. The nucleosynthesis signature from a spinning black hole (BH) can only overlap with that from an HMNS of moderate lifetime (≲60 ms). Finally, we show that angular momentum transport not only contributes with a late-time outflow component, but that it also enhances the neutrino-driven component by moving material to shallower regions of the gravitational potential, in addition to providing additional heating.

accretion, accretion discs, dense matter, gravitational waves, hydrodynamics, neutrinos, nuclear reactions, nucleosynthesis, abundances

1 INTRODUCTION

The astrophysical origin of chemical elements formed through the rapid neutron capture process (r-process) remains an open problem in nuclear astrophysics. Observed abundances in metal-poor Galactic halo stars demand a mechanism that produced a robust abundance pattern – mirroring that in the Solar System – for elements with mass number A ≳ 130 (e.g. Sneden, Cowan & Gallino 2008). This mechanism must also have operated since early on in cosmic history, since r-process elements are found in very old metal-poor stars (e.g. Cowan et al. 1999; Ji et al. 2016). In contrast, observed abundances of light r-process elements in metal-poor stars show deviations from the Solar System pattern relative to heavier elements (e.g. Montes et al. 2007). Meteoritic abundances point to different formation time-scales for light and heavy r-process elements, suggesting that there might be more than one dominant formation site (Wasserburg, Busso & Gallino 1996).

Core-collapse supernovae may be able to produce light r-process elements (A ≲ 130), but recent work indicates that they are most likely not the dominant source of heavy r-process elements (e.g. Fischer et al. 2010; Hüdepohl et al. 2010; Roberts, Woosley & Hoffman 2010; Martínez-Pinedo et al. 2012; Wanajo 2013). On the other hand, mergers of binaries containing two neutron stars (NSNS) or a neutron star and a black hole (NSBH) have long been considered as candidate r-process sites, given the highly neutron-rich conditions achieved in the ejected material (Lattimer & Schramm 1974). The study of NSNS/NSBH mergers has intensified in recent years given that (1) they are likely to be detected in gravitational waves by the Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) and the Virgo interferometer within the next few years (e.g. Abadie et al. 2010; LIGO Scientific Collaboration et al. 2015), (2) recent developments in numerical relativity have enabled merger simulations consistent with Einstein's equations of general relativity (e.g. Lehner & Pretorius 2014; Paschalidis 2017) and (3) the electromagnetic signal from these events can aid with the localization of these sources and provide information complementary to that carried by gravitational waves (e.g. Rosswog 2015; Fernández & Metzger 2016; Tanaka 2016).

Recent work has shown that the dynamical ejecta from NSNS/NSBH mergers can produce a robust solar abundance pattern for A ≳ 130 by virtue of fission cycles (e.g. Goriely et al. 2005), with little sensitivity to binary parameters or the equation of state (e.g. Goriely, Bauswein & Janka 2011; Korobkin et al. 2012; Bauswein, Goriely & Janka 2013), depending instead on nuclear physics properties such as the fission fragment distribution (e.g. Eichler et al. 2015). When neutrino absorption is included in the calculations, a larger fraction of lighter elements (A < 130) can be obtained (e.g. Wanajo et al. 2014; Goriely et al. 2015; Sekiguchi et al. 2015; Radice et al. 2016; Foucart et al. 2016a,b; Roberts et al. 2017).

In addition to the dynamical ejecta, the accretion disc formed in NSNS/NSBH mergers can eject a significant amount of material on time-scales longer than the dynamical time (e.g. Ruffert et al. 1997; Lee, Ramirez-Ruiz & López-Cámara 2009; Metzger, Piro & Quataert 2009b). The neutron-to-seed ratio of this material is lower than that in the dynamical ejecta, because the longer evolutionary time-scales allow weak interactions to modify the composition more significantly (e.g. Dessart et al. 2009; Fernández & Metzger 2013; Perego et al. 2014). Matter can be ejected on the thermal time-scale of the disc (∼30 ms, Section 2.3) by neutrino energy deposition, or on much longer time-scales (∼1 s, Section 2.3) by a combination of angular momentum transport processes and nuclear recombination. The amount of mass ejected by the disc on the longer time-scale can be comparable to that in the dynamical ejecta (e.g. Fernández & Metzger 2016), while the neutrino-driven component is significant only when a hypermassive neutron star (HMNS) phase precedes BH formation (Dessart et al. 2009; Fernández & Metzger 2013; Metzger & Fernández 2014; Just et al. 2015).

Early work on nucleosynthesis from NSNS/NSBH merger discs focused on neutrino-driven outflows and used a parametric treatment to obtain thermodynamic trajectories for composition analysis (e.g. Surman et al. 2008; Wanajo, Janka & Kubono 2011), finding that conditions for both light and heavy r-process elements are possible in these outflows. More recent work has employed tracer particles from long-term hydrodynamic simulations of the disc when a BH is the central object (Just et al. 2015; Wu et al. 2016). These studies have found that the outflow generates a robust abundance of elements around A = 130, with significant production of lighter elements, and a variable yield of heavy r-process elements that depends on binary properties and disc physics. The case of an HMNS at the centre was studied by Martin et al. (2015) using tracer particles from a time-dependent simulation that considered the neutrino-driven wind phase (Perego et al. 2014). The resulting outflow generates primarily nuclei with A < 130, with a significant dependence of the yield on latitude and ejection time. The long-term (viscous) outflow was not captured, however, and the lifetime of the HMNS was accounted for only by looking at subsets of particles that were ejected before a certain time.

In this study, we analyse the nucleosynthesis yields from the long-term outflow generated by an accretion disc around a HMNS of variable lifetime. Following the approach of Metzger & Fernández (2014), we conduct a number of long-term disc simulations in which the HMNS transforms into a BH at different times. This parametrized approach does not require the assumption of a particular equation of state of dense matter. It is also independent of the complex processes that determine transport of angular momentum and cooling, all of which set the survival time of the HMNS (e.g. Paschalidis, Etienne & Shapiro 2012; Kaplan et al. 2014). In order to obtain nucleosynthesis yields, tracer particles are injected in the disc initially, and the resulting thermodynamic trajectories are analysed with the nuclear reaction network SkyNet (Lippuner & Roberts 2017).

This paper is structured as follows. Section 2 describes the numerical method employed and the models studied. Results are presented and discussed in Section 3, and we conclude in Section 4. The Appendix contains further details about our numerical implementation.

2 METHODS

2.1 Disc outflow simulations and thermodynamic trajectories

The long-term evolution of the accretion disc is computed using the approach described in Metzger & Fernández (2014). In short, the equations of Newtonian hydrodynamics and lepton number conservation are solved with the flash code (Fryxell et al. 2000; Dubey et al. 2009) assuming axisymmetry (2D). Source terms include gravity via a pseudo-Newtonian potential, shear viscosity with an α prescription (Shakura & Sunyaev 1973), and contributions due to neutrino absorption and emission to the lepton number and energy equations. A pseudo-Newtonian potential, such as that from Artemova, Bjoernsson & Novikov (1996), approximates the effect of a Schwarzschild or Kerr metric by providing an innermost stable circular orbit (ISCO) in an otherwise Newtonian hydrodynamic simulation. Other aspects, such as the angular dependence of the Kerr metric, are not captured.

Neutrinos are implemented in a leakage scheme for cooling and a disc light-bulb approximation for self-irradiation. Only charged-current interactions are included, since they have the largest cross-section, exchange energy with matter, and drive the evolution of the electron fraction. The HMNS is approximated as a reflecting sphere with a radius of 30 km and a rotation period of 1.5 ms. A time-dependent, isotropic outward neutrino flux is imposed on the HMNS surface with a constant value up to 10 ms and a time dependence t^−1/2 (e.g. Pons et al. 1999) thereafter, with t being the physical simulation time. The electron neutrino and antineutrino luminosities have the same magnitude, which we normalize to 2 × 10⁵² erg s⁻¹ at 30 ms, roughly matching the results of Dessart et al. (2009):

\begin{equation} L_\nu = L_{\bar{\nu }} = 2\sqrt{3}\times 10^{52}\left[\max \left(1,\frac{t}{10\,\rm{ms}}\right)\right]^{-1/2} \rm{erg\,\,s}^{-1}. \end{equation}

(1)

The neutrino and antineutrino temperatures are different, however, chosen to be 4 and 5 MeV, respectively, which roughly corresponds to the values found in proto-neutron stars in core-collapse supernovae (e.g. Janka 2001). Thus, the mean neutrino energies from the HMNS are fixed at 〈E_ν〉 = 12.6 MeV and |$\langle E_{\bar{\nu }}\rangle = 15.8$| MeV (assuming a Fermi–Dirac distribution of neutrinos, the mean energy is given by 〈E_ν〉 = 3.151 k_B T_ν, where k_B is the Boltzmann constant and T_ν is the neutrino temperature). These mean neutrino energies from the HMNS are broadly consistent with results from numerical relativity simulations that include neutrino transport (Foucart et al. 2016a). The HMNS is transformed into a BH by switching the inner radial boundary from reflecting to absorbing and setting the HMNS neutrino luminosities to zero. The viscous stress is also set to zero at this boundary when the BH forms. See Metzger & Fernández (2014) for more details.

The initial condition for most models is an equilibrium torus obtained by solving the Bernoulli equation with constant specific angular momentum and electron fraction Y_e = 0.1 (e.g. Papaloizou & Pringle 1984; Fernández & Metzger 2013). Outside the torus, the computational domain is filled with a low-density ambient medium that follows a power law in radius, with an initial normalization ∼9 orders of magnitude below the maximum torus density (ρ_max ∼ 10¹¹ g cm⁻³). The density floor inside 200 km decreases with time towards a constant asymptotic value of 10 g cm⁻³ over a time-scale of ∼100 ms. In addition, one model initializes the disc using a snapshot from a general-relativistic NSBH simulation reported in Foucart et al. (2015). Details about the mapping procedure can be found in Fernández et al. (2017). Improvements to the code relative to Metzger & Fernández (2014) include the use of separate neutrino temperatures for disc self-irradiation, and a correction to the weak interaction rates; details are provided in the Appendix.

Passive tracer particles record thermodynamic quantities as a function of time, and the resulting information is used as input for the nuclear reaction network calculations (Section 2.2). The initial particle locations are randomly sampled to follow the mass distribution in the disc. We place 10 000 particles initially in all simulations. If there is a BH at the centre, particles can fall into it, which then reduces the total number of tracer particles. Fluid quantities are obtained at each time step from the grid by linear interpolation. When a HMNS is at the centre, the reflecting boundary for particles is placed one cell outside the inner radial boundary to prevent particle trapping in a ‘trench’ of small-magnitude negative radial velocity that forms in the active cells adjacent to this boundary. The small size of this innermost radial cell relative to the inner boundary radius (Δr/r ≃ 1.8 per cent) ensures that the effect on the particle dynamics is minimal considering all other approximations being made.

Most hydrodynamic simulations are evolved up to about 10 s of physical time (Section 2.3), which is sufficient for r-process nucleosynthesis that takes place within a few seconds. To obtain the final abundances, however, the nuclear reaction network needs to be evolved for tens of years, since some of the isotopes produced by the r-process have very long half-lives. One caveat with this post-processing approach is that the energy released by the nuclear reactions does not feed back into the hydrodynamic evolution of the fluid (with the exception of α-recombination, which is accounted for in the hydro evolution). Nuclear heating may slightly change the morphology of the ejecta at late times (see, e.g. fig. 8 of Fernández et al. 2015b) and influence specific features of the nucleosynthesis indirectly via the amount of convection in the outflow (Wu et al. 2016, see also Section 3.4.1). But the majority of heating due to the r-process and nuclear decays occurs after particles have been ejected from the disc, so the dynamics is not significantly altered by not including this heating in the hydrodynamic simulation.

2.2 Nuclear reaction network: skynet

We employ the nuclear reaction network SkyNet for the r-process nucleosynthesis calculations (Lippuner & Roberts 2017). For each thermodynamic trajectory (Section 2.1), we begin the reaction network evolution once the temperature falls below 10 GK or reaches its maximum, if the maximum is less than 10 GK. Since the composition is given by nuclear statistical equilibrium (NSE) at temperatures above 10 GK, there is no need to start the reaction network evolution at higher temperatures.

SkyNet includes the specific viscous heating and neutrino heating/cooling rates recorded by the thermodynamic trajectories, and adds to these a self-consistent calculation of heating from nuclear reactions. In addition to the neutrino heating/cooling rate, the associated rates of neutrino interactions with free neutrons and protons are also given by the thermodynamic trajectories and evolved in SkyNet. The included reactions are electron (anti) neutrino emission (p + e⁻ → n + ν_e, |$\text{n} + e^+ \rightarrow \text{p} + \bar{\nu }_e$|⁠) and neutrino absorption (n + ν_e → p + e⁻, |$\text{p} + \bar{\nu }_e \rightarrow \text{n} + e^+$|⁠).

SkyNet evolves the abundances of 7843 nuclides, ranging from free neutrons and protons to ³³⁷Cn, and includes over 140,000 nuclear reactions. The strong reaction rates are taken from the REACLIB data base provided by the Joint Institute for Nuclear Astrophysics (JINA) (Cyburt et al. 2010), but only the forward rates are used and the inverse rates are computed from detailed balance. Spontaneous and neutron-induced fission rates are taken from Frankel & Metropolis (1947), Mamdouh et al. (2001), Wahl (2002) and Panov et al. (2010). Most of the weak rates come from Fuller, Fowler & Newman (1982), Oda et al. (1994) and Langanke & Martínez-Pinedo (2000) whenever they are available, and otherwise the REACLIB weak rates are used. The nuclear masses and partition functions used in SkyNet are taken from the WebNucleo XML file distributed with REACLIB, which contains experimental data where available and finite-range droplet macroscopic model (FRDM; see e.g. Möller et al. 2016) data otherwise.

2.3 Investigated models

Table 1 summarizes the properties of all investigated models. The main focus of this study is the outflow from a disc around an HMNS of variable lifetime. Our baseline sequence follows the parameter choices of Metzger & Fernández (2014): an HMNS mass of 3 M_⊙ arising from an NSNS merger, with an initial disc mass M_d = 0.03 M_⊙ chosen as a representative case of disc masses obtained in NSNS mergers (e.g. Hotokezaka et al. 2013). We prescribe the lifetime τ of the HMNS to be 0, 10, 30, 100 or 300 ms, after which the HMNS collapses to a non-spinning BH. We also run one case, denoted by τ = ∞, in which the HMNS does not collapse. The HMNS models are denoted by H000, H010, H030, H100, H300 and Hinf, according to their lifetime. Other disc parameters are: density peak radius R_d = 50 km, viscosity parameter α = 0.03, constant initial entropy of 8 k_B baryon⁻¹, constant initial Y_e = 0.1 and maximum evolution time of 8.7 s. These choices are motivated to be broadly compatible with results from dynamical merger simulations (e.g. Ruffert et al. 1997; Oechslin & Janka 2006; Foucart et al. 2016a) and from studies of angular momentum transport in fully-ionized accretion discs (e.g. Davis, Stone & Pessah 2010).

Table 1.

List of investigated models. Columns from left to right show the model name, the compact central object (CCO) type (HMNS or BH), mass M_c of the CCO, lifetime τ of the HMNS, dimensionless spin χ of the BH (the HMNSs all spin at 1.5 ms), radius R_d of the initial disc density peak and viscosity parameter α.

Model	CCO	Mc	τ	χ	Rd	α
		(M⊙)	(ms)		(km)
H000	HMNS	3	0	0	50	0.03
H010			10
H030			30
H100			100
H300			300
Hinf			∞
B070	BH	3	0	0.7	50	0.03
B090				0.9
BF15	BH	8.1	0	0.86	55	0.03
HinfNoVisc	HMNS	3	∞	0	50	0

There are three important time-scales in the problem: the orbital time at the initial disc density peak,

\begin{equation} t_{\rm orb}\simeq 3\left(\frac{R_{\rm d}}{50\rm{\,\,km}} \right)^{3/2} \left(\frac{3\,\,\mathrm{M}_{\odot }}{M_{\rm c}}\right)^{1/2}\rm{ms}, \end{equation}

(2)

where M_c is the mass of the central object (HMNS or BH); the initial thermal time in the disc,

\begin{equation} t_{\rm th}\simeq 30\left(\frac{M_{\rm d}}{0.03\,\,\mathrm{M}_{\odot }}\right)\left(\frac{e_{{\rm i},{\rm d}}}{10^{19}\rm{\,erg\,\,g}^{-1}}\right) \left(\frac{2}{L_{\nu ,52}} \right)\rm{ms}, \end{equation}

(3)

where e_i,d is the initial specific internal energy of the disc (a byproduct of shock heating during the dynamical phase of the merger) and L_ν,52 is a typical neutrino luminosity from the disc in units of 10⁵² erg s⁻¹ (neutrino cooling is in approximate balance with viscous heating at early times); and the initial viscous time of the disc,

\begin{eqnarray} t_{\rm visc}\simeq 200\left(\frac{0.03}{\alpha }\right)\left(\frac{0.3}{H/R}\right)^2 \left(\frac{R_{\rm d}}{50\rm{\,\,km}}\right)^{3/2}\left(\frac{3\,\,\mathrm{M}_{\odot }}{M_{\rm c}}\right)^{1/2}\rm{ms}, \nonumber \\ \end{eqnarray}

(4)

where H/R is the height-to-radius ratio of the disc. Outflows driven primarily by neutrino energy deposition are expected to be launched on the thermal time-scale (equation 3), whereas long-term outflows are launched on the viscous time-scale (equation 4). The latter becomes longer with time, as the disc spreads out and most of the mass in the disc resides at an increasingly larger radius (e.g. Metzger, Piro & Quataert 2009a). Therefore, the length of the simulations has to be several viscous time-scales, which translates into thousands of orbits. We define convergence in mass ejection from the disc as a saturation in the cumulative mass crossing some radius far away from the disc (10⁹ cm). This way, we generally require 3000 t_orb, or about ∼10 s for each simulation.

In order to examine the impact of the central compact object on the nucleosynthesis, we evolve two additional models with spinning BHs at the centre, following the approach of Fernández et al. (2015a). The BH mass is M_c = 3 M_⊙ in both cases and the dimensionless spin parameter |$\chi = Jc/(GM_{\rm c}^2)$| is 0.7 or 0.9, where J is the BH angular momentum, c is the speed of light and G is the gravitational constant. These models are labelled B070 and B090, according to their spin. Other parameters are the same as in the HMNS models. In addition, we investigate the effect of the disc compactness (mass of the central compact object divided by the disc density peak radius, M_c/R_d, which measures the strength of the gravitational field) by evolving a model in which the initial condition is taken from a snapshot of a general-relativistic simulation of an NSBH merger from Foucart et al. (2015), including only the remnant accretion disc. This model has a BH mass of 8.1 M_⊙ and a disc density peak radius of 55 km, resulting in a gravitational potential that is 2.5 times stronger than our fiducial case (which better approximates the remnant of an NSNS merger). The model is denoted by BF15, and the mapping details are described in Fernández et al. (2017).

Finally, we evolve a test model that mirrors the case of an HMNS with τ = ∞, but with the viscosity parameter set to zero, in order to eliminate angular momentum transport and viscous heating. This model, denoted by HinfNoVisc, experiences an outflow driven solely by neutrino heating (possibly aided by nuclear recombination), and is evolved in order to compare results with Martin et al. (2015). The model is evolved for a longer time (14.5 s) since mass ejection converges more slowly with time relative to the viscous case.

3 RESULTS AND DISCUSSION

3.1 Overview of disc evolution

The evolution of the disc and especially the neutrino interactions occurring in the disc set the stage for the outflow and determine its properties. In this section, we present a brief summary of the disc evolution, which was described in detail in Fernández & Metzger (2013) and Metzger & Fernández (2014).

When an HMNS is present, transport of angular momentum causes accreting material to form a boundary layer around the reflecting stellar surface. The outer regions of the disc expand on a thermal time-scale (equation 3) due to energy injection by neutrino heating and viscous heating. Upon collapse of the HMNS to a BH, the boundary layer is swallowed, and a rarefaction wave moves outward from the inner boundary, quenching the thermal outflow (cf. fig. 3 of Metzger & Fernández 2014). The disc re-adjusts on a thermal time-scale, and joins the evolutionary path of a BH accretion disc, which changes on a viscous time-scale (equation 4).

The high densities achieved on the equatorial plane of HMNS discs (∼10¹¹ g cm⁻³) are enough to locally trap neutrinos, resulting in two emission hotspots at mid-latitudes and adjacent to the HMNS surface. The high densities in the mid-plane also result in shadowing of the outer disc, with neutrino irradiation concentrated at latitudes ∼30° away from the equator (cf. fig. 2c of Metzger & Fernández 2014). This general neutrino irradiation geometry is also found when better (Monte Carlo) radiation transport is performed on snapshots of HMNS disc models (Richers et al. 2015).

The high densities near the boundary layer cause the electron fraction in the inner regions of the HMNS disc (r < 10⁷ cm, and within ∼30° from the equatorial plane) to remain low (Y_e ∼ 0.1–0.2) relative to the outer regions. Away from the mid-plane, the weak interaction time-scale becomes shorter, and material that enters the boundary layer at high latitude reaches Y_e ∼ 0.5 within an orbital time-scale, as it transits through the hotspot in neutrino emission. Material ejected within ∼45° of the rotation axis has therefore high electron fraction relative to the rest of the outflow. At t = t_th, most of the disc material within r = 10⁷ cm is close to beta equilibrium, with Y_e ∼ 0.2 − 0.3 in the outer regions. As the disc approaches a viscous time-scale, the density in the inner disc gradually decreases, resulting in two effects: (1) the equilibrium Y_e of the inner disc increases as the degeneracy of the material decreases, and (2) the strength of the weak interactions decreases, until weak interactions become slow relative to the viscous time (e.g. Metzger et al. 2009a). These two effects combine to leave material close to the HMNS with Y_e ∼ 0.4. As long as the HMNS does not collapse, most of this material will eventually be ejected.

The BH accretion disc is such that the inner regions (r < 10⁷ cm) also approach beta equilibrium, but most of the material that reaches high Y_e is accreted on to the BH (cf. fig. 6 of Fernández & Metzger 2013). The fraction of the high-Y_e material that is ejected is a function of the spin of the BH, since a smaller ISCO slows down accretion and allows material to reach regions where weak interactions are stronger, before being ejected (Fernández et al. 2015a). Upon collapse of the HMNS into a BH, a low density funnel of width ∼45° around the rotation axis is created as material that has not yet been unbound is swallowed. Boundary layer material reaches the highest Y_e as it crosses the emission hotspots of the HMNS disc before being ejected, hence formation of the BH precludes further ejection of material with Y_e ∼ 0.5. For the disc masses simulated, the material becomes optically thin within a few orbits if a BH is at the centre, therefore hotspots also disappear upon collapse of the HMNS.

3.2 Neutrino emission

The total electron neutrino and antineutrino luminosities from the disc and the HMNS (if present) as a function of time are shown in the left-hand panel of Fig. 1. The disc luminosities are very similar to the luminosity from the HMNS as long as the latter does not collapse, because the disc absorbs and re-radiates the neutrinos emitted by the HMNS. When the HMNS collapses, the disc luminosities undergo a sudden decrease, and then follow a broken power law as the disc re-adjusts and continues to evolve on the viscous time-scale. Thereafter, |$L_{\bar{\nu }_e}$| starts to dominate over |$L_{\nu _e}$| in all models because the initial neutron-rich (Y_e = 0.1) composition of the disc causes weak interactions to leptonize it, with more positron captures than electron captures (the gradual drop in density over the viscous time-scale increases the equilibrium Y_e of the disc; Section 3.1). As long as the HMNS is present, the disc neutrino luminosities are approximately the same in all cases, consistent with re-radiation of the HMNS luminosities, and the different models are separated by a shift in time depending on the lifetime of the HMNS.

Figure 1.

Left: Neutrino luminosities as a function of time for most models studied in this paper (see Table 1 for a summary). Shown are electron neutrinos (solid lines) and electron antineutrinos (dashed lines) emitted by the disc, as well as the imposed neutrino/antineutrino emission from the HMNS surface (dotted lines, equation 1). Sharp drops in the dotted lines mark the collapse of the HMNS into a BH. Right: Mean energies of electron neutrinos (solid lines) and electron antineutrinos (dashed lines) emitted by the disc. The mean energies of neutrinos and antineutrinos emitted by the HMNS are fixed at 12.6 and 15.8 MeV, respectively.

The mean neutrino energies emitted by the disc are shown in the right-hand panel of Fig. 1. Given the low initial abundance of protons, the optical depth for antineutrinos is initially low, and antineutrino energies are higher by ∼10 per cent than neutrino energies until t ∼ t_visc, by which time the decrease in disc density has gradually lifted the degeneracy and weak interactions have driven the inner disc close to Y_e ∼ 0.4. The neutrino/antineutrino energies are the same between HMNS models until the HMNS collapses, each following the same general evolution pattern but shifted in time. Since the disc mostly re-radiates neutrinos from the HMNS, the mean neutrino energies drop sharply when the HMNS collapses.

In the models that start out with a BH at the centre (B070, B090, BF15 and H000, which has χ = 0), a more rapidly spinning BH results in higher neutrino luminosities and mean neutrino energies. This occurs because larger spins are associated with smaller ISCO radii. Hence the disc material can convert more gravitational energy into thermal energy, resulting in more intense neutrino emission with higher mean energies. In model BF15, the central BH is more massive (8.1 M_⊙ compared to 3 M_⊙ in H000, B070 and B090), which results in a larger ISCO radius. The disc has nearly the same density peak radius as the other models, but the initial condition is not in equilibrium. Therefore, accretion proceeds more intensely at early times than in the other BH models, speeding up the disc evolution despite its slightly longer initial viscous time (smaller H/R in equation 4, which overcomes the effect of a larger BH mass).

3.3 Ejecta properties

In order to associate the thermodynamic properties of the disc with the nucleosynthesis outcome from each trajectory, we use the value of the electron fraction (Y_e,5GK) and specific entropy (s_5GK) at the last time when the temperature of the tracer particle drops below 5 GK. For the rare cases in which the temperature of the trajectory is always below 5 GK, we use the initial values of Y_e and s for Y_e,5GK and s_5GK, respectively. Once the temperature drops below approximately 5 GK, the composition moves out of NSE and a full network evolution is required to evolve the abundances. Therefore Y_e,5GK and s_5GK are the initial conditions for nucleosynthesis. Note that the reaction network evolution starts when the temperature drops below 10 GK, but SkyNet can also evolve the composition while NSE holds.

The distribution of mass ejected as a function of Y_e,5GK is shown in Fig. 2 for all HMNS models with non-zero viscosity. A trajectory is considered to have been ejected when it crosses the surface r = 10⁹ cm. There is a strong correlation between the HMNS lifetime and both the amount of mass ejected and the mean Y_e,5GK of the distribution (Metzger & Fernández 2014). The disc ejecta ranges from 6 to nearly 100 per cent of the initial disc mass. The longer the HMNS lives, the longer the disc material is subject to strong neutrino heating, which combines with viscous heating and nuclear recombination to eject material on the viscous time-scale. A longer HMNS lifetime also allows more material from the inner disc to be ejected instead of being swallowed by the BH. That material from the inner disc reaches beta equilibrium and hence its ejection results in a higher mean electron fraction. Table 2 shows the number of ejected particles (out of the initial 10 000 particles in each model) and the total ejected mass. Also shown is the amount of mass ejected with Y_e,5GK ≤ 0.25, which is neutron-rich enough to robustly make the full r-process (see the next section and e.g. Lippuner & Roberts 2015).

Figure 2.

Mass distribution of the ejecta electron fraction for HMNS models at the time when the temperature is ≥5 GK for the last time for each tracer particle.

Table 2.

Summary of nucleosynthesis results. N_ej is the number of tracer particles that reach a radius r = 10⁹ cm by the end of the hydrodynamic simulation (every simulation starts with 10 000 particles); M_ej is the total ejected mass in 10⁻³ M_⊙ at the same radius; M(Y_e ≤ 0.25) is the ejected mass with Y_e,5GK ≤ 0.25; M_ν-driv is the amount of ejected mass that is driven by neutrino interactions; [Y_1st/Y_2nd] is the log₁₀ of the ratio of the first r-process peak to the second r-process peak, normalized to the solar value, see equation (5) for details; [Y_RE/Y_2nd] and [Y_3rd/Y_2nd] are the same quantities for the rare-earth and third peaks; 〈X_La〉 and 〈X_Ac〉 are the lanthanide and actinide mass fractions averaged over all ejecta particles; and ε_tot 1 d and ε_tot 7 d are the total heating rates of the entire ejecta at 1 and 7 d, respectively.

Model	Nej	Mej	M(Ye ≤ 0.25)	Mν-driv	[Y1st/Y2nd]	[YRE/Y2nd]	[Y3rd/Y2nd]	〈XLa〉	〈XAc〉	εtot 1 d	εtot 7 d
		(M−3⊙)	(M−3⊙)	(M−3⊙)						(erg s−1)	(erg s−1)
H000	527	1.8	1.4	0.013	−1.3	−0.28	−0.30	4.6 × 10−2	6.4 × 10−3	9.0 × 1040	1.4 × 1040
H010	557	1.9	1.1	0.027	−1.1	−0.40	−0.50	3.3 × 10−2	2.0 × 10−3	8.9 × 1040	1.3 × 1040
H030	989	3.3	0.83	0.20	−0.64	−1.0	−1.2	5.1 × 10−3	2.9 × 10−4	1.2 × 1041	1.5 × 1040
H100	2408	7.8	0.52	1.3	−0.0053	−0.91	−1.2	2.1 × 10−3	1.7 × 10−4	1.8 × 1041	1.4 × 1040
H300	5610	18	0.67	6.4	+0.25	−0.88	−1.2	1.1 × 10−3	6.7 × 10−5	3.3 × 1041	3.2 × 1040
Hinf	9587	30	0.69	28*	+0.41	−0.86	−1.1	7.1 × 10−4	4.2 × 10−5	5.2 × 1041	5.4 × 1040
B070	1465	5.4	1.8	0.022	−0.73	−0.65	−0.67	1.3 × 10−2	1.6 × 10−3	2.0 × 1041	2.4 × 1040
B090	2363	7.9	1.6	0.070	−0.54	−0.77	−0.80	7.6 × 10−3	9.7 × 10−4	2.7 × 1041	2.6 × 1040
BF15	910	4.9	0.011	0.022	−0.26	−1.4	−1.3	1.4 × 10−3	7.7 × 10−5	7.8 × 1040	6.7 × 1039

*Since the HMNS persists forever in model Hinf, virtually all trajectories experience significant neutrino interactions and our method becomes inadequate to isolate the component of the ejecta that is driven by neutrino interactions.

We note that the amount of mass ejected with Y_e,5GK ≤ 0.25 is roughly constant between (0.5–0.8) × 10⁻³ M_⊙ once the HMNS lives for 30 ms or longer, despite the total ejecta mass differing by an order of magnitude. This is the result of two competing effects: a longer HMNS lifetime increases the total ejecta mass, but it also increases the average electron fraction of the ejecta, thus reducing the fraction of the ejected mass that has Y_e,5GK ≤ 0.25. These two effects counteract each other, leaving the absolute amount of mass ejected with Y_e,5GK ≤ 0.25 roughly constant.

The thermodynamic properties of the ejecta for model H300 (HMNS with lifetime τ = 300 ms) are illustrated in Fig. 3 through the electron fraction Y_e,5GK, specific entropy s_5GK, final velocity v_final, ejecta mass, and the time t_5GK when the temperature is 5 GK for the last time. Two ejecta components stand out from the scatter plot of Y_e,5GK versus s_5GK. The larger component is ejected before the HMNS collapses, i.e. t_5GK ≲ 300 ms, and it exhibits a tight correlation between the entropy and electron fraction up to Y_e,5GK ∼ 0.5. This is indicative of a neutrino-driven wind, and indeed we would expect the asymptotic Y_e to be ∼0.55 based on the neutrino properties shown in Fig. 1 (Qian & Woosley 1996). Note that the vast majority of the ejecta with Y_e,5GK ≳ 0.4 or v_final ≳ 0.03 c is part of this early wind-like ejecta, with a much smaller group of particles extending to low velocities and low electron fractions.

Figure 3.

Properties of the ejecta of model H300 (central HMNS that collapses after 300 ms). Y_e,5GK and s_5GK are the electron fraction and specific entropy at the time t_5GK, when the temperature of the trajectory is 5 GK for the last time, while v_final is the final velocity of the trajectory. The scatter plot points are colour-coded by t_5GK, which marks the nucleosynthesis starting time. Blue and grey points (t_5GK ≲ 300 ms) start their nucleosynthesis while the HMNS is present, and are thus influenced by the neutrino irradiation from the HMNS. The strong correlation between electron fraction and entropy is similar to what is obtained in a neutrino-driven outflow. The red points (t_5GK ≳ 300 ms) start nucleosynthesis after the HMNS has collapsed to a BH. Hence this component is subject primarily to the action of viscous processes and nuclear recombination in the disc. See the text for details. M_ej,−3⊙ means the amount of ejecta mass in units of 10⁻³ M_⊙.

The second component is ejected after the HMNS has collapsed to a BH, i.e. t_5GK ≳ 300 ms, and it has only a weak correlation between Y_e,5GK and s_5GK. This component is associated with mass ejection as the disc reaches the advective state (very weak or no neutrino cooling/heating) and is driven primarily by heat injection from angular momentum transport processes and nuclear recombination (Metzger & Fernández 2014). Mass ejection in this state is accompanied by vigorous convective activity in the disc.

In order to quantitatively disentangle the wind-like ejecta component seen in Fig. 3 from the other component, we compute the contributions to s_5GK arising from neutrino and viscous heating. While all trajectories experience some degree of viscous heating, only a subset of tracer particles experience an entropy change due to neutrino heating that is larger than 0.1 k_B baryon⁻¹ (ignoring any neutrino cooling). And that subset exactly exhibits the tight, boomerang-shaped correlation between Y_e,5GK and s_5GK that the early, wind-like ejecta exhibits (cf. the blue and grey dots in the lower left panel of Fig. 3). The mass of this neutrino-driven ejecta component is shown in Table 2 under M_ν-driv. This component is absent for HMNS lifetimes shorter than 30 ms and also in the BH models, but once the HMNS lifetime becomes longer, the neutrino-driven component becomes as large or even the dominant fraction of the total ejecta (such as in model Hinf, in which almost the entire ejecta is neutrino-driven). For very long-lived HMNSs all trajectories experience some degree of neutrino interactions, and it becomes difficult to distinguish the wind-like from the viscous component. We emphasize that even the neutrino-driven component experiences significant viscous heating and the entropy change due to viscous heating can be of the same order as that due to neutrino absorption.

One interesting feature of the late-time component is the sharp upper limit of Y_e,5GK ∼ 0.38–0.40, regardless of entropy. These trajectories experience late-time fall back due to large convective eddies in the disc. They get sucked deep inside the disc where the density is much higher than in the outflow, and then they are ejected again almost immediately. This creates a spike in their density profile that results in significant heating, as evidenced by the fact that they all have T ≥ 5GK at t ∼ 2 s. However, before this late-time heating occurs, r-process nucleosynthesis has already taken place in these trajectories, and all free neutrons have been captured on to seed nuclei. Thus, the composition before the heating spike consists of heavy elements with β-decay half-lives of milliseconds to seconds. These elements decay and raise the overall electron fraction of the material to Y_e ∼ 0.38–0.40, which is the characteristic Y_e at 1–3 seconds after neutron exhaustion for the r-process, for a wide range of initial Y_e. The late-time heating then simply pushes the material back into NSE, but the electron fraction remains unchanged. The resulting entropy depends on the amount of heating received by each trajectory, as determined by how far the material falls back into the disc. This class of trajectories therefore ends up with electron fractions Y_e,5GK ∼ 0.38 − 0.40 and nucleosynthesis start times of t_5GK ∼ 2 s, with uncorrelated entropies.

3.4 Nucleosynthesis

3.4.1 Final abundances

The mass-averaged composition of the ejecta for all models with non-zero viscosity is shown in Fig. 4. The abundances are multiplied by the total ejecta mass to emphasize their relative contributions to the different r-process regions. Models H000 and H010 (prompt non-spinning BH and shortest-lived HMNS, respectively) agree most closely with the Solar System r-process abundances (Arnould et al. 2007), which have been scaled to match the second peak at A = 130 (the abundances from our models have not been scaled). The abundances around the third r-process peak in these two models approach the solar values, whereas in all other models production of the third peak is too low compared to solar. H000 and H010 also have the best agreement with the solar rare-earth peak around A ∼ 165. While these two models underproduce the first r-process peak (A ∼ 80), they agree rather well with the feature around A ∼ 100, in contrast to all other models which overproduce it.

Figure 4.

Final trajectory-averaged abundances as a function of mass number, scaled by the total ejecta mass, for all models with non-zero viscosity. The observed solar r-process abundances (Arnould, Goriely & Takahashi 2007) are scaled to match the second peak of the HMNS models at A = 130 (none of the abundances from our models have been scaled).

While the good agreement between models H000/H010 and the solar r-process abundances could be taken as an indication of short HMNSs lifetimes being more common, one has to keep in mind that Fig. 4 assumes that the entire second solar r-process peak is due to the disc outflow. Other sources such as the dynamical ejecta from NSNS/NSBH mergers and core-collapse supernovae can also produce significant amounts of r-process elements. The expected abundance patterns are weighted towards the third peak for the dynamical ejecta (e.g. Goriely et al. 2011; Wanajo et al. 2014; Roberts et al. 2017) and towards the first peak for core-collapse supernovae (e.g. Wanajo 2013; Shibagaki et al. 2016; Vlasov et al. 2017). The solar r-process abundance is thus the outcome of the contribution from each source weighted by their rate and yield per event.

In all models, the third peak is shifted to slightly higher mass numbers, which is a well-known shortcoming of the FRDM mass model (e.g. Mendoza-Temis et al. 2015; Mumpower et al. 2016). We also see an abundance spike at A = 132 in all models. This spike is due to some trajectories experiencing late-time heating that photodissociates neutrons from synthesized heavy elements. This results in additional neutron capture and a pile up of material at the doubly magic nucleus ¹³²Sn (N = 82 and Z = 50). Wu et al. (2016) also observed this phenomenon and described it in detail.

The models with longer HMNS lifetimes have less neutron-rich ejecta (Fig. 2) and hence synthesize a greater fraction of first peak material. Once the HMNS lifetime is longer than 100 ms, the first peak (70 ≤ A ≤ 90) is overproduced with respect to the solar values, when the abundances are normalized to the second peak. Again, we emphasize that the r-process yield from disc outflows is complementary to that from the dynamical ejecta, which tends to produce more neutron-rich nuclei.

We quantify the relative contribution of each model to the different regions of the r-process distribution by computing average abundances around the peaks and normalizing them to the solar values. The abundance of the second peak Y_2nd is computed as the sum of the abundances in the range 125 ≤ A ≤ 135, excluding A = 132 to avoid the spike at that mass number. For the first peak abundance Y_1st, we use the sum of abundances in the range 70 ≤ A ≤ 90. For the rare-earth peak Y_RE, we use 160 ≤ A ≤ 166 and for the third peak we use 186 ≤ A ≤ 203. The quantity [Y_1st/Y_2nd] shown in Table 2 is defined as

\begin{eqnarray} {\,\,}[Y_{\rm 1st}/Y_{\rm 2nd}] = \log _{10}\frac{Y_{{\rm 1st}}}{Y_{{\rm 2nd}}} - \log _{10}\frac{Y_{{\rm 1st},{\odot }}}{Y_{{\rm 2nd},{\odot }}}, \end{eqnarray}

(5)

where Y_{1st, ⊙} and Y_{2nd, ⊙} are the abundances of the third and second peak as observed in the Solar System, respectively. The same procedure is used to compute [Y_RE/Y_2nd] and [Y_1st/Y_2nd]. Using the solar r-process abundances from Arnould et al. (2007), we find log Y_{1st, ⊙}/Y_{2nd, ⊙} = +1.3, log Y_{RE, ⊙}/Y_{2nd, ⊙} = −1.1, and log Y_{3rd, ⊙}/Y_{2nd, ⊙} = −0.42, which we use to normalize the values shown in Table 2.

The different peak ratios shown in Table 2 quantify the trends apparent in Fig. 4. For models H000 and H010, the rare-earth and third peaks are underproduced by only a third to one half of an order of magnitude. But the first peak is underproduced by slightly more than an order of magnitude compared to the second peak in those models. As we go to longer HMNS lifetimes, the rare-earth and third peaks are underproduced by about an order of magnitude regardless of the HMNS lifetime. At the same time, the first peak increases from an under-production of 2/3 of an order of magnitude at τ = 30 ms to an over-production of a factor of 2.6 at τ = ∞.

The spinning BH models underproduce the first, rare-earth, and third peaks in roughly the same amounts, namely between about half to 3/4 of an order of magnitude compared to the second peak. The compact disc model BF15 underproduces the third peak by a similar amount as the HMNS models with long lifetimes. But it also has the lowest rare-earth peak abundance relative to solar, with an underproduction factor of about 1.5 orders of magnitude. The first peak, on the other hand, is only underproduced by a factor of ∼2.

3.4.2 BH spin mimicking HMNS lifetime

Metzger & Fernández (2014) proposed using the relative amount of blue optical emission in the kilonova as an observational test of the HMNS lifetime, given that the latter correlates with the amount of high-Y_e material ejected in the disc outflow. Given that high-Y_e material also correlates with BH spin (Fernández et al. 2015a), the kilonova signature of the two types of central objects can overlap. Here we investigate whether the nucleosynthesis signature offers additional information that can break the degeneracy. We consider the BH models H000, B070 and B090, which have the same mass and varying spins χ = 0, 0.7 and 0.9, respectively, as well as model BF15, which has a larger mass and χ = 0.86. Fig. 5 shows the electron fraction distributions of these BH models.

Figure 5.

Mass histograms of electron fraction for the ejecta from BH models. We include the HMNS model H100 for comparison.

While a larger BH spin has a similar overall effect on the disc ejecta composition as a longer HMNS lifetime, even a spin χ = 0.9 does not reach the same amount of ejecta mass and average value of Y_e,5GK as an HMNS with a lifetime τ = 100 ms (cf. Fig. 5). At best, a rapidly spinning BH can mimic an HMNS of modest lifetime. This can also be seen in Fig. 4 and in the peak ratio values in Table 2. Models B070 and B090 have values of [Y_1st/Y_2nd] that are similar or slightly smaller than in model H030, while the values of [Y_RE/Y_2nd] and [Y_3rd/Y_2nd] from models B070 and B090 fall between those from models H010 and H030. Thus, a BH with a spin χ = 0.7–0.9 produces a disc outflow with a similar final abundance pattern as a HMNS with a lifetime of τ ∼ 15–20 ms. Therefore, the fact that we assume a non-spinning BH after collapse of the HMNS only has a very small impact on the nucleosynthesis. Model BF15, on the other hand, looks more like a HMNS with a lifetime τ ∼ 60 ms, judging from the underproduction of the first and third peaks relative to the second peak. None the less, the rare-earth peak is underproduced by almost half an order of magnitude more in the disc ejecta by BF15 than by the HMNS models with long lifetimes.

3.4.3 Lanthanides, actinides and heating rates

The energy released by the radioactive decay of heavy elements synthesized by the r-process can power an optical or infrared transient called a kilonova (also called macronova in the literature; Li & Paczyński 1998; Kulkarni 2005; Metzger et al. 2010; Roberts et al. 2011; Kasen, Badnell & Barnes 2013; Tanaka & Hotokezaka 2013; Grossman et al. 2014; Metzger 2017; Wollaeger et al. 2017). The two most important microphysical components that determine the light curve and spectrum of kilonovae are the opacity of the material and the radioactive heating rate. Lanthanides (58 ≤ Z ≤ 71) and actinides (90 ≤ Z ≤ 103) have open f-shells, which gives them a very complex atomic line structure that leads to broad-band opacities that are more than an order of magnitude larger than opacities from iron-group elements (Kasen et al. 2013; Tanaka & Hotokezaka 2013; Fontes et al. 2015).

The trajectory-averaged lanthanide 〈X_La〉 and actinide 〈X_Ac〉 mass fractions of the ejecta for all models are summarized in Table 2. For HMNS lifetimes τ ≤ 10 ms, the lanthanide mass fraction is a few times 10⁻², which will result in opacities about an order of magnitude larger than iron group opacities (see fig. 10 in Kasen et al. 2013). Note that the actinide mass fractions are on the order of 10⁻³, which is still a significant contribution to the opacity. Once the HMNS lifetime is longer than about 10 ms, the lanthanide fraction 〈X_La〉 ∼ 10⁻³ monotonically decreases as the lifetime increases. In these models, the lanthanides and actinides will increase the opacity only by a factor of a few relative to iron group opacities. In the BH models, increasing the spin from 0.7 to 0.9 reduces the lanthanide and actinide mass fractions by roughly a factor of 2 from around 10⁻². The compact disc model BF15 has very similar lanthanide and actinide mass fractions as H300. These results apply to the disc outflow component alone; the colour of a kilonova depends also on the spatial distribution and composition of the dynamical ejecta, which tends to be more neutron-rich and hence may contain significant amounts of lanthanides and actinides.

Fig. 6 shows the combined lanthanide and actinide mass fraction of each ejected trajectory as a function of Y_e,5GK for models H010, Hinf, B090 and BF15. For most trajectories, the lanthanide and actinide fraction plummets as the initial electron fraction increases from 0.2 to 0.25 (see also Kasen, Fernández & Metzger 2015; Lippuner & Roberts 2015). Particles with Y_e,5GK ≤ 0.25 have low entropies, because a higher entropy requires either significant neutrino heating, which increases Y_e, or viscous heating over time-scales comparable to the thermal time, which gives weak interactions enough time to increase Y_e towards its equilibrium value. This correlation between Y_e,5GK and s_5GK at Y_e,5GK ≲ 0.3 can be seen in Fig. 3.

Figure 6.

Top: Scatter plots of final lanthanide and actinide mass fraction in each trajectory as a function of the electron fraction Y_e,5GK, for selected models. Points are colour-coded by the initial entropy s_5GK. Bottom: Radioactive heating rate ε at one day as a function of Y_e,5GK, also colour-coded by s_5GK.

Fig. 6 also shows that there is a separate population of particles that contain significant amounts of lanthanides and actinides for electron fractions in the range 0.25 ≲ Y_e,5GK ≲ 0.3, particularly in models B090, BF15, and to a lesser extent in model H010. This group of trajectories has higher entropies (s_5GK ≳ 30 k_B baryon⁻¹) relative to those which do not make lanthanides or actinides for Y_e,5GK ≥ 0.25. We can understand this population by considering the strong dependence of the lanthanide/actinide abundance on the neutron-to-seed ratio Y_n/Y_seed, where Y_seed = ∑_{A ≥ 12}Y_A. In general, a neutron-to-seed ratio Y_n/Y_seed ∼ 40 is necessary to produce a significant amount of lanthanides and actinides: as Y_n/Y_seed increases from about 35 to about 45, X_La + X_Ac increases from ∼10⁻⁵ to ∼10⁻². This effect is illustrated in Fig. 7 for model B090, where the different particle populations are shown in different colours. The trajectories in the range Y_e,5GK = 0.25–0.30 lie on both sides of the Y_n/Y_seed = 35 boundary for the range of entropies encountered in the disc ejecta, with a critical entropy for lanthanide production s_5GK ∼ 30 − 40 k_B baryon⁻¹. Outside of this range of Y_e, lanthanide production is insensitive to s_5GK. The neutron-to-seed ratio increases with increasing entropy because a higher entropy prefers a larger number of particles, and thus the composition contains more lighter particles such as free neutrons.

Figure 7.

Neutron-to-seed ratio at the time when the temperature T ≥ 5 GK for the last time, versus the specific entropy at that time, for ejected trajectories from model B090. Nuclides with A ≥ 12 are counted as seeds. The grey dashed line indicates the minimum value of Y_n/Y_seed = 35 that is required to make lanthanides and actinides. All trajectories with Y_e,5GK < 0.25 (black dots) are above the minimum neutron-to-seed ratio regardless of initial entropy, while almost all trajectories with Y_e,5GK > 0.30 (blue dots) are below it. The remaining trajectories (0.25 ≤ Y_e,5GK ≤ 0.30, red dots) can produce lanthanides or actinides if the initial entropy is ≳30–40 k_B baryon⁻¹.

The population of trajectories with Y_e,5GK = 0.25–0.30 and s_5GK ≳ 30 k_B baryon⁻¹ is absent in model Hinf (cf. Fig. 6). In order to achieve such entropies at modest electron fraction, neutrino irradiation cannot be too strong. In model Hinf, the weak interaction time-scale is comparable to the viscous heating time-scale, raising Y_e above 0.3. Since the HMNS is the strongest source of neutrinos while it persists, this population of trajectories can only exist if there is either no HMNS or if the HMNS collapses to a BH quickly. The eight trajectories in model Hinf that produce lanthanides and actinides at Y_e,5GK ≳ 0.4 are an extreme case. They all have s_5GK > 200 k_B baryon⁻¹, which allows the neutron-to-seed ratio to be high enough to make lanthanides and actinides even at Y_e > 0.4. These particles attain this high entropy because of late-time fallback into (and then rapid ejection from) the innermost part of the disc, where they experience significant heating past 5 GK.

Table 2 gives the total radioactive heating rate in the ejecta ε_tot at 1 d and at 7 d for all models. The quantity scales with the total ejecta mass. When considering the contribution of the disc outflow to a radioactively-powered kilonova, we need to keep in mind the contribution of the dynamical ejecta, which is not simply additive as in the case of nucleosynthesis. The larger fraction of lanthanides and actinides expected for the dynamical ejecta means that its optical opacity is likely to be larger than that of the disc outflow. The geometry of the dynamical ejecta depends on the type of merger involved: quasi-spherical for NSNS mergers (e.g. Hotokezaka et al. 2013), or confined to the equatorial plane for NSBH mergers (e.g. Kawaguchi et al. 2015; Foucart et al. 2017). In the former case, the disc ejecta can be obstructed in all directions by high-opacity material, whereas in the latter a long-lasting blue optical component from the disc ejecta can be detectable from some directions (Kasen et al. 2015). A short-lived blue optical component should be detectable in most cases since the lanthanide-rich material goes through a high-temperature phase and its outer layers let photons escape more freely (Barnes & Kasen 2013; Fernández et al. 2017).

If we only consider the heating rates and lanthanide/actinide content from the disc outflow, we expect models H000 and H010 to make dim infrared kilonovae, while all other models should make bright blue optical kilonovae. The exception is B070, which has about 1.5 per cent lanthanides and actinides by mass and a fairly large heating rate, so this model could make a brighter infrared kilonova. The heating rates reported here are upper limits to the bolometric luminosities of kilonovae from the disc outflow, since the conversion of energy from radioactive decay products into thermal energy of the ejecta has a limited efficiency (e.g. Metzger et al. 2010; Hotokezaka et al. 2016; Barnes et al. 2016).

The bottom row of Fig. 6 shows the heating rates of the individual trajectories at t = 1  d as a function of Y_e,5GK. The heating rate decreases slowly with increasing electron fraction for Y_e,5GK ≲ 0.3. For higher initial Y_e, the decrease steepens until large oscillations start around Y_e,5GK ∼ 0.4. This general behaviour is consistent with the findings of the parametrized r-process nucleosynthesis study of Lippuner & Roberts (2015). The oscillations are due to the initial NSE composition, which can be dominated by a small number of individual nuclei that match the electron fraction of the material. If there is a nuclide with a matching Y_e that has a half-life of about a day, there will be strong heating at t ∼ 1  d as this nuclide decays, but if there is no such nuclide, then there will be significantly less heating. The peak in the heating rate at Y_e,5GK ∼ 0.425, visible in all but the first column of Fig. 6, is caused by ⁶⁶Ni, which has Y_e = 28/66 ∼ 0.424 and is a dominant nuclide in the initial NSE composition. ⁶⁶Ni decays to ⁶⁶Cu with a half-life of 55 h, which then decays to stable ⁶⁶Zn with a half-life of 5 min. Around Y_e,5GK ∼ 0.45, there is a large dip in the heating rate because the initial composition is dominated by ⁶²Ni and ⁶⁶Zn, both of which are stable and have Y_e = 0.45. The little neutron capture that takes place mainly builds up ⁸⁸Sr, which is also stable and has Y_e = 0.43. At Y_e,5GK ∼ 0.49 there is another peak in the heating rate due to ⁶²Zn, ⁶¹Cu and ⁵⁷Ni, which are all unstable and very abundant since they have electron fractions between 0.48 and 0.49. At higher entropies, the initial neutron to seed ratio is enhanced, which will generally produce unstable nuclei and thus possibly break the oscillatory pattern of heating rate versus Y_e,5GK.

The morphology of the ejecta for all models with non-zero viscosity is shown in Fig. 8, with the top row showing final lanthanide and actinide mass fractions. Since the material is close to homologous expansion, Fig. 8 essentially shows the velocity space distribution of the ejected particles. In all models, most of the particles are concentrated in a central blob. When the HMNS persists for 100 ms or more, some particles attain higher velocities (v ∼ 0.1 c) and organize themselves into multiple shells (which arise from shock trains in the gas; e.g. Matsuo, Miyazato & Kim 1999). Since particles are accelerated to high velocities predominantly by neutrino irradiation, there are more high-velocity particles the longer the HMNS persists.

Figure 8.

Spatial distribution of ejected trajectories at time t = 1  d, for all HMNS and BH models with non-zero viscosity. Positions are computed by assuming that the velocity of each particle at the end of the simulation (∼10 s) remains constant at later times. Aside from smoothing of sharp features by radioactive heating at late times (Rosswog et al. 2014; Fernández et al. 2015b), this extrapolation of trajectories yields a good first approximation to the morphology of the homologously expanding disc wind ejecta. Top: Final lanthanide and actinide mass fraction. Because the high-lanthanide component of the ejecta is very small in most models, the points with X_La + X_Ac ≥ 10⁻³ are plotted above all other points with lower mass fractions (except for models H000 and H010). The lanthanide and actinide mass fractions stop changing a few seconds after neutrons are exhausted. Therefore, the final mass fractions are the same as the ones prevailing during the kilonova emission phase. Bottom: Radioactive heating rate at one day. In this row, the points for the different trajectories are drawn in a random order for all models.

In some models (H000, H010 and B070) there appears to be little structure in the distribution of lanthanides and actinides. In others there is a clear preference for the high-lanthanide trajectories to be at low velocities and clustered in the equatorial plane. Model BF15 is an exception in that the few particles that make a significant amount of lanthanides and actinides achieve the highest velocities, because they experienced a significant amount of heating (cf. Fig. 6). The few high-entropy trajectories in model Hinf that make lanthanides and actinides at Y_e,5GK ≳ 0.4 are the high-velocity points moving predominantly in the polar direction, i.e. they have a large absolute y coordinate and a small x coordinate in Fig. 8. For the heating rates, on the other hand, there is no strong spatial correlation between heating magnitude and trajectory location in any of the models. Models H300 and Hinf are the only ones which have a significant component of low heating rate trajectories, although these particles follow the same velocity distribution as those with high heating rates. Note that the low-heating particles are easier to see in the outer part of the ejecta because the lower particle density, but they are also present in the central blob.

3.4.4 Comparison with r-process abundances in metal-poor halo stars

Metal-poor Galactic halo stars that show r-process elements in their spectra are thought to have been enriched by a single or few r-process nucleosynthesis events (Cowan et al. 1999). Fig. 9 shows the final elemental abundances (at t ∼ 10 Myr) of our investigated models, along with the solar r-process abundances from Arnould et al. (2007) and from four metal-poor halo stars as reported in Westin et al. (2000) and Roederer et al. (2012). Abundances are scaled to give the best possible match to the solar values in the range 56 ≤ Z ≤ 75. There is little variation between the different models in this range, and a generally good match to solar and metal-poor star abundances. While Eu (Z = 63) is underproduced by a factor of a few with respect to solar in our models, the metal-poor halo stars also exhibit slightly sub-solar Eu abundances.

Figure 9.

Top: Final elemental abundances for all models with non-zero viscosity. In each case, abundances are scaled so that ∑(log Y_Z/Y_{Z, ⊙})² is minimized for 55 ≤ Z ≤ 75. Black dots show Solar System r-process abundances (Arnould et al. 2007). Bottom: Ratio of model abundances to Solar System r-process abundances (diamonds). Also shown are the observed abundances of four metal-poor halo stars (shown as stars, from Westin et al. 2000; Roederer et al. 2012), scaled to the Solar System values in the same way as model abundances.

Regarding elements in the range 44 ≤ Z ≤ 55, Fig. 9 shows that they are overproduced relative to the solar values by all models except H000. In the range Z = 30–42, our disc outflow models do not match the solar abundances, but none the less agree with the overall increase in abundances with Z (compared to solar) exhibited by metal-poor halo stars. From Z = 44 to 47, the metal-poor halo stars have a declining abundance trend compared to solar, which none of our models reproduce. In fact, the disc outflow generates an increasing trend from Z = 44 to 47. Thus, we can conclude that both the solar and metal-poor halo stars’ abundance patterns are inconsistent with pure disc outflow nucleosynthesis. This means that the nucleosynthesis event (or events) that enriched these metal-poor halo stars in heavy r-process elements could not have been the r-process in merger disc outflow alone. There must have been contributions from additional types of ejecta, or perhaps a different kind of nucleosynthesis event altogether. Contributions from very neutron-rich dynamical ejecta that produces mainly second to third peak material (Z ≳ 55) could make the combined final abundance pattern consistent with the metal-poor halo star observations (cf. Just et al. 2015). For Z ≲ 42, we expect supernovae to contribute to the nucleosynthesis, thus making it difficult to draw any conclusions from the disc outflow nucleosynthesis alone.

3.5 Impact of angular momentum transport

Neutrino irradiation from an HMNS is strong enough that significant amounts of material can in principle be ejected by neutrino energy deposition alone (Ruffert et al. 1997). Therefore, it is useful to clarify the contribution of angular momentum transport processes to the overall mass ejection from the disc. Transport of angular momentum modifies the evolution of the disc in two ways: (1) it causes part of the disc to accrete and part to spread outward as the local contribution of centrifugal forces to hydrostatic balance is modified, and (2) it dissipates energy in the form of heat, increasing the local entropy of the gas. While the detailed form of these two effects depends on the way the process is modelled (i.e. α viscosity), the overall contribution to the disc evolution is generic (e.g. magnetohydrodynamic turbulence also heats up the gas, but with a different spatial distribution than shear viscosity; e.g. Hirose, Krolik & Stone 2006).

We focus here on two models of an HMNS which does not collapse: one including angular momentum transport through α viscosity like all other models (Hinf), and another in which the viscosity is set to zero (HinfNoVisc). Fig. 10 shows the distribution of the ejected particles as a function of initial electron fraction Y_e,5GK, initial specific entropy s_5GK and nucleosynthesis start time (t_5KG).

Figure 10.

Properties of the disc wind ejecta from an HMNS that does not collapse, including and excluding the effect of viscosity on angular momentum transport and heating (models Hinf and HinfNoVisc, respectively). Top: Histograms of the ejected mass in units of 10⁻³ M_⊙ as a function of Y_e,5GK (left) and as a function of the time t_5GK when each particle has cooled to 5 GK (right). Bottom: Scatter plot of Y_e,5GK versus s_5GK, colour-coded by t_5GK for model Hinf (left) and HinfNoVisc (right).

A non-zero viscosity increases the total ejecta mass from 5.5 × 10⁻³ M_⊙ (HinfNoVisc) to 29.6 × 10⁻³ M_⊙ (Hinf), which corresponds to 18 per cent and 99 per cent of the initial disc mass, respectively. Not all of this additional mass is directly ejected by angular momentum transport, however. The scatter plot of Y_e,5GK versus s_5GK in Fig. 10 shows a clear neutrino-driven wind pattern for part of the ejecta from model Hinf (see also Section 3.3 and Fig. 3). This pattern is associated primarily with the early ejecta (t_5GK ≲ 0.2 s), which covers the range Y_e,5GK = 0.2 − 0.54. Some particles that start nucleosynthesis later (t_5GK ∼ 0.5 s), with electron fractions in the range Y_e,5GK = 0.4–0.54, also exhibit strong correlations between Y_e,5GK and s_5GK. This later ejecta could thus also be neutrino-driven, but the different thermodynamic conditions prevailing at later times alter the exact relationship between Y_e,5GK and s_5GK. Ejecta that begins nucleosynthesis at later times (t_5GK ≳ 0.8 s), with Y_e,5GK = 0.3–0.4, shows a much larger dispersion in s_5GK for a given Y_e,5GK. This dispersion is associated with convective motions in the advective phase of the disc, in which mass ejection is driven exclusively by viscous heating and nuclear recombination (Section 3.3).

In contrast, Fig. 10 shows that the model with zero viscosity ejects almost exclusively material with Y_e,5GK = 0.47–0.52, with ejection happening predominantly at late times (t_5GK ≳ 1 s). All of the trajectories in model HinfNoVisc exhibit the characteristic neutrino-driven wind correlation between electron fraction and entropy. Given that the model with viscosity has more neutrino-driven wind ejecta than model HinfNoVisc, we infer that angular momentum transport not only adds an additional late-time ejecta component, but it also enhances the neutrino-driven wind itself. This enhancement arises from two effects. First, spreading of the outer regions of the disk due to angular momentum transport moves material into shallower regions of the gravitational potential, where thermal unbinding requires less energy injection. Second, viscous heating also acts on the neutrino-driven wind ejecta, increasing the rate of internal energy gain. This also explains why the neutrino-driven component is delayed in the case without viscosity: particles are more tightly bound, requiring energy deposition for a longer time in order to be ejected.

Fig. 11 shows the final abundances for models Hinf and HinfNoVisc. Without viscosity, most of the ejecta has Y_e,5GK ∼ 0.5, and so the final abundance is dominated by the iron peak (A ∼ 56) and ⁴He. Note that HinfNoVisc produces the same iron peak and ⁴He abundance as Hinf despite its total ejecta mass being more than a factor of 5 lower. For elements significantly heavier than the iron peak, the final abundance of model HinfNoVisc is one to two orders of magnitude lower than in model Hinf. While model HinfNoVisc can still make the third r-process peak, the material is produced by just a handful of particles that have Y_e,5GK ≤ 0.25.

Figure 11.

Final abundances of non-collapsing HMNS models with and without viscosity (Hinf and HinfNoVisc, respectively). Model abundances are scaled by the total ejecta mass of each model, while the Solar System r-process abundances (Arnould et al. 2007) are scaled to match the second peak of model Hinf.

The nucleosynthesis in the disc outflow from an HMNS has also been studied by Martin et al. (2015), who focused on the neutrino-driven wind without the contribution from angular momentum transport. In their study, they use the ejection time of a particle as a proxy for the HMNS lifetime. Martin et al. find that the total ejecta mass increases with increasing HMNS ejection time, and that the ejecta also becomes less neutron-rich at later times. These results are broadly consistent with what we find in model HinfNoVis.

A more detailed comparison between the results from Martin et al. and ours shows that the initial disc mass can have a significant impact on the properties of this ejecta component. Martin et al. find most of the ejecta having Y_e ∼ 0.3–0.4, and a majority of the particles ejected on a time-scale of 100 ms, whereas model HinfNoVisc has Y_e ∼ 0.5 and most trajectories begin nucleosynthesis after 1 s. The initial condition in Martin et al. is the output of a three-dimensional neutron star merger simulation by Perego et al. (2014), with a disc that is more than six times as massive as in our model (0.19 M_⊙ compared to 0.03 M_⊙). Different physical and thermodynamic conditions in the disc can lead to different results in the nucleosynthesis (e.g. compare the output of our models B070 with BF15). Another source for the differences in the ejecta properties can be the different hydrodynamic methods used, which can have different amounts of numerical viscosity. Finally, one similarity between the results of Martin et al. and ours is the characteristic values for the entropy in the neutrino-driven wind (∼20 k_B baryon⁻¹; compare the bottom right panel of Fig. 10 and fig. 5 in Martin et al. 2015).

4 CONCLUSIONS

We have performed nucleosynthesis calculations in the outflow from a neutron star merger accretion disc when the central object is an HMNS or a BH. We used long-term hydrodynamic simulations of accretion discs to model the ejecta, and detailed nucleosynthesis calculations were carried out on tracer particles using a nuclear reaction network. We have systematically varied the lifetime of the central HMNS to study its impact on the nucleosynthesis. Our simulations continued after the HMNS collapses to a BH, thus allowing us to investigate the long-term effects on the disc and nucleosynthesis. We have also performed some simulations that start with central BHs with different spins, to explore similarities and differences with the HMNS case.

Our results are consistent with previous findings regarding the monotonic increase in mass ejection and mean electron fraction of the disc outflow for longer HMNS lifetimes (Fig. 2; see also Metzger & Fernández 2014). This correlation results in the amount of ejecta that has initial Y_e ≤ 0.25 being almost constant once the HMNS lifetime is τ ≳ 30 ms. The final abundance pattern at large mass numbers (A ≳ 100), normalized by the total ejecta mass, is thus independent of the HMNS lifetime, because only material with Y_e ≲ 0.25 can make those heavy elements (Fig. 4). For very short HMNS lifetimes (τ ≲ 10 ms), there is more neutron-rich ejecta and thus the rare-earth and third r-process abundance peaks are about half an order of magnitude larger. For these short lifetimes, the disc outflow abundances alone are almost consistent with the solar r-process abundances.

For other cases, the inconsistency between the final abundances from the disc outflow and the solar values is not in itself problematic. In most neutron star mergers, we also expect a dynamical ejecta component that tends to be more neutron-rich. This means it could easily synthesize the heavy r-process material that the disc ejecta may be underproducing. If the dynamical ejecta is consistent with the solar r-process abundances between the second and third peak, then our results indicate a preference for short HMNS lifetimes, since abundances from these models are broadly consistent with solar r-process abundances between the second and third peak. If, on the other hand, the dynamical ejecta overproduces the third peak relative to the second peak, then we would require longer HMNS lifetimes, which result in underproduction of the third peak compared to the second, in order to make the combined abundance pattern consistent with solar. We draw similar conclusions from comparing the nucleosynthesis in our disc models to the observed abundances in metal-poor halo stars (Fig. 9).

If the HMNS lifetime is τ ≥ 30 ms, most of the ejected material has Y_e ≳ 0.25 and so only a small amount of lanthanides and actinides are synthesized. Thus we expect the disc outflow to produce a kilonova that peaks in the optical band on a time-scale of a day. However, this optical kilonova component may be obscured by lanthanides and actinides that were synthesized in the more neutron-rich dynamical ejecta, particularly because this ejecta component is quasi-spherical in NSNS mergers (e.g. Bauswein et al. 2013). For a BHNS merger on the other hand, viewing-angle effects are expected to be very important given that the dynamical ejecta should lie primarily in the equatorial plane, allowing emission from lanthanide-free material to more easily escape in the polar direction (e.g. Kasen et al. 2015). However, BHNS mergers do not make HMNSs, so visibility of the HMNS component relies on the dynamical ejecta from the NSNS binary having a small mass. If the HMNS lifetime is 10 ms or less, then the disc outflow contains substantial amounts of lanthanides and actinides. This would also produce an infrared kilonova on a time-scale of a week and could be indistinguishable from the dynamical ejecta contribution. Due to the possible masking of optical emission by the dynamical ejecta, the absence of optical emission could only be an indication for a short-lived HMNS if the viewing-angle can be tightly constrained to be face-on, e.g. by a coincident short gamma ray burst (sGRB) detection. The properties of the radioactive heating rates from our tracer particles (Fig. 6) are consistent with the findings of Lippuner & Roberts (2015). We find no strong spatial correlation between radioactive heating and particle location in the ejecta (Fig. 8).

By comparing the disc outflow abundances when a spinning BH or a HMNS is the central object, we find that a BH of spin 0.7–0.9 can mimic a HMNS with lifetime τ ∼ 15 ms. If the central BH is more massive (and the disc is compact, such as in a NSBH remnant), then it can mimic a longer-lived HMNS with τ ∼ 60 ms. For longer HMNS lifetimes (τ ≳ 100 ms), we find no possible overlap between the amounts of mass ejected or the nucleosynthesis signatures. Determining whether this difference translates into the kilonova signature requires more detailed radiation transport calculations using realistic initial conditions for the central object, disc and dynamical ejecta.

Regarding the disc outflow dynamics, we find that when a HMNS is present, two types of ejecta are clearly distinguishable: one driven primarily by neutrino energy deposition, showing a clear correlation between electron fraction and entropy, and another driven primarily by viscous heating and nuclear recombination, in which no such correlation exists (Fig. 3). We evolved a test model with no explicit angular momentum transport, and found that the contribution of viscosity is not limited to the late-time neutrino-independent outflow. Angular momentum transport also aids the ejection of the neutrino-driven component by moving material to shallower regions of the gravitational potential, and by accelerating thermal ejection through additional heating of the gas (Fig. 10). Thus, excluding angular momentum transport in the HMNS disc evolution can severely underestimate the amount of mass ejected.

While our hydrodynamics simulation includes the most relevant physics (neutrino transport, viscosity and general relativity), all of these are currently implemented in an approximate way. Neutrino interactions in the disc play a crucial role in setting the electron fraction and specific entropy distributions of the outflow. Likewise, the dominance of the HMNS irradiation in driving the early outflow requires a realistic treatment of the merger remnant instead of its approximation as a reflecting boundary. Hence our approximate treatment of the HMNS surface and neutrino interactions are the most important limitations of this study. We believe that our approximate neutrino treatment captures the most important aspects of the neutrino interactions in the disc and so we expect that our results exhibit the correct general trends. However, due to the importance of neutrinos in driving the disc outflow, it is clear that more sophisticated neutrino radiation transport methods are needed in future work. Incorporating an α viscosity prescription into our disc simulations is a significant step forward from simulations with no physical viscosity. Ultimately, however, accretion disc simulations with more realistic physical viscosity (e.g. provided by magnetoturbulence driven by the magnetorotational instability) need to be performed. Furthermore, due to the extreme computational complexity of fully general-relativistic hydrodynamics simulations, no self-consistent simulations of neutron star mergers with subsequent long-term accretion disc evolution have been carried out to date. Our approach of assuming an initial equilibrium torus is thus a necessary approximation in order to investigate nucleosynthesis in merger accretion discs. Mapping the early accretion disc structure found in merger simulations into another code for the disc evolution (as we have done for model BF15) is a more accurate approach, but it is still not fully self-consistent. With next generation general-relativistic hydrodynamics codes, we hope to be able to perform completely self-consistent simulations of neutron star mergers and accretion disc outflows in the future. Finally, we only use the FRDM nuclear mass model for the nucleosynthesis calculations in this study. In future work, we plan to investigate and compare other nuclear mass models for r-process nucleosynthesis in disc outflows.

Acknowledgements

JL and CDO acknowledge support from U.S. National Science Foundation (NSF) grant AST-1333520. RF acknowledges support from the University of California Office of the President, from NSF grant AST-1206097, and from the Faculty of Science at the University of Alberta. DK is supported in part by a U.S. Department of Energy (DOE) Office of Nuclear Physics Early Career Award, and by the Director, Office of Energy Research, Office of High Energy and Nuclear Physics, Divisions of Nuclear Physics, of the DOE under Contract No. DE-AC02-05CH11231. The software used in this work was in part developed by the Flash Center at the University of Chicago funded by the DEO Office of Advanced Scientific Computing Research (OASCR) and the Advanced Simulation and Computing (ASC) program of the DOE's National Nuclear Security Administration (NNSA). This research used resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Some computations were performed on the Edison compute cluster (repository m2058). Some of the calculations were performed on the Zwicky compute cluster at Caltech, supported by NSF under the Major Research Instrumentation (MRI) award PHY-0960291 and by the Sherman Fairchild Foundation. This work was supported in part by NSF grant PHY-1430152 (JINA Center for the Evolution of the Elements). CDO thanks the Yukawa Institute for Theoretical Physics for support and hospitality. This article has been assigned Yukawa Institute report number YITP-17-26.

APPENDIX: Improvements to the hydrodynamic disc simulations

We have corrected an error in the weak interaction rates used since Fernández & Metzger (2013). The error involved the absence of the neutron-proton mass difference in the argument of the Fermi–Dirac integrals for all the electron neutrino rates (but not in the antineutrino rates). After correction of this error, the electron neutrino emission rates decrease by a factor of up to 2 in regions with temperatures ≳1 MeV. The most significant change in the result is an increase in the fraction of matter ejected with Y_e > 0.25 for the non-spinning black hole (BH) case from a few per cent to ∼30 per cent of the outflow, resulting in an increase of ∼10 per cent in the average Y_e of the outflow. A similar increase in the mean Y_e of the wind is observed for the long-lived hypermassive neutron star (HMNS) case. In terms of the amount of mass ejected, the effects are strongest for a non-spinning BH, which ejects up to ∼10 per cent more mass than before the correction is applied (i.e. an additional ∼0.5 per cent of the initial disc mass). For a long-lived HMNS, the mass ejection increases by ∼0.2 per cent.

Relative to Metzger & Fernández (2014), we also account for separate neutrino and antineutrino temperatures for absorption. The mean energy for each neutrino species is calculated by taking the ratio of the globally-integrated energy to number emission rates, as in Ruffert, Janka & Schaefer (1996). The temperatures are then obtained through the relation kT_{ν, i} = F₄(0)/F₅(0)〈ε_{ν, i}〉, where F_i are the Fermi–Dirac functions of integer argument, F₅(0)/F₄(0) ≃ 5.065, and 〈ε_{ν, i}〉 is the neutrino mean energy. When using the updated weak interaction tables, this change translates into an increase of ∼10 per cent in mass ejection for a non-spinning BH, and a decrease of 0.4 per cent in mass ejection for a long-lived HMNS.

To put the effect of these modifications in perspective for the non-spinning BH case, doubling the resolution of the simulations in each dimension can change mass ejection by ∼10 per cent (Fernández et al. 2015b), and a change in the α viscosity parameter leads to an almost directly proportional change in mass ejection (i.e. a factor 10 when going from α = 0.01 to α = 0.1, Wu et al. 2016).