# GRMHD simulations of neutron-star mergers with weak interactions:
r-process nucleosynthesis and electromagnetic signatures of dynamical ejecta

Luciano Combi CITA National Fellow  Instituto Argentino de Radioastronomía
(IAR, CCT La Plata, CONICET/CIC)
C.C.5, (1984) Villa Elisa, Buenos Aires, Argentina  Perimeter Institute for
Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada  Department of Physics,
University of Guelph, Guelph, Ontario N1G 2W1, Canada  Daniel M. Siegel
Institute of Physics, University of Greifswald, D-17489 Greifswald, Germany
Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada
Department of Physics, University of Guelph, Guelph, Ontario N1G 2W1, Canada

###### Abstract

Fast neutron-rich material ejected dynamically over $\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}\\!10$ ms during the
merger of a binary neutron-star (BNS) system can give rise to distinctive
electromagnetic counterparts to the system’s gravitational-wave emission that
can serve as a ‘smoking gun’ to distinguish between a BNS and a NS—black-hole
merger. We present novel ab-initio modeling of the associated kilonova
precursor and kilonova afterglow based on three-dimensional general-
relativistic magneto-hydrodynamic (GRMHD) simulations of BNS mergers with
tabulated, composition-dependent, finite-temperature equations of state
(EOSs), weak interactions, and approximate neutrino transport. We analyze
dynamical mass ejection from 1.35–1.35 $M_{\odot}$ binaries, typical of the
observed Galactic double–NS systems and consistent with inferred properties of
the first observed BNS merger GW170817, using three nuclear EOSs that span the
range of allowed compactness of $1.35\,M_{\odot}$-neutron stars. Nuclear
reaction network calculations yield a robust 2nd-to-3rd-peak r-process. We
find $\text{few}\times 10^{-6}M_{\odot}$ of fast ($v>0.6c$) ejecta that give
rise to broad-band synchrotron emission on $\sim\\!\text{yr}$ timescales,
consistent with recent tentative evidence for excess X-ray/radio emission
following GW170817. We find $\approx\\!2\times 10^{-5}\,M_{\odot}$ of free
neutrons that power a kilonova precursor on $\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}\\!\text{h}$
timescale. A boost in early UV/optical brightness by a factor of a few due to
previously neglected relativistic effects, with appreciable enhancements up to
$\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}\\!10$ h post-merger, provides promising
prospects for future detection with UV/optical telescopes such as _Swift_ or
ULTRASAT out to $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}\\!250\,\text{Mpc}$. We find that a
recently predicted opacity boost due to highly ionized lanthanides at
$\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}\\!70000$ K is unlikely to affect the
early kilonova lightcurve based on the obtained ejecta structures. Azimuthal
inhomogeneities in dynamical ejecta composition for soft EOSs found here
(“lanthanide/actinide pockets”) may have observable consequences for both
early kilonova and late-time nebular emission.

††software: The Einstein Toolkit (Loffler et al. 2012;
http://einsteintoolkit.org), LORENE (Gourgoulhon et al., 2001), GRMHD_con2prim
(Siegel et al. 2018, Siegel & Mösta 2018), Skynet (Lippuner & Roberts, 2017),
Paramo (Rueda-Becerril, 2020), afterglowpy (Ryan et al. 2020), WhiskyTHC
(Radice & Rezzolla 2012, Radice et al. 2014a, Radice et al. 2014b), PyCactus
(Kastaun & Ohme 2021, https://github.com/wokast/PyCactus), Matplotlib (Hunter,
2007), NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020), and hdf5
(The HDF Group, 2000-2010).

[1] #1

##  1 Introduction

The first detection of gravitational waves from the merger of a binary neutron
star (BNS) system, called GW170817 (Abbott et al., 2017a), was associated with
fireworks of electromagnetic (EM) counterparts across the electromagnetic
spectrum. The largest EM follow-up campaign ever conducted (Abbott et al.,
2017b) revealed an associated short gamma-ray burst (GRB) GRB170817, followed
by broad-band GRB afterglow emission, as well as the quasi-thermal counterpart
AT2017gfo—the first unambiguous detection of a kilonova (Li & Paczynski, 1998;
Metzger et al., 2010; Barnes & Kasen, 2013; Metzger, 2020), powered by
radioactive heating from the production of heavy elements via the rapid
neutron-capture process (r-process; Burbidge et al. 1957; Cameron 1957) in
dense, neutron-rich plasma ejected from the merger site.

The GRB afterglow of GW170817 showed several interesting and unusual
properties and has been followed up for several years until today. The
faintness of the prompt GRB gamma-ray emission given the proximity of only 40
Mpc, the unusual rising X-ray flux appearing after nine days (Troja et al.,
2017), peaking at $\approx\\!160$ days after the merger, followed by a steep
decay (Haggard et al., 2017; Hallinan et al., 2017; D’Avanzo et al., 2018;
Lyman et al., 2018; Margutti et al., 2017; Mooley et al., 2018; Troja et al.,
2018; Ruan et al., 2018), as well as the observation of superluminal motion of
the radio centroid, are interpreted as a structured jet expanding into the
interstellar medium (ISM) viewed at an angle of $\theta_{\rm obs}\sim
20^{o}-30^{o}$ from the jet core (Lamb & Kobayashi, 2017; Lamb et al., 2018;
Alexander et al., 2017; Hotokezaka et al., 2018; Wu & MacFadyen, 2019; Fong et
al., 2017; Ghirlanda et al., 2019; Hajela et al., 2019; Lamb et al., 2019;
Troja et al., 2020). After several years, the broad-band synchrotron afterglow
originating in the decelerating external shock driven into the ISM has faded
to levels at which other possible emission components may be revealed, with
first observational indications (Hajela et al., 2022; Troja et al., 2022). One
robust additional emission component expected on timescales of years after the
event is the broad-band synchrotron emission from a transrelativistic shock of
fast dynamical ejecta expanding into the ISM (e.g., Nakar & Piran 2011;
Hotokezaka & Piran 2015; Hotokezaka et al. 2018; Nedora et al. 2021a;
Balasubramanian et al. 2022), which we refer to here as the ‘kilonova
afterglow’. Late-time fall-back accretion onto a remnant black hole offers a
possible alternative explanation to a rebrightening (Metzger & Fernández,
2021). The existence of a rebrightening in GW170817 is still a matter of
debate (Troja et al., 2022; O’Connor & Troja, 2022; Balasubramanian et al.,
2022).

The astrophysical sites that give rise to the synthesis of heavy elements via
the r-process remain a topic of active debate (see Horowitz et al. 2019; Cowan
et al. 2021; Siegel 2022 for recent reviews). Several lines of evidence,
including measurements of radioactive isotopes in deposits on the deep
seafloor (e.g., Wallner et al. 2015; Hotokezaka et al. 2015), abundances of
metal-poor stars formed in the smallest dwarf galaxies (e.g., Ji et al. 2016;
Tsujimoto et al. 2017), in the halo of the Milky Way (Macias & Ramirez-Ruiz,
2018), as well as in globular clusters (Zevin et al., 2019; Kirby et al.,
2020), point to a low-rate, high-yield source, much rarer than ordinary core-
collapse supernovae both in early and recent Galactic history (see Fig. 3 of
Siegel 2019 for a compilation of observational constraints). From the most
promising candidates, which include neutron-star mergers (Lattimer & Schramm,
1974; Symbalisty & Schramm, 1982), rare types of core-collapse supernovae
producing rapidly spinning, strongly magnetized proto-neutron stars
(magnetorotational supernovae; Winteler et al. 2012; Symbalisty et al. 1985;
Nishimura et al. 2006), and hyper-accreting black holes (“collapsars”; Pruet
et al. 2003; Surman et al. 2006; Siegel et al. 2019, 2021; see also Grichener
& Soker 2019), definitive evidence for the production of r-process elements
has only been obtained in the case of neutron-star mergers through the
GW170817 kilonova so far.

Figure 1: Volume renderings of rest-mass density layers from our APR4 BNS
simulation illustrating the different merger stages and associated ejecta
components. From left to right: (i) The two neutron stars reach high angular
speeds before the collision and tidally deform; (ii) At merger, the neutron
star cores bounce off of each other, generating shock-heated material that is
ejected at high speeds, represented here as a red density layer; (iii) Due to
tidal forces during merger, ‘tidal tails’ develop that eject cold and highly
neutron-rich material, represented here by a blue density layer; (iv) After
merger much of the debris material circularizes around the remnant neutron
star to form an accretion disk. Winds launched from this disk eject a
significant fraction of the total unbound material of this BNS merger. The hot
remnant neutron star additionally launches a neutrino and magnetically driven
wind in the polar direction, highlighted here by a green density layer.

Several mechanisms can give rise to the ejection of neutron-rich material
conducive to undergoing r-process nucleosynthesis and powering kilonovae in
BNS mergers111Provided the total mass of the BNS system is below a certain
threshold, typically $\approx\\!1.3-1.7\,M_{\rm TOV}$, where $M_{\rm TOV}$ is
the maximum mass for non-rotating NSs (e.g., Bauswein et al. 2013a). Above
this threshold the merging BNS system promptly collapses into a black hole,
resulting in negligible amounts of ejecta.. In the dynamical merger phase
(typically lasting $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}\\!10$ ms), unbound matter is ejected from
the system by tidal forces (sometimes forming ‘tidal tails’; Rosswog et al.
1999; Rosswog 2013; Korobkin et al. 2012; Goriely et al. 2011), by shocks
generated at the collision interface of the two stars (Hotokezaka et al.,
2013a; Oechslin et al., 2007; Sekiguchi et al., 2016), and by a sequence of
quasi-radial oscillations (‘bounces’; Bauswein et al. 2013b; Radice et al.
2018) of a double-core structure of the remnant that forms immediately upon
merger (see Fig. 1; ‘dynamical ejecta’).

After merger, a remnant NS may be formed, temporarily supported from
collapsing into a black hole by solid-body or differential rotation (Duez et
al., 2006a; Kaplan et al., 2014; Siegel et al., 2014), which can unbind
material in a combination of neutrino-driven (Dessart et al., 2009; Perego et
al., 2014; Desai et al., 2022) and magnetically driven winds (Siegel et al.,
2014; Ciolfi & Kalinani, 2020; Metzger et al., 2018; Mösta et al., 2020;
Curtis et al., 2021) over the lifetime and Kelvin-Helmholtz cooling time of
the remnant ($\sim 10\,\text{ms}-1\,\text{s}$). Bound material surrounding the
remnant circularizes and forms a neutrino-cooled accretion disk (Fernández &
Metzger, 2013; Just et al., 2015; Siegel & Metzger, 2017; Lippuner et al.,
2017; Fujibayashi et al., 2018; Fernández et al., 2019; Miller et al., 2019a;
Nedora et al., 2021b) with typical masses of $0.01–0.3\,{M_{\odot}}$. A
combination of heating-cooling imbalance in the disk corona in the presence of
MHD turbulence (Siegel & Metzger, 2017, 2018) at early times and viscously-
driven outflows at later evolution stages (Fernández et al., 2019), together
with nuclear binding energy release as seed particles for the r-process form,
can eject up to $\sim 30-40\%$ of the initial disk mass.

The GW170817 kilonova was likely generated by a combination of
$\approx\\!0.01\,M_{\odot}$ of dynamical ejecta and winds from the remnant NS
for the early blue emission, and by $\approx\\!0.05\,M_{\odot}$ post-merger
disk ejecta for the red emission (see, e.g., Siegel 2019; Metzger 2020; Radice
et al. 2020; Margutti & Chornock 2021; Pian 2021 for reviews overseeing the
interpretation of the event). The likely dominance of disk outflows over other
ejecta channels may extend to the population of merging BNS systems in general
if these systems follow a distribution similar to the Galactic double-neutron
star systems observable as radio pulsars (Siegel, 2022). The population of
merging BNS systems is a matter of active debate (e.g., Landry et al. 2020)
and future combined gravitational wave and EM observations will be necessary
to infer the observed population. In a population sense, the contribution of
neutron-star–black-hole (NS–BH) systems to Galactic r-process nucleosynthesis
is likely by far subdominant relative to BNS mergers (Chen et al., 2021).

Although likely subdominant with respect to other ejecta channels in many
scenarios, dynamical ejecta can give rise to distinctive EM signals associated
with its most energetic component. The outermost $\sim\\!10^{-8}\,M_{\odot}$
of dynamical ejecta might reach ultra-relativistic velocities and serve as a
breakout medium for a relativistic jet and give rise to prompt gamma-ray
emission similar to GRB 170817A as envisioned by Beloborodov et al. (2020).
The transrelativistic tail of ejecta with velocity $\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}(0.5-0.6)c$ drives a
shock into the ISM, giving rise to the aforementioned broad-band ‘kilonova
afterglow’. Furthermore, this transrelativistic component typically has a
sufficiently rapid expansion timescale ($<5$ ms) to produce free neutrons,
preventing rapid neutron capture to occur in these fastest layers. Free-
neutron decay on a timescale of $\approx\\!15$ min acts with a specific
heating rate roughly an order of magnitude larger than typical r-process
heating at a similar epoch in the outermost layers of the ejecta, giving rise
to a short-lived thermal UV/optical kilonova “precursor” transient on a
timescale of $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}\\!1$ hr after merger (Kulkarni, 2005;
Metzger et al., 2015; Ishii et al., 2018). Such early blue emission might
overlap (and be confused) with additional heating of the outermost ejecta
layers by a cocoon of hot shock-heated material surrounding a jet propagating
into the ejecta envelope (e.g., Kasliwal et al. 2017; Gottlieb et al. 2018).
The aforementioned signatures, if detected early on after merger in future
gravitational-wave events, will help to distinguish between a BNS and an NS–BH
merger, a distinction that is presently not possible based on gravitational-
wave data alone (except when exclusions can be made based on plausible ranges
for inferred component masses). This motivates a detailed investigation of the
physical mechanisms related to dynamical ejecta, a characterization of its
properties, and self-consistent modeling of its unique EM signatures based on
first-principle numerical simulations, as attempted here.

In this paper, we present a set of general-relativistic magnetohydrodynamic
(GRMHD) simulations of equal-mass BNSs typical of the Galactic double-neutron
star systems and consistent with the inferred source parameters of GW170817.
These simulations self-consistently combine several physical ingredients from
inspiral to post-merger, including magnetic fields, tabulated (finite
temperature, composition-dependent) equations of state (EOS), weak
interactions, and approximate neutrino transport via a ray-by-ray scheme. We
employ these simulations to explore in detail the physical ejection mechanisms
of dynamical ejecta, the nucleosynthesis it gives rise to, and its distinctive
EM signatures including neutron precursor emission and the kilonova afterglow.
We focus here on building novel, self-consistent models for the EM signatures
that are directly based on ab-initio numerical simulations of BNS mergers.

The paper is organized as follows. In Sec. 2, we describe our simulation
setup, including numerical methods, initial data, and simulation diagnostics.
In Sec. 3, we discuss simulation results on dynamical ejecta, provide a
comprehensive discussion of ejection mechanisms, and conclude with a brief
description of the post-merger phase. Section 4 reflects on nucleosynthesis
results of the dynamical ejecta. In Sec. 5, we present models of the kilonova
emission including the contribution of free-neutron heating (Secs. 5.1 and
5.2), of the non-thermal kilonova afterglow (Sec. 5.5), and discuss their
application in the context of GW170817 (Secs. 5.2.2 and 5.5.3). We also
discuss the role of exceptionally high opacities associated with high
ionization states of lanthanides recently calculated by Banerjee et al. (2022)
with respect to early kilonova and precursor emission (Sec. 5.3). We also
compute high-energy gamma-ray emission from inverse-Compton and synchrotron
self-Compton processes associated with the kilonova afterglow (Sec. 5.5.4).
Conclusions are presented in Sec. 6. Finally, two appendices elaborate on
numerical details regarding convergence, ejecta analysis, and the use of
passive tracer particles.

##  2 Simulation details

