Modelling the spectra of the kilonova AT2017gfo – II: Beyond the photospheric epochs

Binary neutron star mergers are the ﬁrst conﬁrmed site of rapid neutron capture ( 𝑟 -process) element nucleosynthesis. The kilonova AT2017gfo is the only electromagnetic counterpart of a neutron star merger spectroscopically observed. We analyse the entire spectral sequence of AT2017gfo (from merger to +10.4 days) and identify seven emission-like features. We conﬁrm that the prominent 1.08 𝜇 m feature can be explained by the Sr ii near-infrared triplet evolving from a P-Cygni proﬁle through to pure emission. We calculate the expected strength of the [Sr ii] doublet and show that its absence requires highly clumped ejecta. Near-infrared features at 1.58 and 2.07 𝜇 m emerge after three days and become more prominent as the spectra evolve. We model these as optically thick P-Cygni proﬁles and alternatively as pure emission features (with FWHM (cid:39) 35600 ± 6600 km s − 1 ), and favour the latter interpretation. The proﬁle of the strong 2.07 𝜇 m emission feature is best reproduced with two lines, centred at 2.059 and 2.135 𝜇 m. We search for candidate ions for all prominent features in the spectra. Strong, permitted transitions of La iii, Ce iii, Gd iii, Ra ii and Ac i are plausible candidates for the emission features. If any of these features are produced by intrinsically weak, forbidden transitions, we highlight candidate ions spanning the three 𝑟 -process peaks. The second 𝑟 -process peak elements Te and I have plausible matches to multiple features. We highlight the need for more detailed and quantitative atomic line transition data.


