© The Authors 2023

1. Introduction

Neutron star mergers are a significant source of r-process elements, as suggested by Lattimer et al. (1977) and Symbalisty & Schramm (1982) and indicated by observations of kilonovae (e.g., Tanvir et al. 2013; Coulter et al. 2017), and with confirmation provided by spectroscopy (Watson et al. 2019).

The most detailed constraints on any kilonova so far are provided by the series of X-shooter spectra for AT2017gfo (Smartt et al. 2017; Pian et al. 2017). These daily spectra taken from 1.4 days and onwards after the event provide the temporal evolution of the transient in continuum, emission, and absorption features from the ultraviolet to the near-infrared. It was these data that yielded the first observational identification of a freshly synthesised r-process element in an astrophysical site (Watson et al. 2019). This identification came in the form of resonance lines of Sr⁺ and demonstrated that the early spectra of AT2017gfo can be explained by lighter r-process material. These findings were reproduced and extended with systematic analyses using radiative transfer codes with improved atomic line data (Domoto et al. 2021; Gillanders et al. 2022; Vieira et al. 2023). Another proposed interpretation is that this P Cygni feature originates from the helium 1083.3 nm line, which, as shown in Perego et al. (2022) and Tarumi et al. (2023), in non-local-thermodynamic-equilibrium (NLTE) conditions could reproduce the observed feature. In addition to Sr⁺, the potential presence of Zr⁺ (Gillanders et al. 2022), Y⁺, and Zr⁺ (Vieira et al. 2023) has been inferred at days 1–2 post-merger using radiative transfer modelling, principally based on absorption features at wavelengths λ ≤ 500 nm. The identification of features at λ ≈ 1200–1400 nm with absorption due to La²⁺ and Ce²⁺ has been proposed (Domoto et al. 2022). On the other hand, Watson et al. (2019) noted several near-infrared emission features that emerge most strongly several days post-merger, without suggesting an identification. These works hint at a wealth of information still concealed in these spectra.

Modelling of the X-shooter data has also revealed the spherical nature of the geometry of the early kilonova ejecta, contrary to expectations from most hydrodynamic models (Sneppen et al. 2023b). Here, the fact that the 1 μm feature is a P Cygni was a key discovery, as it allowed the geometry of the kilonova outflow to be assessed.

The identification of another P Cygni feature is therefore important as it can be used to verify the 1 μm line and, because of the anticipated extreme wavelength-dependence of the opacity in r-process-dominated ejecta, as a further probe of the geometry and abundance stratification of the kilonova. In this letter, we revisit these X-shooter spectra and identify a new spectral component, a P Cygni feature at 760 nm rest wavelength, which emerges most clearly several days post-merger. We investigate the wavelength stratification of the line and show that it provides an independent estimate of the expansion velocity of the ejecta, which is in good agreement with the constraints from the 1 μm feature. We also investigate whether this feature can be explained by the expected lines of yttrium at the temperature of the continuum fit and find that 4d²–4d5p transitions of Y⁺ match the feature well. That is, in terms of central wavelength, velocity stratification, and the weaker prominence of the feature compared to its stronger Sr⁺ counterpart, this feature is well-described as a P Cygni profile caused by transition lines from singly ionised yttrium. Furthermore, the fact that the P Cygni feature emerges 3–4 days post-merger is a property of the yttrium identification, because at velocities of 0.2c − 0.3c, the combined interference of different lines of Y⁺ conceals the characteristic P Cygni spectral shape.

2. The 760 nm line detection

For this analysis, we use the X-shooter spectra, which were taken daily beginning 1.4 days after the merger, and follow the data reduction presented in Watson et al. (2019). As illustrated in Fig. 1, the optical spectra from epochs 3–6 (3.4–6.4 days post-merger) contain a clear excess of emission around 760 nm with a corresponding absorption feature at about 650 nm. The overall shape of the spectral component is reminiscent of P Cygni features, which are characteristic of expanding envelopes where the same spectral line yields both an emission peak near the rest wavelength and a blueshifted absorption feature. While the feature is also there in epochs 3 and 6, we focus in this analysis on epochs 4 and 5, where the feature stands out most clearly and its properties can be most readily measured.

Fig. 1. X-shooter spectra for epochs 3–6 (3.4–6.4 days post-merger) with blackbody fit and emphasis on the P Cygni profiles of Sr+ (blue) and an unidentified spectral line at 760 nm (red). For clarity, and to show them on the same plot, the spectra are offset as indicated by the dotted black lines. The feature is clearly noticeable from epoch 3 and onwards. However, the line is poorly constrained in epoch 3 as it less prominent relative to the continuum and overlaps with the most blueshifted parts of the Sr+ P Cygni. It is also poorly constrained in epoch 6 and later as the continuum is noisy and less well fit by a blackbody. Epochs 4 and 5 therefore provide the best constraints on the line properties.