In the following subsections, we briefly discuss the analytic foundations
(Secs. 2.1 and 2.2), numerical methods (Sec. 2.3), initial data (Sec. 2.4) as
well as simulation dignostics (Sec. 2.5) employed here for simulating BNS
mergers. Our methods closely follow Siegel & Metzger (2018) for the GRMHD part
with weak interactions and Radice et al. (2018) for modeling neutrino
transport via a one-moment closure scheme (‘M0 scheme’).

###  2.1 GRMHD with weak interactions

We model BNS mergers within ideal GRMHD coupled to weak interactions. The
equations of ideal GRMHD comprise energy-momentum conservation, baryon and
lepton number conservation, Maxwell’s equations in the ideal MHD regime, and
Einstein’s equations:

$\displaystyle\nabla_{\mu}T^{\mu\nu}$ $\displaystyle=$ $\displaystyle
Qu^{\nu},$ (1) $\displaystyle\nabla_{\mu}(n_{\mathrm{b}}u^{\mu})$
$\displaystyle=$ $\displaystyle 0,$ (2)
$\displaystyle\nabla_{\mu}(n_{\mathrm{e}}u^{\mu})$ $\displaystyle=$
$\displaystyle R,$ (3) $\displaystyle\nabla_{\nu}F^{*\mu\nu}$ $\displaystyle=$
$\displaystyle 0,$ (4) $\displaystyle\quad u_{\mu}F^{\mu\nu}$ $\displaystyle=$
$\displaystyle 0,$ (5)

and

$G_{\mu\nu}=8\pi T_{\mu\nu}.$ (6)

Here, $G_{\mu\nu}$ is the Einstein tensor, $u^{\mu}$ denotes the four-velocity
of the fluid, $n_{\mathrm{b}}$, and $n_{\mathrm{e}}$ are the baryon and
electron number density, $F^{*\mu\nu}$ is the dual of the Faraday
electromagnetic tensor, and the quantities $Q$ and $R$ represent source terms
due to weak interaction and neutrino transport (neutrino absorption and
emission). The energy-momentum tensor is given by

$T^{\mu\nu}=\left(\rho
h+b^{2}\right)u^{\mu}u^{\nu}+\left(p+\frac{b^{2}}{2}\right)g^{\mu\nu}-b^{\mu}b^{\nu},$
(7)

where $\rho=n_{\mathrm{b}}m_{\mathrm{b}}$ is the rest-mass density and
$m_{\mathrm{b}}$ is the baryon mass, $p$ is the pressure,
$h=1+\epsilon+p/\rho$ denotes the specific enthalpy, $\epsilon$ is the
specific internal energy, $b^{\mu}\equiv(4\pi)^{-1/2}F^{*\mu\nu}u_{\nu}$ is
the magnetic field vector in the frame comoving with the fluid, $b^{2}/2\equiv
b^{\mu}b_{\mu}/2$ is the energy density of the magnetic field, and
$g_{\mu\nu}$ is the space-time metric. To close the equations, we assume a
three-parameter, finite-temperature, composition-dependent EOS, i.e., we
assume that every dependent thermodynamic variable $\mathcal{P}$ is a function
of mass density, temperature, and electron fraction,
$\mathcal{P}\equiv\mathcal{P}(\rho,T,Y_{\mathrm{e}})$. The electron or proton
fraction is defined by $Y_{\mathrm{e}}=n_{\rm p}/(n_{\rm p}+n_{\rm n})$, where
$n_{\rm p}$ and $n_{\rm n}$ denote the proton and neutron number density,
respectively. In accordance with the tabulated EOS employed here (Sec. 2.4),
we set the baryon mass $m_{\rm b}$ to the free neutron mass for $\mathtt{APR}$
and $\mathtt{LS220}$, and to the atomic mass unit in the case of the
$\mathtt{SFHo}$ EOS (Sec. 2.4).

Using a 3+1 decomposition of spacetime into non-intersecting space-like
hypersurfaces of constant coordinate time $t$ and timelike unit normal
$n^{\mu}$ (Lichnerowicz, 1944; Arnowitt et al., 2008) , we can write Eqs.
(1)–(5) in conservative form as

$\partial_{t}(\sqrt{\gamma}\mathbf{q})+\partial_{i}[\alpha\sqrt{\gamma}\mathbf{f}^{(i)}(\mathbf{p},\mathbf{q})]=\alpha\sqrt{\gamma}\mathbf{s}(\mathbf{p}),$
(8)

where $\mathbf{f}^{(i)}(\mathbf{p},\mathbf{q})$ and $\mathbf{s}(\mathbf{p})$
represent the fluxes and the source terms, respectively (see Eqs. (18) and
(19) in Siegel & Metzger 2018), $\alpha$ denotes the lapse function, $\gamma$
is the determinant of the spatial metric $\gamma_{ij}$ on the hypersurfaces,

$\mathbf{q}\equiv[D,S_{i},\tau,B^{i},DY_{\mathrm{e}}]$ (9)

denotes the vector of conserved (evolved) variables, and

$\mathbf{p}\equiv[\rho,v^{i},\epsilon,B^{i},Y_{\mathrm{e}}]$ (10)

is the vector of primitive (physical) variables. The conserved vector is
composed of the conserved mass density $D=\rho W$, the conserved momenta
$S_{i}=-n_{\mu}T^{\mu}_{\phantom{\mu}i}$, the conserved energy
$\tau=n_{\mu}n_{\nu}T^{\mu\nu}-D$, the three-vector components of the magnetic
field $B^{\mu}=(4\pi)^{-1/2}F^{*\mu\nu}n_{\nu}$, as well as the conserved
electron fraction $DY_{\mathrm{e}}$—as measured by the Eulerian observer,
defined as the normal observer moving with four velocity $n^{\mu}$
perpendicular to the spatial hypersurfaces. This observer measures a Lorentz
factor of the fluid of $W=-u^{\mu}n_{\mu}$, which moves with three-velocity
$v_{i}=u_{i}/W$ on the hypersurfaces.

###  2.2 Neutrino leakage scheme and neutrino transport

The composition of the fluid as traced by $Y_{\mathrm{e}}$, as well as its
internal energy and momentum, can change due to weak interactions (i.e., due
to the emission and absorption of neutrinos). In the evolution equations, the
effect of these interactions is represented by the source terms $Q$ and $R$ of
the right-hand side of Eqs. (1) and (3), which correspond to the net energy
loss/deposition and net lepton number emission/absorption rate per unit
volume, respectively. We use an energy-averaged leakage scheme to treat
neutrino emission at finite optical depth as described in Siegel & Metzger
(2018), closely following the methods of Galeazzi et al. (2013) and Radice et
al. (2016b), which are, in turn, based on Ruffert et al. (1996).

In more detail, we track the reactions involving electron neutrinos,
${\nu_{\mathrm{e}}}$, electron anti-neutrinos, ${\bar{\nu}_{\mathrm{e}}}$, and
the heavy-lepton neutrinos $\nu_{\mu}$ and $\nu_{\tau}$ summarized as
${\nu_{x}}$. In the rest frame of the fluid, the net neutrino heating/cooling
rate per unit volume and the net lepton emission/absorption rate per unit
volume are given as a local balance of absorption and emission of free-
streaming neutrinos:

$Q=\sum_{\nu_{i}}\kappa_{\nu_{i}}n_{\nu_{i}}E_{\nu_{i}}-\sum_{\nu_{i}}Q_{\nu_{i}}^{\mathrm{eff}}$
(11)

and

$R=\sum_{{\nu_{i}}}\kappa_{\nu_{i}}n_{\nu_{i}}-(R_{\nu_{\mathrm{e}}}^{\mathrm{eff}}-R_{\bar{\nu}_{\mathrm{e}}}^{\mathrm{eff}}).$
(12)

Here, ${\nu_{i}}=\\{{\nu_{\mathrm{e}}},{\bar{\nu}_{\mathrm{e}}},{\nu_{x}}\\}$
and $\kappa_{\nu_{i}}$, $n_{\nu_{i}}$, and $E_{\nu_{i}}$ denote the
corresponding absorption opacities, number densities, and mean energies of the
free-streaming neutrinos in the rest frame of the fluid, respectively. The
effective emission rates $R^{\mathrm{eff}}_{\nu_{i}}$ and
$Q^{\mathrm{eff}}_{\nu_{i}}$ are calculated from intrinsic (free) emission
rates $R_{\nu_{i}}$ and $Q_{\nu_{i}}$ by taking into account the finite
optical depth due to the surrounding medium:

$R^{\mathrm{eff}}_{\nu_{i}}=\frac{R_{\nu_{i}}}{1+t_{\mathrm{diff},{\nu_{i}}}/t^{\mathrm{em,R}}_{\nu_{i}}},\mskip
10.0muQ^{\mathrm{eff}}_{\nu_{i}}=\frac{Q_{\nu_{i}}}{1+t_{\mathrm{diff},{\nu_{i}}}/t^{\mathrm{em,Q}}_{\nu_{i}}}.$
(13)

For a given neutrino species ${\nu_{i}}$, we define
$t_{\mathrm{diff},{\nu_{i}}}=D_{\mathrm{diff}}\tau^{2}_{\nu_{i}}k^{-1}_{\nu_{i}}$
as the local diffusion timescale, where $\tau_{\nu_{i}}$ is the optical depth,
and $D_{\mathrm{diff}}\equiv 6$ is the diffusion normalization factor
(O’Connor & Ott, 2010). We define the local neutrino number and energy
emission timescales as

$t^{\mathrm{em},R}_{\nu_{i}}=\frac{R_{\nu_{i}}}{n^{\mathrm{eq}}_{\nu_{i}}},\mskip
20.0mut^{\mathrm{em},Q}_{\nu_{i}}=\frac{Q_{\nu_{i}}}{e_{\nu_{i}}},$ (14)

where $n^{\mathrm{eq}}_{\nu_{i}}$ is the neutrino number density in chemical
equilibrium, and $e_{\nu_{i}}$ denotes the corresponding neutrino energy
densities. The intrinsic emission rates take into account charged-current
$\beta-$processes, electron-position pair annihilation, and plasmon decay (see
Eqs. (25)–(32) in Siegel & Metzger 2018). Neutrino opacities
$\kappa_{\nu_{i}}$ include contributions from absorption of electron and anti-
electron neutrinos and from scattering on heavy nuclei $A$ and on free
nucleons (see Eqs. (33)–(39) in Siegel & Metzger 2018 and Eq. (8) of Radice et
al. 2018). We compute optical depths $\tau_{\nu_{i}}$ using the effective
local approach of Neilsen et al. (2014).

Finally, using a zeroth-moment (M0) scheme to approximate the general-
relativistic Boltzmann transport equation (Radice et al., 2016b), we evolve
the number density $n_{\nu_{i}}$ and average energy $E_{\nu_{i}}$ of free-
streaming neutrinos. In this scheme, only one moment of the neutrino
distribution functions is used, the neutrino number currents
$J^{\mu}_{\nu_{i}}$. Neglecting scattering, the first moment of the Boltzmann
equation yields an equation for lepton number conservation,

$\nabla_{\mu}J^{\mu}_{\nu_{i}}=R_{\nu_{i}}^{\rm
eff}-\kappa_{{\nu_{i}}}n_{\nu_{i}},$ (15)

which is closed by assuming that free-streaming neutrinos propagate on null
radial rays, $J^{\mu}_{\nu_{i}}=n_{\nu_{i}}k^{\mu}$, with $k^{\mu}$ 11a null
four-vector proportional to the radial coordinate of the grid, and normalized
such that $k^{\mu}u_{\mu}=-1$. In a stationary spacetime with time-like
Killing vector $\partial_{t}$ and in the absence of interactions with the
fluid, the mean neutrino energy $\mathcal{E}_{\nu_{i}}\equiv-
p_{\mu}\partial_{t}^{\mu}$ as seen by an observer with four-velocity
$\partial_{t}^{\mu}$ is conserved along the neutrino worldline (parametrized
by the affine parameter $s$),
$d\mathcal{E}_{\nu_{i}}/ds=d(-p_{\mu}\partial_{t}^{\mu})/ds=0$. Here,
$p_{\mu}=E_{\nu_{i}}k_{\mu}$ is the average neutrino four-momentum. In the
presence of heating and cooling, this becomes

$\frac{d\mathcal{E}_{\nu_{i}}}{ds}=\mathcal{Q}_{{\nu_{i}},{\rm
h}}-\mathcal{Q}_{{\nu_{i}},{\rm c}},$ (16)

where $\mathcal{Q}_{{\nu_{i}},{\rm
h}}=-k_{\mu}\partial_{t}^{\mu}(Q^{\mathrm{eff}}_{\nu_{i}}/n_{\nu_{i}})$ is the
local effective average energy deposition rate per neutrino due to neutrino
absorptions and $\mathcal{Q}_{{\nu_{i}},{\rm
c}}=\mathcal{E}_{{\nu_{i}}}(R^{\mathrm{eff}}_{\nu_{i}}/n_{\nu_{i}})$ is the
average energy loss rate per neutrino due to neutrino emission. The assumption
of a stationary spacetime is approximately satisfied once the neutron stars
merge and the spacetime becomes approximately stationary. We switch on this
transport scheme only as the initially cold neutron stars start to merge and
matter is heated up to at least several MeV to copiously produce neutrinos.
The neutrino number densities $n_{\nu_{i}}$ and mean energies $E_{\nu_{i}}$
are obtained by evolving Eqs. (15) and (16) on radial rays (radial worldlines)
using a separate spherical grid that requires interpolation at each timestep
onto the Cartesian grid on which the GRMHD equations are evolved.

This combined neutrino leakage and M0 transport approach include general-
relativistic and Doppler kinematic effects and it is computationally far
cheaper than more complex schemes such as M1 schemes in GRMHD (Foucart, 2018;
Li & Siegel, 2021; Radice et al., 2022) or Monte-Carlo based radiation
transport schemes (Miller et al., 2019a; Foucart et al., 2020) and thus allow
for longer evolution timescales and wider parameter space studies. While
quantitative differences exist between these methods, it is reassuring that
quantitative comparisons between the methods (Radice et al., 2022; Foucart et
al., 2020) reveal remarkable agreement for most physical observables, in
particular with respect to ejecta mass, velocity, and nucleosynthesis.
Deviations in individual physical quantities are typically on the
$\sim\\!10\%$ level (except for $\nu_{x}$ neutrinos, which, however, are not
the focus here), within the error budget of other assumptions and
approximations.

###  2.3 Numerical set-up

We evolve Einstein’s equations coupled to the general-relativistic MHD
equations in conservative form using the open-source framework of the
EinsteinToolkit222http://einsteintoolkit.org (Goodale et al., 2003; Schnetter
et al., 2004; Thornburg, 2004; Loffler et al., 2012; Babiuc-Hamilton et al.,
2019). Regarding the matter part, we use an enhanced version of the public
GRMHD code GRHydro (Baiotti et al., 2003; Mösta et al., 2013) presented in
Siegel & Metzger (2018) and Siegel et al. (2018). We implement a finite-volume
scheme using piecewise parabolic reconstruction (Colella & Woodward, 1984) and
the approximate HLLE Rieman solver (Harten et al., 1987). The magnetic field
is evolved using the FluxCT method (Tóth, 2000) to maintain the solenoidal
constraint during evolution. The recovery of primitive variables is carried
out using the framework presented in Siegel et al. (2018), which provides
support for general finite-temperature, composition-dependent, three-parameter
EOS. We set a static atmosphere at a rest-mass density of $5\%$ above the
minimum density of the tabulated EOS. For most of our models, this is
$\rho_{\rm atmo}\simeq 6\times 10^{2}$ g cm-3, which is sufficiently small to
capture the most tenuous mass outflows from the system. Indeed, at a typical
extraction radius of $r\approx 450$ km, the total mass of the atmosphere
enclosed at that radius is $(4/3)\pi r^{3}\rho_{\rm atmo}\approx 1\times
10^{-7}~{}{M_{\odot}}$, which is an order of magnitude smaller than the total
amount of material in the fast ejecta tails ($v>0.6c$; see below).

