# Spectroscopic $r$-Process Abundance Retrieval for Kilonovae II: Lanthanides
in the Inferred Abundance Patterns of Multi-Component Ejecta from the GW170817
Kilonova

Nicholas Vieira Trottier Space Institute at McGill and Department of Physics,
McGill University, 3600 rue University, Montreal, Québec, H3A 2T8, Canada
John J. Ruan Department of Physics and Astronomy, Bishop’s University, 2600
rue College, Sherbrooke, Québec, J1M 1Z7, Canada  Daryl Haggard Trottier Space
Institute at McGill and Department of Physics, McGill University, 3600 rue
University, Montreal, Québec, H3A 2T8, Canada  Nicole M. Ford Trottier Space
Institute at McGill and Department of Physics, McGill University, 3600 rue
University, Montreal, Québec, H3A 2T8, Canada  Maria R. Drout David A. Dunlap
Department of Astronomy and Astrophysics, University of Toronto, 50 St. George
St., Toronto, Ontario, M5S 3H4, Canada  Rodrigo Fernández Department of
Physics, University of Alberta, Edmonton, Alberta, T6G 2E1, Canada

###### Abstract

In kilonovae, freshly-synthesized $r$-process elements imprint features on
optical spectra, as observed in AT2017gfo, the counterpart to the GW170817
binary neutron star merger. However, measuring the $r$-process compositions of
the merger ejecta is computationally challenging. Vieira et al. (2023)
introduced Spectroscopic $r$-Process Abundance Retrieval for Kilonovae
(SPARK), a software tool to infer elemental abundance patterns of the ejecta,
and associate spectral features with particular species. Previously, we
applied SPARK to the 1.4 day spectrum of AT2017gfo and inferred its abundance
pattern for the first time, characterized by electron fraction $Y_{e}=0.31$, a
substantial abundance of strontium, and a dearth of lanthanides and heavier
elements. This ejecta is consistent with wind from a remnant hypermassive
neutron star and/or accretion disk. We now extend our inference to spectra at
2.4 and 3.4 days, and test the need for multi-component ejecta, where we
stratify the ejecta in composition. The ejecta at 1.4 and 2.4 days is
described by the same single blue component. At 3.4 days, a new redder
component with lower $Y_{e}=0.16$ and a significant abundance of lanthanides
emerges. This new redder component is consistent with dynamical ejecta and/or
neutron-rich ejecta from a magnetized accretion disk. As expected from
photometric modelling, this component emerges as the ejecta expands, the
photosphere recedes, and the earlier bluer component dims. At 3.4 days, we
find an ensemble of lanthanides, with the presence of cerium most concrete.
This presence of lanthanides has important implications for the contribution
of kilonovae to the $r$-process abundances observed in the Universe.

Nuclear abundances (1128) — R-process (1324) — Radiative transfer simulations
(1967) — Spectral line identification (2073)

††journal: ApJ††software:  approxposterior: Fleming & VanderPlas (2018);
astropy: Astropy Collaboration et al. (2018); cmasher: van der Velden (2020);
corner: Foreman-Mackey (2016); dynesty: Speagle (2020); george: Ambikasaran et
al. (2015); TARDIS: Kerzendorf & Sim (2014); Kerzendorf et al. (2023);
UltraNest: Buchner (2021)

##  1 Introduction

Approximately half of the elements in the Universe heavier than iron are
synthesized by rapid neutron capture nucleosynthesis: the $r$-process (see
Cowan et al. 2021 for a review). Extreme astrophysical environments—namely,
mergers of neutron stars (NS-NS) or a neutron star and black hole (NS-BH) and
other proposed sources like collapsars or magnetorotational supernovae—offer
leading candidate sites for this $r$-process nucleosynthesis due to their
exceptionally high densities of free neutrons (Lattimer & Schramm 1974;
Symbalisty & Schramm 1982; Eichler et al. 1989; Freiburghaus et al. 1999;
Goriely et al. 2011; Korobkin et al. 2012; Bauswein et al. 2013). However, it
is still unclear which of these channels dominates. The NS-NS merger GW170817,
first detected in gravitational waves and then across the electromagnetic
spectrum (Abbott et al. 2017a, b), has provided some insight. Both photometry
(Andreoni et al. 2017; Arcavi et al. 2017; Coulter et al. 2017; Díaz et al.
2017; Drout et al. 2017; Evans et al. 2017; Hu et al. 2017; Kasliwal et al.
2017; Lipunov et al. 2017; Tanvir et al. 2017; Troja et al. 2017; Utsumi et
al. 2017; Valenti et al. 2017)111See Villar et al. (2017) for a compilation of
this photometry considering inter-instrument variation. and spectroscopy
(Chornock et al. 2017; Kasen et al. 2017; Pian et al. 2017; Shappee et al.
2017; Smartt et al. 2017) of the optical/near-infrared counterpart, AT2017gfo,
matched theoretical expectations for a kilonova: an explosive transient event
powered by radioactive decay of freshly-synthesized $r$-process elements.
However, we do not yet know the precise abundance pattern of the $r$-process
elements in the ejecta, nor whether GW170817-like kilonovae could yield the
$r$-process abundances seen across the Universe (Ji et al. 2019; Cowan et al.
2021).

The spectra of kilonovae in particular are marked by absorption and emission
features from a suite of $r$-process elements, and are key for determining the
detailed composition of the merger ejecta. Insights gained from spectra are
independent of and complementary to light curve modelling, which has served
mostly to infer macroscopic properties of the kilonova such as the heating
rates, total ejecta masses, average ejecta velocities, and temperatures of the
ejecta (e.g., Villar et al. 2017; Almualla et al. 2021; Breschi et al. 2021;
Ristic et al. 2022). By modelling the spectra, we can directly infer the
abundances and the conditions of $r$-process nucleosynthesis. These insights
further allow us to assess the importance of mergers versus other proposed
sites as the source of these elements, and, the physical ejection mechanisms
at play during these mergers.

Spectral modelling has already provided evidence that the GW170817 kilonova
ejecta contained $r$-process elements, and, has enabled associations of
certain absorption and/or emission features with individual species. Watson et
al. (2019) analyzed the early time, optically thick spectra of AT2017gfo, and
found the imprint of a P Cygni feature arising from Sr II (strontium, 38Sr) at
$\sim$8000 Å, at 1.4, 2.4, and 3.4 days post-merger. Domoto et al. (2021,
2022) similarly ascribe this feature to Sr II, and find tentative evidence for
doubly-ionized lanthanides La III and Ce III (lanthanum, 57La and cerium,
58Ce) in the near-infrared ($\sim$12,000-14,000 Å). Similarly, Gillanders et
al. (2022) argues for the presence of Sr II, but also ions of adjacent first
$r$-process peak elements Y II and Zr II (yttrium, 39Y and zirconium, 39Zr) at
wavelengths $\lesssim$ 6000 Å. Sneppen & Watson (2023) find that Y II is
present at 4.4 and 5.4 days, producing a P Cygni feature at $\sim$7600 Å.
Gillanders et al. (2022) also suggests the presence of a modest amount of
lanthanide material at these times. At later times, when the ejecta enters an
optically thin regime, the spectrum is dominated by emission features.
Gillanders et al. (2023) find that Ce III may produce emission at $\sim$15,800
Å and $\sim$20,700 Å beyond 3.4 days, up to 10.4 days. At $\gtrsim$ 7 days,
these may instead arise from intrinsically weak lines, with ions Te III and Te
I (tellurium, 52Te) and I II (iodine, 53I) given as the most likely candidates
(Gillanders et al. 2023; Hotokezaka et al. 2023).

While these studies have shed valuable light on some of the species present in
the ejecta of AT2017gfo, the abundances of all elements in the ejecta are not
known. In Vieira et al. (2023) (hereafter V23), we fit the spectrum of
AT2017gfo at 1.4 days post-merger, using our inference approach and software
tool Spectroscopic $r$-Process Abundance Retrieval for Kilonovae (SPARK). With
SPARK, we inferred the complete elemental abundance pattern of the ejecta. We
found that the ejecta was dominated by lighter $r$-process elements,
generating a bluer (relatively lower opacity in the UV and optical) kilonova.
This ejecta had electron fraction $Y_{e}=0.311^{+0.013}_{-0.011}$ and specific
entropy per nucleon $s/k_{\mathrm{B}}=13.6^{+4.1}_{-3.0}$, which yielded an
extremely low lanthanide fraction
$\log_{10}X_{\mathrm{lan}}={-7.03}^{+0.46}_{-0.47}$. This lanthanide fraction
is inconsistent with the $r$-process abundance pattern seen in the Solar
system and beyond (Ji et al. 2019).

We have not yet inferred the abundances at later epochs: 2.4 days, 3.4 days,
and beyond. At later epochs, as the kilonova ejecta expands and becomes more
optically thin, we expect that the photosphere recedes into the ejecta. This
may uncover additional components which power the kilonova at later times,
after being physically hidden underneath the photosphere at early times and/or
out-shined by the early blue emission. The emergence of new components at
later times would be consistent with results from light curve modeling (e.g.,
Villar et al. 2017), which has indicated the presence of multiple ejecta
components; in particular, a redder component emerging at $\sim$3 days post-
merger. These different components may originate from different ejection
mechanisms during the merger, and are characterized by different masses,
velocities, and heating rates as a function of time. In a NS-NS merger, we
expect a redder component mostly confined to the plane of the initial binary
from the tidal ejecta, with little neutrino reprocessing, a bluer squeezed
polar component from the collisional interface between the two NSs, and a more
isotropic red/blue disk wind from an accretion disk around a merger remnant.
Our inferred $Y_{e}$ and $s$ at 1.4 days are consistent with an outflow
produced over timescales longer than the dynamical time which has been
substantially reprocessed by neutrinos, e.g., winds from a hypermassive
neutron star remnant or an accretion disk. At later times, the spectrum might
be dominated or better described by a different ejecta or a multi-component
ejecta configuration.

Here, we explore the time evolution of the inferred abundance pattern and the
need for multi-component ejecta models to fit the spectra of AT2017gfo, at
1.4, 2.4, and 3.4 days post-merger. We first extend our single-component
fitting to 2.4 and 3.4 days to obtain the best single-component models for the
ejecta. We then develop a model in the radiative transfer code where the
ejecta is stratified and multi-component. We compare our best multi-component
fits at 1.4, 2.4, and 3.4 days to their single-component equivalents. The
abundances of these single- and multi-component models are then examined to
paint a picture of the inferred abundance pattern as a function of time.

This paper is organized as follows. In Section 2, we briefly review SPARK and
present the upgrades which enable fitting the later epochs of AT2017gfo with
both single- and multi-component models. In Section 3, we present our fits. In
Section 4, we explore the time evolution of the inferred abundance pattern,
the species present in the ejecta, and the physical origin of different
components. We briefly conclude in Section 5.

##  2 Methods

###  2.1 Spectroscopic $r$-Process Abundance Retrieval for Kilonovae (SPARK)

We briefly describe our tool, SPARK, and refer the reader to V23 for more
detail. SPARK (Spectroscopic $r$-Process Abundance Retrieval for Kilonovae) is
designed as a modular inference engine for extracting key kilonova parameters
from optical spectra, determining the element-by-element abundance pattern of
the ejecta, and associating absorption features in the spectra with particular
species.

