© The Authors 2024

1. Introduction

The detailed spectra of the gravitational-wave detected kilonova GW 170818/AT2017gfo (Abbott et al. 2017; Coulter et al. 2017) revealed the first spectroscopic identification of freshly synthesised r-process material in a binary neutron star (BNS) merger Sr II (Watson et al. 2019). Several r-process line identifications have since been proposed in addition to Sr II: Y II (Sneppen & Watson 2023), Te III (Hotokezaka et al. 2023), La III, and Ce III (Domoto et al. 2022). The appearance time of these features can provide strong tests for line identifications because the rapidly cooling ejecta undergoes quick transitions between ion states, at least under local thermal equilibrium (LTE) conditions (Sneppen et al. 2023a). Tracking the subsequent evolution of the different spectral features also provides strong constraints for ejecta models, including information on geometry and spatial abundance structures, potentially allowing us to ‘dissect’ the kilonova. Spectra obtained with an adequate temporal cadence can constrain different parts of the evolution at various times and wavelength ranges (e.g. Shappee et al. 2017; Andreoni et al. 2017; Nicholl et al. 2017; Chornock et al. 2017; McCully et al. 2017; Pian et al. 2017; Smartt et al. 2017; Tanvir et al. 2017).

To examine both the emergence and evolution of features requires a detailed early- to late-stage spectral series with the broadest possible wavelength coverage, which no individual dataset contains. In this paper, we therefore track the spectral evolution of AT2017gfo by compiling and analysing the early spectra prior to 1.4 days (Shappee et al. 2017; Andreoni et al. 2017; Buckley et al. 2018), the subsequent daily ultraviolet (UV) to near-infrared (NIR) coverage using the VLT/X-shooter spectrograph (Pian et al. 2017; Smartt et al. 2017), as well as the space-based coverage from the Hubble Space Telescope (HST), which provide the only reliable data in the NIR telluric-dominated regions of this event (Tanvir et al. 2017). The combination of these datasets helps elucidate the physics of the spectral features in a way that no individual dataset has been able to before.

We present a summary of the datasets and consider the complementary information they contain in Sects. 2 and 3. We analyse this compilation of all early published spectra of AT2017gfo in a consistent way to provide the evolution of the blackbody temperature (Sect. 4) and a description that is as accurate as possible of the first appearance of the line features (Sect. 5). We analyse the emergence time of the 760 nm P Cygni feature and its identification with Y II in Sect. 5.1. Similarly, we analyse the emergence of the strongest feature (at ∼1 μm) and verify that it is consistent with the Sr II interpretation in Sect. 5.2, and with a possible proposed association with He I in a companion paper (Sneppen et al. in preparation). The NIR spectral features are also analysed, where we show that the 1.4 μm feature in particular displays an apparent P Cygni-like spectral shape. We discuss the possible different origin of the NIR features compared to the optical lines (Sects. 5.3–5.5) and how this relates to the complementary nature of their velocity structure.

2. Overview of the data

On 2017 August 17 at 12:41:04 UTC, the merging neutron star gravitational wave event GW170817 was detected. It was accompanied 1.74 seconds later by a short gamma-ray burst, GRB170817A (Abbott et al. 2017). Following these triggers, a multitude of international observing campaigns searched for the optical counterpart, which was ultimately found in the host galaxy NGC 4993 at 0.45 days post-merger (Coulter et al. 2017). Over the subsequent hours, days, and weeks, a series of detailed spectra recorded the evolution of this rapidly evolving optical transient. In Fig. 1, we show the spectra selected for this analysis because of their 1) unique temporal coverage and then 2) their high signal-to-noise ratios and broad wavelength coverage. The times noted are the time post-merger in the rest frame of the host galaxy NGC 4993, z_cosmic ≈ 0.01 (Hjorth et al. 2017; Howlett & Davis 2020; Nicolaou et al. 2020).

Fig. 1. Spectral series of the kilonova AT2017gfo showing the evolution over the first 10 days post-merger. This spectral compilation of data includes the early spectra taken with the Magellan telescopes (Shappee et al. 2017), ANU 2.3 m telescope (Andreoni et al. 2017) and SALT (Buckley et al. 2018). From 1.4 days post-merger, the X-shooter spectrograph at the European Southern Observatory’s VLT detailed the UV to NIR structure with a daily cadence (Pian et al. 2017; Smartt et al. 2017). To constrain the rapid evolution at early times, we also include the 1.5 day Gemini FLAMINGOS-2 and SOAR spectra (Nicholl et al. 2017; Chornock et al. 2017). Towards intermediate and late-times, HST (Tanvir et al. 2017) infrared spectroscopy additionally constrains the spectral structure in the telluric region around 1400 nm. The spectral continuum is well described as a blackbody in all photospheric epochs and even into nebular-phase epochs with a blackbody temperature evolution that is detailed in Sect. 4. In the zoom-in sub-panels, individual spectral lines are highlighted at different epochs, including (1) the 760 nm (2) 1.0 μm, and (3) 1.4 μm P Cygni features, and (4) the NIR line-emission features.

0.48–0.53 days post-merger. The first two spectra of AT2017gfo were obtained at 0.48 and 0.53 days post-merger by the Magellan telescopes (with the Low Dispersion Survey Spectrograph, LDSS-3, and MagE spectrographs) at Las Campanas Observatory, Chile (Shappee et al. 2017). These spectra display a featureless blue continuum that is well modelled as a blackbody. The spectral coverage (LDSS-3: 425–1000 nm, MagE: 360–700 nm) does not include the blackbody peak at this early epoch, so that the observed temperature (inferred from the spectral slope) is not as well constrained as in later epochs, but is likely around 10 000 K (Shappee et al. 2017; Waxman et al. 2018). AT2017gfo faded measurably in the blue in the short interval covered between these initial spectra, which highlights the rapid cooling of the transient.

0.92 days post-merger. The third independent spectrum of AT2017gfo was taken at 0.92 days with the Australian National University 2.3 m telescope at Siding Spring Observatory (ANU, Andreoni et al. 2017). The spectrum still appears largely featureless at this epoch. The optical wavelength coverage (320–980 nm) now nearly covers the observed blackbody peak, indicating a significantly cooler transient, T ≃ 6000–7000 K. We note that the flux calibration for the SED is somewhat uncertain, and the exact blackbody temperature inferred from the spectral slope is therefore not very tightly constrained. However, the drastic decrease in temperature in the short time span since the earliest spectra is supported by photometric observations (e.g. Drout et al. 2017; Evans et al. 2017; Villar et al. 2017), where interpolating the temperature at 0.92 days would suggest T = 6800 ± 100 K.

1.17 days post-merger. The fourth spectrum in our series of AT2017gfo was taken at 1.17 days with the Southern African Large Telescope (SALT, Buckley et al. 2018). The optical wavelength coverage of SALT (380–900 nm) reveals the first epoch, where a pure blackbody continuum model poorly reproduces the data. Instead, sizeable absorption is seen redward of 650 nm, which we argue in Sect. 5.2 is the first appearance of the 1 μm P Cygni feature. Due to the importance of the precise emergence time of the ∼1 μm feature, we re-reduced the SALT spectrum using different standard stars in order to confirm that this observed feature is consistent and robustly constrained. We describe this re-reduction in Appendix A.

1.4–10.6 days post-merger. From 1.4 days post-merger, AT2017gfo was monitored with a nightly cadence with the X-shooter spectrograph mounted on the European Southern Observatory’s Very Large Telescope (Pian et al. 2017; Smartt et al. 2017). The spectral range (330–2250 nm) provides the first broad spectral coverage encompassing the ultraviolet (UV) through near-infrared (NIR). The earliest X-shooter epoch (1.4 days post-merger) is remarkably well modelled as a blackbody continuum in the UV, optical, and NIR (Sneppen 2023) modulated by spectral features from strong r-process lines (Watson et al. 2019; Gillanders et al. 2022). Over the subsequent days and weeks, the continuum fades, the blackbody approximation becomes increasingly worse and more dominated by the NIR emission, where a series of spectral components emerges more prominently. Many spectra were taken near-contemporaneously with the X-shooter datasets, for example NTT (Smartt et al. 2017; Valenti et al. 2017), Magellan (Shappee et al. 2017; Nicholl et al. 2017), VLT/MUSE (Tanvir et al. 2017), Gemini-South (Chornock et al. 2017), and SOAR (Nicholl et al. 2017). The various spectra are useful for validating the robustness of the spectral reductions, but the uniquely broad wavelength coverage and the high signal-to-noise ratio of the 8.2 meter VLT with the excellent throughput of X-shooter means that the X-shooter data provide the best constraints overall. A comparison across the approximately one-hour difference in observing time between these observations can help us to constrain the reverberation wave of recombining species, as we detail in Sect. 5.2.1. We therefore used the SOAR and Gemini-South spectra (Nicholl et al. 2017; Chornock et al. 2017) from 1.47 days in addition to X-shooter as a complimentary data set at this epoch. The data delivered by HST (Tanvir et al. 2017), which provides a uniquely complementary spectral constraint in the telluric band around 1.4 μm (particularly using the Wide-Field Camera 3 Infrared channel with the grism G141). This region is difficult to confidently constrain from ground-based observations, and therefore, we include the 4.9, 9.4 and 10.6 day post-merger HST spectra in our data series.

