11institutetext:  Cosmic Dawn Center (DAWN) 22institutetext: Niels Bohr
Institute, University of Copenhagen, Jagtvej 128, DK-2200, Copenhagen N,
Denmark 33institutetext: GSI Helmholtzzentrum für Schwerionenforschung,
Planckstraße 1, 64291 Darmstadt, Germany 44institutetext: School of
Mathematics and Physics, Astrophysics Research Centre, Queen’s University
Belfast, Belfast, United Kingdom

# No early helium features in the kilonova AT2017gfo; upper bounds from
spectro-temporal NLTE modelling

Albert Sneppen No early helium features in the kilonova AT2017gfo; upper
bounds from spectro-temporal NLTE modellingNo early helium features in the
kilonova AT2017gfo; upper bounds from spectro-temporal NLTE modellingNo early
helium features in the kilonova AT2017gfo; upper bounds from spectro-temporal
NLTE modellingNo early helium features in the kilonova AT2017gfo; upper bounds
from spectro-temporal NLTE modelling 0000-0002-5460-6126 Rasmus Damgaard No
early helium features in the kilonova AT2017gfo; upper bounds from spectro-
temporal NLTE modellingNo early helium features in the kilonova AT2017gfo;
upper bounds from spectro-temporal NLTE modellingNo early helium features in
the kilonova AT2017gfo; upper bounds from spectro-temporal NLTE modellingNo
early helium features in the kilonova AT2017gfo; upper bounds from spectro-
temporal NLTE modelling 0009-0002-5765-4601 Darach Watson No early helium
features in the kilonova AT2017gfo; upper bounds from spectro-temporal NLTE
modellingNo early helium features in the kilonova AT2017gfo; upper bounds from
spectro-temporal NLTE modellingNo early helium features in the kilonova
AT2017gfo; upper bounds from spectro-temporal NLTE modellingNo early helium
features in the kilonova AT2017gfo; upper bounds from spectro-temporal NLTE
modelling 0000-0002-4465-8264 Christine E. Collins No early helium features in
the kilonova AT2017gfo; upper bounds from spectro-temporal NLTE modellingNo
early helium features in the kilonova AT2017gfo; upper bounds from spectro-
temporal NLTE modelling 0000-0002-0313-7817 Luke Shingles No early helium
features in the kilonova AT2017gfo; upper bounds from spectro-temporal NLTE
modellingNo early helium features in the kilonova AT2017gfo; upper bounds from
spectro-temporal NLTE modelling 0000-0002-5738-1612 Stuart A. Sim No early
helium features in the kilonova AT2017gfo; upper bounds from spectro-temporal
NLTE modellingNo early helium features in the kilonova AT2017gfo; upper bounds
from spectro-temporal NLTE modelling 0000-0002-9774-1192

(Received date / Accepted date )

The spectral features observed in kilonovae (KNe) reveal the elemental
composition and the velocity structures of matter ejected from neutron star
mergers. In the spectra of the kilonova AT2017gfo, a P Cygni line at about
$1\,\mu$m has been linked to Sr ii, providing the first direct evidence of
freshly synthesised r-process material. An alternative explanation to Sr ii
was proposed – He i $\lambda 1083.3$ nm under certain non-local-thermodynamic-
equilibrium (NLTE) conditions. A key way to robustly discriminate between
these identifications, and indeed other proposed identifications, is to
analyse the temporal emergence and evolution of the feature. In this analysis
we trace the earliest appearance of the observed feature and compare with a
full collisional-radiative model of helium. We show that the $1\,\mu$m P Cygni
line is inconsistent with a He i interpretation both in emergence time and in
subsequent spectral evolution, where self-consistent helium masses cannot
reproduced the observed feature.

##  1 Introduction

Observations of the gravitational-wave detected kilonova GW 170818/AT2017gfo
(Abbott2017b; Coulter2017) revealed the first spectroscopic identification of
freshly synthesised r-process material in binary neutron star (BNS) mergers
(Watson2019). Studying such r-process features not only constrains the KN
elemental abundances (e.g. Gillanders2022; Vieira2024), but also places
detailed constraints on the geometry and velocity-structures of the ejecta
(Sneppen2023).

Specifically, the 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 Watson2019; Domoto2022) produce the
most prominent absorption, emission or P Cygni features on the kilonova
continuum in the photospheric epochs. The such feature to be robustly
established was a P Cygni feature at 700–1200 nm (Pian2017; Smartt2017),
interpreted as due to Sr ii – first identified in (Watson2019) and later
corroborated with additional radiative transfer modelling (Domoto2022;
Gillanders2022; Vieira2023; Shingles2023) and in non-local-thermodynamic-
equilibrium (NLTE) models (Pognan2023). The strontium mass required to produce
the feature is compatible with that produced in the dynamical ejecta from
merger simulations (Perego2022). Independent evidence for the Sr ii
identification is that strong UV-absorption is observed at $\lambda\leq 400$
nm around 1.4 days after merger, which could be explained by lines from co-
produced elements: Y ii (Gillanders2022) and Y ii and Zr ii (Vieira2023).
Further corroboration of the Sr ii interpretation is the existence of a 760 nm
P Cygni feature interpreted as due to multiple lines from Y ii (Sneppen2023b)
and the timing of reverberation recombination-waves moving across the ejecta
(Sneppen2024). For completeness we also note that second r-process peak-
elements have been suggested, viz. La iii and Ce iii, in order to explain NIR
features at intermediate times (Domoto2022). It is worth noting that the
spectral features observed are of r-process species that are likely populated
under LTE conditions at the observed blackbody temperature.

Figure 1: Evolution of the 1 $\mu$m P Cygni feature over the first 5 days
post-merger from the unified spectral compilation (Sneppen2024). The spectra
were taken post-merger plus 0.92 day (ANU, Andreoni2017), 1.17 day (SALT,
Buckley2018), 1.4, 2.4, 3.4, and 4.4 day (X-shooter, Pian2017; Smartt2017).
The spectra prior to 1 day contain no strong deviation from the blackbody
continuum (see also Shappee2017). However, between 0.92 and 1.17 days, strong
absorption appears, while an emission feature begins to emerge with a delay of
around 10 hours after this. We note that these are the predicted emergence
time-frames for both absorption and emission for a feature produced by LTE Sr
ii due to a recombination wave passing through the ejecta (Sneppen2023_rapid).
We show in this analysis that this emergence is inconsistent with a He i
$\lambda 1083.0$ nm interpretation of the feature.

