# Relativistic Linear Response in Quantum-Electrodynamical Density Functional
Theory

Lukas Konecny<EMAIL_ADDRESS>Hylleraas Centre for Quantum Molecular
Sciences, Department of Chemistry, UiT The Arctic University of Norway, N-9037
Tromsø, Norway  Max Planck Institute for the Structure and Dynamics of Matter,
Center for Free Electron Laser Science, Luruper Chaussee 149, 22761 Hamburg,
Germany  Valeriia P. Kosheleva<EMAIL_ADDRESS>Max Planck
Institute for the Structure and Dynamics of Matter, Center for Free Electron
Laser Science, Luruper Chaussee 149, 22761 Hamburg, Germany  Heiko Appel
<EMAIL_ADDRESS>Max Planck Institute for the Structure and Dynamics of
Matter, Center for Free Electron Laser Science, Luruper Chaussee 149, 22761
Hamburg, Germany  Michael Ruggenthaler<EMAIL_ADDRESS>Max
Planck Institute for the Structure and Dynamics of Matter, Center for Free
Electron Laser Science, Luruper Chaussee 149, 22761 Hamburg, Germany  Angel
Rubio<EMAIL_ADDRESS>Max Planck Institute for the Structure and
Dynamics of Matter, Center for Free Electron Laser Science, Luruper Chaussee
149, 22761 Hamburg, Germany  Center for Computational Quantum Physics (CCQ),
The Flatiron Institute, 162 Fifth Avenue, New York, New York 10010, USA

###### Abstract

We present the theoretical derivation and numerical implementation of the
linear response equations for relativistic quantum electrodynamical density
functional theory (QEDFT). In contrast to previous works based on the
Pauli–Fierz Hamiltonian, our approach describes electrons interacting with
photonic cavity modes at the four-component Dirac-Kohn-Sham level, derived
from fully relativistic QED through a series of established approximations.
Moreover, we show that a new type of spin–orbit-like (SO) cavity-mediated
interaction appears under the relativistic description of the coupling of
matter with quantized cavity modes. Benchmark calculations performed for atoms
of group 12 elements (Zn, Cd, Hg) demonstrate how a relativistic treatment
enables the description of exciton polaritons which arise from the
hybridization of formally forbidden singlet–triplet transitions with cavity
modes. For atoms in cavities tuned on resonance with a singlet–triplet
transition we discover a significant interplay between SO effects and coupling
to an off-resonant intense singlet-singlet transition. This dynamic
relationship highlights the crucial role of ab initio approaches in
understanding cavity quantum electrodynamics. Finally, using the mercury
porphyrin complex as an example, we show that relativistic linear response
QEDFT provides computationally feasible first-principles calculations of
polaritonic states in large heavy element-containing molecules of chemical
interest.

##  I Introduction

The strength of light–matter coupling, normally weak when a material interacts
with vacuum modes of free space, can be greatly enhanced by placing the
material into an optical cavity. If a cavity mode is in resonance with an
excitation of the material, new hybrid light–matter states called polaritons
emerge. [1, 2, 3, 4] Depending on the nature of the excitation, different
types of polaritonic states exist such exciton–polaritons and vibrational
polaritons resulting from the coupling of light to electronic and vibrational
transitions, respectively. As a result a new field of cavity quantum
electrodynamical (QED) materials [5, 1, 2, 6, 7] has emerged bringing together
various platforms for manipulating and engineering quantum materials with
electromagnetic fields. It encompasses a wide spectrum of different fields
including quantum optics [8], polaritonic chemistry [1, 9, 2, 10, 11, 12, 13,
7], the generation of light-induced states of matter through classical fields
[14, 15] or quantum fields emanating from a cavity [16, 17, 18, 19]. This
rapidly growing field has attracted significant attention from both theorists
and experimentalists, leading to the exploration of various effects under
strong light–matter coupling. Strong light–matter interaction in a cavity
enables the modification of chemical reaction landscapes [20, 21, 22, 23, 24,
25] and ground states of matter [26, 27, 28, 29], enhances charge and energy
transfer [30, 31, 32, 33, 34, 35], and allows for selective manipulation of
electronic states [35], the modification of electron–phonon coupling and
superconductivity [36, 37, 38, 39], controlling of excitons [40, 41, 42] while
achieving exciton-polariton condensation [43, 44]. In recent studies, Quantum
Hall systems, whether operating in the integer regime [45, 46, 47, 48, 49] or
the fractional regime [50, 51], have exhibited remarkable phenomena when
confined within optical cavities. These systems have shown what is known as
ultra-strong coupling with the light field, resulting in significant changes
in their transport properties [52]. Furthermore, there is a growing interest
in exploring the implications of coupling with chiral electromagnetic fields,
and this topic is currently being actively investigated [53, 54, 55, 56]. This
regime of the strong coupling between matter and light can be achieved due to
collective effects, when a large number of molecules is placed in a cavity and
coherently interact with its modes, or for micro- and nano-cavities that
confine light to very small length scales [1, 57, 58, 3, 7].

To describe the interaction of the material with vacuum fluctuations i.e.
cavity modes, the quantum description of light is important making QED the
theory of choice. The standard practical approach is provided by models of
cavity QED, which describe a system of photons strongly coupled to a material
approximated as a few-level system, e.g. the Jaynes–Cummings [59] and the
Dicke [60] model. However, if we want to capture subtle and complex changes of
the material inside a cavity, a first principles description of the matter
subsystem becomes necessary. [2].

To bridge the worlds of quantum chemistry and quantum optics, new
computational methods that simultaneously allow ab initio description of
matter while treating transverse photons as dynamical variables, were
developed [2, 12, 61, 7]. The first such method was quantum electrodynamical
density functional theory (QEDFT) proposed for time-dependent [62, 63, 64] and
ground-state/coupled-vacuum properties [65, 66]. Besides non-interacting
auxiliary systems in the Kohn-Sham formulation of QEDFT also alternative,
explicitly correlated approaches were introduced [67]. The first computational
implementations of QEDFT were reported by Flick et al. for time-dependent
situations [5, 68], for ground-state properties [26] and finally as a linear
response theory of coupled matter–photon system and formulated on a real-space
grid. [69] It motivated subsequent implementations based on Gaussian orbitals
[70, 71] as well as a development of a real-space Sternheimer formalism [72].
Moreover, the work of Yang et al. [70] also considers various simplified cases
based on the neglect of the dipole self-energy terms, on the Tamm–Dancoff
approximation (TDA) for electrons [73] as well as on an analogous
approximation for photons, together with the different combinations of these
approximations. In addition to QEDFT, ab initio theories combining QED with
wave function based quantum chemical methods were developed starting with
quantum electrodynamical Hartree–Fock (QED-HF) theory [74, 75, 76], and
ranging to quantum electrodynamical coupled clusters (QED-CC) [75, 77] and
quantum electrodynamical configuration interaction singles (QED-CIS) [78], as
well as Green function-based approaches [79] and reduced density-matrix
theories [74, 76].

Most of these realizations of QEDFT were based on the non-relativistic
Pauli–Fierz Hamiltonian. However, there exists a plethora of phenomena called
relativistic effects that are not covered by the non-relativistic, Schrödinger
or Pauli equation-based treatment and are described correctly only within
relativistic theories based on the Dirac equation. [80] The most prominent
examples include the color of gold [81, 82], the liquid state of mercury [83],
and the voltage of the lead–acid battery [84]. Further manifestations of
relativistic effects in atomic, molecular, and solid state systems involve
bound-electron $g$ factor [85, 86], electron absorption and X-ray spectra [87,
88, 89, 90, 91], NMR [92, 93, 94, 95], pNMR [96, 97, 98, 99, 100] and EPR
spectroscopies [101, 102, 103, 104, 105, 100], bond lengths, [106] reaction
mechanisms [107], phosphorescence and decay pathways [108], properties of
superheavy elements [109, 110], band structures [111], band gap opening, [112,
113, 114] and the stability of structural phases [115]. Relativistic quantum
chemistry [116, 117] addresses the relativistic effects by describing
electrons at the level of the Dirac equation. The gold standard is represented
by the four-component (4c) Dirac–Coulomb Hamiltonian that combines the one-
electron Dirac equation with the instantaneous Coulomb interaction among
electrons and nuclei and electrons with each other. [118] This Hamiltonian
contains relativistic kinetic energy as well as one-electron spin–orbit
coupling and is sometimes extended to include Breit terms. [118, 119, 120,
121] Moreover, there is an ongoing research effort focused on including the
effects of relativistic QED in quantum chemistry to account for radiative
corrections [122, 123, 124, 125, 126, 127, 128, 129]. At the same time,
researchers are developing approximate two-component Hamiltonians that would
capture the most important relativistic effects at a lower computational cost.
[130, 131, 132, 133]

QEDFT with relativistic Dirac equation-based treatment of electronic structure
thus enables the correct description of molecules containing any element
across the periodic table when strongly coupled to cavity modes (under certain
stability conditions detailed in Sec. II.2). Moreover, the 4c description
includes SOC variationally thus allowing first-principles access to
singlet–triplet transitions that are forbidden according to non-relativistic
theories. While spin–orbit coupling can be added perturbatively, usually up to
the first order in $\frac{Z}{c^{2}}$ ($Z$ being the atomic number and $c$ the
speed of light), the fully relativistic treatment accounts for the higher
order effect that can become prominent for high-$Z$ systems. [134] In
addition, the linear response properties such as the coupling strengths
between the singlet and triplet states or the phosphorescence parameters are
correctly obtained from the first-order perturbation theory [135] instead of
requiring the more demanding quadratic response as they would with
perturbative SOC. [136] Another case for relativistic QEDFT is provided by
cavity engineering of the energy ratio between singlet–singlet and
singlet–triplet excited states. This was suggested to improve the yields of
singlet fission [137] and intersystem crossing [35, 138], particularly for
heavy element-containing systems where the rate of these processes is enhanced
by the increased spin–orbit coupling. [139, 140] Additionally, the inclusion
of SOC in QEDFT theory will enable to compute spin–orbit related phenomena in
complex nanostructures and solids. [141, 142]

At the same time, relativistic QEDFT calculations close the gap between the
stable and non-perturbative low-energy sector of QED described by the
Pauli–Fierz quantum field theory [143, 144, 7] and a fully second-quantized
QED formulation in various flavours [145, 146, 147]. Even without cavities and
their changes in the modes of the electromagentic field, there are many
fundamental questions still to be answered within QED as a general framework
to describe on the most fundamental level the interaction between charged
particles and light. [148]

In this work, we introduce a relativistic quantum-electrodynamical density
functional theory (QEDFT) for molecules that combines four-component
Dirac–Kohn–Sham treatment of electrons with quantized description of photonic
modes. We first present the theoretical foundation of the developed method
where we derive the relativistic Hamiltonian for the coupled light–matter
system in the long-wavelength approximation. We demonstrate the emergence of a
novel spin-orbit-like cavity-mediated interaction within the framework of
relativistic coupling between matter and quantized cavity modes. Then in Sec.
II.4 we formulate a linear response equation in the framework of QEDFT for its
excitation energies that allows the first-principles calculation of
polaritonically modified spectra and present its numerical implementation. As
an example, we calculate the absorption spectra of Zn, Cd, and Hg atoms
focusing on the singlet–triplet transitions in a cavity. Via careful analysis
of the absorption spectra we show that there is a significant interplay
between spin-orbit effects and off-resonant coupling to the intense singlet-
singlet transition. This interplay highlights the crucial role of the ab
initio approach in understanding cavity quantum electrodynamics. Finally, we
considered a large heavy element complex to demonstrate the applicability of
our method to systems of chemical interest.

##  II Theory

The SI system of units is used throughout the paper unless specified
otherwise. For additional notations and conventions see Appendix A.

###  II.1 Light-matter coupling in fully-relativistic limit: total Hamiltonian
of a system

We consider a system of electrons minimally coupled to photons in the presence
of external classical fields. The full QED Hamiltonian of such a system in
Coulomb gauge is formally given as [149, 63]

$\displaystyle\begin{split}\hat{H}_{\mathrm{QED}}(t)&=\int d{\bf
r}:\hat{\bar{\psi}}({\bf r})(-i\hbar
c\bm{\gamma}\cdot\bm{\nabla}+\alpha^{0}mc^{2})\hat{\psi}({\bf
r}):+\frac{1}{c}\int d{\bf r}:\hat{j}^{0}({\bf r}):a_{0}^{\mathrm{ext}}({\bf
r},t)\\\ &+\frac{1}{c^{2}}\int d{\bf r}d{\bf
r}^{\prime}\frac{j_{0}^{\mathrm{ext}}({\bf r},t):\hat{j}^{0}({\bf
r}^{\prime}):}{4\pi\varepsilon_{0}|{\bf r}-{\bf
r^{\prime}}|}+\frac{1}{2c^{2}}\int d{\bf r}d{\bf
r}^{\prime}\frac{:\hat{j}_{0}({\bf r})\hat{j}^{0}({\bf
r}^{\prime}):}{4\pi\varepsilon_{0}|{\bf r}-{\bf r^{\prime}}|}\\\
&-\frac{1}{c}\int d{\bf r}:\hat{\mathbf{j}}({\bf
r})(\mathbf{a}^{\mathrm{ext}}({\bf r},t)+\hat{\mathbf{A}}({\bf
r})):-\frac{1}{c}\int d\mathbf{r}\mathbf{j}^{\text{ ext }}({\bf
r},t)\hat{\mathbf{A}}({\bf r})\\\ &+\frac{\varepsilon_{0}}{2}\int d{\bf
r}:(\hat{\mathbf{E}}_{\perp}^{2}({\bf r})+c^{2}\hat{\mathbf{B}}^{2}({\bf
r})):.\end{split}$ (1)

Here $c$ is the speed of the light in free-space vacuum, $\hbar$ is the
reduced Planck constant, $\varepsilon_{0}$ is the vacuum permittivity, $m$ is
the bare mass of the electron, $\alpha^{\mu}=\gamma^{0}\gamma^{\mu}$ are the
Dirac matrices in the standard representation (see Appendix A). Here all the
operators are in the Schrödinger picture and $::$ denotes their normal
ordering.

The electron-positron field operator $\hat{\psi}({\bf r})$ is given by [150,
149]

$\displaystyle\begin{split}\hat{\psi}({\bf r})&=\sum_{\mu=\pm\frac{1}{2}}\int
d{\bf p}\left(\hat{a}_{{\bf p}\mu}\psi^{(+)}_{{\bf p}\mu}({\bf
r})+\hat{b}^{\dagger}_{{\bf p}\mu}\psi^{(-)}_{{\bf p}\mu}({\bf
r})\right),\end{split}$ (2)

where $\hat{a}_{{\bf p}\mu}$ is an annihilation operator of the electron with
momentum ${\bf p}=(p_{x},p_{y},p_{z})$ and helicity $\mu$, operator
$\hat{b}^{\dagger}_{{\bf p}\mu}$ is a creation operator for the positron with
the same momenta and helicity. Here $\psi^{(+)}_{{\bf p}\mu}({\bf r})$ and
$\psi^{(-)}_{{\bf p}\mu}({\bf r})$ are the solutions of the free Dirac
equation for the electron and positron, respectively. The electron-positron
field operator $\hat{\psi}({\bf r})$ formally obeys the usual Heisenberg
equation of motion (here we indicate with label $\mathrm{H}$ the Heisenberg
picture)

$\frac{d}{dt}\hat{\psi}_{\mathrm{H}}({\bf
r},t)=\frac{i}{\hbar}\left[\hat{H}_{\mathrm{QED,H}}(t),\hat{\psi}_{\mathrm{H}}({\bf
r},t)\right].$ (3)

The Dirac four-current operator $\hat{J}({\bf r})$ is constructed as
$\hat{J}({\bf r})=\left(\hat{j}_{0}({\bf r}),\hat{\mathbf{j}}({\bf r})\right)$
with the components having the following form

$\displaystyle\hat{j}_{0}({\bf r})$ $\displaystyle=ec:\hat{\bar{\psi}}({\bf
r})\gamma_{0}\hat{\psi}({\bf r}):,$ (4) $\displaystyle\hat{\mathbf{j}}({\bf
r})$ $\displaystyle=ec:\hat{\bar{\psi}}({\bf r})\bm{\gamma}\hat{\psi}({\bf
r}):,$ (5)

where $\hat{\bar{\psi}}({\bf r})\equiv\hat{\psi}^{\dagger}({\bf
r})\gamma^{0}$.
The quantized vector potential $\hat{\mathbf{A}}({\bf r})$ in Coulomb gauge is
given by [149] 111We note that in comparison with the usual definition of a
vector potential, we here have an extra multiplication by $c$ in order to have
the units of $A_{\mu}$ in Volts.

$\hat{\mathbf{A}}({\bf
r})=\sqrt{\frac{c^{2}\hbar}{\varepsilon_{0}(2\pi)^{3}}}\int\frac{\mathrm{d}\mathbf{k}}{\sqrt{2\omega_{\mathbf{k}}}}\sum_{\lambda=1}^{2}\left[\hat{a}(\mathbf{k},\lambda)\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}}\bm{\varepsilon}(\mathbf{k},\lambda)+\hat{a}^{\dagger}(\mathbf{k},\lambda)\mathrm{e}^{-\mathrm{i}\mathbf{k}\cdot\mathbf{r}}\bm{\varepsilon}^{*}(\mathbf{k},\lambda)\right],$
(6)

where $\hat{a}(\mathbf{k},\lambda)$ and
$\hat{a}^{\dagger}(\mathbf{k},\lambda)$ are the creation and annihilation
operators of a photon with momentum $\mathbf{k}$ and helicity $\lambda$,
respectively; $\omega_{\mathbf{k}}=c|\mathbf{k}|$ is the energy of the photon
and $\bm{\varepsilon}(\mathbf{k},\lambda)$ is the transverse dimensionless
polarization vector that obeys [149]

$\mathbf{k}\cdot\bm{\varepsilon}(\mathbf{k},\lambda)=\bm{\varepsilon}^{*}(\mathbf{k},1)\cdot\bm{\varepsilon}(\mathbf{k},2)=0.$
(7)

We note that in a Coulomb gauge the field $\hat{\mathbf{A}}({\bf r})$ is fully
transverse, i. e. $\bm{\nabla}\cdot\hat{\mathbf{A}}({\bf r})=0$.
The vector potential $\hat{\mathbf{A}}({\bf r})$ formally obeys the following
equation of motion

$\frac{d}{dt}\hat{\mathbf{A}}_{\mathrm{H}}({\bf
r},t)=\frac{i}{\hbar}\left[\hat{H}_{\mathrm{QED,H}}(t),\hat{\mathbf{A}}_{\mathrm{H}}({\bf
r},t)\right].$ (8)

The last term in Eq. (1) is the free photon Hamiltonian
$\hat{H}_{\text{Ph,free }}$ which creates and annihilates only transverse
photons. Using Eq. (6) one can rewrite it as

$\hat{H}_{\text{Ph,free }}=\frac{\varepsilon_{0}}{2}\int d{\bf
r}:(\hat{\mathbf{E}}_{\perp}^{2}({\bf r})+c^{2}\hat{\mathbf{B}}^{2}({\bf
r})):=\sum_{\lambda=1}^{2}\int\mathrm{d}\mathbf{k}\hbar\omega_{\mathbf{k}}\hat{a}^{\dagger}(\mathbf{k},\lambda)\hat{a}(\mathbf{k},\lambda).$
(9)