In SPARK, we use the 1D TARDIS (Kerzendorf & Sim 2014; Kerzendorf et al. 2023)
radiative transfer code to generate a set of synthetic spectra. In TARDIS,
photon packets are propagated through shell(s) of plasma, where they may
undergo either bound-bound processes or electron scattering. To handle bound-
bound (matter-radiation) interactions in the ejecta, we require a list of
lines for the species in the ejecta. We use a line list of observed lines,
obtained through the Vienna Atomic Line Database (VALD; Ryabchikova et al.
2015; Pakhomov et al. 2019). Each spectrum is parameterized by a set
$\theta_{i}$ of parameters, including a luminosity, density, inner/outer
computational boundary velocities, and three key parameters which set the
abundances in the ejecta: the electron fraction $Y_{e}$, expansion velocity
$v_{\mathrm{exp}}$, and specific entropy per nucleon $s/k_{\mathrm{B}}$. These
latter three parameters describe different abundance patterns output by the
nuclear reaction network calculations of Wanajo (2018), which use parametric
outflow trajectories, allowing us to infer the abundance pattern of the
ejecta.

To perform this inference, we express our likelihood function using the full
formalism of Czekala et al. (2015) for likelihoods involving spectroscopic
data. However, because of the considerable computational cost of spectral
synthesis with TARDIS, we do not use more common methods such as Markov chain
Monte Carlo (MCMC) or nested sampling for inference—rather, we couple TARDIS
to the approximate posterior estimation scheme of approxposterior (Fleming &
VanderPlas 2018; Fleming et al. 2020). In this scheme, we introduce a Gaussian
Process (GP) surrogate for the posterior $L_{p}(\theta)$ and employ Bayesian
Active Posterior Estimation (BAPE; Kandasamy et al. 2017). BAPE is a form of
active learning in which we maximize an acquisition function with terms
including both the mean $\mu(\theta)$ and variance $\sigma^{2}(\theta)$ of the
GP. This acquisition function thus balances exploration (of the parameter
space) and exploitation (sampling around the peak(s) of the posterior). The GP
is iteratively retrained as new points $(\theta,L_{p}(\theta))$ are added to a
training set, and this GP converges to an approximation of the posterior.

In all, inference is dramatically accelerated, and we obtain (among the other
parameters) the $Y_{e}$, $v_{\mathrm{exp}}$, and $s/k_{\mathrm{B}}$ which best
describe the ejecta, with relatively few forward model evaluations. In V23, we
fit the VLT/X-shooter spectrum of AT2017fgfo (Pian et al. 2017; Smartt et al.
2017) at 1.4 days with a base set of 1500 Latin Hypercube samples + 1140 BAPE
active learning samples. This is a factor of $\sim 10^{3}$ fewer samples than
might be required with a standard MCMC for a similar 6-dimensional fit.

###  2.2 Multi-component, stratified ejecta with TARDIS

In V23, we modeled the kilonova ejecta as a single shell with a uniform
abundance pattern. The plasma in this shell is also described by a single
temperature and mass/electron density. This configuration is fully described
by the luminosity at the outer boundary $L_{\mathrm{outer}}$, (which in fact
sets the initial guess for the temperature at the inner boundary), the
normalization in the density power law $\rho_{0}$, the inner and outer
boundary velocities $v_{\mathrm{inner}}$ and $v_{\mathrm{outer}}$222TARDIS
assumes homologous expansion, in which $v\propto r$, such that velocity can be
understood as a spatial coordinate. Given a density profile, the velocity
range sets the mass of the ejecta., and three parameters which set the
abundance pattern: electron fraction $Y_{e}$, expansion velocity
$v_{\mathrm{exp}}$, and specific entropy $s/k_{\mathrm{B}}$. A single-
component fit is thus 7-dimensional, unless one or more of the parameters are
fixed. For example, we fixed $v_{\mathrm{outer}}=0.35c$ in our fit to the 1.4
day spectrum. This setup can describe a single ejecta component such as
dynamical ejecta or some outflow. It may also describe a kilonova in which one
component significantly dominates (by mass or by the strength of the
absorption/emission features) over the other(s).

Here, we implement multi-component ejecta. TARDIS allows for ejecta composed
of stratified radial shells, each with a specific temperature, density, plasma
conditions, and composition. In this configuration, each shell can have a
specific abundance pattern. In single-shell runs, we compute the plasma state
only once given $L_{\mathrm{outer}}$. For these multi-shell runs, we must
employ multiple TARDIS iterations when generating synthetic spectra. At each
iteration, the plasma conditions are updated, and converge towards an ejecta
where the luminosity emitted at the outer boundary matches the user-requested
$L_{\mathrm{outer}}$. We use 10 such shells and 30 such iterations in all
runs.

We begin with a simple two-component ejecta. As with our single-component
model, this two-component model has an outer boundary luminosity
$L_{\mathrm{outer}}$ and density $\rho_{0}$. Each of the components is then
described by inner and outer boundary velocity ($v_{\mathrm{inner}}$ and
$v_{\mathrm{outer}}$) and a $Y_{e}$, $v_{\mathrm{exp}}$, and
$s/k_{\mathrm{B}}$. Two-component fits are thus $2+5+5=12$-dimensional. The
two components necessarily overlap in physical space because TARDIS cannot
simulate a gap between them. The abundance in each shell is then determined by
the component(s) which are in a given shell. For shells where there is overlap
between two components, the abundance is taken as a sum of the two abundance
patterns and renormalized to unity.

Multiple components allow for additional complexity in the spectral
synthesis.333See Kawaguchi et al. (2020) for an exploration of the diversity
of kilonovae which may be produced when multiple components are present. In
particular, we can produce the effect of reprocessing, where emission from one
component is absorbed and re-emitted/scattered by another. We can also produce
the effect of lanthanide curtaining, in which some outer lower-$Y_{e}$ ejecta
containing the lanthanides masks an inner, bluer, lighter-element ejecta,
generating a redder kilonova due to the considerable opacity of the
lanthanides in the near-UV and optical. Some kilonova spectra may be better-
described by these multi-component ejecta models. In V23, we find that the 1.4
day of AT2017gfo is well-described by a single-component ejecta with
$Y_{e}=0.31$. However, at later epochs, as the ejecta expands and becomes more
optically thin, the photosphere recedes into the ejecta and we may unmask
additional components which were hidden or out-shined at early times. Light
curve modelling (e.g., Villar et al. 2017) has shown that the kilonova may
indeed be better-described by multiple components of different opacities, and
some spectral modelling (e.g., Kasen et al. 2017) also invokes multiple
components. This motivates our introduction of multi-component ejecta into
SPARK.

###  2.3 Inference setup

All approxposterior/BAPE hyperparameters and optimizers used in this work are
the same as those used in V23. We again produce a base set of $m_{0}=1500$
Latin Hypercube sampled points at the beginning of each SPARK run. However,
the parameter space allowed by our priors differs for the 1.4, 2.4, and 3.4
day fits. Table 1 includes our (all uniform) priors for each of our single-
component fits. The bounds on the density $\rho_{0}$ and abundance-setting
parameters $Y_{e}$, $v_{\mathrm{exp}}$, and $s$ are the same at all epochs.
The priors differ in the bounds on the luminosity, as expected given the
cooling of the ejecta over time as it expands. We further allow for wider
priors on the inner and outer boundary velocities for the fits at later
epochs. In V23, we fixed $v_{\mathrm{outer}}=0.35c$ during our 1.4 day fit but
noted that we observed similar results for $v_{\mathrm{outer}}$ in the range
$0.35-0.38c$. Here, we allow for greater flexibility in $v_{\mathrm{outer}}$
later epochs.

Our priors for our multi-component fits are given in Table 2. These are broad,
and identical, at 1.4 and 2.4 days. Aside from requiring
$v_{\mathrm{inner,i}}<v_{\mathrm{outer,i}}$, our priors also require some
overlap between the two components, i.e.,
$v_{\mathrm{outer},2}>v_{\mathrm{inner},1}$ or
$v_{\mathrm{outer},1}>v_{\mathrm{inner},2}$. To allow for even sampling of
parameter space for the $m_{0}=1500$ points in the base training set with
these conditional constraints, we use constrained Latin Hypercube Sampling
(Petelet et al. 2009). At 3.4 days, we encounter challenges with the
convergence of the posterior for these broad, conditional priors in
velocities. We thus use tighter priors on the velocities at this epoch, but
all other priors are the same as those at 1.4 and 2.4 days.

When attempting to fit the full 3.4 day spectrum, we find that the fit
converges to a blackbody with little to no absorption to fit the continuum of
the spectrum, especially at the shortest wavelengths $\lesssim$6400 Å. Since
we are interested in determining the species involved in the absorption, and
especially the prominent feature at $\sim$8000 Å, we prioritize fitting the
$\geqslant$6400 Å region of the spectrum by excluding shorter wavelengths in
the computation of the likelihood. We perform this exclusion for both single-
and multi-component fits.

Finally, when expressing the likelihood using the Czekala et al. (2015)
formalism, we use a global covariance term with amplitude
$a_{\mathrm{G}}=10^{-34}~{}(\mathrm{erg~{}s^{-1}~{}cm^{-2}}$ Å${}^{-1})^{2}$
and a correlation length scale $\ell=0.025c$ for all 1.4 and 2.4 day fits. For
3.4 day fits, after some trial and error, we find that a smaller
$a_{\mathrm{G}}=10^{-35}~{}(\mathrm{erg~{}s^{-1}~{}cm^{-2}}$ Å${}^{-1})^{2}$
better matches the uncertainties on the observed spectrum at this epoch.

Table 1: Uniform priors for the single-component fits at 1.4, 2.4, and 3.4 days. $v_{\mathrm{outer}}$ is fixed to $0.35c$ in the 1.4 day fit (see V23). parameter | 1.4 days | 2.4, 3.4 days
---|---|---
$\log_{10}(L_{\mathrm{outer}}/L_{\odot})$ | $[7.6,8.0]$ | $[7.2,7.8]$
$\log_{10}(\rho_{0}/\mathrm{g~{}cm^{-3}})$ | $[-16.0,-14.0]$ | same
$v_{\mathrm{inner}}/c$ | $[0.250,0.340]$ | $[0.100,0.275]$
$v_{\mathrm{outer}}/c$ | $0.35c$ | $[0.280,0.400]$
$v_{\mathrm{exp}}/c$ | $[0.05,0.30]$ | same
$Y_{e}$ | $[0.01,0.40]$ | same
$s~{}[k_{\mathrm{B}}/\mathrm{nucleon}]$ | $[10,35]$ | same
Table 2: Uniform priors for the multi-component fits at 1.4, 2.4, and 3.4 days. The priors on the velocities are tighter at 3.4 days. parameter | 1.4, 2.4 days | 3.4 days
---|---|---
$\log_{10}(L_{\mathrm{outer}}/L_{\odot})$ | $[7.0,8.0]$ | same
$\log_{10}(\rho_{0}/\mathrm{g~{}cm^{-3}})$ | $[-16.0,-14.0]$ | same
$v_{\mathrm{inner,1}}/c$ | $[0.10,0.35]$ | $[0.10,0.29]$
$v_{\mathrm{outer,1}}/c$ | $[0.25,0.40]$ | $[0.30,0.40]$
$v_{\mathrm{exp,1}}/c$ | $[0.05,0.30]$ | same
$Y_{e,1}$ | $[0.01,0.50]$ | same
$s_{1}~{}[k_{\mathrm{B}}/\mathrm{nucleon}]$ | $[10,35]$ | same
$v_{\mathrm{inner,2}}/c$ | $[0.10,0.35]$ | $[0.20,0.30]$
$v_{\mathrm{outer,2}}/c$ | $[0.25,0.40]$ | $[0.31,0.40]$
$v_{\mathrm{exp,2}}/c$ | $[0.05,0.30]$ | same
$Y_{e,2}$ | $[0.01,0.50]$ | same
$s_{2}~{}[k_{\mathrm{B}}/\mathrm{nucleon}]$ | $[10,35]$ | same

