Astronomy & Astrophysics manuscript no.

main ©ESO 2023

December 6, 2023

Kilonova evolution — the rapid emergence of spectral features

Albert Sneppen1, 2 , Darach Watson1, 2 and James H. Gillanders3

1
Cosmic Dawn Center (DAWN)
2
Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, København 2100, Denmark
3
Astrophysics sub-Department, Department of Physics, University of Oxford, Keble Road, Oxford, OX1 3RH, UK

Received date / Accepted date

ABSTRACT
arXiv:2312.02258v1 [astro-ph.HE] 4 Dec 2023

Kilonovae (KNe) are one of the fastest types of optical transients known, cooling rapidly in the first few days following their
neutron-star merger origin. We show here that KN spectral features go through rapid recombination transitions, with features due to
elements in the new ionisation state emerging quickly. Due to time-delay effects of the rapidly-expanding KN, a ‘wave’ of these new
features passing though the ejecta should be a detectable phenomenon. In particular, isolated line features will emerge as blueshifted
absorption features first, gradually evolving into more pronounced absorption/emission P Cygni features and then pure emission
features. In this analysis we present the evolution of the individual exposures of the KN AT2017gfo observed with VLT/X-shooter
that together comprise X-shooter’s first epoch spectrum (1.43 days post-merger). We show that the spectra of these ‘sub-epochs’
show a significant evolution across the roughly one hour of observations, including a decrease of the blackbody temperature and
photospheric velocity. The early cooling is even more rapid than that inferred from later photospheric epochs, and suggest a fixed
power-law relation between temperature and time cannot capture the data. The cooling constrains the recombination-wave, where a
Sr ii interpretation of the AT2017gfo ∼ 1 µm feature predicts both a specific timing for the feature emergence and its early spectral
shape, including the very weak emission component observed at about 1.43 days. This reverberation analysis suggests that temporal
modelling is important for interpreting individual spectra and that higher cadence spectral series, especially when concentrated at
specific times, can provide strong constraints on KN line identifications. Given the use of such short-timescale information, we lay
out improved observing strategies for future KN monitoring.

1. Introduction porally important earlier spectra were obtained at telescopes in

Chile, Australia and South Africa (Andreoni et al. 2017; Mc-
The temporal data on the kilonova AT2017gfo has provided Cully et al. 2017; Shappee et al. 2017; Buckley et al. 2018),
several key insights into its evolution and composition. First, though not with the spectral coverage of the X-shooter spec-
the bolometric luminosity follows a powerlaw-like decay, ex- tra. While we return to these earlier spectra in a companion
pected for a large ensemble of rapid neutron-capture process paper (Sneppen et al. in preparation), each of these ‘daily’ X-
(r-process) isotopes (e.g. Metzger et al. 2010; Wu et al. 2019). shooter spectra is in fact a composite of several individual ex-
Second, follow-up analyses of the spectroscopic data have iden- posures taken over approximately an hour. For standard reduc-
tified P Cygni features proposed to be formed by r-process el- tions these exposures are generally analysed together, but this
ements, the first and most obvious being Sr ii (Watson et al. is not required, given the brightness of the early kilonova, and
2019), but also potentially La iii and Ce iii (Domoto et al. 2022), the mildly relativistic velocities which smear the spectral line
as well as Y ii (Sneppen & Watson 2023). These P Cygni pro- data over tens of thousands of km/s and obviate the need for
files change across epochs because they trace the receding pho- high S/N at high spectral resolution. Conversely, given the rapid
tospheric velocity. They also potentially probe differing ionisa- kilonova evolution it is informative to examine these exposures
tion states and abundances through the ejecta. This means that, individually. Therefore, in Sect. 2, we re-reduce each exposure
more completely, radiative transfer modelling of the spectrum of the VLT, epoch 1 X-shooter spectrum of AT2017gfo (taken
can constrain changes in abundances and velocities (Gillanders on 18–19 August 2017, 1.4 days post-merger). The reduced sub-
et al. 2022; Vieira et al. 2023a,b). Lastly, the continuum flux epoch spectra are made publicly available and are intended as
(across early epochs) is remarkably similar to a blackbody in an extension and addition to the standard reduction presented in
terms of both normalisation and spectral shape (e.g. Sneppen Pian et al. (2017). In Sect. 3, we examine the physical informa-
2023) with especially the NIR continuum being challenging to tion concealed in the differences between these sub-epochs. In
reproduce for current numerical simulations of mergers (e.g. Do- Sect. 4.1, we show that on timescales of order several hours, re-
moto et al. 2022; Collins et al. 2023). verberation effects dictate that KN spectral features will emerge
The daily Very Large Telescope (VLT) spectra taken with rapidly and as blueshifted absorption components first. Lastly,
the X-shooter spectrograph of the kilonova AT2017gfo associ- in Sect. 4.2 we discuss the broader implications for modelling
ated with the gravitational wave event, GW170817, provide the of including reverberation effects and for the optimal observing
highest-quality information yet available for studying the rapidly cadence and timing for future KNe.
evolving electromagnetic output from the merger of neutron-
stars (Pian et al. 2017; Smartt et al. 2017). Due to observational
constraints, spectra were mostly available with a cadence around
24 hours, limiting the temporal analysis to a time-step which is
comparable to the time since merger for the early epochs. Tem-
Article number, page 1 of 8
 A&A proofs: manuscript no. main

