# Frequency-dependent Sternheimer linear-response formalism for strongly
coupled light-matter systems

Davis M. Welakuh<EMAIL_ADDRESS>Max Planck Institute for the
Structure and Dynamics of Matter and Center for Free-Electron Laser Science &
Department of Physics, Luruper Chaussee 149, 22761 Hamburg, Germany  Harvard
John A. Paulson School Of Engineering And Applied Sciences, Harvard
University, Cambridge, Massachusetts 02138, USA  Johannes Flick
<EMAIL_ADDRESS>Center for Computational Quantum Physics,
Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA  Michael
Ruggenthaler<EMAIL_ADDRESS>Max Planck Institute for the
Structure and Dynamics of Matter and Center for Free-Electron Laser Science &
Department of Physics, Luruper Chaussee 149, 22761 Hamburg, Germany  The
Hamburg Center for Ultrafast Imaging, Luruper Chaussee 149, 22761 Hamburg,
Germany.  Heiko Appel<EMAIL_ADDRESS>Max Planck Institute for the
Structure and Dynamics of Matter and Center for Free-Electron Laser Science &
Department of Physics, Luruper Chaussee 149, 22761 Hamburg, Germany  Angel
Rubio<EMAIL_ADDRESS>Max Planck Institute for the Structure and
Dynamics of Matter and Center for Free-Electron Laser Science & Department of
Physics, Luruper Chaussee 149, 22761 Hamburg, Germany  Center for
Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York,
NY 10010, USA  The Hamburg Center for Ultrafast Imaging, Luruper Chaussee 149,
22761 Hamburg, Germany.

###### Abstract

The rapid progress in quantum-optical experiments especially in the field of
cavity quantum electrodynamics and nanoplasmonics, allows to substantially
modify and control chemical and physical properties of atoms, molecules and
solids by strongly coupling to the quantized field. Alongside such
experimental advances has been the recent development of ab-initio approaches
such as quantum electrodynamical density-functional theory (QEDFT) that is
capable of describing these strongly coupled systems from first-principles. To
investigate response properties of relatively large systems coupled to a wide
range of photon modes, ab-initio methods that scale well with system size
become relevant. In light of this, we extend the linear-response Sternheimer
approach within the framework of QEDFT to efficiently compute excited-state
properties of strongly coupled light-matter systems. Using this method, we
capture features of strong light-matter coupling both in the dispersion and
absorption properties of a molecular system strongly coupled to the modes of a
cavity. We exemplify the efficiency of the Sternheimer approach by coupling
the matter system to the continuum of an electromagnetic field. We observe
changes in the spectral features of the coupled system as Lorentzian line
shapes turn into Fano resonances when the molecule interacts strongly with the
continuum of modes. This work provides an alternative approach for computing
efficiently excited-state properties of large molecular systems interacting
with the quantized electromagnetic field.

##  I Introduction

Strong interactions between light and matter within enhanced photonic
environments such as optical cavities and plasmonic devices have attracted
great attention in recent years. The signature of such strong interactions is
the formation of new hybrid light-matter states (polaritons) which are
manifested by a Rabi splitting in the spectrum of the coupled system. These
new states of matter can be used to control, for instance, chemical reactions
Hutchison  _et al._ (2012); Thomas  _et al._ (2016); Schäfer  _et al._ (2021),
enhance charge and energy transport Coles  _et al._ (2014); Zhong  _et al._
(2016); Feist and Garcia-Vidal (2015), selectively manipulate electronic
excited-states Stranius  _et al._ (2018), to name a few examples. Such coupled
light-matter systems have the tendency to exhibit significantly different
properties than the uncoupled subsystems even at ambient conditions, which
suggests various interesting applications in chemistry and material science
Chikkaraddy  _et al._ (2016); Bellessa  _et al._ (2004); Ebbesen (2016);
Santhosh  _et al._ (2016). These intriguing effects caused by the emergence of
polaritons manifest strongly in the excited-state properties of the coupled
systems, for example, in the absorption or emission spectra Chikkaraddy  _et
al._ (2016); Bellessa  _et al._ (2004); Ebbesen (2016); Stranius  _et al._
(2018).

The occurrence of different effects due to the emergence of polaritons shows
the complexity that arises when light and matter strongly mix. Due to this
inherent complexity of the coupled light-matter system, the theoretical
description of these effects are non-trivial. Quite often, the coupled system
is studied with quantum optical models that potentially over-simplify the
matter subsystem. One of such simplifications selects just a few energy-levels
of an atomic or molecular system and couples it to the photon modes of an
optical cavity Dicke (1954); Jaynes and Cummings (1963); Garraway (2011).
Another common simplification is on the photon side where a realistic cavity
which is normally open is reduced to just a few modes that cannot account for
the finite lifetime of excitations. However, in many cases these
phenomenological models are not sufficient to capture important details of the
coupled system, for instance, the emergence of bound states in the continuum
Sidler  _et al._ (2020), how collective strong coupling leads to local
modification of chemical properties Sidler  _et al._ (2021a) and in cavity-
modified chemistry where the reaction-rate is reduced under cavity induced
resonant vibrational strong coupling Thomas  _et al._ (2016); Schäfer  _et
al._ (2020). This calls for ab initio methods which allow to treat from first-
principles complex matter systems interacting with the quantized field
Ruggenthaler  _et al._ (2014); Flick  _et al._ (2015); Haugland  _et al._
(2020) within non-relativistic quantum electrodynamics (QED) Spohn (2004);
Ruggenthaler  _et al._ (2018); Flick  _et al._ (2018a). Amongst the existing
first-principles methods for treating strongly coupled light-matter systems,
quantum electrodynamical density-functional theory (QEDFT) has become a
valuable approach for describing ground- and excited-state properties of
complex matter systems coupled to a photonic environment Flick  _et al._
(2018b, 2019). The Casida-like approach Casida (1996) common to molecular and
quantum chemistry was recently extended to a matter-photon description within
the linear-response framework of QEDFT Flick  _et al._ (2019). The feasibility
of treating the excited-state properties of a single molecule and an ensemble
of molecules coupled to a realistic description of a cavity has been
demonstrated Flick and Narang (2020); Flick  _et al._ (2019); Wang  _et al._
(2021); Sidler  _et al._ (2021a). A different approach within QEDFT is to
solve the time-dependent Kohn-Sham equations coupled to the Maxwell equations
in real time Flick and Narang (2018); Welakuh (2021); Jestädt  _et al._
(2020); Schäfer  _et al._ (2021). One major advantage the real-time approach
has is that it scales favorably with the system size as it involves only
occupied Kohn-Sham orbitals but to obtain a converged response spectrum
requires a long time-propagation, which is not favorable for larger systems.
On the other hand, the Casida approach requires both occupied and unoccupied
orbitals and it also scales with the number of photon modes considered.

In addition to these methods, there is another successful scheme that can
combine the strengths of the previously mentioned methods known as the
Sternheimer approach Sternheimer (1951). The Sternheimer approach has been
used for a long time in the context of density-functional perturbation theory,
for instance, for calculating phonon spectra Baroni  _et al._ (1987). Recent
applications of the Sternheimer equation have also been used to compute the
frequency-dependent electronic response Hübener and Giustino (2014a, b);
Giustino  _et al._ (2010); Andrade  _et al._ (2007). The Sternheimer equation
has been formulated within the framework of time-dependent density-functional
theory (TDDFT) that allows to study dynamic response of much larger complex
systems as it includes only occupied orbitals Andrade  _et al._ (2007);
Hofmann  _et al._ (2018). One of a few advantages the Sternheimer approach has
over real-time TDDFT is that it is formulated in the frequency space and the
responses at different frequencies can be computed independently of each other
allowing for the use of parallelization schemes that speed up the computation.
Another advantage is that since the responses at different frequencies can be
treated independently, we can compute any part of the spectrum without
necessarily starting from the zero-frequency. In this work, we extend the
frequency-dependent Sternheimer approach of TDDFT to the framework of QEDFT.
An advantage the electron-photon Sternheimer approach has over the Casida
approach is that it scales well with the system size since only occupied
orbitals are treated explicitly and the arbitrarily many but finite photon
modes that can be included does not add to this scaling. This approach becomes
useful for investigations in polaritonic chemistry or materials in nano-
plasmonic cavities that usually consider a large number of particles
interacting with the electromagnetic field. We start by showing the
applicability of the method in capturing not only the absorption properties of
strongly coupled light-matter system but also the modified dispersion
properties of the coupled system for the case of an azulene molecule. In
addition, we show the spectra of the photon field which capture similar
features of strong light-matter coupled systems indicating how the hybrid
characteristics can be viewed from either of the subsystems at the same time
highlighting the cross-talk between the subsystems. To show the efficiency of
the Sternheimer approach, we study the absorption spectrum of a lithium
hydride (LiH) molecule coupled to a continuum of photon modes. For the case of
coupling the molecule weakly to half a million photon modes, we recover the
spectrum of the free space case. By effectively enhancing the coupling of the
bath modes to the molecule, we observe changes in the spectrum as the
Lorentzian line shape turn into Fano resonances.

This article is structured as follows. First, we present the physical setting
of a many-electron system coupled to photons in non-relativistic QED and
subsequently present the linear-response setting of this framework. Secondly,
we present in Sec. IV a derivation of the frequency-dependent Sternheimer
approach for electron-photon coupled systems in the linear-response regime
formulated within the framework of QEDFT and discuss numerical details of the
Sternheimer scheme. In the next section, we investigate the complex
polarizability of a molecular system coupled to a high-Q optical cavity and
highlight how the absorption and dispersion properties get modified due to
strong light-matter coupling. Also, we show for the same molecular system the
polaritonic features that arise in the spectra of the photon field.
Furthermore, we demonstrate the efficient and low computational cost of the
Sternheimer approach by coupling a LiH molecule to a (discretized) continuum
of states of photon modes and show the physical effects that arises. Finally
we present a conclusion and an outlook.

##  II From microscopic fields to the quantum description of light-matter
interaction

We are interested in the dynamics of matter interacting with the
electromagnetic field within the setting of non-relativistic QED where both
constituents of the coupled system are treated on an equal quantized footing.
While this setting of slowly moving charged particles can be deduced from QED,
concepts from classical electrodynamics are equally instructive to arrive at
this description. In this regard, we lay emphasizes on the full description of
the electromagnetic field that couples to a matter system. Our starting point
is the inhomogeneous microscopic Maxwell equations for the transverse part of
the electromagnetic field Loudon (2000):

$\displaystyle\boldsymbol{\nabla}\times\textbf{E}(\textbf{r},t)$
$\displaystyle=-\frac{\partial}{\partial t}\textbf{B}(\textbf{r},t)\,,$
$\displaystyle\boldsymbol{\nabla}\times\textbf{B}(\textbf{r},t)$
$\displaystyle=\frac{1}{c^{2}}\left[\frac{\partial}{\partial
t}\textbf{E}(\textbf{r},t)+\mu_{0}c^{2}\textbf{j}(\textbf{r},t)\right]\,,$

where $\textbf{E}(\textbf{r},t)$ and $\textbf{B}(\textbf{r},t)$ are the
classical electric and magnetic fields, respectively. The transverse charge
current $\textbf{j}(\textbf{r},t)$ represent both the free and bound current.
If we consider $\textbf{j}(\textbf{r},t)$ to represent only the bound current,
then it is related to the polarization $\textbf{P}(\textbf{r},t)$ of the
matter as $\textbf{j}(\textbf{r},t)=\frac{\partial}{\partial
t}\textbf{P}(\textbf{r},t)$. The Maxwell’s equations takes into account the
back-reaction of the matter on the electromagnetic field. For a quantum
mechanical description, the field variables are promoted to field operators in
the Heisenberg picture. In this representation, the energy of the transverse
electromagnetic field can be expressed as Abedi  _et al._ (2018):