Our neutrino leakage and M0 transport scheme is based on the implementation in
WhiskyTHC (Radice & Rezzolla, 2012; Radice et al., 2016b). For the M0
transport, we use a uniform spherical grid of $2048$ radial rays extending up
to $\simeq 765$ km, and a radial resolution of $\Delta r=240$ m. We activate
neutrino heating when the center of mass of the stars is at $\simeq 12$ km,
i.e. just after the NSs touch.

Spacetime is evolved using the BSSNOK formulation (Nakamura et al., 1987;
Shibata & Nakamura, 1995; Baumgarte et al., 1999) as implemented in the
McLachlan code (Brown et al., 2009). We use the $1+\log$ gauge and the
$\Gamma$-driver shift condition with the damping parameter set to $\eta=0.75$.
We use a fourth-order finite-difference scheme with fifth-order Kreiss-Oliger
dissipation and dissipation parameter $\epsilon=0.1$. At the time of black-
hole formation, we start to track the apparent horizon with the methods of the
AHFinderDirect code (Thornburg, 2004) and we set the hydrodynamic variables to
atmosphere floor values well inside the horizon. We also excise the M0
transport grid within a radius close to the apparent horizon and we keep the
origin of the spherical transport grid centered on the black hole and move it
with the black hole if the latter receives a ‘kick’ after merger due to non-
axisymmetric effects.

Time integration of the coupled Einstein and GRMHD equations is performed with
the method of lines, using a third-order Runge-Kutta scheme (Gottlieb et al.,
2009). The Courant-Friedrichs-Lewy (CFL) condition is set using a CFL factor
of 0.25.

In all simulations, we use a fixed Cartesian grid hierarchy composed of six
nested Berger-Oliger refinement boxes, doubling the resolution on each level,
as provided by the Carpet infrastructure (Schnetter et al., 2004). The finest
mesh grid covers a radius of $\simeq 38$ km during the entire evolution and
initially contains both NSs. We use fixed refinement boxes and overlap zones
between levels to keep violations of the magnetic solenoidal constraints under
control (Mösta et al., 2014) and confined to the refinement boundaries without
impacting the dynamics. The outer boundary is placed at a radius of $\simeq
1528$ km. For each EOS (see below), we perform two sets of simulations: $(i)$
a set of high-resolution simulations, setting the finest spatial resolution to
$\Delta x=180$ m and imposing reflection symmetry across the equatorial
(orbital) plane for computational efficiency, and $(ii)$ a set of mid-
resolution simulations with finest spatial resolution of $\Delta x=220$ m
without imposing any symmetries. This allows us to test how ejecta properties
change with resolution and to analyze the influence of reflection symmetry
with regard to the properties of the magnetic field.

###  2.4 Models: EOS and initial data

The properties of BNS mergers can be very sensitive to the choice of the EOS
(e.g., Bauswein et al. 2013b; Radice et al. 2016b, 2018; Hotokezaka et al.
2013b; Nedora et al. 2021b). The compactness of the star, as well as finite-
temperature and composition-dependent effects, impact the lifetime of a merger
remnant, the properties of matter outflows, the emitted gravitational waves,
etc. We use three different nuclear EOSs in our set of simulations as
tabulated by O’Connor & Ott (2010) and Schneider et al. (2017), available at
stellarcollapse.org. These tables tabulate all necessary thermodynamic
quantities as a function of the independent variables
$(\rho,T,Y_{\mathrm{e}})$:

  * •

LS220: We use the liquid-drop model with Skyrme interactions as implemented in
Schneider et al. (2017), using a parametrization corresponding to Lattimer &
Swesty (1991). At low densities, the EOS transitions from a single-nucleus
approximation to nuclear statistical equilibrium (NSE) using 3335 nuclides
(Schneider et al., 2017).

  * •

SFHo: This EOS parametrizes nucleonic matter with a relativistic mean field
model (Steiner et al., 2013). The EOS describes nuclei and nonuniform nuclear
matter using the statistical model of Hempel & Schaffner-Bielich (2010).

  * •

APR: This EOS is based on the potential model of Akmal et al. (1998) built to
fit results from nucleon-nucleon scattering and properties of light nuclei. We
use the new 3D, finite-temperature implementation presented in Schneider et
al. (2019).

All three EOS transition to a nuclear statistical equilibrium including heavy
nuclei in a thermodynamically consistent way, thus including the release of
binding energy as alpha-particles and heavier nuclei form. This recombination
energy of $\sim\\!7$ MeV per baryon can be an important source of energy to
accelerate outflows. This effect is important, in particular, in the context
of post-merger accretion disks (Fernández & Metzger, 2013; Siegel & Metzger,
2017, 2018).

The full 3D tabulated LS220 and SFHo EOSs have been used in general-
relativistic hydrodynamic (GRHD) simulations of BNS mergers (Sekiguchi et al.,
2015; Nedora et al., 2021b; Radice et al., 2018) and GRMHD simulations of the
post-merger phase (Miller et al., 2019a; Li & Siegel, 2021; Mösta et al.,
2020), including weak interactions. APR-type EOSs have been used in both GRHD
and GRMHD simulations but mostly using an ideal-gas approximation for the
finite-temperature part (Ciolfi et al., 2017) or without weak interactions
(Hammond et al., 2021). In this work, we present the first simulations of BNS
mergers and their post-merger evolution with the APR EOS using the full 3D
tabulated EOS and weak interactions.

The choice of EOSs roughly spans the current range of allowed neutron-star
radii of $\approx\\!11-13$ km (De et al., 2018; Miller et al., 2019b; Riley et
al., 2019; Capano et al., 2020; Landry et al., 2020; Dietrich et al., 2020)
for a $1.4M_{\odot}$ neutron star (Fig. 2), thus covering the range of
different allowed compactness scenarios at given mass. The radii of a non-
rotating, cold $1.4M_{\odot}$ neutron star in $\beta$-equilibrium are 11.6,
km, 11.9 km, and 12.7 km for APR, SFHo, and LS220, respectively; the maximum
masses of non-rotating neutron stars for these EOS are 2.2 ${M_{\odot}}$, 2.06
${M_{\odot}}$ and 2.04 ${M_{\odot}}$, respectively. With differences in
neutron-star compactness and maximum mass, these three EOSs give rise to a
range of merger and post-merger phenomena. First, the ‘softer’ APR and SFHo
EOSs are expected to produce more violent collisions, owing to the faster
relative velocities of the stars at merger (Bauswein et al., 2013b), and thus
faster ejecta emanating from the collision interface than in the case of the
‘stiffer’ LS220 EOS. Furthermore, the APR EOS leads to long-lived remnants for
a typical $1.35M_{\odot}-1.35M_{\odot}$ binary (Ciolfi et al., 2019; Ciolfi &
Kalinani, 2020), owing to the relatively high maximum TOV mass, while, for
similar parameters, the SFHo and LS220 EOSs lead to short-lived remnants with
light and heavy post-merger disks, respectively (Nedora et al., 2021b; Mösta
et al., 2020).

Figure 2: Gravitational mass as a function of stellar radius of non-rotating
neutron stars in $\beta$-equilibrium at zero temperature for different EOSs.
Dots mark models with $1.4\,M_{\odot}$ and corresponding radii 11.6,km, 11.9
km, and 12.7 km for APR, SFHo, and LS220, respectively. The grey shaded area
represents mass-radius constraints (90% credible intervals) from current non-
parametric inferences using heavy pulsar measurements, gravitational-wave
observations of GW170817 and GW190425, and NICER data on J0030$+$0451 (Landry
et al., 2020). The blue shaded area indicates the 68% credible mass estimates
from the heaviest pulsars measured so far (Cromartie et al., 2020; Demorest et
al., 2010)

We simulate a set of equal-mass binary neutron star systems on quasi-circular
orbits starting at initial separations of $40$ km. 11At this separation, the
stars inspiral for ${\approx}3$ orbits before they merge. The relatively small
separation introduces a small eccentricity to the quasi-circular orbit,
affecting the GW waveforms; this, however, has a small impact on the
hydrodynamics of the merger, compared to other properties, such as initial
compactness of the stars (which is the focus of the present paper). Individual
stars have ADM masses of $M_{\infty}=1.35{M_{\odot}}$ (at infinite
separation), typical of galactic merging double neutron-star systems (Özel &
Freire, 2016). Our selection of binary models are compatible with the inferred
source parameters of GW170817, with chirp mass of
$M_{c}=1.188^{+0.004}_{-0.002}{M_{\odot}}$ and total gravitational mass of
$2.74^{+0.04}_{-0.01}$ (low-spin prior) and mass ratio of 0.7–1.0 at 90%
confidence (Abbott et al., 2017a). Moreover, our simulations do not result in
prompt collapse upon merger, which is strongly disfavored in the case of
GW170817 by electromagnetic observations of the GW170817 kilonova (e.g.,
Margalit & Metzger 2017; Bauswein et al. 2017; Siegel 2019).

The initial data is built using the open-source code LORENE (Gourgoulhon et
al., 2001)333We use the secant fix suggested in
https://ccrgpages.rit.edu/~jfaber/BNSID/ assuming that the stars are at zero
temperature, $\beta$-equilibrated, and non-rotating. For this construction, we
slice the 3D EOS table assuming the aforementioned conditions, subtracting the
pressure contribution of photons (Radice et al., 2016b), and generate a high-
resolution table as a function of rest-mass density (i.e. an effectively
barotropic table). The initial data is imported into the evolution code
adapting the methods in WhiskyTHC444https://bitbucket.org/FreeTHC.

After setting the hydrodynamical variables, we initialize the magnetic field
as a (dynamically) weak poloidal seed field buried inside the stars. For this
purpose, we specify the vector potential at each star according to

$A_{\phi}=A_{b}\>r^{2}_{\mathrm{cyl}}\left(1-\frac{\rho}{\rho_{\rm
max}}\right)^{n_{p}}\mathrm{max}(p-p_{\rm cut},0)^{n_{s}},$ (17)

where $r_{\mathrm{cyl}}$ is the cylindrical radius relative to the star’s
center, $\rho_{\rm max}$ is the initial maximum rest-mass density, and
$\\{A_{b},n_{p},n_{s},p_{\rm cut}\\}$ are free parameters (Liu et al., 2008).
We choose $n_{p}=0$ to set the maximum strength of the field at the center of
the star, $n_{s}=3$, and $p_{\rm cut}$ is set to $p_{\rm cut}=0.4p_{\rm max}$
in order to keep the field within the stellar interior and well off the
stellar surface. Here, $p_{\rm max}$ is the initial maximum pressure within
the star. Finally, we set $A_{b}$ such that the initial maximum field strength
is $B_{\rm max}=5\times 10^{15}$ G for all BNS models. This high value of the
magnetic field strength is unlikely present in typical BNS, in which surface
values of the poloidal component are expected to be closer to $B_{\rm
pole}\sim 10^{12}$ G, as observed in radio pulsars (Tauris et al., 2017),
while in our case, the strength of the magnetic field near the surface (where
$p=0.4p_{\rm max}$) is $\approx\\!3\times 10^{14}$ G.