INTRODUCTION
In Gillanders et al. (2022), hereafter referred to as Paper I, we presented analysis of the early, photospheric spectra of the kilonova (KN) AT2017gfo, the electromagnetic (EM) counterpart to a binary neutron star merger (Abbott et al. 2017;Andreoni et al. 2017;Arcavi et al. 2017;Chornock et al. 2017;Coulter et al. 2017;Cowperthwaite et al. 2017;Drout et al. 2017;Evans et al. 2017;Kasliwal et al. 2017;Lipunov et al. 2017;Nicholl et al. 2017;Pian et al. 2017;Shappee et al. 2017;Soares-Santos et al. 2017;Smartt et al. 2017;Tanvir et al. 2017;Troja et al. 2017;Utsumi et al. 2017;Valenti et al. 2017).KNe are thought to be ideal locations for the synthesis of material by the -process, and many theoretical simulations have corroborated this (Lattimer & Schramm 1974;Eichler et al. 1989;Freiburghaus et al. 1999;Rosswog et al. 1999;Goriely et al. 2011Goriely et al. , 2013Goriely et al. , 2015;;Korobkin et al. 2012;Perego et al. 2014;Wanajo et al. 2014;Just et al. 2015).In Paper I, we focussed on identifying features that we can attribute to a single transition (or a set of transitions) belonging to a specific species.Through the use of a mix of publicly available atomic ★ E-mail: james.gillanders@physics.ox.ac.uk data sets, and tardis (a 1D Monte Carlo radiative transfer spectral synthesis code; Kerzendorf & Sim 2014), we were able to generate a set of models that self-consistently matched the observations of AT2017gfo, from +0.5 d, through to ∼ 1 week post-explosion.At this stage, the observed X-shooter spectra of AT2017gfo appear to transition into an optically thin regime, and the spectral features evolve into broad emission profiles.The photospheric approximation within tardis thus breaks down, and so tardis modelling is not likely to be valid beyond ∼ 6 − 7 days post-merger.
Our focus in Paper I was on studying specific species (e.g.Sr ii, through the identification of the near-infrared triplet producing the broad absorption feature between ∼ 0.7 − 1.0 m), or on identifying sets of elements (e.g. the necessary presence of the lanthanides in our models, combining to suppress the flux in the near-UV and optical parts of the spectra ≳ 2 days).In this paper, our work from Paper I is extended, but now we predominantly focus on the later phases.The project has two main goals.First, we examine the P-Cygni feature that is present between ∼ 0.7−1.2m.In early-phase modelling, this feature has been associated with Sr ii, and it remains persistent in the observed spectra until at least the +7.4 d spectrum.We aim to determine whether this feature is consistent with Sr ii producing an early P-Cygni feature, which evolves through to later phases, and into a pure emission feature, or whether we require contribution from additional species beyond the first few days.Second, we analyse the other spectral features, with the goal of determining ejecta properties, and constraining the composition of the ejecta.We model the two broad near-infrared (NIR) features, present in the intermediate-and latephase X-shooter spectra of AT2017gfo (+3.4 days onward), which we did not successfully reproduce with our tardis modelling in Paper I. We also perform empirical fitting of the late-time (+7.4 days onward) features to deduce ejecta properties, before performing a line identification study to shortlist the species that may be contributing to the observed data.Positive identification would allow us to constrain the composition of the ejecta.
The work presented here is a continuation of that presented in Paper I, and so the reader is directed to that manuscript for further details.Here we use the same atomic data set as Paper I, which made use of data from the Chianti (Dere et al. 1997;Del Zanna et al. 2021), Kurucz (Kurucz 2017, http://kurucz.harvard.edu/atoms.html),and dream data bases (Biémont et al. 1999;Quinet & Palmeri 2020), as well as the Pt and Au data presented by Gillanders et al. (2021) and McCann et al. (2022).In addition, we also make use of collisional data for Sr ii, presented by Bautista et al. (2002).We also use the same spectral data set, which we compiled from Pian et al. (2017) and Smartt et al. (2017), with several additions.Tanvir et al. (2017) and Troja et al. (2017) present Hubble Space Telescope (HST) spectra of AT2017gfo, at four separate epochs (+4.9, +7.3, +9.4 and +10.6 days).The spectra obtained at +4.9 and +9.4 d contain a broad, emission-like feature at ∼ 1.40 m, which corresponds to a telluric region.Hence, this feature is not seen in any of our X-shooter spectra at similar epochs.In the case of the +9.4 d spectrum, the epoch of observation agrees exactly with one of the X-shooter spectra, and so we merged this HST spectrum with the X-shooter spectrum at this epoch, replacing the pixels of X-shooter with the flux-calibrated HST pixels, within the telluric region.The +0.5 day spectrum from Shappee et al. (2017) is included, but it is blue and featureless, and so does not contribute to our analysis in this paper; however, we include it to illustrate the lack of prominent spectral features at very early times.
The manuscript is structured in the following manner.In Section 2 we briefly discuss the different features that are present in the entire observed spectral sequence of AT2017gfo, and their evolution.Section 3 contains our tardis analysis of the two NIR features that are prominent in the intermediate-phase (+3.4 d onward) X-shooter spectra.In Section 4, we investigate the presence of Sr ii at late times, and explore the possibility of it being responsible for the strong emission feature we observe at ∼ 1.08 m in the late-phase spectra.In Section 5 we model the most prominent of the features in the late-time spectra with Gaussians.In Sections 6 and 7, we perform searches for potential candidates that may be responsible for these strong emission features at late times.Finally, we summarise our results, and conclude in Section 8.Here we note that, unless otherwise stated in the text, all wavelengths are given as in vacuum.

SPECTROSCOPIC EVOLUTION OF AT2017GFO
The spectra of AT2017gfo, across all epochs observed (+0.5−10.4d) exhibit remarkably rapid evolution.Rapid cooling is observed and the centroids of the prominent features shift with time.This evolution is unrivalled by any other spectroscopically-observed extragalactic transient.At early times, the majority of the ejecta material is optically thick, and the spectra are dominated by contributions from the small mass of rapidly-expanding ejecta material that is present above the effective photosphere (see Paper I).As the KN ejecta expand and cool however, the outer material becomes transparent and so we can see deeper into the ejecta, meaning our spectra can now contain contributions from a broad range of ejecta material, spanning a large radial (or velocity) region.Given the inferred high velocities and low ejecta mass, it is of interest to investigate if/when we begin to observe evidence of a nebular spectrum.Recently, Pognan et al. (2023) proposed that KN ejecta material does not become completely optically thin very quickly.The results of their simulations show that while the rapidly-expanding outer ejecta quickly do become optically thin, it is possible for the innermost ejecta to remain optically thick out to +20 days.
Here we perform a simple exercise of determining where the spectral features peak, and how they evolve across multiple epochs.We can then attempt identification of the transitions that are responsible and determine information about the distribution of this material within the ejecta.The evolution of the peak wavelengths of the features can give us information about the velocity of the ejecta material contributing to the feature.
In Figure 1, we plot the sequence of AT2017gfo spectra we study, in scaled flux.This is similar to the spectral data set presented in Paper I, with some additions.First we plot the ten X-shooter spectra originally presented by Pian et al. (2017) and Smartt et al. (2017), and the early Magellan spectrum from Shappee et al. (2017).To this previous set, we have added the four HST spectra from Tanvir et al. (2017).Overlaid are lines to highlight the positions of the emission peaks in the spectra.Where features are present across multiple epochs, the lines are joined to illustrate how these positions evolve with time.Table 1 provides the approximate peak positions for the features highlighted in Figure 1.These approximate peak positions have been deduced by visual inspection.In Figure 2, we plot the peak positions of the features relative to their final peak position, to illustrate the evolution of the features and to visualise the data in Table 1.In general, the features shift redward with time, as we expect if the later spectra are dominated by contributions from material deeper within the ejecta (which would typically be expanding more slowly than the material observed at earlier times).
Previous works have drawn attention to some of these features.Smartt et al. (2017), Tanvir et al. (2017), Watson et al. (2019), Domoto et al. (2021Domoto et al. ( , 2022)), Perego et al. (2022), Sneppen et al. (2023), Sneppen & Watson (2023) and Tarumi et al. (2023) present analyses for some of them (and we highlight each work's specific contributions to the identification and study of these features below).However, no work to date has analysed all of the spectral emission-like features, and so here we present our own analysis and interpretation.
We identify six prominent spectral features, at ∼ 0.79 m (+2.4 − 8.4 d), ∼ 1.08 m (+2.4 − 7.4 d), ∼ 1.23 m (+4.4 − 7.4 d), ∼ 1.40 m (+4.9 and +9.4 d), ∼ 1.58 m (+3.4 − 10.4 d), and ∼ 2.07 m (+3.4 − 10.4 d).All features are present across a sequence of consecutive epochs, apart from the ∼ 1.40 m feature, which we only observe at +4.9 and +9.4 d.It is likely that this feature exists at intermediate epochs too, but since this region of our X-shooter spectra is affected by telluric absorption, we can only see it in the two HST spectra.
First we look at the evolution of the 0.79 m feature, which appears in the +2.4 − 8.4 d spectra (previously noted by Tanvir et al. 2017 andSneppen et al. 2023).We consider the most likely explanation of this feature as the edge of the blue wing of the strong Sr ii NIR triplet P-Cygni feature in the early phases (identified by Watson et al. 2019, see also Paper I).This point represents where the Sr ii feature blends into the continuum, but since the continuum The early Magellan spectrum is plotted (green), as are the spectral sequences obtained with X-shooter (black) and HST (blue).Each observation has been annotated with its phase, relative to the gravitational wave (GW) trigger.The peaks of the most prominent emission-like features have been estimated, and joined (red lines) across epochs, to highlight their evolution.The vertical grey bands correspond to telluric regions.
exhibits suppressed flux blueward of this, the location at which it joins the continuum becomes pronounced.Alternatively, it is possible that there is some as-yet unidentified species (unrelated to Sr ii) producing a pronounced emission feature (see Sneppen & Watson 2023, where they propose the feature is produced by Y ii).The peak of this feature shifts redward with time, indicating that the spectral features are being produced by sequentially slower-expanding material.If the material is expanding homologously, then this is simply a result of the outer regions of faster expanding ejecta becoming optically thin, allowing us to observe the effects of the inner regions of slower-expanding ejecta on the emergent spectrum.This leads to the blue wing not extending as far from the Sr ii NIR triplet rest wavelengths, which agrees with our model results presented in Paper I. Alternatively, the observed redward-shifting of the feature could be explained by an unidentified emission line or blend, which evolves redward with time.
Next we analyse the 1.08 m feature, which spans from +2.4 − 7.4 d (noted by Tanvir et al. 2017;analysed by Smartt et al. 2017;Watson et al. 2019;Domoto et al. 2021Domoto et al. , 2022;;Perego et al. 2022;Sneppen et al. 2023;Tarumi et al. 2023).Watson et al. (2019) proposed that this feature, at early times, can be readily explained by emission from the Sr ii NIR triplet.In Paper I, we confirmed  1).The horizontal grey dashed line marks a wavelength ratio of unity -the ratios for all features should tend towards this value (by design).
that a composition extracted from a realistic hydrodynamical simulation of a binary neutron star (BNS) merger (Goriely et al. 2011(Goriely et al. , 2013(Goriely et al. , 2015;;Bauswein et al. 2013) contained sufficient strontium for this feature to be produced by Sr ii.As the ejecta material evolves from its photospheric phase, through to a quasi-nebular regime, the peak of the feature emission should broadly line up with that of the rest wavelength of the transition (or set of transitions) that produce it.The weighted average of the expected centroid of the Sr ii NIR triplet blend (≃ 1.047 m), is blueward of the observed peak wavelength of the emission feature at late times (a point we will return to in Section 4).We note that in our tardis models presented in Paper I, we find generally good agreement between the peak posi-tions of the emission components of the Sr ii P-Cygni profiles produced by our sequence of models, and the observed emission peaks.Perego et al. (2022) and Tarumi et al. (2023) have suggested the He i  air = 10830 Å line as a possible alternative to the Sr ii NIR triplet (a point we return to in Section 4.4).We identify a weak emission component at ∼ 1.23 m, which appears at +4.4 days, and disappears beyond +7.4 days (also noted by Tanvir et al. 2017 andSneppen et al. 2023).It shifts redward with time, as seen for the other features.It appears as a weak feature, that just rises above the continuum redward of the broad and pronounced P-Cygni feature previously discussed.
As noted above, we find evidence for an emission feature at 1.40 m, but due to the telluric region in the X-shooter spectra, we cannot constrain its presence at any epochs apart from +4.9 and +9.4 d, where we are able to supplement our X-shooter observations with the HST spectra presented by Tanvir et al. (2017).This emission feature has been previously identified by Tanvir et al. (2017) and Sneppen et al. (2023).Domoto et al. (2022) model its absorption component (present at early times), and propose that this feature is produced by lanthanide absorption -specifically, La iii (a point we return to in Section 6).
Next we highlight the two emission features, present at ∼ 1.58 and ∼ 2.07 m.These appear in the +3.4 d spectrum, with both persisting until +10.4 d.These features were noted in Paper I, as they rose significantly above the NIR continuum, although our modelling did not explore their properties (we return to the nature of these two features in Section 3).The 1.58 m feature has been previously noted by Smartt et al. (2017) and Tanvir et al. (2017).Watson et al. (2019) and Sneppen et al. (2023) note the presence of both the 1.58 and 2.07 m features.None of these works identify the cause of these features.Domoto et al. (2022) model the +1.4 − 3.4 day X-shooter spectra, and show that a Ce iii P-Cygni feature can appear prominently at ∼ 14000 − 16000 Å in their models at these phases (see Section 6).These two emission features again follow the trend of shifting towards redder wavelengths as time evolves, for the reasons discussed for the previous features.
One of the key outstanding issues with interpreting the sources of these emission features is whether they originate from an optically thick or thin regime.To enable us to better understand these features, we consider the case that the emission features we have identified above are all the emission components of P-Cygni profiles, as are routinely produced from photospheric regimes (e.g. the Sr ii NIR triplet feature at ∼ 0.7−1.2m).If this is indeed the case, then we can reasonably expect to see evidence for absorption troughs blueward of the respective emission components.We can make predictions for the wavelength ranges these proposed absorption features span by considering the Sr ii feature.
If the ∼ 1.08 m feature is produced by the Sr ii NIR triplet (which has a weighted rest wavelength of ≃ 1.047 m), and the ∼ 0.79 m feature corresponds to the edge of our blueshifted absorption feature from this Sr ii P-Cygni profile, we are able to estimate the maximum ejecta velocity of the Sr ii material that is having an observable impact on the spectrum.Assuming that Sr ii is uniformly distributed throughout the ejecta, we use the inferred velocities extracted from the Sr ii P-Cygni feature as representative velocities for the ejecta (at least the ejecta within the line-forming region).We then use these velocities to constrain the positions of the absorption features for the other proposed P-Cygni profiles at each epoch the Sr ii blue wing feature is prominent (+2.4 − 8.4 d).From this, we can determine the expected positions of the absorption features (blue wings) of the other features in the spectra.For this calculation, we assume rest wavelengths of 1.24, 1.58 and 2.059 m for the features at ∼ 1.23, 1.58 and 2.07 m, respectively (from Tables 1 and 5).We present the results of this analysis in Figure 3.
The purple shaded region in Figure 3 corresponds exactly to the proposed absorption region of the Sr ii NIR triplet P-Cygni profile (by design), and acts to highlight the spectral shape we should see for the other features, if they are also P-Cygni features.The brown shaded region (∼ 1.23 m feature) does not show any evidence for an absorption trough, while the green and orange shaded regions (∼ 1.58 and 2.07 m features, respectively) are somewhat ambiguous.There does appear to be some flux deficit in these regions, but it is unclear whether this is a result of absorption, or whether these regions simply lie between prominent emission features, resulting in them appearing to have suppressed flux beneath some continuum.
One explanation for there being no convincing absorption components for these features is that emission from one (or many) other species is 'filling in', or replenishing, the flux in these wavelength regions.Alternatively, it may indicate that these features are not associated with P-Cygni profiles, and are indeed pure emission features, likely originating from some distinct component of optically thin ejecta material.Indeed, we proposed this may be the case in Paper I, since the shapes of our tardis model continua seemed to reasonably match the slopes of the observed NIR spectral regions (aside from the prominent emission components at ∼ 1.58 and 2.07 m).From this analysis, it appears we cannot rule either scenario out (optically thick versus optically thin), and so we explore both cases in greater detail in Sections 3 and 5, to better constrain which scenario best describes the data.

MODELLING THE EARLY PHASE FEATURES
There are two prominent emission-like features present in the NIR region, at  ≃ 1.58 and 2.07 m (see Section 2).It is unclear what produces these features, and whether each is produced by a single transition, or a set of transitions, belonging to one or multiple species.
Given the possibility that the 1.08 m feature might be (partly) due to the He i 2p 3 P -2s 3 S  air = 10830 Å transition (see Section 4), it is particularly interesting to consider whether the observed 2.07 m feature might be associated with the He i 2p 1 P -2s 1 S  air = 20581 Å transition (see also Tarumi et al. 2023, where they explore the presence of the 10830 and 20581 Å He i transitions in KNe).Such an identification does provide a rather good wavelength match with the spectral feature (within 10 Å of the peak position of the 2.059 m blue component of the blended 2.07 m feature; see Section 5), making it an intriguing possibility.However, the apparent relative strengths of the observed 1.08 and 2.059 m features does not seem to be compatible with expectations for these being produced consistently by He ii recombination.In particular, at +7.4 days, the approximate ratio of feature luminosities is ∼ 0.4 (see Table 5), while estimates for recombination line strengths suggest that the emissivity of the 10830 Å line should be at least an order of magnitude higher than the 20581 Å line (Benjamin et al. 1999).Particularly in view of the fact that the observed 2.059 m feature persists to later times after the 1.08 m feature has faded, it therefore seems unlikely that the 2.059 m spectral feature could be dominated by He i.
A critical factor in identifying these NIR features is the lack of line information at these long wavelengths.As noted in Paper I, our atomic data set is incomplete for species with  ≳ 40.Additionally, the heavy species for which we do have data typically only include We plot the early Magellan spectrum (green), as well as the sequences of X-shooter (black) and HST (blue) spectra.All observations are labelled with their phase, relative to the GW trigger.We also plot shaded regions corresponding to the locations where we should see any absorption components of P-Cygni features, with the intensity of the shading increasing with blueshift.The purple, brown, green and orange shaded regions correspond to the locations of proposed absorption for the ∼ 1.08, 1.23, 1.58 and 2.07 m features, respectively.The proposed rest wavelengths for these features are also plotted, as vertical dashed red lines.The vertical grey bands correspond to telluric regions.
transitions with wavelengths ≲ 1 m. 1 Therefore, it is likely that the currently available atomic data is insufficient to model these features.In Sections 6 and 7, we carry out a speculative search for candidate transitions -however, here we first attempt to establish the best framework to use to interpret the data; i.e. can they be explained by optically thick material producing P-Cygni features (similar to the 1.08 m feature)?If so, what constraints can we place on ejecta properties (e.g.distribution of the species involved)?

Motivation of parameters
We investigate whether these features can be formed in the same way as the Sr ii feature at the same epochs; i.e. whether these features are 1 For example, only three of the species that we included from the dream data base (Biémont et al. 1999;Quinet & Palmeri 2020, see Section 3 and Table 1 in Paper I for more on our atomic data set) have data for transitions at  > 1 m (La i, Pr iv and Yb ii).
produced by line scattering in the same region of the ejecta as Sr ii.
To constrain the nature of these features, we generate atomic data for two synthetic species.Both species are treated in a two-level atom approximation; i.e. we include a single excited level in each, where the energy gap between this level and the ground state corresponds to the wavelength of the observed features ( = 1.58 and 2.07 m).
The two-level approximation means that tardis will treat these as pure resonance scattering lines.We dub our mystery species X and Z, respectively. 2e adopt Einstein -values of 10 6 s −1 for each of these features, which is characteristic of strong dipole transitions at these wavelengths.We then incorporated these synthetic species into our atomic data, using the carsus package,3 and generated tardis models to determine if we could reproduce these NIR features in this way.We note that the work presented in this section is based on the tardis modelling we presented in Paper I, and where we refer to 'our bestfitting tardis models', we are referring to the models presented in Paper I that most closely matched the observational data.
We generate models for the +3.4−7.4 d AT2017gfo spectra, which correspond to the epochs modelled in Paper I that have prominent emission.The results of our modelling of these NIR features are shown in Figure 4. We plot the observed spectra, as well as the continua from our best-fitting tardis models presented in Paper I, for comparison.As noted in Paper I, our best-fitting models do not reproduce the flux level of the NIR continuum well.As a result, our model continua do not exactly match the observations at wavelengths ≳ 1.2 m.As there is some systematic uncertainty in the continuum position, we give ourselves freedom to shift the continuum of our model to improve agreement between the feature profiles from our mystery species X and Z, and the observations.
The mass fractions of the two-level atoms X and Z needed to produce features with strength similar to observations are presented in Table 2. 4 Interpretation of these mass fractions requires care, owing to the simplicity of the two-level atom approach used.We summarise the dominant factors that are relevant with: . (1) Here, MF tardis is the mass fraction of species X/Z in our tardis models, whereas MF real represents the mass fraction of the elements that produce X/Z in a real astrophysical explosion. ion is the fraction of the element that is ionised as species X/Z. l is the fractional population of the ion in the lower level for the transition. tardis is a variable that encapsulates the level populations as computed in our tardis models (in our two-level atoms, typically  tardis ∼ 1, and in all our calculations  tardis > 0.86; thus  tardis has only a modest impact on any interpretation of results). ion is the atomic mass number of species X/Z (the calculations are performed assuming a mass number  ion = 100, characteristic of elements between the first and second -process peaks).Although all the parameters in this expression may vary, it is reasonable to expect that the largest potential deviations are in the ion fraction ( ion ) and lower level population fraction ( l ).Our two-level atomic models effectively assume both of these parameters are (close to) 1, but both can be significantly smaller (e.g. the transition of interest may originate from a sub-dominant ion, or from a low-population excited state).red).Also plotted are the best-fitting tardis models with our mystery species X (blue dashed lines) and Z (green dash-dotted lines) included.All spectra have been plotted in scaled flux, and are offset for clarity (2.4,1.8, 1.2 and 0.6 for the third, fourth, fifth and sixth epoch spectra, respectively).The models with our mystery species X and Z have in some instances an additional additive offset applied, to better match the local continua around the observed features (see Table 2 for these model offsets).Telluric regions have been marked by grey shaded regions, and the observed flux in these regions are uncertain.The epochs have been annotated with their phase relative to GW detection.
Thus we expect that MF tardis represents a minimum mass fraction for the relevant species, since our two-level atom approach represents a best-case scenario for line formation.

Comparison between mystery species X and observations
First, we look at the feature at ∼ 1.58 m, which is produced by our mystery species X.For the models at both +3.4 and +4.4 days, we find that no continuum offset is necessary to match the observations.The models have some small component of absorption blueward of the emission feature, that mostly lies in the telluric region of the observational data. 5This prevents any definitive conclusion to be drawn with regards to the presence of a P-Cygni feature, versus a pure emission feature, although we do note that the spectrum immediately blueward 5 Here we include the data in the telluric regions, but stress that the flux within these regions is uncertain, and so should not be used for any inferences as to model agreement.of the emission feature (at ∼ 1.45 m) appears to rise to the blue, in apparent contradiction with our models, which have absorption features at this wavelength.Our model emission components appear to be much broader than the observed features, and so we do not get good agreement between our models and the observed spectra.This indicates that the emission feature in the observed spectrum is produced by material that is expanding more slowly than the velocities required to reproduce the optical region of the spectra (the Sr ii feature in particular).Domoto et al. (2022) analysed the +1.4 − 3.4 d spectra, assuming an optically thick photosphere, with a hybrid line list.They identify the same feature in the data, although they focus on modelling the wavelength position of the associated absorption component.Domoto et al. (2022) propose this feature is a blend of Ce iii lines, and their synthetic spectra produce a P-Cygni profile similar to ours.However, it appears to peak at a redder wavelength, and is broader than the observed feature.The discrepancy in feature width between the synthetic and observed spectra is similar to what we find; hence, we reserve judgement here on whether or not mystery species X can be definitively explained by these Ce iii transitions.
Looking at later times, we find reasonable agreement with our model and the ∼ 1.58 m feature at +5.4 days.Again we require no continuum offset, but we do require a lower mass fraction of X than in the previous two epochs (4 × 10 −5 at +5.4 days, versus 1.2 × 10 −4 for +3.4 and +4.4 d).Again our emission feature is broader than the observation, but we get good agreement with the observed data at ∼ 1.45 m, directly between the emission component and the telluric region.Our +6.4 d model requires a small continuum offset to better match observation, but the same mass fraction of X as needed at +5.4 d.Again, we reasonably match the flux between the telluric and the emission feature.Our emission feature is still broader than the observations.At +7.4 d, our model does not agree with the observations.Our model absorption feature has shifted to redder wavelengths (due to the receding photosphere within our tardis models, leading to higher densities at lower velocities, and therefore most of our mystery species in this model has slower velocity than previous phases), and now does not agree with the observed spectrum.We apply an offset to the model to better match the continuum redward of the emission feature, but overall the fit is poor.
Our model emission features being broader than the observed features indicates our model expansion velocities are too high.This implies that our mystery species X is not present in the fastestexpanding material, which corresponds to the regions of ejecta that Sr ii is present in (i.e. the regions that we modelled with tardis in Paper I, where we obtained good agreement to the observations, up to ≲ 1 week).This can be due to (i) ionisation effects, where at higher velocities, our mystery species is not the dominant ion, or (ii) ejecta stratification from different ejection mechanisms.

Comparison between mystery species Z and observations
Now we summarise the results of modelling the second emission feature at ∼ 2.07 m, with our mystery species Z. Across all epochs, our model emission components are broader than the observations, similar to the behaviour of X.Our model for +3.4 d requires a small offset to match the continuum blueward of the ∼ 2.07 m feature, but this leaves us with poor agreement to the data redward of the emission feature, which is perhaps evidence for needing a steeper continuum than our tardis model produces.We see similar behaviour at +4.4 d, where again we match the continuum at blue wavelengths, but not at redder wavelengths.The model emission component is still too broad to match the data.
At +5.4 d, we get good agreement to the data, with a reasonable match to the continuum either side of the observed emission feature.Additionally, the shape of our emission feature appears to broadly agree with the observed one.Our model predicts a reasonably strong absorption feature, which cannot be compared with observation since it lies in a telluric region.For the final two phases (+6.4 and +7.4 days), the quality of our fits degrades.Now we cannot obtain a good match to the continuum at both sides of the emission feature.At +6.4 d, the spectrum blueward of the emission feature appears to rise in flux (to the blue), at odds with our model, which predicts an absorption component at this position (∼ 1.95 m).To reduce this discrepancy, we lowered our Z mass fraction (from 3 × 10 −4 at +3.4 − 5.4 d, to 1.8 × 10 −4 ), but still cannot match the observed data well.At +7.4 d, we match the continuum at wavelengths blueward of the emission feature, but not redward, and our model now requires a much higher mass fraction of Z (6 × 10 −4 ) to produce an emission feature comparable in strength to the data.As a result, we now have a very pronounced absorption feature in our model that does not match the blue wing of the observed emission feature well.The fact our model requires such a large increase in mass fraction (≳ 3× increase) to reproduce the observed feature strength indicates that some additional species may be contributing to this feature at this epoch (see Section 5, although there we only see strong evidence for an additional emission component at +8.4 d and later).
We have shown with our tardis models that we cannot get good agreement with these features overall.Our model features are much broader than the equivalent features in the observed data, indicating that this ∼ 2.07 m feature originates from ejecta moving more slowly than that needed to reproduce the Sr ii feature.This could be due to ionisation or stratification effects, where Z is only present at lower velocities.

Summary
The analysis above demonstrates that we are unable to convincingly reproduce the observations of the ∼ 1.58 and 2.07 m features from +3.4 − 7.4 d under the assumption that these are produced by scattering lines forming in the same ejecta regions that are described by our tardis models for the Sr ii feature.Our model P-Cygni profiles generally do not match the observations, especially the absorption components.Therefore, we ultimately disfavour this scenario (i.e.optically thick ejecta producing P-Cygni features), and conclude that it is more likely that these features are in fact pure emission features, that are produced by either collisional excitation, recombination or fluorescence from some region of the ejecta, typically at lower velocity.In Section 5, we parameterise such a scenario by simple Gaussian fitting of the emission features.

ANALYSING THE PRESENCE OF STRONTIUM AT
LATE TIMES Watson et al. (2019) first presented evidence for the presence of Sr in the early spectra of AT2017gfo (later corroborated by Domoto et al. 2021, Perego et al. 2022, andPaper I).Specifically, they identified two features which they suggest are produced by Sr ii absorption.These features lie at ∼ 3500 Å (produced by the Sr ii resonance lines from the ground state, with wavelengths,  air = 4077.7 and 4215.5 Å) and ∼ 7000 − 10000 Å (produced by the Sr ii transitions from the metastable 4p6 4d levels, with wavelengths,  air = 10036.7,10327.3 and 10914.9Å).However, they did not present detailed predictions for how the Sr ii features may be expected to evolve at later times, beyond the early, photospheric phase.
In Paper I, we quantitatively confirmed that the Sr ii NIR triplet can reproduce the strong 1.08 m feature, but found that the effect of absorption from the resonance doublet was hidden by the strong line-blanketing in the blue ends of the spectra, dominated by Y ii and Zr ii in our early models (< 2 days), and the lanthanides in our later models (> 2 days).Here we analyse whether the spectral feature at ∼ 1.08 m (prominent in the +7.4 d X-shooter spectrum) is consistent with the Sr ii NIR triplet evolving from a P-Cygni profile in the early spectra, to a pure emission feature at this position after ∼ 7 days.
The Sr ii NIR triplet transitions, and their connection to the ground state, are analogous to the well-known Ca ii NIR triplet.The Ca ii NIR triplet is commonly observed to transition from a P-Cygni profile to a pure emission feature in supernova spectra (see e.g.Kumar et al. 2022, for a recent example of such behaviour for SN 2017iro, a well-sampled SN Ib).As the density drops, the [Ca ii] doublet (at  air = 7291.5 and 7323.9Å) becomes prominent as a strong emission feature through de-excitation processes.As a verification of our method (presented below), we are able to qualitatively reproduce this Ca ii evolutionary behaviour.If the 1.08 m feature observed in the photospheric spectra of AT2017gfo is indeed produced by Sr ii, as the density drops we expect to see this P-Cygni profile transition into emission and for the equivalent [Sr ii] doublet to appear, as we typically see for Ca ii.The laboratory-confirmed line centroids for the [Sr ii] doublet lie at wavelengths,  air = 6738.4and 6868.2Å.Here we analyse the +7.4−10.4d X-shooter spectra, and investigate if there is evidence for the Sr ii NIR triplet evolving into a pure emission feature, and for the appearance of the expected [Sr ii] emission.

Method
Once the ejecta become optically thin, the photospheric approximation used for modelling the early spectra breaks down, and can no longer be used.Additionally, we can no longer assume level populations remain close to their local thermal equilibrium (LTE) populations due to the fact that, in diffuse material, radiative effects will be comparable in strength to (or even dominate over) collisional effects.For example, Hotokezaka et al. (2021) and Pognan et al. (2022a,b) demonstrate that non-LTE effects may become important for determining accurate level populations after just a few days in KN ejecta.If the ejecta material is optically thin, then it will be dominated by emission, which may explain the presence of emission-like features in the spectra of AT2017gfo taken after ∼ 1 week.
For this study, we use the collisional data for Sr ii, as presented by Bautista et al. (2002). 6With these data, we are able to estimate the level populations for the lowest few levels of Sr ii, based on both radiative and collisional effects, which is more physically motivated than a simple LTE approximation.We do not discuss the details of calculating collision strengths or excitation and deexcitation rates here (readers are instead directed to Burgess & Tully 1992, for further information on the topic).However, in brief, the thermally averaged collision strengths, Υ, are used to derive the excitation ( lu ) and de-excitation ( ul ) rates for a level by collisions with free electrons. lu and  ul both depend on temperature ().For the following calculations, we assume a range of plausible temperatures,  ∈ [2500, 3500, 4500] K, for the late-time ejecta.From Paper I, our tardis inner boundary temperatures and velocities are  = 3200, 3100 and 2900 K, and  min = 0.09, 0.07 and 0.05 , for the +5.4,+6.4 and +7.4 day spectra, respectively.Accounting for Doppler corrections, these correspond to ejecta temperatures of ∼ 3000−3500 K.We propose that this range of temperatures explored by our models ( ∈ [2500, 3500, 4500] K) correspond to reasonable expectations for the ejecta temperature after ∼ 1 week, as well as sampling slightly cooler or hotter ejecta than expected.
We also require an estimate for the electron density.Jerkstrand (2017) presents a review of supernovae in the nebular phase, which includes a useful parameterisation of the electron density of latephase SN ejecta from a spherically symmetric distribution, given as: where  is the mean atomic weight,   is the proton mass,  is the ejecta mass of the system,  is the typical velocity of the ejecta,   is the electron fraction (  =   / nuclei ),  is the expansion time, or time since explosion (assuming homologous expansion), and  is the so-called filling factor (see the discussion below).
Here, we use this equation to estimate the electron density in KN ejecta at late times.Assuming  = 100,  = 10 −2 M ⊙ ,  = 0.1  (see Paper I, or any of Fernández & Metzger 2013;Metzger & Fernández 2014;Perego et al. 2014;Just et al. 2015;Wu et al. 2016;Siegel & Metzger 2017;Fujibayashi et al. 2018;Miller et al. 2019;Curtis et al. 2023;Just et al. 2022) and that, on average, each atom present in the ejecta has been singly ionised (  = 1), we obtain: Suitable values for  in KN ejecta are unknown, and so we assume two different scenarios -uniform density material (  = 1), and highly clumped material (  = 10 −2 ).We stress that when we refer to clumped material, we are not referring to large distinct clumps of ejecta with different compositions.Instead, we are referring to small ejecta clumps, which do not impact the column density along our line of sight, known as microclumping.These clumps however do lead to regions of over-and under-dense ejecta compared to a homogeneous distribution, which impact the relative importance of collisional to radiative processes.Microclumping has been invoked to explain observations of stellar winds (see e.g.Hamann & Koesterke 1998;Hillier & Miller 1999;Puls et al. 2006), and it has also been invoked empirically in the context of supernova ejecta (see e.g.Dessart et al. 2018;Wilk et al. 2020).We therefore argue that it is reasonable to consider the potential impact of microclumping in the kilonova Table 3. Computed electron densities at times corresponding to the phases of the late-time X-shooter spectra of AT2017gfo.
The thermally averaged collision strengths (Υ), combined with the three  values (2500, 3500 and 4500 K) enable us to determine  lu and  ul values (Burgess & Tully 1992).With these, the two sets of   estimates (Table 3), and the Einstein -values for all relevant transitions, we determine the level populations for the lowest five levels of Sr ii (as in, e.g.Dere et al. 1997).With quantitative estimates for level populations computed, we can then calculate the number of escaping photons, and the luminosities that we would expect for the emission lines ( em ).
We need to account for the possibility of these emitted photons being scattered or re-absorbed by another Sr ii ion before escaping.Therefore, we calculate the escape probability of the emitted photons,  esc , using: where  s is the Sobolev optical depth of the transition (see e.g.Noebauer & Sim 2019).We calculate  s from the level populations computed by assuming a uniform density sphere, with a radius derived from a constant expansion velocity,  = 0.1  (assuming homologous expansion), and with the temperatures and masses presented in Table 4.  esc decreases with increasing  s , and so optically thick material will prevent photons from escaping.Since the escape probability of some of the transitions are not of order unity, we need to take this effect into consideration.9To do this, we computed effective Einstein -values for all transitions ( esc ×  ul ).We then re-calculated the level populations with these effective -values, to obtain a truer representation of the level populations we would expect for Sr ii.Line luminosities and profiles are calculated (as in Gillanders et al. 2021) using: where  u is the number of atoms or ions in the excited state, and  is the wavelength corresponding to the transition.Gaussian line profiles for each Sr ii transition were simulated with these  em values and a full-width, half-maximum (FWHM) velocity of 0.1  (see Section 5).We scale the mass of Sr ii in our level population calculations, Table 4. Masses of Sr ii included in our synthetic emission spectra calculations, for both our uniform density (  = 1) and clumped (  = 10 −2 ) ejecta material.
such that the strongest emission features in our +7.4d models have luminosities comparable to the observed features in the spectrum of AT2017gfo at this epoch.These model masses remained fixed for the subsequent epochs, and are provided in Table 4.To put these masses in context, our best-fitting models from Paper I require a range of Sr ii masses (above the photosphere), from  Sr ii = 10 −7 M ⊙ for our +1.4d model, to  Sr ii = 2.25 × 10 −5 M ⊙ for the model evolved to +7.4 d.We expect our  Sr ii values here to be larger than these, since we are sensitive to the entire ejecta in this regime (as opposed to just the line-forming region in the tardis models).We impose an upper mass limit of  Sr ii = 10 −3 M ⊙ , to prevent our analysis exploring unphysically large mass estimates.
The individual Gaussian profiles were co-added to create a composite emission spectrum.Figure 5 shows the emission spectra we generated for our range of temperatures and electron densities.In the upper panel, we plot the emission spectra we generated, assuming uniform density ejecta material (  = 1), whereas the lower panel contains the models we generated assuming clumped ejecta material (  = 10 −2 ).The models in both panels have been over-plotted onto the observed late-phase spectra of AT2017gfo (+7.4 − 10.4 d), such that the models and observations that have the same phase are overlaid, to aid comparison.

Uniform density ejecta material ( 𝒇 = 1)
Our 3500 and 4500 K models, presented in the upper panel of Figure 5, produce a strong emission feature at ∼ 6800 Å that is not present in the observed spectra.This feature is produced by the [Sr ii] 6738.4,6868.2Å lines.Although it may be possible for lineblanketing from other species to suppress some of the flux from this feature, the likelihood that a feature of this strength could be so heavily negated, such that it has no observational trace, seems low.Alternatively, it may be possible to explain the lack of this feature by invoking most of the Sr ii material in the ejecta remains beneath the photosphere.However, this explanation becomes weaker for the later epochs, by which time the photosphere will have receded (or disappeared completely).Our model continues to predict this feature is prominent at these late times, and so this solution cannot reasonably explain the discrepancy between the models and the data.At these temperatures, electron densities and Sr ii masses, the NIR triplet feature (which we conclude produces the ∼ 0.7−1.2m P-Cygni feature at early times; see Paper I), is much weaker than the data necessitate.Our 2500 K model does not produce a strong [Sr ii] doublet, but also does not reproduce the NIR triplet.In fact, at this temperature, no prominent features are produced, even though we have included our maximum Sr ii mass ( Sr ii = 10 −3 M ⊙ ).As the models evolve with time, the strong [Sr ii] feature weakens, but still remains discrepant with the observations.
Our results indicate that, for the parameter space we have explored with this set of models, we cannot reproduce the observed feature at ∼ 1.08 m with Sr ii.We would need significantly more Sr ii material to increase the luminosity of the NIR triplet feature, but this Here we present our uniform density material (  = 1) models, which are also offset to match the observations, to aid with comparisons.The different colours correspond to different temperatures (2500, 3500 and 4500 K spectra are red, orange and blue, respectively).Lower panel: Same as the upper panel, but with the synthetic emission spectra computed assuming clumped ejecta material (  = 10 −2 ).
would also lead to an increase in the luminosity of the [Sr ii] doublet (and exceed our Sr ii ejecta mass upper limit).Given there is already a discrepancy between our models and the observations at this wavelength, this suggests that a larger mass would not improve the agreement to the data (as well as being unphysical).The strong [Sr ii] feature, for which we see no observational signature, strongly indicates that these models already possess too much Sr ii material.
If the ejecta of AT2017gfo is well-represented by the range of temperatures and electron densities encompassed by this (i.e. = 1) set of models, then the feature at ∼ 1.08 m cannot be produced by instantaneous Sr ii emission.

Clumped ejecta material ( 𝒇
For the set of models generated assuming  = 10 −2 , we find that the 2500 K model is too cool to produce any observable features, even when we include our maximum Sr ii mass ( Sr ii = 10 −3 M ⊙ ).
For the hotter models (3500 and 4500 K), we find that we need less Sr ii mass overall (relative to the  = 1 models previously presented in Section 4.2) to produce features with luminosities comparable to the observations, since for these models we have higher   , which leads to more highly populated excited levels, which in turn leads to stronger emitted luminosities.Another effect of this increase in   is reflected in the relative strengths of the emission features present in the synthetic spectra we generate.In the previous set of models, the dominant feature was produced by the [Sr ii] doublet.However, for  = 10 −2 , the luminosity of the NIR triplet is comparable to (or exceeds) the luminosity of the [Sr ii] doublet.For the 3500 and 4500 K models, we have emission features that are comparable in strength to the observed ∼ 1.08 m feature.
However, these  = 10 −2 models highlight a concerning discrepancy between the synthetic spectra and observations.Here, the emission feature produced by the NIR triplet is too blue to reproduce the observations.The NIR triplet emission feature peaks at ∼ 10450 Å, whereas the observed emission feature peaks at ∼ 10790 Å (from our analysis of the late-time emission features; see Section 5).This corresponds to an offset of ∼ 340 Å.To reconcile these features, we would require the bulk of the ejecta to be redshifted by ∼ 9800 km s −1 .This velocity is easily reconcilable with typical KN ejecta velocities (since we see evidence for material at velocities ≳ 0.1 ; see Paper I and Section 5), but the notion that the bulk of the ejecta will be redshifted by this amount is unexpected.
Although a bulk offset of ∼ 10000 km s −1 seems unusual, without other identifiable features to independently constrain ejecta properties, it is impossible to conclude if this is indeed redshifted Sr ii emission, or some other effect.An alternative explanation for this discrepancy is that the feature is being produced by a P-Cygni profile, with an absorption feature to the blue that affects the perceived peak wavelength of the emission.However, there exists no evidence for such absorption in the observed spectrum, although we note that this could be concealed if there is some other species 'filling in' the absorption trough with its own emission.Another explanation is that there is some other contributing species redward of the Sr ii NIR triplet, which blends to produce a single emission feature that peaks redward of the weighted peak wavelength of the Sr ii NIR triplet.Other explanations include asymmetrical ejecta structures (which could lead to a disproportionate amount of Sr ii material moving away from us), or observational effects arising from -dependent arrival times from a rapidly-expanding medium.Finally, it is also possible that this feature is in fact not produced by Sr ii at all, although this seems unlikely, given the agreement with Sr ii at early phases (see Watson et al. 2019, Perego et al. 2022 andPaper I).Tarumi et al. (2023) have proposed that the He i  air = 10830 Å line could be the dominant source, which is a much closer match to the line centroid of the observed emission feature at 10790 Å (see Section 4.4 for more details).
The [Sr ii] feature is also present in these models, although not as prominently as in the  = 1 models.These model features again do not agree with the observations.As previously discussed, it may be possible that another species is contributing to the same region of the spectrum, and its absorption is removing any trace of this emission feature.A similar scenario may also explain the discrepancy between the 4500 K model spectra and the observed data at ∼ 4150 Å, which is produced by the Sr ii  air = 4077.7,4215.5 Å resonance lines.Without some mechanism to destroy the emission from these features (both the resonance lines and the [Sr ii] doublet feature for the 4500 K model, or just the [Sr ii] doublet for the 3500 K model), these models cannot reproduce the data.Therefore, we either need to invoke the presence of line-blanketing in the spectrum at wavelengths ≲ 7500 Å, or the ejecta of AT2017gfo do not contain this amount of Sr ii (5 × 10 −4 M ⊙ for the 3500 K model, and 10 −4 M ⊙ for the 4500 K model).
So far, we have not considered how our clumped ejecta models agree with the observations beyond the +7.4 d spectrum.Our good models (3500 and 4500 K) become less luminous with time, as the expanding material leads to lower electron densities.The strong emission feature at ∼ 1.08 m in the observed spectra also decays with time, and the rate at which it fades seems to roughly match the fading of the model features at the same position (assuming our continuum is placed correctly at +8.4 d, which is not clear).The resonance feature and the forbidden doublet remain reasonably persistent across all epochs, and so if it is the case that line-blanketing (or some other effect) is acting to mask their presence in the observed data, it must be the case that it continues to do so across all epochs that we model here (+7.4 − 10.4 d).Alternatively, these features could be removed at later epochs by recombination (or ionisation) effects; i.e. most of the Sr ii recombines to Sr i (or ionises to Sr iii), which would also accelerate the rate of fading of these Sr ii features, improving the agreement to the data.
In Paper I, we presented the inferred masses of Sr ii needed to model the observed sequence of X-shooter spectra of AT2017gfo, up to +7.4 d.Due to the approximations made within tardis, we argued that our models are capable of accurately representing the ejecta properties of AT2017gfo at +0.5 − 4.4 days.The models for the +5.4 and +6.4 d spectra do an adequate job of reproducing the observations, with some caveats, as the observed spectra appear to be transitioning away from a purely photospheric regime.The +7.4 d model is a poor reproduction of the data as, at this phase, tardis is no longer suitable for modelling the ejecta properties.Therefore, we expect our model parameters (including inferred Sr ii mass) to be most reliable up to +4.4 days, be reasonably reliable at +5.4 and +6.4 days, and not reliable at +7.4 days and beyond.In Figure 6, we present the total Sr ii mass in each of our tardis models.The points are qualitatively colour-coded to indicate the reliability of the tardis Sr ii masses presented at each phase (with the green, orange and red points indicating good, medium and poor reliability, respectively).Since tardis only simulates the line-forming region of the ejecta, our tardis masses are actually lower limits for the total Sr ii mass at each phase (hence they are marked as lower limits on the plot; see Paper I for more on these models).Also plotted in Figure 6 are the inferred Sr ii masses from our modelling efforts here, invoking  = 3500 and 4500 K. Clearly, we need more Sr ii material to model the +7.4 d spectrum in this way than we did using tardis for the same epoch.This is not too informative, as the reliability of this tardis model is poor.Note that the total Sr mass included in our tardis model compositions (from +2.4 d onward) does not change, but since the inner boundary is receding with time, we are modelling a greater portion of the total system, and so we infer a larger Sr ii mass.Therefore, we also expect our modelling work here to require more Sr ii mass than our tardis models, which they do.This is consistent with most (or all) of the ejecta becoming optically thin, allowing us to see the effects of all of the Sr ii material in the ejecta.We note that the inferred masses for Sr ii from this analysis are not unphysically high either.Our best-fitting composition profile for the early spectra (+2.4 d onward) of AT2017gfo contains a Sr mass fraction of 0.0325 (see Paper I).Assuming a total ejecta mass for the system,  ≃ 10 −2 M ⊙ , this implies a total Sr mass,  Sr = 3.3 × 10 −4 M ⊙ , which is close to the inferred mass of Sr ii for our 3500 K model (5 × 10 −4 M ⊙ ), and .The reliability of these masses is indicated by marker colour, with the green, orange and red triangles indicating good, medium and poor reliability, respectively.Also plotted are the Sr ii masses inferred from our modelling here, with the points for the  = 3500 and 4500 K models included.We also plot the evolution of the inner boundary velocities from our sequence of tardis models, to indicate how much of the total velocity range of the system our tardis lower limit masses represent at each epoch.
comfortably higher than the Sr ii mass inferred for the 4500 K model (10 −4 M ⊙ ). 10  Watson et al. ( 2019) invoke Sr masses of 1 − 5 × 10 −5 M ⊙ (in the line-forming region) to reproduce the feature between 7000−10000 Å in the +1.4 − 4.4 d spectra, while Perego et al. (2022) suggest that ∼ 8 − 16 × 10 −5 M ⊙ of Sr was ejected by AT2017gfo.Both of these ranges are typically larger than the range of values we present for Sr ii masses from our tardis modelling, which is somewhat expected, since we are comparing Sr masses with Sr ii masses.

Summary and discussion
Although this modelling approach seems quite simplistic, we can deduce some useful constraints for the ejecta of AT2017gfo.First, we find that the identification of the Sr ii NIR triplet transitions as the dominant source of the ∼ 1.08 m feature is problematic for uniform density material, with reasonable estimates for ejecta temperature, velocity and ionisation.This is because the [Sr ii] doublet transitions are the dominant transitions, leading to strong emission at ∼ 6800 Å, for which we see no observational trace (see the upper panel of Figure 5).Therefore, if the ejecta is smooth (i.e. = 1), then the ∼ 1.08 m feature cannot be produced by instantaneous Sr ii emission.
Second, we show that, by invoking significant clumping of the ejecta material, we can reach regimes of high enough electron density such that the NIR triplet emission is stronger (or at least comparably as strong) as the [Sr ii] doublet (see the lower panel of Figure 5).Clumped ejecta material has been invoked to interpret SN ejecta (see e.g.Dessart et al. 2018;Wilk et al. 2020), and so it seems reasonable 10 Note that we have not made any assumptions regarding the ionisation state of Sr at late times, and so this mass comparison between Sr from the early phases, and Sr ii from the present work, although informative, will be somewhat uncertain.A more complete understanding of both non-LTE and non-thermal effects is needed to enable a sensible comparison.
that we invoke it here.If the ∼ 1.08 m feature is produced by Sr ii, then we require clumped ejecta material to boost the relative strength of the NIR triplet feature over the [Sr ii] doublet and the resonance lines.Either the ejecta material was always clumped, or became clumped as it expanded and cooled.
Third, even accounting for clumped ejecta material, which enables our models to produce a feature comparable in strength to observations at ∼ 1.08 m, we still have an issue with feature offset.Our model emission peaks at ∼ 10450 Å, which is ∼ 340 Å blueward of the observed feature peak at +7.4 d, corresponding to a ∼ 9800 km s −1 bulk redshift in the observed feature.The source of this discrepancy may be the result of a number of factors, the most likely of which is our assumption that these late-time spectra originate from a fully nebular regime, which likely does not completely encapsulate the properties of the system.Alternatively, the difference could be evidence of a contribution from some (one or multiple) unidentified species that contributes to (or completely dominates) the emission immediately redward of the Sr ii NIR feature.Other sources for the discrepancy could be some asymmetrical ejecta structure, or observational effects arising from -dependent arrival times from a rapidly-expanding medium.
In light of these issues, it is worth considering the alternative possibility as investigated by Perego et al. (2022) and Tarumi et al. (2023), namely that He i can produce the ∼ 1.08 m feature.As noted above, the He i  air = 10830 Å line closely matches the peak wavelength of the observed feature at late times.There are no other contaminating lines that obviously do not match the data, which could be used to rule out the presence of He i.We note in particular that the suggestion by Tarumi et al. (2023) -that the feature could be dominated by Sr ii at early times, with He i becoming dominant at later times -could alleviate the issues identified above.Specifically, late-time contribution from He i would help explain the apparent wavelength shift, as well as being able to reproduce the 1.08 m feature without the need to invoke a large Sr ii mass that leads to strong [Sr ii] emission (in the  = 1 case).
The case for He i is also strengthened when considering nucleosynthesis arguments.He is among the most abundant elements in many of the composition profiles we present in Paper I. These composition profiles were extracted from nuclear network calculations based on a realistic hydrodynamical simulation of a BNS merger (Goriely et al. 2011(Goriely et al. , 2013(Goriely et al. , 2015;;Bauswein et al. 2013), with some modification.The data have a prescribed distribution of   , from which we extract 'representative' composition profiles, effectively binned by   .This approach was taken to obtain a diverse range of feasible composition profiles for the ejecta from BNS mergers.These composition profiles span a range of   values (0.44 − 0.05), which correspond to compositions dominated by the light iron-group and first -process peak elements (  = 0.44), through to compositions dominated by the heavy lanthanides and third -process peak elements (  = 0.05).Consulting these composition profiles, we see that He makes up 15 per cent of all material (by mass) in our   −0.21b composition profile (see Paper I for more details).This is in agreement with other studies; for example, Just et al. (2023) show BNS mergers are expected to produce significant quantities of He (mass fractions of ∼ 10 per cent).As the analysis by Tarumi et al. (2023) and the calculations presented here show, a full understanding of this (and other) features will depend on further studies, in which full non-thermal and non-LTE excitation and ionisation effects are incorporated for realistic ejecta models.The vertical blue lines mark where we estimate the emission features blend into the continuum.As the position of the continuum is quite uncertain, we do not make predictions for these positions at the blue and red wings of all features.The vertical grey bands correspond to telluric regions.

MODELLING THE LATE PHASE FEATURES
As previously highlighted in Paper I, the evolution of the spectral energy distribution (SED) of AT2017gfo is unlike any previously observed extragalactic transient.After ∼ 1 week, the ejecta has expanded and cooled sufficiently such that the spectra appear to be originating from an optically thin regime.As such, the photospheric approximation used with our tardis modelling in Paper I can no longer be applied to model the spectra.Hence, we need an alternative method to analyse the late-time observations.If the spectra are truly originating from an optically thin region, then the features observed will be produced predominantly by collisionally-excited emission and/or recombination.Here we fit Gaussian profiles to the most prominent emission-like features in the late-time spectra to constrain the peak positions, which we then use in our attempts to identify candidate species for each of the features (see Sections 6 and 7).In addition, measured line widths and centroids can be used to search for evolution of the ejecta properties, and fluxes can be used to constrain the properties of the emitting lines.Note that this approach represents an alternative interpretation for the features at ∼ 1.58 and 2.07 m than that presented in Section 3, and we note that both features become more prominent in the +7.4 − 10.4 d spectra.
We focussed on fitting the most prominent of these emission features -specifically, we fit the ∼ 1.08, 1.40, 1.58 and 2.07 m features.We used a -squared analysis to determine the best-fitting Gaussians for these features, and these are shown in Figure 7.The parameters corresponding to the fits are presented in Table 5.
The 1.08 m feature is only prominent in the +7.4 d spectrum, and so we were only able to fit this feature at this single epoch.This is the feature that we have concluded is likely produced by Sr ii at early times (see Paper I), and late times (see Section 4).
The telluric region at 1.40 m prevents identification of any emission feature at this wavelength in the X-shooter data.However, the +9.4 d spectrum, with the merged HST spectrum covering this wavelength range, clearly shows an emission feature.The presence of this feature in the +9.4 d spectrum should not be completely unexpected, since we can see it is prominent in the earlier +4.9 d HST spectrum (see Figure 1).
The feature at 1.58 m is prominent in the +7.4 − 10.4 d spectra (although by +10.4 d it has weakened significantly, and appears to barely rise above the continuum to the red).The blue wing of this feature in the +9.4 d spectrum contains data from the HST spectrum that has been merged into the X-shooter spectrum here to demonstrate the emission feature at 1.40 m.This join in the two spectra may have introduced some error in this region of the spectrum.For these reasons, we trust the results of our fit to the earlier two epochs at this wavelength range more.Hence, we freely fit the 1.58 m feature in the +7.4 and +8.4 d spectra, and then took an average of these values and refit.We found that the average of these fits was able to broadly reproduce the feature across both epochs, thus motivating no additional analysis; i.e. the data can be adequately modelled with a single emission feature with identical width and peak wavelength across both epochs.We then took the peak wavelength and width from the +7.4 and +8.4 d fits, and imposed those on the +9.4 and +10.4 d spectra, allowing only the line luminosity to vary.For +9.4 d, this fit does not replicate the blue wing of the feature, but this is not surprising, based on the discussion above.For +10.4 d, this fit matches the data reasonably well, although we again note the continuum to the red is higher than the blue, which our Gaussian does not match.
The 2.07 m feature exhibits the most evolution across these four epochs.In the +7.4 d spectrum, it peaks at ∼ 2.06 m, but by +10.4 d, this has shifted by about 400 Å, peaking at ∼ 2.10 m (see Section 2).We find that the evolution of this feature can be easily modelled by the presence of two distinct emission features,11 the relative strengths of which evolve across the multiple observations.By assuming that the feature in the +7.4 d spectrum is composed primarily of a single emission feature, and by fixing the properties of this best-fitting Gaussian for all subsequent epochs, we are able to satisfactorily reproduce the emission feature across all subsequent features by introducing a second Gaussian emission feature redward of the first.We also fixed the properties of this Gaussian across the +8.4 − 10.4 d spectra.The only changes to either of these Gaussians are the strengths of the features -both absolute luminosity from the transitions, and the relative strengths of the features.By transitioning from a regime dominated by the bluer Gaussian, to a blend between the two, with increasing contribution from the redder Gaussian as time evolves, we can reproduce the feature across all four spectra.We find that this feature can be well-reproduced by these two Gaussians, and there is no requirement from the data to have different FWHM velocities, indicating they are both originating from similar regions of the ejecta, which seems to be a physically motivated explanation for their source.These individual Gaussians are plotted in orange and purple in Figure 7, while the co-added resultant Gaussian is plotted in red.
The maximum velocity of the emitting region can be estimated by considering where the emission features blend into the continuum ( edge ).However, this can be difficult to determine accurately due to the uncertainty in continuum location.Nevertheless, we manage to estimate some of these positions (see Table 5 and Figure 7).
For cases when identifying where the feature blends into the con- tinuum is more difficult, we adopt the FWHM velocity as a proxy for the velocity extent of the emitting region, noting that FWHM estimates are argued to be comparable to half-width zero-intensity measurements in a study of nebular phase supernova spectra (Nicholl et al. 2019).Our FWHM velocities are generally comparable to our  edge values, indicating the validity of this approach.From our best-fitting Gaussians, we conclude that a FWHM ejecta velocity,  fwhm = 35600 ± 6600 km s −1 , is capable of reproducing all of the prominent emission features in the late-time spectra of AT2017gfo.Hence, for simplicity in all of our analyses (e.g.Section 4, Appendix C2, and previously in Gillanders et al. 2021), we adopt a FWHM velocity of 0.1  as a characteristic ejecta velocity.We find that all emission features across all late phases (+7.4 − 10.4 days) can be well-matched by invoking a velocity of ∼ 0.1 .We find no evidence for continued feature peak evolution at these late phases (i.e. the peak positions remain constant across the +7.4 − 10.4 d spectra).We also found that the 2.07 m feature evolution is best explained as a blend of two prominent components, which have fixed peak positions, but whose relative strengths vary.Finally, our feature luminosity measurements here can be used to constrain the masses needed to reproduce the observations.These luminosity values enable us to make quantitative statements regarding species we favour or disfavour, based on the validity of their derived mass estimates (see e.g.Section 3 and Appendix C2).

SEARCH FOR POTENTIAL CANDIDATE SPECIES (STRONG TRANSITIONS)
With estimates for the peak position of each emission feature from Sections 2 and 5, we can now attempt to identify the transitions that are responsible for producing these features.Additionally, the luminosity required to produce these features can be used to infer the mass of the emitting species (if we know the intrinsic line strength).
In this section, we consider permitted line transitions (which may produce P-Cygni features as explored in Section 3).In Section 7 we extend the search to weaker candidate transitions.

Method
For our search for permitted transitions coincident with the prominent emission features in the spectra of AT2017gfo, we generate synthetic line lists from the level information that is available for the heavy elements.Specifically, we extracted the level information from the National Institute of Standards and Technology Atomic Spectra Data base (NIST ASD; Kramida et al. 2020) for all neutral, singly and doubly ionised species, with 38 ≤  ≤ 92. 12 Tanaka et al. (2020) find that for KN ejecta at times,  ≳ 1 d, and temperatures,  ≲ 20000 K, the dominant ionisation stages of heavy elements are neutral up to triply ionised (i − iv).Since our tardis models for the third epoch (and all later epochs) have characteristic temperatures much lower than this (3400 K for the +3.4 d best-fitting tardis model), we assume that the dominant ionisation stages will be neutral up to doubly ionised (i − iii).We can validate this choice by analysing the ionisation of the species in our best-fitting tardis models, presented in Paper I. As an example, for each of the lanthanides, the mass fraction of triply ionised material,  X iv ≲ 10 −14 , and so they can be safely excluded from consideration.Therefore, we focus only on extracting the first three ion stages from the NIST ASD.With this level information, we perform a rudimentary set of calculations, where we determine the energy differences between all levels.This gives us the wavelengths for any theoretical transitions between all known levels within the NIST ASD.We note that all energy levels within the NIST ASD are critically evaluated, and so our transition wavelengths that we compute are sufficiently accurate for line identification studies.However, this approach lacks Einstein -coefficients, and so we use a different method to estimate line strength (see below).Since we are considering permitted transitions as the cause of the features, we can make some further cuts on our line list.First, permitted, or electric dipole, transitions (traditionally labelled E1) require a parity change between levels.Therefore, we can discard any lines from our list that do not satisfy this clause.Additionally, E1 transitions can only have variations in the quantum J number, |ΔJ| = 0 or 1, so we excluded all lines from our list that do not satisfy this rule too.For the instances where the configurations and terms for levels can be easily expressed by the LS-coupling formalism, we also eliminate lines that do not obey the LS-coupling rules. 13Finally, we also discard any lines from these lists that have an upper level with a higher excitation energy than the ionisation energy for the species under consideration (also extracted from the NIST ASD).
We lack any measure of intrinsic line strengths from this approach, and so, in an effort to begin to constrain whether the transitions we propose are likely to be prominent, we calculate relative level populations for the upper levels of each of the transitions.We assume LTE, and so the level populations will be governed simply by the Boltzmann equation: where  u is the number of atoms or ions in the upper level,  t is the total number of atoms or ions,  u is the statistical weight of the upper level,  is the LTE partition function,  u is the energy of the upper level, and  is the temperature.We assume that LTE level populations should be a reasonable approximation, at least for the earlier epochs.We cannot compute total numbers of populated levels without an accurate constraint on the ejecta composition and mass, but we can compute relative level populations using the statistical weights and energies of the levels, and some estimate for the ejecta temperature.From Paper I, our tardis models for the +3.4 − 7.4 d spectra have ejecta temperatures ranging from 3400 K at +3.4 d, to 2900 K at +7.4 d.For simplicity, we assume  = 3000 K for the calculations here.This approach allows us to better understand the relative strengths of these transitions, since, if we assume they are all equally strong lines (identical -values), then larger level populations directly correspond to more flux.Hence, we use upper level population as a proxy for line strength, and subsequently filter our data, to identify what we consider to be the most prominent, and therefore most relevant, lines in our analysis.
We normalise these level populations to the upper level that has the strongest potential E1 transition with a wavelength that lies within the range covered by the X-shooter spectra of AT2017gfo (0.5 ≤  ≤ 2.5 m).As an example for how to interpret these relative upper level population values, consider La iii (see Table A1).This ion has two candidate transitions for the 1.40 m feature.Based on the relative level population calculation, the 13898 Å transition comes from an excited level that is ∼ 1.5 times more populous than the excited level that produces the 14100 Å transition.This indicates that, if all other parameters remain identical, the number of photons emitted through the 13898 Å candidate transition will be approximately 1.5 times that emitted through the 14100 Å candidate transition.
With our line list and associated estimates for line strengths, we make comparisons to the observed spectra of AT2017gfo.Species that do not produce any strong emission features coincident with the observed features of interest (∼ 0.79, 1.08, 1.23, 1.40, 1.58 and 2.07 m) are discarded.Additionally, species that produce features coincident with one (or more) emission feature, but also produce comparably prominent emission at wavelengths where we can quantitatively rule out significant flux are also discarded.However, we do not rule out species that only have a few lines in the observable wavelength range of AT2017gfo on this basis.Due to uncertainties associated with our estimates, we retain the cases with a small number of lines, as the contaminant line may be weaker than our analysis predicts, and thus the model could potentially match the observations.When we consider species with many lines with broad flux contamination, we discard the species, since even if the relative strengths are uncertain, we still expect the effect of many blended lines on the spectra to be prominent.
We are interested in the lines that may be producing the observed features, but the locations of these peak positions are somewhat uncertain.Therefore, we adopted an arbitrary wavelength range that corresponds to a 0.05  Doppler shift tolerance for the peak wavelength of each feature, within which lines are considered viable candidate transitions.

Results
We recover two shortlists of candidate species that contain transitions that match the observational data: (i) a viable list of candidate ions that match the data with no perceived limitations, and (ii) a list of ions that may contribute to the data, with some caveat. 14Our list of viable candidate lines for each of the observed features is given in Table A1. 15Our viable list of candidate species include Sr ii, La iii, Ce iii, Gd iii, Ra ii and Ac i (see Table 6).
We find candidate species for all prominent emission features, but some are less plausible than others.Additional points of consideration are the number of total contaminating lines, and the upper level energy. 16Note that where we refer to the number of lines obtained from our analysis (to indicate the numbers of contaminant lines), we are referring to all lines that lie within the wavelength range, 0.5 ≤  ≤ 2.5 m, that have strengths > 10 −3 the strength of the strongest line for that species in our analysis (within the same wavelength range).Our most viable candidate species are presented below.
• Sr ii: In our analysis, this ion yields two lines coincident with the 1.08 m feature, providing a consistency check for our method (see Appendix A2).
14 An example of a caveat is an ion that possesses a line in agreement with an observed feature that is expected to have another comparably strong line at a position where we do not see any observational evidence for emission (see Section 6.1). 15For completeness, we present the wavelengths of the candidate transitions in Table A1 both in vacuum and air. 16Our analysis computes a relative level population, so if only highly excited transitions are available, they will appear prominently in our analysis, despite the fact they originate from energy levels beyond what we expect to be significantly populated in a 3000 K plasma (ignoring any effects that can alter the level populations significantly from LTE).Candidate lines originating from highly-excited levels ( u > 4 eV) are discarded.
• La iii: This species returned only three candidate lines in our wavelength range of interest.Two of these lie in the wavelength range of the 1.40 m feature (13898 and 14100 Å).The third lies in a telluric region (17883 Å), and so we cannot determine whether it is present in the observational data.
• Ce iii: This species returned 30 lines in our analysis, the strongest of which are coincident with many of the observed features.For instance, there is good agreement between our shortlist and the observed features, with the 10723 Å (∼ 1.08 m), 12760 and 12825 Å (∼ 1.23 m), 14659 Å (∼ 1.40 m), 15720, 15852, 15961 and 16133 Å (∼ 1.58 m), and 20691 Å (coincident with the blended ∼ 2.07 m feature) lines all matching the observations.
• Gd iii: Our analysis returned nine lines in our wavelength range of interest.The strongest of these is coincident with the ∼ 1.40 m feature (14337 Å).There are other non-negligible lines that lie at wavelengths ≳ 17000 Å.Some of these are coincident with the blended ∼ 2.07 m feature (20002 and 21265 Å), while there is no observational evidence for the others.Gd iii is a viable candidate for the 1.40 m feature, but to be able to reproduce the blended 2.07 m feature, the relative strengths of the contaminant lines need to be weaker than what we have estimated here.
• Ra ii: There are only three lines extracted by our analysis; two of these are coincident with the 0.79 and 1.08 m features (8021.9 and 10791 Å).
• Ac i: This species returned 27 lines in our analysis.The strongest line is coincident with the 1.40 m feature (13374 Å), although it lies at the extreme blue end of our permitted wavelength range.Ac i also has a strong line coincident with the 0.79 m emission feature (8145.6Å).
Apart from the above six species, the full list of candidate species that may contribute to the data, albeit with some caveats, are summarised and presented in Appendix A1 -these include In i, Ba i, Eu ii, Gd iii, Tb ii, Ho i, Er i−ii, Tm i−ii, Yb i, Lu i, Hf i−ii, Fr i and Th iii (see also Table 6).

Discussion
We present a summary of our shortlisted viable and potential contributing species for each of the observed features in Table 6 (for individual discussion on each shortlisted species, see Section 6.2 and Appendix A1).We also present some verification of our analysis by comparing with previous works in Appendix A2.
We note that many of our shortlisted species that have satisfied our various cuts belong to the lanthanides.This is not too unexpected, since the valence  -shell electron structure leads to many low-lying levels, between which many lines exist, some of which we expect to coincidentally line up with our observed features.From Paper I, lanthanides are expected to make up around 5 per cent (by mass) of the total ejecta material ( ln = 0.05 +0.05 −0.02 ).If the assumption that the observed spectral features are produced by permitted transitions is valid, then we can speculate that our mystery species X (see the analysis in Section 3) is most likely to be Ce iii, but could possibly be Eu ii, Tb ii, Yb i or Th iii.Mystery species Z is also most likely to be Ce iii, but could also be any of Ba i, Gd iii, Yb i, Hf i and Th iii.For both cases, the most likely candidate is Ce iii, a lanthanide.It is also worth noting again that these features may not be made up of emission from a unique species -at the high velocities characteristic for KN ejecta, line blending is also plausible.If this is indeed the case, we would not be able to obtain a unique identification for either mystery species X or Z.
Identification of features produced by multiple ions of the same element are intriguing, as it would lend additional credence to the presence of that element in the ejecta.We find multiple candidate species for the same element: Er i−ii, Tm i−ii and Hf i− ii, although we note that none of these are a convincing viable candidate species for any of the features under investigation.Further analysis of these elements in particular would be worthwhile to see if they do contribute significantly to the spectra -we note that having positive identifications of multiple ions for the same element in the ejecta would be extremely useful for constraining the ionisation.Finally, the specific cases of In i, Ce iii, Eu ii, Gd iii, Ho i, Er i, Yb i, Hf i−ii, Ra ii, Ac i and Th iii are intriguing, since we find each of these species are capable of contributing to multiple observed features.In particular, Ce iii shows remarkable agreement with the observations, as it is a viable candidate for five of the emission features observed at late times (∼ 1.08, 1.23, 1.40, 1.58 and 2.07 m).Not only is Ce iii the most promising candidate, but it is also supported by nucleosynthesis arguments, since it is among the most abundant lanthanides produced in our realistic composition profiles (see Paper I).We note however, that since the lanthanide elements are expected to be co-produced, it may be difficult to disentangle evidence for a single lanthanide species from the contribution expected from all the other lanthanides (although particularly strong lines belonging to the more abundant lanthanide species may be identifiable).
Our rudimentary line search here successfully recovers the candidate La iii and Ce iii lines presented by Domoto et al. (2022) (see Appendix A2).Specifically, they propose that La iii is responsible for absorption at ∼ 12000 Å, and that Ce iii is responsible for absorption at ∼ 14000 Å (these absorption features are linked to the emission features we analyse here, at ∼ 1.40 and 1.58 m, respectively).While the synthetic KN spectral models presented by Domoto et al. (2022) show that Ce iii is capable of producing a prominent absorption feature at ∼ 1.40 m, their models do not produce any other prominent absorption due to Ce iii.The absence of absorption for the other Ce iii transitions we shortlist here as candidates for the other observed features does not rule them out from contributing to the observed emission.Since we are exploring the emission components, we are affected by the upper level populations, whereas the study of absorption features will be dependent on lower level populations.Additionally, the upper level populations can be significantly affected by both non-LTE and non-thermal effects, and so these could act to boost the level populations needed for our shortlisted transitions to produce prominent emission.
We do not recover Y ii as a candidate for the 0.79 m feature, as proposed by Sneppen & Watson (2023).Their models assume a single blackbody fit to the data, and P-Cygni profiles for the strongest Y ii lines.They show that strong Y ii transitions can produce a prominent P-Cygni profile that matches the data at intermediate phases.The presence of strong absorption from Y ii has been noted in other works (e.g.Paper I, Shingles et al. 2023).Those studies predict broad absorption due to Y ii across blue wavelengths, but do not produce a prominent and detectable feature due to the Y ii 4d5p -4d 2 family of transitions (with  air = 7264.2− 7881.9Å), as proposed by Sneppen & Watson (2023).However, the strength of this feature may be dependent on details of the velocity evolution in the models, and further work is required.For our work in this section, we are considering the source of the emission at later times (> 7 days), and so we are focussed on species that can produce prominent emission.In our analysis for Y ii, we recover many contaminant lines that we also expect to be producing prominent emission, for which we see no evidence.Hence, we rule it out as a viable candidate for dominating the 0.79 m emission feature (at +7.4 days and beyond).Pognan et al. (2023) found that the Y ii candidate lines from Sneppen & Watson (2023) are expected to be weak in their KN ejecta simulations, and instead propose Rb i as an alternative candidate for the same 0.79 m P-Cygni feature.
We stress that without estimates for line strengths (i.e.-values) for the lines, it is impossible to definitively determine whether any of our shortlisted lines are the cause of the features we observe in the AT2017gfo spectra.However, this type of analysis does highlight which atoms or ions could be likely contributors, considering our atomic data restrictions.Despite the uncertainties associated with our analysis, we have obtained a list of candidate transitions for the observed features (see Table A1).We present a summary of our shortlisted species for each prominent feature in the observed spectra in Table 6.This table contains a summary of all shortlisted species discussed above, and we argue that this shortlist encapsulates the species that most pertinently need better atomic data.This list can act as a 'priority list' for future atomic data studies.

SEARCH FOR POTENTIAL CANDIDATE SPECIES (WEAK TRANSITIONS)
In Section 5 we showed that the observed features in the late-time spectra of AT2017gfo can be reproduced by pure emission features.We proposed that these may be produced by intrinsically weak transitions originating from some optically thin region of ejecta.Therefore, here we perform a similar search as in Section 6, but this time searching for intrinsically weak transitions that we expect to be prominent.Hotokezaka et al. (2021) and Pognan et al. (2022a,b) suggest that forbidden lines are important for understanding the late phases of KNe, and so here we consider that they may have an effect on the late-phase spectra of AT2017gfo.

Method
Here we follow broadly the same approach as outlined in Section 6, but now consider only forbidden magnetic dipole and electric quadrupole transitions, typically denoted as M1 and E2, respectively.The main steps of the analysis are summarised here, with emphasis on the differences with the approach presented in Section 6.
As in Section 6, we ingest the level information from the NIST ASD for all species with  ≥ 38.We again focus our efforts on the lowest three ion stages of these elements (for the same reason discussed in Section 6).We calculate energy differences between all the levels, which gives us a list of all possible transitions.As highlighted in Section 6, the levels in the NIST ASD (Kramida et al. 2020) have been critically evaluated, and so these computed transition wavelengths will be accurate and reliable for line identification studies.For both M1 and E2 transitions, parity must be conserved.Additionally, only specific variations in the quantum J number are allowed.For M1 transitions, |ΔJ| = 0 or 1, whereas for E2, |ΔJ| = 0, 1 or 2. As before, where the levels have configurations and terms that can be easily expressed by the LS-coupling formalism, we consider the LS-coupling rules.Consulting all of the above, we discard any lines that do not satisfy these rules.Finally, we truncate all lines that have upper levels above the ionisation threshold of the species under investigation.
In the earlier analysis in Section 6, we considered a wavelength range corresponding to a 0.05  Doppler shift, to set our tolerated wavelength range for each of the emission features explored.We deliberately set a broad tolerable threshold for the optically thick case, as we wanted to accommodate the effect that any potential P-Cygni absorption might have on the inferred emission feature peak location.However, this effect is not of concern in this optically thin regime, and so we set a smaller tolerable range (Doppler velocity of 0.02 ).
We again compute the upper level population, and use these values as a proxy for intrinsic line strength (as in Section 6), still assuming  = 3000 K, for convenience.These are normalised to the upper level with the strongest potential M1 or E2 transition within some observable wavelength range (0.5 ≤  ≤ 2.5 m).
In addition to the observed spectra of AT2017gfo presented throughout this work, there are other observations that may help to constrain potential candidate species.Villar et al. (2018) and Kasliwal et al. (2022) present Spitzer Space Telescope (Spitzer) observations of AT2017gfo at +43 and +74 days.There is a significant flux detection at 4.5 m (in the IRAC band, which covers the 3.96 − 5.02 m wavelength range), whereas there is no such signal present at 3.6 m (which samples the 3.18 − 3.96 m wavelength range).We can use the detections from Spitzer to additionally favour or disfavour shortlisted candidate ions, although we note that these observations are much later than the spectra we analyse, and so there may be significant evolution in the ejecta.
We perform a complementary analysis to that presented in this section using the Kurucz atomic data (Kurucz 2017, see Appendix C2).This analysis includes E1, M1 and E2 transitions, and possesses estimates for intrinsic line strengths.However, the Kurucz data base is incomplete for heavy elements.However, utilising this data set can still provide useful results, and as we highlight below, we recover [Y ii] as a viable candidate ion from our analysis with the Kurucz data (see Appendix C2 for details).

Results
Our list of forbidden candidate lines that we extract from our analysis are presented in Table B1. 17 6).Our analysis yields multiple candidate transitions for each of the features identified in the late-time spectra of AT2017gfo.As in Section 6, we discard any candidate species that have highly excited upper levels.We also disfavour species with many lines, since these produce an intolerable amount of flux, for which we have no observational evidence.We summarise our most viable candidate species below.

• [Y ii]:
We recover 24 lines in our analysis.The strongest lines (12990,13355,13711,13731,13961,14372,14539 and 15259 Å) produce a candidate blend, which we expect to produce an emission feature that peaks slightly blueward of the 1.40 m feature (note that many of these lines are the same lines identified in our line identification search using the Kurucz data; see Appendix C2).There are other emission features produced, but at  = 3000 K, these are subdominant.This approach recovers the same strong lines we found in Appendix C2, although we do not obtain the same relative strengths for the features.Note that only one of these strong lines satisfies our tolerable wavelength range for the 1.40 m feature (13961 Å).The recovery of the lines shortlisted from our more detailed analysis in Appendix C2 acts to verify the validity of this rudimentary analysis (see Appendix B2.) • [Rh ii]: Our analysis recovers six lines for this species.Of these, two are coincident with the 1.23 m feature (12248 and 12325 Å).
Table 6.Summary of the candidate species we identify in our analysis for each of the prominent emission features in the late-time spectra of AT2017gfo.The ions we refer to as viable candidate species are those which we propose are capable of reproducing the observations, with no perceived limitations.The ions we propose as potential candidate species are those that may contribute to the observations, but with some caveat (discussed in the main text). Approx.

𝜆
• [Pd iii]: Our analysis recovers seven lines, the strongest of which (21338 Å) is coincident with the 2.135 m feature.The next three strongest lines lie at 9775.2, 14284 and 18039 Å.The 14284 Å line is a potential candidate for the 1.40 m feature, if it is stronger than these two contaminant lines.
• [Ag iii]: Only one line is predicted from our analysis, which is coincident with the 2.135 m feature (21696 Å).
• [Te i]: Our analysis recovers two lines, both of which coincide with the 2.135 m feature (21049 and 21247 Å).
• [Te iii]: Our analysis returns three lines, the strongest of which is coincident with the 2.135 m feature (21050 Å).Another strong line is predicted at 12248 Å, which is coincident with the small emission feature at 1.23 m.
• [I ii]: There are two lines recovered by our analysis.These are coincident with the 1.40 and 1.58 m features (14111 and 15509 Å), although the redder of these is at the extreme blue edge of our wavelength cut.
• [I iii]: We recover five lines from our analysis.The two most prominent lie at 7943.9 and 10641 Å, which are coincident with the 0.79 m and 1.08 m features (see Section 2).
• [Ba ii]: Only two lines recovered by our analysis -the strongest two (17622 and 20518 Å) are coincident with a telluric region and the 2.059 m feature.
• [Tb iii]: Our analysis recovers 56 lines, but the three strongest coincide with the 2.059 and 2.135 m features (20760 and 21121 Å), with the third lying just blueward of the 1.58 m feature (15456 Å).Tb iii has a transition between two levels in the ground configuration that coincides exactly with the 3.6 m Spitzer band non-detection (35655 Å).We predict this to be our strongest line, and so some explanation for not seeing this line at later times (e.g.recombination/ionisation effects, intrinsically weak transition, etc.) is needed.
• [Dy i]: Our analysis returns 140 lines, but we predict most of these are weak, ineffectual transitions.The strongest line within our wavelength range is coincident with the 1.40 m feature (14183 Å).

• [Ir i]:
There are 52 lines returned by our analysis, three of which we predict to be prominent.Two of these strongest lines are coincident with the 1.40 m (14071 Å) and 1.58 m (15813 Å) features, while we have no observational evidence for the third (17287 Å).If this contaminant line is weaker than we estimate, then we have very good agreement with the other two strong lines.
• [Ir ii]: In our analysis, this ion yields 30 lines, with the strongest within our observable wavelength range coinciding with the 2.059 m feature (20886 Å).There is another weaker (but still detectable) line at 12215 Å, which is coincident with the 1.23 m feature.
• [Au i]: This species returned two lines in our analysis, only one of which is expected to be prominent.This line is coincident with the 1.08 m feature (10916 Å), and is the same line as previously identified by Gillanders et al. (2021).
• [Ac i]: Our analysis identifies 50 lines.The strongest line in our analysis lies at 44814 Å, which is exactly coincident with the 4.5 m Spitzer detection.The strongest lines within our observable wavelength range lie at 20837 and 21458 Å, coincident with the 2.059 and 2.135 m features, respectively.
• [Ac ii]: This species has 23 lines in our analysis.The strongest lie at 11004, 37218 and 46310 Å, which correspond to the 1.08 m feature, and the 3.6 and 4.5 m Spitzer bands, respectively.Although we do not see strong evidence for flux in the 3.6 m Spitzer band, we propose [Ac ii] is a viable candidate for the epochs of the X-shooter spectra.
Our candidate species that may contribute emission to the observed spectra, albeit with some caveats, are summarised and presented in Appendix B1 -these include [Pt i] and [Au ii] (see Table 6).

Discussion
We present a summary of our shortlisted viable and potential contributing species for each of the features in Table 6.For individual discussion on each of the shortlisted species, the reader is referred to Section 7.2 and Appendix B1.
For our modelling at early times (see Paper I), we favoured a composition with little material heavier than the lanthanides.However, as we probe later times, and become sensitive to different regions of ejecta, it is possible that the composition deviates from that inferred at earlier times, and so we cannot rule out the possibility of these features being produced by trans-lanthanide species, such as We note that many of our candidate species are lanthanides, as was also the case in Section 6.This is again a result of the lanthanide species having valence  -shell electrons, and so there exist many low-lying levels, between which there are many lines.
We expect some of these candidate species to have coincidental agreement with at least one of the observed features in the spectra.As before, the best way to constrain the presence of a particular species is to search for evidence of it producing (or at least contributing to) a number of features.From our list, [Pd iii], [Te iii], [I ii], [I iii], [Ce iii] (although this is due to a single line that lies within the tolerable wavelength range of both the 2.059 and 2.135 m features), [Tb iii], [Er ii], [Ta i], [Os iii], [Ir i], [Ir ii] and [Ac i] all potentially contribute to more than one emission feature.
Elements for which we shortlist multiple ion stages include Positive identification of any of these multiple ion stages would be extremely useful for constraining the ionisation of the ejecta material.
Both I and Te possess multiple candidate transitions that match the observed emission features, across multiple ion stages.In addition, they are both predicted to be abundant by nucleosynthesis calculations (see Paper I).I (which has atomic mass,  = 127 and 129) corresponds to the top of the second -process peak, and is expected to be quite abundant, 18 and so I having a detectable effect on the observed spectra seems plausible.We find agreement with four of the emission features in the observed spectra (0.79, 1.08, 1.40 and 1.58 m) with two ion stages of I ([I ii] & [I iii]) -see Table 6.
For the specific case of Te, we shortlist both [Te i] and [Te iii] as viable candidates for the 2.135 m feature, as well as finding that [Te ii] is capable of producing a prominent feature coincident with the 4.5 m Spitzer detection.Also, [Te ii] and [Te iii] have lines that coincide with the emission feature at 1.23 m.Te (which has atomic mass,  = 128 and 130) corresponds to the top of the second -process peak, due to the  = 82 magic number.This implies that we can expect Te to be abundant in many different -process scenarios.In fact, Te is among the most abundant elements in many of our composition profiles presented in Paper I (and has a mass fraction of ≈ 10 per cent in our best-fitting   −0.29a composition profile).
We conclude that there are a reasonable number of potential candidate ions for the features in the late phase spectra of AT2017gfo (see Table B1), and we propose that our rudimentary study here has produced a reasonable shortlist of ions that warrant further analysis.As in Section 6, we again present a summary of the candidate species we identify from our analysis in Table 6.This table contains a summary of all shortlisted species discussed above, and we again argue that this shortlist contains the species that most pertinently need better atomic data.
Although our approach here is quite rudimentary, we are performing as detailed an analysis as is currently possible with the existing atomic data.Despite the obvious limitations encountered by not having access to complete line lists, we still manage to identify candidate transitions for each of the observed emission features observed in the spectra of AT2017gfo.While our shortlist is likely incomplete, we propose that it is still useful, as it demonstrates the species that we think are most likely to contribute to the observed data.These species (see Table 6) are ideal candidates for future atomic data studies, which would allow us to build on our conclusions here.Improved atomic data would enable us to quantitatively determine whether our proposed candidate ions contribute to the observations of AT2017gfo.We note that it is important that future atomic data studies not only expand to heavier elements and to longer wavelengths, but also to encompass line information for forbidden transitions, since these are key to interpreting the spectra of kilonovae at late times.

CONCLUSIONS
The spectra of the kilonova AT2017gfo exhibit remarkable evolution.Due to the high expansion velocities and low ejecta mass, the transient evolves extremely rapidly, as is evidenced by the extreme spectral evolution.The spectra quickly reach a quasi-nebular regime, with prominent emission features rising above the continuum.To date, no study has modelled the late-phase spectra in detail, and so here we present an empirical approach to aid interpretation of these spectra.First we show that the spectral feature peaks evolve redward with time.We find that there are seven perceived features in the ejecta of AT2017gfo (at ∼ 0.79, 1.08, 1.23, 1.40, 1.58, 2.059 and 2.135 m), most of which appear prominently across multiple epochs.We have mapped the evolution of these features and performed a line identification study.
We modelled the 1.58 and 2.07 m NIR features at intermediate epochs (+3.4 − 7.4 days) with tardis as single strong permitted transitions produced by some mystery species X and Z.We find that such an interpretation can only achieve agreement with the data at a few epochs (+5.4 and +6.4 d for X, and +5.4 d for Z), and that generally these models produce features that are broader than the observations.This indicates that mystery species X and Z exist at lower velocities than the Sr ii material producing the ∼ 1.08 m P-Cygni feature at these intermediate phases.This could be a result of ionisation (where X and Z are not the dominant ions at higher velocities), or ejecta stratification arising from different ejection mechanisms.
Next we analysed the ∼ 1.08 m feature evolution.The P-Cygni feature at ∼ 0.7 − 1.2 m is thought to be produced by the Sr ii NIR triplet in the first few days.The Sr ii NIR triplet feature plausibly evolves into a pure emission feature at later times (+7.4 d), but the emission component for Sr ii would require a ∼ 9800 km s −1 bulk redshift to match observations.It is possible that this feature is not originating from a fully nebular regime at +7.4 d, or it may be contaminated by another emerging species at redder wavelengths.Our calculations of the emission expected from Sr ii show that if this emission profile is due to the NIR triplet, the [Sr ii] doublet should be detectable.There is no detection of this feature in the observed data, implying the ejecta is highly clumped (  = 0.01).It may also cast doubt on the Sr ii NIR triplet emission persisting to the later epochs.We explore the possibility of He i emission contributing to the data.
While the He i  air = 10830 and 20581 Å lines are remarkably good wavelength matches for the observed 1.08 and 2.059 m features, we conclude that He i emission cannot reproduce the observed relative strengths of these features.While He i could contribute to the 1.08 m feature, it cannot be responsible for the 2.059 m feature.
We next modelled all emission features present in the late-time (+7.4 days onward) spectra as Gaussian-shaped emission lines arising from a nebular regime, with characteristic FWHM ejecta velocity,  fwhm = 35600 ± 6600 km s −1 .We favour this scenario for line formation and perform a search for candidate permitted transitions of ions within the NIST ASD.We find that almost all of our candidate species are lanthanides, and that the 1.58 and 2.07 m features (mystery species X and Z) are both likely to be Ce iii.Although Sr ii is still the favoured explanation for the 1.08 m emission profile, our search returned alternative identifications including Ce iii or Ra ii.While we propose the observed feature at 0.79 m likely corresponds to the point at which the Sr ii NIR triplet P-Cygni rejoins the continuum, our search returned Ra ii and Ac i as alternative explanations.We also shortlist Ce iii as the species most likely responsible for the 1.23 m feature.The shortlisted species responsible for the ∼ 1.40 m feature include La iii, Ce iii, Gd iii and Ac i.We note it is likely that lineblending is commonplace in KN ejecta -even at later times -and so it is possible that none of these emission features have unique and identifiable contributions from a single species.
An alternative interpretation of the emission features in the spectra ≳ 7 days is that they are produced by intrinsically weak lines in optically thin emitting material, and so we perform a similar search for possible forbidden transitions.Candidate species include neutral, singly and doubly ionised ions from the light -process elements (Y -Ba), lanthanides (Tb & Dy), and very heavy trans-lanthanide elements (Ir -Ac).We propose that Te and I are the most promising elements that warrant further investigation, as they possess multiple candidate transitions that match the observed emission features, across multiple ion stages, in addition to lying close to one of the -process abundance peaks.However, with the existing atomic data available, we cannot make unique and definitive identifications of forbidden lines for the observed emission profiles.
The major limitation throughout this work is a lack of atomic data.We lack estimates for intrinsic line strengths for many of our shortlisted species for the observed features, which prevents a more in-depth analysis.Not only does the community need access to an accurate and experimentally-calibrated atomic line list, but it needs atomic data that also encompasses the NIR (and MIR), that also contains intrinsically weak (i.e.forbidden) transitions.Having these would make line modelling and identification studies much more accessible.Finally, although line blending makes unique identifications difficult, future observations extending into the MIR will be invaluable, since the density of transitions decreases with increasing wavelength, making individual line identifications more feasible.

APPENDIX A: SEARCH FOR POTENTIAL CANDIDATE SPECIES (STRONG TRANSITIONS) A1 Shortlisted species
Our candidate species that may contribute to the data, albeit with some caveats, are summarised and presented below (see Table A1).
• In i: Our analysis recovers 15 lines.Two of the strongest lines are coincident with the 1.23 and 1.40 m features (12916 and 13434 Å).Although the upper level energies are quite high for these transitions, they do lie below the upper level energy limit imposed for this analysis.
• Ba i: We obtain 35 lines for this species in our analysis.The three strongest lines lie at 22318, 23260 and 25522 Å.The bluest of these lies at the red end of our tolerable wavelength range for the 2.135 m feature.Although this line is likely too red to be the main source of the flux for the 2.135 m feature, it is capable of contributing to the data (although we would need the two contaminating lines to be weaker than what we have estimated).
• Eu ii: Our analysis returns 81 lines, although the majority of these are clustered at the blue end of the observations (< 8000 Å).Of the redder lines in our analysis, many are coincident with the 1.08 m (10312, 10739 and 10907 Å), 1.40 m (13611, 13882 and 14150 Å) and 1.58 m (15075 and 15322 Å) features.The only contaminant lines at wavelengths > 8000 Å lie immediately blueward of the 1.08 m feature, which produce a candidate blend for a feature at ∼ 1 m, for which we see no evidence.We need line blanketing for wavelengths < 8000 Å, as well as the contaminant lines at ∼ 1 m to be relatively weaker than what we have estimated here for Eu ii to be a viable candidate.
• Gd iii: See the main text (Section 6.2).
• Tb ii: There are 162 lines in our wavelength range of interest.Some of the strongest lines are coincident with the 1.58 m feature (15967, 16068 and 16204 Å).However, there are multiple weaker (but non-negligible) lines in our observed wavelength range, for which we see no evidence.
• Ho i: Our analysis recovered 262 lines for this species.Of these, many are expected to be weak.Many of the most prominent transitions are coincident with the observed features; specifically, the 0.79 m (8100.7 and 8104.4Å), 1.08 m (10932 Å), 1.23 m (11935 and 11866 Å) and 1.40 m (14441 and 14452 Å) features.The relative strengths of these lines do not match the observations, and so we need these strongest transitions to have different relative strengths for Ho i to be considered a viable species.
• Er i: This ion returns 161 lines in our analysis, although many of these are negligibly weak.The strongest two lines are coincident with the 1.23 and 1.40 m features (12992 and 13934 Å).Er i also has a small cluster of weak, but still detectable lines at ∼ 8500 Å, which we see no strong evidence for, and so we need these lines to be weaker than what we estimate (or we need to invoke line-blanketing) for Er i to be a viable candidate.
• Er ii: There are 96 lines in our wavelength range of interest.Er ii has its most prominent line coincident with the 1.40 m feature, although we note that it lies at the very edge of our permitted wavelength range for this feature (14653 Å).There are many weaker (but non-negligible) lines between ∼ 7500 − 10000 Å, which may blend to produce some detectable flux that we do not see evidence for in the data.
• Tm i: Our analysis returned 242 lines for this species.Of these, the strongest transition lies at 7622.2 Å, coincident with the 0.79 m feature.We need estimates for the strengths of the many other lines in our analysis to verify that they are much weaker than this prominent line, for Tm i to be a viable species.
• Tm ii: This species returned 168 lines in our analysis.We find very strong transitions in the near-UV and blue end of the visible spectrum ( ≤ 6000 Å), which corresponds with the heavily lineblanketed region of the AT2017gfo spectra.Two of its strongest lines at redder wavelengths are coincident with the 1.08 m feature (10714 and 11090 Å).
• Yb i: Our analysis returned 60 lines for this species.Many of the strongest lines are coincident with the 1.40 m (13888 Å), 1.58 m (15391 Å) and blended 2.07 m feature (19835, 20926, 22202 and 22276 Å).We note the presence of two comparably strong lines, which lie at 17984 and 14793 Å -these lines are coincident with a telluric region, and the gap between the 1.40 and 1.58 m features, respectively.While telluric absorption can explain away the contaminant 17984 Å line, we see no evidence for the 14793 Å line in the observed data.We need the relative strength of this line to be much weaker than our estimate for Yb i to be a viable candidate.
• Lu i: Our analysis shortlists only three transitions (13375, 18240 and 24177 Å).The strongest is coincident with the 1.40 m feature (although it lies at the extreme blue end of permitted wavelengths), one lies in a telluric region, and the third lies just beyond the wavelengths covered by our X-shooter spectra.The 13375 Å line is too blue to be the sole producer of the 1.40 m feature, but it may contribute.
• Hf i: Our analysis finds 68 lines within our wavelength range of interest.Many of these lines lie at the blue end of the spectrum (≲ 7500 Å), where the data may be significantly impacted by line blanketing effects.Two of the strongest lines recovered are coincident with the blended 2.07 m feature (19865 and 20533 Å).Other weaker, but still prominent lines are coincident with the 0.79 m feature (8206.8 and 8279.2Å).There are other detectable lines across the spectral range, with a clustering around ∼ 1 m, which are too blue to reproduce the 1.08 m feature.Some explanation for not seeing evidence for prominent flux from Hf i at wavelengths ≲ 7500 Å is needed before Hf i can be considered a viable candidate.
• Hf ii: Our analysis returns 63 lines for this ion, many of which lie at wavelengths ≲ 7000 Å.There are a small number of strong lines at redder wavelengths -10787 and 10904 Å (which are coincident with the 1.08 m feature), 12604 and 12875 Å (which are coincident with the emission feature at ∼ 1.23 m), and a small grouping of lines at ∼ 9500 Å, for which we see no evidence in the observations.
• Fr i: Our analysis recovered 19 lines, two of which are expected to be much more prominent than the rest (7181.8 and 8171.7 Å).The redder of these two lines is coincident with the 0.79 m feature, while there is no observational evidence for the bluer line -if this bluer transition is relatively weaker than we predict here, then Fr i can be a viable candidate for the 0.79 m feature.
• Th iii: We obtained 63 lines from our analysis.The strongest lines generally agree with the observed peak positions, with the 13446 Å (1.40 m), 16064 Å (1.58 𝜇m), 19948, 20011, 20306, 20993 and 21510 Å (2.07 m) lines all broadly agreeing with observation.However, there are a number of weaker (but non-negligible) linesthe most prominent of these lie at 12726, 14767 and 14781 Å.We see no evidence for any of these three lines in the observed data.

A2 Verification of analysis
Here we present some discussion to highlight the validity of our rudimentary analysis.By comparing to previous results, we can verify whether our analysis is returning sensible results.
As a sanity check, and proof of the ability to recover sensible candidate transitions from this analysis, we note that we recover two candidate transitions from Sr ii for the 1.08 m feature.These are two of the three Sr ii lines that are responsible for producing this feature at early times (i.e. the Sr ii NIR triplet; see Paper I).It is worth noting however that only the reddest two of these three transitions satisfies our wavelength cut, indicating again that this feature does indeed lie redward of the weighted average of the expected centroid of the Sr ii NIR triplet blend (see Section 4).
As discussed in Section 3.2, Domoto et al. (2022) focus on identifying the NIR features present in the same +1.4−3.4 d X-shooter spectra of AT2017gfo as we analyse here.They assume the ejecta is optically thick, and attempt to identify the cause of the inferred absorption features that have troughs at ∼ 12000 and 14000 Å. Domoto et al. (2022) propose La iii as the cause of absorption at ∼ 1.20 m, and Ce iii for the absorption at ∼ 1.40 m.Specifically, they identify two strong transitions for La iii (with log [  ] > −3) at 13898 and 14100 Å, as the cause of the observed absorption trough at ∼ 1.20 m.These are the same transitions as we find in our search for candidates for the ∼ 1.40 m emission feature (see Table A1).Hence it is plausible that these strong La iii lines cause both shallow absorption in the early spectra and evolve into pure emission in the later spectra.Domoto et al. (2022) further identify 38 strong transitions for Ce iii at NIR wavelengths, and find that the strongest of these could produce a P-Cygni feature in the first ∼ 3 days, with emission peaking around 1.59 m.We discussed this in Section 3.2 and could not quantitatively confirm that our synthetic, photospheric phase spectra matched the observations for these lines.However, these strong, permitted Ce iii lines could potentially be capable of reproducing the 1.58 m emission feature that develops later.In our line search, we recover many of the same lines as Domoto et al. (2022).As well as being able to match the emission feature at ∼ 1.58 m, many of our strongest Ce iii lines also broadly match the emission features at ∼ 1.08, 1.23, 1.40 and 2.07 m.
Our agreement with many of the candidate lines proposed by Domoto et al. (2022), belonging to both La iii and Ce iii, demonstrates that we are able to recover sensible shortlists of candidate transitions with our approach here.nent line we find in Gillanders et al. (2021).For the [Au i] case, we recover the same three most prominent transitions shortlisted in Gillanders et al. (2021), albeit with different relative strengths (due to the effects of consulting -values).We note that for the [Pt i] case, we fail to recover some of the weaker (but non-negligible) lines that we shortlist in Gillanders et al. (2021) in this analysis.This is due to our analysis discarding these lines as possibilities, due to them not satisfying our imposed LS-coupling rules.It is known that the LS-coupling formalism does not apply well to heavy ions, and so this result is somewhat expected.This effect will be more pronounced for heavier ions, and so we propose the possibility of missing transitions that are relevant to our inferences in this work is smaller for the lighter ions.
While expanding our analysis to include the intermediate coupling rules as well as the LS-coupling rules would assist our recovery of the missing [Pt i] lines, we opted not to pursue a study of these additional lines, for the following reason.Throughout our analysis, we have found that there is no convincing candidate ion or element that can explain all (or many) of the observed features.We are already finding many possible candidate ions (even after enforcing the strict LS-coupling rules), and so adding additional lines to this analysis does nothing to improve our ability to shortlist a small number of candidate ions -instead, it will likely only increase the numbers of possible line identifications. 19If all of the spectral features arise from different species, or if they arise from blended emission from multiple lines, then we really cannot advance with this type of line identification study without possessing Einstein -values.
Although this comparison highlights the uncertainty of our approach, it also demonstrates that level population as a proxy for line strength is a useful approximation, since we can use this approach to find the strongest few lines -although the relative strengths of the strongest lines appear to be uncertain relative to one another (which is something we have considered while compiling our shortlist of viable candidates).

APPENDIX C: SEARCH FOR POTENTIAL CANDIDATE SPECIES (KURUCZ)
The Kurucz data base (Kurucz 2017), from which we exported atomic data to include in our tardis modelling (see Paper I), contains atomic data for most of the lowest few ionisation stages of species between  = 38 − 46, as well as data for  = 56.Specifically, it contains line lists that include forbidden M1 and E2 transitions in addition to the permitted E1 transitions that we needed for our earlier tardis modelling.The Kurucz line lists contain a mix of theoretical and experimentally measured lines, but we note that our analysis favours the transitions between low-lying levels, which tend to be the ones that have been experimentally verified.With these line lists, we search for candidate transitions for each of the emission features present in the late-time spectra of AT2017gfo.
Note that the Kurucz atomic data base presents all transition wave-spectra of AT2017gfo.We also present the lines that are coincident with the 3.6 and 4.5 m Spitzer bands -see Table C1.
From our analysis of the species in the Kurucz data base, we find that [Y ii], [Zr ii], Nb i, Mo i, [Rh ii] and [Rh iii] all produce features coincident with the observed data.The lines that are prominent in our models, and potentially produce the observed emission features are presented in Table C1.Details of each of these candidate species, and their viability, are discussed below.We note that for all candidate ions presented here, the inferred masses needed to match the observed feature strengths are typically higher than those supported by nucleosynthesis arguments.22 • [Y ii]: This ion produces two emission features, centred at ∼ 7800 and 13750 Å.These are coincident with the emission features we observe at ∼ 0.79 m (the point at which we proposed the Sr ii NIR triplet absorption rejoins the continuum; see Section 2), and ∼ 1.40 m.It is possible that the ∼ 0.79 m region of the spectra has contribution from some emission feature, at least at late times.One issue with Y ii as the main source of flux for this 1.40 m feature is the larger than expected mass needed to power the observed emission feature.We require model masses,  Yii ∼ 8.0 × 10 −2 , 9.2 × 10 −3 , and 4.2 × 10 −3 M ⊙ for our 2000, 3500 and 5000 K models, respectively, to reproduce the measured flux of the observed 1.40 m feature at +9.4 d (see Section 5).Note that higher temperature lowers the required mass, but increases the relative strength of the 7800 Å feature (which is not supported by the data, since the 1.40 m feature is more prominent than the 0.79 m feature).For Y ii to be producing the ∼ 1.40 m feature, we need even higher  (or some alternative non-thermal mechanism to populate the relevant levels).
• [Zr ii]: This species produces an emission feature at 10800 Å, exactly coincident with one of our emission features.However, there is significant contribution from many weaker, blended lines.This emission feature becomes more prominent (relative to these other weaker lines) at higher .We require large masses to reproduce the flux needed to power the 1.08 m feature at +7.4 d ( Zr ii ∼ 3.3 × 10 −2 M ⊙ at 5000 K).
• Nb i: This ion produces a broad emission feature at 10720 Å, coincident with the 1.08 m feature.However, it produces comparably strong emission features at 7800 Å (coincident with the 0.79 m feature), and 9160 and 16760 Å, which do not match the observed spectra.Variations in temperature across the range explored does not improve the relative strengths, and so we disfavour this ion as a prominent source of emission in the AT2017gfo spectra.This candidate ion has (intrinsically weak) E1 lines as the source of its most prominent emission in this regime.We note the presence of mid-infrared (MIR) lines that are coincident with the 3.6 and 4.5 m Spitzer bands (35636 and 46410 Å; see Table C1).We find the 3.6 m line is stronger than the 4.5 m line across all model temperatures, and so it seems unlikely there will be much Nb i present (at least at the times of the Spitzer observations).
• Mo i: This species produces a prominent feature (from E1 transitions) at 20800 Å, coincident with our blended 2.07 m feature.This species also produces strong features at 5470 and 6700 Å, due to M1 and E2 transitions, respectively.These features are predicted to be stronger than our favourable 20800 Å feature.We do not see either of these features in the spectra, although they lie at wavelengths that suffer from line-blanketing (at least at early times).We need  Mo i ∼ 6.1 × 10 −2 M ⊙ (at  = 5000 K) to match the flux producing the 2.07 m feature at +7.4 d.
• [Rh ii]: This species produces a number of broad emission features, centred at ∼ 6730, 8020, 10820, 17350, 21810 and 41650 Å.The emission feature at 1.08 m is dominant at higher temperatures, and there is reasonable agreement between the model emission features (at ∼ 8020 and 21810 Å) and the observed 0.79 and 2.07 m features.The features at ∼ 6730 and 17350 Å are also expected to be reasonably strong, but we see no evidence for them in the observed spectra.Therefore, we ultimately disfavour this ion as a significant contributor to the emission components observed in the late-phase spectra of AT2017gfo.We note the presence of the MIR feature at 41650 Å, which lies within the wavelength range of the 4.5 m Spitzer band, making it potentially detectable if present in significant quantities.We need masses,  Rh ii ∼ 6.0 × 10 −2 and 0.13 M ⊙ (at  = 5000 K), to produce the flux needed to power the +7.4 d 1.08 m and +8.4 d 2.135 m features, respectively.
• [Rh iii]: This ion produces a prominent M1 feature at 7690 Å, which is semi-coincident with the observed ∼ 0.79 m emission feature.We also find a strong emission feature at 46550 Å, which is coincident with the 4.5 m Spitzer band, making this species detectable in the MIR.

C3 Discussion
From this analysis, we find no convincing viable candidate species for any of the observed emission features, assuming they originate from an optically thin medium that is dominated by intrinsically weak lines (see Table C1 for a summary of all candidate transitions).From this analysis, the only possible identification comes from [Y ii], although this ion has the issue of requiring more mass than is plausible for the AT2017gfo ejecta.In fact, we find that for this approach, we need more mass for all potentially viable candidates than we expect is reasonable, based on our composition profiles from modelling the early phases, and the inferred ejecta mass for AT2017gfo (see Paper I).The fact that the inferred masses for many of the ions discussed above are higher than what we can tolerate in the ejecta of AT2017gfo is perhaps indicative of some error with our approach here; i.e. perhaps LTE level populations are not sufficiently accurate to reproduce reasonable estimates for the emission of the relevant lines.Alternatively, it may be evidence of intrinsically weak lines not being capable of reaching the observed luminosities of the features.If this is the case, then it indicates that intrinsically weak transitions cannot be responsible for the observed emission features in the AT2017gfo spectra at late times.This would suggest the correct interpretation is that these features are produced by strong transitions (e.g. the lines presented in Section 6).To quantify this issue, a more accurate calculation including non-thermal and non-LTE effects is necessary (such a study is not currently possible with the available data).
Ultimately, we lack extensive atomic data to explore the species heavier than  = 46, and so we cannot constrain any more viable candidate ions than [Y ii] from this analysis.We need to either rule out other species (for which we do not yet have complete atomic data), or identify other features that are also possibly produced by [Y ii], which would help to corroborate its presence in the ejecta (although then we would have to reconcile the required mass with our expected composition).
Table A1.Candidate E1 transitions for the ∼ 0.79, 1. 08, 1.23, 1.40, 1.58 and 2.07 m features in the spectra of AT2017gfo.Only transitions with a relative level population ≥ 0.1 are included, for brevity.For each feature, two sets of transitions are given, separated by a vertical space.The upper set of transitions are those that belong to our most viable transitions, while the lower set are the potential candidate transitions.The species, wavelength (both in vacuum and air, for completeness), level energies, quantum J numbers, and the relative level population for the upper level are presented for each transition.The configurations and terms are displayed (from the NIST ASD; Kramida et al. 2020)

Figure 1 .
Figure1.Sequence of AT2017gfo spectra, scaled and offset for clarity.The early Magellan spectrum is plotted (green), as are the spectral sequences obtained with X-shooter (black) and HST (blue).Each observation has been annotated with its phase, relative to the gravitational wave (GW) trigger.The peaks of the most prominent emission-like features have been estimated, and joined (red lines) across epochs, to highlight their evolution.The vertical grey bands correspond to telluric regions.

Figure 3 .
Figure 3. Sequence of AT2017gfo spectra, scaled and offset for clarity.We plot the early Magellan spectrum (green), as well as the sequences of X-shooter (black) and HST (blue) spectra.All observations are labelled with their phase, relative to the GW trigger.We also plot shaded regions corresponding to the locations where we should see any absorption components of P-Cygni features, with the intensity of the shading increasing with blueshift.The purple, brown, green and orange shaded regions correspond to the locations of proposed absorption for the ∼ 1.08, 1.23, 1.58 and 2.07 m features, respectively.The proposed rest wavelengths for these features are also plotted, as vertical dashed red lines.The vertical grey bands correspond to telluric regions.

Figure 4 .
Figure 4. Sequence of AT2017gfo spectra, from +3.4 − 7.4 d, plotted alongside the continua from our best-fitting tardis models (presented in Paper I;red).Also plotted are the best-fitting tardis models with our mystery species X (blue dashed lines) and Z (green dash-dotted lines) included.All spectra have been plotted in scaled flux, and are offset for clarity (2.4,1.8, 1.2 and 0.6 for the third, fourth, fifth and sixth epoch spectra, respectively).The models with our mystery species X and Z have in some instances an additional additive offset applied, to better match the local continua around the observed features (see Table2for these model offsets).Telluric regions have been marked by grey shaded regions, and the observed flux in these regions are uncertain.The epochs have been annotated with their phase relative to GW detection.

Figure 6 .
Figure6.Total Sr ii masses inferred from our tardis modelling in Paper I (lower limits; triangles).The reliability of these masses is indicated by marker colour, with the green, orange and red triangles indicating good, medium and poor reliability, respectively.Also plotted are the Sr ii masses inferred from our modelling here, with the points for the  = 3500 and 4500 K models included.We also plot the evolution of the inner boundary velocities from our sequence of tardis models, to indicate how much of the total velocity range of the system our tardis lower limit masses represent at each epoch.

Figure 7 .
Figure7.Best-fitting Gaussians overlaid on the late-time spectra of AT2017gfo.The observed spectra are offset for clarity (4.5 × 10 36 , 3.0 × 10 36 and 1.5 × 10 36 erg s −1 Å −1 , for +7.4,+8.4 and +9.4 d, respectively), and are annotated with their phase, relative to the GW trigger.The solid red lines represent our best-fitting Gaussians for the most prominent emission features.The feature at 2.07 m in the +8.4 − 10.4 d spectra can be reproduced well by the composite of two Gaussians, which are plotted in orange and purple, while the co-added resultant Gaussian is plotted in red.The vertical blue lines mark where we estimate the emission features blend into the continuum.As the position of the continuum is quite uncertain, we do not make predictions for these positions at the blue and red wings of all features.The vertical grey bands correspond to telluric regions.

Table 1 .
Approximate peak wavelengths of the prominent emission-like features in the observed X-shooter spectra of AT2017gfo.The peak locations of the features shift over time, indicating feature evolution.Although this feature is only detectable in the HST spectrum, we include it here since this HST spectrum was merged with the X-shooter spectrum taken at the same epoch.Note.The emission feature with an approximate peak position of 1.40 m peaks at 1.38 m in the +4.9 d HST spectrum (see Figure1).Ratio of the peak positions of the observed emission features and their final peak position.This is a visualisation of the data in Table1, and helps illustrate the redward-shifting nature of the observed emission features with time.The 1.40 m data point at +4.9 d is included, although it does not appear in the main part of Table1(but is clearly visible in the HST spectrum presented in Figure *

Table 2 .
tardis mass fractions (MF tardis ) of mystery species X and Z needed to reproduce the features in the observed spectra of AT2017gfo.These offset values correspond to the additive offsets that have been applied to our models to better match the local continua around the observed features (and are quoted in scaled flux).

Table 5 .
Parameters of the best-fitting Gaussians presented in this work.
* These  edge values possess some associated error from spectral splicing, and so we expect them to be very uncertain. peak(m)

Table A1 .
Birch & Downs 1994;Morton 2000;Ryabchikova et al. 2015)eme.(continued)Note.airvalues have been computed assuming the standard vacuum-to-air conversion from VALD3 (seeBirch & Downs 1994;Morton 2000;Ryabchikova et al. 2015).This is a candidate transition for the bluer component of the feature at ∼ 2.07 m (2.059 m; see Section 5).b This is a candidate transition for the redder component of the feature at ∼ 2.07 m (2.135 m; see Section 5).
o Denotes an odd parity.a