Here $\hat{\mathbf{E}}_{\perp}({\bf
r})=\left[-\frac{\partial\hat{\mathbf{A}}_{\rm H}({\bf r},t)}{c\partial
t}\right]_{\rm S}$ with the subindex S denoting the Schrödinger picture, and
$\hat{\mathbf{B}}({\bf r})=\frac{1}{c}\bm{\nabla}\times\hat{\mathbf{A}}({\bf
r})$, are quantized transverse electric and magnetic fields, respectively
[149], that expressed via photonic creation and annihilation operators assume
the form

$\displaystyle\hat{\mathbf{E}}_{\perp}({\bf r})$
$\displaystyle=\sqrt{\frac{\hbar}{\varepsilon_{0}(2\pi)^{3}}}\int\frac{\mathrm{d}\mathbf{k}i\omega_{\mathbf{k}}}{\sqrt{2\omega_{\mathbf{k}}}}\sum_{\lambda=1}^{2}\left[\hat{a}(\mathbf{k},\lambda)\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}}\bm{\varepsilon}(\mathbf{k},\lambda)-\hat{a}^{\dagger}(\mathbf{k},\lambda)\mathrm{e}^{-\mathrm{i}\mathbf{k}\cdot\mathbf{r}}\bm{\varepsilon}^{*}(\mathbf{k},\lambda)\right],$
(10) $\displaystyle\hat{\mathbf{B}}({\bf r})$
$\displaystyle=\sqrt{\frac{\hbar}{\varepsilon_{0}(2\pi)^{3}}}\int\frac{\mathrm{d}\mathbf{k}}{\sqrt{2\omega_{\mathbf{k}}}}\sum_{\lambda=1}^{2}\left[\hat{a}(\mathbf{k},\lambda)\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}}[\mathbf{k}\times\bm{\varepsilon}(\mathbf{k},\lambda)]+\hat{a}^{\dagger}(\mathbf{k},\lambda)\mathrm{e}^{-\mathrm{i}\mathbf{k}\cdot\mathbf{r}}[\mathbf{k}\times\bm{\varepsilon}^{*}(\mathbf{k},\lambda)]\right].$
(11)

In our system, we also introduce external (with respect to the system’s
electrons plus photons) classical magnetic vector potential
$\mathbf{a}^{\mathrm{ext}}({\bf r},t)$ and scalar potential
$a^{\mathrm{ext}}_{0}({\bf r},t)$ as well as classical external four-currents
$J^{\mathrm{ext}}({\bf r},t)=\left(j_{0}^{\mathrm{ext}}({\bf
r},t)=c\rho^{\mathrm{ext}}({\bf r},t),\mathbf{j}^{\mathrm{ext}}({\bf
r},t)\right)$ with $\rho^{\mathrm{ext}}({\bf r},t)$ being a charge density.
Using Eq. (6) one can rewrite the energy due to coupling to a classical
external charge current $\mathbf{j}^{\mathrm{ext}}({\bf r},t)$ as [65]

$\frac{1}{c}\int d\mathbf{r}\mathbf{j}^{\text{ ext }}({\bf
r},t)\hat{\mathbf{A}}({\bf r})=\int
d\mathbf{k}\hbar\omega_{\mathbf{k}}\left(\hat{a}(\mathbf{k},\lambda)[j^{\text{
ext }}]^{*}(\mathbf{k},\lambda,t)+j^{\text{ ext
}}(\mathbf{k},\lambda,t)\hat{a}^{\dagger}(\mathbf{k},\lambda)\right),$ (12)

with the expansion coefficients $j^{\text{ ext
}}(\mathbf{k},\lambda,t)=[j^{\text{ ext
}}(-\mathbf{k},\lambda,t)]^{*}=\left(\frac{1}{2}\omega_{\mathbf{k}}^{3}\varepsilon_{0}\hbar(2\pi)^{3}\right)^{-1/2}\int
d\mathbf{r}\bm{\varepsilon}^{*}(\mathbf{k},\lambda)\cdot\mathbf{j}^{\text{ ext
}}(\mathbf{r},t)\mathrm{e}^{-\mathrm{i}\mathbf{k}\cdot\mathbf{r}}$. We note,
that with these expansion coefficients, one can recover the transverse part of
the external current since $\bm{\varepsilon}(\mathbf{k},\lambda)$ is a
transverse vector

$\mathbf{j}^{\text{ ext }}_{\perp}({\bf
r},t)=\sqrt{\frac{\hbar\varepsilon_{0}}{(2\pi)^{3}}}\int\frac{\mathrm{d}\mathbf{k}\omega_{\mathbf{k}}^{2}}{\sqrt{2\omega_{\mathbf{k}}}}\sum_{\lambda=1}^{2}\left[j^{\text{
ext
}}(\mathbf{k},\lambda,t)\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}}\bm{\varepsilon}(\mathbf{k},\lambda)+[j^{\text{
ext
}}(\mathbf{k},\lambda,t)]^{*}\mathrm{e}^{-\mathrm{i}\mathbf{k}\cdot\mathbf{r}}\bm{\varepsilon}^{*}(\mathbf{k},\lambda)\right].$
(13)

###  II.2 Approximations

We note that so far, we have chosen implicitly circularly polarised vectors
$\bm{\varepsilon}(\mathbf{k},\lambda)$. However, for the quantization of the
electromagnetic field, one can choose any pair of polarization vectors obeying
(7), for example, linearly polarised ones. We also note that circular and
linear representations are connected by a canonical transformation that
preserves the form of the equations. Therefore, without loss of generality, we
restrict ourselves to the case of linear polarisation, i.e., real transverse
polarization vectors. Then the coupling term from Eq. (12) reads as

$\frac{1}{c}\int d\mathbf{r}\mathbf{j}^{\text{ ext }}({\bf
r},t)\hat{\mathbf{A}}({\bf r})=\int
d\mathbf{k}\hbar\omega_{\mathbf{k}}j^{\text{ ext
}}(\mathbf{k},\lambda,t)\left(\hat{a}(\mathbf{k},\lambda)+\hat{a}^{\dagger}(\mathbf{k},\lambda)\right),$
(14)

with the real expansion coefficients

$j^{\text{ ext }}(\mathbf{k},\lambda,t)=[j^{\text{ ext
}}(\mathbf{k},\lambda,t)]^{*}=\left(\frac{1}{2}\omega_{\mathbf{k}}^{3}\varepsilon_{0}\hbar(2\pi)^{3}\right)^{-1/2}\int
d\mathbf{r}\bm{\varepsilon}(\mathbf{k},\lambda)\cdot\mathbf{j}^{\text{ ext
}}(\mathbf{r},t)\exp^{-\mathrm{i}\mathbf{k}\cdot\mathbf{r}}.$ (15)

Next, we would like to work with a stable vacuum, i.e. the energy gap between
the positive and negative energy solutions of the Dirac equation is nonzero.
If the intensity of the external fields is small compared to the Schwinger
limit, and if one considers the long-wavelength approximation then no pair
production occurs [151, 152, 153]. The energy gap then substitutes the
boundedness-from-below (minimal-energy state) of non-relativistic theories.
Since no pair creation is present, in addition to charge conservation,
particle-number conservation is satisfied, allowing us to proceed to a first-
quantized formulation for electrons. Specifically due to the long-wavelength
approximation, Eq. (6) can be written as

$\hat{\mathbf{A}}=\sqrt{\frac{\hbar
c^{2}}{\varepsilon_{0}(2\pi)^{3}}}\int\frac{\mathrm{d}\mathbf{k}}{\sqrt{2\omega_{\mathbf{k}}}}\sum_{\lambda=1}^{2}\bm{\varepsilon}(\mathbf{k},\lambda)\left[\hat{a}(\mathbf{k},\lambda)+\hat{a}^{\dagger}(\mathbf{k},\lambda)\right],$
(16)

and then the Hamiltonian of Eq. (1) takes a form

$\displaystyle\begin{split}\hat{H}_{\mathrm{QED}}^{\mathrm{Dip}}(t)&\approxeq\sum_{l=1}^{N}[-i\hbar
c\bm{\alpha}_{l}\cdot\bm{\nabla}_{l}+\beta_{l}mc^{2}+ea_{0}^{\mathrm{ext}}({\bf
r}_{l},t)]+\frac{e}{c}\sum_{l=1}^{N}\int d{\bf
r}\frac{j_{0}^{\mathrm{ext}}({\bf r},t)}{4\pi\varepsilon_{0}|{\bf r}-{\bf
r}_{l}|}\\\ &+\frac{1}{2}\sum_{l\neq
m}^{N}\frac{e^{2}}{4\pi\varepsilon_{0}\left|{\bf r}_{l}-{\bf
r}_{m}\right|}-e\sum_{l=1}^{N}(\bm{\alpha}_{l}\cdot\mathbf{a}^{\mathrm{ext}}(t)+\bm{\alpha}_{l}\cdot\hat{\mathbf{A}})\\\
&+\sum_{\lambda=1}^{2}\int\mathrm{d}\mathbf{k}\hbar\omega_{\mathbf{k}}[\hat{a}^{\dagger}(\mathbf{k},\lambda)\hat{a}(\mathbf{k},\lambda)-j^{\text{
ext
}}(\mathbf{k},\lambda,t)\left(\hat{a}(\mathbf{k},\lambda)+\hat{a}^{\dagger}(\mathbf{k},\lambda)\right)],\end{split}$
(17)

where $N$ is the number of electrons in the system. Assume we are dealing with
a cavity, in which case the free-space mode expansion (16) can be converted
into the corresponding expansion for the photonic structure e. g. Fabry–Pérot
cavity [154]. Since we are working in a dipole approximation, i.e. no momentum
is transferred in a cavity, the actual spatial mode structure and the momentum
matching are no longer important. As a result, one can modify free-space
frequencies, coupling strength, and polarizations to fit a given cavity
structure without breaking fundamental symmetries [7]

$\hat{\mathbf{A}}=\sqrt{\hbar
c^{2}}\sum_{\alpha}^{M_{p}}g_{\alpha}\bm{\varepsilon}_{\alpha}\hat{Q}_{\alpha}.$
(18)

Here we have assumed a discretized continuum of $M_{p}$ modes labeled by
$\alpha$, with each $\alpha$ corresponding to a different frequency
$\omega_{\alpha}$, coupling strength $g_{\alpha}$, and polarization
$\bm{\varepsilon}_{\alpha}$. The modes are given in terms of generalized
coordinates $\hat{Q}_{\alpha}$ such as

$\displaystyle\hat{a}_{\alpha}=\sqrt{\frac{\omega_{\alpha}}{2}}\left(\hat{Q}_{\alpha}+\frac{1}{\omega_{\alpha}}\frac{\partial}{\partial\hat{Q}_{\alpha}}\right),$
(19)
$\displaystyle\hat{a}_{\alpha}^{\dagger}=\sqrt{\frac{\omega_{\alpha}}{2}}\left(\hat{Q}_{\alpha}-\frac{1}{\omega_{\alpha}}\frac{\partial}{\partial\hat{Q}_{\alpha}}\right),$
(20)
$\displaystyle\left[\hat{Q}_{\alpha},-i\frac{\partial}{\partial\hat{Q}_{\alpha^{\prime}}}\right]=i\delta_{\alpha\alpha^{\prime}},$
(21)
$\displaystyle\left[\hat{a}_{\alpha},\hat{a}_{\alpha^{\prime}}^{\dagger}\right]=\delta_{\alpha\alpha^{\prime}},$
(22)

where $-i\frac{\partial}{\partial\hat{Q}_{\alpha}}$ is conjugate momenta of a
generalized coordinate, $\hat{a}_{\alpha}$ and $\hat{a}_{\alpha}^{\dagger}$
are annihilation and creation of the photon in a mode $\alpha$, respectively.
We note that in the case of free space with a quantization volume $l^{3}$, the
quantities from Eq. (18) correspond to
$\mathbf{k}_{\mathbf{m}}=2\pi\mathbf{m}/l$,
$\alpha\equiv\left(\mathbf{k}_{\mathbf{m}},\lambda\right),\omega_{\alpha}=c\left|\mathbf{k}_{\mathbf{m}}\right|$
and $g_{\alpha}=\sqrt{1/\varepsilon_{0}l^{3}}$, where
$\mathbf{m}\in\mathbb{Z}^{3}$. In the case of general photonic structures the
respective quantities can be determined from e.g. macroscopic QED [155].

Figure 1: Schematic representation of a Fabry–Pérot cavity. We assume that
molecules of interest are localized around a cavity’s center, far away from
the cavity mirrors. Then one can approximate the main frequency due to the
distance between mirrors $L$ by $\omega_{n}=c\frac{\pi}{L}|\bf{n}|$ with
$\bf{n}$ parallel to L. The influence of the continuum free-space modes
(perpendicular to L) is incorporated in the observed mass of the particles.
The coupling strength $g_{\alpha\equiv(n,\lambda)}$ for the two independent
polarization directions $\lambda$ of the polarization vector
$\bf{\varepsilon}_{\alpha\equiv(n,\lambda)}$ increases proportionally to the
$\sqrt{\frac{1}{V}}$, where $V$ is an effective volume of a cavity.

Using Eqs. (18)-(22) one can rewrite Eq. (17) as

$\displaystyle\begin{split}\hat{H}_{\mathrm{QED,cavity}}^{\mathrm{Dip}}(t)&\approxeq\sum_{l=1}^{N}[-i\hbar
c\bm{\alpha}_{l}\cdot\bm{\nabla}_{l}+\beta_{l}mc^{2}+ea_{0}^{\mathrm{ext}}({\bf
r}_{l},t)]\\\ &+\frac{1}{2}\sum_{l\neq
m}^{N}\frac{e^{2}}{4\pi\varepsilon_{0}\left|{\bf r}_{l}-{\bf
r}_{m}\right|}-e\sum_{l=1}^{N}(\bm{\alpha}_{l}\cdot\mathbf{a}^{\mathrm{ext}}(t)+\bm{\alpha}_{l}\cdot\hat{\mathbf{A}})\\\
&+\sum_{\alpha=1}^{M_{p}}\left(-\frac{\hbar}{2}\frac{\partial^{2}}{\partial\hat{Q}_{\alpha}^{2}}+\frac{\hbar\omega_{\alpha}^{2}}{2}\hat{Q}_{\alpha}^{2}-\sqrt{2\hbar^{2}\omega_{\alpha}^{3}}j_{\alpha}^{\mathrm{ext}}(t)\hat{Q}_{\alpha}\right),\end{split}$
(23)

where

$j_{\alpha}^{\mathrm{ext}}(t)=[j_{\alpha}^{\mathrm{ext}}(t)]^{*}=\left(\frac{1}{2}\omega_{\mathbf{k}}^{3}\hbar\right)^{-1/2}\frac{1}{\varepsilon_{0}g_{\alpha}}\bm{\varepsilon}_{\alpha}\cdot\mathbf{j}^{\text{
ext }}(t).$ (24)

Note, that we probe our system with an external perturbation
$\mathbf{a}^{\mathrm{ext}}({\bf r}_{l},t)$ that does not generate longitudinal
currents, i.e. $j_{0}^{\mathrm{ext}}({\bf r},t)=0$.
In the dipole approximation, it is possible to eliminate the term
$e\sum_{l=1}^{N}(\bm{\alpha}_{l}\cdot\mathbf{a}^{\mathrm{ext}}(t)+\bm{\alpha}_{l}\cdot\hat{\mathbf{A}})$
from Eq. (23) by performing the unitary length-gauge transformation which is
defined as [154, 64, 156]

$\hat{H}_{\mathrm{L}}=\hat{U}^{\dagger}\hat{H}\hat{U}-i\hbar\hat{U}^{\dagger}(\partial_{t}\hat{U}),\quad\hat{U}=e^{\frac{\mathrm{i}}{\hbar
c}\left[\hat{\mathbf{A}}+\mathbf{a}^{\mathrm{ext}}(t)\right]\cdot\bm{\mu}},$
(25)

where $\bm{\mu}=\sum_{i=1}^{N}e\mathbf{r}_{i}$ is the total dipole operator of
the electrons in a system, and index $\mathrm{L}$ indicates the length gauge.
Using Eq. (25) one can rewrite the Hamiltonian in length gauge as

$\displaystyle\begin{split}\hat{U}^{\dagger}\hat{H}_{\mathrm{QED,cavity}}^{\mathrm{Dip}}(t)\hat{U}=\hat{H}_{\mathrm{L}}(t)&=\sum_{l=1}^{N}[-i\hbar
c\bm{\alpha}_{l}\cdot\bm{\nabla}_{l}+\beta_{l}mc^{2}+ea_{0}^{\mathrm{ext}}({\bf
r}_{l},t)]\\\ &+\frac{1}{2}\sum_{l\neq
m}^{N}\frac{e^{2}}{4\pi\varepsilon_{0}\left|{\bf r}_{l}-{\bf
r}_{m}\right|}-\sum_{l=1}^{N}v^{\text{ext}}({\bf r}_{l},t)\\\
&+\sum_{\alpha=1}^{M_{p}}\left(\frac{\hbar}{2}\left(\omega_{\alpha}\hat{q}_{\alpha}-\frac{g_{\alpha}\bm{\varepsilon}_{\alpha}\cdot\bm{\mu}}{\sqrt{\hbar}}\right)^{2}-\frac{\hbar}{2}\frac{\partial^{2}}{\partial\hat{q}_{\alpha}^{2}}+i\sqrt{2\hbar^{2}\omega_{\alpha}}j_{\alpha}^{\mathrm{ext}}(t)\frac{\partial}{\partial\hat{q}_{\alpha}}\right),\end{split}$
(26)

where $v^{\text{ext}}({\bf r}_{l},t)=e\mathbf{E}^{\mathrm{ext}}(t)\cdot{\bf
r}_{l}$ is the external potential generated by the external electric field
$\mathbf{E}^{\mathrm{ext}}(t)=-\tfrac{\partial}{c\partial
t}\mathbf{a}^{\mathrm{ext}}(t)$. Here we have performed the additional
canonical transformations
$-i\frac{\partial}{\partial_{\hat{Q}_{\alpha}}}\mapsto-\omega_{\alpha}\hat{q}_{\alpha},\hat{Q}_{\alpha}\mapsto-i\frac{\partial}{\omega_{\alpha}\partial_{\hat{q}_{\alpha}}}$,
which leave the commutation relations unchanged, i. e.
$\left[\hat{Q}_{\alpha},-i\frac{\partial}{\partial\hat{Q}_{\alpha^{\prime}}}\right]=\left[-i\frac{\partial}{\omega_{\alpha}\partial_{\hat{q}_{\alpha}}},-\omega_{\alpha}\hat{q}_{\alpha}\right]=i\delta_{\alpha\alpha^{\prime}}$.
We want to note that $\hat{q}_{\alpha}$ does not correspond to a pure photonic
quantity when we couple the modes to the electronic system. It is connected to
the auxiliary displacement field of the macroscopic Maxwell equations and is
usually referred to as a displacement coordinate [154, 156].
We are interested in studying the electronic structure properties of
molecules/atoms coupled to photon modes. Then $ea_{0}^{\mathrm{ext}}({\bf
r}_{l},t)\rightarrow v_{\mathrm{nucl}}({\bf
r}_{l})=-\frac{1}{4\pi\varepsilon_{0}}\sum_{j=1}^{N_{\mathrm{n}}}\frac{Z_{j}e^{2}}{\left|\mathbf{r}_{l}-\mathbf{R}_{j}\right|}$
is a Coulomb potential describing the interaction of $N_{\mathrm{n}}$ nuclei
with the bound electrons. Yet we perform another time-dependent unitary
transformation
$\hat{U}(t)=e^{\left[\frac{\mathrm{i}}{\hbar}\sum_{\alpha=1}^{M_{p}}\hat{q}_{\alpha}\sqrt{2\hbar^{2}\omega_{\alpha}}j_{\alpha}^{\mathrm{ext}}(t)\right]}$
aiming to eliminate the last term in Eq. (26). Then the Hamiltonian (26) takes
the final form