For consistency in the comparison of the spectra, we applied a single Galactic dust extinction correction, using the reddening law of Fitzpatrick (2004), with R_V = 3.1 and E(B − V) = 0.11. This E(B − V) is the uncertainty-weighted average of the Green et al. (2018) and Planck Collaboration XLVIII (2016) constraints on the Galactic dust extinction along the line of sight. Similar to all previous reductions, the dust extinction at the location of AT2017gfo in the host galaxy was assumed to be negligible (for detailed arguments on the dust extinction, see Sneppen et al. 2023c).

3. One-day transition and importance of early spectra

To predict the expected emergence times of features, we show the average ionisation states for different r-process elements as a function of time under LTE conditions in Fig. 2. The elements include the r-process elements on the left side of the periodic table (i.e. groups I–III of periods 5 and 6, including the lanthanides). We include ₃₇Rb and ₅₅Cs for comparison. Except for these two, these are the elements that produce the strongest lines because they have only a few valence electrons, low-lying energy states, and small partition functions (i.e. Sr, Y, Zr, and Ba), and the elements with a high optical/NIR opacity with many low-lying states (i.e. the lanthanides). These latter two sets of elements display a rapid change in ionisation around one day post-merger.

Fig. 2. Average ionisation inferred for the r-process elements in groups I–III. These elements are likely to produce the strongest lines under LTE conditions. We show them as a function of the observed time post-merger. The ionisation states were calculated using the Saha equation, where we assume delectron density, ne(t) = 108 cm−3(t/1.5 days)−3 and the blackbody temperature evolution found for AT2017gfo (Sneppen et al. 2023b). The vertical blue and red vertical regions indicate the observations with X-shooter and other instruments (see the spectral series in Fig. 1). At certain temporal intervals, particularly at around 1 day post-merger, the r-process elements that dominate the opacity transition rapidly between ionisation states. For t ⪆ 5 days post-merger, NLTE-effects are likely to become important, implying that neutral, single-, and double-ionised species of these elements may be present (Pognan et al. 2022).

For Fig. 2, we assumed the temperature evolution inferred from the Doppler-corrected blackbody continuum (Sneppen et al. 2023b) with a homologously expanding (and thus diluting) electron density, n_e(t) = 10⁸ cm⁻³(t/1.5days)⁻³ (Sneppen et al. 2023a). At early times, the temperature cools dramatically with T ∝ t^−α with α ∈ [0.5; 0.8] (see Table 1, Drout et al. 2017; Waxman et al. 2018; Sneppen et al. 2023a), and the abundance of lowly ionised species therefore rapidly transitions from being a negligible population to the dominant species (and thus entering the domain of new prominent spectral features, such as discussed in the models in Tak et al. 2024). Specifically, this sensitivity in LTE to temperature (or equivalently time) follows from the Saha equation, where the relative density of ions n_i + 1/n_i (with the subscript indicating the state of ionisation) is given by the electron density, n_e, temperature, T, the ionisation energy, 𝒳_i, the partition functions for each ion, Z_i + 1 and Z_i, and the thermal electron de Broglie wavelength, λ_dB,

$$\begin{array}{c}\hfill \frac{n_{\mathrm{i}+1}}{n_{\mathrm{i}}}=\frac{2}{{\lambda }_{\mathit{dB}}^3}\frac{Z_{\mathrm{i}+1}}{Z_{\mathrm{i}}}\frac{1}{n_{\mathrm{e}}}{e}^{-\frac{{\mathcal{X}}_{\mathrm{i}}}{k_{\mathrm{B}}T}}.\end{array}$$(1)

The transition between ion states is especially sensitive to the temperature because of its exponential dependency. Due to the rapid drop in temperature over the first hours and days (see 4), the ionisation transition is expected to be abrupt under LTE conditions. For Sr II at 1.4 days (where the $n_{\mathrm{Sr}\mathrm{I}\mathrm{I}\mathrm{I}}/n_{\mathrm{Sr}\mathrm{I}\mathrm{I}}$ fraction is near unity) the ionisation energy is 𝒳₁ ≈ 30k_BT, which implies that a subtle decrease in temperature results in a strong change in the Sr II number density. As the ion recombines, the column density, and hence the optical depth, τ, of the lines produced by the relevant species, will change by orders of magnitude. This implies a drastic change in the transmission (I = I₀e^−τ) and thus a very rapid appearance or disappearance of features in LTE.

Thus, the timing of the early spectral data from 0.5–1.5 days post-merger (see Fig. 2) allows us to probe the ejecta before, during, and after this rapid recombination wave passes through the ejecta, which transition from double- to single-ionised states. The subsequent spectra constrain the evolving photospheric phase, the recession of the photospheric surface, and the transition to the nebular regime. The step change between the dominant ionisation states under LTE followed by longer periods of relatively little evolution is a key observational difference to NLTE features. NLTE features evolve (by their nature) more smoothly and continuously during their (comparatively) slower formation time and their subsequent evolution.

Fig. 3. Spectral series of the kilonova AT2017gfo. The best-fit blackbodies are overlaid. The temporal span here covers the photospheric epochs (i.e. spectra up to 5.4 days post-merger), where a blackbody model approximates the continuum well. The blackbody temperatures and the spectral features are listed in Table 1 and in Sect. 5, respectively.

4. Temperature evolution of the blackbody

The photometric (e.g. Fig. 3 in Drout et al. 2017) and spectroscopic (e.g. Pian et al. 2017) observations of AT2017gfo are remarkably well described by a blackbody framework for the observed epochs. For instance, the first X-shooter spectrum from 1.43 days displays a temperature consistent at the 2% level when inferred from the UV, optical, or NIR separately (i.e. from wavelengths blueward and redward of the spectral peak), which leaves little territory for perturbations away from the blackbody (Sneppen 2023). Epoch-by-epoch estimates of the best-fit blackbody temperature that is required to match photometry and spectra were provided by Drout et al. (2017) and Waxman et al. (2018). However, as these analyses predated the discovery and identification of the various spectral features, and their modelling constituted a simple blackbody fit to the whole spectral range. Not accounting for the modulations from the spectral feature may have biased the inferred blackbody temperatures. This bias can push the best-fit blackbody to colder temperatures (when absorption in the UV/optical or NIR emission is not modelled in early spectra) or hotter temperatures (when the prevalence of emission structures in the nebular phase from one week post-merger is neglected).

Therefore, we provide an updated summary of the best-fit epoch-by-epoch temperatures in Table 1, which we visualise in Fig. 3, where the spectral modelling now includes additional parameters for the features discussed in Sect. 5. We summarise the main uncertainties in estimating the blackbody temperature in the various observed periods below.

In the first day post-merger, the flux peaks blueward of the optical spectroscopic coverage (λ ≲ 400 nm). This means that the temperature is only directly constrained from the spectral slope on the blue side of the blackbody in the earliest spectra. Fortunately, photometry from Swift satellite (probing 200–400 nm) at 0.6 and 1.1 days constrains the UV flux, showing that it peaks at 300–400 nm. Thus, while the ANU spectrum at 0.92 days for instance does not independently provide strong constraints on the blackbody temperature due to its limited wavelength coverage, the preceding, contemporaneous, and subsequent photometry places strong bounds on the radiation field of the blackbody.

When the spectral peak is constrained and during the photospheric epochs (1 days ≲ t ≲ 5–6 days), the main uncertainty in the observed blackbody temperature lies in the exact Galactic dust extinction in the direction of the host galaxy (see details of dust-extinction estimates for AT2017gfo in Sneppen & Watson 2023). Based on this work, we used E(B − V) = 0.11, which is the weighted average of Green et al. (2018) and Planck Collaboration XLVIII (2016), and we conservatively used a 15% uncertainty that includes the absolute uncertainty in this direction. In addition, we also adopted a 10% uncertainty on the R_V, based on the standard deviation of R_V values found along different Milky Way sight-lines, which we propagated throughout the uncertainty estimates of this analysis. These total dust-extinction uncertainties imply an inherent uncertainty of as much as 5% in the relative flux at 500 nm relative to 1μm, which corresponds to an uncertainty of 1–2% in the temperature inferred solely from the spectral shape and location of the spectral peak. We note that as the dust correction is the same for all epochs, the blackbody temperature uncertainty is correlated between epochs.