Our initial magnetic fields anticipate field strengths similar to those
expected promptly after merger without the need to resolve the computationally
extremely costly amplification process. Magnetic fields are amplified during
the merger process to an equipartition level of $\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}\\!10^{15}-10^{16}$ G
via the Kevin-Helmholtz instability (e.g., Price & Rosswog 2006; Kiuchi et al.
2014; Anderson et al. 2008; Zrake & MacFadyen 2013 and the magnetorotational
instability (e.g., Balbus & Hawley 1991; Duez et al. 2006b; Siegel et al.
2013) on ms timescales. Specifying an initial poloidal seed field is an
acceptable assumption motivated by the fact that the topology is ‘reset’
during merger due to global dynamical effects (see below) as well as by the
aforementioned small-scale amplification effects (Palenzuela et al., 2021;
Aguilera-Miret et al., 2020). Although large at absolute value, the initial
magnetic fields are sufficiently weak to be dynamically insignificant and they
do not alter the inspiral dynamics of the stars.

11Finally, our initial magnetic field configuration neglects the external
magnetosphere and large-scale fields likely present in the inspiral phase of a
BNS system. Simulating regions of extremely high magnetic pressure in ideal
MHD is challenging; various strategies have been pursued, such as imposing
high-density atmospheres (Paschalidis et al., 2015; Ruiz et al., 2016) to
reduce the plasma $\beta$ parameter and to maintain stability, improving the
conservative-to-primitive scheme (Kastaun et al., 2021), or using non-ideal
(resistitive) MHD approximations (Andersson et al., 2022). Our conservative-
to-primitive scheme has already been significantly optimized to handle large
plasma-$\beta$. However, in order to handle magnetospheric field strengths of
$\sim 10^{14}$ G using the first approach, one would still need to increase
the atmospheric density to levels at which it absorbs much of the high-
velocity ejecta. This would thus preclude a detailed analysis of dynamical
ejecta, which is one of the foci of the present paper. While the large-scale
structure of the pre-merger field could have an important role in launching a
short GRB after merger (Mösta et al., 2020; Ruiz et al., 2016), the properties
of the magnetosphere in BNS are still largely unknown. For the purposes of
analyzing the dynamical ejecta, we expect magnetospheric effects to be at
least of second-order importance, compared to other hydrodynamical effects.

###  2.5 Simulation diagnostics and analysis

In the following, we briefly describe some key simulation diagnostic tools.

####  2.5.1 Ejecta properties

This work focuses on dynamical ejecta from the merger process itself and the
electromagnetic radiation it can give rise to. There are various proposed
criteria to determine whether a fluid element becomes unbound from the merging
system. Assuming a fluid element moves along a geodesic in a stationary,
asymptotically flat spacetime, neglecting fluid pressure and self-gravity, a
fluid element can reach spatial infinity if the specific kinetic energy at
spatial infinity is $E_{\infty}:=-u_{t}+1>0$, where $-u_{t}$ is conserved
along the geodesic. The fluid then has an asymptotic escape Lorentz factor
given by

$W_{\infty}=-u_{t}$ (18)

and the escape velocity is $v_{\infty}=\sqrt{1-1/W_{\infty}^{2}}$. This is
known as the geodesic criterion. Because conversion of thermal energy into
kinetic energy by pressure gradients can occur in outflows, one guaranteed
source of internal energy being the recombination of nucleons into alpha-
particles and heavier nuclei, this criterion imposes a lower limit on the
total unbound mass.

Assuming a stationary fluid, a stationary spacetime, and an asymptotic value
of the specific enthalpy of $h_{\infty}$, an unbound fluid element fulfills
the Bernoulli criterion, if the specific energy at infinity satisfies
$-hu_{t}>h_{\infty}$. Typically, the asymptotic enthalpy is $h_{\infty}=1$, if
the EOS is independent of composition (e.g., for a polytropic EOS or ideal
gas). In general, however, the asymptotic enthalpy depends on $Y_{\rm e}$. The
asymptotic value of $Y_{\rm e}$ varies from fluid element to fluid element if
one takes into account neutrino interactions and r-process nucleosynthesis.
The former necessarily plays an important role for most ejecta elements due to
the high ($\sim\\!10$ MeV) temperatures reached during the merger process, and
the latter necessarily sets in as the ejected fluid element decompresses
and/or moves out of nuclear statistical equilibrium (NSE) as it is being
ejected from the merger site. We calculate $h_{\infty}$ taking into account
the contribution of r-process heating as $h_{\infty}=1+\epsilon_{\infty}$,
where $\epsilon_{\infty}$ is the average binding energy of the nuclei formed
by the r-process (see also Foucart et al. 2021; Fujibayashi et al. 2020). This
number depends on how the binding energy is defined in the EOS table, and it
is approximately independent of $Y_{\rm e}$ if we ignore neutrino cooling
during the r-process. For our SFHo table, where the reference mass is the
atomic mass unit, we have $h_{\infty}\simeq 1$, while LS220 and APR use the
free neutron mass, which corresponds to $h_{\infty}\simeq 0.992$. Although
this difference is small, it has a non-negligible effect on the measured
kinetic energy of the ejecta. Using this criterion, the asymptotic Lorentz
factor is:

$W_{\infty}=-\frac{hu_{t}}{h_{\infty}}.$ (19)

We note that the Bernoulli condition might overestimate the total unbound mass
if the system is not stationary (Kastaun & Galeazzi, 2015).

We analyze the physical properties of the ejecta using two methods:
discretizing the outflow on spherical detector surfaces and computing surface
integrals, and by injecting passive tracer particles into the simulation
domain that record plasma properties along their Lagrangian trajectories. We
place different spherical surface detectors at radii $R_{\rm
D}\approx\\{300,450,600,750,900,1050,1350\\}$ km, and extract various
quantities of the fluid on spherical surface grids with resolution
$(n_{\theta},n_{\phi})=(56,96)$. The cumulative ejected mass across a detector
(coordinate) sphere with radius $R_{\rm D}$ is given by

${M_{\rm ej}}=\int dt\oint_{R_{\rm D}}\alpha\sqrt{\gamma}\mathcal{W}\rho
u^{r}dA,$ (20)

where $dA=R_{\rm D}^{2}\sin\theta\,d\theta d\phi$, and
$\mathcal{W}(\theta,\phi)$ is a function that is one (zero) if the fluid at a
given grid point is unbound (bound). We also calculate mass-averaged
properties of the fluid (such as $Y_{\mathrm{e}}$ and specific entropy)
through the detectors as

$\langle\mathcal{P}\rangle\equiv\frac{1}{T}\sum_{i}dt_{i}\frac{1}{M_{i}}\sum_{j}dM_{i,j}\mathcal{P}(j),$
(21)

where $dM_{i,j}$ is the amount of mass passing through the detector during a
time interval $dt_{i}$ for the fluid component with property
$\mathcal{P}\in[\mathcal{P}(j)-\delta,\mathcal{P}(j)+\delta]$, where $\delta$
refers to an appropriately chosen bin width; $M_{i}$ is the total mass
crossing the detector surface in $dt_{i}$, and $T$ is the total time interval
considered. When using surface integration to calculate outflow properties,
one ideally uses large surface radii for extraction to capture the entire
unbound outflow component and to ensure that the fluid is approximately
stationary. This requires, in particular, to evolve the system for
sufficiently long timescales for outflows to cross the relevant detector
spheres. Here, we evolve the systems for more than $30$ ms post-merger, which
is sufficient to yield approximate convergence of the unbound component at
large radii using the geodesic criterion (see Appendix A).

####  2.5.2 Tracers

Ejected material from BNS mergers undergoes rapid neutron-capture
nucleosynthesis and produces heavy elements. We compute nucleosynthesis
abundances arising from the simulation ejecta using a nuclear reaction
network. Specifically, we sample thermodynamic properties such as
$\rho,T,Y_{\mathrm{e}},$ and $W$ using several families of passive tracer
particles injected into the simulation domain and input the recorded
Lagrangian fluid trajectory histories into the nuclear reaction network in a
post-processing step (Sec. 2.5.3).

We place a total of $\sim\\!10^{5}$ tracer particles into the crust of each
star using a probability function proportional to rest-mass density. At the
time of merger, but before ejection of material commences, we allocate $1/4$
of all tracers into a second family of tracers and resample regions of high
specific entropy ($s>30k_{b}$) using the same probability function, to ensure
sufficient sampling of the shock-heated ejecta component that originates in
the collision interface. Finally, after merger, when the system has settled
into a quasi-stationary state ($\sim 20$ ms after merger), we allocate 1/4 of
all tracers into a third family to adequately track the outflows from the
post-merger accretion disks. These tracers are placed within a radius of
$r\in(20,100)$ km and a polar angle of $\theta\in(40,140)$ using the same
probability function as before. This three-stage resampling method allows us
to accurately identify, track, characterize, and distinguish between different
components of outflows, which are ejected by different physical mechanisms.

In order to compute mass distributions and mass-weighted averages of ejecta
properties based on tracers, we (re)assign mass to each tracer once it becomes
part of an outflow. For a set of tracers crossing a sphere of fixed radius
($R_{\rm D}=440$ km) at a given time $t$, we define the mass of those tracers
integrating the mass-flux using the rest-mass density and velocity recorded by
the tracer, $m\equiv\sum_{i}\alpha_{i}\sqrt{\gamma_{i}}\rho_{i}u_{i}^{r}(4\pi
R_{\rm D}^{2})\Delta t$, where $\Delta t$ is the output frequency of the
tracer data (see also Bovard & Rezzolla 2017). We find that this method of
assigning mass very well reproduces the mass distribution of outflows as
extracted by the grid-based spherical surface outflow detectors. In Appendix
B, we present a comparison between the $Y_{\rm e}$ and $v_{\infty}$
distributions as obtained by the grid-based and tracer methods, respectively.

####  2.5.3 Nucleosynthesis setup

Using the full set of unbound tracers, we perform nuclear reaction network
calculations with the open-source reaction network
SkyNet555https://bitbucket.org/jlippuner/skynet (Lippuner & Roberts, 2017) in
post-processing. SkyNet includes 7852 isotopes from free neutrons to 336Cn,
including 140000 nuclear reactions. Nuclear masses (experimental data where
available, data from the finite range droplet macroscopic model (FRDM; Möller
et al. 2016) otherwise), partition functions, and forward strong reaction
rates are taken from the JINA REACLIB database (Cyburt et al., 2010). Inverse
strong reaction rates are derived from detailed balance. Where available, weak
reaction rates are taken from Fuller et al. (1982), Oda et al. (1994),
Langanke & Martınez-Pinedo (2000), and from REACLIB otherwise. Spontaneous and
neutron-induced fission rates are obtained from Frankel & Metropolis (1947),
Mamdouh et al. (2001), Panov et al. (2010), and Wahl (2002).

The Lagrangian trajectories of unbound tracers are followed through the grid
until the end of each simulation or until they exit the grid. At that point,
we smoothly transition and extrapolate their trajectories assuming homologous
expansion ($\rho\propto t^{-3}$). We start evolving the composition with
SkyNet at $t=t_{7\rm GK}$, once the temperature drops below $7$ GK. At this
point, a fluid element represented by a tracer particle is still in NSE, and
SkyNet initially evolves the composition in NSE accordingly until the
temperature drops further to $5$ GK. At $t=t_{5\rm GK}$ NSE starts to break
down and SkyNet automatically switches over to a full network evolution. As a
result the $Y_{\rm e}$ of a fluid element then decouples from the $Y_{\rm e}$
evolution on the simulation grid. Reaction network calculations are performed
up to $10^{9}$ s, when a ‘stable’ final abundance pattern has emerged. Using
the full distribution of tracer particles instead of, e.g., feeding the
outflow history recorded by spherical detectors into the nuclear reaction
network (e.g., Radice et al. 2018; Nedora et al. 2021b) allows us to
distinguish between different ejecta components and to determine the initial
conditions for nucleosynthesis self-consistently for each ejecta component and
fluid element of the flow.

Table 1: Properties of our suite of high-resolution BNS simulations. Here, $M_{1,2}$ denotes the gravitational mass of the neutron stars, $R_{\rm NS}$ their radius, $d_{12}$ is the initial binary separation, and $B_{\rm max}(0)$ is the initial maximum magnetic field strength; $t_{\rm BH}$ is the time of black-hole formation, $M_{\rm ej}$ is the total dynamical ejecta mass, $M_{\rm ej}(v>0.6)$ is the amount of fast moving dynamical ejecta, and ${M_{{\rm ej},n}}$ is the amount of ejecta mass in free neutrons. In brackets, we show several mass and time-averaged ejecta properties, such as the electron fraction ($Y_{\rm e}$), the specific entropy ($s$), the asymptotic velocity ($v_{\infty}$), the angular direction of the ejecta measured from the pole ($\theta$), and the electron fraction extracted at $T=5$ GK by tracer particles ($Y_{\rm e}(T=5\,{\rm GK})$). Simulation | APR | SFHo | LS220
---|---|---|---
$M_{1}$-$M_{2}$ [${M_{\odot}}$] | $1.35-1.35$ |  |
$R_{\rm NS}$ [km] | 11.6 | 11.9 | 12.7
$d_{12}$ [km] | 40 |  |
$|B_{\rm max}(0)|$ [$10^{15}\>$ G] | $5$ |  |
$t_{\rm BH}-t_{\rm merger}$ [ms] | $>50$ | $21$ | $15$
${M_{\rm ej}}\>[10^{-3}{M_{\odot}}]$ | $2.29$ | $3.26$ | $0.92$
${M_{\rm ej}}(v>0.6c)\>[10^{-6}{M_{\odot}}]$ | $4.9$ | $3.8$ | $0.6$
${M_{{\rm ej},n}}\>[10^{-5}{M_{\odot}}]$ | $2.03$ | $1.57$ | $2.44$
$\langle Y_{\mathrm{e}}\rangle$ | $0.25$ | $0.28$ | $0.24$
$\langle s\rangle$ [$k_{\rm b}/\rm b$] | $13.1$ | $13.4$ | $15$
$\langle v_{\infty}\rangle$ [c] | $0.206$ | $0.217$ | $0.195$
$\langle\theta\rangle$ [rad] | $64.8$ | $63$ | $64.6$
$\langle Y_{\mathrm{e}}\>(T=5\textrm{GK})\rangle$ | $0.17$ | $0.24$ | $0.23$

##  3 Simulation results

###  3.1 Overview of simulations

Figure 3: Key dynamical ejecta properties as measured at a radius of 440 km
according to the geodesic criterion: histograms of estimated asymptotic ejecta
velocity ($v_{\infty}$), electron fraction ($Y_{e}$), and specific entropy
($s$), for each high-resolution simulation. We choose the geodesic criterion
here to largely exclude secular ejecta from the remnant NS (not of interest
here) and to focus on dynamical ejecta only. The high-velocity tails of the
ejecta distributions that give rise to free-neutron decay and associated
kilonova precursor emission (Secs. 4.2 and 5.2) are indicated as color-shaded
areas.

In each BNS simulation, the NSs inspiral for approximately three orbits,
emitting GWs, before merging into a massive remnant neutron star. All remnants
created in these runs are gravitationally unstable and eventually collapse
into a black hole upon removal or redistribution of sufficient angular
momentum due to magnetic stresses, spanning a wide range of lifetimes (cf.
Tab. 1). The remnant in the $\mathtt{APR}$ case is a long-lived star with mass
in the supramassive666Supramassive neutron stars refer to stars with mass
above the maximum mass for nonrotating configurations and below the maximum
mass for uniformly rotating configurations. regime and a lifetime of likely
more than a few hundred milliseconds (Ciolfi et al., 2019). In contrast, the
$\mathtt{LS220}$ and $\mathtt{SFHo}$ binaries considered here lead to stars in
the hypermassive777Configurations above the maximum mass for uniformly
rotating neutron stars are referred to as hypermassive neutron stars, which
can be temporarily stabilized against gravitational collapse by differential
rotation. regime with short lifetimes of $\approx\\!15$ ms and $\approx\\!21$
ms, respectively.

Our simulations self-consistently incorporate weak interactions and
approximate neutrino transport, which is pivotal for accurately modeling
ejecta properties, the compositional distribution represented by
$Y_{\mathrm{e}}$, in particular. Furthermore, our simulations include magnetic
fields, which play a key role in generating angular moment transport and
outflows in the post-merger phase. In our setup, magnetic fields are
initialized well inside the stars (cf. Sec. 2.4) and only ‘leak’ out of the
stars during the inspiral in an insignificant way. At merger,
$\beta^{-1}\equiv b^{2}/p$ remains small, and the ‘buried’ fields do not
influence the ejection of dynamical ejecta. In this early stage of the merger
process, our results resemble closely purely hydrodynamic simulations that
include weak interactions and approximate neutrino absorption, but neglect
magnetic fields (e.g., Sekiguchi et al. 2016; Radice et al. 2018).

In this paper, we focus on ejection mechanisms, ejecta properties, and
observables of material ejected during the dynamical phase of the merger
itself. We consider dynamical ejecta only, defined as material that is unbound
by global dynamical processes. Table 1 provides an overview of the mass-
averaged properties of the dynamical ejecta. Corresponding distributions of
ejecta mass relevant for observables according to composition
($Y_{\mathrm{e}}$), asymptotic escape speed ($v_{\infty}$), and specific
entropy ($s_{\infty}$) are summarized in Fig. 3. We extract physical
quantities at a radius of $R=300\,M\simeq 440$ km, where $M$ is the total
binary mass, using the geodesic criterion. We mainly focus on the geodesic
criterion here, since at close separations of 440 km it is somewhat
insensitive to secular outflows such as neutrino-driven winds from the merger
remnant (not of interest for the present study) and it thus acts as a filter
for dynamical ejecta. The merger process during which dynamical ejecta is
generated according to the geodesic criterion lasts approximately $10$ ms in
all our simulations (Sec. 3.2 and Fig. 6). We turn to a discussion of the
details of mass ejection in the following subsections.

Figure 4: Evolution of the electron fraction (upper half-plane) and specific
entropy (lower half-plane) in the equatorial plane during merger for the
$\mathtt{SFHo}$ model. White and gold lines indicate rest-mass density
contours at $\rho=10^{13.5},10^{14.8}\,\rm g\,cm^{-3}$, respectively. After
the first bounce (Sec. 3.2, Fig. 6), two neutron-rich tidal tails form (first
panel). After the second quasi-radial bounce, a shock with associated heating
and neutrino emission reprocesses part of the tidal tail to a higher electron
fraction in a preferential direction defined by the phase of rotation (second
panel). Finally, the remnant’s oscillations are heavily damped, the ambient
flow circularizes, and an accretion disk with $Y_{e}\approx 0.15-0.25$ and
$s\approx 8-10\,k_{B}$ per baryon forms (third panel).

###  3.2 Ejecta dynamics and fast outflow

Two types of ejecta can be distinguished at merger: tidal and shock-heated
ejecta. Tidal torques extract material from the surface of the stars during
the final inspiral and merger process, creating spiral arms that expand into
the orbital plane as they transport angular momentum outwards, expelling cold,
neutron-rich material ($Y_{\mathrm{e}}\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}0.1$) into the
interstellar medium (see Fig. 4, first panel). Because neutron stars are more
compact in general relativity compared to Newtonian gravity, these tidal tails
are not as prominent here as in Newtonian simulations (e.g., Rosswog et al.
1999; Korobkin et al. 2012; Rosswog 2013). Furthermore, for equal-mass
binaries one expects a minimum of tidal ejecta: For a given EOS, tuning the
binary mass ratio away from unity generally enhances the tidal torque on the
lighter companion. This leads to increased tidal ejecta while reducing the
shock-heated component originating in the collision interface (Hotokezaka et
al., 2013a; Bauswein et al., 2013b; Dietrich et al., 2015; Lehner et al.,
2016a; Sekiguchi et al., 2016). This is because the less massive companion
becomes tidally elongated and seeks to ‘avoid’ a (radial) collision. Finally,
for a given binary mass ratio, changing the EOS from stiff (large NS radii) to
soft (small NS radii) one expects the shock-heated component to be enhanced
while reducing the tidal component (Hotokezaka et al., 2013a; Bauswein et al.,
2013b; Dietrich et al., 2015; Lehner et al., 2016a; Sekiguchi et al., 2016;
Palenzuela et al., 2015). This is because tidal forces are smaller for less
extended objects, and NSs with smaller radii approach closer before merger,
reaching higher orbital velocities at the collision, thus enhancing the shock
power and associated ejecta mass.