An alternative interpretation for the 1.0 $\upmu$m P Cygni feature has been
proposed – that it could originate from the neutral He i $\lambda$1083.3 nm
line (1s2s-1s2p transition) which in NLTE conditions can produce a feature at
a similar wavelength (see discussions in Perego2022; Tarumi2023), while Sr ii
may potentially be ionised away by a population of non-thermal electrons from
radioactive heating. It is expected that at later times, e.g. ($t\gtrapprox 5$
days post-merger, NLTE ionisation modelling should become important
(Pognan2022) – and several NLTE features have been suggested for these epochs,
including Rb i (Pognan2023) and Te iii (Gillanders2023; Hotokezaka2023).
Nevertheless, the validity of LTE, alongside the transition point to when NLTE
modelling is required, has so far been difficult to determine observationally.
In this context, helium is a particularly interesting NLTE species to consider
as the required atomic data is largely known and has been applied in
astrophysical context for a wealth of supernovae studies (e.g. Lucy2002).
Furthermore, as helium can be prolifically produced in the nucleosynthethic
limits associated with high-entropies and/or high electron-fraction (such as
the ejecta from neutrino-driven winds expelled during the hypermassive
neutron-star phase, see Just2023; Kawaguchi2023), this implies a helium
feature could constrain the ejection mechanism.

In this paper, we therefore re-visit the NLTE helium proposition to determine
if the observed feature can be self-consistently modelled by NLTE He i in the
early epochs of AT2017gfo (from 0.92–4.4 days post-merger). In Sect. 2, we
briefly summarise the observed evolution of the 1 $\upmu$m feature in the
early epochs. In Sect. 3, we discuss the assumptions required to produce an
observable He i $\lambda$1083.3 nm feature. In Sect. 4, we enumerate the
various relevant transitions between states and implement a full collisional-
radiative model for helium. In Sect. 5, we summarise the predicted evolution
of a helium NLTE feature and show that it is inconsistent with AT2017gfo’s
spectral series. Lastly, in Sect. 6, we discuss the broader lessons this
reveals on the radioactive ionisation and the LTE-likeness of the ejecta.

##  2 The observed evolution of the 1 $\upmu$m feature

To provide the strongest constraints for modelling the 1 $\upmu$m feature, we
need to examine the first appearance and subsequent evolution of the feature.
In Sneppen2024, we have compiled the early spectra of AT2017gfo from Magellan,
ANU, SALT, VLT/X-shooter, and Gemini to provide a higher cadence analysis at
early times than any individual data-series provides (see Fig. 1). This
revealed that the first detected appearance of the 1 $\upmu$m P Cygni feature
is from SALT (1.17 days post-merger, Buckley2018) with an absorption component
at around 800 nm. The spectrum taken 6 hours earlier from ANU (0.92 days post-
merger, Andreoni2017) displays no clear deviation from an underlying blackbody
continuum. Thus the feature emerges rapidly over a relatively short timescale
of a few hours. From 1.4 days post-merger, AT2017gfo was monitored at a
nightly cadence with the X-shooter spectrograph mounted on the European
Southern Observatory’s Very Large Telescope (Pian2017; Smartt2017). The
spectral range (330–2250 nm) provides the first broad spectral coverage
encompassing the near-infrared (NIR) emission peak of the 1 $\upmu$m P Cygni
feature. At 1.4 days this feature has not yet formed in emission likely due to
light-travel time-delay and reverberation effects if the feature is due to Sr
ii (Sneppen2023_rapid). Particularly, the more distant ejecta (where less time
has passed post-merger) has not yet cooled sufficiently to form singly ionised
strontium. Indeed, over the subsequent hours and days, the feature forms in
emission as well as absorption producing a full P Cygni profile. Towards the
later spectra, the characteristic velocities decrease as the photosphere
recedes deeper into the ejecta and the Sobolev optical depth of the absorption
decreases. Post 5 days the feature becomes more emission dominated, which
would be produced by reverberation effects (see Sneppen2023_rapid), but may
also be contributed to by the decreasing densities reducing the optical depth
in the line(s) and allowing radiative de-excitation to form pure emission
features (see Gillanders2023).

For quantitative comparisons of He i $\lambda$1083.3 nm model with the
observed 1 $\upmu$m feature (as undertaken in Sect. 5), we will use the P
Cygni implementation in the Elementary Supernova model (see Jeffery1990)111We
adapt Ulrich Noebauer’s pcygni_profile.py in
https://github.com/unoebauer/public-astro-tools. Here the profile is expressed
in terms of the rest wavelength $\lambda_{0}$, the line optical depth
$\tau(v)=\tau(v_{\rm ph})\cdot e^{-(v-v_{\rm ph})/v_{\rm e}}$, with a scaling
velocity $v_{\rm e}$, a photospheric velocity $v_{\rm ph}$, and a maximum
ejecta velocity $v_{\rm max}$. These parameters can all be optimised to fit
the data with the exception of the rest wavelength of $1083.0$ nm. We note,
this central wavelength is slightly redward (0.02-0.06c) of the observed
emission peak in all photospheric epochs. The parameters related to velocity
can be chosen to mimic the observed feature velocity-scale as the ejecta-
structure is a priori (largely) unknown. In contrast, the optical depth of a
line is a well-defined physical prediction given the oscillator strength and
the density of the relevant level. The observed 1$\upmu$m feature is not
detected at 0.92 days indicating it is very optically thin, is
marginally/fully optically thick from 1.17-3.4 day spectra and then becomes
(increasingly) optically thin from the 4.4 day spectrum and onwards.

The underlying continuum is well-modelled as a Planck function, $B(\lambda,T)$
with temperature, $T$, at wavelength $\lambda$, emitting from a sharply
defined photosphere. At early times, a blackbody provides excellent fit to the
data with percent-level consistency in the required temperature from the UV
through the NIR (Sneppen2023_bb). Towards the later photospheric epochs (see
Fig. 1), NIR emission features emerge at 1.6 and 2.0$\upmu$m, which if not
modelled will bias a simple best-fit blackbody towards lower temperatures than
inferred from the UV-side. As the object is dramatically cooling in early
epochs, the observed feature persists over a broad range of temperatures from
4900-2900 K (1.17-3.4 days), which over the subsequent days then fades away,
while the temperature cools a relatively minor few hundred Kelvin.

##  3 He i $\lambda$1083.3 nm line

The critical aspect in whether He i $\lambda$1083.3 nm (1s2s 3S–1s2p 3P) can
explain the observed feature relies on the population of the lower energy-
level of the transition, which is the lowest-lying level (pseudo ground state)
of the triplet states of neutral helium (see Fig. 2). There are three key
conditions required for producing a He i feature:

  1. 1.

Naturally one needs an appropriate abundance of helium, which, in contrast to
most r-process elements, varies by as much as four orders of magnitude between
different nucleosynthesis simulations. Minor helium abundances can be produced
when modelling dynamical ejecta (such as Perego2022, where $X_{\rm He}\approx
0.01\%$), while high helium abundances are produced in the ejecta material
from the HMNS phase (such as Kawaguchi2022; Just2023, where $X_{\rm He}\approx
1-50\%$).

  2. 2.

As the triplet ground-state is 19.8 eV above the He i ground state (compared
with $k_{B}T\approx 0.25$–$0.4$ eV for early epochs in AT2017gfo) the triplet
population is not efficiently populated by excitation due to thermal
particles. The triplet states are instead populated by recombination from He
ii, which is itself populated by ionisation of He i by radioactive decay
particles. This means that most effective populating of the He i triplet
ground state will occur when the majority of helium is singly ionised.

  3. 3.

The triplet levels are also highly sensitive to photoionisation from near-UV
photons due to their proximity to the ionisation energy ($\lesssim 5$ eV). The
most detailed NLTE helium modelling in the KN context up to this point
(Tarumi2023) only includes ionisation from the triplet ground state, 1s2s 3S
(e.g. photons with energies $>4.8$ eV, $\lambda<260$ nm). However, in this
analysis, we find that ionisation from the 1s2p 3P-state ($\lambda\approx 340$
nm) must be included because the photoionisation and natural decay pathway
leaving the triplet states in all early epochs is _not_ dominated by the 1s2s
3S-state, but by the 1s2p 3P-state, as we will discuss in Sec. 4.3.

We focus in this modelling on the last issue, because this allows quantitative
predictions to be made for the temporal evolution of the He i feature. This
predicted evolution can be compared and contrasted with the predicted
evolution of features due to strontium, and with the 1 $\upmu$m feature as it
is observed. In particular, this aspect is most strongly constrained
observationally in the earliest appearance of the feature, when comparing with
the UV flux constraints from Swift-UVOT (see Sec. 3.3).

###  3.1 Helium yields from r-process nucleosynthesis

By computing synthetic spectra, Perego2022 found that detecting helium
spectral features in kilonovae would be very unlikely for their expected
helium masses (mass fraction, $X_{\rm He}\approx 0.01\%$;
$10^{-5}$–$10^{-6}\,M_{\odot}$) and kilonova luminosities, even when including
NLTE effects. Conversely, Tarumi2023 argue that the 1 $\upmu$m feature could,
in their NLTE modelling, be produced by He i $\lambda 1083.3$ nm (see Fig. 2)
with a relatively modest helium mass fraction of $X_{\rm He}\approx 0.2\%$.
For neutron-rich ejecta, $\alpha$-particles can be produced by radioactive
decay of heavy nuclei ($209<A<250$, Pb–U). However, such pathways can only
produce limited helium with a mass fraction of a few percent. For the ejecta
masses calculated for AT2017gfo (Siegel2019) this corresponds to $M_{\rm
He}\approx 10^{-3}$–$10^{-4}M_{\odot}$. For high electron-fraction or high-
entropy ejecta, $Y_{e}>0.45$ or $s>40{\rm\,k_{B}\,baryon^{-1}}$, the helium
abundances can be increased further through $\alpha$-rich freeze out
(Fernandez2013), although such conditions would quickly imply different
abundances than typically predicted, with, for instance negligible abundances
of elements from the second r-process peak or heavier (e.g. see Fig. 2 in
Vieira2023). However, as we are observing AT2017gfo through the polar ejecta,
which is commonly associated with the high-$Y_{e}$ component (and thus highest
helium abundance) it is indeed possible we are encountering such a
nucleosynthesis regime. Lastly, we note the dynamical ejecta in Just2023
provides orders-of-magnitude smaller helium mass-fractions than found from
post-merger ejecta. These end-to-end NS merger models highlight, that the
helium is predominantly produced in the hyper-massive neutron-star phase (with
its high-$Y_{e}$ ejecta) if the collapse to the black hole takes more than
tens of milliseconds.

While the exact helium mass is thus not known a priori with different ejection
mechanisms producing very different yields of helium, a longer lifetime of the
hyper-massive neutron-star (HMNS) is required for a very high helium
production ($X_{\rm He}\gtrapprox 0.05$). Conversely, a robust upper-limit on
the helium-mass can be used to constrain the lifetime of the HMNS (ie. the
time-scale before black-hole formation).

###  3.2 Radioactive ionisation

The particles released through radioactive decay gradually deposit energy
until they establish equilibrium with the thermal reservoir of the ejecta. The
time-dependent thermalisation of radioactive decay products is well-studied
and can be analytically approximated (e.g. see Barnes2016; Kasen2019;
Hotokezaka2020; Shenhar2024). The ionisation-state of helium provides a
relatively well-defined case-study with experimentally-tested atomic data
(e.g. Nahar2010). In contrast, the ionisation-state of r-process elements in
KNe environments still remain largely unknown due to limited atomic data on
their cross-section for photoionisation, cross sections for non-thermal
collisional ionisation and recombination-rates. In Sect. 6, we return to the
observational constraints from AT2017gfo on the ionisation-state of the
ejecta.

Figure 2: Grotrian diagram showing the He i energy levels. The 1s2s 3S–1s2p 3P
1083 nm transition could produce a 1 $\upmu$m P Cygni feature in NLTE models
(e.g. Tarumi2023), because 1s2s 3S can be difficult to ionise (e.g. requiring
UV-photons $<260\,$nm) and the natural transition to the ground-state is very
slow (e.g. characteristic time-scale of hours). However, as shown here, the
other well-populated triplet state, 1s2p 3P, is several orders of magnitude
more susceptible to ionisation and naturally decays to the ground state
$10^{6}$ times faster, on a timescale of 0.01 seconds. Figure 3: SALT and
VLT/X-shooter spectra with _Swift_ -UVOT photometric constraints on the early
UV flux. Photons with energies sufficient to ionise the population at 1s2s 3P
($\lambda<340$ nm rest-frame or equivalently $\lambda<280$ nm corrected for
ejecta expanding at $0.28c$) is shown by the red region. The UV flux detected
by _Swift_ -UVOT in Evans2017 decreases with time, falling below the 3$\sigma$
(2$\sigma$) detection threshold after 1.1 days (3.0 days). For comparison with
the _Swift_ observations the maximum UV flux allowed so that He i could still
produce the 1 $\upmu$m feature is shown in black – the _Swift_ flux is around
three orders of magnitude larger than permitted for He i to produce a feature
within the Tarumi2023 NLTE helium model.

###  3.3 UV sensitivity of helium triplet states

The He i $\lambda 1083.3$ nm interpretation requires strong limits on the UV
flux. That is, the relevant energy levels (1s2s3S at 19.8 eV and 1s2p3P at
21.0 eV above the ground state) are highly sensitive to photoionisation to He
ii (at 24.6 eV above the ground state). Thus, UV photons (with $\lambda<260$
nm for 1s2s and $\lambda<340$ nm for 1s2p) would quickly depopulate these
triplet states of helium through photoionisation. Tarumi2023 note that 0.3% of
the flux inferred from the Wien tail of a 4400 K blackbody (e.g. approximately
the Doppler-corrected 1.4 day blackbody) is sufficient to depopulate the
triplet states in their models, but given the lack of UV observations it is
tentatively suggested that line-blanketing could suppress the UV flux to this
level (and thus permit a NLTE population in the 1s2s 3S state).

However, near-contemporaneous with the SALT spectrum (and the first
chronological observation of the 1 $\upmu$m feature), there are observational
constraints on the UV-flux taken with _Swift_ -UVOT (Evans2017; Nicholl2017).
We show these constraints in Fig. 3, which suggest i) the _Swift_ data
naturally extends the optical blackbody into the UV (showing no wholesale UV-
blanketing) and ii) the photometric observations are far in excess of the
model-allowed UV flux upper limit, which would thus photoionise and hence,
depopulate, the triplet states. Given the observed decrease of the UV-flux
from 0.6 to 1.1 days (and that we are three orders of magnitude above the UV
upper limit), the UV flux is likely to have been to large to permit a
sufficient He i triplet population to form a strong feature – thus excluding
the He i interpretation – in the first days post-merger. In the subsequent
days, a lower UV flux is tentatively detected – with the last detection (at
$2\sigma$ significance) at 3.0 days post-merger in the _Swift_ UVW2 band.
While extrapolating the _Swift_ data beyond three days would be quite
uncertain, the consistency of the UVOT photometry with a single temperature
blackbody inferred from the X-shooter data in the first few days with no
significant evidence for line-blanketing, gives some confidence that
extrapolating the UV flux in wavelength from blackbody fits to the X-shooter
UVB arm at 4–5 days may be reasonable. Such spectral extrapolation strongly
suggests the UV flux still exceeds the upper limit until at least 4–5 days
post-merger. A _HST_ STIS UV spectrum provides a 3$\sigma$ upper limit of
$F_{UV}\lessapprox 1.5\times 10^{-18}$ erg s-1 cm-2 Å-1 at 250 nm at 5.5 days
(Nicholl2017), which ultimately is not below the limits of the UV-flux
expected from extrapolations of the X-shooter SED.