During the nebular epochs (from 5–6 days post-merger), the systematic uncertainties associated with identifying what constitutes the continuum and what constitutes nebular emission features will likely come to dominate. By 10.4 days post-merger, comparable flux is contributed by the emission features as from the underlying continuum, which implies that any spectral continuum modelling is fraught with uncertainty. Further complicating the constraints, the inferred spectral peak has shifted redward of 2.25 μm by this time and is thus beyond the spectral regime constrained by X-shooter. It is noteworthy that NIR follow-up from the JWST at 1–5 μm of the recent GRB230307A (interpreted as being due to a BNS merger) is well modelled as a blackbody continuum at ∼670 K at 29 days post-merger (Levan et al. 2024). This would tentatively suggest that the blackbody continuum in kilonovae may even persist in some form at late times.

5. Observed spectral features

In the following, we consider each of the main spectral features, their emergence time, velocity, and evolution. The elements associated with the proposed line identifications span the first and second r-process peaks and could potentially originate from different mass ejection channels and be distributed spatially inhomogeneously, both azimuthally and radially. This would be encoded in the emergence time, the line profiles, and in different velocities for each individual observed feature. In Sneppen et al. (2023b), we showed that the line-of-sight velocity from the Sr II P Cygni profile is highly consistent with the velocity inferred from the blackbody normalisation. This indicates that the outer ejecta (which defines the early photosphere) is highly symmetric, and it potentially shows that the light r-process element distribution may be as well. However, it is still unclear how constraining this velocity consistency is for the underlying ejecta structures, with the simulation presented in Collins et al. (2024) showing that highly symmetric ejecta is not necessarily required to produce near-spherical line-forming regions. In Sneppen & Watson (2023), we extended this analysis to the newly discovered 760 nm yttrium P Cygni profile, showing that it traces the velocity structure and evolution of its strontium counterpart. In addition to these light r-process elements, the series of NIR features noted in previous work (Pian et al. 2017; Smartt et al. 2017; Watson et al. 2019) has been linked to the second r-process peak elements (Domoto et al. 2022), but their velocity structure and their evolution has not yet been analysed within this framework. In this paper, we therefore used our spectral compilation to extend the velocity constraints to a broader temporal span, including the NIR features.

5.1. P Cygni feature at 760 nm – Y II?

The 760 nm P Cygni feature emerges 3–4 days post-merger, and in terms of central wavelength (constrained by the emission peak) and the prominence of the feature (relative to its stronger Sr II counterpart), it is well reproduced by the 4d²–4d5p transitions of Y II (Sneppen & Watson 2023). At earlier times, the higher velocities mean that the various lines of Y II blend together and thereby conceal the characteristic P Cygni spectral shape. Modelling with the TARDIS code in Sneppen & Watson (2023) showed that a prediction of the yttrium identification is that the P Cygni structure emerges clearly when the velocity drops to about 0.15–0.2c, which is 3–4 days after the merger. Tak et al. (2024) have now also identified the 4d²–4d5p Y II transitions as likely to strongly affect the line signatures around 760 nm within their models. Other transitions to the 4d5p energy level of Y II have been linked to the strong absorption observed at λ ≤ 400 nm at early times (Gillanders et al. 2022; Vieira et al. 2023). Conversely, Pognan et al. (2023) reported that the Y II transitions within their models are too weak to account for the observed 760 nm features, and they instead proposed that in NLTE conditions Rb I may be able to produce a feature at roughly the correct wavelength. However, as they used an uncalibrated theoretical line list for yttrium that has incorrect wavelengths, placing these transitions at 1.1 μm (below the stronger Sr II lines), it is unclear how constraining their conclusions on yttrium would be in a self-consistent model with the features placed at the correct wavelengths. As regards Rb I, no explicit fit of the observed feature was presented, and their NLTE models cannot extend to the early times (2.4, 3.4, and 4.4 days post-merger) where the observed feature is first revealed and rises to prominence. Finally, the Rb I transition rest wavelength (around 790 nm) is 0.04–0.05c redward from the observed emission peak at early/intermediate epochs (in contrast to all other proposed line IDs, Sr II, Ce III, La III, (and indeed Y II), which are consistent with the observed emission peaks at 1.0, 1.4, and 1.6 (and 0.76) μm). This means that the Rb I transition wavelength would be at odds with the spectral location and absorption/emission P Cygni shape (assuming Sr II, Ce III, or La III are reasonable proxies for the other observed features). In the later epochs, the Rb I transition wavelength becomes more consistent with the observed feature, which 6–8 days post-merger has evolved into what appears as a pure emission-like component.

The velocity stratification of the feature 4.4–5.4 days post-merger (see Fig. 1, first sub-panel) suggests a photospheric velocity around 0.15c. This agrees with the velocity inferred from the 1.0 μm P Cygni feature at this time. At later times, the feature evolves analogously to the 1.0 μm feature, reverberating into pure emission and fading away towards the more nebular epochs. Reverberation here means that the photons from different wavelengths of the P Cygni profile are scattered at different physical times due to light travel-time effects. Thus, as the line fades away, it does so first in absorption and subsequently in emission. We discuss this later-time evolution further in the context of the more prominent and, by extension, better constrained 1.0 μm feature.

5.2. P Cygni feature at 1.0 μm – Sr II

The Sr II 4p⁶4d–4p⁶5p triplet of strong lines producing the 1.0 μm P Cygni feature was first proposed (Watson et al. 2019) and later reproduced and extended with systematic analyses using different radiative transfer codes (Domoto et al. 2021; Gillanders et al. 2022; Vieira et al. 2023) and with full 3D line-by-line opacity modelling (Shingles et al. 2023). The strontium mass required to produce the feature is compatible with what is produced in the dynamical and spiral-wave wind ejecta from merger simulations (Perego et al. 2022). Independent evidence for the Sr II 4p⁶4d–4p⁶5p identification is that a large UV absorption is observed at λ ≤ 400 nm around 1.4 days, which could be explained by Sr II transitions from the ground state to 4p⁶5p (Watson et al. 2019), or similar lines from nearby and co-produced elements: Zr II (Gillanders et al. 2022) and Y II (Vieira et al. 2023). Further corroboration of the Sr II interpretation is the existence of a 760 nm P Cygni feature due to yttrium, as discussed in the previous section (Sneppen & Watson 2023).

An alternative interpretation has been proposed: The 1.0 μm P Cygni feature could originate from the neutral He Iλ1083.3 nm line (1s2s-1s2p transition), which in NLTE conditions, could produce a feature at a similar wavelength (see Perego et al. 2022; Tarumi et al. 2023), while Sr II would be depopulated to higher states of ionisation. In this case, the apparently corroborating features of Sr II, Y II, and Zr II at other wavelengths would need to be explained through different line identifications with other transitions. However, as we show in Sneppen et al. (2024), a He I interpretation is not compatible with the emergence time and subsequent spectral evolution of the observed 1 μm feature. This is in contrast to the situation for Sr II, as we discuss below.

5.2.1. Emergence of the 1.0 μm P Cygni feature

A key way to robustly distinguish these and indeed other proposed identifications is to analyse the temporal properties and evolution of the feature. That is, in LTE, a Sr II feature is predicted to begin its emergence only at a specific time, as the temperature falls rapidly, and thus, the dominant ionisation level will also change quickly (see Fig. 2). As mentioned above, this is in contrast to a helium interpretation, for example, which contradicts the first appearance and the subsequent evolution of the feature. In the following, we show that the early evolution of the Sr II feature (as per the predictions in Sneppen et al. 2023a) agrees with the observations. Previous analyses have focused on the X-shooter spectra, but the earlier spectra taken with the Magellan telescopes (Shappee et al. 2017), ANU (Andreoni et al. 2017), and SALT (Buckley et al. 2018) allow us to constrain the earliest appearance of the observed feature.

Given the blackbody temperature recovered in the X-shooter epoch 1 spectra (∼1.4 days post-merger), the expected ejecta electron densities, and the observed cooling rate, we argued that a Sr II P Cygni feature should first start to form in absorption around 0.8–1.0 days post-merger, with a significant absorption present by around 1.2–1.3 days post-merger (Sneppen et al. 2023a).