2.1. Statistical significance

To determine the significance of the feature, we compare the χ² with and without fitting this spectral feature. For the null model, we fit the continuum as a blackbody as empirically proposed in many previous works (e.g., Malesani et al. 2017; Drout et al. 2017; Shappee et al. 2017; Waxman et al. 2018; Watson et al. 2019) and as theoretically motivated in Sneppen (2023). However, to account for the effects of spectral features on the continuum fit, we include parameters for the observed Gaussian emission lines at 1.5 μm and 2.0 μm and the Sr⁺ P Cygni profile following the prescription in Watson et al. (2019). The analysis is robust regardless of whether these spectral features are parameterised in the fit as we do here or the analysis is limited to wavelengths without these spectral components (Sneppen et al. 2023b). To complete the full model, in addition to continuum and nuisance spectral components, we add a P Cygni prescription to model the 760 nm line.

For this analysis, we use the implementation of the P Cygni profile in the Elementary Supernova Model¹, where the profile is set by several properties of the kilonova atmosphere. The profile is expressed in terms of the rest wavelength, λ₀, and the line optical depth, τ, with a velocity stratification parameterised by a scaling velocity, v_e, a photospheric velocity, v_phot, and a maximum ejecta velocity, v_max. This parameterisation assumes an optical depth with an exponential decay in velocity as presented in Thomas et al. (2011), but one can equivalently assume a power-law decay with the resulting constraints being indistinguishable. In short, the optical depth determines the strength of absorption and emission, while the velocity of the ejecta sets the wavelength stratification.

For epochs 4 and 5, no blackbody continuum model can fit the observed profile well from 600 to 800 nm (as seen in Fig. 2), with the reduced chi-squared of ${\chi }_{\nu ,\text{4th}}^2=1.65$ and ${\chi }_{\nu ,\text{5th}}^2=1.88$. In contrast, fitting the residual structure with the additional P Cygni model instead yields ${\chi }_{\nu ,\text{4th}}^2=1.17$ and ${\chi }_{\nu ,\text{5th}}^2=1.21$ over the same wavelengths (see the posterior landscape in Figs. A.1 and A.2). The fact that ${\chi }_{\nu }^2>1$ for the full model may indicate that the errors are underestimated, that a single P Cygni prescription is an over-simplification, or that there are unknown spectral substructures.

Fig. 2. Spectra for epochs 4 and 5 divided by the model continuum flux, with single-line P Cygni profiles fit around 760 nm overlaid. For both epochs, the existence of the feature is highly significant.

Nevertheless, the substantial improvement in χ² while only requiring the five additional parameters of the P Cygni suggests that the feature is highly significant and we investigate this further here. To examine the significance using more observationally related uncertainties, we artificially inflate the errors of the full model so that ${\chi }_{\nu }^2=1$. In this case, the simpler continuum model without the 760 nm P Cygni provides a poorer fit to the data, with χ² values going from 14 321 to 10 102 and 15 686 to 10 096 for epochs 4 and 5, respectively. This results in Δχ² = 4219 and 5590 for five additional degrees of freedom. This corresponds to a null-hypothesis probability rejected at 26σ and 33σ, again for epochs 4 and 5, respectively. One could also apply the deviance information criterion (DIC; see e.g., Liddle 2007) here, which results in a ΔDIC < −4000 for both epochs corresponding to decisive evidence on the Jeffreys scale (e.g. Jeffreys 1939; Kass & Raftery 1995).

2.2. The central wavelength

The best-fit central wavelength of the single P Cygni line fit is 750.5 ± 1 nm and 769.1 ± 1.2 nm for epochs 4 and 5, respectively. The unweighted average of these would indicate a central wavelength of roughly 760 ± 10 nm. While the statistical uncertainties of fitted parameters from each epochs are very small and is even at a level where the central wavelength value is inconsistent between the two epochs, we emphasise that there are systematic effects that are expected to bias the central value significantly in one or both of the epochs. First, without a prior spectral identification, we assume the origin of the line to be a single transition, when a blending of nearby lines may produce a similar profile – indeed, we think multiple transitions from related energy levels are highly likely (see Sect. 3). Second, the most redshifted part of the 760 nm feature is coincident in wavelength with the most blueshifted absorption from Sr⁺, and so the exact location of the emission peak and the maximum ejecta velocity of the Sr line-forming region are degenerate. Third, the relativistic velocity of the ejecta moving at ≈0.15c suggest that there are significant time delays, which implies the different wavelengths of the profile may be probing different times, that is, at different temperatures and optical depths.

2.3. Velocity stratification

