# HARM3D+NUC: A new method for simulating the post-merger phase of binary
neutron star mergers with GRMHD, tabulated EOS and neutrino leakage

Ariadna Murguia-Berthier  Department of Astronomy and Astrophysics, University
of California, Santa Cruz, CA 95064, USA  DARK, Niels Bohr Institute,
University of Copenhagen, Blegdamsvej 17, 2100 Copenhagen, Denmark  Scott C.
Noble  Gravitational Astrophysics Lab, NASA Goddard Space Flight Center,
Greenbelt, MD 20771, USA  Luke F. Roberts  NSCL, Michigan State University,
East Lansing, MI 48824, USA  Enrico Ramirez-Ruiz  Department of Astronomy and
Astrophysics, University of California, Santa Cruz, CA 95064, USA  DARK, Niels
Bohr Institute, University of Copenhagen, Blegdamsvej 17, 2100 Copenhagen,
Denmark  Leonardo R. Werneck  Center for Gravitational Waves and Cosmology,
West Virginia University, Chestnut Ridge Research Building, Morgantown, WV
26505  Michael Kolacki  Center for Computational Relativity, Rochester
Institute of Technology, Rochester, New York 14623, USA  School of Physics and
Astronomy, Rochester Institute of Technology, Rochester, New York 14623, USA
Zachariah B. Etienne  Department of Physics and Astronomy, West Virginia
University, Morgantown, WV 26506  Center for Gravitational Waves and
Cosmology, West Virginia University, Chestnut Ridge Research Building,
Morgantown, WV 26505  Mark Avara  Center for Computational Relativity and
Gravitation, Rochester Institute of Technology, 85 Lomb Memorial Drive,
Rochester, New York 14623, USA  Manuela Campanelli  Center for Computational
Relativity and Gravitation, Rochester Institute of Technology, 85 Lomb
Memorial Drive, Rochester, New York 14623, USA  School of Mathematical
Sciences, Rochester Institute of Technology, Rochester, New York 14623, USA
School of Physics and Astronomy, Rochester Institute of Technology, Rochester,
New York 14623, USA  Riccardo Ciolfi  INAF, Osservatorio Astronomico di
Padova, Vicolo dell’Osservatorio 5, I-35122 Padova, Italy  INFN, Sezione di
Padova, Via Francesco Marzolo 8, I-35131 Padova, Italy  Federico Cipolletta
Leonardo Corporate LABS - via Raffaele Pieragostini 80, 16149 Genova GE -
Italy  Brendan Drachler  Center for Computational Relativity and Gravitation,
Rochester Institute of Technology, Rochester, New York 14623, USA  School of
Physics and Astronomy, Rochester Institute of Technology, Rochester, New York
14623, USA  Lorenzo Ennoggi  Center for Computational Relativity and
Gravitation, Rochester Institute of Technology, 85 Lomb Memorial Drive,
Rochester, New York 14623, USA  Joshua Faber  Center for Computational
Relativity and Gravitation, Rochester Institute of Technology, 85 Lomb
Memorial Drive, Rochester, New York 14623, USA  School of Mathematical
Sciences, Rochester Institute of Technology, Rochester, New York 14623, USA
School of Physics and Astronomy, Rochester Institute of Technology, Rochester,
New York 14623, USA  Grace Fiacco  Center for Computational Relativity and
Gravitation, Rochester Institute of Technology, Rochester, New York 14623, USA
School of Physics and Astronomy, Rochester Institute of Technology, Rochester,
New York 14623, USA  Bruno Giacomazzo  Universitá degli Studi di Milano -
Bicocca, Dipartimento di Fisica G. Occhialini, Piazza della Scienza 3, I-20126
Milano, Italy  INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3,
I-20126 Milano, Italy  INAF, Osservatorio Astronomico di Brera, via E. Bianchi
46, I-23807 Merate (LC), Italy  Tanmayee Gupte  Center for Computational
Relativity and Gravitation, Rochester Institute of Technology, Rochester, New
York 14623, USA  School of Physics and Astronomy, Rochester Institute of
Technology, Rochester, New York 14623, USA  Trung Ha  Center for Computational
Relativity and Gravitation, Rochester Institute of Technology, Rochester, New
York 14623, USA  Bernard J. Kelly  Department of Physics, University of
Maryland Baltimore County, 1000 Hilltop Circle Baltimore, MD 21250, USA
Gravitational Astrophysics Lab, NASA Goddard Space Flight Center, Greenbelt,
MD 20771, USA  Center for Research and Exploration in Space Science and
Technology, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA  Julian
H. Krolik  Physics and Astronomy Department, Johns Hopkins University,
Baltimore, MD 21218, USA  Federico G. Lopez Armengol  Center for Computational
Relativity and Gravitation, Rochester Institute of Technology, Rochester, New
York 14623, USA  Ben Margalit  Astronomy Department and Theoretical
Astrophysics Center, University of California, Berkeley, Berkeley, CA 94720,
USA  Tim Moon  Center for Computational Relativity and Gravitation, Rochester
Institute of Technology, 85 Lomb Memorial Drive, Rochester, New York 14623,
USA  School of Mathematical Sciences, Rochester Institute of Technology,
Rochester, New York 14623, USA  Richard O’Shaughnessy  Center for
Computational Relativity and Gravitation, Rochester Institute of Technology,
85 Lomb Memorial Drive, Rochester, New York 14623, USA  School of Mathematical
Sciences, Rochester Institute of Technology, Rochester, New York 14623, USA
School of Physics and Astronomy, Rochester Institute of Technology, Rochester,
New York 14623, USA  Jesús M. Rueda-Becerril  Center for Computational
Relativity and Gravitation, Rochester Institute of Technology, 85 Lomb
Memorial Drive, Rochester, New York 14623, USA  Jeremy Schnittman
Gravitational Astrophysics Lab, NASA Goddard Space Flight Center, Greenbelt,
MD 20771, USA  Yossef Zenati  Physics and Astronomy Department, Johns Hopkins
University, Baltimore, MD 21218, USA  Yosef Zlochower  Center for
Computational Relativity and Gravitation, Rochester Institute of Technology,
85 Lomb Memorial Drive, Rochester, New York 14623, USA  School of Mathematical
Sciences, Rochester Institute of Technology, Rochester, New York 14623, USA
School of Physics and Astronomy, Rochester Institute of Technology, Rochester,
New York 14623, USA

###### Abstract

The first binary neutron star merger has already been detected in
gravitational waves. The signal was accompanied by an electromagnetic
counterpart including a kilonova component powered by the decay of radioactive
nuclei, as well as a short $\gamma$-ray burst. In order to understand the
radioactively-powered signal, it is necessary to simulate the outflows and
their nucleosynthesis from the post-merger disk. Simulating the disk and
predicting the composition of the outflows requires general relativistic
magnetohydrodynamical (GRMHD) simulations that include a realistic, finite-
temperature equation of state (EOS) and self-consistently calculating the
impact of neutrinos. In this work, we detail the implementation of a finite-
temperature EOS and the treatment of neutrinos in the GRMHD code HARM3D+NUC,
based on HARM3D. We include formal tests of both the finite-temperature EOS
and the neutrino leakage scheme. We further test the code by showing that,
given conditions similar to those of published remnant disks following neutron
star mergers, it reproduces both recombination of free nucleons to a neutron-
rich composition and excitation of a thermal wind.

accretion disks- general relativistic magnetohydrodynamical simulations,
neutrino leakage

##  1 Introduction

On August 17, 2017, the LIGO/VIRGO collaboration detected the first
gravitational wave signal arising from the merger of two neutron stars (Abbott
et al., 2017a). This signal was accompanied by a counterpart observed all over
the electromagnetic spectrum (Abbott et al., 2017b; Murguia-Berthier et al.,
2017; Coulter et al., 2017; Shappee et al., 2017). This event, named GW170817,
gave credence to the idea that at least a subset of neutron star mergers give
rise to short $\gamma$-ray bursts (sGRBs; Eichler et al., 1989; Narayan et
al., 1992; Lee & Ramirez-Ruiz, 2007; Nakar, 2007).

In order to understand the electromagnetic emission, we need to study the
properties of merger. After the two neutron stars merge, the fate of the
remnant depends on the final mass of the resulting object. If the final mass
is less than the mass allowed for an object with rigid rotation, then the
remnant will be a stable neutron star. On the other hand, if the final mass is
larger, then it can result in a hot hyper-massive neutron star (HMNS),
supported by differential rotation, or it can promptly collapse to a black
hole (Shibata & Taniguchi, 2006; Baiotti et al., 2008; Ravi & Lasky, 2014). In
both cases, the compact object will be surrounded by an accretion disk
(Eichler et al., 1989; Baiotti et al., 2008). If the result is an HMNS, there
will be transport of mass and angular momentum from the inner edge to the
outer edge that will drive the HMNS to rigid rotation, where it can either
remain stable, or undergo a delayed collapse to a black hole (BH; see Nakar,
2019, for a recent review). It is widely believed that GW170817 resulted in a
delayed collapse to a black hole (Margalit & Metzger, 2017). In any case, the
compact object is left surrounded by an accretion disk containing highly
neutron-rich material (Lee & Ramirez-Ruiz, 2007).

The post-merger accretion disk will be entirely opaque to photons (Popham et
al., 1999; Narayan et al., 2001; Lee et al., 2004, 2005, 2009). As we go
deeper in the disk, due to the high density and temperature, neutrinos (and
anti-neutrinos) will be created via the charged $\beta$-process, electron-
positron annihilation, and plasmon decay (Narayan et al., 2001; Di Matteo et
al., 2002; Chen & Beloborodov, 2007). In the region where the neutrinos are
created, matter will be optically thin to neutrinos. In even deeper regions,
matter will be optically thick to neutrinos. In the optically thin region,
free neutrinos will carry energy away, and cool the disk, making it
geometrically thinner (Chevalier, 1989; Houck & Chevalier, 1991).

Further out in the the disk, where neutrinos are no longer created in
substantial numbers, free nucleons will recombine into $\alpha$-particles. The
photons will still be trapped in the disk, therefore the disk will be thicker
and radiatively inefficient (Popham et al., 1999; Narayan et al., 2001; Lee et
al., 2004, 2005, 2009). An outflow arises due to instabilites in the accretion
disk from its magnetic field (Balbus & Hawley, 1998). The instabilities will
transport angular momentum at significant rates, dissipating energy and
driving a high velocity outflow. In addition, the recombination of free
nucleons into $\alpha-$particles is capable of unbinding part of the material
from the disk (Lee et al., 2009; Fernández & Metzger, 2013a).