Figure 5: Total unbound mass flux according to the geodesic criterion through
a coordinate sphere of radius $440$ km for the $\mathtt{APR}$ (upper panel),
$\mathtt{SFHo}$ (middle panel), and $\mathtt{LS220}$ (lower panel) models as a
function of time after the onset of the merger. We distinguish between unbound
ejecta with high specific entropy ($s>20\>k_{\rm B}$ baryon -1), fast
velocities (asymptotic speed $v_{\infty}>0.6c$), as well as between polar
(within $45^{\circ}$ from the polar axis) and equatorial (within $45^{\circ}$
from the equatorial plane) ejecta. We choose the geodesic criterion here to
largely exclude secular wind ejecta from the remnant NS (not of interest here)
and focus on dynamical ejecta only. The total amount of dynamical ejecta
converges around 14 ms at this distance (see also Fig. 22) and the mass flux
drops.

With our NSs spanning the compactness range of currently allowed EOSs for
typical galactic double neutron star masses, we find our runs span a range of
dynamical mass ejection phenomena. A detailed analysis shows (see below) that
for all systems considered here, by far most of the ejecta is expelled by
shock waves produced in quasi-radial bounces of an oscillating double-core
remnant structure that forms after the onset of the merger, with only a
negligible amount of material being ejected by tidal tails (see Fig. 5 and
below; Sec. 3.3). This ejecta material is, in general, faster and more proton-
rich than tidal ejecta. A large fraction of the material released in such
waves has been heated considerably due to hydrodynamical shocks at the
collision interface during the merger process and is further heated as it
shocks into slower surrounding merger debris. Associated neutrino emission in
such a neutron-rich environment favors electron captures onto neutrons over
positron captures on protons, raising the electron fraction and causing the
electron fraction to widen from the initially cold and neutron-rich conditions
of the individual NSs ($Y_{\mathrm{e}}\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}\\!0.1$) into broad
distributions with prominent high-$Y_{\mathrm{e}}$ tails (Fig. 3). The
stronger the heating, as indicated by the specific entropy distributions (Fig.
3), the higher the resulting final $Y_{\mathrm{e}}$ of the ejecta. Neutrino
reabsorption, which also contributes to increasing $Y_{\mathrm{e}}$ in a
neutron-rich environment, only plays a minor role in setting $Y_{\mathrm{e}}$
at this early dynamical stage. Overall, the more compact NS binaries
($\mathtt{APR}$ and $\mathtt{SFHo}$) eject a factor of $\simeq 2-3$ more
ejecta than the less compact $\mathtt{LS220}$ binary (cf. Tab. 1), since tidal
ejecta is subdominant and post-merger bounces are more energetic and lead to
more violent shocks unbinding more material for more compact NSs.

Figure 6 summarizes a more detailed analysis of the matter ejection mechanism
post-merger, correlating the generation of unbound material with radial
oscillations of the remnant. Unbound mass, as sampled by passive tracers
according to the Bernoulli criterion, initially increases in ‘steps’ when the
maximum density reaches a minimum during oscillations.

Figure 6: Upper panel: Fraction of unbound ejected mass sampled by passive
tracer particles according to the Bernoulli criterion as a function of time
since the onset of the merger. Lower panel: evolution of maximum density of
the merger remnant relative to its initial value at merger, showing
perturbations that result from radial oscillations of a double-core structure
within the remnant.

At this point of maximum decompression and least compactness, some amount of
material becomes unbound and is ejected, the remnant undergoes a bounce, and
the maximum density starts to increase again. In the softest EOS simulations,
nearly $75\%$ of the outflow becomes unbound after the first three bounces,
followed by a secular growth that lasts up to 10 ms after merger, reflecting
the integrated action of several weaker bounces. In the stiffer
$\mathtt{LS220}$ case, only two individual mass ejection episodes can be
discerned by this tracer method ($\simeq 30\%$ of the total ejected mass).

The Bernoulli criterion used here is valid for stationary fluid flows only and
might overestimate the total ejecta mass near shocks (see Kastaun & Galeazzi
2015 for a discussion). However, confidence in the present interpretation
comes from the fact that all ejecta curves are monotonically increasing over
time, together with the fact that the amount of unbound material for both the
geodesic and the Bernoulli criteria at large distances agree up to saturation
of the dynamical ejecta as measured by the geodesic criterion (cf. Fig. 22).
The importance of post-merger bounces for mass ejection has been noted
previously in smoothed-particle hydrodynamics (SPH) simulations (Bauswein et
al., 2013b) and more recently in grid-based GRHD simulations (Radice et al.,
2018; Nedora et al., 2021a).

The three (two) individual mass ejection episodes discussed above for
$\mathtt{APR}$ and $\mathtt{SFHo}$ ($\mathtt{LS220}$) drive shock waves into
the surrounding medium (see Fig. 4) and manifest themselves in step-wise
increases of high specific entropy mass flux ($s>20k_{\rm B}$ baryon-1) at
large distances (see Fig. 5).

Figure 5 also indicates that the first high-entropy mass flux events for
$\mathtt{LS220}$ and $\mathtt{APR}$ are associated with predominantly
equatorial outflows (within $45^{\circ}$ from the plane), while the first
event for $\mathtt{SFHo}$ also contains a significant fraction of high-entropy
material in the polar direction (within $45^{\circ}$ from the polar axis). In
general, mass ejection in equatorial directions dominates over polar
directions, as shown by Figs. 5 and 8. However, the relative ratio of
equatorial to polar ejecta is larger for the softer EOSs ($\mathtt{SFHo}$ and
$\mathtt{APR}$), owing to stronger bounces and associated shock waves, cf.
Figs. 5 and 8).

Figure 5 also illustrates that the first two (one) mass ejection episodes
discussed above are (is) associated with the ejection of fast material with
asymptotic speed $v_{\infty}>0.6c$. We find a fast-moving ejecta component
with $v_{\infty}>0.6c$ of mass $\approx\\!(4-5)\times 10^{-6}\,{M_{\odot}}$
for the $\mathtt{APR}$ and $\mathtt{SFHo}$ runs (cf. Tab. 1), which is more
than an order of magnitude larger than the corresponding swept-up atmosphere
mass and is thus well captured; for the stiff $\mathtt{LS220}$ case, we also
find a fast component of mass $\approx\\!6\times 10^{-7}\,{M_{\odot}}$, which
is still larger than the swept-up atmosphere mass by a factor
$\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}7$. Measuring the amount of fast ejecta at
a closer distance of $R_{D}=300$ km, we find that the total mass of this fast
tail is larger by about $5\%$ for both $\mathtt{APR}$ and $\mathtt{SFHo}$ and
by approximately $30\%$ for $\mathtt{LS220}$. The extracted values for the
maximum velocity remain essentially insensitive to the detector radius (see
also Fig. 22 for a convergence test).

This fast outflow component is ejected by the first two bounces after merger,
consistent with other grid-based GRHD simulations (Nedora et al., 2021a). For
$\mathtt{SFHo}$, we observe in Fig. 5 that these fast outflows are ejected in
two impulsive trains of similar strength, while for the even softer
$\mathtt{APR}$ EOS, the first ejection event is much stronger. For the stiffer
$\mathtt{LS220}$ EOS, there is only a single fast ejection event because the
bounces are generally weaker.

The angular distribution of cumulative kinetic energy $E_{\rm
K}=(\Gamma-1)M_{\rm ej}$ of the ejecta as a function of $\beta\Gamma$ is shown
in Fig. 7. Here, $\beta$ denotes the ejecta speed in units of the speed of
light and $\Gamma$ the corresponding Lorentz factor. We observe that the most
energetic outflows for $\mathtt{SFHo}$ are ejected at latitudes close to the
pole (polar angles $\approx\\!0^{\circ}-35^{\circ}$), consistent with previous
findings (Radice et al., 2018). For the softest EOS $\mathtt{APR}$, the
energetics of the ejecta are dominated by near-equatorial outflows (polar
angles $\approx\\!60^{\circ}-90^{\circ}$). For the stiffer EOS
$\mathtt{LS220}$, we obtain milder energetics (mostly $\Gamma\beta<1$), and a
fairly homogenous angular distribution, with more energetic outflows near the
equatorial plane.

Figure 7: Isotropic equivalent cumulative kinetic energy distribution of the
ejecta as a function of ejecta velocity $\beta$ in units of the speed of light
times Lorentz factor $\Gamma$ for six different polar angle bins equidistant
in $\cos\theta$ and different EOS.

The total amount of fast ejecta measured in simulations is sensitive to
different factors such as (a) grid resolution or sampling of high-velocity
ejecta tails with a small number of particles in SPH simulations, (b) the grid
topology, (c) atmosphere levels and magnetosphere treatment, and (d)
definition of asymptotic velocity. In our present study (a) is fixed by
constraints arising from computational cost. High-resolution Newtonian
convergence studies (Dean et al., 2021) suggest that at our resolution of 180
m convergence of the high-velocity ejecta ($v>0.6c$) may, in principle, be
achieved within measurement uncertainties of a factor of $\sim\\!2$. Aspect
(b) is fixed by the Cartesian grid paradigm our code is based on. We elaborate
on aspects (c) and (d) below.

In this work, we use lower floor/atmosphere levels than previously reported
BNS simulations (e.g., smaller by a factor of 10 compared to Radice et al.
2018) to minimize the effects of an artificial atmosphere on the properties of
the fast ejecta (at relevant extraction radii of $\approx\\!400$ km).888We
note that our simulation set-up does not use refluxing through mesh refinement
boundaries; we are thus not ensuring mass conservation at round-off level for
the outflows. However, the other mentioned effects are likely more important
than this one (Reisswig et al., 2013).

The asymptotic velocity $v_{\infty}$ of unbound outflows are computed using
the Lorentz factor at infinity $W_{\infty}=-u_{t}$, where
$W_{\infty}=(1-v^{2}_{\infty})^{-1/2}$. Using the Bernoulli criterion instead,
according to which $W_{\infty}=-hu_{t}$, we find a total amount of fast ejecta
$15\%$ larger, which is expected as thermal effects can still accelerate the
ejecta at this detector distance. There are other techniques used in the
literature to extract outflows, e.g. volume integrations (Kastaun et al.,
2013; Ciolfi et al., 2017; Hotokezaka et al., 2013a), and other definitions of
asymptotic velocity that explicitly include the gravitational potential
(Shibata & Hotokezaka, 2019). Nedora et al. (2021a, b) compute asymptotic
velocities using the Newtonian expression $v_{\rm
N,\infty}=\sqrt{2E_{\infty}}$, where $E_{\infty}$ is the specific energy at
infinity. We find this expression overestimates the speed and amount of fast
ejecta: using the latter expression, we find a larger amount (by almost an
order of magnitude) of fast ejecta in our simulations, comparable to that of
Nedora et al. (2021b).

Previous analyses of fast-ejecta with GRHD simulations using neutrino
transport (Radice et al., 2018; Nedora et al., 2021b) and Newtonian
simulations at ultra-high resolution in axisymmetry (Dean et al., 2021) report
fast-ejecta masses within factors of $5-10$ higher than found here. SPH
simulations find even higher masses than grid-based simulations, of the order
of $10^{-4}-10^{-5}{M_{\odot}}$ (Metzger et al., 2015; Kullmann et al., 2022;
Rosswog et al., 2022). Our default mass estimates are similar to the ones
obtained in the grid-based ultra-high resolution simulations (purely
hydrodynamic, without weak interactions) of Kiuchi et al. (2017),
specifically, their HB EOS case, which leads to a similar NS compactness to
$\mathtt{APR}$ (see Figure 7 in Hotokezaka et al. (2018)). Differences to
other grid-based merger simulations (Nedora et al., 2021b) can be understood
in terms of the definition of asymptotic velocity (see above), while the
higher level of ejecta in grid-based convergence studies (Dean et al., 2021)
may be attributable to differences in the overall setup (such as, e.g.,
Newtonian vs. general-relativistic dynamics and geometric effects due to
assumed symmetries), since ejecta masses are expected to be roughly converged
according to these studies at our present resolution (at least up to a factor
of a few). Fast ejecta tails are resolved by $\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}\\!\text{tens}$ of
SPH particles in present SPH simulations (Kullmann et al., 2022; Rosswog et
al., 2022), and details about the dependence of fast ejecta masses with
resolution (number of particles) have been reported only recently (Rosswog et
al., 2022). Rosswog et al. (2022) show that the mass of fast ejecta decreases
by a factor of $\approx\\!10$ when the resolution is increased from $10^{6}$
to $5\times 10^{6}$ SPH particles used in the simulation. The total fast
ejecta mass in the latter paper—estimated using the relativistic asymptotic
velocity similar to our work—is still a factor of $\approx\\!2$ higher than
the value we find for a BNS merger with a similar set-up and using the APR3
EOS, but much closer than previous SPH simulations.

Introducing a dedicated tracer family to track the shock-heated fast ejecta
(cf. Sec. 2.5.2) allows us to sample the fast ejecta tail with velocities
$v>0.6c$ by typically $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}\\!500$ tracer particles. Appendix B
provides more details on tracer sampling of the fast tail of the ejecta
distribution in our simulations.

###  3.3 Electron fraction distribution and shock reprocessing

The $Y_{\mathrm{e}}$ distribution of the ejecta extracted at $440$ km spans a
wide range of values from $0.05$ to $0.5$ and is centered around
$\approx\\!0.3$ for all simulations, with a more pronounced tail to higher
$Y_{\rm e}$ in the case of the softer $\mathtt{APR}$ and $\mathtt{SFHo}$ EOS
(Fig. 3; Tab. 1). The angular dependence of the mass-weighted $Y_{\mathrm{e}}$
is shown in Fig. 8. Near the pole, we observe high values $>0.3$ of the
electron fraction for all EOSs. This can be attributed mostly to the dominance
of shock-heated ejecta in polar regions as well as to neutrino irradiation of
the polar ejected material by the newly formed remnant. For $\mathtt{APR}$ and
$\mathtt{SFHo}$, the NS remnant reaches higher temperatures due to larger
compressions and oscillations after collision (owing to higher compactness of
the individual stars), whereas the remnant in $\mathtt{LS220}$ is cooler. As a
result, since the heating rate of material due to neutrino irradiation is
roughly the same in all remnants, the entropy of the polar outflows in
$\mathtt{LS220}$ is slightly higher.

Figure 8: Sky map (Moillede) projection of the mass-weighted electron fraction
(north) and of the ejected mass (south) for the total time-integrated
dynamical BNS ejecta measured at $440$ km and the three different EOS
considered here. Equatorial outflows with low-to-medium electron fraction
dominate the ejecta. Shock-reprocessing of tidal ejecta leads to azimuthal
asymmetries in $\langle Y_{\rm e}\rangle$ for soft EOS (‘lanthanide/actinide
pockets’).

Near the equatorial plane, shocks reprocess part of the neutron-rich matter to
values of $Y_{\rm e}\geq 0.1$. For the soft EOSs simulations, we observe a
clear azimuthal asymmetry in the $Y_{\mathrm{e}}$ angular distribution, while
in $\mathtt{LS220}$ case the azimuthal distribution is nearly homogeneous
(Fig. 8). This azimuthal asymmetry for the $\mathtt{APR}$ and $\mathtt{SFHo}$
runs occurs because the first violent bounces of the remnant produce shock
waves in a preferential direction dictated by the phase of orbital rotation.

When the two NSs merge, a rotating remnant with a bar-like ($m=2$ mode) shape
is formed, which undergoes quasi-radial oscillations of its two-core structure
(cf. Sec. 3.2, Fig. 6). The first quasi-radial bounce generates a fast shock-
wave expelling material that propagates outwards without encountering
significant resistance by the ambient medium. The system then develops two
spiral arms due to enhanced tidal torques at maximum decompression, polluting
the environment with cold neutron-rich material (see Figure 4, first panel).
When the remnant bounces again, it generates another violent shock wave
directed along the axis formed by the double-core structure (the $x$-axis in
Figure 4, second panel). This fast shock wave reprocesses the previously
ejected tidal tail material to high $Y_{\rm e}$ through shock-heating and
associated neutrino emission in a relatively wide azimuthal wedge. The
resulting angular distribution shows pockets of low-$Y_{\rm e}$ material
around the equatorial plane (in directions orthogonal to the propagation of
the second shock wave; Fig. 8). These ‘lanthanide/actinide pockets’ possess
much higher opacities than the surrounding shock-reprocessed material (owing
to the lower electron fraction, which is conducive to forming high-opacity
lanthanide and actinide elements). Observable imprints of these pockets in the
kilonova emission can be explored with multi-dimensional kilonova radiation
transport simulations, which are, however, beyond the scope of the present
paper.