In addition to the central wavelength, the P Cygni profile also provides information on the properties of the expanding ejecta. The spherical nature of the early ejecta of AT2017gfo was first shown by Sneppen et al. (2023b) from the high degree of consistency between the line-of-sight velocity, v_∥, from the Sr⁺ lines and the cross-sectional velocity, v_⊥, derived from the continuum flux. Here, the cross-sectional velocity (or equivalently area) is also determined by comparing the observed flux with the intrinsic blackbody luminosity, while assuming Planck cosmology (Planck Collaboration VI 2020) and the Hubble flow recession velocity of the host galaxy in Mukherjee et al. (2021). In addition, the line fits for Sr⁺ suggest the photospheric velocity is v_∥(Sr⁺)≈0.17c and 0.15c in epochs 4 and 5, respectively (Sneppen et al. 2023b). This allows the degree of sphericity to be parameterised in terms of the asymmetry index:

$$\begin{array}{c}\hfill \mathrm{{\rm Y}}({\text{Sr}}^{+})=\frac{v_{\perp }-v_{\Vert }({\text{Sr}}^{+})}{v_{\perp }+v_{\Vert }({\text{Sr}}^{+})}.\end{array}$$(1)

Similarly, the new line profile also contains spectral constraints on the velocity stratification of the expanding ejecta, which is shown in the posterior probability distributions of Fig. 3. As seen in Table 1, the photospheric velocity derived from the line, v_∥(760 nm) ≈ 0.16c, is larger than the cross-sectional velocity of the continuum and broadly similar to that derived from the Sr⁺ lines. Therefore, this new line seems to confirm an increasing discrepancy in later epochs between the velocity indicators with v_⊥ < v_∥ as was hinted at by the Sr⁺ P Cygni analysis (Sneppen et al. 2023b). However, this does not necessarily imply an increasingly prolate geometry; for example, the blackbody assumption is a worsening approximation in later epochs, resulting in a substantial increase in systematic uncertainty.

Fig. 3. Corner plots showing the posterior probability distribution of key parameters for epochs 4 (left) and 5 (right) fitting as a P Cygni feature from the set of Y+ lines from the 4d5p − 4d2 transitions with λ ≥ 700 nm. The fit parameters are (1) the cross-sectional velocity, v⊥, derived from the blackbody normalisation fit assuming Planck cosmology (Planck Collaboration VI 2020), (2) the photospheric velocity, vphot, (3) the maximum outer ejecta velocity, vmax, (4) the velocity range, Δv = vmax − vphot, and (5) the average optical depth of the Y+ lines weighted by their relative strength, τavg. Velocities are given in units of the speed of light, c. The central wavelength of a line, λ0, and its relative strength is set by the 4d5p to 4d2 transitions in Y+ (see Biémont et al. 2011 and Table 2) in LTE at about 0.4 eV. These constraints show only the statistical uncertainty in the fitting model.

We note that the fits to the line for epoch 5 appear saturated. However, there is a strong degeneracy between the optical depth and the thickness of the line-forming region (e.g., Fig. 3). Therefore, the apparent saturation could be caused by a thin line-forming region in the fit (i.e., a small difference between v_phot and v_max) producing a weak feature, which in the best fit is compensated by increasing the optical depth. However, it is likely that line blending, reverberation effects, or other systematic effects may change the velocity stratification and, by extension, the optical depth. Additionally, we note a high-opacity line is problematic for any direct comparisons between the photospheric velocity and the cross-sectional velocity, as a high opacity may detach the line-forming region from the continuum, inducing a bias in the inferred photospheric velocity (Sim 2017).

3. Yttrium identification

The new P Cygni line detection is useful, not only because it allows comparison with the velocity analysis of Sr⁺, but also because it provides a tantalising potential identification, which may subsequently improve the constraint on the composition and r-process synthesis of kilonovae. However, any identification is difficult, as current atomic line lists are incomplete and/or suffer from significant inaccuracies in wavelength and/or line strength (e.g., VALD, Kurucz; Ryabchikova et al. 2015; Kurucz 2018, respectively), especially for elements beyond the first r-process peak. Due to the difficulty in generating this information for heavier elements, despite significant recent efforts (e.g., Tanaka et al. 2020), line lists with reliable wavelengths and transition strengths remain mostly incomplete beyond the first peak, and especially longward of 1 μm, which hampers modelling efforts.

We examined the line lists of the most obvious spectral candidates, that is, those with strong transitions due to a small number of valence electrons and low-lying energy levels (Watson et al. 2019; Domoto et al. 2022). First, the elements near the first r-process peak, which are abundant and contain the only known prior identification, namely of Sr. Of these elements, neither Ru, Sr, or Zr have strong lines at the relevant wavelengths for their likely dominant ionisation states. There are several noteworthy strong lines from r-process second peak elements (notably Ba⁺ and La⁺), but these are located at considerably shorter (or longer) wavelengths.

