# Self-consistent 3D radiative transfer for kilonovae: directional spectra
from merger simulations

Luke J. Shingles GSI Helmholtzzentrum für Schwerionenforschung, Planckstraße
1, 64291 Darmstadt, Germany  Christine E. Collins GSI Helmholtzzentrum für
Schwerionenforschung, Planckstraße 1, 64291 Darmstadt, Germany  Vimal Vijayan
GSI Helmholtzzentrum für Schwerionenforschung, Planckstraße 1, 64291
Darmstadt, Germany  Department of Physics and Astronomy, Ruprecht-Karls-
Universität Heidelberg, Im Neuenheimer feld 226, 69120 Heidelberg, Germany
Andreas Flörs GSI Helmholtzzentrum für Schwerionenforschung, Planckstraße 1,
64291 Darmstadt, Germany  Oliver Just GSI Helmholtzzentrum für
Schwerionenforschung, Planckstraße 1, 64291 Darmstadt, Germany  Astrophysical
Big Bang Laboratory, RIKEN Cluster for Pioneering Research, 2-1 Hirosawa,
Wako, Saitama 351-0198, Japan  Gerrit Leck GSI Helmholtzzentrum für
Schwerionenforschung, Planckstraße 1, 64291 Darmstadt, Germany  Institut für
Kernphysik (Theoriezentrum), Technische Universität Darmstadt,
Schlossgartenstraße 2, D-64289 Darmstadt, Germany  Zewei Xiong GSI
Helmholtzzentrum für Schwerionenforschung, Planckstraße 1, 64291 Darmstadt,
Germany  Andreas Bauswein GSI Helmholtzzentrum für Schwerionenforschung,
Planckstraße 1, 64291 Darmstadt, Germany  Gabriel Martínez-Pinedo GSI
Helmholtzzentrum für Schwerionenforschung, Planckstraße 1, 64291 Darmstadt,
Germany  Institut für Kernphysik (Theoriezentrum), Technische Universität
Darmstadt, Schlossgartenstraße 2, D-64289 Darmstadt, Germany  Stuart A. Sim
Astrophysics Research Centre, School of Mathematics and Physics, Queens
University Belfast, Belfast BT7 1NN, UK

###### Abstract

We present three-dimensional radiative transfer calculations for the ejecta
from a neutron star merger that include line-by-line opacities for tens of
millions of bound-bound transitions, composition from an r-process nuclear
network, and time-dependent thermalization of decay products from individual
$\alpha$ and $\beta^{-}$ decay reactions. In contrast to expansion opacities
and other wavelength-binned treatments, a line-by-line treatment enables us to
include fluorescence effects and associate spectral features with the emitting
and absorbing lines of individual elements. We find variations in the
synthetic observables with both the polar and azimuthal viewing angles. The
spectra exhibit blended features with strong interactions by Ce III, Sr II, Y
II, and Zr II that vary with time and viewing direction. We demonstrate the
importance of wavelength-calibration of atomic data using a model with
calibrated Sr, Y, and Zr data, and find major differences in the resulting
spectra, including a better agreement with AT2017gfo. The synthetic spectra
for near-polar inclination show a feature at around 8000 Å, similar to
AT2017gfo. However, they evolve on a more rapid timescale, likely due to the
low ejecta mass (0.005 ${\rm M}_{\sun}$) as we take into account only the
early ejecta. The comparatively featureless spectra for equatorial observers
gives a tentative prediction that future observations of edge-on kilonovae
will appear substantially different from AT2017gfo. We also show that 1D
models obtained by spherically averaging the 3D ejecta lead to dramatically
different direction-integrated luminosities and spectra compared to full 3D
calculations.

Neutron stars — nuclear astrophysics — r-process — transient sources —
gravitational wave astronomy — radiative transfer simulations

††software:  NumPy (Oliphant, 2007), Matplotlib (Hunter, 2007),
artis666https://github.com/artis-mcrt/artis (Sim, 2007; Kromer & Sim, 2009;
Shingles et al., 2020), and artistools777https://github.com/artis-
mcrt/artistools (Shingles et al., 2023)

##  1 Introduction

The gravitational wave event GW170817 from a neutron star merger (NSM) was
complemented by a rich set of electromagnetic signals (Abbott et al., 2017)
including the associated kilonova AT2017gfo (Smartt et al., 2017; Villar et
al., 2017). Simple analysis of the luminosity decline rate of AT2017gfo shows
consistency with radioactively decaying material that has undergone the rapid
neutron capture process (r-process), responsible for synthesizing many of the
elements heavier than Fe (Metzger et al., 2010). Another indication that very
heavy elements have been produced is the inferred high opacity of the ejecta
(Kasen et al., 2017), which is considered to be a signature for the presence
of lanthanides or actinides (Kasen et al., 2013).

The observed time series of spectra and the luminosity evolution of AT2017gfo
have provided a powerful set of constraints for testing theoretical NSM
models, the high-density Equation of State (EoS), and r-process
nucleosynthesis. However, linking the electromagnetic observations back to the
underlying physical conditions requires accurate models of radiative transfer
in the ejecta, which depend on accurate atomic data.

The availability of atomic data from experiments is too limited to constitute
a dataset for radiative transfer calculations. Instead, systematic atomic
structure calculations for trans-Fe elements (including lanthanides) have been
performed by Fontes et al. (2020) and Tanaka et al. (2020), and both provide
tables of wavelength-binned opacities (i.e., opacities obtained as averages
over discrete wavelength intervals). In addition, the energy levels and
transition data of Tanaka et al. (2020) are available via the Japan-Lithuania
Opacity Database for Kilonova111http://dpc.nifs.ac.jp/DB/Opacity-Database/
(Kato et al., 2021), which enables the calculation of individual line
opacities. However, the energy levels of Tanaka et al. (2020) have not been
calibrated to be consistent with observed transition lines (as noted by Domoto
et al., 2021), which can severely affect the locations of spectral features
and hinder the comparison between synthetic and observed spectra. For some
elements and ionization stages, calibrated atomic data is available (e.g.,
Kurucz, 2018), but the coverage of the heavy elements is incomplete.