Since bounces are much less violent in the $\mathtt{LS220}$ run, azimuthal
inhomogeneities in $Y_{\rm e}$ as described above are much less prominent.
This is because a series of weaker shock waves run into previously (tidally)
ejected material, which leads to an overall more homogeneous shock heating of
the ejecta across all azimuthal directions and thus to a more homogeneous
azimuthal $Y_{\rm e}$ profile. After the first bounces, when most of the
dynamical ejecta is expelled, neutrino-driven outflows and tidal torques help
launch additional material from the remnant, some of which join previously
launched, circularizing debris material in forming a disk with an initial
$Y_{\rm e}\approx 0.15-0.25$ and specific entropy 8-10 $k_{B}$ per baryon (see
cf. Fig. 4, third panel). These entropies that arise from a complicated
process of shock reprocessing as described above are almost identical to the
initial entropies of $8\,k_{B}$ per baryon assumed in previous post-merger
disk simulations (e.g., Fernández & Metzger 2013; Siegel & Metzger 2017).

Figure 9: Snapshots of the electron fraction (upper part) and magnetic field
(lower part) in the meridional plane at different times after merger for the
$\mathtt{APR}$ EOS run. White contours indicate rest-mass densities of
$\rho=\\{10^{9},10^{10},10^{10.5},10^{11},10^{14}\\}\rm~{}g\,cm^{-3}$.

###  3.4 Post-merger phase

In this section, we provide a brief overview of the post-merger evolution for
each simulation. We defer a more detailed discussion of the post-merger stage
to a forthcoming paper.

After merger, the deformed remnant keeps shedding material to the environment
through spiral arms that contribute to forming a disk. The remnant NS cores of
the $\mathtt{LS220}$ and $\mathtt{SFHo}$ binaries keep oscillating until they
merge into a single core and eventually collapse to a BH after $15$ ms and
$25$ ms post-merger, respectively. In contrast, oscillations in the
$\mathtt{APR}$ case damp on a much shorter (few ms) timescale (Fig. 6). The
time of collapse of the remnant to a black hole depends critically on the
angular momentum transport mechanisms present in the system. In this set of
simulations, we incorporate magnetic fields and neutrino transport self-
consistently, which play an important role in the evolution of the angular
momentum of the remnant and its accretion flow.

Angular momentum in the disk is transported outwards by spiral waves, powered
by the $m=2$ azimuthal mode of the remnant, which then develops an $m=1$ mode
through hydrodynamical instabilities (see also Lehner et al. 2016b; Radice et
al. 2016a). The disk increases in mass until $\sim\\!20$ ms post-merger. At
this point, the average magnetic field has increased by two orders of
magnitude from its initial value and MRI-driven turbulence starts contributing
to angular momentum transport in the disk, opposing outward angular momentum
transport by the hydrodynamical spiral waves in the disk. As the MRI fully
develops, the bulk of the disk becomes more turbulent, as can be seen
qualitatively in Fig. 9 (third panel). As a result of the MRI-mediated
angular-momentum transport, the accretion rate increases well above the
ignition threshold for neutrino cooling in the disk (Chen & Beloborodov, 2007;
Metzger et al., 2008a, b; De & Siegel, 2021). Thus cooling becomes
energetically significant, the disk height decreases, and the disk neutronizes
due to electrons becoming degenerate, keeping the disk midplane at mild
electron degeneracy and $Y_{\rm e}\sim 0.1-0.15$ in a self-regulated process
(Chen & Beloborodov, 2007; Siegel & Metzger, 2017, 2018). This can be seen in
Fig. 9 (second and third panel), where the disk $Y_{\rm e}$ evolves from
$\sim\,0.25$ to $\sim\,0.15$, the magnetic field strength increases throughout
the disk, and the disk becomes less ‘puffy’.

Material from the disk winds and merger debris initially pollutes the polar
regions right after merger, suppressing outflows driven by neutrino absorption
(Perego et al., 2014; Desai et al., 2022). After $\sim\\!15-20$ ms post-
merger, a neutrino-driven wind emerges, ejecting only mildly neutron-rich
material (Fig. 9, second and third panel). Such a wind does not develop in the
$\mathtt{LS220}$ case, in which the remnant collapses to a black hole too
early for a wind to develop. In the $\mathtt{SFHo}$ simulation, a coherent
neutrino-driven wind is launched, but the remnant collapses shortly
thereafter. The long-lived remnant formed in the $\mathtt{APR}$ case does
develop a steady neutrino-driven wind, which is additionally aided by magnetic
fields. We defer a more detailed discussion to a forthcoming paper.

##  4 Nucleosynthesis

###  4.1 r-process abundance pattern

Detailed nucleosynthesis analyses of the ejecta using a nuclear reaction
network are performed in a post-processing step (Sec. 2.5.3) based on passive
tracer particles that record various thermodynamic quantities of the flow
along their respective Lagrangian trajectories (Sec. 2.5.2). In each
simulation, $>\\!10\%$ of the initially placed tracers are unbound during the
dynamical phase of the merger, which constitutes approximately $5\times
10^{3}-8\times 10^{3}$ tracers in total.

We find that the injected tracer particles sample the outflow properties very
well compared to those of the mass flux recorded by detector spheres on the
Cartesian grid of the Eulerian observer. Figure 24 shows the ejecta mass
distribution in $Y_{\rm e}$ for the $\mathtt{SFHo}$ run as recorded by a
spherical outflow detector and as recorded by tracers passing through the same
detector, as well as the total number of tracers per $Y_{\rm e}$ bin. The mass
of each tracer is assigned as described in Sec. 2.5.2. Similar good agreement
is obtained in all other simulations and across other quantities, such as
entropy, velocity, etc.; this is of pivotal importance to ensure accurate
nucleosynthesis analyses. Merely the mass distribution in velocity at high
velocities $(v\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}0.6)$ is typically slightly oversampled
(by factors of 2–3) with respect to grid detectors (see Appendix B). This is a
result of placing a large number of tracers in initially highly tenuous plasma
(Sec. 2.5.2). Since only a relatively tiny fraction of mass resides at these
high velocities, such oversampling does not influence the final abundance
patterns (cf. Appendix B); it is beneficial for resolving the properties of
and nucleosynthetic processes within the fast outflows.

Figure 10:  Top: Mass-weighted final nuclear abundances (arbitrary units) as a
function of mass number $A$ for each simulation, based on all tracers sampling
the dynamical ejecta of the respective run. A robust 2nd-to-3rd peak r-process
independent of the EOS is obtained. For comparison, solar r-process abundances
from Sneden et al. (2008) are shown as black dots. Bottom: Mass distribution
as a function of $Y_{\rm e}$ extracted at $T=5\textrm{GK}$ using tracer
particles, showing the collective composition of dynamical ejecta at the onset
of the r-process.

The final mass-averaged abundances for the dynamical ejecta of each simulation
are shown in the top panel of Fig. 10, with solar r-process abundances added
for comparison as black dots. The second (mass number $A\sim 130)$ and third
($A\sim 195)$ r-process peak are well-reproduced in all simulations—a robust
2nd-to-3rd peak r-process is obtained irrespective of the EOS. This is
because, as far as dynamical ejecta is concerned, weak interactions involving
both emission and absorption of neutrinos (chiefly via $e^{+}+n\leftrightarrow
p+\bar{\nu}_{e}$ and $e^{-}+p\leftrightarrow n+\nu_{e}$) only have a finite
($\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}$ few ms; Fig. 6) amount of time to
increase the electron fraction of merger material from the cold, highly
neutron-rich conditions $Y_{\rm e}\sim 0.05-0.15$ of the colliding NSs via
direct shock heating at the collision interface, via reheating of neutron-rich
tidal material by shock waves from the oscillating merger remnant (Sec. 3.2),
and via absorption of strong neutrino radiation from the hot remnant NS that
is being formed. Provided a dominant fraction of the ejecta remains at $Y_{\rm
e}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}0.25$ around 5 GK when NSE breaks down, as
satisfied here (cf. Fig. 10, bottom panel), a sufficiently high neutron-to-
seed ratio can be achieved, such that a pile-up of material in the fission
region at freeze-out occurs. The subsequent decay via fission then guarantees,
depending on the fission model, a robust r-process pattern in the region
$A\approx 120-180$ (Mendoza-Temis et al., 2015). The robustness of the 2nd-
to-3rd peak r-process abundance pattern largely independent of the EOS is in
agreement with other recent BNS simulations including weak interactions and
approximate neutrino transport (Radice et al., 2018; Kullmann et al., 2022).
The robustness of the pattern also holds when extended to non-equal mass
mergers, which suppress the amount of high-$Y_{\rm e}$ material and increase
the neutron-rich tidal ejecta component (and thus the nucleon-to-seed ratio).
For all runs, we find actinide abundances at the level of uranium similar to
solar abundances.

Some deviations from the solar abundance pattern visible in Fig. 10 are likely
the result of nuclear input data for the nucleosynthesis calculations. Since
we employ the FRDM mass model for nucleosynthesis, the third r-process peak is
systematically shifted to the right (slightly higher mass numbers) for all
simulations, which has been attributed to neutron captures after freeze-out
combined with relatively slower $\beta$-decays of third-peak nuclei in the
FRDM model (Eichler et al., 2015; Mendoza-Temis et al., 2015; Caballero et
al., 2014). The trough in abundances between $A\sim 140-170$ relative to solar
is likely due to the fission fragment distribution employed here, as pointed
out by r-process sensitivity studies (e.g., Eichler et al. 2015; Mendoza-Temis
et al. 2015).

The first r-process peak is under-produced in all models, which is due to only
partially reprocessed ejecta material with respect to the original cold,
neutron-rich matter ($Y_{\rm e}\sim 0.05-0.1$) of the individual NSs (see
above). Reproducing the first solar r-process abundances requires a $Y_{\rm
e}$-distribution that extends well above 0.25 for a significant fraction of
the ejecta (e.g., Lippuner & Roberts 2015). Given that the mean of the $Y_{\rm
e}$ distribution is slightly higher for $\mathtt{SFHo}$ and $\mathtt{APR}$,
these models have a larger fraction of first peak material, which is however
still under-produced with respect to solar values. Light r-process elements in
the first-to-second r-process peak region are preferentially synthesized in
the post-merger phase in winds launched from a remnant neutron star and from
the accretion disk around the remnant, which can give rise to broad $Y_{\rm
e}$-distributions (e.g., Perego et al. 2014; Lippuner et al. 2017; Siegel &
Metzger 2018; De & Siegel 2021). We defer a more detailed discussion on
nucleosynthesis including post-merger ejecta to a forthcoming paper.

Figure 11:  Evolution of the specific heating rate calculated with SkyNet for
all unbound tracers of the dynamical ejecta and $\mathtt{APR}$ (upper panel;
color-coded by asymptotic velocity) as well as $\mathtt{LS220}$ (lower panel;
color-coded by electron fraction). Thick lines represent mass averages over
all tracers, which closely follow the expected $\propto t^{-1.3}$ power-law
for r-process heating. Black dashed-lines correspond to an analytic
approximation to the heating rate of free-neutron decay (Kulkarni, 2005).
Heating due to free neutron decay is present over a wide range of EOS
softness/stiffness and originates in high-velocity
($v_{\infty}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}0.5c$) and low-to-moderate $Y_{\rm
e}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}0.2$ outflows.

###  4.2 Fast ejecta and free neutrons

The possibility of tracing individual fluid elements allows us to evaluate the
radioactive heating rate within the ejecta in more detail. In particular, we
are interested in the fast portion of the ejecta that generates free neutrons,
which then $\beta$-decay with a half-life of $\approx\\!10$ min and provide
additional heating of the material at timescales of up to hours relative to
what would be expected from pure r-process heating. Such early excess heating
can power bright UV emission known as a kilonova precursor (Metzger et al.,
2015). The amount of free neutrons generated by the ejecta sensitively depends
on the expansion timescale (i.e., on the ejecta velocity) and the proton
fraction $Y_{\rm e}$.

Figure 11 shows the specific heating rate as recorded by each unbound tracer
that samples the dynamical ejecta. Irrespective of the stiffness/softness of
the EOS, we find a fast $v_{\infty}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}0.5c$ and neutron-rich $Y_{\rm
e}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}0.2$ component of tracers that generate
excess heating on a $\sim\\!10$ min timescale (as expected for free-neutron
decay with half-life $\approx\\!10$ min), relative to the standard
$\dot{Q}\propto t^{-1.3}$ heating rate of a pure r-process. This excess self-
consistently obtained from nuclear reaction network calculations is in good
agreement with analytic predictions for free-neutron decay (dashed lines in
Fig. 11; Kulkarni 2005). Although excess heating from this light, fast ejecta
component does not significantly alter the mass-averaged heating rate (thick
solid lines in Fig. 11) with respect to its power-law behavior, this excess
heating occurs in the outermost layers of the ejecta where even a small number
of free neutrons can give rise to observable emission on timescales of hours
(Sec. 5.2).

Figure 12 shows the detailed mass distributions of free neutrons according to
the asymptotic speed of ejecta. We find that almost all material with
velocities faster than $\approx\\!0.6c$ produces free neutrons. Slower parts
of the ejecta, however, also contribute. These slower layers reside deeper
within the ejecta and thus contribute to the kilonova precursor at a higher
initial optical depth. Their contribution to the total luminosity is thus
dimmer and peaks at slightly later times relative to injection in the
outermost (fastest) layers (Sec. 5.2). The total mass of free neutrons is
$\approx\\!2\times 10^{-5}\,{M_{\odot}}$ in all runs (Tab. 1).

Figure 12: Distribution of dynamical ejecta mass ${M_{{\rm ej},n}}$ (cf. Tab.
1) that results in free neutrons as a function of asymptotic speed and as
sampled by tracer particles (solid lines). Also shown are the corresponding
distributions of total dynamical ejecta mass (dashed lines).

We find that the stiffer EOS simulation $\mathtt{LS220}$ provides a slightly
higher heating rate at $t\approx 10^{-2}$ days than the softer EOSs employed
here (Fig. 13), which we mostly attribute to free neutron heating. Although
the fraction of fast ejecta producing free neutrons in the former model is
smaller (and slower) than those of the latter (Tab. 1, Fig. 12), it is more
neutron-rich owing to less shock heating ($Y_{\rm e}\sim 0.2$ for the fast
ejecta in $\mathtt{LS220}$, compared to $Y_{\rm e}\sim 0.3$ for
$\mathtt{SFHo}$ and $\mathtt{APR}$, see Fig. 11). Overall, this results in
more free neutrons nevertheless. This conclusion is expected from parametric
explorations of r-process nucleosynthesis (Lippuner et al., 2017) (see their
Fig. 3) and also supported by recent parametric studies in merger ejecta using
SPH simulations (Kullmann et al. 2022, their Fig. 12).

Figure 13:  Mass-weighted specific heating rate of all tracer particles
sampling the dynamical ejecta for each simulation. The inset at $t=10^{-2}$
days = $900$ s shows the differences of the mass-weighted heating rates due to
free neutron decay.

##  5 Electromagnetic signatures