Figure 4: Network Grotrian diagram indicating the dominant pathways entering
and leaving any level. The arrow-width shows the transition rate on a
logarithmic scale with the colour indicating the dominant mechanism. At early
times the dominant pathways leaving triplet He i is photoionisation. At later
times (as the upper triplet levels are less populated), collisional pathways
out of the triplet become increasingly important.

##  4  Collisional-radiative NLTE model

To accurately model NLTE helium at various epochs, we have implemented a full
collisional-radiative model. In Fig. 4, we show an energy-level diagram and
highlight the various pathways and rates connecting the He i levels. To
produce a strong He i $\lambda 1083.3$ nm (1s2s 3S-1s2p 3P transition) feature
requires the triplet (orthohelium) population i) to be decoupled from the
singlet (parahelium) states, as these rapidly radiatively de-excite to the
(isolated) ground state and ii) to have only weak pathways for photoionisation
to He ii. Crucially, accurately modelling the 1s2s 3S level population must
include higher energy-levels and their transitions, as these contribute
significantly to the pathways depopulating orthohelium. We have collisional
data between all levels with $n\leq 4$, so this is the fiducial model
presented here. We find that including further energy-levels is of lesser
importance given their small NLTE populations, with a less than a 10% decrease
in the triplet population when including all levels up to $n\leq 8$ for all
epochs.

###  4.1 The transition channels

In the following, we enumerate the transition channels between levels. The
Tarumi2023 model includes collisional deexcitation, bound-bound radiative
excitation and photoionisation directly leaving the 1s2s 3S state. The
difference between the work of Tarumi2023 and the model presented here is
solely that we include all the direct and indirect pathways that connect
triplet and singlet states, that is, the inclusion of collisional, radiative,
and photoionisation pathways from the triplet states above 1s2s 3S. As 1s2s 3S
is strongly coupled to these higher triplet levels, the pathways from higher
energy levels are important for accurate calculations (see Fig. 4).

In our model, for ease of comparison, we assume an identical ejecta structure
to Tarumi2023, i.e. a spherically symmetric and homogeneous ejecta with a mass
of $0.04M_{\odot}$, distributed as $\rho=\rho_{0}v^{-5}t^{-3}$ between
$v=0.1c$ and $v=0.5c$, with the value of $\rho_{0}$ set by the total helium
mass in the ejecta. The electron density has the same spatial distribution as
$\rho$, as the free electrons trace the low ionisation r-process elements –
for this we follow Tarumi2023 and assume $n_{e}=1.5\times 10^{8}\,$cm-3 (at
$v=0.3c$, $t=1$ day) as our default model. We have also explored using other
ejecta-structures with the ultimate conclusions on the self-consistency of the
model remaining unchanged (see Sect. 5.1).

The line-forming region is illuminated by a blackbody photon background from
the underlying photosphere. Specifically, we assume a geometrically diluted
blackbody radiation field which is inferred from the Doppler-corrected
continuum. In practice, this means the best-fit blackbody to the continuum is
assumed to set the local radiation field (needed for modelling radiative
transitions) and the electron temperature (needed for modelling recombination
and collisional transitions). The geometric dilution,
$W(v)=0.5\cdot\left(1-\sqrt{1-(v_{phot}/v)^{2}}\right)$, here represents that
the radiation field decreases for velocity shells increasingly distant from
the photosphere. We use the Sobolov escape probability,
$P_{esc}(\tau)=(1+\exp{-\tau})/\tau$, and the blue-wing intensity from the
blackbody model to describe the mean intensity in the line (e.g following the
procedure for large velocity fields in Lucy2002).

####  4.1.1 Radiative transitions (bound-bound)

At all times investigated, radiative transitions within triplet levels
dominate over all other processes and thus set the relative level population.
We model the three different types of radiative transitions, namely
absorption, and spontaneous and stimulated emission.

We have used $A$-values for all transitions extracted from NIST (Kramida2023).
For multiple lines between different fine structure levels, we sum over the
transition’s lower levels and average over all upper states involved
(Axner2004). Notably, the transition strengths between triplet and singlet
states systematically increase for higher energy levels with particularly the
1s4d3D–1s4f1F and 1s4f3F–1s4d1D transitions being very strong, $A\sim
10^{8}\,$s-1 (Drake2006).

We have not included any additional radiative decay mechanisms such as two-
photon emission, since these occur on much longer timescales. In the limit of
detailed balance between two levels, the level population ratios is Boltzmann-
distributed with a geometric dilution factor of 0.5 (as exemplified in Eq. 2,
below). Due to the strong transitions between triplet states this provides a
good approximation of the relative orthohelium populations at all times
investigated.

####  4.1.2 Collisional transitions

In the low-temperature limit transitions from triplet to singlet states are
driven by collisions with thermal electrons. We use the thermally-averaged
transition rates from Berrington1987. When testing the effect of including
higher energy levels ($5\leq n\leq 8$), we have assumed effective collisional
strengths of unity between unlisted levels. However, in practice, collisional
transition rates from high-lying levels, due to the their smaller population,
are negligible when compared to low-lying triplet levels. The rate of
ionisation by thermal electrons is negligible compared to the collisional de-
excitation rate, similar to the findings in previous works (e.g. Tarumi2023).

####  4.1.3 Photoionisation

The wavelength-dependent photoionisation cross-section, $\sigma_{\rm PI}$, for
different levels is taken from Nahar2010. In the first days post-merger,
photoionisation is a significant source of depopulation. For higher levels the
photoionisation cross-section increases and more of the spectral flux can
ionise as the required photon energy decreases. The 1s2s 3S level has the
smallest photoionisation cross-section of any triplet species and only extends
up to 260 nm. 1s2p 3P-level has a photoionisation cross-section that is
$\sim$5 times greater and extends up to 340 nm. The levels from $n=3$ and
above can be photoionised from photons near the blackbody peak, where the
photon density is at its highest. This means that photoionisation of $n=3$
levels is unaffected by UV flux modelling. While 1s2s 3S may have a larger
population then higher energy levels, this level is never a dominant
photoionisation pathway for any observed KN temperature. Photoionisation does
not affect He ii (given the blackbody radiation field observed) due to its
high ionisation potential of $53$ eV.

####  4.1.4 Recombination

For recombination, Nahar2010 provides the temperature-dependent state-specific
recombination coefficients, $\alpha$. Due to the natural transition cascades
from excited levels the excited electron quickly decay to the triplet and
singlet ground-state. For all relevant temperatures this implies approximately
half of recombinations effectively end up populating 1s2 1S and 1s2s 1S, while
the remaining half populates 1s2s 3S.

####  4.1.5 Non-thermal electron ionisation

Nuclear decay creates non-thermal electrons that can ionise helium, thus
producing the He ii required for recombination to populate triplet helium. We
model this process identically to Tarumi2023 and have not revisited their
approximations for evaluating the Spencer-Fano calculations. Specifically, for
ionisation we assume a non-thermal electron flux from Hotokezaka2020 of $1$ eV
s${}^{-1}t^{-1.3}$, with $t$ in days and a work per ion of 593 eV for He i and
3076 eV for He ii (from Tarumi2023 who refer to work in preparation by
Hotokezaka et al.). Non-thermal electrons can ionise any level of He i, but it
is only the limiting rate-step for the ground-state. For all energy-levels
above the ground-state, the radiative and thermal collisional pathways will
dominate. While we thus model it identically for every He i level, it is only
important from the ground-state. For helium the main uncertainty in non-
thermal ionisation lies in the deposition rate, as (to first order) the
ionisation rate coefficients scale linearly with the deposition rate. The
deposition rates are similar across models within a factor of a few, which
implies the dominant uncertainty of the helium ionisation-state comes from the
recombination rate (ie. the exact electron-density).

It is worth noting that the non-thermal modelling leads to the most generous
conditions for the triplet-helium population. At all times considered, He ii
is a dominant species near the photosphere, which maximises the recombination
rate and thus the possible triplet density. If one were to significantly
increase or decrease the radioactive ionisation, He ii could become sub-
dominant to either He i or He iii.

Figure 5:  The fraction of helium in its various ion-states as a function of
electron density (assuming the recombination-rate from (Nahar2010) and the
non-thermal electron ionisation at 4.4 days, see Sect. 4.1.5). The dashed
orange region indicates the range in $n_{e}$, where
$X_{\rm{He\textsc{ii}}}>0.1X_{\rm He}$. For the electron densities expected
near the photosphere ($n_{e}\approx 10^{7}-10^{8}{\rm cm}^{-3}$), He ii should
constitute a major ionisation-state. To deplete He ii would likely require a
several-order-of-magnitude shift from fiducial values of deposition-rates
and/or electron-density. For comparison, we show the corresponding Sr ii
fraction assuming a work-per-ion of 300 eV and a recombination-rate,
$\alpha=3\times 10^{-12}\,{\rm s}^{-1}\,{\rm cm}^{3}$ at $T=3000$ K (C.
Ballance, personal communication). While Sr ii is a sub-dominant ion-state, we
note only minor amounts of Sr ii is needed to produce the observed feature
with $M_{\rm{Sr\textsc{ii}}}\lessapprox 10^{-5}M_{\odot}$ (Watson2019;
Gillanders2022).

###  4.2 Inferred He mass

Using this NLTE population modelling we can estimate and compare the helium
mass needed to produce the observed feature. Specifically, the NLTE model
determines what fraction of helium is in the 1s2s 3S level for any given
epoch, where we infer the radiation field (needed for radiative transition)
and the particle temperature (needed for collision-rates and recombination-
rates) from the Doppler-corrected observed blackbody continuum in each epoch.
Given a total helium-density and time post-merger, the NLTE thus calculates
the 1s2s 3S-density across the line-forming region. To produce a Sobolev
optical depth, $\tau$, given the oscillator strength of the He i $\lambda
1083.3$ nm feature, will require a density n(1s2s 3S)$\approx 7.4\tau
t_{d}^{-1}\,$cm-3 where $t_{d}$ is the time in days (Tarumi2023). The Sobolev
optical depth (at the photosphere) inferred from a best-fit P Cygni to the
absorption feature thus defines the 1s2s 3S-density needed to match the
observed feature. Given the 1s2s 3S-density and in extension (by the
collisional-radiative NLTE model) the total helium density, we can deduce the
total mass by integrating over each velocity-shell:

$M_{\rm He}=\int 4\pi(vt)^{2}\rho\ d(vt),$ (1)

from the inner photospheric shell and outwards. Given the rapid decline of
density with velocity, $\rho\propto v^{-5}$, our mass is dominated by the
ejecta near the photosphere. The resulting masses are detailed in Fig. 6.
Within this framework the dominant uncertainty in the required helium mass
comes from the exact 1s2s 3S level-population relative to the total helium
population – which is dominated by the uncertainty in the ejecta temperature.
We emphasise this ‘helium mass’ is tightly coupled to solely the He ii mass as
this is a dominant ion (see Fig. 5).

In Fig. 7, we show how the He i $\lambda$1083.3 nm feature would evolve across
epochs given various helium mass fractions. The continuum is modelled by the
best-fit blackbody overlayed with the P Cygni parameterised from the
Elementary Supernova model (deliberated in Sect. 2), which has an optical
depth set by the 1s2s 3S-density. The reference total ejecta mass of the high-
velocity ejecta ($v\gtrapprox 0.2c$) is around 0.01M⊙ (Siegel2019). It is
likely that there is even less mass in the line-forming region from 0.92-2.4
days, because the velocity $v\gtrapprox 0.25-0.3c$ (inferred from the
1$\upmu$m line or the blackbody continuum, see detailed velocity-constraints
in Sneppen2023), but as the photosphere recedes towards $0.2c$ this estimate
likely becomes representative of the total ejecta mass around 3-5 days post-
merger.

###  4.3 NLTE helium modelling – sensitivity to 1s2p 3P

This NLTE modelling of helium, assuming a geometrically diluted blackbody
radiation field, indicates that the relative level-populations are well-
described as Boltzmann distributed with a dilution-factor $W\approx 0.5$ for
1s2p 3P (see Bhatia1986; Tarumi2023),

$\frac{n({\rm 1s2p\,^{3}P})}{n({\rm 1s2s\,^{3}S})}=\frac{3}{1}{\rm
exp}\left[\frac{E({\rm 1s2p\,^{3}P})-E({\rm 1s2s\,^{3}S})}{k_{B}T}\right]W$
(2)

For the dereddened blackbody temperatures $T$ for the early 1 $\upmu$m
feature, i.e. 5000 K and 4200 K at 1.17 and 1.43 days, this would imply 10%
and 6% respectively of orthohelium is in the 1s2p 3P-level. We again note that
we are equating in this model the temperature of the local blackbody radiation
field with the observed Doppler-corrected blackbody temperature. However, even
perturbing this assumption, we are almost certainly in the regime that a few
percent of the population is in higher triplet levels.

This subtle minority of the triplet population is important in two respects.
Firstly, the photoionisation cross-section for 1s2p 3P is $\sim$5 times larger
than for 1s2s 3S (Nahar2010) and it can photoionise with redder photons, i.e.
$\lambda\approx 340$ nm, where the flux-level is at least 2–20 times larger,
as we are getting closer to the spectral peak. For the relevant temperatures
of AT2017gfo where the 1 $\upmu$m feature is present, the photoionisation rate
from 1s2p 3P is around 1–2 orders of magnitude larger than from 1s2s 3S. 1s2s
3S and 1s2p 3P are photoionised within respectively $\sim 0.01$ s and $\sim
0.0005$ s by 1.17 days and $\sim 300$ s and $\sim 1$ s by 4.4 days. Secondly,
the natural decay to the ground level from 1s2p 3P (having a lifetime of $\sim
0.01$ s) is a million times faster than 1s2s 3S (which has a lifetime of $\sim
2$ hours). It is thus only at late times, that natural decay becomes important
for estimating the triplet population.

Including even higher triplet energy-levels (with quantum numbers n=3, n=4,
etc.) adds to the depopulation channels, but only significantly at early
times. This is because these higher energy-levels are very weakly populated in
both LTE and NLTE ($\approx 10^{-4}$ less populated than 1s2s 3S at 1-2 days
post-merger and rapidly dropping with time). However, the depopulation rate
from these higher states due to photoionisation is also $O(10^{4})$ times
larger, as optical photons near the spectral peak can ionise and the cross-
section for photoionisation is larger. Additionally, some of their radiative
decays to the singlet states are very large, with $A$-values of order
$10^{6}\,$s-1. Thus, these higher levels are actually comparable in importance
to 1s2p 3P at early times.

Tarumi2023 concluded that a He i feature could be relatively stable when
examining the pathways leaving 1s2s 3S in the limit of no UV-flux. By
including all relevant energy levels, we however find that the He i feature’s
strength is temperature sensitive and will evolve significantly between
epochs.

##  5 Evolution of a NLTE He i feature

Ultimately, a NLTE He i interpretation in early spectra is observationally
contradicted in four respects.

First, the time of the feature’s first appearance is not consistent with He i.
Forming a feature at these epochs requires a helium mass comparable to the
entire mass of the KN higher-velocity ejecta (i.e. the line-forming region).
The required helium mass is lower towards the end of the photospheric epochs
and begins to enter the range of nucleosynthesis calculations (see Fig. 6 and
Fig. 7).

Second, the feature forms too suddenly for He i. Modelling the 1 $\upmu$m
feature as He i at 1.17 days, such a helium-rich model would also create a
feature in the earlier spectra at 0.92 days, which is not observed.
Specifically, a helium mass $M_{He}>0.003M_{\odot}$ is ruled out at $5\sigma$
by the non-existence of a feature in the ANU spectrum from 0.92 days, while a
helium mass of $M_{He}=0.013\pm 0.005M_{\odot}$ is required to fit the
observed feature at 1.17 days (see Fig. 7)).

Figure 6: Estimated helium mass required to produce the absorption of the 1
$\upmu$m feature in early spectra. These estimates are derived using the NLTE
helium model described in Sec. 4. The high velocity ejecta component
constitutes around 0.01–0.02 M☉ for AT2017gfo (Siegel2019). The required
helium mass changes drastically between epochs and for any epoch is largely
inconsistent with all remaining epochs. At early times, the helium mass
required for a feature is comparable to the total mass of the early, higher-
velocity ($v\approx 0.2$-$0.3c$) ‘blue’ ejecta. The absence of a helium
feature at 0.92 days provides a $5\sigma$ upper limit on the helium mass
(reported as solely the statistical uncertainty from the observed spectra).
Figure 7:  The evolution of a NLTE He i $\lambda$1083.3 nm feature with the
vertical line indicating the transition’s rest-wavelength. The helium model of
the 1 $\upmu$m feature is a development of the model of Tarumi2023 with
additions to the original model including all important transitions between
triplet and singlet energy levels (e.g. including radiative absorption and
emission, photoionisation, and electron collision of higher energy triplet
levels, see Sec. 4.3). A NLTE He i feature will become increasingly strong
with time as the dominant pathways depleting the triplet states are from
excited triplet levels which are less populated at lower temperatures. This
highlights the differences in the temporal trends, which show a very weak
initial, but then increasingly strong He line in contrast to the observed
feature that emerges suddenly and then reduces over time.

Third, the higher-triplet level (1s2p 3P) being significantly populated at
early times implies the other triplet He i lines at 587.7 and 706.7 nm, seen
in type-Ib supernovae, should also produce weaker features if a He i 1083.3 nm
line was present. Especially the He i $\lambda 587.7$ nm should be pronounced
in early epoch, as evaluated from the level-population and $A$-values, it
should have around a tenth of the strength of the 1083.3 nm line at 1.17 and
1.43 days. This would fade away to a few percentages by the later epochs. This
is around the typical range of line-ratios seen for type Ib SNe (Dessart2020).
A larger continuum opacity in the optical could suppress the lines, but as
other features are not obscured around these wavelengths (e.g. the absorption
below 400 nm at 1.4 days (Gillanders2022; Vieira2023) and the 760 nm P Cygni
(Sneppen2023b)), such a suppression is unlikely to kill the corroborating
features. Lastly, we note in passing that the 2.06 $\upmu$m line from singlet
helium (1s2s 1S–1s2p 1P) would _not_ be observable due to the rapid radiative
transition to the ground state-level from 1s2p 1P, which keeps these levels
depopulated (Tarumi2023; Gillanders2023).

Fourth, the observed feature strength evolves in the opposite direction
expected in a helium interpretation (see Fig. 7). The triplet 1s2s-level is
increasingly easier to populate at lower temperatures (as the other triplet
states are becoming less populated and these are more sensitive to
photoionisation and radiative decay), so a He i feature should grow in optical
depth with time. In contrast, the observed feature decreases in optical depth.
This also means the He i mass required to produce a feature in any early epoch
is inconsistent with the 1 $\upmu$m feature for any other epoch. The helium
mass required in any early epoch over-predicts the strength in all the
subsequent epochs by 1–2 orders of magnitude – and the lower helium mass
required in later epochs cannot explain the prominence of the feature observed
at earlier times. We note that sizeable radial composition gradients can be
produced in hydrodynamical simulations, but crucially the line-forming region
is overlapping between epochs – making the helium abundance inconsistent
within the same ejecta regions.

To describe the observed feature, we instead require a line which suddenly
appears, remains optically thick and then subsequently gradually becomes
weaker in absorption with time. As showed by Sneppen2023_rapid; Sneppen2024,
all these aspects are predictions of a LTE Sr ii interpretation of the
feature. Specifically, Sr ii in LTE conditions has to i) become rapidly
optically thick at 1.0$\pm$0.1 days in absorption, ii) remain optically thick
for several days, iii) transition to being optically thin at around 3-4 days
post-merger, where-after the feature should decay away. One could potentially
argue, as mentioned in Tarumi2023, that a blend of Sr ii and He i combine to
produce the observed feature with the latter only contributing at the later
epochs. However, such a blend is not observationally required by the data; the
decreasing Sobolev optical depth due to the increasingly dilute Sr ii self-
consistently captures the decreasing strength of the observed absorption. As
such, the results of this temporal analysis naturally favours the conclusion
found in the models of (Perego2022), where strontium rather than helium
produces the observed 1 $\upmu$m feature.

Figure 8:  The helium mass fraction at 1 day from nucleosynthesis calculation
as a function of electron-fractions, $Y_{e}$, and entropy, $s$, given an
expansion time-scale of dynamical ejecta (1 ms, left panel) and neutrino-
driven winds (10 ms, right panel). The nucleosynthesis calculation is from
Gross2023 using the parametric density profiles in Lippuner2015 except that
the full nuclear network is started at 8 GK, instead of 6 GK to account for
the quasi-equilibrium evolution (Meyer1998). The non-existence of a helium
feature suggests low entropy ($s\lessapprox 40\,{\rm k_{B}\,Baryon^{-1}}$) and
provides upper bounds on the electron-fraction ($Y_{e}\lessapprox 0.45$).
These limits become increasingly constraining towards smaller expansion time-
scales.

