11institutetext: Cosmic Dawn Center (DAWN) 22institutetext: Niels Bohr
Institute, University of Copenhagen, Jagtvej 128, DK-2200, Copenhagen N,
Denmark 33institutetext: Trottier Space Institute at McGill and Department of
Physics, McGill University, 3600 rue University, Montreal, Québec, H3A 2T8,
Canada 44institutetext: South African Astronomical Observatory, P.O. Box 9,
Observatory 7935, Cape Town, South Africa 55institutetext: Finnish Centre for
Astronomy with ESO, FINCA, University of Turku, Turku, FI-20014, Finland

# Emergence hour-by-hour of r-process features in the kilonova AT2017gfo

Albert Sneppen Emergence hour-by-hour of r-process features in the kilonova
AT2017gfoEmergence hour-by-hour of r-process features in the kilonova
AT2017gfoEmergence hour-by-hour of r-process features in the kilonova
AT2017gfoEmergence hour-by-hour of r-process features in the kilonova
AT2017gfo 0000-0002-5460-6126 Darach Watson Emergence hour-by-hour of
r-process features in the kilonova AT2017gfoEmergence hour-by-hour of
r-process features in the kilonova AT2017gfoEmergence hour-by-hour of
r-process features in the kilonova AT2017gfoEmergence hour-by-hour of
r-process features in the kilonova AT2017gfo 0000-0002-4465-8264 Rasmus
Damgaard Emergence hour-by-hour of r-process features in the kilonova
AT2017gfoEmergence hour-by-hour of r-process features in the kilonova
AT2017gfoEmergence hour-by-hour of r-process features in the kilonova
AT2017gfoEmergence hour-by-hour of r-process features in the kilonova
AT2017gfo 0009-0002-5765-4601 Kasper E. Heintz Emergence hour-by-hour of
r-process features in the kilonova AT2017gfoEmergence hour-by-hour of
r-process features in the kilonova AT2017gfoEmergence hour-by-hour of
r-process features in the kilonova AT2017gfoEmergence hour-by-hour of
r-process features in the kilonova AT2017gfo 0000-0002-9389-7413 Nicholas
Vieira Emergence hour-by-hour of r-process features in the kilonova
AT2017gfoEmergence hour-by-hour of r-process features in the kilonova
AT2017gfo 0000-0001-7815-7604 Petri Väisänen and Antoine Mahoro Emergence
hour-by-hour of r-process features in the kilonova AT2017gfoEmergence hour-by-
hour of r-process features in the kilonova AT2017gfoEmergence hour-by-hour of
r-process features in the kilonova AT2017gfoEmergence hour-by-hour of
r-process features in the kilonova AT2017gfo 0000-0001-7673-4850 Emergence
hour-by-hour of r-process features in the kilonova AT2017gfoEmergence hour-by-
hour of r-process features in the kilonova AT2017gfo 0000-0002-6518-781X

(Received date / Accepted date )

The spectral features in the optical/near-infrared counterparts of neutron
star mergers (kilonovae, KNe), evolve dramatically on hour timescales. To
examine the spectral evolution we compile a temporal series complete at all
observed epochs from 0.5 to 9.4 days of the best optical/near-infrared (NIR)
spectra of the gravitational-wave detected kilonova AT2017gfo. Using our
analysis of this spectral series, we show that the emergence times of spectral
features place strong constraints on line identifications and ejecta
properties, while their subsequent evolution probes the structure of the
ejecta. We find that the most prominent spectral feature, the 1 $\upmu$m P
Cygni line, appears suddenly, with the earliest detection at 1.17 days. We
find evidence in this earliest feature for the fastest kilonova ejecta
component yet discovered, at 0.40–0.45$c$; while across the observed epochs
and wavelengths, the velocities of the line-forming regions span nearly an
order of magnitude, down to as low as 0.04–0.07$c$. The time of emergence
closely follows the predictions for Sr ii, due to the rapid recombination of
Sr iii under local thermal equilibrium (LTE) conditions. The time of
transition between the doubly and singly ionised states provides the first
direct measurement of the ionisation temperature, This temperature is highly
consistent, at the level of a few percent, with the temperature of the emitted
blackbody radiation field. Further, we find the KN to be isotropic in
temperature, i.e. the polar and equatorial ejecta differ by less than a few
hundred Kelvin or $\lesssim 5$%, in the first few days post-merger, based on
measurements of the reverberation time-delay effect. This suggests that a
model with very simple assumptions, with single-temperature LTE conditions,
reproduces the early kilonova properties surprisingly well.