This time frame is observationally probed by the early spectra (see Fig. 4). The early spectra at 0.49, 0.53, and 0.92 days all display a near featureless blackbody-like continuum. In contrast, the 1.17 day SALT spectrum continuum is depressed by absorption in the red part of the optical (λ ⪆ 650 nm for a blackbody with a temperature T_obs ≈ 6400 K inferred from the location of the spectral peak). This is the first appearance of the 1.0 μm P Cygni feature, consistent with the expectation from recombination of Sr III to Sr II (Fig. 2). The observed absorption is a robust feature located in the middle of the wavelength coverage of SALT’s prime focus Robert Stobie spectrograph. The feature is also highly significant statistically: A blackbody continuum model, which adequately describes the spectrum below 640 nm, describes the spectral region at 640–680 nm poorly with a χ² null hypothesis probability of ∼10⁻²⁶ when this region is added to the fit, and it can thus be rejected at 10σ. Interestingly, the absorption reaches blueward to around λ ≈ 640 nm, which indicates that the line-forming region must extend out to quite high velocities: 0.40–0.45c (given absorption from ejecta expanding along the line of sight). This is the fastest kilonova ejecta component observed so far. It is 0.05c faster than the fastest component ejecta found from the 1.0 μm P Cygni feature at 1.43 days post-merger (Sneppen et al. 2023b; Vieira et al. 2024). While rapidly receding, this ejecta velocity is still within the velocity range predicted from numerical simulations of kilonova dynamical ejecta (e.g. Bauswein et al. 2013).

Fig. 4. Emergence of the 1 μm P Cygni feature in time. The selected spectra are from 0.48 days (LDSS, Shappee et al. 2017), 0.92 days (ANU, Andreoni et al. 2017), 1.17 days (SALT, Buckley et al. 2018), 1.43 days (X-shooter, Pian et al. 2017), 1.47 days (SOAR, Nicholl et al. 2017 and Gemini-South, Chornock et al. 2017), and 2.42 days (X-shooter, Smartt et al. 2017). The spectra before 1 day do not deviate strongly from the blackbody continuum. However, between 0.92 and 1.17 days, strong absorption appears, while the bulk emission begins to emerge with a delay of around 10 hours. This is the predicted emergence time-frame for both emission and absorption (and the early spectral shape) for a feature produced by LTE Sr II due to a recombination wave passing through the ejecta (Sneppen et al. 2023a), and it is inconsistent with a He Iλ1083.0 nm interpretation (Sneppen et al. 2024).

Due to the long time delays across the ejecta, we expect to first observe the blueshifted outflowing material undergo an ionisation change, after which the more distant matter, producing the emission component of the P Cygni feature, will be observed to undergo change in ionisation state with a delay of about 10 hours. Thus, the earliest feature should form in very blueshifted absorption first, gradually evolving into a more pronounced absorption/emission P Cygni feature. In reality, the first spectra with coverage of the Sr II emission peak (X-shooter epoch 1) at 1.44 days display only very weak emission, as expected due to this reverberation effect. At 1.47 days (just 1 hour after the epoch 1, X-shooter spectrum), the Gemini spectrum shows weak but slightly stronger emission that begins to form (Sneppen et al. 2023a). Consistent with this idea, by epoch 2 (2.4 days post-merger), the feature has evolved into a clear P Cygni shape with a pronounced emission component (Fig. 4).

That the observed temperature evolution of the blackbody spectra successfully predicts the line formation time indicates that the ejecta are very near to LTE conditions and that the blackbody temperature is a good probe of the ionisation/excitation temperature of the ejecta (we revisit this further in Sect. 6.1).

5.2.2. Evolution of the 1.0 μm P Cygni feature

The 1.0 μm P Cygni feature, as the first clearly detected spectral feature, allows us to trace the evolution of the ejecta from 1.17 days post-merger. The feature evolves in several key respects (Fig. 1, second sub-panel).

First, the outer and inner (photospheric) velocities of the P Cygni feature decrease with time (see Fig. 1 or Watson et al. 2019) because the line-forming region recedes deeper into the ejecta as the outer layers become optically thin. The photospheric velocity decreases from 0.28c (at 1.4 days) to around 0.15c (at 5.4 days). The optical depth of the absorption remains borderline optically thick throughout the early photospheric epochs with τ ∈ [1; 3] from 1.17–3.4 days, after which it becomes optically thin, and it has ultimately faded away by 8 days post-merger.

Another important change in the feature is the increasingly prominent nature of the emission. This is caused by the time-delay effect of the line noted above and by the rapidly cooling photosphere, but also by the rapidly fading continuum; these effects create spectral features that are quite sensitive to reverberation effects (Sneppen et al. 2023a). This means that the amplitude of the P Cygni absorption/emission is modulated not only by the changing ionisation state of the gas, but also by the changing strength of the underlying continuum, which itself rapidly fades with time. Modelling the changing continuum only produces a minor effect at times t < 2 days because of the slower flux evolution on the Rayleigh-Jeans side of the blackbody. However, this becomes pronounced at later times (t > 3 days), when the blackbody peak has shifted redward and into the feature, because the wavelength-dependent fading evolves rapidly near the blackbody peak (Sneppen et al. 2023a). This effect is investigated in more detail in a future work (McNeill et al., in prep.). In the early nebular phase epochs (5.4–8.4 days post-merger), the 1.0 μm feature fades away first in absorption and subsequently in emission.

5.3. P Cygni feature at 1.4 μm – La III?

A feature around 1.4 μm has been suggested to be due to the 5p⁶5d–5p⁶4f transitions in La III (Domoto et al. 2022). The feature is poorly constrained from the ground-based spectra because most of the profile resides within the telluric region. However, at 4.9 days, HST spectra were taken of AT2017gfo (Tanvir et al. 2017), which revealed a P Cygni-like feature with an emission peak around 1400 nm and a blueshifted absorption around 1300 nm (see the third sub-panel in Fig. 1 and Fig. 5). It is likely that this feature is also present in other epochs, potentially complicating the modelling of the 1.6 μm feature. Significant absorption from the continuum is noticeable around 1250 nm as early as 2.4–3.4 days post-merger (see Fig. 1, main panel). Notably, this is so late that the low blackbody radiation temperature would suggest that La III should be almost entirely absent under LTE conditions. This identification therefore seems to require either significant NLTE effects at a few days post-merger or that the line-forming region for this feature originates from an inner and potentially hotter ejecta that allows for doubly ionised lanthanides. We discuss this prospect further in Sect. 6.3. In the subsequent HST spectra at 9.4 and 10.6 days, a feature around 1.4 μm is still detected, but on a continuum that is 5–10 times fainter (Tanvir et al. 2017).

Fig. 5. HST spectrum of AT2017gfo at 4.9 days post-merger from Tanvir et al. (2017). The grey shaded region indicates the continuum constraints from the interpolation of the X-shooter spectra taken at 4.4 and 5.4 days post-merger. Regardless of the exact continuum modelling, a P Cygni absorption/emission feature around 1400 nm is clearly visible. Domoto et al. (2022) suggested that this feature could be due to La III lines (here indicated with dashed blue lines). The characteristic velocity of this feature is relatively low, vph = 0.07 ± 0.01c, but may be systematically biased by the strong nearby features (including the prominent reverberation emission peak around 1.1 μm).

Given the relatively limited wavelength range of the HST spectra, the continuum flux around this P Cygni feature is difficult to ascertain directly. However, the spectral evolution is relatively minor between the temporally adjacent X-shooter spectra (at 4.4 and 5.4 days post-merger), which allow strong modelling constraints on the continuum. Thus, within the range of allowable blackbody models at 4.4 and 5.4 days, we can model the continuum, and by extension, probe the nature of the 1.4 μm P Cygni feature (see Fig. 5). Fitting the P Cygni profile with this range of continuum models would suggest that the photospheric velocity is v_ph = 0.07 ± 0.01c. This is a relatively low velocity that is around a factor of two lower than the contemporaneous light r-process P Cygni features at 760 nm (Y II, Sneppen & Watson 2023) and at 1 μm (Sr II, Watson et al. 2019). This may suggest that the optically thick region at these NIR wavelengths resides deeper within the ejecta (where the expansion velocity is lower). However, the velocity we inferred may also be biased by nearby features, for instance, by the feature around 1.6 μm and the prominent reverberation emission peak around 1.1 μm, which could create a bias towards lower velocities. In addition, there appears to be a marginal feature at around 1.2 μm that is unaccounted for so far.

5.4. Feature at 1.6 μm – Ce III?

The 1.6 μm feature (see Fig. 1, fourth sub-panel) was interpreted by Domoto et al. (2022) as potentially being due to a combination of lines from Ce III, (predominantly 4f²–4f5d transitions), which also produces several strong lines in the corresponding wavelength range of HR 465 (the binary GY And), which is a metal-enriched star (Tanaka et al. 2023). Conversely, Spectroscopic r-Process Abundance Retrieval for Kilonovae (SPARK) models in Vieira et al. (2024) showed that Ce II can produce a spectral feature at a similar location in the NIR, but likely of lower prominence. Gillanders et al. (2023) additionally noted the presence of nearby intrinsically weak [I II] lines that might contribute significantly to a feature at late times. No statistically significant feature is present at this wavelength in the 1.4 day X-shooter spectrum. Curiously, the first weak emission has already begun to emerge around 1.47 days in the Gemini spectrum, apparently in the narrow time-frame associated with the transition from double- to single-ionised r-process line-forming species. Over the subsequent days, the emission feature slowly becomes more prominent than the continuum. As in the case of La III, the blackbody temperature and LTE conditions would suggest that no significant Ce III population is present in the photosphere from 2.4 days onwards (see Fig. 2). Again, this would require significant NLTE effects at this time. Alternatively, the NIR line-forming region might originates from hotter ejecta.

