The landscape of disc outflows from black hole–neutron star mergers

ABSTRACT

We investigate mass ejection from accretion discs formed in mergers of black holes (BHs) and neutron stars (NSs). The third observing run of the LIGO/Virgo interferometers provided BH–NS candidate events that yielded no electromagnetic (EM) counterparts. The broad range of disc configurations expected from BH–NS mergers motivates a thorough exploration of parameter space to improve EM signal predictions. Here we conduct 27 high-resolution, axisymmetric, long-term hydrodynamic simulations of the viscous evolution of BH accretion discs that include neutrino emission/absorption effects and post-processing with a nuclear reaction network. In the absence of magnetic fields, these simulations provide a lower limit to the fraction of the initial disc mass ejected. We find a nearly linear inverse dependence of this fraction on disc compactness (BH mass over initial disc radius). The dependence is related to the fraction of the disc mass accreted before the ouflow is launched, which depends on the disc position relative to the innermost stable circular orbit. We also characterize a trend of decreasing ejected fraction and decreasing lanthanide/actinide content with increasing disc mass at fixed BH mass. This trend results from a longer time to reach weak freezout and an increasingly dominant role of neutrino absorption at higher disc masses. We estimate the radioactive luminosity from the disc outflow alone available to power kilonovae over the range of configurations studied, finding a spread of two orders of magnitude. For most of the BH–NS parameter space, the disc outflow contribution is well below the kilonova mass upper limits for GW190814.

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

1 INTRODUCTION

The Advanced LIGO interferometer has completed three observing runs – with Advanced Virgo joining part of the way – resulting in the official detection of 11 binary black hole (BH) mergers and two neutron star (NS) mergers (Abbott et al. 2019, 2020a; LIGO Scientific Collaboration & Virgo Collaboration 2020), with many more in candidate status at the time of this writing. The increased sensitivity of the third observing run also yielded an event that can be either a BH–BH or a BH–NS merger (GW190814) (e.g. Abbott et al. 2020b; Ackley et al. 2020; Andreoni et al. 2020; Coughlin et al. 2020; Thakur et al. 2020; Vieira et al. 2020). Only one of these events (GW170817), however, has had electromagnetic (EM) counterparts detected (Abbott et al. 2017). While multiple reasons can account for the lack of an EM detection (such as a large distance, large localization area, galactic extinction, or Sun constraints; e.g. Foley et al. 2020), the possibility remains that these sources were intrinsically fainter than the kilonova from GW170817.

BH–NS mergers can lead to a wide range of ejected masses depending on whether the NS is tidally disrupted by the BH. The outcome depends on the masses of the ingoing BH and NS, as well as on the spin of the BH and the compactness of the NS (e.g. Foucart, Hinderer & Nissanke 2018). The dynamical ejecta emerges as a very neutron-rich, equatorial tidal tail that quickly leaves the system. The nucleosynthesis properties and contribution to the kilonova transient are mostly set at the time of ejection (e.g. Roberts et al. 2017), and its properties can be parametrized by direct comparison with dynamical merger simulations (e.g. Kawaguchi et al. 2016; Krüger & Foucart 2020).

The accretion disc, on the other hand, ejects mass on a longer time-scale, as angular momentum is transported initially by gravitational torques and later by magnetohydrodynamic (MHD) turbulence (see e.g. Fernández & Metzger 2016; Shibata & Hotokezaka 2019 for reviews). This longer evolutionary time-scale allows weak interactions to modify the composition, resulting in a different r-process yield (and possibly a different kilonova colour) than the dynamical ejecta. The complexity of this evolution makes it very expensive to realistically model the disc, however. Only a handful of 3D general-relativistic (GR) MHD simulations of discs around BHs have been carried out including at least some important microphysics or neutrino effects (Siegel & Metzger 2017, 2018; Hossein Nouri et al. 2018; Christie et al. 2019; Fernández et al. 2019; Miller et al. 2019). Furthermore, all of these simulations either focus on a narrow subset of parameter space and/or do not evolve the system for long enough to achieve completion of mass ejection.

More extensive studies of BH accretion discs have been carried out using axisymmetric hydrodynamic simulations with a wide variety of approximations to the physics (Fernández & Metzger 2013; Fernández et al. 2015a; Just et al. 2015; Fujibayashi et al. 2020). None of these studies covers a significant fraction of all the possible BH accretion disc configurations, however.

Despite missing the magnetic field, hydrodynamic simulations can provide a good description of the late-time thermal component of the outflow that arises when weak interactions freeze out (Metzger, Piro & Quataert 2009) and heating of the disc by viscous stresses (in lieu of MHD turbulence) is unbalanced. Close comparison between GRMHD and hydrodynamic simulations show that the latter provide a lower limit to the fraction of the disc ejected, with magnetic enhancements dependent on the strength of the initial poloidal field in the disc (Christie et al. 2019; Fernández et al. 2019).

Here we carry out an extensive set of long-term hydrodynamic simulations of accretion discs around BH remnants, with the aim of sampling the entirety of the parameter space resulting from BH–NS mergers, and thereby improving parameter estimation models that take disc outflow properties as input (e.g. Barbieri et al. 2019; Hinderer et al. 2019; Coughlin et al. 2020). Along the way, a broad probe of parameter space allows us to identify trends in the disc ejection physics, helping to focus on areas where improvements in the physics (e.g. neutrino transport) have the most impact. Finally, by providing a lower limit to the disc mass ejected, we are also estimating the lower limit to the raw radioactive heating available to power kilonova transients.

The structure of the paper is as follows. Section 2 presents our methods, including our choice of initial conditions from the plausible parameter space of BH–NS mergers. Section 3 presents our results, divided into mass ejection, outflow composition, and implications for EM counterparts. We close with a summary and discussion in Section 4.

2 METHODS

2.1 Initial conditions

In order to sample representative initial disc masses for our simulations, we map the parameter space of BH–NS merger remnants using analytical formulae that are calibrated to numerical relativity simulations. For a given ingoing BH mass M_bh(in) and NS mass M_ns, we uniformly sample the range 9–13 km for NS radii and 0–0.7 for the ingoing BH spin, and compute distributions of (1) the remnant baryon mass left outside the BH using the formula of Foucart et al. (2018), (2) the disc mass using the formula of Krüger & Foucart (2020) for the dynamical ejecta, and (3) the post-merger BH mass M_bh(out) and its spin using the formulae of Pannarale (2014) and the output from steps (1)–(2). This approach is intended to be agnostic about the properties of the EOS of dense matter and initial BH spin distribution, within plausible limits.

The resulting cumulative distributions of initial disc masses are shown in Fig. 1 for two NS masses |$\lbrace 1.35,1.45\rbrace \, \mathrm{M}_\odot$| and ingoing BH masses such that the median value of M_bh(out) is |$\lbrace 3,5,8,10,15\rbrace \, \mathrm{M}_\odot$|⁠. The lowest post-merger BH mass is chosen to explore the hypothetical case of a very low-mass ingoing BH, or the prompt collapse of an NS–NS system. The highest BH mass is chosen such that at least |$10{{\ \rm per\ cent}}$| of mergers result in disruption. Median disc masses range from |$0.1 \, \mathrm{M}_\odot$| for the |$3 \, \mathrm{M}_\odot$| post-merger BH, to |$0.02 \, \mathrm{M}_\odot$| for the |$15\, \mathrm{M}_\odot$| BH. In most cases, sensitivity to the specific choice of NS mass does alter the shape of the histogram but not the extreme values. The post-merger BH spin distributions have medians in the range 0.85–0.9, except for the lowest mass BHs considered, which have spins >0.95.

Figure 1.

Cumulative distribution of post-merger disc masses, obtained by using analytic formulae for the post-merger BH mass and spin (Pannarale 2014), remnant mass outside the BH (Foucart et al. 2018), and mass in dynamical ejecta (Krüger & Foucart 2020). For fixed ingoing NS mass M_ns and BH mass M_bh(in), the intervals 9–13 km for NS radii and 0–0.7 for ingoing BH spins are uniformly sampled. The uncertainty range in the post-merger BH mass M_bh(out) indicates the range in median values obtained by using the NS masses shown, with the central value corresponding to |$M_{\rm ns}=1.4\, \mathrm{M}_\odot$|⁠. The triangles indicate the initial disc masses selected for our hydrodynamic simulations (cf. Table 1). The fraction of mergers that result in NS disruption for each {M_bh(in), M_ns} pair is shown in purple.