Fig. 1. Top panel: X-shooter spectra of the first (blue) and second (red) pair of exposures from epoch 1. The timeline of exposures is shown in the
the top-right inset). Corresponding best-fit models are shown with dashed lines. The spectral shapes are similar, but the later exposures are less
luminous with a lower blackbody temperature. Shaded bars indicate regions with strong tellurics and poorly-constrained or noisy regions at the
edge of the UVB, VIS and NIR arms of the spectrograph. Lower panel: The flux ratio as a function of wavelength (with the binned weighted mean
overlaid), which shows a large fractional change in flux for the UV-arm, a gradual evolution across the VIS-arm and minor and somewhat noisy
evolution in the NIR-arm. The cooling blackbody prediction is taken from the epoch 1–2 evolution from Sneppen et al. (2023a). The data match
the model reasonably well overall. However, the data seem to require a stronger temperature decline at this early time than the model prediction
based on the overall epoch 1–2 evolution.

2. Sub-epoch reduction and predicted variation EFOSC2 (Smartt et al. 2017); DECAM (Cowperthwaite et al.
2017); LDSS (Shappee et al. 2017); Swope (Coulter et al. 2017)
In the following, we outline the data reduction and observational and VIRCAM (Tanvir et al. 2017). Therefore we follow the con-
difference between the individual exposures of AT2017gfo made vention in Pian et al. (2017), Watson et al. (2019) and Snep-
during epoch 1 (1.43 days post-merger). For the first observ- pen et al. (2023b) and do not artificially scale the spectra to any
ing night, the medium-resolution, ultraviolet (320 nm) to near- specific photometric data as these have significant internal scat-
infrared (2480 nm) spectrograph X-shooter took four exposures ter. Note, the “contemporaneous” photometric fluxes commonly
in nodding mode (Pian et al. 2017), which we label epochs 1.1, used for calibration are actually taken over a timescale of two
1.2, 1.3 and 1.4. With 10 minute exposure-times and additional hours, which should be considered a sizeable systematic effect
telescope overhead, these exposures monitor the evolution of the in the reductions, given the rapid nature of the kilonova evolu-
transient over about one hour. We focus on this early epoch, be- tion at this time.
cause it has the highest ratio of exposure timescale to the to-
tal time passed since the merger. The four exposures are re-
duced in nodding mode for the first two and latter two expo- 2.1. Atmospheric observing conditions
sures (tdiff ≈ 25 minutes, see Fig. 1, top panel). We also reduced On the timescale of the sub-epochs, the atmospheric observing
all four exposures individually in stare-mode, which shows the conditions can vary and affect the reduction output in three no-
same temporal trend as the nodding-mode reduction. However, table ways.
as the NIR background subtraction is difficult in stare-mode, we First, the transmission curves and, by extension, the required
here focus on the nodding mode constraints. telluric corrections, may change. Wavelengths with low atmo-
The spectral reduction pipeline follows those previously de- spheric transmission, such as large parts of the NIR, are most
tailed in earlier publications (Pian et al. 2017; Smartt et al. 2017), susceptible to these changes, and a minor increase in transmis-
where the X-Shooter pipeline supplied by ESO is used for flat- sion, for example, could lead to prominent spikes of the telluric-
fielding, order tracing, rectification, and wavelength calibration. corrected flux. These spikes will be narrow (like the typical
The flux calibration of the X-shooter spectra is broadly con- width of atmospheric absorption lines) and statistically weak, as
sistent with contemporaneous photometry from GROND and the actual observed flux-count is small. To determine the impact
Article number, page 2 of 8
 Sneppen et al.: Kilonova evolution — the rapid emergence of spectral features