The feature and its mostly emission-like nature have been discussed previously (Smartt et al. 2017; Watson et al. 2019; Sneppen et al. 2023b; Gillanders et al. 2023). The central wavelength of this feature shifts redward with time, starting at 1.52 μm (2.4 days post-merger) and shifting to 1.58 μm by 7–10 days post-merger (Gillanders et al. 2023). The velocity width of this feature does not change substantially over time when modelled as a Gaussian emission feature, with a σ width of only around 0.04–0.05c. This is very narrow compared to the velocity measurements of other features, including even the redder 2.1 μm emission feature, which at face value appears to be broader, and which we consider next. We note that this might be a full (and broader) P Cygni absorption/emission feature modified by possible surrounding features. However, there is no very clear signature for this hypothesis because the corresponding absorption would reside within the telluric regions and would be highly degenerate with the 1.4 μm P Cygni feature emission peak.

5.5. Feature at 2.1 μm – [Te III]?

The 2.1 μm feature (see Fig. 1, fourth sub-panel) is well modelled as pure emission and was interpreted as a [Te III] fine structure line by Hotokezaka et al. (2023) when they modelled the spectra from 7.4–10.3 days post-merger. Gillanders et al. (2023) also identified [Te III] as a plausible line-forming species, but further argued that due to the width of the line, the feature would likely have to be produced by two separate lines to match the narrowness of the width of the 1.6 μm feature. This feature (or a similar feature at more or less the same wavelengths) is also statistically significant in earlier epochs. It is relatively minor in the first spectra with NIR coverage (X-shooter epoch 1, 1.4 days post-merger), grows in prominence with time (equalling the continuum at 5.4 days post-merger), and in the later, probably more nebular phase, becomes the most prominent feature. The central wavelength drifts redward with time, starting at 2.0 μm at 1.4 days, but shifting to 2.1 μm by 10 days post-merger.

A feature at a similar wavelength has also been detected in late JWST spectroscopy associated with a long-duration GRB and was argued to be a kilonova from a merger (Levan et al. 2024; Yang et al. 2024). This feature was also tentatively associated with [Te III] (Levan et al. 2024; Gillanders et al. 2023), and may possibly have substructure, which would be indicative of more than a single line contributing.

Despite the considerable evolution in prominence of this feature in AT2017gfo, the velocity width of the features remains largely constant at around 0.07c throughout all epochs from 1.4–10.3 days post-merger. If the origin is several lines at different wavelengths, this would naturally suggest an even lower expansion velocity, which, as argued in Gillanders et al. (2023), could lower the velocity to be consistent with the measured width of the 1.6 μm feature. This velocity is particularly striking in two respects. First, in early epochs, the velocity width of the 2.1 μm feature is significantly smaller than the velocity inferred from the continuum, the 760 nm P Cygni feature or the 1 μm P Cygni feature, but it is consistent with the 1.4 μm feature (the reddest of the P Cygni features) from the HST spectra. This suggests that the r-process elements creating this feature may not be co-located spatially with the light r-process elements producing the optical P Cygni features. This either requires these lines to be located spatially near the equatorial plane (where the line-of-sight velocity is low due to projection effects for a near-polar observer) or deeper within the ejecta (visible due to the lower opacity at longer wavelengths). Second, the velocity width of the feature does not evolve significantly. If the elements producing the feature were confined to the equatorial plane of the merger (which could explain the lower projected velocity compared the optical P Cygni features), then the relative velocity should still decrease as the line-forming region recedes inwards. If the line-forming region is located radially deeper in the ejecta, this non-evolution still requires a constant velocity width of the line-forming region. At the same time, the observed drift of the central wavelength requires that the line-forming region is initially offset from the merger centre or perhaps that reverberation effects are significant. Ultimately, modelling both a drifting central velocity and an invariant velocity width seems difficult to reproduce, and this requires further study.

5.6. Minor and tentative features

In the high S/N X-shooter spectrum from 1.4 days, a few smaller, narrower spectral wiggles are noticeable in the absorption of the P Cygni feature. The most noticeable peaks lie at around 780 nm and 860 nm. These spectral features are persistent and appear in every single exposure taken throughout epoch 1 (Sneppen et al. 2023a). Similar wiggles are seen in coincident MUSE (Tanvir et al. 2017) and NTT (Smartt et al. 2017) spectra, suggesting that they may not be artefacts of the processing pipeline, but may be the contribution of less prominent transitions that are yet to be identified. These spectral wiggles with their smaller inferred (projected) velocity scales could potentially follow from ejecta that expand equatorially given the polar inclination angle.

Alternatively, the positive and negative interference of various lines centred at different wavelength could in conjunction produce spectral wiggles with a somewhat smaller velocity/wavelength width than that of any individual line. These minor wiggles indeed have counterparts at similar wavelengths in the subsequent spectra. The 780 nm wiggle has a clear spectral counterpart in the Y II 4d²–4d5p transitions, which produces the 760 nm P Cygni feature from 3.4–6.4 days when the characteristic velocities decrease and the lines have deblended. A minor feature around 860 nm may also be noticeable in intermediate epochs (4.4–7.4 days post-merger) as a small peak resides between the 0.76 and 1 μm features (see Fig. 1, main panel). This weak feature is particularly poorly constrained across the photospheric epochs (especially in epoch 2) as it is located within the absorption valley of the 1 μm P Cygni feature. Coincidentally, 860 nm is quite near the rest-wavelength of the strong NIR Ca II triplet. Even a minor Ca abundance, Ca/Sr ratio < 0.002 in mass fraction, could produce a feature of this prominence, which would remain within the calcium upper bounds for AT2017gfo argued in Domoto et al. (2021). However, the width of the feature does not fit the lines identified with Sr and Y well, and if the feature really were to be due to Ca II, the modelling is complex and beyond the scope of this paper.

At 1.2 μm, a weak spectral wiggle is present in the spectra from 4–8 days post-merger (as previously noted by Tanvir et al. 2017; Sneppen et al. 2023b; Gillanders et al. 2023). This feature is particularly poorly constrained because it resides between the prominent 1 μm reverberated emission peak and the absorption of the 1.4 μm P Cygni feature. It is worth noting that in the metal-enhanced star HR 465, which Tanaka et al. (2023) analysed for its lanthanide features, a few prominent (currently unidentified) lines are present around 1.2 μm, although they did not associate those lines with this feature.

6. Implications for the ionisation structure

Having considered the various features and their observed structure, we now explore the information revealed when we combined the various spectral identifications and the complete dataset. First, in Sect. 6.1, we discuss how the emergence of the 1 μm feature constrains the ejecta to be close to LTE conditions. Second, in Sect. 6.3, we discuss the complementary nature of comparing the optical and NIR velocity structures. Lastly, in Sect. 6.4, we discuss how the spatial distribution of elements may be constrained from the various observed P Cygni features.

6.1. Quantifying the validity of LTE at early times

Perhaps the most striking aspect of the many early spectra taken of AT2017gfo is their relative simplicity. A simple blackbody empirically describes the continuum very well at early times (e.g. the first X-shooter spectrum at 1.43 days has percent-level consistency in the inferred temperature from the UV through the NIR, see Sneppen 2023). The spectral perturbations away from the blackbody are produced in the surrounding atmosphere by the LTE r-process species with the strongest individual transitions (e.g. elements from the left side of the periodic table with low-lying energy-levels and small partition functions, see Watson et al. 2019; Domoto et al. 2022). While both the continuum and the spectral features are thus well described by LTE populations at early times, it is still unclear how strong NLTE effects are permitted by observation.

By studying the emergence of features, we can now begin to constrain the allowed level of deviations from LTE conditions at early times observationally. That is, the blackbody temperature not only predicts the dominant ionisation level (and thus, which spectral features are present), but the temperature evolution accurately predicts the transition time between ionisation states (and thus, when features form), as shown in Fig. 6. In Sect. 5.2.1, the spectral series allowed us to trace the formation of the 1 μm feature on an hourly basis by first tracing the formation of the absorption component at around 1.0 day and subsequently (with a 10 hour time delay) the emission component, which started to form from around 1.4 days. Both of these time frames were predictions in Sneppen et al. (2023a) due to an LTE recombination wave from Sr III to Sr II. This empirically validates that the radiation temperature of the blackbody is well coupled to the ejecta temperature because for every 100 K of difference invoked between these temperatures, the emergence time would have to change by around one hour. When we assume that the electron density is unknown to an order of magnitude (i.e. propagate a total uncertainty in n_e of 1.0 dex), this corresponds to an uncertainty of 2 hours in emergence time towards either earlier or later formation. Within this level of uncertainty, the emergence time in absorption and emission is as predicted for LTE, implying that the difference between the excitation/ionisation and blackbody radiation temperatures must be within ∼200 K at these times, that is, it must be consistent at a level of some percent.