###  5.1 Generalisations of the ejecta structure

This tension between the absorption of the observed 1$\upmu$m feature’s
evolution (ie. fading after the first few days) and that predicted for a He i
interpretation (ie. generically growing stronger in time) is not strongly
sensitive to the specific ejecta-structure modelled here. We have explored
density-structures $\rho\propto v^{-\alpha}$ with various powerlaw-slopes
$\alpha\in[0;10]$. Furthermore, the baseline model assumes a higher
ionisation-state than predicted for LTE, so we have also investigated a
broader range with $n_{e}=(0.6-1.2)\cdot 10^{8}\,$cm-3 (again at $v=0.3c$,
$t=1$ day), as would be suggested from single-ionised nuclei with $A\sim 100$
and an ejecta mass of $0.04-0.08M_{\odot}$ (Siegel2019). Varying the electron
density or density-structure ultimately does not alleviate the inconsistency
with the observed feature, as this is a consequence of the diminishing
strength of radiative pathways leaving triplet helium.

However, towards the later epoch where recombination and electron collision
begin to dominate the pathways leaving and entering triplet helium, any He ii
mass constraints (as discussed next in Sect. 3.1) are affected by the
electron-density. The highest upper-limit on the He ii mass (ie. the weakest
observational constraints) are provided by the lowest $n_{e}$, because this
weakens the recombination-rate flux entering the triplet He i state. To derive
robust upper limits on He ii-masses, we will assume the lowest electron
density allowed (ie. the remaining species are single-ionised ejecta as is
suggested by LTE for these temperatures). However, NLTE modelling would
suggest that the typical r-process species ionisation-state likely increases
towards later times (Pognan2023), which may suggest larger electron densities,
higher recombination-rates and thus that the He ii mass upper-limits derived
in this analysis are conservative.

###  5.2 Helium mass upper limit and potential implications for
nucleosynthesis

Notably as we approach 5 days post-merger, the blackbody temperature begins to
plateau at around 2500–2800 K (Drout2017; Waxman2018) where the pathways
leaving the triplet states are dominated by electron collisions from 1s2s 3S.
In this limit, the relative triplet population is no longer strongly time-
dependent, so a helium model would be able to produce a self-consistent
optical depth across multiple days. However, at this relatively late time-
frame the observed feature has begun to become optically thin and is in the
process of fading away. As it is observationally unclear what contribution, if
any, the He i $\lambda$1083.3 nm line provides to the 1 $\upmu$m feature, the
helium mass allowed to produce the observed feature may be treated as one form
of upper bound. Any larger He ii abundance would produce a stronger absorption
feature than observed at this time (see Fig. 7).

It is in these later epochs, when the upper limit of the allowed He ii mass is
at its lowest, that a helium interpretation provides the strongest constraints
on the ejecta properties. Specifically, $M_{\rm{He\textsc{ii}}}\lesssim 4\cdot
10^{-5}$ M☉ from the optical depth of the absorption at 4.4 days (when we are
looking at a photospheric velocity of $\sim 0.15c$) would suggest
$X_{\rm{He\textsc{ii}}}\lesssim 0.004$ if one compares with the high-velocity
ejecta mass estimate, $M\approx 0.01\,$M☉. It would be
$X_{\rm{He\textsc{ii}}}\lesssim 0.004-0.001$ if we assumed a total ejecta mass
of $M=0.06\,$M☉ with a density structure from $\rho\propto v^{-7}$ to
$\rho\propto v^{-3}$. To turn this into a total helium mass limit, we need to
know the fractional distribution of helium ionisation-states. The baseline-
model (see Sect. 4) would suggest a maximum total helium mass-fraction,
$X_{\rm He}\lesssim 0.025$. Across a broad range in electron-density,
$n_{e}/{\rm cm}^{-3}\in[10^{6};10^{9}]$, more than 10 % of the helium should
be in He ii, which would indicate constraints $X_{\rm He}\lesssim 0.04$.

Similar to the analysis in Domoto2022, where the non-detection of Ca ii lines
(and thus its relative mass to Strontium) provided nucleosynthesis constraints
for the kilonova, these He ii upper-limits can inform us on the
nucleosynthesis conditions (under the assumption He ii is a significant
ionisation-state). Specifically, these conditions are within nucleosynthetic
yield predictions from heavy element radioactive decay but does not suggest
substantial $\alpha$-rich freeze-out (see Fig. 8). These helium mass upper
bounds are consistent with an electron fraction $Y_{e}<0.45$ and $s<50{\rm
k_{B}\,baryon^{-1}}$ even for the typical expansion-timescale of neutrino-
driven winds ($\sim$10 ms). If single-ionised helium is more heavily
suppressed than generically suggested from the radioactive ionisation
modelling, these constraints can be weakened. For $X_{\rm{He\textsc{ii}}}\sim
0.05X_{\rm He}$, this would suggest looser constraint $Y_{e}<0.47$ and
$s<100{\rm k_{B}\,baryon^{-1}}$. We emphasise that these upper-bounds on the
helium mass only directly constrain the abundance from absorption. The near-
polar viewing angle of the system means we are in effect only constraining the
helium abundance in the regions close to the pole from the absorption.
However, it is also the polar ejecta, with its high $Y_{e}$ and high-entropy
tail, that is commonly interpreted to have the highest helium compositions
(e.g. Fernandez2013; Just2015; Fujibayashi2018).

##  6 Is NLTE needed for modelling KNe at early times?

Perhaps the most striking aspect of the early spectra taken of AT2017gfo is
their apparent simplicity. A simple blackbody empirically describes the
continuum very well at early times (e.g. the first X-shooter spectrum at 1.4
days has percent-level consistency in the inferred temperature from the UV
through the NIR, see Sneppen2023_bb). The spectral perturbations away from the
blackbody are produced in the KN atmosphere by the r-process species with the
strongest individual transitions - with LTE appearing to be a good assumption
for the relevant ionisation-state of those species (Watson2019; Domoto2022).
Indeed, the emergence time and early spectral shape of features in AT2017gfo
follow the predictions of a recombination wave between ionisation states in
LTE, which constrains the blackbody temperature of the radiation field to be
within 5% of the electron temperature (Sneppen2024). While both the continuum
and the spectral features are thus well-described by LTE conditions and LTE
populations at early times, it remains unclear to what extent NLTE effects are
required and indeed, permitted, by the early data.

We note further complexities could be explored in future models, such as
variety of ejecta-structures, including more refined radioactive ionisation
modelling, and allowing different excitation temperatures beyond the observed
colour temperature. However, as noted above, the NLTE-model developed by
Tarumi2023 and extended here already assumes optimistic conditions for forming
a He feature. That it still cannot describe the observed feature is because
one would need to produce a NLTE model that populates the 1s2s 3S-level
significantly without populating any other excited levels of He i. This seems
difficult, since one would have to decouple 1s2s 3S from 1s2p 3P, and these
are the relevant lower and upper levels of the 1083 nm transition
respectively.

One cannot alleviate the unphysical helium mass required at early times by
introducing stronger (or weaker) radioactive ionisation because the models
presented by Tarumi2023 and here already assume the majority of helium is
singly-ionised and thus already assume the maximum recombination rate
possible. Indeed, any significant change to the radioactive ionisation
modelling would only decrease He ii, thus decreasing the recombination rate to
He i and increasing the helium mass tension. Tarumi2023 note that non-thermal
electrons should deplete Sr ii within their framework in favour of more
highly-ionised species with time. Conversely, the NLTE modelling in Pognan2023
still finds substantial contribution from Sr ii even for late-time spectra.
Ultimately, current ionisation-state modelling is highly uncertain for
r-process elements. Due to the lack of photoionisation cross sections for
r-process elements, the recombination rates are simply scaled hydrogenic
approximations in Tarumi2023, and non-thermal collisions are commonly
approximated with generic formulae for ionisation cross sections (see
Pognan2022). Improvements to the atomic data through detailed theoretical
calculations and experiments will be instrumental to permitting stronger
statements based on the observed data.

Nevertheless, it is worth considering what the observed consistency of the 1
$\upmu$m feature (and its continued relation to Sr ii) suggests for modelling
radioactive ionisation, because models which predict that the Sr ii population
should substantially diminish in time (due to radioactive ionisation) do not
describe the observed feature’s evolution well. For context, we note that
detonation models of type Ia SNe predict that non-thermal ionisation leads to
the presence of multiply-ionised species and the depletion of low-ionisation
species, which again is in contrast to the observed emission lines from SN Ia,
which are dominated by singly- and doubly-ionised Fe-group nuclei (Wilk2018;
Shingles2022). Such discrepancies may potentially be alleviated by introducing
clumping (Wilk2020), while it may also hint at broader tensions between
current non-thermal ionisation modelling and the observed spectra.

##  7 Conclusion

In this analysis, we have revisited whether He i $\lambda 1083.3$ nm in NLTE
can explain the observed $1\,\mu$m P Cygni feature in the kilonova AT2017gfo.
Employing a full collisional-radiative model, we found that the density of
triplet helium (and thus any He i $\lambda 1083.3$ nm feature) is highly
temperature sensitive. Due to the decreasing importance of radiative
transitions a helium feature will generically become stronger in time, which
is in contrast to the decreasing Sobolev optical depth of the observed
feature. Thus, the helium mass required to produce the $1\,\mu$m feature is
inconsistent between each observed epoch.

The earliest appearance of the feature is in the SALT spectrum at 1.17 days,
where the helium mass required in the line-forming region, $M_{\rm He}\sim
10^{-2}\,{\rm M}_{\odot}$, is comparable to the entire high-velocity ejecta
component (i.e. $X_{\rm He}\sim 1$), which effectively excludes a He i
interpretation from a nucleosynthesis perspective on mass grounds alone at
early times. This interpretation is further disfavoured by the non-detection
of a corroborating features (particularly He i $\lambda 587.7$ nm), the
emergence-time being more rapid than suggested by NLTE feature’s formation and
the subsequent fading in the line. While the Sr ii model can reproduce the
data at all times, the increased strength of a He i line at later epochs
allows us to use this as a limit on the maximum He ii abundance in the
kilonova. By 4.4 day spectra the required He ii mass is
$M_{\rm{He\textsc{ii}}}\lessapprox 4\cdot 10^{-5}\,{\rm M}_{\odot}$
($X_{\rm{He\textsc{ii}}}\lesssim 0.004$). If sizeable helium is single ionised
($X_{\rm{He\textsc{ii}}}\gtrapprox 0.05X_{\rm He}$), this put strong limits on
the mass of the high-$Y_{e}$ component.

## Data availability