$\displaystyle\begin{split}\hat{H}_{\mathrm{L}}(t)&=\sum_{l=1}^{N}[-i\hbar
c\bm{\alpha}_{l}\cdot\bm{\nabla}_{l}+\beta_{l}mc^{2}+v_{\mathrm{nucl}}({\bf
r}_{l})]\\\ &+\frac{1}{2}\sum_{m\neq
l}^{N}\sum_{l=1}^{N}\frac{e^{2}}{4\pi\varepsilon_{0}\left|{\bf r}_{l}-{\bf
r}_{m}\right|}-\sum_{l=1}^{N}v^{\text{ext}}({\bf r}_{l},t)\\\
&+\sum_{\alpha=1}^{M_{p}}\left(\frac{\hbar}{2}\left(\omega_{\alpha}\hat{q}_{\alpha}-\frac{g_{\alpha}\bm{\varepsilon}_{\alpha}\cdot\bm{\mu}}{\sqrt{\hbar}}\right)^{2}-\frac{\hbar}{2}\frac{\partial^{2}}{\partial\hat{q}_{\alpha}^{2}}+\sqrt{2\hbar^{2}\omega_{\alpha}}\left[\partial_{t}j_{\alpha}^{\mathrm{ext}}(t)\right]\hat{q}_{\alpha}\right),\end{split}$
(27)

where we disregarded the physically irrelevant time-dependent energy/phase
shift $-\hbar\left[j_{\alpha}^{\mathrm{ext}}(t)\right]^{2}$.

A particular consequence of the Dirac Hamiltonian is the emergence of the
spin–orbit coupling (SOC). While SOC can be added ad hoc into a non-
relativistic theory based on the Schrödinger equation, it arises naturally
only in the Dirac equation due to its multi-component structure. The SOC can
be made explicit when expressing the Dirac equation in a two-component form,
or when performing expansion in different orders of the speed of light [116].
There are formally different ways of obtaining such an expansion [157, 158].
In our case, the most convenient form is from the relativistic energy
correction, such that the potential from the nuclei gives rise to one-electron
SOC contributions

$\Delta
E_{\mathrm{1e\,SOC}}=\frac{1}{4m^{2}c^{2}}\left\langle\left(\sum_{l=1}^{N}\bm{\sigma}_{l}\cdot\hat{\mathbf{p}}_{l}\right)\Psi^{L}\right|\sum_{m=1}^{N}v_{\mathrm{nucl}}({\bf
r}_{m})\left.\left(\sum_{l^{\prime}=1}^{N}\bm{\sigma}_{l^{\prime}}\cdot\hat{\mathbf{p}}_{l^{\prime}}\right)\Psi^{L}\right\rangle,$
(28)

while the electron–electron Coulomb interactions leads to two-electron SOC
contributions

$\Delta
E_{\mathrm{2e\,SOC}}=\frac{1}{4m^{2}c^{2}}\left\langle\left(\sum_{l=1}^{N}\bm{\sigma}_{l}\cdot\hat{\mathbf{p}}_{l}\right)\Psi^{L}\right|\frac{1}{2}\sum_{m^{\prime}\neq
m}^{N}\sum_{m=1}^{N}\frac{e^{2}}{4\pi\varepsilon_{0}\left|{\bf r}_{m}-{\bf
r}_{m^{\prime}}\right|}\left.\left(\sum_{l^{\prime}=1}^{N}\bm{\sigma}_{l^{\prime}}\cdot\hat{\mathbf{p}}_{l^{\prime}}\right)\Psi^{L}\right\rangle,$
(29)

where $\hat{\mathbf{p}}=-i\hbar\bm{\nabla}$ and $\Psi^{L}$ is the respective
“large” component of the Dirac wavefunction. In the presence of SOC the spin
operator no longer commutes with the Hamiltonian meaning that the eigenstates
of the Hamiltonian are not eigenstates of the spin operator, referred to as
spin no longer being a good quantum number. This results for example in non-
zero transition amplitudes between states of different spin symmetry.
Following the same strategy as done for the standard terms of Eqs. (28) and
(29) also for the new cavity-coupling term, we obtain cavity-mediated SOC
contributions from the relativistic energy correction as

$\displaystyle\begin{split}\Delta
E&=\frac{1}{4m^{2}c^{2}}\left\langle\left(\sum_{l=1}^{N}\bm{\sigma}_{l}\cdot\hat{\mathbf{p}}_{l}\right)\Psi^{L}\left|\sum_{\alpha=1}^{M_{p}}\frac{\hbar}{2}\left(-\frac{2\omega_{\alpha}\hat{q}_{\alpha}g_{\alpha}\bm{\varepsilon}_{\alpha}\cdot\bm{\mu}}{\sqrt{\hbar}}+\frac{g_{\alpha}^{2}(\bm{\varepsilon}_{\alpha}\cdot\bm{\mu})^{2}}{\hbar}\right)\left(\sum_{l^{\prime}=1}^{N}\bm{\sigma}_{l^{\prime}}\cdot\hat{\mathbf{p}}_{l^{\prime}}\right)\Psi^{L}\right.\right\rangle\\\
&=-\left\langle\Psi^{L}\left|\sum_{l=1}^{N}\sum_{\alpha=1}^{M_{p}}\frac{\hat{q}_{\alpha}\sqrt{\hbar^{3}}\omega_{\alpha}g_{\alpha}e}{8m^{2}c^{2}}\left(\bm{\varepsilon}_{\alpha}\cdot\left(\hat{\mathbf{p}}_{l}\times\bm{\sigma}_{l}\right)+\left(\hat{\mathbf{p}}_{l}\times\bm{\sigma}_{l}\right)\cdot\bm{\varepsilon}_{\alpha}\right)\right|\Psi^{L}\right\rangle+\\\
&+\left\langle\Psi^{L}\left|\sum_{l=1}^{N}\sum_{k=1}^{N}\sum_{\alpha=1}^{M_{p}}\frac{g_{\alpha}^{2}\hbar
e^{2}}{16m^{2}c^{2}}\left(\bm{\varepsilon}_{\alpha}\cdot\mathbf{r}_{k}\right)\left[\bm{\varepsilon}_{\alpha}\cdot\left(\hat{\mathbf{p}}_{l}\times\bm{\sigma}_{l}\right)+\left(\hat{\mathbf{p}}_{l}\times\bm{\sigma}_{l}\right)\cdot\bm{\varepsilon}_{\alpha}\right]\right|\Psi^{L}\right\rangle+\ldots\end{split}$
(30)

where the ellipsis stands for other terms that contain a correction to kinetic
energy and a Darwin-like term. Among others, we find a cavity-induced SOC term
of the form

$\displaystyle-\frac{e\hbar}{8m^{2}c^{2}}\sum_{l=1}^{N}\left(\hat{\mathbf{p}}_{l}\cdot(\bm{\sigma}_{l}\times\hat{\mathbf{E}})+(\bm{\sigma}_{l}\times\hat{\mathbf{E}})\cdot\hat{\mathbf{p}}_{l}\right),$
(31)

and the electric field operator of the cavity in the length gauge is given by

$\displaystyle\hat{\mathbf{E}}=\sum_{\alpha=1}^{M_{p}}\left(\sqrt{\hbar}g_{\alpha}\omega_{\alpha}\hat{q}_{\alpha}-\frac{g_{\alpha}^{2}}{2}\bm{\mu}\cdot\bm{\varepsilon}_{\alpha}\right)\bm{\varepsilon}_{\alpha}.$
(32)

Note that the electric dipole moment operator $\bm{\mu}$ contains a sum over
all electrons making the second term a collective effect. Here the the first
part of the electric field is proportional to the coupling strength
$g_{\alpha}$ and usually dominates for single-molecule situations in the weak
to strong coupling regime [2, 3]. However, for collective-coupling situation,
where it is the response of the total ensemble of molecules, the second term
proportional to $g_{\alpha}^{2}$ can become dominant [159, 160]. One can
therefore imagine that even a collectively-coupled ensemble of molecules could
show enhanced SOC effects. As a final remark, in a similar way to Eq. (31) the
term $v^{\text{ext}}({\bf r}_{l},t)$ in Eq. (27) can transiently lead to a
light-induced SOC due to external classical fields.

###  II.3 Kohn–Sham QEDFT equations of motion

After specifying the Hamiltonian one can find the time evolution of a system
from a given initial state
$\Psi_{0}\left(x_{1}\ldots,x_{N};q_{1},\ldots,q_{M_{p}}\right)$ by solving the
coupled electron–photon equation of motion

$i\hbar\partial_{t}\Psi\left(x_{1}\ldots,x_{N};q_{1},\ldots,q_{M_{p}};t\right)=\hat{H}_{\text{L}}(t)\Psi\left(x_{1}\ldots,x_{N};q_{1},\ldots,q_{M_{p}};t\right),$
(33)

where $x_{i}\equiv\left(\mathbf{r}_{i},\tau_{i}\right)$ denotes spatial
coordinate $\mathbf{r}_{i}$ and the four spin-components $\tau_{i}$ of $i$-th
Dirac electron. While solutions of Eq. (33) contain all the information about
the given system, it becomes numerically intractable for realistic many-
electron systems and therefore approximate methods have to be employed. In
this work, our method of choice is density functional theory in the form of
QEDFT and in notation we follow the formalism outlined in Ref. 69. The
electron–photon system is described by the electronic density
$n\left(\mathbf{r},t\right)=N\sum_{\tau,\tau_{2},...,\tau_{N}}\int|\Psi(\mathbf{r},\tau,x_{2},...,x_{N})|^{2}\mathrm{d}\mathbf{r}_{2}...\mathrm{d}\mathbf{r}_{N}$
and the expectation value of photonic coordinates $q_{\alpha}\left(t\right)$
as basic variables [64, 72, 69, 63] which uniquely determine its state. We
work in the Kohn-Sham (KS) formalism, where we consider an auxiliary system of
independent particles that reproduces the density of the real system and is
described by molecular orbitals $\varphi_{j}(x,t)$, with index $j$ enumerating
the $N_{\mathrm{occ}}$ occupied orbitals. The electron density is then
calculated as [161]

$n(\mathbf{r},t)=\sum_{j=1}^{N_{\mathrm{occ}}}\sum_{\tau=1}^{4}\left|\varphi_{j}(\mathbf{r},\tau,t)\right|^{2}.$
(34)

Then the equations of motion for the basic variables in terms of the auxiliary
system read as [69, 63]

$\displaystyle i\hbar\partial_{t}\varphi_{j}(x,t)$
$\displaystyle=\hat{H}^{\text{KS}}(t)\varphi_{j}(x,t),$ (35)
$\displaystyle\left(\frac{\partial^{2}}{\partial
t^{2}}+\omega_{\alpha}^{2}\right)q_{\alpha}(t)$
$\displaystyle=-\sqrt{2}\partial_{t}j_{\alpha,\mathrm{KS}}\left(\left[n,q_{\alpha},t\right]\right),$
(36)

with the single-particle KS Hamiltonian $\hat{H}^{\text{KS}}(t)$

$\hat{H}^{\text{KS}}(t)=-i\hbar c\bm{\alpha}\cdot\bm{\nabla}+\beta
mc^{2}+v_{\mathrm{KS}}\left(\left[v,n,q\right];{\bf r},t\right),$ (37)

and the Kohn–Sham current
$j_{\alpha,\mathrm{KS}}\left(\left[n,q_{\alpha},t\right]\right)$ defined below
in Eq. (44). The Kohn–Sham potential
$v_{\mathrm{KS}}\left(\left[v,n,q\right];{\bf r},t\right)$ in Eq. (37) is
defined as

$v_{\mathrm{KS}}\left(\left[v,n,q\right];{\bf r},t\right)=v({\bf
r},t)+v_{\mathrm{Mxc}}\left(\left[n,q\right];{\bf r},t\right).$ (38)

Here the external potential $v({\bf r},t)$ according to Eq. (27) is given by

$v({\bf r},t)=v_{\mathrm{nucl}}({\bf r})+v^{\text{ext}}({\bf r},t),$ (39)

and $v_{\mathrm{Mxc}}\left(\left[n,q\right];{\bf r},t\right)$ is the mean
field exchange-correlation potential

$v_{\mathrm{Mxc}}\left(\left[n,q\right];{\bf
r},t\right)=v_{\mathrm{M}}\left(\left[n,q\right];{\bf
r},t\right)+v_{\mathrm{xc}}\left(\left[n,q\right];{\bf r},t\right).$ (40)

Since we consider a system of electrons coupled to photonic modes, one can
split the mean-field potential $v_{\mathrm{M}}\left(\left[n,q\right];{\bf
r},t\right)$ into a contribution which comes solely from the field generated
by electrons, the Hartree potential $v_{\mathrm{H}}([n];{\bf r},t)$, and mean
field generated by photons $v_{\mathrm{P}}\left(\left[n,q\right];{\bf
r},t\right)$

$v_{\mathrm{M}}\left(\left[n,q\right];{\bf r},t\right)=v_{\mathrm{H}}([n];{\bf
r},t)+v_{\mathrm{P}}\left(\left[n,q\right];{\bf r},t\right),$ (41)

where

$v_{\mathrm{P}}\left(\left[n,q\right];{\bf
r},t\right)=\sum_{\alpha=1}^{M_{p}}\left(eg_{\alpha}\int d{\bf
r}^{\prime}n\left({\bf r}^{\prime},t\right)\bm{\varepsilon}_{\alpha}\cdot{\bf
r}^{\prime}-\sqrt{\hbar}\omega_{\alpha}q_{\alpha}(t)\right)eg_{\alpha}\bm{\varepsilon}_{\alpha}\cdot{\bf
r}.$ (42)

Similarly, one can split the exchange-correlation potential
$v_{\mathrm{xc}}\left(\left[n,q\right];{\bf r},t\right)$ into a sum of
electron–electron $v_{\mathrm{xc}}^{e-e}([n];{\bf r},t)$ and photon–electron
$v_{\mathrm{xc}}^{\mathrm{e-P}}\left(\left[n,q\right];{\bf r},t\right)$
exchange-correlation potentials,

$v_{\mathrm{xc}}\left(\left[n,q\right];{\bf
r},t\right)=v_{\mathrm{xc}}^{e-e}([n];{\bf
r},t)+v_{\mathrm{xc}}^{\mathrm{e-P}}\left(\left[n,q\right];{\bf r},t\right),$
(43)

respectively. Note that in density functional approximations employed in
practical calculations, the potentials $v_{\mathrm{xc}}^{e-e}([n];{\bf r},t)$
and $v_{\mathrm{xc}}^{\mathrm{e-P}}([n,q];{\bf r},t)$ can further dependent on
the gradient of the density or on spin densities, as well as include an
admixture of the exact exchange. Finally, we introduce KS current
$j_{\alpha,\mathrm{KS}}\left(\left[n,q_{\alpha},t\right]\right)$ in Eq. (36)

$\displaystyle\begin{split}\partial_{t}j_{\alpha,\mathrm{KS}}\left(\left[n,q_{\alpha},\partial_{t}j_{\alpha}^{\text{ext}}\right];t\right)&=\partial_{t}j_{\alpha,Mxc}\left(\left[n,q_{\alpha}\right];t\right)+\partial_{t}j^{\text{ext}}_{\alpha}(t)\\\
&=\partial_{t}j_{\alpha,M}\left(\left[n\right];t\right)+\partial_{t}j_{\alpha,xc}\left(\left[n,q_{\alpha}\right];t\right)+\partial_{t}j^{\text{ext}}_{\alpha}(t)\\\
&=-\frac{\omega_{\alpha}e}{\sqrt{2\hbar}}\int
d\mathbf{r}g_{\alpha}\bm{\varepsilon}_{\alpha}\cdot\mathbf{r}n(\mathbf{r})+\partial_{t}j^{\text{ext}}_{\alpha}(t),\end{split}$
(44)

where the exchange-correlation current $\partial_{t}j_{\alpha,xc}(t)$ is zero
in the case of dipole coupling and only the mean-field Kohn-Sham current
$\partial_{t}j_{\alpha,M}\left(\left[n\right];t\right)$ contributes [5, 69].

###  II.4 Linear response equations for relativistic QEDFT

In the present work, we assume that the external probe potential
$v^{\text{ext}}({\bf r},t)$ and current $j_{\alpha}^{\text{ext}}(t)$ are
sufficiently weak to allow the use of time-dependent perturbation theory for
the description of their effect on a molecule coupled to photon modes. Our
goal is to develop a solution of Eqs. (35) and (36) to the first order of
perturbation theory with the external potential and current acting as the
perturbation. Note that we do not assume that the coupling to the cavity modes
is weak, only that the external perturbation to the coupled cavity-matter
system is small. We will proceed in the spirit of standard time-dependent
(Rayleigh–Schrödinger) perturbation theory following our earlier work [162].
We express the time-dependent spin-orbitals as a perturbation series

$\varphi_{i}(x,t)=\varphi_{i}^{(0)}(x,t)+\lambda\varphi_{i}^{(1)}(x,t)+\lambda^{2}\varphi_{i}^{(2)}(x,t)+\ldots,$
(45)

with $\lambda$ being the perturbation parameter and $\varphi_{i}^{(k)}(x,t)$
the $k$-th order correction to the $i$-th molecular spinor. Furthermore, the
$k$-th order corrections are expanded in the basis of ground-state (canonical)
molecular spinorbitals

$\varphi_{i}^{(k)}(x,t)=d^{(k)}_{pi}(t)\varphi_{p}(x)e^{-i\varepsilon_{i}t},$
(46)

where $\varphi_{p}(x)$ are the canonical spin-orbitals obtained as the
solution of time-independent SCF procedure for $\hat{H}^{\text{KS}}(t)$ from
Eq. (37), and $d^{(k)}_{pi}(t)$ are the time-dependent expansion coefficients.
The exponential $e^{-i\varepsilon_{i}t}$ results from the solution for free
evolution $\varphi_{i}(x)e^{-i\varepsilon_{i}t}$ where the weak perturbation
then modifies $\varphi_{i}(x)$ to $\sum_{k}d^{(k)}_{pi}(t)\varphi_{p}(x)$.

Since the Kohn–Sham Hamiltonian $H^{\mathrm{KS}}(t)$ depends on the electron
density and in turn on the molecular orbitals $\varphi_{i}({\bf r},t)$, as
well as on the photonic coordinate $q_{\alpha}(t)$, it is also modified by the
perturbation and can be expanded in the orders of $\lambda$

$H^{\mathrm{KS}}=H^{\mathrm{KS,(0)}}+\lambda
H^{\mathrm{KS,(1)}}+\lambda^{2}H^{\mathrm{KS,(2)}}+\ldots$ (47)