$\displaystyle\hat{H}_{\text{EM}}$ $\displaystyle=\frac{\epsilon_{0}}{2}\int
d^{3}\textbf{r}\left[\hat{\textbf{E}}^{2}(\textbf{r})+c^{2}\hat{\textbf{B}}^{2}(\textbf{r})\right]$
$\displaystyle=\frac{\epsilon_{0}}{2}\int
d^{3}\textbf{r}\left[\frac{1}{\epsilon_{0}^{2}}\left(\hat{\textbf{D}}(\textbf{r})-\hat{\textbf{P}}(\textbf{r})\right)^{2}+c^{2}\hat{\textbf{B}}^{2}(\textbf{r})\right],$
(1)

where we introduced the displacement field
$\hat{\textbf{D}}(\textbf{r})=\epsilon_{0}\hat{\textbf{E}}(\textbf{r})+\hat{\textbf{P}}(\textbf{r})$.
Equation (1) will differ from other forms only in the choice of canonical
variables and here we impose a commutation relation between $\hat{\textbf{B}}$
and $\hat{\textbf{D}}$ to be
$\left[\hat{\textbf{B}}^{i}(\textbf{r}),\hat{\textbf{D}}^{j}(\textbf{r}^{\prime})\right]=-i\hbar\varepsilon^{ijk}{\partial_{k}}\delta(\textbf{r}-\textbf{r}^{\prime})$.
For any photonic environment of varying geometry, the fields in Eq. (1) can be
expanded in the modes Abedi  _et al._ (2018)

$\displaystyle\hat{\textbf{D}}(\textbf{r})$
$\displaystyle=\sum_{\alpha}\textbf{S}_{\alpha}(\textbf{r})\hat{d}_{\alpha}\,,$
(2) $\displaystyle\hat{\textbf{B}}(\textbf{r})$
$\displaystyle=\sum_{\alpha}\frac{1}{\omega_{\alpha}}\hat{b}_{\alpha}\boldsymbol{\nabla}\times\textbf{S}_{\alpha}(\textbf{r})\,,$
(3) $\displaystyle\hat{\textbf{P}}(\textbf{r})$
$\displaystyle=\sum_{\alpha}\textbf{S}_{\alpha}(\textbf{r})\left(\textbf{S}_{\alpha}(\textbf{r}_{0})\cdot\hat{\textbf{R}}\right)\,.$
(4)

The expansion coefficients $\hat{d}_{\alpha}$ and $\hat{b}_{\alpha}$ are
respectively the field amplitudes of the electric displacement and the
magnetic field, $\textbf{S}_{\alpha}(\textbf{r})$ are the mode functions and
$\hat{\textbf{R}}=\sum_{i}e\,\hat{\textbf{r}}_{i}$ is the total electronic
dipole operator. In Eq. (4), we employed the dipole approximation when the
electromagnetic field interacts with the matter system via the electronic
dipole. Making a substitution of Eqs. (2)-(4) into (1) results to the
following expression of the electromagnetic energy Abedi  _et al._ (2018):

$\displaystyle\hat{H}_{\text{EM}}$
$\displaystyle=\sum_{\alpha}\frac{1}{2}\left[\hat{p}^{2}_{\alpha}+\omega^{2}_{\alpha}\left(\hat{q}_{\alpha}\\!-\\!\frac{\boldsymbol{\lambda}_{\alpha}}{\omega_{\alpha}}\cdot\hat{\textbf{R}}\right)^{2}\right].$
(5)

The displacement coordinate $\hat{q}_{\alpha}$ and momentum operator
$\hat{p}_{\alpha}$ are related to the amplitudes as
$\hat{d}_{\alpha}=\sqrt{\epsilon_{0}}\omega_{\alpha}\hat{q}_{\alpha}$ and
$\hat{b}_{\alpha}=\sqrt{1/\epsilon_{0}}\hat{p}_{\alpha}$ where they satisfy
the commutation relation
$\left[\hat{q}_{\alpha},\hat{p}_{\alpha^{\prime}}\right]=i\hbar\,\delta_{\alpha,\alpha^{\prime}}$.
Equation (5) tells us that what would normally be the purely photonic
Hamiltonian is now a mixture of matter and photon degrees. The term
$\boldsymbol{\lambda}_{\alpha}$ in Eq. (5) is the light-matter coupling
strength given as

$\displaystyle\boldsymbol{\lambda}_{\alpha}=\frac{1}{\sqrt{\epsilon_{0}}}\,\textbf{S}_{\alpha}(\textbf{r}_{0})\,,$
(6)

where the mode function is evaluated at the center of charge Milonni (1993).
In deriving Eq. (5) we have assumed a finite photonic environment with
appropriate boundary conditions. For real systems we, however, have usually a
continuum of modes, i.e., the cavity geometry is open to free infinite space.
We can approximate this situation by extending the quantization volume of the
electromagnetic field beyond the photonic environment and thus work with a
discretized continuum. By making this discretization finer and finer, i.e.
take the quantization volume to infinity, we can approximate the open-cavity
situation arbitrarily well. The discrete continuum description of the photon
field has the advantage that it accounts for the emission or absorption of a
photon in real space Flick  _et al._ (2017) and allows for modeling an open
photonic environment Welakuh  _et al._ (2021). Together with the Hamiltonian
representing the bound charged particles, i.e. the kinetic energy, binding and
interaction potentials, equation (5) constitutes the so-called Pauli-Fierz
Hamiltonian in the length gauge Rokaj  _et al._ (2018); Schäfer  _et al._
(2020). In the case where we include time-dependent external perturbations,
the length gauge Hamiltonian is given by

$\displaystyle\hat{H}(t)$
$\displaystyle=\sum\limits_{i=1}^{N}\left(\frac{\hat{\textbf{p}}_{i}^{2}}{2m}+v_{\textrm{ext}}(\hat{\textbf{r}}_{i},t)\right)+\sum\limits_{i>j}^{N}w(|\hat{\textbf{r}}_{i}-\hat{\textbf{r}}_{j}|)$
$\displaystyle+\sum_{\alpha=1}^{M}\frac{1}{2}\left[\hat{p}^{2}_{\alpha}+\omega^{2}_{\alpha}\left(\hat{q}_{\alpha}\\!-\\!\frac{\boldsymbol{\lambda}_{\alpha}}{\omega_{\alpha}}\cdot\hat{\textbf{R}}\right)^{2}\right]\\!+\\!\sum_{\alpha=1}^{M}\frac{j_{\textrm{ext}}^{(\alpha)}(t)}{\omega_{\alpha}}\hat{q}_{\alpha}.$
(7)

Here the $N$ electrons are described by the electronic coordinates
$\hat{\textbf{r}}_{i}$ and the momentum operator $\hat{\textbf{p}}_{i}$, which
satisfy the commutation relation
$\left[\hat{\textbf{r}}_{i},\hat{\textbf{p}}_{j}\right]=i\hbar\,\delta_{ij}$.
The interaction due to the longitudinal part of the photon field
$w(|\hat{\textbf{r}}_{i}-\hat{\textbf{r}}_{j}|)$ can be written as a mode-
expansion in Coulomb gauge which for the free-space case results to the
standard Coulomb interaction

$\displaystyle
w(|\hat{\textbf{r}}-\hat{\textbf{r}}^{\prime}|)=\sum_{n}\frac{e^{2}}{\textbf{k}_{n}}\frac{e^{i\textbf{k}_{n}(\hat{\textbf{r}}-\hat{\textbf{r}}^{\prime})}}{\epsilon_{0}L^{3}}\xrightarrow[]{L\rightarrow\infty}\frac{e^{2}}{4\pi\epsilon_{0}}\frac{1}{\left|\hat{\textbf{r}}-\hat{\textbf{r}}^{\prime}\right|}\,,$

where $\textbf{k}_{n}=2\pi\textbf{n}/L$ are the allowed wave vectors of the
photon field for an arbitrarily large but finite box of length $L$ Greiner and
Reinhardt (1996). For the transverse field, we consider an arbitrarily large
but finite number of photon modes $M$. It is important to note that when we
sample a large number of modes to describe the photon continuum, we might need
to use the bare mass of the electrons instead of the renormalized physical
mass Spohn (2004); Rokaj  _et al._ (2022). In Sec. V.2 of this work, we make
the common assumption that only the sampled continuum due to a cavity or
photonic nanostructure is changed with respect to the free space case. The
rest of the continuum of modes not affected by the cavity is subsumed in the
already renormalized physical mass of the electrons. The time-dependent
external potential and current in Eq. (7) can be split into

$\displaystyle v_{\textrm{ext}}(\textbf{r},t)$
$\displaystyle=v(\textbf{r})+\delta v(\textbf{r},t),$ $\displaystyle
j_{\textrm{ext}}^{(\alpha)}(t)$ $\displaystyle=j_{\alpha}+\delta
j_{\alpha}(t),$

where $v(\textbf{r})$ describes the attractive potentials of the nuclei and
$\delta v(\textbf{r},t)$ a classical external probe field that couples to the
electronic subsystem. Also, the static part $j_{\alpha}$ merely polarizes the
vacuum and the time-dependent part $\delta j_{\alpha}(t)$ then generates
photons in the mode $\alpha$. Either of these perturbations can be used to
probe the coupled light-matter system.

##  III Linear Response formulation in the length gauge

To characterize the properties of a system, one can investigate its response
to an external perturbation. In the case of a weak external perturbation, we
have access to linear response properties of the system such as its
polarizability which gives access to its excitation energies and oscillator
strengths. Usually this requires knowledge of the linear density response and
in the case of a coupled matter-photon system, we have access to the
displacement field Flick  _et al._ (2019). In the length gauge, the linear
response of the electron density $n(\textbf{r},t)$ to the external potential
$\delta v(\textbf{r},t)$ and charge current $\delta j_{\alpha}(t)$ yields the
response equation Flick  _et al._ (2019)

$\displaystyle\delta n(\textbf{r},t)$ $\displaystyle=\int dt^{\prime}\int
d^{3}\textbf{r}^{\prime}\chi_{n}^{n}(\textbf{r},t;\textbf{r}^{\prime},t^{\prime})\delta
v(\textbf{r}^{\prime},t^{\prime})$ $\displaystyle\quad+\sum_{\alpha=1}^{M}\int
dt^{\prime}\chi^{n}_{q_{\alpha}}(\textbf{r},t;t^{\prime})\delta
j_{\alpha}(t^{\prime}).$ (8)

Due to the coupling between light and matter, we can equally compute the
linear response of the photon displacement coordinate $q_{\alpha}(t)$ due to
the external potential $\delta v(\textbf{r},t)$ and current $\delta
j_{\alpha}(t)$ that results to the response equation Flick  _et al._ (2019)

$\displaystyle\delta q_{\alpha}(t)$ $\displaystyle=\int dt^{\prime}\int
d^{3}\textbf{r}^{\prime}\chi^{q_{\alpha}}_{n}(t;\textbf{r}^{\prime},t^{\prime})\delta
v(\textbf{r}^{\prime},t^{\prime})$
$\displaystyle\quad+\sum_{\alpha^{\prime}=1}^{M}\int
dt^{\prime}\chi^{q_{\alpha}}_{q_{\alpha^{\prime}}}(t,t^{\prime})\delta
j_{\alpha^{\prime}}(t^{\prime}).$ (9)