Fig. 6. Ionisation degree of Sr in LTE as a function of temperature and electron density. We overlay the approximate locations in parameter landscape of AT2017gfo with time. The lines indicate the evolution modelled as constraints from reverberation effects within a single epoch. The temperature error bars are for the part of the ejecta with line-of-sight expansion (Sneppen 2023), and the density assumes a spherically symmetric, homogeneous, and homologously expanding high-velocity (0.2c ≤ v ≤ 0.3c) ejecta component of 0.01 M⊙. The electron density can only be estimated to an order of magnitude, but regardless of the exact ne, a rapid step change in ionisation state occurs over a narrow interval in temperature around 4500 K. The 1 μm feature is not detected in observations prior to 1 day (marked in red), but from 1.17 days and forward, the feature is present (marked in blue).

These strong constraints on the temperature are attainable because the timing of ionisation transitions depends so strongly on the temperature and is only relatively weakly affected by the electron density, as illustrated in Fig. 7. Specifically, this shows the recombination transition to Sr II – and, by extension, the formation time of a Sr II feature, given Eq. (1), typical KN electron-densities, and for various relations between electron temperature and the blackbody radiation temperature. For comparison, we calculated the fraction of strontium required as Sr II to produce the observed Sobolov optical depth τ = 1, where the strontium density was computed using the total strontium mass estimate M_Sr = 6 × 10⁻⁶ M_⊙ (as inferred from in the line-forming region at 2.4 days, Gillanders et al. 2023) homogeneously distributed between 0.2c < v < 0.4c. We used the mass estimate at 2.4 day as a benchmark as this is the first epoch where Sr II is the dominant ionisation state, but we note that there are large uncertainties in the exact strontium densities and indeed in its radial structure. The y-axis error bars in Fig. 7 thus solely represent the statistical uncertainty in the optical depth of the observed feature from each spectrum as modelled by the Sr II 4p⁶4d–4p⁶5p lines given LTE level populations.

Fig. 7. Fraction of Sr II under LTE conditions as a function of time post-merger. On the upper x-axis, we illustrate the physical time of scattering probed by the 1 μm feature in absorption and emission by the various spectra from ANU, SALT, VLT/X-shooter, and Gemini FLAMINGOS-2. The dotted lines indicate that no feature has clearly emerged (prior to ∼1.4 days in both emission and absorption). We indicate the required fraction of Sr II needed to produce the optical depths observed in each spectrum (assuming MSr = 6 × 10−6 M⊙ homogeneously distributed across 0.2c < v < 0.4c; see main text). The ANU (X-shooter) and SALT (Gemini) spectra indicate a rapid transition in absorption (emission) from an optically thin to optically thick regime around 1.3–1.5 days. For the feature to emerge around 1.4 days post-merger, the ionisation temperature is required to be nearly identical to the radiation temperature inferred from the blackbody emission. For the feature to emerge at near-contemporaneous physical time in the polar and equatorial ejecta, a near-isotropic temperature is required.

However, especially at early times, the evolution of τ(t) is very constrained under LTE. In particular, we must i) form a feature on a timescale of hours around 1.3 days, ii) maintain an optically thick line on a timescale of a few days, and iii) subsequently become optically thin and fade away at late times. This first aspect is observationally constrained in the polar ejecta by the ANU and SALT spectra, which are marked by no clear absorption (τ < 0.08; 3σ upper limit) and a strong absorption (τ = 2.2 ± 0.6), respectively. In physical time, this represents an increase of an order of magnitude in the optical depth of the line from 1.3 to 1.5 days. The optical depth of the emission feature can similarly be estimated using the P Cygni framework under the assumption of pure scattering. In this case, the optical depth increases from τ < 0.09 (X-shooter; 3σ upper limit) to τ = 0.8 ± 0.2 (Gemini FLAMINGOS-2) from 1.43 days to 1.47 days. This highlights a similarly rapid line formation in the equatorial ejecta, but now confined to the timescale of an hour due to the higher cadence around this period.

It is interesting to note that the X-shooter NIR-arm data constrain the scattering of light (i.e. emission) at an earlier physical time than the optical light from the SALT spectrum (i.e. the absorption), although the former spectra were obtained 6 hours after the latter.

6.2. The temperature isotropy of the line-forming region

We note that the fact that we can successfully predict both the formation time in absorption and (in the delayed) emission by accounting only for light travel time effects indicates that there are no large temperature gradients from the polar to equatorial ejecta. Specifically, it begins to form in absorption (i.e. the polar ejecta) between ANU and SALT spectra (indicating a formation around 1.3–1.5 days in physical time post-merger) and is beginning to form in emission around the time of X-shooter and Gemini (1.43–1.5 days). This near consistency in emergence to within 2–3 hours suggests that the polar and equatorial line-forming region must reach the temperature for which Sr II can first exist (i.e. T ∼ 4500–5000 K) at nearly the same physical time. This requires a consistency in ejecta temperature between the polar and equatorial ejecta within 200–300 K, that is, an angular consistency in temperature at a level of a few percent.

This again illustrates that the remarkably simple framework of a single-temperature blackbody (as suggested by the spectral shape) may surprisingly capture the complexities of the ejecta. For comparison, in the most recent 3D models presented in Shingles et al. (2023), large temperature anisotropies exist between the equatorial and polar ejecta at a level of up to 50%. This class of models cannot produce the emergence of the 1 μm features observed in AT2017gfo. Thus, the discovery of these new observables measuring the ionisation temperature and the temperature anisotropy around ionisation transitions can inform and constrain future modelling of BNS mergers, specifically, the thermalisation of decay products and the radiative transfer in KNe.

At intermediate times, complexities emerge as we may begin to see both single- and double-ionised species producing spectral features. This is the focus of our next discussion.

6.3. Complementary velocity and ionisation constraints

Sr II and Y II can only co-exist with the higher ionised species La III and Ce III over a small range of temperatures at LTE because the first and second ionisation energies are similar (see Fig. 2). It is therefore unlikely that these species are co-spatial, with their features appearing together across the wide temperature range probed at the different photospheric epochs. As argued above, if both optical and near-infrared identifications are correct, their co-existence in LTE would require their line-forming regions to be spatially separated by ejecta with different physical temperatures. The more highly ionised species would only exist, for example, in deeper (and thus, likely hotter) parts of the ejecta. These outer and inner ejecta could plausibly be simultaneously probed by different wavelengths due to the expected decline in overall opacity towards the NIR. This spatial separation between the species would map onto different expansion velocities, such that the velocity difference between the optical P Cygni line profiles (from single-ionised species) and the near-infrared features (double-ionised species) could be due to different ejecta temperatures.

This interpretation can be tested for its self-consistency as the temperature must both be consistent with 1) the ionisation level of the proposed species and 2) the continuum flux given the emitting area set by the photospheric velocity of the line and free expansion. In the UV-optical, the temperature and velocity are individually well constrained by the spectral shape and location of the blackbody peak, but in the Rayleigh-Jeans tail toward the NIR, the temperature and emitting area can be varied along a line of degeneracy to obtain a similar continuum (Fig. 8). To produce a flux at one specific wavelength, there is a complete degeneracy between emitting area and temperature, but the NIR spectral slope is sensitive to the temperature. For example, the NIR temperature is consistent with the UV/optical temperature at the percent level at 1.4 days (Sneppen 2023). For double-ionised species to be present at typical electron densities (at 4 days, this is around 10⁷ cm⁻³), the electron temperature would need to be around 4000 K or hotter. This is much hotter than the optical-UV temperature at this epoch. However, this ejecta temperature is broadly consistent with the NIR continuum flux for the observed velocity inferred from the 1.4 μm line. It is worth noting, however, that these NIR-continuum temperature constraints for any fixed velocity are likely upper bounds because they neglect modelling the contribution of fluorescence from all layers at greater radii, which requires detailed radiative transfer modelling. The actual radiation temperature could therefore be lower.

Fig. 8. Degeneracy between the blackbody temperature and the inferred emitting area (or velocity) for a given NIR continuum. The X-shooter epoch 4 spectrum is overlaid with blackbodies of various temperatures and areas that produce similar NIR flux levels. The observed NIR features can be modelled as P Cygni or pure emission components, but make the exact underlying continuum difficult to estimate precisely. If the near-observed infrared features are produced by La III and Ce III for LTE conditions, these species must reside in relatively hot ejecta. Hotter ejecta produces a higher NIR continuum flux, which would conversely require a smaller emitting area or lower velocity to match the observed continuum. While the optical P Cygni features are consistent with the observed blackbody temperature and velocity, these NIR identifications would need higher temperatures and be spatially distinct from their optical counterparts.

