# 3D radiative transfer kilonova modelling for binary neutron star merger
simulations

Christine E. Collins,1 Andreas Bauswein,1 Stuart A. Sim,2 Vimal Vijayan,1,3
Gabriel Martínez-Pinedo,1,4 Oliver Just1,5, Luke J. Shingles1 and Markus
Kromer6
1GSI Helmholtzzentrum für Schwerionenforschung, Planckstraße 1, 64291
Darmstadt, Germany
2Astrophysics Research Center, School of Mathematics and Physics, Queen’s
University Belfast, Belfast BT7 1NN, Northern Ireland, UK
3Department of Physics and Astronomy, Ruprecht-Karls-Universität Heidelberg,
Im Neuenheimer feld 226, 69120 Heidelberg, Germany
4Institut für Kernphysik (Theoriezentrum), Fachbereich Physik, Technische
Universität Darmstadt, Schlossgartenstraße 2, 64298 Darmstadt, German
5Astrophysical Big Bang Laboratory, RIKEN Cluster for Pioneering Research, 2-1
Hirosawa, Wako, Saitama 351-0198, Japan
6Heidelberger Institut für Theoretische Studien (HITS), Schloss-
Wolfsbrunnenweg 35, 69118 Heidelberg E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The detection of GW170817 and the accompanying electromagnetic counterpart,
AT2017gfo, have provided an important set of observational constraints for
theoretical models of neutron star mergers, nucleosynthesis, and radiative
transfer for kilonovae. We apply the 3D Monte Carlo radiative transfer code
ARTIS to produce synthetic light curves of the dynamical ejecta from a neutron
star merger, which has been modelled with 3D smooth-particle hydrodynamics
(SPH) and included neutrino interactions. Nucleosynthesis calculations provide
the energy released from radioactive decays of r-process nuclei, and radiation
transport is performed using grey opacities given as functions of the electron
fraction. We present line-of-sight dependent bolometric light curves, and find
the emission along polar lines of sight to be up to a factor of $\sim 2$
brighter than along equatorial lines of sight. Instead of a distinct emission
peak, our bolometric light curve exhibits a monotonic decline, characterised
by a shoulder at the time when the bulk ejecta becomes optically thin. We show
approximate band light curves based on radiation temperatures and compare
these to the observations of AT2017gfo. We find that the rapidly declining
temperatures lead to a blue to red colour evolution similar to that shown by
AT2017gfo. We also investigate the impact of an additional, spherically
symmetric secular ejecta component, and we find that the early light curve
remains nearly unaffected, while after about $1\,$ day the emission is
strongly enhanced and dominated by the secular ejecta, leading to the shift of
the shoulder from 1-2 to 6-10 days.

###### keywords:

Radiative transfer – (Transients:) neutron star mergers – Methods: numerical

††pubyear: 2022††pagerange: 3D radiative transfer kilonova modelling for
binary neutron star merger simulations–3D radiative transfer kilonova
modelling for binary neutron star merger simulations

##  1 Introduction

The detection of GW170817 (Abbott et al., 2017) and its optical counterpart
AT2017gfo (e.g. Smartt et al., 2017, see Villar et al. 2017 and references
therein) confirmed the prediction that a kilonova would accompany the merging
of binary neutron stars (BNS) (Li & Paczyński 1998; Metzger et al. 2010, see
e.g. Metzger 2019 for a recent review). The observations of a power-law
luminosity decline consistent with r-process material and a sufficient event
rate are consistent with binary neutron star mergers being the dominant site
of r-process element production in the Universe (Kasen et al., 2017; Drout et
al., 2017), which had previously been proposed by some authors on theoretical
grounds (Lattimer & Schramm, 1976; Eichler et al., 1989).

Simulations of BNS mergers predict complex ejecta structures. The merger
ejecta is comprised of material expelled on dynamical timescales of tens of
milliseconds, which is followed by secular ejecta expelled on timescales of
seconds. The dynamical ejecta typically have mass of $10^{-4}$ – $10^{-2}$
${\mathrm{M}}_{\odot}$, and high ejecta velocities, typically 0.2 – 0.3c (e.g.
Bauswein et al., 2013; Hotokezaka et al., 2013; Tanaka & Hotokezaka, 2013;
Krüger & Foucart, 2020; Radice et al., 2018; Ardevol-Pulpillo et al., 2019).
The secular ejecta, i.e. the material unbound on longer time scales, are
typically more massive ($10^{-2}$ – $10^{-1}$ ${\mathrm{M}}_{\odot}$) and
ejected with lower velocities of $\sim$0.1c (e.g. Fernández & Metzger, 2013;
Perego et al., 2014; Just et al., 2015; Fujibayashi et al., 2018; Siegel &
Metzger, 2018).

The ejecta are predicted to have low lanthanide fraction material in the polar
directions, where there is a lower neutron abundance, corresponding to
electron fractions of $Y_{e}\gtrsim 0.3$, which would lead to low opacities
and therefore ‘blue’ colours in the predicted kilonova. Higher lanthanide
fractions are expected in the equatorial directions, where there is a higher
neutron abundance ($Y_{e}$ $\lesssim 0.3$) which would lead to higher
opacities and ‘red’ colours in the predicted kilonova (e.g. Metzger &
Fernández, 2014; Sekiguchi et al., 2015; Just et al., 2015; Perego et al.,
2014; Foucart et al., 2020; Radice et al., 2022).

The observations of the kilonova AT2017gfo initially showed blue colours with
the spectra peaking in the UV/blue, which rapidly evolved to redder colours
with spectra peaking in the near-infrared (see e.g. Smartt et al., 2017;
Villar et al., 2017). It has been suggested that to explain this rapid colour
evolution, the outermost layers of the merger ejecta must be composed of low
opacity, high $Y_{e}$ material (which would lead to low lanthanide fraction
material), while the inner ejecta layers require higher opacity, low $Y_{e}$
material (high lanthanide fraction) to reproduce the red colours (e.g.
Cowperthwaite et al., 2017; Villar et al., 2017, based on the analytical two
compnent model introduced by Metzger 2017).

Initial progress in this field has relied on the use of analytic ejecta models
that are described by e.g., power law density structures or idealised
geometries (e.g., Metzger et al., 2010; Barnes & Kasen, 2013; Watson et al.,
2019; Banerjee et al., 2020; Even et al., 2020; Domoto et al., 2021; Heinzel
et al., 2021; Korobkin et al., 2021; Wollaeger et al., 2021; Pognan et al.,
2022). Kilonova studies based on parameterised ejecta configurations avoid the
complexities connected to hydrodynamics modelling, while other aspects, such
as thermalisation (e.g. Barnes et al. 2016) can be studied in more detail.
Ultimately, however, a reliable interpretation of future multi-messenger
observations of BNS will require kilonova models to be based on self-
consistent simulations of the merger and its ejecta.

