# Bound Debris Expulsion from Neutron Star Merger Remnants

Yossef Zenati ISEF International Fellowship  Physics and Astronomy Department,
Johns Hopkins University, Baltimore, MD 21218, USA  Julian H. Krolik Physics
and Astronomy Department, Johns Hopkins University, Baltimore, MD 21218, USA
Leonardo R. Werneck Department of Physics, University of Idaho, Moscow, ID
83843, USA  Ariadna Murguia-Berthier NASA Einstein fellow  Center for
Interdisciplinary Exploration and Research in Astrophysics (CIERA), 1800
Sherman Ave., Evanston, IL 60201, USA  Zachariah B. Etienne Department of
Physics, University of Idaho, Moscow, ID 83843, USA  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  Scott C. Noble Gravitational
Astrophysics Lab, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
Tsvi Piran Racah Institute of Physics, Hebrew University, Jerusalem, 91904,
Israel Yossef Zenati<EMAIL_ADDRESS>
(March 2023)

###### Abstract

Many studies have found that neutron star mergers leave a fraction of the
stars’ mass in bound orbits surrounding the resulting massive neutron star or
black hole. This mass is a site of $r-$ process nucleosynthesis and can
generate a wind that contributes to a kilonova. However, comparatively little
is known about the dynamics determining its mass or initial structure. Here we
begin to investigate these questions, starting with the origin of the disk
mass. Using tracer particle as well as discretized fluid data from numerical
simulations, we identify where in the neutron stars the debris came from, the
paths it takes in order to escape from the neutron stars’ interiors, and the
times and locations at which its orbital properties diverge from those of
neighboring fluid elements that end up remaining in the merged neutron star.

Neutron Star — General relativity — Accretion disk — black hole

††software: astropy (Astropy Collaboration et al., 2018), IllinoisGRMHD
(Etienne et al., 2015), Numpy (Harris et al., 2020), and Matplotlib (Hunter,
2007).

##  1 Introduction

On August 17, 2017, the LIGO/Virgo collaboration (LVC) detected the first
gravitational wave (GW) signal from a binary neutron star (BNS) merger, the
GW170817 event (Abbott et al., 2017a; Villar et al., 2017). This was also the
first multi-messenger detection of such an event, as it was seen across the
electromagnetic (EM) spectrum (Fermi, Chandra, UV, optical, IR, radio).

From GW170817 and similar BNS mergers, we can learn about many fundamental
questions in general relativity, cosmology, and astrophysics (Coulter et al.,
2017; Guidorzi et al., 2017; Haggard et al., 2017; Margutti et al., 2017). For
example, the distance to this event inferred from GWs was consistent with
entirely independent measurements of the Hubble constant (Abbott et al.,
2017b). Analysis of the GW170817 GW form led to important constraints on the
nuclear equation of state (EOS) (Damour et al., 2012; Radice et al., 2018b;
Bernuzzi, 2020). Its optical/IR emission (optical/IR analysis papers) gave
support to the idea that such mergers could be a main channel of rapid
neutron-capture ($r$-process) nucleosynthesis (e.g, Lattimer & Schramm, 1974;
Eichler et al., 1989; Arcones & Thielemann, 2013; Cowan et al., 2021; Pian,
2023). In addition, it also demonstrated that, as predicted many years before,
neutron star (NS) mergers can create short $\gamma$-ray bursts (see, Eichler
et al., 1989; Rosswog & Ramirez-Ruiz, 2002; Nakar, 2007; Berger, 2014;
Goldstein et al., 2017)

There is a strong theoretical consensus that a portion of the neutron stars’
original mass is not held in the remnant, whether it is a massive neutron star
or a black hole. Mass can be separated during numerous stages of a BNS merger:
the late inspiral phase, the merger stage, and the subsequent evolution of the
initially bound debris. Multiple mechanisms can contribute, including tidal
gravity, shock heating, neutrino heating, nuclear heating, magneto-centrifugal
winds (from both a debris disk and highly-magnetized massive neutron star
merger remnant), and magnetohydrodynamic (MHD) turbulent dissipation (see,
Metzger et al., 2010; Hotokezaka et al., 2013; Fernández et al., 2015; Radice
et al., 2020; Margutti & Chornock, 2021; Sarin & Lasky, 2021; Most & Quataert,
2023; Combi & Siegel, 2023; Shibata et al., 2023).

Some of the matter that escapes from the remnant is ejected from the system
altogether during the late inspiral and the merger itself. This part is called
“dynamical ejecta”, and it is generally ascribed to a combination of tidal
forces and shock heating when the two stars come into direct contact
(Hotokezaka et al., 2013). Its mass is rather small, ${\sim}10^{-3}M_{\odot}$
(Radice et al., 2016; Siegel & Metzger, 2017; Fujibayashi et al., 2018;
Shibata et al., 2023). The dynamical ejecta have received much attention for
their possible role in producing a kilonova (KN) event, particularly those
features indicating especially high outflow velocities (Metzger et al., 2010;
Barnes & Kasen, 2013; Radice et al., 2018a; Domoto et al., 2021; Curtis et
al., 2023).

However, the same studies predict that the rest of the matter left outside the
remnant is, at least temporarily, bound; its mass is rather larger,
${\sim}$0.01–0.1$M_{\odot}$. Here we call this portion the “bound debris”. It,
too, likely plays a major role in neutron star merger phenomenology. The
magnetic flux and confinement pressure needed by relativistic jets may be
supplied by this debris. In addition, on rather longer timescales than
dynamical ejection, the heat released inside the bound debris by
nucleosynthesis and dissipation of MHD turbulence may drive a wind carrying
away neutron-enriched isotopes. The photosphere of such a wind could also
contribute to kilonova radiation (e.g, Rosswog & Ramirez-Ruiz, 2002; Fernández
et al., 2015; Just et al., 2015; Korobkin et al., 2020).

Despite its likely importance to the phenomenology of neutron star mergers,
there has been little effort to understand the mechanisms governing how much
mass is in bound orbits or the nature of its dynamical state, i.e., its
distribution of angular momentum and energy: Most et al. (2021) studied the
correlations between the angular momentum distribution of debris and other
properties, but only for NS-BH mergers. The distribution of mass with angular
momentum has special interest because it determines the initial surface
density profile of the bound debris, which is almost certainly far from inflow
equilibrium. The initial angular momentum distribution, therefore, governs the
bound debris’ subsequent evolution.

In this paper, we will begin the study of these issues, starting with the most
basic questions: From which part(s) of the NSs is the bound debris drawn? How,
and at what point in the merger, is the bound debris separated from the
majority of the merged NS mass, the part that remains in the newly-formed NS
or is carried into the black hole (BH) if one is created?

Figure 1: Schematic presentation of the different stages and outcomes of BNS
mergers: from left to right, inspiral, merger, post-merger remnant plus debris
disk; and, possibly, a stable NS. Timescales for each stage are shown below
the corresponding sketches. The bottom line lists the GW frequency and the
principal physical mechanisms of each stage.

##  2 Calculation

###  2.1 Overview

