# Igniting weak interactions in neutron star postmerger accretion disks

Soumi De Computational Physics and Methods, Los Alamos National Laboratory,
Los Alamos, NM 87545, USA  Department of Physics, Syracuse University,
Syracuse, NY 13244, USA  Kavli Institute for Theoretical Physics, University
of California, Santa Barbara, CA 93106, USA  Daniel M. Siegel Perimeter
Institute for Theoretical Physics, Waterloo, Ontario, Canada, N2L 2Y5
Department of Physics, University of Guelph, Guelph, Ontario, Canada, N1G 2W1
Kavli Institute for Theoretical Physics, University of California, Santa
Barbara, CA 93106, USA

###### Abstract

The merger of two neutron stars or of a neutron star and a black hole
typically results in the formation of a post-merger accretion disk. Outflows
from disks may dominate the overall ejecta from mergers and be a major source
of r-process nuclei in our universe. We explore the parameter space of such
disks, their outflows, and r-process yields by performing three-dimensional
general-relativistic magnetohydrodynamic (GRMHD) simulations with weak
interactions and approximate neutrino transport. We discuss the mapping
between initial binary parameters and the parameter space of resulting disks,
chiefly characterized by their initial accretion rate. We demonstrate the
existence of an ignition threshold for weak interactions at around $\sim
10^{-3}M_{\odot}\mathrm{s}^{-1}$ for typical parameters by means of analytic
calculations and numerical simulations. We find a degenerate, self-regulated,
neutrino-cooled regime above the threshold and an advection dominated regime
below the threshold. Excess heating in the absence of neutrino cooling below
the threshold leads to $\gtrsim 60\%$ of the initial disk mass being ejected
in outflows, with typical velocities $\sim(0.1-0.2)c$, compared to $\lesssim
40\%$ at $\sim(0.1-0.15)c$ above the threshold. While disks below the
threshold show suppressed production of light r-process elements, disks above
the threshold can produce the entire range of r-process elements in good
agreement with the observed solar system abundances. Disks below the ignition
threshold may produce an overabundance of actinides seen in actinide-boost
stars. As gravitational-wave detectors start to sample the neutron star merger
parameter space, different disk realizations may be observable via their
associated kilonova emission.

accretion, accretion disks – magnetohydrodynamics – gravitation – neutrinos –
nuclear reactions: nucleosynthesis, abundances – stars: neutron – stars: black
holes

##  1 Introduction

The gravitational-wave observatories advanced LIGO and Virgo, recently joined
by KAGRA, are starting to make routine detections of compact binary mergers
(Abbott et al., 2019, 2020b; Nitz et al., 2019a,
b)111https://gracedb.ligo.org/superevents/public/O3/. The first neutron-star
(NS) merger GW170817 (Abbott et al., 2017a) and a second NS merger candidate
GW190425 (Abbott et al., 2020a) already probe drastically different parts of
the NS binary parameter space. While the component masses of GW170817 are
typical of previously known Galactic double NS systems, the total mass of
GW190425 of $\sim\\!3.4\,M_{\odot}$—if indeed a double NS system—represents an
outlier by $5\sigma$ with respect to the known Galactic distribution of binary
NSs (Abbott et al., 2020a) (or may represent a peculiar neutron star - black
hole or binary black hole system otherwise). Furthermore, while GW170817 was
followed by electromagnetic counterparts across the entire electromagnetic
spectrum (Abbott et al., 2017b), including the first unambiguous detection of
a ‘kilonova’ (Li & Paczyński, 1998; Kulkarni, 2005; Metzger et al., 2010;
Metzger, 2019), GW190425 did not lead to the detection of electromagnetic
counterparts, which may be due to intrinsically dim emission, the large
distance of $\sim\\!160$ Mpc to the binary, and poor sky localization of
$\sim\\!8,300\,\mathrm{deg}^{2}$ (Abbott et al., 2020a; Antier et al., 2020;
Lundquist et al., 2019; Hosseinzadeh et al., 2019; Coughlin et al., 2020).

Six decades after Burbidge et al. (1957) and Cameron (1957) realized that
about half of the cosmic abundances of nuclei heavier than iron are created by
rapid neutron capture onto light seed nuclei (‘the r-process’), GW170817
provided the first direct observation of cosmic synthesis of such elements.
The quasi-thermal emission in the ultraviolet, optical, and near-infrared was
consistent with a kilonova, i.e., with being powered by the radioactive decay
of r-process nuclei synthesized in the merger ejecta (see, e.g., Metzger 2019;
Siegel 2019 for overiews discussing the interpretation of the event).

Binary neutron star (BNS) or NSBH mergers that lead to tidal disruption of the
NS outside the innermost stable circular orbit can give rise to neutron-rich
ejecta material conducive to the r-process in a number of ways. This includes
dynamical ejecta with tidal and shock-heated components (Ruffert et al., 1997;
Rosswog et al., 1999; Oechslin et al., 2007; Hotokezaka et al., 2013a, b),
neutrino-driven and magnetically driven winds from a (meta-)stable remnant
(Dessart et al., 2009; Siegel et al., 2014; Ciolfi et al., 2017, 2019), and
outflows from a post-merger accretion disk (Metzger et al., 2009, 2008a;
Fernández & Metzger, 2013; Just et al., 2015; Siegel & Metzger, 2017). The
details and relative importance of these ejecta components depend on binary
parameters and the unknown equation of state (EOS) of nuclear matter at
supranuclear densities (see Sec. 2.1 for more details). In particular, the
bulk of the GW170817 ejecta, specifically the material giving rise to the
‘red’ lanthanide-bearing kilonova component, is most naturally explained by
outflows from a post-merger accretion disk (Kasen et al., 2017; Siegel &
Metzger, 2017), while the origin of the ‘blue’ emission in GW170817 may be due
to a different source or combination of sources (Siegel & Metzger 2018a;
Fernández et al. 2019; Miller et al. 2019b; Nedora et al. 2019; Metzger et al.
2018; Ciolfi & Vijay Kalinani 2020; see, e.g., Metzger 2019 and Siegel 2019
for more discussion). Due to issues with chemical evolution (see Siegel 2019
for a brief summary) and the possibility of other sources such as magneto-
rotational supernovae (Winteler et al., 2012; Halevi & Mösta, 2018) and
collapsars (Siegel et al., 2019) contributing significantly, it still remains
an open question whether NS mergers are the dominant source of r-process
elements.

Post-merger accretion disks form as a significant amount of merger debris
circularized around the remnant. Such disks also provide a promising central
engine to generate collimated relativistic jets needed to generate short
gamma-ray bursts (Aloy et al., 2005; Shibata & Taniguchi, 2006; Paschalidis et
al., 2015; Ruiz et al., 2016). Numerical studies of the evolution of such
accretion disks exist with various levels of approximation and computational
complexity (Fernández & Metzger, 2013; Just et al., 2015; Siegel & Metzger,
2017; Fernández et al., 2019, 2020; Christie et al., 2019a; Miller et al.,
2019b; Fujibayashi et al., 2017, 2020; Nedora et al., 2019; Just et al., 2021;
Li & Siegel, 2021). Recent studies indicate that about 20–40% of the disk
material may be unbound into powerful neutron-rich outflows, which makes them
a strong source of kilonova emission and a potentially dominant source of
r-process ejecta across a wide region in NS binary parameter space (see Sec.
2.1 for a discussion). However, due to the computational complexity and cost,
to date there exist only a few specific simulations that take all necessary
physical ingredients—general relativity, magnetic fields, weak interactions,
neutrino transport—into account (Siegel & Metzger, 2017; Siegel & Metzger,
2018a; Fernández et al., 2019; Christie et al., 2019a; Miller et al., 2019b;
Li & Siegel, 2021) (see Secs. 2.1 and Sec. 3 for more details), and most of
the post-merger parameter space and associated physics remains largely
unexplored.222During the writing of this paper, Fernández et al. (2020)
published a first survey of disk models throughout the parameter space
employing two-dimensional Newtonian disk simulations with an
$\alpha$-viscosity to mimic angular momentum transport in the absence of
magnetic fields.

This paper presents the first exploration of the parameter space of neutrino-
cooled accretion disks across two orders of magnitude in accretion rates and
disk masses by means of self-consistent three-dimensional general-relativistic
magnetohydrodynamic (GRMHD) simulations with weak interactions. This study is
conducted in anticipation of future detections of binary mergers by LIGO,
Virgo, and Kagra, which will soon sample the NS merger parameter space with
many more detections. This study focuses on the transition across an ignition
threshold for weak interactions, which distinguishes qualitatively distinct
states and properties of such accretion disks and related parts of the neutron
star binary parameter space. We begin by elaborating on this ignition
threshold (Sec. 2.1) and by relating post-merger disks to NS binary parameters
and future detections (Sec. 2.2). A brief overview of numerical methods is
provided in Sec. 3. Section 4 summarizes our results, including global and
local disk properties as well as r-process nucleosynthesis. Discussion and
conclusions are presented in Sec. 5.

##  2 Physical model

###  2.1 NS post-merger disks: ignition threshold

Compact accretion disks while optically thick to photons may be cooled via
neutrino emission (Popham et al., 1999; Narayan et al., 2001; Di Matteo et
al., 2002; Beloborodov, 2003; Kohri et al., 2005; Chen & Beloborodov, 2007;
Kawanaka & Mineshige, 2007). At sufficiently high midplane density and
temperature, weak interaction rates become high relative to the rate of radial
advection of thermal energy. This gives rise to two limiting states of such
disks: $(i)$ weak interactions are important and the disk is neutrino-cooled
(predominantly via electron and positron capture: $e^{-}+p\rightarrow
n+\nu_{e}$, $e^{+}+n\rightarrow p+\bar{\nu}_{e}$); $(ii)$ weak interactions
and neutrino cooling are negligible. In addition to changing the
thermodynamics, weak interactions also change the lepton number and thus the
composition of the disk and its outflows. The composition in stationary state
as parametrized by the electron fraction $Y_{e}=n_{\rm p}/n_{\rm b}$, with
$n_{\rm p}$ and $n_{\rm b}$ denoting the proton and total baryon number
densities, is determined by the degree of electron/positron degeneracy
(Beloborodov, 2003; Chen & Beloborodov, 2007; Siegel & Metzger, 2017; Siegel &
Metzger, 2018a), which we explore further in this paper (Sec. 4.2). This has
important consequences for r-process nucleosynthesis and the nature of the
associated kilonova emission (Sec. 4.3).

Weak interactions are expected to become important above a certain ‘ignition’
threshold on the accretion rate (Metzger et al. 2008b; Chen & Beloborodov
2007; Metzger et al. 2008a; Sec. 2.3),

$\displaystyle\dot{M}_{\rm ign}$ $\displaystyle\approx$
$\displaystyle\dot{\mathcal{M}}_{\rm ign}(M_{\rm BH},\chi_{\rm BH})\alpha_{\rm
vis}^{\frac{5}{3}}$ (1) $\displaystyle\approx$ $\displaystyle 2\times
10^{-3}M_{\odot}\,\text{s}^{-1}\left(\frac{M_{\rm
BH}}{3M_{\odot}}\right)^{\frac{4}{3}}\left(\frac{\alpha_{\rm
vis}}{0.02}\right)^{\frac{5}{3}}.$ (2)

In the second step, we have evaluated the expression for the regime of post-
merger disks, assuming a black-hole of mass $M_{\rm BH}=3M_{\odot}$ and
dimensionless spin of $\chi_{\rm BH}\approx 0.8$ (see Sec. 3) and normalizing
to a dimensionless Shakura-Sunyaev viscosity coefficient $\alpha_{\rm
vis}=0.02$ (see Sec. 4.1.1). While this relation has been found numerically
for 1D disk solutions in Kerr spacetime (Chen & Beloborodov, 2007), we show
here that the scaling $\dot{M}\propto M_{\rm BH}^{4/3}\alpha_{\rm vis}^{5/3}$
can be obtained analytically (see Sec. 2.3). Essentially, the ignition
threshold can be written as a condition on the accretion rate as a result of
the fact that viscous heating, neutrino cooling, and accretion rate scale with
the midplane density (see Sec. 2.3).

At even higher accretion rates of $\dot{M}\gtrsim
0.1\,M_{\odot}\,\text{s}^{-1}$ the disk is expected to become opaque to
neutrinos and at even higher rates of $\dot{M}\gtrsim
1\,M_{\odot}\,\text{s}^{-1}$, neutrinos start to become trapped (Chen &
Beloborodov, 2007). Over a viscous timescale (Eq. (27)), disks initially in
such a state may evolve into an optically thin state through viscous
spreading. These more extreme initial regimes are not the focus of this paper.

Many previous simulations have been performed in hydrodynamics adopting
$\alpha_{\rm vis}$ as a parameter (Fernández & Metzger, 2013; Fernández et
al., 2015; Just et al., 2015; Fujibayashi et al., 2017, 2020), some of which
include or model general-relativistic effects by using a pseudo-Newtonian
potential. While such $\alpha$-disk models are able to qualitatively capture
the evolution of disk density and angular momentum, the nature of turbulence
(convection) is fundamentally different from self-consistent
magnetohydrodynamic turbulence driven by the magnetorotational instability
(MRI; Hawley & Balbus 1992; Balbus & Hawley 2002). MHD disks self-consistently
set an effective $\alpha_{\rm vis}$ and thus self-consistently control the
relative importance of weak interactions to viscous energy transport333This is
only true in three spatial dimensions, as the anti-dynamo theorem in
axisymmetry (Cowling, 1933) does not allow for a steady turbulent state.;
this, in turn, sets the composition of disk and outflow material and thus
determines the nucleosynthetic r-process yields and kilonova colors such disks
give rise to. Furthermore, while $\alpha$-disks dissipate heat generated by
viscosity locally and predominantly in the disk midplane (proportional to the
gas density), MHD disks dissipate a significant fraction non-locally via
reconnection in low-density regions of a disk corona (Jiang et al., 2014;
Siegel & Metzger, 2018a). This difference is crucial in launching outflows,
which originate from this ‘hot’ corona with the additional help of free nuclei
recombining into $\alpha$ particles (Siegel & Metzger, 2018a). Indeed, a
comparison between a post-merger GRMHD and $\alpha$-disk simulation shows that
MHD disks are much more effective in evaporating material early on, giving
rise to most of their ejecta mass during the first few hundred milliseconds,
before viscous spreading (as in an $\alpha$-disk) takes over (Fernández et
al., 2019). This led to the preliminary conclusion that MHD disks can eject up
to 30-40% of their initial disk mass (Siegel & Metzger, 2018a; Fernández et
al., 2019). We demonstrate in this paper (Sec. 4.1.4) that the ejecta masses
depend on the significance of weak interactions (cooling) and that higher
relative mass loss occurs below the ignition threshold.