To date, most radiative transfer simulations for kilonovae use simplified
ejecta models, such as with analytical density structures (Metzger et al.,
2010; Kasen et al., 2013; Bulla, 2019; Even et al., 2020; Banerjee et al.,
2020; Korobkin et al., 2021; Domoto et al., 2021; Wollaeger et al., 2021;
Gillanders et al., 2022; Pognan et al., 2022a, b). With recent claims of
observational constraints on the ejecta geometry (e.g., Sneppen et al., 2023),
it is especially important to understand whether multi-dimensional simulations
can be reconciled with observations. While there are studies that consider NSM
ejecta distributions directly from simulations in two dimensions (Kasen et
al., 2015; Kawaguchi et al., 2018; Korobkin et al., 2021; Kawaguchi et al.,
2021; Bulla, 2023; Just et al., 2022; Kawaguchi et al., 2022) or three
dimensions (Tanaka & Hotokezaka, 2013; Collins et al., 2023; Neuweiler et al.,
2023), these calculations use some form of simplified opacity treatment,
generally by combining line transitions into wavelength bins rather than
calculating individual line-by-line opacities.

While line-by-line opacities are widely used for supernova calculations, the
same method applied to kilonovae becomes computationally more expensive due to
the substantial abundances of trans-Fe heavy elements, which have complex
electronic structures that lead to vast numbers of energy levels and
transition lines (Kasen et al., 2013). Published kilonova calculations with
line-by-line opacities are limited to simplified ejecta models and often a
small selection of elements. For example, Pognan et al. (2022a, b) apply line-
by-line opacities in their non-LTE investigation of Te, Ce, Pt, and Th with a
spherically-symmetric ejecta model but do not predict synthetic spectra.
Gillanders et al. (2022) present synthetic spectra calculated with line-by-
line opacities for an extensive set of atomic data from H to U, but assume
spherically symmetric ejecta with an inner photospheric boundary.

Using wavelength-binned opacities instead of line-by-line opacities has
consequences beyond merely decreasing the wavelength resolution of the
resulting synthetic spectra. With binned opacities, the absorbing atomic
transitions are unknown in the radiative transfer calculation, and re-
radiation following absorption is typically handled either by assuming perfect
scattering (no wavelength change) or perfect thermalization (emission
according to the Planck function). In order to avoid this simplification and
apply a detailed treatment of the subsequent de-excitation radiation
(including fluorescence, Lucy, 2002), it is necessary to associate each bound-
bound absorption with a specific transition. When emission and absorption are
associated with specific transitions, features in the synthetic spectra can be
directly linked with particular elements, which would otherwise (i.e., with a
wavelength-binned scheme) have to be inferred based on estimated line
strengths or a comparison of several calculations with different compositions
(e.g., Domoto et al., 2021; Sneppen & Watson, 2023).

The only direct comparison between synthetic spectra calculated using line-by-
line opacities, line-binned opacities, and expansion opacities is made by
Fontes et al. (2020) for a one-dimensional test problem with a pure-Nd
composition, and neglecting fluorescence effects. Even this simple model leads
to differences in resulting spectra between line-by-line and binned opacities
(see Fontes et al., 2020, figure 13).

In this work, we present 3D radiative-transfer calculations that include, for
the first time, a line-by-line Sobolev treatment for tens of millions of lines
based on three-dimensional ejecta self-consistently obtained from NSM
simulations with a nuclear network. We show that our forward modeling approach
can produce a sequence of spectra in the polar direction that look similar to
AT2017gfo, and predict relatively featureless spectra for equatorial
observers. We also associate spectral features with particular elements in the
ejecta composition and demonstrate the importance of using multidimensional
models.

##  2 Method

###  2.1 NSM ejecta model and r-process calculation

Our simulation of two merging 1.35 M⊙ neutron stars adopts the SFHo equation
of state (Steiner et al., 2013) and has been calculated with a 3D general
relativistic smoothed-particle hydrodynamics (SPH) code (Oechslin et al.,
2002; Bauswein et al., 2010) employing the conformal flatness condition
(Isenberg & Nester, 1980; Wilson et al., 1996) and including a recently-
developed leakage-plus-absorption neutrino treatment (ILEAS, Ardevol-Pulpillo
et al., 2019).

Since the SPH simulation terminates at 20 ms after the merger, we only capture
the early, dynamical ejecta, which in our case has a mass of 0.005 ${\rm
M}_{\sun}$ and is made up of $\sim$2000 SPH particles. The same NSM simulation
was used in three-dimensional radiative transfer calculations with wavelength-
independent (gray) opacities by Collins et al. (2023).

For performing 3D frequency-dependent radiative transfer calculations with
artis 222https://github.com/artis-mcrt/artis/
The C++ source code is provided as is with no support or documentation under
the BSD 3-Clause License. (Sim, 2007; Kromer & Sim, 2009; Shingles et al.,
2020), we require a snapshot of the density, composition, and thermal energy
of each cell in a Cartesian velocity grid. We use the same SPH mapping
technique as Collins et al. (2023), in which particles are propagated with
constant velocity for 0.5 s after the SPH simulation terminates, except that
our grid size is $50^{3}$ (instead of $128^{3}$). The SPH densities are then
mapped onto the Cartesian velocity grid spanning -0.48c to 0.48c along each
axis, from which homologous expansion is assumed (i.e., $\rho\propto t^{-3}$,
see Neuweiler et al. 2023 for a study on the homologous assumption for
dynamical NSM ejecta).

The first few minutes after a NSM involve a complex set of r-process
reactions, so we use a nuclear network calculation (as employed previously by
Collins et al., 2023) to model the early heating and composition evolution.
The network calculation is run separately along the trajectory of each SPH
particle, starting once the temperature falls below 10 GK and continuing
beyond the end of the SPH simulation by assuming homologous expansion.

Synthetic light curves and spectra may depend sensitively on the assumptions
and approximations used for the heating and thermalization treatment (e.g.,
Barnes et al., 2016; Bulla, 2023). Although the network calculation gives
summarized energy production rates that can be used for approximate heating
(e.g., Collins et al., 2023), a detailed treatment of thermalization requires
complete information about the energies of emitted $\gamma$-rays and charged
particles, so we extend artis to treat individual $\alpha$ and $\beta^{-}$
decay reactions. We take a snapshot of the network calculations at 0.1 days to
ensure that capture reactions have slowed sufficiently such that only $\alpha$
and $\beta^{-}$ decay reactions are significant (see Appendix B for comparison
with the full network calculations). We then map the partial density of each
nuclide from the particles to the Cartesian grid in a similar way as the total
mass density. The network calculations also give the integrated energy
released in each cell (excluding neutron losses and assuming adiabatic
expansion) from all reactions prior to the snapshot time, which artis uses to
seed thermal energy packets and set the initial cell temperatures.