If these identifications of double-ionised species are correct, this diversity in the ionisation levels of species emerging at intermediate epochs would constitute the first observational indications of a multi-temperature ejecta/blackbody. That the more highly ionised species first appear (and increase in prominence) after the observed optical-UV blackbody temperature cools away from the dominant temperature regime of the species (T ≲ 4000 K) might be considered surprising at first, but it may be explained by the relevant species initially being hidden within optically thick inner regions (e.g. Vieira et al. 2024). On the other hand, the appearance of more highly ionised species in the NIR could be a distinctly NLTE phenomenon. The gradual emergence and continuously increasing prominence are indeed key observational properties of NLTE features (in contrast to the LTE evolution). Towards late times (t ≳ 5 days), NLTE effects may boost the abundance of higher-ionisation species (Pognan et al. 2022). Finally, the identifications with more highly ionised species may simply not be correct, and we await the correct identifications of these lines. However, in this case, the apparently low velocities of these lines would still need to be explained.

6.4. Angular distribution of line-forming species

The observation of several P Cygni structures in the line-forming regions (see Fig. 1, right panels), requiring both absorbing and emitting components provides indications of the spatial geometric distribution of the elements. As we observed this object with a near-polar inclination angle θ_inc = 22°   ±  3° (Mooley et al. 2018, 2022), the blueshifted absorption and the rest-wavelength emission here closely correlate with the more polar and equatorial components of the ejecta, respectively.

With the emergence of the 1 μm feature as detailed in this compiled spectral series, we further tightened the connection between Sr II and the observed emission of the 1 μm P Cygni feature. This requires a sizeable component of strontium (i.e. a first-peak r-process element) to be located in the equatorial plane. Furthermore, the emission part of the 0.76 μm P Cygni feature, again interpreted as due to first-peak r-process species, Y II (Sneppen & Watson 2023), supports the interpretation of light r-process elements in the equatorial plane in addition to the polar (absorption) component. Equatorially distributed light r-process elements like this are certainly not a general expectation from BNS merger simulations, which typically produce light r-process ejecta principally in the polar direction (e.g. Kasen et al. 2015; Collins et al. 2023; Just et al. 2023). Lastly, if the 1.4 μm feature is P Cygni-like, as argued in Fig. 5, and if the La III interpretation suggested in Domoto et al. (2022) is correct, this would similarly suggest that second-peak r-process elements are also spatially distributed in a manner that allows both absorption and emission. This would apparently suggest a similar abundance of lanthanides at and above the NIR photosphere in the polar ejecta as in the more generally expected equatorial (emission) component, which is often found in simulations (e.g. Just et al. 2023).

Thus, it is likely that the near-spherical distribution of Sr II in AT2017gfo found in Sneppen et al. (2023b), which contrasts with predictions from simulations, is not exclusively an issue in strontium. However, the full landscape of these geometrical constraints has yet to be explored. The development of full 3D radiative transfer modelling coupled self-consistently to hydrodynamically simulated kilonovae, such as that pioneered in Shingles et al. (2023), Collins et al. (2024), may provide an avenue to bridge the interpretability of BNS merger simulations and the observed spectra of AT2017gfo and future KNe. These models currently do not predict or reproduce P Cygni structures due to the large spatial asymmetries of the various line-forming species from hydrodynamic simulations.

7. Conclusion

We have compiled a detailed early to late spectral series of AT2017gfo to investigate the emergence and evolution of spectral features from freshly synthesised r-process elements. Combining these data has revealed the points listed below.

1. The earliest chronological appearance of the 1 μm P Cygni feature (in the SALT spectrum from 1.17 days) reveals the fastest kilonova ejecta component yet discovered at 0.40 − −0.45c.

2. The 1 μm emerges as blueshifted absorption first, gradually evolving into a more pronounced P Cygni-like feature. This follows the predictions in Sneppen et al. (2023a) for the exact emergence time and early spectral shape from strontium undergoing a III–II recombination wave in LTE. The transition time between ionisation states empirically shows that the radiation temperature of the blackbody is well coupled to the ionisation/excitation temperature of the particles themselves, and it observationally constrains the ejecta to be close to LTE.

3. The observed spectral features at 0.76, 1.0, 1.4, 1.6, and 2.1 μm were described and compared in their velocity structures. The two optical P Cygni features are broadly similar in velocity over time, while their NIR counterparts appear to have lower velocities. This difference in velocities may relate to a difference in the proposed identifications for the optical (Sr II, Y II) and NIR features (La III, Ce III, Te III), that is, they are first versus second r-process-peak elements and single- versus double-ionised. This may indicate a possible spatial separation.

4. The abundance of P Cygni structures in the line-forming regions, requiring both absorbing and emitting components, indicates the spatial geometric distribution of the elements. In particular, the Sr II and Y II emissions requires a sizeable component of first-peak elements to be located in the equatorial plane. This element distribution is in contrast with the typical spatial predictions from simulations. However, this is likely no longer exclusively an issue in Sr II (or Y II), as these dual emission/absorption components are also seen for the 1.4 μm P Cygni feature. Furthermore, that the recombination of Sr III to Sr II in absorption and emission occurs at the same physical time (i.e. 1.3–1.5 days post-merger) indicates a temperature in the line-forming-region around this time that is isotropic at the level of 5% from the polar to the equatorial ejecta.

Data availability

The newly reduced spectra are available at https://github.com/Sneppen/Kilonova-analysis

Acknowledgments

The authors would like to thank Igor Andreoni, Steven Crawford, Jonathan Selsing, and Jesse Palmerio for sharing, introducing and clarifying the various datasets and data-reductions. We would further like to express our gratitude to Stuart Sim, Christine Collins, Luke Shingles, and Fiona McNeill for discussions on the spectral modelling. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant DNRF140. AS, DW, RD, and KEH are co-funded by the European Union (ERC, HEAVYMETAL, 101071865). Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. PV and AM acknowledge support from the National Research Foundation of South Africa, and AM financial support from the Swedish International Development Cooperation Agency (SIDA) through the International Science Programme (ISP) – Uppsala University to the University of Rwanda through the Rwanda Astrophysics, Space and Climate Science Research Group (RASCSRG). Data and code: As discussed, the spectral series presented in this paper is composed of a series of different observing programmes at various telescopes. We request that any use of the data, whether compiled or re-reduced for this analysis, includes appropriate citation to the original papers. Spectra from the 6.5 m Magellan Telescopes (0.49 and 0.53 days) located at Las Campanas Observatory, Chile, can be found at: https://arxiv.org/abs/1710.05432. The 0.92 day spectra were obtained with the Australian National University (ANU) the 2.3 m telescope located at Siding Spring Observatory. The 1.17 day data obtained with the Southern African Large Telescope (SALT) under the Director’s Discretionary Time programme 2017–1-DDT-009, are available at https://ssda.saao.ac.za with the newly reduced spectra (e.g. with improved flux-calibration, see Appendix A) now available at https://github.com/Sneppen/Kilonova-analysis. X-shooter data from European Space Observatory (ESO) telescopes at the Paranal Observatory under programmes 099.D-0382 (principal investigator [PI]: E. Pian), 099.D-0622 (PI: P. D’Avanzo), 099.D-0376 (PI: S. J. Smartt), which are available at http://archive.eso.org and WISeREP (https://wiserep.weizmann.ac.il/). The re-reduced X-shooter spectra examining evolution in sub-epoch exposures at 1.4 days (derived in Sneppen et al. (2023a)) are available from https://github.com/Sneppen/Kilonova-analysis. Gemini-south data were obtained at the Gemini Observatory (Program IDs GS-2017B-Q-8 and GS-2017B-DD-4; PI: Chornock) and are together with SOAR spectra available at https://kilonova.org. HST observations were obtained using programs GO 14771 (PI: N. Tanvir) and GO 14804 (PI: A. Levan). We used the implementation of the P Cygni profile in the elementary supernova from https://github.com/unoebauer/public-astro-tools with generalisations to account for reverberation, time-delay effects, and special relativistic corrections (see https://github.com/Sneppen/Kilonova-analysis).

References

Appendix A: Appendix 1.1. SALT Reduction

As reported in Buckley et al. (2018), the SALT spectrum was observed in the evening twilight of August 18, 2017, with the Robert Stobie Spectrograph (Burgh et al. 2003) in long-slit mode, using a 2″ slit-width with a low-resolution PG300 grating giving access to the whole visible wavelength range, with a mean R ≈ 380 spectral resolution. For this work, we used the originally-reduced data and 1D extraction, and used a wavelength range of approximately 380 to 900 nm.