By inserting the expansion for the Hamiltonian and the molecular orbitals into
the time-dependent Kohn–Sham equation, we obtain differential equations for
the expansion coefficients of all orders. Similarly for photons, we consider
the perturbation expansion of the displacement coordinate

$q_{\alpha}(t)=q^{(0)}_{\alpha}(t)+\lambda
q^{(1)}_{\alpha}(t)+\lambda^{2}q^{(2)}_{\alpha}(t)+\ldots$ (48)

and formulate and solve EOMs for the different orders in the expansion.
In Eqs. (45), (48) the zeroth order refers to the ground state solution, i.e.
the external probe potential $v^{\text{ext}}({\bf r},t)$ and current
$j^{\text{ext}}_{\alpha}(t)$ are absent. The stationary solution of Eq. (36)
has the form

$q_{\alpha}^{(0)}=\frac{eg_{\alpha}}{\omega_{\alpha}\sqrt{\hbar}}\int
d\mathbf{r}\bm{\varepsilon}_{\alpha}\cdot\mathbf{r}n(\mathbf{r}),$ (49)

which causes the photonic potential $v_{\mathrm{p}}$ to vanish (as can be
verified by direct substitution of Eq. (49) into Eq. (42)). The coupling to
photon modes is still present via the photon–electron XC potential
$v_{\mathrm{xc}}^{\mathrm{e-P}}\left(\left[n,q\right];{\bf r}_{l},t\right)$,
see Eq. (43), which is, however, part of the unperturbed Hamiltonian whose
eigenfunctions are the canonical orbitals. Therefore, the zeroth order
electronic EOM has the usual form

$i\partial_{t}d^{(0)}_{pi}(t)=H^{\mathrm{KS,(0)}}_{pq}d^{(0)}_{qi}(t),$ (50)

with corresponding solutions

$d^{(0)}_{pi}(t)=\delta_{pi}.$ (51)

Here, $\delta_{pi}$ is the Kronecker delta resulting from the fact that the
zeroth order KS Hamiltonian is diagonal in the basis of ground state molecular
orbitals with molecular orbital energies on the diagonal. Note that in
practical calculations we employ the pRPA which neglects
$v_{\mathrm{xc}}^{\mathrm{e-P}}\left(\left[n,q\right];{\bf r}_{l},t\right)$
after Eq. (43) and, in addition, start from the uncoupled ground state.
However, the derivation presented in the following text is general and holds
also with the full photon–electron exchange–correlation potential included.

In first order, the electronic EOM takes the form

$\displaystyle\begin{split}i\partial_{t}d^{(1)}_{ai}(t)=&A_{ai,bj}(\omega_{\mathrm{ext}})\,d^{(1)}_{bj}(t)+B_{ai,bj}(\omega_{\mathrm{ext}})\,d^{(1)*}_{bj}(t)\\\
&+\frac{e^{2}g_{\alpha}^{2}}{2}(\mathbf{r}_{ai}\cdot\bm{\varepsilon}_{\alpha})(\mathbf{r}_{jb}\cdot\bm{\varepsilon}_{\alpha})\,d^{(1)}_{bj}(t)+\frac{e^{2}g_{\alpha}^{2}}{2}(\mathbf{r}_{ai}\cdot\bm{\varepsilon}_{\alpha})(\mathbf{r}_{bj}\cdot\bm{\varepsilon}_{\alpha})\,d^{(1)*}_{bj}(t)\\\
&-eg_{\alpha}\sqrt{\hbar}\omega_{\alpha}(\mathbf{r}_{ai}\cdot\bm{\varepsilon}_{\alpha})q^{(1)}_{\alpha}(t)\\\
&+P_{\mathrm{ai}}e^{-i\omega_{\mathrm{ext}}t}+P_{\mathrm{ai}}^{*}e^{i\omega_{\mathrm{ext}}t},\end{split}$
(52)

where the right-hand side contains the electronic kernel, the interaction with
the cavity modes, and the interaction with the external field. Specifically,
in Eq. (52), the first line on the right-hand side contains the usual Coulomb,
exchange, and exchange–correlation coupling matrices

$\displaystyle A_{ai,bj}(\omega_{\mathrm{ext}})$
$\displaystyle=(\varepsilon_{a}-\varepsilon_{i})\delta_{ab}\delta_{ij}+K_{\mu\nu\lambda\tau}(\omega_{\mathrm{ext}})C_{\mu
a}^{*}C_{\nu i}C_{\lambda j}^{*}C_{\tau b},$ (53a) $\displaystyle
B_{ai,bj}(\omega_{\mathrm{ext}})$
$\displaystyle=K_{\mu\nu\lambda\tau}(\omega_{\mathrm{ext}})C_{\mu a}^{*}C_{\nu
i}C_{\lambda b}^{*}C_{\tau j},$ (53b) $\displaystyle
K_{\mu\nu\lambda\tau}(\omega_{\mathrm{ext}})$
$\displaystyle=K_{\mu\nu\lambda\tau}^{\mathrm{HF}}(\xi)+K_{\mu\nu\lambda\tau}^{\mathrm{XC}}(\xi)(\omega_{\mathrm{ext}})+K_{\mu\nu\lambda\tau}^{\mathrm{e-P}}(\omega_{\mathrm{ext}}),$
(53c)

with molecular orbital energies $\varepsilon_{p}\ (p=a,i)$ and coefficients
$C_{\mu p}$ obtained from the solution of the ground-state SCF procedure for
the time-independent Dirac–Coulomb equation, and $K_{\mu\nu\lambda\tau}$ being
matrix elements of the response kernel with HF, XC, e–P denoting its
Hartree–Fock, electron–electron exchange–correlation, and photon–electron
parts, respectively, in hybrid functionals mixed with weight $\xi$. We note
that, strictly speaking, in the case of a hybrid auxiliary systems the
respective exchange-correlation contributions of the electron-electron and
electron-photon parts get modified, but we will not indicate this for
notational simplicity [163]. The XC kernels can be in general non-adiabatic,
denoted by their dependence on $\omega_{\mathrm{ext}}$. The second and the
third line contain terms arising from the mean-field coupling to the photon
field, with the dipole self-energy terms in the second line and in the third
one the interaction with the displacement field of the cavity modes
$q^{(1)}_{\alpha}(t)$. Finally, the last line contains the coupling to the
external classical electromagnetic field representing the probe field used to
perform absorption spectroscopy on the molecule, i.e. the perturbation.

For photons, the EOM in the first-order reads

$\left(\frac{\partial^{2}}{\partial
t^{2}}+\omega_{\alpha}^{2}\right)q^{(1)}_{\alpha}(t)=-\sqrt{2}\left[(\partial_{t}j^{(1)}_{\alpha,\mathrm{KS}}\left(\left[n,q_{\alpha},t\right]\right)+(\partial_{t}j^{\mathrm{ext}}(t))\right],$
(54)

where the first-order contribution to the current comes only from its
dependence on the (first-order) electron density in the form

$\partial_{t}j^{(1)}_{\alpha,\mathrm{KS}}\left(\left[n,t\right]\right)=-\frac{g_{\alpha}\omega_{\alpha}e\int
d^{3}rn^{(1)}(r,t)(\mathbf{r}\cdot\bm{\varepsilon}_{\alpha})}{\sqrt{2\hbar}}.$
(55)

By expressing the first-order electron density via the canonical molecular
orbitals

$n^{(1)}(r,t)=\varphi_{j}^{\dagger}(x)d^{(1)}_{bj}(t)\varphi_{b}(x)+d^{(1)*}_{bj}(t)\varphi_{b}^{\dagger}(x)\varphi_{j}(x),$
(56)

we obtain the EOM for $q^{(1)}_{\alpha}(t)$ coupled to electronic coefficients
$d^{(1)}_{ai}(t)$ and $d^{(1)*}_{ai}(t)$ which reads

$\left(\frac{\partial^{2}}{\partial
t^{2}}+\omega_{\alpha}^{2}\right)q^{(1)}_{\alpha}(t)=\frac{g_{\alpha}\omega_{\alpha}e}{\sqrt{\hbar}}\left[(\mathbf{r}_{jb}\cdot\bm{\varepsilon}_{\alpha})\,d^{(1)}_{bj}(t)+(\mathbf{r}_{bj}\cdot\mathbf{\varepsilon}_{\alpha})\,d^{(1)*}_{bj}(t)\right]-\sqrt{2}\partial_{t}j^{\mathrm{ext}}(t).$
(57)

Eq. (57) together with Eq. (52) and its complex conjugate constitute a system
of coupled differential equations for variables $q^{(1)}_{\alpha}(t)$,
$d^{(1)}_{ai}(t)$, and $d^{(1)*}_{ai}(t)$. Note that since $q_{\alpha}(t)$ is
real, we do not need a separate equation for $q^{(1)*}_{\alpha}(t)$. However,
while Eq. (52) is a first-order differential equation, Eq. (57) is a second-
order equation. We can transform it to a first-order equation by defining a
new variable $p^{(1)}_{\alpha}(t)$ as

$p^{(1)}_{\alpha}(t)=-i\frac{\partial}{\partial t}q^{(1)}_{\alpha}(t),$ (58)

which obeys a first-order differential equation

$i\frac{\partial}{\partial
t}p^{(1)}_{\alpha}(t)=-\omega_{\alpha}^{2}q^{(1)}_{\alpha}(t)+\frac{g_{\alpha}e\omega_{\alpha}}{\sqrt{\hbar}}\left[(\mathbf{r}_{jb}\cdot\bm{\varepsilon}_{\alpha})\,d^{(1)}_{bj}(t)+(\mathbf{r}_{bj}\cdot\bm{\varepsilon}_{\alpha})\,d^{(1)*}_{bj}(t)\right]-\sqrt{2}\partial_{t}j^{\mathrm{ext}}(t).$
(59)

Thus we have replaced a second-order differential equation (57) with two
first-order differential equations (58) and (59) The imaginary unit in Eq.
(58) was added for convenience in order to obtain a first-order EOM containing
$i\partial_{t}$ in accord with Eq. (52).

The final complete system of coupled differential equations for first-order
corrections to the electronic and photonic degrees of freedom is

$\displaystyle i\partial_{t}d^{(1)}_{ai}(t)$
$\displaystyle=(A_{ai,bj}+\Delta_{ai,bj})d^{(1)}_{bj}(t)+(B_{ai,bj}+\Delta^{\prime}_{ai,bj})d^{(1)*}_{bj}(t)-L_{ai,\alpha}q^{(1)}_{\alpha}(t)+P_{\mathrm{ai}}e^{-i\omega_{\mathrm{ext}}t}+P_{\mathrm{ai}}^{*}e^{i\omega_{\mathrm{ext}}t},$
(60a) $\displaystyle-i\partial_{t}d^{(1)*}_{ai}(t)$
$\displaystyle=(A^{*}_{ai,bj}+\Delta^{*}_{ai,bj})d^{(1)*}_{bj}(t)+(B^{*}_{ai,bj}+\Delta^{\prime*}_{ai,bj})d^{(1)}_{bj}(t)-L^{*}_{ai,\alpha}q^{(1)}_{\alpha}(t)+P^{*}_{\mathrm{ai}}e^{i\omega_{\mathrm{ext}}t}+P_{\mathrm{ai}}e^{-i\omega_{\mathrm{ext}}t},$
(60b) $\displaystyle i\frac{\partial}{\partial t}p^{(1)}_{\alpha}(t)$
$\displaystyle=+Q_{\alpha,bj}^{\prime*}d^{(1)}_{bj}(t)+Q^{\prime}_{\alpha,bj}d^{(1)*}_{bj}(t)-\omega_{\alpha}^{2}q^{(1)}_{\alpha}(t)-Je^{-i\omega_{\mathrm{ext}}t}-J^{*}e^{i\omega_{\mathrm{ext}}t},$
(60c) $\displaystyle-i\frac{\partial}{\partial t}q^{(1)}_{\alpha}(t)$
$\displaystyle=p^{(1)}_{\alpha}(t),$ (60d)

where we introduced the notation

$\displaystyle\Delta_{ai,bj}$
$\displaystyle=\frac{e^{2}g_{\alpha}^{2}}{2}(\mathbf{r}_{ai}\cdot\bm{\varepsilon}_{\alpha})(\mathbf{r}_{jb}\cdot\bm{\varepsilon}_{\alpha}),$
(61a) $\displaystyle\Delta^{\prime}_{ai,bj}$
$\displaystyle=\frac{e^{2}g_{\alpha}^{2}}{2}(\mathbf{r}_{ai}\cdot\bm{\varepsilon}_{\alpha})(\mathbf{r}_{bj}\cdot\bm{\varepsilon}_{\alpha}),$
(61b) $\displaystyle L_{ai,\alpha}$
$\displaystyle=-eg_{\alpha}\sqrt{\hbar}\omega_{\alpha}(\mathbf{r}_{ai}\cdot\bm{\varepsilon}_{\alpha}),$
(61c) $\displaystyle Q^{\prime}_{ai,\alpha}$
$\displaystyle=\frac{g_{\alpha}e\omega_{\alpha}}{\sqrt{\hbar}}(\mathbf{r}_{jb}\cdot\bm{\varepsilon}_{\alpha}),$
(61d)

and where we considered the external current to have a time dependence given
by
$\sqrt{2}\partial_{t}j^{\mathrm{ext}}(t)=Je^{-i\omega_{\mathrm{ext}}t}+J^{*}e^{i\omega_{\mathrm{ext}}t}$.
The explicit dependence of coupling matrices on $\omega_{\mathrm{ext}}$
resulting from such dependence in the exchange–correlation kernel was omitted
for brevity. The first-order equations of motion for the coupled
electron–photon system, Eqs. (60), can be solved by different approaches such
as the method of undetermined coefficients shown in Appendix B, or by the
direct solution of the system of equations as a matrix differential equation,
followed by a similarity transformation shown in Appendix C. Both of these
approaches lead to an algebraic equation corresponding to an eigenvalue
problem whose eigenvalues are excitation energies and eigenvectors transition
vectors of the coupled light–matter system. Alternatively, one can arrive at
the final equation for excitation energies by considering density–density
linear response function as presented at the non-relativistic level of theory
in Ref. 69.

The linear response equation obtained from Eqs. (60) has the form of the
Casida equation [164, 165, 161], of linear-response time dependent density
functional theory (TDDFT) extended by blocks describing photons

$\begin{pmatrix}\mathbf{A}+\mathbf{\Delta}&\mathbf{B}+\mathbf{\Delta^{\prime}}&-\mathbf{L}&-\mathbf{L}\\\
\mathbf{B}^{*}+\mathbf{\Delta^{\prime}}^{*}&\mathbf{A}^{*}+\mathbf{\Delta}^{*}&-\mathbf{L}^{*}&-\mathbf{L}^{*}\\\
-\mathbf{Q}&-\mathbf{Q}^{*}&\bm{\omega}&\mathbf{0}\\\
-\mathbf{Q}&-\mathbf{Q}^{*}&\mathbf{0}&\bm{\omega}\\\
\end{pmatrix}\begin{pmatrix}\mathbf{X}_{n}\\\ \mathbf{Y}_{n}\\\
\mathbf{M}_{n}\\\ \mathbf{N}_{n}\\\
\end{pmatrix}=\Omega_{n}\begin{pmatrix}\mathbf{1}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&-\mathbf{1}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{1}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}&-\mathbf{1}\\\
\end{pmatrix}\begin{pmatrix}\mathbf{X}_{n}\\\ \mathbf{Y}_{n}\\\
\mathbf{M}_{n}\\\ \mathbf{N}_{n}\\\ \end{pmatrix},$ (62)

where the terms $\mathbf{A}$ and $\mathbf{B}$ are the electronic coupling
matrices introduced in Eqs. (53), the electron–photon coupling and self-energy
terms $\mathbf{L}$,
$Q_{jb,\alpha}=\frac{1}{2\omega_{\alpha}}Q^{\prime}_{jb,\alpha}=\frac{1}{2}\frac{g_{\alpha}e}{\sqrt{\hbar}}(\mathbf{r}_{jb}\cdot\bm{\varepsilon}_{\alpha})$,
$\mathbf{\Delta}$, and $\mathbf{\Delta^{\prime}}$ are defined in Eqs. (53).
The solutions of Eq. (62) are the eigenvalues $\Omega_{n}$ corresponding to
excitation energies of the coupled light-matter system and the eigenvectors
that contain the electronic excitation and deexcitation amplitudes
$\mathbf{X}_{n}$ and $\mathbf{Y}_{n}$ parametrizing $d^{(1)}_{ai}$ and
$d^{(1)*}_{ai}$, together with the photonic creation and annihilation
amplitudes $\mathbf{M}_{n}$ and $\mathbf{N}_{n}$, parametrizing
$q^{(1)}_{\alpha}$ and $p^{(1)}_{\alpha}$ (for this interpretation see
Appendix C). The linear response equation (62) is formulated in the basis of
the ground state molecular orbitals where the dimensions of the terms in the
equation expressed via the number of occupied orbitals $N_{\mathrm{occ}}$, the
number of virtual orbitals $N_{\mathrm{vir}}$, and the number of photon modes
$N_{\mathrm{ph}}$ are

  * •

$\mathbf{A}$, $\mathbf{B}$, $\mathbf{\Delta}$, $\mathbf{\Delta^{\prime}}$:
matrices of dimension $N_{\mathrm{vir}}N_{\mathrm{occ}}\times
N_{\mathrm{vir}}N_{\mathrm{occ}}$,

  * •

$\mathbf{L}$, $\mathbf{Q}$: matrices of dimension
$N_{\mathrm{vir}}N_{\mathrm{occ}}\times N_{\mathrm{ph}}$,

  * •

$\bm{\omega}$: diagonal matrix of dimension $N_{\mathrm{ph}}\times
N_{\mathrm{ph}}$,

  * •

$\mathbf{X}_{n}$, $\mathbf{Y}_{n}$: vectors of dimension
$N_{\mathrm{vir}}N_{\mathrm{occ}}$,

  * •

$\mathbf{M}_{n}$, $\mathbf{N}_{n}$: vectors of dimension $N_{\mathrm{ph}}$.

Note that while the transition vectors $\mathbf{X}_{n}$, $\mathbf{Y}_{n}$ are
originally defined as matrices of dimension $N_{\mathrm{vir}}\times
N_{\mathrm{occ}}$ ($X_{ai}$, $Y_{ai}$ where $a\in\mathrm{vir}$,
$i\in\mathrm{occ}$), for the purposes of solution and computer implementation
of Eq. (62), they are unrolled into one-dimensional vectors ($X_{\kappa}$,
$Y_{\kappa}$, where $\kappa$ runs over all products $ai$). Similarly, the
terms $\mathbf{A}$, $\mathbf{B}$, $\mathbf{\Delta}$,
$\mathbf{\Delta^{\prime}}$ originally defined as rank 4 tensors, and the terms
$\mathbf{L}$, $\mathbf{Q}$ defined as rank 3 tensors, became matrices by
considering virtual–occupied pair $ai$ to label a single dimension. The matrix
on the left-hand side of Eq. (62) is not symmetric which has consequences for
the calculation of spectra as detailed in Section II.5 and Appendix D.
Previous works on linear response QEDFT either transformed Eq. (62) into an
equation of half the original dimension with eigenvalue $\Omega_{n}^{2}$ made
possible due to real-valued orbitals used in the non-relativistic theory, [69]
or worked with a symmetric equation obtained due to approximations made in the
electron–photon Hamiltonian [70, 71]. In Appendix C we show that Eq. (62) can
be symmetrized by a similarity transformation to a complex version of this
symmetric equation, due to equivalent approximations made in the KS potential
and current in the practical implementation.