Of the studies that have considered simulated ejecta from BNS mergers, these
have mostly carried out radiative transfer simulations in 1D (Curtis et al.,
2021; Gillanders et al., 2022; Wu et al., 2022) and 2D (Kasen et al., 2015;
Kawaguchi et al., 2021; Just et al., 2022; Kawaguchi et al., 2022; Klion et
al., 2022), although we note some of these were based on 3D merger simulations
(e.g. Kawaguchi et al., 2021; Just et al., 2022). 3D radiative transfer
calculations have been carried out by Tanaka & Hotokezaka (2013), but for
axisymmetric merger data, and by Neuweiler et al. (2022). Bulla et al. (2021)
and Nativi et al. (2021) have carried out 3D radiative transfer simulations
for the 3D merger simulation of a BH-NS system, and for a 3D simulation of
neutrino-driven winds from a neutron star merger remnant, respectively.

In this paper, we carry out radiative transfer simulations for merger ejecta
extracted directly from a BNS simulation. Our hydrodynamic simulation does not
only self-consistently provide the density distribution of the dynamical
ejecta, but also – thanks to its sophisticated neutrino treatment – the
electron fraction, which determines the r-process composition and therefore
the energy powering the kilonova, as well as the opacities.

##  2 Methods

###  2.1 Hydrodynamical Model

The model we consider in this paper is of the dynamical ejecta produced from
the merger of binary neutron stars. The merger simulation was carried out
using a 3D general relativistic smooth-particle hydrodynamics (SPH) code
(Oechslin et al., 2002; Bauswein et al., 2013), and an advanced neutrino
leakage treatment, ILEAS (Improved leakage-equilibration-absorption scheme
Ardevol-Pulpillo et al., 2019). As a first step, we consider equal mass
neutron stars, each of 1.35 M⊙ (gravitational mass for an orbit with infinite
separation), and each with 150000 SPH particles, using the SFHo equation of
state (Steiner et al., 2013). Hence, the chosen setup represents a system well
compatible with GW170817 considering its mass and equation of state. The stars
were placed on an orbit with a separation of $\sim$38 km and we assumed no
intrinsic spin of the stars. The stellar matter was initially cold and in
neutrinoless beta-equilibrium. The simulation was evolved until 0.02 seconds
after the time when both stars first touched. The calculation thus only covers
the early dynamical mass ejection, while considerable amounts of ejecta may
still be produced on longer time scales. The system did not collapse to a
black hole within the timescale of the simulation. For the given mass and
equation of state a delayed gravitational collapse may be expected.

####  2.1.1 Mapping SPH particles to a grid

The SPH particles from the merger simulation must be mapped onto a Cartesian
grid for the radiative transfer simulation. We first propagate the particles
based on the velocities they had at the end of the SPH simulation for a
further 0.5 seconds. After this time we assume homologous expansion and map
the particle positions onto a 1283 Cartesian grid spanned in velocity space.
We set the maximum absolute velocity of the velocity grid (i.e. the grid
boundary) equal to 0.5c (along each axis, noting that the corners of the
Cartesian grid extend to velocities up to 0.87c). The bulk of the ejecta mass
is well below this velocity. Few SPH particles are at velocities higher than
0.5c, corresponding to 5% of the total ejecta mass, and we do not expect these
high velocity particles at very low densities to have a strong influence on
the predicted light curves at the times we consider in our simulation. We map
the density and $Y_{e}$ of the SPH particles to the Cartesian grid. The
mapping of a quantity $A$ is performed via

$A(r)=\frac{1}{N}\sum_{a}A_{a}\frac{m_{a}}{\rho^{*}_{a}}W(|r-r_{a}|;h)$ (1)

with the sum running over all particles $a$, whose kernel function $W$ is
finite at the position $r$. Here $A_{a}$ is the value of A associated with
particle $a$, $m_{a}$ is the particle mass, $\rho^{*}_{a}$ is the conserved
rest-mass density, $r_{a}$ is the position of particle $a$, $W$ is the SPH
kernel function with smoothing length $h$. $N$ is introduced for normalization
and is defined by $N=\sum_{a}\frac{m_{a}}{\rho^{*}_{a}}W(|r-r_{a}|;h)$. For
the density we do not use the normalization but compute
$\rho(r)=\sum_{a}m_{a}W(|r-r_{a}|;h)$ (this expression includes an additional
approximation, namely equating $\rho$ and $\rho^{*}$). For the mapping, we
adopt the smoothing length and smoothing kernel from the SPH simulation. In
order to reduce the clumpiness of the outermost ejecta (likely related to the
limited particle resolution used in the underlying SPH simulation) we increase
the smoothing length by 50% for particles with absolute velocities $v>0.5$c.

####  2.1.2 Dynamical ejecta

Figure 1: 3D rendering of the dynamical ejecta mapped to a Cartesian grid. The
vz axis corresponds to the rotation axis (or polar axis) of the system and
v${}_{\rm z}=0$ represents the equatorial plane. For visual purposes, in this
figure we do not show grid cells containing density contributions from less
than 5 SPH particles. The colour scale is set by isosurfaces showing the
average $Y_{e}$ in a grid cell. Material near the equator generally has a
lower $Y_{e}$ than in the polar directions.  Figure 2: Distribution of ejecta
mass (of all SPH trajectories) into electron fraction bins
($\Delta$$Y_{e}=0.04$) (blue), and separated into polar (red) and equatorial
regions (green). The electron fraction near the poles is generally higher than
near the equator.  Figure 3: Electron fraction ($Y_{e}$) in each model grid
cell on the Cartesian grid, plotted against absolute velocity. We mark the
average $Y_{e}$ within 10 uniform angle bins to indicate angle dependence. We
note that due to the low resolution of the outer ejecta, at high velocities we
find that one SPH particle can influence the $Y_{e}$ of many grid cells (shown
in this plot by a line at the same $Y_{e}$ for a range of velocities).  Figure
4: Cell mass densities (on the Cartesian grid) versus absolute radial velocity
at 0.52 seconds after the merger, corresponding to the time at which SPH
particles are mapped to the Cartesian grid and homologous expansion is
assumed. The average densities in ten polar angle bins have been plotted in
colour to indicate the density range along different lines of sight. The outer
ejecta are less well resolved numerically, leading to scatter in the densities
at high velocities.

The ejecta is resolved by $\sim$2000 SPH particles, which become unbound
during the simulation (see Oechslin et al. 2002; Bauswein et al. 2013 for
ejecta criterion). The total mass of dynamical ejecta mapped onto the
radiative transfer grid was 0.0051 ${\mathrm{M}}_{\odot}$ (although we note
that mass ejection through different mechanisms still continues after the end
of the merger simulation). In Figure 1 we show a 3D rendering of the dynamical
ejecta once it has been mapped to the grid, indicating the 3-dimensional
structure produced by the merger simulation. The colour scale indicates the
$Y_{e}$ of the ejecta (where $Y_{e}$ is defined at the end time of the merger
simulation). Generally, lower $Y_{e}$ material is found near the equator,
while higher $Y_{e}$ material is found in the polar directions. This is
demonstrated by Figure 2, showing the $Y_{e}$ of the trajectories at the end
of the SPH simulation. The trend of drecreasing $Ye$ from the pole to the
equator is consistent with previous works, e.g. Radice et al. (2018); Foucart
et al. (2020); Kullmann et al. (2022). However, we also find that lower
$Y_{e}$ material is mixed with the higher $Y_{e}$ material in the poles (i.e.
the polar directions are not composed solely of high $Y_{e}$ material). A
similar effect was reported by Just et al. (2022). As a result we expect lower
opacities in the polar directions, due to the lower lanthanide fraction
expected to be synthesised for the lower $Y_{e}$ material. The $Y_{e}$ mapped
to the Cartesian grid is shown in Figure 3. To analyse the 3D density
structure of the ejecta, we plot the density of each model grid cell in Figure
4. To indicate the angular dependence, we divide the ejecta into 10 uniform
solid-angle bins in the polar direction relative to the positive z-axis and
plot the average density within each angle bin. Due to the limited resolution
of the underlying SPH simulation, the outermost velocities may suffer from
purely numerical particle noise. The impact of this particle noise will need
to be assessed in future studies using a higher resolution. While we would
expect there to be some difference in the densities in each polar direction
due to statistical fluctuations, we can not rule out that the relatively large
difference in density between the poles in this model could be numerical.