The spectral series used in this paper is composed from a series of different
observing programmes at various telescopes. 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, see Sneppen2024)
now publicly 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 compiled data-sets (and re-processed
when necessary, see Sneppen2024) can be found at
https://github.com/Sneppen/Kilonova-analysis.

## Acknowledgements

The authors would like to thank Kenta Hotekezaka for clarifying and discussing
the NLTE helium modelling in previous literature.

The Cosmic Dawn Center (DAWN) is funded by the Danish National Research
Foundation under grant DNRF140. AS, RD, DW, CEC and SAS are funded/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. LJS acknowledge support by the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program (ERC Advanced
Grant KILONOVA No. 885281). LJS acknowledges support by Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) - Project-ID
279384907 - SFB 1245 and MA 4248/3-1.


# A Lindblad master equation capable of describing hybrid quantum systems in
the ultra-strong coupling regime

Maksim Lednev  Francisco J. García-Vidal  Johannes Feist  Departamento de
Física Teórica de la Materia Condensada and Condensed Matter Physics Center
(IFIMAC), Universidad Autónoma de Madrid, E-28049 Madrid, Spain

###### Abstract

Despite significant theoretical efforts devoted to studying the interaction
between quantized light modes and matter, the so-called ultra-strong coupling
regime still presents significant challenges for theoretical treatments and
prevents the use of many common approximations. Here we demonstrate an
approach that can describe the dynamics of hybrid quantum systems in any
regime of interaction for an arbitrary electromagnetic (EM) environment. We
extend a previous method developed for few-mode quantization of arbitrary
systems to the case of ultrastrong light-matter coupling, and show that even
such systems can be treated using a Lindblad master equation where decay
operators act only on the photonic modes by ensuring that the effective
spectral density of the EM environment is sufficiently suppressed at negative
frequencies. We demonstrate the validity of our framework and show that it
outperforms current state-of-the-art master equations for a simple model
system, and then study a realistic nanoplasmonic setup where existing
approaches cannot be applied.

Light-matter interaction in the strong coupling regime in which matter and
electromagnetic (EM) modes hybridize has enabled the manipulation of the
physical and chemical properties of hybrid light-matter systems at a quantum
level [1, 2, 3, 4, 5, 6, 7]. Typically, these are dissipative systems in which
the interaction with an external reservoir introduces irreversible dynamics
such as decay of excitations by photon emission. Such effects are customarily
treated through Lindblad master equations, in which baths are represented by
Lindblad dissipation terms [8, 9, 10]. These terms typically act on the
uncoupled components to, e.g., represent losses of a cavity mode due to
leakage through the mirrors. While it has long been known that decay operators
derived for an uncoupled system can lead to unphysical effects in the coupled
system [11], their use often remains a reasonable approximation. However, this
fails in the ultra-strong coupling (USC) regime where the coupling strength
becomes a significant fraction of the system transition frequencies, which has
been achieved in many physical systems ranging from organic molecules in
Fabry-Pérot cavities to superconducting qubit-oscillator circuits [12, 13, 14,
15, 16, 17, 18, 19, 20, 21, 22]. In this regime, the commonly used rotating-
wave approximation for light-matter interaction breaks down, leading to
entangled ground states with virtual excitations and opening new opportunities
for nonlinear optics [23, 24]. Decay operators in the uncoupled basis then
introduce unphysical effects such as artificial emission from the ground state
since they act on the virtual excitations [25, 26, 27].

In order to mitigate these problems, decay operators acting in the coupled or
dressed basis have been derived [26, 28, 29, 30]. However, the price to be
paid is that the system Hamiltonian has to be diagonalized and that the decay
operators become significantly more complex. The use of such approaches is
thus restricted to simple cases where the dynamics is dominated by a single
lossy cavity mode and the emitter has limited structure, while at the same
time the use of few-mode and few-state approximations in the ultrastrong
coupling regime becomes questionable [31, 32, 33]. Finally, these methods
still employ an underlying Born-Markov approximation for the coupling of the
system (emitter and discrete mode) to the outside environment and are in this
respect similar to the standard Bloch-Redfield (BR) master equation of open
quantum systems [8, 9].

Figure 1: A realistic spectral density at zero temperature and its Lorentzian
approximation, corresponding to a single-mode Lindblad master equation. The
tail at negative frequencies causes artificial pumping of energy into the
emitter, leading to incorrect dynamics.

In this Letter, we introduce a method for treating hybrid quantum systems in
the USC regime that works for arbitrary EM environments, is independent of the
properties of the quantum emitter, and does not make any Markovian
approximation, while still maintaining the simple form of a Lindblad master
equation in which the dissipators are simple decay terms for the “cavity”
modes. We achieve this by extending a recently developed approach [34, 35] to
the case of extreme coupling between constituents or with external baths. The
crucial step relies on the fact that the failure of the conventional approach
can be understood from the perspective of open quantum systems theory [9, 36,
37] by realizing that a single lossy cavity mode corresponds to an effective
environment that contains negative-frequency components (see Fig. 1). We
demonstrate that for an EM environment consisting of interacting modes [34,
35], their interference in the coupling to the emitter can suppress these
negative-frequency components, making it possible to avoid the artificial
effects inherent in the standard Lindblad master equation. We consider several
limiting cases, benchmarking our approach by comparing with numerically exact
solutions, and compare it to current state-of-the-art methods. Finally, we
illustrate the capability of our method to go beyond simple ad-hoc model
systems by considering the dynamics of an emitter placed in the hot spot of a
nanoplasmonic structure formed by a dimer of silver spheres.

We start from the Hamiltonian describing the interaction of an emitter with
its EM environment (with $\hbar=1$ here and below) [38, 39]

$H=\frac{\omega_{e}}{2}\sigma_{z}+\int\limits_{-\infty}^{\infty}(\omega
a_{\omega}^{\dagger}a_{\omega}+\sqrt{J(\omega)}(a_{\omega}^{\dagger}+a_{\omega})\sigma_{x})\mathrm{d}\omega,$
(1)

where for simplicity we use a two-level emitter described by Pauli operators
$\sigma_{x}$ and $\sigma_{z}$, with transition frequency $\omega_{e}$. The
operators $a_{\omega}$ ($a_{\omega}^{\dagger}$) are bosonic annihilation
(creation) operators for photon modes at frequency $\omega$. The spectral
density $J(\omega)$ encodes the full information about the EM environment and
its interaction with the emitter [36, 37, 39].

We apply a recently developed few-mode quantization method based on replacing
the EM system with an equivalent model consisting of $N$ lossy and interacting
modes [34, 35]. Its dynamics are given by

$\displaystyle H_{\mathrm{mod}}$
$\displaystyle=\frac{\omega_{e}}{2}\sigma_{z}+\sum_{i,j}\omega_{ij}a_{i}^{\dagger}a_{j}+\sum_{i}g_{i}(a_{i}^{\dagger}+a_{i})\sigma_{x}$
(2a) $\displaystyle\dot{\rho}$
$\displaystyle=-i[H_{\mathrm{mod}},\rho]+\sum_{i}\kappa_{i}L_{a_{i}}[\rho],$
(2b)

where $L_{O}[\rho]=O\rho O^{\dagger}-\frac{1}{2}\\{O^{\dagger}O,\rho\\}$ is a
Lindblad dissipator, $\omega_{ij}$ encodes the mode energies and couplings,
$\kappa_{i}$ their decay rates, and $g_{i}$ their coupling to the emitter.
This model is exactly equivalent to Eq. (missing)1 with spectral density
$J_{\mathrm{mod}}(\omega)=\frac{1}{\pi}\vec{g}\cdot\operatorname{Im}[(\tilde{\textbf{H}}-\omega)^{-1}]\cdot\vec{g}$,
where $\tilde{\textbf{H}}_{ij}=\omega_{ij}-\frac{i}{2}\delta_{ij}\kappa_{i}$
[34, 35]. By varying parameters $\omega_{ij}$, $\kappa_{i}$, and $g_{i}$, the
model can be adjusted to obtain a spectral density as close as desired to the
original problem by performing a nonlinear fit of $J_{\mathrm{mod}}(\omega)$
to the physical spectral density $J(\omega)$ for a sufficient number of modes
$N$.

The single-mode case $N=1$ of Eq. (missing)2 is the conventional quantum Rabi
model [40, 41], and thus corresponds to a Lorentzian spectral density
$J_{\mathrm{lor}}(\omega)=\frac{g^{2}}{\pi}\frac{\kappa/2}{(\omega_{c}-\omega)^{2}+\kappa^{2}/4}$
[42, 43, 34], where $\omega_{c}$ is the mode frequency and $\kappa$ its decay
rate. In contrast to physical spectral densities, $J_{\mathrm{lor}}(\omega)$
is non-zero along the whole real axis, including for negative frequencies.
Since counterrotating coupling terms $a\sigma^{-}$, $a^{\dagger}\sigma^{+}$
are resonant at negative frequencies (where $\sigma^{\pm}$ are Pauli jump
operators, with $\sigma_{x}=\sigma^{+}+\sigma^{-}$), and the emission of
negative-energy quanta corresponds to the absorption of energy by the system,
this viewpoint provides a simple intuitive explication for the artificial
pumping observed when cavity decay is described with a Lindblad term and the
rotating-wave approximation is not performed. At the same time, it provides a
natural recipe for preventing such effects: Ensuring that the spectral density
at negative frequencies is sufficiently small. While this is not possible for
noninteracting modes (corresponding to a sum of Lorentzians), we show below
that interactions can enable destructive interference between the modes that
allows this to be achieved even with a relatively small number of modes.

As a first test case, we study a single-mode setup corresponding to a
physically allowed extension of the quantum Rabi model. This is obtained by
coupling a single mode to an Ohmic “background” bath. The effective spectral
density of the full EM environment can then be analytically obtained as

$J_{\mathrm{sm}}(\omega)=\theta(\omega)\frac{2g^{2}}{\pi}\frac{\kappa\omega_{c}\omega}{(\omega_{c}^{2}-\omega^{2})^{2}+\kappa^{2}\omega^{2}},$
(3)