For each post-merger BH mass, Fig. 1 shows the disc masses sampled in our study. The lowest disc mass in all cases is taken to be |$0.01\, \mathrm{M}_\odot$|⁠, which is optically thin to neutrinos, while the largest disc mass is such that it is at the uppermost end of plausible values. In all of our models, we take the spin of the post-merger BH to be 0.8, which while somewhat lower than the median values of the distributions obtained, lies within the range of plausible values (except for |$M_{\rm bh(out)}=3\, \mathrm{M}_\odot$|⁠) and does not demand a prohibitively small time-step. Furthermore, the effect of BH spin on the disc outflow properties is known, with more mass with higher average electron fraction being ejected for higher BH spins (Fernández et al. 2015a; Fujibayashi et al. 2020).

Our simulations start from idealized equilibrium tori (Section 2.2) which approach a Keplerian angular velocity distribution after a few orbits. This equilibrium configuration requires more parameters in addition to the disc mass: a radius of maximum density, an entropy (internal energy content), and an electron fraction (composition). While these parameter choices would not be necessary if we mapped the disc directly from a merger simulation, our broad coverage of possible BH–NS combinations would be limited by accessible merger simulation data. The remaining initial disc parameters are therefore chosen by inspecting the output of numerical relativity simulations of BH–NS mergers and making educated guesses about these values in regimes not covered by simulations.

We adopt initial disc radii that roughly follow the location of density maxima outside of the BH in published BH–NS simulations (Etienne et al. 2009; Kyutoku, Shibata & Taniguchi 2010; Foucart et al. 2011, 2013, 2014; Kawaguchi et al. 2015; Kyutoku et al. 2015; Foucart et al. 2017; Brege et al. 2018). This radius of maximum density is not well defined, however, since it depends on (1) the time at which it is measured, (2) the metric, and (3) on whether the local density, surface density, or enclosed mass are reported. Since there is no consistency across the literature for this quantity, we adopt a fiducial set of initial disc radii R_d = {50, 50, 60, 90, 120} km for post-merger BH masses |$\lbrace 3,5,8,10,15\rbrace \, \mathrm{M}_\odot$|⁠, respectively.

The entropy of all discs is taken to be |$8\, k_{\rm B}$| per baryon, which results in ratios of isothermal sound speed to orbital speed |${\sim}10{-}30{{\ \rm per\ cent}}$| at the point of maximum density. While this choice has some effect on the amount of mass ejected, we use a constant value for uniformity. Finally, the default electron fraction of the discs is set to |$Y_{e,\rm ini} = 0.2$|⁠, although we vary this parameter in our simulations given that it is dependent on the quality of the neutrino transport implementation in the merger simulation, and has a non-negligible effect on the disc outflow composition.

2.2 Hydrodynamic simulations

We perform time-dependent hydrodynamic simulations with flash version 3.2 (Fryxell et al. 2000; Dubey et al. 2009), with the modifications described in Fernández & Metzger (2013), Metzger & Fernández (2014), Fernández et al. (2015a), and Lippuner et al. (2017). The code solves the Euler equations in axisymmetric spherical polar coordinates (r, θ), subject to source terms that include the pseudo-Newtonian gravitational potential of a spinning BH (Artemova, Bjoernsson & Novikov 1996) without disc self-gravity, shear viscosity with an α parametrization (Shakura & Sunyaev 1973), and a leakage scheme for neutrino emission, with absorption included as a disc like light bulb (Fernández & Metzger 2013; Metzger & Fernández 2014). We only include electron type neutrinos/antineutrinos interacting with nucleons via charged-current weak interactions. The code employs the equation of state (EOS) of Timmes & Swesty (2000) with the abundances of neutrons, protons, and alpha particles in nuclear statistical equilibrium (NSE) above a temperature T = 5 × 10⁹ K and accounting for the nuclear binding energy of these particles. The electron–positron quantities are extended above the high-density limit of the table using analytic expressions (Bludman & van Riper 1978; Bethe, Applegate & Brown 1980).

The initial condition is an equilibrium torus with constant angular momentum, entropy, and electron fraction, with mass fractions assumed to be in NSE (e.g. Fernández & Metzger 2013). Parameters are chosen according to Section 2.1. The floor of density is set to 10 g cm⁻³ at r = 4R_d, and has an initial radial dependence r⁻². For r ≤ 4R_d, the radial exponent of the floor is smoothly brought to zero on a time-scale of 40 orbital times at r = R_d, reaching a flat floor in this region (Fernández et al. 2019, see also Just et al. 2015). The initial ambient density is set at 1.1 times the floor.

The computational domain extends from an inner radius r_in midway between the radius of the innermost stable circular orbit (ISCO) r_isco and the BH horizon, to an outer radius r_out = 10⁴r_in, with the polar angle spanning the range [0, π]. The grid is discretized logarithmically in radius, using 128 cells per decade, and a polar grid equi-spaced in cos θ using 112 cells. On the equatorial plane, this results in a spacing |$\Delta r /r \simeq 1.8{{\ \rm per\ cent}} \simeq 1^\circ \simeq \Delta \theta$|⁠. This resolution is double that of the models in Fernández et al. (2017), equivalent to that of the high-resolution models of Fernández & Metzger (2013) and Fernández et al. (2015a), and the same as in Fahlman & Fernández (2018) and the hydrodynamic models of Fernández et al. (2019). The boundary conditions are set to outflow in radius and reflecting in θ.

2.3 Nuclear reaction network post-processing

Passive tracer particles are initially placed in the disc following the density distribution. For each hydrodynamic simulation we employ 10⁴ particles, each representing an equal amount of mass. Particles are advected with the flow and record various kinematic and thermodynamic quantities as a function of time. Particles that are ejected with positive Bernoulli parameter beyond a radius of 10⁹ cm by the end of the simulation are considered to be part of the disc outflow.

Outflow trajectories are post-processed with the nuclear reaction network SkyNet (Lippuner & Roberts 2017), using the same settings as in Lippuner et al. (2017). The network employs 7843 nuclides and more than 1.4 × 10⁵ reactions, including strong forward reaction rates from the REACLIB data base (Cyburt et al. 2010) with inverse rates computed from detailed balance; spontaneous and neutron-induced fission rates from Frankel & Metropolis (1947), Mamdouh et al. (2001), Wahl (2002), and Panov et al. (2010); weak rates from Fuller, Fowler & Newman (1982), Oda et al. (1994), Langanke & Martínez-Pinedo (2000), and the REACLIB data base; and nuclear masses from the REACLIB data base, which includes experimental values were available, or otherwise theoretical masses from the finite-range droplet macroscopic model (FRDM) of Möller et al. (2016).

The rates of electron neutrino/antineutrino absorption/emission recorded by the trajectory are included in the evolution of the proton and neutron fraction. Likewise, the temperature and entropy are evolved self-consistently by accounting for nuclear heating from the network, as well as viscous heating and neutrino heating/cooling in the hydrodynamic simulation as recorded by the trajectory.

Processing begins when the trajectory reaches 10 GK for the last time, or when the temperature is maximal if lower than 10 GK at all times. For the portion of the evolution in which the temperature is higher than 7 GK, abundances are evolved in NSE, subject to neutrino interactions, while full network integration is carried out at lower temperatures. Trajectories are extended beyond the end of the simulation (12–25 s) by assuming that the density decays with time as t⁻³, to allow r-process nuclei with long half-lives to decay. Since the r-process is complete by the time this transition is made, most of the nuclear heating has already been deposited and the exact time dependence of the density decay is not important. While trajectories are evolved until 30 yr, information is extracted at t = 1 d and t = 1 week to estimate the properties of the kilonova at peak.

2.4 Models evolved

Table 1 shows all of the hydrodynamic models we evolve. As a baseline set, we take discs with initial conditions as described in Section 2.1: black hole and disc masses as in Fig. 1, initial entropy 8k_B per baryon, and initial electron fraction Y_e,ini = 0.2. Model names follow the convention bXXdYY, where XX and YY refer to the BH mass and disc mass, respectively. Tori are constructed as an equilibrium solution to the momentum equation with constant Bernoulli parameter, constant angular momentum, and the pseudo-Newtonian potential of the BH (Fernández & Metzger 2013). The torus shape is controlled by a dimensionless distortion parameter d (Stone, Pringle & Begelman 1999) which is solved for by fixing the entropy, Y_e, and torus mass. The distortion parameter is related to the torus initial Bernoulli parameter b_ini, black hole mass M_bh ≡ M_bh(out), and radius of initial density peak R_d by