##  3 Fitting Later Epochs and Multi-Component Ejecta

Figure 1: Compilation of all best fit models to the spectrum of the GW170817 kilonova from this work: 1.4, 2.4, and 3.4 days post-merger, for both single- and multi-component models. The preferred models are highlighted; disfavored are more translucent. The single-component models are favored at 1.4 and 2.4 days (solid blue and purple lines). At 3.4 days, we require a multi-component ejecta to adequately fit the spectrum (dashed red line). The abundance patterns corresponding to these preferred models are included in Figure 2. Zoomed in versions of these best fits are included in Figures 6 and 7 (Appendix A). Best fits are obtained as the median of the posterior, or the median of some mode of the posterior; these full posteriors are also included in Appendix A. Table 3: Best fit parameters for single-component fits to the GW170817 kilonova at $1.4$, $2.4$, and $3.4$ days. The $1.4$ day fit is the “purple + warm” model of V23. $v_{\mathrm{outer}}$ is fixed to $0.35c$ in the 1.4 day fit. The inner boundary temperature $T_{\mathrm{inner}}$ and lanthanide mass fraction $X_{\mathrm{lan}}$ are derived parameters, i.e., they are not dimensions in $\theta$-space. parameter | 1.4 days | 2.4 days | 3.4 days
---|---|---|---
$\log_{10}(\frac{L_{\mathrm{outer}}}{L_{\odot}})$ | $7.782^{+0.013}_{-0.014}$ | $7.594^{+0.040}_{-0.061}$ | $7.531^{+0.017}_{-0.016}$
$\log_{10}(\frac{\rho_{0}}{\mathrm{g~{}cm^{-3}}})$ | $-15.016^{+0.320}_{-0.316}$ | $-15.443^{+0.742}_{-0.463}$ | $-14.586^{+0.313}_{-0.384}$
$v_{\mathrm{inner}}/c$ | $0.313^{+0.013}_{-0.014}$ | $0.249^{+0.017}_{-0.032}$ | $0.253^{+0.011}_{-0.021}$
$v_{\mathrm{outer}}/c$ | $0.35$ | $0.342^{+0.047}_{-0.050}$ | $0.309^{+0.032}_{-0.023}$
$v_{\mathrm{exp}}/c$ | $0.240^{+0.055}_{-0.082}$ | $0.172^{+0.107}_{-0.101}$ | $0.132^{+0.094}_{-0.056}$
$Y_{e}$ | $0.311^{+0.013}_{-0.011}$ | $0.306^{+0.055}_{-0.204}$ | $0.226^{+0.062}_{-0.067}$
$s~{}[k_{\mathrm{B}}/\mathrm{nuc}]$ | $13.6^{+4.1}_{-3.0}$ | $17.6^{+7.1}_{-6.3}$ | $15.4^{+7.1}_{-6.3}$
$T_{\mathrm{inner}}~{}[\mathrm{K}]$ | $3958^{+87}_{-94}$ | $3050^{+126}_{-223}$ | $2455^{+59}_{-104}$
$\log_{10}X_{\mathrm{lan}}$ | ${-7.03}^{+0.46}_{-0.47}$ | $-5.64^{+1.05}_{-1.31}$ | $-2.10^{+1.28}_{-6.20}$
Table 4: Best fit parameters for multi-component fits to the GW170817 kilonova at $1.4$, $2.4$, and $3.4$ days. The ejecta mass above the photosphere $M_{\mathrm{phot},1;2}$ (not the total ejecta mass), lanthanide mass fractions $X_{\mathrm{lan,1;2}}$, and the total $X_{\mathrm{lan,total}}$ are derived parameters. parameter | 1.4 days | 2.4 days | 3.4 days
---|---|---|---
$\log_{10}(\frac{L_{\mathrm{outer}}}{L_{\odot}})$ | $7.854^{+0.012}_{-0.017}$ | $7.700^{+0.065}_{-0.073}$ | $7.605^{+0.049}_{-0.040}$
$\log_{10}(\frac{\rho_{0}}{\mathrm{g~{}cm^{-3}}})$ | $-15.095^{+0.127}_{-0.401}$ | $-15.440^{+0.813}_{-0.467}$ | $-14.505^{+0.323}_{-0.372}$
$v_{\mathrm{inner,1}}/c$ | $0.323^{+0.007}_{-0.020}$ | $0.260^{+0.061}_{-0.039}$ | $0.213^{+0.056}_{-0.035}$
$v_{\mathrm{outer,1}}/c$ | $0.326^{+0.009}_{-0.012}$ | $0.335^{+0.054}_{-0.065}$ | $0.344^{+0.035}_{-0.038}$
$v_{\mathrm{exp,1}}/c$ | $0.118^{+0.029}_{-0.025}$ | $0.170^{+0.108}_{-0.100}$ | $0.198^{+0.092}_{-0.102}$
$Y_{e,1}$ | $0.139^{+0.028}_{-0.114}$ | $0.288^{+0.129}_{-0.187}$ | $0.228^{+0.073}_{-0.088}$
$s_{1}~{}[k_{\mathrm{B}}/\mathrm{nuc}]$ | $24.6^{+3.4}_{-5.1}$ | $21.1^{+10.5}_{-9.5}$ | $14.9^{+8.0}_{-4.5}$
$v_{\mathrm{inner,2}}/c$ | $0.281^{+0.011}_{-0.013}$ | $0.268^{+0.049}_{-0.045}$ | $0.232^{+0.038}_{-0.027}$
$v_{\mathrm{outer,2}}/c$ | $0.355^{+0.015}_{-0.015}$ | $0.335^{+0.045}_{-0.070}$ | $0.334^{+0.037}_{-0.022}$
$v_{\mathrm{exp,2}}/c$ | $0.195^{+0.034}_{-0.036}$ | $0.183^{+0.089}_{-0.110}$ | $0.134^{+0.091}_{-0.069}$
$Y_{e,2}$ | $0.340^{+0.022}_{-0.021}$ | $0.261^{+0.163}_{-0.124}$ | $0.161^{+0.149}_{-0.104}$
$s_{2}~{}[k_{\mathrm{B}}/\mathrm{nuc}]$ | $20.1^{+3.8}_{-3.6}$ | $22.0^{+10.0}_{-10.1}$ | $21.5^{+8.3}_{-9.5}$
$M_{\mathrm{phot},1}~{}[M_{\odot}]$ | $6.9^{+22.1}_{-6.0}\times 10^{-7}$ | $2.3^{+12.5}_{-0.7}\times 10^{-5}$ | $5.3^{+1.8}_{-1.3}~{}\times 10^{-4}$
$\log_{10}X_{\mathrm{lan,1}}$ | $-1.08^{+0.34}_{-0.34}$ | $-0.63^{+0.26}_{-11.1}$ | $-0.97^{+0.23}_{-0.32}$
$M_{\mathrm{phot},2}~{}[M_{\odot}]$ | $7.0^{+0.8}_{-1.0}\times 10^{-5}$ | $4.0^{+4.8}_{-0.9}\times 10^{-5}$ | $3.8^{+1.1}_{-0.9}~{}\times 10^{-4}$
$\log_{10}X_{\mathrm{lan,2}}$ | ${-12.99}^{+3.15}_{-3.45}$ | $-4.53^{+2.61}_{-3.53}$ | $-0.88^{+0.32}_{-0.55}$
$\log_{10}X_{\mathrm{lan,total}}$ | $-12.99^{+8.82}_{-3.45}$ | $-3.39^{+1.34}_{-1.75}$ | $-0.97^{+0.23}_{-0.32}$
Figure 2: Best fit abundance patterns at 1.4, 2.4, and 3.4 days. The
abundances at 1.4 and 2.4 days are taken from the preferred single-component
fits, while the abundance at 3.4 days is taken from the preferred multi-
component fit (see Figure 1). The overall abundance of the multi-component
ejecta at 3.4 days is the mass-weighted sum over all 10 shells in the
stratified ejecta. A new, redder kilonova component emerges at 3.4 days post-
merger. Uncertainty bands are obtained by taking additional samples from the
the posterior, effectively propagating the uncertainty on $Y_{e}$,
$v_{\mathrm{exp}}$, and $s$ into the abundances. We also show the Solar
$r$-process pattern, computed using data from Lodders et al. (2009) subtracted
by the $s$-process residuals of Bisterzo et al. (2014). The best fit
abundances are evidently non-Solar at the first two epochs, but are closer to
Solar at the third.

1.4 days

2.4 days

3.4 days

Figure 3: Inferred physical ejecta configuration, and the composition-setting
parameters, for the multi-component fits. Top to bottom: 1.4, 2.4, and 3.4
days. At 1.4 days, a thin red component containing $\sim 1$% of the total mass
is inferred. This component has a negligible impact on the emergent spectrum,
indicating that this model is effectively single-component. At 2.4 days, two
components with highly similar compositions completely overlap in space, also
effectively indicating a single-component ejecta. At 3.4 days, two components
with different compositions are present; the ejecta is multi-component. The
1.4 and 2.4 day fits are disfavored compared to the single-component
equivalents, while the multi-component fit shown here is favored for 3.4 days.

We model the 1.4, 2.4, and 3.4 day spectra of AT2017gfo with both single- and
multi-component (2-component) models, and compare the results at each epoch.
In Figure 1, we compile our best fit spectra, and highlight our favored
models. The parameters for our best fit single- and multi-component models are
included in Tables 3 and 4, respectively. We find that a single-component blue
ejecta model best describes the 1.4 and 2.4 day spectra. At 3.4 days, we
require an additional red, lanthanide-rich component in a multi-component
ejecta configuration. In the following sections, we present and examine each
of the best fits in turn, at each epoch.

Appendix A contains compilations of all $m_{0}+m_{\mathrm{active}}$ spectra
generated in a given SPARK run, for the 1.4, 2.4, and 3.4 day single- and
multi-component runs. Posteriors generated by the single- and multi-component
runs are also included. Unless stated otherwise, we use the median and 2.5% /
97.5% quantiles of the posterior, or the median and these quantiles of some
mode of a multi-modal posterior, to determine our best fit and uncertainties.
All posteriors converge with sufficient $m_{\mathrm{active}}$ active learning
samples, according to the $z$ convergence diagnostic given in Fleming et al.
(2020).

###  3.1 Single-and multi-component modelling at 1.4 days

In V23, we fit the 1.4 day spectrum of AT2017gfo with a single-component
ejecta model. This fit required $m_{0}+m_{\mathrm{active}}=1500+1140$ points
to obtain a good fit and converged posterior (see Figure 6). We describe the
best fit “purple + warm” model here, and refer the reader to V23 for details.
The purple + warm model is characterized by a photospheric velocity of
$v_{\mathrm{inner}}/c=0.313^{+0.013}_{-0.014}$; while the outer boundary
velocity was fixed to $v_{\mathrm{outer}}=0.35c$. For abundance-setting
parameters, we found an electron fraction $Y_{e}=0.311^{+0.013}_{-0.011}$,
expansion velocity of $v_{\mathrm{exp}}/c=0.240^{+0.055}_{-0.082}$, and
entropy of $s/k_{\mathrm{B}}=13.6^{+4.1}_{-3.0}$. In V23, we referred to this
model as “purple + warm”, but emphasize that this is nonetheless substantially
blue ejecta within the broader context of kilonova components. These
parameters correspond to an abundance pattern that is poor in lanthanides, but
rich in strontium (38Sr), allowing for a good fit to the prominent $\sim$8000
Å feature. The fit struggled at shorter wavelengths: over $3500-4500$ Å, we
overestimated the absorption from yttrium (39Y) and zirconium (40Zr), while at
$\lesssim$3500 Å, we underestimated. Nonetheless, the single-component model
remains our preferred model at 1.4 days, as we elaborate on below.