The most promising candidate, yttrium, should – at these temperatures and densities – principally be in the singly ionised state, Y⁺, which has strong line transitions from 4d5p to 4d². These transitions have an LTE-weighted mean wavelength of around 760–770 nm with the most prominent line being at 788.19 nm (see Table 2 and Fig. 4). Given the strong transition, its relatively high expected abundance, and its atomic number proximity to Sr, yttrium seems a likely candidate. Interestingly, a specific prediction of the yttrium identification is that the feature should first emerge clearly 3–4 days post merger, as is indeed observed in the spectra. This is because in earlier epochs with greater characteristic velocities, the absorption from the lines around 760 nm will coincide in wavelength with the emission of the 670 nm Y⁺ lines, and act to suppress the prominence of the feature, as discussed further in Sect. 4.

Fig. 4. Grotrian diagram showing the lowest energy levels in Y+. The transitions between energy levels are indicated, with the line transparency set by the transition strength, log gf. The 4d5p–4d2 (28 000 cm−1–14 000 cm−1) transitions are shown with thick lines in blue. Radiative transfer modelling analysis in Vieira et al. (2023) tentatively links the 4d5p–4d2 (28 000 cm−1–8000 cm−1) transitions to absorption at λ < 450 nm in the epoch 1 spectrum; other notable transitions are around 300 nm, too short to be observed by X-shooter (4d5p–4d5s, 28 000 cm−1–1000 cm−1), and signatures of the 5s5p–4d2 transitions (24 000 cm−1–14 000 cm−1) at 1 μm, which are much weaker than the Sr+ line transitions at similar wavelengths. The 4d5p–4d2 transitions have lower oscillator strength than several of the Y+ lines at shorter wavelengths, but due to the rapid decline of opacity with wavelength for typical r-process elements, these transitions have higher opacity compared to the continuum at 760 nm. When no lines are drawn between levels, the transition is negligible with log gf < −2.

The ionisation energies for Y and Y⁺ are 6.2 and 12.2 eV, respectively. Using the Saha equation and the range of electron densities expected for the ejecta (n_e ≈ 10⁶ − 10⁸ cm³) suggests characteristic temperatures for the first and second ionisation of 2200–2400 K and 4100–4700 K, respectively. For the observed temperature in the analysed epochs, T ≈ 3000 K, the fraction of yttrium in the singly ionised state should be near unity. We note that for earlier epochs (for instance 1.4 days post-merger), the inferred temperature is near the second ionisation threshold, which suggests a substantial fraction of yttrium at those epochs may be present as Y⁺⁺. The growing prominence of this feature over the epochs could in part be the result of the decreasing ionisation of yttrium in LTE. However, given the presence of a strong Sr⁺ absorption in the first epoch spectrum (Watson et al. 2019), and the lower first and second ionisation energies of Sr, decreasing ionisation alone is not a very satisfactory explanation for the increased prominence of the 760 nm feature at epochs 3–6.

3.1. P Cygni analysis of the yttrium feature

In Fig. 5, we show the continuum and P Cygni fit for the Y⁺ line transitions with λ ≥ 700 nm. We note that these transitions provide a good fit to the observed absorption trough and emission peak. This is despite the fact that, for this fit, we fixed the central wavelength and the relative strengths of the lines according to an LTE calculation based on the atomic line data of Biémont et al. (2011). The multi-line nature of the fit even allows the spectral substructures to be fit, that is, tracing both the overall line shape but also the small-scale wiggles observed, for example, in the epoch 5 fit. Interestingly, this shows that characteristic wiggles may be produced by the combination of nearby lines with varying transition strengths, which in turn suggests that spectral fingerprints of specific r-process elements could be hiding within other such wiggles of the data. We note that including shorter wavelength transitions generally leads to smaller fitted velocities (as these transitions naturally extend the profile to slightly shorter wavelengths), but the exact lower wavelength cutoff below which lines are not considered in the fit is unimportant in terms of the goodness-of-fit. Without substantially changing the best-fit profile-shape, we can include all 4d5p − 4d² transitions, but as the ∼670 nm lines are not the driving lines behind the observed feature and there are strong indications that these shorter wavelength features may be suppressed by line blanketing (see Sect. 3.2), for the presented fit we emphasise the dominant higher-wavelength transitions. In Sect. 3.2, we use the full line list and explore the relation between the ∼670 nm and ∼760 nm transitions.

Fig. 5. X-shooter spectra with model fit including Y+ lines with λ ≥ 700 nm for epochs 4 and 5. The blackbody continuum model is shown as a dashed black line, the blackbody plus the Sr+ P Cygni lines as a dotted purple line, and the full model including both Sr+ and Y+ as a solid red line. The partially transparent purple and red regions indicate the successive contributions to these models of the Sr+ and Y+ P Cygni lines, respectively. The grey hatched region indicates the overlapping noisy regions between the UVB and VIS arms of the spectrograph. The Y+ lines (from the 4d2–4d5p transitions) have fixed central wavelengths and their relative strengths are set by the LTE condition.