Aside from material ejected from the disk, there are other outflows from the
binary merger that will significantly contribute to the electromagnetic
emission, including a dynamical ejecta (see, for example Rosswog et al., 1999;
Fernández et al., 2015; Radice et al., 2016), and a neutrino-driven wind
(Dessart et al., 2009; Fernández & Metzger, 2013a; Perego et al., 2014; Kasen
et al., 2017; Fernández et al., 2017). As the different outflows expand and
cool down, heavy elements are synthesised via the rapid neutron capture
process ($r$-process) (Freiburghaus et al., 1999; Kulkarni, 2005; Fernández &
Metzger, 2013b; Lippuner & Roberts, 2015; Palenzuela et al., 2015; Radice et
al., 2016; Roberts et al., 2017; Fernández et al., 2017; Lippuner et al.,
2017; Radice et al., 2018; Zenati et al., 2019; Radice et al., 2020). After
neutrons are exhausted, elements will radioactively decay and heat the
surrounded material, which will thermally emit in the optical/IR bands (Li &
Paczyński, 1998; Metzger et al., 2010; Roberts et al., 2011; Kasen et al.,
2013; Barnes & Kasen, 2013; Tanaka & Hotokezaka, 2013; Grossman et al., 2014;
Kasen et al., 2015; Barnes et al., 2016; Rosswog et al., 2017; Kasen & Barnes,
2019; Siegel, 2019); in particular, see Metzger (2019) and references within.
This emission, called a _kilonova_ , was detected for GW170817 (Drout et al.,
2017; Kilpatrick et al., 2017; Soares-Santos et al., 2017; Tanvir et al.,
2017; Smartt et al., 2017; Nicholl et al., 2017; Cowperthwaite et al., 2017;
Villar et al., 2017; Kasen et al., 2017; Pian et al., 2017; Kasliwal et al.,
2019). It is predicted that if the composition of the ejecta includes
lanthanides, the emission tends to be more red and peak at later times,
whereas if there are no third peak elements, the emission tends to be bluer
and peaks earlier (Barnes & Kasen, 2013; Tanaka & Hotokezaka, 2013).
Understanding the nucleosynthesis, and the amount of mass ejected is therefore
important when deciding the best strategy to observe and perform surveys for
kilonovae. This paper will focus on the disk ejecta.

The key parameter in determining the rate of nucleosynthesis, and in
particular whether third peak $r$-process elements (including the lanthanides)
are created in the disk ejecta, is the electron fraction of the ejected
material (Kasen et al., 2013; Lippuner & Roberts, 2015; Roberts et al., 2017;
Lippuner et al., 2017; Kasen et al., 2017; Just et al., 2021). The problem is
that the composition of these ejecta varies between different simulations with
results ranging from compositions dominated by iron peak elements to ejecta
dominated by lanthanides (e.g., Janiuk, 2014; Fernández et al., 2015; Foucart
et al., 2018; Janiuk, 2019; Siegel & Metzger, 2018; Miller et al., 2019b). One
of the significant differences between the simulations is in the neutrino
treatment. Neutrinos carry away energy and lepton number, altering the
electron fraction and the final ejecta mass and they significantly alter the
composition of the ejected material. Thus, simulations need to model the
composition and thermodynamic state of the ejecta as realistically as possible
to understand and model the kilonova emission.

In order to model the post-merger disk, we need to self-consistently include
multiple relevant physical processes. Due to the compact nature of the BH, we
need to consider general relativity (GR). Due to the importance of the
magnetic stresses we need to include magneto-hydrodynamics (MHD).
Additionally, to self-consistently include the addition of neutrinos and
recombination energy, we need both a realistic equation of state (EOS) and a
way in which to consider the impact of neutrinos in the optically thick and
thin regions.

There have been many previous efforts to simulate a black hole surrounded by
an accretion disk in the context of a binary neutron star. A brief (and
certainly incomplete) summary of the numerical efforts is below.

Numerical simulations initially added neutrino physics by adding pressure
terms in the EOS and adding emission and heating/cooling terms from weak
reactions in hydrodynamical simulations (Popham et al., 1999; Narayan et al.,
2001; Di Matteo et al., 2002; Kohri & Mineshige, 2002; Lee et al., 2004, 2005;
Metzger et al., 2008; Zalamea & Beloborodov, 2011). There have been general
relativistic magnetohydrodynamical (GRMHD) simulations in 2d with analytical
terms for the neutrino pressure with approximations by Di Matteo et al. (2002)
that also include nuclear reactions using the GRMHD code HARM2D (Janiuk et
al., 2013; Janiuk, 2014, 2019). There have been efforts performing simulations
of binary neutron stars, or a hyper-massive NS, with an accretion disk in 3d
with GRMHD but without neutrinos (for example, Siegel et al., 2014; Kiuchi et
al., 2014, 2015; Dionysopoulou et al., 2015; Ruiz et al., 2016; Ciolfi et al.,
2017; Kiuchi et al., 2018; Ruiz et al., 2018). Also, groups simulated disks
after the merger of binary NS including GR with some kind of neutrino
transport but including no magnetic fields (Foucart et al., 2016; Fujibayashi
et al., 2017; Nedora et al., 2021). Other groups performed hydrodynamical
calculations with neutrino physics, including neutrino leakage schemes and a
transport scheme but no magnetic fields (Ruffert et al., 1996; Rosswog &
Liebendörfer, 2003; Metzger & Fernández, 2014; Perego et al., 2014; Martin et
al., 2015; Fernández et al., 2015; Just et al., 2015).

Foucart et al. (2015, 2018) performed general relativistic hydrodynamical
(GRHD) simulations and compared different neutrino treatments, including
neutrino transport and leakage schemes. Additionally, Siegel & Metzger (2018)
and De & Siegel (2020) performed GRMHD simulations of a magnetized torus with
a neutrino leakage scheme and the Helmholtz equation of state. Hossein Nouri
et al. (2018) compared 3d simulations of magnetized and unmagnetized accretion
disks with GRMHD including a neutrino leakage scheme. Li & Siegel (2021)
performed an M1 scheme with neutrino conversions. There have also been GRMHD
simulations that included a tabulated EOS with neutrino transport using Monte-
Carlo methods (Miller et al., 2019a, b).

In this paper, we present simulations using HARM3D+NUC, based on HARM3D,
considering the impact of neutrinos through a leakage scheme and a multi-
component, finite-temperature EOS. HARM3D is a versatile GRMHD code that has
been well tested and used in many astrophysical scenarios. It uses arbitrary
coordinates, allowing for a more accurate conservation of angular momentum.
Additionally, it has copious analysis tools developed over the years. The
addition of a neutrino leakage scheme and tabulated EOS into HARM3D+NUC is a
stepping stone that allows for further advances. The paper is structured as
follows: in Section 2 we discuss how we implemented the realistic EOS and the
leakage scheme. In Section 3 we describe the tests we performed to validate
the implementation of the tabulated EOS including a torus in hydrostatic
equilibrium. In Section 4 we describe the tests we performed to validate the
leakage scheme, and in Section 5 we use both the tabulated EOS and leakage
scheme to better simulate a torus with a magnetic field.

##  2 Methods

In order to accurately simulate accretion disks, we need the ability to solve
the general relativistic magnetohydrodynamics (GRMHD) equations with a
realistic equation of state (EOS) and a way to account for the effect
neutrinos and anti-neutrinos have on the material’s energy and electron
fraction. In this section, we explain how we added a tabulated EOS and
neutrino leakage scheme to HARM3D, in a new code called HARM3D+NUC.

###  2.1 HARM3D+NUC

HARM3D (Gammie et al., 2003; Noble et al., 2006, 2009) solves the GRMHD
equations in conservative form. HARM3D is a well tested code that can handle
arbitrary coordinate systems, which allows for less numerical diffusion and
better conservation of angular momentum when using coordinate systems that
more closely conform to local symmetries of the problem (Zilhão & Noble,
2014). Below we set $G=c=1$. The GRMHD equations of motion include the baryon
conservation equation,

$\nabla_{\mu}\left(n_{b}u^{\mu}\right)=0\quad,$ (1)

the energy-momentum conservation equations (with a heating/cooling source,
neglecting momentum transfer)

$\nabla_{\mu}{T^{\mu}}_{\nu}=\mathcal{Q}u_{\nu}\quad,$ (2)

and Maxwell’s equations

$\nabla_{\nu}{{}^{{}^{*}}\\!\\!F}^{\mu\nu}=0\quad,$ (3)
$\nabla_{\nu}F^{\mu\nu}=J^{\mu}\quad,$ (4)

where $u^{\mu}$ is the 4-velocity of the fluid, $\mathcal{Q}$ is the energy
change rate per volume in the comoving fluid frame (due to neutrino
heating/cooling), $n_{b}$ is the number density of baryons, $F^{\mu\nu}$ is
the Faraday tensor times $1/\sqrt{4\pi}$, ${{}^{{}^{*}}\\!\\!F}^{\mu\nu}$ is
the dual of this tensor or the Maxwell tensor times $1/\sqrt{4\pi}$, and
$J^{\mu}$ is the 4-current111We follow Gammie et al. (2003) in our definition
of the electromagnetic field tensor and magnetic field variables.. In
practice, we don’t use Eq. (4), since we work in the limit of ideal MHD. The
change in the conservation of lepton number is

$\nabla_{\mu}\left(n_{e}u^{\mu}\right)=\mathcal{R}/m_{b}\quad,$ (5)

where $n_{e}$ is the number density of electrons,
$\mathcal{R}=-\mathcal{R}_{\nu_{e}}+\mathcal{R}_{{\bar{\nu}_{e}}}$ is the
difference in the net rate of neutrino and anti-neutrino number per volume in
the comoving fluid frame.

Note that the rest-mass density of the gas (mass per unit volume) is dominated
by the baryon mass, $\rho\approx m_{b}n_{b}$, where $m_{b}$ is the baryon
mass. The baryon number conservation equation can then be replaced by the
regular continuity equation:

$0=m_{b}\nabla_{\mu}\left(n_{b}u^{\mu}\right)=\nabla_{\mu}\left(m_{b}n_{b}u^{\mu}\right)=\nabla_{\mu}\left(\rho
u^{\mu}\right)\quad.$ (6)

Instead of using $n_{e}$ and $n_{b}$, we may use the fluid density $\rho$ and
the electron fraction $Y_{e}$:

$Y_{e}\equiv\frac{n_{e}}{n_{b}}=\frac{n_{e}}{\rho/m_{b}}=\frac{m_{b}n_{e}}{\rho}$
(7)

or $Y_{e}\rho=m_{b}n_{e}$ and we can therefore multiply Eq. (5) by $m_{b}$ to
yield the electron fraction equation:

$\nabla_{\mu}\left(\rho Y_{e}u^{\mu}\right)=\mathcal{R}\quad.$ (8)

The total stress-energy tensor is the sum of the fluid part,

$T^{\mu\nu}_{\mathrm{fluid}}=\rho hu^{\mu}u^{\nu}+Pg^{\mu\nu},$ (9)

and the electromagnetic part

$\displaystyle T^{\mu\nu}_{\mathrm{EM}}$ $\displaystyle=$ $\displaystyle
F^{\mu\lambda}{F^{\nu}}_{\lambda}-\frac{1}{4}g^{\mu\nu}F^{\lambda\kappa}F_{\lambda\kappa}$
(10) $\displaystyle=$
$\displaystyle{{||b||^{2}}}u^{\mu}u^{\nu}+\frac{1}{2}{{||b||^{2}}}g^{\mu\nu}-b^{\mu}b^{\nu}\quad,$
(11)

where we adopt the ideal MHD condition

$u_{\lambda}F^{\lambda\kappa}=0\quad,$ (12)