$$\begin{eqnarray*} b_{\rm ini} = -\frac{1}{2d}\frac{GM_{\rm bh}}{R_d}. \end{eqnarray*}$$

(1)

The baseline set is evolved with a viscosity parameter α = 0.03

Table 1.

Hydrodynamic models evolved and input parameters. Columns from left to right show model name, black hole mass, disc mass, radius of initial disc density peak, initial electron fraction, torus distortion parameter, viscosity parameter, and maximum evolution time.

Model	Mbh	Md	Rd	|$Y_{e,\rm ini}$|	d	α	tmax
	(M⊙)	(M⊙)	(km)				(s)
b03d01	3	0.01	50	0.20	1.52	0.03	16.5
b03d03		0.03			1.82
b03d10		0.10			2.40
b03d30		0.30			3.37
b05d01	5	0.01			1.39		12.2
b05d03		0.03			1.59
b05d10		0.10			1.97
b05d30		0.30			2.57
b08d01	8	0.01	60		1.30		12.1
b08d03		0.03			1.43
b08d10		0.10			1.67
b08d20		0.20			1.87
b10d01	10	0.01	90		1.20		20.6
b10d03		0.03			1.30
b10d10		0.10			1.46
b10d20		0.20			1.59
b15d01	15	0.01	120		1.15		25.4
b15d03		0.03			1.22
b15d10		0.10			1.32
b03d01-y10	3	0.01	50	0.10	1.50		16.5
b03d30-y10		0.30			3.30
b08d03-y10	8	0.03	60		1.42		12.1
b03d01-y15	3	0.01	50	0.15	1.51		16.5
b03d30-y15		0.30			3.32
b08d03-y15	8	0.03	60		1.43		12.1
b03d01-v10	3	0.01	50	0.20	1.52	0.10	16.5
b08d03-v10	8	0.03	60		1.43		12.1

To assess the effects of initial composition, we evolve a few models with lower initial electron fraction than the baseline set Y_e,ini = {0.10, 0.15}. Likewise, we evolve two models with higher viscosity parameter, α = 0.1.

All hydrodynamic models are evolved for 5000 orbits at the initial density peak radius, which corresponds to ≃12–25 s of physical time (Table 1). This time is chosen such that the mass ejection from the disc is mostly complete. Tracer particles from each simulation are then post-processed with the nuclear reaction network as described in Section 2.3

3 RESULTS

The overall evolution of neutrino cooled accretion disks follows well-known stages (e.g. Popham, Woosley & Fryer 1999; Ruffert & Janka 1999; Di Matteo, Perna & Narayan 2002; Setiawan, Ruffert & Janka 2006; Chen & Beloborodov 2007; Lee, Ramirez-Ruiz & López-Cámara 2009). Depending on the initial disc mass, neutrinos can be trapped or escape freely. In the former case, an initial optically thick phase ensues until the density has decreased sufficiently for transparency. Thereafter, neutrino cooling is important compared to viscous heating, the inner disc is not too thick vertically, and accretion proceeds efficiently. After about a viscous time |$R_{\rm d}^2/(\alpha c_i^2/\Omega _{\rm K})\sim$| few 100 ms (with c_i the isothermal sound speed and Ω_K the orbital frequency; Shakura & Sunyaev 1973), the density becomes low enough that weak interactions freeze out, shutting down cooling (e.g. Metzger et al. 2009). At this point, the disc is radiatively inefficient, with viscous heating and nuclear recombination of α particles being unbalanced by cooling, thus an outflow is launched until no more mass is available to be ejected.

In the absence of magnetic fields, this thermal outflow is the only relevant mass ejection channel when a BH sits at the centre, as neutrino-driven winds are weak given that self-irradiation is not efficient (e.g. Just et al. 2015). A comparison with long-term GRMHD simulations shows that the outflow from hydrodynamic simulations is of similar quantity and has similar velocities as the analogue process occurring due to dissipation of MHD turbulence (Fernández et al. 2019). The magnetic field provides for additional, faster components that can eject a comparable amount of mass than the thermal outflow.

In the following, we discuss mass ejection properties across the range of models we evolve, the composition of these ouflows, and the implications for EM counterparts of BH–NS mergers.

3.1 Mass ejection

The disc mass ejection rate in all directions |$\dot{M}_{\rm out}$| is measured at a radius r_out = 10⁹ cm, far enough outside the disc that mostly complete ejection is achieved before the disc viscously spreads to that radius. Material is considered to be unbound when its Bernoulli parameter

$$\begin{eqnarray*} b = \frac{1}{2}\left[v_r^2 + v_\theta ^2 + \left(\frac{j}{r\sin \theta }\right)^2 \right] + e_{\rm int} + \frac{p}{\rho } + \Phi \end{eqnarray*}$$

(2)

is positive. In equation (2), v_r and v_θ are the radial and meridional velocities, e_int is the specific internal energy, p is the pressure, ρ is the density, j is the specific angular momentum, and Φ is the gravitational potential. Table 2 shows the ejected outflow mass M_ej – the time integral of |$\dot{M}_{\rm out}$| over the simulation – for all models.

Table 2.

Summary of results. Columns from left to right show model name, disc compactness (equation 3), ejected mass, fraction of initial disc mass ejected M_ej/M_d, average ouflow velocity, average outflow electron fraction, ejecta mass with X_{La + Ac} < 10⁻⁴ (M_blue), ejecta mass with X_{La + Ac} > {10⁻⁴, 10⁻³, 10⁻²} ({M₋₄, M₋₃, M₋₂}, respectively), and radioactive heating power (in units of 10⁴⁰ erg s⁻¹) at 1 and 7 d, ignoring thermalization efficiency.

Model	Cd	Mej	Mej/Md	〈v/c〉	〈Ye〉	Mblue/Mej	M−4/Mej	M−3/Mej	M−2/Mej	L40,1d	L40,1w
		(⁠|$10^{-2}\, \mathrm{M}_\odot$|⁠)		(10−2)
b03d01	0.60	0.21	0.21	3.4	0.28	0.75	0.25	0.18	0.10	8.9	0.99
b03d03		0.57	0.19	3.3	0.28	0.84	0.16	0.09	0.03	21	2.2
b03d10		1.8	0.18	3.4	0.29	0.93	0.07	0.03	0.01	60	5.9
b03d30		4.8	0.16	3.1	0.29	0.86	0.14	0.08	0.05	150	17
b05d01	1.00	0.11	0.11	3.4	0.29	0.72	0.28	0.20	0.09	4.5	0.51
b05d03		0.32	0.11	3.4	0.30	0.87	0.13	0.04	0.02	12	1.2
b05d10		0.98	0.09	3.5	0.31	0.98	0.02	0.01	0.01	30	2.8
b05d30		2.7	0.09	3.1	0.31	0.98	0.02	0.02	0.01	72	6.3
b08d01	1.33	0.05	0.05	3.3	0.28	0.66	0.34	0.25	0.12	2.3	0.25
b08d03		0.15	0.05	3.9	0.30	0.82	0.18	0.08	0.05	5.5	0.58
b08d10		0.49	0.05	3.6	0.33	0.99	0.01	0.01	0.00	13	1.1
b08d20		0.87	0.04	3.9	0.33	1.00	0.00	0.00	0.00	19	1.6
b10d01	1.11	0.09	0.09	3.4	0.29	0.72	0.28	0.20	0.12	3.9	0.47
b10d03		0.27	0.09	3.4	0.31	0.86	0.14	0.06	0.02	8.6	0.99
b10d10		0.90	0.09	3.3	0.33	0.98	0.02	0.02	0.01	23	2.5
b10d20		1.7	0.08	3.5	0.34	0.99	0.01	0.01	0.00	40	4.3
b15d01	1.25	0.06	0.06	3.7	0.27	0.63	0.37	0.27	0.14	2.9	0.35
b15d03		0.19	0.06	3.5	0.30	0.79	0.21	0.12	0.04	6.8	0.80
b15d10		0.62	0.06	3.5	0.34	0.98	0.02	0.01	0.01	16	1.6
b03d01-y10	0.60	0.21	0.21	3.4	0.26	0.62	0.38	0.32	0.25	9.0	1.1
b03d30-y10		4.1	0.14	3.4	0.28	0.86	0.14	0.10	0.08	130	13
b08d03-y10	1.33	0.16	0.05	3.7	0.29	0.75	0.25	0.20	0.16	6.0	0.71
b03d01-y15	0.60	0.20	0.20	3.3	0.27	0.69	0.31	0.25	0.19	8.9	1.1
b03d30-y15		4.3	0.14	2.8	0.28	0.85	0.15	0.12	0.09	120	13
b08d03-y15	1.33	0.15	0.05	3.7	0.30	0.76	0.24	0.18	0.10	5.5	0.64
b03d03-v10	0.60	0.28	0.28	5.6	0.27	0.55	0.45	0.41	0.34	13	1.2
b08d03-v10	1.33	0.31	0.10	5.5	0.31	0.77	0.23	0.20	0.15	12	0.83