###  2.2 NS post-merger disks: relation to binary parameters

Figure 1: Approximate mapping between binary parameters and accretion disk as
well as ejecta masses for binary neutron star systems (based on fitting
formulae to numerical relativity simulations as considered by Krüger & Foucart
2020). Left: Disk mass as a function of mass and radius of the secondary
(less-massive) neutron star. Middle: Ratio of disk to dynamical ejecta mass as
a function of the binary mass ratio ($Q=M_{1}/M_{2}$) and radius of the
component neutron stars (assuming both stars have roughly the same radius),
with the mass of the secondary fixed to $M_{2}=1.2~{}M_{\odot}$ and assuming
disk ejecta to be 30% of the initial disk mass. Right: Same as the middle
panel, but for $M_{2}=1.4~{}M_{\odot}$. The dot-dashed lines (in white)
delineate the disk masses $M_{\rm disk}=0.03~{}M_{\odot}$,
$0.002~{}M_{\odot}$, $0.0002~{}M_{\odot}$ considered for post-merger
simulations in this work. Dotted lines (in cyan) show the lower and upper
limits of neutron star radii allowed by recent observations (De et al., 2018;
Capano et al., 2020; Miller et al., 2019c; Riley et al., 2019; Landry et al.,
2020). Figure 2: Approximate mapping between binary parameters and accretion
disk as well as ejecta masses for neutron star–black hole binary systems
(based on fitting formulae to numerical relativity simulations as considered
by Foucart et al. 2018 and Krüger & Foucart 2020). Shown are post-merger disk
masses and ratios of disk to dynamical ejecta (assuming disk ejecta to be 30%
of the initial disk masses) as a function of the NS radius and the binary mass
ratio $Q=M_{\rm BH}/M_{\rm NS}$ for different scenarios, including light and
heavier NSs ($M_{\rm NS}=1.2M_{\odot},\,1.4M_{\odot}$) and slowly or rapidly
spinning BHs (dimensionless spin projected on the binary orbital angular
momentum, $\chi_{1z}$, 0.4 or 0.8). The dot-dashed lines (in white) indicate
the disk masses $M_{\rm disk}=0.03~{}M_{\odot}$, $0.002~{}M_{\odot}$,
$0.0002~{}M_{\odot}$ considered for the simulations in this work. Dotted lines
(in cyan) show the lower and upper limits of neutron star radii allowed by
observations as in Fig. 1.

Figures 1 and 2 show mappings between binary parameters, accretion disk
masses, and ejecta masses (dynamical and disk ejecta) for plausible binary
neutron star (BNS) and neutron star black hole (NSBH) systems. Mappings from
binary parameters to BNS disk masses, BNS dynamical ejecta, and NSBH dynamical
ejecta are based on fitting formulae to numerical relativity simulations as
considered by Krüger & Foucart (2020) for the respective cases. Mappings from
binary parameters to NSBH disk masses are performed by calculating the mass
outside the remnant black hole post merger, $M_{\rm rem}$, using the fitting
formula provided by Foucart et al. (2018), and subtracting the NSBH dynamical
ejecta mass from $M_{\rm rem}$. White lines in Figs. 1 and 2 indicate the disk
models simulated here, highlighting the regime of parameter space we focus on
in this work.

For BNS systems, Krüger & Foucart (2020) find that disk masses extracted from
existing BNS simulations can be fit to $\sim 35\%$ accuracy by a formula of
the type $M_{\rm disk}=M_{2}[aC_{2}+c]^{d}$, which is effectively insensitive
to mass ratio to leading order. Disk masses scale with the mass and
compactness of the secondary (lighter) neutron star
$C_{2}=GM_{2}/(R_{2}c^{2})$, increasing with stiffer EOSs and larger secondary
component masses. For small total mass BNS systems (middle panel of Fig. 1),
mergers give rise to both dynamical ejecta and disk ejecta, irrespective of
the stiffness of the EOS, while for high total mass systems (right panel of
Fig. 1) both disk masses and the amount of dynamical ejecta are reduced or
even non-existent due to the fact that the BNS quickly collapse to a black
hole after merger and little to no material is left outside the event horizon.

Both the binary mass ratio and the EOS determine the dominant source of ejecta
in BNS mergers. Disk ejecta dominates across most of the parameter space for
small total-mass systems (middle panel of Fig. 1), except for very small NS
radii $\lesssim\\!10.5$ km and large mass ratios $Q\gtrsim 1.6$, while
dynamical ejecta dominates in high total-mass systems (right panel of Fig. 1),
except for larger NS radii $\gtrsim\\!12$ km and small-to-medium mass ratios
$Q\lesssim 1.4$. For low-mass BNS systems, the disks simulated here cover the
parameter space for small NS radii $\lesssim\\!11$ km, while for high-mass
systems they span a wide range of $11\lesssim R_{\rm NS}\lesssim 12.5$ km. We
note that the results of our highest-$M_{\rm disk}$ run can be qualitatively
extrapolated to larger disk masses and thus cover most of the remaining
parameter space; however, for very massive disks effects of self-irradiation
of the outflows become increasingly significant, leading to a more pronounced
tail of ‘blue’ ejecta (Miller et al., 2019b). Depending on the secondary mass,
the disk outflows of our simulated models may or may not dominate the total
ejecta of the system.

In NSBH systems, dynamical ejecta and post-merger accretion disks only arise
if the NS is tidally disrupted by the BH in the binary. This disruption
process requires the tidal disruption radius to reside outside of the
innermost stable circular orbit of the BH, thus depending on the spin of the
BH projected on the binary orbital angular momentum, $\chi_{1z}$ and its mass
(and thus on the mass ratio $Q=M_{\rm BH}/M_{\rm NS}$ of the binary). Outflows
from accretion disks typically dominate across most of the parameter space
(cf. Fig. 2), except for light NSs ($\sim\\!1.2M_{\odot}$) with large NS radii
$\gtrsim 11.5$ km and medium-to-large mass ratios $Q\gtrsim 3-4$, somewhat
dependent on the BH spin.

In the NSBH parameter space, our models simulated here reside in the unequal
mass ratio, low BH-spin regimes. In the regime spanned by our low-mass disk
model, disk ejecta is dominant and dynamical ejecta is absent. In the regime
spanned by the highest-mass disk model in this work, dynamical ejecta
dominates only for light NSs ($\sim\\!1.2M_{\odot}$) and large NS radii
$\gtrsim 11.5$ km.

###  2.3 Ignition threshold and simplified 1D disk model

In this section, the basic scaling $\dot{M}_{\rm ign}\propto\alpha_{\rm
vis}^{5/3}$ (Eq. (1)) of the accretion rate at the ignition threshold for weak
interactions in an accretion disk with the dimensionless, effective Shakura-
Sunyaev viscosity coefficient $\alpha_{\rm vis}$ is derived analytically. In
doing so, we obtain a simplified one-dimensional accretion disk model that
proves useful for further simulation analysis in Sec. 4.

We work within a height-integrated, one-dimensional model for accretion disks
in Kerr spacetime with metric $g_{\mu\nu}$ using Boyer-Lindquist coordinates
$x^{i}=(t,r,\theta,\phi)$ (e.g., Beloborodov 1998; Gammie & Popham 1998; Chen
& Beloborodov 2007). In such 1D models, angular averages are performed by
approximating $\int{\rm d}\theta{\rm d}\phi\sqrt{g}q\simeq 4\pi
Hq(\theta=\pi/2)$, where $q$ represents any physical quantity, $H(r)$ is the
characteristic angular half-thickness of the disk, and $g={\rm
det}(g_{\mu\nu})$. We work in the “thin-disk” approximation (Shakura &
Sunyaev, 1973; Bardeen et al., 1972; Novikov & Thorne, 1973), well justified
if neutrino cooling is significant, in which terms in $(H/r)^{2}$ are
neglected and gas pressure is negligible. This results in the assumption that
the fluid orbits with four-velocity $u^{\mu}=(u^{t},u^{r},0,u^{\phi})$ and
Keplerian angular velocity $\Omega\equiv u^{\phi}/u^{t}=c/(r\sqrt{2r/r_{\rm
g}}+\chi_{\rm BH}r_{\rm g}/2)$. Here, $r_{\rm g}=2GM_{\rm BH}/c^{2}$ is the
gravitational radius, and $M_{\rm BH}$ and $\chi_{\rm BH}$ denote, as before,
the mass and the dimensionless spin of the black hole, respectively. In this
thin-disk limit, the equation of vertical hydrodynamical equilibrium for the
disk can be written as (cf. Abramowicz et al. 1997; Beloborodov 1998)

$\left(\frac{H}{r}\right)^{2}=\frac{2r}{c^{2}r_{\rm g}J(\chi_{\rm
BH},r)}\frac{p}{\rho},\mskip 40.0mu$ (3)

or

$c_{\rm s}=c\left(\frac{J(\chi_{\rm BH},r)r_{\rm
g}}{2r}\right)^{\frac{1}{2}}\frac{H}{r},$ (4)

where $c_{\rm s}=\sqrt{p/\rho}$ is the isothermal sound speed, and

$J(\chi_{\rm BH},r)\equiv\frac{2\left(r^{2}-\chi_{\rm BH}r_{\rm
g}\sqrt{2r_{\rm g}r}+3\chi_{\rm BH}^{2}r_{\rm g}^{2}/4\right)}{2r^{2}-3r_{\rm
g}r+\chi_{\rm BH}r_{\rm g}\sqrt{2r_{\rm g}r}}.$ (5)

The equation of baryon number conservation reads

$\dot{M}=-2\pi rcu^{r}\Sigma,$ (6)

where $\Sigma=2H\rho$ is the disk surface density. From the equations of
energy and angular momentum conservation, one can derive the identity (Page &
Thorne, 1974; Chen & Beloborodov, 2007)

$2\nu\Sigma r\sigma^{r}_{\phi}=-\frac{cr_{\rm g}}{4\pi}\dot{M}F(x,\chi_{\rm
BH}),$ (7)

where $\nu$ is the kinematic viscosity,

$\sigma^{r}_{\phi}=(1/2)cg^{rr}g_{\phi\phi}\sqrt{-g^{tt}}\gamma^{3}({\rm
d}\Omega/{\rm d}r)$ (8)

is the shear, with $\gamma=u^{t}/\sqrt{-g^{tt}}$ being the Lorentz factor of
the fluid measured by a zero-angular momentum observer, and $F(x,\chi_{\rm
BH})$ is a function of dimensionless radius $x\equiv\sqrt{2r/r_{\rm g}}$ and
black-hole spin given in Appendix A for completeness.

Substituting Eq. (7) into Eq. (6) and using Eq. (4) one finds

$u^{r}=\frac{8}{3}\left(\frac{J(r,\chi_{\rm BH})r}{2c^{2}r_{\rm
g}}\right)^{\frac{1}{2}}\frac{\sigma^{r}_{\phi}}{F(x,\chi_{\rm
BH})}\alpha_{\rm vis}\left(\frac{H}{r}\right)^{2},$ (9)

where we have adopted the Shakura-Sunyaev parametrization
$\nu=\frac{2}{3}\alpha_{\rm vis}c_{\rm s}H$ with viscosity coefficient
$\alpha_{\rm vis}$. Substituting this back into Eq. (6), we obtain

$\displaystyle\dot{M}$ $\displaystyle=$ $\displaystyle-\frac{32\pi r_{\rm
g}^{2}}{3\sqrt{2}}J^{\frac{1}{2}}(r,\chi_{\rm
BH})\mskip-5.0mu\left(\frac{r}{r_{\rm
g}}\right)^{\frac{5}{2}}\mskip-15.0mu\frac{\sigma^{r}_{\phi}}{F(x,\chi_{\rm
BH})}\alpha_{\rm vis}\rho\mskip-5.0mu\left(\frac{H}{r}\right)^{3}$ (10)
$\displaystyle=$ $\displaystyle 2\sqrt{2}\pi cr_{\rm
g}^{2}S^{-1}J^{\frac{1}{2}}\left(\frac{r}{r_{\rm
g}}\right)^{\frac{3}{2}}\alpha_{\rm vis}\rho\left(\frac{H}{r}\right)^{3},$
(11)

where we have introduced the function $S(r,\chi_{\rm BH})\equiv-(3/8)c(r_{\rm
g}/r)[F(x,\chi_{\rm BH})/\sigma^{r}_{\phi}]$, which varies between zero at the
marginally stable orbit and unity at $r\rightarrow\infty$.

The viscous heating rate per unit area of the disk is given by
$Q^{+}=2\nu\Sigma hc\sigma^{r}_{\phi}\sqrt{-g^{tt}}\gamma({\rm d}\Omega/{\rm
d}r)$, where $h\approx 1$ is the specific enthalpy of the fluid in the thin
disk (Beloborodov, 1998). Using the identity (7) one can rewrite this as

$Q^{+}=\frac{c^{2}}{4\pi}\left(\frac{r}{r_{\rm g}}\right)^{-1}F(x,\chi_{\rm
BH})\mathcal{Q}(r,M_{\rm BH},\chi_{\rm BH})\dot{M},$ (12)