and where $g_{\mu\nu}$ is the metric, $h=\left(1+\epsilon+P/\rho\right)$ is
the specific enthalpy, $P$ is the pressure, $\epsilon$ is the specific
internal energy density, $b^{\mu}={{}^{{}^{*}}\\!\\!F}^{\nu\mu}u_{\nu}$ is the
magnetic field 4-vector, and ${{||b||^{2}}}\equiv b^{\mu}b_{\mu}$ is twice the
magnetic pressure $P_{m}$.

Equations (2-6) can be expressed in flux conservative form

${\partial}_{t}{\bf U}\left({{\mathbf{P}}}\right)=-{\partial}_{i}{\bf
F}^{i}\left({{\mathbf{P}}}\right)+\mathbf{S}\left({{\mathbf{P}}}\right)\,$
(13)

where ${\bf U}$ is a vector of “conserved” variables, ${\bf F}^{i}$ are the
fluxes, $\mathbf{S}$ is a vector of source terms, and ${{\mathbf{P}}}$ is the
vector of primitive variables. Explicitly, these are

${{\mathbf{P}}}=\left[\rho,{{\mathcal{B}}}^{k},\tilde{u}^{i},Y_{e},T\right]^{T}$
(14) ${\bf U}\left({{\mathbf{P}}}\right)=\sqrt{-g}\left[\rho
u^{t},{T^{t}}_{t}+\rho u^{t},{T^{t}}_{j},B^{k},\rho Y_{e}u^{t}\right]^{T}$
(15) ${\bf F}^{i}\left({{\mathbf{P}}}\right)=\sqrt{-g}\left[\rho
u^{i},{T^{i}}_{t}+\rho
u^{i},{T^{i}}_{j},\left(b^{i}u^{k}-b^{k}u^{i}\right),\rho
Y_{e}u^{i}\right]^{T}$ (16)
$\mathbf{S}\left({{\mathbf{P}}}\right)=\sqrt{-g}\left[0,{T^{\kappa}}_{\lambda}{\Gamma^{\lambda}}_{t\kappa}+\mathcal{Q}u_{t},{T^{\kappa}}_{\lambda}{\Gamma^{\lambda}}_{j\kappa}+\mathcal{Q}u_{i},0,\mathcal{R}\right]^{T}\,,$
(17)

where $g$ is the determinant of the metric, ${\Gamma^{\lambda}}_{\mu\kappa}$,
is the metric’s affine connection, $T$ is the temperature, and
$B^{i}={{\mathcal{B}}}^{i}/\alpha={{}^{{}^{*}}\\!\\!F}^{it}$ is the magnetic
field.

The primitive velocity is the flow’s 4-velocity projected into a frame moving
orthogonal to the space-like hypersurface:

$\tilde{u}^{\mu}=\left({\delta^{\mu}}_{\nu}+n^{\mu}n_{\nu}\right)u^{\nu}$ (18)

which only has spatial coefficients

$\tilde{u}^{i}=u^{i}+\alpha\gamma g^{ti}\quad,$ (19)

where $\alpha=1/\sqrt{-g^{tt}}$ is the lapse function,
$\beta^{i}=-g^{ti}/g^{tt}$ is the shift function, $\gamma=\alpha u^{t}$ is the
Lorentz factor, and $n^{\mu}$ is the 4-velocity of the orthogonal frame:
$n_{\mu}=[-\alpha,0,0,0]$ and $n^{\mu}=[1/\alpha,-\beta^{i}/\alpha]^{T}$.
Defining a fluid three-velocity $v^{i}=\tilde{u}^{i}/\gamma$, it can be shown
that $\gamma=1/\sqrt{1-v^{2}}$, where $v^{2}=v_{i}v^{i}$.

###  2.2 Implementation of a tabulated EOS in HARM3D+NUC

In the following section, we describe the implementation of a tabulated EOS in
HARM3D+NUC.

The tables and routines for interpolating tabulated quantities are provided
by222The link to the tabulated EOS is the following:
https://stellarcollapse.org/SROEOS, and the link to the interpolation routines
is: https://bitbucket.org/zelmani/eosdrivercxx/src O’Connor & Ott (2010) and
Schneider et al. (2017). The finite-temperature tables give thermodynamic
variables, including, for example, the sound speed, and the chemical
potentials of the nucleons, electrons/positrons and neutrinos/anti-neutrinos,
as a function of the temperature ($T$), the electron fraction ($Y_{e}$), and
the rest-mass density ($\rho$). The linear interpolation routines are provided
by O’Connor & Ott (2010) and Schneider et al. (2017). The interpolation is
done in $\log T$, $\log\rho$, $Y_{e}$ space for $\log\epsilon$, $\log P$, and
the rest of the thermodynamical variables.

The tables consider an interpolation between a single nucleus approximation
(SNA) in the high density regime and nuclear statistical equilibrium (NSE) of
several nucleides in the low density regime. The SNA is composed of free
nucleons, electrons, positrons, $\alpha-$particles, and photons. In the high
density regime, nuclei are included using the liquid drop model. The regimes
are smoothly interpolated. Using the tables, we have the advantage that the
nuclear binding energy release due to recombination energy from the
$\alpha$-particles is included.

There are three main calls to the EOS in HARM3D+NUC:

  * •

We call the EOS when setting the characteristic velocity in order to solve the
Riemann problem (Gammie et al., 2003). The wave velocities depend on the
relativistic sound speed (Gammie et al., 2003), which can be interpolated
directly from the tables.

  * •

We replaced the primitive variable $u=\rho\epsilon$ with the temperature as a
reconstructed variable, which makes the interpolation of the pressure faster
as all independent variables are known and can be used to perform the
interpolation immediately.. This means that we call the EOS to obtain the
primitive energy density $u$ after we update $\rho$, $T$, and $Y_{e}$ from the
conservation equations.

  * •

We call the EOS repeatedly when converting from conserved variables to
primitive variables.

Our implementation of a tabulated EOS into the conserved to primitive
variables routine in HARM3D+NUC follows Siegel et al. (2018).

####  2.2.1 Primary recovery: 3d routine

The primary recovery routine follows a 3-parameter root-finding method similar
to ones implemented in Cerdá-Durán et al. (2008); Siegel et al. (2018). We
call this routine the ‘3d’ routine. For this routine, we reduce the GRMHD
equations into three equations that have three unknowns, allowing us to solve
the following system:

$\tilde{Q}^{2}=\bigg{(}1-\frac{1}{\gamma^{2}}\bigg{)}(\mathcal{B}^{2}+W)^{2}-\frac{(Q_{\mu}\mathcal{B}^{\mu})^{2}(\mathcal{B}^{2}+2W)}{W^{2}}$
(20)
$Q_{\mu}n^{\mu}=-\frac{{{\mathcal{B}}}^{2}}{2}\bigg{(}2-\frac{1}{\gamma^{2}}\bigg{)}+\frac{(Q_{\mu}{{\mathcal{B}}}^{\mu})^{2}}{2W^{2}}-W+P(\rho,Y_{e},T)$
(21) $\epsilon=\epsilon(\rho,Y_{e},T)\quad.$ (22)

Using these equations, we perform Newton-Raphson iterations until we obtain
sufficiently accurate values for the independent variables $\gamma$, $T$ and
$W$. Here $Q_{\mu}=-n_{\nu}T^{\nu}_{\mu}=\alpha T^{t}_{\mu}$, $W$ is related
to the specific enthalpy through $W=h\rho\gamma^{2}$,

$\tilde{Q}^{\mu}={j^{\mu}}_{\nu}Q^{\nu}$,
$j_{\mu\nu}=g_{\mu\nu}+n_{\mu}n_{\nu}$, and $P$ is the pressure interpolated
from tables.

####  2.2.2 Backup recovery 1: 2d routine

We also implemented backup routines that recover the conserved variables. One
of them follows an optimized version of the ”2d” method of Noble et al.
(2006). We call this routine the ‘2d’ routine. In this routine, the
independent variables are $W$ and $v^{2}$, found using equations (20-21). The
previous time step’s set of primitive variables are used as initial guesses to
the Newton-Raphson procedure. As was done in Siegel et al. (2018), we obtain
the pressure and the temperature for each $W$ and $v^{2}$. This is done by
first constructing the specific enthalpy: $h=W/(\gamma^{2}\rho)$, which can
also be constructed with quantities from the EOS tables: $h(\rho,T,Y_{e})$.
Then, with the density, the electron fraction and the specific enthalpy, we
perform a Newton-Raphson method to obtain the temperature from the tables,
solving the equation: $h=h(\rho,T,Y_{e})$. Note that this inversion is time
expensive, which is why this routine is slower than the 3d routine.

####  2.2.3 Backup recovery 2: 2d ’safe-guess’ routine

If there is non-convergence for this backup routine, we include an initial
‘safe guess’ as described in Cerdá-Durán et al. (2008). We call this routine
the ‘2d safe guess’ routine. In this scenario, we use the upper limits of the
EOS table to obtain the maximum thermodynamical quantities:

$\rho_{\rm max}=D,$ (23) $T_{\rm max}=T_{\rm max,tables},$ (24) $P_{\rm
max}=P(\rho_{\rm max},Y_{e},T_{\rm max})\quad.$ (25)

Were $D$ is the density measured in the orthogonal frame:

$D\equiv-\rho n_{\mu}u^{\mu}=\gamma\rho\quad.$ (26)

Then we can estimate the initial ‘safe guess’ for the root-finding procedure:

$\gamma_{\rm guess}=\gamma_{\rm max}=50,$ (27) $W_{\rm
guess}=Q_{\mu}n^{\mu}+P_{\rm max}-\frac{{{\mathcal{B}}}^{2}}{2}\quad.$ (28)

####  2.2.4 Backup recovery 3: 2d dog leg routine

If the ’safe guess’ option does not converge, this routine includes a backup
root-finding method: a trust-region, dog leg routine that is more robust than
a Newton-Raphson (Press et al., 1992; Powell, 1968). We call this routine the
‘2d dog leg’ routine.

####  2.2.5 Backup recovery 4: ’Palenzuela’ routine

If all else fails, we use the routine described in Palenzuela et al. (2015).
This routine solves a 1d equation using the Brent method. In this routine,
called ’Palenzuela’, the independent variable is a rescaled variable

$x_{\rm pal}\equiv\frac{\rho h\gamma^{2}}{\rho\gamma}\ .$ (29)

We use the auxiliary rescaled variables:

$\displaystyle q_{\rm pal}\equiv\frac{-(Q_{\mu}n^{\mu}+D)}{D}$ ,
$\displaystyle\ r_{\rm pal}\equiv\frac{\tilde{Q}^{2}}{D^{2}}\ ,\ $ (30)
$\displaystyle s_{\rm pal}\equiv\frac{\mathcal{B}^{2}}{D}$ , $\displaystyle\
t_{\rm pal}\equiv\frac{Q_{\mu}{{\mathcal{B}}}^{\mu}}{D^{3/2}}\ .$ (31)

The independent variable should be bracketed between:

$1+q_{\rm pal}-s_{\rm pal}>x_{\rm pal}>2+2q_{\rm pal}-s_{\rm pal}\ .$ (32)

