11institutetext: Cosmic Dawn Center (DAWN) 22institutetext: Niels Bohr
Institute, University of Copenhagen, Blegdamsvej 17, København 2100, Denmark

#  Discovery of a 760 nm P Cygni line in AT2017gfo: Identification of yttrium
in the kilonova photosphere

Albert Sneppen Corresponding author email<EMAIL_ADDRESS>Discovery
of a 760 nm P Cygni line in AT2017gfo: Identification of yttrium in the
kilonova photosphereDiscovery of a 760 nm P Cygni line in AT2017gfo:
Identification of yttrium in the kilonova photosphereDiscovery of a 760 nm P
Cygni line in AT2017gfo: Identification of yttrium in the kilonova
photosphereDiscovery of a 760 nm P Cygni line in AT2017gfo: Identification of
yttrium in the kilonova photosphere Darach Watson Discovery of a 760 nm P
Cygni line in AT2017gfo: Identification of yttrium in the kilonova
photosphereDiscovery of a 760 nm P Cygni line in AT2017gfo: Identification of
yttrium in the kilonova photosphereDiscovery of a 760 nm P Cygni line in
AT2017gfo: Identification of yttrium in the kilonova photosphereDiscovery of a
760 nm P Cygni line in AT2017gfo: Identification of yttrium in the kilonova
photosphere

(Received 15 March 2023 / Accepted 3 June 2023 )

Neutron star mergers are believed to be a major cosmological source of rapid
neutron-capture elements. The kilonovae associated with neutron star mergers
have to date yielded only a single well-identified spectral signature: the P
Cygni line of Sr+ at about $1\,\mu$m in the spectra of the optical transient,
AT2017gfo. Such P Cygni lines are important, because they provide significant
information not just potentially on the elemental composition of the merger
ejecta, but also on the velocity, geometry, and abundance stratification of
the explosion. In this paper we show evidence for a previously unrecognised P
Cygni line in the spectra of AT2017gfo that emerges several days after the
explosion, located at $\lambda\approx 760\,$nm. We show that the feature is
well-reproduced by 4d2–4d5p transitions of Y+, which have a weighted mean
wavelength around 760–770 nm, with the most prominent line at 788.19 nm. While
the observed line is weaker than the Sr+ feature, the velocity stratification
of the new line provides an independent constraint on the expansion rate of
the ejecta similar to the constraints from Sr+.

##  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 $\lambda\leq 500$ nm.
The identification of features at $\lambda\approx$ 1200–1400 nm with
absorption due to La2+ and Ce2+ 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\,\mu$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\,\mu$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 paper, 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\,\mu$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 4d2–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.

Figure 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. Figure 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. Table 1: Epochs 4 and 5 best-fit parameters for single-line P Cygni (with free-to-vary central wavelength) and for the Yttrium P Cygni model. This includes (1) the best-fit cross-sectional velocity inferred from the blackbody radii assuming Planck cosmology, $v_{\perp}$, and several parameters of the P Cygni feature, including (2) the photospheric velocity, $v_{\rm phot}$, (3) the maximum outer ejecta velocity, $v_{\rm max}$, (4) the e-folding scale of the optical depth, $v_{e}$, (5) the optical depth of the line at the photospheric velocity, $\tau$, and (6) the central wavelength, $\lambda_{0}$. The uncertainties are derived from the 16th and 84th percentile and only include the statistical uncertainties in the fitting framework; i.e. they do not include the potentially significant systematic effects from line-blending and reverberation effects (see Sneppen et al. 2023a). All velocities are given in units of the speed of light, $c$.  | $v_{\rm\perp}$ | $v_{\rm phot}$ | $v_{\rm max}$ | $v_{\rm e}$ | $\tau$ |  $\lambda_{0}$ [$\mu$m]
---|---|---|---|---|---|---
Single P Cygni |  |  |  |  |
Epoch 4 | $0.129_{-0.007}^{+0.007}$ | $0.156_{-0.001}^{+0.001}$ | $0.183_{-0.001}^{+0.001}$ | $0.5_{-0.2}^{+0.3}$ | $0.54_{-0.02}^{+0.02}$ | $0.751_{-0.001}^{+0.001}$
Epoch 5 | $0.112_{-0.006}^{+0.006}$ | $0.191_{-0.004}^{+0.005}$ | $0.210_{-0.004}^{+0.005}$ | $0.7_{-0.2}^{+0.2}$ | $4.7_{-1.3}^{+3.1}$ | $0.769_{-0.001}^{+0.001}$
| $v_{\rm\perp}$ | $v_{\rm phot}$ | $v_{\rm max}$ | $v_{\rm e}$ | $\tau_{\rm avg}$ |
Y+ P Cygni from 4d5p–4d2 line-transitions with $\lambda\in[726\,{\rm nm};788\,{\rm nm}]$ |  |
Epoch 4 | $0.129_{-0.007}^{+0.007}$ | $0.168_{-0.003}^{+0.004}$ | $0.220_{-0.002}^{+0.002}$ | $0.4_{-0.1}^{+0.1}$ | $0.161_{-0.011}^{+0.013}$ |
Epoch 5 | $0.110_{-0.006}^{+0.006}$ | $0.176_{-0.002}^{+0.003}$ | $0.191_{-0.003}^{+0.003}$ | $0.6_{-0.2}^{+0.2}$ | $1.74_{-0.18}^{+0.21}$ |