##  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 time of appearance 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 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, 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 UV-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.

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 as
accurate a description as possible of when line features first appear. We
analyse the emergence time of the 760 nm P Cygni feature and its
identification with Y ii in Sect. 4.1. Similarly, we analyse the emergence of
the strongest feature (at $\sim 1$ $\upmu$m) and check its consistency with
the Sr ii interpretation in Sect. 4.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 particularly
the 1.4 $\upmu$m feature displays an apparent P Cygni-like spectral shape. We
discuss the possible different origin of the NIR features compared to the
optical lines (Sect. 4.3-4.5) and how this relates to the complementary nature
of their velocity structure.

Figure 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),
Australian National University 2.3-m telescope (ANU, Andreoni et al. (2017)),
Southern African Large Telescope (SALT, Buckley et al. (2018)). From 1.4 days
post-merger, the X-shooter spectrograph at the European Southern Observatory’s
Very Large Telescope details 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 days 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. In the zoom-in
sub-panels, individual spectral lines are highlighted at different epochs,
including (1) the 760 nm (2) 1.0 $\upmu$m, and (3) 1.4 $\upmu$m P Cygni
features, and (4) the NIR line-emission features.

##  2 Overview of the data

On 2017 August 17 at 12:41:04 UTC, the merging neutron star gravitational wave
event GW170817 was detected and 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 restframe of the host galaxy NGC 4993, $z_{\rm
cosmic}\approx 0.01$ (Hjorth et al. 2017; Howlett & Davis 2020; Nicolaou et
al. 2020).

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 the Las
Campanas Observatory, Chile (Shappee et al. 2017). These spectra display a
featureless blue continuum 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 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, illustrating 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 the
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\simeq$ 6000–7000 K. We note the flux-
calibration for the SED is somewhat uncertain, so the exact blackbody
temperature inferred from the spectral slope is 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\pm 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. 4.2 is the first appearance of the 1 $\upmu$m P Cygni feature. Due to
the importance of the precise emergence time of the $\sim 1$ $\upmu$m feature,
we have 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 Sect. Appendix 1.1.

1.4–9.4 days post-merger: From 1.4 days post-merger, AT2017gfo was monitored
with 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, with the
blackbody approximation increasingly poor, and becoming more dominated by the
NIR emission, where a series of spectral components emerge more prominently.
Many spectra were taken near-contemporaneously with the X-shooter datasets,
e.g. 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 overall provides the best constraints. Comparing across the
approx. one hour difference in observing time between these observations can
help in constraining the reverberation wave of recombining species as will
detail in Sect. 4.2.1. We therefore use 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 compliment 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 $\upmu$m. This region is difficult to confidently
constrain from ground-based observations and therefore we include the 4.9 day
post-merger _HST_ spectrum in our spectral set.

For consistency in comparing spectra, we apply 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 et al. (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 is assumed to be negligible. For detailed arguments on the
dust extinction see Sneppen et al. (2023c).

Figure 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 observed time post-merger. The
ionisation states are calculated using the Saha equation, where we assume
electron density, $n_{e}(t)=10^{8}$ cm${}^{-3}(t/1.5\,\textrm{days})^{-3}$ and
the blackbody temperature evolution found for AT2017gfo (Sneppen et al.
2023b). The blue and red vertical regions indicate respectively 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 dominating the opacity make rapid transitions
between ionisation states. For $t\gtrapprox 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).

##  3 The one day transition and the 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 37Rb and 55Cs for comparison, otherwise
these are the elements which produce the strongest lines due to their few
valence electrons, low-lying energy-states and small partition functions (i.e.
Sr, Y, Zr, Ba) and the elements with high optical/NIR opacity with many low-
lying states (i.e. the lanthanides). Both of these latter sets of elements
display a rapid change in ionisation around one day post-merger.

For Fig. 2, we assume 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^{8}$
cm${}^{-3}(t/1.5\ {\rm days})^{-3}$ (Sneppen et al. 2023a). At early times the
temperature is cooling dramatically with $T\propto t^{-\alpha}$ with
$\alpha\in[0.5;0.8]$ (Drout et al. 2017; Waxman et al. 2018; Sneppen et al.
2023a), which is why the abundance of lowly ionised species rapidly transition
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/time
follows from the Saha equation, where the relative density of ions $n_{\rm
i+1}/n_{\rm i}$ (with the subscript indicating the state of ionisation) is
given by the electron-density, $n_{e}$, temperature, $T$, the ionisation
energy, $\mathcal{X}_{\rm i}$, the partition functions for each ion, $Z_{\rm
i+1}$ & $Z_{\rm i}$, and the thermal electron de-Broglie wavelength,
$\lambda_{\rm dB}$:

$\frac{n_{\rm i+1}}{n_{\rm i}}=\frac{2}{\lambda_{dB}^{3}}\frac{Z_{\rm
i+1}}{Z_{\rm i}}\frac{1}{n_{e}}e^{-\frac{\mathcal{X}_{\rm i}}{k_{B}T}}$ (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, the ionisation transition is expected to be
abrupt under LTE conditions. For Sr ii at 1.4 days (where the
$n_{\rm{Sr\textsc{iii}}}/n_{\rm{Sr\textsc{ii}}}$ fraction is near-unity) the
ionisation energy is $\mathcal{X}_{\rm 1}\approx 30k_{B}T$, which implies a
subtle decrease in temperature results in a large change in the Sr ii number
density. As the ion recombines, the column density, and hence the optical
depth, $\tau$, of the lines produced by the relevant species will change by
orders of magnitude. This implies a drastic change in the transmission
($I=I_{0}e^{-\tau}$) 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 – transitioning from doubly to
singly-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
are (by their nature) more smoothly and continuously evolving both during
their (comparatively) slower formation-time and their subsequent evolution.

##  4 Observed spectral features

In the following, we consider each of the major 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, potentially, 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 use our spectral compilation to extend
the velocity constraints to a broader temporal span including the NIR
features.

Figure 3: Emergence of the 1 $\mu$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 contain no strong deviation 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. We
note that 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 is inconsistent with a He i $\lambda 1083.0$ nm interpretation (see
Sneppen et al. in preparation). That the temperature-evolution constrained
from the blackbody spectra predicts the formation time, indicates that the
ejecta is close to LTE conditions and that the radiation blackbody temperature
is a good probe of the ejecta ionisation temperature.

###  4.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), is well reproduced
by the 4d2–4d5p transitions of Y ii (Sneppen & Watson 2023). At earlier times
the larger velocities mean that the various lines of Y ii blend together
thereby concealing the characteristic P Cygni spectral shape. Modelling with
the TARDIS code in Sneppen & Watson (2023) shows that a prediction of the
yttrium identification is that the P Cygni structure emerges clearly when the
velocity drops to about 0.15–0.2$c$, which is 3–4 days after the merger. Tak
et al. (2024) have now also identified the 4d2–4d5p Y ii transitions as likely
to strongly affect 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 $\lambda\leq 400$ nm at early times (Gillanders et al.
2022; Vieira et al. 2023b). Conversely, Pognan et al. (2023) find that the Y
ii transitions within their models are too weak to account for the observed
760 nm features, and instead propose that in NLTE conditions Rb i may be able
to produce a feature at roughly the right wavelength. However, as they use an
uncalibrated theoretical line-list for yttrium which has incorrect
wavelengths, placing these transitions at 1.1 $\upmu$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 is
presented and their NLTE models cannot extend to the early times (3.4 & 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) $\upmu$m). This means 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 is
in agreement with the velocity inferred from the 1.0 $\upmu$m P Cygni feature
at this time. At later times, the feature evolves analogously to the 1.0
$\upmu$m feature, reverberating into pure emission and fading away towards the
more nebular epochs. Reverberation here means in the sense 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 will discuss
this later-time evolution further in the context of the more prominent and, by
extension, better constrained 1.0 $\upmu$m feature.

###  4.2 P Cygni feature at 1.0 $\mu m$ – Sr ii

The Sr ii 4p64d—4p65p triplet of strong lines producing the 1.0 $\upmu$m P
Cygni feature were 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. 2023b) 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 4p64d—4p65p identification is
that a large UV-absorption is observed at $\lambda\leq 400$ nm around 1.4
days, which could be explained by Sr ii transitions from the ground-state to
4p65p (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. 2023b).
Further corroboration of the Sr ii interpretation is the existence of a 760 nm
P Cygni feature due to yttrium, discussed in the previous section (Sneppen &
Watson 2023).

An alternative interpretation has been proposed, that the 1.0 $\upmu$m P Cygni
feature could originate from the neutral He i $\lambda$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 that 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. (in
preparation), a He i interpretation is not compatible with the emergence-time
and subsequent spectral evolution of the observed 1 $\upmu$m feature. This is
contrary to the situation for Sr ii, as we discuss below.