where $\theta(\omega)$ is the Heaviside theta function. This expression
fulfills the physical constraints for EM spectral densities: It only contains
positive-frequency components and it tends to zero for $\omega\to 0$. We note
that it can also be obtained as an antisymmetrized extension of the Lorentzian
spectral density, i.e.,
$J_{\mathrm{sm}}(\omega)=\theta(\omega)(J_{\mathrm{lor}}(\omega)-J_{\mathrm{lor}}(-\omega))$
(with a small renormalization of the parameters). We choose parameters typical
for Landau polaritons formed in semiconductor quantum wells in the USC regime
[24, 21, 22], with values $\omega_{c}=\omega_{e}=0.58$ meV, $g=0.25$ meV, and
$\kappa=0.1$ meV. The corresponding spectral density is shown in Fig. 2(a),
together with fits using either a single Lorentzian or 10 interacting modes
and including the additional suppression at negative frequencies. We compare
various approaches to calculate the population dynamics of an initially
excited emitter in Fig. 2(b): (i) The single-mode Lindblad master equation,
(ii) a BR master equation (for which the Ohmic bath coupling to the “main”
mode is treated perturbatively), (iii) the generalized master equation (GME)
introduced in Ref. [30], and (iv) our approach with a collection of
interacting modes. As expected, the single-mode Lindblad equation
significantly overestimates the steady-state population due to the presence of
artificial pumping. Comparison with “exact” (i.e., numerically converged)
results obtained through direct discretization demonstrates that our approach
produces converged results within the linewidth of the figure. Both the GME
and BR approaches produce similar results that are relatively close to the
exact ones, but show significant deviations at later times. We attribute this
to the Markov approximation for the coupling between the main cavity mode and
its background bath inherent in these approaches.

Figure 2: (a) The spectral density of the system under study. Black dashed
line represents the preset threshold value for model spectral density at
negative frequencies. Grey lines show real part of the complex resonances of
the model spectral density, and the brown dotted line indicates the emitter
transition frequency. (b) The population of the two-level emitter interacting
with photonic environment in the USC regime for different methods. Figure 3:
Average relative error of the emitter population dynamics obtained with our
approach for different number of modes involved into fitting. Red and green
arrows show the error of GME and BR master equations. In inset: time
dependence of relative error. Orange, red and green lines represent result for
our approach (with 10 modes at the threshold value = $10^{-8}$ meV), GME and
BR master equation, respectively.

In order to quantify the accuracy of the different methods, in Fig. 3 we show
the relative error in the emitter population (compared to the “exact”
numerical solution) for the different approaches. The BR and GME methods
contain no free parameters and thus provide no way to systematically improve
the approximation. Their time-dependent relative error in the emitter
population is shown in the inset of Fig. 3, with its value averaged over the
propagation time indicated by the arrows in the main figure. In contrast, the
fitting-based approach introduced here permits to choose a tradeoff between
accuracy and complexity of the model by choosing the number of modes as well
as the threshold value below which we require the spectral density at negative
frequencies to be suppressed. The behavior of the method is demonstrated for
the cases of $N=3$, $N=5$, and $N=10$ modes. When the threshold is set
relatively high, the error is dominated by the artificial pumping and somewhat
independent of the number of modes used for fitting, although we note that
even for the largest threshold values considered here, the errors in our
approach are comparable to the ones of the GME and BR methods. Decreasing the
threshold initially leads to a significant reduction in the error for any
number of modes. However, the error increases again at some point as the
threshold is decreased, essentially because the fit quality in the positive-
frequency components cannot be maintained when the constraints at negative
frequencies are too stringent. This can be mitigated by increasing the number
of modes in the fit, with $N=10$ modes providing enough flexibility in the
current case to maintain accuracy even at the smallest threshold values we set
($J(\omega)\leq 10^{-8}$ meV for $\omega<0$). We note that even for 3 modes,
the optimal accuracy of our approach is considerably higher than that of the
BR and GME methods, while providing the additional advantage that the
dissipators in the Lindblad master equation act only on the cavity modes and
no diagonalization of the full Hamiltonian is required.

Next, we demonstrate the implementation of our model in a nanophotonic
structure with a more complex spectral density: a dimer of silver spheres
(radius 15 nm) with a 1 nm gap between them, embedded in a GaP matrix
($\epsilon=9$). We study the excited state population of a two-level quantum
emitter placed in the center of the gap between the spheres. The emitter has a
transition frequency $\omega_{e}=2.4$ eV and dipole moment $\mu=20$ D oriented
along the line connecting the centers of the spheres, which are typical
parameters for organic molecules or quantum dots. Using the boundary element
method implemented in SCUFF-EM [44], we computed the dyadic EM Green function
of the structure, which determines its spectral density [39].

Figure 4: (a) The spectral density of the considered system. Black dashed line
represents the preset threshold value for model spectral density at negative
frequencies and the brown dotted line indicates the emitter transition
frequency. (b) Emitter population as a function of time for numerically exact
and our approach applied for narrow ($(0.7;4)$ eV) and broad ($(-5;5)$ eV)
spectral windows. In inset: The total number of photons accumulated in the
system and bath modes.

The spectral density, shown in Fig. 4(a), consists of several localized
surface plasmon resonances of various multipolar orders. This already implies
that the number of required modes for an accurate fit will be significantly
larger than in the simplified models treated up to now. The high-permittivity
dielectric background red-shifts these resonances compared to the free-space
situation, making the fitting procedure for our approach more challenging
since the target spectral density shows a significant gradient close to zero
frequency. In order to obtain an accurate fit for the spectral density with
sufficient suppression at negative frequencies, we used 28 modes. Since in
this case, it is not possible to compare with simplified master equations or
single-mode models, we instead compare with a fit where the negative frequency
range is ignored and only the resonance peaks (in the frequency range between
$0.7$ eV and $4$ eV) are fitted. For this simpler situation, the fit is
already well-converged with 12 modes.

The emitter frequency $\omega_{e}$ is resonant with the pseudomode formed by
overlapping high-order surface plasmon modes [45, 46], i.e., close to the peak
of the spectral density. We show the emitter population after initial
excitation in Fig. 4(b). After several oscillations due to interaction with
the complex EM environment, it gradually decreases to zero while undergoing
several further oscillations. Our method again provides essentially perfect
agreement with the numerically exact discretization method. Moreover, we note
that it is essential to include the suppression of the spectral density at
negative frequencies to obtain converged results when counter-rotating
coupling terms are included. The fit taking into account only the positive-
frequency region from $0.7$ to $4$ eV leads to considerable deviations in the
dynamics. As in the previous cases, the differences between the results of our
approach for broad and narrow spectral ranges are caused by artificial pumping
and, notably, are most pronounced at long times. In contrast, the oscillatory
dynamics at short times is represented correctly in both cases since it is
principally determined by the interaction between the emitter and the resonant
region of the spectral density.

In addition to emitter population, in the inset of Fig. 4(b), we show the
time-dependent total number of photons in the system and bath modes. The exact
dynamics demonstrates that the steady-state photon number slightly exceeds
unity. This is a consequence of the importance of the counter-rotating terms,
which lead to the (virtual) excitation of the cavity modes in the coupled
system even in the ground state. Since the decay operators in our approach are
of a Lindblad form, the population of the bath modes can be easily tracked as
$P_{bath}(t)=\sum_{i=1}^{N}\int_{0}^{t}\kappa_{i}\langle
a_{i}^{\dagger}a_{i}\rangle(t^{\prime})dt^{\prime}$. If the fit is performed
only at positive frequencies without taking the suppression at negative
frequencies into account, this results in a linear increase in the photon
number at long times due to the continued artificial pumping through emission
of (negative-frequency) photons from the ground state. In contrast, the result
of our approach for the extended spectral region perfectly follows the trend
of the numerically exact results. This shows that even though the discrete
modes in the fit are in some sense arbitrary, their combination not only
provides the correct emitter dynamics, but also correctly reproduces the
temporal behavior of the photon number.

We furthermore note that the steady state reached by the hybrid system becomes
essentially pure when the artificial pumping is sufficiently suppressed. In
the case of the extended spectral range, the steady-state density matrix
$\rho_{s}$ is highly pure, $1-\Tr(\rho_{s}^{2})<10^{-3}$, meaning that
$\rho_{s}$ can be expressed to a good approximation as describing a pure
quantum state $\rho_{s}\approx\ket{\psi_{0}}\bra{\psi_{0}}$. Furthermore,
$\ket{\psi_{0}}$ corresponds to the eigenstate of the Hamiltonian with the
smallest number of excitations. This is a consequence of the fact that the
steady state is simply the ground state of the full system, which is
unaffected by the decay terms in the master equation. The steady state could
thus be obtained from just the Hamiltonian without having to invoke the
Lindblad master equation.

To conclude, we have demonstrated a powerful method able to reproduce dynamics
of a quantum system interacting with an arbitrary EM environment in any
coupling regime by exploiting the mapping between nanophotonic spectral
densities and Lindblad-form master equations [34, 35]. By suppressing the
spectral density at negative frequencies, artificial effects inherent to
standard Lindblad forms can be removed. For simple model systems, we showed
that the accuracy of our approach considerably exceeds the one of state-of-the
art master equations. Furthermore, the approach can deal with realistic
nanophotonic systems where simple models are not available. It thus offers a
straightforward way for investigating extreme regimes of light-matter
interaction while taking into account the full mode spectrum of complex EM
environments.

###### Acknowledgements.

We acknowledge support by the European Research Council through Grant
ERC-2016-STG-714870, by the Spanish Ministry for Science, Innovation, and
Universities – Agencia Estatal de Investigación (AEI) through Grants
PID2021-125894NB-I00 and CEX2018-000805-M (through the María de Maeztu program
for Units of Excellence in Research and Development), by the Comunidad de
Madrid through Proyecto Sinérgico CAM 2020 Y2020/TCS-6545 (NanoQuCo-CM) and by
the “(MAD2D-CM)-UAM7” project funded by the Comunidad de Madrid, by the
Recovery, Transformation and Resilience Plan from Spain, and by
NextGenerationEU from the European Union.