###  2.2 Nucleosynthesis Calculations

To obtain the energy released by the merger simulation, time dependent
nucleosynthesis calculations are carried out for each SPH particle trajectory.
The simulation provided thermodynamical histories for all the SPH particles
only up to $t_{0}\sim 20$ ms after merger. For later times, we extrapolate the
density evolution assuming homologous expansion, i.e.
$\rho(t)r(t)^{3}=\rho_{0}r_{0}^{3}$ and $r(t)=r_{0}+v_{0}(t-t_{0})$ with the
subscript $0$ denoting the end of the simulation data. This corresponds to a
density evolution as:

$\rho(t)=\rho_{0}\left(\frac{\Delta+t_{0}}{\Delta+t}\right)^{3},\quad\Delta=\frac{r_{0}}{v_{0}}-t_{0}.$
(2)

We start the nucleosynthesis calculations at a temperature of $T=10$ GK or the
lowest temperature reached in the SPH simulation if this is higher. Under
these conditions our initial composition is well described by nuclear
statistical equilibrium. During the nucleosynthesis calculations, the density
is evolved based on the SPH simulation data together with equation (2). We use
the same nuclear reaction network as in Mendoza-Temis et al. (2015) together
with the set of nuclear reactions labelled “FRDM”. Briefly, it consists of
neutron-capture and photodissociation rates computed within the statistical
model using the FRDM masses (Möller et al., 1995). For nuclei with
experimentally unknown $\beta$-decay rates, we use the compilation of Möller
et al. (2003). Fission rates (Panov et al., 2010) have been computed based on
the Thomas-Fermi fission barriers of Myers & Świaţecki (1999). Finally,
$\alpha$-decay rates are computed using a Viola-Seaborg formula (Dong & Ren,
2005) for those nuclei without experimental values. We use the $Y_{e}$ value
at the end of the SPH simulation to determine the composition at the beginning
of the network calculations.

In Figure 5 we show the specific heating rate, $\dot{\mathrm{Q}}$, averaged
over all trajectories. $\dot{\mathrm{Q}}$ represents the total rate of energy
released into the ejecta from nuclear reactions (including the contribution
from neutrinos). Also plotted are the average energies released from
$\beta$-decays and $\alpha$-decays. Fission also contributes to the total
energy released, although the contribution is orders of magnitudes less than
$\beta$-decays and $\alpha$-decays by $10^{-3}$ seconds, and is most
significant at very early times. Energy from $\beta$-decays accounts for most
of the total heating rate in the time range considered here.

Figure 5: Heating rates from nucleosynthesis calculations. The total energy
released (on average for all of the trajectories) is shown by the dashed line.
We also show the total energy from $\beta$-decays and from $\alpha$-particles.

###  2.3 Radiative transfer

We use the time-dependent, multi-dimensional Monte Carlo radiative transfer
code, artis (Sim, 2007; Kromer & Sim, 2009, based on the methods of Lucy 2002,
2003, 2005) to predict the kilonova from the neutron star merger simulation
described in Section 2.1. We propagate $1.15\times 10^{8}$ Monte Carlo packets
between 0.02 and 120 days after the merger. We produce viewing-angle dependent
light curves by assigning escaping UVOIR-packets into time and directional
bins. The bin sizes have been chosen to keep Monte Carlo noise at a low level
for the observables we present here. We define 100 uniform solid-angle bins in
the polar and azimuthal directions.

####  2.3.1 Opacity treatment

Table 1: Mass absorption cross sections ($\kappa$) adopted in this study for each electron fraction $Y_{e}$ range. The values are motivated by the Planck-mean opacities calculated by Tanaka et al. (2020, table 1). These calculations assumed a temperature range of 5000 - 10000 K and a density of $\rho=1\times 10^{-13}$ g cm-3. *The lowest $Y_{e}$ opacity is underestimated due to lack of complete atomic data for actinides (Tanaka et al., 2020).  $Y_{e}$ | $\kappa$
---|---
| cm2 g-1
$Y_{e}$ $\leq 0.1$ | 19.5*
$0.1<$ $Y_{e}$ $\leq 0.15$ | 32.2
$0.15<$ $Y_{e}$ $\leq 0.2$ | 22.3
$0.2<$ $Y_{e}$ $\leq 0.25$ | 5.60
$0.25<$ $Y_{e}$ $\leq 0.3$ | 5.36
$0.3<$ $Y_{e}$ $\leq 0.35$ | 3.30
$Y_{e}$ > 0.35 | 0.96

We adopt a temperature independent grey-opacity treatment for the propagation
of all radiation (described by Sim 2007). We use $Y_{e}$ dependent grey
absorption cross-sections, based on the Planck mean opacities listed in table
1 of Tanaka et al. (2020). These mean opacities were calculated for
temperatures of 5000 - 10000 K (with $\rho=1\times 10^{-13}$ g cm-3 at a time
of 1 day). Specifically, we use the opacities listed in Table 1 for model
cells with an average $Y_{e}$ within the ranges in Table 1. We note that the
$Y_{e}$ in each grid cell is based on the $Y_{e}$ of the SPH particles at the
end of the merger simulation (not the $Y_{e}$ evolved in the nuclear network
calculation).

####  2.3.2 Energy from nucleosynthesis calculations

Figure 6: Total energy deposited between 0.02 and 120 days versus radial
velocity. Top: Deposited energy per unit mass. Bottom: Deposited energy per
cubic grid cell (constant velocity extent in x-y-z). Mean values within polar
angle bins are indicated in colour. We note that due to the low resolution of
the outer ejecta, at high velocities we find that one SPH particle can
influence many grid cells (shown in the upper plot by a line at the same
energy for a range of velocities).

In artis, the total energy in the simulation is defined at the start of the
simulation. We then create Monte Carlo packets of equal energy which are
propagated through the simulation grid. We define what fraction of the total
energy will be deposited in each model grid cell, and place packets according
to this distribution.

We define how energy is distributed in the model grid cells according to the
total energy released by the SPH particle trajectories that contribute to that
model grid cell. The total amount of energy released along a trajectory over
time is mapped to the 3D Cartesian grid using the same method as described in
Section 2.1.1. These energies are shown in Figure 6. This energy (in erg/g)
and the total mass in the model grid cell are used to set the total energy
generated in each model grid cell, and in the framework of the code, the
fraction of energy which will be deposited in the cell (also shown in Figure
6).