The column density inferred from fitting the P Cygni Y⁺ 4d5p–4d² transition fits suggests total yttrium masses of 0.6 × 10⁻⁴ M_⊙ and 2.4 × 10⁻⁴ M_⊙ would be required to produce these features for epochs 4 and 5, respectively, assuming spherical symmetry. The derived masses should be treated with caution as there is no correction for light-travel-time effects or additional blending lines, which could also potentially bias inferred properties. These are lower limits on the total amount of material ejected, as they trace only the matter between the photospheric front and the outer atmosphere. This is especially apparent for the best-fit in epoch 5, where the line profile is saturated and the line-forming region may therefore be detached from the inner emitting mass. Given the total mass within the atmosphere ejected inferred from light-curve modelling, M_ej ≈ 3 − 5 × 10⁻² M_⊙ (Kasen et al. 2017), this would suggest a lower bound on the mass fraction of yttrium of around 0.5%–0.8% from epoch 5. This is within the yield estimated from numerical simulations of neutron-star mergers (Perego et al. 2021), and is consistent with estimates of the r-process fraction in the Sun (Bisterzo et al. 2014; Lodders et al. 2009).

3.2. TARDIS analysis of the yttrium feature

In addition to the P Cygni fitting framework, we also test the prominence of features due to the proposed Y⁺ transitions using the Monte Carlo radiative-transfer spectral synthesis code TARDIS (Kerzendorf & Sim 2014). In Fig. 6, we show the spectra resulting from a TARDIS simulation with a solar r-process element abundance distribution for the r-process elements below the second peak (34 ≤ Z ≤ 50). Our solar r-process distribution is based on the solar abundance (Lodders et al. 2009) with the s-process abundance distribution subtracted (Bisterzo et al. 2014). In the same figure, we also show the spectrum resulting from a simulation with the same parameters but without yttrium, as well as the difference between the models with and without yttrium. When yttrium is included, a P Cygni component is visible at 760 nm. A similar plot, but this time including all r-process elements (34 ≤ Z ≤ 92), is also shown in Fig. 6. In both scenarios, that is, with all elements or with only the lighter ones, the P Cygni feature at about 760 nm is still produced by the 4d²–4d5p transitions of Y⁺.

Fig. 6. Example TARDIS model spectra including yttrium (blue line) and excluding yttrium (red dashed line) but otherwise composed of r-process elements with solar r-process abundance ratios (Lodders et al. 2009; Bisterzo et al. 2014). The left panel is limited to the first r-process peak elements (34 ≤ Z ≤ 50), while the right panel shows a model with all r-process elements (33 ≤ Z ≤ 92). The spectrum at the photosphere is modelled as a blackbody, which is then propagated through a relativistically expanding atmosphere (0.15c ≤ v ≤ 0.2c) using the radiative transfer equations with a radially uniform composition using the Kurucz (2018) line list. The red region indicates the 760 nm feature, which is clearly associated with the presence of yttrium. The yttrium abundance in the model constitutes 0.8% of the total mass. The 0.8–1 μm feature is the Sr+ P Cygni feature (Watson et al. 2019).

In addition to the ∼760 nm feature that we observe, the presence of yttrium also results in (1) strong absorption of flux at short wavelengths ≤450 nm, (2) a small P Cygni feature at about 650 nm in emission (the presence or absence of which is hard to constrain observationally from the X-shooter spectra as the absorption sits in the noisy region between UVB and VIS arms of the spectrograph), and (3) a near negligible contribution to the 1 μm feature from the 5s5p–4d² transitions. The relative strengths of these features are affected by the opacity from other lines in the same region of the spectrum. This is seen most clearly in the difference between the left and right panels of Fig. 6, where the more blueward features are significantly less prominent with the higher blue-wavelength opacity introduced to the radiative transfer model by the r-process second peak elements. The relative prominence of the ∼760 nm feature compared to its lower-wavelength counterparts could therefore be considered tentative evidence of a significantly higher opacity below ∼700 nm, which is consistent with the arguments in Gillanders et al. (2022) that for synthetic models to match the observed flux distribution for λ < 750 nm in later epochs requires the inclusion of lanthanides. However, as we cannot confidently exclude a feature at λ ∼ 600 nm, and because the line lists we currently have are incomplete – a more complete line list could well result in greater blue-wavelength opacity from the lighter elements (with 34 ≤ Z ≤ 50) –, we do not regard the dominance of the 760 nm feature as definitive evidence of lanthanides; nevertheless, in agreement with Gillanders et al. (2022), we do find it suggestive. Regardless, it is interesting to note the unique prominence of the P Cygni feature from Y⁺ in these models as well as the relative strength compared with the Sr⁺ feature. We emphasise that the potential identification of a 760 nm P Cygni feature from Y⁺ noted here is complementary to, and independent of, a deficiency of flux noted at short wavelengths in early epochs, which has been tentatively suggested to be due to Y⁺ and, to a lesser degree, to Zr⁺ (Vieira et al. 2023). We note that Vieira et al. (2023) explored the spectral signatures of Y⁺ and did not find this feature. However, the 4d5p to 4d² transitions were mistakenly omitted from the compiled line list in this latter analysis (N. Vieira, priv. comm.).