The method uses an initial guess for $x_{\rm pal}$ from the previous time
step, and gets approximate quantities. Using them, it updates $x_{\rm pal}$
and iterates again until convergence is reached. The method is the following
(where approximate quantities will be denoted by a hat):

We obtain an approximate Lorentz factor $\hat{\gamma}^{-2}$:

$\hat{\gamma}^{-2}=1-\frac{x_{\rm pal}^{2}r_{\rm pal}+(2x_{\rm pal}+s_{\rm
pal})t_{\rm pal}^{2}}{x_{\rm pal}^{2}(x_{\rm pal}+s_{\rm pal})^{2}}\ .$ (33)

With that, we can estimate:

$\hat{\rho}=\frac{D}{\hat{\gamma}}$ (34)

and an approximate specific energy:

$\hat{\epsilon}=\hat{\gamma}-1+\frac{x_{\rm
pal}}{\hat{\gamma}}(1-\hat{\gamma}^{2})+\hat{\gamma}\bigg{(}q_{\rm pal}-s_{\rm
pal}+\frac{t_{\rm pal}^{2}}{2x_{\rm pal}^{2}}+\frac{s_{\rm
pal}}{2\hat{\gamma}^{2}}\bigg{)}\ .$ (35)

A call to the EOS will give the pressure
$\hat{P}(\hat{\rho},\hat{\epsilon},Y_{e})$, and with all those approximate
quantities, we can solve for $x_{\rm pal}$ using the Brent method by solving:

$0=f(x_{\rm pal})=x_{\rm
pal}-\hat{\gamma}\bigg{(}1+\hat{\epsilon}+\frac{\hat{P}}{\hat{\rho}}\bigg{)}\
.$ (36)

We repeat the estimation of all the hat quantities until the solution for
$x_{\rm pal}$ converges.

Figure 1: Relative error comparing primitive variables created from a grid of
density and temperature after we performed the conversion from conserved
variables to primitive variables. The primitive variables were created with
$Y_{e}=0.1$, the Lorentz factor $\gamma=2$, $\log{\frac{P_{\rm mag}}{P_{\rm
gas}}}=-5$, and a Minkowski metric. We perturbed them by $5\%$ and then
recovered them using our conserved to primitive routines. The error is
calculated by summing over the relative error of each primitive variable
compared to the original. We did this for $2^{14}$ points in the shown range.
The 2d routines failed only once, the 3d routines failed 11 times, and the
Palenzuela routine did not fail in this range. Density is in units of ${\rm
g/cm^{3}}$, and temperature is in units of Kelvin. Here we compare different
routines, described in the text.

###  2.3 Neutrino leakage scheme

In the following section, we describe how we implemented a leakage scheme that
takes into account the heating/cooling due to neutrinos, as well as how their
emission and absorption affect the electron fraction. This leakage scheme is
suited to describe the contribution of neutrinos to the composition, and
energy.

####  2.3.1 Rates

The scheme calculates the absorption/emission rate as well as the energy loss
rates due to neutrinos. We use these rates in the source terms of Eq. (2) and
Eq. (5). The scheme uses energy-averaged quantities.

Like Ruffert et al. (1996); Galeazzi et al. (2013); Siegel & Metzger (2018),
we consider the following neutrino reactions, each with their own
absorption/emission rate (which has units of $\rm{cm}^{-3}{\rm s}^{-1}$) and
the energy loss rate rate due to neutrinos (with units of
$\rm{erg}\,\rm{cm}^{-3}{\rm s}^{-1}$):

  * •

Charged $\beta$-process with $\mathcal{R}^{\beta}_{\nu_{i}}$ and
$\mathcal{Q}^{\beta}_{\nu_{i}}$:

$e^{-}+p\rightarrow n+\nu_{e}$ (37) $e^{+}+n\rightarrow p+\bar{\nu}_{e}$ (38)

  * •

Plasmon decay with $\mathcal{R}^{\gamma}_{\nu_{i}}$ and
$\mathcal{Q}^{\gamma}_{\nu_{i}}$:

$\gamma\rightarrow\nu_{e}+\bar{\nu}_{e}$ (39)
$\gamma\rightarrow\nu_{x}+\bar{\nu}_{x}$ (40)

where x is the muon and tauon, and in this case, $\gamma$ corresponds to a
photon.

  * •

Electron-positron pair annihilation with $\mathcal{R}^{ee}_{\nu_{i}}$ and
$\mathcal{Q}^{ee}_{\nu_{i}}$

$e^{-}+e^{+}\rightarrow\nu_{e}+\bar{\nu}_{e}$ (41)
$e^{-}+e^{+}\rightarrow\nu_{x}+\bar{\nu}_{x}\ .$ (42)

Using the above reactions, we calculate the total number emission in the
optically thin regime from species $i$ as (Ruffert et al., 1996):

$\mathcal{R}_{\nu_{i}}=\mathcal{R}^{\beta}_{\nu_{i}}+\mathcal{R}^{\gamma}_{\nu_{i}}+\mathcal{R}^{ee}_{\nu_{i}}$
(43)

and the total energy loss rate rate in the optically thin regime is:

$\mathcal{Q}_{\nu_{i}}=\mathcal{Q}^{\beta}_{\nu_{i}}+\mathcal{Q}^{\gamma}_{\nu_{i}}+\mathcal{Q}^{ee}_{\nu_{i}}\quad,$
(44)

where ”$i$” denotes the different neutrino/anti-neutrino flavors: electron, or
muon and tauon.

The total emission/absorption rates and the energy loss rates are given by an
interpolation between the diffusive optically thick regime and the transparent
optically thin regime (Ruffert et al., 1996):

$\mathcal{R}_{\nu_{i}}^{\rm eff}=\mathcal{R}_{\nu_{i}}\left(1+\frac{t_{\rm
diff}}{t_{\rm emission,\mathcal{R}}}\right)^{-1}$ (45)
$\mathcal{Q}_{\nu_{i}}^{\rm eff}=\mathcal{Q}_{\nu_{i}}\bigg{(}{1+\frac{t_{\rm
diff}}{t_{\rm emission,\mathcal{Q}}}}\bigg{)}^{-1}.\ .$ (46)

Here the diffusion timescale is given by:

$t_{\rm diff}=\frac{D_{\rm diff}\tau^{2}}{c\kappa_{\nu_{i}}}\ ,$ (47)

where $D_{\rm diff}=6$ (Rosswog & Liebendörfer, 2003; O’Connor & Ott, 2010;
Siegel & Metzger, 2018) and $\tau$ is the optical depth, and
$\kappa_{\nu_{i}}$ the energy averaged opacity (in units of
$\mathrm{cm}^{-1}$) of $\nu_{i}$. The absorption/emission and energy loss
timescales are $t_{\rm
emission,\mathcal{R}}=\mathcal{R}_{\nu_{i}}/n_{\nu_{i}}$, with $n_{\nu_{i}}$
being the neutrino number density (at chemical equilibrium), and $t_{\rm
emission,\mathcal{Q}}=\mathcal{Q}_{\nu_{i}}/\varepsilon_{\nu_{i}}$, with
$\varepsilon_{\nu_{i}}$ being the neutrino energy density. In the optically
thick regime, the neutrino loss rate is less than the diffusion time, which
results in $\mathcal{R}_{\nu_{i}}^{\rm eff}=n_{\nu_{i}}/t_{\rm diff}$ and
$\mathcal{Q}_{\nu_{i}}^{\rm eff}=\varepsilon_{\nu_{i}}/t_{\rm diff}$, whereas
in the optically thin regime, we recover the rates from equation (43) and
(44). The rates for the muon and tauon neutrinos/anti-neutrinos estimated in
Ruffert et al. (1996) take into account all four of those species. We also
note that several quantities, including the chemical potentials, are obtained
from EOS table interpolation.

####  2.3.2 Optical depth

The transition between the two regimes will be set by the optical depth
$\tau_{\nu_{i}}$, which is also needed to obtain the diffusion timescale. In
order to get the optical depth, we consider the following reactions as the
source of neutrino opacity:

$\nu_{e}+n\rightarrow p+e^{-}$ (48) $\bar{\nu}_{e}+p\rightarrow n+e^{+}\ $
(49) $\nu_{i}+p\rightarrow\nu_{i}+p$ (50) $\nu_{i}+n\rightarrow\nu_{i}+n.$
(51)

The opacities are obtained from Ruffert et al. (1996). Electron scattering is
neglected.

The usual global approach to calculate the optical depth of a point in the
flow would be to integrate the opacity over all directions and determine the
path of minimal absorption. This approach assumes that the neutrino will
follow a straight path. However, we follow Neilsen et al. (2014); Siegel &
Metzger (2018), where a local, iterative approach is used instead of a global
calculation, and where crooked minimal paths are acceptable. The optical depth
is calculated by obtaining the shortest path of the neutrino out of the star
using its neighbors. For the first timestep, we begin by initializing the
optical depth grid to zero. Next, we perform the first iteration, where we
estimate the optical depth at each cell as the minimum of the optical depth of
its neighbor ($\tau_{{\nu_{i}},\rm{neighbor}}$) plus the optical depth needed
for the neutrino to reach that neighbor
($\bar{\kappa}_{\nu_{i}}{(g_{{k}j}dx^{{k}}dx^{j})}^{1/2}$):

$\tau_{\nu_{i}}=\min\big{(}\tau_{{\nu_{i}},\rm{neighbor}}+\bar{\kappa}_{\nu_{i}}{(g_{{k}j}dx^{{k}}dx^{j})}^{1/2}\big{)}$
(52)

where $\tau_{{\nu_{i}},\rm{neighbor}}$ is the optical depth of the neighboring
cell, $\bar{\kappa}_{\nu_{i}}$ is the average opacity between the cell and its
neighbor, and ${(g_{{k}j}dx^{{k}}dx^{j})}^{1/2}$ is the distance to the
neighboring cell calculated by taking the average value of $g_{kj}$ between
the local and neighboring cells. We minimize over all neighbors.

This essentially traces the path of least resistance of the neutrino to a
neighbor. We update the entire grid, and perform the next iteration, where
again, we minimize over all the adjacent neighbors. The next iteration will
show the path to the neighbor two cells away. As we do more iterations, we
trace the path of least resistance that the neutrinos will take out of the
star. This will lead us to the final optical depth. During the first timestep,
we initialize the optical depth by by performing $20N_{\rm max}$ iterations,
where $N_{\rm max}$ is the maximum number of cells in each direction,
independent of resolution. This is done to trace a path to the edge of the
domain initially. After the initial calculation, which has a fixed number of
iterations, we continue to do iterations to obtain the final optical depth,
however we impose a convergence criterion in order to minimize the number of
iterations. In order to converge, we set conditions on the difference between
iteration ${k}-1$ and ${k}$:

$R_{\rm
change,\tau}({k})\equiv\frac{|\sum\tau_{{k}-1}-\sum\tau_{{k}}|}{\sum\tau_{{k}-1}}<\epsilon_{1}$
(53)

or

$\frac{|R_{\rm change,\tau}({k}-1)-R_{\rm change,\tau}({k})|}{R_{\rm
change,\tau}({k}-1)}<\epsilon_{2}$ (54)