At the start of the simulation we also define the rate at which energy will be
deposited (the time at which a Monte Carlo packet will be placed in a model
grid cell). We assume that the energy in all model grid cells is deposited at
the same rate, which we take as the average $\dot{\mathrm{Q}}$ of all unbound
SPH particle trajectories (plotted in Figure 5). We consider energy released
from 8 seconds after the merger ($\ll$ start time of radiative transfer
simulation) until 120 days after the merger. Energy released before the start
of the simulation is accounted for by placing Monte Carlo packets in the
ejecta which are advected with the homologous flow. The energy of these
packets is reduced to account for adiabatic losses before the beginning of the
simulation.

Since $\beta$-decays are the dominant source of energy at the times we are
considering, we assume that all of the energy in our simulations comes from
$\beta$-decays. Following the results of Barnes et al. (2016), we assume that
35% of the $\beta$-decay energy emerges as neutrinos, 20% as $\beta$-particles
and 45% as $\gamma$-rays. The energy from the neutrinos will never thermalise,
so this energy is ignored in the simulations.

Barnes et al. (2016) find that low energy $\beta$-particles effectively
thermalise on kilonova timescales, and that even slightly tangled magnetic
fields are effective at trapping high energy $\beta$-particles. We assume that
all of the energy from the $\beta$-particles will thermalise locally and
instantaneously. In reality, the $\beta$-particles will continuously deposit
energy over some distance through Coulomb interactions with thermal electrons,
by ionising or exciting bound atomic electrons or by Bremsstrahlung emission.
The thermalisation of $\beta$-particles is dependent on the magnetic field
structure in the ejecta, which we neglect.

We include $\gamma$-ray transport in our simulations, as described by Sim
(2007). We approximate the emission energies of the $\gamma$-rays guided by
the early time $\gamma$-ray emission spectra obtained by Barnes et al. (2016)
(see their figure 3), which peaks at several hundred keV. The exact energies
sampled and the probability of injecting a gamma-packet at each energy are
shown in Figure 7.

Figure 7: Approximated probability distribution of $\gamma$-ray energies for
radioactive decays of r-process material, guided by the early time
$\gamma$-ray emission spectrum by Barnes et al. (2016). This includes 115
energy bins, which are sampled to obtain $\gamma$ energies used in this work.

##  3 Results

###  3.1 Dynamical ejecta

####  3.1.1 Bolometric light curves

Figure 8: Line of sight bolometric UVOIR light curves (upper panel) and
$\gamma$-ray light curves (lower panel). Each coloured line is the azimuthally
averaged light curve (averaged over each angle bin in cos($\theta$)), while
all 100 uniformly spaced viewing angle bins are plotted in light grey. When
presenting the angle-dependent light curves, we plot them as equivalent
isotropic luminosities (i.e. from the simulation we record the energy emitted
per second into each sold angle bin to obtain light curves in erg/s/sr for
each orientation; we then scale these to an equivalent isotropic luminosity by
multiplication by $4\pi\cdot$sr). Also marked is the total energy released (on
average) by nuclear reactions, $\dot{\mathrm{Q}}_{\gamma+\beta}$, excluding
the energy assumed to be lost to neutrino emission (35%). In the upper panel
we additionally plot the energy due to $\beta$-particles
($\dot{\mathrm{Q}}_{\beta}$) and in the lower panel the heating due to
$\gamma$-rays ($\dot{\mathrm{Q}}_{\gamma}$). The bolometric light curve of
AT2017gfo from Waxman et al. (2018) is plotted with the bolometric light
curves for reference.  Figure 9: Bolometric (UVOIR) equivalent isotropic
luminosities in 100 uniform viewing angle bins at 1 day after the merger.

The bolometric light curves for the dynamical ejecta are shown in Figure 8.
From very early times photons are able to escape the ejecta, hence we find
that the light curves begin declining within a fraction of a day, and do not
show a significant rising phase in our simulations. According to the
nucleosynthesis calculations, significant numbers of $\beta$-decays will occur
in the outer ejecta layers, hence, in our simulations, energy is deposited in
low density, low opacity outer regions, where it thermalises, and is emitted
as optical radiation from very early times. Since we assume that all
$\beta$-particle energy thermalises, this may overestimate the true amount of
energy that would thermalise, however, at such times the thermalisation
efficiency is likely to be high (Barnes et al., 2016). The lack of rise to
peak in bolometric light curves has previously been found, e.g. by Banerjee et
al. (2020), Klion et al. (2022) and Kawaguchi et al. (2022). The light curves
do, however, show a ‘shoulder’ at $\sim 1$ day, particularly in the polar
directions. At this time, the ejecta are becoming optically thin, and energy
stored in the ejecta is able to escape. Photons are preferentially emitted in
the polar directions due to the lower optical depths in these lines of sight
(see Figure 8). At times around 1 day after the merger, we find the strongest
angle variation in the light curves (see Figure 9), however this decreases
over time as the ejecta become optically thin. As seen in Figure 9, the light
curve also exhibits $\sim 30\%$ variations in azimuthal viewing angle, likely
reflecting that the ejecta do not show perfect cylindrical symmetry.

The heating rate marked in Figure 8 shows the angle-averaged, total amount of
energy available for heating the ejecta over time. Since we assume all energy
comes from $\beta$-decays in our simulations, the heating rate represents the
total energy from $\gamma$-rays plus $\beta$-particles from $\beta$-decays.
The energy lost to neutrino emission is excluded from the heating rate shown
in Figure 8. Since $\gamma$-rays thermalise inefficiently (see Section 3.1.2)
it is predominantly energy from $\beta$-particles powering the light curve. By
$\sim 5$ days the ejecta has become optically thin, and the light curve is
equal to the $\beta$-particle heating rate. At this time $\beta$-particles are
no longer expected to thermalise efficiently (e.g. Barnes et al., 2016,
although we note this was for a simplified ejecta structure) and so the
$\beta$-particle heating rate is likely overestimated here. However, energy
from $\alpha$-particles becomes more significant at these times. Since the
late time light curve is dependent on the assumed heating rate, this
highlights the importance of calculating the fractions of $\beta$-decay energy
going into neutrinos, $\gamma$-rays and $\beta$-particles, by following the
radioactive decays of r-process elements in the radiative transfer
simulations, as well as calculating how much of that energy will thermalise.
Accounting for $\alpha$-particles, and the corresponding thermalisation
efficiency of the decay products is likely also important at late times.

For reference, the bolometric light curve of AT2017gfo, constructed by Waxman
et al. (2018), is also plotted in Figure 8. The mass of the dynamical ejecta
in our model is relatively small ($\sim$10 times less than e.g. the mass
inferred for AT2017gfo by Smartt et al. 2017, using an Arnett-type model, of
$0.04\pm 0.01$ ${\mathrm{M}}_{\odot}$) and therefore the total amount of
energy produced in our model is lower, leading to fainter, faster evolving
light curves.

####  3.1.2 $\gamma$-ray light curves

We also show the emerging $\gamma$-ray light curves in each viewing angle bin
in Figure 8. For all except the earliest times, the $\gamma$-rays do not
thermalise and are able to free-stream out of the ejecta. By $\sim 0.1$ days
the emerging $\gamma$-ray energy is equal to the total $\gamma$-ray energy
rate, marked in Figure 8. The $\gamma$-ray light curves do not show a viewing
angle dependence. To observe the peak of the $\gamma$-ray light curve,
observations would need to be within the first hour after the merger.