of this effect, we find both the typical telluric correction esti- powerlaw relation (e.g. Drout et al. 2017; Waxman et al. 2018):
mated from the standard star (taken at comparable airmass im- !−α
t
mediately following the last exposure) and derive telluric correc- T (t) = T ref (1)
tions for each sub-epoch individually using Molecfit (Smette tref
et al. 2015; Kausch et al. 2015). Molecfit was run using the Here, α is the cooling powerlaw index and T ref is the refer-
standard, well-tested configuration and wavelength-range for fit- ence temperature at some reference time, tref . Fitting blackbody
ting described in Kausch et al. (2015) – albeit not including the models to photometry and spectra from 1 to 6 days post-merger
wavelength-band around 2.35µm in the fit due to the weak signal suggests the observed temperature has a power-law decline with
at these wavelengths. Because the standard star was only taken at α = 0.54 ± 0.02 (e.g. Waxman et al. 2018). However, the tem-
the end of the observing sequence, only Molecfit can constrain perature evolution may be more complex than a single powerlaw
the changes in telluric corrections over the individual exposures. prescription. Indeed, more advanced spectral modelling which,
Ultimately the difference between the two methods proves minor along with a blackbody, also takes account of the 1 µm P Cygni
for the continuum, although fewer artificial and narrow peaks are and the observed NIR emission lines, found a single powerlaw
seen in the intermittent telluric forests redward of 1300 nm when does not accurately describe the temperature evolution (Snep-
using Molecfit. pen et al. 2023a). It is cooling more rapid initially (α ≈ 0.62) at
Second, atmospheric dispersion for objects at low inclina- epochs 1–2 (1.4–2.4 days post-merger), and slower (α ≈ 0.4) at
tion is sizeable and chromatic. However, given the atmospheric epochs 2–3 (2.4–3.4 days post-merger).
dispersion correctors on X-shooter and the parallactic angle of The emitting area evolves subtly between the exposures. The
the slit, we can verify that the spectral trace shows no skew with atmosphere expands radially with R(t) ∝ t for homologous ex-
wavelength in any of the exposures. pansion. However, as the outer surface becomes increasingly op-
Lastly, the seeing conditions were quite favourable across tically thin, the photospheric boundary will recede deeper into
the first epoch with a seeing of 0.58′′ − 0.74′′ at 500 nm. How- the ejecta, so Rph (t) ∝ tβ with β < 1. Fitting blackbody mod-
ever, if the PSF width becomes comparable to the slitwidth, slit- els from 1–6 days post-merger suggests β = 0.61 ± 0.05 (Wax-
loss becomes significant and must be accounted for, otherwise it man et al. 2018). However, the 1 µm P Cygni velocity and the
would lead to an underestimation of the normalisation. Due to Doppler-corrected blackbody-velocity from epochs 1–2 indicate
the wavelength-dependence of seeing, there would also be an in- a smaller recession effect with β ≈ 0.75 (Sneppen et al. 2023a).
creasing slitloss towards shorter wavelengths (which would lead This estimate improves on the previous blackbody velocity by
to underestimation of the temperature). To establish an accurate accounting for the relativistic corrections (see Sneppen 2023)
flux calibration, slitloss corrections were calculated using the and by modelling the observed features in the spectrum (includ-
theoretical wavelength-dependence on seeing with the average ing 1 µm P Cygni and NIR features, see Watson et al. 2019).
seeing FWHM in each exposure (Fried 1966), identical to the Changing the emitting area predominantly produces an achro-
previous reductions (e.g. Pian et al. 2017). However, the seeing matic shift in the normalisation of the spectrum. Additionally,
in epoch 1 is good enough that not accounting for slitloss correc- higher-order chromatic (albeit for this analysis negligible) ef-
tions provides a similar wavelength-dependent spectral evolution fects may follow from a change in the velocity, such as different
between sub-epochs, indicating the robustness of the reduction light-travel time-delays and relativistic corrections (as quantified
to this systematic effect. in Sneppen 2023).
Ultimately, the evolution in epochs 1–2 provides the tempo-
2.2. Predicted kilonova cooling
rally nearest constraints on the predicted spectral evolution over
the first epoch exposures. We will therefore define this model
In a KN atmosphere, the temporal evolution of temperature is as the cooling blackbody prediction, though, a priori, one might
set by the competition between adiabatic cooling, radiative cool- expect a slightly faster evolution during epoch 1 itself than the
ing and heating from the decay of r-process elements. The de- transition from epoch 1–2 would give.
cays from an ensemble of isotopes is predicted to provide a
heating rate well-approximated by a single power-law (Met- 3. Observed sub-epoch evolution
zger et al. 2010), which is observationally supported by the
AT2017gfo lightcurve, where the bolometric luminosity from In Fig. 1, we show the epoch 1 spectra from the nodding-mode
1 to 6 days post-merger follows a temporal (t) powerlaw-like reduction, which clearly shows that the transient is subtly, but
decay, Lbol ∝ t−0.95±0.06 (Waxman et al. 2018). The spectra of statistically significantly, evolving between the exposures. The
AT2017gfo in early epochs is observed to be very similar to drop in flux is particularly pronounced near the blackbody peak
a blackbody as quantified in Sneppen (2023). In this case, for and in the UV, while the two spectra converge towards the redder
a spherical photosphere the wavelength-specific luminosity is wavelengths of the optical.
LλBB (t) = 4πRph (t)2 πB(λ, T (t)), where Rph (t) is the photospheric To constrain and quantify the spectral evolution, we model
radius, B(λ, T ) is the Planck function with temperature, T , at the spectra as a blackbody perturbed with a P Cygni profile from
wavelength λ. The ratio of fluxes at different times for a simple the three strong resonance lines of Sr+ at 1.0037, 1.0327, and
blackbody cooling model is thus merely the ratio of 1) Planck 1.0915 µm following the framework detailed in Watson et al.
functions of different temperatures and 2) emitting areas. We de- (2019) and Sneppen et al. (2023a,b). The flux drop close to the
liberate on each of these in the following. blue edge of the spectrograph at around 350–400 nm is likely
The spectral energy distribution ratio is quite sensitive to the due to Y+ and Zr+ absorption (Gillanders et al. 2022; Shingles
evolution in temperature since we sample the Rayleigh-Jeans tail et al. 2023; Sneppen & Watson 2023; Vieira et al. 2023b). Hence
(i.e. Bλ≫hc/(kb T ) ∝ T ) which has a linear scaling of flux with the spectral shape in this region of the spectrum is degenerate
temperature, the peak (where the flux-ratio for blackbodies is with the modelling of the UV line-opacity. Therefore we only
proportional to T 5 ), and the exponentially growing discrepancy fit the wavelengths λ ≥ 400 nm, where the spectrum is better
between the Wien tails. A commonly used prescription for the constrained. In the following, we consider the evolution of the
temperature, T , as a function of the time t since the merger is a continuum and spectral features.
Article number, page 3 of 8
 A&A proofs: manuscript no. main