####  4.2.1 Emergence of the 1.0 $\upmu$m P Cygni feature

A key way to robustly discriminate between 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 is falling rapidly and thus the
dominant ionisation level will also quickly change (see Fig. 2). As mentioned
above, this is in contrast to a helium interpretation, for example, which
conflicts with the first appearance and the subsequent evolution of the
feature. In the following, we show that the Sr ii feature’s early evolution
(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 observed feature’s earliest appearance.

Given the blackbody temperature recovered in the X-shooter epoch 1 spectra
($\sim 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 forming 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 timeframe is observationally probed by the early spectra (see Fig. 3).
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 ($\lambda\gtrapprox
650\,$nm for a blackbody with a temperature $T_{\rm obs}\approx 6400$ K
inferred from the location of the spectral peak). This is the first appearance
of the 1.0 $\upmu$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 statistically
significant: a blackbody continuum model, which adequately describes the
spectrum below 640 nm, describes the spectral region at 640–680 nm poorly with
a $\chi^{2}$ null hypothesis probability of $\sim 10^{-26}$ when this region
is added to the fit, and can thus be rejected at 10$\sigma$. Interestingly,
that the absorption reaches blueward to around $\lambda\approx 640\,$nm, which
indicates that the line-forming region must extend out to quite high
velocities: 0.40–0.45$c$ (given absorption from ejecta expanding along the
line-of-sight). This is the fastest kilonova ejecta component observed so far
– $0.05c$ faster than the fastest component ejecta found from the 1.0 $\upmu$m
P Cygni feature at 1.43 days post-merger (Sneppen et al. 2023b; Vieira et al.
2023a). 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).

Due to the large 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, does
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 beginning to form (Sneppen
et al. 2023a). And 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. 3).

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
ejecta’s ionisation/excitation temperature (we will revisit this further in
Sect. 5.1).

####  4.2.2 Evolution of the 1.0 $\upmu$m P Cygni feature

The 1.0 $\upmu$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
$\tau\in[1;3]$ from 1.17-3.4 days, where-after it becomes optically thin and
ultimately has faded away by 8 days post-merger.

Another important change of the feature is the increasingly prominent nature
of the emission. This is caused by the line’s time delay effect noted above
and the rapidly cooling photosphere, but also the rapidly fading continuum;
these effects creates 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 $\upmu$m feature fades away first in absorption and
subsequently in emission.

Figure 4:  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 this feature could be due to La iii lines (here indicated
with blue dashed lines). The characteristic velocity of this feature is
relatively low, $v_{\rm ph}=0.07\pm 0.01c$, but may be systematically biased
by the strong nearby features (including the prominent reverberation emission
peak around 1.1 $\upmu$m).

###  4.3 P Cygni feature at 1.4 $\mu m$ – La iii?

A feature around 1.4 $\upmu$m has been suggested as being due to the
5p65d—5p64f transitions in La iii (Domoto et al. 2022). The feature is poorly-
constrained from the ground-based spectra as 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
Fig. 4). It is likely that this feature is also present in other epochs,
potentially complicating the modelling of the 1.6 $\upmu$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). Notably, this is so late that
the low blackbody radiation temperature would suggest 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 which allows for doubly ionised lanthanides. We discuss this
prospect further in Sect. 6).

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 $\upmu$m P Cygni feature (see
Fig. 4). Fitting the P Cygni profile with this range of continuum models would
suggest the photospheric velocity is $v_{\rm ph}=0.07\pm 0.01c$. This is a
relatively low velocity – 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 $\upmu$m (Sr ii, Watson et al. 2019). This may suggest
the optically thick region at these NIR wavelengths resides deeper within the
ejecta (where the expansion velocity is smaller). However, the velocity
inferred may also be biased by nearby features, e.g. the feature around 1.6
$\upmu$m and the prominent reverberation emission peak around 1.1 $\upmu$m,
which could create a bias towards lower velocities. In addition, there appears
to be a marginal feature at around 1.2 $\upmu$m which is so far unaccounted
for.

###  4.4 Feature at 1.6 $\mu m$ – Ce iii?