4. Discussion

The 760 nm P Cygni feature is less prominent than its strontium P Cygni counterpart, but shows a similar width. Furthermore, we show that it is well-fit by 4d5p to 4d² transitions in Y⁺. If this identification is genuine, it raises several prospects in terms of r-process nucleosynthesis, velocity constraints, and the geometry of the ejecta.

First, the two elements identified so far, strontium and yttrium, are both near the first r-process-peak, which is consistent with analyses of the observed early colours of AT2017gfo (Waxman et al. 2018) and arguments in Gillanders et al. (2022) and Vieira et al. (2023), which allow a negligible to modest lanthanide abundance to reproduce the observed optical and UV flux. However, the lack of clear spectral evidence of heavier r-process elements is somewhat surprising given the expected abundances produced in binary neutron-star mergers from typical numerical simulations (e.g., Ekanger et al. 2023). This is especially notable considering the variations in abundances from different mass-ejection channels, which are believed to contribute and dominate in different epochs and depths. If both the Sr and Y line identifications are genuine, these findings reinforce the hypothesis of a light r-process dominance despite a large radial change across epochs, which would allow for the possibility that light r-process elements dominated the nucleosynthetic output of AT2017gfo.

Second, the 760 nm P Cygni feature is about an order of magnitude weaker than the strontium feature in emission, though this difference is much less in absorption, while it weakens more slowly than the strontium feature at 5–6 days post-merger compared to its strontium counterpart. These differences are surprising given that Sr and Y are immediate neighbours in the periodic table and their abundance and radial distribution should be closely linked. Furthermore, the yttrium and strontium ionisation states do not vary significantly over the temperatures inferred for these epochs, with the vast majority being in the singly ionised state. However, given the significant metastable states involved in these transitions for both Sr and Y, it seems likely that NLTE effects change the Sr⁺ and Y⁺ populations over time as the ejecta expands. Such effects would need to be modelled accurately in order to infer the relative Sr to Y ratio.

Importantly, a specific prediction of the yttrium identification is that the 760 nm feature should only clearly emerge when the velocity drops to about 0.15–0.2c, which is several days after the merger, which is precisely what observations reveal. This is not because yttrium is not in the outer layers, but because the characteristic velocities are so large that the absorption of the 760 nm P Cygni is coincident in wavelength with the emission peak predominantly due to the 661.4 nm line from the 4d5p a³P₂–4d²z³D₃° transition. At velocities of 0.3c, these features are strongly broadened and the characteristic P Cygni feature is effectively cancelled by the summed contribution of these transitions. We illustrate this effect using a TARDIS model for three different characteristic ejecta velocities; see Fig. 7. On a related note, the strong flux deficit shortward of 450 nm in the first epoch noted in Watson et al. (2019) – and suggested there to be related to Sr⁺ – is mostly ascribed to Y⁺ lines in radiative transfer modelling (Gillanders et al. 2022; Domoto et al. 2022; Vieira et al. 2023). A feature at ≈760 nm in the first epochs could bias the inferred ejecta properties from the Sr⁺ lines. However, as argued in Sneppen et al. (2023b), this is unlikely to dramatically shift the constraints given the relative weakness of such features.

Fig. 7. Example TARDIS model spectra including yttrium (blue) and excluding yttrium (red) but otherwise composed of all r-process elements (33 ≤ Z ≤ 92) with solar r-process abundance ratios (Lodders et al. 2009; Bisterzo et al. 2014) for three different ejecta velocities. The three bottom panels show the difference in the spectral shape between the yttrium-free and yttrium-containing model for each characteristic velocity. Only when the characteristic ejecta velocity has decreased below ∼0.2c (at about 3–4 days post merger) does the P Cygni feature become apparent, because it no longer merges with lower-wavelength yttrium features.

Lastly, we emphasise that while the statistical uncertainties are small for any specific fitting model, there are larger systematic uncertainties. This is illustrated by the relatively large variations in v_phot or λ₀ for different fitting models and epochs, respectively. For the constraints on the line shape, these systematic uncertainties may include reverberation effects, potential unmodelled blending lines, the simplicity of the assumed parameterisation of the optical depth, or even limitations in inferring the continuum flux. For the comparison with the cross-sectional velocity, v_⊥, the blackbody framework becomes less robust in later epochs due to emerging lines and the observed increased complexity of the spectrum. In contrast, the geometrical constraints in early epochs are derived from the more prominent Sr⁺ lines, which are less sensitive to the modelling of other features (Sneppen et al. 2023b), and the continuum is well-fit by a blackbody (Sneppen 2023). Therefore, while the Y⁺ feature provides an interesting approximate cross-check on the Sr⁺ constraints, we caution that the statistical uncertainties derived here from this 760 nm feature do not represent the full uncertainty.