Figure 10: Mean ejecta velocity at which Monte Carlo packets were last emitted
within a given angle bin, indicating the ejecta velocity where packets last
interacted.

(a) 0.2 days

(b) 0.5 days

(c) 1 day

(d) 2 days

Figure 11: The spatial distributions (slices at v${}_{\rm y}=0$) of last
interactions of escaping radiation at three epochs. Note packets can escape in
any direction – we make no selection based on escaping angle. Packets that
escaped representing $\gamma$-rays are excluded from this.

####  3.1.3 Ejecta emission velocities

As discussed in Section 3.1.1, initially most of the photon emission occurs in
the outer ejecta layers, hence we find that the bolometric light curves begin
declining from very early times and do not rise to a peak within the time
frame of our simulation. Figure 10 shows the mean ejecta velocity at which
escaping Monte Carlo packets underwent their last interaction in the
simulation – this gives an indication of the regions of the ejecta that are
contributing to the kilonova emission. The outer ejecta layers quickly become
optically thin, due to the high expansion velocities. As a result, the ejecta
velocities from which packets are emitted decrease rapidly within the first
day. We show the mean ejecta velocities from which packets are emitted in
Figure 10.

We also show slices of the ejecta indicating where the Monte Carlo packets
escape from at a given time in Figure 11. Initially, all packets are emitted
from the outermost (i.e. fastest) ejecta. Until around 1 day, this resembles a
photosphere, beneath which radiation does not escape the ejecta.

####  3.1.4 Temperature

Figure 12:  Average radiation temperature at the last interaction locations of
escaping packets, indicating the ejecta temperatures where radiation is
escaping from. We also mark the inferred temperatures from the spectra of
AT2017gfo by Smartt et al. (2017).  Figure 13: Light curves estimated from
radiation temperature. Also plotted are the light curves for AT2017gfo from
Villar et al. (2017), assuming a distance of 40 Mpc. We note that the light
curves show some Monte Carlo noise.

The extremely high expansion velocities cause the temperature of the ejecta to
drop rapidly. We estimate the temperature in each model grid cell by matching
the energy density of the Monte Carlo radiation field in the cell during each
timestep to that of a black body radiation field: we denote this equivalent
black body temperature as TJ. We show the average temperatures of the model
grid cells from which packets are escaping over time in Figure 12.
Specifically, we record the cell temperature where a packet last interacted
before escaping and bin the packets in time. The mean temperature in each time
bin is plotted.

The opacities we adopt from Tanaka et al. 2020 (see Table 1) were calculated
for temperatures of 5000 - 10000 K. In our simulations the temperatures are
approximately in this range, at least in the regions from which packets
escape, between 0.1 – 1 days.

At late times it is likely that we underestimate the temperature, due to the
grey and local thermodynamic equilibrium (LTE) approximations assumed here. In
the nebular phase, non-thermal processes would provide heating to the ejecta
(Hotokezaka et al., 2021; Pognan et al., 2022) and likely prevent cooling to
such low temperatures as suggested by our TJ estimate.

We also plot the mean ejecta temperatures at the location of the last
interaction of escaping radiation (specifically, the temperature of the model
grid cell where a packet last interacted before escaping), binned by time in
Figure 12, and we show how the temperatures compare to the inferred
temperatures of AT2017gfo by Smartt et al. (2017). Although the model
temperatures are cooler than AT2017gfo by the time of the observations, the
evolution of the temperature is similar. The cooler temperatures are likely
due to the lower mass of the ejecta model, and therefore the lower amount of
energy generated in the ejecta, as well the LTE assumptions made in our
simulation.

####  3.1.5 Colour curves from radiation temperature

Figure 14: Angle-averaged colour evolution, estimated from radiation
temperature. Also plotted are the colours of AT2017gfo for points where
observations were taken at the same time (i.e. we do not interpolate between
light curve points). Note that Monte Carlo noise can be seen here.

Our simulations use only grey opacities, so we do not calculate any frequency
dependence to produce spectra or colour information from the Monte Carlo
packets directly. However, we record the temperature of the grid cell where
each packet was last emitted from before escaping the ejecta, i.e. the point
at which the packet last underwent an interaction (see Figure 12). From this,
we assign an approximate frequency to the Monte Carlo packet by making the
simplistic assumption that it is governed by the equilibrium radiation
distribution (a black body function) at the local radiation temperature, TR
(TR = TJ under the LTE assumption). We do this by randomly sampling a
frequency from a black body at the temperature in the model grid cell from
which the packet escaped. From this we obtain an approximate spectral energy
distribution, and can generate band-limited light curves, which we show in
Figure 13.

The elemental compositions of the kilonova ejecta are dominated by lanthanides
and actinides, which are very effective at line-blocking blue wavelengths, and
therefore we would expect that the true spectra would be redder than a pure
black body, particularly at later times (e.g. Gillanders et al., 2022).
However, we would also expect the temperature evolution to change in non-grey
and non-LTE simulations, which would also affect the colour evolution. To
accurately determine band-limited light curves, frequency dependent
simulations are required, which goes beyond the scope of this work.

The light curves show an angle variation in magnitude of up to $\sim$1 mag.
The emission from polar directions is brighter than from equatorial
directions. As discussed in Section 2.1, we have a higher $Y_{e}$, and
therefore lower lanthanide fraction in the polar directions, which represents
a ‘blue’ component, and a lower $Y_{e}$ in the direction of the equator, which
would lead to a higher lanthanide fraction, representing a ‘red’ component.
However, this does not lead to significantly redder colours in the equatorial
direction than the polar directions. We show the angle dependent colour
evolution in Figure 14, which does not show a significant angle variation.

Due to the evolution of the temperature alone, we find that the colours show a
rapid evolution from blue to red over time in all lines of sight, which can be
seen in Figure 14.

####  3.1.6 Comparison to AT2017gfo

We plot the light curves of AT2017gfo (from Villar et al., 2017) in Figure 13
for reference, although again we note that the model we consider here is much
less massive than what was inferred for AT2017gfo, and therefore we do not
expect to match the brightness. The light curves have been corrected for
foreground reddening, assuming an extinction of E(B-V) = 0.11 mag (Smartt et
al., 2017).

Interestingly, the model light curves show a similar evolution to AT2017gfo,
although the model light curves are fainter and evolve faster. This could
suggest that in future events where less mass is ejected than in the case of
AT2017gfo, the light curves would show a faster decline and may have already
faded in the bluer bands by the time the first detections of AT2017gfo were
made. The red bands, however, remain brighter for longer, suggesting that
searches for the electromagnetic counterparts of future kilonova events should
focus on these bands. We note that the model H and K band light curves are
likely too red, since these are similar in brightness to AT2017gfo. As
discussed in Section 3.1.4, due to our grey approximation we likely
underestimate the temperature at later times, which could be responsible for
the very red late-time colours.