where $\mathcal{Q}(r,M_{\rm BH},\chi_{\rm BH})\equiv-u^{t}({\rm d}\Omega/{\rm
d}r)$.

We assume that cooling of the disk is dominated by electron and positron
capture (URCA cooling). Ignoring final state blocking in the neutrino phase
space, the cooling rate per unit area of the disk is then approximately given
by (Tubbs & Schramm, 1975; Bruenn, 1985; Qian & Woosley, 1996; Popham et al.,
1999)

$Q^{-}_{\nu}=2H\mathcal{C}_{\nu}\rho T^{6}.$ (13)

Here, $\mathcal{C}_{\nu}$ is a constant times the mass fraction of nucleons
$X_{\rm nuc}$, which is roughly unity in the inner parts of the accretion disk
where photodisintegration breaks down nuclei into neutrons and protons once
$T\sim 10^{10}$ K. At sufficiently small $\dot{M}$ (low midplane density), the
pressure is dominated by radiation pressure, $p=\frac{11}{12}a_{\rm SB}T^{4}$,
where contributions of relativistic electron-positron pairs have been included
(Popham et al., 1999). Substituting into Eq. (4), this yields the disk
midplane temperature

$T=\left(\frac{6}{11}\frac{c^{2}}{a_{\rm
SB}}\right)^{\frac{1}{4}}J^{\frac{1}{4}}(r,\chi_{\rm BH})\left(\frac{r}{r_{\rm
g}}\right)^{-\frac{1}{4}}\left(\frac{H}{r}\right)^{\frac{1}{2}}\rho^{\frac{1}{4}}.$
(14)

Using this relation in Eq. (13) together with Eq. (11), one obtains the
following expression for neutrino cooling:

$\displaystyle Q^{-}_{\nu}$ $\displaystyle=$ $\displaystyle
2\mathcal{C}_{\nu}\left(\frac{1}{2\sqrt{2}\pi
c}\right)^{\frac{5}{2}}\left(\frac{6}{11}\frac{c^{2}}{a_{\rm
SB}}\right)^{\frac{3}{2}}$ (16) $\displaystyle\times r_{\rm
g}^{-4}S^{\frac{5}{2}}J^{\frac{1}{4}}\left(\frac{r}{r_{\rm
g}}\right)^{-\frac{17}{4}}\left(\frac{H}{r}\right)^{-\frac{7}{2}}\alpha_{\rm
vis}^{-\frac{5}{2}}\dot{M}^{\frac{5}{2}}.$

Adopting the condition $Q^{-}_{\nu}/Q^{+}=1/2$ for weak interactions to become
energetically significant, one can employ the expressions (12) and (16) to
formulate this as a condition on the accretion rate:

$\displaystyle\dot{M}_{\rm ign}$ $\displaystyle=$
$\displaystyle\frac{11}{12}\sqrt{2}^{\frac{5}{3}}\pi a_{\rm
SB}\mathcal{C}_{\nu}^{-\frac{2}{3}}$ (18) $\displaystyle\times r_{\rm
g}^{\frac{8}{3}}S^{-\frac{5}{3}}J^{-\frac{1}{6}}F^{\frac{2}{3}}\mathcal{Q}^{\frac{2}{3}}\mskip-5.0mu\left(\frac{r}{r_{\rm
g}}\right)^{\frac{13}{6}}\mskip-5.0mu\left(\frac{H}{r}\right)^{\frac{7}{3}}\mskip-5.0mu\alpha_{\rm
vis}^{\frac{5}{3}}$ $\displaystyle\equiv$ $\displaystyle\dot{\mathcal{M}}_{\rm
ign}(r,M_{\rm BH},\chi_{\rm BH})\alpha_{\rm vis}^{\frac{5}{3}}.$ (19)

where $\dot{\mathcal{M}}_{\rm ign}(r,M_{\rm BH},\chi_{\rm BH})$ is a function
that depends on the black-hole parameters. Apart from $S(r,\chi_{\rm BH})$,
which must be calculated numerically (but may be approximated analytically),
$\dot{\mathcal{M}}_{\rm ign}(r,M_{\rm BH},\chi_{\rm BH})$ can be analytically
evaluated on a horizon-scale $r\sim r_{\rm g}$ to provide the characteristic
accretion rate onto the black hole. Noting that at $r\sim r_{\rm g}$,
$\mathcal{Q}\propto r_{\rm g}^{-2}$, we find that the critical accretion rate
onto the black hole scales as

$\dot{M}_{\rm ign}\propto M_{\rm BH}^{\frac{4}{3}}\alpha_{\rm
vis}^{\frac{5}{3}},$ (20)

with a prefactor of order unity that depends on the black hole spin $\chi_{\rm
BH}$.

##  3 Numerical Methods

###  3.1 Simulation setup

We perform simulations in ideal GRMHD and full 3D with a fixed background
spacetime for computational efficiency using the code and numerical setup
described in Siegel & Metzger (2018a). The code is based on GRHydro (Mösta et
al., 2014) and makes use of the Einstein Toolkit444http://einsteintoolkit.org
(Babiuc-Hamilton et al., 2019; Löffler et al., 2012; Schnetter et al., 2004;
Goodale et al., 2003; Thornburg, 2004), with neutrino interactions implemented
via a leakage scheme based on Bruenn (1985) and Ruffert et al. (1996), and
follows the implementation of Galeazzi et al. (2013) and Radice et al. (2016).
Thermodynamic properties of matter are based on the Helmholtz EOS (Timmes &
Arnett, 1999; Timmes & Swesty, 2000), and we compute abundances of nuclei at a
given density, temperature, and electron fraction $Y_{e}$, assuming nuclear
statistical equilibrium.

The simulations include a Kerr black hole of mass $3.0M_{\odot}$ and
dimensionless spin $\chi_{\rm BH}=0.8$, initially surrounded by a torus of
constant specific angular momentum, small constant specific entropy of
$8~{}k_{\rm B}$ per baryon, and initial electron fraction $Y_{\rm e}=0.1$; we
refer to Tab. 1 for a summary of initial torus properties of the simulation
runs. Figure 3 shows disk densities for the initial configurations. Run MD_M03
has been discussed before (Siegel & Metzger, 2017; Siegel & Metzger, 2018a),
and is further elaborated on here by comparing it to the two new runs MD_M002
and MD_M0002, which represent lighter accretion disks. The BH-disk problem is
formulated in Cartesian, horizon-penetrating Kerr-Schild coordinates. The
black-hole mass and spin reflect typical NS merger scenarios (see Sec. 2.2).
The black hole spins in the case of prompt black hole formation from BNS
mergers are typically not larger than $\chi_{\rm BH}\approx 0.8$ (Kiuchi et
al., 2009; Rezzolla et al., 2010; Bernuzzi et al., 2014; Kastaun et al.,
2013), and black hole spins in case of delayed black hole formation are
$\chi_{\rm BH}\lesssim 0.7$ (Sekiguchi et al., 2016). Furthermore, $\chi_{\rm
BH}\sim 0.8$ is significant enough to disrupt the NS and lead to a post-merger
accretion disk in NSBH mergers across a wide range in mass ratio (Foucart
2012; see the discussion in Sec. 2.2). Our initial tori masses are chosen to
reflect a mass range covering the ignition threshold for weak interactions in
post-merger disks (cf. Sec. 2.1) and is typical both for BNS and NSBH
scenarios (Sec. 2.2).

Table 1: Initial configurations of the accretion disks before relaxation. From left to right: black-hole mass and dimensionless spin, disk mass, inner and outer radius of the disk, radius at maximum density, specific entropy, electron fraction, and maximum magnetic-to-fluid pressure ratio. Run | $M_{\rm BH}$ | $\chi_{\rm BH}$ | $M_{\rm d,0}$ | $R_{\rm in,0}$ | $R_{\rm out,0}$ | $R_{0}$ | $s_{0}$ | $Y_{\rm e,0}$ | $p_{b}/p_{\rm f}$
---|---|---|---|---|---|---|---|---|---
| [$M_{\odot}$] |  | [$M_{\odot}$] | [km] | [km] | [km] | [$k_{B}$/b] |  |
MD_03 | 3 | 0.8 | 0.03 | 17.71 | 106.27 | 30 | 8 | 0.1 | $<5\times 10^{-3}$
MD_002 | 3 | 0.8 | 0.002 | 31 | 88.56 | 45.61 | 8 | 0.1 | $<5\times 10^{-3}$
MD_0002 | 3 | 0.8 | 0.0002 | 53.14 | 88.56 | 66.27 | 8 | 0.1 | $<5\times 10^{-3}$
Figure 3: Initial disks in the three simulations performed in this work. The
top and middle panels show density slices along the equatorial and meridional
plane, respectively, at $t=0$ ms. The bottom panels show radial profiles of
density.

The tori are initialized with weak poloidal magnetic seed fields, confined to
the interior of the tori and defined by the magnetic vector potential
$A^{r}=A^{\theta}=0$ and $A^{\phi}=A_{b}$ max$\\{p-p_{\rm cut},0\\}$. Here,
$p$ denotes the fluid pressure, $p_{\rm cut}$ is the pressure below which the
magnetic field is set to zero, and $A_{b}$ sets the initial field strength.
Here, $p_{\rm cut}\approx 1\times 10^{-2}p_{\rm max}$ in all cases, where
$p_{\rm max}$ is the pressure at maximum density in the torus. We choose
$p_{\rm cut}$ such that the magnetic field covers the bulk volume of the
torus, while preventing it from becoming buoyant in the outermost layers and
violently breaking out of the torus at the start of the simulation. We adjust
$A_{b}$ such that the magnetic-to-fluid pressure ratio in the torus is a small
value, $p_{\rm B}/p_{\rm f}<5\times 10^{-3}$. This ratio provides relatively
‘weak’ (dynamically unimportant) initial magnetic field strengths with maximum
values of $\approx\\!3\times 10^{14}$ G for MD_M03, $\approx\\!6.4\times
10^{13}$ G for MD_M002, and $\approx\\!1.3\times 10^{13}$ G for MD_M0002.

The initial torus is embedded in a tenuous atmosphere with $T=10^{5}$ K,
$Y_{\rm e}=1$, and $\rho\approx$ 37 g cm-3, $\rho\approx$ 3.7 g cm-3, and
$\rho\approx$ 0.37 g cm-3 for runs MD_M03, MD_M002, and MD_M0002,
respectively. The density and temperature are chosen such that they are
sufficiently low to neither impact the dynamics nor the composition of the
disk outflows. The atmosphere densities are set to approximately scale with
the maximum density of the accretion disk during the evolution (cf. Tab. 2).
The total atmosphere mass of the entire computational domain is $3.8\times
10^{-5}M_{\odot}$ for MD_M03, $3\times 10^{-7}M_{\odot}$ for MD_M002, and
$3\times 10^{-8}$ for MD_M0002, orders of magnitude smaller than the disk
ejecta (cf. Tab. 2); over a volume of radius 1000 km, which we consider as the
minimum radius for outflow material to be unbound from the BH-disk system, the
corresponding atmosphere masses are $1.8\times 10^{-8}M_{\odot}$ for MD_M03,
$1.8\times 10^{-9}M_{\odot}$ for MD_M002, and $1.8\times 10^{-8}M_{\odot}$ for
MD_M0002. At the chosen atmosphere temperature of $T=10^{5}$ K weak
interations are frozen out.

The computational domain represents a Cartesian grid hierarchy centered around
the black hole. For MD_M03, the grid has eight refinement levels, with an
extent in each coordinate direction of $1.53\times 10^{4}$ km. For MD_M002 and
MD_M0002, the grid has seven refinement levels, with an extent in each
coordinate direction of $1.14\times 10^{4}$ km. The initial tori have
diameters of 240 km, 206 km, and 206 km for simulations MD_M03, MD_M002, and
MD_M0002, respectively. The initial tori are encompassed by the finest
refinement level of the corresponding grid hierarchy. Following previous work
(Siegel et al., 2013; Siegel & Metzger, 2018a; Siegel et al., 2019), the
finest resolution is $\Delta_{xyz}\approx 850$ m for all simulations, chosen
such that the MRI is well resolved in the stationary turbulent state of the
disk (typically by at least ten grid points per fastest-growing MRI mode),
which ensures convergence of global observables (see, e.g., Siegel et al.
2019).

Figure 4: Slices of rest-mass density and optical depth for $\bar{\nu}_{e}$
neutrinos in the meridional plane at 30 ms shortly after the disks have
relaxed into a quasi-stationary state. This state represents the actual
initial data for our simulations. All disks are optically thin to neutrinos.

Angular momentum transport is mediated by magnetic turbulence, driven by the
MRI, in our setups. The initial tori undergo a relaxation phase with self-
consistent magnetic-field amplification, and settle into a quasi-stationary
phase at $\sim\\!20$ ms. We consider the relaxed state of the disks at this
time as the actual initial data for our simulations, and exclude the early
transient phase from our analysis—in particular, all material accreted onto
the black hole or ejected via outflows during this phase. Figure 4 shows
snapshots of the disks just after having settled into the quasi-stationary
phase, indicating that all configurations are optically thin to neutrinos.

###  3.2 Diagnostics

In order to monitor certain physical quantities in the disk, we compute radial
profiles of a quantity $\chi(\varpi)$ by performing azimuthal, density-
weighted averages, integrating up to one scale height of the disk:

$\langle\chi\rangle_{{\rm
azim},z_{H}}=\frac{\int^{z_{H}}_{-z_{H}}\int_{0}^{2\pi}\chi\hat{D}\varpi d\phi
dz}{\int^{z_{H}}_{-z_{H}}\int_{0}^{2\pi}\hat{D}\varpi d\phi dz}.$ (21)