###  2.1 Statistical significance

To determine the significance of the feature, we compare the $\chi^{2}$ 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 and physically 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\mu$m and $2.0\mu$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 Model111The adopted implementation is Ulrich Noebauer’s
pcygni_profile.py https://github.com/unoebauer/public-astro-tools., where the
profile is set by several properties of the kilonova atmosphere. The profile
is expressed in terms of the rest wavelength, $\lambda_{0}$, and the line
optical depth, $\tau$, with a velocity stratification parameterised by a
scaling velocity, $v_{e}$, a photospheric velocity, $v_{\rm phot}$, and a
maximum ejecta velocity, $v_{\rm 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_{\rm\nu,4th}^{2}=1.65$ and $\chi_{\rm\nu,5th}^{2}=1.88$. In contrast,
fitting the residual structure with the additional P Cygni model instead
yields $\chi_{\rm\nu,4th}^{2}=1.17$ and $\chi_{\rm\nu,5th}^{2}=1.21$ over the
same wavelengths (see the posterior landscape in Appendix Fig. 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.

Nevertheless, the substantial improvement in $\chi^{2}$ 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^{2}_{\nu}=1$. In this case,
the simpler continuum model without the 760 nm P Cygni provides a poorer fit
to the data, with $\chi^{2}$ values going from 14 321 to 10 102 and 15 686 to
10 096 for epochs 4 and 5, respectively. This results in $\Delta\chi^{2}=4219$
and 5590 for five additional degrees of freedom. This corresponds to a null-
hypothesis probability rejected at 26$\sigma$ and 33$\sigma$, 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
$\Delta\mathrm{DIC}<-4000$ for both epochs corresponding to decisive evidence
on the Jeffreys scale (e.g. Jeffreys 1939; Kass & Raftery 1995).

Figure 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 $4{\rm d}5{\rm p}$-$4{\rm d}^{2}$
transitions with $\lambda\geq 700$ nm. The fit parameters are (1) the cross-
sectional velocity, $v_{\perp}$, derived from the blackbody normalisation fit
assuming Planck cosmology (Planck Collaboration et al. 2020), (2) the
photospheric velocity, $v_{\rm phot}$, (3) the maximum outer ejecta velocity,
$v_{\rm max}$, (4) the velocity range, $\Delta v=v_{\rm max}-v_{\rm phot}$,
and (5) the average optical depth of the Y+ lines weighted by their relative
strength, $\tau_{avg}$. Velocities are given in units of the speed of light,
$c$. The central wavelength of a line, $\lambda_{0}$, and its relative
strength is set by the $4$d$5$p to $4$d2 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.

###  2.2 The central wavelength

The best-fit central wavelength of the single P Cygni line fit is $750.5\pm 1$
nm and $769.1\pm 1.2$ nm for epochs 4 and 5, respectively. The unweighted
average of these would indicate a central wavelength of roughly $760\pm 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
$\approx 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_{\parallel}$, from the Sr+ lines and the cross-sectional velocity,
$v_{\perp}$, 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 et al. 2020) and the Hubble flow recession velocity of
the host galaxy in Mukherjee et al. (2021). In addition, the line fits for
$\textrm{Sr}^{+}$ suggest the photospheric velocity is
$v_{\parallel}(\textrm{Sr}^{+})\approx 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:

$\Upsilon(\textrm{Sr}^{+})=\frac{v_{\perp}-v_{\parallel}(\textrm{Sr}^{+})}{v_{\perp}+v_{\parallel}(\textrm{Sr}^{+})}.$
(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_{\parallel}{\rm(760\,nm)}\approx 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_{\perp}<v_{\parallel}$ 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.

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_{\rm phot}$ and $v_{\rm 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\,\mu$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.

Figure 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–14000 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 $\lambda<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 (24000 cm-1–14000 cm-1) at 1$\mu$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$.  Table 2: 4d5p to 4d2 transitions for Y+. All transitions are used for the tardis modelling, while we tested the Y+ P Cygni framework using only the transitions listed here with $\lambda>700$ nm and all the transitions in this table. Vacuum wavelengths ($\lambda$), energy levels ($E$), and oscillator strengths ($\log{gf}$) for the Y+ 4d5p to 4d2 lines are from Biémont et al. (2011). $\lambda$ | Lower level | Upper level | $\log{gf}$
---|---|---|---
(nm) | $4{\rm d}^{2}$ |  $E$ (cm-1) | $4{\rm d}5{\rm p}$ |  $E$ (cm-1) |
788.19 | $b^{1}D_{2}$ | 14 832 | $z^{1}P_{1}^{\circ}$ | 27 517 | $-0.60$
745.03 | $a^{3}P_{2}$ | 14 098 | $z^{1}P_{1}^{\circ}$ | 27 517 | $-1.44$
733.29 | $a^{3}P_{0}$ | 13 883 | $z^{1}P_{1}^{\circ}$ | 27 517 | $-1.98$
726.42 | $b^{1}D_{2}$ | 14 832 | $z^{3}D_{1}^{\circ}$ | 28 595 | $-1.11$
689.60 | $a^{3}P_{2}$ | 14 098 | $z^{3}D_{1}^{\circ}$ | 28 595 | $-1.67$
685.82 | $a^{3}P_{1}$ | 14 018 | $z^{3}D_{1}^{\circ}$ | 28 595 | $-1.69$
683.24 | $a^{3}P_{2}$ | 14 098 | $z^{3}D_{2}^{\circ}$ | 28 730 | $-1.86$
679.54 | $a^{3}P_{0}$ | 13 883 | $z^{3}D_{1}^{\circ}$ | 28 595 | $-1.54$
679.53 | $a^{3}P_{1}$ | 14 018 | $z^{3}D_{2}^{\circ}$ | 28 730 | $-1.03$
661.37 | $a^{3}P_{2}$ | 14 098 | $z^{3}D_{3}^{\circ}$ | 29 213 | $-0.83$

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 $4$d$5$p to $4$d2. 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.

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}\approx 10^{6}-10^{8}\,\rm{cm}^{3}$) 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\approx 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.

Figure 5: X-shooter spectra with model fit including Y+ lines with
$\lambda\geq 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.

###  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 $\lambda\geq 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 $4{\rm d}5{\rm
p}$-$4{\rm d}^{2}$ transitions, but as the $\sim 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 Sec. 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 $\sim 670$ nm and $\sim 760$ nm transitions.