5. Conclusion

We analysed the X-shooter spectra of the kilonova AT2017gfo associated with GW170817, focussing on epochs 4 and 5, that is, 4.4 and 5.4 days post-merger, respectively. We find strong evidence for a broad P Cygni feature with the emission component centred around 760 nm. The velocity structure of this line is similar to that of the ∼1 μm line previously associated with Sr⁺ (Watson et al. 2019). We propose a possible identification of this feature as the sum of the 4d²-4d5p transitions of Y⁺ between 726.42 nm and 788.19 nm. We find this to be reasonable interpretation because a r-process abundance of yttrium yields a P Cygni feature with a central wavelength, a velocity stratification, a prominence, and a timing that match those of the observed feature. This is the second clear identification of a particular feature with an element in a kilonova and would strengthen the case for AT2017gfo being light r-process-element dominated.

Acknowledgments

We thank Stuart Sim and Nicholas Vieira for useful correspondence, discussions and feedback. The Cosmic Dawn Center is funded by the Danish National Research Foundation under grant number 140.

References

Appendix A: Additional figures

We provide the posterior distribution of the single P Cygni line fit in Figs. A.1 and A.2. We emphasise that these plots show only the statistical uncertainty in the fitting model and not the (potentially sizable) systematic effects, such as those from line blending and reverberation effects.

Fig. A.1. Corner plots showing the posterior probability distribution of key parameters for epoch 4 with a single P Cygni line fit. The fit parameters are (1) the cross-sectional velocity, v⊥, derived from the blackbody normalisation fit assuming Planck cosmology (Planck Collaboration VI 2020), and several parameters of the 760 nm P Cygni feature, including (2) the photospheric velocity, vphot, (3) the maximum outer ejecta velocity, vmax, (4) the velocity range, Δv = vmax − vphot, (5) the optical depth of the line, τ, and (6) the central wavelength, λ0, in μm. Velocities are indicated in units of the speed of light, c. These constraints show only the statistical uncertainty in the fitting model.

Fig. A.2. Corner plots showing the posterior probability distribution of key parameters for epoch 5 with a single P Cygni line fit. The fit parameters are (1) the cross-sectional velocity, v⊥, derived from the blackbody normalisation fit assuming Planck cosmology (Planck Collaboration VI 2020), and several parameters of the 760 nm P Cygni feature, including (2) the photospheric velocity, vphot, (3) the maximum outer ejecta velocity, vmax, (4) the velocity range, Δv = vmax − vphot, (5) the optical depth of the line, τ, and (6) the central wavelength, λ0, in μm. Velocities are indicated in units of the speed of light, c. These constraints show only the statistical uncertainty in the fitting model.

All Tables

All Figures

Fig. 1. X-shooter spectra for epochs 3–6 (3.4–6.4 days post-merger) with blackbody fit and emphasis on the P Cygni profiles of Sr+ (blue) and an unidentified spectral line at 760 nm (red). For clarity, and to show them on the same plot, the spectra are offset as indicated by the dotted black lines. The feature is clearly noticeable from epoch 3 and onwards. However, the line is poorly constrained in epoch 3 as it less prominent relative to the continuum and overlaps with the most blueshifted parts of the Sr+ P Cygni. It is also poorly constrained in epoch 6 and later as the continuum is noisy and less well fit by a blackbody. Epochs 4 and 5 therefore provide the best constraints on the line properties.

In the text

	Fig. 2. Spectra for epochs 4 and 5 divided by the model continuum flux, with single-line P Cygni profiles fit around 760 nm overlaid. For both epochs, the existence of the feature is highly significant.
In the text

Fig. 3. Corner plots showing the posterior probability distribution of key parameters for epochs 4 (left) and 5 (right) fitting as a P Cygni feature from the set of Y+ lines from the 4d5p − 4d2 transitions with λ ≥ 700 nm. The fit parameters are (1) the cross-sectional velocity, v⊥, derived from the blackbody normalisation fit assuming Planck cosmology (Planck Collaboration VI 2020), (2) the photospheric velocity, vphot, (3) the maximum outer ejecta velocity, vmax, (4) the velocity range, Δv = vmax − vphot, and (5) the average optical depth of the Y+ lines weighted by their relative strength, τavg. Velocities are given in units of the speed of light, c. The central wavelength of a line, λ0, and its relative strength is set by the 4d5p to 4d2 transitions in Y+ (see Biémont et al. 2011 and Table 2) in LTE at about 0.4 eV. These constraints show only the statistical uncertainty in the fitting model.