Here, $\varpi=\sqrt{x^{2}+y^{2}}$ is the cylindrical radius,
$\hat{D}=\sqrt{\gamma}\rho W$ is the conserved rest-mass density as seen by
the Eulerian observer moving normal to the spatial hypersurfaces in the 3+1
split of Kerr-Schild spacetime, $\gamma$ is the determinant of the spatial
metric $\gamma_{ij}$ in 3+1 split, and $W$ the Lorentz factor of the fluid.
The density scale height $z_{H}$ is defined as

$z_{H}(\varpi)=\frac{\int\int_{0}^{2\pi}|z|\hat{D}\varpi d\phi
dz}{\int\int_{0}^{2\pi}\hat{D}\varpi d\phi dz}.$ (22)

By integrating up to the local density scale height, we exclude the disk
corona and winds from the calculation. For some quantities, temporal averages
$\langle\cdot\rangle_{t}$ (cf., e.g., Eqs. (26) and (25)) over a specified
time window are taken, in order to reduce the effect of temporal fluctuations
of a turbulent medium.

For some disk quantities, it is useful to compute the rest-mass density
average and study its evolution over time. We define this averaging as

$\langle\chi\rangle_{\hat{D}}=\frac{\int\chi\hat{D}d^{3}x}{\int\hat{D}d^{3}x}.$
(23)

In some cases, this average is calculated by integrating only up to the local
scale height,

$\langle\chi\rangle_{\hat{D},z_{H}}=\frac{\int^{z_{H}}_{-z_{H}}\chi\hat{D}\varpi
dz}{\int^{z_{H}}_{-z_{H}}\hat{D}\varpi dz},$ (24)

which allows us to explicitly exclude the disk corona and disk wind regions.

##  4 Numerical Results

We proceed by first discussing our results for global disk properties, such as
MHD mediated angular momentum transport (Sec. 4.1.1), accretion (Sec. 4.1.2),
neutrino emission (Sec. 4.1.3), disk ejecta (Sec. 4.1.4), before discussing
disk evolution locally in terms of compositional changes and weak interactions
(Sec. 4.2), and nucleosynthesis from the disk outflows (Sec. 4.3).

###  4.1 Global properties

Table 2: Various properties of the accretion disks simulated here. From left to right: i) average density of the inner accretion disk, ii) disk mass after relaxation, iii) accretion rate after relaxation, iv) estimated total unbound disk outflow mass ($M_{\rm ej}$; ‘ejecta’) in units of the effective initial disk mass after relaxation, v) with mean electron fraction, vi) maximum total neutrino luminosity in electron and anti-electron neutrinos, vii) effective $\alpha$-viscosity parameter of the inner accretion disk, viii) viscous timescale and ix) total simulated time for the respective runs. The averages $\bar{\rho}_{\rm d}$ and $\bar{\alpha}_{\mathrm{vis}}$ with the quoted uncertainties represent density-weighted average, maximum, and minimum values around a 10 ms window centered at $t=30$ ms (see the text for details). The values for $\dot{M}_{d}$ with quoted uncertainties represent the average, maximum, and minimum values over the same time window at $t=30$ ms. The values for $L_{\nu,\mathrm{max}}$ represent the average value extracted from the neutrino luminosity over a 10 ms window centered at the global peak luminosity of anti-electron neutrinos.  Run | $\bar{\rho}_{\rm d}$ | $M_{\rm d,20}$ | $\dot{M}_{d}$ | $M_{\rm ej}$ | $\bar{Y}_{e}$ | $L_{\nu,\mathrm{max}}$ | $\bar{\alpha}_{\mathrm{vis}}$ | $t_{\rm vis}$ | $t_{\rm sim}$
---|---|---|---|---|---|---|---|---|---
| [$\mathrm{g}\,\mathrm{cm}^{-3}$] | [$10^{-2}M_{\odot}$] | [$M_{\odot}\,\mathrm{s}^{-1}$] | [$M_{\rm d,20}$] |  | [$\mathrm{erg}\,\mathrm{s}^{-1}$] | $\times 10^{-2}$ | [ms] | [ms]
MD_03 | $5.53_{-1.4}^{+2.2}\times 10^{10}$ | 1.89 | $2.6_{-0.7}^{+1.4}\times 10^{-1}$ | 0.17 | 0.174 | 2.1$\times 10^{52}$ | $1.76_{-1.7}^{+2.8}$ | 279 | 381
MD_002 | $6.6^{+0.5}_{-0.9}\times 10^{9}$ | 0.184 | $1.14_{-0.4}^{+0.4}\times 10^{-2}$ | 0.22 | 0.114 | 14$\times 10^{50}$ | $1.01_{-0.6}^{+1.5}$ | 874 | 309
MD_0002 | $3.4^{+0.7}_{-0.6}\times 10^{8}$ | 0.0199 | $3.79_{-1.7}^{+2.5}\times 10^{-4}$ | 0.30 | 0.101 | 2.1$\times 10^{48}$ | $1.14_{-1.1}^{+2.0}$ | 1731 | 294

####  4.1.1 MHD turbulence & effective viscosity

Soon after the start of the simulations, our accretion disks show vigorous
magnetic turbulence, triggered by the MRI, a local fluid instability developed
in differentially rotating magnetized fluids (Velikhov, 1959; Chandrasekhar,
1960; Balbus & Hawley, 1991, 1998; Balbus, 2003; Armitage, 2011). The initial
weak magnetic field in the simulations, which is dynamically not important and
purely poloidal, is amplified by the MRI at an exponential rate. A toroidal
field component is initially generated and amplified by magnetic winding,
before the toroidal field becomes susceptible to the MRI as well. This
combination of the MRI and magnetic winding causes an overall increase of the
maximum total magnetic field strength by one to two orders of magnitude from
its initial value to saturation. At $t\approx 20$ ms, a steady turbulent state
is achieved by the disk and the magnetic field, self-consistently amplified by
the MRI, reaches saturation, losing memory of the initial magnetic field
configuration. We refer to Siegel & Metzger (2018a) for a more detailed
discussion on how this turbulent state arises. The resolution employed here
(see Sec. 3.1) guarantees a converged saturation level for the magnetic field
strength of each run; thus, the efficiency of MHD mediated angular momentum
transport does not depend on the particular resolution employed here (Siegel
et al., 2019). Figure 5 shows the number of grid points per fastest growing
MRI wavelength, $\lambda_{\mathrm{MRI}}/\Delta{x}$, in the three accretion
disk models, shortly after the disks have reached a quasi-stationary state.
The MRI can be considered to be well resolved when
$\lambda_{\mathrm{MRI}}/\Delta{x}\gtrsim 10$. Figure 5 shows that at $t=30$
ms, $\lambda_{\mathrm{MRI}}/\Delta{x}$ is at least $\sim\\!10$ in all the
models, and it continues to be so through the rest of the simulations. This
condition can be violated locally in space and time due to the turbulent
nature of the accretion flow.

Figure 5: Meridional slices through the simulated disks, showing the number of
grid points per fastest growing MRI wavelength at $t=30$ ms. Overlaid are
contours of density at $\rho=10^{7},10^{8},10^{9},$ and $10^{10}$ g cm-3.

MHD turbulence in the disk operates as large-scale viscosity, which can be
parametrized by an effective Shakura-Sunyaev viscosity $\alpha_{\rm vis}$. We
define this quantity as the ratio of total stress to fluid pressure,

$\alpha_{\rm vis}(\varpi)=\frac{\langle\langle|T^{r,\phi}|\rangle_{{\rm
azim},z_{H}}\rangle_{t}}{\langle\langle p\rangle_{{\rm
azim},z_{H}}\rangle_{t}},$ (25)

where $T^{r,\phi}$ is the $r-\phi$ component of the stress-energy tensor in
the frame comoving with the fluid, and $p$ the fluid pressure. Here, time
averages are taken over a few neighboring data snapshots.

Figure 6: Radial profiles of the effective Shakura-Sunyaev viscosity parameter
$\alpha_{\rm vis}$ defined in Eq. (25) at different times during the evolution
for the three simulation runs in this work.

Figure 6 shows radial profiles of $\alpha_{\rm vis}$, computed using Eq. (21),
at different times during the disks’ evolution, for all three simulation runs.
We find that $\alpha_{\rm vis}$ is roughly constant among the simulations,
independent of the initial disk mass. Table 2 reports global values for
$\alpha_{\rm vis}$ for each simulation. The latter values describe the
accretion state of the inner disk. We compute these values as an average over
a 10 ms time-window centered at $t=30$ ms, extracted from the time evolution
of the absolute value of the rest-mass density average of $\alpha_{\rm vis}$,

$\alpha_{\rm
vis}=\frac{\langle\langle|T^{r,\phi}|\rangle_{\hat{D},z_{H}}\rangle_{t}}{\langle\langle
p\rangle_{\hat{D},z_{H}}\rangle_{t}}.$ (26)

The rest-mass density average is calculated following the procedure in
Equation (24), restricting to the region between 30-175 km from the centers of
the black holes. By applying this restriction along the radial extent, only
the inner disk regions that can be directly related to the accretion rate onto
the black hole are used in the calculation; the data of the outer parts of the
disks that start to viscously spread are excluded from the calculation. Also
excluded are the innermost regions that show a rise of
$\alpha_{\mathrm{visc}}(\varpi)$ toward the innermost stable circular orbit
and the black-hole horizon due to increasing mean magnetic field strengths
(Penna et al. 2013; cf. Fig. 6). As indicated by the radial profiles in Fig.
6, the averages for $\alpha_{\rm vis}$ in the inner disks as reported in Tab.
2 are roughly constant among the different accretion disks explored here.

The extraction of effective $\alpha$-viscosities allows us to compute
approximate viscous evolution timescales for our disks,

$\displaystyle t_{\rm vis}$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha}\left(\frac{R_{0}^{3}}{GM_{\rm
BH}}\right)^{1/2}\left(\frac{z_{H}}{\varpi}\right)^{-2}$ (27)
$\displaystyle\approx$ $\displaystyle\mskip-15.0mu650\,{\rm
ms}\left(\frac{\alpha}{0.01}\right)^{-1}\mskip-5.0mu\left(\frac{R_{0}}{30\,{\rm
km}}\right)^{\frac{3}{2}}\mskip-5.0mu\left(\frac{M_{\rm
BH}}{3M_{\odot}}\right)^{-\frac{1}{2}}\mskip-5.0mu\left(\frac{z_{H}/\varpi}{0.2}\right)^{-2}\mskip-5.0mu,$

where we have normalized to typical values extracted from our simulations in
the second line. This estimate is remarkably similar to the duration of most
short gamma-ray bursts, the bulk of which have typical durations of a few
hundred milliseconds (as measured by their $T_{90}$, the time interval in
which 5% to 95% of the total fluence is received by the detector; e.g., Berger
2014; Kumar & Zhang 2015). The approximate viscous timescales for our
simulation runs are listed in Tab. 2; they span the range of short gamma-ray
burst $T_{90}$ durations up to $T_{90}\lesssim 2\,{\rm s}$.

####  4.1.2 Accretion

Figure 7: Top: Black hole accretion rates as a function of time for the three
simulation runs in this work. Bottom: Density-averaged cylindrical radius of
matter as a function of time, indicating viscous spreading over time.  Figure
8: Total accreted mass onto the black hole as a function of time for the three
simulation runs in this work. As a result of viscous spreading, the total
cumulative accreted mass essentially converges over the duration of the
simulation runs, and all remaining material will be eventually unbound from
the system.

Magnetic turbulence in the disks generated by the MRI drives accretion onto
the black hole at the center, and outward transport of angular momentum in the
disks. Figure 7 shows the evolution of accretion rates $\dot{M}$ with time for
all three simulations. The accretion rates for models MD_M03, MD_M002, and
MD_M0002, time-averaged over a 10 ms window around $t=30$ ms are
$2.6_{-0.7}^{+1.4}\times 10^{-1}~{}M_{\odot}s^{-1}$, $1.14_{-0.4}^{+0.4}\times
10^{-2}~{}M_{\odot}s^{-1}$, and $3.79_{-1.7}^{+2.5}\times
10^{-4}~{}M_{\odot}s^{-1}$ respectively; the accretion rate changes by a
similar order of magnitude as the disk masses among the three simulation runs.
This is expected from one-dimensional disk models, for which we have (see Eq.
(11))

$\dot{M}\propto\alpha_{\rm vis}\bar{\rho}_{\rm d}M_{\rm
BH}^{2}\left(\frac{z_{H}}{\varpi}\right)^{3},$ (28)

where $\bar{\rho}_{\rm d}$ is the disk midplane density. We compute a
representative value for $\bar{\rho}_{\rm d}$ in the inner accretion disk by
calculating a three-dimensional spatial average of the disk following Eq.
(24), time-averaged over a 10 ms window around $t=30$ ms (see Tab. 2). The
average, minimum, and maximum values for the scale height $z_{H}/\varpi$ are
extracted in an analogous way as for $\alpha_{\rm vis}$, restricting to the
inner disk region between 30-175 km in a 10 ms window centered at $t=30$ ms.
We find $0.23_{-0.1}^{+0.1}$ for MD_M03, $0.23_{-0.1}^{+0.2}$ for MD_M002, and
$0.21_{-0.1}^{+0.2}$ for MD_M0002. We note that the quantities, self-
consistently set by MHD turbulence in our simulations, roughly satisfy
relation (28), within the uncertainties in extracting these numbers from the
3D turbulent evolution of the disks.

As a result of accretion onto the black hole, the outer parts of the disk are
forced to viscously spread. The bottom panel of Fig. 7 shows the early
evolution of the density-averaged cylindrical radius of matter
$\langle\varpi\rangle_{{\rm azim},z_{H}}$ in the simulations, integrated up to
the local density scale height. This parameter can serve as a rough indicator
of radial viscous spreading; such spreading is indeed evident from Fig. 7. The
smaller the viscous timescale (cf. Tab. 2) the faster the viscous spreading
proceeds initially.

