Spectroscopic r-process Abundance Retrieval for Kilonovae. II. Lanthanides in the Inferred Abundance Patterns of Multicomponent Ejecta from the GW170817 Kilonova

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. 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 multicomponent 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 modeling, 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.


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).
The spectra of kilonovae in particular are marked by absorption and emission features from a suite of rprocess 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 modeling, 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 modeling 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 as well as the physical ejection mechanisms at play during these mergers.
Spectral modeling 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, 38 Sr) at ∼8000 Å, at 1.4, 2.4, and 3.4 days post-merger.Domoto et al. (2021Domoto et al. ( , 2022) ) similarly ascribe this feature to Sr II, and find tentative evidence for doubly ionized lanthanides La III and Ce III (lanthanum, 57 La and cerium, that Y II is present at 4.4 and 5.4 days, producing a P Cygni feature at ∼7600 Å. Gillanders et al. (2022) also suggest 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 ∼15,800 Å and ∼20,700 Å beyond 3.4 days, up to 10.4 days.At ≳ 7 days, these may instead arise from intrinsically weak lines, with ions Te III and Te I (tellurium, 52 Te) and I II (iodine, 53 I) 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 B = 13.6 +4.1 −3.0 , which yielded an extremely low lanthanide fraction log 10 X lan = −6.74+1.15  −1.85 .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 outshined 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 ∼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 an 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 that has been substantially reprocessed by neutrinos, e.g., winds from a hypermassive NS remnant or an accretion disk.At later times, the spectrum might be dominated or better described by a different ejecta or a multicomponent ejecta configuration.
Here we explore the time evolution of the inferred abundance pattern and the need for multicomponent ejecta models to fit the spectra of AT2017gfo, at 1.4, 2.4, and 3.4 days post-merger.We first extend our singlecomponent 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 multicomponent.We compare our best multicomponent fits at 1.4, 2.4, and 3.4 days to their single-component equivalents.The abundances of these single-and multicomponent 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 that enable fitting the later epochs of AT2017gfo with both single-and multicomponent 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.

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 (matterradiation) 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).This line list is identical to that used in V23, and includes only empirical lines, i.e., no lines from theoretical atomic structure calculations.Importantly, this line list is highly incomplete owing to the lack of empirical measurements of transitions and energy levels for many of the heavier elements.In associating any features in the spectrum with elements in our line list, we can only claim that a given element can produce a given feature, but not that it is the only element that could produce this feature.
Each spectrum is parameterized by a set θ i of parameters, including a luminosity, density, inner/outer computational boundary velocities, and three key parameters that set the abundances in the ejecta: the electron fraction Y e , expansion velocity v exp , and specific entropy per nucleon s/k B .The last 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, the total time for inference with more common methods such as Markov Chain Monte Carlo (MCMC) or nested sampling would be prohibitively long, due to the large numbers of samples required by these methods.Instead, we couple TARDIS to the approximate posterior estimation scheme of approxposterior (Fleming & VanderPlas 2018;Fleming et al. 2020).In this scheme, rather than directly evaluating the likelihood function, we introduce a Gaussian Process (GP) surrogate for the posterior L p (θ) and employ Bayesian Active Posterior Estimation (BAPE; Kandasamy et al. 2017) to intelligently sample parameter space.BAPE is a form of active learning in which we maximize an acquisition function with terms including both the mean µ(θ) and variance σ 2 (θ) 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 (θ, L p (θ)) 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 exp , and s/k B that best describe the ejecta, with relatively few forward model evaluations.We begin with a total of m 0 baseline samples over all parameter space obtained via Latin Hypercube Sampling (LHS), then follow with m active active learning samples obtained with BAPE.In V23, we fit the VLT/X-shooter spectrum of AT2017fgfo (Pian et al. 2017;Smartt et al. 2017) at 1.4 days with m 0 = 1500 and m active = 1140 samples.This is a factor of ∼ 103 fewer samples than might be required with a standard MCMC for a similar 6-dimensional fit.