In our 1.4 day multi-component SPARK run, we require
$m_{0}+m_{\mathrm{active}}=1500+1440$ points to obtain a good fit and
converged posterior (see Figure 7). More active learning samples are required
than for single-component due to the increased dimensionality of the multi-
component fits. The resultant fit generally captures the shape of the
spectrum, and the broad absorption at $\sim$8000 Å, like the single-component
equivalent. However, the fit does not capture the continuum at wavelengths
$\gtrsim$10,500 Å nor $\sim$ 3500 - 5000 Å. The fit is marginally closer to
the observed spectrum at the shortest wavelengths ($\lesssim$ 3500 Å) but this
region lies at the edge of the X-shooter spectrograph sensitivity, and overall
the single-component model is favored by visual inspection.

We infer an outer boundary luminosity of
$\log_{10}(L_{\mathrm{outer}}/L_{\odot})=7.854^{+0.012}_{-0.017}$. This is
brighter than is inferred in the single-component equivalent. This is as
expected: due to the use of 30 TARDIS iterations (rather than 1),
$L_{\mathrm{outer}}$ is updated at each iteration. Hence, in our multi-
component models, the meaning of $L_{\mathrm{outer}}$ is different and cannot
be mapped directly to $T_{\mathrm{inner}}$. ($L_{\mathrm{outer}}$ is also
higher in the 2.4 and 3.4 day multi-component models, compared to the single-
component equivalents).

This multi-component run yields two components with substantially different
abundance patterns. The first component has a lower electron fraction
$Y_{e,1}=0.139^{+0.028}_{-0.114}$ (and specific entropy
$s_{1}/k_{\mathrm{B}}=24.6^{+3.4}_{-5.1}$), yielding. The second has higher
$Y_{e,2}=0.340^{+0.022}_{-0.021}$ (and
$s_{2}/k_{\mathrm{B}}=20.1^{+3.8}_{-3.6}$), distinctly bluer. Indeed, the
$Y_{e,2}$ of this second, bluer component is consistent with the single-
component fit at 1.4 days. As in the single-component fit, expansion
velocities ($v_{\mathrm{exp},1}$ and $v_{\mathrm{exp},2}$) are poorly
constrained. Interestingly, the better-constrained entropies ($s_{1}$ and
$s_{2}$) are both higher than and inconsistent with the single-component
equivalent.

However, despite the substantial differences between the two components, one
component dominates over the other in producing the emergent spectrum. A
physical picture for the components is included in Figure 3. We see that the
red component is confined to only 2 shells, while the blue component extends
over a greater radius. The mass contained in these two components is
$M_{\mathrm{phot},1}/M_{\odot}=6.9^{+22.1}_{-6.0}\times 10^{-7}$ and
$M_{\mathrm{phot},2}/M_{\odot}=7.0^{+0.8}_{-1.0}\times 10^{-5}$, respectively.
We emphasize that this is the mass above the photosphere, and not the total
ejecta mass. Nonetheless, the bluer component contains $\sim$100$\times$ as
much mass as the redder, and dominates the absorption/emission in the
spectrum. In this sense, despite the different compositions of the two
components (Figure 13, Appendix B), this multi-component SPARK run has
effectively converged to a single-component model.

###  3.2 Single-and multi-component modelling at 2.4 days

In our 2.4 day single-component SPARK run, We require
$m_{0}+m_{\mathrm{active}}=1500+600$ points to obtain a good fit and converged
posterior (see Figure 6). The resultant posterior shows some bimodality, most
evident in the $s/k_{\mathrm{B}}$ dimension. Discarding samples from the
posterior with $s/k_{\mathrm{B}}\geqslant 25.0$, we obtain a superior fit, and
use this as our preferred model. This fit captures most of the continuum of
the observed spectrum, and most of the $\sim$8000 Å absorption. If this
absorption belongs to a P Cygni feature (as has been argued in Watson et al.
2019; Sneppen et al. 2023), we partially miss the red wing of this feature.
Furthermore, we overestimate the absorption, or underestimate the continuum,
at wavelengths $\lesssim$4500 Å.

We infer an $\log_{10}(L_{\mathrm{outer}}/L_{\odot})=7.594^{+0.040}_{-0.061}$.
In the case of single-shell, single-iteration TARDIS runs, we can use these
values to determine an inner boundary temperature of
$T_{\mathrm{inner}}=3050^{+126}_{-223}~{}\mathrm{K}$. As expected, the
kilonova has dimmed and cooled since the previous epoch. The inferred
velocities also match expectations: we obtain inner and outer boundary
velocities $v_{\mathrm{inner}}/c=0.249^{+0.017}_{-0.032}$ and
$v_{\mathrm{outer}}=0.342^{+0.047}_{-0.050}$. This $v_{\mathrm{outer}}$ is
consistent with the fixed $v_{\mathrm{outer}}=0.35c$ from 1.4 days, while
$v_{\mathrm{inner}}$ (effectively the photospheric velocity) has receded into
the ejecta, compared to $0.31c$ at 1.4 days. This is evidence for the ejecta
expanding, cooling, and become optically thinner, as expected.

Finally, we obtain $Y_{e}=0.306^{+0.055}_{-0.204}$ and
$s/k_{\mathrm{B}}=17.6^{+7.1}_{-6.3}$. Similar to our results at 1.4 days,
$v_{\mathrm{exp}}$ is poorly constrained, while $s$ is better-constrained.
This $s$ is also consistent with the previous epoch. The posterior
distribution in all dimensions is wider at 2.4 days; in particular, $Y_{e}$
exhibits a tail extending to smaller electron fractions. The posterior is
nonetheless peaked at $Y_{e}=0.306^{+0.055}_{-0.204}$, which is slightly lower
than at 1.4 days, but consistent to within the uncertainties.

In our 2.4 day, multi-component SPARK run, we require
$m_{0}+m_{\mathrm{active}}=1500+1020$ points to obtain a good fit and
converged posterior (see Figure 7). The resultant fit once again captures the
broad shape and $\sim$8000 Å absorption, but the depth of the $\sim$8000 Å
absorption feature is overestimated. The multi-component fit does, however,
achieve a better fit to the continuum at wavelengths $\lesssim$7000 Å.

In contrast with the 1.4 day, multi-component run, we see that the two
inferred components are remarkably similar. The first component is described
by $Y_{e,1}=0.288^{+0.129}_{-0.187}$ and
$s_{1}/k_{\mathrm{B}}=21.1^{+10.5}_{-9.5}$, while the latter has
$Y_{e,2}=0.261^{+0.163}_{-0.124}$ and
$s_{2}/k_{\mathrm{B}}=22.0^{+10.0}_{-10.1}$. Again, expansion velocities are
poorly constrained, while the better-constrained entropies are higher than and
inconsistent with the single-component equivalent. Both entropies are
consistent with the entropies inferred in the 1.4 day, multi-component model.
Finally, $Y_{e,1}$ and $Y_{e,2}$ are both lower than, but still consistent
with, the dominant blue component in the 1.4 day, multi-component model and
the 1.4 day, single-component (purple+warm) model.

Our inferred geometry for the two components is shown in Figure 3. As
expected, the photosphere recedes into the ejecta. The two components almost
completely overlap in physical space; indeed, the mass contained in these two
components is roughly equally split into
$M_{\mathrm{phot},1}/M_{\odot}=2.3^{+12.5}_{-0.7}\times 10^{-5}$ and
$M_{\mathrm{phot},2}/M_{\odot}=4.0^{+4.8}_{-0.9}\times 10^{-5}$, respectively.
Furthermore, the similarity between $(Y_{e,1},~{}v_{\mathrm{exp},1},~{}s_{1})$
and $(Y_{e,2},~{}v_{\mathrm{exp},2},~{}s_{2})$ yields similar compositions for
both components (Figure 13, Appendix B). Given the two components’ similar
compositions and equal contributions to the emergent spectrum, this multi-
SPARK run has also converged to an effectively single-component model.

###  3.3 Single-and multi-component modelling at 3.4 days

Our best single-component model for 3.4 days is in general a poor fit to
AT2017gfo. Nonetheless, we review it here for completeness. As at 2.4 days, we
require $m_{0}+m_{\mathrm{active}}=1500+600$ points (see Figure 6) to obtain a
converged posterior. The resultant posterior shows a high degree of multi-
modality, most evident in the $\rho_{0}$, $Y_{e}$, and $s/k_{\mathrm{B}}$
dimensions. Of the multiple modes, we obtain our best (but still poor) fit
with a higher-$\rho_{0}$, mid-$Y_{e}$, lower-$s$ mode. This fit produces some,
but not all, of the dominant absorption feature at $\sim$8000 Å. This fit also
overestimates the absorption, or underestimates the continuum, for wavelengths
$\lesssim$5500 Å.

We find $\log_{10}(\rho_{0}/\mathrm{g~{}cm^{-3}})=-14.586^{+0.313}_{-0.384}$.
This density is larger than that inferred at earlier epochs. We also infer
$v_{\mathrm{outer}}=0.309^{+0.032}_{-0.023}$. This $v_{\mathrm{outer}}$ is
much smaller than the fixed value at 1.4 days or the inferred at 2.4 days.
This may explain the insufficient absorption at $\sim$8000 Å, if the photon
packets are not interacting with enough matter to be absorbed in large numbers
before escaping the ejecta. Finally, we find $Y_{e}=0.226^{+0.062}_{-0.067}$
and $s/k_{\mathrm{B}}=15.4^{+7.2}_{-3.1}$ at this epoch. $Y_{e}$ is lower and
inconsistent with both previous epochs. Interestingly, despite the lower
$Y_{e}$, the inferred abundance of Sr is within $1$% of that of our single-
component model at 1.4 days, and in fact $\sim$25% larger than that of our
single-component model at 2.4 days. Nonetheless, Sr does not yield the
expected prominent absorption feature at $\sim$8000 Å. This suggests that the
shortcomings of this model may not be due to the inferred abundance pattern,
but rather the density, velocity structure, and/or some inherent shortcoming
of the single-component model.

For our preferred, multi-component model, we require
$m_{0}+m_{\mathrm{active}}=1500+1140$ points to obtain a good fit and
converged posterior at 3.4 days (see Figure 7). The resultant fit is better
than the single-component equivalent. This fit better captures the $\sim$8000
Å absorption, and the red wing of the nominal P Cygni feature. However, this
fit also overestimates the absorption, or underestimates the continuum, for
wavelengths $\lesssim$5500 Å. This over/underestimation is more severe than in
the single-component model.

Interestingly, we find
$\log_{10}(\rho_{0}/\mathrm{g~{}cm^{-3}})=-14.505^{+0.323}_{-0.372}$, which is
higher than, and inconsistent with, the density inferred at 1.4 and 2.4 days.
This inconsistency (also present in the disfavored single-component model at
3.4 days) may arise from the emergence of a new component, with a different
density profile (either in functional form or in normalization), at this
epoch.