For comparison with the full 3D model, we also construct a 1D ejecta model by
spherically-averaging and binning the 3D ejecta into 25 radial velocity shells
(uniformly spaced from $v_{rad}=0$ to $v_{rad}=0.48c$) with densities chosen
such that the mass, as well as the partial masses of individual nuclides, are
preserved within each shell.

###  2.2 Thermalization of decay products

The tracking of individual decay reactions in artis enables the use of per-
decay tabulated $\gamma$-ray spectra (in contrast to the single averaged
$\gamma$ emission spectrum used by Collins et al. 2023), as well as a new
time-dependent treatment of $\alpha$ and $\beta^{-}$ particle thermalization
that uses emitted particle energies adapted to the local conditions in the
ejecta. The models presented here include 2,591 nuclides with $\alpha$ and
$\beta^{-}$ reaction data ($\gamma$-ray spectra and mean emitted particle
energies) from ENDV/B-VII.1 (Chadwick et al., 2011)333As provided by
https://github.com/hotokezaka/HeatingRate.

For $\gamma$-rays, we use the same frequency-dependent transport scheme
previously applied to Type Ia calculations (Lucy, 2005; Sim, 2007). For
$\alpha$ and $\beta^{-}$ particles, we assume that energy deposition takes
place at the location of the particle emission, which is a plausible
assumption given that magnetic fields may prevent particles from crossing a
significant fraction of the ejecta before losing their energy (Barnes et al.,
2016).

Radioactively emitted particles lose energy over time to the ejecta until they
reach equilibrium with the thermal pool. We adopt the approximation by Barnes
et al. (2016) for the energy loss rates of $\alpha$\- and $\beta$-particles of
$\dot{E}_{\alpha}=5\times 10^{11}\rho$ and $\dot{E}_{\beta}=4\times
10^{10}\rho$ [MeV/s], where $\rho$ is the density in units of g/cm3. Although
this loss rate function is independent of energy, the loss rate is defined per
particle, and therefore the total energy deposition rate depends on the number
of particles yet to thermalize, which requires knowledge of individual
radioactive decay rates and their associated emitted particle energies. Even
for the same total $\beta^{-}$ particle energy emission rate, energy from a
large number of reactions that emit low-energy particles will thermalize on a
shorter timescale than a slower rate of reactions that emit high-energy
particles. This treatment captures the local dependence of thermalization on
both the 3D density structure and the composition (which determines the decay
rates and energy spectrum of emitted particles) that are neglected in
analytical formulae provided by Barnes et al. (2016), which are based on a
uniform density model and used by the majority of existing kilonova models. In
Appendix C, we show a comparison between our calculated thermalization
efficiencies and the Barnes et al. (2016) analytical formulae.

###  2.3 Radiative transfer and atomic data

We consider two different sets of atomic data, AD1 and AD2. The full atomic
dataset for the AD1 models includes 71 elements, 283 ions, 194,035 energy
levels and a total of 43.6 million bound-bound transitions. Levels and
transitions for elements (C, O, Ne, Mg, Si, S, Ar, Ca, Fe, Co, Ni) are sourced
from the CMFGEN compilation444Available at
http://kookaburra.phyast.pitt.edu/hillier/web/CMFGEN.htm (Hillier, 1990;
Hillier & Miller, 1998), similar to the dataset used for the nebular Type Ia
models by Shingles et al. (2020, 2022) except that photoionization is not
included. For elements Cu to Ra, we incorporate all ions and levels of the
Japan-Lithuania database (Tanaka et al., 2020), which covers ionization stages
from neutral to triply-ionized. We do not include any atomic data for
actinides, which we will investigate in a future study.

Since the level energies in the Japan-Lithuania database are not calibrated to
known transition lines, the low-lying permitted Sr II 5p $\rightarrow$ 4d
transitions (likely identified in AT2017gfo by Watson et al. 2019, but see
Tarumi et al. 2023 for an alternative explanation for the same feature with
He) occur at rest wavelengths of 17903, 18190, and 20291 Å instead of 10036,
10327, and 10914 Å (Sansonetti, 2012). To address this, we have separately
produced calculations with a different atomic dataset called AD2 that is the
same as AD1, except that the data for Sr, Y, and Zr are from the Kurucz (2018)
extended line list, which includes calibrated wavelengths where available.

We simulate a time range of 0.1 to 80 days for 3D AD2 (for the other models,
0.1 to 20 days), although the validity of our synthetic observables (light
curves and spectra) is limited by light travel time effects to between 0.18
and 13.6 days (0.18 to 3.4 days for other models, see Appendix A for details).

The radiative transfer calculations assume local thermodynamic equilibrium
(LTE) with Boltzmann level populations and ionization balance according to the
Saha equation. The LTE assumption should be reasonably accurate within the
first few days (Pognan et al. 2022a, but see Tarumi et al. 2023 for possible
NLTE effects on Sr and He) and beyond this time, we only consider the
luminosity to be reliable, since it simply converges to the deposition rate in
optically thin ejecta and no longer has a dependence on the excitation or
ionization state.

##  3 Results

###  3.1 Luminosity evolution

Figure 1: Direction-integrated luminosity versus time for the models 3D AD1,
3D AD2, 1D AD1, 1D AD2, the 3D gray opacity model of Collins et al. (2023),
and inferred bolometric luminosity of AT2017gfo (Smartt et al., 2017).

Figure 1 shows the spherically-averaged light curves for all four models with
3D or 1D ejecta and AD1 or AD2 atomic datasets, the Collins et al. (2023) 3D
gray opacity calculations for the same NSM ejecta, as well as the inferred
bolometric (all wavelengths) flux for AT2017gfo (Smartt et al., 2017).

Compared to AT2017gfo, the bolometric luminosities of all of the models are
lower by about a factor of ten, mainly because our model mass includes only
early ejecta with a mass of 0.005 M⊙, which is about ten times lower than the
mass estimated for AT2017gfo (Smartt et al., 2017).