The response functions $\chi_{n}^{n}$, $\chi^{n}_{q_{\alpha}}$,
$\chi^{q_{\alpha}}_{n}$ and $\chi^{q_{\alpha}}_{q_{\alpha^{\prime}}}$ are
intrinsic properties of the electron-photon coupled system which can be
computed to obtain excited-state properties of the system. However, computing
these response equations or the response functions directly is usually very
challenging even for the electron-only system. One possible way to do this
efficiently is to reformulate the response equations using the Maxwell-Kohn-
Sham system of QEDFT that reproduces the same response of the density and
photon coordinate Flick  _et al._ (2015); Ruggenthaler  _et al._ (2014); Flick
_et al._ (2019). In such a setting, the response functions can be computed
approximately giving access to, for instance, excitation energies and
oscillator strengths. We recently extended the Casida equation Casida (1996)
within the framework of QEDFT to treat electron-photon coupled systems Flick
_et al._ (2019). This approach computes the excitation energies and oscillator
strengths of either of the coupled response functions
$\left\\{\chi_{n}^{n},\chi^{q_{\alpha}}_{n}\right\\}$ and
$\left\\{\chi^{n}_{q_{\alpha}},\chi^{q_{\alpha}}_{q_{\alpha^{\prime}}}\right\\}$
by diagonalizing a pseudo-eigenvalue equation Flick  _et al._ (2019). The
Casida approach which requires both occupied and unoccupied Kohn-Sham orbitals
in addition with the sampled photon modes scales favorably for small coupled
systems Sidler  _et al._ (2021a). However, for larger electronic systems
coupled to many photon modes, this leads to a rapid increase in computational
effort in the Casida approach as the Casida matrix equation increases in size.

An alternative approach which rather computes the response equations instead
of the response function is the frequency-dependent Sternheimer equation
Sternheimer (1951). Formulated within the framework of TDDFT, this method
computes the linear density response due to an external weak perturbation
Andrade  _et al._ (2007); Hofmann  _et al._ (2018); Ullrich (2011) as well as
non-linear responses Marquesa  _et al._ (2012); Andrade  _et al._ (2007). This
approach has several advantages the main one being that it relies only on the
occupied Kohn-Sham orbitals thereby restricting the computation complexity for
very large systems. In the following we extend this approach to treat
electron-photon coupled system within the framework of QEDFT.

##  IV The Sternheimer approach for electron-photon coupled systems

Practical ab-initio methods for computing optical excitation spectra are
usually achieved by applying many-body methods that solve the correlated
problem in an approximate way. One of such many-body methods is TDDFT which is
considered a very promising methodology since it provides a good balance
between accuracy and computational cost. Within the context of TDDFT there
exist different formalism for computing optical excitation spectra Marquesa
_et al._ (2012). The frequency-dependent Sternheimer method formulated within
TDDFT is a perturbative approach on the Kohn-Sham orbitals that computes the
density response without relying explicitly on unoccupied Kohn-Sham orbitals
Andrade  _et al._ (2007); Hofmann  _et al._ (2018). Based on this advantage,
an extension of this approach to the setting of QEDFT to treat complex atomic
and molecular systems coupled to arbitrary large but finite number of photon
modes is an important alternative method to existing QEDFT methods Flick  _et
al._ (2019); Flick and Narang (2018). The derivation presented here is solely
in the frequency space following that of Ref. Ullrich (2011). In an electron-
only description, the Sternheimer approach obtains only electronic observables
such as the electron density response. However, when this method is formulated
within the QEDFT framework we have in addition to the density response, the
response of the photon displacement coordinate (field). The mode-resolved
response of the field gives access to physical processes such as the
absorption or emission process. Starting with the reformulation of the density
and photon displacement coordinate responses in the QEDFT framework, the
coupled response due to a weak external potential $\delta
v(\textbf{r},\omega)$ are Flick  _et al._ (2019):

$\displaystyle\delta n(\textbf{r},\omega)$ $\displaystyle=\int
d^{3}\textbf{r}^{\prime}\chi_{n,s}^{n}(\textbf{r},\textbf{r}^{\prime},\omega)\delta
v_{\text{KS}}(\textbf{r}^{\prime},\omega)\,,$ (10) $\displaystyle\delta
q_{\alpha}(\omega)$
$\displaystyle=\chi_{q_{\alpha,s}}^{q_{\alpha}}(\omega)\delta
j_{\alpha,\text{KS}}(\omega)\,,$ (11)

where the first-order Kohn-Sham potential and currents in Eqs. (10) and (11)
are given in terms of the interacting density and photon coordinate responses
including mean-field exchange-correlation kernels:

$\displaystyle\delta v_{\text{KS}}(\textbf{r}^{\prime},\omega)$
$\displaystyle=\delta v(\textbf{r}^{\prime},\omega)+\int
d^{3}\textbf{y}f_{\text{Mxc}}^{n}(\textbf{r}^{\prime},\textbf{y},\omega)\delta
n(\textbf{y},\omega)$
$\displaystyle\quad+\sum_{\alpha}f_{\text{Mxc}}^{q_{\alpha}}(\textbf{r}^{\prime},\omega)\delta
q_{\alpha}(\omega)\,,$ (12) $\displaystyle\delta j_{\alpha,\text{KS}}(\omega)$
$\displaystyle=\int
d^{3}\textbf{y}\,g_{\text{M}}^{n_{\alpha}}(\omega,\textbf{y})\delta
n(\textbf{y},\omega)\,.$ (13)

The mean-field exchange-correlation kernels $f_{\text{Mxc}}^{n}$ and
$f_{\text{Mxc}}^{q_{\alpha}}$ are defined to be the variation of the mean-
field exchange-correlation potential with respect to the density and photon
coordinate, respectively while the mean-field kernel
$g_{\text{M}}^{n_{\alpha}}$ is the variation of the mean-field current with
respect to the electron density Flick  _et al._ (2019). These kernels account
for the correlations in the Kohn-Sham setting of QEDFT in linear response. The
non-interacting response functions of the decoupled electronic and photonic
subsystems of Eq. (10) and (11) are given explicitly as Flick  _et al._
(2019):

$\displaystyle\chi_{n,s}^{n}(\textbf{r},\textbf{r}^{\prime},\omega)$
$\displaystyle=\sum_{k=1}^{N_{v}}\sum_{j=1}^{\infty}\left[\frac{\varphi_{j}(\textbf{r})\varphi_{k}(\textbf{r}^{\prime})\varphi_{k}^{*}(\textbf{r})\varphi_{j}^{*}(\textbf{r}^{\prime})}{\omega-(\epsilon_{j}-\epsilon_{k})+i\eta}\right.$
$\displaystyle\qquad\qquad\quad\left.-\frac{\varphi_{k}(\textbf{r})\varphi_{j}(\textbf{r}^{\prime})\varphi_{j}^{*}(\textbf{r})\varphi_{k}^{*}(\textbf{r}^{\prime})}{\omega+(\epsilon_{j}-\epsilon_{k})+i\eta}\right],$
(14) $\displaystyle\chi_{q_{\alpha,s}}^{q_{\alpha}}(\omega)$
$\displaystyle=\frac{1}{2\omega_{\alpha}^{2}}\left(\frac{1}{\omega-\omega_{\alpha}+i\eta^{\prime}}-\frac{1}{\omega+\omega_{\alpha}+i\eta^{\prime}}\right).$
(15)

Here, $\epsilon_{k}$ and $\varphi_{k}(\textbf{r})$ are the ground-state
energies and orbitals of the Kohn-Sham system and $\omega_{\alpha}$ is the
frequency of the $\alpha$ mode. The parameters $\eta$ and $\eta^{\prime}$
shifts the poles (excitation energies) of Eqs. (14) and (15) to the lower half
of the complex plane and are in general not equal in both uncoupled systems.

Since the Sternheimer method is a perturbative approach on the Kohn-Sham
orbitals, we start by describing the unperturbed equilibrium setting of the
coupled electron-photon system as this corresponds to the zeroth-order of a
perturbation expansion, for example, of the density. For this case, we start
by describing the static Kohn-Sham system of ground-state QEDFT Ruggenthaler
(2015) where we have to solve the coupled Kohn-Sham equations

$\displaystyle\hat{h}\,\varphi_{k}(\textbf{r})$
$\displaystyle=\left[\frac{\hat{\textbf{p}}^{2}}{2m}+v(\textbf{r})+v_{\textrm{Mxc}}\left(\left[n,q_{\alpha}\right];\textbf{r}\right)\right]\varphi_{k}(\textbf{r})$
$\displaystyle=\epsilon_{k}\,\varphi_{k}(\textbf{r})\,,$ (16)
$\displaystyle\omega_{\alpha}^{2}q_{\alpha}$
$\displaystyle=-\frac{1}{\omega_{\alpha}}\left(j_{\alpha}-\omega_{\alpha}^{2}\int
d^{3}\textbf{r}\,\boldsymbol{\lambda}_{\alpha}\cdot\textbf{r}\,n(\textbf{r})\right),$
(17)

where $\hat{h}=\hat{h}\left(\left[v,n,q_{\alpha}\right]\right)$ is the ground-
state Kohn-Sham Hamiltonian. The ground-state density can be obtained from the
Kohn-Sham orbitals as $n(\textbf{r})=\sum_{k}|\varphi_{k}(\textbf{r})|^{2}$
and the photon coordinate from Eq. (17). The mean-field exchange-correlation
potential $v_{\textrm{Mxc}}(\textbf{r})$ represents the longitudinal
interactions between the electrons as well as all the transversal interactions
of the electrons with the photon field.

To solve for the linear density response and photon displacement coordinate of
Eq. (10) and (11), we first start by substituting Eq. (14) into the density
response $n(\textbf{r},\omega)$ of Eq. (10). The density response can now be
written in a form that includes a sum over only occupied orbitals as

$\displaystyle\delta n(\textbf{r},\omega)$
$\displaystyle=\sum_{k=1}^{N_{v}}\left[\varphi_{k}^{*}(\textbf{r})\varphi_{k}^{(+)}(\textbf{r},\omega)+\varphi_{k}(\textbf{r})\left[\varphi_{k}^{(-)}(\textbf{r},\omega)\right]^{*}\right],$
(18)

where the first-order response of the Kohn-Sham orbitals
$\varphi_{k}^{(\pm)}(\textbf{r},\omega)$ in Eq. (18) are given by

$\displaystyle\varphi_{k}^{(+)}(\textbf{r},\omega)$ $\displaystyle=\int
d^{3}\textbf{r}^{\prime}\sum_{j=1}^{\infty}\frac{\varphi_{j}(\textbf{r})\varphi_{j}^{*}(\textbf{r}^{\prime})\varphi_{k}(\textbf{r}^{\prime})}{\omega-(\epsilon_{j}-\epsilon_{k})+i\eta}\delta
v_{\textrm{KS}}(\textbf{r}^{\prime},\omega),$ (19)
$\displaystyle\varphi_{k}^{(-)}(\textbf{r},\omega)$ $\displaystyle=-\int
d^{3}\textbf{r}^{\prime}\sum_{j=1}^{\infty}\frac{\varphi_{j}(\textbf{r})\varphi_{j}^{*}(\textbf{r}^{\prime})\varphi_{k}^{*}(\textbf{r}^{\prime})}{\omega+(\epsilon_{j}-\epsilon_{k})+i\eta}\delta
v_{\textrm{KS}}(\textbf{r}^{\prime},\omega).$ (20)