The fraction of the initial disc mass ejected is shown in Fig. 2 for the baseline model sequence. This fraction ranges from |$4{{\ \rm per\ cent}}$| for the heaviest disc with |$M_{\rm bh}=8\, \mathrm{M}_\odot$| to |$21{{\ \rm per\ cent}}$| for the lightest disc around the |$M_{\rm bh}=3\, \mathrm{M}_\odot$| BH. The most important trend in mass ejection apparent from Table 2 is a monotonically decreasing ejected fraction with strength of gravity at the disc, ∝M_bh/R_d. For convenience, we define a disc compactness parameter as

$$\begin{eqnarray*} C_{\rm d} = \left(\frac{M_{\rm bh}}{5\, \mathrm{M}_\odot }\right)\left(\frac{50\, \textrm {km}}{R_{\rm d}}\right). \end{eqnarray*}$$

(3)

Fig. 2 also shows that the dependence of the ejected fraction with compactness is approximately linear, although at low BH masses the disc mass becomes an additional factor. This dependence on the strength of gravity has previously been documented in Fernández & Metzger (2013), and is also apparent from the results of Just et al. (2015), Fernández et al. (2017), and Fujibayashi et al. (2020), although this is the first time that it is sampled over an extended region of parameter space.

Figure 2.

Left-hand panel: Fraction of the initial disc mass ejected with positive Bernoulli parameter (equation 2) as a function of BH mass. Different symbols and colours correspond to different disc masses as labelled. The grey number above each symbol column corresponds to the disc compactness parameter C_d (equation 3). Right-hand panel: Fraction of the initial disc mass ejected as a function of disc compactness parameter, using the same colour and symbol coding as in the left-hand panel. The grey numbers above each symbol column denote the corresponding BH mass. The red dotted line is a linear fit to the ejected fraction for discs with |$M_{\rm d}=0.03\, \mathrm{M}_\odot$|⁠. GRMHD effects can enhance the ejected fraction and average velocity by up to a factor of ∼2 relative to hydrodynamic models (Fernández et al. 2019).

The second mass ejection trend in all models of the baseline sequence is a monotonic decrease in the ejected fraction with increasing disc mass at constant compactness. Fig. 2 shows that the strength of this dependence on disc mass is itself a function of disc compactness, with low-compactness discs being the most sensitive to the initial disc mass, while in high compactness systems this property has a smaller impact on the ejected fraction.

The physical origin of these trends in mass ejection can be traced back to the nature of the ejection mechanism. Most of the outflow is launched once weak interactions freeze out in the disc, removing the source of cooling. The fraction of the disc mass available to be ejected depends on how much has already been lost to accretion on to the BH by the time freezout occurs. This interplay is illustrated in Fig. 3, which shows the evolution of the mass accretion rate (⁠|$\dot{M}_{\rm isco}$|⁠) and cumulative mass accreted at the ISCO (M_acc), mass ejection rate at large radii (⁠|$\dot{M}_{\rm out}$|⁠), and the electron neutrino luminosity.

Figure 3.

Top: Mass accretion rate at the ISCO (⁠|$\dot{M}_{\rm isco}$|⁠, dotted lines), cumulative accreted mass at the ISCO (M_acc, dashes lines), and mass outflow rate in unbound material at r = 10⁹ cm (⁠|$\dot{M}_{\rm out}$|⁠, solid lines) as a function of time for models b08d03 (high-compactness), b03d01 (low compactness, low disc mass), and b03d30 (low compactness, high disc mass), as labelled. To facilitate comparison, the data from models b08d03 and b03d30 has been normalized to a disc mass of |$0.01\, \mathrm{M}_\odot$| (as in model b03d01). Bottom: Electron neutrino luminosity for models b08d03, b03d01, and b03d30, as labelled. The fraction of the disc ejected is related to the fraction of the disc accreted at the time when weak interactions freeze out.

In the model with the highest compactness (b08d03), accretion starts much earlier when measured in orbital times than in the lower compactness models. By the time weak interactions freeze out (steep plummet in neutrino luminosity at about 100 orbits) a significant fraction of the disc (⁠|$85{{\ \rm per\ cent}}$|⁠) has already been accreted to the BH. This earlier onset of accretion, despite having the same viscosity parameter, is due to the disc being closer to the ISCO. In terms of dimensionless numbers: r_isco/R_d = {0.26, 0.57} for C_d = {0.60, 1.33}, respectively.

For discs of the same compactness, the evolution of the accretion rate is very similar. Fig. 3 shows that mass ejection begins later in the more massive disc, which also takes longer time to reach freezout of weak interactions. More massive discs are more optically thick to neutrinos owing to their higher initial density, and for the same strength of viscosity, it takes more orbits for the density to decrease to a level where neutrino processes are no longer effective in cooling the disc. At the time when the neutrino luminosity reaches 10⁵⁰ erg s⁻¹, models b03d01 and b03d30 have accreted |$60{{\ \rm per\ cent}}$| and |$78{{\ \rm per\ cent}}$| of their initial disc masses, respectively.

The density dependence of the neutrino optical depth is not the only factor influencing the freezout time. Table 1 and equation (1) show that to keep the entropy constant, initial equilibrium discs with higher masses also have a higher internal energy content and therefore higher temperatures. This effect is stronger for lower compactness configurations. While disc material is more weakly bound, it can also radiate neutrinos at relevant levels for a longer time, and so this acts in the direction of delaying the onset of mass ejection. The post-merger entropy of the disc is thus an important parameter to keep track of in dynamical merger simulations given its effect on mass ejection efficiency.

Using the outflow trajectories from hydrodynamic simulations we can further analyse the energetics of mass ejection. In particular, we can quantify the strength and importance of different processes that change the heat content of the fluid: viscous heating, neutrino heating/cooling, and nuclear recombination of alpha particles. For each outflow trajectory, we compute the time integral of the local energy source term, yielding the heat gained or lost by the fluid element per unit mass:

$$\begin{eqnarray*} \Delta q_{\rm visc} = \int _0^{t_{\rm max,p}} \dot{q}_{\rm visc} \mathrm{ d}t \end{eqnarray*}$$

(4)

$$\begin{eqnarray*} \Delta q_\nu = \int _0^{t_{\rm max,p}} \dot{q}_{\rm \nu } \mathrm{ d}t \end{eqnarray*}$$

(5)

$$\begin{eqnarray*} \Delta q_\alpha = \frac{B_\alpha }{m_\alpha }\left[X_{\rm alpha}(t_{\rm max,p}) - X_\alpha (0) \right] \end{eqnarray*}$$

(6)

$$\begin{eqnarray*} \Delta q_{\rm net} = \Delta q_{\rm visc} + \Delta q_\nu + \Delta q_\alpha , \end{eqnarray*}$$

(7)

where |$\lbrace \dot{q}_{\rm visc},\dot{q}_{\rm \nu }\rbrace$| stand for viscous heating and neutrino heating/cooling, respectively, B_α/m_α ≃ 6.8 × 10¹⁸ erg g⁻¹ is the specific nuclear binding energy of alpha particles, X_α is the mass fraction of alpha particles, t_max,p is the maximum time for the particle evolution, either when it leaves the outer boundary of the computational domain or when the simulation ends, whichever is shorter, and Δq_i is the heat gained or lost from process i. The resulting quantities are shown in Table 3 for selected models.

Table 3.