Both the AD1 and AD2 atomic datasets lead to approximately equal total
bolometric luminosities, which indicates that wavelength calibration of Sr, Y,
and Zr does not significantly affect the overall optical depth of the ejecta.

The luminosities of the 3D models are higher than for the 1D models until
about two days. This is a consequence of the strong anisotropies in our NSM
model. A complex 3D ejecta structure with local opacity variations exhibits
pathways of lower optical depth and therefore allows light to effectively
escape faster than in a 1D configuration. After a couple of days, when the
expanding ejecta have become optically thin, the luminosities of the 1D and 3D
models converge to the rate of thermalized beta particle energy.

The luminosities of the 3D AD1 and 3D AD2 models agree reasonably well with
the 3D gray opacity calculation of Collins et al. (2023) for the first couple
of days. This shows that the $Y_{e}$-dependent gray opacity tables of Tanaka
et al. (2020) are a good approximation to the opacities calculated line-by-
line with the time- and location-dependent temperatures and composition of our
model. At times later than about two days, the thermalization efficiency for
$\beta$ particles begins to decrease significantly and the luminosity of our
model decreases more rapidly than the Collins et al. (2023) model, which
assumes fully-efficient $\beta^{-}$ thermalization.

Overall, we find that the use of 3D ejecta is critical for modeling the early
(<2 days, while optically thick) luminosity evolution, and the opacity
treatment is comparatively less important. This is also approximately the time
range within which the assumption of fully-efficient thermalization of $\beta$
particles is valid.

###  3.2 Observing direction maps

Figure 2: Mollweide projections of direction-dependent quantities for 3D AD2
UVOIR packets arriving at the observer between 1.3 and 1.7d: radiant intensity
times 4$\pi$ solid angle, mean temperature at last interaction, and line of
sight velocity at last interaction. For these figures, we use 32x32 direction
bins, uniformly spaced in azimuthal angle (horizontal) and cosine of the polar
angle (vertical) to give the same solid angle in each bin.

With 3D models, we are able to explore the direction dependence of the
synthetic observables. Figure 2 shows Mollweide projections for the isotropic-
equivalent luminosity (radiant intensity times $4\pi$), the mean temperature
of the emission locations, and the mean line-of-sight velocity of emitting
locations for light arriving at the observer between 1.3 and 1.7 days.

The faintest lines of sight are those viewing toward the equator (merger
plane), where the model has the lowest $Y_{e}$ and a more lanthanide-rich
composition. The polar-angle variation of the composition (e.g., Kullmann et
al., 2022; Collins et al., 2023) leads to a higher mean last-interaction
ejecta velocity for equatorial observers, where the ejecta density and
temperature are lower than lower-velocity central regions. Compared to polar
observers, equatorial observers receive emission originating from regions with
a higher average line-of-sight velocity (0.16–0.17c versus 0.10–0.13c), and
lower temperature (2500–2800 K versus 3200-3400 K).

Apart from the variation with polar angle, which was already found in previous
studies (e.g., Just et al., 2022; Neuweiler et al., 2023), we also see
substantial variation in the emission as a function of the azimuthal angle.
Since most existing multi-dimensional kilonova studies assume 2D axisymmetry,
they cannot capture the dependence on azimuthal angle. Our results suggest
that this dependence can be important as well, at least for merger models
similar to the one considered here.

###  3.3 Spectra toward pole and equator

Figure 3: Spectra for polar and equatorial viewing directions for the 3D AD1
and 3D AD2 models at 0.8 days. The height of each wavelength point is colored
according to the emitting species of the last interactions of the emerging
radiation packets. The area under the horizontal axis shows the distribution
of frequencies (colored by absorbing/scattering ion) just prior to the last
interactions of the emerging packets. The 11 most-significant ions are
separately colored, while the ”Other” group combines many smaller
contributions from other ions.

Figure 3 shows spectra at 0.8 days for the 3D AD1 and 3D AD2 models for polar
and equatorial inclinations (averaged over azimuthal angle to reduce Monte
Carlo noise). At each wavelength, we divide the flux into the set of emitting
ions based on the last interactions of the emerging radiation packets. Below
the horizontal axis, we show the distribution of absorption frequencies (in
the observer frame) at the last interaction and associate this flux with the
set of absorbing/scattering ions. This visualization gives some indication of
the ions that are important for shaping the spectrum, although it does not
consider the full history of interactions experienced by the emerging packets.

We notice substantial variation in the spectra with both the polar viewing
angle as well as the adopted atomic dataset. For both 3D models, the
equatorial spectra are fainter and show comparatively few pronounced features
compared to the spectra for the near-polar inclination. In contrast to
supernova spectra, for which individual spectral features can often be
attributed to the transitions of one or two ions, our model results suggest
that the spectral formation in kilonovae is much more complex. Across the
wavelength range, the emerging packets of radiation have interacted with
multiple ions. Fluorescence by several elements (especially Sr, Y, Zr, and Ce)
is highly effective at redistributing flux bluewards of $\sim$4000 Å to longer
(infrared) wavelengths.

The effect of wavelength calibration is dramatic for the Sr II contributions
(see subsection 2.3 for details of the rest wavelength differences): In the 3D
AD1 spectra, an extremely broad Sr II triplet feature exists around
15000-20000 Å, while in the 3D AD2 model, the same feature occurs at around
7500-12500 Å, blue-shifted from the rest wavelengths of around 10000 Å. In the
3D AD2 case, the wavelength-calibrated Sr II transitions occur much closer to
the peak of the continuum, and therefore contribute more strongly to the
overall spectral distribution. Absorption by Sr II (below the axis) is further
blue-shifted than the Sr II emission, as would be expected from a P Cygni
feature (Watson et al., 2019). However, we see that Sr II also emits by
fluorescence following the absorption of shorter-wavelength ($\sim$4000 Å)
radiation.

The complex interplay of ion interactions is also indicated by the degree to
which the calibration of Sr, Y, and Zr has consequences for the interactions
of other ions (e.g., Pr III) that have identical atomic data in the AD1 and
AD2 models.