The 1.6 $\upmu$m feature was interpreted by Domoto et al. (2022) as
potentially being due to a combination of lines from Ce iii, (predominantly
4f2–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. (2023a) show
that Ce ii can produce a spectral feature at a similar location in the NIR,
though likely of lesser prominence. Gillanders et al. (2023) additionally
notes the presence of nearby intrinsically weak [I ii] lines which could
contribute significantly to a feature towards late times. There is no
statistically significant feature 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 – seemingly in the narrow time-
frame associated with the transition from doubly to singly ionised r-process
line-forming species. Over the subsequent days the emission feature slowly
becomes more prominent relative to the continuum. As in the case of La iii,
the blackbody temperature and LTE conditions would suggest no significant Ce
iii population would be present in the photosphere from 2.4 days onwards (see
Fig. 2). Again, this would require significant NLTE effects at this time, or
it may be possible that the NIR line-forming region 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 $\upmu$m (2.4 days post-merger) and shifting to 1.58 $\upmu$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 $\sigma$-width of only around 0.04–0.05$c$. This is
very narrow compared to the velocity measurements of other features, including
even the redder 2.1 $\upmu$m emission feature, which at face value appears to
be broader, and which we tackle next. We note that it may be possible this is
a full (and broader) P Cygni absorption/emission feature modified by possible
surrounding features. However there is no very clear signature of such a
hypothesis, as the corresponding absorption would reside within the telluric
regions and be highly degenerate with the 1.4 $\upmu$m P Cygni feature
emission peak.

###  4.5 Feature at 2.1 $\mu m$ – [Te iii]?

The 2.1 $\upmu$m feature is well-modelled as pure emission and was interpreted
as a [Te iii] fine structure line by Hotokezaka et al. (2023) when modelling
the spectra from 7.5–10.5 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 1.6 $\upmu$m
feature width. 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 $\upmu$m at 1.4 days but shifting to 2.1 $\upmu$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 argued to be a kilonova
from a merger (Levan et al. 2024; Yang et al. 2024). That feature was also
tentatively associated with [Te iii] (Levan et al. 2024; Gillanders et al.
2023), and may possibly have substructure, indicative of more than a single
line.

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.4 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 $\upmu$m
feature. This velocity is particularly striking in two respects. First, in
early epochs the velocity-width of the 2.1 $\upmu$m feature is significantly
smaller than the velocity inferred from the continuum, the 760 nm P Cygni
feature, or the 1 $\upmu$m P Cygni feature, while it is consistent with the
1.4 $\upmu$m feature (the reddest of the P Cygni features) from the _HST_
spectra. This suggests 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 small due to
projection-effect for a near-polar observer) or deeper within the ejecta
(visible due to the lower opacity at longer wavelengths). Second, the
feature’s velocity-width does not evolve significantly. If the elements
producing the feature were confined to the merger’s equatorial plane (which
could explain the smaller 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 requires further study.

###  4.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
with the most noticeable peaks 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 they may
not be artefacts of the processing pipeline, but potentially the contribution
of less prominent transitions yet to be identified. Such spectral “wiggles”
with their smaller inferred (projected) velocity-scales could potentially
follow from ejecta expanding 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.
Indeed, these minor wiggles are also marginally detected in later epochs – the
780 nm wiggle has a clear spectral counterpart in the Y ii 4d2–4d5p
transitions, which produces the 760 nm P Cygni feature at later times when the
characteristic velocities decrease and the lines have deblended. The minor
feature around 860 nm is also noticeable in intermediate epochs (3–7 days
post-merger), but 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 $\upmu$m P Cygni feature. Coincidentally, 860 nm is
quite near the rest-wavelength of the strong NIR Ca ii doublet. 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 well with the lines identified with Sr and Y, and if the feature
really were to be due to Ca ii, the modelling is complex and outside the scope
of this paper.

At 1.2 $\upmu$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, as it resides between the prominent 1 $\upmu$m reverberated
emission peak and the absorption of the 1.4 $\upmu$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 $\upmu$m, though they did not
associate those lines with this feature.

##  5 Implications for the ionisation-structure

Having considered the various features and their observed structure, we now
explore the information revealed when combining the various spectral
identifications and the complete dataset. First, in Sect. 5.1, we discuss how
the 1 $\upmu$m feature’s emergence constrains the ejecta to be close to LTE
conditions. Second, in Sect. 5.2, we discuss the complementary nature of
comparing the optical and NIR velocity-structures. Lastly, in Sect. 5.3, we
discuss how the spatial distribution of elements may be constrained from the
various P Cygni features observed.

Figure 5: 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, with lines indicating the evolution modelled as
constraints from reverberation effects within a single epoch; the temperature
errorbars 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\leq v\leq 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 $n_{e}$ a rapid step-change in
ionisation-state occurs over a narrow interval in temperature around 4500 K.
The 1 $\upmu$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).  Figure 6: The 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 $\upmu$m feature in absorption and emission by the
various spectra from ANU, SALT, VLT/X-shooter and Gemini FLAMINGOS-2. Dotted
lines indicate that no feature has clearly emerged (prior to $\sim$1.4 days in
both emission and absorption). We indicates the required fraction of Sr ii
needed to produce the optical depths observed in each spectrum (assuming
$M_{\rm Sr}=6\cdot 10^{-6}$ M☉ homogeneously distributed across $0.2c<v<0.4c$,
see main text). The ANU (X-shooter) and SALT (Gemini) spectrum 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 requires the ionisation temperature 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 requires a near-isotropic temperature.  Figure 7: 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 which 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, such species must reside in relatively hot ejecta.
Hotter ejecta produces a higher NIR continuum flux, which would conversely
require a smaller emitting area/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.