Average time-integrated heat gain or loss per unit mass (equations 4–7), and change in Y_e (equations 10–17), due to various processes acting on tracer particles during the hydrodynamic evolution of selected models. The specific heat gain due to process i is normalized as |$\bar{\Delta } q_{i,19} = \bar{\Delta } q_i/(10^{19}\, \rm {erg\,g}^{-1})$|⁠, and |$\bar{\Delta }$| denotes average over all outflow particles. Each quantity is separately rounded for clarity, net sums match when including all significant digits.

Model	|$\bar{\Delta }q_{\rm net}/|b_{\rm ini}|$|	|$\bar{\Delta }q_{\rm net,19}$|	|$\bar{\Delta }q_{\rm visc,19}$|	|$\bar{\Delta }q_{\nu ,19}$|	|$\bar{\Delta }q_{\alpha ,19}$|	|$\bar{\Delta }Y_e^{\rm net}$|	|$\bar{\Delta }Y_e^{\rm em,\nu _e}$|	|$\bar{\Delta }Y_e^{\rm em,\bar{\nu }_e}$|	|$\bar{\Delta }Y_e^{\rm abs,\nu _e}$|	|$\bar{\Delta }Y_e^{\rm abs,\bar{\nu }_e}$|	|$\bar{\Delta }Y_e^{\rm em}$|	|$\bar{\Delta }Y_e^{\rm abs}$|
b03d01	0.81	2.1	2.6	−0.8	0.3	0.08	0.42	0.45	0.07	0.01	0.03	0.05
b03d03	0.79	1.7	2.5	−1.1	0.3	0.08	0.64	0.61	0.15	0.04	−0.03	0.12
b03d10	0.91	1.5	2.8	−1.6	0.3	0.09	1.00	0.91	0.25	0.07	−0.09	0.18
b03d30	1.12	1.3	2.9	−1.9	0.3	0.09	1.12	1.04	0.26	0.08	−0.08	0.18
b08d01	0.51	3.5	3.9	−0.7	0.3	0.09	0.34	0.41	0.03	0.01	0.07	0.02
b08d03	0.55	3.4	4.2	−1.1	0.4	0.10	0.53	0.58	0.07	0.02	0.05	0.05
b08d10	0.60	3.2	4.5	−1.7	0.4	0.11	0.81	0.83	0.13	0.04	0.02	0.09
b08d20	0.66	3.1	4.9	−2.2	0.4	0.13	1.10	1.08	0.21	0.06	−0.02	0.15
b03d01-y10	0.79	2.1	2.6	−0.8	0.3	0.16	0.37	0.48	0.07	0.01	0.11	0.05
b03d01-y15	0.78	2.1	2.6	−0.8	0.3	0.12	0.40	0.47	0.07	0.01	0.07	0.05
b03d01-v10	0.94	2.5	2.7	−0.5	0.3	0.07	0.22	0.25	0.05	0.02	0.04	0.03
b03d30-y10	1.16	1.4	3.1	−2.1	0.3	0.19	1.16	1.16	0.27	0.08	0.00	0.19
b03d30-y15	1.05	1.3	2.7	−1.8	0.3	0.13	1.05	1.01	0.25	0.07	−0.04	0.17
b08d03-y10	0.53	3.3	4.0	−1.1	0.3	0.19	0.46	0.60	0.07	0.02	0.14	0.05
b08d03-y15	0.55	3.4	4.3	−1.2	0.3	0.15	0.55	0.65	0.07	0.02	0.10	0.05
b08d03-v10	0.67	4.2	4.9	−1.0	0.4	0.12	0.40	0.48	0.06	0.02	0.08	0.04

Fig. 4 shows the distribution of net heat gained Δq_net (equation 7) by outflow particles for the same set of simulations shown in Fig. 3. Specific energies are shown in units of the initial Bernoulli parameter of the disc (equation 1). In all three models, the distribution shows a peak at around |$60{{\ \rm per\ cent}}$| of |b_ini|, with a tail towards high gain that is more extended for less compact models and higher disc masses. For a large fraction of the outflow, unbinding is not completely achieved by absorption of heat alone, which means that a significant part of the energy gain comes from adiabatic work. Discs with higher compactness are also less effective at absorbing heat, despite the fact that the absolute value of the net heat gain is higher (Table 3).

Figure 4.

Distribution of net specific heat gained by tracer particles due to source terms (equation 7) in the hydrodynamic evolution of models b08d03, b03d01, and b03d30, as labelled. The heat gain is normalized to the initial Bernoulli parameter in the disc (equation 1).

It is worth emphasizing that while discs with a higher compactness are more gravitationally bound in an absolute sense (b_ini ∝ M_bh/R_d in equation 1), they also undergo more net heating than lower compactness models (Table 3). Thus the absolute value of the gravitational potential (in a Newtonian sense), while correlating negatively with the fractional amount of heat absorbed, does not by itself explain the efficiency of mass ejection without accounting for how close the disc is to the ISCO radius.

For individual models, viscous heating and nuclear recombination contribute with net heating, while neutrinos primarily cool the disc (Table 3). This is illustrated in Fig. 5, which shows the distribution of the individual heating/cooling terms for the outflow from model b08d03. While there is a non-negligible fraction of particles for which neutrinos provide net heating, the magnitude of this heating (∼10¹⁷ erg g⁻¹) is dynamically negligible when compared to the dominant source terms (∼10¹⁹ erg g⁻¹).

Figure 5.

Distribution of specific heat gained or lost by tracer particles due to source terms (equations 4–6), as labelled, in the hydrodynamic evolution of model b08d03. The inset shows a snapshot of the initial positions of outflow particles in the disc, with grey/black particles corresponding to the shaded/dotted subsets of the neutrino and viscous heat gain histograms, respectively (see also Wu et al. 2016 for a larger plot of the initial particle distribution).

The integrated heat gain due to α particle recombination in the hydrodynamic simulation is subdominant compared to that from viscous heating. Despite its low global value, however, the heating due to α recombination is deposited over a short amount of time as the outflow is launched, and it can become comparable or even exceed the rate of viscous heating in this phase. The bulk of neutrino heating and cooling takes place before weak freezout. Fig. 5 shows that the nuclear recombination gain is narrowly distributed around the mean value 3–4 × 10¹⁸ erg g⁻¹ (Table 3). This value can be understood from the fact that upon expansion and cooling, all fluid elements achieve the maximal alpha particle mass fraction set by charge conservation, |$X_{\alpha ,\rm max} = 2Y_e$|⁠. The average electron fraction of the outflow 〈Y_e〉 ≃ 0.3 (Table 2), then sets the average amount of energy gained.

Fig. 5 also shows that the tracer particles in model b08d03 follow bimodal distributions of viscous heating and neutrino cooling. This bimodality can be traced back to the initial positions of the particles in the disc. These particles originate from two regions: (1) the equatorial plane of the disc, and (2) regions above the equatorial plane around the initial density peak (see also Wu et al. 2016). The first group of particles experiences stronger viscous heating and neutrino cooling, while the second group experiences heating or cooling with less intensity. The latter group includes all the particles that experience net neutrino heating. The initial position of the particles is related to the way the disc overturns in the poloidal direction due to viscous heating, and may differ from that obtained when MHD turbulence transports angular momentum.

The fraction of the disc mass ejected and average outflow velocity are nearly insensitive to the initial Y_e of the disc except for very massive discs in low-compactness systems, where differences of a few per cent of the disc mass can arise (Table 3). Changes in the viscosity parameter, on the other hand, result in important changes to both ejected fraction and outflow velocity. The average outflow radial velocity is in the range 0.03–0.04 c for all models that use α = 0.03, while this average increases to 0.05–0.06 c for the models with α = 0.1. The fraction of the disc mass ejected increases by |${\sim}30{{\ \rm per\ cent}}$| for the model with low compactness (b03d01-v10) and by a factor of 2 for the high compactness model (b08d03-v10). Table 3 shows that the increase in the viscosity parameter results in more viscous heating in the high-compactness model b08d03-v10 and less neutrino cooling in the least compact model b03d01-v10, in both cases increasing the net heat gain of the outflow.

3.2 Outflow composition