We also see that Ce III is often the last interaction of packets emerging
around 10000-12000 Å, while Domoto et al. (2022) have shown that calibration
would substantially weaken Ce III lines in this range. However, despite
frequent interactions with Ce III, a comparison with a separate simulation
performed without Ce (but otherwise identical to 3D AD2) suggests that the Ce
only has a modest impact on the shape of the spectral distribution, and
interacts mainly by scattering, with some fluorescence of optical flux to
infrared wavelengths.

###  3.4 Spectral evolution

Figure 4: Time series of spectra in the polar direction of the 3D AD2 model
compared to reddening and redshift corrected spectra of AT2017gfo (Pian et
al., 2017; Smartt et al., 2017). The area under the spectra have been coloured
by the emitting species of the last interactions of the emerging packets. The
times of the artis and AT2017gfo spectra intentionally do not match.

Figure 4 shows a time series of spectra for the polar direction of the 3D AD2
model and a series of AT2017gfo observations555Engrave data release available
at http://www.engrave-eso.org/AT2017gfo-Data-Release/ (Pian et al., 2017;
Smartt et al., 2017). We select a viewing angle bin closest to the +z pole
($\pi$-$\theta\leq 37$°, where $\theta$ is the angle from the positive z-axis,
as defined by Collins et al. 2023) containing the estimated viewing angle to
AT2017gfo of between 14° and 28° (Mooley et al., 2022).

Although the model evolves about 3-4 times quicker than AT2017gfo, the 3D AD2
spectra pass through phases that are remarkably similar to the observed
spectra. The evolution being too rapid could be a consequence of the low-mass
of our model that includes only the early ejecta, as a larger mass would
increase the diffusion timescale.

The 3D AD2 polar spectrum at 0.4 days has a single broad bump at around
4000-7000 Å, which is similar to AT2017gfo at 1.4 days. In our model, the
broad peak is a blend, where the Monte Carlo packets have had notable
interactions with Zr II, Y II, Nd III, and Pr III, as well as interactions
with many other ions.

At 0.7 days, the 3D AD2 spectrum begins to exhibit more pronounced spectral
features, with similarities to the AT2017gfo spectrum at 2.42 days. However,
the spectral features in the model occur at shorter wavelengths (by up to
$\sim$1000 Å) which might suggest that the expansion velocities in the ejecta
model are too high.

At 0.8 days, the distribution has shifted further to the red, with stronger
features and a similar spectrum to AT2017gfo at 3.41 days. There are spectral
features in AT2017gfo redward of 12000 Å that are not reproduced by our model
(possibly due to our use of uncalibrated La III and Ce III lines, see Domoto
et al. 2022 and Tanaka et al. 2023), but the approximate agreement with the
polar spectrum is interesting when contrasted with the predicted spectra in
the equatorial direction. The equatorial spectrum of 3D AD2 at 0.8 days (shown
in Figure 3) is comparatively lacking in pronounced spectral features and is
much less similar to the AT2017gfo spectrum. The closer spectral match in the
polar direction of our model independently supports the polar inclination of
AT2017gfo, which has been previously inferred with other methods that involve
different assumptions about the merger (Mooley et al., 2022).

At 0.9 days, the 3D AD2 spectral energy distribution shifts further toward
redder wavelengths, with the overall spectrum appearing similar to AT2017gfo
at 4.40 days.

These results show that forward modeling of a merger simulation that has not
been tuned to match AT2017gfo, nevertheless shows remarkable similarities with
observations when viewed in the polar direction.

###  3.5 Spherically-averaged ejecta models

Figure 5: Spherically-averaged spectra at 0.8 days for the 3D AD1 (solid
blue), 3D AD2 (solid orange), 1D AD1 (dashed blue), 1D AD2 (dashed orange)
models.

Figure 5 shows the spectra at 0.8 days averaged over all viewing directions
for the models calculated in either full 3D, or with 1D spherically-averaged
ejecta with either the AD1 or AD2 atomic datasets. Even when the same atomic
dataset is used, the 1D calculations predict very different spectral energy
distributions compared to the 3D calculations. In general, the 3D models
produce substantially more flux shortward of $\sim$8000 Å, which highlights
the need for multidimensional models when predicting spectral energy
distributions.

For the AD1 case, the 1D and 3D models produce very different sets of spectral
features. Indeed, the 3D AD1 model shows a flux maximum at $\sim$7000 Å, while
1D AD1 has a maximum at around 12000 Å. For the AD2 case, the 1D and 3D models
each produce a spectral maximum at around 10000 Å, although the 1D model also
shows a second distinct feature at around 12000–13000 Å that is not present
for the 3D model.

##  4 Discussion and Conclusions

We have presented the so-far most self-consistent radiative-transfer
calculations of NSM ejecta, taking into account the full 3D spatial and time
dependence of the density, composition, heating rate, and thermalization
efficiency, and using for the first time, a line-by-line opacity treatment.
The line-by-line treatment enables the handling of fluorescence and allows the
association of spectral features with individual atomic transitions. We have
shown that the Sr II triplet (5p $\rightarrow$ 4d) plays an important role in
shaping the spectrum around 8000 to 11000 Å, as has been suggested on the
basis of simplified models (e.g., Watson et al., 2019), which is an important
confirmation given the substantial differences between 1D and 3D model results
presented in subsection 3.5. However, we show that important contributions to
the spectra are also made by Ce III (although our Ce III lines are
uncalibrated, see Domoto et al., 2022), Zr I, Zr II, and Y II (as suggested
from observations by Sneppen & Watson, 2023).

We have produced a time series of viewing-angle-dependent spectra that, in the
polar direction, are similar to AT2017gfo, and which show distinctly more
features compared to equatorial spectra. These results independently support
the inferred polar inclination of AT2017gfo and suggest the possibility that
distinctly different spectra will be obtained from future kilonovae observed
at different inclination angles. While we can only broadly constrain the
inclination with our current model, future theoretical models offer the
possibility of further constraining the inclination, which is highly desirable
for gravitational-wave data analysis, distance estimates, Hubble constant
inference, and the interpretation of NSM observations in general.

We obtain dramatically different spectra between the 1D and the 3D ejecta
models derived from the same NSM model, despite keeping all other aspects of
the radiative transfer method identical. In particular, there are strong
observer-orientation effects (see Figure 3), and moreover, the spectra from a
spherically-averaged 1D model do not resemble those of the underlying 3D model
in any viewing angle. We contend that this demonstrates the need for
multidimensional models of kilonovae, even when the direction-dependence of
the observables is not required. This has important implications for
interpreting analyses based on 1D empirical models (e.g., Gillanders et al.,
2022).