Most significantly, we infer two components with substantially different
properties. The first component is described by
$Y_{e,1}=0.228^{+0.073}_{-0.088}$ and
$s_{1}/k_{\mathrm{B}}=14.9^{+8.0}_{-4.5}$, while the latter has
$Y_{e,2}=0.161^{+0.149}_{-0.104}$ and
$s_{2}/k_{\mathrm{B}}=21.5^{+8.3}_{-9.5}$. Again, expansion velocities are
poorly constrained. Both entropies are better-constrained and consistent with
the entropies inferred in the single-component (and multi-component) models at
1.4 and 2.4 days. The $Y_{e,1}$ and $Y_{e,2}$ are both substantially lower
than, and inconsistent with, those of the single- and multi-component fits at
1.4 and 2.4 days. This leads to an ejecta that is substantially more
lanthanide-rich: The two components have
$\log_{10}X_{\mathrm{lan},1}=-0.97^{+0.23}_{-0.32}$ and
$\log_{10}X_{\mathrm{lan},2}=-0.88^{+0.32}_{-0.55}$, respectively. This yields
a total lanthanide fraction of
$\log_{10}X_{\mathrm{lan,total}}=-0.97^{+0.23}_{-0.32}$ (the same as
$\log_{10}X_{\mathrm{lan},1}$ after rounding). Figure 13 (Appendix B) shows
the abundances of the two components, with uncertainties. Both components
contain some substantial abundance of lanthanides, but the higher-$Y_{e}$ of
the two has $\sim$10$\times$ more of the crucial element Sr.

In Figure 3, we show a physical picture of the inferred ejecta. The
photosphere further recedes compared to 1.4 and 2.4 days. Both components
mostly overlap in space and the mass is roughly equally split: they have
masses $M_{\mathrm{phot},1}/M_{\odot}=5.3^{+1.8}_{-1.3}\times 10^{-4}$ and
$M_{\mathrm{phot},2}/M_{\odot}=3.8^{+1.1}_{-0.9}\times 10^{-4}$, respectively.
Given the roughly equal contributions from a component rich in lanthanides and
another with a similar abundance of lanthanides but $\sim$10$\times$ as much
Sr, and the fact that the multi-component model is a substantially better fit
than the single-component at 3.4 days, we interpret the ejecta as being
genuinely two-component at this epoch. As we will see, the equal contribution
from two components with distinct properties shapes the emergent spectrum.

###  3.4 Favored models

Figure 4: Hypothetical models at 1.4, 2.4, and 3.4 days, where the
higher-$Y_{e}$ component of each multi-component fit has been replaced by that
of the single-component, 1.4 days (purple + warm) model of V23. Here, we take
the best fit parameters from our multi-component runs (Table 4) and, after
fitting, swap out the higher-$Y_{e}$ component with the purple + warm
$(Y_{e}=0.311,~{}v_{\mathrm{exp}}/c=0.240,~{}s/k_{\mathrm{B}}=13.6)$.
$L_{\mathrm{outer}}$, $\rho_{0}$, and the velocities are left unchanged, and
thus the mass in each component is also unchanged. Comparing to the original
multi-component models (dashed lines, Figure 1), the differences are
negligible at all epochs.

Thus far, we have selected our favored models by visual inspection, seeing
which produces a better fit to the observed spectrum at a given epoch. Here,
we explore the physical consistency of our favored models.

We have found that the 1.4 and 2.4 day spectra are well-fit by a single bluer
component. Even when allowed to yield multiple components, the multi-component
SPARK runs yield effectively single-component best fits at 1.4 and 2.4 days.
At 1.4 days, the bluer component contains $\sim$100$\times$ the mass of the
redder, and the presence of the red component is thus inconsequential. At 2.4
days, the two components are roughly equal in mass and characterized by
similar abundance patterns; this fit is thus also effectively single-
component.

In contrast, at 3.4 days, the two components in the multi-component fit have
different abundance patterns. In particular, an additional lower-$Y_{e}$
component rich in lanthanides is required to fit the observed spectrum. The
lower-$Y_{e}$ component contains as much as $\sim$10$\times$ more lanthanides
than the higher-$Y_{e}$ component, given the broadness of the posterior. The
remaining higher-$Y_{e}$ component, however, contains $\sim$10$\times$ as much
Sr. Given our atomic data, line list, and (most importantly) abundance
patterns from the reaction network calculations of Wanajo (2018), no single
component is able to simultaneously yield the abundance of Sr and lanthanides
needed to reproduce the observed 3.4 day spectrum of AT2017gfo. As we will see
in the following sections, the role of the higher-$Y_{e}$ component is to
provide enough Sr to yield the $\sim$8000 Å absorption, while the
lower-$Y_{e}$ yields most of the lanthanides needed to produce the short-
wavelength $\lesssim$7500 Å absorption.

Aside from asking whether a single- or multi-component ejecta is individually
preferred at each epoch, we also ask: does the addition of a new component to
the original, purple + warm ejecta at 1.4 days improve the fit? To test this,
we replace the higher-$Y_{e}$ components of the multi-component models at 1.4,
2.4, and 3.4 days with this purple + warm component. (These higher-$Y_{e}$
components have $Y_{e}=$ $0.340^{+0.022}_{-0.021}$, $0.288^{+0.129}_{-0.187}$,
and $0.228^{+0.073}_{-0.088}$ at 1.4, 2.4, and 3.4 days, respectively). Figure
4 shows the spectra which result from replacing these components with the
parameters
$(Y_{e}=0.311,~{}v_{\mathrm{exp}}/c=0.240,~{}s/k_{\mathrm{B}}=13.6)$ at each
epoch. At all epochs, there is minimal impact on the resultant spectrum and we
retain a good fit. This is not surprising at 1.4 and 2.4 days, where the
inferred $Y_{e}$ is consistent with the purple + warm model. It is, however,
surprising at 3.4 days. Crucially, replacing the higher-$Y_{e}$ component at
3.4 days changes the overall abundance pattern, but the overall abundance of
Sr remains unchanged due to the presence of substantial Sr in the purple +
warm ejecta. This further suggests that the role of the higher-$Y_{e}$
component at 3.4 days is primarily to provide enough Sr to produce the
$\sim$8000 Å absorption, with the lower-$Y_{e}$ providing sufficient
lanthanides.

Overall, we infer that a new, redder, lanthanide-rich component emerges at 3.4
days. This new red ejecta component has photospheric velocity $0.21c$ and
extends out to $0.35c$, suggesting that this component was in fact present at
earlier epochs, but was partially buried underneath the photosphere and/or
out-shined by the brighter blue component.

##  4 Discussion

###  4.1 Abundances of the favored models

Given our favored inferred models, and their best fit ($Y_{e}$,
$v_{\mathrm{exp}}$, $s$), we study the inferred abundance pattern of the
ejecta as a function of time. We favor the single-component models at 1.4 and
2.4 days, and multi-component model at 3.4 days. We show the abundance
patterns of these favored models in Figure 2. The overall multi-component
abundance pattern at 3.4 days is computed as the mass-weighted sum over all
shells, including both components.

The 1.4 and 2.4 day abundance patterns are remarkably similar, but the 2.4 day
model has greater uncertainties due to the broader overall posterior. These
2.4 day abundances show evidence for $\sim 10^{-9}-10^{-6}$ mass fractions for
some of the lighter lanthanides. However, as we see in the following section,
lanthanides do not leave any clear imprint on the spectrum at 2.4 days.
Overall, the 1.4 and 2.4 day spectra are well-described by a single component
dominated by light $r$-process elements around the first $r$-process peak.

In contrast, 3.4 day abundance pattern indicates the clear presence of heavier
elements and especially lanthanides in the ejecta. The total lanthanide
fraction is $\log_{10}X_{\mathrm{lan,total}}=-0.97^{+0.23}_{-0.32}$, several
orders of magnitude larger than inferred at 1.4 and 2.4 days. The
higher-$Y_{e}$ component of this multi-component ejecta then provides Sr, the
role of which is crucial, as we will see in the the following section. Indeed,
we infer abundances of $Y_{\mathrm{Sr},1.4}=0.069^{+0.021}_{-0.010}$ and
$Y_{\mathrm{Sr},2.4}=0.087^{+0.060}_{-0.020}$ at 1.4 and 2.4 days, and
$Y_{\mathrm{Sr},3.4}=0.067^{+0.035}_{-0.002}$ (mass-weighted sum over all
shells, including both components) at 3.4 days. These are remarkably
consistent with each other.

The inferred abundance patterns at 1.4 and 2.4 days are markedly non-Solar in
Figure 2, where the Solar $r$-process abundance pattern is taken using the
data of Lodders et al. (2009) with the $s$-process residuals subtracted
according to Bisterzo et al. (2014). If all NS-NS/BH mergers produced such
blue ejecta, this would point to the inability of these mergers to yield the
$r$-process abundances seen in several astrophysical settings. The 3.4 day
abundance pattern, in contrast, is much closer to Solar up to $Z\sim 80$. In
particular, the presence of lanthanides improves this agreement.

Our abundance results are in broad agreement with inference on the light
curves and spectral fitting of AT2017gfo (see Ji et al. 2019, and references
therein). While light curve modelling typically uses grey (wavelength-
independent) opacities, this modelling generally recovers the emergence of
multiple components, including a higher-opacity redder component at $\gtrsim$3
days, as we find here. Our results for lanthanide fractions deviate somewhat
from previous works. In particular, previous works find lanthanide fractions
$\log_{10}X_{\mathrm{lan}}\sim-6$ to $\sim-4$ for a blue component, and
$\sim-2$ for a red component. We find a more lanthanide-poor bluer component,
and a marginally more lanthanide-rich ($\log_{10}X_{\mathrm{lan}}\sim-1$)
redder component. This is notable, as Ji et al. (2019) find that the
$\log_{10}X_{\mathrm{lan}}\sim-2$ inferred in previous works falls on the low
end of the lanthanide fractions measured in metal-poor stars. If our higher
inferred lanthanide fraction is correct, we relieve some of the tension
between the lanthanide fractions of the AT2017gfo ejecta and metal-poor stars.
However, we caution that we do not probe all of the ejecta, and in particular
are not sensitive to the ejecta below the photosphere. This ejecta could be
more or less lanthanide-rich, complicating these comparisons.

###  4.2 Elements & ions present in the favored models

Figure 5: Leave-one-out spectra for the favored models: single-component for
1.4 (top) and 2.4 (center) days, and multi-component for 3.4 days (bottom).
All models show clear absorption from strontium (38Sr) at $\sim$8000 Å. At 3.4
days, when the ejecta is richer in heavier elements, we also see absorption
(at $\sim$7000 Å and $\sim$12,000 Å) and emission (at $\sim$8000 Å) from the
lanthanide cerium (58Ce). However, we also see over-absorption from barium
(56Ba) at $\sim$4500Å. Spectral DEComposition (SDEC) plots in Figure 14
(Appendix C), provide a complementary view of the dominant species.

To assess the impact of different elements on the best fits, we generate
leave-one-out spectra. In each, we iteratively leave out a single element by
setting its abundance to 0 and transferring that original abundance to a
filler element. We use helium (He) as our filler element, since it should not
have a marked impact on the emergent spectrum when the ejecta remains
optically thick and local thermodynamic equilibrium (LTE) is a valid
approximation (Perego et al. 2022; Tarumi et al. 2023).