Table 2 shows that the average electron fraction of the disc outflow is in the range 0.25–0.35 for all simulated models. These values are just above the critical transition at which the nucleosynthesis changes from rich to poor in elements with mass numbers A > 130 (e.g. Kasen, Fernández & Metzger 2015; Lippuner & Roberts 2015). For the purposes of predicting kilonova properties, the mass fraction in lanthanides (57 ≤ Z ≤ 72, with Z the atomic number) and actinides (89 ≤ Z ≤ 104) is the most important, since these species have an outsize influence on the ejecta opacity – and therefore on the kilonova colour, luminosity, and duration – given their atomic complexity (Barnes & Kasen 2013; Kasen, Badnell & Barnes 2013; Tanaka & Hotokezaka 2013; Fontes et al. 2015; Tanaka et al. 2020).

We therefore refine our diagnostic of the outflow composition by analysing the output of post-processed tracer particles with SkyNet (Section 2.3). For each simulation, we report in Table 2 the fraction of the outflow particles with lanthanide and actinide mass fraction X_La + X_Ac < 10⁻⁴ and define it as the ‘blue mass’

$$\begin{eqnarray*} M_{\rm blue} = \int _{0}^{10^{-4}}\frac{\mathrm{ d}M_{\rm ej}}{\mathrm{ d}(X_{\rm La}+X_{\rm Ac})}\, \mathrm{ d}(X_{\rm La}+X_{\rm Ac}). \end{eqnarray*}$$

(8)

The value of 10⁻⁴ is small enough that the outflow opacity is indistinguishable from that dominated by iron-group elements (Kasen et al. 2013). Likewise, we define three ‘red’ masses

$$\begin{eqnarray*} M_{-k} = \int _{10^{-k}}^1 \frac{\mathrm{ d}M_{\rm ej}}{\mathrm{ d}(X_{\rm La}+X_{\rm Ac})}\, \mathrm{ d}(X_{\rm La}+X_{\rm Ac}). \end{eqnarray*}$$

(9)

using −k = {−4, −3, −2}. These four numbers provide a description of the incidence of heavy r-process elements in the ejecta that, while coarser than detailed abundances, is more informative than just the electron fraction. Fig. 6 shows the distribution of mass fractions at 1 d in the outflow particles from models b08d01-b08d20 (highest compactness), along with the regions encompassed by M_blue and M_−k, and the isotopic abundances for comparison. At 1 d, the mass fractions of actinides are in general 10 times smaller than that of lanthanides in our models, and thus we lump them together when discussing composition effects (for more detailed studies on actinide production; see Eichler et al. 2019 or Holmbeck et al. 2019).

Figure 6.

Top: Distribution of lanthanide and actinide mass fractions at 1 d, for nuclear-network-processed particles from models b08d01-b08d20, as labelled. Since each particle represents an equal mass element in the disc, a sum over the bins yields the fraction of the mass with a given Lanthanide and Actinide fraction (M_blue and |$\lbrace M_{-4},M_{-3},M_{-2}\rbrace$|⁠, cf. equations (8)–(9) and Table 2). The histograms continue to mass fractions lower than 10⁻¹² with similar slope and were truncated for clarity. Bottom: Isotopic abundances at 1 d for nuclear-network-processed particles from models b08d01-b08d20 as labelled. Abundances are normalized such that their mass fractions add up to unity. The Solar system r-process abundances from Goriely (1999) are normalized to model b08d03 at A = 130. Note that the dynamical ejecta is rich in elements with A > 130, and in combination with the disc outflow can supply the entire range of r-process elements (e.g. Just et al. 2015).

The most robust composition trends from Table 2 are that (1) the majority of the disc outflow mass is lanthanide poor (M_blue), and that (2) the fraction of lanthanide poor material is a monotonically increasing function of the disc mass, for constant BH mass. The dependence on disc mass is clearly illustrated in Fig. 6: the fraction of particles with high lanthanide abundance is a steep function of the disc mass, as is the abundance of elements with A > 130 (this dependence has also been reported by Just et al. 2015 and Fujibayashi et al. 2020). While there is some dependence on the compactness C_d at constant disc mass, this dependence is not fully monotonic, and is weaker than the dependence on disc mass at fixed compactness. Therefore, this variation with compactness is more likely to be dependent on the details of how the disc evolution is modelled.

The trend of more lanthanide poor ejecta with increasing disc mass is apparent from the distribution of electron fraction (Fig. 7). Despite having a very similar average value, the Y_e distribution of model b08d20 extends to significantly higher values than in model b08d01. Naively, one would expect that higher electron fraction is associated with a higher equilibrium Y_e arising from a higher abundance of positrons given higher entropies (e.g. Beloborodov 2003). However, the entropy distributions of models b08d01 and b08d20 show that the average entropy decreases with higher disc masses. The lower entropy can be understood from the fact that the temperature is primarily set by the strength of the gravitational potential once accretion is established, and is similar in both models. On the other hand, the densities are higher at any given time for a higher disc mass, with a correspondingly lower entropy than at lower disc masses.

Figure 7.

Mass histograms of electron fraction (left-hand panel) and entropy (right-hand panel) for models b08d01 (low disc mass), b08d20 (high disc mass), and b08d03-v10 (high viscosity), obtained by the end of the hydrodynamic simulation at r = 10⁹ cm. To facilitate comparison, histogram masses have been normalized to that of model b08d01. The hatched region is the subset of the histogram of model b08d03-v10 for times t < 1 s .

We can further analyse the evolution of the electron fraction distribution of the outflow by computing the time-integrated contribution of the different neutrino processes that change Y_e. The net rate of change of Y_e due to neutrino emission and absorption arises from the following reactions:

$$\begin{eqnarray*} \Gamma ^{\rm em, \nu _e}: e^- + p \rightarrow n + \nu _e \end{eqnarray*}$$

(10)

$$\begin{eqnarray*} \Gamma ^{\rm em, \bar{\nu }_e}: e^+ + n \rightarrow p + \bar{\nu }_e \end{eqnarray*}$$

(11)

$$\begin{eqnarray*} \Gamma ^{\rm abs, \nu _e}: \nu _e + n \rightarrow e^- + p \end{eqnarray*}$$

(12)

$$\begin{eqnarray*} \Gamma ^{\rm abs, \bar{\nu }_e}: \bar{\nu }_e + p \rightarrow e^+ + n \end{eqnarray*}$$

(13)

$$\begin{eqnarray*} \Gamma ^{\rm em} = \Gamma ^{\rm em, \bar{\nu }_e} - \Gamma ^{\rm em, \nu _e} \end{eqnarray*}$$

(14)

$$\begin{eqnarray*} \Gamma ^{\rm abs} = \Gamma ^{\rm abs, \nu _e} - \Gamma ^{\rm abs, \bar{\nu }_e} \end{eqnarray*}$$

(15)

$$\begin{eqnarray*} \Gamma ^{\rm net} = \Gamma ^{\rm em} + \Gamma ^{\rm abs} \end{eqnarray*}$$

(16)

with the net rate setting the overall evolution of the electron fraction in the hydrodynamic simulation: |$\partial Y_e/\partial t + \mathbf {v}\cdot \nabla Y_e = \Gamma ^{\rm net}$|⁠. For each tracer particle, we compute a separate time integral for each of the rates above, obtaining a contribution to the change in electron fraction:

$$\begin{eqnarray*} \Delta Y_e^i = \int _0^{t_{\rm max,p}} \Gamma ^i \mathrm{ d}t, \end{eqnarray*}$$

(17)

where again t_max,p is the maximum time of the particle in the simulation.

Table 3 shows that the rates of neutrino emission dominate over neutrino absorption for all models, in line with the overall dominance of neutrino cooling over neutrino heating (cf. Fig. 5). However, the change in Y_e is set by differences in the rates of neutrino/antineutrino emission and absorption. If the two emission rates are closer in magnitude than the two absorption rates, the latter can dominate the change in Y_e despite being smaller in magnitude than the former.

Fig. 8 shows the average change in Y_e for tracer particles as a function of disc mass in the baseline sequence with |$M_{\rm bh}=8\, \mathrm{M}_\odot$| along with the breakdown of this change between emission and absorption of electron neutrinos/antineutrinos. At low disc masses, emission processes dominate the change in Y_e, with decreasing relative importance with increasing disc mass. Emission processes act towards driving Y_e to its local equilibrium value set by the entropy, and this equilibrium is lower at higher disc masses given the lower entropy (Fig. 7). At the highest disc mass (model b08d20), this change in Y_e due to emission is even negative.