###  II.5 Calculation of spectra

The frequency-dependent polarizability [166] that defines the absorption
spectrum is given by (compare with Eq. (75))

$\bm{\alpha}(\omega_{\mathrm{ext}})=\begin{pmatrix}\mathbf{P}^{*}&\mathbf{P}&\mathbf{0}&\mathbf{0}\end{pmatrix}\begin{bmatrix}\begin{pmatrix}\mathbf{A}+\mathbf{\Delta}&\mathbf{B}+\mathbf{\Delta}^{\prime}&-\mathbf{L}&-\mathbf{L}\\\
\mathbf{A}^{*}+\mathbf{\Delta}^{*}&\mathbf{B}^{*}+\mathbf{\Delta}^{\prime*}&-\mathbf{L}^{*}&-\mathbf{L}^{*}\\\
-\mathbf{Q}&-\mathbf{Q}^{*}&\bm{\omega}&\mathbf{0}\\\
-\mathbf{Q}&-\mathbf{Q}^{*}&\mathbf{0}&\bm{\omega}\\\
\end{pmatrix}-\omega_{\mathrm{ext}}\begin{pmatrix}\mathbf{1}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&-\mathbf{1}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{1}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}&-\mathbf{1}\\\
\end{pmatrix}\end{bmatrix}^{-1}\begin{pmatrix}\mathbf{P}\\\ \mathbf{P}^{*}\\\
\mathbf{0}\\\ \mathbf{0}\end{pmatrix},$ (63)

where $\mathbf{P}$ is the matrix representation of the electric dipole moment
operator. We can replace the matrix $[\ldots]^{-1}$ in Eq. (63) by its
spectral decomposition from eigenvectors obtained by solving Eq. (62).
However, the electron–photon Casida equation, Eq. (62), corresponds to a
generalized eigenvalue problem of a non-Hermitian non-symmetric matrix (unlike
the electronic Casida equation containing only the Hermitian electron–electron
block $\begin{pmatrix}\mathbf{A}&\mathbf{B}\\\
\mathbf{B}^{*}&\mathbf{A}^{*}\end{pmatrix}$ ). Therefore, the eigenproblem has
different left and right eigenvectors, a fact that has to be taken into
account when performing the spectral decomposition. In this case, the
expression for the frequency-dependent linear polarizability tensor derived in
Appendix D has the from

$\bm{\alpha}(\omega)=\sum_{n}\left[\frac{t^{\mathrm{R*}}_{n}t^{\mathrm{L}}_{n}}{\Omega_{n}+\omega}-\frac{t^{\mathrm{R}}_{n}t^{\mathrm{L*}}_{n}}{\omega-\Omega_{n}}\right],$
(64)

where index $n$ runs over positive energy eigenvalues and the left and right
transition dipole moments are

$\displaystyle t^{\mathrm{R}}_{n}$
$\displaystyle=\mathbf{P}^{\dagger}\mathbf{X}^{\mathrm{R}}_{n}+\mathbf{P}\mathbf{Y}^{\mathrm{R}}_{n},$
(65a) $\displaystyle t^{\mathrm{L}}_{n}$
$\displaystyle=\mathbf{P}^{\dagger}\mathbf{X}^{\mathrm{L}}_{n}-\mathbf{P}\mathbf{Y}^{\mathrm{L}}_{n},$
(65b)

with $\mathbf{X}^{\mathrm{R/L}}_{n}$ and $\mathbf{Y}^{\mathrm{R/N}}_{n}$
corresponding to electronic parts of the right/left $n$-th eigenvector,
respectively. An empirical broadening parameter $\gamma$ can be added to
obtain spectra with finite peak widths by replacing $\omega$ with
$\omega+i\gamma$ in Eq. (64). Note, however, that the empirical broadening
parameter is only a convenience used when representing a cavity by a single or
few modes. Physical radiative broadening can be described in linear response
QEDFT by considering coupling to a dense spectrum of photon modes [69]. The
absorption spectrum of a molecule in a cavity is then calculated using the
usual formula for the dipole strength function

$S(\omega)=\frac{4\pi\omega}{3c}\Im\,\textrm{Tr}\left[\bm{\alpha}(\omega)\right],$
(66)

where $c$ is the speed of light, $\Im$ denotes the imaginary part, and Tr is
the trace over the Cartesian components. Note that while this is the formula
used in QEDFT literature, the issue of calculating a proper rotationally
averaged spectral function for the molecule-in-cavity set-up comparable to
experiments is more subtle and requires further work [167, 160].

###  II.6 Numerical solver

Eq. (62) is thus an eigenproblem for a large matrix whose construction and
solution comprises the bulk of the computational effort of relativistic linear
response QEDFT. The equation was implemented into the relativistic quantum
chemistry DFT package ReSpect [168] and is built upon an earlier four-
component molecular linear response TDDFT code [169] that solves the Casida
equations for electrons. The linear response QEDFT equation, Eq. (62), is
solved using an iterative Davidson–Olsen solver [170, 171] as is common in
TDDFT implementations, here extended to consider the photon part as well. The
solver expresses the solution vector $\mathbf{Z}=(\mathbf{X}\ \mathbf{Y}\
\mathbf{M}\ \mathbf{N})^{\mathrm{T}}$ as a linear combination of so-called
trial vectors, i.e. $\mathbf{Z}=\sum_{\tau}c_{\tau}\mathbf{t}_{\tau}$.
Moreover, the solvers developed for electronic relativistic linear response
TDDFT employ a paired structure of the trial vectors [172, 173, 169], where
each trial vector $\mathbf{t}_{\tau}=(\mathbf{t}^{X}_{\tau}\
\mathbf{t}^{Y}_{\tau})^{\mathrm{T}}$, a pair vector
$\tilde{\mathbf{t}}_{\tau}=(\mathbf{t}^{Y*}_{\tau}\
\mathbf{t}^{X*}_{\tau})^{\mathrm{T}}$ is considered. However, since the
calculation of polarizability requires the left eigenvectors, these also have
to be correctly described by the trial basis. To this end, trial vectors with
a special internal structure informed by the analytic relationship between the
left and right eigenvectors have to be added as well. Such an extension was
first suggested by Komorovsky, Cherry, and Repisky [169] for open-shell TDDFT,
and includes adding a trial vector
$\mathbf{t}_{\tau}^{\mathrm{L}}=(\mathbf{t}^{X}_{\tau}\
-\mathbf{t}^{Y}_{\tau})^{\mathrm{T}}$ for each vector
$\mathbf{t}_{\tau}\equiv\mathbf{t}_{\tau}^{\mathrm{R}}$ as well as a vector
$\tilde{\mathbf{t}}_{\tau}^{\mathrm{L}}=(\mathbf{t}^{Y*}_{\tau}\
-\mathbf{t}^{X*}_{\tau})^{\mathrm{T}}$ for its pair
$\tilde{\mathbf{t}}_{\tau}\equiv\tilde{\mathbf{t}}_{\tau}^{\mathrm{R}}$. In
the case of linear response QEDFT describing a coupled electron–photon system,
the special structure of paired and left trial vectors is translated to the
photon part as well giving

$\begin{pmatrix}\mathbf{X}\\\ \mathbf{Y}\\\ \mathbf{M}\\\ \mathbf{N}\\\
\end{pmatrix}=\sum_{\tau}c_{\tau}^{\mathrm{R}}\begin{pmatrix}\mathbf{t}_{\tau}^{X}\\\
\mathbf{t}_{\tau}^{Y}\\\ \mathbf{t}_{\tau}^{M}\\\ \mathbf{t}_{\tau}^{N}\\\
\end{pmatrix}+\sum_{\tau}\tilde{c}_{\tau}^{\mathrm{R}}\begin{pmatrix}\mathbf{t}_{\tau}^{Y*}\\\
\mathbf{t}_{\tau}^{X*}\\\ \mathbf{t}_{\tau}^{N*}\\\ \mathbf{t}_{\tau}^{M*}\\\
\end{pmatrix}+\sum_{\tau}c_{\tau}^{\mathrm{L}}\begin{pmatrix}\mathbf{t}_{\tau}^{X}\\\
-\mathbf{t}_{\tau}^{Y}\\\ \mathbf{t}_{\tau}^{M}\\\ -\mathbf{t}_{\tau}^{N}\\\
\end{pmatrix}+\sum_{\tau}\tilde{c}_{\tau}^{\mathrm{L}}\begin{pmatrix}\mathbf{t}_{\tau}^{Y*}\\\
-\mathbf{t}_{\tau}^{X*}\\\ \mathbf{t}_{\tau}^{N*}\\\
-\mathbf{t}_{\tau}^{M*}\\\ \end{pmatrix},$ (67)

as the ansatz for the solution. The expansion coefficients $c_{\tau}$ and
$\tilde{c}_{\tau}$ are determined by solving Eq. (62) projected onto the
subspace defined by the trial vectors. New trial vectors are generated in each
iteration by preconditioning the residue vectors where each desired eigenvalue
can contribute up to two trial vectors ($\mathbf{t}$ and its pair
$\tilde{\mathbf{t}}$) in each iteration. The subspace is expanded until
convergence determined by the norm of the residue vector for the particular
frequency is reached. The solutions for the desired number of eigenvalues are
sought in a common basis composed of trial vectors contributed by all the
eigenvalues. New trial vectors are orthogonalized i) with respect to the
previously added trial vectors using Gram–Schmidt orthogonalization and
Kahan–Partlett “twice is enough” algorithm; and ii) within the
$(\mathbf{t},\tilde{\mathbf{t}})$ pair using an analytic formula. The vectors
that are reduced to a length below a set threshold during the
orthogonalization are discarded from the trial basis. This protocol achieves
robust convergence of the solver across different molecular systems and the
number of desired eigenvalues. To cross-validate our results, we implemented
the symmetrized version of the equation as well and confirmed that these two
equations lead to the same spectra.

##  III Computational details

In the benchmark calculations presented below we consider the series of group
12 atoms (Zn, Cd, Hg), and a large mercury porphyrin (Hg@porphyrin). The
geometry of the Hg porphyrin was taken from Ref. 135 and corresponds to a non-
planar configuration of C_4v symmetry.

All spectra were calculated using the relativistic spectroscopy DFT program
ReSpect. [168] The reference orbitals for the linear response calculations
were obtained from the ground state uncoupled to photons. The relativistic
four-component molecular spinorbitals for electrons were described Gaussian-
type orbital (GTO) basis sets where the scalar GTOs were used for the large
component and the basis for the small-component basis was generated by
imposing the restricted kinetically balanced (RKB) relation. [174] The
selected scalar basis sets were the uncontracted Dyall’s VDZ basis sets for
metals (Zn, Cd, Hg) [175, 176, 177, 178] and the uncontracted Dunning’s cc-
pVDZ basis sets [179] for light elements. and the B3LYP XC approximation of
the density functional. [180, 181, 182, 183, 184] The numerical integration of
the noncollinear exchange–correlation potential and kernel was done with an
adaptive molecular grid of medium size (program default). In the 4c
calculations, atomic nuclei of finite size were approximated by a Gaussian
charge distribution model.

The electron–photon correlation is treated under the photon random phase
approximation (pRPA) [69, 72] which amounts to disregarding the
photon–electron exchange–correlation potential
$v_{\mathrm{xc}}^{\mathrm{e-P}}\left(\left[n,q\right];{\bf r}_{l},t\right)$
and kernel $K_{\mu\nu\lambda\tau}^{\mathrm{e-P}}$ in Eqs. (43) and (53c),
respectively. Moreover, we consider the adiabatic approximation in the
construction of the XC kernels, i.e. all dependencies on
$\omega_{\mathrm{ext}}$ in Eqs. (53) are dropped. At present, we restrict
ourselves to closed-shell systems for which the adiabatic noncollinear XC
kernel has the form given in Eqs. 24-26c in Ref. 162. The eigenvalue equation
was solved iteratively for the first 50 or 100 excitation energies for systems
coupled to a single cavity mode. The spectra were subsequently evaluated from
the excitation energies and transition moments obtained from the left and
right eigenvectors, in the case of coupling to a single cavity mode, empirical
Lorentzian broadening was used to obtain finite-width peaks.

##  IV Results and discussion

###  IV.1 Atoms with singlet–triplet transitions

First, we consider the case of atoms of group 12 metals, Zn, Cd, and Hg whose
spectra were previously studied at the level of four-component Dirac–Coulomb
Hamiltonian using linear response [185, 173] and real-time TDDFT [186]. The
first two excited states correspond to the triplet state ^3P and the singlet
state ^1P, where the transition from the ^1S0 ground state to the excited
triplet state is forbidden in theories neglecting spin–orbit coupling (SOC).
With SOC considered, the transition to ^3P1 becomes allowed while the
transitions to the other two triplet components (^3P0, ^3P2) remain forbidden
due to the $\Delta J$ selection rule within the dipole approximation. The
excitation energies for the allowed ^1S0 $\rightarrow$ ^3P1 and ^1S0
$\rightarrow$ ^1P1 transitions are summarized in Table S1 in the Supplemental
Material while the calculated values reproduce reference experimental data
[187] reasonably well given the functional and basis set size. Both the ^1S0
$\rightarrow$ ^3P1 and ^1S0 $\rightarrow$ ^1P1 transitions are composed of
three degenerate spectral lines. In the orbital picture, the highest occupied
shell corresponds to $n$s orbitals and the lowest unoccupied to $n$p while
after SOC is taken into account the p shell splits into $n$p_1/2 and $n$p_3/2
levels. There are in total 12 possible single-electron excitations between
$n$s and $n$p_1/2 and $n$p_3/2 orbitals and these constitute the dominant
contributions to TDDFT transition vectors. However, these 12 simple
configurations cannot be straightforwardly identified with the components of
the singlet and triplet excited states as their transition vectors contain
more than one significant orbital pair contribution.

To investigate the effect of the strong coupling to photons on these atomic
spectra, we solve the linear response QEDFT equation for a single atom coupled
to a single effective cavity mode with the cavity frequency scanning the
region around the singlet–triplet (S–T) and singlet–singlet (S–S) transitions.
By focusing on a single efficient cavity mode, it is implied that the
remaining spectrum of modes is accounted for by working with the renormalized
(observable) masses of the electrons [7, 188].

The results of series of calculations with different cavity frequencies are
excitation energies and corresponding transition moments from which it is
possible to calculate absorption spectra. Such spectra are depicted for Hg
atom in Figure 2(c) showing both non-relativistic one-component (1c) and fully
relativistic four-component (4c) calculations (results for Zn and Cd atoms are
in the Supplemental Material, Figure S1). The $x$-axis shows the cavity
frequency, the $y$-axis corresponds to excitation energy while the magnitude
of the spectral function is depicted as the color gradient in logarithmic
($\log_{10}$) scale to allow comparison between the S–T and S–S transitions of
vastly different magnitudes. The spectra were artificially broadened using
Lorentzian broadening with broadening parameter $\gamma=0.00002\,\mathrm{au}$
(see text after Eqs. (65)). Only excited states with nonzero transition dipole
moments are visible as signals in the spectra, i.e. ^1P1 and ^3P1, while the
other components of the triplet state with total angular momentum $J=0$ and
$J=2$ are dark.

The plots depict a common scenario wherein light-matter mixing increases as
the cavity frequency approaches a resonance with an excitation energy. Both
the ^3P1 and ^1P1 excited states are originally triply degenerate with the
degenerate levels having perpendicular transition dipole moments. Therefore,
the photon state mixes with only one of the three states – the one with the
transition dipole moment parallel with the photon polarization – and the
uncoupled excitation energies are still visible due to the two remaining
levels of both ^3P1 and ^1P1. When compared with the non-relativistic picture,
the relativistic description causes a shift in the energy of the excited
singlet state – the non-relativistic spectra were artificially shifted so that
the S–S excitation energies would be aligned. The specific shifts were
$0.16\,\mathrm{eV}$, $0.43\,\mathrm{eV}$, and $1.57\,\mathrm{eV}$, for Zn, Cd,
and Hg, respectively. Furthermore, the inclusion of relativity introduces a
transition to the excited triplet state. The intensity of the S–T transition
as well as the spectral shifts are more prominent for heavier atoms. As the
transition to the excited triplet state becomes allowed, a mixing with the
cavity mode is observed near the resonant frequency. Since the coupling to the
cavity mode is proportional to the transition dipole moment, which for the S–T
transition increases with atomic number, the Rabi splitting of the polaritons
originating from the triplet state is the largest for mercury. From the
zoomed-in insets in the dispersion plots in Figure 2(c) we see that the
maximum (50:50) light–matter mixing does not occur at exactly the S–T resonant
frequency, but rather at a higher frequency. The maximum light-matter mixing
around the S–S resonance also occurs at a frequency shifted from the resonance
by approximately the same amount as the around the S–T transition.

We calculated absorption spectra with cavities tuned to these maximum mixing
frequencies, i.e. $4.0204\,\mathrm{eV}$, $3.744\,\mathrm{eV}$, and
$4.724\,\mathrm{eV}$, for Zn, Cd, and Hg, respectively. The spectra for Zn and
Hg are depicted in Figure 3 while the spectrum for Cd is available in the
Supplemental Material, Figure S2. As discussed above, two S–T transitions with
perpendicular transition dipole are still present as excitation energies and
visible in the spectra of Cd and Hg as a line with 2/3 of the intensity seen
in the uncoupled system (for Zn, the extremely weak line corresponding to the
uncoupled S–T transitions is obscured by peak broadening). In the spectral
region around the S–T transition, the spectra contain very intense polaritonic
peaks with the increase of the transition moment with respect to the original
S–T transition ranging from a striking 50-fold increase for Zn to about 20%
for Hg (compared to the intensity of one of the three degenerate states
meaning that the peak composed of two degenerate uncoupled S–T states is more
intense than the polaritonic peak in Hg).

While it is tempting to suggest an interpretation based on cavity enhancement
of the S–T transition, it is important to note that the peak appears also in
the one-component non-relativistic (1c) spectra. Since the 1c spectra are
evaluated from the non-relativistic linear response calculations that do not
contain the S–T transition, the peak must correspond to the lower polariton
originating from the mixing of cavity mode with S–S exciton. Even if the
cavity frequency is far from the resonance with the S–S transition and the
polaritonic state is light-like ($>90\%$), its intensity is sufficient to be
comparable with the S–T transition in the relativistic (4c) calculations. This
is in agreement with the 2D spectra in Figure 2(c) where the intensity of the
lower polaritonic branch in the 1c calculations is more intense (Zn) or
comparable (Hg, also Cd) with the S–T transition in the 4c spectra.