Figure 8 shows the evolution of the accreted mass by the black hole as
measured by the mass flux through a spherical coordinate detector surface
placed at a radius of 12 km over the duration of the simulations. The total
accreted mass by the black hole is $\approx$ 46%, $\approx$ 60%, and $\approx$
31% of the disk mass at 20 ms for MD_M03, MD_M002, and MD_M0002 respectively.
Most of the accretion completes within the first $\sim$ 50-100 ms, after which
the accretion rate starts to drop rapidly, with a comparatively negligible
amount of mass projected to be accreted onto the black hole past the end of
the simulation. We ascribe this effect to viscous spreading of the disks,
which forces the disk material remaining at the end of the simulations to be
unbound from the system. Properties of ejecta are discussed in Sec. 4.1.4.

####  4.1.3 Neutrino emission

Figure 9: Characteristics of neutrino emission for the disk models explored in
this work. Top: total neutrino luminosity. Bottom: mean neutrino emission
temperature. The simulations end at $t\sim\\!290$ ms after which the
quantities are extrapolated by power law fits to late times. Figure 10:
Radiative efficiency as a function of accretion rate for accretion onto a
black hole with mass $3~{}M_{\odot}$ and dimensionless spin 0.8, as derived
from the simulations in this work. The efficiency rises with increasing
$\dot{M}$, reaches a maximum just above the ignition threshold $\dot{M}_{\rm
ign}\sim 1\times 10^{-3}\,M_{\odot}\,{\rm s}^{-1}$ (cf. Eq. (2)), then
decreases with increasing $\dot{M}$. The black dashed line shows the maximum
radiative efficiency possible in such a system (see the text for details).

Properties of neutrino radiation are reported in Figure 9 and Tab. 2. As in
Siegel & Metzger (2018a), we define the total neutrino luminosity $L_{\nu}$ of
a given species and the corresponding mean neutrino emission temperature
$\bar{T}_{\nu}$ as

$L_{\nu_{i}}=\int\alpha WQ_{\nu_{i}}^{\rm eff}\alpha\sqrt{\gamma}d^{3}x,$ (29)

and

$\bar{T}_{\nu_{i}}=\frac{\int TQ_{\nu_{i}}^{\rm
eff}W\alpha\sqrt{\gamma}d^{3}x}{\int Q_{\nu_{i}}^{\rm
eff}W\alpha\sqrt{\gamma}d^{3}x},$ (30)

respectively. Here, $Q_{\nu_{i}}^{\rm eff}$ is the effective neutrino
emissivity, $W$ is the Lorentz factor, and $\alpha$ is the lapse function.
Neutrino species are labelled by $\nu_{i}\in{\nu_{e},\bar{\nu}_{e},\nu_{x}}$,
where $\nu_{x}$ represents all heavier neutrino species collectively. We fit
power laws to the late time simulation data to extrapolate these quantities
beyond the time range modeled in the simulations.

Figure 9 shows that the electron and anti-electron neutrino luminosities are
at least an order of magnitude larger than the heavier neutrino luminosities.
Emission channels for the heavier species are relatively suppressed at the
comparatively low densities and temperatures of such accretion disks. For
accretion disks above or close to the ignition threshold (see below and Sec.
2.1), luminosities reach their maximum initially when the disks are still
compact and in their high-density and high-temperature regime, before starting
to viscously spread (cf. MD_M03 and MD_M002). For disks below the ignition
threshold, neutrino luminosities may peak at a later time, once energy
released through viscous accretion further heats up the gas and thus increases
the neutrino emission (cf. MD_M0002, right panels of Fig. 9). Table 2 reports
the total $L_{\nu}$-value for $\nu_{e}$ and $\bar{\nu}_{e}$ for each
simulation run, extracted as an average over a $10$ ms time-window around the
peak luminosities of the respective runs. These range between
$\approx\\!1\times 10^{52}$ erg s-1 (MD_M03) and $\approx\\!1\times 10^{48}$
erg s-1 (MD_M0002). At later simulation times, as the accretion disks spread
radially and become less compact, luminosities start to quickly fade over the
timescales of the simulations. This indicates that neutrino self-irradiation
of outflows in the context of r-process nucleosynthesis (see Sec. 4.3) is
likely only important initially, and much less so for the less luminous disks
MD_M002 and MD_M0002.

The qualitatively different behavior of our disk models in terms of weak
interactions is captured by the differences in radiative efficiency
$L_{\nu}/\dot{M}c^{2}$ among the simulations. Figure 10 shows the variation of
the radiative efficiency of the disks as a function of their accretion rate
$\dot{M}$. The ratio $L_{\nu}/\dot{M}c^{2}$ represents the amount of accreted
rest-mass energy that is turned into radiation per unit time. In order to
assess radiative efficiency, for each simulation, we extract $L_{\nu}$ and
$\dot{M}$ as mean values over the time range $t=25-35$ ms. Fig. 10 shows the
resulting efficiencies of $5.61^{+0.9}_{-0.9}\times 10^{-2}$,
$6.73^{+3.7}_{-1.8}\times 10^{-2}$, and $9.01^{+6.0}_{-2.5}\times 10^{-4}$ for
MD_M03, MD_M002, and MD_M0002, respectively, compared to the maximum possible
radiative efficiency. The latter is a fundamental limit on the amount of
energy that can be extracted from a black hole accretion flow, determined by
the available binding energy (Thorne, 1974),

$[L_{\nu}/\dot{M}c^{2}]_{\rm max}=1-E_{\rm ms},$ (31)

where

$E_{\rm ms}=\frac{1-2M_{\rm BH}}{3r_{\rm ms}}$ (32)

is the specific energy at the marginally stable circular orbit of a Kerr black
hole (Bardeen et al., 1972),

$r_{\rm ms}=M_{\rm BH}\\{3+Z_{2}\mp[(3-Z_{1})(3+Z_{1}+2Z_{2})]^{1/2}\\},$ (33)

with

$\displaystyle Z_{1}$ $\displaystyle=$ $\displaystyle 1+(1-\chi_{\rm
BH}^{2}/M_{\rm BH}^{2})^{1/3}[(1+\chi_{\rm BH}/M_{\rm BH})^{1/3}$ (34)
$\displaystyle+(1-\chi_{\rm BH}/M_{\rm BH})^{1/3}],$ $\displaystyle Z_{2}$
$\displaystyle=$ $\displaystyle(3\chi_{\rm BH}^{2}/M_{\rm
BH}^{2}+Z_{1}^{2})^{1/2}.$ (35)

The efficiency in realistic scenarios, as also seen from the simulation data
in Fig. 10, is smaller than the maximum theoretical value, as a fraction of
the binding energy is stored in the disk and radially advected into the black
hole in the form of heat. The amount of heat radiated away thus depends on the
efficiency of radiative processes with respect to radial advection of energy
(cf. Sec. 2.3). At low $\dot{M}$ values, the low midplane densities and
temperatures in the inner disk suppress neutrino emission, resulting in a low
radiative efficiency. As $\dot{M}$ increases, higher midplane densities and
temperatures enhance neutrino emission, causing the radiative efficiency to
rise and to reach a maximum just above the ignition threshold (Secs. 2.1 and
2.3). As $\dot{M}$ increases further, neutrino cooling initially becomes more
effective; eventually, however, the accretion timescale becomes shorter than
the neutrino cooling timescale due to an increase in optical depth: the very
high $\dot{M}$ values lead to high midplane densities (cf. Eq. (28)) and thus
to an increase in the optical depth close to the midplane, eventually trapping
neutrinos and reducing the cooling volume. One therefore expects enhanced
radial advection of energy and a decrease in radiative efficiency with
increasing accretion rate.

This behavior is evident from Fig 10, which shows a stark rise in radiative
efficiency between the runs MD_M0002 and MD_M002 around an ignition threshold
of $\dot{M}_{\rm ign}\sim 1\times 10^{-3}\,M_{\odot}\,{\rm s}^{-1}$ as
predicted by Eq. (2). For even more massive accretion disks, a strong decline
in the radiative efficiency would be expected as argued above (see also Chen &
Beloborodov 2007), but such disks are beyond the scope of the present study.
We note that our results are qualitatively consistent with previous studies
from one-dimensional disk models (Chen & Beloborodov, 2007).

####  4.1.4 Ejecta

Figure 11: Outflow properties of the three accretion disks simulated here.
Top: Mass distribution of unbound and total outflows as recorded by tracer
particles in terms of their cylindrical radius $\varpi_{\rm ej}$ at which they
get ejected from the disk corona. Bottom: Cumulative distribution of the
kinetic energy of unbound outflows (normalized by the corresponding total
kinetic energy) as a function of asymptotic escape velocity $v_{\infty}$ as
estimated by a spherical detector surface placed at 1000 km from the black
hole and integrated up to different times.

Outflows from the disks originate in high specific entropy (‘hot’) disk
coronae as a result of an imbalance between heating and cooling at high
latitudes. Viscous heating as well as dissipation of magnetic energy far from
the midplane is not offset by neutrino cooling, which becomes energetically
subdominant at the comparatively low densities off the midplane. The precise
amount of material unbound from the disk depends on the details of the self-
consistent heating-cooling imbalance generated by MHD turbulence. We refer to
Siegel & Metzger (2018a) for a detailed discussion on the emergence of these
outflows in the presence of MHD turbulence.

We track properties of these outflows using $10^{4}$ tracer particles of equal
mass, placed throughout the disk initially, with probability density
proportional to the conserved rest mass density $\hat{D}=\sqrt{\gamma}\rho W$.
The top panels of Fig. 11 show the mass distribution of total and unbound disk
outflow (ejecta) as recorded by the tracer particles in terms of the
cylindrical radius at which the tracer particles are ejected from the disk
corona. This ejection radius is defined as the cylindrical radius after which
their radial coordinate position only increases with time. We define unbound
material as having reached a coordinate radius of $10^{3}$ km and additionally
having positive specific energy at infinity, $-hu_{0}>0$, where $h$ denotes
the specific enthalpy and $u_{0}$ is the covariant 0-component of the fluid
four-velocity. For all simulations, we find that the total outflow material is
almost evenly distributed over a broad range of radii, whereas the majority of
the ejecta mass originates from the inner accretion disk ($\varpi\lesssim
100-250$ km), where most of the binding energy is released through viscous
heating. The total amount of material ejected by the disks is estimated to be
$3.179\times 10^{-3}~{}M_{\odot}(\approx 17\%$ of the initial disk mass) for
MD_M03, $4.03\times 10^{-4}~{}M_{\odot}(\approx 22\%$ of the initial disk
mass) for MD_M002, and $6.046\times 10^{-5}~{}M_{\odot}(\approx 30\%$ of the
initial disk mass) for MD_M0002 by the end of the corresponding simulations.

The disks still generate steady winds at the end of our simulations, which
indicates that the remaining disks would continue evaporating themselves if
the simulations were evolved over longer timescales. Given that the total
cumulative accreted masses have already converged over the duration of the
simulation runs (cf. Sec. 4.1.2), we predict that the remaining disk will
generate outflow material corresponding to $\approx\\!37\%$ of the initial
disk mass for MD_M03, $\approx\\!18\%$ of the initial disk mass for MD_M002,
and $\approx\\!38\%$ of the initial disk mass for MD_M0002. This material will
be ejected over longer timescales by a combination of the MHD-driven outflows
as discussed here and viscous spreading of the disk over several viscous
timescales (see Fernández et al. 2019 for a discussion of late-time viscous
outflows). Ejecta percentages are uncertain by at least $\pm 10\%$ due to the
ambiguities and approximations in defining and measuring ejected mass.

While we find disks above the ignition threshold (cf. Secs. 2.1 and 4.1.3)
eject roughly 30–50% of their initial disk mass after relaxation, disks below
the ignition threshold even eject more $\gtrsim 60\%$. We ascribe such
increased evaporation to the absence of a significant cooling mechanism for
the accretion flow not just off the midplane, but also in the midplane and
thus throughout the disk. Strongly enhanced heating is evident both in the
midplane (cf. the large midplane entropy for MD_M0002 compared to MD_M03 and
MD_M002 in Fig. 13) as well as in the ‘corona’ (cf. the increase in specific
entropy at high latitudes between MD_M0002 and MD_M002 in Fig. 14).

The bottom panels of Figure 11 show the distribution of kinetic energy of the
disk ejecta in terms of the asymptotic escape speed $v_{\infty}$ as measured
by a spherical detector surface placed at $\approx\\!1000$ km from the black
hole. The asymptotic velocities $v_{\infty}$ are computed from the asymptotic
Lorentz factor $W_{\infty}=-hu_{0}$. In all three simulations, we find the
emergence of a high-velocity tail, with the bulk ejecta residing around
$v_{\infty}\approx 0.1c$. We find kinetic-energy weighted mean escape
velocities of $\bar{v}_{\infty}\approx 0.14c,\,0.21c,\,0.17c$ for MD_M03,
MD_M002, and MD_M0002, respectively. This velocity scale is naturally
explained as a combination of moderate outflow velocities from the hot corona
of typically $(0.03-0.1)c$ together with the energy released by recombination
of individual nucleons into $\alpha$-particles ($\approx\\!7$ MeV per baryon
per $\alpha$-particle formed). Our velocity distributions are qualitatively
similar to Fernández et al. (2019) and Christie et al. (2019a), albeit their
GRMHD run—similar to our MD_M03 case—shows more outflow mass at higher
velocities; this has been commented on in Fernández et al. (2019) and can be
attributed in part to their initial magnetic field configuration, which is
optimized for fast magnetically-dominated outflows. The high-velocity tail for
the low-$\dot{M}$ disks MD_M002 and MD_M0002 are more pronounced, which may be
ascribed to more violent viscous heating as neutrino cooling becomes less
important at low accretion rates (cf. increase in specific entropy in Figs. 13
and 14).

###  4.2 Weak interactions and disk composition