Here, solving for the Kohn-Sham orbital responses is highly involved since we
need to first determine infinitely many Kohn-Sham orbitals and evaluate an
infinite sum over all these orbitals. However, this can be circumvented by
acting with $\left(\omega-\hat{h}+\epsilon_{k}+i\eta\right)$ and
$\left(\omega+\hat{h}-\epsilon_{k}+i\eta\right)$ on Eqs. (19) and (20),
respectively and using the static Kohn-Sham equation of Eq. (16) in concert
with the completeness of the infinite set of ground-state Kohn-Sham orbitals,
i.e.,
$\sum_{j=1}^{\infty}\varphi_{j}(\textbf{r}^{\prime})\varphi_{j}^{*}(\textbf{r})=\delta(\textbf{r}^{\prime}-\textbf{r})$.
After some algebra Welakuh (2021) this leads to the frequency-dependent
Sternheimer equations of the following form

$\displaystyle\left(\omega-\hat{h}+\epsilon_{k}+i\eta\right)\varphi_{k}^{(+)}(\textbf{r},\omega)$
$\displaystyle=\delta
v_{\textrm{KS}}(\textbf{r},\omega)\varphi_{k}(\textbf{r}),$ (21)
$\displaystyle\left(\omega+\hat{h}-\epsilon_{k}+i\eta\right)\varphi_{k}^{(-)}(\textbf{r},\omega)$
$\displaystyle=-\varphi_{k}^{*}(\textbf{r})\delta
v_{\textrm{KS}}(\textbf{r},\omega),$ (22)

where the first-order Kohn-Sham potential $\delta
v_{\textrm{KS}}(\textbf{r},\omega)$ is given by

$\displaystyle\delta v_{\textrm{KS}}(\textbf{r},\omega)$ $\displaystyle=\delta
v(\textbf{r},\omega)+\int
d^{3}\textbf{r}^{\prime}f_{\text{Mxc}}^{n}(\textbf{r},\textbf{r}^{\prime},\omega)\delta
n(\textbf{r}^{\prime},\omega)$
$\displaystyle\quad+\sum_{\alpha}f_{\text{Mxc}}^{q_{\alpha}}(\textbf{r},\omega)\delta
q_{\alpha}(\omega).$ (23)

The response of the photon coordinate $\delta q_{\alpha}(\omega)$ in Eq. (10)
to the external potential $\delta v(\textbf{r},\omega)$ can be expressed in
the following form

$\displaystyle\delta q_{\alpha}(\omega)$ $\displaystyle=\delta
q_{\alpha}^{(+)}(\omega)+\delta q_{\alpha}^{(-)}(\omega),$ (24)

where we substituted Eq. (15) into (10). The first-order responses of the
photon coordinates $\delta q_{\alpha}^{(+)}(\omega)$ and $\delta
q_{\alpha}^{(-)}(\omega)$ are given explicitly as

$\displaystyle\delta q_{\alpha}^{(+)}(\omega)$
$\displaystyle=\frac{1}{2\omega_{\alpha}^{2}}\left(\frac{1}{\omega-\omega_{\alpha}+i\eta^{\prime}}\right)\int
d^{3}\textbf{r}^{\prime}g_{\text{M}}^{n_{\alpha}}(\textbf{r}^{\prime})\delta
n(\textbf{r}^{\prime},\omega),$ (25) $\displaystyle\delta
q_{\alpha}^{(-)}(\omega)$
$\displaystyle=-\frac{1}{2\omega_{\alpha}^{2}}\left(\frac{1}{\omega+\omega_{\alpha}+i\eta^{\prime}}\right)\int
d^{3}\textbf{r}^{\prime}g_{\text{M}}^{n_{\alpha}}(\textbf{r}^{\prime})\delta
n(\textbf{r}^{\prime},\omega).$ (26)

To obtain the response of the density and photon coordinate of Eqs. (18) and
(24), we have to solve Eqs. (21)-(26) self-consistently. The self-consistency
in solving these equations becomes evident by noting that the right-hand side
of the Sternheimer equations (21) and (22) depends on the solution through
$\delta v_{\textrm{KS}}(\textbf{r},\omega)$ which in turn depends on $\delta
n(\textbf{r},\omega)$ and $\delta q_{\alpha}(\omega)$. These two quantities
depends on the first-order perturbed Kohn-Sham orbitals
$\varphi_{k}^{(\pm)}(\textbf{r},\omega)$ and photon responses $\delta
q_{\alpha}^{(\pm)}(\omega)$. It is important to note that the first-order
response of the Kohn-Sham orbitals must satisfy the orthogonality condition
with the ground-state Kohn-Sham orbitals Andrade  _et al._ (2007); Hofmann
_et al._ (2018):

$\displaystyle\int
d^{3}\textbf{r}\,\varphi_{k}^{(*)}(\textbf{r})\,\varphi_{k}^{(\pm)}(\textbf{r},\omega)=0\,.$
(27)

From solving the self-consistent Sternheimer equations we can compute the
dynamic polarizability of the coupled system which is given in terms of the
variation of the density

$\displaystyle\alpha_{\mu\nu}(\omega)=\int d^{3}\textbf{r}\,\delta
n_{\nu}(\textbf{r},\omega)\,\textbf{r}_{\mu}\,,$ (28)

and is related to the photo-absorption cross-section as
$\sigma(\omega)=(4\pi\omega/3c)\,\textrm{Tr}\,\bar{\alpha}_{\mu\mu}$ Ullrich
(2011). Since the solutions $\varphi_{k}^{\pm}(\textbf{r},\omega)$ of Eqs.
(26) and (25) are complex-valued, the density response of Eq. (18) become
complex as well. This gives rise to the polarizability
$\alpha_{\mu\nu}(\omega)$ having real and imaginary parts. The imaginary part
of polarizability describes the absorption of radiation, and the real part
defines the refraction properties of the matter system due to a perturbation
from an external electromagnetic field Boyd (1992).

In the decoupling limit of light and matter when
$|\boldsymbol{\lambda}_{\alpha}|\rightarrow 0$, the Sternheimer equations (21)
and (22) still retain the same form, however, the potential $\delta
v_{\textrm{KS}}(\textbf{r},\omega)$ simplifies to that of an electron-only
interacting system as $f_{\text{Mxc}}^{n}\rightarrow f_{\text{Hxc}}^{n}$ and
$f_{\text{Mxc}}^{q_{\alpha}}\rightarrow 0$. Also, the ground-state Kohn-Sham
Hamiltonian in Eqs. (21) and (22) reduces to $\hat{h}=\hat{h}([v,n])$ as
$v_{\textrm{Mxc}}([n,q_{\alpha}];\textbf{r})\rightarrow
v_{\textrm{Hxc}}([n];\textbf{r})$, thus, decoupling the photon contribution of
Eq. (17). The derivation of the Sternheimer scheme for the electron density
and photon displacement coordinate responses in the QEDFT framework due to a
weak external charge current $\delta j_{\alpha}(\omega)$ follows the same
steps as above Welakuh (2021).

Details about the numerical treatment of Eqs. (21) and (22) have been
discussed in the TDDFT framework of the frequency-dependent Sternheimer method
Andrade  _et al._ (2007); Hofmann  _et al._ (2018). Therefore, we only
summarize features in the numerical application of these equations. First, the
positive infinitesimal parameter $\eta$ is required for numerical stability
for the solution of the Sternheimer equations close to the resonance
frequencies as it removes the divergences. It is also necessary to obtain the
imaginary part of the polarizability. In addition, this parameter accounts for
the artificial linewidth that represents the finite lifetimes of the
excitations. Our extension of the Sternheimer method to treat electron-photon
coupled systems introduced the small positive infinitesimal parameter
$\eta^{\prime}$ that enters the self-consistent Sternheimer equations as in
Eqs. (25) and (26). This parameter is necessary to ensure that the poles at
$\omega_{\alpha}$ are finite. In our simulations we found that
$\hbar\eta^{\prime}=0.001$ eV is the ideal value to obtain converged results
and we used $\hbar\eta=0.1$ eV.

For the electron-photon Casida approach, the resulting dimension of the
coupled but truncated matrix is
$\left((N_{v}*N_{c}+M)\times(N_{v}*N_{c}+M)\right)$ where $N_{v}$ and $N_{c}$
denote the number of occupied and unoccupied Kohn-Sham orbitals, respectively
Flick  _et al._ (2019) and $M$ describes the number of photon modes. The
dimensionality of the matrix increases with $N_{c}$ and $M$-photon modes. We
have been so far able to treat a finite matter system coupled to $150,000$
modes with an efficient massive parallel implementation of the Casida equation
Flick  _et al._ (2019); Welakuh (2021). In terms of scaling with system size,
the electron-photon Sternheimer approach is better when compared to the Casida
approach since it still scales the same as the electron-only Sternheimer case
Andrade  _et al._ (2007); Hofmann  _et al._ (2018); Hübener and Giustino
(2014a). This is evident since we can substitute Eqs. (24)-(26) into (23) such
that the complexity rests in solving the Sternheimer equations (21) and (22).
We have implemented the linear-response frequency-dependent Sternheimer
equations (18) and (21)-(26) into the real-space code OCTOPUS Marques  _et
al._ (2003); Andrade  _et al._ (2007).

##  V Applications of the Frequency-dependent Sternheimer approach

Figure 1: Schematic setup of an azulene molecule confined within a high-Q
optical cavity. The cavity field is polarized along the $x$-axis with mode
coupling $\boldsymbol{\lambda}_{\alpha}$ and the photon propagation vector is
along the cavity axis of length $L$ in the $z$-direction. The frequency of the
photon mode is $\omega_{\alpha}$.

Figure 2: Spectrum of an azulene molecule in free space (i.e. $\lambda=0$) and
coupled to a high-Q optical cavity (i.e. $\lambda>0$) showing the line shapes
characteristic of the real and imaginary parts of the polarizability near the
$\pi-\pi^{*}$ resonance at $4.825$ eV. Panel (a) shows the region near the
resonance where the Re $\left\\{\alpha_{xx}(\omega)\right\\}$ is asymmetric
about the resonance while the Im $\left\\{\alpha_{xx}(\omega)\right\\}$ is
symmetric about the resonance. Coupling the cavity mode resonantly to the
$\pi-\pi^{*}$ transition and increasing the coupling strength continuously as
in (b) to (d) result to a Rabi splitting into lower and upper polariton
branches which each have an asymmetric line shape for the different Re
$\left\\{\alpha_{xx}(\omega)\right\\}$.

Figure 3: (a) The real part of the polarizability of azulene showing the
change in the anomalous dispersion in free space ($\lambda_{\alpha}=0$) and
when coupled to a cavity mode ($\lambda_{\alpha}>0$). (b) The analogous
anomalous dispersion in the photon spectrum occurs only when both subsystems
are coupled. Panels (a) and (b) show clearly how this feature can be
controlled by coupling to a cavity mode.

Figure 4: Spectrum of the photon displacement coordinate of an azulene
molecule in free space (i.e. $\lambda=0$) and coupled to a high-Q optical
cavity (i.e. $\lambda>0$). Panel (a) has no response as the photons are
decoupled. Coupling the cavity mode resonantly to the $\pi-\pi^{*}$ transition
and increasing the coupling strength lead to a splitting into lower and upper
polaritonic branches in the photonic spectrum as shown in (b) to (d) for Im
$\left\\{q_{\alpha}(\omega)\right\\}$. The Re
$\left\\{q_{\alpha}(\omega)\right\\}$ for these cases show an antisymmetric
lineshape opposite to Re $\left\\{\alpha_{xx}(\omega)\right\\}$ in Fig. 2.

In this section, we now apply the introduced electron-photon frequency-
dependent Sternheimer approach for studying excited-state properties of
molecular systems coupled to a photon mode or a continuum of modes. This
approach has been validated by comparing the optical absorption spectrum of a
single benzene ring coupled to photons to that obtained using the electron-
photon Casida and time-propagation methods of QEDFT Flick and Narang (2018);
Welakuh (2021). This makes the frequency-dependent Sternheimer method of QEDFT
a valid alternative for studying excited states properties of strongly coupled
light-matter systems.