Although dynamical ejecta are expected to be subdominant with respect to post-
merger ejecta for typical BNS systems such as those considered here, but also
more generally looking at the expected BNS population as a whole (Siegel,
2019, 2022), they give rise to a number of distinctive electromagnetic
transients that we intend to model here in a self-consistent way based on ab-
initio GRMHD simulations. Based on characteristics of dynamical ejecta
extracted from our GRMHD simulations (Sec. 3) and detailed nucleosynthesis
analysis (Sec. 4), we present a model to calculate the corresponding early
kilonova lightcurves taking into account the contribution of free-neutron
decay (the kilonova precursor), as well as relativistic effects associated
with the fast-moving outflow (Sec. 5.1). We also present a model to calculate
the non-thermal afterglow emission related to the mildly relativistic shock
wave that the fast ejecta gives rise to as it expands into the interstellar
medium (Sec. 5.5). We build a multi-angle model that solves a 1D hydrodynamic
shock-wave problem in multiple directions (Sec. 5.5.1) and computes the
associated radiation of a non-thermal distribution of electrons accelerated
within the shock, solving the time-dependent Fokker-Planck equation (Sec.
5.5.2). Results of both models based on our simulations runs and in the
context of GW170817 are discussed in Secs. 5.2 and 5.5.3, respectively. We
also investigate the potential impact of recently calculated extremely large
opacities of highly ionized lanthanides on the early kilonova emission (Sec.
5.3) as well as signatures of high-energy gamma-ray emission due to inverse
Compton and synchrotron self-Compton processes in the kilonova afterglow shock
(Sec. 5.5.4).

###  5.1 Kilonova Semi-Analytical Model

We calculate kilonova light curves generated by the dynamical ejecta of our
simulations using a new semi-analytical model. We formulate this model based
on the approach of Hotokezaka & Nakar (2020), which, in turn, closely follows
the work of Waxman et al. (2019) and Kasen & Barnes (2019). This model allows
us to compute the heating rates of r-process elements using experimentally
measured values for the injection energies of nuclear decay chains. Our code
builds on the public version of Hotokezaka & Nakar
(2020)999https://github.com/hotokezaka/HeatingRate, which we optimized in
terms of computational efficiency. We further developed the framework to
additionally include (a) arbitrary mass distributions from numerical
simulations, (b) heating rates due to free neutrons, (c) arbitrary opacities,
(d) flux calculation in different wavelength bands, and (c) relativistic
effects. Furthermore, we embed this 1D framework in a generalized
infrastructure for computing a multi-angle kilonova signal following Perego et
al. (2017). For this work, however, we restrict ourselves to a spherically
symmetric (angular-averaged) approximation and defer a multi-angle approach
including post-merger ejecta components to forthcoming work. In the following,
we briefly outline the main features of the model.

Given an angular-averaged mass distribution $M(v)$ extracted from a BNS merger
simulation, we assume homologous expansion of the outflow and we discretize it
into velocity shells $v_{i}$. Each velocity shell is characterized by a mass
$M(v_{i})$ and grey opacity $\kappa(v_{i})$. We typically use at least 30
velocity bins, which we find is sufficient to obtain converged light curves.

The energy of each ejecta shell is determined by losses due to pdV work
associated with the adiabatic expansion of the fluid ($E_{i}/t$), radiative
loses due to photons escaping the flow ($L_{\rm rad,i}$), and heating from
nuclear decay products that thermalize in the ejecta ($\dot{Q}_{i}$). From
energy conservation, the equation for each shell in the comoving frame reads

$\frac{dE_{i}}{dt}=-\frac{E_{i}}{t}+\dot{Q}_{i}-L_{\rm rad,i},$ (22)

which we integrate using a fourth-order Runge-Kutta algorithm. This requires
the computation of the heating term $\dot{Q}_{i}$ and of the radiative
luminosity $L_{{\rm rad},i}$

The ejected outflow contains a variety of heavy neutron-rich isotopes
synthesized by the r-process as well as free neutrons (Sec. 4). Unstable
lighter nuclei will mainly $\beta$-decay, while translead nuclei will also
generate $\alpha$-decays and undergo spontaneous fission. The products of
these processes, mainly fast-moving electrons, $\alpha$-particles, and gamma-
rays, will in part be absorbed by the outflow. The deposition of energy from
charged particles in radioactive decay occurs through collisional ionization
and Coulomb collisions. It is determined by the energy-loss rate of fast
particles $\beta K_{\rm st}$, where $K_{\rm st}$ is the stopping cross-section
per unit mass in units of MeV cm2 g-1 and $\beta$ is the velocity of the
particle in untis of the speed of light $c$. As in Hotokezaka & Nakar (2020),
we take the experimentally measured values of the stopping cross-sections from
ESTAR and ASTAR of the NIST database101010http://physics.nist.gov/Star.

If the timescales of energy loss of the charged particles created in
radioactive decays are much shorter than the dynamical time of the flow, their
energy is entirely deposited in the ejecta. At later times, thermalization is
not as efficient anymore, and we need to follow the energy loss of the
particles. Let us consider a nuclear element $j$ that injects an electron with
energy $e_{0}$ at time $t_{0}$. Taking into account adiabatic losses, the
energy $e_{j}$ of the electron follows (Kasen & Barnes, 2019)

$\frac{de_{j}}{dt}=-3(\gamma(p)-1)\frac{e_{j}}{t}-\beta K_{\rm
st}(e_{j})\rho(t)c,$ (23)

where $\gamma(p)$ is the adiabatic index of charged particles (Nakar, 2020)
and $\rho$ is the density of a mass shell. Once we have solved for
$e_{j}(t;e_{0},t_{0})$, the deposited energy into the mass shell is obtained
by integrating the ionization losses as a function of time, $\beta K_{\rm
st}(t;e_{0},t_{0})$. In order to calculate the total heating rate of charged
decay products at time $t$, one also needs to take into account particles
injected at earlier times as they may not have previously lost all of their
energy. Summing over all nuclear species, the total heating rate per unit mass
is

$\dot{q}_{\rm th,e}(t)=\sum_{j}\int^{t}_{t0_{j}}\beta K_{\rm
st}(t;e_{0},t^{\prime})\rho(t)c\frac{N_{j}(t^{\prime})}{\tau_{j}}dt^{\prime},$
(24)

where $N_{j}(t^{\prime})$ is the number of a radioactive element $j$ per unit
mass, $\tau_{j}$ the mean lifetime of that element, and $t0_{j}$ the injection
time of the oldest non-thermal particle surviving at $t$. To compute the
initial injection energies and mean lifetimes we use the Evaluated Nuclear
Data File library (ENDF/B-VII.1; Chadwick et al. 2011). We use abundance
distributions of nuclei with $A\geq 70$ directly extracted from our
nucleosynthesis calculations (Sec. 4; Fig. 10). We compute Eq. (24) separately
for $\beta$-decay, $\alpha$-decay, and fission, using the corresponding
stopping cross-sections for each process.

Thermalization of gamma-rays produced in radioactive decays occurs through
Compton scattering, absorption, and pair creation. The fraction of energy that
is deposited by gamma-rays is calculated as $f_{\gamma}=1-\exp(-\tau_{\rm
eff})$, where the effective optical depth is defined as $\tau_{\rm
eff}=\kappa_{\rm eff}\Sigma(t)$, with $\kappa_{\rm eff}$ the absorptive
opacity, and $\Sigma(t)$ the mass weighted column density of the mass shell
(see Sec. 3.2 in Hotokezaka & Nakar 2020 for more details). The specific
heating rate is then $\dot{q}_{\rm th,\gamma}=f_{\gamma}(t)\dot{q}_{\gamma}$,
where, again, we use data from ENDF/B-VII.1 to calculate $\dot{q}_{\gamma}$.

Summing the contributions to the heating rates from $\alpha$ and
$\beta$-decays as well and from fission by both charged particles and gamma-
rays, we obtain $\dot{q}_{\rm r}$, the total net heating rate per unit mass of
the synthesized r-process elements. Upon adding the contribution of free
neutrons that undergo $\beta$-decay, the total heating rate of each shell
$\dot{Q}_{i}$ in Eq. (22) is obtained as

$\dot{Q}_{i}=M_{i}(1-X_{\rm fn,i})\,\dot{q}_{\rm r}(t)+M_{i}X_{\rm
fn,i}\,\dot{q}_{\rm fn}(t),$ (25)

where $X_{\rm fn,i}$ is the mass fraction of free neutrons in each velocity
shell as extracted from the nuclear reaction network calculations based on the
simulation runs; for the specific heating rate $\dot{q}_{\rm fn}(t)$ of free
neutrons we set (Kulkarni, 2005)

$\dot{q}_{\rm fn}(t)=3.2\times 10^{14}{\rm
erg\,g^{-1}\,s^{-1}}\,\exp{(-t/\tau_{\rm N})},$ (26)

which provides a good approximation to the free-neutron decay heating rate in
our reaction network calculations (Sec. 4, Fig. 11). Here, $\tau_{\rm N}\simeq
880$ s is the mean lifetime of a free neutron.

Having calculated the effective heating rate for all mass shells, we need to
estimate how photons escape from each shell. In general, the frequency-
integrated (grey) opacity of r-process ejecta sensitively depends on the mass
fraction of lanthanides and actinides (owing to their valence f-shells, which
result in an increase in line expansion opacities by orders of magnitude with
respect to light r-process elements) as well as on the temperature of the
outflow. Since we focus on the early light curve evolution from the outermost
(dynamical) ejecta, the photospheric temperatures satisfy $T>4000-5000$ K, at
which opacities are largely insensitive to temperature (e.g., Kasen et al.
2013; Tanaka et al. 2020; Banerjee et al. 2020). Recent opacity calculations
by Banerjee et al. (2022) find a strong enhancement in the opacity of
lanthanide elements when highly ionized at temperatures $\approx 8\times
10^{4}\,\text{K}$. We discuss in Sec. 5.3 that this temperature regime is not
of concern for the ejecta structure of our binary simulations presented here,
as the photosphere never reaches this temperature regime. The opacity is
largely controlled by the lanthanide content, which, in our scenario, mainly
depends on the electron fraction of the ejecta. We thus ignore the temperature
dependence and use a 1D-parametrization of the opacity as a function of
$Y_{\rm e}$ proposed by Wu et al. (2022),

$\kappa(Y_{\rm e})=\frac{9}{1+(4Y_{\rm e})^{12}}\quad\rm cm^{2}\rm g^{-1},$
(27)

which models the expected range of grey opacity, with a characteristic rapid
drop around $Y_{\rm e}\approx 0.25$ (above which lanthanide production is
suppressed; Lippuner & Roberts 2015). For each velocity shell $v_{i}$, we
calculate a mass-averaged electron fraction $Y_{\rm e}$ to obtain
$\kappa(v_{i})$ according to Eq. (27). Although there is some (residual)
mixing of material (and thus of different $Y_{\rm e}$ components) within and
in-between neighboring velocity shells, we find this mapping is robust when
applied at different extraction radii.

The radiative luminosity of each shell, $L_{\rm rad,i}$, is calculated by
taking into account, approximately, the trapping, diffusive, and free-
streaming radiative regimes. Generalizing the single-zone model of Piro &
Nakar (2013), each velocity shell is assigned an energy escape fraction as

$L_{\rm rad,i}=\frac{f_{\rm Esc,i}E_{i}}{t_{\rm Esc,i}},$ (28)

where $f_{\textrm{Esc},i}=\textrm{erfc}(\sqrt{t_{\rm diff,i}/t})$ is the
energy escape fraction from each shell. Here, erfc denotes the complementary
error function, $t_{\textrm{diff},i}=\tau_{i}v_{i}t/c$ is the diffusion time,
and the escape time $t_{\textrm{Esc},i}$ of the shell is defined by
$t_{\textrm{Esc},i}=\text{min}(t_{\rm diff},t)+v_{i}t/c$ (see also Eqs.
(30)–(32) in Hotokezaka & Nakar 2020 for more details). The optical depth is
determined by $\tau_{i}(t)=\int^{\infty}_{v_{i}t}\kappa(r)\rho(r)dr$, where
$\rho(r)$ is the density of the shell. The radiative luminosity is assigned a
blackbody spectrum $L_{{\rm rad},i,\nu}$, defined by the effective temperature
that is determined by the shell’s radius, $T_{{\rm eff},i}=[L_{{\rm
rad},i}/(4\pi\sigma_{\rm SB}v_{i}^{2}t^{2})]^{1/4}$, where $\sigma_{\rm SB}$
is the Stefan-Boltzmann constant. The total luminosity in the comoving frame
is obtained by adding the radiative luminosity of all shells, $L_{\rm
bol}=\sum_{i}\int\mathrm{d}\nu\,L_{{\rm rad},i,\nu}$. For reference, we also
define the location of the photosphere as the shell $i$ with radius $R_{\rm
ph}$ where $\tau_{i}=1$. The corresponding photospheric temperature is then
defined as $T_{\textrm{ph}}\equiv[L_{\rm bol}/(4\pi\sigma_{\rm SB}R_{\rm
ph}^{2}))]^{1/4}$.

Finally, since we are dealing with mildly relativistic outflows, we need to
incorporate relativistic effects. Since the observed luminosity strongly
scales with the Doppler factor, even mildly relativistic ejecta shell
velocities can have an appreciable impact on the observed lightcurves.
Luminosity and fluxes in the observer frame are computed taking into account
all relevant relativistic effects (the relativistic Doppler effect, the time-
of-flight effect, and relativistic beaming) following the reconstruction
method based on an ‘energy package’ approach presented by Siegel & Ciolfi
(2016). The photosphere at a given time $t$ is discretized in angular bins and
the emergent radiation from each patch of the surface is further discretized
into frequency bins. The resulting energy packages $\Delta
E(\theta_{k},\nu_{l})$ are emitted from the instantaneous photosphere in the
comoving frame and are received by the observer at their respective arrival
times. The observer light curve is reconstructed by binning these packages
into observer time and frequency bins accounting for relativistic effects (see
Sec. 5.7 in Siegel & Ciolfi 2016 for details).

###  5.2 Kilonova Light Curves

####  5.2.1 Bolometric light curves

Figure 14 presents the bolometric luminosity of the kilonova emission and
neutron precursor emission resulting from the dynamical ejecta of all our runs
as computed with the model discussed in Sec. 5.1. Shown are separate
calculations that distinguish the effects of heating due to free neutron decay
and relativistic effects from pure r-process heating only.

Heating from free-neutron decay enhances the total luminosity at early times
($\sim\\!0.5$ h after merger) by a factor of 8–10. Differences to Metzger et
al. (2015), who find an enhancement of the order of 15–20, are likely due to
their mass of free neutron material being an order of magnitude larger and
distributed at larger velocities ($v\sim 0.8c$) than what we find in our
present simulations (see the discussion in Sec. 3.2). We expect relativistic
effects (chiefly Doppler boosting) to enhance the observed luminosity roughly
by the Doppler factor to the third power (cf. Siegel & Ciolfi 2016).111111This
follows from the fact that $I/\nu^{3}$ is a Lorentz invariant, where $I$ is
the specific intensity and $\nu$ the frequency of emitted radiation. Since,
$\nu^{\prime}=D\nu$ in another inertial frame, where $D$ is the Doppler
factor, one has $I^{\prime}=ID^{3}$. We crudely estimate the angular averaged
enhancement factor to

$D^{3}\approx\left(\frac{\sqrt{1-\beta^{2}}}{1-\beta\langle\cos\theta\rangle}\right)^{4}\approx
3.5,$ (29)

where we used a typical (solid angle averaged) value of
$\langle\cos\theta\rangle=\pi/4$. This estimate is consistent with Fig. 14
which shows an additional enhancement by a factor of 3–4. In total, at early
times ($\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}1\,\text{h}$), we find the bolometric
luminosity of this kilonova with an amount of
$\approx\\!2\times\\!10^{-5}{M_{\odot}}$ of free-neutron producing material
and typical velocity of $\langle v\rangle\sim 0.6c$ (cf. Sec. 4.2; Tab. 1) is
a factor $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}\\!10-20$ higher than a non-relativistic
kilonova powered by freshly synthesized r-process elements only. Relativistic
effects subside as the photosphere recedes to lower expansion speeds further
within the ejecta; we find noticeable enhancements in observed luminosity up
to $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}10$ h.

The neutron precursor signal is expected to peak in the source frame when the
free-neutron producing ejecta becomes transparent to photons, $t_{\rm
diff}\approx t$, which, for typical values in our simulations, can be roughly
estimated to