Figure 8.

Average electron fraction change in particles due to neutrino emission or absorption (equations 10–16 and 17, Table 3), during the hydrodynamic evolution of models b08d01-b08d20, as labelled. The black horizontal lines denote the net change in Y_e when including all processes. The relative importance of neutrino absorption increases as the disc mass increases, all else being equal, accounting for the higher Y_e of the outflow with lower overall entropies (Fig. 7).

Absorption, on the other hand, is set by the ambient flux of incident neutrinos and the mass fractions of neutrons and protons. At higher disc masses, neutrino/antineutrino luminosities are higher and stay high for a longer time (Fig. 3) thus increasing the ambient neutrino flux and thus the magnitude of absorption terms. The asymmetry in the neutron–proton mass fraction then results in different absorption efficiencies and a net change in Y_e which always acts in the direction of increasing it (because the reaction in equation 12 occurs more frequently). Fig. 8 shows that at high disc masses, absorption dominates the evolution of Y_e in the high-compactness models shown. Table 3 shows that this trend is also present in models with the lowest compactness (b03d01-b03d30), with an even stronger effect of absorption on Y_e than in the high compactness sequence.

Models with lower initial Y_e = {0.10, 015} eject a higher proportion of lanthanides and actinides, as expected, although the change is at a |${\lesssim}10{{\ \rm per\ cent}}$| level by mass relative to our baseline parameters (Fig. 9). While the Y_e distributions extend to slightly lower electron fractions than the default models, the bulk of the outflow has Y_e > 0.2 in all cases. Table 3 shows that these changes are driven primarily by neutrino emission processes, which adjust to provide the required change in Y_e towards equilibrium.

Figure 9.

Distribution of lanthanide and actinide mass fractions at 1 d, for nuclear-network-processed particles from models b08d03 (baseline), b08d03-y10 and b08d03-y15 (varying initial Y_e), and b08d03-v10 (high viscosity), as labelled. The bin size is the same is in the top panel of Fig. 6. The lowest bin contains all the particles with X_La + X_ac < 10^−14.5.

Models with higher viscosity have a similar average Y_e than their low-viscosity counterparts, but a significantly higher fraction of material rich in lanthanide and actinides. The electron fraction distribution of model b08d03-v10 has a tail to low Y_e that extends well below the lower limit of the distribution of model b08d01. Fig. 7 shows that the material with the lowest Y_e is ejected at the earliest times in the high-viscosity model, in line with the general expectation that the faster the evolution of a disc, the stronger the sensitivity of the outflow composition to initial conditions. This is consistent with the results of GRMHD simulations, which show even stronger memory of initial conditions given their faster evolution (Fernández et al. 2019).

3.3 Implications for EM counterpart searches

The key question we are interested in is how does the disc outflow contribute to a kilonova transient. The answer depends on the amount of mass ejected and its velocity, its composition, as well as how efficiently does the radioactive heating thermalize (e.g. Metzger 2019).

The mass ejected, and the composition to a lesser extent, determines how much radioactive power is available for a kilonova. Fig. 10 shows the total radioactive heating luminosity at 1 d as a function of compactness parameter C_d and initial disc mass. For fixed compactness, the total radioactive heating is proportional to the ejected mass, since the average radioactive heating per unit mass is close to 2 × 10¹⁰ erg g⁻¹ s⁻¹ for most models, given their similar composition. The dependence of ejected fraction on compactness results in an additional variation of a factor ∼5 between the least and most compact models, for fixed disc mass.

Figure 10.

Total radioactive heating luminosity at 1 d as a function of disc compactness (equation 3), for various disc masses, as labelled. The BH mass is shown in grey under each symbol column. The crosses denote the heating rates interpolated to the median disc masses from Fig. 1. The total radioactive heating rate is an upper limit to the bolometric luminosity of the kilonova, being subject to thermalization efficiency and radiative transfer effects.

The raw radioactive heating at 1 d ranges from 2 × 10⁴⁰ erg s⁻¹ for the lightest and most compact disc (b08d01) to 1.5 × 10⁴² erg s⁻¹ for the heaviest disc in the least compact configuration (b03d30). Thermalization efficiency can result in a decrease in these values by a factor ∼2 (Barnes et al. 2016; Kasen & Barnes 2019; Waxman, Ofek & Kushnir 2019; Hotokezaka & Nakar 2020), while radiative transfer effects (dependent on the opacity and thus on composition) set whether this power can readily escape the ejecta or is trapped until later times. Table 2 shows that the radioactive power at 1 week is about 10 times smaller than that at 1 d for most models.

If we take the median disc masses from Fig. 1 as representative values for each BH mass and interpolate the ejected masses from Table 2, we obtain disc outflow masses |$\lbrace 2.00, 0.70, 0.22, 0.25, 0.14\rbrace \times 10^{-2}\, \mathrm{M}_\odot$| for BH masses |$\lbrace 3, 5, 8, 10, 15\rbrace \, \mathrm{M}_\odot$|⁠, respectively. These values are subject to enhancement by a factor ≲2 if GRMHD effects were to be included (Christie et al. 2019; Fernández et al. 2019).

In the case of GW190425, for which our lowest BH mass model would be applicable, the median disc outflow mass would be similar to the total ejecta from GW170817 within a factor of two, and hence it would have been detected with good sky coverage (Kyutoku et al. 2020). A BH–NS merger with a low-mass BH and high-mass NS is most efficient at tidally disrupting the NS and most inefficient at producing dynamical ejecta, with most mass ejection coming from the disc (Foucart et al. 2019). In contrast, two massive NS that collapse promptly to a BH are the least efficient configuration for ejecting mass and forming a disc (e.g. Radice et al. 2018, but see Kiuchi et al. 2019; Bernuzzi et al. 2020 for the case of an asymmetric mass ratio NS–NS merger generating a more massive BH accretion disc than a symmetric binary).

Regarding the BH–NS merger candidate GW190814, which had a much smaller localization area and deeper EM coverage relative to other events, constraints on the total mass ejection are in the range |$0.02-0.1\, \mathrm{M}_\odot$| depending on viewing angle, composition, and distance (e.g. Andreoni et al. 2020; Kawaguchi, Shibata & Tanaka 2020b; Morgan et al. 2020; Vieira et al. 2020). Our results indicate that, with the exception of a very low-mass BH and/or very high disc masses, most BH–NS merger systems would not have generated sufficient disc outflow for a detectable kilonova.

An additional factor influencing the kilonova appearance is the relative masses of the disc and the dynamical ejecta. In a BH–NS merger, the dynamical ejecta is produced along the equatorial plane, and it blocks only part of the viewing angles. It is expected to be very rich in lanthanides and thus have a high opacity that blocks the light from the disc towards most equatorial directions. The results of Fernández et al. (2015b, 2017) – who studied the combined long-term evolution of these two components mapped from dynamical merger simulations – show that the bulk of fallback accretion mixes in with the disc before weak freezout occurs and the disc outflow is launched. The net outflow has a very similar composition as the dynamical ejecta and disc when evolved separately, but with an added component that has intermediate electron fraction values. The expected kilonova colour has an important dependence on viewing angle, as the bluer disc emission would only be detectable from directions not obscured by the dynamical ejecta (see, e.g. Bulla et al. 2019; Darbha & Kasen 2020; Kawaguchi, Shibata & Tanaka 2020a;Korobkin et al. 2020 for more recent work on viewing angle dependencies of different ejecta configurations).

In hydrodynamic simulations of BH accretion discs, the spatial dependence of the composition is quite generic, with the highest Y_e material being ejected first along intermediate latitudes, and then wrapping around the outermost edge of the disc outflow (Fernández et al. 2015a). This configuration would guarantee the existence of a blue spike at early times (albeit a faint one) if the merger remnant is viewable from polar latitudes. All of our models are such that at least |$50{{\ \rm per\ cent}}$| of the disc outflow is lanthanide-poor, and in some cases this fraction can reach even |$100{{\ \rm per\ cent}}$| (Table 2). GRMHD models show, however, that magnetic fields begin to eject matter much earlier than weak freezout, thus adding material that has not been sufficiently processed by weak interactions from the initial post-merger composition (e.g. Siegel & Metzger 2017; Christie et al. 2019; Fernández et al. 2019; Miller et al. 2019). This will likely increase the lanthanide-rich fraction at the leading edge of the disc outflow and thus modify the colour of the disc-contributed kilonova.