In the following, we first investigate the physical cavity QED setup in which
a single molecule is strongly coupled to a photon mode of a high-Q cavity
where we expect to capture the hallmark of strong light-matter coupling (Rabi
splitting). In the next setup, we include a large but finite number of photon
modes that simulates the electromagnetic vacuum and investigate situations
where a molecular system couples weakly and strongly to the continuum.

###  V.1 Single molecule strong coupling

The first example studies intrinsic properties of a strongly coupled light-
matter system that is commonly not considered, for instance, the real part of
the polarizability (in Fig. (2)) and the photon displacement field (in Eq.
(4)). These quantities are particularly interesting as they give insight into
the dispersive properties of the coupled system (for the real part of the
polarizability) and how energy is exchanged between the electron-photon system
(for the photon displacement field).

The molecular system considered here is an azulene
($\text{C}_{10}\text{H}_{8}$) molecule which is a bicyclic, nonbenzenoid
aromatic hydrocarbon studied in Ref. Flick  _et al._ (2018b). We describe in
detail how we compute the electronic structure of azulene in App. A. Before
looking at how these observables get modified due to strong light-matter
coupling, we will first present the absorption spectra (obtained from the
imaginary part of the polarizability) of the molecular system strongly coupled
to photons that captures the Rabi splitting between polaritonic peaks Flick
_et al._ (2019); Chikkaraddy  _et al._ (2016).

To study the spectral properties of the coupled system we now confine the
azulene molecule inside an optical high-Q cavity that couples to a photon mode
with increased strength. The cavity field is polarized along the $x$-direction
with a coupling strength $\boldsymbol{\lambda}_{\alpha}$ as shown in Fig. (1).
The optical absorption spectra of the azulene molecule has been computed with
TDDFT which captures the $\pi-\pi^{*}$ transition occurring at around $4.825$
eV Malloci  _et al._ (2007, 2005). In Fig. (2 a.), we show the $x$-component
of the polarizability of the uncoupled azulene molecule. The imaginary part of
the polarizability captures a sharp peak occurring at $4.825$ eV due to the
$\pi-\pi^{*}$ excitation. Based on the Kramers-Kronig relations, an absorption
usually occurs simultaneously with an anomalous dispersion Boyd (1992). The
anomalous dispersion describes a sudden change in the material’s dispersion
spectrum in the vicinity of a resonant absorption. We also find in real part
of the polarizability an anomalous dispersion around the $\pi-\pi^{*}$
excitation which shows how its dispersive properties decreases when the
excitation energy increases. This is characterized by the asymmetric line
shape about this resonance while the imaginary part is symmetric as usually
observed Bonin and Kresin (1997). We now couple a single cavity mode in
resonance to the $\pi-\pi^{*}$ excitation and tune the coupling
$\lambda=|\boldsymbol{\lambda}|$ as in Fig. (2 b-d). For the Im
$\left\\{\alpha_{xx}(\omega)\right\\}$, an increasing coupling strength
results to an increased Rabi splitting of the $\pi-\pi^{*}$ peak into lower
and upper polaritonic branches where the lower branch has more intensity,
compared to the upper polaritonic peak as measured in experiments Chikkaraddy
_et al._ (2016) and not captured by common phenomenological models such as the
Jaynes-Cumming model Wang  _et al._ (2021). This splitting which is a
characteristic of strong light-matter coupling shows how excited-state
properties of matter get modified when strongly coupled to a cavity mode. For
the Re $\left\\{\alpha_{xx}(\omega)\right\\}$, we find for each of the lower
and upper polariton peaks for different $\lambda$, asymmetric line shapes
about their respective excitation energies indicating anomalous dispersion
usually occurs simultaneously with absorption even for strongly coupled
systems. In addition, the anomalous dispersion can be controlled for strongly
coupled systems by varying the coupling strength. This is clearly shown in
Fig. (3 a) where the anomalous dispersion (in particular for the lower
polariton) is smaller for the coupled case when compared to the uncoupled
result. The emergence of polaritonic features in the Re
$\left\\{\alpha_{xx}(\omega)\right\\}$ highlights that the dispersion
properties of the matter system becomes modified due to strong light-matter
coupling. The modification of dispersion properties for strongly coupled
light-matter systems has potential in controlling optical dipole traps. This
can be made clear by considering the interaction potential of the induced
dipole moment normally expressed as
$U=-\frac{1}{2\varepsilon_{0}c}\text{Re}\left\\{\alpha_{xx}(\omega)\right\\}I$,
where $I$ is the field intensity Grimm  _et al._ (2000). The standard approach
for realizing optical dipole traps is by laser detuning from a specific
resonance of the bare matter system, for instance, laser detuning from an
atomic resonance such that the dipole potential minima occur at regions with
maximum intensity for red-detuned traps Grimm  _et al._ (2000). For
polaritonic resonances that emerge in strongly coupled light-matter systems,
the optical dipole traps that can be realized by detuning the external field
from these polaritonic resonances can be controlled by strongly coupling to
the photon field. This is evident in Fig. (2) where the Re
$\left\\{\alpha_{xx}(\omega)\right\\}$ is modified under strong coupling and
highlights a new perspective with potential applications in engineering
optical dipole traps for neutral atoms or molecules.

Next, we study the spectral properties of the photon field when we probe the
matter subsystem. This observable $\delta q_{\alpha}(\omega)$ is now
accessible since we treat the photon field as a dynamical part of the coupled
light-matter system. We note that the displacement field in this case
represents a mixed (matter and photon) spectroscopic observable since its
response function $\chi_{n}^{q_{\alpha}}$ is a commutator between photonic and
electronic quantities Flick  _et al._ (2019). The observable $\delta
q_{\alpha}(\omega)$ indicates how the photon field reacts in a standard
absorption or emission measurement when the system is probed by an external
field represented by the potential $\delta v(\textbf{r},\omega)$. In Fig. (4),
we show the spectrum of the photon displacement coordinate in free space (when
$\lambda=0$) and coupled to a cavity mode (when $\lambda>0$). As expected the
free space case has no response since light and matter decouple and we have
access only to the observables in Fig. (2 a). However, coupling to the photon
mode and increasing the coupling strength $\lambda>0$ we observe in the
imaginary part of $\delta q_{\alpha}(\omega)$ a Rabi splitting into lower and
upper polaritons peaks. The polaritonic peaks are asymmetric about the
$\pi-\pi^{*}$ excitation energy to which the mode was initially coupled to and
the lower polariton peaks are negative with more intensity compared to the
upper polarition. Physically, this result highlights that excitations due to
an external perturbation from $\delta v(\textbf{r},\omega)$ can be exchanged
between the coupled subsystems and the hybrid light-matter features occur not
only in the matter subsystem but also in the photon subsystem due to the self-
consistent interaction. For the Re $\left\\{q_{\alpha}(\omega)\right\\}$, we
also find for each of the lower and upper polariton branches an asymmetric
line shape about the energies of the respective polariton peaks with varying
strengths for different $\lambda$. In analogy to the Re
$\left\\{\alpha_{xx}(\omega)\right\\}$ where the anomalous dispersion get
modified due to strong light-matter coupling, the same holds true for the
anomalous region in the spectrum of Re $\left\\{q_{\alpha}(\omega)\right\\}$
as shown in Fig. (3 b). Due to the self-consistent back-reaction between
subsystems, we expect that the Re $\left\\{q_{\alpha}(\omega)\right\\}$ can be
made to influence the optical dipole potential thereby controlling how the
matter subsystem is trapped in the field. It is important to note that for the
responses of the subsystems, the excitation energies of the strongly coupled
system are the same but with differing oscillator strengths (see App. A). The
results presented here demonstrates that the electron-photon Sternheimer
formalism is able to describe excited-state properties of strong light-matter
coupled system.

###  V.2 Changes in the matter spectral features

In this section, we consider the case where a molecular system is coupled
explicitly to a wide range of photon modes and show how spectral features of
the system change when we effectively increase its coupling to the continuum
of the electromagnetic field. This computation will at the same time show the
advantages the Sternheimer approach has over the Casida approach in terms of
scaling with the number of photon modes.

Figure 5: (a) Photo-absorption cross-section of a LiH molecule coupled to
500,000 photon modes (blue dashed) in a quasi one-dimensional cavity and its
comparison to the uncoupled case (red solid). (b) Zoom-in view to the
$X^{\,1}\Sigma^{+}\rightarrow A^{\,1}\Sigma^{+}$ transition around $3.2$ eV
where we observe a slight deviation in the peak amplitude between the
uncoupled and the case coupled to the continuum.

Figure 6: (a) Photo-absorption spectrum of a LiH molecule coupled to a
continuum of photon modes where the coupling is effectively enhanced by
changing the cavity volume via the area $L_{y}L_{z}$. The Lorentzian line
shapes turn into Fano line shapes for increasing effective coupling strength.
(b) Zoom-in view to the $X^{\,1}\Sigma^{+}\rightarrow A^{\,1}\Sigma^{+}$
transition around $3.2$ eV where we observe clearly the asymmetry of the Fano
resonances when compared to Lorentzian lineshape (blue solid line).

We now consider as matter system a lithium hydride (LiH) molecule coupled to a
wide range of photon modes that sample densely the electromagnetic vacuum.
Since the Sternheimer approach for electron-only system is known to scale
favorably with the system size Andrade  _et al._ (2007); Hofmann  _et al._
(2018), the focus here will be to demonstrate that the photon modes do not add
to this scaling. Here we sample modes of a quasi one-dimensional mode space by
employing the coupling
$\boldsymbol{\lambda}_{\alpha}=\sqrt{\frac{2}{\epsilon_{0}L_{x}L_{y}L_{z}}}\text{sin}(\omega_{\alpha}/c\,x_{0})\textbf{e}_{x}$,
where $x_{0}=L_{x}/2$ is the position of the molecule in $x$-direction and
$\omega_{\alpha}=\alpha c\pi/L_{x}$ are the frequencies of the modes Flick
_et al._ (2019). The volume $V=L_{x}L_{y}L_{z}$ with $L_{x}=3250$ $\mu$m,
$L_{y}=10.58$ Å and $L_{z}=2.65$ Å are chosen such that the sampled modes
couple weakly to the molecular system and we assume a constant mode function
in the $y$ and $z$-directions.

In this first example, we couple the molecule to 500,000 photon modes of a
one-dimensional mode space with an energy cut-off of 190.74 eV and a spacing
between modes of 0.38 meV. Sampling the continuum of modes serve to constitute
the line width of the excitations and also represent dissipation channels in
the coupled system Welakuh  _et al._ (2021); Flick  _et al._ (2019). The one-
dimensional sampling of mode frequencies that couple weakly to the matter
subsystem will not capture the actual three-dimensional lifetimes. In the
matter-only (uncoupled) case, we use a broadening $\hbar\eta=0.1$ eV (as in
Sec. V.1) to account for the finite lifetime of the excited states. When the
molecule is coupled weakly to the photon continuum, we obtain the same
spectral broadening as the uncoupled case. The results of this calculation is
shown in Fig. (5) where we compare the photo-absorption cross-section of the
uncoupled LiH molecule and the case when it is weakly coupled to 500,000
photon modes. We find that the two results are qualitatively the same which is
evident for the lowest electronic transition $X^{\,1}\Sigma^{+}\rightarrow
A^{\,1}\Sigma^{+}$ around $3.2$ eV that corresponds to an electronic
transition from the bonding to the antibonding $\sigma$-orbital Villaume  _et
al._ (2010); Stwalley and Zemke (1993). This result shows that the weak
coupling of the molecule to the continuum of modes reproduces the results of
the matter-only case. We note that obtaining this result using the electron-
photon Casida approach will increase the computational cost drastically even
for the case of coupling to 100,000 photon modes. Computationally, this result
demonstrates that the electron-photon Sternheimer method scales favorably not
only with system size but also with the number of photon modes.