###  5.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 observationally
constrain the allowed level of deviations from LTE conditions at early times.
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. 5. In Sect. 4.2.1 we saw how,
on an hour-by-hour basis, the spectral series allowed us to trace the
formation of the 1 $\upmu$m feature, 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 timeframes 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. If we
assume we don’t know the electron density 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 the difference between the
excitation/ionisation and blackbody radiation temperatures must be within
$\sim$200 K at these times, i.e. consistent at the level of some percent.

These strong constraints on the temperature are attainable because the timing
of ionisation transitions is so strongly dependent on temperature and only
weakly affected by the exact electron density, as illustrated in Fig. 6.
Specifically, this shows the recombination-transition to Sr ii – and in
extension the formation-time of a Sr ii feature – given Eq. 1, typical KN
electron-densities and for various relationships between electron temperature
and the blackbody radiation temperature. For comparison, we calculate the
fraction of strontium required as Sr ii to produce the observed Sobolov
optical depth $\tau=1$, where the strontium-density is computed using the
total strontium mass-estimate from later epochs ($M_{\rm Sr}=6\cdot 10^{-6}$
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 have 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 note there are large uncertainties in the
exact strontium densities or indeed its radial structure. However, especially
at early times the evolution of $\tau(t)$ is very-constrained under LTE.
Particularly, we must i) form a feature on a time-scale 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 respectively no clear absorption
($\tau<0.08$; 3$\sigma$ upper limit) and a strong absorption ($\tau=2.2\pm
0.6$). In physical time, this represents an order-of-magnitude increase 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 $\tau<0.09$ (X-shooter; 3$\sigma$ upper limit) to $\tau=0.8\pm 0.2$
(Gemini FLAMINGOS-2) from 1.43 days to 1.47 days. This highlights a similar
rapid line-formation in the equatorial ejecta, but now confined to the time-
scale of an hour due to the higher cadence around this period.

Curiously, the X-shooter’s NIR-arm actually constrains the scattering of light
(ie. emission) at an earlier physical time than the SALT’s optical scattered
light (ie. the absorption), despite the former spectra being obtained 6 hours
after the latter.

####  5.1.1 The temperature isotropy of the line-forming region