4 SUMMARY AND DISCUSSION

We have performed axisymmetric hydrodynamic simulations of the viscous evolution of accretion discs formed in BH–NS mergers. Our models include the effects of neutrino emission and absorption on the outflow composition, and the results are post-processed with a nuclear reaction network for a more precise diagnostic of the nucleosynthesis yield. These hydrodynamic models provide a lower limit to the amount of mass in the disc outflow relative to 3D GRMHD models, and hence provide a lower limit to the contribution of disc outflows to the kilonova transient. Our simulations cover, for the first time, a large fraction of the plausible parameter space of BH and disc masses expected from these mergers (Fig. 1). Our main results are the following:

1. The fraction of the initial disc mass ejected as an unbound outflow has an approximately inverse linear scaling with the initial compactness of the disc, and can vary by a factor of ∼4 (Fig. 2 and Table 2). While this dependence on compactness was implicit in previous work, this is the first time that it is systematically probed over a wide parameter space.

   The origin of this dependence can be traced back to the earlier onset of accretion in more compact discs, as they are located closer to the ISCO. Compared to lower compactness discs, a higher fraction of the initial mass is accreted by the time weak interactions freeze out and the outflow is launched (Fig. 3).

2. At constant compactness, the fraction of the disc mass ejected decreases for higher disc masses. This effect is related to the longer time to reach weak freezout in more massive discs, which delays the onset of mass ejection to a time when more mass has been accreted to the BH (Fig. 3). The dependence on initial disc mass is weaker for higher compactness systems (Fig. 2 and Table 2). The initial density and entropy of the disc are the key variables regulating this mass dependence of the ejection efficiency (Section 3.1).

3. The disc outflow is more lanthanide/actinide-poor for higher disc masses (Fig. 6) at constant compactness (this trend has also been found in previous studies). This effect can be traced back to neutrino absorption becoming more important relative to emission in increasing Y_e for more massive discs (Fig. 8 and Table 3). Stronger absorption counteracts the action of neutrino emission in lowering Y_e given the lower entropy of the outflow from more massive discs (Fig. 7).

   While our disc outflows are |$50{-}100{{\ \rm per\ cent}}$| lanthanide-poor by mass, magnetically driven contributions – not included here – can add a significant amount of lanthanide-rich matter, hence the net composition of disc outflows requires simulations with more complete physics for reliable predictions.

4. The ejected fraction and final composition are sensitive to the viscosity parameter of the disc, as known from previous work, with more mass ejected, at higher velocities, and with a higher lanthanide-rich fraction for higher viscosity parameter (Table 2). The most neutron-rich material is produced at the earliest times in the simulation (Fig. 7) and is thus related to the shorter evolution time of these discs.

5. In most cases, the initial Y_e of the disc has a negligible effect on the amount of mass ejected, with the exception of massive discs in low-compactness systems, where the effect modifies the ejection efficiency by a few per cent of the disc mass (similar to numerical resolution). Hence, the ejected fraction is mostly robust to variations in the initial composition (Table 2). The final composition does depend on the initial conditions, with corrections of the order of |$10{{\ \rm per\ cent}}$| to the fraction of the outflow mass that is lanthanide-rich (Fig. 9)

6. The range of ejecta masses from the disc outflow can result in a range of two orders of magnitude in raw radioactive luminosity over the BH–NS parameter space probed (Fig. 10). Except for very low-mass BHs and/or very high disc masses, most BH–NS mergers should generate disc outflows that are below constraints for the total ejecta mass from the BH–NS merger candidate GW190814 (Section 3.3).

Our results are consistent with previous hydrodynamic models of BH accretion discs. Model t-a80-hr from Fernández et al. (2015a) is equivalent to our model b03d01 but evolved until 3000 orbits. By that time, the total and unbound (positive energy) mass ejected in their model are |$19{{\ \rm per\ cent}}$| and |$12{{\ \rm per\ cent}}$|⁠, respectively, while here we obtain |$21{{\ \rm per\ cent}}$| and |$14{{\ \rm per\ cent}}$|⁠, respectively. The |${\sim}10{{\ \rm per\ cent}}$| difference can be attributed to the lower ambient and floor of density used here, and on improvements in the neutrino implementation as reported in Lippuner et al. (2017). Similarly, model Fdisk of Fernández et al. (2017) falls in between models b08d03 and b08d10 in terms of compactness, but has a slightly higher BH spin (0.86). The higher ejected fraction in their model (⁠|$8{{\ \rm per\ cent}}$|⁠) can be partially accounted for by the higher value of the BH spin, and also by a more extended initial density distribution with radius in the torus mapped from the merger simulation (covering a wider range in compactness)

We have made specific choices for the disc entropy (Section 2.1), which can have implications for the sensitivity of the ejected fraction to initial disc mass. Similarly, the choice of viscosity parameter is on the low end of values that bracket the amount of angular momentum transport seen in GRMHD simulations of equivalent accretion discs (Fernández et al. 2019). More realistic values for these two parameters must come from direct mapping of the outcome of GRMHD simulations of BH–NS mergers that include neutrino transport, where the magnetic field geometry and strength replaces the viscosity parameter. A direct mapping would also avoid the need to make well-motivated but ultimately arbitrary choices for an initial torus radius, which is required by an equilibrium initial condition and directly enters the compactness parameter (equation 3). Mapping from merger simulations would also inform the initial Y_e distribution which, while not crucial for determining the amount of mass ejected, has implications for the outflow composition at the level of detail needed for accurate predictions of kilonovae and nucleosynthesis yields. Compact object merger simulations that account for both MHD and neutrino effects are few and only implement the latter via leakage schemes (e.g. Neilsen et al. 2014; Most et al. 2019). Simulations with more advanced (e.g. M1) neutrino transport do not yet include MHD effects.

Our results suggest that GRMHD simulations carried out with no neutrino absorption (e.g. Siegel & Metzger 2018; Christie et al. 2019; Fernández et al. 2019) can give reasonable Y_e distributions only for very low mass discs |${\lt}0.01\, \mathrm{M}_\odot$|⁠, for which the absorption contribution to Y_e is subdominant. For the configuration studied by Siegel & Metzger (2018), Fernández et al. (2019), and Christie et al. (2019) (⁠|$M_{\rm bh}=3\, \mathrm{M}_\odot$| and |$M_{\rm d}=0.03\, \mathrm{M}_\odot$|⁠), we find that absorption is already more important than emission in setting the Y_e distribution (Table 3), in line with the results of Miller et al. (2019) who showed a significant increase in the Y_e of the outflow when including neutrino absorption.

Further GRMHD studies of BH accretion discs over a wide region of parameter space are needed to determine whether the ejected fraction has the same dependence on compactness as in pure hydrodynamic models, and whether the fraction of the outflow that is lanthanide-rich is significantly larger than what we find here. Both of these questions are crucial to improve predictions of EM counterparts to BH–NS and NS–NS sources, and require (1) models evolved for long time-scales, (2) realistic initial field strengths, geometries, entropies, and electron fractions, and (3) inclusion of neutrino absorption. Such simulations remain challenging given current algorithms and computational resources.

ACKNOWLEDGEMENTS

We thank the anonymous referee for helpful comments that improved the presentation of the paper. RF acknowledges support from the National Sciences and Engineering Research Council (NSERC) of Canada through Discovery Grant RGPIN-2017-04286, and from the Faculty of Science at the University of Alberta. FF gratefully acknowledges support from the U.S. National Science Foundation (NSF) through grant PHY-1806278, from the National Aeronautics and Space Administration (NASA) through grant80NSSC18K0565, and from the U.S. Department of Energy CAREER grant DE-SC0020435. This work was supported by the U.S. Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001). This work was assigned report number LA-UR-20-23877. The software used in this work was in part developed by the U.S. Department of Energy NNSA-ASC OASCR Flash Center at the University of Chicago. This research was enabled in part by support provided by WestGrid (www.westgrid.ca), the Shared Hierarchical Academic Research Computing Network (SHARCNET, www.sharcnet.ca), Calcul Québec (calculquebec.ca), and Compute Canada (www.computecanada.ca). Computations were performed on Graham, Cedar, and Béluga. This research also used storage resources of the U.S. 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 (repository m2058). Graphics were developed with matplotlib (Hunter 2007).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.