In the text

Fig. 4. Grotrian diagram showing the lowest energy levels in Y+. The transitions between energy levels are indicated, with the line transparency set by the transition strength, log gf. The 4d5p–4d2 (28 000 cm−1–14 000 cm−1) transitions are shown with thick lines in blue. Radiative transfer modelling analysis in Vieira et al. (2023) tentatively links the 4d5p–4d2 (28 000 cm−1–8000 cm−1) transitions to absorption at λ < 450 nm in the epoch 1 spectrum; other notable transitions are around 300 nm, too short to be observed by X-shooter (4d5p–4d5s, 28 000 cm−1–1000 cm−1), and signatures of the 5s5p–4d2 transitions (24 000 cm−1–14 000 cm−1) at 1 μm, which are much weaker than the Sr+ line transitions at similar wavelengths. The 4d5p–4d2 transitions have lower oscillator strength than several of the Y+ lines at shorter wavelengths, but due to the rapid decline of opacity with wavelength for typical r-process elements, these transitions have higher opacity compared to the continuum at 760 nm. When no lines are drawn between levels, the transition is negligible with log gf < −2.

In the text

Fig. 5. X-shooter spectra with model fit including Y+ lines with λ ≥ 700 nm for epochs 4 and 5. The blackbody continuum model is shown as a dashed black line, the blackbody plus the Sr+ P Cygni lines as a dotted purple line, and the full model including both Sr+ and Y+ as a solid red line. The partially transparent purple and red regions indicate the successive contributions to these models of the Sr+ and Y+ P Cygni lines, respectively. The grey hatched region indicates the overlapping noisy regions between the UVB and VIS arms of the spectrograph. The Y+ lines (from the 4d2–4d5p transitions) have fixed central wavelengths and their relative strengths are set by the LTE condition.

In the text

Fig. 6. Example TARDIS model spectra including yttrium (blue line) and excluding yttrium (red dashed line) but otherwise composed of r-process elements with solar r-process abundance ratios (Lodders et al. 2009; Bisterzo et al. 2014). The left panel is limited to the first r-process peak elements (34 ≤ Z ≤ 50), while the right panel shows a model with all r-process elements (33 ≤ Z ≤ 92). The spectrum at the photosphere is modelled as a blackbody, which is then propagated through a relativistically expanding atmosphere (0.15c ≤ v ≤ 0.2c) using the radiative transfer equations with a radially uniform composition using the Kurucz (2018) line list. The red region indicates the 760 nm feature, which is clearly associated with the presence of yttrium. The yttrium abundance in the model constitutes 0.8% of the total mass. The 0.8–1 μm feature is the Sr+ P Cygni feature (Watson et al. 2019).

In the text

Fig. 7. Example TARDIS model spectra including yttrium (blue) and excluding yttrium (red) but otherwise composed of all r-process elements (33 ≤ Z ≤ 92) with solar r-process abundance ratios (Lodders et al. 2009; Bisterzo et al. 2014) for three different ejecta velocities. The three bottom panels show the difference in the spectral shape between the yttrium-free and yttrium-containing model for each characteristic velocity. Only when the characteristic ejecta velocity has decreased below ∼0.2c (at about 3–4 days post merger) does the P Cygni feature become apparent, because it no longer merges with lower-wavelength yttrium features.

In the text

Fig. A.1. Corner plots showing the posterior probability distribution of key parameters for epoch 4 with a single P Cygni line fit. The fit parameters are (1) the cross-sectional velocity, v⊥, derived from the blackbody normalisation fit assuming Planck cosmology (Planck Collaboration VI 2020), and several parameters of the 760 nm P Cygni feature, including (2) the photospheric velocity, vphot, (3) the maximum outer ejecta velocity, vmax, (4) the velocity range, Δv = vmax − vphot, (5) the optical depth of the line, τ, and (6) the central wavelength, λ0, in μm. Velocities are indicated in units of the speed of light, c. These constraints show only the statistical uncertainty in the fitting model.

In the text

Fig. A.2. Corner plots showing the posterior probability distribution of key parameters for epoch 5 with a single P Cygni line fit. The fit parameters are (1) the cross-sectional velocity, v⊥, derived from the blackbody normalisation fit assuming Planck cosmology (Planck Collaboration VI 2020), and several parameters of the 760 nm P Cygni feature, including (2) the photospheric velocity, vphot, (3) the maximum outer ejecta velocity, vmax, (4) the velocity range, Δv = vmax − vphot, (5) the optical depth of the line, τ, and (6) the central wavelength, λ0, in μm. Velocities are indicated in units of the speed of light, c. These constraints show only the statistical uncertainty in the fitting model.

In the text