Fig. 2. Corner-plot indicating the posterior probability distributions of key parameters from fitting the spectra in epoch 1.1+1.2 (blue) and epoch
1.3+1.4 (red). The fit parameters are (1) the best-fit observed blackbody temperature, T obs , (2) the cross-sectional velocity, v⊥ , derived from the
blackbody normalisation, (3) the photospheric velocity of the 1 µm P Cygni feature, and (4) the outer ejecta velocity inferred from the 1µm P Cygni
feature, vmax . The grey dashed line in the v⊥ vs. v∥ plot indicates the line of equality, i.e. a spherically symmetric velocity surface. Computing the
cross-sectional velocity requires an estimate of the distance to the host, where we here assume DL = 44.2 ± 2.3 Mpc derived from the Planck
cosmology (Planck Collaboration et al. 2020) and the cosmological recession velocity of the host galaxy, zcosmic = 0.00986 ± 0.00049 (Mukherjee
et al. 2021; Sneppen et al. 2023b). We do not propagate the 5% uncertainty in the absolute luminosity distance into v⊥ , as we are interested in the
relative change and uncertainties in fitting the spectra themselves. Additional systematic effects like line-blending could shift the posteriors, but
are likely to affect sub-epochs in a similar fashion.

3.1. Continuum evolution The constraints on the temperature evolution are further
quantified in Fig. 3. The inferred α = −0.77 ± 0.06 from the
In Fig. 1, lower panel, we plot the ratio of fluxes for different sub-epoch evolution suggests a rapid cooling, which is faster
wavelength bins. The wavelength-dependent decline in flux pre- than the 1–2 epoch evolution (α = −0.62 ± 0.01, Sneppen
dicted from a blackbody cooling in time (indicated with grey et al. 2023a) and the average 1–6 days post-merger evolution
shading) traces the UV and optical decline fairly well. We em- (α = −0.54 ± 0.02, Waxman et al. 2018). This suggests that the
phasise this is not a fit to the data, but a prediction of a black- early temperature evolution is steeper than inferred from com-
body with a powerlaw decline of temperature across time de- paring with later epochs. Naturally, the statistical constraining
rived from comparing epochs 1–2 (see Sec. 2.2). On these short power is lower due to the smaller temporal range, but it is pre-
timescales the evolution in flux from 400 nm to 1000 nm broadly cisely this temporal locality which makes this estimate com-
follows the predictions of blackbody cooling. Fitting the black- plementary in nature to the cross-epoch (but temporally sparse)
body temperature and normalisation for each spectrum indepen- analysis. That the powerlaw-slope, α, varies with time is a strong
dently highlights the evolution in ejecta properties (see Fig. 2). indication that a single powerlaw prescription for the tempera-
The best-fit temperature clearly cools between the exposures and ture is an oversimplification.
the emitting area inferred from the blackbody normalisation sug-
gests that v⊥ has decreased.
Article number, page 4 of 8
 Sneppen et al.: Kilonova evolution — the rapid emergence of spectral features