In the 4c calculations, the effects of the coupling to the low-intensity S–T
transition and the intense but off-resonant S–S transition combine. Besides
the spectral signature in the S–T region discussed above, in the region around
the S–S transition, one of the originally three degenerate spectral lines is
blueshifted and its response vector contains a small but non-negligible
($10^{-3}$) photonic contribution. Most of the intensity increase of the
polaritonic peaks arising from the S–T transition is accounted for by the
corresponding decrease in intensity of the shifted S–S transition. The
decrease in intensity of the shifted S–S line is in fact larger than the gain
of the S–T transition, with the gain being about 80% of the decrease. The
Thomas–Reiche–Kuhn sum rule is then preserved by states higher in the spectrum
when one state of three degenerate states is also slightly shifted in energy
and its intensity is now higher. This repeats for all triply degenerate states
at higher energies. The sum of photonic contributions from the lower and upper
S–T polariton and the shifted S–S line is 1 with the accuracy of
$10^{-4}$–$10^{-5}$ and lower for the higher lying states. Note that in pRPA
QEDFT with single cavity mode, the total sum of photon contributions should be
1 and only one extra peak can appear in the spectrum of an atom or molecule
coupled to the cavity.

These results show the importance of the off-resonant coupling where the
photon mode is coupled not only to the resonant S–T exciton but also to the
off-resonant yet high-transition-moment S–S exciton as well as _all higher-
lying states_. In turn, this demonstrates the importance of ab initio QEDFT
treatment: a simple two-state model built from the S–T excitation energy and
corresponding transition dipole moment would produce two polaritonic states
whose intensity would sum to the original S–T intensity. Another effect beyond
the two-state model is the shift of S–T resonance (50:50 light-matter mixing)
toward higher frequencies. This becomes more pronounced as the coupling
strength increases as shown in Figure 4(c). Here, the results of QEDFT
calculations are compared with Jaynes–Cummings (JC) energies
$E_{\pm}=(\hbar\omega_{\alpha}-\frac{1}{2}\hbar\Delta)\pm\frac{1}{2}\hbar\sqrt{\Delta^{2}+\Omega^{2}}$,
where $\Delta=\omega_{\alpha}-\omega_{\mathrm{S-T}}$ and
$\Omega=2g_{\alpha}\sqrt{\frac{\omega_{\mathrm{\alpha}}}{2\hbar}}\left\langle{}^{3}\mathrm{P}_{1}|\bm{\mu}\cdot\bm{\varepsilon}_{\alpha}|^{1}\mathrm{S}_{0}\right\rangle$
with S–T excitation energy $\omega_{\mathrm{S-T}}$ and transition dipole
moment
$\left\langle{}^{3}\mathrm{P}_{1}|\bm{\mu}\cdot\bm{\varepsilon}_{\alpha}|^{1}\mathrm{S}_{0}\right\rangle$
taken from a TDDFT calculation without coupling to cavity.

While models can be extended to account for more than two states, it is
impossible to a priori say how many excited states should be included,
particularly for systems with dense spectra such as molecules. Therefore, ab
initio calculations are an indispensable tool for qualitative and quantitative
predictions of properties of strongly bound light–matter systems, as well as
for guiding the development of models. The quest to extend their capabilities
to include further effects, such as relativity and spin–orbit coupling, is
thus an important direction in polaritonic chemistry.

Figure 2: Absorption spectra of Hg atom in a cavity strongly coupled
($g_{\alpha}=0.01\,\mathrm{au}$) to photonic modes, non-relativistic one-
component (1c) vs relativistic four-component (4c) calculations.

(a) Hg 1c

(b) Hg 4c

(c) Hg 4c zoom

Figure 3: Absorption spectra of Zn and Hg atoms in cavities strongly coupled
($g_{\alpha}=0.01\,\mathrm{au}$) to a cavity set to effective resonance with
the singlet–triplet (S–T) transition defined by 50:50 light–matter mixing of
the lower polariton rather than by the numerical value compared to reference
spectra of free atoms without cavities. The region around the low-intensity
S–T transition is magnified to ease reading by a factor specified in each
figure.

(a) Zn

(b) Hg

Figure 4: Absorption spectra of Hg atom in a cavity with different coupling
strength $g_{\alpha}$ compared to a two-level Jaynes–Cummings model (black
lines) considering only the S–T transition.

(a) $g_{\alpha}=0.02$ au

(b) $g_{\alpha}=0.05$ au

(c) $g_{\alpha}=0.08$ au

###  IV.2 Large molecular metal complex

To demonstrate that our implementation is not restricted to atoms, but
applicable also to large molecules containing heavy elements, we calculate the
electron absorption spectrum of a mercury porphyrin inside a Fabry–Pérot
cavity. Porphyrins are a class of molecules with important biological roles as
well as technological applications and there is a large number of both
experimental [189] and theoretical [190, 191] works devoted to their study
including the modification of their properties in optical cavities. [139, 192,
193] Moreover, the absorption spectra of porphyrin complexes of group 12 atoms
including the Hg porphyrin considered here, were studied before at the four-
component relativistic level of theory by Fransson et al. [135] Note that
unlike the reference work [135] we do not explicitly consider point group
symmetry in our calculations thus performing an all-atom all-electron
calculation in the ground state.

The spectrum of the free Hg-substituted porphyrin is dominated by the strong
B-band at $3.39\,\mathrm{eV}$, while other notable features are the N- and
L-bands at 3.76 and $4.12\,\mathrm{eV}$, respectively with the Q-band being
very weak and not visible in this spectrum. All of these spectral lines arise
from doubly degenerate excited states with perpendicular transition dipole
moments, i.e. $y$ and $z$ for the molecule oriented such that its C4 symmetry
axis coincides with the $x$ axis. In the calculations of polaritonically
modified spectra we consider the Hg porphyrin complex embedded in a
Fabry–Pérot cavity with an effective single mode of frequency that is in
resonance with the B band excitation energy with progressively stronger
coupling constant $g_{\alpha}$. The resulting spectra together with the
spectrum of the uncoupled system are depicted in Figure 5. While for the
weakest coupling $g_{\alpha}=0.005$ (dark blue line), the dominating
modification of the spectrum due to the cavity is the splitting of one of the
two degenerate B lines into an upper and lower polariton, for the stringer
values of the coupling, other lines in the spectrum start to be affected as
well. Particularly for $g_{\alpha}=0.02$ we see breaking of the degeneracy in
the form of frequency and intensity shifts of other spectral lines including
the N and L lines. Therefore, besides demonstrating the readiness of our
method for the treatment of chemically relevant molecules, these results also
stress the importance of ab initio QED methods.

Figure 5: Comparison of relativistic (4c) electron absorption spectra (both
band spectra, i.e. spectral function $S$ on the left vertical axis, and line
spectra as oscillator strength $f$ on the right vertical axis, only lines with
$f>10^{-5}$ are plotted.) of Hg porphyrin complex in free space and in cavity
set to resonance with the porphyrin B line with different coupling strengths
$g_{\alpha}$ [au]. The Lorentzian spectra were obtained using broadening
parameter $\gamma=0.002$. Inset: the structure of the molecule.

##  V Conclusions

In this article we presented the first implementation of a relativistic ab
initio quantum electrodynamical method combining relativistic description of
electronic structure with the treatment of transverse photons as dynamical
variables in the form of linear response QEDFT based on the four-component
Dirac–Kohn–Sham Hamiltonian. We presented a detailed theory derivation of the
final electron–photon Hamiltonian starting from relativistic QED while keeping
track of all approximations made on the way. For the linear response equations
we showed how a spectral function can be calculated from a non-hermitian
Casida-like equation by considering both the left and right eigenvectors, as
well as showed how this equation can be transformed to the symmetric equation
often considered in ab initio QEDFT literature.

We applied the developed method to the calculation of polaritonically modified
spectra of heavy-element atoms with singlet–triplet transitions and
demonstrated that the intensity of the transition can be enhanced by the
coupling to the cavity as well showed the coupling between the singlet–singlet
and singlet–triplet transitions via the cavity modes that results in non-
resonant splitting of spectral lines as well the transfer of intensity between
them. Of specific interest is that the cavity even in the long-wavelength
limit can influence the SOC due to a fully relativistic (four-component)
treatment of the electrons. Our results showing the interplay of off-resonant
coupling and spin–orbit driven transition demonstrate the need for accurate ab
initio methods such as QEDFT to describe cavity modification of molecular
properties as well as help to guide the formulation of approximate models – in
the cases presented a minimal two-state model would be insufficient in
describing the cavity modified of spectra around the singlet–triplet
resonance. Moreover, the linear response QEDFT methodology developed with the
primary motivation of describing polaritonic chemistry in cavities allows to
address other situations when electrons are coupled to photons as well. Such
an example is the coupling to vacuum modes that presents a way of calculating
radiative lifetimes from first principles thereby resolving one of the
shortcomings of electronic-structure linear-response theories that add this
variable as an empirical parameter. Finally by addressing the spectrum of
mercury porphyrin in the cavity we demonstrated that our program is capable of
handling even large molecules of chemical interest.

Ab initio QED methods have emerged recently as a first principles approach for
the prediction and interpretation of phenomena occurring under strong coupling
of electronic and vibrational degrees of freedom of a material to the photonic
modes of an optical cavity. In this work, we expanded the toolbox of ab initio
QED methods further by bringing them into the relativistic domain. Considering
the advances in understanding also collective coupling effects [167, 159, 160]
and the here presented changes in the SOC, we think that this toolbox can help
to design cavity setups such that SOC and other relativistic effects can be
enhanced. This also conceptually is very interesting since it highlights that
the transverse photons, usually discarded when deriving relativistic quantum
chemistry equations from full QED, can be decisive for real atomic and
molecular systems. For instance, taking into account (approximately) the
continuum of modes can be used to study non-perturbatively very fundamental
effects of quantum field theories, such as the mass renormalization of the
electrons [7]. We also showed how relativistic description of electronic
structure of a material inside a photonic structure reveals the presence of a
new type of spin–orbit interaction mediated by cavity photons. This opens a
new avenue in the fields of cavity QED and polaritonic chemistry with a
potential for follow-up works focused on singlet–triplet transitions in
singlet fission systems, the exploration of the new SOC term in solid state
systems, or the development of effective Hamiltonians.

##  Appendix A Notations and conventions

We use the following definitions: $x^{\mu}$ for contravariant four-vector,
$x_{\mu}=g_{\mu\nu}x^{\nu}$ for covariant four-vector, $g_{\mu\nu}$ is the
metric tensor

$g_{\mu\nu}=g^{\mu\nu}=\left(\begin{array}[]{rrrr}1&0&0&0\\\ 0&-1&0&0\\\
0&0&-1&0\\\ 0&0&0&-1\end{array}\right),$

and $\alpha^{\mu}=\gamma^{0}\gamma^{\mu}$ are the Dirac matrices in their
standard representation

$\alpha^{0}=\left(\begin{array}[]{cc}\mathds{1}_{2\times 2}&0\\\
0&\mathds{1}_{2\times
2}\end{array}\right),\quad\bm{\alpha}=\left(\begin{array}[]{cc}0&\bm{\sigma}\\\
\bm{\sigma}&0\end{array}\right),$
$\gamma^{0}\equiv\beta=\left(\begin{array}[]{cc}\mathds{1}_{2\times 2}&0\\\
0&-\mathds{1}_{2\times
2}\end{array}\right),\quad\bm{\gamma}=\left(\begin{array}[]{cc}0&\bm{\sigma}\\\
-\bm{\sigma}&0\end{array}\right),$

and the $2\times 2$ Pauli matrices

$\quad\sigma_{x}=\left(\begin{array}[]{cc}0&1\\\
1&0\end{array}\right),\quad\sigma_{y}=\left(\begin{array}[]{cc}0&-i\\\
i&0\end{array}\right),\quad\sigma_{z}=\left(\begin{array}[]{rr}1&0\\\
0&-1\end{array}\right).$

Here $\mathds{1}_{2\times 2}$ is a unit $2\times 2$ matrix

$\mathds{1}_{2\times 2}=\left(\begin{array}[]{ll}1&0\\\
0&1\end{array}\right).$

We use the italic style $(r)$ for scalars, boldface $(\bf{r})$ for three
vectors, and Roman style $(\mathrm{r})$ for four-vectors and their components.
Four vectors have the form $\mathrm{p}=\left(p_{0},\bf{p}\right)$ and
coordinates in four space are also denoted as
$\mathrm{r}=\left(t,\bf{r}\right)$. The scalar product of four-vectors is
defined as $\mathrm{k}\mathrm{p}=k^{\mu}p_{\mu}=k^{0}p_{0}-\bf{k\cdot p}$. The
summing over repeated indices is assumed by default. Fock space vectors and
matrices are also denoted by Roman boldface with $\mathbf{1}$ and $\mathbf{0}$
denoting the unit and zero matrix, respectively.

##  Appendix B Solution of first order EOMs by undetermined coefficients

In the method of undetermined coefficients, we use the following ansatz for
the electronic equation

$\displaystyle d^{(1)}_{ai}(t)$
$\displaystyle=X_{ai}e^{-i\omega_{\mathrm{ext}}t}+Y^{*}_{ai}e^{i\omega_{\mathrm{ext}}t},$
(68a) $\displaystyle d^{(1)*}_{ai}(t)$
$\displaystyle=X^{*}_{ai}e^{i\omega_{\mathrm{ext}}t}+Y_{ai}e^{-i\omega_{\mathrm{ext}}t},$
(68b)

motivated by the inhomogeneity in Eq. (52). Similarly, we formulate an ansatz
for the photonic displacement coordinate in the form

$q^{(1)}_{\alpha}(t)=\tilde{M}_{\alpha}e^{-i\omega_{\mathrm{ext}}t}+\tilde{M}^{*}_{\alpha}e^{i\omega_{\mathrm{ext}}t},$
(69)

from which follows

$p^{(1)}_{\alpha}(t)=-\omega_{\mathrm{ext}}\tilde{M}_{\alpha}e^{-i\omega_{\mathrm{ext}}t}+\omega_{\mathrm{ext}}\tilde{M}^{*}_{\alpha}e^{i\omega_{\mathrm{ext}}t},$
(70)

and where the real nature of $q^{(1)}_{\alpha}(t)$ was taken into account.

We consider only a single inhomogeneity in the whole system of equations at a
time, i.e., the external field or the external current. The ansätze will then
contain its respective frequency (i.e. $\omega_{\mathrm{ext}}$ or
$\omega^{\prime}_{\mathrm{ext}}$). By inserting expressions from Eqs.
(68)-(70) into Eqs. (60) and selecting only the terms with time dependence
$\exp(-i\omega_{\mathrm{ext}}t)$ we obtain a set of coupled algebraic
equations

$\displaystyle\omega_{\mathrm{ext}}X_{ai}$
$\displaystyle=(A_{ai,bj}+\Delta_{ai,bj})X_{bj}+(B_{ai,bj}+\Delta^{\prime}_{ai,bj})Y_{bj}-L_{ai,\alpha}\tilde{M}_{\alpha}+P_{\mathrm{ai}},$
(71a) $\displaystyle-\omega_{\mathrm{ext}}Y_{ai}$
$\displaystyle=(A^{*}_{ai,bj}+\Delta^{*}_{ai,bj})Y_{bj}+(B^{*}_{ai,bj}+\Delta^{\prime*}_{ai,bj})X_{bj}-L^{*}_{ai,\alpha}\tilde{M}_{\alpha}+P^{*}_{\mathrm{ai}},$
(71b)
$\displaystyle(-\omega_{\mathrm{ext}}^{2}+\omega_{\alpha}^{2})\tilde{M}_{\alpha}$
$\displaystyle=Q^{\prime}_{jb,\alpha}X_{bj}+Q^{\prime*}_{jb,\alpha}Y_{bj}.$
(71c)

Unfortunately, the equation for $\tilde{M}_{\alpha}$ contains
$\omega_{\mathrm{ext}}^{2}$, a recurrence of the second order differential
equation (57) as we would obtain the system of equations (71) also by
considering Eq. (57) instead of Eqs. (58) and (59) for the photonic degrees of
freedom. Since we want to get an equation linear in $\omega_{\mathrm{ext}}$,
we remove the second power by employing the identity

$\frac{1}{\omega_{\mathrm{ext}}^{2}-\omega_{\alpha}^{2}}=\frac{1}{2\omega_{\alpha}}\left(\frac{1}{\omega_{\mathrm{ext}}+\omega_{\alpha}}-\frac{1}{\omega_{\mathrm{ext}}-\omega_{\alpha}}\right),$
(72)

which transforms Eq. (71c) into

$\tilde{M}_{\alpha}=\frac{1}{2\omega_{\alpha}}\left(\frac{Q^{\prime}_{jb,\alpha}X_{bj}+Q^{\prime*}_{jb,\alpha}Y_{bj}}{\omega_{\mathrm{ext}}+\omega_{\alpha}}-\frac{Q^{\prime}_{jb,\alpha}X_{bj}+Q^{\prime*}_{jb,\alpha}Y_{bj}}{\omega_{\mathrm{ext}}-\omega_{\alpha}}\right).$
(73)

We can now divide $\tilde{M}_{\alpha}$ into two components
$\tilde{M}_{\alpha}=M_{\alpha}+N_{\alpha}$ that satisfy Eq. (73) term by term

$\displaystyle M_{\alpha}$
$\displaystyle=-\frac{1}{2\omega_{\alpha}}\frac{Q^{\prime}_{jb,\alpha}X_{bj}+Q^{\prime*}_{jb,\alpha}Y_{bj}}{\omega_{\mathrm{ext}}-\omega_{\alpha}},$
(74a) $\displaystyle N_{\alpha}$
$\displaystyle=\frac{1}{2\omega_{\alpha}}\frac{Q^{\prime}_{jb,\alpha}X_{bj}+Q^{\prime*}_{jb,\alpha}Y_{bj}}{\omega_{\mathrm{ext}}+\omega_{\alpha}},$
(74b)

and both contain $\omega_{\mathrm{ext}}$ only in the first power. This way, we
can write four coupled algebraic equations linear in $\omega_{\mathrm{ext}}$
for the variables $X_{ai}$, $Y_{ai}$, $M_{\alpha}$, and $N_{\alpha}$ that
collected in a matrix give

$\begin{bmatrix}\begin{pmatrix}\mathbf{A}+\mathbf{\Delta}&\mathbf{B}+\mathbf{\Delta}^{\prime}&-\mathbf{L}&-\mathbf{L}\\\
\mathbf{A}^{*}+\mathbf{\Delta}^{*}&\mathbf{B}^{*}+\mathbf{\Delta}^{\prime*}&-\mathbf{L}^{*}&-\mathbf{L}^{*}\\\
-\mathbf{Q}&-\mathbf{Q}^{*}&\bm{\omega}&\mathbf{0}\\\
-\mathbf{Q}&-\mathbf{Q}^{*}&\mathbf{0}&\bm{\omega}\\\
\end{pmatrix}-\omega_{\mathrm{ext}}\begin{pmatrix}\mathbf{1}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&-\mathbf{1}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{1}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}&-\mathbf{1}\\\
\end{pmatrix}\end{bmatrix}\begin{pmatrix}\mathbf{X}\\\ \mathbf{Y}\\\
\mathbf{M}\\\ \mathbf{N}\\\ \end{pmatrix}=\begin{pmatrix}\mathbf{P}\\\
\mathbf{P}^{*}\\\ \mathbf{0}\\\ \mathbf{0}\end{pmatrix},$ (75)

where we defined