Figure 12: Snapshots of the electron fraction $Y_{e}$ and electron degeneracy
$\eta=\mu_{e}/k_{\rm B}T$ over time for the three simulation runs performed
here. Plotted are slices of these quantities in the disk midplane ($x$-$y$
plane) at times $t=42.8$ ms, $t=130$ ms, and $t=250$ ms. Above the ignition
threshold for weak interactions (Eq. (2)), the high accretion rates of runs
MD_M03 and MD_M002 cause the inner disk to remain strongly neutron-rich over
time, whereas below the ignition threshold the low accretion rate of run
MD_M0002 leads to accelerated protonization in the inner part.  Figure 13:
Evolution of the ratio of anti-electron to electron neutrino number emission
rates $R_{\bar{\nu}_{\rm e}}^{\rm eff}/R_{\nu_{\rm e}}^{\rm eff}$ and specific
entropy $s$ for the three simulation runs performed here. Plotted are
snapshots of these quantities in the disk midplane ($x$-$y$ plane) at times
$t=42.8$ ms, $t=130$ ms, and $t=250$ ms. Figure 14: Space-time diagram of
specific entropy in accretion disks above (top) and below (bottom) the
ignition threshold. For each disk, variation of specific entropy with height
(in the $z$ direction) measured from the disk midplane ($z=0$) is shown as a
function of time; at each epoch, specific entropy values are extracted from
the $x-z$ plane of the disk with a radial averaging performed between 45 km
and 70 km. Figure 15: Snapshots of the electron fraction $Y_{e}$ at $t=250$ ms
in the meridional plane ($x$-$z$ plane) for the three simulation runs
performed here, indicating lateral as well as radial composition gradients and
significant mixing.

Accretion disks below and above the ignition threshold for weak interactions
give rise to qualitatively different evolution of their composition as we
shall illustrate in this section. Discussing weak interactions in the disk in
more detail also benefits the interpretation of the behavior of some global
disk quantities (Sec. 4.1) as well as of nucleosynthesis in disk outflows
(Sec. 4.3).

Figures 12 and 13 show various quantities pertaining to the evolution of disk
composition, including snapshots of the electron fraction, electron degeneracy
$\eta=\mu_{e}/k_{\rm B}T$, where $\mu_{e}$ is the electron chemical potential,
the ratio of $\bar{\nu}_{e}$ to $\nu_{e}$ number emission rates
($R_{\bar{\nu}_{\rm e}}^{\rm eff}/R_{\nu_{\rm e}}^{\rm eff}$), and specific
entropy $s$ over the course of the three simulation runs.

Neglecting absorption of neutrinos, roughly appropriate for the disk midplane
in the regimes probed here, the weak interactions controlling disk composition
are the charged-current $\beta$-processes,

$\displaystyle e^{-}+p$ $\displaystyle\rightarrow$ $\displaystyle n+\nu_{\rm
e},$ (36) $\displaystyle e^{+}+n$ $\displaystyle\rightarrow$ $\displaystyle
p+\bar{\nu}_{\rm e}.$ (37)

Starting off from neutron-rich initial conditions $Y_{e}\approx 0.1$ (cf. Sec.
3), one expects the disks to protonize over time due to the dominant reaction
Eq. (37). This is evident from $R_{\bar{\nu}_{\rm e}}^{\rm eff}/R_{\nu_{\rm
e}}^{\rm eff}>0$ and the gradual increase of $Y_{e}$ over time in the outer
parts of the accretion disks in Figs. 13 and 12. However, less massive disks
with lower $\dot{M}$ such as MD_M002 and MD_M0002 take longer times compared
to MD_M03 to raise their electron fraction from its initial value; this can be
understood from an analysis of the disk protonization timescale.

From the equations of lepton number and baryon number conservation, we compute

$R=\nabla_{\mu}(n_{e}u^{\mu})=\nabla_{\mu}(n_{\rm b}u^{\mu}Y_{e})=n_{\rm
b}u^{\mu}\nabla_{\mu}Y_{e},$ (38)

and thus

$u^{\mu}\nabla_{\mu}Y_{e}=R\frac{m_{\rm b}}{\rho},$ (39)

where $u^{\mu}$ is the four-velocity of the fluid and $R$ denotes a source
term to account for the change in lepton number due to weak interactions.
Equation (39) shows that along a fluid trajectory, i.e., in the comoving frame
of the fluid, the electron fraction changes by a rate of $Rm_{\rm b}/\rho$. We
can therefore define the characteristic timescale for the composition of
matter to change by $t_{Y_{e}}\equiv Y_{e}\rho/(Rm_{\mathrm{b}})$. Neglecting
neutrino absorption, appropriate for the disk midplane, one has

$t_{Y_{e}}=\frac{Y_{e}}{R_{\bar{\nu_{e}}}^{\mathrm{eff}}-R_{\nu_{e}}^{\mathrm{eff}}}\frac{\rho}{m_{\rm
b}}.$ (40)

Ignoring final state blocking in the neutrino phase space, the effective
emission rates scale as (Tubbs & Schramm, 1975; Bruenn, 1985)
$R_{\bar{\nu}_{e},\nu_{e}}^{\mathrm{eff}}\propto\rho
T^{5}\propto\dot{M}^{9/4}$, where we have used the approximate expression Eq.
(14) for the midplane temperature in the second step. Furthermore, since
$\dot{M}\propto\rho$ (cf. Eq. (11)), we deduce that, approximately,

$t_{Y_{e}}\propto\dot{M}^{-5/4},$ (41)

which explains the decrease in protonization timescale with increasing
accretion rate seen in the outer accretion disks in Figs. 12 and 13.

The inner accretion disk shows qualitatively different behavior depending on
whether the accretion disk resides in a state above or below the ignition
threshold (Eq. (2), Sec. 4.1.3). We find that above the ignition threshold,
the disk midplane density $\rho\propto\dot{M}$ (cf. Eq. (11)) is sufficiently
large that electrons become degenerate, as shown in Fig. 12. This, in turn,
suppresses positron creation via $\gamma\rightarrow e^{+}+e^{-}$ and thus Eq.
(37) relative to (36), and hence leads to slight neutronization even in very
neutron-rich initial conditions (see MD_M03 in Fig. 13). This neutronization
mechanism is also why collapsar accretion disks starting with much higher
initial $Y_{e}$ are able to generate neutron-rich outflows and synthesize
r-process elements (Siegel et al., 2019). However, the inner part of the
accretion disks cannot become arbitrarily degenerate and neutron rich, thanks
to a self-regulation mechanism discussed in Siegel & Metzger (2018a) and first
pointed out by Beloborodov (2003) (see also Chen & Beloborodov (2007)). As
shown in Fig. 12, disk self-regulation leads to a heating-cooling balance that
results in moderate electron degeneracy $\eta\sim 1$ and corresponding neutron
richness of $Y_{e}\approx 0.1$. While the disk is in such a self-regulated,
moderately degenerate phase, the inner part of the accretion disk feeds highly
neutron-rich material into the outflows. At lower $\dot{M}$, close to the
ignition threshold, degeneracy and self-regulation become somewhat weaker and
less pronounced; however, MD_M002 still qualitatively shows the same behavior
as MD_M03.

As a result of viscous spreading, the disk density decreases over time, which
may suppress degeneracy eventually late in the disk’s evolution. The onset of
a decrease in degeneracy is noticeable from the snapshots in Fig. 12. However,
by the time massive disks such as MD_M03 reach this break-down of degeneracy,
most of the ejecta material has already been launched and the change in disk
composition has little effect on the overall composition of ejecta.

Below the ignition threshold, we find an inverted scenario in terms of
midplane composition in the inner parts of the accretion disk. The timescale
for protonization in the innermost part of the accretion disk is significantly
smaller than in the outer parts, resulting in high-$Y_{e}$ material
surrounding the black hole already at $t\sim 100$ ms (cf. MD_M0002 in Fig.
12). We attribute this to excess viscous heating in the absence of
energetically significant neutrino cooling. The resulting high-entropy
environment (cf. Figs. 13 and 14) leads to a prolific generation of positrons
via pair production, $\gamma\rightarrow e^{+}+e^{-}$, and thus decreases the
timescale for protonization via Eq. (37). This innermost region in MD_M0002 is
too small, however, to feed material into the outflows sufficient enough to
form a pronounced high-$Y_{e}$ tail of the ejecta (cf. Sec. 4.3).

Figure 15 shows a snapshot of $Y_{e}$ in the meridional plane on large spatial
scales $\lesssim 700$ km. In general, the proton fraction is enhanced in polar
outflows, in agreement with enhanced pair production in the low-density polar
funnels (Fig. 4). A radial $Y_{e}$ gradient is also present, which results
from the gradual protonization of the outer disks over time according to Eq.
(41). As a result of a long ‘engine lifetime’ on a timescale $0.1-1$ s,
continuous ejection of wind material leads to significant mixing of outflow
material before it enters a homologous expansion phase. Three-dimensional
radiation transport calculations will be required to compute detailed
predictions for kilonova lightcurves and spectra.

###  4.3 Nucleosynthesis

Figure 16: Final elemental abundances at $10^{9}$ s for the accretion disk
ejecta of our three simulation runs. Shown are mean abundances of all unbound
tracer particles (excluding the disk relaxation phase). Also plotted for
reference are the observed solar system abundances from Arnould et al. (2007),
scaled to match the MD_M03 abundance at $A=195$.  Figure 17: Mass fractions of
the nucleosynthetic yields in the first (summed over $A=70-90$), second
(summed over $A=125-135$), and third (summed over $A=186-203$) r-process peaks
for the disk ejecta of simulation runs MD_M03, MD_M002, and MD_M0002 presented
here. A qualitative change in nucleosynthesis products across the ignition
threshold for weak interactions (Eq. (2)) is evident.  Figure 18: Mass
distribution of unbound outflows from our simulated accretion disks according
to the electron fraction $Y_{\rm e}$, extracted from tracer particles at 5 GK.
A significant high-$Y_{e}$ tail is established across the ignition threshold.

The neutron-rich outflows generated by our post-merger accretion disks are
sites for r-process nucleosynthesis. We record thermodynamic properties of the
ejecta by tracer particles distributed throughout the disks initially (cf.
Sec. 4.1.4), and use such tracer profiles as input to nuclear reaction network
calculations. The initial conditions that determine the outcome of the
r-process are mostly set by the electron fraction, specific entropy, and
expansion timescale of the flow at $T\approx 5$ GK (Lippuner & Roberts, 2015),
which is the characteristic temperature at which nuclear statistical
equilibrium breaks down and neutron-capture reactions set in.

Neutrino absorption is taken into account approximately by a ring-like ‘light-
bulb’ scheme in post-processing (Fernández & Metzger, 2013), which irradiates
the ejecta with neutrino luminosities as computed in Sec. 4.1.3. We refer to
Fernández & Metzger (2013) and Siegel & Metzger (2018a) for more details on
this approach.

Starting at 10 GK, we perform full nuclear reaction network calculations on
the tracer particles using the reaction network SkyNet (Lippuner & Roberts,
2017), in order to track nuclear abundances as the outflows undergo r-process
nucleosynthesis. Figure 16 shows the final total abundance yields from all
tracers after $10^{9}$ s for all three disk simulation runs. These abundances
are the result of integrated ejecta that are unbound over the course of the
simulation runs, most of which originate from the inner accretion disk region
with radii between $\approx(100-250)$ km (cf. Fig. 11).

Final r-process abundance patterns qualitatively change across the ignition
threshold for weak interactions in post-merger disks (see Fig. 16). At high
accretion rates, such as run MD_M03, we find abundances in good agreement with
residual solar r-process abundances (Arnould et al., 2007) across the entire
range of mass numbers from the first ($A\approx 80$) to the third ($A\approx
195$) r-process peak (see Siegel & Metzger (2018a) for a more detailed
discussion on this case). While across all accretion rates we find good
agreement with the solar abundance pattern between the second peak ($A\approx
130$) and the third r-process peak, run MD_M002 (close to the ignition
threshold) already shows somewhat suppressed first-to-second peak elements,
and light r-process elements are strongly suppressed below the ignition
threshold (run MD_M0002).

The qualitative change in nucleosynthesis patterns across the ignition
threshold is further illustrated by Fig. 17, which shows the mass fractions of
nuclei grouped into first-, second-, and third-peak elements as a function of
accretion rate. While 2nd-peak elements are synthesized in roughly constant
amounts across the $\dot{M}$ range investigated here, abundances of light
r-process nuclei rise strongly toward and around the ignition threshold, while
3rd-peak nuclei are somewhat overproduced in that regime.

This behavior is a result of a qualitative change in the distribution of
$Y_{e}$ of the ejecta across the ignition threshold as evident from Fig. 18.
While at low $\dot{M}$, below the ignition threshold, the $Y_{e}$ distribution
of the unbound outflows at 5 GK is still centered around the initial value of
$Y_{e}=0.1$ (run MD_M0002), run MD_M002 close to the ignition threshold has
already built up a significant tail toward higher electron fraction; finally,
run MD_M03 shows a broad distribution of $Y_{e}$ at least up to $\lesssim
0.4$. The emerging high $Y_{e}$ tail with increasing disk masses across the
ignition threshold was also reported by Metzger et al. (2009). This effect can
be ascribed in part to the strongly changing protonization timescale as a
function of accretion rate across the ignition threshold (cf. Eq. (41); Sec.
4.2). There is a possibility of our tracer particles slightly under-resolving
the high-$Y_{e}$ tail—a significantly larger number of tracers, not considered
here due to computational cost and feasibility, may better resolve the effects
of, e.g., an inner high-$Y_{e}$ part as in the MD_M0002 disk, but is unlikely
to qualitatively change the conclusions reached here. Toward high $\dot{M}$,
neutrino absorption by the ejecta contributes to flattening the $Y_{e}$
distribution as well. This effect is already noticeable for MD_M03 and is
expected to become more important at even higher-$\dot{M}$ disks (Miller et
al., 2019b), which are beyond the scope of the present paper.