In Figure 5, we show leave-one-out spectra for the new favored models at 2.4
and 3.4 days, and the 1.4 day model from V23. We see the clear imprint of Sr
at $\sim$8000 Å at all epochs. This is consistent with Watson et al. (2019);
Gillanders et al. (2022); Domoto et al. (2021, 2022); Sneppen et al. (2023);
Sneppen & Watson (2023), which all argue for the importance of Sr in the
spectra. In our new model for 2.4 days, we also see some light evidence for
absorption from an adjacent first $r$-process peak element, yttrium (39Y), at
short wavelengths $\lesssim$4500 Å. Considering our previous identification of
Y at 1.4 days in V23, this identification at 2.4 days strengthens our previous
claim of the importance of Y, and is in agreement with Gillanders et al.
(2022). Interestingly, Y does not have the same pronounced impact at 3.4 days.
This is in agreement with Sneppen & Watson (2023), which finds that absorption
from Y is less prominent at 3.4 days, before a P Cygni feature emerges at
$\sim$7600 Å at 4.4 and 5.4 days. In our favored multi-component model, the
absorption from the lanthanides is much stronger at these (and shorter)
wavelengths.

At 3.4 days, we identify two new, heavier elements: barium (56Ba), an open
$s$-shell element in the same periodic table group as Sr, and cerium (58Ce), a
lanthanide. Ba is actually responsible for some of the overestimated
absorption at $\lesssim$5000 Å. However, recall that we neglect wavelengths
$\leqslant$6400 Å in our computation of the likelihood at 3.4 days to obtain a
fit which accurately captures the Sr absorption. Nonetheless, the omission of
Ba improves the fit, suggesting that we overestimate the abundance of Ba in
our best fit abundance pattern. The other new element, Ce, produces absorption
at $\sim$7000 Å which is reprocessed and emitted at $\sim$8000 Å, and its
omission worsens the fit. Ce also introduces some absorption at $\sim$12,000
Å. Interestingly, this absorption may originate from astrophysically-measured
Ce II lines with rest wavelengths in the range of $\sim$$15,200$ $-$ $16,700$
Å (Cunha et al. 2017; Majewski et al. 2017), blueshifted due to the expansion
of the ejecta. Domoto et al. (2021) similarly noted that these Ce II lines may
be prominent in kilonova spectra. The observation of Ce is also consistent
with the kilonova parameter space clustering analysis of Ford et al. (2023),
which finds that Ce II may broadly be an important ion in kilonova spectra
from 1.4 to 3.4 days. Ce III may also be important at later epochs (Gillanders
et al. 2023).

We complement these leave-one-out analyses using the Spectral DEComposition
(SDEC) tool of TARDIS. SDEC allows us to measure which elements or ions absorb
and/or emit the greatest luminosity during a given TARDIS run. All SDEC plots,
for favored and disfavored models, are compiled in Figure 14 (Appendix C). In
both models at 1.4 days, we see that the absorption of photon packets is
dominated by just three singly-ionized species: Sr II, Y II, and Zr II.
Indeed, 98.9% of all absorbed luminosity is absorbed by just these three
species in the favored single-component model. These ions remain important at
2.4 days, though the impact of Zr II is less clear. We also see some minor
absorption from Ba II at 2.4 days, at the shortest wavelengths.

At 3.4 days, while Sr II remains important, the SDEC plots reveal the presence
of several new heavier ions. In the favored multi-component model, we see
clearer evidence for Ba II, with the caveat that the abundance of Ba is likely
overestimated. More interesting, we see absorption from four singly-ionized
lanthanides: Ce II, Nd II, Sm II, and Eu II (cerium, 58Ce; neodymium, 60Nd;
samarium, 62Sm; europium, 63Eu). The absorption from Ce is strongest.
Altogether, Ce II is responsible for 27.5% of the absorbed luminosity at this
epoch, compared to 20.4% from Nd II + Sm II + Eu II, and 30.3% from Sr II.
This is a clear indication of the presence of lanthanides in the ejecta.
Interestingly, Eu is a “pure” $r$-process element, in that it is produced in
negligible quantities by the $s$-process (e.g., Bisterzo et al. 2014). The
presence of Eu is thus further proof for the operation of the $r$-process in
NS-NS merger ejecta. Finally, at 3.4 days, we also see some light (5.2% of the
total) absorption from an ion of bismuth (83Bi), Bi II. This is the heaviest
element detected in any of our models, but we caution that its detection is
marginal (we measure $Y_{\mathrm{Bi},3.4}=7.84^{+24.6}_{-7.43}\times
10^{-6}$), and it does not leave a clear imprint on the leave-one-out spectra.

Considering both the leave-one-out and SDEC analyses, we confidently identify
Sr at all epochs. We solidify our identification of Y at 1.4 and 2.4 days. Ba
may also be present in the ejecta, but its abundance is likely overestimated
at 3.4 days, complicating this claim. Finally, at 3.4 days, we infer the
presence of singly-ionized species of the lanthanides Ce, Nd, Sm, and Eu, with
the detection of Ce being most concrete. This ensemble of lanthanides, which
was not present at 1.4 and 2.4 days, emerges at 3.4 days to significantly
shape the emergent spectrum.

###  4.3 Physical origin of the ejecta components

Inferring the parameters of the ejecta component(s) also allows us to map
these components to different ejection mechanisms in the NS-NS merger which
generated the kilonova. Given the inferred $Y_{e}$ and $s$ (we neglect
$v_{\mathrm{exp}}$, since it is poorly constrained), the ejecta at 1.4 and 2.4
days can be attributed to matter that has undergone strong neutrino
reprocessing. This ejecta can arise from the neutrino-reprocessed wind of a
hypermassive neutron star remnant plus disk, either magnetized (e.g., Combi &
Siegel 2023; Curtis et al. 2023b) or unmagnetized (e.g., Fujibayashi et al.
2023; Just et al. 2023). Alternatively, accretion disks around black
holes—evolved in MHD—also yield ejecta with a broad distribution of electron
fractions (e.g., Siegel & Metzger 2018; Fernández et al. 2019; Christie et al.
2019; Miller et al. 2019; Just et al. 2022; Fahlman & Fernández 2022; Hayashi
et al. 2023; Curtis et al. 2023a). An outflow from such an accretion disk
could in principle provide both the ejecta at 1.4 and 2.4 days, and, the
neutron-rich component which emerges at 3.4 days.

The lanthanide-bearing component at 3.4 days
($Y_{e,2}=0.161^{+0.149}_{-0.104}$, $s_{2}/k_{\mathrm{B}}=21.5^{+8.3}_{-9.5}$)
can also be accounted for by dynamical ejecta. Numerical relativity
simulations that include neutrino absorption typically find $Y_{e}\sim
0.15-0.25$ and $s/k_{\mathrm{B}}\sim 20$ (e.g., Zappa et al. 2023), consistent
with our inferred parameters. The slightly higher-$Y_{e}$ component at 3.4
days ($Y_{e,1}=0.228^{+0.073}_{-0.088}$,
$s_{1}/k_{\mathrm{B}}=14.9^{+8.0}_{-4.5}$) is also consistent with such
dynamical ejecta. However, the dynamical ejecta alone is unlikely to account
for most of the mass generating the kilonova, as numerical relativity
simulations that employ parameters consistent with GW170817 produce $\lesssim
10^{-2}M_{\odot}$ in this dynamical component (e.g., Shibata et al. 2017; Most
et al. 2019; Nedora et al. 2021). Our inferred masses $M_{\mathrm{phot},1}\sim
M_{\mathrm{phot},2}\sim 10^{-4}M_{\odot}$ are consistent with this bound, but
recall that we are only sensitive to the mass above the photosphere. Moreover,
light curve modelling has found that more mass is contained in the red than
the blue component, which is challenging to accomplish if our red component is
indeed dynamical ejecta. Modelling the later spectra of AT2017gfo, when we are
sensitive to more of the ejecta mass, will be important for understanding
whether this new, redder, lanthanide-rich component is more consistent with
dynamical ejecta or an outflow from an accretion disk.

##  5 Conclusions

We fit single- and multi-component ejecta models to the spectra of the
GW170817 kilonova, AT2017gfo, during the early, optically thick phase at 1.4,
2.4, and 3.4 days post-merger. With these fits, we infer the element-by-
element abundance patterns at each of these epochs. We find that a single-
component model is favored at 1.4 and 2.4 days, while a multi-component model
is favored at 3.4 days.

This single component at 1.4 and 2.4 days is characterized by a high electron
fraction $Y_{e}\sim 0.3$ and moderate specific entropy $s/k_{\mathrm{B}}\sim
13-18$, yielding an ejecta dominated by lighter $r$-process elements and a
blue kilonova. This ejecta is consistent with material which has undergone
substantial neutrino reprocessing, e.g., winds from a remnant hypermassive
neutron star and/or an accretion disk. The multi-component ejecta at 3.4 days
contains a higher $Y_{e}\sim 0.23$ component and lower $Y_{e}\sim 0.16$
component, with entropies in the range $s/k_{\mathrm{B}}\sim 15-22$. These new
components contain heavier elements, and especially the lanthanides, yielding
a redder kilonova. The most substantial contribution of lanthanides comes from
the lower-$Y_{e}$ component, which is consistent with either dynamical ejecta
or a neutron-rich outflow from a remnant accretion disk.

The emergence of a new red component at 3.4 days is broadly in agreement with
modelling of the light curves of AT2017gfo. Physically, we infer that the
photosphere recedes into the ejecta over time: from $0.31c$ at 1.4 days, to
$0.25c$ at 2.4 days, to $0.21c$ at 3.4 days. This recession of the
photosphere, and the dimming of the earlier, blue component, reveals this new
red component that was not inferred at earlier epochs.

Using both leave-one-out and Spectral DEComposition analyses, we assess the
contributions of individual elements to the emergent spectra. We find that
strontium Sr produces the $\sim$8000 Å absorption at each epoch. We also
strengthen our identification of yttrium Y at short wavelengths $\lesssim$4500
Å at 1.4 and 2.4 days. At 3.4 days, this absorption from Y is less clear, as
absorption from an ensemble of lanthanides (cerium Ce, neodymium Nd, samarium
Sm, and europium Eu) dominates the absorption at these short wavelengths. The
identification of Ce is most concrete, producing absorption at $\sim$7000 Å
and $\sim$12,000 Å.

The abundance patterns at 1.4 and 2.4 days show a dearth of lanthanides and
heavier elements which makes these inconsistent with the Solar $r$-process
abundance pattern. However, at 3.4 days, the emergence of ejecta with
lanthanide fraction $X_{\mathrm{lan}}\sim-1$ substantially improves the
agreement between the inferred abundance pattern and the Solar $r$-process,
and the distribution of lanthanide fractions measured in metal-poor stars. The
better agreement between the Solar $r$-process/metal-poor stars and inferred
abundance pattern at 3.4 days, and the possible presence of multiple
components (remnant/disk wind and dynamical ejecta), lends more credence to
the ability of NS-NS/BH mergers to dominate the $r$-process in the Universe.