While the use of advanced treatments of radioactive decay and thermalization,
atomic absorption and emission help to reduce systematic uncertainties in the
radiative transfer calculations, the results of our study also highlight the
importance of using accurate atomic data. We have shown the importance of
calibrating energy levels in atomic structure calculations to observed
transition wavelengths, with major differences in the resulting synthetic
spectra being produced when calibrated atomic data are used for Sr, Y, and Zr.
There are many heavy ions for which no calibrated atomic data are published,
while work in this direction is ongoing (e.g., Flörs et al., 2023). Future
applications of calibrated data to radiative transfer calculations are likely
to help in explaining additional features of kilonova spectra and correlating
these with merger dynamics and remnant properties.

The authors thank the anonymous referee for providing helpful comments, as
well as Brian Metzger for providing useful feedback. AF, GL, GMP, LJS, and ZX
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). AB, OJ, CEC, and VV 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, CEC, AF, OJ, GL, GMP,
LJS, and ZX acknowledge support by Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) - Project-ID 279384907 - SFB 1245 and MA 4248/3-1.
AB and VV acknowledge support by DFG - Project-ID 138713538 - SFB 881 (“The
Milky Way System”, subproject A10). AB, CEC, AF, GL, GMP, OJ, LJS, and ZX
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, ST/X00094X/1]. 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. The authors gratefully acknowledge the
Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this
project by providing computing time through the John von Neumann Institute for
Computing (NIC) on the GCS Supercomputer JUWELS (Jülich Supercomputing Centre,
2021) at Jülich Supercomputing Centre (JSC) and the VIRGO cluster at GSI for
computational support.

##  Appendix A ARTIS calculation and time range of synthetic observables

We use artis version v2023.5.16.3 (Sim et al., 2023). The 3D AD2 model (0.1 to
80 days) with $10^{8}$ Monte Carlo energy packets used 114 hours of wallclock
time on 40 nodes with 48 Intel Skylake CPU cores (220 kilo-core-hours). The
memory usage was about 2.8 GB per core. The 3D AD1 model (0.1 to 20 days) with
$5\times 10^{7}$ packets on the same node configuration used 128 kilo-core-
hours.

We simulate a time range of 0.1 to 80 days for 3D AD2 (for the other models,
0.1 to 20 days), although the range of validity of our synthetic observables
is limited by light travel time to between 0.18 and 13.6 days in 3D AD2 (0.18
and 3.4 days for other models).

Packets of radiation can escape from the simulation volume up to the end of
the simulated time range, but the escape surface at the boundary of the
simulation expands with constant velocity ($0.48c$ per axis, $0.83c$ at the
corners), and packets might have already traveled some distance toward the
observer before escaping. The light-crossing time of the component of distance
along the line of sight to the observer ($\vec{r}\cdot\hat{v}/c$) must be
subtracted from the escape time to get a consistent time of arrival for the
observer that is independent of the arbitrary escape surface.

The greatest time correction will occur for a packet moving in a radial
direction and escaping at a corner of the simulation at the end time
($|\vec{r}_{\mathrm{corner}}|=0.83c\times 80$ days). In this case, the
observer arrival time is less than the escape time by
$(\vec{r}\cdot\hat{v})/c=0.83c\times 80$ days $=66.4$ days. The latest
observer time for which we have complete coverage of escaping packets is then
$80-66.4=13.6$ days for 3D AD2 (3.4 days for other models).

For the earliest observation times, a packet emitted radially inwards from a
simulation corner must have time to reach the origin (beyond which further
travel time will be subtracted), which is $0.8c\times 0.1$ days, or 0.08 days
after the 0.1 day simulation start time. Therefore the synthetic observables
are valid from 0.18 to 13.6 days (or to 3.4 days), while internal quantities
such as deposition and temperatures are defined for the entire simulation time
range.

##  Appendix B Validation of the ARTIS nuclear decay treatment

The artis treatment of $\alpha$ and $\beta^{-}$ decays uses a fast analytical
calculation. With the assumption of pure decays (i.e., no capture reactions),
there are no loops in the nuclear network, and therefore we can enumerate each
of the finite number of paths that begin at a nuclide in the snapshot
composition. Given the snapshot abundances at time $t_{s}$, a modified Bateman
(1910) equation (incorporating branching factors for each decay type) gives
each path’s contribution to the abundances of the nuclide at the end of the
path at any time $t>t_{s}$. The abundance of any nuclide is then given by the
sum of the contributions from all associated decay paths.

Here, we show that (with our chosen snapshot time of $t_{s}=0.1$ days) the
decay treatment in artis maintains close agreement with the iterative matrix-
based calculation that includes capture reactions and spontaneous fission.

Figure 6: Comparison of energy release rates of the artis treatment in 3D AD2
versus the output of the full nuclear network (GSI network). Neglecting
fission reactions has only a negligible effect on the total energy release
rates.

Figure 6 shows the global energy release rates as a function of time for both
artis and the full nuclear network calculation. The artis treatment is able to
match the evolution of $\dot{Q}_{\alpha}$ and $\dot{Q}_{\beta}$ closely, which
are the total (including neutrino) energy release rates for the $\alpha$ and
$\beta^{-}$ decays, respectively. The total $\dot{Q}$ is not shown because it
is extremely close to the $\dot{Q}_{\beta}$ value.

Figure 7: Comparison between the nuclear mass fractions from the artis decay-
only treatment and the full nuclear network calculation for a cell at $z=0$
and radial velocity $v_{r}=0.24c$ with $Y_{e}=0.13$ .

Figure 7 shows a comparison of the abundance evolution for selected nuclides
in both the artis decay-only treatment and the full nuclear network
calculation. The abundances maintain close agreement throughout the time
range, and in particular, the 4He abundance increases due to many $\alpha$
decays.

##  Appendix C Deposition and thermalization

Figure 8: Left panel: synthetic direction-integrated bolometric luminosity
(all wavelengths excluding $\gamma$-rays), global deposition rate for
$\gamma$-rays, and emission and deposition rates for $\beta^{-}$ particles for
the 3D AD2 model. Right panel: global thermalization ratios (power deposited
divided by emitted) versus time for $\gamma$-rays, $\beta^{-}$, and $\alpha$
particles in the 3D AD2 artis model. We also show the Barnes et al. (2016)
analytical fits (Equations C1 and C2) for comparison.