Multicomponent, 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 outer , (which in fact sets the initial guess for the temperature at the inner boundary), the normalization in the density power law ρ 0 , the inner and outer boundary velocities v inner and v outer2 , and three parameters that set the abundance pattern: electron fraction Y e , expansion velocity v exp , and specific entropy s/k B .A single-component fit is thus 7-dimensional, unless one or more of the parameters are fixed.For example, we fixed v 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 multicomponent 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 outer .For these multi-shell runs, we must employ multiple TARDIS iterations when generating synthetic spectra.Beginning with a guess for the inner boundary temperature based on the user-requested value for L outer , at each iteration the plasma conditions (temperatures and dilution factors in each shell) are perturbed, and electron densities and level populations are recomputed.Photon packets are then propagated through all shells, and an estimate of the emergent spectrum is obtained.With sufficient iterations, the plasma converges to a state where the luminosity emitted at the outer boundary matches the user-requested L outer .TARDIS employs 20 such shells and 20 such iterations, by default.We use 10 shells to decrease computation times; in our tests, we find no marked difference in the resultant synthetic spectra with this smaller number of shells.We use 30 iterations to guarantee convergence of the plasma to a state described by the input parameters.
We begin with a simple two-component ejecta.As with our single-component model, the ejecta has an outer boundary luminosity L outer and power-law density profile with normalization ρ 0 .L outer describes the luminosity at the outer boundary of the overall ejecta, where photon packets are binned into a synthetic spectrum.The ejecta has an inner bound that is the lesser of v inner,1 , v inner,2 and an outer bound that is the greater of v outer,1 , v outer,2 .Each of two components is thus described by two velocities and three abundance-settingparameters Y e , v exp , and s/k B , i.e., five parameters each.Two-component fits are thus 2 + 5 + 5 = 12dimensional.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) that 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. 3In particular, we can produce the effect of reprocessing, where emission from one component is absorbed and reemitted/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 owing to the considerable opacity of the lanthanides in the near-UV and optical.Some kilonova spectra may be better described by these multicomponent ejecta models.In V23, we find that the 1.4-day spectrum of AT2017gfo is well-described by a single-component ejecta with Y e = 0.311 +0.013  −0.011 , where 0.311 is the median of our posterior distribution for Y e and the upper/lower bounds are the 97.5th and 2.5th percentiles, respectively.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 that were hidden or outshined at early times.Light-curve modeling (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 modeling (e.g., Kasen et al. 2017) also invokes multiple components.This motivates our introduction of multicomponent ejecta into SPARK.

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 ρ 0 and abundance-setting parameters Y e , v 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.However, all priors on luminosity are quite broad and based on the estimated bolometric luminosity of AT2017gfo (e.g., Wu et al. 2019) and our inference at 1.4 days.We further allow for wider priors on the inner and outer boundary velocities for the fits at later epochs.In V23, we fixed v outer = 0.35c during our 1.4-day fit but noted that we observed similar results for v outer in the range 0.35c − 0.38c.Here we allow for greater flexibility in v outer at later epochs.
Our priors for our multicomponent fits are given in Table 2.These are broad, and identical, at 1.4 and 2.4 days.Aside from requiring v inner,i < v outer,i , our priors also require some overlap between the two components, i.e., v outer,2 > v inner,1 or v outer,1 > v 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 LHS (Petelet et al. 2009).At 3.4 days, we encounter challenges with the convergence of the posterior for these broad, conditional priors in velocities.These broad, complicated priors result in a highly multimodal posterior, and the GP struggles to capture the relative amplitudes of the different modes, leading to inferred parameters that evolve as more samples are added to the GP.This challenge with GPs has been documented in El Gammal et al. (2023).In an effort to obtain a less complex posterior, we use tighter, nonconditional priors on the velocities and find that these lead to a convergent posterior.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 ≲6400 Å.Since we are interested in determining the species involved in the absorption, and especially the prominent feature at ∼8000 Å, we prioritize fitting the ⩾6400 Å region of the spectrum by excluding shorter wavelengths in the computation of the likelihood.We perform this exclusion for both single-and multicomponent fits.Finally, when expressing the likelihood using the Czekala et al. (2015) formalism, we introduce a global covariance term in the likelihood that introduces pixelto-pixel (i.e., wavelength bin-to-wavelength bin) correlations into the likelihood, beyond the standard individual uncorrelated uncertainties that would yield a simple χ 2 likelihood.This global covariance is described by a Matérn-3/2 kernel with some hyperparameters: an amplitude a G and correlation length scale ℓ.We use a global covariance term with amplitude a G = 10 −34 (erg s −1 cm −2 Å −1 ) 2 and a correlation length scale ℓ = 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 G = 10 −35 (erg s −1 cm −2 Å −1 ) 2 better matches the uncertainties on the observed spectrum at this epoch.

FITTING LATER EPOCHS AND MULTICOMPONENT EJECTA
We model the 1.4-, 2.4-, and 3.4-day spectra of AT2017gfo with both single-and multicomponent (twocomponent) models, and we 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 multicomponent models are included in Tables 3 and 4, respectively.In all cases, we use the median of the posterior as our best-fit parameter and the 2.5th / 97.5th quantiles of the posterior as the lower/upper bounds, respectively.For multimodal posteriors where we have selected some mode as our best fit, we use median and percentiles of that mode.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 multicomponent ejecta configuration.The element-byelement abundances with uncertainties for the best-fit models are included in Figure 2, and a detailed list of the mass fractions of each element is included in Appendix B. In the following sections, we present and examine each of the best fits in turn, at each epoch.
In all SPARK runs, we obtain a sufficient number of m 0 baseline samples and m active BAPE active learning samples such that the posterior converges, according to the z convergence diagnostic given in Fleming et al. (2020).Appendix A contains compilations of all m 0 + m active spectra (generated by LHS and BAPE active learning sampling, respectively) in a given SPARK run, for the 1.4-, 2.4-, and 3.4-day single-and multicomponent runs.
Posteriors generated by the single-and multicomponent runs are also included in Appendix A.
Among the parameters that set the abundances, we find that electron fraction Y e and entropy s are the most tightly constrained, while expansion velocity v exp is poorly constrained at all epochs.Given the importance of Y e and s in setting the abundances, we highlight the posteriors for (Y e , s) at all epochs for the favored models in Figure 3.We discuss the inferred Y e , s in greater detail in the following sections.We emphasize that the posterior distributions of the ejecta parameters should not be understood as the distribution of the ejecta parameters themselves, but rather a probabilistic distribution of the unique best-fit parameters for the ejecta.Our model still assumes a single Y e (or two, in the case of two components); the posterior for Y e is not the distribution of Y e in the ejecta itself, but rather quantifies the confidence of our best-fit parameters.

Single-and multicomponent modeling 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 active = 1500 + 1140 points to obtain a good fit and converged posterior.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 inner /c = 0.313 +0.013 −0.014 ; while the outer boundary velocity was fixed to v outer = 0.35c.For abundance-setting parameters, we found an electron fraction Y e = 0.311 +0.013 −0.011 , expansion velocity of v exp /c = 0.240 +0.055 −0.082 , and entropy of s/k B = 13.6 +4.1 −3.0 .In V23, we referred to this model as "purple + warm", but we 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 ( 38 Sr), allowing for a good fit to the prominent ∼8000 Å feature.The fit struggled at shorter wavelengths: over 3500 − 4500 Å, we overestimated the absorption from yttrium ( 39 Y) and zirconium ( 40 Zr), while at ≲3500 Å, we underestimated.Nonetheless, the singlecomponent model remains our preferred model at 1.4 days, as we elaborate on below.
In our 1.4-day multicomponent SPARK run, we require m 0 + m active = 1500 + 1440 points to obtain a good fit and converged posterior.More active learning samples are required than for the single-component model owing to the increased dimensionality of the multicomponent fits.The resultant fit generally captures the shape of the spectrum and the broad absorption at ∼8000 Å, like the single-component equivalent.However, the fit does not capture the continuum at wavelengths ≳10,500 Å nor 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 multicomponent models.The preferred models are highlighted: opaque lines are the preferred fits.We justify their preference in Section 3 (especially 3.4).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 multicomponent 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 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. The preferred fits capture the shape of the continuum and the dominant absorption feature at ∼8000 Å, but struggle to capture the emission bump of this feature in the interpretation that it is a P Cygni feature.
∼ 3500 -5000 Å.The fit is marginally closer to the observed spectrum at the shortest wavelengths (≲ 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 outer /L ⊙ ) = 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 outer is updated at each iteration.Hence, in our multicomponent models, the meaning of L outer is different and cannot be mapped directly to T inner .(L outer is also higher in the 2.4-and 3.4-day multicomponent models, compared to the single-component equivalents).
This multicomponent run yields two components with substantially different abundance patterns.A comparison of the abundance patterns is given in Figure 4.The first component has a lower electron fraction Y e,1 = 0.139 +0.028 −0.114 (and specific entropy s 1 /k B = 24.6 +3.4 −5.1 ) and could be described as red.The second has higher Y e,2 = 0.340 +0.022 −0.021 (and s 2 /k B = 20.1 +3.8 −3.6 ), distinctly bluer.Indeed, the Y e,2 of this second, bluer component is consistent 4 with the single-component fit at 1.4 days.As in the single-component fit, expansion velocities (v exp,1 and v 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 5.We see that the red component is confined to only two shells, while the blue component extends over a greater radius.The mass contained in these two components is M phot,1 /M ⊙ = 6.9 +22.1 −6.0 × 10 −7 and M phot,2 /M ⊙ = 7.0 +0.8 −1.0 × 10 −5 , respectively.We emphasize that this is the mass above the photosphere, and not the total ejecta mass.Nonetheless, the bluer component contains ∼100× as much mass as the redder component and dominates the absorption/emission in the spectrum.Indeed, while the redder component has a lanthanide mass fraction of log 10 X lan,1 = −1.24+0.68  −0.20 , the total lanthanide mass fraction of both components, computed as the mass-weighted sum over all shells, is log 10 X lan,total = −13.05+9.88  −3.32 due to the much larger mass of the bluer component.This small presence of lanthanides in the redder component might mark the photometric evolution of the kilonova, but it does not shape the spectrum.In this sense, despite the different compositions of the two components, this multicomponent SPARK run has effectively converged to a singlecomponent model.

Single-and multicomponent modeling at 2.4 days
In our 2.4-day single-component SPARK run, we require m 0 +m active = 1500+600 points to obtain a good fit and converged posterior.The resultant posterior shows some bimodality, most evident in the s/k B dimension.If we select the higher-entropy mode by selecting all samples with s/k B ⩾ 25.0, the resultant posterior has median Y e = 0.25 and v outer /c = 0.33.The synthetic spectrum generated at the median of this higher-entropy mode resembles a blackbody with some absorption at ∼5000 Å, which is not seen in the observed spectrum, and is a poor fit.Instead selecting samples from the posterior with s/k B ⩽ 25.0 (selecting the lower-entropy mode), we obtain a superior fit, and we use this as our preferred model.There is some arbitrariness in the cutoff s/k B ; we obtain similar results for cutoff s/k B in the range 22 − 28.This preferred fit captures most of the continuum of the observed spectrum and most of the ∼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 ≲4500 Å.
We infer n log 10 (L outer /L ⊙ ) = 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 −0.032 and v outer = 0.342 +0.047 −0.050 .This v outer is consistent with the fixed v outer = 0.35c from 1.4 days, while v 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 becoming optically thinner, as expected.
Finally, we obtain Y e = 0.306 +0.055 −0.204 and s/k B = 17.6 +7.1 −6.3 .Similar to our results at 1.4 days, v 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, multicomponent SPARK run, we require m 0 + m active = 1500 + 1020 points to obtain a good fit and converged posterior.The resultant fit once again captures the broad shape and ∼8000 Å absorption, but the depth of the ∼8000 Å absorption feature is overestimated.The multicomponent fit does, however, achieve a better fit to the continuum at wavelengths ≲7000 Å.
In contrast with the 1.4-day, multicomponent run, we see that the two inferred components are remarkably similar.Figure 4  Our inferred geometry for the two components is shown in Figure 5.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 phot,1 /M ⊙ = 2.3 +12.5 −0.7 × 10 −5 and M phot,2 /M ⊙ = 4.0 +4.8 −0.9 × 10 −5 , respectively.Given the two components' similar compositions and equal contributions to the emergent spectrum, this multicomponent SPARK run has also converged to an effectively single-component model.

Single-and multicomponent modeling 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 active = 1500 + 600 points to obtain a converged posterior.The resultant posterior shows a high degree of multimodality, most evident in the ρ 0 , Y e , and s/k B dimensions.Of the multiple modes, we obtain our best (but still poor) fit with a higher-ρ 0 , mid-Y e , lowers mode.This fit produces some, but not all, of the dominant absorption feature at ∼8000 Å.This fit also overestimates the absorption, or underestimates the continuum, for wavelengths ≲5500 Å.
We find log 10 (ρ 0 /g cm −3 ) = −14.586+0.313 −0.384 .This density is larger than that inferred at earlier epochs.We also infer v outer = 0.309 +0.032 −0.023 .This v outer is much smaller than the fixed value at 1.4 days or the inferred value at 2.4 days.This may explain the insufficient absorption at ∼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 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 ∼25% larger than that of our single-component model at 2.4 days.Nonetheless, Sr does not yield the expected prominent absorption feature at ∼8000 Å.This suggests that the shortcomings of this model may be due not to the inferred abundance pattern, but rather to the density, velocity structure, and/or some inherent shortcoming of the singlecomponent model.
For our preferred, multicomponent model, we require m 0 + m active = 1500 + 1140 points to obtain a good fit and converged posterior at 3.4 days.The resultant fit is better than the single-component equivalent.This fit better captures the ∼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 ≲5500 Å.This over/underestimation is more severe than in the singlecomponent model.
Interestingly, we see evolution in the inferred densities.Indeed, we find log 10 (ρ 0 /g cm −3 ) of −15.016 +0.320  −0.316 , −15.443 +0.742  −0.463 in our preferred models at 1.4 and 2.4 days, respectively, and −14.505 +0.323  −0.372 in this new multicomponent fit at 3.4 days.While this density at 3.4 days is higher than that inferred at earlier epochs, which may hint at the emergence of a new mass component, these densities are consistent with each other, within uncertainties.Note that the inferred total mass above the photosphere increases with time.However, this is due not to a lack of mass preservation, but rather to the recession of the photosphere into the ejecta.
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 B = 14.9 +8.0 −4.5 , while the latter has Y e,2 = 0.161 +0.149  −0.104 and s 2 /k 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 singlecomponent (and multicomponent) models at 1.4 and 2.4 days.The Y e,1 and Y e,2 are both substantially lower than those of the single-and multicomponent 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 lan,1 = −1.78+1.11  −3.69 and log 10 X lan,2 = −0.53+0.08 −2.86 , respectively.This yields a total lanthanide fraction of log 10 X lan,total = −1.58+0.86  −1.46 .Figure 4 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 ∼10× more of the crucial element Sr.
In Figure 5, we show a physical picture of the inferred ejecta.The photosphere further recedes com-  In contrast, the abundance patterns for the two components are similar at 2.4 days.At 3.4 days, there is a lower-Ye component with more lanthanides and a relative dearth of Sr, while the higher-Ye component contains ∼10× as much Sr.The components are approximately equipartitioned in mass at 2.4 and 3.4 days, such that both contribute to the emergent spectrum.Single-component models are favored over the multicomponent models shown here at 1.4 and 2.4 days, while the multicomponent model is favored at 3.4 days.pared to 1.4 and 2.4 days.Both components mostly overlap in space, and the mass is roughly equally split: they have masses M phot,1 /M ⊙ = 5.3 +1.8 −1.3 × 10 −4 and M phot,2 /M ⊙ = 3.8 +1.1 −0.9 × 10 −4 , respectively.Given the roughly equal contributions from a component rich in lanthanides and another with a similar abundance of lanthanides but ∼10× as much Sr, as well as the fact that the multicomponent model is a substantially better fit than the single-component model at 3.4 days, we interpret the ejecta as being genuinely composed of two unique components at this epoch.As we will see, the equal contribution from two components with distinct properties shapes the emergent spectrum.Here we take the best-fit parameters from our multicomponent runs (Table 4) and, after fitting, swap out the higher-Ye component with the purple + warm component (Ye = 0.311, vexp/c = 0.240, s/k B = 13.6).The legend indicates the parameters that have been changed to those of the purple + warm model.Louter, ρ 0 , and the velocities are left unchanged, and thus the mass in each component is also unchanged.Comparing to the original multicomponent models (dashed lines, Figure 1), the differences are negligible at all epochs.

Favored models
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 al-lowed to yield multiple components, the multicomponent SPARK runs yield effectively single-component best fits at 1.4 and 2.4 days.At 1.4 days, the bluer component contains ∼100× the mass of the redder one, and the presence of the red component is thus inconsequential for the spectrum.At 2.4 days, the two components are roughly equal in mass and characterized by similar abundance patterns; this fit is thus also effectively a single-component fit.
In contrast, at 3.4 days, the two components in the multicomponent 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 ∼10× more lanthanides than the higher-Y e component, given the broadness of the posterior.The remaining higher-Y e component, however, contains ∼10× 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 ∼8000 Å absorption, while the lower-Y e component yields most of the lanthanides needed to produce the short-wavelength ≲7500 Å absorption.
Aside from asking whether a single-or multicomponent ejecta is individually preferred at each epoch, we also ask, is the higher-Y e component of the multicomponent models interchangeable with the purple + warm ejecta model obtained at 1.4 days?To test this, we replace the higher-Y e components of the multicomponent models at 1.4, 2.4, and 3.4 days with this purple + warm component, while leaving the lower-Y e component intact.(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 6 shows the spectra that result from replacing the higher-Y e components of each multicomponent fit with the parameters of the purple + warm model (Y e = 0.311, v exp /c = 0.240, s/k B = 13.6) at each epoch.All other parameters are unchanged, i.e., we only change the abundance pattern of the higher-Y e component.This amounts to testing whether we produce a good fit with some hypothetical two-component model that includes the purple + warm component plus the lower-Y e component of each respective fit.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 higher 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 relatively unchanged owing 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 ∼8000 Å absorption, while the lower-Y e component provides sufficient lanthanides.
Overall, we infer that a new, redder, lanthanide-rich component emerges in addition to the previously inferred bluer ejecta 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 outshined by the brighter blue component.

Abundances of the favored models
Given our favored inferred models and their best fit (Y e , v exp , s), we study the inferred abundance pattern of the ejecta as a function of time.We favor the singlecomponent models at 1.4 and 2.4 days and the multicomponent model at 3.4 days.We show the abundance patterns of these favored models in Figure 2. The overall multicomponent 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 owing to the broader overall posterior.These 2.4-day abundances show evidence for ∼10 −9 to 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 rprocess elements around the first r-process peak.
In contrast, the 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 lan,total = −1.58+0.86  −1.46 , several orders of magnitude larger than inferred at 1.4 and 2.4 days.The higher-Y e component of this multicomponent ejecta then provides Sr, the role of which is crucial, as we will see in the following section.Indeed, we infer abundances of Y Sr,1.4 = 0.069 +0.053  −0.022 and Y Sr,2.4 = 0.087 +0.035 −0.059 at 1.4 and 2.4 days and Y Sr,3.4 = 0.067 +0.017 −0.046 (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 nonsolar in Figure 2, where the solar rprocess 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 rprocess abundances seen in several astrophysical settings.The 3.4-day abundance pattern, in contrast, is much closer to solar up to Z ∼ 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 modeling typically uses gray (wavelength-independent) opacities, this modeling generally recovers the emergence of multiple components, including a higher-opacity redder component at ≳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 lan ∼ −6 to ∼−4 for a blue component and ∼−2 for a red component (e.g., Chornock et al. 2017;Kasen et al. 2017).We find a more lanthanide-poor bluer component (log 10 X lan ∼ −7 as per the inference at 1.4 and 2.4 days) and a more lanthanide-rich (log 10 X lan ∼ −1.6, as per the inference at 3.4 days) redder component.Notably, Ji et al. (2019) compile the inferred lanthanide fractions of AT2017gfo from several studies and find that the log 10 X lan ∼ −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.

Elements and ions present in the favored models
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 7, we show leave-one-out spectra for the new favored models at 2.4 and 3.4 days and the 1.4- single component for 1.4 (top) and 2.4 (center) days, and multicomponent for 3.4 days (bottom).All models show clear absorption from strontium ( 38 Sr) at ∼8000 Å.At 3.4 days, when the ejecta is richer in heavier elements, we also see absorption (at ∼7000 and ∼12,000 Å) and emission (at ∼8000 Å) from the lanthanide cerium ( 58 Ce).However, we also see overabsorption from barium ( 56 Ba) at ∼4500Å.SDEC plots in Figure 8 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 ( 39 Y), at short wavelengths ≲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 ∼7600 Å at 4.4 and 5.4 days.In our favored multicomponent 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 ( 56 Ba), an open s-shell element in the same periodic table group as Sr, and cerium ( 58 Ce), a lanthanide.Ba is actually responsible for some of the overestimated absorption at ≲5000 Å.However, recall that we neglect wavelengths ⩽6400 Å in our computation of the likelihood at 3.4 days to obtain a fit that 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 ∼7000 Å that is reprocessed and emitted at ∼8000 Å, and its omission worsens the fit.Ce also introduces some absorption at ∼12,000 Å.Interestingly, this absorption may originate from astrophysically measured Ce II lines with rest wavelengths in the range of ∼15, 200 − 16, 700 Å (Cunha et al. 2017;Majewski et al. 2017), blueshifted owing 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 8.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 multicomponent model, we see clearer evidence for Ba II, with the caveat that the abundance of Ba is likely overestimated.More interestingly, we see absorption from four singly ionized lanthanides: Ce II, Nd II, Sm II, and Eu II (cerium, 58 Ce; neodymium, 60 Nd; samarium, 62 Sm; europium, 63 Eu).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: at the epoch of solar system formation, Eu was likely 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 ( 83 Bi), Bi II.This is the heaviest element detected in any of our models, but we caution that its detection is marginal (we measure Y Bi,3.4 ∼ 7.8 × 10 −6 ), and it does not leave a clear imprint on the leave-one-out spectra.Recall that all lines in our line list are empirical.
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.

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 that generated the kilonova.Given the inferred Y e and s see Figure 3; we neglect v exp here 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 NS remnant plus disk, either magnetized (e.g., Combi & Siegel 2023;Curtis et al. 2023b;Fahlman et al. 2023;Kiuchi et al. 2023) or unmagnetized (e.g., Fujibayashi et al. 2023;Just et al. 2023).Alternatively, accretion disks around BHs-evolved in MHDalso 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 that 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 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 ∼ 0.15−0.25 and s/k B ∼ 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 B = 14.9 +8.0 −4.5 ) may be consistent with a wind from an accretion disk around a BH or longer-lived NS, but it is also consistent with 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 as little as ≲ 10 −2 M ⊙ in this dynamical component and more than this amount in post-merger ejecta (e.g., Shibata et al. 2017;Most et al. 2019;Nedora et al. 2021).Our inferred masses M phot,1 ∼ M phot,2 ∼ 10 −4 M ⊙ are consistent with this bound, but recall that we are only sensitive to the mass above the photosphere.Moreover, some twocomponent light-curve modeling has found that more mass is contained in the red than the blue component (e.g., Chornock et al. 2017;Cowperthwaite et al. 2017;Villar et al. 2017), which is challenging to accomplish if our red component is indeed dynamical ejecta.Modeling 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.

Impact of nuclear and atomic physics
All of our results rely on (1) the use of a particular list of elements' energy levels and lines, i.e., a particular atomic physics, and (2) the results of the Wanajo (2018) nuclear network calculations, i.e., a particular nuclear physics.Changing this nuclear or atomic physics might alter the inferred parameters.As highlighted in Barnes et al. (2021), the photometric evolution of kilonovae is sensitive to the employed nuclear physics.Whether the nuclear physics inputs also impact spectral evolution remains to be seen, but is likely.In particular, the assumed nuclear mass models, decay rates, prescription for fission, and treatment of neutrino transport (if included) all enter into the nuclear reaction network calculations.These inputs in turn affect the mapping between elemental abundances and Y e , v exp , and s.
With respect to the atomic physics, the incompleteness of our line list, which uses only empirical lines as obtained from VALD, may bias the inference.Injecting theoretical lines into our list might worsen this bias, as it has been observed that different theoretical atomic structure calculations generate a large spread in line transition wavelengths and energies given the same initial configurations (e.g., Tanaka et al. 2020;Flörs et al. 2023).In practice, the use of a larger line list also significantly increases the runtime for spectral syntheses.Nevertheless, given our use of an incomplete but empirical line list, we can claim that a particular element can produce a particular feature in the observed spectra.We cannot claim that a particular element is the only element that can produce said feature.

CONCLUSIONS
We fit single-and multicomponent 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 elementby-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 multicomponent 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 ∼ 0.3 and moderate specific entropy s/k B ∼ 13−18, yielding an ejecta dominated by lighter r-process elements and a blue kilonova.This ejecta is consistent with material that has undergone substantial neutrino reprocessing, e.g., winds from a remnant hypermassive NS and/or an accretion disk.The multicomponent ejecta at 3.4 days contains a higher Y e ∼ 0.23 component and lower Y e ∼ 0.16 component, with entropies in the range s/k B ∼ 15 − 22.These new components contain heavier elements, 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 modeling 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, as well as 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 SDEC analyses, we assess the contributions of individual elements to the emergent spectra.We find that strontium Sr produces the ∼8000 Å absorption at each epoch.We also strengthen our identification of yttrium Y at short wavelengths ≲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 ∼7000 Å and ∼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 log 10 X lan ∼ −1.6 substantially improves the agreement between the inferred abundance pattern and the solar r-process, as well as the distribution of lanthanide fractions measured in metalpoor stars.The better agreement between the solar rprocess/metal-poor stars and inferred abundance pattern at 3.4 days, as well as 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.
All of our inference relies on the assumption of a particular atomic physics (line list of transitions and energy levels for ions) and nuclear physics (reaction network calculations).SPARK has been designed in a modular way such that either the atomic physics or nuclear physics, could be swapped out for another to quantify how the different physics impact the inference.A further simplification in our model is the use of single ejecta parameters (especially single Y e ) in generating synthetic spectra.In reality, the ejecta is described by a distribution of such parameters.In future models, it could be possible to introduce a toy distribution in the ejecta parameters (based on, e.g., the results of GRMHD simulations) and fit for the parameters of this distribution.
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 nonnegligible.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 kilonovae 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.versity of Vienna.We thank Nikolai Piskunov and Eric Stempels for help in obtaining the VALD data.
This research also made use of TARDIS, a communitydeveloped 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's (ESA) 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.Figure 9 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 10 shows the same for the multicomponent equivalents.These include the m 0 = 1500 Latin hypercube samples, m active active learning samples, and the best-fit spectrum for each SPARK run.
Posterior probabilities are computed for each of these spectra, and these are used by the GP to construct a surrogate posterior.All posteriors are obtained through dynamic nested sampling of the surrogate GP using dynesty (Speagle 2020).Figures 11 and 12 show the full posteriors for the new single-component fits to the 2.4-and 3.4-day spectra, respectively.Both posteriors are multimodal.We highlight the modes that correspond to the best fits in Table 3. Figures 13, 14, and 15 show the full posteriors for the new multicomponent 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.

B. ELEMENTAL MASS FRACTIONS
In Table 5, we compile the mass fractions of all elements and their uncertainties, for the favored models: single component at 1.4 and 2.4 days, and multicomponent at 3.4 days.We separately show the abundance of the individual components at 3.4 days and the massweighted sum of the two components.We omit elements' mass fractions when these are below 10 −9 .For elements for which the relative uncertainty on the mass fraction is > 100% the best-fit parameter, we use a tilde (∼) to indicate the large uncertainty.For elements for which the uncertainty is small (< 1%), we report just the best-fit value (i.e., posterior median).
We have used boldface for certain elements of interest, namely, Sr, Y, Zr, Ba, and Ce.We caution that we overestimate the absorption of Ba in our 3.4-day, multicomponent fit.The mass fraction of Ba should thus be taken conservatively as an upper limit.We select samples with s/k B ⩽ 25.0 to select the lower-entropy mode of the posterior, which is a better fit and is our preferred model in Table 3.The parameters of this mode are given in Table 3.  evident in the v 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 15.Posterior from the multicomponent SPARK run at 3.4 days.As at 2.4 days, the posterior contains one lower-probability mode, most evident in the v 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.

Figure 2 .
Figure2.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 multicomponent fit (see Figure1).The overall abundance of the multicomponent 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 posterior, effectively propagating the uncertainty on Ye, vexp, and s into the abundances.We also show the solar r-process pattern, computed using data fromLodders et al. (2009) subtracted by the s-process residuals ofBisterzo et al. (2014).The best-fit abundances are evidently nonsolar at the first two epochs but are closer to solar at the third.temperature of T inner = 3050 +126 −223 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 inner /c = 0.249 +0.017 −0.032 and v outer = 0.342 +0.047 −0.050 .This v outer is consistent with the fixed v outer = 0.35c from 1.4 days, while v 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 becoming optically thinner, as expected.Finally, we obtain Y e = 0.306 +0.055 −0.204 and s/k B = 17.6 +7.1 −6.3 .Similar to our results at 1.4 days, v 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, multicomponent SPARK run, we require m 0 + m active = 1500 + 1020 points to obtain a good fit and converged posterior.The resultant fit once again captures the broad shape and ∼8000 Å absorption, but the depth of the ∼8000 Å absorption feature is overestimated.The multicomponent fit does, however, achieve a better fit to the continuum at wavelengths ≲7000 Å.
compares the abundance patterns of both individual components.The first component is described by Y e,1 = 0.288 +0.129 −0.187 and s 1 /k B = 21.1 +10.5 −9.5 , while the latter has Y e,2 = 0.261 +0.163 −0.124 and s 2 /k 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 singlecomponent equivalent.Both entropies are consistent with the entropies inferred in the 1.4-day, multicomponent 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, multicomponent model and the 1.4-day, single-component (purple + warm) model.

Figure 3 .
Figure 3.Comparison of the posterior distributions of (Ye, s) for the favored models: single component at 1.4 and 2.4 days, and multicomponent at 3.4 days.The lower-Ye mode at 1.4 days is the preferred "purple + warm" model of V23.The two contours at 3.4 days describe the two different components, one of which has a higher Ye.The posteriors broaden at later epochs, but nonetheless, one of the components at 3.4 days has a substantially lower Ye than at any earlier epoch.Posteriors for all other ejecta parameters are included in Appendix A.

Figure 4 .
Figure 4. Comparison of the abundances of different components for the multicomponent models.Top to bottom: 1.4, 2.4, and 3.4 days.Importantly, these plots show abundances, and not mass fractions.The red component at 1.4 days has ∼100× less mass than the blue one and its contribution to the emergent spectrum is negligible, despite having a substantially different abundance pattern.In contrast, the abundance patterns for the two components are similar at 2.4 days.At 3.4 days, there is a lower-Ye component with more lanthanides and a relative dearth of Sr, while the higher-Ye component contains ∼10× as much Sr.The components are approximately equipartitioned in mass at 2.4 and 3.4 days, such that both contribute to the emergent spectrum.Single-component models are favored over the multicomponent models shown here at 1.4 and 2.4 days, while the multicomponent model is favored at 3.4 days.

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

Figure 6 .
Figure 6.Hypothetical models at 1.4, 2.4, and 3.4 days, where the higher-Ye component of each multicomponent fit has been replaced by that of the single-component, 1.4-day (purple + warm) model of V23.Here we take the best-fit parameters from our multicomponent runs (Table4) and, after fitting, swap out the higher-Ye component with the purple + warm component (Ye = 0.311, vexp/c = 0.240, s/k B = 13.6).The legend indicates the parameters that have been changed to those of the purple + warm model.Louter, ρ 0 , and the velocities are left unchanged, and thus the mass in each component is also unchanged.Comparing to the original multicomponent models (dashed lines, Figure1), the differences are negligible at all epochs.

Figure 7 .
Figure 7. Leave-one-out spectra for the favored models:

Figure 8 .
Figure8.SDEC plots for all best fits, at 1.4, 2.4, and 3.4 days, for single-and multicomponent models.The left column contains single-component models for 1.4, 2.4, and 3.4 days; the right contains multicomponent equivalents.The single-component models are favored at 1.4 and 2.4 days; the multicomponent is favored at 3.4 days.Above the horizontal axis, the height of the histogram indicates the total emitted luminosity in some wavelength bin, across different species, including electron scattering, and including photon packets that escape without interaction.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 reemitted in another.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-day multicomponent model, which is favored over the single-component model, 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 multicomponent model.Electron scattering is negligible compared to bound-bound interactions, in all models.

Figure 9 .
Figure 9.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 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 as the top left panel, but for m 0 + m active = 1500 + 600 spectra for the 2.4-day fit.Bottom center: the same as the top left panel, but for the 3.4-day fit.The 1.4-and 2.4-day single-component models (shown here) are favored over the multicomponent models; the 3.4-day single-component model is disfavored with respect to the multicomponent model shown in Figure 10.

Figure 10 .
Figure 10.All spectra in the training set, and the best fits, for the multicomponent 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 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 as the top left panel, but for the 2.4-day fit, where m active = 1020 active learning samples are required.Bottom center: the same as the top left panel, but for the 3.4-day fit, where m active = 1140 active learning samples are required.The 1.4-and 2.4-day multicomponent models shown here are disfavored with respect to the single-component models; the 3.4-day multicomponent model shown here is favored over the single-component model.

Figure 11 .
Figure 11.Posterior from the single-component SPARK run at 2.4 days.All posteriors are obtained through dynamic nested sampling of the surrogate GP using dynesty (Speagle 2020).The posterior presents some bimodality, most evident in the s/k B dimension.We select samples with s/k B ⩽ 25.0 to select the lower-entropy mode of the posterior, which is a better fit and is our preferred model in Table3.

Figure 12 .
Figure 12.Posterior from the single-component SPARK run at 3.4 days.The posterior at this epoch is highly multimodal, as most evident in the ρ 0 , Ye, and s/k B dimensions.Although none of these modes in the posterior correspond to a good fit and the multicomponent model is favored at 3.4 days, the higher-ρ 0 , mid-Ye, low-s mode produces the best possible fit with single-component ejecta.We select this mode of the posterior by selecting only samples in the range log 10 (ρ 0 /g cm −3 ) ∈ [−15.0,−14.0] ∪ Ye ∈ [0.13, 0.29] ∪ s/k B ∈ [10.0, 25.0].The parameters of this mode are given in Table3.

Figure 13 .
Figure13.Posterior from the multicomponent 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 14 .
Figure 14.Posterior from the multicomponent SPARK run at 2.4 days.The posterior contains one lower-probability mode, most

Table 1 .
Uniform priors for the single-component fits at 1.4, 2.4, and 3.4 days.vouter is fixed to 0.35c in the 1.4 day fit (see V23).

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.vouter is fixed to 0.35c in the 1.4-day fit.The inner boundary temperature Tinner and lanthanide mass fraction X lan are derived parameters, i.e., they are not dimensions in θ-space.

Table 4 .
Best-fit parameters for multicomponent fits to the GW170817 kilonova at 1.4, 2.4, and 3.4 days.The ejecta mass above the photosphere M phot,1;2 (not the total ejecta mass), lanthanide mass fractions X lan,1;2 , and the total X lan,total are derived parameters.

Table 5 .
Mass fractions Xi for the best-fit single-component models at 1.4 and 2.4 days, the two individual components at 3.4 days, and the mass-weighted sum of both components at 3.4 days.