$Q_{jb,\alpha}=\frac{1}{2\omega_{\alpha}}Q^{\prime}_{jb,\alpha}=\frac{g_{\alpha}e\omega_{\alpha}}{\sqrt{\hbar}}\frac{1}{2\omega_{\alpha}}(\mathbf{r}_{jb}\cdot\bm{\varepsilon}_{\alpha})=\frac{1}{2}\frac{g_{\alpha}e}{\sqrt{\hbar}}(\mathbf{r}_{jb}\cdot\bm{\varepsilon}_{\alpha}),$
(76)

and $\bm{\omega}$ is a diagonal matrix with elements $\omega_{\alpha}$ on the
diagonal. An analogous equation can be derived for the case of external
current with the right-hand side having the form $(\mathbf{0}\ \mathbf{0}\
\mathbf{J}\ \mathbf{J}^{*})^{\mathrm{T}}$. This is the Sternheimer equation
describing the response of a molecular system in a cavity to an external field
of frequency $\omega_{\mathrm{ext}}$ coupled to the electronic subsystem via
operator $\mathbf{P}$. The algebraic form of the homogeneous equation to the
system (60) can now be obtained from Eq. (75) by discarding the right-hand
side with the frequency of the external field gaining the meaning of a
resonance frequency of the system, i.e. $n$-th excitation energy $\Omega_{n}$,
leading to Eq. (62).

##  Appendix C Solution of the matrix differential equation and the
interpretation of photon terms

It is interesting to note that an alternative approach to Eq. (62) exists
based on a direct solution of the coupled system of differential equations
(60). This system corresponds to a matrix differential equation in the form
$i\dot{\mathbf{x}}(t)=\bm{\mathcal{M}}\mathbf{x}(t)+\mathbf{b}(t)$, where

$\displaystyle\mathbf{x}(t)$ $\displaystyle=(\mathbf{d}^{(1)}(t),\
\mathbf{d}^{(1)*}(t),\ \mathbf{p}^{(1)}(t),\
\mathbf{q}^{(1)}(t))^{\mathrm{T}},$ (77a) $\displaystyle\bm{\mathcal{M}}$
$\displaystyle=\begin{pmatrix}\mathbf{A}+\mathbf{\Delta}&\mathbf{B}+\mathbf{\Delta}^{\prime}&\mathbf{0}&-\mathbf{L}\\\
-\mathbf{B}^{*}-\mathbf{\Delta}^{\prime*}&-\mathbf{A}^{*}-\mathbf{\Delta}^{*}&\mathbf{0}&\mathbf{L}^{*}\\\
\mathbf{Q}^{\prime*}&\mathbf{Q}^{\prime}&\mathbf{0}&-\bm{\omega}^{2}\\\
\mathbf{0}&\mathbf{0}&-\mathbf{1}&\mathbf{0}\end{pmatrix},$ (77b)
$\displaystyle\mathbf{b}(t)$
$\displaystyle=(\mathbf{P}e^{-i\omega_{\mathrm{ext}}t}+\mathbf{P}^{*}e^{i\omega_{\mathrm{ext}}t},\
-\mathbf{P}^{*}e^{i\omega_{\mathrm{ext}}t}-\mathbf{P}e^{-i\omega_{\mathrm{ext}}t},\
\mathbf{J}e^{-i\omega^{\prime}_{\mathrm{ext}}t}-\mathbf{J}^{*}e^{i\omega^{\prime}_{\mathrm{ext}}t},\
\mathbf{0})^{\mathrm{T}}.$ (77c)

The homogeneous equation takes the form
$i\dot{\mathbf{x}}(t)=\bm{\mathcal{M}}\mathbf{x}(t)$ and its general solution
is $\mathbf{x}(t)=\sum_{n}c_{n}e^{-i\Omega_{n}t}\mathbf{u}_{n}$, where
$\mathbf{u}_{n}$ are the eigenvectors of matrix $\bm{\mathcal{M}}$,
$\Omega_{n}$ their corresponding eigenvalues and $c_{n}$ constant
coefficients. The solution of the homogeneous equation thus boils down to
finding the eigenvectors and eigenvalues of the matrix $\bm{\mathcal{M}}$
defined in Eq. (77b). Since two matrices related by a similarity
transformation have the same set of eigenvalues, one can alternatively solve
the eigenproblem for its similar matrix
$\bm{\mathcal{M}}^{\prime}=\mathbf{S}\bm{\mathcal{M}}\mathbf{S}^{-1}$ with
eigenvectors $\mathbf{u}^{\prime}=\mathbf{S}\mathbf{u}$, if the new matrix has
some more favorable properties. For example, matrix $\bm{\mathcal{M}}$ defined
in Eq. (77b) can be transformed into a diagonally dominant matrix by

$\begin{pmatrix}\mathbf{A}+\mathbf{\Delta}&\mathbf{B}+\mathbf{\Delta}^{\prime}&-\mathbf{L}&-\mathbf{L}\\\
-\mathbf{B}^{*}-\mathbf{\Delta}^{\prime*}&-\mathbf{A}^{*}-\mathbf{\Delta}^{*}&\mathbf{L}^{*}&\mathbf{L}^{*}\\\
-\mathbf{Q}^{*}&-\mathbf{Q}&\bm{\omega}&\mathbf{0}\\\
\mathbf{Q}^{*}&\mathbf{Q}&\mathbf{0}&-\bm{\omega}\\\
\end{pmatrix}=\begin{pmatrix}\mathbf{1}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{1}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&-\bm{\frac{1}{2\omega}}&\mathbf{\frac{1}{2}}\\\
\mathbf{0}&\mathbf{0}&\bm{\frac{1}{2\omega}}&\mathbf{\frac{1}{2}}\\\
\end{pmatrix}\begin{pmatrix}\mathbf{A}+\mathbf{\Delta}&\mathbf{B}+\mathbf{\Delta}^{\prime}&\mathbf{0}&-\mathbf{L}\\\
-\mathbf{B}^{*}-\mathbf{\Delta}^{\prime*}&-\mathbf{A}^{*}-\mathbf{\Delta}^{*}&\mathbf{0}&\mathbf{L}^{*}\\\
\mathbf{Q}^{\prime*}&\mathbf{Q}^{\prime}&\mathbf{0}&-\bm{\omega}^{2}\\\
\mathbf{0}&\mathbf{0}&\mathbf{-1}&\mathbf{0}\\\
\end{pmatrix}\begin{pmatrix}\mathbf{1}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{1}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&-\bm{\omega}&\bm{\omega}\\\
\mathbf{0}&\mathbf{0}&\mathbf{1}&\mathbf{1}\\\ \end{pmatrix},$ (78)

with the same definition of $\mathbf{Q}$ as above and matrices
$\pm\bm{\frac{1}{2\omega}}$, $\mathbf{\frac{1}{2}}$, and $\pm\bm{\omega}$
corresponding to diagonal matrices with $\pm\frac{1}{2\omega_{\alpha}}$,
$\frac{1}{2}$, and $\pm\omega_{\alpha}$, respectively, on the diagonal. This
transformation recovers the matrix in Eq. (62) (after multiplication with the
matrix $\mathrm{diag}(\mathbf{1},-\mathbf{1},\mathbf{1},-\mathbf{1})$ on the
right-hand side). Moreover, the similarity transformation also shows a
relationship between the new photonic variables $M_{\alpha}$ and $N_{\alpha}$,
and the old $\mathbf{p}^{(1)}(t)$ and $\mathbf{q}^{(1)}(t)$ in the form

$\displaystyle M_{\alpha}$
$\displaystyle=-\frac{1}{2\omega_{\alpha}}p^{(1)}_{\alpha}+\frac{1}{2}q^{(1)}_{\alpha}$
(79a) $\displaystyle N_{\alpha}$
$\displaystyle=\frac{1}{2\omega_{\alpha}}p^{(1)}_{\alpha}+\frac{1}{2}q^{(1)}_{\alpha},$
(79b)

which is analogous but not yet identical to the definition of creation and
annihilation operators in Eqs. (20) and (19), suggesting an interpretation of
$M_{\alpha}$ as a creation and $N_{\alpha}$ as an annihilation amplitude,
respectively.

Furthermore, we can perform yet another similarity transformation to
symmetrize the matrix $\bm{\mathcal{M}}^{\prime}$ defined in Eq. (78),

$\begin{pmatrix}\mathbf{A}+\mathbf{\Delta}&\mathbf{B}+\mathbf{\Delta}^{\prime}&-\mathbf{g}&-\mathbf{g}\\\
-\mathbf{B}^{*}-\mathbf{\Delta}^{\prime*}&-\mathbf{A}^{*}-\mathbf{\Delta}^{*}&\mathbf{g}^{*}&\mathbf{g}^{*}\\\
-\mathbf{g}^{*}&-\mathbf{g}&\bm{\omega}&\mathbf{0}\\\
\mathbf{g}^{*}&\mathbf{g}&\mathbf{0}&-\bm{\omega}\\\
\end{pmatrix}=\begin{pmatrix}\mathbf{1}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{1}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\bm{\sqrt{2\omega}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}&\bm{\sqrt{2\omega}}\\\
\end{pmatrix}\begin{pmatrix}\mathbf{A}+\mathbf{\Delta}&\mathbf{B}+\mathbf{\Delta}^{\prime}&-\mathbf{L}&-\mathbf{L}\\\
-\mathbf{B}^{*}-\mathbf{\Delta}^{\prime*}&-\mathbf{A}^{*}-\mathbf{\Delta}^{*}&\mathbf{L}^{*}&\mathbf{L}^{*}\\\
-\mathbf{Q}^{*}&-\mathbf{Q}&\bm{\omega}&\mathbf{0}\\\
\mathbf{Q}^{*}&\mathbf{Q}&\mathbf{0}&-\bm{\omega}\\\
\end{pmatrix}\begin{pmatrix}\mathbf{1}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{1}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\bm{\frac{1}{\sqrt{2\omega}}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}&\bm{\frac{1}{\sqrt{2\omega}}}\\\
\end{pmatrix},$ (80)

where the photon–electron and the electron–photon coupling has the form

$g_{ai,\alpha}=\sqrt{\frac{\hbar\omega_{\alpha}}{2}}\,g_{\alpha}\,\vec{\mu}_{ai}\cdot\vec{\varepsilon}_{\alpha},$
(81)

and matrices $\bm{\sqrt{2\omega}}$ and $\bm{\frac{1}{\sqrt{2\omega}}}$ are
diagonal matrices with $\sqrt{2\omega_{\alpha}}$ and
$\frac{1}{\sqrt{2\omega_{\alpha}}}$, respectively, on the diagonal.
Symmetrization in Eq. (80) recovered the matrix appearing in the Casida
equation solved in some other works on linear response QEDFT. [70, 71]
Furthermore, the symmetrizing similarity transform gives

$\displaystyle M_{\alpha}$
$\displaystyle=\sqrt{\frac{\omega_{\alpha}}{2}}\left(-\frac{1}{\omega_{\alpha}}p^{(1)}_{\alpha}+q^{(1)}_{\alpha}\right)$
(82a) $\displaystyle N_{\alpha}$
$\displaystyle=\sqrt{\frac{\omega_{\alpha}}{2}}\left(\frac{1}{\omega_{\alpha}}p^{(1)}_{\alpha}+q^{(1)}_{\alpha}\right),$
(82b)

which exactly corresponds to the definition of the creation and annihilation
operators in Eqs. (20) and (19) giving $M_{\alpha}$ and $N_{\alpha}$ the
meaning of creation and annihilation amplitudes, respectively.

##  Appendix D Derivation of the expression for polarizability using left and
right eigenvectors

To compactify the following discussion, let us adopt the notation common in
quantum chemistry where the Sternheimer equation (75) is written term-by-term
as
$[\mathbf{E}^{[2]}-\omega_{\mathrm{ext}}\mathbf{S}^{[2]}]\mathbf{Z}=\mathbf{G}$.
and Eq. (63) takes the form
$\bm{\alpha}(\omega_{\mathrm{ext}})=-\mathbf{G}^{\dagger}\left[\mathbf{E}^{[2]}-\omega_{\mathrm{ext}}\mathbf{S}^{[2]}\right]^{-1}\mathbf{G}$.
The right and left eigenvectors are defined as

$\displaystyle\mathbf{S}^{[2]}\mathbf{E}^{[2]}\mathbf{Z}^{\mathrm{R}}$
$\displaystyle=\mathbf{Z}^{\mathrm{R}}\mathbf{\Omega},$ (83a)
$\displaystyle\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{S}^{[2]}\mathbf{E}^{[2]}$
$\displaystyle=\mathbf{\Omega}\mathbf{Z}^{\mathrm{L\dagger}},$ (83b)

where matrix $\mathbf{\Omega}$ is a diagonal matrix containing all the
eigenvalues and the columns of matrices $\mathbf{Z}^{\mathrm{R}}$ and
$\mathbf{Z}^{\mathrm{L}}$ contain all the right and left eigenvectors,
respectively. Moreover, we made use of the fact that
$\mathbf{S}^{[2],-1}=\mathbf{S}^{[2]}$. Left and right eigenvectors
corresponding to different eigenvalues are orthogonal, and can be transformed
to be so for degenerate subspaces, meaning that their matrices satisfy
$\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{Z}^{\mathrm{R}}=\mathbf{1}$. Let us
note that the definition of $\mathbf{Z}^{\mathrm{L\dagger}}$ as a left
eigenvector of matrix $\mathbf{S}^{[2]}\mathbf{E}^{[2]}$ is due to the
computer implementation that solves the eigenproblem for matrix
$\mathbf{S}^{[2]}\mathbf{E}^{[2]}$ instead of the generalized eigenproblem.
The alternative formulation of left eigenvectors stemming from the generalized
eigenproblem
$\tilde{\mathbf{Z}}^{\mathrm{L\dagger}}\mathbf{E}^{[2]}=\mathbf{\Omega}\tilde{\mathbf{Z}}^{\mathrm{L\dagger}}\mathbf{S}^{[2]}$
is related to ours via
$\tilde{\mathbf{Z}}^{\mathrm{L\dagger}}\mathbf{S}^{[2]}=\mathbf{Z}^{\mathrm{L\dagger}}$
and the eigenvectors satisfy the condition
$\tilde{\mathbf{Z}}^{\mathrm{L\dagger}}\mathbf{S}^{[2]}\mathbf{Z}^{\mathrm{R}}=\mathbf{1}$.
From Eqs. (83) it also follows that
$\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{E}^{[2]}\mathbf{Z}^{\mathrm{R}}=\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{S}^{[2]}\mathbf{Z}^{\mathrm{R}}\mathbf{\Omega}$.

The expression for the polarizability can be worked out as

$\displaystyle\begin{split}\bm{\alpha}(\omega)&=-\mathbf{G}^{\dagger}\left[\mathbf{E}^{[2]}-\omega\mathbf{S}^{[2]}\right]^{-1}\mathbf{G}\\\
&=-\mathbf{G}^{\dagger}\mathbf{Z}^{\mathrm{R}}\mathbf{Z}^{\mathrm{L\dagger}}\left[\mathbf{E}^{[2]}-\omega\mathbf{S}^{[2]}\right]^{-1}\mathbf{Z}^{\mathrm{R}}\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{G}\\\
&=-\left(\mathbf{G}^{\dagger}\mathbf{Z}^{\mathrm{R}}\right)\left(\mathbf{Z}^{\mathrm{L\dagger}}\left[\mathbf{E}^{[2]}-\omega\mathbf{S}^{[2]}\right]^{-1}\mathbf{Z}^{\mathrm{R}}\right)\left(\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{G}\right)\\\
&=-\left(\mathbf{G}^{\dagger}\mathbf{Z}^{\mathrm{R}}\right)\left(\mathbf{Z}^{\mathrm{R,-1}}\left[\mathbf{E}^{[2]}-\omega\mathbf{S}^{[2]}\right]\mathbf{Z}^{\mathrm{L\dagger,-1}}\right)^{-1}\left(\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{G}\right)\\\
&=-\left(\mathbf{G}^{\dagger}\mathbf{Z}^{\mathrm{R}}\right)\left(\mathbf{Z}^{\mathrm{L\dagger}}\left[\mathbf{E}^{[2]}-\omega\mathbf{S}^{[2]}\right]\mathbf{Z}^{\mathrm{R}}\right)^{-1}\left(\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{G}\right)\\\
&=-\left(\mathbf{G}^{\dagger}\mathbf{Z}^{\mathrm{R}}\right)\left(\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{E}^{[2]}\mathbf{Z}^{\mathrm{R}}-\omega\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{S}^{[2]}\mathbf{Z}^{\mathrm{R}}\right)^{-1}\left(\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{G}\right)\\\
&=-\left(\mathbf{G}^{\dagger}\mathbf{Z}^{\mathrm{R}}\right)\left(\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{S}^{[2]}\mathbf{Z}^{\mathrm{R}}\mathbf{\Omega}-\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{S}^{[2]}\mathbf{Z}^{\mathrm{R}}\omega\right)^{-1}\left(\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{G}\right)\\\
&=-\left(\mathbf{G}^{\dagger}\mathbf{Z}^{\mathrm{R}}\right)\left(\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{S}^{[2]}\mathbf{Z}^{\mathrm{R}}\left[\mathbf{\Omega}-\omega\mathbf{1}\right]\right)^{-1}\left(\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{G}\right)\\\
&=-\left(\mathbf{G}^{\dagger}\mathbf{Z}^{\mathrm{R}}\right)\left(\left[\mathbf{\Omega}-\omega\mathbf{1}\right]^{-1}\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{S}^{[2]}\mathbf{Z}^{\mathrm{R}}\right)\left(\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{G}\right)\\\
&=-\left(\mathbf{G}^{\dagger}\mathbf{Z}^{\mathrm{R}}\right)\left[\mathbf{\Omega}-\omega\mathbf{1}\right]^{-1}\left(\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{S}^{[2]}\mathbf{G}\right),\\\
\end{split}$ (84)

where we renamed $\omega_{\mathrm{ext}}$ to $\omega$ both for clarity as well
as to indicate its meaning as the frequency at which the frequency-dependent
polarizability is evaluated. The matrix $\mathbf{\Omega}-\omega\mathbf{1}$ is
diagonal and can be easily inverted. However, the formula can be further
simplified by taking into account the fact that the eigenvalues come in
positive–negative pairs, i.e. the matrix $\mathbf{\Omega}$ has the structure

$\mathbf{\Omega}=\begin{pmatrix}\mathbf{\Omega}^{+}&\mathbf{0}\\\
\mathbf{0}&-\mathbf{\Omega}^{+}\end{pmatrix},$ (85)

where $\mathbf{\Omega}^{+}$ is a diagonal matrix of the positive eigenvalues
labeled from 1 to $n$. Therefore, the expression for the polarizability tensor
becomes

$\bm{\alpha}(\omega)=-\sum_{n}\left[\frac{t^{\mathrm{R}}_{n}t^{\mathrm{L*}}_{n}}{\Omega_{n}-\omega}-\frac{t^{\mathrm{R}}_{-n}t^{\mathrm{L*}}_{-n}}{-\Omega_{n}-\omega}\right],$
(86)

where variables labeled by $n$ and $-n$ belong to excitation energies
$\Omega_{n}$ and $-\Omega_{n}$, respectively, and the right and left
transition dipole moments are defined as