The column density inferred from fitting the P Cygni Y+ 4d5p-4d2 transition
fits suggests total yttrium masses of $0.6\times 10^{-4}\,\rm{M}_{\odot}$ and
$2.4\times 10^{-4}\,\rm{M}_{\odot}$ 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_{\rm ej}\approx 3-5\times 10^{-2}\,\rm{M}_{\odot}$ (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).

Figure 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\leq Z\leq 50$), while the right panel shows a model with all r-process
elements ($33\leq Z\leq 92$). The spectrum at the photosphere is modelled as a
blackbody, which is then propagated through a relativistically expanding
atmosphere ($0.15c\leq v\leq 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$\mu$m feature is the Sr+ P Cygni feature (Watson et
al. 2019).

###  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\leq Z\leq 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\leq Z\leq
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 4d2–4d5p transitions of Y+.

In addition to the $\sim 760$ nm feature that we observe, the presence of
yttrium also results in (1) strong absorption of flux at short wavelengths
$\leq 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$\mu$m feature from the $5$s$5$p-$4$d2 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 $\sim 760$ nm feature compared to its lower-wavelength
counterparts could therefore be considered tentative evidence of a
significantly higher opacity below $\sim 700$ nm, which is consistent with the
arguments in Gillanders et al. (2022) that for synthetic models to match the
observed flux distribution for $\lambda<750$ nm in later epochs requires the
inclusion of lanthanides. However, as we cannot confidently exclude a feature
at $\lambda\sim 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\leq Z\leq 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 $4{\rm d}5{\rm p}$ to $4{\rm d}^{2}$ transitions were
mistakenly omitted from the compiled line list in this analysis (N. Vieira,
personal communication).

Figure 7: Example tardis model spectra including yttrium (blue) and excluding
yttrium (red) but otherwise composed of all r-process elements ($33\leq Z\leq
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 $\sim 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.

##  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 $4$d$5$p to $4$d2 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.2$c$, 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 $4{\rm d}5{\rm
p}\leavevmode\nobreak\ {\rm a}^{3}{\rm P}_{2}$–$4{\rm
d}^{2}\leavevmode\nobreak\ {\rm z}^{3}{\rm D}_{3}^{\circ}$ 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 $\approx 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.

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_{\rm phot}$ or
$\lambda_{0}$ 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_{\perp}$, 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 $\sim 1\,\mu$m line previously associated with Sr+
(Watson et al. 2019). We propose a possible identification of this feature as
the sum of the 4d2-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.

##  6 Acknowledgements

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

  * Biémont et al. (2011) Biémont É., et al., 2011, MNRAS, 414, 3350
  * Bisterzo et al. (2014) Bisterzo S., Travaglio C., Gallino R., Wiescher M., Käppeler F., 2014, ApJ, 787, 10
  * Coulter et al. (2017) Coulter D. A., et al., 2017, Science, 358, 1556
  * Domoto et al. (2021) Domoto N., Tanaka M., Wanajo S., Kawaguchi K., 2021, ApJ, 913, 26
  * Domoto et al. (2022) Domoto N., Tanaka M., Kato D., Kawaguchi K., Hotokezaka K., Wanajo S., 2022, ApJ, 939, 8
  * Drout et al. (2017) Drout M. R., et al., 2017, Science, 358, 1570
  * Ekanger et al. (2023) Ekanger N., Bhattacharya M., Horiuchi S., 2023, arXiv e-prints, p. arXiv:2303.00765
  * Gillanders et al. (2022) Gillanders J. H., Smartt S. J., Sim S. A., Bauswein A., Goriely S., 2022, MNRAS, 515, 631
  * Jeffreys (1939) Jeffreys H., 1939, Theory of Probability
  * Kasen et al. (2017) Kasen D., Metzger B., Barnes J., Quataert E., Ramirez-Ruiz E., 2017, Nature, 551, 80
  * Kass & Raftery (1995) Kass R. E., Raftery A. E., 1995, Journal of the American Statistical Association, 90, 773
  * Kerzendorf & Sim (2014) Kerzendorf W. E., Sim S. A., 2014, MNRAS, 440, 387
  * Kurucz (2018) Kurucz R. L., 2018, in Workshop on Astrophysical Opacities. p. 47
  * Lattimer et al. (1977) Lattimer J. M., Mackie F., Ravenhall D. G., Schramm D. N., 1977, ApJ, 213, 225
  * Liddle (2007) Liddle A. R., 2007, MNRAS, 377, L74
  * Lodders et al. (2009) Lodders K., Palme H., Gail H. P., 2009, Landolt Börnstein, 4B, 712
  * Malesani et al. (2017) Malesani D., Watson D., Hjorth J., 2017, GRB Coordinates Network, 21577, 1
  * Mukherjee et al. (2021) Mukherjee S., Lavaux G., Bouchet F. R., Jasche J., Wandelt B. D., Nissanke S., Leclercq F., Hotokezaka K., 2021, A&A, 646, A65
  * Perego et al. (2021) Perego A., Thielemann F. K., Cescutti G., 2021, in , Handbook of Gravitational Wave Astronomy. p. 13, doi:10.1007/978-981-15-4702-7˙13-1
  * Perego et al. (2022) Perego A., et al., 2022, ApJ, 925, 22
  * Pian et al. (2017) Pian E., et al., 2017, Nature, 551, 67
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Ryabchikova et al. (2015) Ryabchikova T., Piskunov N., Kurucz R. L., Stempels H. C., Heiter U., Pakhomov Y., Barklem P. S., 2015, Phys. Scr, 90, 054005
  * Shappee et al. (2017) Shappee B. J., et al., 2017, Science, 358, 1574
  * Sim (2017) Sim S. A., 2017, in Alsabti A. W., Murdin P., eds, , Handbook of Supernovae. p. 769, doi:10.1007/978-3-319-21846-5˙28
  * Smartt et al. (2017) Smartt S. J., et al., 2017, Nature, 551, 75
  * Sneppen (2023) Sneppen A., 2023, arXiv e-prints, p. arXiv:2306.05452
  * Sneppen et al. (2023a) Sneppen A., Watson D., Poznanski D., Just O., Bauswein A., Wojtak R., 2023a, arXiv e-prints, p. arXiv:2306.12468
  * Sneppen et al. (2023b) Sneppen A., Watson D., Bauswein A., Just O., Kotak R., Nakar E., Poznanski D., Sim S., 2023b, Nature, 614, 436
  * Symbalisty & Schramm (1982) Symbalisty E., Schramm D. N., 1982, Astrophys. Lett., 22, 143
  * Tanaka et al. (2020) Tanaka M., Kato D., Gaigalas G., Kawaguchi K., 2020, MNRAS, 496, 1369
  * Tanvir et al. (2013) Tanvir N. R., Levan A. J., Fruchter A. S., Hjorth J., Hounsell R. A., Wiersema K., Tunnicliffe R. L., 2013, Nature, 500, 547
  * Tarumi et al. (2023) Tarumi Y., Hotokezaka K., Domoto N., Tanaka M., 2023, arXiv e-prints, p. arXiv:2302.13061
  * Thomas et al. (2011) Thomas R. C., Nugent P. E., Meza J. C., 2011, PASP, 123, 237
  * Vieira et al. (2023) Vieira N., Ruan J. J., Haggard D., Ford N., Drout M. R., Fernández R., Badnell N. R., 2023, ApJ, 944, 123
  * Watson et al. (2019) Watson D., et al., 2019, Nature, 574, 497
  * Waxman et al. (2018) Waxman E., Ofek E. O., Kushnir D., Gal-Yam A., 2018, MNRAS, 481, 3423

## 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.

Figure 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_{\perp}$, derived from the blackbody
normalisation fit assuming Planck cosmology (Planck Collaboration et al.
2020), and several parameters of the 760 nm P Cygni feature, including (2) the
photospheric velocity, $v_{\rm phot}$, (3) the maximum outer ejecta velocity,
$v_{\rm max}$, (4) the velocity range, $\Delta v=v_{\rm max}-v_{\rm phot}$,
(5) the optical depth of the line, $\tau$, and (6) the central wavelength,
$\lambda_{0}$, in $\mu$m. Velocities are indicated in units of the speed of
light, $c$. These constraints show only the statistical uncertainty in the
fitting model.  Figure 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_{\perp}$, derived from
the blackbody normalisation fit assuming Planck cosmology (Planck
Collaboration et al. 2020), and several parameters of the 760 nm P Cygni
feature, including (2) the photospheric velocity, $v_{\rm phot}$, (3) the
maximum outer ejecta velocity, $v_{\rm max}$, (4) the velocity range, $\Delta
v=v_{\rm max}-v_{\rm phot}$, (5) the optical depth of the line, $\tau$, and
(6) the central wavelength, $\lambda_{0}$, in $\mu$m. Velocities are indicated
in units of the speed of light, $c$. These constraints show only the
statistical uncertainty in the fitting model.