##  5 Discussion & Conclusions

Our simulation results on post-merger accretion disks address a number of key
questions regarding r-process nucleosynthesis and the electromagnetic
counterpart observables of LIGO-Virgo’s BNS and NSBH merger source
populations. In view of future gravitational-wave detections, we explore the
parameter space of BNS and NSBH systems by presenting mappings between binary
parameters and accretion disk as well as ejecta masses, based on current fits
to numerical relativity simulations (Sec. 2.2). We estimate that disk ejecta
may dominate over dynamical ejecta across most of the parameter space for NSBH
systems and BNS systems with a light-to-intermediate secondary neutron star.
Dynamical ejecta may dominate in BNS systems with either very high mass ratio
or a heavy secondary star.

We argue that to lowest order, the nucleosynthesis and kilonova signatures of
post-merger accretion disks are determined by a single parameter, the
(initial) accretion rate $\dot{M}$, which is largely determined by the initial
mass of the disk. We demonstrate the existence of an ignition threshold
$\dot{M}_{\rm ign}$ for weak interactions by employing a simplified 1D
analytical disk model and self-consistent 3D GRMHD simulations with weak
interactions. This threshold divides the 1D parameter space of post-merger
accretion disks into two distinct regions (radiative efficient versus
advection dominated) that lead to qualitatively different behavior with regard
to nucleosynthesis and kilonova emission. We explore this change in behavior
by investigating three disk models across this threshold in detail. Some key
findings from these simulations include:

  * •

We find that the radiative efficiency—the fraction of accreted rest mass
energy turned into radiation per unit time—is very small below the ignition
threshold ($\lesssim\\!10^{-3}$). It rises to almost the theoretical maximum
($\gtrsim 0.1$ for the configurations considered here) just above the ignition
threshold $\dot{M}_{\rm ign}\sim 1\times 10^{-3}M_{\odot}s^{-1}$, and finally
decreases with increasing $\dot{M}$ (or initial disk mass) as a result of
increasing optical depth in the accretion flow that starts to trap neutrinos
(Fig. 10).

  * •

Protonization timescales $t_{Y_{e}}$ for the outer (non-degenerate) parts of
our accretion flows vary strongly, with approximate scaling
$t_{Y_{e}}\propto\dot{M}^{-5/4}$ (Figs. 12 and 13). This leads to massive
accretion disks acquiring a high-$Y_{e}$ tail of ejecta much faster than
lighter accretion disks (Fig. 18). For the inner disks, we report a
qualitative difference in compositional evolution across the ignition
threshold. For high-mass disks above the threshold, high accretion rates cause
dense accretion flows and force electrons to become degenerate. Weak
interactions then lead to a self-regulation mechanism that keeps the disk at
moderate degeneracy $\eta\sim 1$ and high neutron richness $Y_{e}\approx 0.1$.
For the low-mass disk below the ignition threshold studied here, negligible
neutrino cooling causes a ‘pile-up’ of heat in the inner disks, which triggers
$e^{+}e^{-}$ pair production. The injection of additional $e^{+}$ strongly
accelerates protonization of the neutron-rich disk material in the vicinity of
the black hole (Figs. 12 and 13).

  * •

All three disk models investigated here produce large amounts of heavy
$r$-process elements beyond the second r-process peak—corresponding to red
kilonova components (cf. Figs. 16 and 17). Our highest-$\dot{M}$ model is able
to produce light $r$-process elements as well, with an overall r-process
pattern in good agreement with solar abundances. This is the result of $(i)$
rapid protonization of the disks at high accretion rates and $(ii)$ self-
irradiation of the outflows due to strong neutrino radiation. Our results
indicate that disks close to or below the ignition threshold show heavily
suppressed light r-process patterns. These disks largely keep their initial
composition—weak interactions are suppressed to a level that protonization
occurs on timescales longer than ejecta generation. Blue kilonova components
are thus only expected for massive disks well above the ignition threshold.

  * •

The effect of an ignition threshold is also imprinted on the bulk properties
of the disk ejecta. We find that disks above the ignition threshold eject
$\sim 30-50\%$ of their initial disk mass, whereas disks below the ignition
threshold can eject $\gtrsim 60\%$ (estimates are to be taken with at least
10% uncertainty). This increase is a result of the strong heating-cooling
imbalance throughout the disk and its corona in absence of energetically
significant cooling (Figs. 13 and 14). The increase of viscous heating with
decreasing disk mass also enhances the high-velocity tail of the disk ejecta
(Fig. 11).

  * •

Outflows driven by MHD turbulence, with an effective $\alpha$ viscosity
parameter and thus accretion rate set self-consistently (cf. Sec. 4.1.1), are
generated early on during the disk evolution. We note that $\sim\\!50\%$ of
the total ejecta has already left the detection sphere at 1000 km radius
within the simulated time frame, a fraction of the effective viscous timescale
of the disks (cf. Tab. 2, Sec. 4.1.4). This is in contrast to
$\alpha$-viscosity disks, in which ejecta is generated on much longer
(viscous) timescales (Fernández et al., 2019, 2020; Kyutoku et al., 2020).
This has important consequences for ejecta composition and kilonova colors:
for MHD disks there is less time for ejecta material to protonize and produce
light r-process elements and blue kilonova components.

  * •

_Actinide-boost stars:_ One immediate conclusion from our simulations is that
disks below the ignition threshold retain their original composition (which
may be arbitrarily neutron-rich, depending on the initial cold NS matter) and
can thus contribute to _actinide-boost events_ —ejecta with an overabundance
of actinides relative to lanthanides. Environments polluted by such ejecta
could be conducive to forming actinide-boost stars. Roughly $\sim\\!30\%$ of
all r-process enhanced stars show such an overabundance (Mashonkina et al.,
2014). This has so far been attributed to NS merger events with dominating and
very neutron-rich tidal ejecta (e.g., Holmbeck et al. 2019). Here we
demonstrate that, depending on the EOS, post-merger accretion disks below the
ignition threshold could contribute to actinide overabundance to a similar or
even higher level than dynamical ejecta (cf. right panel of Fig. 1, and
considering that such disks may evaporate $>60\%$ of disk material into
outflows). Above the ignition threshold, the lanthanide-heavy disk ejecta will
tend to ‘wash out’ possible actinide overabundance from dynamical ejecta,
unless lanthanide production is severely suppressed at very high accretion
rates (Miller et al., 2019a; Li & Siegel, 2021). Furthermore, for NSBH
systems, our parameter study suggests that actinide overabundance may mainly
result from post-merger accretion disks, rather than dynamical ejecta (cf.
Fig. 2). This is because disk ejecta may dominate over most of the parameter
space. Actinide overabundance can thus be realized for disks below the
ignition threshold, while above the ignition threshold potential actinide
overabundance from dynamical ejecta may be (over)compensated by lanthanide-
heavy ejecta from the disk (with the same caveat at high accretion rates
mentioned for BNS systems above). This is also relevant to inferences about
metal-poor star populations (e.g., Holmbeck et al. 2020).

Finally, we note a few caveats and limitations of the present study:

  * •

The late-time viscous regime of the outflows is not simulated here for reasons
of computational cost. These late-time outflows may moderately enhance the
high-$Y_{e}$ tail of the ejecta distribution (Fernández et al., 2019, 2020;
Metzger et al., 2009, 2008a). However, given that the disks have already
undergone significant viscous spreading (cf. Fig. 7), we expect weak
interactions to be close to freeze-out already, so that compositional changes
may be moderate (in particular for disks below the ignition threshold).

  * •

Additional variation in the total mass of MHD outflows on the order of
$\sim\\!10\%$ may result from variation of initial magnetic field
configurations (Christie et al., 2019a), which we do not explore here, but
rather defer to future work.

  * •

At accretion rates higher than simulated here, more accurate neutrino
transport is required to reliably capture the effect of neutrino self-
irradiation (Miller et al., 2019a; Li & Siegel, 2021).

The current suite of simulations represents a first step in exploring the
parameter space of post-merger accretion disks in GRMHD. Future studies along
the lines of Fernández et al. (2020) are required to investigate the
sensitivity of ejecta properties to black-hole mass and spin, initial
compactness of the disk, initial magnetic field configurations, etc., to a
level that is not yet encoded in the accretion rate parameter, and to probe a
wider range in $\dot{M}$ systematically.

The authors thank Brian Metzger, Jonah Miller, Duncan Brown, Eric Coughlin,
and Stefan Ballmer for discussions. This work was performed in part at the
Kavli Institute for Theoretical Physics; SD and DMS thank the organizers and
participants of the KITP program “The New Era of Gravitational-Wave Physics
and Astrophysics” for stimulating discussions and KITP for hospitality. This
research was supported in part by the National Science Foundation under Grant
No. NSF PHY-1748958. Resources supporting this work were provided by the NASA
High-End Computing (HEC) Program through the NASA Advanced Supercomputing
(NAS) Division at Ames Research Center and by Syracuse University. This
research was enabled inpart by support provided by SciNet (www.scinethpc.ca)
and Compute Canada (www.computecanada.ca). SD acknowledges support from the
Research Excellence Doctoral Fellowship at Syracuse University and the
Director’s Postdoctoral Fellowship at Los Alamos National Laboratory—with the
Laboratory Directed Research and Development Program project number
20200687PRD2. Los Alamos National Laboratory is operated by Triad National
Security, LLC, for the National Nuclear Security Administration of U.S.
Department of Energy under Contract No. 89233218CNA000001. DMS acknowledges
the support of the Natural Sciences and Engineering Research Council of Canada
(NSERC), funding reference number RGPIN-2019-04684. Research at Perimeter
Institute is supported in part by the Government of Canada through the
Department of Innovation, Science and Economic Development Canada and by the
Province of Ontario through the Ministry of Colleges and Universities.

##  Appendix A Ignition threshold and simplified 1D disk model

For completeness, we provide the expression for $F(x,\chi_{\rm BH})$ that
appears in the simplified one-dimensional $\alpha$-disk model described in
Sec. 2.3:

$\displaystyle F(x,\chi_{\rm BH})$ $\displaystyle\equiv$
$\displaystyle\frac{x^{3}+\chi_{\rm BH}}{(x^{3}-3x+2\chi_{\rm
BH})^{1/2}x^{3/2}}\Bigg{[}(x-x_{0})-\frac{3}{2}\chi_{\rm
BH}\ln\frac{x}{x_{0}}-\frac{3(x_{1}-\chi_{\rm
BH})^{2}}{x_{1}(x_{1}-x_{2})(x_{1}-x_{3})}\ln\left(\frac{x-x_{1}}{x_{0}-x_{1}}\right)$
(A1) $\displaystyle-\frac{3(x_{2}-\chi_{\rm
BH})^{2}}{x_{2}(x_{2}-x_{1})(x_{2}-x_{3})}\ln\left(\frac{x-x_{2}}{x_{0}-x_{2}}\right)-\frac{3(x_{3}-\chi_{\rm
BH})^{2}}{x_{3}(x_{3}-x_{1})(x_{3}-x_{2})}\ln\left(\frac{x-x_{3}}{x_{0}-x_{3}}\right)\Bigg{]}.$

Here, $x\equiv\sqrt{2r/r_{\rm g}}$, $x_{0}$ corresponds to the location of the
marginally stable orbit, and $x_{1}$, $x_{2}$, $x_{3}$ are the roots of
$x^{3}-3x+2\chi_{\rm BH}=0$. Explicitly,
$x_{1}=2\cos(\frac{1}{3}\cos^{-1}\chi_{\rm BH}-\pi/3)$,
$x_{2}=2\cos(\frac{1}{3}\cos^{-1}\chi_{\rm BH}+\pi/3)$, and
$x_{3}=-2\cos(\frac{1}{3}\cos^{-1}\chi_{\rm BH})$. See also Page & Thorne
(1974) and Chen & Beloborodov (2007).