We note that the fact that one can successfully predict both the formation
time in absorption and (in the delayed) emission, accounting only for light
travel time effects, indicates that there are not large temperature gradients
from the polar to equatorial ejecta. Specifically it begins forming in
absorption (i.e. the polar ejecta) between ANU and SALT spectra (indicating
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 the polar and
equatorial line-forming region must reach the temperature for which Sr ii can
first exist (ie. T$\sim$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, i.e. an angular consistency in temperature at the
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 state-of-the-
art 3D models presented in Shingles et al. (2023) large temperature
anisotropies exist between the equatorial and polar ejecta – at the level of
up to 50%. Such class of models cannot produce the emergence of the 1$\upmu$m
features observed in AT2017gfo. Thus, the discovery of these new observables
measuring the ionisation temperature and the temperature anisotropy around
ionisation transitions can thus help to inform and constrain future modelling
of BNS mergers, specifically the thermalization of decay products and the
radiative transfer in KNe.

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

###  5.2  Velocity/ionisation complementarity

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 due to similar ionisation-
energies (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. This 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 species would map onto
different expansion velocities, such that the velocity difference between the
optical P Cygni line profiles (from singly-ionised species) and the near-
infrared features (doubly-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 – however in the Rayleigh-Jeans tail toward the
NIR the temperature and emitting area can be varied along a line-of-degeneracy
to get a similar continuum (Fig. 7).To produce a flux at one specific
wavelength there is a complete degeneracy between emitting area and
temperature, but the NIR spectral slope is still sensitive to the temperature
- with notably the NIR temperature being consistent with the UV/optical-
temperature at the percent level at 1.4 days (Sneppen 2023). For doubly-
ionised species to be present at typical electron densities (at 4 days around
$10^{7}$ cm-3) one would require the electron temperature to be around 4000 K
or hotter (much hotter than the optical-UV temperature at this epoch).
However, such an ejecta temperature is broadly consistent with the NIR
continuum flux for the observed velocity inferred from the 1.4 $\upmu$m line.
But, it is worth noting that these NIR-continuum temperature constraints for
any fixed velocity are likely upper bounds, as they neglect modelling the
contribution of fluorescence from all layers at greater radii, which requires
detailed radiative transfer modelling, so the actual radiation temperature
could be lower.

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

###  5.3 The 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 hints of the spatial geometric distribution of the elements. As we
are observing this object with a near-polar inclination angle
$\theta_{inc}=22\degree\pm 3\degree$ (Mooley et al. 2018, 2022), the
blueshifted absorption and the rest-wavelength emission here closely correlate
respectively with the more polar and equatorial components of the ejecta.

With the emergence of the 1 $\upmu$m feature as detailed in this compiled
spectral series, we have further tightened the connection between Sr ii and
the observed emission of the 1 $\upmu$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. Further, the emission-part of the 0.76 micron
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. Such equatorially-distributed light r-process elements 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 $\upmu$m feature is P Cygni-like as argued in Fig. 4 and 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 which allows both absorption and emission. This would
seemingly suggest a moderate abundance of lanthanides at and above the NIR
photosphere in the polar ejecta – in addition to the more generally expected
equatorial (emission) component found often in simulations (e.g. Just et al.
2023).

Thus, it is likely that the spherical-like 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 are 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
Cygnii structures due to the large spatial asymmetries of the various line-
forming species from hydrodynamic simulations.

##  6 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:

  1. 1.

The earliest chronological appearance of the 1 $\upmu$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. 2.

The 1 $\upmu$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) on the exact emergence-time and early spectral shape
from strontium undergoing a 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 observationally constraining the ejecta to be close
to LTE.

  3. 3.

The observed spectral features at 0.76, 1.0, 1.4, 1.6 and 2.1 $\upmu$m have
been 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), i.e. being respectively
first vs. second r-process-peak elements and single vs. doubly ionised and
hints at a possible spatial separation.

  4. 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. Particularly, the Sr ii and Y ii emissions
requires a sizeable component of first peak elements to be located in the
equatorial plane. Such an 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 $\upmu$m P Cygni feature. Furthermore,
that the recombination of Sr iii to Sr ii in absorption and emission occurs at
the same physical time (ie. 1.3-1.5 days post-merger) indicates a temperature
in the line-forming-region around this time, which is isotropic at the level
of 5% from the polar to the equatorial ejecta.

## Data and code availability

As discussed, the spectral series presented in this paper is composed from 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 meter 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 meter 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, seeSect. Appendix 1.1)
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 use 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.

## Acknowledgements

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
for 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).