$\displaystyle t^{\mathrm{R}}_{n}$
$\displaystyle=\mathbf{G}^{\dagger}\mathbf{Z}^{\mathrm{R}}$
$\displaystyle=\mathbf{P}^{\dagger}\mathbf{X}^{\mathrm{R}}_{n}+\mathbf{P}\mathbf{Y}^{\mathrm{R}}_{n},$
(87a) $\displaystyle t^{\mathrm{L}}_{n}$
$\displaystyle=\mathbf{Z}^{\mathrm{L\dagger}}\mathbf{S}^{[2]}\mathbf{G}$
$\displaystyle=\mathbf{P}^{\dagger}\mathbf{X}^{\mathrm{R}}_{n}-\mathbf{P}\mathbf{Y}^{\mathrm{R}}_{n}.$
(87b)

The fact that only the electronic parts $\mathbf{X}^{\mathrm{R/L}}_{n}$ and
$\mathbf{Y}^{\mathrm{R/L}}_{n}$ of the transition vector
$\mathbf{Z}^{\mathrm{R/L}}_{n}=(\mathbf{X}^{\mathrm{R/L}}_{n},\
\mathbf{Y}^{\mathrm{R/L}}_{n},\ \mathbf{M}^{\mathrm{R/L}}_{n},\
\mathbf{N}^{\mathrm{R/L}}_{n})^{\mathrm{T}}$ are used to calculate the
transition dipole moment in Eqs. (65) is the consequence of the fact that the
electric dipole moment is an electron-only property described by
$\mathbf{G}=(\mathbf{P},\ \mathbf{P}^{*},\,\mathbf{0},\
\mathbf{0})^{\mathrm{T}}$. The eigenvector corresponding to an eigenvalue with
the opposite sign has
$\mathbf{X}^{\mathrm{R/L}}_{-n}=\mathbf{Y}^{\mathrm{R/L*}}_{n}$ and
$\mathbf{Y}^{\mathrm{R/L}}_{-n}=\mathbf{X}^{\mathrm{R/L*}}_{n}$, therefore
$t^{\mathrm{R}}_{-n}=t^{\mathrm{R*}}_{n}$ and
$t^{\mathrm{L}}_{-n}=t^{\mathrm{L*}}_{n}$. After inserting these expressions
into Eq. (86), we obtain the final expression for the frequency-dependent
polarizability calculated from the distinct right and left eigenvectors

$\bm{\alpha}(\omega)=\sum_{n}\left[\frac{t^{\mathrm{R*}}_{n}t^{\mathrm{L}}_{n}}{\Omega_{n}+\omega}-\frac{t^{\mathrm{R}}_{n}t^{\mathrm{L*}}_{n}}{\omega-\Omega_{n}}\right].$
(88)

## Supplemental Material

Additional spectra, JC model definition, geometry of the Hg@porphyrin.

## ACKNOWLEDGMENTS

LK acknowledges support by the Research Council of Norway through its Centres
of Excellence scheme, project no. 262695, and its Mobility Grant scheme,
project no. 314814. LK would like to thank Stanislav Komorovsky, Johannes
Flick, and Davis Welakuh for helpful discussions. VK would like to thank
Simone Latini for the help in the preparation figure of the cavity. The
calculations were performed on resources provided by Sigma2—the National
Infrastructure for High Performance Computing and Data Storage in Norway,
Grant Nos. NN4654K and NN14654K.

## References

  * [1] Thomas W. Ebbesen.  Hybrid light–matter states in a molecular and material science perspective.  Acc. Chem. Res., 49(11):2403–2412, 2016.  PMID: 27779846.
  * [2] Michael Ruggenthaler, Nicolas Tancogne-Dejean, Johannes Flick, Heiko Appel, and Angel Rubio.  From a quantum-electrodynamical light–matter description to novel spectroscopies.  Nat. Rev. Chem., 2(3):1–16, 2018.
  * [3] Anton Frisk Kockum, Adam Miranowicz, Simone De Liberato, Salvatore Savasta, and Franco Nori.  Ultrastrong coupling between light and matter.  Nature Reviews Physics, 1(1):19–40, 2019.
  * [4] Dmitri N Basov, Ana Asenjo-Garcia, P James Schuck, Xiaoyang Zhu, and Angel Rubio.  Polariton panorama.  Nanophotonics, 10(1):549–577, 2020.
  * [5] 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. USA, 112:15285, 2015.
  * [6] F. Schlawin, D. M. Kennes, and M. A. Sentef.  Cavity quantum materials.  Appl. Phys. Rev., 9(1):011312, 02 2022.
  * [7] Michael Ruggenthaler, Dominik Sidler, and Angel Rubio.  Understanding polaritonic chemistry from ab initio quantum electrodynamics.  Chemical Reviews, 0(0):null, 0.  PMID: 37729114.
  * [8] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg.  Photons and Atoms-Introduction to Quantum Electrodynamics.  Wiley-VCH, Weinheim, 1997.
  * [9] Johannes Feist, Javier Galego, and Francisco J. Garcia-Vidal.  Polaritonic chemistry with organic molecules.  ACS Photonics, 5(1):205–216, 2018.
  * [10] Tao E. Li, Abraham Nitzan, and Joseph E. Subotnik.  On the origin of ground-state vacuum-field catalysis: Equilibrium consideration.  J. Chem. Phys., 152(23):234107, 06 2020.
  * [11] 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.  J. Chem. Phys., 156(23):230901, 06 2022.
  * [12] Arkajit Mandal, Michael AD Taylor, Braden M Weight, Eric R Koessler, Xinyang Li, and Pengfei Huo.  Theoretical advances in polariton chemistry and molecular cavity quantum electrodynamics.  Chemical Reviews, 123(16):9786–9879, 2023.
  * [13] Blake S Simpkins, Adam D Dunkelberger, and Igor Vurgaftman.  Control, modulation, and analytical descriptions of vibrational strong coupling.  Chem. Rev., 123(8):5020–5048, 2023.
  * [14] D. Basov, R. Averitt, and D. Hseeh.  Towards properties on demand in quantum materials.  Nat. Mater., 16:1077–1088, 2017.
  * [15] M. Buzzi, M. Först, R. Mankowsky, and A Cavalleri.  Probing dynamics in quantum materials with femtosecond X-rays.  Nat. Rev. Mater., 3:299–311, 2018.
  * [16] Xiao Wang, Enrico Ronca, and Michael A. Sentef.  Cavity quantum electrodynamical chern insulator: Towards light-induced quantized anomalous hall effect in graphene.  Phys. Rev. B, 99:235156, Jun 2019.
  * [17] Martin Kiffner, Jonathan R. Coulthard, Frank Schlawin, Arzhang Ardavan, and Dieter Jaksch.  Manipulating quantum materials with quantum light.  Phys. Rev. B, 99:085116, Feb 2019.
  * [18] Jiajun Li, Denis Golez, Giacomo Mazza, Andrew J. Millis, Antoine Georges, and Martin Eckstein.  Electromagnetic coupling in tight-binding models for strongly correlated light and matter.  Phys. Rev. B, 101:205140, May 2020.
  * [19] Yuto Ashida, Ataç İmamoğlu, and Eugene Demler.  Cavity Quantum Electrodynamics at Arbitrary Light-Matter Coupling Strengths.  Phys. Rev. Lett., 126:153603, Apr 2021.
  * [20] James A Hutchison, Tal Schwartz, Cyriaque Genet, Eloïse Devaux, and Thomas W Ebbesen.  Modifying chemical landscapes by coupling to vacuum fields.  Angew. Chem. Int. Ed., 51(7):1592–1596, 2012.
  * [21] 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.  Angew. Chem. Int. Ed., 55(38):11462–11466, 2016.
  * [22] Xinyang Li, Arkajit Mandal, and Pengfei Huo.  Cavity frequency-dependent theory for vibrational polariton chemistry.  Nat. Commun., 12(1):1315, 2021.
  * [23] Francisco J Garcia-Vidal, Cristiano Ciuti, and Thomas W Ebbesen.  Manipulating matter by strong coupling to vacuum fields.  Science, 373(6551):eabd0336, 2021.
  * [24] 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.  Nat. Commun., 13:7817, 2022.
  * [25] Wonmi Ahn, Johan F Triana, Felipe Recabal, Felipe Herrera, and Blake S Simpkins.  Modification of ground-state chemical reactivity via light–matter coherence in infrared cavities.  Science, 380(6650):1165–1168, 2023.
  * [26] 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(3):992–1005, 2018.
  * [27] Simone Latini, Dongbin Shin, Shunsuke A. Sato, Christian Schäfer, Umberto De Giovannini, Hannes Hübener, and Angel Rubio.  The ferroelectric photo ground state of SrTiO3: Cavity materials engineering.  Proc. Natl. Acad. Sci. U.S.A., 118(31):e2105618118, 2021.
  * [28] Vasil Rokaj, Simeon I. Mistakidis, and H. R. Sadeghpour.  Cavity induced collective behavior in the polaritonic ground state.  SciPost Phys., 14:167, 2023.
  * [29] Emil Vinas Boström, Adithya Sriram, Martin Claassen, and Angel Rubio.  Controlling the magnetic state of the proximate quantum spin liquid $\alpha$-RuCl3 with an optical cavity.  arXiv:2211.07247, 2022.
  * [30] 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.  Nat. Commun., 5:5561, 2014.
  * [31] E. Orgiu, J. George, J. A. Hutchison, E. Devaux, J. F. Dayen, B. Doudin, F. Stellacci, C. Genet, J. Schachenmayer, C. Genes, G. Pupillo, P. Samori, and T. W. Ebbesen.  Conductivity in organic semiconductors hybridized with the vacuum field.  Nat. Mater, 14:1123–1129, 2015.
  * [32] Johannes Schachenmayer, Claudiu Genes, Edoardo Tignone, and Guido Pupillo.  Cavity-enhanced transport of excitons.  Phys. Rev. Lett., 114(19):196403, 2015.
  * [33] 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.  Angew. Chem. Int. Ed., 55(21):6202–6206, 2016.
  * [34] Johannes Feist and Francisco J. Garcia-Vidal.  Extraordinary exciton conductance induced by strong coupling.  Phys. Rev. Lett., 114:196402, May 2015.
  * [35] Kati Stranius, Manuel Hertzog, and Karl Börjesson.  Selective manipulation of electronically excited states through strong light–matter interactions.  Nat. Commun., 9(1):2273, 2018.
  * [36] Ovidiu Cotleţ, Sina Zeytinoǧlu, Manfred Sigrist, Eugene Demler, and Ata ç Imamoǧlu.  Superconductivity and other collective phenomena in a hybrid bose-fermi mixture formed by a polariton condensate and an electron system in two dimensions.  Phys. Rev. B, 93:054510, Feb 2016.
  * [37] M. A. Sentef, M. Ruggenthaler, and Rubio A.  Cavity quantum-electrodynamical polaritonically enhanced electron- phonon coupling and its influence on superconductivity.  Sci. Adv., 4:eaau6969, 2018.
  * [38] Jonathan B. Curtis, Zachary M. Raines, Andrew A. Allocca, Mohammad Hafezi, and Victor M. Galitski.  Cavity quantum eliashberg enhancement of superconductivity.  Phys. Rev. Lett., 122:167002, Apr 2019.
  * [39] Frank Schlawin, Andrea Cavalleri, and Dieter Jaksch.  Cavity-mediated electron-photon superconductivity.  Phys. Rev. Lett., 122:133602, Apr 2019.
  * [40] Simone Latini, Enrico Ronca, Umberto De Giovannini, Hannes Hübener, and Angel Rubio.  Cavity control of excitons in two-dimensional materials.  Nano Lett., 19(6):3473–3479, 2019.  PMID: 31046291.
  * [41] Jesper Levinsen, Guangyao Li, and Meera M. Parish.  Microscopic description of exciton-polaritons in microcavities.  Phys. Rev. Res., 1:033120, Nov 2019.
  * [42] Michael Förg, Léo Colombier, Robin K. Patel, Jessica Lindlau, Aditya D. Mohite, Hisato Yamaguchi, Mikhail M. Glazov, David Hunger, and Alexander Högele.  Cavity-control of interlayer excitons in van der waals heterostructures.  Nat. Commun., 10(1):3697, Aug 2019.
  * [43] 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 Le Si Dang.  Bose–einstein condensation of exciton polaritons.  Nature, 443(7110):409–414, Sep 2006.
  * [44] Jonathan Keeling and Stéphane Kéna-Cohen.  Bose–Einstein Condensation of Exciton-Polaritons in Organic Microcavities.  Annu. Rev. Phys. Chem., 71(1):435–459, 2020.  PMID: 32126177.
  * [45] David Hagenmüller, Simone De Liberato, and Cristiano Ciuti.  Ultrastrong coupling between a cavity resonator and the cyclotron transition of a two-dimensional electron gas in the case of an integer filling factor.  Phys. Rev. B, 81:235303, Jun 2010.
  * [46] G. Scalari, C. Maissen, D. Turčinková, D. Hagenmüller, S. De Liberato, C. Ciuti, C. Reichl, D. Schuh, W. Wegscheider, M. Beck, and J. Faist.  Ultrastrong coupling of the cyclotron transition of a 2d electron gas to a thz metamaterial.  Science, 335(6074):1323–1326, 2012.
  * [47] Xinwei Li, Motoaki Bamba, Qi Zhang, Saeed Fallahi, Geoff C. Gardner, Weilu Gao, Minhan Lou, Katsumasa Yoshioka, Michael J. Manfra, and Junichiro Kono.  Vacuum bloch–siegert shift in landau polaritons with ultra-high cooperativity.  Nat. Photon., 12(6):324–329, Jun 2018.
  * [48] Janine Keller, Giacomo Scalari, Felice Appugliese, Shima Rajabali, Mattias Beck, Johannes Haase, Christian A. Lehner, Werner Wegscheider, Michele Failla, Maksym Myronov, David R. Leadley, James Lloyd-Hughes, Pierre Nataf, and Jérôme Faist.  Landau polaritons in highly nonparabolic two-dimensional gases in the ultrastrong coupling regime.  Phys. Rev. B, 101:075301, Feb 2020.
  * [49] Vasil Rokaj, Markus Penz, Michael A Sentef, Michael Ruggenthaler, and Angel Rubio.  Polaritonic Hofstadter butterfly and cavity control of the quantized Hall conductance.  Phys. Rev. B, 105(20):205424, 2022.
  * [50] Sylvain Ravets, Patrick Knüppel, Stefan Faelt, Ovidiu Cotlet, Martin Kroner, Werner Wegscheider, and Atac Imamoglu.  Polaron polaritons in the integer and fractional quantum hall regimes.  Phys. Rev. Lett., 120:057401, Jan 2018.
  * [51] Stephan Smolka, Wolf Wuester, Florian Haupt, Stefan Faelt, Werner Wegscheider, and Ataç Imamoglu.  Cavity quantum electrodynamics with many-body states of a two-dimensional electron gas.  Science, 346(6207):332–335, 2014.
  * [52] Gian L. Paravicini-Bagliani, Felice Appugliese, Eli Richter, Federico Valmorra, Janine Keller, Mattias Beck, Nicola Bartolo, Clemens Rössler, Thomas Ihn, Klaus Ensslin, Cristiano Ciuti, Giacomo Scalari, and Jérôme Faist.  Magneto-transport controlled by landau polariton states.  Nat. Phys., 15(2):186–190, Feb 2019.
  * [53] Jan Petersen, Jürgen Volz, and Arno Rauschenbeutel.  Chiral nanophotonic waveguide interface based on spin-orbit interaction of light.  Science, 346(6205):67–71, 2014.
  * [54] Peter Lodahl, Sahand Mahmoodian, Søren Stobbe, Arno Rauschenbeutel, Philipp Schneeweiss, Jürgen Volz, Hannes Pichler, and Peter Zoller.  Chiral quantum optics.  Nature, 541(7638):473–480, Jan 2017.
  * [55] Fan Zhang, Juanjuan Ren, Lingxiao Shan, Xueke Duan, Yan Li, Tiancai Zhang, Qihuang Gong, and Ying Gu.  Chiral cavity quantum electrodynamics with coupled nanophotonic structures.  Phys. Rev. A, 100:053841, Nov 2019.
  * [56] Hannes Hübener, Umberto De Giovannini, Christian Schäfer, Johan Andberger, Michael Ruggenthaler, Jerome Faist, and Angel Rubio.  Engineering quantum materials with chiral optical cavities.  Nat. Mater., 20(4):438–442, Apr 2021.
  * [57] 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(7610):127–130, 2016.
  * [58] Alexey V Kavokin, Jeremy J Baumberg, Guillaume Malpuech, and Fabrice P Laussy.  Microcavities, volume 21.  Oxford university press, 2017.
  * [59] Edwin T Jaynes and Frederick W Cummings.  Comparison of quantum and semiclassical radiation theories with application to the beam maser.  Proc. IEEE, 51(1):89–109, 1963.
  * [60] Robert H Dicke.  Coherence in spontaneous radiation processes.  Phys. Rev., 93(1):99, 1954.
  * [61] Jonathan J Foley IV, Jonathan F McTague, and A Eugene DePrince III.  Ab initio methods for polariton chemistry.  arXiv:2307.04881, 2023.
  * [62] Michael Ruggenthaler, F Mackenroth, and Dieter Bauer.  Time-dependent kohn-sham approach to quantum electrodynamics.  Phys. Rev. A, 84(4):042107, 2011.
  * [63] 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.
  * [64] I. V. Tokatly.  Time-dependent density functional theory for many-electron systems interacting with cavity photons.  Phys. Rev. Lett., 110:233001, 2013.
  * [65] Michael Ruggenthaler.  Ground-state quantum-electrodynamical density-functional theory.  arXiv:1509.01417, 2015.
  * [66] Markus Penz, Erik I. Tellgren, Mihály A. Csirik, Michael Ruggenthaler, and Andre Laestadius.  The Structure of the Density-Potential Mapping. Part II: Including Magnetic Fields.  ACS Phys. Chem Au, 0(0):null, 0.
  * [67] Soeren Ersbak Bang Nielsen, Christian Schäfer, Michael Ruggenthaler, and Angel Rubio.  Dressed-orbital approach to cavity quantum electrodynamics and beyond.  arXiv preprint arXiv:1812.00388, 2018.
  * [68] Johannes Flick, Michael Ruggenthaler, Heiko Appel, and Angel Rubio.  Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry.  Proc. Natl. Acad. Sci. USA, 114:3026, 2017.
  * [69] Johannes Flick, Davis M Welakuh, Michael Ruggenthaler, Heiko Appel, and Angel Rubio.  Light–matter response in nonrelativistic quantum electrodynamics.  ACS Photonics, 6(11):2757–2778, 2019.
  * [70] Junjie Yang, Qi Ou, Zheng Pei, Hua Wang, Binbin Weng, Zhigang Shuai, Kieran Mullen, and Yihan Shao.  Quantum-electrodynamical time-dependent density functional theory within Gaussian atomic basis.  J. Chem. Phys., 155(6):064107, 2021.
  * [71] Marcus D Liebenthal, Nam Vu, and A Eugene DePrince III.  Mean-field cavity effects in quantum electrodynamics density functional and coupled-cluster theories.  arXiv:2303.10821, 2023.
  * [72] Davis M Welakuh, Johannes Flick, Michael Ruggenthaler, Heiko Appel, and Angel Rubio.  Frequency-Dependent Sternheimer Linear-Response Formalism for Strongly Coupled Light–Matter Systems.  J. Chem. Theory Comput., 18(7):4354–4365, 2022.