To place our calculation in context, we first review the stages of a BNS
merger (see Fig. 1). Once formed, a NS binary can evolve only by emission of
gravitational radiation. During very nearly the entire span of this period,
the stars’ separation is so large that there is no direct interaction between
them, and even gravitational tidal forces are negligible. When, however, their
separation shrinks to only a few times their radii, tidal forces begin to
distort them; this is the period we call “late inspiral”. Eventually, they
touch, beginning the “merger proper”, which creates a very massive NS (MNS).
Depending upon the mass and rotational state of this newly-formed star, it may
collapse to a BH (almost immediately or after some delay) or live for an
indefinitely long time. In the course of these events, a small fraction of the
initial total mass is expelled.

When the mass distribution is substantially asymmetric (the inspiral, merger,
and possible collapse phases), there is no simple description of the
spacetime.

Consequently, it is necessary to employ numerical relativity methods in order
to follow gravitational dynamics. On the other hand, once the mass
distribution relaxes to axisymmetry, either as a stable NS or a Kerr BH, there
is no need for a further solution of the Einstein field equations. For this
reason, we divide our simulation of such events into two pieces. The first
part is performed with the numerical relativity plus MHD code IllinoisGRMHD
(Etienne et al., 2015, 2020; Werneck et al., 2023); the second part uses
HARM3D+NUC (Murguia-Berthier et al., 2021). A “hand-off” procedure assures
consistency between the two in terms of the spacetime and the fluid properties
(Lopez Armengol et al., 2022). In the present paper, we make use of data taken
only from the first part; future work in this project will employ data from
the second part for studies of the bound debris’ evolution, including jet-
launching, nucleosynthesis, and wind-driving.

###  2.2 Physics treated and equations solved

IllinoisGRMHD solves the equations of MHD,

$\displaystyle\nabla_{\mu}\left(n_{\rm b}u^{\mu}\right)$ $\displaystyle=0\;,$
(1) $\displaystyle\nabla_{\mu}\left(n_{\rm e}u^{\mu}\right)$
$\displaystyle={\cal R}\;,$ (2) $\displaystyle\nabla_{\mu}T^{\mu\nu}$
$\displaystyle=\mathcal{Q}u^{\nu}\;,$ (3)
$\displaystyle\nabla_{\mu}{}^{*}\\!F^{\mu\nu}$ $\displaystyle=0\;,$ (4)

corresponding to the conservation of the baryon number, conservation of Lepton
number, conservation of energy-momentum, and the two homogeneous Maxwell’s
equations, respectively. In coordination with the solution of these equations,
the BSSN evolution thorn Baikal solves the Einstein Field Equations.

In these equations, $n_{\rm b}$ ($n_{\rm e}$) is the baryon (lepton) number
density and $u^{\mu}$ is the fluid four-velocity. The net rate of change in
the lepton number ${\cal R}$ is the difference between the creation rates of
electron anti-neutrinos and neutrinos. The cooling rate ${\cal Q}$ is the rate
at which neutrinos (and anti-neutrinos) carry energy out of the gas. It is
calculated from a local emission model and a leakage rate formalism (see
Murguia-Berthier et al., 2021; Werneck et al., 2023, for details about the
calculation of ${\cal R}$ and ${\cal Q}$).

The energy-momentum tensor is assumed to be that of a perfect fluid plus an EM
contribution,

$T^{\mu\nu}=\left(\rho_{\rm
b}h+b^{2}\right)u^{\mu}u^{\nu}+\left(P+\frac{b^{2}}{2}\right)g^{\mu\nu}-b^{\mu}b^{\nu},$
(5)

where $\rho_{\rm b}=m_{\rm b}n_{\rm b}$ is the baryon density, $m_{\rm b}$ is
the mean baryon mass,
${}^{*}\\!F^{\mu\nu}=(1/2)\tilde{\epsilon}^{\mu\nu\rho\sigma}F_{\rho\sigma}$
is the dual of the Faraday tensor $F^{\mu\nu}$, and
$\tilde{\epsilon}^{\mu\nu\rho\sigma}$ is the Levi-Civita tensor. The parameter
$h=1+\epsilon+P/\rho_{\rm b}$ is the specific enthalpy, $\epsilon$ is the
specific internal energy, and $P$ is the fluid pressure, which is given
through a tabulated form of the SFHo EOS for hot degenerate matter (O’Connor &
Ott, 2010).111The SFHo EOS table was downloaded from
http://stellarcollapse.org.

Further, $b^{\mu}=(4\pi)^{-1/2}B^{\mu}_{(u)}$ is the rescaled magnetic
4-vector in the fluid frame, where

$\displaystyle B^{0}_{(u)}$ $\displaystyle=u_{i}B^{i}/\alpha\;,$ (6)
$\displaystyle B^{i}_{(u)}$
$\displaystyle=\bigl{(}B^{i}/\alpha+B^{0}_{(u)}u^{i}\bigr{)}/u^{0}\;.$ (7)

Here $B^{i}$ is the magnetic field in the frame normal to the hypersurface,
$b^{2}\equiv b^{\mu}b_{\mu}$ is the magnetic energy density, $g_{\mu\nu}$ is
the spacetime metric, and $\alpha$ is the lapse function. The lapse, the shift
vector, and the metric are all defined on spatial hypersurfaces of constant
coordinate time $t$. Ideal MHD ($u_{\mu}F^{\mu}=0$) is assumed throughout.

In addition to solving these equations, we also follow the locations and
velocities of a large number of tracer particles. The initial positions of the
tracer particles are selected by randomly choosing points $p=(x,y,z)$ within a
radius $R_{\rm seed}=75$ km from the system center of mass and accepting the
particle if

$\frac{\rho_{\rm b}(x,y,z)}{\max\rho_{\rm b}}>\zeta\;,$ (8)

where $\zeta\in[0,1]$ is a random number. The process is repeated iteratively
until all $N_{\rm particles}=50,000$ have been seeded. The fact that the NSs
are extremely compact and have central densities many orders of magnitude
larger than that of the surrounding gas ensures that the vast majority of the
tracer particles are seeded within each NS, as shown in Fig.
LABEL:fig:rho_particles.

Once the positions $p_{n}=(x_{n},y_{n},z_{n})$ of the particles are known,
where $n=0,1,\ldots,N_{\rm particles}-1$, we update them in time using

$\frac{dp^{i}_{n}}{dt}=v^{i}(p_{n})\;,$ (9)

where $v^{i}(p_{n})$ is the fluid three-velocity $v^{i}=u^{i}/u^{0}$
interpolated to the particle’s position $p_{n}$. The tracer particle positions
are updated every 16 local time steps, and the three-velocities are obtained
using fourth-order Lagrange interpolation.

###  2.3 Numerical setup

We perform BNS simulations using the Einstein Toolkit (Loffler et al., 2012;
Haas et al., 2022) 222See https://github.com/zachetienne/nrpytutorial. applied
to a Cartesian grid with adaptive mesh refinement (AMR) provided by Carpet
(Schnetter et al., 2004). Initially, the grid contains eight refinement levels
differing by factors of two, such that the resolution at the finest refinement
level is ${\approx}185$ m. Once the minimum value of the lapse function drops
below 0.1, BH formation is assumed to be imminent and two additional
refinement levels are added, such that the resolution of the finest refinement
level becomes ${\approx}46$ m.