The left panel of Figure 8 shows the direction-integrated bolometric
(ultraviolet, optical, and infrared) luminosity, as well as the global rates
of energy deposited by $\gamma$-rays, and the energy released and deposited by
$\beta$ particles for the 3D AD2 model. The right panel of Figure 8 shows the
time evolution of the thermalization ratios for $\gamma$-rays, $\beta^{-}$,
and $\alpha$ particles in artis and for comparison, the Barnes et al. (2016)
analytical fits using the total mass and average velocity of our ejecta model.

The approximate thermalization efficiencies for $\gamma$-rays and massive
particles are given by equations 32 and 33 of Barnes et al. (2016),

$\displaystyle f_{p}(t)$
$\displaystyle=\frac{\ln\left[1+2\left(\frac{t}{t_{\text{ineff,p}}}\right)^{2}\right]}{2\left(\frac{t}{t_{\text{ineff,p}}}\right)^{2}}\quad\text{and}$
(C1) $\displaystyle f_{\gamma}(t)$
$\displaystyle=1-\exp\left[-\left(\frac{t}{t_{\text{ineff},\gamma}}\right)^{-2}\right],$
(C2)

respectively, where $t_{\text{ineff}}$ refers to the time at which
$\gamma$-rays or particles fail to thermalize efficiently (when the
thermalization timescale exceeds the expansion time).

The thermalization timescales are defined by equations 16, 17, 20, 25 of
Barnes et al. (2016):

$\displaystyle t_{\text{ineff},\gamma}$
$\displaystyle=1.4M_{5}^{1/2}v_{2}^{-1}\text{ days}$ (C3) $\displaystyle
t_{\text{ineff},\beta}$ $\displaystyle=7.4\left(\frac{E_{\beta,0}}{0.5\text{
MeV}}\right)^{-1/2}M_{5}^{1/2}v_{2}^{-3/2}\text{ days}$ (C4) $\displaystyle
t_{\text{ineff},\alpha}$ $\displaystyle=4.3\times
1.8\left(\frac{E_{\alpha,0}}{6\text{
MeV}}\right)^{-1/2}M_{5}^{1/2}v_{2}^{-3/2}\text{ days},$ (C5)

where $M_{5}$ is the ejecta mass in units of $10^{-5}$ ${\rm M}_{\sun}$, and
$v_{2}$ is the average velocity computed according to
$v_{2}=\sqrt{\int{v^{2}dm}/\int{dm}}/(0.2c)$. For the analytic values
presented here, we use the suggested initial particle energies of
$E_{\beta,0}=0.5$ MeV and $E_{\alpha,0}=6$ MeV and, although the Barnes
formulae are not intended to be applied to a 3D density structure, we use the
total mass $M_{5}$ and average velocity $v_{2}$ of our simulated ejecta model.

The majority of the energy released after 0.1 days is in the form of
$\gamma$-rays. However, the ejecta trap only 60 percent of this energy at 0.1
days and progressively become optically thin to $\gamma$-rays, with the
thermalization efficiency decreasing to about 10 percent at one day. The form
of Equation C2 requires thermalization ratios to start from 1.0 at time $t=0$.

Since we assume that $\beta$ particles remain locally trapped, they do not
diffuse out from the ejecta, and their thermalization efficiency decreases
more gradually. The efficiency of $\beta$ thermalization does decrease with
the ongoing dilution of the expanding ejecta, causing emitted particles to
deposit their kinetic energy on longer timescales (Kasen & Barnes, 2019).
After a couple of hours, the primary contribution to energy deposition is from
$\beta$ particles ($\alpha$ particles make a much smaller contribution). The
thermalization remains efficient out to timescales of days, only dropping
below 90 percent at around two days.

We find that our thermalization ratios for $\alpha$ and $\beta$ particles are
generally close to the analytical approximations over the first two days (less
than about 5 percent difference) and agree fairly well (up to 20 percent
difference) up to about five days. Compared to the analytical formulae, which
predict similar thermalization ratios of $\alpha$ and $\beta$ particles, we
find that $\beta$ particles are more efficiently thermalized than $\alpha$
particles. This is likely related to differences in the composition-dependent
rates of individual $\alpha$ and $\beta^{-}$ decay reactions, which each emit
particles with different energies and therefore have different thermalization
timescales.