Fig. 3. The sub-epoch evolution of the power-law index of cooling (see Fig. 4. Fraction of Sr in its first ionised state given the time since the
Eq. 2.2) of the X-shooter spectra at 1.43 days post-merger (red his- merger in LTE. The Sr ii is shown as a fraction of total Sr and is tran-
togram). The temperature evolution is highly significantly (> 10σ from sitioning (recombining) from the Sr iii state. The blue lines are based
α = 0). The inferred cooling between exposures, α = −0.77 ± 0.06, on the temperature evolution from the sub-epoch analysis in this paper,
is more rapid than inferred from the average epoch 1–2 evolution while the grey background lines assume the slower cooling in Waxman
(α = −0.62±0.01, Sneppen et al. 2023a), which is itself more rapid than et al. (2018). Dashed and dotted lines show the fraction assuming differ-
the cross-epoch analysis over the period 1–6 days post-merger (Wax- ent electron densities. Regardless of the exact cooling-rate, and across a
man et al. 2018). broad range of electron densities, there is a sharp transition from Sr iii to
Sr ii at 1.4–1.7 days post-merger (or equivalently for an observer the ab-
sorption will form at 1.0–1.2 days). This suggests that 1) spectra taken
That v⊥ decreases may indicate we are observationally prob- just a few hours earlier would have a significantly different Sr ii profile
ing the recession of the thermalisation surface, which drops and 2) the more distant ejecta will have less Sr ii. On the upper axis we
show the time of the X-shooter epochs 1–2 (Pian et al. 2017; Smartt
deeper into the ejecta as the outer layers become increasingly et al. 2017). We note that due to light-travel time delays the more dis-
optically thin. While it is possible that the variable seeing with tant ejecta (near the limb) is observed at a different time to the anterior
non-optimal modelling of the slitloss could produce a similar ejecta as indicated with the grey shaded regions.
shift, nevertheless, we believe the effect we are observing is real,
as there is no other indication of strong slitloss effects and, as we
discuss in Sect. 3.2, the Sr ii line also suggests a recession of the flux over these wavelengths is more pronounced between expo-
photosphere. sures and may indicate increasing opacity at short wavelengths
and thus more reprocessing of the light towards redder colours.
3.2. Spectral features
By later epochs, this absorption feature has disappeared, which
may follow from lanthanide (or d-block element) line-blanketing
The 1 µm line has a P Cygni nature, which allows constraints washing out all spectral features below 600–700 nm (Gillanders
on the velocity of the line-forming region. We show these con- et al. 2022; Sneppen & Watson 2023).
straints as part of the model fit in Fig. 2. The line photospheric The NIR spectra are the most difficult to interpret, partly be-
velocity, v∥ , and the outer velocity of the line-forming region, cause some of the strongest features appear in the NIR at later
vmax , decrease slightly with time. This independently illustrates epochs and may be contributing in epoch 1, and these features
and constrains the photospheric surface recession into the ejecta. are still not well-explained, and partly because the atmospheric
We note here that, similar to the expanding photosphere method conditions, which also change on these timescales, affect the at-
for AT2017gfo presented in Sneppen et al. (2023a), Fig. 2 clearly mospheric transmission curve. However, the decrease in the bulk
shows that the thermalisation radius of the blackbody and the flux is least in the NIR, and seems to be more or less as predicted
photospheric radius of the line decrease coherently and remain from the cooling blackbody framework.
consistent with each other across the exposures within the first
epoch. The constraints presented in Fig. 2, only include the sta-
tistical uncertainty of the fitting model; systematic effects, such 4. Discussion
as line-blending, time-delay or reverberation-effects for v∥ or the
uncertainty in luminosity distance for v⊥ , could potentially shift In this manuscript, we have reduced and compared the sub-
the means of the posteriors, but such effects will likely affect epochs of the first phase of observations obtained with X-shooter
sub-epochs in a similar way, so the temporal evolution of the fit- for AT2017gfo. We have highlighted the physical information
ted parameters still follows the trend suggested from cross-epoch of the kilonova these changes on short-timescale contain, com-
analysis. plementary to the sparse long-time scales between observations.
As mentioned above, the absorption feature below 400 nm These considerations naturally lead to a discussion on the con-
(the UV deficit compared to the blackbody) in the first epochs straints attainable given observations with better cadence. In
has been interpreted as tentative evidence of Y ii and Zr ii lines Sect. 4.1 below we therefore discuss the additional constraints
(Gillanders et al. 2022; Vieira et al. 2023b). While we do not for line-identifications a higher cadence allows and in Sect. 4.2
fit these ranges in the model, we can still see that the change in we discuss how better to optimise future observing strategies.
Article number, page 5 of 8
 A&A proofs: manuscript no. main