At the outer boundary, located at ${\approx}5670\ \rm{km}$, we apply radiation
boundary conditions to the metric quantities using Newrad (Alcubierre et al.,
2003; Loffler et al., 2012) and a simple copy boundary condition for the MHD
fields.

The spacetime is evolved using Baikal (Haas et al., 2022), while the MHD
fields are evolved using a newly developed version of IllinoisGRMHD (Werneck
et al., 2023) that supports finite-temperature, microphysical EOSs, and
neutrino physics via a leakage scheme. In addition to the MHD fields, we
evolve the entropy by assuming that it is conserved (see e.g., Noble et al.,
2009, for a similar strategy), an approximation that, while poor at shocks,
allows us to use the entropy as a backup variable during primitive recovery
(see Werneck et al., 2023, for more details).

Initial data for an equal-mass binary is constructed using LORENE (Gourgoulhon
et al., 2001; Feo et al., 2017; Gourgoulhon et al., Accessed in 2022), with
each NS having a baryonic (gravitational) mass of 1.550$M_{\odot}$
(1.348$M_{\odot}$; see Table 1). The initial separation of the stars is 45 km
and the stellar radii are ${\approx}9.3$ km. Each is given a strong poloidal
magnetic field (see e.g., Appendix A of Etienne et al., 2015) such that
$\max\bigl{(}P_{\rm mag}/P\bigr{)}=\max\bigl{(}b^{2}/2P\bigr{)}=10^{-4}$,
corresponding to $\max\sqrt{b^{2}}=5.05\times 10^{15}~{}G$.

Parameter | Value
---|---
Initial data
EoS | SFHo
NS baryonic mass | 1.550$M_{\odot}$
NS gravitational mass | 1.348$M_{\odot}$
NS radius $R_{*}$ | 9.3 km
Initial separation | 45 km
$\max\sqrt{b^{\mu}b_{\mu}}$ | $5.05\times 10^{15}G$
Number of tracer particles | $50000$
Post-merger
BH dimensionless spin parameter, $\chi=a/M$ | 0.795
BH irreducible mass, $M_{\rm irr}$ | 2.444$M_{\odot}$
BH mass, $M_{\rm BH}$ | 2.726 $M_{\odot}$
Mass-weighted mean $Y_{\rm e}$ $\left(r=60\,km,t\simeq 17msec\right)$ | $0.122$
Table 1: Simulation parameters for an equal-mass, magnetized BNS merger
performed with IllinoisGRMHD using a microphysical, finite-temperature EOS and
a neutrino leakage scheme.

Where polar coordinates are more easily interpretable, we define such a system
through the simple coordinate transformation $r^{2}=(x^{2}+y^{2}+z^{2})$,
$\cos\theta=z/r$, $\tan\phi=y/x$, and the $xy$-plane is identical to the
initial binary orbital plane. In these coordinates, the polar component of
angular momentum is $u_{\phi}\equiv xu_{y}-yu_{x}$ becomes a conserved
quantity when the spacetime is axisymmetric.

##  3 Results

###  3.1 Original location of matter escaping the merger remnant

Figure 2: The final distance $R_{f}$ (from the system center-of-mass) of the
surviving tracer particles plotted against their initial distance $R_{i}$ from
the center of their home neutron star. Red circles distinguish unbound
particles.

The most basic distinction between the tracer particles is their post-merger
fate: unbound, bound, or confined within either a stable NS or a BH. We define
unbound as having $hu_{t}<-1$, where $h$ is the specific enthalpy of the fluid
surrounding the particle and $u_{t}$ is the particle’s conserved energy. A
bound particle is one with $hu_{t}>-1$, but located (in this case) outside the
innermost stable circular orbit (ISCO) of the BH that forms (for this case,
${\simeq}3.7r_{g}=14.8$ km from the center of the BH). We call the latter two
categories “survivors”.

The simplest description of the particles’ post-merger position is their
radial distance from the system center-of-mass, $R_{f}$; here “post-merger”
means at the end of the simulation, ${\sim}18$ ms after its start. One might
expect that particles initially placed near the surface of a NS would be more
likely to “survive” than particles whose starting position was deeper inside
one of the stars. The relation between $R_{f}$ and the distance $R_{i}$
between each survivor particle and the center of its home NS is shown in
Figure 2. Only ${\simeq}950$ of the particles survive, 1.9% of the original
number. Roughly half are found outside $15r_{g}\simeq 60\rm{km}$, and
${\sim}20\%$ outside $100\rm{km}$. Almost all those outside $\rm 200\ km$
(${\simeq}8\%$ of the surviving particles) are unbound. Most importantly for
our present purposes, there is no correlation at all between $R_{i}$ and
$R_{f}$ except for a slight (and statistically weak) tendency for the unbound
particles to originate at larger $R_{i}$ than the bound particles. Note that
both $R_{i}$ and $R_{f}$ are coordinate distances defined by
$\left[(x-x_{c})^{2}+(y-y_{c})^{2}+(z-z_{c})^{2}\right]^{1/2}$, where
$\left(x_{c},y_{c},z_{c}\right)$ are either the coordinates of the center-of-
mass ($x_{c}=y_{c}=z_{c}=0$) or the coordinates of the nearest NS center.

However, there is a sense in which particles whose initial home was at the
surface of a NS are more likely to survive from the remnant. As shown by the
upper panel of Figure 3, the probability that a tracer particle escapes grow
for $R_{i}>8$ km, i.e., within ${\sim}1$ km of the surface. Compared to
particles with smaller $R_{i}$ this probability is greater by
${\sim}3\text{--}4$. The reason why Figure 2 shows no correlation is revealed
by the bottom panel of Figure 3: the probability density for the initial
radius of a surviving particle is considerably larger over the range
$R_{i}\approx 2\text{--}8$ km than for $R_{i}\gtrsim 8$ km. Even if a large
fraction of all the surface particles escapes, while only a small fraction of
interior particles do, there is so little mass near the surface compared to
the mass in the interior that the great majority of surviving particles
originate in the deep interior of a NS.

Figure 3: Upper panel: The probability that a particle with initial radius
$R_{i}$ survives black hole collapse. Lower panel: The probability density (in
arbitrary units) that a surviving particle had initial radius $R_{i}$.

###  3.2 Escape of debris from the NSs’ interiors

Although answering the first question we raised, this result prompts a second:
how does the bulk of the debris find its way from the deep interior of a NS to
an orbit well outside the merged NS or BH? Put another way, how is it that
${\sim}$2% of the matter in the interior of a NS is somehow separated from its
immediate neighbors and given enough energy and angular momentum to be
released from the merged NS?