## References

  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 848, L12
  * 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, The Astrophysical Journal, 901, 29
  * Barnes et al. (2016) Barnes, J., Kasen, D., Wu, M.-R., & Martinez-Pinedo, G. 2016, ApJ, 829, 110
  * Bateman (1910) Bateman, H. 1910, Proc. Cambridge Philos. Soc., 15, 423
  * Bauswein et al. (2010) Bauswein, A., Oechslin, R., & Janka, H. T. 2010, Physical Review D, 81, 024012
  * Bulla (2019) Bulla, M. 2019, MNRAS, 489, 5037
  * Bulla (2023) —. 2023, MNRAS, 520, 2558
  * Chadwick et al. (2011) Chadwick, M. B., Herman, M., Obložinský, P., et al. 2011, Nuclear Data Sheets, 112, 2887
  * Collins et al. (2023) Collins, C. E., Bauswein, A., Sim, S. A., et al. 2023, MNRAS, 521, 1858
  * Domoto et al. (2022) Domoto, N., Tanaka, M., Kato, D., et al. 2022, ApJ, 939
  * Domoto et al. (2021) Domoto, N., Tanaka, M., Wanajo, S., & Kawaguchi, K. 2021, ApJ, 913, 26
  * Even et al. (2020) Even, W., Korobkin, O., Fryer, C. L., et al. 2020, ApJ, 899, 24
  * Flörs et al. (2023) Flörs, A., Silva, R. F., Deprince, J., et al. 2023, arXiv, arXiv:2302.01780
  * Fontes et al. (2020) Fontes, C. J., Fryer, C. L., Hungerford, A. L., Wollaeger, R. T., & Korobkin, O. 2020, MNRAS, 493, 4143
  * Gillanders et al. (2022) Gillanders, J. H., Smartt, S. J., Sim, S. A., Bauswein, A., & Goriely, S. 2022, MNRAS, 515, 631
  * Hillier (1990) Hillier, D. J. 1990, Astronomy and Astrophysics (ISSN 0004-6361), 231, 116
  * Hillier & Miller (1998) Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407
  * Hunter (2007) Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90
  * Isenberg & Nester (1980) Isenberg, J., & Nester, J. 1980, General Relativity and Gravitation. Vol. 1. One hundred years after the birth of Albert Einstein. Edited by A. Held. New York, 1, 23
  * Jülich Supercomputing Centre (2021) Jülich Supercomputing Centre. 2021, Journal of large-scale research facilities, 7
  * Just et al. (2022) Just, O., Kullmann, I., Goriely, S., et al. 2022, MNRAS, 510, 2820
  * Kasen et al. (2013) Kasen, D., Badnell, N. R., & Barnes, J. 2013, ApJ, 774, 25
  * Kasen & Barnes (2019) Kasen, D., & Barnes, J. 2019, ApJ, 876, 128
  * Kasen et al. (2015) Kasen, D., Fernandez, 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, 192, L145
  * Kato et al. (2021) Kato, D., Murakami, I., Tanaka, M., et al. 2021, Japan-Lithuania Opacity Database for Kilonova (version 1.1)
  * Kawaguchi et al. (2022) Kawaguchi, K., Fujibayashi, S., Hotokezaka, K., Shibata, M., & Wanajo, S. 2022, The Astrophysical Journal, 933, 22
  * Kawaguchi et al. (2021) Kawaguchi, K., Fujibayashi, S., Shibata, M., Tanaka, M., & Wanajo, S. 2021, The Astrophysical Journal, 913, 100
  * Kawaguchi et al. (2018) Kawaguchi, K., Shibata, M., & Tanaka, M. 2018, The Astrophysical Journal Letters, 865, L21
  * Korobkin et al. (2021) Korobkin, O., Wollaeger, R. T., Fryer, C. L., et al. 2021, ApJ, 910, 116
  * Kromer & Sim (2009) Kromer, M., & Sim, S. A. 2009, MNRAS, 398, 1809
  * Kullmann et al. (2022) Kullmann, I., Goriely, S., Just, O., et al. 2022, MNRAS, 510, 2804
  * Kurucz (2018) Kurucz, R. L. 2018, Workshop on Astrophysical Opacities. ASP Conference Series, 515, 47
  * Lucy (2002) Lucy, L. B. 2002, A&A, 384, 725
  * Lucy (2005) —. 2005, A&A, 429, 19
  * Metzger et al. (2010) Metzger, B. D., Martínez-Pinedo, G., Darbha, S., et al. 2010, MNRAS, 406, 2650
  * Mooley et al. (2022) Mooley, K. P., Anderson, J., & Lu, W. 2022, Nature, 610, 273
  * Neuweiler et al. (2023) Neuweiler, A., Dietrich, T., Bulla, M., et al. 2023, Physical Review D, 107
  * Oechslin et al. (2002) Oechslin, R., Rosswog, S., & Thielemann, F.-K. 2002, Physical Review D, 65, 103005
  * Oliphant (2007) Oliphant, T. E. 2007, Computing in Science & Engineering, 9, 10
  * Pian et al. (2017) Pian, E., D’Avanzo, P., Benetti, S., et al. 2017, Nature, 213, 225
  * Pognan et al. (2022a) Pognan, Q., Jerkstrand, A., & Grumer, J. 2022a, MNRAS, 513, 5174
  * Pognan et al. (2022b) —. 2022b, MNRAS, 510, 3806
  * Sansonetti (2012) Sansonetti, J. E. 2012, Journal of Physical and Chemical Reference Data, 41, 013102
  * Shingles et al. (2023) Shingles, L. J., Collins, C. E., Holas, A., Callan, F., & Sim, S. 2023, artistools, v2023.5.16.3, doi: 10.5281/zenodo.7942849
  * Shingles et al. (2022) Shingles, L. J., Flörs, A., Sim, S. A., et al. 2022, MNRAS, 512, 6150
  * Shingles et al. (2020) Shingles, L. J., Sim, S. A., Kromer, M., et al. 2020, MNRAS, 492, 2029
  * Sim et al. (2023) Sim, S., Kromer, M., Shingles, L., et al. 2023, artis-mcrt/artis: v2023.5.2, v2023.5.2, Zenodo, doi: 10.5281/zenodo.7981387
  * Sim (2007) Sim, S. A. 2007, MNRAS, 375, 154
  * Smartt et al. (2017) Smartt, S. J., Chen, T. W., Jerkstrand, A., et al. 2017, Nature, 551, 75
  * Sneppen & Watson (2023) Sneppen, A., & Watson, D. 2023, arXiv, arXiv:2306.14942
  * Sneppen et al. (2023) Sneppen, A., Watson, D., Bauswein, A., et al. 2023, Nature, 614, 436
  * Steiner et al. (2013) Steiner, A. W., Hempel, M., & Fischer, T. 2013, The Astrophysical Journal, 774, 17
  * Tanaka & Hotokezaka (2013) Tanaka, M., & Hotokezaka, K. 2013, ApJ, 775, 113
  * Tanaka et al. (2020) Tanaka, M., Kato, D., Gaigalas, G., & Kawaguchi, K. 2020, MNRAS, 496, 1369
  * Tanaka et al. (2023) Tanaka, M., Domoto, N., Aoki, W., et al. 2023, arXiv, arXiv:2306.04697
  * Tarumi et al. (2023) Tarumi, Y., Hotokezaka, K., Domoto, N., & Tanaka, M. 2023, arXiv, arXiv:2302.13061
  * Villar et al. (2017) Villar, V. A., Guillochon, J., Berger, E., et al. 2017, The Astrophysical Journal Letters, 851, L21
  * Watson et al. (2019) Watson, D., Hansen, C. J., Selsing, J., et al. 2019, Nature, 574, 497
  * Wilson et al. (1996) Wilson, J. R., Mathews, G. J., & Marronetti, P. 1996, Physical Review D (Particles, 54, 1317
  * Wollaeger et al. (2021) Wollaeger, R. T., Fryer, C. L., Chase, E. A., et al. 2021, ApJ, 918, 10