Fig. 5. The emergence of a Sr ii P Cygni given different electron densities (left) and time post-merger (right). Left: The X-shooter epoch 1 spectrum
with a model blackbody + Sr ii P Cygni feature overlaid. The wavelength-dependent strength of the Sr ii feature is modulated by the fractional
abundance of single ionised Sr at the time of observation, tobs = 1.43 days (see Fig. 4). There is only a weak dependence on the electron density
as the transition-time is a relatively robust prediction in time. This means that the electron density is not tightly constrained, with the majority of
the feature indicating a moderate-to-high electron density of ne ≈ 108 − 109 cm−3 . Right: given an electron density ne = 109 cm−3 , we show the
corresponding Sr ii P Cygni from 0.9 to 1.8 days post-merger in intervals of 0.1 days. The feature rapidly emerges with the reverberation wave
moving from the blueshifted (and nearer) ejecta towards the red. At the time of observation the absorption is fully formed while the emission is
weak due to the decreased fraction of Sr ii in the emitting region at the time of emission. This means that the Sr ii interpretation makes strong
predictions for the timing of the emergence of the feature and the early spectral shape and its evolution.

4.1. The rapid emergence of spectral lines of magnitude over a few hours. In the radiative transfer equa-
tions, this implies a drastic change in the transmission, I = I0 e−τ ,
Conventionally, line identifications are determined by matching and thus a rapid formation of the line. This immediately suggests
observed features with modelled spectral lines which have con- two temporal predictions of the Sr ii appearance.
sistent transition wavelength and transition probabilities. Thus,
First, spectra taken significantly earlier than about 1.2 day
identifications follow from analysing the wavelength dimension
post-merger should not have a 1 µm P Cygni feature due to Sr ii
of the data, while spectra at different times merely serve as dif-
(as shown in Fig. 4). This fraction of Sr ii is computed assuming
ferent spectral samples for wavelength models. However, for
LTE with the Saha ionisation equation and the temporal depen-
transients with rapid evolution, the timing of the emergence of
dence on temperature in Eq. 2.2 with the cooling-rate, α ≈ 0.77,
spectral features is itself a strong prediction for specific line-
inferred in Sect. 3.1. Applying a slower cooling-rate such as that
identifications both in terms of which ionisation states are domi-
suggested in Waxman et al. (2018) does not significantly change
nant and in terms of being able to deblend lines, as we exemplify
the rapid nature of the recombination1 and the emergence time is
with two of the proposed identifications below.
only relatively weakly dependent on the exact electron density.
In Fig. 5 (right panel), we illustrate how rapidly a Sr ii feature
Sr ii: With a second ionisation energy for Sr of 11.0 eV, typi- would emerge.
cal KN electron densities (i.e. ne ≈ 107 − 109 cm−3 ) and LTE Second, the early spectral shape of Sr ii is very distinct due to
conditions, Sr iii (or even higher-ionised species) will become light travel-time effects. The more distant (earlier and thus hot-
the dominant species for temperatures greater than T ≈ 4300 − ter) parts of the ejecta will have more Sr iii and less Sr ii (Fig. 4).
4700 K. For comparison the emitted blackbody temperature for As the distant ejecta is the origin of the emission peak of the
the ejecta nearest the observer in epoch 1, T (cos(θ) = 1, tobs = P Cygni, and the nearer parts are the origin of the absorption,
1.43 days) = 4150 ± 60 K, while the ejecta further from the ob- this would suggest less emission relative to absorption. That is, at
server (which due to light travel-time is observed at an earlier early times a Sr ii feature would emerge in absorption blueshifted
and hotter time) T (cos(θ) = β, tobs = 1.43 days) = 4900 ± 70 K away from the rest-wavelength of the line, but at increasingly
(Sneppen 2023). Here θ is the angle between the direction of ex- later times as the recombination wave passes through the ejecta,
pansion and the line-of-sight. The transition between ionisation a full P Cygni feature would emerge. This seems to be observa-
states is closely dictated by the temperature because of its ex- tionally suggested by the spectra of epoch 1, where the absorp-
ponential dependency in LTE, while the transition time is much
less sensitive to the exact electron density. As the ionisation-state 1
Or indeed, combination, as it may be the first time these freshly-
transitions, the optical depth, τ, of the line will change by orders formed atoms have ever reached this state.

Article number, page 6 of 8

 Sneppen et al.: Kilonova evolution — the rapid emergence of spectral features

tion dominates and the emission peak is small (see Fig. 5) and the Y ii: There are two series of lines from the 4d2 –4d5p transitions
later epochs where the emission feature grows in prominence. of Y ii, one prominent set of lines centred around 760 nm and
slightly weaker transitions at 670 nm (Biémont et al. 2011). For
In general, given a line-ID (i.e. assuming an ionisation- the early characteristic expansion velocities (i.e. v ≈ 0.3c) the
energy) the timing of the feature emergence allows the elec- features from these components will blend together, where the
tron density to be determined under LTE conditions. Singly- emission and absorption features will produce a series of com-
to doubly-ionised r-process ejecta with Mejecta = 0.05M⊙ , v = plex wiggles in the spectra. However, as velocities decrease these
0.25c and t = 1.43 days would have a mean electron density: lines gradually deblend. When the characteristic velocity recedes
log10 (ne /cm−3 ) = 8.3 − 8.6. Naturally, many systematic effects below 0.2c (around 4 days post-merger) a P Cygni feature should
could shift this estimate substantially, but this density is notably emerge from the 760 nm lines if yttrium is responsible for these
close to that inferred from the 1µm feature given the Sr inter- lines. We emphasise that the time when the feature should first
pretation (see Fig. 5, left panel). Conversely, assuming a broad appear is a prediction specific to the configuration of atomic lev-
range of electron densities one could determine what ionisation els in Y ii. This feature is indeed found to emerge clearly between
energy would be compatible with the line-shape and the forma- epochs 3–4 for AT2017gfo (Sneppen & Watson 2023).
tion time (i.e. putting strong constraints on the spectral line-ID).
The optimal constraints on this recombination wave would be
a series of spectra detailing the feature shift over several hours 4.2. Optimal observing strategy
(see Fig. 5, right panel).