Now, we effectively enhance the coupling strength $|\boldsymbol{\lambda}|$ by
reducing the cavity volume along the $y$\- and $z$-direction. The results are
shown in Fig. 6 where the blue line is the result shown in Fig. 5 that has a
Lorentzian profile. We find that upon reducing the area $L_{y}L_{z}$, this
effectively enhances the coupling to the photon continuum such that the
symmetric Lorentz line shapes turn into asymmetric Fano line shapes. Fano
resonances occur due to the interference of discrete quantum states with a
continuum of states Ott  _et al._ (2013); Limonov  _et al._ (2017). The
asymmetry is characterized as a ratio of the transition amplitude to a given
discrete state and that of a transition to a continuum state Miroshnichenko
_et al._ (2010). As this ratio becomes finite due to strong coupling to the
continuum, this indicates the onset of a competition between constructive and
destructive interference that gives rise to the asymmetric line shape Welakuh
and Narang (2021). Also, the broadening of the spectra (see Fig. (5 b)) and
decrease in amplitude is a consequence of the interference Welakuh and Narang
(2021). These results show the changes in the spectral features of excited-
states of a matter system strongly coupled to the electromagnetic continuum.
Thus, the electron-photon Sternheimer approach is a valid alternative method
for studying excited-state properties of real systems strongly interacting
with the quantized electromagnetic field.

##  VI Conclusion and outlook

In this work we presented a linear-response method that solves the response
equations of non-relativistic QED in the length gauge setting. The approach is
based on the Sternheimer equation formulated within the framework of QEDFT
that is capable of computing excited-state properties of strongly coupled
light-matter systems. This approach serves as an alternative linear-response
method for studying response properties of large systems coupled to the
quantized electromagnetic field since it scales favorably with the system size
as it utilizes only the occupied Kohn-Sham orbitals and it also scales
favorably with the number of photon modes. Using the Sternheimer approach we
computed different observables of strongly coupled systems. These observables
showed how both the dispersion and absorption properties of the matter system
changes with potential applications in modifying and controlling optical
dipole traps. Also, we showed examples where we lift the restriction to one
cavity mode in the dipole approximation and sampled densely the
electromagnetic continuum. In one case we showed that when a LiH molecule is
weakly coupled to the photon continuum, we reproduce the free space absorption
spectrum of the molecule. When the coupling strength between the light and
matter is effectively enhanced, we find changes in the absorption spectrum as
symmetric Lorentzian line shapes turn into asymmetric Fano line shapes.

The electron-photon Sternheimer method presented here is a suitable approach
for studying excited-state properties of large systems coupled to a single
mode or to the electromagnetic continuum. In the fast-growing field of
polaritonic chemistry where there is an ongoing debate about the mesoscopic
scale of quantum-collectively of coupled molecules Sidler  _et al._ (2021a,
b), ab-initio methods such as the electron-photon Sternheimer method become
desirable to capture intricate details of the complex interactions between the
coupled subsystems. Another important property of the Sternheimer approach is
that it can be generalized to higher orders to obtain higher order
polarizabilities by solving a hierarchy of Sternheimer equations Andrade  _et
al._ (2007). For the coupled electron-photon system, this will give access to
higher order polarizabilities with signatures of strong light-matter coupling.

##  VII Acknowledgments

We acknowledge financial support from the European Research Council
(ERC-2015-AdG-694097) and the SFB925 ”Light induced dynamics and control of
correlated quantum systems”. This work was supported by the Excellence Cluster
”CUI: Advanced Imaging of Matter” of the Deutsche Forschungsgemeinschaft
(DFG), EXC 2056, project ID 390715994. The Flatiron Institute is a division of
the Simons Foundation.

##  Appendix A Photo-absorption cross-section of the azulene molecule

Figure 7: Comparison of the photo-absorption cross-section of an azulene
molecule using the electron-photon Sternheimer and Casida methods. (a)
Absorption cross-section in free space (i.e. $\lambda_{\alpha}=0$). Coupling
the cavity mode resonantly to the $\pi-\pi^{*}$ transition and increasing the
coupling strength continuously as in (b) to (d) leads to a Rabi splitting into
lower and upper polariton branches. Both approaches are in good agreement for
the computed spectra.

Figure 8: The spectrum of the imaginary parts of the polarizability and photon
displacement coordinate of an azulene molecule in free space and in an optical
cavity. (a) The free space case where light decouples from matter and there is
no response from the photon subsystem. Panels (b)-(d) shows the light-matter
coupled case where both observables are shown to have the same excitation
energies but different oscillator strengths.