However, given the criticality of the overall spectral shape in this work, we revisited the flux calibration. First of all, we did this independently, and several times, using different parameter choices for the functions and orders, when fitting the spectrophotometric standard star, EG21, observed during the same night. Moreover, to check the robustness of the results, we also used different spectrophotometric standards, observed with exactly the same instrumental configuration as searched and found from the SALT database. We tested the calibration results using ones found closest to the date of observation: CD-32d9927, observed on 28 August, LTT4364 from 8 July, and LTT6248 from 10 September. Our tests showed that all the various independent choices using the original and most appropriate star, EG21, result in the same overall shape, with only very minor differences irrelevant to the final results. As a further demonstration of confidence, while the calibrations done with CD-32d992, LTT4364, and LTT6248 do show some variations in spectral shape¹, corresponding to 3%-level differences in the best-fit temperature, the Sr II absorption break location is similar in all cases beginning abruptly around the same wavelength, λ ∼ 640 − 680 nm. Thus, even the variation between the standard stars does not impact the break, and we conclude that the shape of the SALT spectrum is fully robust for the analysis undertaken in this work. The most appropriate spectrum, that calibrated with EG21, is displayed in all the figures and fits (see e.g. Fig 4).

Finally, we note that using calibrations with LTT4364 and LTT6248, i.e. the ones furthest away from the date of observation, the depth of the absorption feature will become somewhat smaller, translating to a difference in the feature’s optical depth and by extension in the abundance required to fit the feature (see Fig. A.1). However, this is of minor concern to our analysis, as exact mass-estimates already have sizeable systematic uncertainties and as we (for testing line-identifications and constraining the ionisation-transition) are mainly concerned with whether the feature has emerged at this time or not.

Fig. A.1. The spectral shape of the SALT spectrum at 1.17 days using various standard stars for flux calibration. While the overall spectral shape is only weakly affected by the choice of standard star, the inferred depth of the absorption feature is smaller using spectrophotometric standards more than 1–2 weeks from the date of observation.

All Tables

All Figures

Fig. 1. Spectral series of the kilonova AT2017gfo showing the evolution over the first 10 days post-merger. This spectral compilation of data includes the early spectra taken with the Magellan telescopes (Shappee et al. 2017), ANU 2.3 m telescope (Andreoni et al. 2017) and SALT (Buckley et al. 2018). From 1.4 days post-merger, the X-shooter spectrograph at the European Southern Observatory’s VLT detailed the UV to NIR structure with a daily cadence (Pian et al. 2017; Smartt et al. 2017). To constrain the rapid evolution at early times, we also include the 1.5 day Gemini FLAMINGOS-2 and SOAR spectra (Nicholl et al. 2017; Chornock et al. 2017). Towards intermediate and late-times, HST (Tanvir et al. 2017) infrared spectroscopy additionally constrains the spectral structure in the telluric region around 1400 nm. The spectral continuum is well described as a blackbody in all photospheric epochs and even into nebular-phase epochs with a blackbody temperature evolution that is detailed in Sect. 4. In the zoom-in sub-panels, individual spectral lines are highlighted at different epochs, including (1) the 760 nm (2) 1.0 μm, and (3) 1.4 μm P Cygni features, and (4) the NIR line-emission features.

In the text

Fig. 2. Average ionisation inferred for the r-process elements in groups I–III. These elements are likely to produce the strongest lines under LTE conditions. We show them as a function of the observed time post-merger. The ionisation states were calculated using the Saha equation, where we assume delectron density, ne(t) = 108 cm−3(t/1.5 days)−3 and the blackbody temperature evolution found for AT2017gfo (Sneppen et al. 2023b). The vertical blue and red vertical regions indicate the observations with X-shooter and other instruments (see the spectral series in Fig. 1). At certain temporal intervals, particularly at around 1 day post-merger, the r-process elements that dominate the opacity transition rapidly between ionisation states. For t ⪆ 5 days post-merger, NLTE-effects are likely to become important, implying that neutral, single-, and double-ionised species of these elements may be present (Pognan et al. 2022).

In the text

	Fig. 3. Spectral series of the kilonova AT2017gfo. The best-fit blackbodies are overlaid. The temporal span here covers the photospheric epochs (i.e. spectra up to 5.4 days post-merger), where a blackbody model approximates the continuum well. The blackbody temperatures and the spectral features are listed in Table 1 and in Sect. 5, respectively.
In the text

Fig. 4. Emergence of the 1 μm P Cygni feature in time. The selected spectra are from 0.48 days (LDSS, Shappee et al. 2017), 0.92 days (ANU, Andreoni et al. 2017), 1.17 days (SALT, Buckley et al. 2018), 1.43 days (X-shooter, Pian et al. 2017), 1.47 days (SOAR, Nicholl et al. 2017 and Gemini-South, Chornock et al. 2017), and 2.42 days (X-shooter, Smartt et al. 2017). The spectra before 1 day do not deviate strongly from the blackbody continuum. However, between 0.92 and 1.17 days, strong absorption appears, while the bulk emission begins to emerge with a delay of around 10 hours. This is the predicted emergence time-frame for both emission and absorption (and the early spectral shape) for a feature produced by LTE Sr II due to a recombination wave passing through the ejecta (Sneppen et al. 2023a), and it is inconsistent with a He Iλ1083.0 nm interpretation (Sneppen et al. 2024).

In the text

Fig. 5. HST spectrum of AT2017gfo at 4.9 days post-merger from Tanvir et al. (2017). The grey shaded region indicates the continuum constraints from the interpolation of the X-shooter spectra taken at 4.4 and 5.4 days post-merger. Regardless of the exact continuum modelling, a P Cygni absorption/emission feature around 1400 nm is clearly visible. Domoto et al. (2022) suggested that this feature could be due to La III lines (here indicated with dashed blue lines). The characteristic velocity of this feature is relatively low, vph = 0.07 ± 0.01c, but may be systematically biased by the strong nearby features (including the prominent reverberation emission peak around 1.1 μm).

In the text

Fig. 6. Ionisation degree of Sr in LTE as a function of temperature and electron density. We overlay the approximate locations in parameter landscape of AT2017gfo with time. The lines indicate the evolution modelled as constraints from reverberation effects within a single epoch. The temperature error bars are for the part of the ejecta with line-of-sight expansion (Sneppen 2023), and the density assumes a spherically symmetric, homogeneous, and homologously expanding high-velocity (0.2c ≤ v ≤ 0.3c) ejecta component of 0.01 M⊙. The electron density can only be estimated to an order of magnitude, but regardless of the exact ne, a rapid step change in ionisation state occurs over a narrow interval in temperature around 4500 K. The 1 μm feature is not detected in observations prior to 1 day (marked in red), but from 1.17 days and forward, the feature is present (marked in blue).

In the text

Fig. 7. Fraction of Sr II under LTE conditions as a function of time post-merger. On the upper x-axis, we illustrate the physical time of scattering probed by the 1 μm feature in absorption and emission by the various spectra from ANU, SALT, VLT/X-shooter, and Gemini FLAMINGOS-2. The dotted lines indicate that no feature has clearly emerged (prior to ∼1.4 days in both emission and absorption). We indicate the required fraction of Sr II needed to produce the optical depths observed in each spectrum (assuming MSr = 6 × 10−6 M⊙ homogeneously distributed across 0.2c < v < 0.4c; see main text). The ANU (X-shooter) and SALT (Gemini) spectra indicate a rapid transition in absorption (emission) from an optically thin to optically thick regime around 1.3–1.5 days. For the feature to emerge around 1.4 days post-merger, the ionisation temperature is required to be nearly identical to the radiation temperature inferred from the blackbody emission. For the feature to emerge at near-contemporaneous physical time in the polar and equatorial ejecta, a near-isotropic temperature is required.

In the text

Fig. 8. Degeneracy between the blackbody temperature and the inferred emitting area (or velocity) for a given NIR continuum. The X-shooter epoch 4 spectrum is overlaid with blackbodies of various temperatures and areas that produce similar NIR flux levels. The observed NIR features can be modelled as P Cygni or pure emission components, but make the exact underlying continuum difficult to estimate precisely. If the near-observed infrared features are produced by La III and Ce III for LTE conditions, these species must reside in relatively hot ejecta. Hotter ejecta produces a higher NIR continuum flux, which would conversely require a smaller emitting area or lower velocity to match the observed continuum. While the optical P Cygni features are consistent with the observed blackbody temperature and velocity, these NIR identifications would need higher temperatures and be spatially distinct from their optical counterparts.

In the text

	Fig. A.1. The spectral shape of the SALT spectrum at 1.17 days using various standard stars for flux calibration. While the overall spectral shape is only weakly affected by the choice of standard star, the inferred depth of the absorption feature is smaller using spectrophotometric standards more than 1–2 weeks from the date of observation.
In the text