While we have fit spectra of AT2017gfo at 1.4, 2.4, and 3.4 days, the observed
spectra of AT2017gfo extend to 10.4 days. At later times, the ejecta is
expected to leave the photospheric phase and enter the optically-thin nebular
phase, when non-LTE effects become non-negligible. Given the modular nature of
SPARK, we could swap out TARDIS for a code specifically suited to non-LTE
radiative transfer, or, use the non-LTE functionality already available in
TARDIS. Our Bayesian inference framework would also enable quantitative model
comparison between models with and without non-LTE effects at epochs where
their relative importance is not yet understood (e.g., 4.4 and 5.4 days).
Fitting later epochs will also be crucial for inferring the abundances and
masses of all of the ejecta components, including those hidden below the
photosphere at 3.4 days. Beyond the later epochs of AT2017gfo, SPARK can
perform inference on the growing zoo of kilonova yielded by NS-NS/BH binaries
with different parameters, allowing us to understand the connections between
binary parameters, ejection mechanisms, elemental compositions, and
fundamental $r$-process conditions in the ejecta.

NV works in Tiohtiá:ke / Mooniyang, also known as Montréal, which lies on the
unceded land of the Haudenosaunee and Anishinaabeg nations. This work made use
of high-performance computing resources in Tiohtiá:ke / Mooniyang and in
Burnaby, British Columbia, the unceded land of the Coast Salish peoples,
including the Tsleil-Waututh, Kwikwetlem, Squamish, and Musqueam nations. We
acknowledge the ongoing struggles of Indigenous peoples on this land, and
elsewhere on Turtle Island, and hope for a future marked by true
reconciliation. We thank the attendees of an April 2023 workshop at University
of California—Santa Cruz, for fruitful discussions on kilonovae and the
$r$-process. We are also grateful to Shinya Wanajo for kindly sharing their
reaction network calculations. Finally, we thank Jessica Birky and David
Fleming for useful discussions on approximate Bayesian inference and the use
of approxposterior. This work made extensive use of the Narval and Cedar
clusters of the Digital Research Alliance of Canada at the École de
technologie supérieure and Simon Fraser University (with regional partner
WestGrid), respectively. We thank the support staff of Calcul Québec in
particular for their assistance at various steps in this project. This work
made use of the Vienna Atomic Line Database (VALD), operated at Uppsala
University, the Institute of Astronomy RAS in Moscow, and the University of
Vienna. We thank Nikolai Piskunov and Eric Stempels for help in obtaining the
VALD data. This research also made use of TARDIS, a community-developed
software package for spectral synthesis in supernovae (Kerzendorf & Sim 2014).
The development of TARDIS received support from the Google Summer of Code
initiative and from the European Space Agency (ESA)’s Summer of Code in Space
program. TARDIS is a fiscally sponsored project of NumFOCUS. TARDIS makes
extensive use of astropy and PyNE. We thank Andrew Fullard, Wolfgang
Kerzendorf, and the entire TARDIS development team for their assistance and
their commitment to the development and maintenance of the code. N.V.
acknowledges funding from the Natural Sciences and Engineering Research
Council of Canada (NSERC) Canada Graduate Scholarship - Doctoral (CGS-D), the
Murata Family Fellowship, and the the Bob Wares Science Innovation Prospectors
Fund. J.J.R. and D.H. acknowledge support from the Canada Research Chairs
(CRC) program, the NSERC Discovery Grant program, the FRQNT Nouveaux
Chercheurs Grant program, and the Canadian Institute for Advanced Research
(CIFAR). J.J.R. acknowledges funding from the Canada Foundation for Innovation
(CFI), and the Québec Ministère de l’Économie et de l’Innovation. N.M.F.
acknowledges funding from the Fondes de Recherche Nature et Technologies
(FRQNT) Doctoral research scholarship number 328732. M.R.D. acknowledges
support from the NSERC through grant RGPIN-2019-06186, the Canada Research
Chairs Program, and the Dunlap Institute at the University of Toronto. R.F.
acknowledges support from NSERC of Canada through Discovery Grant
RGPIN-2022-03463.