To describe the azulene molecule studied in Sec. V, we compute the electronic
structure of the azulene molecule using the same setup as that in Ref. Flick
_et al._ (2019) for the benzene molecule. That is, a cylindrical real space
grid of $8$ Å in length with the radius of 6 Å in the $x$-$y$ plane and a
spacing $\Delta x=\Delta y=\Delta z=0.22$ Å is used. We treat the core
electrons of the carbon and hydrogen atoms using the local-density
approximation Troullier-Martins pseudopotentials Troullier and Martins (1991).
Since the Sternheimer method considers only occupied orbitals, the $48$
valence electrons that amounts to $24$ doubly occupied orbitals are used in
the computation. When we use the Casida approach, we include 500 unoccupied
states which amounts to $N_{v}*N_{c}=24*500=12000$ pairs of occupied-
unoccupied states. To account for the electron-electron and electron-photon
interaction in the linear-response QEDFT framework Flick  _et al._ (2019), we
apply the adiabatic LDA (ALDA) to the Hartree exchange-correlation kernel
(i.e. $f^{n}_{\text{Hxc}}\rightarrow f^{n}_{\text{Hxc,ALDA}}$) and the photon
random-phase approximation to the photon exchange-correlation kernels (i.e.
$f^{n}_{\text{pxc}},f^{q_{\alpha}}_{\text{pxc}}\rightarrow
f^{n}_{\text{p}},f^{q_{\alpha}}_{\text{p}}$), respectively. We compute the
photo-absorption cross-section of the molecule in free space and inside the
high-Q optical cavity and the results are shown in Fig. (7. This result is
related to the imaginary part of Fig. (2) and it shows more excitations at
higher energies.

In the linear-response setting of non-relativistic QED in the length gauge
Flick  _et al._ (2019), the response quantities that arise have the same
excitation energies but different oscillator strengths as shown in Fig. (8)
for two observables of interest. This is because their respective response
functions have the same excitation energies and different oscillator strengths
Flick  _et al._ (2019).

## References

  * Hutchison  _et al._ (2012) James A. Hutchison, Tal Schwartz, Cyriaque Genet, Eloïse Devaux, and Thomas W. Ebbesen, “Modifying chemical landscapes by coupling to vacuum fields,” Angewandte Chemie International Edition 51, 1592–1596 (2012).
  * Thomas  _et al._ (2016) Anoop Thomas, Jino George, Atef Shalabney, Marian Dryzhakov, Sreejith J. Varma, Joseph Moran, Thibault Chervy, Xiaolan Zhong, Eloïse Devaux, Cyriaque Genet, James A. Hutchison, and Thomas W. Ebbesen, “Ground-state chemical reactivity under vibrational coupling to the vacuum electromagnetic field,” Angewandte Chemie International Edition 55, 11462–11466 (2016).
  * Schäfer  _et al._ (2021) Christian Schäfer, Johannes Flick, Enrico Ronca, Prineha Narang, and Angel Rubio, “Shining light on the microscopic resonant mechanism responsible for cavity-mediated chemical reactivity,” arXiv:2104.12429 (2021).
  * Coles  _et al._ (2014) David M. Coles, Yanshen Yang, Yaya Wang, Richard T. Grant, Robert A. Taylor, Semion K. Saikin, Alán Aspuru-Guzik, David G. Lidzey, Joseph Kuo-Hsiang Tang, and Jason M. Smith, “Strong coupling between chlorosomes of photosynthetic bacteria and a confined optical cavity mode,” Nature Communications 5, 5561 (2014).
  * Zhong  _et al._ (2016) Xiaolan Zhong, Thibault Chervy, Shaojun Wang, Jino George, Anoop Thomas, James A. Hutchison, Eloise Devaux, Cyriaque Genet, and Thomas W. Ebbesen, “Non-radiative energy transfer mediated by hybrid light-matter states,” Angewandte Chemie International Edition 55, 6202–6206 (2016).
  * Feist and Garcia-Vidal (2015) Johannes Feist and Francisco J. Garcia-Vidal, “Extraordinary exciton conductance induced by strong coupling,” Phys. Rev. Lett. 114, 196402 (2015).
  * Stranius  _et al._ (2018) Kati Stranius, Manuel Hertzog, and Karl Börjesson, “Selective manipulation of electronically excited states through strong light-matter interactions,” Nature Communications 9, 2273 (2018).
  * Chikkaraddy  _et al._ (2016) Rohit Chikkaraddy, Bart de Nijs, Felix Benz, Steven J. Barrow, Oren A. Scherman, Edina Rosta, Angela Demetriadou, Peter Fox, Ortwin Hess, and Jeremy J. Baumberg, “Single-molecule strong coupling at room temperature in plasmonic nanocavities,” Nature 535, 127–130 (2016).
  * Bellessa  _et al._ (2004) J. Bellessa, C. Bonnand, J. C. Plenet, and J. Mugnier, “Strong coupling between surface plasmons and excitons in an organic semiconductor,” Phys. Rev. Lett. 93, 036404 (2004).
  * Ebbesen (2016) Thomas W. Ebbesen, “Hybrid light–matter states in a molecular and material science perspective,” Accounts of Chemical Research 49, 2403–2412 (2016).
  * Santhosh  _et al._ (2016) Kotni Santhosh, Ora Bitton, Lev Chuntonov, and Gilad Haran, “Vacuum rabi splitting in a plasmonic cavity at the single quantum emitter limit,” Nat. Commun. 7, 11823 (2016).
  * Dicke (1954) R. H. Dicke, “Coherence in spontaneous radiation processes,” Phys. Rev. 93, 99 (1954).
  * Jaynes and Cummings (1963) E.T. Jaynes and F.W. Cummings, “Comparison of quantum and semiclassical radiation theories with application to the beam maser,” Proceedings of the IEEE 51, 89–109 (1963).
  * Garraway (2011) Barry M. Garraway, “The dicke model in quantum optics: Dicke model revisited,” Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 369, 1137–1155 (2011).
  * Sidler  _et al._ (2020) Dominik Sidler, Michael Ruggenthaler, Heiko Appel, and Angel Rubio, “Chemistry in quantum cavities: Exact results, the impact of thermal velocities, and modified dissociation,” J. Phys. Chem. Lett. 11, 7525–7530 (2020).
  * Sidler  _et al._ (2021a) Dominik Sidler, Christian Schäfer, Michael Ruggenthaler, and Angel Rubio, “Polaritonic chemistry: Collective strong coupling implies strong local modification of chemical properties,” J. Phys. Chem. Lett. 12, 508–516 (2021a).
  * Schäfer  _et al._ (2020) Christian Schäfer, Michael Ruggenthaler, Vasil Rokaj, and Angel Rubio, “Relevance of the quadratic diamagnetic and self-polarization terms in cavity quantum electrodynamics,” ACS Photonics 7, 975–990 (2020).
  * Ruggenthaler  _et al._ (2014) Michael Ruggenthaler, Johannes Flick, Camilla Pellegrini, Heiko Appel, Ilya V. Tokatly, and Angel Rubio, “Quantum-electrodynamical density-functional theory: Bridging quantum optics and electronic-structure theory,” Phys. Rev. A 90, 012508 (2014).
  * Flick  _et al._ (2015) Johannes Flick, Michael Ruggenthaler, Heiko Appel, and Angel Rubio, “Kohn-sham approach to quantum electrodynamical density-functional theory: Exact time-dependent effective potentials in real space,” Proc. Natl. Acad. Sci. U. S. A. 112, 15285–15290 (2015).
  * Haugland  _et al._ (2020) Tor S. Haugland, Enrico Ronca, Eirik F. Kjonstad, Angel Rubio, and Henrik Koch, “Coupled cluster theory for molecular polaritons: Changing ground and excited states,” Phys. Rev. X 10, 041043 (2020).
  * Spohn (2004) Herbert Spohn,  _Dynamics of charged particles and their radiation field_ (Cambridge university press, 2004).
  * Ruggenthaler  _et al._ (2018) Michael Ruggenthaler, Nicolas Tancogne-Dejean, Johannes Flick, Heiko Appel, and Angel Rubio, “From a quantum-electrodynamical light–matter description to novel spectroscopies,” Nature Reviews Chemistry 2, 0118 (2018).
  * Flick  _et al._ (2018a) Johannes Flick, Nicholas Rivera, and Prineha Narang, “Strong light-matter coupling in quantum chemistry and quantum photonics,” Nanophotonics 7, 1479–1501 (2018a).
  * Flick  _et al._ (2018b) Johannes Flick, Christian Schäfer, Michael Ruggenthaler, Heiko Appel, and Angel Rubio, “Ab initio optimized effective potentials for real molecules in optical cavities: Photon contributions to the molecular ground state,” ACS Photonics 5, 992–1005 (2018b).
  * Flick  _et al._ (2019) Johannes Flick, Davis M. Welakuh, Michael Ruggenthaler, Heiko Appel, and Angel Rubio, “Light-matter response in nonrelativistic quantum electrodynamics,” ACS Photonics 6, 2757–2778 (2019).
  * Casida (1996) Mark Casida, “Time-dependent density functional response theory of molecular systems: Theory, computational methods, and functionals,” in  _Recent Developments and Applications of Modern Density Functional Theory_ , edited by J.M. Seminario (Elsevier, Amsterdam, 1996).
  * Flick and Narang (2020) Johannes Flick and Prineha Narang, “Ab initio polaritonic potential-energy surfaces for excited-state nanophotonics and polaritonic chemistry,” J. Chem. Phys. 153, 094116 (2020).
  * Wang  _et al._ (2021) Derek S. Wang, Tomáš Neuman, Johannes Flick, and Prineha Narang, “Light-matter interaction of a molecule in a dissipative cavity from first principles,” J. Chem. Phys. 154, 104109 (2021).
  * Flick and Narang (2018) Johannes Flick and Prineha Narang, “Cavity-correlated electron-nuclear dynamics from first principles,” Phys. Rev. Lett. 121, 113002 (2018).
  * Welakuh (2021) Davis M. Welakuh,  _Ab initio Strong Light-matter interaction in Non-relativistic Quantum Electrodynamics (QED)_ , Ph.D. thesis, Universität Hamburg (2021).
  * Jestädt  _et al._ (2020) René Jestädt, Michael Ruggenthaler, Micael J. T. Oliveira, Angel Rubio, and Heiko Appel, “Light-matter interactions within the ehrenfest-maxwell-pauli-kohn-sham framework: fundamentals, implementation, and nano-optical applications,” Advances in Physics 68, 225–333 (2020).
  * Sternheimer (1951) R. Sternheimer, “On nuclear quadrupole moments,” Phys. Rev. 84, 244–253 (1951).
  * Baroni  _et al._ (1987) Stefano Baroni, Paolo Giannozzi, and Andrea Testa, “Elastic constants of crystals from linear-response theory,” Phys. Rev. Lett. 59, 2662 (1987), 10.1103/PhysRevLett.59.2662.
  * Hübener and Giustino (2014a) Hannes Hübener and Feliciano Giustino, “Time-dependent density-functional theory using atomic orbitals and the self-consistent sternheimer equation,” Phys. Rev. B 89, 085129 (2014a), 10.1103/PhysRevB.89.085129.
  * Hübener and Giustino (2014b) Hannes Hübener and Feliciano Giustino, “Linear optical response of finite systemsusing multishift linear system solvers,” J. Chem. Phys. 141, 044117 (2014b), 10.1063/1.4890736.
  * Giustino  _et al._ (2010) Feliciano Giustino, Marvin L. Cohen, and Steven G. Louie, “Gw method with the self-consistent sternheimer equation,” Phys. Rev. B 81, 115105 (2010).
  * Andrade  _et al._ (2007) Xavier Andrade, Silvana Botti, Miguel A. L. Marques, and Angel Rubio, “Time-dependent density functional theory scheme for efficient calculations of dynamic (hyper)polarizabilities,” J. Chem. Phys. (2007), 10.1063/1.2733666.
  * Hofmann  _et al._ (2018) Fabian Hofmann, Ingo Schelter, and Stephan Kümmel, “Linear response time-dependent density functional theory without unoccupied states: The kohn-sham-sternheimer scheme revisited,” J. Chem. Phys. 149, 024105 (2018).
  * Loudon (2000) Rodney Loudon,  _The Quantum Theory of Light_ (Oxford University Press, 2000).
  * Abedi  _et al._ (2018) Ali Abedi, Elham Khosravi, and Ilya V. Tokatly, “Shedding light on correlated electron-photon states using the exact factorization,” Eur. Phys. J. B 91, 194 (2018).
  * Milonni (1993) Peter W. Milonni,  _The Quantum Vacuum: An Introduction to Quantum Electrodynamics_ (Academic Press, Boston, 1993).
  * Flick  _et al._ (2017) Johannes Flick, Michael Ruggenthaler, Heiko Appel, and Angel Rubio, “Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (qed) chemistry,” Proceedings of the National Academy of Sciences 114, 3026–3034 (2017).
  * Welakuh  _et al._ (2021) Davis M. Welakuh, Michael Ruggenthaler, Mary-Leena M. Tchenkoue, Heiko Appel, and Angel Rubio, “Down-conversion processes in ab-initio non-relativistic quantum electrodynamics,” Phys. Rev. Research 3, 033067 (2021).
  * Rokaj  _et al._ (2018) Vasil Rokaj, Davis M Welakuh, Michael Ruggenthaler, and Angel Rubio, “Light–matter interaction in the long-wavelength limit: no ground-state without dipole self-energy,” Journal of Physics B: Atomic, Molecular and Optical Physics 51, 034005 (2018).
  * Greiner and Reinhardt (1996) Walter Greiner and Joachim Reinhardt,  _Field quantization_ (Springer, 1996).
  * Rokaj  _et al._ (2022) Vasil Rokaj, Michael Ruggenthaler, Florian G. Eich, and Angel Rubio, “The free electron gas in cavity quantum electrodynamics,” Phys. Rev. Research 4, 013012 (2022).
  * Ullrich (2011) Carsten A Ullrich,  _Time-dependent density-functional theory: concepts and applications_ (OUP Oxford, 2011).
  * Marquesa  _et al._ (2012) M. Marquesa, N. Maitra, F. Nogueira, E. Gross, and A. Rubio,  _Fundamentals of Time-Dependent Density Functional Theory_ (Lecture Notes in Physics, Springer Berlin Heidelberg, 2012).
  * Ruggenthaler (2015) Michael Ruggenthaler, “Ground-State Quantum-Electrodynamical Density-Functional Theory,” ArXiv e-prints (2015), arXiv:1509.01417 [quant-ph] .
  * Boyd (1992) Robert W. Boyd,  _Nonlinear Optics_ (Academic, New York, 1992).
  * Marques  _et al._ (2003) Miguel A.L. Marques, Alberto Castro, George F. Bertsch, and Angel Rubio, “octopus: a first-principles tool for excited electron–ion dynamics,” Computer Physics Communications 151, 60 – 78 (2003).
  * Malloci  _et al._ (2007) G. Malloci, C. Joblin, and G. Mulas, “On-line database of the spectral properties of polycyclic aromatic hydrocarbons,” Chem. Phys. 332, 353-359 (2007), 10.1016/j.chemphys.2007.01.001.
  * Malloci  _et al._ (2005) G. Malloci, G. Mulas, G. Cappellini, V. Fiorentini, and I. Porceddu, “Theoretical electron affinities of pahs and electronic absorption spectra of their mono-anions,” Astron. Astrophys. 432, 585-594 (2005), 10.1051/0004-6361:20042246.
  * Bonin and Kresin (1997) Keith D. Bonin and Vitaly V. Kresin,  _Electric-Dipole Polarizabilities of Atoms, Molecules, and Clusters_ (World Scientific, 1997).
  * Grimm  _et al._ (2000) Rudolf Grimm, Matthias Weidemüller, and Yurii B.Ovchinnikov, “Optical dipole traps for neutral atoms,” Advances in atomic, molecular, and optical physics 42, 95–170. (2000).
  * Villaume  _et al._ (2010) Sebastien Villaume, Trond Saue, and Patrick Norman, “Linear complex polarization propagator in a four-component kohn-sham framework,” J. Chem. Phys. 133, 064105 (2010).
  * Stwalley and Zemke (1993) William C. Stwalley and Warren T. Zemke, “Spectroscopy and structure of the lithium hydride diatomic molecules and ions,” J. Phys. Chem. Ref. Data 22, 87 (1993).
  * Ott  _et al._ (2013) Christian Ott, Andreas Kaldun, Philipp Raith, Kristina Meyer, Martin Laux, Jörg Evers, Christoph H. Keitel, Chris H. Greene, and Thomas Pfeifer, “Lorentz meets fano in spectral line shapes: A universal phase and its laser control,” Science 340, 716–720 (2013).
  * Limonov  _et al._ (2017) Mikhail F. Limonov, Mikhail V. Rybin, Alexander N. Poddubny, and Yuri S. Kivshar, “Fano resonances in photonics,” Nature Photonics 11, 543–554 (2017).
  * Miroshnichenko  _et al._ (2010) Andrey E. Miroshnichenko, Sergej Flach, and Yuri S. Kivshar, “Fano resonances in nanoscale structures,” Rev. Mod. Phys. 82, 2257 (2010).
  * Welakuh and Narang (2021) Davis M. Welakuh and Prineha Narang, “Transition from lorentz to fano spectral line shapes in non-relativistic quantum electrodynamics,” arXiv:2112.05114 (2021).
  * Sidler  _et al._ (2021b) Dominik Sidler, Michael Ruggenthaler, Christian Schäfer, Enrico Ronca, and Angel Rubio, “A perspective on ab initio modeling of polaritonic chemistry: The role of non-equilibrium effects and quantum collectivity,” arXiv:2108.12244 (2021b).
  * Troullier and Martins (1991) N. Troullier and Jose Luis Martins, “Efficient pseudopotentials for plane-wave calculations,” Phys. Rev. B 43, 1993–2006 (1991).


* (120) Scaringi, S., Groot, P. J., Knigge, C., Lasota, J. P., de Martino, D., Cavecchi, Y., Buckley, D. A. H., and Camisassa, M. E. (2022c).  Triggering micronovae through magnetically confined accretion flows in accreting white dwarfs.  MNRAS, 514(1):L11–L15.
  * Scaringi et al., (2017) Scaringi, S., Maccarone, T. J., D’Angelo, C., Knigge, C., and Groot, P. J. (2017).  Magnetically gated accretion in an accreting ‘non-magnetic’ white dwarf.  Nature, 552(7684):210–213.
  * Scepi et al., (2019) Scepi, N., Dubus, G., and Lesur, G. (2019).  Magnetic wind-driven accretion in dwarf novae.  A&A, 626:A116.
  * Scepi et al., (2018) Scepi, N., Lesur, G., Dubus, G., and Flock, M. (2018).  Impact of convection and resistivity on angular momentum transport in dwarf novae.  A&A, 609:A77.
  * Schaefer et al., (2022) Schaefer, B. E., Pagnotta, A., and Zoppelt, S. (2022).  Discovery of extreme, roughly daily superflares on the recurrent nova V2487 Oph.  MNRAS, 512(2):1924–1943.
  * Schreiber, (2007) Schreiber, M. R. (2007).  The statistical significance of the superhump signal in U Geminorum.  A&A, 466(3):1025–1030.
  * Schreiber, (2013) Schreiber, M. R. (2013).  One Good Measure.  Science, 340(6135):932–933.
  * Schreiber et al., (2003) Schreiber, M. R., Hameury, J. M., and Lasota, J. P. (2003).  Delays in dwarf novae I: The case of SS Cygni.  A&A, 410:239–252.
  * Schreiber et al., (2004) Schreiber, M. R., Hameury, J. M., and Lasota, J. P. (2004).  Delays in dwarf novae II: VW Hyi, the tidal instability and enhanced mass transfer models.  A&A, 427:621–635.
  * Shakura and Sunyaev, (1973) Shakura, N. I. and Sunyaev, R. A. (1973).  Black holes in binary systems. Observational appearance.  A&A, 24:337–355.
  * Simonsen, (2011) Simonsen, M. (2011).  The Z CamPaign: Year 1.  J. Am. Assoc. Var. Star Obs., 39(1):66.
  * Smak, (1984) Smak, J. (1984).  Accretion in cataclysmic binaries. IV. Accretion disks in dwarf novae.  Ac. Ast., 34:161–189.
  * (132) Smak, J. (2009a).  New Interpretation of Superhumps.  Ac. Ast., 59(1):121–130.
  * (133) Smak, J. (2009b).  On the Origin of Tilted Disks and Negative Superhumps.  Ac. Ast., 59(4):419–430.
  * Smak, (2017) Smak, J. (2017).  DW UMa and the Irradiation Modulated Mass Transfer Model for Superhumps.  Ac. Ast., 67(3):273–280.
  * (135) Smak, J. (2020a).  On the Critical Mass Ratio for the 3:1 Resonance in Accretion Disks.  Ac. Ast., 70(4):313–315.
  * (136) Smak, J. (2020b).  On the Periods and Nature of Superhumps.  Ac. Ast., 70(4):317–321.
  * Smak and Waagen, (2004) Smak, J. and Waagen, E. O. (2004).  The 1985 Superoutburst of U Geminorum. Detection of Superhumps.  Ac. Ast., 54:433–442.
  * Smak, (1991) Smak, J. I. (1991).  On the models for superoutbursts in dwarf novae of the SU UMa type.  Ac. Ast., 41:269–277.
  * Spruit and Taam, (1993) Spruit, H. C. and Taam, R. E. (1993).  An Instability Associated with a Magnetosphere-Disk Interaction.  ApJ, 402:593.
  * Szkody et al., (2013) Szkody, P., Albright, M., Linnell, A. P., Everett, M. E., McMillan, R., Saurage, G., Huehnerhoff, J., Howell, S. B., Simonsen, M., and Hunt-Walker, N. (2013).  A Study of the Unusual Z Cam Systems IW Andromedae and V513 Cassiopeia.  PASP, 125(934):1421.
  * Tetarenko et al., (2020) Tetarenko, B. E., Dubus, G., Marcel, G., Done, C., and Clavel, M. (2020).  Thermally driven disc winds as a mechanism for X-ray irradiation heating in black hole X-ray binaries: the case study of GX339-4.  MNRAS, 495(4):3666–3682.
  * Tetarenko et al., (2018) Tetarenko, B. E., Lasota, J. P., Heinke, C. O., Dubus, G., and Sivakoff, G. R. (2018).  Strong disk winds traced throughout outbursts in black-hole X-ray binaries.  Nature, 554(7690):69–72.
  * Thomas et al., (2022) Thomas, J. K., Charles, P. A., Buckley, D. A. H., Kotze, M. M., Lasota, J.-P., Potter, S. B., Steiner, J. F., and Paice, J. A. (2022).  Large optical modulations during 2018 outburst of MAXI J1820 + 070 reveal evolution of warped accretion disc through X-ray state change.  MNRAS, 509(1):1062–1074.
  * Townsend et al., (2017) Townsend, L. J., Kennea, J. A., Coe, M. J., McBride, V. A., Buckley, D. A. H., Evans, P. A., and Udalski, A. (2017).  The 2016 super-Eddington outburst of SMC X-3: X-ray and optical properties and system parameters.  MNRAS, 471(4):3878–3887.
  * Tsygankov et al., (2017) Tsygankov, S. S., Doroshenko, V., Lutovinov, A. A., Mushtukov, A. A., and Poutanen, J. (2017).  SMC X-3: the closest ultraluminous X-ray source powered by a neutron star with non-dipole magnetic field.  A&A, 605:A39.
  * van Paradijs, (1996) van Paradijs, J. (1996).  On the Accretion Instability in Soft X-Ray Transients.  ApJLett, 464:L139.
  * van Paradijs and McClintock, (1994) van Paradijs, J. and McClintock, J. E. (1994).  Absolute visual magnitudes of low-mass X-ray binaries.  A&A, 290:133–136.
  * Veledina et al., (2023) Veledina, A., Muleri, F., Poutanen, J., Podgorný, J., Dovčiak, M., Capitanio, F., Churazov, E., De Rosa, A., Di Marco, A., Forsblom, S., Kaaret, P., Krawczynski, H., La Monaca, F., Loktev, V., Lutovinov, A. A., Molkov, S. V., Mushtukov, A. A., Ratheesh, A., Rodriguez Cavero, N., Steiner, J. F., Sunyaev, R. A., Tsygankov, S. S., Zdziarski, A. A., Bianchi, S., Bright, J. S., Bursov, N., Costa, E., Egron, E., Garcia, J. A., Green, D. A., Gurwell, M., Ingram, A., Kajava, J. J. E., Kale, R., Kraus, A., Malyshev, D., Marin, F., Matt, G., McCollough, M., Mereminskiy, I. A., Nizhelsky, N., Piano, G., Pilia, M., Pittori, C., Rao, R., Righini, S., Soffitta, P., Shevchenko, A., Svoboda, J., Tombesi, F., Trushkin, S., Tsybulev, P., Ursini, F., Weisskopf, M. C., Wu, K., Agudo, I., Antonelli, L. A., Bachetti, M., Baldini, L., Baumgartner, W. H., Bellazzini, R., Bongiorno, S. D., Bonino, R., Brez, A., Bucciantini, N., Castellano, S., Cavazzuti, E., Chen, C.-T., Ciprini, S., Del Monte, E., Di Gesu, L., Di Lalla, N., Donnarumma, I., Doroshenko, V., Ehlert, S. R., Enoto, T., Evangelista, Y., Fabiani, S., Ferrazzoli, R., Gunji, S., Hayashida, K., Heyl, J., Iwakiri, W., Jorstad, S. G., Karas, V., Kislat, F., Kitaguchi, T., Kolodziejczak, J. J., Latronico, L., Liodakis, I., Maldera, S., Manfreda, A., Marinucci, A., Marscher, A. P., Marshall, H. L., Massaro, F., Mitsuishi, I., Mizuno, T., Negro, M., Ng, C.-Y., O’Dell, S. L., Omodei, N., Oppedisano, C., Papitto, A., Pavlov, G. G., Peirson, A. L., Perri, M., Pesce-Rollins, M., Petrucci, P.-O., Possenti, A., Puccetti, S., Ramsey, B. D., Rankin, J., Roberts, O., Romani, R. W., Sgrò, C., Slane, P., Spandre, G., Swartz, D., Tamagawa, T., Tavecchio, F., Taverna, R., Tawara, Y., Tennant, A. F., Thomas, N. E., Trois, A., Turolla, R., Vink, J., Xie, F., and Zane, S. (2023).  Astronomical puzzle Cyg X-3 is a hidden Galactic ultraluminous X-ray source.  arXiv e-prints, page arXiv:2303.01174.
  * Viallet and Hameury, (2007) Viallet, M. and Hameury, J. M. (2007).  Hydrodynamic simulations of irradiated secondaries in dwarf novae.  A&A, 475(2):597–606.
  * Vincent et al., (2021) Vincent, F. H., Wielgus, M., Abramowicz, M. A., Gourgoulhon, E., Lasota, J. P., Paumard, T., and Perrin, G. (2021).  Geometric modeling of M87* as a Kerr black hole or a non-Kerr compact object.  A&A, 646:A37.
  * Vogt, (1983) Vogt, N. (1983).  VW Hydri revisited : conclusions on dwarf nova outburst models.  A&A, 118:95–101.
  * Weng et al., (2017) Weng, S.-S., Ge, M.-Y., Zhao, H.-H., Wang, W., Zhang, S.-N., Bian, W.-H., and Yuan, Q.-R. (2017).  Swift Observations of SMC X-3 during Its 2016-2017 Super-Eddington Outburst.  ApJ, 843(1):69.
  * Wiktorowicz et al., (2014) Wiktorowicz, G., Belczynski, K., and Maccarone, T. (2014).  Black hole X-ray transients: The formation puzzle.  In de Grijs, R., editor, Binary Systems, their Evolution and Environments, page 37.
  * Wood et al., (1986) Wood, J., Horne, K., Berriman, G., Wade, R., O’Donoghue, D., and Warner, B. (1986).  High-speed photometry of the dwarf nova Z Cha in quiescence.  MNRAS, 219:629–655.
  * Yang et al., (2023) Yang, Q., Green, P. J., MacLeod, C. L., Plotkin, R. M., Anderson, S. F., Bieryla, A., Civano, F., Eracleous, M., Graham, M., Ruan, J. J., Runnoe, J., and Zhao, X. (2023).  Probing the Origin of Changing-look Quasar Transitions with Chandra.  arXiv e-prints, page arXiv:2303.06733.
  * Yungelson and Lasota, (2008) Yungelson, L. R. and Lasota, J. P. (2008).  Evolutionary models of short-period soft X-ray transients: comparison with observations.  A&A, 488(1):257–265.
  * Yungelson et al., (2006) Yungelson, L. R., Lasota, J. P., Nelemans, G., Dubus, G., van den Heuvel, E. P. J., Dewi, J., and Portegies Zwart, S. (2006).  The origin and fate of short-period low-mass black-hole binaries.  A&A, 454(2):559–569.
  * Zdziarski et al., (2022) Zdziarski, A. A., You, B., and Szanecki, M. (2022).  Corrections to Estimated Accretion Disk Size due to Color Correction, Disk Truncation, and Disk Wind.  ApJLett, 939(1):L2.
  * Zel’dovich and Novikov, (1965) Zel’dovich, Y. B. and Novikov, I. D. (1965).  Mass of Quasi-Stellar Objects.  Soviet Physics Doklady, 9:834.
  * Zemko et al., (2014) Zemko, P., Orio, M., Mukai, K., and Shugarov, S. (2014).  X-ray observations of VY Scl-type nova-like binaries in the high and low state.  MNRAS, 445(1):869–880.