The rapid blue to red colour evolution found for our model is similar to that
observed for AT2017gfo (see Figure 14). Since we only include dynamical ejecta
in this simulation, this shows that a late-time, high opacity component is not
necessarily required to explain the colour evolution shown by AT2017gfo, as
has been found by e.g. Cowperthwaite et al. (2017), although we note that we
do not claim there was no secondary component. This suggests that the colour
evolution could be driven by the cooling of the ejecta. Future frequency
dependent simulations will be required to confirm this.

###  3.2 Late-time ejecta

So far, we have only considered the dynamical ejecta, expelled on timescales
of tens of milliseconds. Now we investigate in an approximate manner the
effects of ejecta expelled on longer timescales.

####  3.2.1 Model for late-time ejecta

Figure 15: Same as Figure 4 but including the spherically averaged late-time
ejecta. The densities of the late-time ejecta have been scaled (assuming
homologous expansion) to match the time in Figure 4 of 0.52 seconds.

We base the density structure for the secular ejecta on a simulation of a BH-
accretion torus, which is formed after the collapse of the hypermassive
neutron star. The BH-torus model is similar to the models of Just et al.
(2015), but includes a special relativistic treatment as well as additional,
manually constructed ejecta components. The components represent dynamical
ejecta and a neutrino-driven wind from a hyper-massive neutron star (HMNS) as
described by Ito et al. (2021, see therein for more details). The initial disk
mass was 0.1 M⊙, the mass of the HMNS wind was 0.006 M⊙, the black hole mass
was 2.7 M⊙, and the assumed time of BH formation after the merger was 0.05 s.
The torus and wind components from the long term evolution simulation are
spherically averaged and then mapped to the 3D Cartesian grid (i.e. the
dynamical ejecta component from this simulation is neglected). In each cell of
our Cartesian grid, we add the corresponding mass from the remnant simulation
to the mass of the dynamical ejecta, described in Section 2.1. We do not use
the $Y_{e}$ information from the remnant model, but keep the $Y_{e}$ and
energy (per unit of mass) distributions that were determined for the dynamical
ejecta. Any empty cells in the original dynamical ejecta model without a
$Y_{e}$ (close to the centre where the bound remnant was removed) are given a
$Y_{e}$ of 0.5 and a total heating energy equal to the average heating rate of
the dynamical ejecta trajectories (marked in Figure 5) integrated over time.
The additional mass of the late-time ejecta is 0.019 ${\mathrm{M}}_{\odot}$,
which gives a total ejecta mass of 0.024 ${\mathrm{M}}_{\odot}$.

In Figure 15 we show the model cell densities after including the mass from
the torus and wind ejecta. Due to the very high central densities introduced
by the late-time ejecta model, significantly higher amounts of energy will be
generated in the low velocity, central regions of the ejecta (since energy
released is proportional to mass). We therefore expect the light curve to be
brighter at late times when this energy diffuses out of the ejecta.

####  3.2.2 Effect on light curves of late-time ejecta

Figure 16: Bolometric angle-averaged light curves for the model including the
secular ejecta component and without. When presenting the angle-dependent
light curves, we plot them as equivalent isotropic luminosities (multiplied by
$4\pi\cdot$sr). We note that Monte Carlo noise can be seen in this plot.

We carry out a radiative transfer simulation on the model combining the late-
time ejecta with the dynamical ejecta. This simulation was carried out between
0.05 and 120 days with $2.88\times 10^{7}$ Monte Carlo packets.

To demonstrate the effect of the late-time ejecta on the bolometric light
curve in comparison to the dynamical ejecta alone, we plot the angle-averaged
bolometric light curves of each model in Figure 16. Due to the greater mass in
the center of this model, the increased energy generated in the high density
inner regions leads to a peak (or shoulder) in the bolometric light curve at
$\sim 5$ days after the merger, which was not produced by the dynamical ejecta
alone. The increased energy thermalising at low velocities and high optical
depths takes longer to diffuse out of the ejecta, and we find a later peak in
the bolometric light curve. This highlights the necessity of long term
hydrodynamic simulations, since the late-time light curve is dominated by the
lower velocity, late-time ejecta. At the same time, the small difference
between the two cases at early times ($\lesssim 1$ day) suggests that the
early light curve may be dominated by the dynamical ejecta.

At late times (>10 days) the model light curve shows a slower decline than the
observations of AT2017gfo. This is likely due to our assumption that all
$\beta$-particle energy will thermalise at all times, and suggests that by
these times the actual thermalisation rate is lower than we assume.

The late time light curve is now a similar brightness to AT2017gfo, suggesting
the mass of our model at low velocities may be similar to that of AT2017gfo.
However, the relative faintness at early times indicates that our model would
require more mass at higher velocities to match the brightness of AT2017gfo.

####  3.2.3 Effect on mean emission velocities

Figure 17: Mean angle-averaged ejecta radial velocity ($v=rt$) of last
interaction point of escaping radiation.

We show the mean angle-averaged ejecta velocities at which Monte Carlo packets
are emitted from for the simulations with and without the late time ejecta
component in Figure 17. At early times, the mean velocities are the same in
both models, but at later times the mean velocities are much lower in the
model with combined dynamical and late-time ejecta. At these times the
emission is primarily from the low velocity, high density late-time ejecta
component.

####  3.2.4 Effect on light curve angle-dependence

Figure 18: Line of sight dependent bolometric light curves for the model
including the late-time ejecta component. We note that Monte Carlo noise is
visible in this plot.

We also show the angle-dependent bolometric light curves for the model with
the late-time ejecta in Figure 18. We note that the higher densities, and
therefore larger amount of energy at low velocities, have slightly increased
the predicted variation with observation angle. We find that the brightest
lines of sight are still in the polar directions, which have the lowest
densities and lowest opacities, and therefore we find photons are
preferentially emitted in these directions. The faintest lines of sight are
those near the equator, as these have the highest densities and highest
opacities. The increased viewing angle dependence is due to more packets of
photons being generated at lower velocities in the central regions of the
ejecta. Instead of the packets escaping quickly from the lower density
dynamical ejecta (mostly in the line of sight where the packet was created),
the packets take longer to diffuse out of the higher density material, but
preferentially in the polar directions since these have lower opacities, which
leads to the increased angle dependence. This highlights the need for long-
term evolution simulations, and the need to include ejecta expelled on longer
timescales. Future simulations where the late-time ejecta are self-
consistently calculated are required.

##  4 Discussion and conclusions

In this paper we aimed to predict the kilonova emission from the ejecta
density and composition directly from a neutrino-hydrodynamics simulation of a
binary neutron star merger. Initially, we considered the dynamical ejecta
alone, which are expelled up to $\sim 20$ milliseconds after the merger. We
presented bolometric light curves, and found that we do not expect a rise to
peak within the time frame of our simulation (starting at 0.02 days). This is
due to the distribution of energy generated in the ejecta. Energy is generated
and thermalised at very high velocities and due to low optical depths is able
to immediately diffuse out of the ejecta, leading to a bolometric light curve
that begins declining from very early times after the merger. A similar result
has previously been reported, e.g. by Banerjee et al. (2020), Klion et al.
(2022) and Kawaguchi et al. (2022).