where $\sum\tau_{k}$ is the sum of all the optical depths in the grid at
iteration ${k}$, and $\epsilon_{1}$ and $\epsilon_{2}$ are parameters that we
choose to be $\epsilon_{1}=10^{-4}$ and $\epsilon_{2}=10^{-3}$, respectively.
Only a few iterations are needed for convergence after the initial guess.

##  3 Validation tests for the tabulated EOS

In this section, we describe the tests performed to validate the implemented
EOS tables.

###  3.1 Testing the conserved to primitive variables routine

In order to validate the routines that transform the conserved variables into
primitive variables with tabulated EOS, we created primitive variables out of
a grid of density and temperature values within the EOS table. The magnetic
field was set randomly to be either aligned or anti-aligned with the velocity
vector. The magnitude of the magnetic field was set to be such that:
$b^{2}/2=\left(\frac{P_{\rm mag}}{P_{\rm gas}}\right)P_{\rm gas}$, where
$(P_{\rm mag}/P_{\rm gas})$ is set as a parameter, $P_{\rm mag}$ is the
magnetic pressure, and $P_{\rm gas}$ is the gas pressure. We then obtained a
set of conserved variables based on these primitives. The true primitives were
then varied by randomly adding or subtracting a $5\%$ perturbation to each
primitive. This test is based on Siegel et al. (2018).

We then used these primitives as initial guesses for the various routines that
transform the conserved variables to primitive variables and compared the
resultant solution to the original.

We show the error we obtained for all primitive variables in Figure 1. It can
be seen that the recovery error is low. Additionally, the figure shows that
the 3d method is less robust, but more accurate, which is the reason it is set
as the primary routine. The different 2d methods, and the ’Palenenzuela’
routine are more robust, but less accurate (and slower) than the 3d method, so
they serve better as backup routines.

Figure 2:  Top panel: Evolution of a torus in hydrostatic equilibrium with a
tabulated EOS and no neutrino leakage scheme.. We show the meridional (Left)
and equatorial (Right) cut. The initial conditions are set as described in
section 3.2. Here $x_{1}$, $x_{2}$, $x_{3}$ correspond to the coordinates $x$,
$z$, $y$ respectively. Bottom panel: Density as a function of radius for
different times in the equator. Figure 3: Comparison of the semi-analytical
solution (dotted line) with the simulation (solid line) for the evolution of
the temperature and electron fraction of an isotropic, optically thin gas with
constant density.  Figure 4: Comparison of the analytical solution (dotted
line) of the optical depth with simulations (solid line) for different
resolutions. The labels indicate the number of cells in each direction. The
Top panel is the anti-neutrino optical depth, and the Bottom panel is the
neutrino optical depth.  Figure 5: Comparison of the analytical solution
(dotted line) of the optical depth to neutrinos with simulations (solid line).

###  3.2 Torus in hydrostatic equilibrium

To test the EOS implementation, we simulated a non-magnetized torus that is in
hydrostatic equilibrium with no leakage scheme, following Fishbone & Moncrief
(1976).

Figure 2 shows the 3d hydrodynamical evolution of a torus constructed to be in
hydrostatic equilibrium with a tabulated EOS without neutrino cooling. There
are perturbations particularly near the BH due to accretion onto the BH, but
the density is low in those regions. As can be seen from the figure, the torus
remains in hydrostatic equilibrium throughout the simulation.

####  3.2.1 Initial conditions inside the torus

The specific enthalpy inside the torus is implemented via Equation (3.6) of
Fishbone & Moncrief (1976), but adding $\ln{h_{\rm min}}$ to the integration
constant (see Section 3.2.2). By construction, the torus is in hydrostatic
equilibrium with the ambient atmosphere. We also set the torus to be
isentropic, and have uniform electron fraction. Given a specific entropy
$s_{\rm disk}$, a specific enthalpy given by Fishbone & Moncrief (1976), and
an electron fraction $Y_{e\rm,disk})$, the temperature and density of the disk
are found by solving the following equations:

$s_{\rm disk}=s(\rho,T,Y_{e})$ (55) $h=h(\rho,T,Y_{e})$ (56)

where $s$ is the specific entropy.

####  3.2.2 Atmosphere

In the classical torus, the boundaries of the torus are defined where $h=1$.
In the tabulated EOS, though, negative internal energy densities are allowed
since the internal energy per nucleon is measured relative to the free neutron
rest mass energy. In this case, the minimum specific enthalpy is not
restricted to $1$, but rather it can be $1>h_{\rm min}>0$, where $h_{\rm min}$
is the specific enthalpy from the table given the atmospheric density, and the
disk’s electron fraction and specific entropy. Thus, we set the torus boundary
to be where $h=h_{\rm min}$. For the background atmosphere, we set the minimum
atmospheric density $\rho_{\rm atm}$ as a parameter. Then we find the minimum
specific enthalpy by doing a table inversion and finding $h_{\rm
min}=h(\rho_{\rm atm},s_{\rm disk},Y_{e\rm,disk})$. We also find the
atmospheric temperature by doing a table inversion $T_{\rm atm}=T(\rho_{\rm
atm},s_{\rm disk},Y_{e\rm,disk})$.

The density in the background is set to:

$\max\bigg{(}\rho_{\rm atm},\frac{\rho_{0}}{r^{2}}\bigg{)}$ (57)

Where we set $\rho_{0}$ as a parameter as well. The background atmosphere
temperature is set to:

$\max\bigg{(}T_{\rm atm},\frac{T_{0}}{r}\bigg{)}\ ,$ (58)

where $T_{0}$ is a parameter. The power-law dependence is set to provide the
background atmosphere with more pressure support so that it does not rapidly
accrete onto the BH. This ultimately helps with robustness near the BH, as the
low density and low temperature zones with high velocity are where the
conserved to primitive routines tend to fail. We note that in the region where
there is a power-law dependence, the specific enthalpy is not a constant,
whereas once the background atmosphere is set to be constant, everything is
thermodynamically consistent because it was constructed with the tabulated EOS
tables.We set the electron fraction of the atmosphere to a constant value
found by assuming $\beta$-equilibrium (where the neutrino chemical potential
is zero) at $T_{\rm atm}$ and $\rho_{\rm atm}$.

The units are normalized so that the maximum density in the torus is set to
$\rho_{\rm max}=1$ in code units, which in this case corresponds to $\rho_{\rm
max}=5.4\times 10^{8}{\rm g/cm^{3}}$ in cgs units. In the simulation we
performed, the torus has a constant electron fraction of $Y_{e}=0.1$ and a
specific entropy of $10{\rm k_{B}/baryon}$, where ${\rm K_{B}}$ is Boltzmann’s
constant. The background atmosphere is characterized by $\rho_{\rm
atm}=6000{\rm g/cm^{3}}$, $\rho_{0}=3\times 10^{5}{\rm g/cm^{3}}$,
$T_{0}=0.4$MeV. We used the SLy4 table with NSE from Schneider et al. (2017),
and with that table the minimum specific enthalpy for our parameters is set to
$h_{\rm min}=0.9974$ (in code units), and $T_{\rm atm}=0.0053$MeV. The
electron fraction in the atmosphere, given by $\beta$-equilibrium, is set to
$Y_{e,{\rm atm}}=0.45$. The boundary conditions are outflow in the outer
radial boundary, reflective in the angular coordinate $\theta$, and periodic
in the angular coordinate $\phi$. The metric is Kerr-Schild in spherical
coordinates for a non-spinning BH.

##  4 Validation tests for the leakage scheme

In this subsection, we describe how we tested the leakage scheme in the
optically thin regime and for finite optical depth.

###  4.1 Testing the optically thin regime

Following Miller et al. (2019a), we tested the leakage scheme in an optically
thin regime by considering an isotropic gas of constant density and
temperature such that the gas is optically thin to neutrinos. We tested both
reactions in the charged $\beta-$process separately where we included only
either the neutrinos or the anti-neutrinos.

In this case, the GRMHD equations reduce to:

$\partial_{t}{{T^{t}}_{t}}=\mathcal{Q}.\ ,$ (59)
$\partial_{t}{Y_{e}}=\mathcal{R}/\rho\ ,$ (60)

where $\mathcal{R}$ and $\mathcal{Q}$ are the emission/absorption and energy
loss rates due to neutrinos or anti-neutrinos of the reactions in
$\beta-$process separately. The rates need to be calculated semi-analytically,
since they depend on interpolated quantities, such as the degeneracy
parameters. We can then solve the equations semi-analytically with a set of
initial conditions and compare to simulations. We chose the initial density
and temperature such that the medium is optically thin to neutrinos and anti-
neutrinos.

For the initial conditions, we used an initial density of $617714{\rm
g/cm^{3}}$ and temperature of $1$MeV, chosen so that the medium is optically
thin to neutrinos and anti-neutrinos. We used $Y_{e,0}=0.5$, $Y_{e,0}=0.005$
for the electron neutrino and anti-neutrino tests respectively. We used
$2\times 2\times 1$ number of cells in each direction using a Cartesian grid
with Minkowski metric.

In Figure 3 we show the comparison between the semi-analytical solution and
the simulation for the $\beta$-process both for neutrinos and anti-neutrinos.
We compare the change in the electron fraction due to the absorption/emission
rate, and the change in temperature due to the heating/cooling rate. As can be
seen from the figure, HARM3D+NUC is able to recreate the semi-analytical
solution.

###  4.2 Testing the optically thick regime

####  4.2.1 Constant density circular disk

In order to test the optical depth calculation, we simulated a circular disk
with uniform density and temperature embedded in an optically thin medium of
constant density and temperature. The advantage of this scenario is that we
can calculate the opacity inside the circle and then calculate the optical
depth analytically. This way we can compare to the simulation. The simulations
were performed in 2d, and the domain is $2r_{\rm g}$, where $r_{\rm
g}=GM/c^{2}$ is the gravitational radius. We used a Minkowski metric with
spherical coordinates. There are outflow conditions on the radial boundaries.
The optical depth in the outer radial boundary was set to zero so that the
neutrinos and anti-neutrinos could escape the domain. We simulated an
optically thick circular disk that has a constant density of $9.8\times
10^{13}{\rm g/cm^{3}}$, an electron fraction of $0.1$ and a temperature of
$8$MeV embedded in an optically thin medium, with a density of $6\times
10^{7}{\rm g/cm^{3}}$, an electron fraction of $0.5$ and a temperature of
$0.01$MeV. Figure 4 shows the optical depth for both the electron neutrino and
anti-neutrino for different resolutions. As can be seen from the figure, the
initial guess for the optical depth is accurate and the convergence to the
solution does not change with resolution. At smaller optical depths, the
optical depth is slightly overestimated at lower resolutions, but as the
optical depth increases, the solution doesn’t depend noticeably on resolution.

####  4.2.2 Stripes