In order to approach this question, we begin by identifying the time from the
beginning of the simulation at which this separation occurs. We do so by
choosing a small sample of surviving particles and then pairing each with the
closest non-surviving particle in the initial state of the system. These were
selected by finding survivors for which $R_{*}-2$ km$\ <R_{i}<R_{*}$ in a
single neutron star and were all located in different quadrants. We also
constructed comparison samples varying these choices, but the results were
quite similar.

In Figure 4, we display how the positions of these particles projected onto
the $x$-$y$ (orbital) plane evolve from $t=13.0$ ms, 0.2 ms before the two
stars first touch, to 17.62 ms, 1 ms after the black hole forms. The
oscillations reflect the rotation of the NSs, first as they orbit one another,
and later as they rotate as an asymmetric merged object. As this figure shows,
four of the _six_ pairs begin to separate at ${\simeq}14.5\ \rm{ms}$, while
the other two do not substantially part ways until ${\simeq}15.5$ ms. This
small sample suggests the critical time for distinguishing surviving particles
from captured particles is $\simeq$14.5–15.5 ms; as we will see shortly, this
suggestion is confirmed when we examine the data more closely.

Figure 4: The evolution of the separation between _six_ surviving particles
and their initial close neighbors. The top panel shows how their _y-_
coordinates evolve, the bottom panel shows the evolution of their _x-_
coordinates. Each pair’s starting location $(X,Y)$ is shifted by an arbitrary
amount relative to the other pairs so that they can be visually distinguished.
In some cases (e.g., the bottom two), the members of the pair start out so
close together that their curves are superposed.

To see what is happening during this crucial millisecond, we present images of
the density and $Y_{\rm e}$ maps at that time (Fig. 5). The first of the
density images ($\rm t=14.66$ ms) shows the system ${\simeq}0.5$ ms after the
star’s first touch. Over the course of the following ${\sim}2$ ms, the
configuration of the NSs changes from one in which the two stars remain
partially distinguishable but are in contact over a sizable surface to one in
which most of the stars’ mass has taken on a much rounder shape, but a
noticeable amount of gas has been expelled into spiral arms.

Complementary images of the specific internal energy density are shown in
Figure 6. For these, we begin at a slightly earlier time ($t=14.50$ ms) and
follow it with the same three times shown in Fig. 5. The internal energy
images make shock locations easily visible. In particular, at the two earlier
times, they show clearly the shock at the contact surface, while at the two
later times they who how it weakens. They further point out shocks in
surrounding gas driven by both polar outflow and the rotating asymmetric
merged star (Hotokezaka et al., 2013).

Thus, the moment at which the paths of escaping and captured particles diverge
is the moment when the two stars are in direct contact, but have not yet
combined their masses into a single (in this case, short-lived) NS. During the
first part of this stage, the mass distribution is not only highly asymmetric;
it also rotates around the center of mass, so spacetime is both asymmetric and
time-dependent. However, during the latter part, the mass distribution becomes
significantly more axisymmetric, weakening the azimuthal component of
gravitational acceleration. Throughout this period, as shown by the numerous
narrow regions of elevated internal energy per unit mass in Figure 6, there
are many shocks, where kinetic energy is dissipated into heat, and deflection
can exchange angular momentum between different fluid elements.

|  |
---|---|---
|  |
|  |
(a) $Y_{\rm e}$ | (b) $\log\rho_{\rm xy}$ | (c) $\log\rho_{\rm xz}$
Figure 5: Evolution of the electron fraction, _xy-_ plane mass density, and _xz-_ plane mass density. $Y_{\rm e}$ is shown on a linear color scale (_xy-_ plane), the mass densities on a logarithmic scale. Each quantity is presented at times $\rm 14.66\ ms$, $\rm 15.29\ ms$, $\rm 16.31\ ms$, and $\rm 17.25\ ms$; the last is shortly after the collapse to a BH. Note that the colorscale for $Y_{\rm e}$ at $t=14.66$ ms is different from the colorscale at the other times.  |
---|---
|
|
|
(a) $\log\epsilon_{\rm xy}$ | (b) $\log\epsilon_{\rm xz}$
Figure 6: Evolution of the specific internal energy $\epsilon_{\rm int}$ in
the _xy-_ and _xz-_ planes, shown on a logarithmic scale at times $\rm 14.50\
ms$, $\rm 14.66\ ms$, $\rm 15.37\ ms$; $\rm 16.39\ ms$.

To understand more deeply _why_ the surviving particles diverge from their
neighboring captured particles, we now take a closer look at the _six_
surviving particles in our sample. Figure 7 displays their $x$-$y$ positions
at $t=14.50$ ms in relation to the fluid density in the orbital plane. Four of
the _six_ particles’ locations trace a curve in the _x-y_ plane running along
the outer edge of the merged star; the other two are found along the “seam-
line” between the two merging stars. Their _z-_ coordinates inform us that
these particles are typically $\sim(0.1-0.3)R_{*}$ out of the orbital plane,
where $R_{*}$ is the radius of the original NSs.