We do however, find a ‘shoulder’ in the light curve, particularly in the polar
directions, around the time when the ejecta start to become optically thin.
Just et al. (2022) predicted light curves for a similar model to the dynamical
ejecta model we consider here (they also consider a 1.35-1.35 M⊙ model using
the SFHo equation of state and the neutrino leakage scheme ILEAS). They also
find a ‘shoulder’ feature at a similar time, although their shoulder
luminosity is lower than ours by a factor of a few.

We find that light curves are brighter in the polar directions around the time
of the shoulder compared to the equatorial lines of sight. After the shoulder,
once the ejecta become optically thin we no longer see any angle dependence.
The angle variation is similar to the level found by Just et al. (2022) and
Kawaguchi et al. (2021).

We assumed that all energy came from $\beta$-decays. We found that
$\gamma$-rays are only thermalised within approximately the first hour after
the merger. The late time light curve (after the shoulder) closely follows the
injected heating rate, since the ejecta have become optically thin. In future
simulations it will be important to calculate the fraction of thermalised
heating energy self-consistently.

We showed that in the day after the merger, most Monte Carlo packets of
radiation were emitted from the outer layers of the ejecta, which at early
times resembles a photosphere. This suggests that photospheric approximations
may be reasonable at early times. We note that the angular dependence of this
photospheric-like emission (e.g. 0.5 days) does not necessarily lead to strong
asymmetries in the light curves (i.e. the structure of the photospheric-like
emission does not correspond to structure in the angle variation of the light
curves). Due to the high expansion velocities, the density of the ejecta
quickly decreases, as does the optical depth, leading to lower emission
velocities. By $\sim$ 2 days we find packets are emitted from all velocities
in the ejecta, and the emission no longer resembles a photosphere.

The ejecta temperatures also rapidly decrease, again due to the high expansion
velocities. Although we do not carry out frequency dependent simulations, we
estimate the frequency of escaping packets from the radiation temperature in
the cell from which the packet escaped. We did this by sampling a black body
distribution at the radiation temperature. We generated band light curves
based on this, and showed that due to the temperature evolution we find a
rapid evolution from blue to red colours, which is similar to that observed in
AT2017gfo, however, the evolution is on much shorter timescales. This suggests
that cooling of the ejecta could drive the colour evolution. However,
frequency dependent simulations of models including self-consistently evolved
late time ejecta will be required to investigate this further. Since this
model was not tuned in any way to match the observations of AT2017gfo, it is
promising that we find a similar colour evolution from our model with
relatively simple assumptions.

We also carried out a simulation where we included, in an approximate,
spherically symmetric manner, ejecta expelled from the merger remnant. We
investigated the effect on the predicted light curves, and the extent to which
dynamical ejecta may be responsible for the observed light curve of AT2017gfo.
We found that the light curve at early times does not show a strong
sensitivity with respect to the late-time ejecta component, although it
becomes increasingly important with time. After $\sim$1 day the light curve is
dominated by the late-time ejecta component, due to the higher densities, and
therefore higher energy in the center of the ejecta, which is able to diffuse
out at later times. The late time light curve is of similar brightness to that
of AT2017gfo, however, we still do no match the initial brightness. This
suggests that the higher masses of high velocity ejecta to match the initial
brightness of AT2017gfo.

We also showed that the viewing angle dependence does not disappear, and even
increases, with the addition of a massive, spherically symmetric late-time
ejecta. This suggests that the dynamical ejecta can have a significant impact
on the angle dependence even at late times. This highlights the importance of
long-term evolution simulations for predicting kilonovae.

## Acknowledgements

We thank Ricard Ardevol-Pulpillo and Thomas Janka for providing the ILEAS
scheme. CEC, AB and OJ acknowledge support by the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation program
under grant agreement No. 759253. AB, G.M.-P and OJ acknowledge support by
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Project-ID
279384907 - SFB 1245. AB and VV acknowledge support by DFG - Project-ID
138713538 - SFB 881 (“The Milky Way System”, subproject A10). AB and G.M.-P
acknowledge support by the State of Hesse within the Cluster Project ELEMENTS.
The work of SAS was supported by the Science and Technology Facilities Council
[grant numbers ST/P000312/1, ST/T000198/1]. G.M.-P and LJS acknowledge support
by the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (ERC Advanced Grant KILONOVA No. 885281). This
work was performed using the Cambridge Service for Data Driven Discovery
(CSD3), part of which is operated by the University of Cambridge Research
Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The
DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital
grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1.
DiRAC is part of the National e-Infrastructure. OJ is grateful for
computational support by the HOKUSAI computing facility at RIKEN. CEC, OJ and
VV are grateful for computational support by the VIRGO cluster at GSI. NumPy
and SciPy (Oliphant, 2007), IPython (Pérez & Granger, 2007), Matplotlib
(Hunter, 2007), PyVista (Sullivan & Kaszynski, 2019) and
artistools111https://github.com/artis-mcrt/artistools/ were used for data
processing and plotting.

## Data Availability

Data will be made available upon reasonable request.