Earlier spectra than the VLT/X-shooter epoch 1 were taken While the timescale of evolution for astrophysical objects, even
with Magellan telescopes using the LDSS-3 and MagE spectro- most transients, is typically much longer than the exposure time,
graphs (at 0.49 and 0.53 days post-merger) and indeed show for fast-evolving transients like KNe this is not always the case.
no 0.7 − 1 µm feature (Shappee et al. 2017). A later spectrum The situation is likely to be even more extreme for earlier iden-
taken using the ANU 2.3m, WiFeS instrument (0.93 days, An- tifications post-merger than for AT2017gfo. Furthermore, given
dreoni et al. 2017) displays only weak or non-existent absorp- the emergence of JWST as a viable tool to produce spectra of
tion, while the South African Large Telescope spectrum (1.18 gamma-ray burst–triggered KNe (Rastinejad et al. 2022; Levan
days, Buckley et al. 2018) has a sizeable absorption component et al. 2023) as well as the improvements to the LIGO sensitiv-
as predicted from Fig. 5. However, spectral modelling is diffi- ity and thus the increased distance of most future KNe, reaching
cult as these spectra have relatively low S/N and only X-shooter a similar signal-to-noise ratio (S/N) as obtained for AT2017gfo
(with its UV-NIR coverage) provides strong constraints on the will likely require longer exposure times, at least before the ad-
continuum shape and the 1 µm emission peak. While this emer- vent of the next generation of extremely large ground-based op-
gence time is thus promising for the Sr ii interpretation, we leave tical telescopes. Thus, in the context of the physical information
a detailed comparative analysis for future work (Sneppen et al., contained within this short timescale perspective, it is natural to
in preparation). It remains to be seen whether a He i interpreta- consider what would be an optimal observing strategy for future
tion of the 1 µm feature, as introduced in Perego et al. (2022) KNe.
and argued for in Tarumi et al. (2023), can reproduce both the First, higher cadence can offer remarkable constraining
observed rapid absorption-to-emission evolution and the time of power for KN modelling. As the temperature and density evo-
appearance, both of which seem to be required by the Sr ii iden- lution is rapid, species cross the threshold between their ionisa-
tification. tion states on the timescale of an hour (i.e. their corresponding
spectral lines will emerge or disappear rapidly), the photosphere
At late times, when the Sr ii species becomes depopulated the recedes deeper within the ejecta, and lines deblend to form inter-
feature will recede and again due to reverberation this will first pretable features (see Sect. 4.1). Thus, even if a future KN due to
be in the blueshifted ejecta (i.e. leaving an emission peak which observational limitations is only briefly observable (such as two
then subsequently fades away). The modelling of the reverber- consecutive hours) it is highly constraining to leverage the full
ation discussed here mainly focuses on the changing ionisation. temporal-span to study the evolutionary trend in parameters.
The change in the strength of the observed lines and the fading of
the spectral continuum with time could produce additional rever- Second, observing blocks should prioritise several shorter
beration effects as discussed in Sneppen et al. (2023b). However, exposures rather than few long exposures, particularly within
modelling the changing continuum only produces a minor effect the first day or two post-merger. The overhead associated with
at early times because of the softer flux evolution towards the increased readout noisemay be relatively minor compared to the
Rayleigh-Jeans–tail of the blackbody. Nevertheless, these con- gain in information from probing the KN evolution between in-
tinuum reverberation effects are essential for understanding the dividual exposures. For instance in the case of AT2017gfo, the
0.7 − 1 µm feature for intermediate-to-late time evolution when four exposures in epoch 1 allow for both temporally-resolved
the continuum shifts redward and the feature’s emission peak nodding and stare-mode reductions. Conversely, the two (albeit
is observed to become dominant (see McNeill et al. in prepa- longer) exposures taken in epoch 2 limits the temporal informa-
ration). Lastly, we emphasise that non-LTE effects, a highly in- tion and removes the possibility of robustness-testa of the reduc-
homogeneous temperature distribution (and potentially a com- tion pipeline across sub-epochs.
plex radial/angular density distribution) could make the tran-
sitions between ionisation states more gradual. Full radiative Third, while spectra derived from averaging several expo-
transfer modelling of merger simulations is required to explore sures for optimal S/N are the standard, early sub-epoch spectra
these dependencies. Conversely, strong observational constraints contain additional insights into the time derivative of KN prop-
on the emergence of lines can thus potentially reveal insights erties and therefore should ideally be reduced and published as
into the homogeneity of temperature and the validity of LTE- well, to provide the full synergy of the time-series spectral evo-
approximations in the ejecta. lution.
Article number, page 7 of 8
 A&A proofs: manuscript no. main