## References

  * Abbott et al. (2017a) Abbott, B. P., et al. 2017a, Phys. Rev. Lett., 119, 161101.
  * Abbott et al. (2017b) —. 2017b, ApJ, 848, L12
  * Abbott et al. (2019) —. 2019, PRX, 9, 031040.
  * Abbott et al. (2020a) —. 2020a, ApJ, 892, L3.
  * Abbott et al. (2020b) Abbott, R., et al. 2020b, arXiv:2010.14527
  * Abramowicz et al. (1997) Abramowicz, M. A., Lanza, A., & Percival, M. J. 1997, ApJ, 479, 179.
  * Aloy et al. (2005) Aloy, M. A., Janka, H.-T., & Müller, E. 2005, A&A, 436, 273
  * Antier et al. (2020) Antier, S., Agayeva, S., Aivazyan, V., et al. 2020, MNRAS, 492, 3904.
  * Armitage (2011) Armitage, P. J. 2011, ARA&A, 49, 195
  * Arnould et al. (2007) Arnould, M., Goriely, S., & Takahashi, K. 2007, Phys. Rep., 450, 97
  * Babiuc-Hamilton et al. (2019) Babiuc-Hamilton, M., Brandt, S. R., Diener, P., et al. 2019, The Einstein Toolkit, “Somerville” ET_2015_11 and “Tesla” ET_2018_02 versions, Zenodo doi:10.5281/zenodo.3522086.
  * Balbus (2003) Balbus, S. A. 2003, ARA&A, 41, 555
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
  * Balbus & Hawley (1998) —. 1998, Rev. Mod. Phys. , 70, 1
  * Balbus & Hawley (2002) —. 2002, ApJ, 573, 749
  * Bardeen et al. (1972) Bardeen, J. M., Press, W. H., & Teukolsky, S. A. 1972, ApJ, 178, 347
  * Beloborodov (1998) Beloborodov, A. M. 1998, MNRAS, 297, 739.
  * Beloborodov (2003) —. 2003, ApJ, 588, 931
  * Berger (2014) Berger, E. 2014, Annu. Rev. Astron. Astrophys., 52, 43
  * Bernuzzi et al. (2014) Bernuzzi, S., Dietrich, T., Tichy, W., & Brügmann, B. 2014, Phys. Rev. D, 89, 104021
  * Bruenn (1985) Bruenn, S. W. 1985, ApJS, 58, 771.
  * Burbidge et al. (1957) Burbidge, E. M., Burbidge, G. R., Fowler, W. A., & Hoyle, F. 1957, , 29, 547
  * Cameron (1957) Cameron, A. G. W. 1957, AJ, 62, 9
  * Capano et al. (2020) Capano, C. D., Tews, I., Brown, S. M., et al. 2020, Nature Astronomy, 1, publisher: Nature Publishing Group.
  * Chandrasekhar (1960) Chandrasekhar, S. 1960, Proc. Natl. Acad. Sci., 46, 253
  * Chen & Beloborodov (2007) Chen, W.-X., & Beloborodov, A. M. 2007, ApJ, 657, 383
  * Christie et al. (2019a) Christie, I., Lalakos, A., Tchekhovskoy, A., et al. 2019a, Mon. Not. Roy. Astron. Soc., 490, 4811
  * Ciolfi et al. (2017) Ciolfi, R., Kastaun, W., Giacomazzo, B., et al. 2017, Phys. Rev. D, 95, 063016
  * Ciolfi et al. (2019) Ciolfi, R., Kastaun, W., Kalinani, J. V., & Giacomazzo, B. 2019, Physical Review D, 100, 023005.
  * Ciolfi & Vijay Kalinani (2020) Ciolfi, R., & Vijay Kalinani, J. 2020, arXiv e-prints, 2004, arXiv:2004.11298.
  * Coughlin et al. (2020) Coughlin, M. W., Dietrich, T., Antier, S., et al. 2020, MNRAS, 492, 863.
  * Cowling (1933) Cowling, T. G. 1933, MNRAS, 94, 39.
  * De et al. (2018) De, S., Finstad, D., Lattimer, J. M., et al. 2018, Phys. Rev. Lett., 121, 091102, [Erratum: Phys.Rev.Lett. 121, 259902 (2018)]
  * Dessart et al. (2009) Dessart, L., Ott, C. D., Burrows, A., Rosswog, S., & Livne, E. 2009, ApJ, 690, 1681
  * Di Matteo et al. (2002) Di Matteo, T., Perna, R., & Narayan, R. 2002, ApJ, 579, 706
  * Fernández et al. (2020) Fernández, R., Foucart, F., & Lippuner, J. 2020, MNRAS, 497, 3221.
  * Fernández et al. (2015) Fernández, R., Kasen, D., Metzger, B. D., & Quataert, E. 2015, Mon. Not. R. Astron. Soc., 446, 750
  * Fernández & Metzger (2013) Fernández, R., & Metzger, B. D. 2013, MNRAS, 435, 502
  * Fernández et al. (2019) Fernández, R., Tchekhovskoy, A., Quataert, E., Foucart, F., & Kasen, D. 2019, MNRAS, 482, 3373.
  * Foucart (2012) Foucart, F. 2012, Phys. Rev. D, 86, 124007
  * Foucart et al. (2018) Foucart, F., Hinderer, T., & Nissanke, S. 2018, Physical Review D, 98, 081501.
  * Fujibayashi et al. (2017) Fujibayashi, S., Sekiguchi, Y., Kiuchi, K., & Shibata, M. 2017, ApJ, 846, 114\.
  * Fujibayashi et al. (2020) Fujibayashi, S., Shibata, M., Wanajo, S., et al. 2020, arXiv e-prints, 2001, arXiv:2001.04467.
  * Galeazzi et al. (2013) Galeazzi, F., Kastaun, W., Rezzolla, L., & Font, J. A. 2013, Phys. Rev. D, 88, 064009
  * Gammie & Popham (1998) Gammie, C. F., & Popham, R. 1998, The Astrophysical Journal, 498, 313.
  * Goodale et al. (2003) Goodale, T., Allen, G., Lanfermann, G., et al. 2003, in Vector and Parallel Processing – VECPAR’2002, 5th International Conference, Lecture Notes in Computer Science (Berlin: Springer).  http://edoc.mpg.de/3341
  * Halevi & Mösta (2018) Halevi, G., & Mösta, P. 2018, MNRAS, 477, 2366.
  * Hawley & Balbus (1992) Hawley, J. F., & Balbus, S. A. 1992, ApJ, 400, 595
  * Holmbeck et al. (2020) Holmbeck, E. M., Frebel, A., McLaughlin, G. C., et al. 2020, arXiv e-prints, arXiv:2010.01621
  * Holmbeck et al. (2019) Holmbeck, E. M., Sprouse, T. M., Mumpower, M. R., et al. 2019, ApJ, 870, 23
  * Hosseinzadeh et al. (2019) Hosseinzadeh, G., Cowperthwaite, P. S., Gomez, S., et al. 2019, ApJ, 880, L4.
  * Hotokezaka et al. (2013a) Hotokezaka, K., Kiuchi, K., Kyutoku, K., et al. 2013a, Phys. Rev. D, 87, 024001
  * Hotokezaka et al. (2013b) Hotokezaka, K., Kyutoku, K., Tanaka, M., et al. 2013b, ApJ, 778, L16
  * Jiang et al. (2014) Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2014, ApJ, 796, 106
  * Just et al. (2015) Just, O., Bauswein, A., Pulpillo, R. A., Goriely, S., & Janka, H.-T. 2015, MNRAS, 448, 541
  * Just et al. (2021) Just, O., Goriely, S., Janka, H.-T., Nagataki, S., & Bauswein, A. 2021, arXiv:2102.08387
  * Kasen et al. (2017) Kasen, D., Metzger, B., Barnes, J., Quataert, E., & Ramirez-Ruiz, E. 2017, Nature, 551, 80
  * Kastaun et al. (2013) Kastaun, W., Galeazzi, F., Alic, D., Rezzolla, L., & Font, J. A. 2013, Phys. Rev. D, 88, 021501
  * Kawanaka & Mineshige (2007) Kawanaka, N., & Mineshige, S. 2007, ApJ, 662, 1156.
  * Kiuchi et al. (2009) Kiuchi, K., Sekiguchi, Y., Shibata, M., & Taniguchi, K. 2009, Phys. Rev. D, 80, 064037
  * Kohri et al. (2005) Kohri, K., Narayan, R., & Piran, T. 2005, Astrophys. J., 629, 341
  * Krüger & Foucart (2020) Krüger, C. J., & Foucart, F. 2020, Phys. Rev. D, 101, 103002
  * Kulkarni (2005) Kulkarni, S. R. 2005, arXiv e-prints, astro.
  * Kumar & Zhang (2015) Kumar, P., & Zhang, B. 2015, Phys. Rep., 561, 1
  * Kyutoku et al. (2020) Kyutoku, K., Fujibayashi, S., Hayashi, K., et al. 2020, ArXiv e-prints, arXiv:2001.04474.
  * Landry et al. (2020) Landry, P., Essick, R., & Chatziioannou, K. 2020, arXiv:2003.04880
  * Li & Paczyński (1998) Li, L.-X., & Paczyński, B. 1998, ApJ, 507, L59
  * Li & Siegel (2021) Li, X., & Siegel, D. M. 2021, ArXiv e-prints, arXiv:2103.02616.
  * Lippuner & Roberts (2015) Lippuner, J., & Roberts, L. F. 2015, Astrophys. J., 815, 82
  * Lippuner & Roberts (2017) —. 2017, ApJS, 233, 18
  * Löffler et al. (2012) Löffler, F., Faber, J., Bentivegna, E., et al. 2012, Class. Quant. Grav. , 29, 115001
  * Lundquist et al. (2019) Lundquist, M. J., Paterson, K., Fong, W., et al. 2019, ApJ, 881, L26.
  * Mashonkina et al. (2014) Mashonkina, L., Christlieb, N., & Eriksson, K. 2014, A&A, 569, A43
  * Metzger et al. (2009) Metzger, B., Piro, A., & Quataert, E. 2009, Mon. Not. Roy. Astron. Soc., 396, 304
  * Metzger (2019) Metzger, B. D. 2019, Living Reviews in Relativity, 23, 1.
  * Metzger et al. (2008a) Metzger, B. D., Piro, A. L., & Quataert, E. 2008a, MNRAS, 390, 781\.
  * Metzger et al. (2008b) Metzger, B. D., Thompson, T. A., & Quataert, E. 2008b, Astrophys. J., 676, 1130
  * Metzger et al. (2018) —. 2018, ApJ, 856, 101
  * Metzger et al. (2010) Metzger, B. D., Martínez-Pinedo, G., Darbha, S., et al. 2010, MNRAS, 406, 2650
  * Miller et al. (2019a) Miller, J. M., Sprouse, T. M., Fryer, C. L., et al. 2019a, arXiv e-prints, arXiv:1912.03378.
  * Miller et al. (2019b) Miller, J. M., Ryan, B. R., Dolence, J. C., et al. 2019b, Phys. Rev. D, 100, 023008.
  * Miller et al. (2019c) Miller, M., et al. 2019c, Astrophys. J. Lett., 887, L24
  * Mösta et al. (2014) Mösta, P., Mundim, B. C., Faber, J. A., et al. 2014, Class. Quant. Grav. , 31, 015005
  * Narayan et al. (2001) Narayan, R., Piran, T., & Kumar, P. 2001, Astrophys. J., 557, 949
  * Nedora et al. (2019) Nedora, V., Bernuzzi, S., Radice, D., et al. 2019, ApJ, 886, L30.
  * Nitz et al. (2019a) Nitz, A. H., Capano, C., Nielsen, A. B., et al. 2019a, Astrophys. J., 872, 195
  * Nitz et al. (2019b) Nitz, A. H., Dent, T., Davies, G. S., et al. 2019b, Astrophys. J., 891, 123
  * Novikov & Thorne (1973) Novikov, I. D., & Thorne, K. S. 1973, in Black Holes (Les Astres Occlus), 343–450
  * Oechslin et al. (2007) Oechslin, R., Janka, H.-T., & Marek, A. 2007, A&A, 467, 395
  * Page & Thorne (1974) Page, D. N., & Thorne, K. S. 1974, ApJ, 191, 499.
  * Paschalidis et al. (2015) Paschalidis, V., Ruiz, M., & Shapiro, S. L. 2015, ApJ, 806, L14
  * Penna et al. (2013) Penna, R. F., S¸adowski, A., Kulkarni, A. K., & Narayan, R. 2013, MNRAS, 428, 2255
  * Popham et al. (1999) Popham, R., Woosley, S. E., & Fryer, C. 1999, Astrophys. J., 518, 356
  * Qian & Woosley (1996) Qian, Y.-Z., & Woosley, S. E. 1996, ApJ, 471, 331
  * Radice et al. (2016) Radice, D., Galeazzi, F., Lippuner, J., et al. 2016, MNRAS, 460, 3255
  * Rezzolla et al. (2010) Rezzolla, L., Baiotti, L., Giacomazzo, B., Link, D., & Font, J. A. 2010, Class. Quant. Grav. , 27, 114105
  * Riley et al. (2019) Riley, T. E., et al. 2019, Astrophys. J. Lett., 887, L21
  * Rosswog et al. (1999) Rosswog, S., Liebendörfer, M., Thielemann, F.-K., et al. 1999, A&A, 341, 499\.
  * Ruffert et al. (1996) Ruffert, M., Janka, H.-T., & Schaefer, G. 1996, A&A, 311, 532
  * Ruffert et al. (1997) Ruffert, M., Janka, H.-T., Takahashi, K., & Schaefer, G. 1997, A&A, 319, 122
  * Ruiz et al. (2016) Ruiz, M., Lang, R. N., Paschalidis, V., & Shapiro, S. L. 2016, ApJ, 824, L6.
  * Schnetter et al. (2004) Schnetter, E., Hawley, S. H., & Hawke, I. 2004, Class. Quantum Grav., 21, 1465
  * Sekiguchi et al. (2016) Sekiguchi, Y., Kiuchi, K., Kyutoku, K., Shibata, M., & Taniguchi, K. 2016, Phys. Rev. D, 93, 124046
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, Astron. Astrophys., 24, 337
  * Shibata & Taniguchi (2006) Shibata, M., & Taniguchi, K. 2006, Phys. Rev. D, 73, 064027
  * Siegel (2019) Siegel, D. M. 2019, Eur. Phys. J. A, 55, 203.
  * Siegel et al. (2019) Siegel, D. M., Barnes, J., & Metzger, B. D. 2019, Nature, 569, 241.
  * Siegel et al. (2013) Siegel, D. M., Ciolfi, R., Harte, A. I., & Rezzolla, L. 2013, Phys. Rev. D, 87, 121302
  * Siegel et al. (2014) Siegel, D. M., Ciolfi, R., & Rezzolla, L. 2014, ApJ, 785, L6
  * Siegel & Metzger (2017) Siegel, D. M., & Metzger, B. D. 2017, Phys. Rev. Lett., 119, 231102
  * Siegel & Metzger (2018a) —. 2018a, ApJ, 858, 52
  * Thornburg (2004) Thornburg, J. 2004, Class. Quantum Grav., 21, 743
  * Thorne (1974) Thorne, K. S. 1974, ApJ, 191, 507.
  * Timmes & Arnett (1999) Timmes, F. X., & Arnett, D. 1999, ApJS, 125, 277
  * Timmes & Swesty (2000) Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501
  * Tubbs & Schramm (1975) Tubbs, D. L., & Schramm, D. N. 1975, Astrophysical Journal, 201, 467
  * Velikhov (1959) Velikhov, E. P. 1959, Sov. Phys. JETP, 36, 995
  * Winteler et al. (2012) Winteler, C., Käppeli, R., Perego, A., et al. 2012, Astrophys. J. Lett., 750, L22