## References

  * Abbott et al. (2017) Abbott B. P., et al., 2017, ApJ, 848, L13
  * Andreoni et al. (2017) Andreoni I., et al., 2017, PASA, 34, e069
  * Bauswein et al. (2013) Bauswein A., Goriely S., Janka H. T., 2013, ApJ, 773, 78
  * Buckley et al. (2018) Buckley D. A. H., et al., 2018, MNRAS, 474, L71
  * Burgh et al. (2003) Burgh E. B., Nordsieck K. H., Kobulnicky H. A., Williams T. B., O’Donoghue D., Smith M. P., Percival J. W., 2003, in Iye M., Moorwood A. F. M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes. pp 1463–1471, doi:10.1117/12.460312
  * Chornock et al. (2017) Chornock R., et al., 2017, ApJ, 848, L19
  * Collins et al. (2023) Collins C. E., Bauswein A., Sim S. A., Vijayan V., Martínez-Pinedo G., Just O., Shingles L. J., Kromer M., 2023, MNRAS, 521, 1858
  * Collins et al. (2024) Collins C. E., et al., 2024, MNRAS,
  * 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
  * Evans et al. (2017) Evans P. A., et al., 2017, Science, 358, 1565
  * Fitzpatrick (2004) Fitzpatrick E. L., 2004, in Witt A. N., Clayton G. C., Draine B. T., eds, Astronomical Society of the Pacific Conference Series Vol. 309, Astrophysics of Dust. p. 33 (arXiv:astro-ph/0401344), doi:10.48550/arXiv.astro-ph/0401344
  * Gillanders et al. (2022) Gillanders J. H., Smartt S. J., Sim S. A., Bauswein A., Goriely S., 2022, MNRAS, 515, 631
  * Gillanders et al. (2023) Gillanders J. H., Sim S. A., Smartt S. J., Goriely S., Bauswein A., 2023, MNRAS,
  * Green et al. (2018) Green G. M., et al., 2018, MNRAS, 478, 651
  * Hjorth et al. (2017) Hjorth J., et al., 2017, ApJ, 848, L31
  * Hotokezaka et al. (2023) Hotokezaka K., Tanaka M., Kato D., Gaigalas G., 2023, MNRAS, 526, L155
  * Howlett & Davis (2020) Howlett C., Davis T. M., 2020, MNRAS, 492, 3803
  * Just et al. (2023) Just O., et al., 2023, ApJ, 951, L12
  * Kasen et al. (2015) Kasen D., Fernández R., Metzger B. D., 2015, MNRAS, 450, 1777
  * Levan et al. (2024) Levan A. J., et al., 2024, Nature, 626, 737
  * McCully et al. (2017) McCully C., et al., 2017, ApJ, 848, L32
  * Mooley et al. (2018) Mooley K. P., et al., 2018, Nature, 561, 355
  * Mooley et al. (2022) Mooley K. P., Anderson J., Lu W., 2022, Nature, 610, 273
  * Nicholl et al. (2017) Nicholl M., et al., 2017, ApJ, 848, L18
  * Nicolaou et al. (2020) Nicolaou C., Lahav O., Lemos P., Hartley W., Braden J., 2020, MNRAS, 495, 90
  * 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. (2016) Planck Collaboration et al., 2016, A&A, 596, A109
  * Pognan et al. (2022) Pognan Q., Jerkstrand A., Grumer J., 2022, MNRAS, 513, 5174
  * Pognan et al. (2023) Pognan Q., Grumer J., Jerkstrand A., Wanajo S., 2023, MNRAS, 526, 5220
  * Shappee et al. (2017) Shappee B. J., et al., 2017, Science, 358, 1574
  * Shingles et al. (2023) Shingles L. J., et al., 2023, ApJ, 954, L41
  * Smartt et al. (2017) Smartt S. J., et al., 2017, Nature, 551, 75
  * Sneppen (2023) Sneppen A., 2023, ApJ, 955, 44
  * Sneppen & Watson (2023) Sneppen A., Watson D., 2023, A&A, 675, A194
  * Sneppen et al. (2023a) Sneppen A., Watson D., Gillanders J., 2023a, in prep, 614, 436
  * Sneppen et al. (2023b) Sneppen A., Watson D., Bauswein A., Just O., Kotak R., Nakar E., Poznanski D., Sim S., 2023b, Nature, 614, 436
  * Sneppen et al. (2023c) Sneppen A., Watson D., Poznanski D., Just O., Bauswein A., Wojtak R., 2023c, A&A, 678, A14
  * Tak et al. (2024) Tak D., Uhm Z. L., Gillanders J. H., 2024, arXiv e-prints, p. arXiv:2402.05471
  * Tanaka et al. (2023) Tanaka M., et al., 2023, ApJ, 953, 17
  * Tanvir et al. (2017) Tanvir N. R., et al., 2017, ApJ, 848, L27
  * Tarumi et al. (2023) Tarumi Y., Hotokezaka K., Domoto N., Tanaka M., 2023, arXiv e-prints, p. arXiv:2302.13061
  * Valenti et al. (2017) Valenti S., et al., 2017, ApJ, 848, L24
  * Vieira et al. (2023a) Vieira N., Ruan J. J., Haggard D., Ford N. M., Drout M. R., Fernández R., 2023a, arXiv e-prints, p. arXiv:2308.16796
  * Vieira et al. (2023b) Vieira N., Ruan J. J., Haggard D., Ford N., Drout M. R., Fernández R., Badnell N. R., 2023b, ApJ, 944, 123
  * Villar et al. (2017) Villar V. A., et al., 2017, ApJ, 851, L21
  * 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
  * Yang et al. (2024) Yang Y.-H., et al., 2024, Nature, 626, 742

## Appendix

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

###  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\approx 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
shape111This is not unexpected due to the time difference; we note that
normally SALT Operations do not consider using spectrophotometric standards
taken more than 5–7 days from the date of observations, 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, $\lambda\sim 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 3).

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.