5. Data availability
Work in this paper was based on observations made with Euro-
pean Space Observatory (ESO) telescopes at the Paranal Obser-
vatory under programmes 099.D-0382 (principal investigator E.
Pian), 099.D-0622 (principal investigator P. D’Avanzo), 099.D-
0376 (principal investigator S. J. Smartt). The data are avail-
able at http://archive.eso.org. The re-reduced sub-epoch spec-
tra are made available from: https://github.com/Sneppen/
Kilonova-analysis

6. Acknowledgements
We thank Jonatan Selsing for useful correspondence on the
original reduction pipeline in Pian et al. (2017) and Tom
Reynolds for the implementation of Molecfit. We thank
Stephen Smartt, Stuart Sim, Christine Collins, Luke Shingles
and Kasper Heintz for comments and discussions on the evo-
lution of the recombination-wave. The Cosmic Dawn Center is
funded by the Danish National Research Foundation under grant
number 140. AS and DW are funded by the European Union
(ERC, HEAVYMETAL, 101071865). Views and opinions ex-
pressed are however those of the authors only and do not nec-
essarily reflect those of the European Union or the European Re-
search Council. Neither the European Union nor the granting
authority can be held responsible for them.

References
Andreoni I., et al., 2017, PASA, 34, e069
Biémont É., et al., 2011, MNRAS, 414, 3350
Buckley D. A. H., et al., 2018, MNRAS, 474, L71
Collins C. E., et al., 2023, arXiv e-prints, p. arXiv:2309.05579
Coulter D. A., et al., 2017, Science, 358, 1556
Cowperthwaite P. S., et al., 2017, ApJ, 848, L17
Domoto N., Tanaka M., Kato D., Kawaguchi K., Hotokezaka K., Wanajo S.,
2022, ApJ, 939, 8
Drout M. R., et al., 2017, Science, 358, 1570
Fried D. L., 1966, Journal of the Optical Society of America (1917-1983), 56,
1380
Gillanders J. H., Smartt S. J., Sim S. A., Bauswein A., Goriely S., 2022, MN-
RAS, 515, 631
Kausch W., et al., 2015, A&A, 576, A78
Levan A., et al., 2023, arXiv e-prints, p. arXiv:2307.02098
McCully C., et al., 2017, ApJ, 848, L32
Metzger B. D., et al., 2010, MNRAS, 406, 2650
Mukherjee S., Lavaux G., Bouchet F. R., Jasche J., Wandelt B. D., Nissanke S.,
Leclercq F., Hotokezaka K., 2021, A&A, 646, A65
Perego A., et al., 2022, ApJ, 925, 22
Pian E., et al., 2017, Nature, 551, 67
Planck Collaboration et al., 2020, A&A, 641, A6
Rastinejad J. C., et al., 2022, Nature, 612, 223
Shappee B. J., et al., 2017, Science, 358, 1574
Shingles L. J., et al., 2023, ApJ, 954, L41
Smartt S. J., et al., 2017, Nature, 551, 75
Smette A., et al., 2015, A&A, 576, A77
Sneppen A., 2023, ApJ, 955, 44
Sneppen A., Watson D., 2023, A&A, 675, A194
Sneppen A., Watson D., Bauswein A., Just O., Kotak R., Nakar E., Poznanski
D., Sim S., 2023a, Nature, 614, 436
Sneppen A., Watson D., Poznanski D., Just O., Bauswein A., Wojtak R., 2023b,
A&A, 678, A14
Tanvir N. R., et al., 2017, ApJ, 848, L27
Tarumi Y., Hotokezaka K., Domoto N., Tanaka M., 2023, arXiv e-prints, p.
arXiv:2302.13061
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 N., Ruan J. J., Haggard D., Ford N., Drout M. R., Fernández R., Badnell
N. R., 2023b, ApJ, 944, 123
Watson D., et al., 2019, Nature, 574, 497
Waxman E., Ofek E. O., Kushnir D., Gal-Yam A., 2018, MNRAS, 481, 3423
Wu M.-R., Barnes J., Martínez-Pinedo G., Metzger B. D., 2019, Phys. Rev. Lett.,
122, 062701

Article number, page 8 of 8