$\displaystyle t_{\rm p}$ $\displaystyle\approx$
$\displaystyle\left(\frac{3{M_{{\rm ej},n}}\kappa}{4\pi\alpha
vc}\right)^{1/2}$ $\displaystyle\approx$ $\displaystyle
0.4\,\text{h}\left(\frac{{M_{{\rm ej},n}}}{2\times
10^{-5}M_{\odot}}\right)^{\frac{1}{2}}\mskip-5.0mu\left(\frac{\kappa}{1\,\text{cm}^{2}\text{g}^{-1}}\right)^{\frac{1}{2}}\mskip-5.0mu\left(\frac{v}{0.6c}\right)^{-\frac{1}{2}},$

in good agreement with the peak timescales of the dashed-dotted curves in Fig.
14. Here, we have assumed a steep outermost ejecta profile with $M(v)\propto
v^{-\alpha}$, where $\alpha\approx 10$ (cf. Figs. 19 and 12). Relativistic
effects then decrease that timescale roughly to $t_{\rm
p,obs}=\sqrt{1-\beta^{2}}t_{\rm p}\approx 0.8t_{\rm p}\approx 0.3\,h$ for
$v=\beta c=0.6c$, consistent with the shift evident from Fig. 19. The observed
peak timescale is thus closer to the timescale of free neutron heating of
about $0.25\,\text{h}\approx 900$ s (cf. Eq. (26)).

Softer EOSs show a brighter transient at peak light with a slightly earlier
peak time, while $\mathtt{LS220}$ shows a somewhat broader peak and slightly
delayed peak time. While this behavior is expected from the fact that softer
EOSs lead to larger velocities of the outermost layers of the dynamical ejecta
(cf. the discussion in Sec. 3.2, Figs. 3, 7), the electron fraction of the
outflows complicates the picture. Part of the fastest ejecta is shock-heated
and thus turns less neutron-rich, with softer EOS giving rise to more shock
heating. As free neutrons become increasingly buried within slower ejecta
shells, their heating acts at higher initial optical depths, which implies
longer peak timescales, hence a broader peak, and a dimmer light curve at peak
light (cf. $\mathtt{LS220}$ in Fig. 14). We note that despite the smaller
ejecta velocities, heating due to free neutron decay can still be significant
for stiffer EOS due to a higher neutron-richness in the fast outflows.

Figure 14: Bolometric light curves of kilonovae and neutron precursor emission
produced by the dynamical ejecta only for each simulation, taking into account
free-neutron decay and relativistic effects (thick lines), free-neutron decay
without relativistic effects (dot-dashed lines), and r-process heating only
(dashed lines).

####  5.2.2 Application to GW170817

Given the high photospheric temperatures, the kilonova precursor emission is
most prominent in the U-V/U band. The binary configurations explored here are
compatible with the binary parameters inferred for GW170817 (Sec. 2). Figure
15 shows the predictions of our kilonova model (Sec. 5.1) for the UV/U/B bands
based on the $\mathtt{SFHo}$ and $\mathtt{LS220}$ runs when applied at the
distance of GW170817 of $\approx\\!41$ Mpc, together with kilonova
observations in these bands starting at 11 h after merger as compiled in
Villar et al. (2017). Our results suggest that if earlier observations could
have been conducted, the neutron precursor signal would likely have been
detected. Based on our modeling, we find that the $\approx 2\times
10^{-5}\,M_{\odot}$ neutron precursor generated by the binary systems typical
of Galactic double neutron-star systems considered here can be detected with
planned missions such as the wide-field ULTRASAT satellite (Sagiv et al.,
2014)121212https://www.weizmann.ac.il/ultrasat/ out to $\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}\\!250$ Mpc, assuming
a $5\sigma$ limiting magnitude of $\approx\\!23$ for a 1 h target-of-
opportunity integration at a wavelength of 220–280 nm.

The double-peaked structure in the UV/U/B bands ($\lambda=268,365,445$ nm,
respectively) evident here is due to the neutron precursor and conventional
r-process heating, respectively. Figure 15 illustrates that the relative
strength of the first and second peak (i.e. neutron precursor to dynamical
r-process ejecta) encodes information about the softness of the EOS, which
might prove useful to place additional EOS constraints in future merger
events. As also evident from Fig. 15 (cf. the short peak timescales), the
dynamical ejecta in our runs is not massive enough to explain the blue
kilonova data in GW170817 (Just et al., 2022). This is consistent with
previous conclusions that the blue kilonova emission in GW170817 likely
requires a substantial contribution from post-merger ejecta (Siegel, 2019;
Metzger, 2020), such as a magnetar-accelerated neutrino-driven wind (Metzger
et al., 2018; Mösta et al., 2020; Curtis et al., 2021).

Figure 15: AB magnitude of blue wavelength bands for the $\mathtt{SFHo}$
(solid lines) and the $\mathtt{LS220}$ model (dashed lines). Also shown as
dots with corresponding colors are early kilonova observations of GW170817 as
compiled in Villar et al. (2017).

###  5.3 Impact of highly ionized lanthanides on early kilonova emission

Recent opacity calculations by Banerjee et al. (2022) point out that high
ionization states of lanthanide elements reached at temperatures
$\approx\\!80000$ K can boost the opacity by a factor of 10–100, up to
$\approx 10^{3}\,\text{cm}^{2}\,\text{g}^{-1}$. If this temperature regime is
realized in the early evolution of merger ejecta $\sim\\!1$ h that contains
only a small fraction of lanthanides, Banerjee et al. (2022) show that the
boosted opacity can significantly dim the early $\sim\\!1\,\text{h}$ kilonova
lightcurve. This affecting, in particular, the early UV brightness is a
potential concern for detection prospects of the neutron precursor emission as
discussed in the previous sections.

Figure 16: Gas temperature of our ejecta structure as a function of velocity
coordinate and time for the $\mathtt{SFHo}$ model. The white line represents
the location of the photosphere, indicating that the local gas temperature
stays below the temperature regime of $\approx 80000$ K for an opacity boost
due to highly ionized states of lanthanides (Banerjee et al., 2022).

Figure 16 shows the gas temperature of our ejecta structure in the
$\mathtt{SFHo}$ case as a function of velocity coordinate as well as time.
Similar behavior applies to our other merger simulation runs. We reconstruct
the gas temperature in each ejecta velocity shell $i$ with mass $dM_{i}$ and
volume $V_{i}$ from the kilonova model (Sec. 5.1) assuming contributions of
radiation and ions to the internal energy $E_{i}$,

$E_{i}=aV_{i}T^{4}+\frac{3}{2}\frac{dM_{i}}{\langle A\rangle m_{u}}k_{\rm
B}T,$ (31)

where $a$ is the radiation constant, $\langle A\rangle$ is the average mass
number as obtained from the nucleosynthesis calculations (Sec. 4), and $m_{u}$
is the atomic mass unit (see also Just et al. (2022)). In agreement with the
model calculations of Banerjee et al. (2022), we find a 70000 K temperature
regime at low velocity coordinates $v\sim 0.2$ c and early times
$\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}\\!1$ h, in which high ionization states
of lanthanides can be realized. However, the photosphere in our merger ejecta
resides at much lower gas temperatures at all times, leading us to conclude
that the opacity boost from highly ionized lanthanides is unlikely to
significantly affect the early kilonova emission in practice.

###  5.4 Impact of “lanthanides pockets” on the kilonova and its nebular phase

11At late times, the kilonova reaches a nebular phase, in which the ejecta
material becomes optically thin and the electromagnetic emission driven by
radioactivity is dominated by emission lines (Hotokezaka et al., 2021). The
spectral signatures of this late-time phase can contain unique information
about atomic species synthesized by neutron capture.

11Our simulations show that the azimuthal distribution of the electron
fraction for soft EOSs is non-uniform because of the preferred direction of
core bounces (see Sec. 3.3 and Fig. 8). The ejected material contains
“pockets” of neutron-rich material that synthesizes lanthanides and actinides.
The opacity of these pockets is thus significantly higher than the surrounding
ejecta, which has undergone stronger heating due to shocks during the
dynamical merger phase. The resulting “opacity pockets” are illustrated in the
upper panel of Figure 17, extracted at a distance of ${\approx}440$ km from
the merger site. For the stiff EOS, these lanthanide pockets are absent and
the opacity distribution is more uniform, see Figure 18.

Figure 17: Upper part: Skymap of mass-weighted opacity for unbound dynamical
ejecta in soft EOS runs (here: $\mathtt{SFHo}$), computed using the $Y_{\rm
e}$ distribution extracted from the simulation at a distance of ${\approx}440$
km from the merger site and using Eq. (27). For soft EOS a latitudinally and
azimuthally dependent opacity profile emerges due to inhomogeneous shock
heating/reprocessing of tidal ejecta (Sec. 3.3), giving rise to
‘lanthanide/actinide pockets’, which might lead to interesting observational
consequences both in the early kilonova lightcurve and the late-time nebular
phase. 11Lower part: Corresponding skymap of mass-weighted quotient of radial
and azimuthal velocities. Figure 18: 11Same as Figure 17, but for the stiff
$\mathtt{LS220}$ EOS simulation. The opacity distribution is more uniform in
the azimuthal direction and shows the expected latitudinal dependence.

11 If the structure of these pockets is preserved to larger distances (late
times), this will have an impact on the observable properties of the nebular
emission since the dominant nebular lines are determined by elements that
dominate the cooling of the plasma. We estimate angular spreading of the
ejecta based on the mass-weighted radial ($\langle v_{r}\rangle$) and
azimuthal ($\langle v_{\phi}\rangle$) velocities extracted at ${\approx}440$
km, as shown in the lower panel of Fig. 17 in the case of $\mathtt{SFHo}$. In
general, we observe that radial velocities in these pockets dominate by at
least a factor of ${\approx}3-5$ over azimuthal velocities, even though this
low-$Y_{\rm e}$ component arises mostly from the tidal tails. For surrounding
material this ratio is larger, typically ${\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}}6-10$. As a result
of the ratio $\langle v_{r}\rangle$/$\langle v_{\phi}\rangle$, material from
within the pockets and from the surroundings will azimuthally spread and tend
to disperse over time. A rough estimate for the angular dispersion of opacity
structures in Fig. 17 is given by

$\Delta\phi\approx\frac{d_{\phi}(t)}{r(t)}\sim\frac{\langle
v_{\phi}\rangle}{\langle v_{r}\rangle}\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}\left(\frac{1}{5}-\frac{1}{3}\right)\approx(10-20)^{\circ},$
(32)

where $d_{\phi}(t)$ and $r(t)$ denote azimuthal and radial distance traveled
by the ejecta, respectively. It is therefore plausible that such structures
may persist at timescales of interest for both early photospheric ($\sim$days)
and late-time ($\sim$weeks) nebular emission.

###  5.5 Kilonova afterglows

As the mildly relativistic dynamical ejecta expands into the interstellar
medium (ISM), it sweeps up ISM material and generates a long-lived blast-wave
(Nakar & Piran, 2011; Piran et al., 2013; Margalit & Piran, 2015; Hajela et
al., 2022; Hotokezaka et al., 2018). Within the shock, randomly oriented
magnetic fields are generated, amplified, and particles, mainly electrons, are
accelerated to non-thermal distributions and generate synchrotron emission
(Sari et al., 1998). In this section, we present a model to calculate the
dynamics of the shock propagation and the generation of non-thermal radiation
directly based on the results of numerical relativity simulations.

####  5.5.1 Shock dynamics

To describe the hydrodynamical propagation of the blast wave the ejecta runs
into the ISM, we assume the ejecta has entered homologous, quasi-spherical
expansion with a mass profile $M(v)$ as, e.g., in Fig. 3. We assume this shock
sweeps up ISM material, which remains concentrated in a thin slab close to the
shock front, where most of the electrons are accelerated. If we assume that
the shock is adiabatic, and the EOS of the fluid is trans-relativistic
(Mignone & McKinney, 2007; Nava et al., 2013; van Eerten, 2013), the evolution
of the outermost ejecta with velocity $\beta$ (in units of $c$) and associated
Lorentz factor $\Gamma$ is determined by energy conservation (see, e.g., van
Eerten 2013; Ryan et al. 2020 for more detailed discussion), and the velocity
$\dot{R}$ of the shock at radial position $R$ from the merger site is
determined by the Rankine-Hugoniot jump conditions,

$\frac{\dot{R}}{c}=\frac{4\Gamma u}{4u^{2}+3},$ (33)

where $u=\Gamma\beta$ is the four-velocity.

In a merger event, the outflow has a radial velocity structure (Fig. 3). As
the fastest part of the outflow starts to decelerate because of its
interaction with the ISM, velocity shells deeper within the ejecta inject
kinetic energy into the shock region, leading to a so-called refreshed shock
(Panaitescu et al., 1998; Rees & Meszaros, 1998). In this case, conservation
of energy is expressed by

$E_{\rm K}(>\\!u)=(\Gamma-1)M_{0}+\frac{4}{9}R^{3}\rho_{\rm
ISM}\frac{u^{2}(4u^{2}+3)}{1+u^{2}},$ (34)

where $E_{\textrm{K}}(>\\!u)$ is the kinetic energy of the flow extracted from
simulations (see Fig. 19), $M_{0}=M(u_{\rm max})$ is the mass of the outer-
most part of the fluid (we impose a mass cut-off of
$M_{0}=10^{-8}{M_{\odot}}$), and the second term is the kinetic plus thermal
energy of the shocked ISM material that has been swept up by the blast wave.
Here, $\rho_{\rm ISM}$ is the constant density of the ISM. The fluid velocity
as a function of shock position $R$ can be obtained from the non-linear
algebraic Eq. (34), and then the shock Lorentz factor $\Gamma_{\rm sh}(R)$ can
be found using Eq. (33). We note that in the non-relativistic limit, Eq. (34)
reduces to the solution used by Piran et al. (2013) and Chevalier (1982).

Figure 19: Cumulative kinetic energy of the dynamical ejecta as a function of
four-velocity $u=\Gamma\beta$ for the three BNS merger simulation runs
considered here.

Employing an ISM density of $n=\rho_{\rm ISM}/m_{p}=0.001\,\rm cm^{-3}$, where
$m_{p}$ is the proton mass, we compare the evolution of the shock Lorentz
factor $\Gamma_{\rm sh}=\left[1-(\dot{R}/c)^{2}\right]^{-1/2}$ between the
different simulation runs as a function of $R$ in Fig. 20. The (forward) shock
propagates at mildly relativistic velocities at early times, with differences
of tens of percent between different EOSs and larger Lorentz factors for
softer EOS. The shock waves become Newtonian ($\Gamma_{\rm sh}\approx 1$) once
reaching their respective Sedov radius, defined as the radius $R_{\rm d}$
within which the rest mass energy of the material within the blast wave equals
the energy $E_{\rm K}$ of the explosion,

$\displaystyle R_{\rm d}$ $\displaystyle\approx$
$\displaystyle\left(\frac{E_{k}}{(4\pi/3)nm_{p}c^{2}}\right)^{1/3}$
$\displaystyle\approx$ $\displaystyle 2.5\times
10^{18}\,\text{cm}\,\left(\frac{E_{\rm
K}}{10^{50}\,\text{erg}}\right)^{1/3}\left(\frac{n}{10^{-3}\,\text{g\,cm}^{-3}}\right)^{-1/3},$

in good agreement with the behavior shown in Fig. 20. Here, we have evaluated
the expression for typical isotropic equivalent energies of the dynamical
ejecta of our simulations (Fig. 19). The shock waves decelerate more gradually
than expected from a single-shell scenario $(\Gamma_{\rm sh}\propto R^{-3/2})$
due to the fact that the shock is continuously being refreshed by kinetic
energy of shells deeper within the ejecta that catch up with the decelerating
shock front.

Figure 20: Evolution of the Lorentz factor $\Gamma_{\rm sh}$ of the shock
driven into the ISM (with density $n=0.001\,\text{cm}^{-3}$) by the dynamical
ejecta as a function of radial distance to the merger site and for the three