## References

  * Abbott et al. (2017) Abbott B. P., et al., 2017, Phys. Rev. Lett., 119, 161101
  * Ardevol-Pulpillo et al. (2019) Ardevol-Pulpillo R., Janka H. T., Just O., Bauswein A., 2019, MNRAS, 485, 4754
  * Banerjee et al. (2020) Banerjee S., Tanaka M., Kawaguchi K., Kato D., Gaigalas G., 2020, ApJ, 901, 29
  * Barnes & Kasen (2013) Barnes J., Kasen D., 2013, ApJ, 775, 18
  * Barnes et al. (2016) Barnes J., Kasen D., Wu M.-R., Martínez-Pinedo G., 2016, ApJ, 829, 110
  * Bauswein et al. (2013) Bauswein A., Goriely S., Janka H.-T., 2013, ApJ, 773, 78
  * Bulla et al. (2021) Bulla M., et al., 2021, MNRAS, 501, 1891
  * Cowperthwaite et al. (2017) Cowperthwaite P. S., et al., 2017, ApJ, 848, L17
  * Curtis et al. (2021) Curtis S., Mösta P., Wu Z., Radice D., Roberts L., Ricigliano G., Perego A., 2021, arXiv e-prints, p. arXiv:2112.00772
  * Domoto et al. (2021) Domoto N., Tanaka M., Wanajo S., Kawaguchi K., 2021, ApJ, 913, 26
  * Dong & Ren (2005) Dong T., Ren Z., 2005, doi:10.1140/epja/i2005-10142-y, 26, 69
  * Drout et al. (2017) Drout M. R., et al., 2017, Science, 358, 1570
  * Eichler et al. (1989) Eichler D., Livio M., Piran T., Schramm D. N., 1989, Nature, 340, 126
  * Even et al. (2020) Even W., et al., 2020, ApJ, 899, 24
  * Fernández & Metzger (2013) Fernández R., Metzger B. D., 2013, MNRAS, 435, 502
  * Foucart et al. (2020) Foucart F., Duez M. D., Hebert F., Kidder L. E., Pfeiffer H. P., Scheel M. A., 2020, ApJ, 902, L27
  * Fujibayashi et al. (2018) Fujibayashi S., Kiuchi K., Nishimura N., Sekiguchi Y., Shibata M., 2018, ApJ, 860, 64
  * Gillanders et al. (2022) Gillanders J. H., Smartt S. J., Sim S. A., Bauswein A., Goriely S., 2022, MNRAS, 515, 631
  * Heinzel et al. (2021) Heinzel J., et al., 2021, MNRAS, 502, 3057
  * Hotokezaka et al. (2013) Hotokezaka K., Kiuchi K., Kyutoku K., Okawa H., Sekiguchi Y.-i., Shibata M., Taniguchi K., 2013, Phys. Rev. D, 87, 024001
  * Hotokezaka et al. (2021) Hotokezaka K., Tanaka M., Kato D., Gaigalas G., 2021, MNRAS, 506, 5863
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Ito et al. (2021) Ito H., Just O., Takei Y., Nagataki S., 2021, ApJ, 918, 59
  * Just et al. (2015) Just O., Bauswein A., Ardevol Pulpillo R., Goriely S., Janka H. T., 2015, MNRAS, 448, 541
  * Just et al. (2022) Just O., Kullmann I., Goriely S., Bauswein A., Janka H. T., Collins C. E., 2022, MNRAS, 510, 2820
  * Kasen et al. (2015) Kasen D., Fernández R., Metzger B. D., 2015, MNRAS, 450, 1777
  * Kasen et al. (2017) Kasen D., Metzger B., Barnes J., Quataert E., Ramirez-Ruiz E., 2017, Nature, 551, 80
  * Kawaguchi et al. (2021) Kawaguchi K., Fujibayashi S., Shibata M., Tanaka M., Wanajo S., 2021, ApJ, 913, 100
  * Kawaguchi et al. (2022) Kawaguchi K., Fujibayashi S., Hotokezaka K., Shibata M., Wanajo S., 2022, ApJ, 933, 22
  * Klion et al. (2022) Klion H., Tchekhovskoy A., Kasen D., Kathirgamaraju A., Quataert E., Fernández R., 2022, MNRAS, 510, 2968
  * Korobkin et al. (2021) Korobkin O., et al., 2021, ApJ, 910, 116
  * Kromer & Sim (2009) Kromer M., Sim S. A., 2009, MNRAS, 398, 1809
  * Krüger & Foucart (2020) Krüger C. J., Foucart F., 2020, Phys. Rev. D, 101, 103002
  * Kullmann et al. (2022) Kullmann I., Goriely S., Just O., Ardevol-Pulpillo R., Bauswein A., Janka H. T., 2022, MNRAS, 510, 2804
  * Lattimer & Schramm (1976) Lattimer J. M., Schramm D. N., 1976, ApJ, 210, 549
  * Li & Paczyński (1998) Li L.-X., Paczyński B., 1998, ApJ, 507, L59
  * Lucy (2002) Lucy L. B., 2002, A&A, 384, 725
  * Lucy (2003) Lucy L. B., 2003, A&A, 403, 261
  * Lucy (2005) Lucy L. B., 2005, A&A, 429, 19
  * Mendoza-Temis et al. (2015) Mendoza-Temis J. J., Wu M.-R., Langanke K., Martínez-Pinedo G., Bauswein A., Janka H.-T., 2015, Phys. Rev.~C, 92, 055805
  * Metzger (2017) Metzger B. D., 2017, Living Reviews in Relativity, 20, 3
  * Metzger (2019) Metzger B. D., 2019, Living Reviews in Relativity, 23, 1
  * Metzger & Fernández (2014) Metzger B. D., Fernández R., 2014, MNRAS, 441, 3444
  * Metzger et al. (2010) Metzger B. D., et al., 2010, MNRAS, 406, 2650
  * Möller et al. (1995) Möller P., Nix J. R., Myers W. D., Swiatecki W. J., 1995, doi:10.1006/adnd.1995.1002, 59, 185
  * Möller et al. (2003) Möller P., Pfeiffer B., Kratz K.-L., 2003, Phys. Rev.~C, 67, 055802
  * Myers & Świaţecki (1999) Myers W. D., Świaţecki W. J., 1999, Phys. Rev.~C, 60, 014606
  * Nativi et al. (2021) Nativi L., Bulla M., Rosswog S., Lundman C., Kowal G., Gizzi D., Lamb G. P., Perego A., 2021, MNRAS, 500, 1772
  * Neuweiler et al. (2022) Neuweiler A., Dietrich T., Bulla M., Vivekanandji Chaurasia S., Rosswog S., Ujevic M., 2022, arXiv e-prints, p. arXiv:2208.13460
  * Oechslin et al. (2002) Oechslin R., Rosswog S., Thielemann F.-K., 2002, Phys. Rev. D, 65, 103005
  * Oliphant (2007) Oliphant T. E., 2007, Computing in Science & Engineering, 9, 10
  * Panov et al. (2010) Panov I. V., Korneev I. Y., Rauscher T., Martínez-Pinedo G., Kelić-Heil A., Zinner N. T., Thielemann F., 2010, doi:10.1051/0004-6361/200911967, 513, A61
  * Perego et al. (2014) Perego A., Rosswog S., Cabezón R. M., Korobkin O., Käppeli R., Arcones A., Liebendörfer M., 2014, MNRAS, 443, 3134
  * Pérez & Granger (2007) Pérez F., Granger B. E., 2007, Computing in Science & Engineering, 9, 21
  * Pognan et al. (2022) Pognan Q., Jerkstrand A., Grumer J., 2022, MNRAS, 510, 3806
  * Radice et al. (2018) Radice D., Perego A., Hotokezaka K., Fromm S. A., Bernuzzi S., Roberts L. F., 2018, The Astrophysical Journal, 869, 130
  * Radice et al. (2022) Radice D., Bernuzzi S., Perego A., Haas R., 2022, MNRAS, 512, 1499
  * Sekiguchi et al. (2015) Sekiguchi Y., Kiuchi K., Kyutoku K., Shibata M., 2015, Phys. Rev. D, 91, 064059
  * Siegel & Metzger (2018) Siegel D. M., Metzger B. D., 2018, ApJ, 858, 52
  * Sim (2007) Sim S. A., 2007, MNRAS, 375, 154
  * Smartt et al. (2017) Smartt S. J., et al., 2017, Nature, 551, 75
  * Steiner et al. (2013) Steiner A. W., Hempel M., Fischer T., 2013, ApJ, 774, 17
  * Sullivan & Kaszynski (2019) Sullivan C. B., Kaszynski A., 2019, Journal of Open Source Software, 4, 1450
  * Tanaka & Hotokezaka (2013) Tanaka M., Hotokezaka K., 2013, ApJ, 775, 113
  * Tanaka et al. (2020) Tanaka M., Kato D., Gaigalas G., Kawaguchi K., 2020, Monthly Notices of the Royal Astronomical Society, 496, 1369
  * Villar et al. (2017) Villar V. A., et al., 2017, ApJ, 851, L21
  * Watson et al. (2019) Watson D., et al., 2019, Nature, 574, 497
  * Waxman et al. (2018) Waxman E., Ofek E. O., Kushnir D., Gal-Yam A., 2018, MNRAS, 481, 3423
  * Wollaeger et al. (2021) Wollaeger R. T., et al., 2021, ApJ, 918, 10
  * Wu et al. (2022) Wu Z., Ricigliano G., Kashyap R., Perego A., Radice D., 2022, MNRAS,