We can also test the optical depth algorithm by simulating stripes of high
density material with low density material in between. In this scenario, it is
expected that a neutrino created in the region with high optical depth
material will travel to the region with low optical depth and stream freely
from the surface. For the simulation, we used $4096\times 96\times 1$ cells.
The simulations were performed in 2d, and the domain is $1r_{\rm g}$ large in
radial extent, where $r_{\rm g}=GM/c^{2}$ is the gravitational radius. We used
a Minkowski metric with spherical coordinates and outflow conditions at the
radial boundaries. The optical depth at the outer radial boundary was set to
zero so that the neutrinos and anti-neutrinos could escape the domain. We
simulated three stripes of material with high optical depth: $\rho=9.8\times
10^{13}{\rm g/cm^{3}}$, $Y_{e}=0.1$, $T=8\mathrm{MeV}$. In between the
stripes, the optically thin gas was initialized to $\rho=6\times
10^{7}\mathrm{g/cm^{3}}$, $Y_{e}=0.5$, and $T=0.01\mathrm{MeV}$. The high
opacity stripes start at $r=0r_{\rm g}$, and have a width of $r=0.1r_{\rm g}$.
The next stripes are located in $r=0.2r_{\rm g}$ and $r=0.4r_{\rm g}$.

We show the results from this setup in Figure 5, where we compare the results
from the simulation with the analytical estimate (length units are in $r_{\rm
g}$):