To interpret this 3D geometry, we first plot the $x$-$y$ locations of _all_
the tracer particles at this same time (blue dots in the top panel of Fig. 8).
The four particles in our _six_ -member sample found along the edge of the
merged star are all found in the lower hook of the $\mathbf{``S"}$-curve
traced out by the entire ensemble of particles; the densest parts of both
hooks coincide with the plumes of gas density visible in Figure 7. The other
two are in the central bar of the $\mathbf{``S"}$-curve, which coincides
closely with the seam-line. It, therefore, appears that particles escape when
they move outward through the seam region and then reach the surface. The
curved portions of the $\mathbf{``S"}$-curve result from the time sequence at
which the particles crossed through the surface of the merged NS. Note that
the orbital plane velocities of the survivor particles are primarily
azimuthal, and the radial gradient in their speeds indicates that they acquire
greater angular momentum as they move outward.

Figure 7: The $x$-$y$ locations of the same _six_ surviving particles shown in
Fig. 4 superposed on the $t=14.50\rm{ms}$ density map. Their _z-_ coordinates
are $z=2.558,-2.225,-3.097,2.182,1.738,3.504$ km for the particles shown as a
plus-sign, a square, a circle, a dot, an ’X’, and a caret, respectively.

Figure 8: (Top): Position in the _xy-_ plane of the surviving particles at
$t=14.506\ \rm{msec}$ (blue dots) and at $t=14.705\ \rm{msec}$ (yellow
crosses). The line labeled “$l_{x}$” is fit to the line-like feature in the
middle of the $\mathbf{``S"}$-shape at $t=14.50\ \rm{msec}$. The red arrows
illustrate the _xy-_ velocity of a random sample of surviving particles at
$t=14.50$ ms; the arrows’ lengths are proportional to the speeds of the
particles’ orbital plane motion. (Bottom): Position in the _$l_{x}$ z-_plane
of the surviving particles in the line-like feature at $t=14.50$ ms.

Further information about vertical structure can be derived from examining the
locations of the particles within the bar at $t=14.50$ ms. Their orbital-plane
coordinates all lie along the linear feature in the $\mathbf{``S"}$, so that
their 3D positions all lie in a vertical plane defined by the line feature and
the $z-$axis (the simple geometry of this structure is likely due to our
choice of an equal-mass merger). The bottom panel of Figure 8 shows their
distribution in this plane. All of the particles lie a short distance inside
the outline of the merged star, a figure that is flattened in the vertical
direction (see the images in the third column of Fig. 5). Moreover, they are
weakly concentrated toward the orbital plane and tend to avoid the polar axis.

The evolution of this structure in time is illustrated by the difference
between the figures outlined by the blue and gold dots in the upper panel of
Figure 8. At $t=14.50$ ms, numerous particles can be found in the plane of the
inter-star seam, while the “hooks”, taken together, extend over $\sim 3\pi/2$
in azimuthal angle. Only 0.2 ms later, hardly any particles remain in the
seam, while the hooks now cover all $2\pi$ in azimuth. Thus, the time period
we have singled out as critical in terms of separating surviving particles
from initially neighboring captured particles coincides with the time when the
expulsion of matter from the seam region is completed.

The variety of paths these particles follow in order is illustrated in Figure
9, which shows three particle trajectories during the $\sim 2$ ms spanning the
time of first contact between the stars. Just before contact, the tracer
particles follow their star’s orbit, but moving in and out relative to the
system center-of-mass as the star rotates. Once the two stars come into
contact (at $t=13.2$ ms), the character of their motion changes: now they are
part of the merged NS, partaking in its rotation around the system center of
mass. However, the axial fluid motions within the contact surface permit
particles to travel both outward and inward; the red path in Figure 9 is an
example of the latter option. Large amplitude motions within this zone are
expected because the strong shear at the boundary between the two stars
excites turbulence through the Kelvin-Helmholtz instability (Rasio & Shapiro,
1999; Kiuchi et al., 2014; Aguilera-Miret et al., 2022). By $t=14.7$ ms, only
a few particles remain in the seam-plane. Instead, the particles destined to
survive that had been inside the merger remnant at 14.5 ms are almost all
found in a shell surrounding the remnant. The survivor particles that at 14.5
ms had been in the hooks of the “S” are now in either of two places: roughly
half are in the shell surrounding the remnant, many of them having moved
inward to do so; the other half are outside the boundary of the figure.

Figure 9: Motion in $(x,y)$ of three surviving particles (blue, green, and red
dotted paths) from $t=12$ ms to $t=16.6$ ms, when the star collapses to a BH.
The black circle shows the outline of their parent star at $t=12$ ms. Cusps
are the signature of shock deflection. Figure 10: Time-dependence of the
angular momentum $u_{\phi}$ for our sample of _six_ pairs of particles.
Survivors are shown with solid lines, and captured particles with dashed
lines; matching colors identify the pairs. The time of collapse to a BH,
$16.6$ ms is marked with a gray dashed vertical line; for captured particles
after that time, the physical significance of $u_{\phi}$ is unclear.

###  3.3 Separation of surviving particles from neighboring particles

The surviving particles are $\simeq 2\%$ of the NSs’ mass; as can be seen in
Figure 8, $\simeq 90\%$ have reached positions just outside the star by $\sim
1$ ms after contact. From the data illustrated in Figure 4, it is clear that
the captured particles that began close to the survivors are still close to
them at this time. It immediately follows that the distinction between
survivors and the captured matter is accomplished in two stages. First, a few
percent of the total mass of NSs is conveyed to locations outside the stars by
motion through the contact surface. Second, when this matter reaches locations
immediately outside the merged star, local effects divide this mass into two
subsets, one that escapes, the other remaining with the merged star.

How this distinction is made is portrayed in Figure 10. For the first $\sim 2$
ms shown, the angular momenta of the surviving particles are very similar to
those of their matching captured particles; both oscillate as a result of the
NSs’ rotation. At around $\rm 14-15\ ms$, the angular momentum of surviving
particles begins to rise (as shown by the arrows in Fig. 8), while the angular
momentum of the neighboring particles that remain with the star for the long
run do not change appreciably. This growth in survivor particle angular
momentum finishes shortly before collapse to a BH ($t=16.6$ ms). By then, the
surviving particles’ angular momentum is $\sim 2-10\times$ the angular
momentum of the captured particles that started out from locations very near
them.

That these results are characteristic of _all_ the surviving particles are
shown in Figure 11. A small number of particles have large excursions in
$u_{\phi}$ during the inspiral, but the overwhelming majority of particles
follow the pattern of the small sample we have chosen for a particular study.

The torques creating the change in angular momentum could be the result of
non-axisymmetric gravity (Hotokezaka et al., 2013; Radice et al., 2018b;
Fujibayashi et al., 2018; Shibata et al., 2023), smooth pressure gradients, or
shocks (Hotokezaka et al., 2011, 2013; Sekiguchi et al., 2016). If the
magnetic field is as strong as Kiuchi et al. (2014) found (but not seen here),
magnetic forces might also contribute. Shocks can be recognized by the sharp
steps in angular momentum evident in Figure 10 as well as the cusps seen in
the trajectories of Figure 9, which might be associated with shock features
visible in the later images of Figure 6. Both gravitational acceleration and
smooth pressure gradients could lead to smooth portions of trajectories. The
short timescale smaller fluctuations in $u_{\phi}$ shown in Figure 10 are more
likely associated with pressure gradients than with gravity because there are
many short length-scale density gradients, whereas the gravitational field is
the result of the global stress-energy distribution. The general upward trend
could be associated with non-axisymmetric gravity, and its end at $t\approx
16$ ms roughly coincides with the axisymmetrizing of the merged NS (see Fig.
5). On the other hand, it is unclear how a gravitational effect would
distinguish the particles that ultimately survive from the nearby particles
that are ultimately captured; non-axisymmetric gravity could become important
only after other dynamics place fluid elements in the regions where
gravitational torques exist. In this context, it is also worth noting that the
captured particles show much less short timescale fluctuation in $u_{\phi}$,
suggesting that the regions through which they travel have much smoother
density and pressure distributions. This would, of course, be a natural
consequence of staying closer to the merged star.

Figure 11: Evolution of $u_{\phi}$, the conserved angular momentum component
in an axisymmetric spacetime, for a large sample of the survivor tracer
particles: $220$ particles selected at random. The small number of surviving
particles with higher $u_{\phi}$ and noisy $u_{\phi}$ histories (blue and
green lines) were initially located very close to the surface of one or the
other of the NSs.

##  4 Discussion

###  4.1 Initial electron fraction in the bound debris

The electron fraction $Y_{\rm e}$, the ratio of the number of protons $n_{p}$
to the total number of nucleons, is one of the key quantities for determining
the elements synthesized over longer timescales in the bound debris. Through
this influence, it is also the key to determining the opacity of the
optical/UV photosphere where the KN emission is radiated (see, Barnes & Kasen,
2013; Hotokezaka et al., 2013; Siegel & Metzger, 2017; Radice et al., 2018a).

If there were special locations in the NSs that contributed disproportionately
to the bound debris, the distribution of mass with electron fraction
$dM/dY_{\rm e}$ in its initial state could be quite different from $dM/dY_{\rm
e}$ in the NSs as a whole. However, we have shown that, in rough terms, the
probability that a fluid element inside one of the NSs is expelled from the
remnant is very nearly constant for all radii between $\simeq 1$ and $\simeq
8$ km, a region containing the great majority of mass in the star. There are
only two exceptions to this constant probability. One is a spherical shell
$\simeq 1$ km thick at the outside of the star, where the probability of
escape rises to a few times greater than elsewhere. The other is a sphere of
radius $\sim 1$ km at the center of the star, where the escape probability is
$\lesssim 0.1\times$ the value in the star’s bulk. Consequently, $dM/dY_{\rm
e}$ in the pre-merger NSs should be a very good predictor of $dM/dY_{\rm e}$
in the bound debris immediately after it leaves the merged NS.

On the other hand, it should also be borne in mind that $Y_{\rm e}$ can be
affected by events happening during the merger and while the bound debris
travel from the merged NS to its place in the surrounding disk. Shock heating
in the merger and during the subsequent evolution of the merger remnant can
raise the matter temperature above $\rm 10MeV$. In such a high-temperature
environment, electron-positron pair creation can substantially raise the rate
of the reaction $e^{+}+n\rightarrow p^{+}+\bar{\nu}_{e}$, thereby increasing
$Y_{\rm e}$. Neutrino irradiation, within the merged star, as the gas travels
to the disk, and in the disk can also raise $Y_{\rm e}$ via reactions like
$\nu_{e}+n\rightarrow p^{+}+e^{-}$. The degree to which any of these
mechanisms may or may not alter $Y_{e}$ requires further investigation
(Shibata et al., 2023).

###  4.2 Parameter-dependence

Quantitative details in the single calculation on which we base our results
could change for different values of the principal system parameters, the mass
ratio $q$ and the total mass $M$. In particular, the bilateral symmetry so
evident in the images of Figures 5, 7, and 8 will be broken when $q\neq 1$. It
has also long been known (Oechslin & Janka, 2006; Shibata et al., 2023) that
the total mass removed from the merger remnant increases for smaller values of
$q$, so the total amount of bound debris mass in this simulation is on the
lower end of what might be expected for a merger of this total mass.

The results should also depend on $M$. If it is so large that the merger
results in immediate collapse to a BH, the amount of expelled mass, both in
dynamical ejecta and bound debris, would be much smaller. Moreover, because
the tidal gravity contribution to the expelled mass would be a larger fraction
of the total, the bound fraction would likely be considerably smaller. On the
other hand, if $M$ were small enough for the merged NS to survive for a longer
time, the expelled mass could be rather larger.

The combination of different $M$ and different $q$ can alter the results
quantitatively in several ways. Smaller values of both imply particularly
small mass for the less-massive star, giving it both lower mean density and a
less centrally-concentrated internal density profile (see, Lattimer & Prakash,
2007; Bauswein et al., 2012; Shibata et al., 2021; Raithel et al., 2022).
Relatively large tidal mass-loss would likely follow. Larger $M$ and smaller
$q$ could enhance the contribution of the contact shock if the larger mass
star has $\gtrsim 2M_{\odot}$ because for such large masses, the radius of the
NS diminishes, whereas it hardly changes (for most EOSs) for masses between
$\sim 0.7$ and $\sim 2M_{\odot}$ (Kashyap et al., 2022). The sense in which
these changes tip the balance between bound and unbound debris requires a
detailed calculation to determine.

###  4.3 Uncertainty in the equation of state

Although the gravitational waveform measured in the event GW170817 has helped
prune the long list of proposed nuclear EOSs (see, Hotokezaka et al., 2011;
Bauswein et al., 2013; Radice et al., 2018b), there still remain many
candidates. Our quantitative results could be altered in several ways if a
different EOS applies. First, the internal density profile of an equilibrium
NS is sensitive to the EOS, and that could alter the probability of escape as
a function of radius. Second, when the two stars come into contact, the
pressure of the shocked region between the two depends on the EOS. Stiffer
EOSs support larger NSs that strike each other with lower velocity,
diminishing the post-shock temperature (Hotokezaka et al., 2011); on the other
hand, a larger ratio between pressure and internal energy density would make
the shocked region wider and enlarge the pressure contrast between the shocked
region and the surface of the merged star. Third, the EOS governs both the
critical mass dividing quick collapse from a longer-lived merger NS and the
time to collapse if rapid collapse is, indeed, the result (Bauswein et al.,
2013; Radice et al., 2018b; Margalit et al., 2022). A longer lifetime for the
merged NS likely leads to more escaping mass, and the ratio of bound to
unbound mass could change over this span of time. Thus, different EOSs could
lead to changes in both the total mass removed from the remnant and how it is
divided into bound and unbound portions.

##  5 Conclusions

In much of the literature on NS mergers, when matter that escapes being
confined within the merger remnant is discussed, the focus is solely on the
unbound fraction. This fraction is then further divided into two categories,
“dynamical” and “secular” (sometimes called “post-merger”) (Metzger et al.,
2010; Fernández et al., 2015; Siegel & Metzger, 2017; Fujibayashi et al.,
2018; Radice et al., 2018a). “Dynamical” refers to the material expelled
during the merger itself, whether by tidal forces, shocks heating gas, or
neutrinos that heat or push the gas. “Secular” refers to matter driven off the
disk of bound debris many milliseconds after the merger.

However, the majority of the matter outside the merger remnant remains bound
to the system, forming the disk that is the source of the secular ejecta and
may also support any jet that forms. Hitherto, the mechanisms by which this
material leaves the merging neutron stars and is deposited in the debris disk
have gone largely unexamined. Here we have determined the origin of the bound
debris within the NSs and examined how its fate is distinguished from that of
matter remaining inside the merged star.

In particular, by making use of the complementary information provided by
Eulerian fluid data and Lagrangian tracer particle data, we have uncovered
several surprising facts about the process by which matter is expelled from
merging neutron stars.

  * •

The probability that a given parcel of mass leaves the merger remnant is
almost, but not quite, independent of its location inside the pre-merger NSs.
Matter placed near the surface of a NS is considerably more likely than most
matter to be expelled, but there is so much more matter in the bulk of the
star that, despite its lower escape probability, nearly all the escaping
matter, whether bound or unbound, has its origin deep in a star’s interior.

  * •

Turbulent motions within the “seam” surface formed when the two stars touch
can mix matter from a wide range of initial depths within the stars. These
same motions transport some of this matter through the seam plane to the
surface of the merged star.

  * •

A significant amount of matter is brought out to the surface of the merged NS
but does not escape. Which matter is captured and which escapes depends on
details of the shock structure in the region where the matter that has been
squeezed outward by the shock between the two merging stars reaches the
surface.

  * •

The surviving matter acquires significant additional angular momentum as it
leaves the remnant. It is likely that torques due to the non-axisymmetric
gravity of the merged NS, smooth pressure gradients, and shocks all
contribute.

  * •

All of these processes occur within a few milliseconds after the stars first
touch.

Thus, we have shown that the process by which matter escapes the remnant
created by a neutron star merger is considerably more complex and subtle than
previously envisioned. Because the disk of bound matter circulating around the
merged NS or BH is far from inflow equilibrium, its subsequent evolution
depends on what sort of non-equilibrium structure it is given. Understanding
the mechanisms regulating this structure is, therefore, an essential part of
understanding many of the dramatic observational properties of NS
mergers—jets, kNe, and $r-$process nucleosynthesis.

## Acknowledgements

This work was partially supported by NASA TCAN grant NNH17ZDA001N. YZ thanks
Prof. Will Farr, Prof. Yuri Levin, and Mor Rozner for their helpful
discussions. L.R.W. and Z.B.E. gratefully acknowledge support from NSF awards
PHY-1806596, PHY-2110352, OAC-2004311, as well as NASA award
ISFM-80NSSC18K0538. A.M-B. is supported by NASA through the NASA Hubble
Fellowship grant HST-HF2-51487.001-A awarded by the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., for NASA, under contract NAS5-26555. T.P. was supported in
part by ERC advanced grant “Multijets”. This research made use of Idaho
National Laboratory computing resources which are supported by the Office of
Nuclear Energy of the U.S. Department of Energy and the Nuclear Science User
Facilities under Contract No. DE-AC07-05ID14517.

## References

  * Abbott et al. (2017a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, Nature, 551, 85, doi: 10.1038/nature24471
  * Abbott et al. (2017b) —. 2017b, ApJ, 848, L12, doi: 10.3847/2041-8213/aa91c9
  * Aguilera-Miret et al. (2022) Aguilera-Miret, R., Viganò, D., & Palenzuela, C. 2022, ApJ, 926, L31, doi: 10.3847/2041-8213/ac50a7
  * Alcubierre et al. (2003) Alcubierre, M., Bruegmann, B., Diener, P., et al. 2003, Phys. Rev. D, 67, 084023, doi: 10.1103/PhysRevD.67.084023
  * Arcones & Thielemann (2013) Arcones, A., & Thielemann, F. K. 2013, Journal of Physics G Nuclear Physics, 40, 013201, doi: 10.1088/0954-3899/40/1/013201
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Barnes & Kasen (2013) Barnes, J., & Kasen, D. 2013, ApJ, 775, 18, doi: 10.1088/0004-637X/775/1/18
  * Bauswein et al. (2013) Bauswein, A., Baumgarte, T. W., & Janka, H. T. 2013, Phys. Rev. Lett., 111, 131101, doi: 10.1103/PhysRevLett.111.131101
  * Bauswein et al. (2012) Bauswein, A., Janka, H. T., Hebeler, K., & Schwenk, A. 2012, Phys. Rev. D, 86, 063001, doi: 10.1103/PhysRevD.86.063001
  * Berger (2014) Berger, E. 2014, ARA&A, 52, 43, doi: 10.1146/annurev-astro-081913-035926
  * Bernuzzi (2020) Bernuzzi, S. 2020, General Relativity and Gravitation, 52, 108, doi: 10.1007/s10714-020-02752-5
  * Combi & Siegel (2023) Combi, L., & Siegel, D. M. 2023, arXiv e-prints, arXiv:2303.12284, doi: 10.48550/arXiv.2303.12284
  * Coulter et al. (2017) Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Science, 358, 1556, doi: 10.1126/science.aap9811
  * Cowan et al. (2021) Cowan, J. J., Sneden, C., Lawler, J. E., et al. 2021, Reviews of Modern Physics, 93, 015002, doi: 10.1103/RevModPhys.93.015002
  * Curtis et al. (2023) Curtis, S., Mösta, P., Wu, Z., et al. 2023, MNRAS, 518, 5313, doi: 10.1093/mnras/stac3128
  * Damour et al. (2012) Damour, T., Nagar, A., & Villain, L. 2012, Phys. Rev. D, 85, 123007, doi: 10.1103/PhysRevD.85.123007
  * Domoto et al. (2021) Domoto, N., Tanaka, M., Wanajo, S., & Kawaguchi, K. 2021, ApJ, 913, 26, doi: 10.3847/1538-4357/abf358
  * Eichler et al. (1989) Eichler, D., Livio, M., Piran, T., & Schramm, D. N. 1989, Nature, 340, 126, doi: 10.1038/340126a0
  * Etienne et al. (2020) Etienne, Z. B., Paschalidis, V., Haas, R., Moesta, P., & Shapiro, S. L. 2020, IllinoisGRMHD: GRMHD code for dynamical spacetimes, Astrophysics Source Code Library, record ascl:2004.003.  http://ascl.net/2004.003
  * Etienne et al. (2015) Etienne, Z. B., Paschalidis, V., Haas, R., Mösta, P., & Shapiro, S. L. 2015, Class. Quant. Grav., 32, 175009, doi: 10.1088/0264-9381/32/17/175009
  * Feo et al. (2017) Feo, A., De Pietri, R., Maione, F., & Löffler, F. 2017, Class. Quant. Grav., 34, 034001, doi: 10.1088/1361-6382/aa51fa
  * Fernández et al. (2015) Fernández, R., Kasen, D., Metzger, B. D., & Quataert, E. 2015, MNRAS, 446, 750, doi: 10.1093/mnras/stu2112
  * Fujibayashi et al. (2018) Fujibayashi, S., Kiuchi, K., Nishimura, N., Sekiguchi, Y., & Shibata, M. 2018, ApJ, 860, 64, doi: 10.3847/1538-4357/aabafd
  * Goldstein et al. (2017) Goldstein, A., Veres, P., Burns, E., et al. 2017, ApJ, 848, L14, doi: 10.3847/2041-8213/aa8f41
  * Gourgoulhon et al. (Accessed in 2022) Gourgoulhon, E., Grandclément, P., & Novak, J. Accessed in 2022, LORENE: Langage Objet pour la RElativité NumériquE
  * Gourgoulhon et al. (2001) Gourgoulhon, E., Grandclement, P., Taniguchi, K., Marck, J.-A., & Bonazzola, S. 2001, Phys. Rev. D, 63, 064029, doi: 10.1103/PhysRevD.63.064029
  * Guidorzi et al. (2017) Guidorzi, C., Margutti, R., Brout, D., et al. 2017, ApJ, 851, L36, doi: 10.3847/2041-8213/aaa009
  * Haas et al. (2022) Haas, R., Cheng, C.-H., Diener, P., et al. 2022, The Einstein Toolkit, The ”Sophie Kowalevski” release, ET_2022_11, Zenodo, doi: 10.5281/zenodo.7245853
  * Haggard et al. (2017) Haggard, D., Nynka, M., Ruan, J. J., et al. 2017, ApJ, 848, L25, doi: 10.3847/2041-8213/aa8ede
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357–362, doi: 10.1038/s41586-020-2649-2
  * Hotokezaka et al. (2013) Hotokezaka, K., Kiuchi, K., Kyutoku, K., et al. 2013, Phys. Rev. D, 87, 024001, doi: 10.1103/PhysRevD.87.024001
  * Hotokezaka et al. (2011) Hotokezaka, K., Kyutoku, K., Okawa, H., Shibata, M., & Kiuchi, K. 2011, Phys. Rev. D, 83, 124008, doi: 10.1103/PhysRevD.83.124008
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Just et al. (2015) Just, O., Bauswein, A., Ardevol Pulpillo, R., Goriely, S., & Janka, H. T. 2015, MNRAS, 448, 541, doi: 10.1093/mnras/stv009
  * Kashyap et al. (2022) Kashyap, R., Das, A., Radice, D., et al. 2022, Phys. Rev. D, 105, 103022, doi: 10.1103/PhysRevD.105.103022
  * Kiuchi et al. (2014) Kiuchi, K., Kyutoku, K., Sekiguchi, Y., Shibata, M., & Wada, T. 2014, Phys. Rev. D, 90, 041502, doi: 10.1103/PhysRevD.90.041502
  * Korobkin et al. (2020) Korobkin, O., Wollaeger, R., Fryer, C., et al. 2020, arXiv e-prints, arXiv:2004.00102.  https://arxiv.org/abs/2004.00102
  * Lattimer & Prakash (2007) Lattimer, J. M., & Prakash, M. 2007, Phys. Rep., 442, 109, doi: 10.1016/j.physrep.2007.02.003
  * Lattimer & Schramm (1974) Lattimer, J. M., & Schramm, D. N. 1974, ApJ, 192, L145, doi: 10.1086/181612
  * Loffler et al. (2012) Loffler, F., et al. 2012, Class. Quant. Grav., 29, 115001, doi: 10.1088/0264-9381/29/11/115001
  * Lopez Armengol et al. (2022) Lopez Armengol, F. G., Etienne, Z. B., Noble, S. C., et al. 2022, Phys. Rev. D, 106, 083015, doi: 10.1103/PhysRevD.106.083015
  * Margalit et al. (2022) Margalit, B., Jermyn, A. S., Metzger, B. D., Roberts, L. F., & Quataert, E. 2022, ApJ, 939, 51, doi: 10.3847/1538-4357/ac8b01
  * Margutti & Chornock (2021) Margutti, R., & Chornock, R. 2021, ARA&A, 59, 155, doi: 10.1146/annurev-astro-112420-030742
  * Margutti et al. (2017) Margutti, R., Berger, E., Fong, W., et al. 2017, ApJ, 848, L20, doi: 10.3847/2041-8213/aa9057
  * Metzger et al. (2010) Metzger, B. D., Martínez-Pinedo, G., Darbha, S., et al. 2010, MNRAS, 406, 2650, doi: 10.1111/j.1365-2966.2010.16864.x
  * Most et al. (2021) Most, E. R., Papenfort, L. J., Tootle, S. D., & Rezzolla, L. 2021, MNRAS, 506, 3511, doi: 10.1093/mnras/stab1824
  * Most & Quataert (2023) Most, E. R., & Quataert, E. 2023, ApJ, 947, L15, doi: 10.3847/2041-8213/acca84
  * Murguia-Berthier et al. (2021) Murguia-Berthier, A., Noble, S. C., Roberts, L. F., et al. 2021, ApJ, 919, 95, doi: 10.3847/1538-4357/ac1119
  * Nakar (2007) Nakar, E. 2007, Phys. Rep., 442, 166, doi: 10.1016/j.physrep.2007.02.005
  * Noble et al. (2009) Noble, S. C., Krolik, J. H., & Hawley, J. F. 2009, Astrophys. J., 692, 411, doi: 10.1088/0004-637X/692/1/411
  * O’Connor & Ott (2010) O’Connor, E., & Ott, C. D. 2010, Class. Quant. Grav., 27, 114103, doi: 10.1088/0264-9381/27/11/114103
  * Oechslin & Janka (2006) Oechslin, R., & Janka, H. T. 2006, MNRAS, 368, 1489, doi: 10.1111/j.1365-2966.2006.10238.x
  * Pian (2023) Pian, E. 2023, Universe, 9, 105, doi: 10.3390/universe9020105
  * Radice et al. (2020) Radice, D., Bernuzzi, S., & Perego, A. 2020, Annual Review of Nuclear and Particle Science, 70, 95, doi: 10.1146/annurev-nucl-013120-114541
  * Radice et al. (2016) Radice, D., Galeazzi, F., Lippuner, J., et al. 2016, MNRAS, 460, 3255, doi: 10.1093/mnras/stw1227
  * Radice et al. (2018a) Radice, D., Perego, A., Hotokezaka, K., et al. 2018a, ApJ, 869, 130, doi: 10.3847/1538-4357/aaf054
  * Radice et al. (2018b) Radice, D., Perego, A., Zappa, F., & Bernuzzi, S. 2018b, ApJ, 852, L29, doi: 10.3847/2041-8213/aaa402
  * Raithel et al. (2022) Raithel, C. A., Espino, P., & Paschalidis, V. 2022, MNRAS, 516, 4792, doi: 10.1093/mnras/stac2450
  * Rasio & Shapiro (1999) Rasio, F. A., & Shapiro, S. L. 1999, Classical and Quantum Gravity, 16, R1, doi: 10.1088/0264-9381/16/6/201
  * Rosswog & Ramirez-Ruiz (2002) Rosswog, S., & Ramirez-Ruiz, E. 2002, MNRAS, 336, L7, doi: 10.1046/j.1365-8711.2002.05898.x
  * Sarin & Lasky (2021) Sarin, N., & Lasky, P. D. 2021, General Relativity and Gravitation, 53, 59, doi: 10.1007/s10714-021-02831-1
  * Schnetter et al. (2004) Schnetter, E., Hawley, S. H., & Hawke, I. 2004, Class. Quant. Grav., 21, 1465, doi: 10.1088/0264-9381/21/6/014
  * Sekiguchi et al. (2016) Sekiguchi, Y., Kiuchi, K., Kyutoku, K., Shibata, M., & Taniguchi, K. 2016, Phys. Rev. D, 93, 124046, doi: 10.1103/PhysRevD.93.124046
  * Shibata et al. (2023) Shibata, M., Fujibayashi, S., Hayashi, K., Kiuchi, K., & Wanajo, S. 2023, IAU Symposium, 362, 190, doi: 10.1017/S1743921322001351
  * Shibata et al. (2021) Shibata, M., Fujibayashi, S., & Sekiguchi, Y. 2021, Phys. Rev. D, 104, 063026, doi: 10.1103/PhysRevD.104.063026
  * Siegel & Metzger (2017) Siegel, D. M., & Metzger, B. D. 2017, Phys. Rev. Lett., 119, 231102, doi: 10.1103/PhysRevLett.119.231102
  * Villar et al. (2017) Villar, V. A., Guillochon, J., Berger, E., et al. 2017, ApJ, 851, L21, doi: 10.3847/2041-8213/aa9c84
  * Werneck et al. (2023) Werneck, L. R., Etienne, Z. B., Murguia-Berthier, A., et al. 2023, Phys. Rev. D, 107, 044037, doi: 10.1103/PhysRevD.107.044037