## References

  * Garcia-Vidal  _et al._ [2021] F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, Manipulating Matter by Strong Coupling to Vacuum Fields, Science 373, eabd0336 (2021).
  * Kasprzak  _et al._ [2006] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and L. S. Dang, Bose-Einstein Condensation of Exciton Polaritons, Nature 443, 409 (2006).
  * Zhong  _et al._ [2017] X. Zhong, T. Chervy, L. Zhang, A. Thomas, J. George, C. Genet, J. A. Hutchison, and T. W. Ebbesen, Energy Transfer between Spatially Separated Entangled Molecules, Angew. Chem. Int. Ed. 56, 9034 (2017).
  * Chervy  _et al._ [2016] T. Chervy, J. Xu, Y. Duan, C. Wang, L. Mager, M. Frerejean, J. A. W. Münninghoff, P. Tinnemans, J. A. Hutchison, C. Genet, A. E. Rowan, T. Rasing, and T. W. Ebbesen, High-Efficiency Second-Harmonic Generation from Hybrid Light-Matter States, Nano Lett. 16, 7352 (2016).
  * Zasedatelev  _et al._ [2021] A. V. Zasedatelev, A. V. Baranikov, D. Sannikov, D. Urbonas, F. Scafirimuto, V. Y. Shishkov, E. S. Andrianov, Y. E. Lozovik, U. Scherf, T. Stöferle, R. F. Mahrt, and P. G. Lagoudakis, Single-Photon Nonlinearity at Room Temperature, Nature 597, 493 (2021).
  * Feist and Garcia-Vidal [2015] J. Feist and F. J. Garcia-Vidal, Extraordinary Exciton Conductance Induced by Strong Coupling, Phys. Rev. Lett. 114, 196402 (2015).
  * Galego  _et al._ [2016] J. Galego, F. J. Garcia-Vidal, and J. Feist, Suppressing Photochemical Reactions with Quantized Light Fields, Nat. Commun. 7, 13841 (2016).
  * Carmichael [1993] H. Carmichael,  _An Open Systems Approach to Quantum Optics_, Lecture Notes in Physics Monographs, Vol. 18 (Springer, Berlin, Heidelberg, 1993).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2007).
  * Manzano [2020] D. Manzano, A Short Introduction to the Lindblad Master Equation, AIP Advances 10, 025106 (2020).
  * Carmichael and Walls [1973] H. J. Carmichael and D. F. Walls, Master Equation for Strongly Interacting Systems, J. Phys. Math. Nucl. Gen. 6, 1552 (1973).
  * Askenazi  _et al._ [2017] B. Askenazi, A. Vasanelli, Y. Todorov, E. Sakat, J.-J. Greffet, G. Beaudoin, I. Sagnes, and C. Sirtori, Midinfrared Ultrastrong Light–Matter Coupling for THz Thermal Emission, ACS Photonics 4, 2550 (2017).
  * Yoshihara  _et al._ [2017a] F. Yoshihara, T. Fuse, S. Ashhab, K. Kakuyanagi, S. Saito, and K. Semba, Superconducting Qubit–Oscillator Circuit beyond the Ultrastrong-Coupling Regime, Nat. Phys. 13, 44 (2017a).
  * Yoshihara  _et al._ [2017b] F. Yoshihara, T. Fuse, S. Ashhab, K. Kakuyanagi, S. Saito, and K. Semba, Characteristic Spectra of Circuit Quantum Electrodynamics Systems from the Ultrastrong- to the Deep-Strong-Coupling Regime, Phys. Rev. A 95, 053824 (2017b).
  * Yoshihara  _et al._ [2018] F. Yoshihara, T. Fuse, Z. Ao, S. Ashhab, K. Kakuyanagi, S. Saito, T. Aoki, K. Koshino, and K. Semba, Inversion of Qubit Energy Levels in Qubit-Oscillator Circuits in the Deep-Strong-Coupling Regime, Phys. Rev. Lett. 120, 183601 (2018).
  * Bayer  _et al._ [2017] A. Bayer, M. Pozimski, S. Schambeck, D. Schuh, R. Huber, D. Bougeard, and C. Lange, Terahertz Light–Matter Interaction beyond Unity Coupling Strength, Nano Lett. 17, 6340 (2017).
  * Maissen  _et al._ [2014] C. Maissen, G. Scalari, F. Valmorra, M. Beck, J. Faist, S. Cibella, R. Leoni, C. Reichl, C. Charpentier, and W. Wegscheider, Ultrastrong Coupling in the near Field of Complementary Split-Ring Resonators, Phys. Rev. B 90, 205309 (2014).
  * Benz  _et al._ [2016] F. Benz, M. K. Schmidt, A. Dreismann, R. Chikkaraddy, Y. Zhang, A. Demetriadou, C. Carnegie, H. Ohadi, B. de Nijs, R. Esteban, J. Aizpurua, and J. J. Baumberg, Single-Molecule Optomechanics in “Picocavities”, Science 354, 726 (2016).
  * Barachati  _et al._ [2018] F. Barachati, J. Simon, Y. A. Getmanenko, S. Barlow, S. R. Marder, and S. Kéna-Cohen, Tunable Third-Harmonic Generation from Polaritons in the Ultrastrong Coupling Regime, ACS Photonics 5, 119 (2018).
  * Gambino  _et al._ [2014] S. Gambino, M. Mazzeo, A. Genco, O. Di Stefano, S. Savasta, S. Patanè, D. Ballarini, F. Mangione, G. Lerario, D. Sanvitto, and G. Gigli, Exploring Light–Matter Interaction Phenomena under Ultrastrong Coupling Regime, ACS Photonics 1, 1042 (2014).
  * Paravicini-Bagliani  _et al._ [2019] G. L. Paravicini-Bagliani, F. Appugliese, E. Richter, F. Valmorra, J. Keller, M. Beck, N. Bartolo, C. Rössler, T. Ihn, K. Ensslin, C. Ciuti, G. Scalari, and J. Faist, Magneto-Transport Controlled by Landau Polariton States, Nat. Phys. 15, 186 (2019).
  * Keller  _et al._ [2020] J. Keller, G. Scalari, F. Appugliese, S. Rajabali, M. Beck, J. Haase, C. A. Lehner, W. Wegscheider, M. Failla, M. Myronov, D. R. Leadley, J. Lloyd-Hughes, P. Nataf, and J. Faist, Landau Polaritons in Highly Nonparabolic Two-Dimensional Gases in the Ultrastrong Coupling Regime, Phys. Rev. B 101, 075301 (2020).
  * Kockum  _et al._ [2019] A. F. Kockum, A. Miranowicz, S. D. Liberato, S. Savasta, and F. Nori, Ultrastrong Coupling between Light and Matter, Nat. Rev. Phys. 1, 19 (2019).
  * Forn-Díaz  _et al._ [2019] P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Ultrastrong Coupling Regimes of Light-Matter Interaction, Rev. Mod. Phys. 91, 025005 (2019).
  * Ciuti and Carusotto [2006] C. Ciuti and I. Carusotto, Input-Output Theory of Cavities in the Ultrastrong Coupling Regime: The Case of Time-Independent Cavity Parameters, Phys. Rev. A 74, 033811 (2006).
  * Beaudoin  _et al._ [2011] F. Beaudoin, J. M. Gambetta, and A. Blais, Dissipation and Ultrastrong Coupling in Circuit QED, Phys. Rev. A 84, 043832 (2011).
  * Salmon  _et al._ [2022] W. Salmon, C. Gustin, A. Settineri, O. D. Stefano, D. Zueco, S. Savasta, F. Nori, and S. Hughes, Gauge-Independent Emission Spectra and Quantum Correlations in the Ultrastrong Coupling Regime of Open System Cavity-QED, Nanophotonics 11, 1573 (2022).
  * Ma and Law [2015] K. K. W. Ma and C. K. Law, Three-Photon Resonance and Adiabatic Passage in the Large-Detuning Rabi Model, Phys. Rev. A 92, 023842 (2015).
  * Bamba and Ogawa [2014] M. Bamba and T. Ogawa, Recipe for the Hamiltonian of System-Environment Coupling Applicable to the Ultrastrong-Light-Matter-Interaction Regime, Phys. Rev. A 89, 023817 (2014).
  * Settineri  _et al._ [2018] A. Settineri, V. Macrí, A. Ridolfo, O. Di Stefano, A. F. Kockum, F. Nori, and S. Savasta, Dissipation and Thermal Noise in Hybrid Quantum Systems in the Ultrastrong-Coupling Regime, Phys. Rev. A 98, 053834 (2018).
  * De Bernardis  _et al._ [2018] D. De Bernardis, P. Pilar, T. Jaako, S. De Liberato, and P. Rabl, Breakdown of Gauge Invariance in Ultrastrong-Coupling Cavity QED, Phys. Rev. A 98, 053819 (2018).
  * Di Stefano  _et al._ [2019] O. Di Stefano, A. Settineri, V. Macrì, L. Garziano, R. Stassi, S. Savasta, and F. Nori, Resolution of Gauge Ambiguities in Ultrastrong-Coupling Cavity Quantum Electrodynamics, Nat. Phys. 15, 803 (2019).
  * Sáez-Blázquez  _et al._ [2022] R. Sáez-Blázquez, D. de Bernardis, J. Feist, and P. Rabl, Can We Observe Non-Perturbative Vacuum Shifts in Cavity QED?, arxiv:2212.08675 (2022).
  * Medina  _et al._ [2021] I. Medina, F. J. García-Vidal, A. I. Fernández-Domínguez, and J. Feist, Few-Mode Field Quantization of Arbitrary Electromagnetic Spectral Densities, Phys. Rev. Lett. 126, 093601 (2021).
  * Sánchez-Barquilla  _et al._ [2022] M. Sánchez-Barquilla, F. J. García-Vidal, A. I. Fernández-Domínguez, and J. Feist, Few-Mode Field Quantization for Multiple Emitters, Nanophotonics 11, 4363 (2022).
  * Weiss [2012] U. Weiss,  _Quantum Dissipative Systems_, 4th ed. (World Scientific, 2012).
  * de Vega and Alonso [2017] I. de Vega and D. Alonso, Dynamics of Non-Markovian Open Quantum Systems, Rev. Mod. Phys. 89, 015001 (2017).
  * Buhmann and Welsch [2008] S. Y. Buhmann and D.-G. Welsch, Casimir-Polder Forces on Excited Atoms in the Strong Atom-Field Coupling Regime, Phys. Rev. A 77, 012110 (2008).
  * Feist  _et al._ [2020] J. Feist, A. I. Fernández-Domínguez, and F. J. García-Vidal, Macroscopic QED for Quantum Nanophotonics: Emitter-Centered Modes as a Minimal Basis for Multiemitter Problems, Nanophotonics 10, 477 (2020).
  * Rabi [1937] I. I. Rabi, Space Quantization in a Gyrating Magnetic Field, Phys. Rev. 51, 652 (1937).
  * Braak [2011] D. Braak, Integrability of the Rabi Model, Phys. Rev. Lett. 107, 100401 (2011).
  * Imamoğlu [1994] A. Imamoğlu, Stochastic Wave-Function Approach to Non-Markovian Systems, Phys. Rev. A 50, 3650 (1994).
  * Tamascelli  _et al._ [2018] D. Tamascelli, A. Smirne, S. F. Huelga, and M. B. Plenio, Nonperturbative Treatment of Non-Markovian Dynamics of Open Quantum Systems, Phys. Rev. Lett. 120, 030402 (2018).
  * Reid and Johnson [2015] M. T. H. Reid and S. G. Johnson, Efficient Computation of Power, Force, and Torque in BEM Scattering Calculations, IEEE Trans. Antennas Propag. 63, 3588 (2015).
  * Delga  _et al._ [2014] A. Delga, J. Feist, J. Bravo-Abad, and F. J. Garcia-Vidal, Quantum Emitters Near a Metal Nanoparticle: Strong Coupling and Quenching, Phys. Rev. Lett. 112, 253601 (2014).
  * Li  _et al._ [2016] R.-Q. Li, D. Hernángomez-Pérez, F. J. García-Vidal, and A. I. Fernández-Domínguez, Transformation Optics Approach to Plasmon-Exciton Strong Coupling in Nanocavities, Phys. Rev. Lett. 117, 107401 (2016).