$\tau_{\rm analytical}=\left\\{\begin{array}[]{lllll}\int_{0}^{0.2}\kappa
dr&r\leq 0.2\\\ \int_{0.2}^{0.25}\kappa dr&0.2\leq r\leq 0.25\\\
\int_{0.25}^{0.35}\kappa dr&0.25\leq r\leq 0.35\\\ \int_{0.4}^{0.45}\kappa
dr&0.4\leq r\leq 0.45\\\ \int_{0.45}^{0.55}\kappa dr&0.45\leq r\leq
0.55\end{array}\right.$ (61)

Figure 6: Shown is a meridional cut of the MRI quality factors $Q_{\rm mri,2}$
and $Q_{\rm mri,3}$ at 114ms, where the subscripts for $Q_{\rm mri,2}$,
$Q_{\rm mri,3}$ correspond to the coordinates $\theta,\phi$ respectively.
Figure 7:  Top and Middle panel: Shell-integrated mass-weighted quality
factors as a function of radius, averaged over different epochs of time.
Bottom panel: Mass-weighted quality factors integrated over angles and radii
that are less than 150km: $\int^{150\rm km}_{0\rm km}\int\int Q_{\rm
mri}\rho\sqrt{-g}drd\phi d\theta/\int^{150\rm km}_{0\rm
km}\int\int\sqrt{-g}\rho drd\phi d\theta$.  Figure 8: Density of a magnetized
torus including the impact of neutrinos at 114ms. Shown is an equatorial cut
(top panel) and a meridional cut (bottom panel).  Figure 9: Electron fraction
of a magnetized torus including the impact of neutrinos at 114ms. Shown is an
equatorial cut (top panel) and a meridional cut (bottom panel). Figure 10:
Zoomed in version of the density of a magnetized torus including the impact of
neutrinos at 114ms. Shown is an equatorial cut (top panel) and a meridional
cut (bottom panel).  Figure 11: Zoomed in version of the electron fraction of
a magnetized torus including the impact of neutrinos at 114ms. Shown is an
equatorial cut (top panel) and a meridional cut (bottom panel).

##  5 Magnetized disk

In this section, we apply our new code HARM3D+NUC to a magnetized torus in 3d
that approximates a post-merger disk. We use both the tabulated EOS and the
leakage scheme in this test.

###  5.1 Initial conditions

The initial conditions inside the torus follow a similar setup to that of
section 3.2.1, but with the addition of a poloidal magnetic field. In order to
start with a magnetic field devoid of magnetic monopoles, we first set the
vector potential to a prescribed distribution and calculate its curl using a
finite difference operator compatible with our constrained transport method
(see Zilhão & Noble, 2014, for further details). Our poloidal magnetic field
distribution results from a vector potential with only one non-zero component:

$A_{\phi}=\max\left(\overline{\rho}/\rho_{\rm max}-\rho_{0,\rm mag}\
,0\right)$ (62)

where $\overline{\rho}$ is the average density at that position, and
$\rho_{\rm max}=1.66\times 10^{11}{\rm g/cm^{3}}$ is the maximum density of
the torus. We set $\rho_{0,\rm mag}=0.2$ in code units, which corresponds to
$\rho_{0,\rm mag}=3.33\times 10^{10}{\rm g/cm^{3}}$. Then we build the
magnetic field with the vector potential and normalize its magnitude such that
the ratio of the integrated gas pressure to integrated magnetic pressure is
100. Inside the disk, the matter is set to be neutron rich, $Y_{e}=0.1$. The
treatment of the atmosphere is the same as in section 3.2.2, except the
density scales as $r^{-3/2}$. In the atmosphere, the electron fraction is set
to its value in $\beta$-equilibrium, where the chemical potential of the
neutrinos is set to zero. We show the parameters used in Table 1.

The simulations were performed in 3d on a grid designed to focus more cells
about the equator and towards the black hole horizon. We use the same grid as
defined in Noble et al. (2010) but with different parameters. The azimuthal
grid spacing is uniform. The logarithmic radial grid is such that $\Delta r/r$
is fixed and the $i^{\mathrm{th}}$ cell center is located at:

$r_{i}=r_{\mathrm{min}}\exp{\left[\left(i+1/2\right)\log_{\mathrm{10}}{(r_{\mathrm{max}}/r_{\mathrm{min}})}/N_{r}\right]}\
,$ (63)

with $r_{\mathrm{min}}=1.303r_{\rm g}$, $r_{\mathrm{max}}=2000r_{\rm g}$, and
$i\in\left[0,N_{r}-1\right]$. The $\theta$ grid uses a high-order polynomial
function to provide a nearly uniform grid spacing spacing near the equator:

$\theta_{j}=\frac{\pi}{2}\left[1+\left(1-\xi\right)\left(2x^{(2)}_{j}-1\right)+\left(\xi-\frac{2\theta_{c}}{\pi}\right)\left(2x^{(2)}_{j}-1\right)^{n}\right]\
,$ (64)

where $\xi$ is a parameter controlling the severity of the focusing, $n$ is
the order of polynomial used in the transformation, $\theta_{c}$ is the
opening angle of the polar regions we excise,
$x^{(2)}_{j}\equiv\left(j+1/2\right)/N_{\theta}$, and
$j\in\left[0,N_{\theta}-1\right]$ . In our run, we used $\theta_{c}=\pi
10^{-14}$, $\xi=0.65$, and $n=7$. The number of cells per dimension used was
$N_{r}\times N_{\theta}\times N_{\phi}=1024\times 160\times 256$.

Parameter | Value
---|---
Disk radius of maximum pressure | $9r_{\rm g}$
Disk inner radius | $4r_{\rm g}$
Mass of disk | $0.03M_{\odot}$
$Y_{e}$ in the disk | 0.1
Specific entropy in the disk |  $7\ k_{\rm b}$/baryon
$\beta$ | 100
BH spin | 0.9375
BH mass | $3M_{\odot}$
Specific enthalpy at boundary | 0.9977 [code units]
Temperature at radius of maximum pressure | 4.4 MeV
Table 1: Parameters used in the simulation.

###  5.2 Scaling tests

We performed scaling tests for this run for 3 different number of processors:
5120, 2560 and 1280 processors. For this setup, the number of time steps in
the code per second per processor were: 0.000723, 0.000781, 0.000868
respectively. The difference between 5120 and 1280 processors is around
$17\%$. If we don’t include the neutrino leakage scheme but include only a
tabulated EOS, for 2560 processors, the number of steps per second per
processor is 0.001328, which makes the leakage $58\%$ slower than only
considering the tabulated EOS.

###  5.3 Magnetic turbulence

In order to confirm that we are adequately resolving magnetic turbulence, we
display in Figure 6 the number of grid cells per wavelength of the fastest
growing mode of the magneto-rotational instability (MRI), defined as (Noble et
al., 2010; Hawley et al., 2011; Sorathia et al., 2012; Hawley et al., 2013):

$Q_{\rm mri,x}=\frac{\lambda_{x,\rm mri}}{\Delta_{x}}$ (65)

where $x=\theta,\phi$, $\Delta_{x}$ is the cell size, and the wavelength of
the fastest MRI growing mode is:

$\lambda_{x,\rm mri}=\frac{2\pi}{\Omega}\frac{|b^{x}|}{\sqrt{\rho h+b^{2}}}\
.$ (66)

As can be seen in Fig 7, our grid satisfies the criterion of Sano et al.
(2004) everywhere except for later times within $r<50$km. While our results
fail to meet criteria for asymptotic MRI convergence set forth in Hawley et
al. (2011), our disk does satisfy $Q_{\rm mri,3}>10$ everywhere, and $Q_{\rm
mri,2}>6$ for $r\gtrsim 50$km for most of the run. One reason why our
simulation may not reach larger $Q_{\rm mri}$ values is because we used the
same random perturbations across all MPI processes in the initial conditions.
Because we used 16 subdomains in the azimuthal dimension, this means that the
simulation is periodic over $\Delta\phi=\pi/8$, and the azimuthal modes with
$m<8$ start off significantly weaker as they are seeded with perturbations at
only the round-off error level.

Figure 12: Shown is the geometrical thickness ($H/r$) of the disk, as a
function of radius. The thickness is averaged between the indicated time in
the legend. Figure 13: Top panel: Mass accretion rate onto the innermost
stable circular orbit (ISCO) of the BH as a function of time. Bottom panel:
Average mass accretion rate as a function of radius. We averaged the mass
accretion rate between the times indicated in the legend. Figure 14:
Luminosity due to the different neutrino species as a function of time.

###  5.4 Impact of neutrinos and EOS

Magnetic stresses will transport angular momentum in the disk, heating the
gas, which will produce a high velocity outflow (Fernández & Metzger, 2013a;
Siegel & Metzger, 2018). This outflow will be affected by the addition of
neutrinos formed through weak reactions. In the midplane, neutrinos will carry
significant amounts of energy, which will cool and make the disk geometrically
thinner. Another outflow is also expected to occur in the outer regions of the
disk due to the release in nuclear binding energy when there is recombination
of free nucleons into $\alpha-$particles, which produces enthalpy and unbinds
material (Lee et al., 2009; Fernández & Metzger, 2013a). In this subsection we
show the impact of both the emission of neutrinos and the recombination of
free nucleons.

In Figures 8 and 11 we display the outflows that results from our simulations
of a neutrino-cooled magnetized disk at 114 ms. In Figures 9 and 10 we plot
the electron fraction and the density, respectively, at $t={114}\mathrm{ms}$.
Neutrino cooling is expected to happen in diffusion timescales, which is on
the order of milliseconds, much shorter than our evolution timescale. The
inner regions of the disk are very neutron rich, confirming the self-
regulating phase found in Siegel & Metzger (2017, 2018). In this phase, there
is a balance between the neutrino cooling and the heating driven by MHD that
self-regulates the electron degeneracy parameter, and the final state is a
neutron rich disk (Siegel & Metzger, 2017, 2018). We note that, although this
new code does not include neutrino absorption in the ejecta, the absorption
should modify the electron fraction in the outflow (Just et al., 2021).

In the top panel of Figure 13 we show the mass accretion rate through the
innermost stable circular orbit (ISCO) as a function of time, and show the
accretion rate as a function of radius in the bottom panel. The outflow can be
clearly seen as a negative mass accretion rate at larger radii, as well as a
settling of the mass accretion rate as time passes.

In Figure 12, we plot the geometrical thickness of the disk, or $H/r$. We
estimated this thickness using the scale height $H$ following Noble et al.
(2012):

$H=\frac{\langle\rho\sqrt{g_{\theta\theta}}\left|\theta-\pi/2\right|\rangle}{\langle\rho\rangle}$
(67)

where $\langle X\rangle$ is the average of the quantity $X$ over a spherical
shell:

$\langle X\rangle=\frac{\int X\sqrt{-g}d\theta d\phi}{\int\sqrt{-g}d\theta
d\phi}\ .$ (68)

In the deepest regions of the disk, the heating due to MHD turbulence helps
create neutrinos/anti-neutrinos, which escape, remove energy, and
geometrically thin the disk. Recombination of free nucleons into
$\alpha$-particles releases binding energy, effectively increasing the
enthalpy and unbinds material. The effect of recombination is less severe than
the geometrically thinning due to neutrino/anti-neutrino losses. This
transition can be seen at around $150\mathrm{km}$.

We may obtain the amount of energy radiated by each species of neutrino and
anti-neutrino as was done in Siegel & Metzger (2018):

$L_{\nu_{i}}=\int\alpha\gamma\mathcal{Q}_{\nu_{i}}^{\rm eff}\sqrt{-g}d^{3}x\
.$ (69)

In Figure 14, we show the luminosity for each species. It can be seen that the
electron neutrino (and anti-neutrino) dominate the emission over all of the
other species of neutrino. The luminosity roughly follows the mass accretion
rate as seen in Figure 13, as heating from the magnetic stresses ignite the
creation of neutrinos/anti-neutrinos. This suggests the radiative efficiency
of neutrino/anti-neutrinos emission remains relatively steady.

Our initial conditions are similar (although not identical) to the initial
conditions in Siegel & Metzger (2018). They performed 3d simulations of a
post-merger accretion disk with a relatively higher specific entropy and lower
spin than this simulation. They used Cartesian coordinates, a Helmholtz EOS
for relatively low densities, and a neutrino leakage scheme. They evolved the
disk for longer times (380 ms). Even though we use a different EOS (Sly4), the
disk thickness is qualitatively similar. At the inner regions of the disk,
neutrino cooling dominates, whereas at outer regions (at radius higher than
around 100km), recombination is responsible for making the disk geometrically
thicker. The neutrino/anti-neutrino luminosities are comparable, Siegel &
Metzger (2018) has a higher luminosity, but that could be attributed to the
difference in the initial disk specific entropy.

As the outflow expands, it will cool, and heavy elements will be created via
the $r-$process. We will explore this nucleosynthesis in a future paper.

##  6 Summary

GRMHD simulations of post-merger accretion disks have advanced over the last
few years with better treatment of neutrinos and a more realistic EOS. In this
paper we present the addition of a neutrino leakage scheme and a tabulated EOS
into the computationally efficient, versatile GRMHD code HARM3D. This new
addition to HARM3D, called HARM3D+NUC, has the potential to be used in a range
of simulations where neutrinos are present. In the paper, we use the new code
HARM3D+NUC to simulate an accretion disk resembling the post-merger phase of a
binary neutron star, though other applications include collapsars (e.g.,
Siegel et al., 2019; Miller et al., 2020).

The paper shows how we implemented the tabulated EOS in the conserved variable
to primitive variable routines, and the different methods we implemented and
tested for performing this inversion. We show that using the 3d primary
recovery method is the most accurate and efficient, but least robust, choice
which is why we also employ several 2d and 1d backup routines. The leakage
scheme is implemented by adding the neutrino/anti-neutrino heating/cooling and
emission/absorption terms as source terms in the equations of motion. We
describe in detail an approach to obtain the optical depth locally and how we
can use a convergence criterion to get the optical depth after a few
iterations once the initial guess is made.

We show several tests for our new code. The tabulated EOS is tested by
determining the relative error between original primitive variables and the
recovered primitive variables. We also test the EOS by performing a simulation
of a torus in hydrostatic equilibrium, showing that it stays in hydrostatic
equilibrium throughout the entire simulation. We test the neutrino leakage
scheme in the optically thin regime by investigating the $\beta-$process in a
constant density gas. We test the optical depth algorithm in a constant
density circular disk and a stripes setup.

With our new machinery, we simulate a magnetized high-density torus, which
serves as an approximation to the accretion flow after the merger of two
neutron stars. Magnetic stresses transport angular momentum from the disk,
driving a high velocity outflow. The outflow is affected by both the addition
of neutrinos and the nuclear binding energy released from the recombination of
nucleons to $\alpha-$particles, which acts to geometrically thicken the disk.
Neutrinos will alter the electron fraction of the ejecta especially in the
inner regions of the disk, whereas the recombination of nucleons is more
prominent in the outer regions of the disk. This highlights the importance of
modeling the accretion disk including neutrinos and an EOS that considers this
extra unbinding of material due to recombination.

We plan to use the new code to do long-term evolutions of binary neutron star
mergers starting from before the neutron stars merge to the evolution of the
outflow. Heavy elements should be created via the $r-$process in this outflow
as it expands and cools. We plan to use different codes and methods to treat
the initial data, pre-merger/merger, and post-merger phases. The initial data
for the neutron stars will be constructed using a modified version of LORENE
(Gourgoulhon et al., 2016) we have developed. Binaries will be evolved until
they merge and eventually form a black hole surrounded by an accretion disk
using two GRMHD codes: IllinoisGRMHD (Etienne et al., 2015), and Spritz
(Cipolletta et al., 2020). After the remnant has collapsed to a BH and the
numerical metric has stabilized, we will interpolate the MHD primitives and
numerical metric into the grid of HARM3D+NUC (López Armengol et al. in prep).
After doing the appropriate tensorial transformations from the Cartesian base
to the coordinate base of HARM3D+NUC, we will continue the post-merger
evolution with HARM3D+NUC.

We thank the anonymous referee, R. Foley, B. Villaseñor, D. Radice, T. Piran,
V. Mewes, D. Siegel, S. Rosswog, A. Janiuk, J. Miller, J. Dolence, M. C.
Miller, P. Mösta, N. M. Lloyd-Ronning, A. Batta, G. Koenigsberger, D. Kasen
for valuable conversations. A.M-B and E. R-R are supported by the Heising-
Simons Foundation, the Danish National Research Foundation (DNRF132), NSF
(AST-1911206 and AST-1852393). A.M-B is supported by the UCMEXUS-CONACYT
Doctoral Fellowship. B.J.K was supported by the NASA Goddard Center for
Research and Exploration in Space Science and Technology (CRESST) II
Cooperative Agreement under award number 80GSFC17M0002. R.OS was supported by
NSF PHY-2012057 and AST-1909534. This work was made possible by the NASA TCAN
award TCAN-80NSSC18K1488. Computational resources were provided by the NCSA’s
Blue Waters sustained-petascale computing NSF projects OAC-1811228 and
OAC-1516125, by the TACC’s Frontera NSF projects PHY20010 and AST20021, and
the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315.

## References

  * Abbott et al. (2017a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, Physical Review Letters, 119, 161101
  * Abbott et al. (2017b) —. 2017b, ApJ, 848, L12
  * Baiotti et al. (2008) Baiotti, L., Giacomazzo, B., & Rezzolla, L. 2008, Phys. Rev. D, 78, 084033
  * Balbus & Hawley (1998) Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1
  * Barnes & Kasen (2013) Barnes, J., & Kasen, D. 2013, ApJ, 775, 18
  * Barnes et al. (2016) Barnes, J., Kasen, D., Wu, M.-R., & Martínez-Pinedo, G. 2016, ApJ, 829, 110
  * Cerdá-Durán et al. (2008) Cerdá-Durán, P., Font, J. A., Antón, L., & Müller, E. 2008, A&A, 492, 937
  * Chen & Beloborodov (2007) Chen, W.-X., & Beloborodov, A. M. 2007, ApJ, 657, 383
  * Chevalier (1989) Chevalier, R. A. 1989, ApJ, 346, 847
  * Ciolfi et al. (2017) Ciolfi, R., Kastaun, W., Giacomazzo, B., et al. 2017, Phys. Rev. D, 95, 063016
  * Cipolletta et al. (2020) Cipolletta, F., Kalinani, J. V., Giacomazzo, B., & Ciolfi, R. 2020, Classical and Quantum Gravity, 37, 135010
  * Coulter et al. (2017) Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Science, 358, 1556
  * Cowperthwaite et al. (2017) Cowperthwaite, P. S., Berger, E., Villar, V. A., et al. 2017, ApJ, 848, L17
  * De & Siegel (2020) De, S., & Siegel, D. 2020, arXiv e-prints, arXiv:2011.07176
  * 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
  * Dionysopoulou et al. (2015) Dionysopoulou, K., Alic, D., & Rezzolla, L. 2015, Phys. Rev. D, 92, 084064
  * Drout et al. (2017) Drout, M. R., Piro, A. L., Shappee, B. J., et al. 2017, Science, 358, 1570
  * Eichler et al. (1989) Eichler, D., Livio, M., Piran, T., & Schramm, D. N. 1989, Nature, 340, 126
  * Etienne et al. (2015) Etienne, Z. B., Paschalidis, V., Haas, R., Mösta, P., & Shapiro, S. L. 2015, Classical and Quantum Gravity, 32, 175009
  * Fernández et al. (2017) Fernández, R., Foucart, F., Kasen, D., et al. 2017, Classical and Quantum Gravity, 34, 154001
  * Fernández & Metzger (2013a) Fernández, R., & Metzger, B. D. 2013a, MNRAS, 435, 502
  * Fernández & Metzger (2013b) —. 2013b, ApJ, 763, 108
  * Fernández et al. (2015) Fernández, R., Quataert, E., Schwab, J., Kasen, D., & Rosswog, S. 2015, MNRAS, 449, 390
  * Fishbone & Moncrief (1976) Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962
  * Foucart et al. (2018) Foucart, F., Duez, M. D., Kidder, L. E., et al. 2018, Phys. Rev. D, 98, 063007
  * Foucart et al. (2016) Foucart, F., O’Connor, E., Roberts, L., et al. 2016, Phys. Rev. D, 94, 123016
  * Foucart et al. (2015) —. 2015, Phys. Rev. D, 91, 124021
  * Freiburghaus et al. (1999) Freiburghaus, C., Rosswog, S., & Thielemann, F. K. 1999, ApJ, 525, L121
  * Fujibayashi et al. (2017) Fujibayashi, S., Sekiguchi, Y., Kiuchi, K., & Shibata, M. 2017, ApJ, 846, 114
  * Galeazzi et al. (2013) Galeazzi, F., Kastaun, W., Rezzolla, L., & Font, J. A. 2013, Phys. Rev. D, 88, 064009
  * Gammie et al. (2003) Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444
  * Gourgoulhon et al. (2016) Gourgoulhon, E., Grandclément, P., Marck, J.-A., Novak, J., & Taniguchi, K. 2016, LORENE: Spectral methods differential equations solver, , , ascl:1608.018
  * Grossman et al. (2014) Grossman, D., Korobkin, O., Rosswog, S., & Piran, T. 2014, MNRAS, 439, 757
  * Hawley et al. (2011) Hawley, J. F., Guan, X., & Krolik, J. H. 2011, ApJ, 738, 84
  * Hawley et al. (2013) Hawley, J. F., Richers, S. A., Guan, X., & Krolik, J. H. 2013, ApJ, 772, 102
  * Hossein Nouri et al. (2018) Hossein Nouri, F., Duez, M. D., Foucart, F., et al. 2018, Phys. Rev. D, 97, 083014
  * Houck & Chevalier (1991) Houck, J. C., & Chevalier, R. A. 1991, ApJ, 376, 234
  * Janiuk (2014) Janiuk, A. 2014, A&A, 568, A105
  * Janiuk (2019) —. 2019, ApJ, 882, 163
  * Janiuk et al. (2013) Janiuk, A., Mioduszewski, P., & Moscibrodzka, M. 2013, ApJ, 776, 105
  * Just et al. (2021) Just, O., Goriely, S., Janka, H.-T., Nagataki, S., & Bauswein, A. 2021, arXiv e-prints, arXiv:2102.08387
  * Just et al. (2015) Just, O., Obergaulinger, M., & Janka, H. T. 2015, MNRAS, 453, 3386
  * Kasen et al. (2013) Kasen, D., Badnell, N. R., & Barnes, J. 2013, ApJ, 774, 25
  * Kasen & Barnes (2019) Kasen, D., & Barnes, J. 2019, ApJ, 876, 128
  * Kasen et al. (2015) Kasen, D., Fernández, R., & Metzger, B. D. 2015, MNRAS, 450, 1777
  * Kasen et al. (2017) Kasen, D., Metzger, B., Barnes, J., Quataert, E., & Ramirez-Ruiz, E. 2017, Nature, 551, 80
  * Kasliwal et al. (2019) Kasliwal, M. M., Kasen, D., Lau, R. M., et al. 2019, MNRAS, L14
  * Kilpatrick et al. (2017) Kilpatrick, C. D., Foley, R. J., Kasen, D., et al. 2017, Science, 358, 1583
  * Kiuchi et al. (2018) Kiuchi, K., Kyutoku, K., Sekiguchi, Y., & Shibata, M. 2018, Phys. Rev. D, 97, 124039
  * Kiuchi et al. (2014) Kiuchi, K., Kyutoku, K., Sekiguchi, Y., Shibata, M., & Wada, T. 2014, Phys. Rev. D, 90, 041502
  * Kiuchi et al. (2015) Kiuchi, K., Sekiguchi, Y., Kyutoku, K., et al. 2015, Phys. Rev. D, 92, 064034
  * Kohri & Mineshige (2002) Kohri, K., & Mineshige, S. 2002, ApJ, 577, 311
  * Kulkarni (2005) Kulkarni, S. R. 2005, arXiv e-prints, astro
  * Lee & Ramirez-Ruiz (2007) Lee, W. H., & Ramirez-Ruiz, E. 2007, New Journal of Physics, 9, 17
  * Lee et al. (2009) Lee, W. H., Ramirez-Ruiz, E., & López-Cámara, D. 2009, ApJ, 699, L93
  * Lee et al. (2004) Lee, W. H., Ramirez-Ruiz, E., & Page, D. 2004, ApJ, 608, L5
  * Lee et al. (2005) —. 2005, ApJ, 632, 421
  * 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 et al. (2017) Lippuner, J., Fernández, R., Roberts, L. F., et al. 2017, MNRAS, 472, 904
  * Lippuner & Roberts (2015) Lippuner, J., & Roberts, L. F. 2015, ApJ, 815, 82
  * Margalit & Metzger (2017) Margalit, B., & Metzger, B. D. 2017, ApJ, 850, L19
  * Martin et al. (2015) Martin, D., Perego, A., Arcones, A., et al. 2015, ApJ, 813, 2
  * Metzger (2019) Metzger, B. D. 2019, Living Reviews in Relativity, 23, 1
  * Metzger et al. (2010) Metzger, B. D., Arcones, A., Quataert, E., & Martínez-Pinedo, G. 2010, MNRAS, 402, 2771
  * Metzger & Fernández (2014) Metzger, B. D., & Fernández, R. 2014, MNRAS, 441, 3444
  * Metzger et al. (2008) Metzger, B. D., Piro, A. L., & Quataert, E. 2008, MNRAS, 390, 781
  * Miller et al. (2019a) Miller, J. M., Ryan, B. R., & Dolence, J. C. 2019a, ApJS, 241, 30
  * Miller et al. (2020) Miller, J. M., Sprouse, T. M., Fryer, C. L., et al. 2020, ApJ, 902, 66
  * Miller et al. (2019b) Miller, J. M., Ryan, B. R., Dolence, J. C., et al. 2019b, Phys. Rev. D, 100, 023008
  * Murguia-Berthier et al. (2017) Murguia-Berthier, A., Ramirez-Ruiz, E., Kilpatrick, C. D., et al. 2017, ApJ, 848, L34
  * Nakar (2007) Nakar, E. 2007, Phys. Rep., 442, 166
  * Nakar (2019) —. 2019, arXiv e-prints, arXiv:1912.05659
  * Narayan et al. (1992) Narayan, R., Paczynski, B., & Piran, T. 1992, ApJ, 395, L83
  * Narayan et al. (2001) Narayan, R., Piran, T., & Kumar, P. 2001, ApJ, 557, 949
  * Nedora et al. (2021) Nedora, V., Bernuzzi, S., Radice, D., et al. 2021, ApJ, 906, 98
  * Neilsen et al. (2014) Neilsen, D., Liebling, S. L., Anderson, M., et al. 2014, Phys. Rev. D, 89, 104029
  * Nicholl et al. (2017) Nicholl, M., Berger, E., Kasen, D., et al. 2017, ApJ, 848, L18
  * Noble et al. (2006) Noble, S. C., Gammie, C. F., McKinney, J. C., & Del Zanna, L. 2006, ApJ, 641, 626
  * Noble et al. (2009) Noble, S. C., Krolik, J. H., & Hawley, J. F. 2009, ApJ, 692, 411
  * Noble et al. (2010) —. 2010, ApJ, 711, 959
  * Noble et al. (2012) Noble, S. C., Mundim, B. C., Nakano, H., et al. 2012, ApJ, 755, 51
  * O’Connor & Ott (2010) O’Connor, E., & Ott, C. D. 2010, Classical and Quantum Gravity, 27, 114103
  * Palenzuela et al. (2015) Palenzuela, C., Liebling, S. L., Neilsen, D., et al. 2015, Phys. Rev. D, 92, 044045
  * Perego et al. (2014) Perego, A., Rosswog, S., Cabezón, R. M., et al. 2014, MNRAS, 443, 3134
  * Pian et al. (2017) Pian, E., D’Avanzo, P., Benetti, S., et al. 2017, Nature, 551, 67
  * Popham et al. (1999) Popham, R., Woosley, S. E., & Fryer, C. 1999, ApJ, 518, 356
  * Powell (1968) Powell, M. J. 1968, A Fortran subroutine for solving systems of nonlinear algebraic equations, Tech. rep., United Kingdom Atomic Energy Authority
  * Press et al. (1992) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in C, 2nd edn. (Cambridge, USA: Cambridge University Press)
  * Radice et al. (2020) Radice, D., Bernuzzi, S., & Perego, A. 2020, Annual Review of Nuclear and Particle Science, 70, 95
  * Radice et al. (2016) Radice, D., Galeazzi, F., Lippuner, J., et al. 2016, MNRAS, 460, 3255
  * Radice et al. (2018) Radice, D., Perego, A., Hotokezaka, K., et al. 2018, ApJ, 869, 130
  * Ravi & Lasky (2014) Ravi, V., & Lasky, P. D. 2014, MNRAS, 441, 2433
  * Roberts et al. (2011) Roberts, L. F., Kasen, D., Lee, W. H., & Ramirez-Ruiz, E. 2011, ApJ, 736, L21
  * Roberts et al. (2017) Roberts, L. F., Lippuner, J., Duez, M. D., et al. 2017, MNRAS, 464, 3907
  * Rosswog et al. (2017) Rosswog, S., Feindt, U., Korobkin, O., et al. 2017, Classical and Quantum Gravity, 34, 104001
  * Rosswog & Liebendörfer (2003) Rosswog, S., & Liebendörfer, M. 2003, MNRAS, 342, 673
  * 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
  * Ruiz et al. (2016) Ruiz, M., Lang, R. N., Paschalidis, V., & Shapiro, S. L. 2016, ApJ, 824, L6
  * Ruiz et al. (2018) Ruiz, M., Shapiro, S. L., & Tsokaros, A. 2018, Phys. Rev. D, 97, 021501
  * Sano et al. (2004) Sano, T., Inutsuka, S.-i., Turner, N. J., & Stone, J. M. 2004, ApJ, 605, 321
  * Schneider et al. (2017) Schneider, A. S., Roberts, L. F., & Ott, C. D. 2017, Phys. Rev. C, 96, 065802
  * Shappee et al. (2017) Shappee, B. J., Simon, J. D., Drout, M. R., et al. 2017, Science, 358, 1574
  * Shibata & Taniguchi (2006) Shibata, M., & Taniguchi, K. 2006, Phys. Rev. D, 73, 064027
  * Siegel (2019) Siegel, D. M. 2019, European Physical Journal A, 55, 203
  * Siegel et al. (2019) Siegel, D. M., Barnes, J., & Metzger, B. D. 2019, Nature, 569, 241
  * 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 (2018) —. 2018, ApJ, 858, 52
  * Siegel et al. (2018) Siegel, D. M., Mösta, P., Desai, D., & Wu, S. 2018, ApJ, 859, 71
  * Smartt et al. (2017) Smartt, S. J., Chen, T. W., Jerkstrand, A., et al. 2017, Nature, 551, 75
  * Soares-Santos et al. (2017) Soares-Santos, M., Holz, D. E., Annis, J., et al. 2017, ApJ, 848, L16
  * Sorathia et al. (2012) Sorathia, K. A., Reynolds, C. S., Stone, J. M., & Beckwith, K. 2012, ApJ, 749, 189
  * Tanaka & Hotokezaka (2013) Tanaka, M., & Hotokezaka, K. 2013, ApJ, 775, 113
  * Tanvir et al. (2017) Tanvir, N. R., Levan, A. J., González-Fernández, C., et al. 2017, ApJ, 848, L27
  * Villar et al. (2017) Villar, V. A., Guillochon, J., Berger, E., et al. 2017, ApJ, 851, L21
  * Zalamea & Beloborodov (2011) Zalamea, I., & Beloborodov, A. M. 2011, MNRAS, 410, 2302
  * Zenati et al. (2019) Zenati, Y., Perets, H. B., & Toonen, S. 2019, MNRAS, 486, 1805
  * Zilhão & Noble (2014) Zilhão, M., & Noble, S. C. 2014, Classical and Quantum Gravity, 31, 065013