## References

  * Abbott et al. (2017a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, AJ, 848, L12
  * Abbott et al. (2017b) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, ApJ, 848, L13
  * Almualla et al. (2021) Almualla, M., Ning, Y., Bulla, M., et al. 2021, arXiv:2112.1547
  * Ambikasaran et al. (2015) Ambikasaran, S., Foreman-Mackey, D., Greengard, L., et al. 2015, IEEE Transactions on Pattern Analysis and Machine Intelligence, 38, 252
  * Andreoni et al. (2017) Andreoni, I., Ackley, K., Cooke, J., et al. 2017, PASA, 34, e069
  * Arcavi et al. (2017) Arcavi, I., Hosseinzadeh, G., Howell, D. A., et al. 2017, Nature, 551, 64
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Bauswein et al. (2013) Bauswein, A., Goriely, S., & Janka, H.-T. 2013, ApJ, 773, 78
  * Bisterzo et al. (2014) Bisterzo, S., Travaglio, C., Gallino, R., et al. 2014, ApJ, 787, 10
  * Breschi et al. (2021) Breschi, M., Perego, A., Bernuzzi, S., et al. 2021, MNRAS, 505, 1661
  * Buchner (2021) Buchner, J. 2021, The Journal of Open Source Software, 6, 3001
  * Bulla (2019) Bulla, M. 2019, MNRAS, 489, 5037
  * Chornock et al. (2017) Chornock, R., Berger, E., Kasen, D., et al. 2017, ApJ, 848, L19
  * Christie et al. (2019) Christie, I. M., Lalakos, A., Tchekhovskoy, A., Fernández, R., Foucart, F., Quataert, E., Kasen, D. 2019, MNRAS, 490, 4811
  * Combi & Siegel (2023) Combi, L., Siegel, D. M. 2023, arXiv:2303.12284
  * Coulter et al. (2017) Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Science, 358, 1556
  * Cowan et al. (2021) Cowan, J. J., Sneden, C., Lawler, J. E., et al. 2021, Reviews of Modern Physics, 93, 015002
  * Cunha et al. (2017) Cunha, K., Smith, V. V., Hasselquist, S., et al. 2017, ApJ, 844, 145
  * Curtis et al. (2023a) Curtis, S., Miller, J. M., Fröhlich, C., Sprouse, T., Lloyd-Ronning, N., Mumpower, M. 2023, ApJ, 945, L13
  * Curtis et al. (2023b) Curtis, S., Bosch, P., Mösta, P., et al. 2023, arXiv:2305.07738
  * Czekala et al. (2015) Czekala, I., Andrews, S. M., Mandel, K. S., et al. 2015, ApJ, 812, 128
  * Díaz et al. (2017) Díaz, M. C., Macri, L. M., Garcia Lambas, D., et al. 2017, ApJ, 848, L29
  * Domoto et al. (2021) Domoto, N., Tanaka, M., Wanajo, S., et al. 2021, ApJ, 913, 26
  * Domoto et al. (2022) Domoto, N., Tanaka, M., Kato, D., et al. 2022, ApJ, 939, 8
  * Drout et al. (2017) Drout, M. R., Piro, A. L., Shappee, B. J., et al. 2017, Science, 358, 6370, 1570-1574
  * Eichler et al. (1989) Eichler, D., Livio, M., Piran, T., et al. 1989, Nature, 340, 126
  * Evans et al. (2017) Evans, P. A., Cenko, S. B., Kennea, J. A., et al. 2017, Science, 358, 1565
  * Fahlman & Fernández (2022) Fahlman, S., & Ferández, R. 2022, MNRAS, 513, 2689
  * Fernández et al. (2019) Fernández, R., Tchekhovskoy, A., Quataert, E., et al. 2019, MNRAS, 482, 3373
  * Fleming & VanderPlas (2018) Fleming, D. P., & VanderPlas, J. 2018, The Journal of Open Source Software, 3, 781
  * Fleming et al. (2020) Fleming, D. P., Barnes, R., Luger, R., et al. 2020, ApJ, 891, 155
  * Ford et al. (2023) Ford, N. M., Vieira, N., Ruan, J. J., et al. 2023, arXiv:2308.15657
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24
  * Freiburghaus et al. (1999) Freiburghaus, C., Rosswog, S., & Thielemann, F.-K. 1999, ApJ, 525, L121
  * Fujibayashi et al. (2023) Fujibayashi, S., Kiuchi, K., Wanajo, S., et al. 2023, ApJ, 942, 39
  * Gillanders et al. (2022) Gillanders, J. H., Smartt, S. J., Sim, S. A., et al. 2022, MNRAS, 515, 631
  * Gillanders et al. (2023) Gillanders, J. H., Sim, S. A., Smartt, S. J., et al. 2023, arXiv:2306.15055
  * Goriely et al. (2011) Goriely, S., Bauswein, A., & Janka, H.-T. 2011, ApJ, 738, L32
  * Hayashi et al. (2023) Hayashi, K., Kiuchi, K., Kyutoku, K., Sekiguchi, Y., Shibata, M. 2023, Phys. Rev. D, 107, 123001
  * Hotokezaka et al. (2023) Hotokezaka, K., Tanaka, M., Kato, D., et al. 2023, arXiv:2307.00988
  * Hu et al. (2017) Hu, L., Wu, X., Andreoni, I., et al. 2017, Science Bulletin, 62, 1433
  * Ji et al. (2019) Ji, A. P., Drout, M. R., & Hansen, T. T. 2019, ApJ, 882, 40
  * Just et al. (2022) Just, O., Goriely, S., Janka, H.-T., et al. 2022, MNRAS, 509, 1377
  * Just et al. (2023) Just, O., Vijayan, V., Xiong, Z., Goriely, S., Soultanis, T., Bauswein, A., Guilet, J., Janka, H. -Th., Martínez-Pinedo, G. 2023, ApJ, 951, L12
  * Kandasamy et al. (2017) Kandasamy, K., Schneider, J., Póczos, B. 2017, Artificial Intelligence, 243
  * Kasen et al. (2017) Kasen, D., Metzger, B., Barnes, J., et al. 2017, Nature, 551, 80
  * Kasliwal et al. (2017) Kasliwal, M. M., Nakar, E., Singer, L. P., et al. 2017, Science, 358, 1559
  * Kawaguchi et al. (2020) Kawaguchi, K., Shibata, M., & Tanaka, M. 2020, ApJ, 889, 171
  * Kerzendorf & Sim (2014) Kerzendorf, W. E., & Sim, S. A. 2014, MNRAS, 440, 387
  * Kerzendorf et al. (2023) Kerzendorf, W., Sim, S., Vogl, C., et al. 2023, Zenodo
  * Korobkin et al. (2012) Korobkin, O., Rosswog, S., Arcones, A., Winteler, C., et al. 2012, MNRAS, 426, 3, 1940-1949
  * Lattimer & Schramm (1974) Lattimer, J. M., & Schramm, D. N. 1974, ApJ, 192, L145
  * Lipunov et al. (2017) Lipunov, V. M., Gorbovskoy, E., Kornilov, V. G., et al. 2017, ApJ, 850, L1
  * Lodders et al. (2009) Lodders, K., Palme, H., & Gail, H.-P. 2009, Landolt-Börnstein, 4B, 712
  * Majewski et al. (2017) Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94
  * Miller et al. (2019) Miller, J. M., Ryan, B. R., Dolence, J. C., et al. 2019, Phys. Rev. D, 100, 023008
  * Most et al. (2019) Most, E. R., Papenfort, L. J., Tsokaros, A., et al. 2019, ApJ, 884, 40
  * Nedora et al. (2021) Nedora, V., Bernuzzi, S., Radice, D., et al. 2021, ApJ, 906, 98
  * Nelder & Mead (1965) Nelder, J. A., Mead, R. 1965, Computer Journal, 7, 308
  * Pakhomov et al. (2019) Pakhomov, Y. V., Ryabchikova, T. A., & Piskunov, N. E. 2019, Astronomy Reports, 63, 1010
  * Perego et al. (2022) Perego, A., Vescovi, D., Fiore, A., et al. 2022, ApJ, 925, 22
  * Petelet et al. (2009) Petelet, M., Iooss, B., Asserin, O., et al. 2009, arXiv:0909.0329
  * Pian et al. (2017) Pian, E., D’Avanzo, P., Benetti, S., et al. 2017, Nature, 551, 67
  * Pozanenko et al. (2018) Pozanenko, A. S., Barkov, M. V., Minaev, P. Y., et al. 2018, ApJ, 852, L30
  * Ristic et al. (2022) Ristic, M., Champion, E., O’Shaughnessy, R., et al. 2022, Physical Review Research, 4, 013046
  * Ryabchikova et al. (2015) Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Phys. Scr, 90, 054005
  * Savitzky & Golay (1964) Savitzky, A. & Golay, M. J. E. 1964, Analytical Chemistry, 36, 1627
  * Shappee et al. (2017) Shappee, B. J., Simon, J. D., Drout, M. R., et al. 2017, Science, 358, 1574
  * Shibata et al. (2017) Shibata, M., Fujibayashi, S., Hotokezaka, K., Kiuchi, K., Kyutoku, K., Sekiguchi, Y., Tanaka, M. 2017, Phys. Rev. D, 96, 123012
  * Siegel & Metzger (2018) Siegel, D. M. & Metzger, B. D. 2018, ApJ, 858, 52
  * Smartt et al. (2017) Smartt, S. J., Chen, T.-W., Jerkstrand, A., et al. 2017, Nature, 551, 75
  * Sneppen et al. (2023) Sneppen, A., Watson, D., Bauswein, A., et al. 2023, Nature, 614, 436
  * Sneppen & Watson (2023) Sneppen, A. & Watson, D. 2023, arXiv:2306.14942
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132
  * Symbalisty & Schramm (1982) Symbalisty, E. & Schramm, D. N. 1982, Astrophys. Lett., 22, 143
  * Tanvir et al. (2017) Tanvir, N. R., Levan, A. J., González-Fernández, C., et al. 2017, ApJ, 848, L27
  * Tarumi et al. (2023) Tarumi, Y., Hotokezaka, K., Domoto, N., et al. 2023, arXiv:2302.13061. doi:10.48550/arXiv.2302.13061
  * Troja et al. (2017) Troja, E., Piro, L., van Eerten, H., et al. 2017, Nature, 551, 71
  * Utsumi et al. (2017) Utsumi, Y., Tanaka, M., Tominaga, N., et al. 2017, PASJ, 69, 101
  * Valenti et al. (2017) Valenti, S., Sand, D. J., Yang, S., et al. 2017, ApJ, 848, L24
  * van der Velden (2020) van der Velden, E. 2020, The Journal of Open Source Software, 5, 2004
  * Vieira et al. (2023) Vieira, N., Ruan, J. J., Haggard, D., et al. 2023, ApJ, 944, 123
  * Villar et al. (2017) Villar, V. A., Guillochon, J., Berger, E., et al. 2017 AJ, 851, L21
  * Wanajo (2018) Wanajo, S. 2018, ApJ, 868, 65
  * Watson et al. (2019) Watson, D., Hansen, C. J., Selsing, J., et al. 2019, Nature, 574, 497
  * Zappa et al. (2023) Zappa, F., Bernuzzi, S., Radice, D., et al. 2023, MNRAS, 520, 1481

##  Appendix A All Training Spectra & Posteriors

Figure 6 shows all of the spectra used to construct the training sets in the
single-component SPARK runs at 1.4, 2.4, and 3.4 days post-merger. Figure 7
shows the same for the multi-component equivalents. These include the
$m_{0}=1500$ Latin Hypercube Samples, $m_{\mathrm{active}}$ active learning
samples, and the best fit spectrum, for each.

Posterior probabilities are computed for each of these spectra, and these are
used by the Gaussian Process to construct a surrogate posterior. All
posteriors are obtained through dynamic nested sampling of the surrogate
Gaussian Process using dynesty (Speagle 2020). Figures 8 and 9 show the full
posteriors for the new single-component fits to the 2.4 and 3.4 day spectra,
respectively. Both posteriors are multi-modal. We highlight the modes which
correspond to the best fits in Table 3. Figures 10, 11, and 12 show the full
posteriors for the new multi-component fits to the 1.4, 2.4, and 3.4 day
spectra, respectively. The medians of these posteriors are given as the best
fit parameters in Table 4.

Figure 6: All spectra in the training set, and the best fits, for the single-
component SPARK runs at 1.4, 2.4, and 3.4 days. Top left: Spectra
corresponding to the $m_{0}=1500$ initial Latin Hypercube samples to explore
parameter space and $m_{\mathrm{active}}=1140$ additional active learning
points selected using BAPE, for the fit to the 1.4 day spectrum of AT2017gfo,
taken from V23. We also show the best fit, as obtained from the median of the
posteriors. Top right: The same, for $m_{0}+m_{\mathrm{active}}=1500+600$
spectra for the 2.4 day fit. Bottom center: The same, for the 3.4 day fit. The
1.4 and 2.4 day single-component models (shown here) are favored over the
multi-component models; the 3.4 day single-component model is disfavored with
respect to the multi-component model shown in Figure 7.

Figure 7: All spectra in the training set, and the best fits, for the multi-
component SPARK runs at 1.4, 2.4, and 3.4 days. Top left: Spectra
corresponding to the $m_{0}=1500$ initial Latin Hypercube samples to explore
parameter space and $m_{\mathrm{active}}=1440$ additional active learning
points selected using BAPE, for the fit to the 1.4 day spectrum of AT2017gfo.
We also show the best fit, as obtained from the median of the posteriors. Top
right: The same, for the 2.4 day fit, where $m_{\mathrm{active}}=1020$ active
learning samples are required. Bottom center: The same for the 3.4 day fit,
where $m_{\mathrm{active}}=1140$ active learning samples are required. The 1.4
and 2.4 day multi-component models shown here are disfavored with respect to
the single-component models; the 3.4 day multi-component model shown here is
favored over the single-component. Figure 8: Posterior from the single-
component SPARK run at 2.4 days. All posteriors are obtained through dynamic
nested sampling of the surrogate Gaussian Process using dynesty (Speagle
2020). The posterior presents some bimodality, most evident in the
$s/k_{\mathrm{B}}$ dimension. We discard samples with
$s/k_{\mathrm{B}}\geqslant 25.0$ to select the lower-entropy mode of the
posterior, which is a better fit and is our preferred model in Table 3. Figure
9: Posterior from the single-component SPARK run at 3.4 days. The posterior at
this epoch is highly multi-modal, as most evident in the
$\rho_{0},~{}Y_{\mathrm{e}},~{}\mathrm{and}~{}s/k_{\mathrm{B}}$ dimensions.
Although none of these modes in the posterior correspond to a good fit, and
the multi-component model is favored at 3.4 days, the higher-$\rho_{0}$,
mid-$Y_{\mathrm{e}}$, low-$s$ mode produces the best possible fit with single-
component ejecta. We select this mode of the posterior by discarding samples
outside the range
$\log_{10}(\rho_{0}/\mathrm{g~{}cm^{-3}})\in[-15.0,-14.0]~{}\cup~{}Y_{\mathrm{e}}\in[0.13,0.29]~{}\cup~{}s/k_{\mathrm{B}}\in[10.0,25.0]$.
The parameters of this mode are given in Table 3. Figure 10: Posterior from
the multi-component SPARK run at 1.4 days. The region of parameter space
sampled using active learning is substantially larger than the higher-
probability region of the posterior, which is relatively narrow. Figure 11:
Posterior from the multi-component SPARK run at 2.4 days. The posterior
contains one lower-probability mode, most evident in the
$v_{\mathrm{inner},1}$ dimension. We ignore this mode, as the median of the
higher-probability mode yields our best fit, as given in Table 4. Figure 12:
Posterior from the multi-component SPARK run at 3.4 days. As at 2.4 days, the
posterior contains one lower-probability mode, most evident in the
$v_{\mathrm{inner},1}$ dimension. We once again ignore this mode, as the
median of the higher-probability mode yields our best fit, as given in Table
4.

##  Appendix B Comparison of multi-component abundances

Figure 13: Comparison of the abundances of different components for the multi-
component models at 1.4, 2.4, and 3.4 days. The abundance patterns for the two
components are similar at 2.4 days. At 3.4 days, there is a lower-$Y_{e}$
component with more lanthanides and a relative dearth of Sr, while the
higher-$Y_{e}$ component contains $\sim$10$\times$ as much Sr. Importantly,
these plots show abundances, and not mass fractions. The red component at 1.4
days has $\sim$100$\times$ less mass than the blue and its contribution to the
emergent spectrum is negligible, despite having a substantially different
abundance pattern. The components are approximately equipartitioned in mass at
2.4 and 3.4 days, such that both contribute to the emergent spectrum.

Figure 13 compares the abundances of each component in the multi-component
models at 1.4, 2.4, and 3.4 days. The abundance patterns for the two
components are highly similar at 2.4 days. At 3.4 days, there is a
lower-$Y_{e}$ component with more lanthanides and a relative dearth of Sr,
while the higher-$Y_{e}$ component contains $\sim$10$\times$ as much Sr.
Crucially, these are abundances, not mass fractions: while the components are
approximately equipartitioned in mass at 2.4 and 3.4 days, at 1.4 days, the
blue has $\sim$100$\times$ the mass of the red such that the red has
negligible impact on the emergent spectrum.

##  Appendix C Spectral DEComposition Plots

Figure 14 shows the Spectral DEComposition (SDEC) plots for all 1.4, 2.4, and
3.4 day single- and multi-component fits. These SDEC plots show the species
(elements or ions) which dominate the absorption (measured by absorbed
luminosity) and emission (measured by emitted luminosity) in the spectrum.
Above the horizontal axis, the height of the histogram indicates the total
emitted luminosity in some wavelength bin, across different species. Below the
horizontal axis, the depth of the histogram indicates the total absorbed
luminosity. These plots also show how emission may be reprocessed, i.e.,
absorbed in some wavelength range and re-emitted in another. These are
complementary to leave-one-out plots (Figure 5), which highlight the elements
that have the largest impact on the shape of the spectrum, in terms of either
absorption or emission.

single-component | multi-component |
---|---|---
|  |  1.4 days
|  |  2.4 days
|  |  3.4 days
Figure 14: Spectral DEComposition (SDEC) plots for all best fits, at 1.4, 2.4,
and 3.4 days, for single- and multi-component models. The left column contains
single-component models for 1.4, 2.4, and 3.4 days; the right contains multi-
component equivalents. The ions Sr II, Y II, and Zr II produce significant
absorption in all 1.4 day and 2.4 day models. The 2.4 day models also show
absorption from Ba II. The 3.4 day single-component model (which is a poor fit
to the observed spectrum, but nonetheless the best fit possible with a single
component) exhibits absorption from all four of these ions, except Zr I in the
place of Zr II, due to the cooler ejecta. The spectrum also shows moderate
absorption from ions of three lanthanides: La II, Ce II, and Nd II. The 3.4
multi-component model, which is favored over the single-component, shows
absorption from two open $s$-shell elements in the same periodic table group
(Sr II and Ba II), four lanthanide ions (Ce II, Nd II, Sm II, and Eu II), and
one ion of an especially heavy element, Bi II. Furthermore, unlike at earlier
epochs, ions of Y and Zr do not dominate absorption in the 3.4 day multi-
component model.