# Polaritonic response theory for exact and approximate wave functions

Matteo Castagnola  Rosario Roberto Riso  Department of Chemistry, Norwegian
University of Science and Technology, 7491 Trondheim, Norway  Alberto Barlini
Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy  Enrico
Ronca<EMAIL_ADDRESS>Dipartimento di Chimica, Biologia e
Biotecnologie, Università degli Studi di Perugia, Via Elce di Sotto, 8,06123,
Perugia, Italy  Henrik Koch<EMAIL_ADDRESS>Department of Chemistry,
Norwegian University of Science and Technology, 7491 Trondheim, Norway  Scuola
Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy

###### Abstract

Polaritonic chemistry is an interdisciplinary emerging field that presents
several challenges and opportunities in chemistry, physics and engineering. A
systematic review of polaritonic response theory is presented, following a
chemical perspective based on molecular response theory. We provide the reader
with a general strategy for developing response theory for ab initio cavity
quantum electrodynamics (QED) methods and critically emphasize details that
still need clarification and require cooperation between the physical and
chemistry communities. We show that several well-established results can be
applied to strong coupling light-matter systems, leading to novel perspectives
on the computation of matter and photonic properties. The application of the
Pauli-Fierz Hamiltonian to polaritons is discussed, focusing on the effects of
describing operators in different mathematical representations. We thoroughly
examine the most common approximations employed in ab initio QED, such as the
dipole approximation. We introduce the polaritonic response equations for
recently developed ab initio QED Hartree-Fock and QED coupled cluster methods.
The discussion focuses on the similarities and differences from standard
quantum chemistry methods, providing practical equations for computing the
polaritonic properties.

##  1 Introduction

Hybrid light-matter states, referred to as polaritons, were described in the
1950s from the interaction of optical modes and lattice excitons in
crystals.1, 2, 3 In 1963, Jaynes and Cummings described with a
phenomenological quantum picture the formation of polaritons from the
interaction of photons and atomic states.4 Their predictions were later
confirmed experimentally for Rydberg atoms in a quantum microcavity.5, 6 The
use of quantum optical devices to confine the fields, such as Fabry-Pérot
cavities or metal nanostructures,7, 8, 9, 10, 11, 12 increase the coupling
between matter and light, and polaritons have been observed since the 1970s
for inorganic semiconductors13, 14 and from the 1990s for organic molecules in
optical cavities.15, 16 Polaritonics has recently attracted the interest of
chemists since the pioneering work of Ebbesen and coworkers, who showed that
the formation of polaritons can influence photochemical17 and ground state
reactivity. 18, 19, 20 Several experiments have now shown that electromagnetic
confinement leads to essential modifications of processes such as chemical
reactions,17, 18, 19, 20, 21, 22 singlet fission,23, 24, 25 intersystem
crossing,26, 27, 28 crystallization and assembly,29, 30 as well as optical
properties like absorption, scattering and emission,31, 32, 33, 34, 35, 36,
37, 38, 39, 40, 41, 42, 43, 44, 45, 46 although the reproducibility of the
chemical reactivity modifications is not always straightforward.47 Crucial
experimental aspects are still to be clarified, and cooperation between
experimentalists and theoreticians can be of great assistance in unravelling
the conditions and the consequences of vibrational and electronic strong
coupling. Since both the photon and the molecular degrees of freedom are
important in such systems, their theoretical modelling is a challenging
interdisciplinary task.48, 49, 50, 51, 52, 53, 54, 55 Following well-
established electronic structure theory, several ab initio quantum
electrodynamics (QED) approaches have been developed recently. These include
QED density functional theory (QEDFT),56 reduced density matrix approaches for
matter-photon systems,57, 58 polaritonic coupled cluster,59 QED Hartree-Fock
(QED-HF),60, 61 QED coupled cluster (QED-CC),60, 61 QED full configuration
interaction (QED-FCI),60, 61 polarized Fock states approach (PFSs),62
polaritonic unitary coupled cluster (QED-UCC),63 strong coupling (SC)-QED-
HF,64 QED Møller-Plesset perturbation theory of second-order (QED-MP2),65 and
second-order QED algebraic diagrammatic construction scheme for the
polarization propagator (QED-ADC(2)).65 In addition, several approaches have
been investigated combining phenomenological models from the quantum optics
community with computational chemistry.66, 4, 67, 68, 69, 70, 71, 72, 73, 52,
74, 75, 76, 77, 78, 79, 80, 81, 82, 2, 83

The response theory formalism is a very successful and well-known framework
employed in quantum chemistry to calculate approximate ground and excited
state properties.84, 85, 86, 87, 88, 89, 90, 91, 92, 93 While several of these
results are still valid for QED systems, essential and subtle differences
arise due to the explicit modelling of the photon field, highlighting
additional challenges and opportunities. Although some response schemes for
QED methods have been proposed,94, 95, 96, 97, 98, 99, 100, 101 a systematic
discussion of polaritonic response theory and its connections to molecular
response theory is still not available.

This paper aims to fill this gap by proposing a systematic discussion of wave
function based methods for computing polaritonic properties. We will review
several results from molecular response theory that can be applied to
polaritonic wave functions and critically discuss the limitations and future
challenges such methods have to tackle. In section (2), we briefly introduce
the most common optical devices employed to achieve strong coupling. In
section (3), we present the Pauli-Fierz Hamiltonian and focus on how this
framework is used to model the strong coupling regime. We review the most
common approximations for ab initio QED in section (4). In section (5), we
discuss polaritonic response theory for exact states, focusing on similarities
and differences with molecular response theory, and provide several examples
of polaritonic response properties. The derived equations are equivalent to
QED-FCI response theory. In section (6), we discuss how to develop response
theory for approximate models and derive linear and quadratic response
equations for QED-HF and QED-CC. Finally, section (7) recaps the main results
described in this paper and contains some concluding remarks and future
perspectives.

##  2 Optical cavities

This section provides a brief introduction to quantum optical devices. The
discussion, although far from being exhaustive, aims at reviewing the most
common experimental setups employed in polaritonic chemistry and provides the
reader with phenomenological basics of light-matter strong coupling.

The use of quantum optical devices provides a non-invasive way to engineer a
material’s properties by exploiting the confinement of the electromagnetic
field. The geometry and the material of the device define the resonator
eigenmodes. In a real cavity, the confinement is always defective and depends
specifically on the frequency of the field and the cavity material. The so-
called $Q$-factor measures the quality of a cavity102

$Q(\omega_{0})=\omega_{0}\frac{\textit{energy stored}}{\textit{power loss}},$
(1)

where $\omega_{0}$ is the resonance frequency of the device. This equation
implies that the electric field will exhibit an exponential decay

$E(t)=E_{0}e^{-i\omega_{0}t}e^{-\frac{\omega_{0}}{2Q}t},$ (2)

such that the energy spectrum will show a Lorenzian lineshape

$|E(\omega)|^{2}\propto\frac{1}{(\omega-\omega_{0}^{2})+(\omega_{0}/2Q)^{2}}.$
(3)

From Eq. (3), we have a more practical definition of the $Q$-factor as the
ratio between the frequency of the hosted field mode $\omega_{0}$ and the full
width at half maximum (FWHM) $\Delta\omega$ of the cavity spectrum

$Q(\omega_{0})=\frac{\omega_{0}}{\Delta\omega}.$ (4)

The coupling $\lambda$ between matter and the confined cavity modes is
proportional to the inverse square root of the effective confinement volume
$V_{eff}$

$\lambda\propto\frac{1}{\sqrt{V_{eff}}}.$ (5)

The efficiency of a device in enhancing the interaction between light and
matter can then be measured by the $Q/V_{eff}$ ratio. Strong coupling
conditions are obtained once the coupling strength exceeds the material and
cavity losses (i.e. the coherent energy exchange rate between the
electromagnetic field and the molecule exceeds the dissipation processes).
Decoherence interactions with the environment are suppressed, and hybrid
states between the cavity and the molecular eigenstates naturally describe the
system. In this regime, the molecular properties are therefore affected even
in the absence of an external photon pumping of the device.

The most straightforward optical device to confine light is a Fabry-Pérot
cavity, composed of two planar parallel metal plates kept at a fixed
distance.103, 104, 105, 106, 107, 108 These resonators are often produced by
sputtering gold or silver on a substrate (often $SiO_{2}$, $BaF_{2}$ or
$ZnSe$), and direct contact with the sample is prevented by using a thin
polymer film. The sample is then coated by this multilayer structure or
injected into the hollow cavity. The cavity eigenmodes are linearly-polarized
standing waves whose wavelength $\lambda_{n}$ depends on the distance $L$
between the mirrors

$\lambda_{n}=\;\frac{2L}{n}\qquad n=1,2,\;\dots$ (6)

The spatial distribution of the optical mode can be observed by changing the
position of a thin sample slab inside the cavity. The strong coupling features
of this system are a function of the slab position, which gives us the field
strength distribution.46, 109 Changing the form of the plates, e.g. using
curved mirrors similar to lenses, can lead to modifications of the field shape
and possibly an increase in the coupling constant.110, 111, 112, 113, 114,
115, 116 The coupling strength $\lambda$ in Fabry-Pérot cavities is limited by
the diffraction limit, which restrains the field confinement. Devices which
bound the fields on a scale smaller than the diffraction limit are called
subwavelength cavities and can facilitate ultrastrong coupling. To this end,
using nanostructured metamaterials in a suitable geometric arrangement can
allow efficient confinement and overcome the diffraction limit, leading to
efficient subwavelength cavities. In Fig. (1a), we show an example of a
subwavelength "dogbone" metamaterial resonator, which can confine light on a
scale of about one-tenth of the mode wavelength.117, 118, 119, 120, 121, 122
Metallic nanostructures can be used to confine electromagnetic fields through
surface plasmon resonances, allowing for remarkably low-volume confinement. It
is also possible to employ non-metallic plasmonic structures, such as
graphene.123, 124, 125, 126, 127, 128 Coupling two or more plasmonic
nanoparticles can lead to extreme field localization in hot spots between
them, as schematically shown in Fig. (1b) for two nano prisms on a silicon
substrate.129, 130, 131, 132, 133 However, these systems can be very lossy,
i.e. they have a low $Q$-factor, which effectively limits their use for strong
coupling. Alternatively, dielectric materials can also be used to confine the
electromagnetic field. In Fig. (1c), we show a cavity obtained by two
distributed Bragg reflectors (DBRs) composed of several alternate dielectric
slabs. By carefully choosing the number of slabs and their material, DBRs can
achieve low-volume confinement and high $Q$-factors.134, 135, 136, 137, 138
Transparent high-reflection index dielectrics in curved or polygonal shapes
can confine electromagnetic waves by total internal reflection. The optical
modes of such resonators, called whispering-gallery modes (WGMs), can achieve
extreme field localization and quality factors.139, 140, 141, 142, 143, 144,
145, 146
The devices described so far cannot, in general, host any chiral eigenmodes.
The helicity of a circularly polarized wave changes sign when reflected by a
mirror surface, creating a mode with zero global helicity. In recent years,
different helicity-preserving optical devices (chiral cavities) have been
proposed based on the introduction of a chiral element in the resonator.147,
148, 149, 150, 115, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161 In
Fig. (1d), we show a chiral Fabry-Pérot cavity formed by the insertion of a 2D
chiral polymer slab.147 The 2D chirality is achieved by applying torsional
sheer stress, and chiral eigenmodes are sustained in the device by a circular
conversion dichroism mechanism.147 Other chiral optical resonators have been
proposed using chiral metasurfaces115, 151, 149, 150 and chiral
plasmonics.161, 160, 159, 158, 157

Figure 1: (a) Metallic nanostructures can be used to confine electromagnetic
fields through surface plasmon resonances. A "dogbone" metamaterial can be
employed as a subwavelength resonator, confining the electromagnetic field on
a scale of a tenth of the mode wavelength.119 Subwavelength optical devices
can overcome the diffraction limit by employing suitable geometries and
metamaterials, facilitating the fulfilment of strong coupling conditions.117,
118, 119, 120, 121, 122 (b) Two metal nanoparticles are coupled on a crystal
substrate, increasing the electric field’s intensity in a hotspot between
them.129, 130, 131, 132, 133 Metal nanoparticles can be used to achieve very
low-volume confinement due to surface plasmon resonances, although dissipation
processes are often significant. (c) Distributed Bragg reflectors (DBRs) are
composed of alternate slabs of dielectric materials. Carefully choosing the
number of slabs and their material, it is possible to achieve efficient
confinement of the electromagnetic fields between two DBRs.134, 135, 136, 137,
138 (d) A 2D chiral polymer film, obtained by torsional sheer stress, inserted
inside a Fabry-Pérot resonator endows the system with chiral eigenmodes.147

##  3 Pauli-Fierz Hamiltonian

The starting point of molecular modelling is the nonrelativistic molecular
Hamiltonian

$\displaystyle
H=\sum_{M}\frac{1}{2m_{M}}\mathbf{p}_{M}^{2}+\sum_{i}\frac{1}{2}\mathbf{p}_{i}^{2}$
$\displaystyle+\frac{1}{2}\sum_{i\neq
j}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}-\sum_{i,M}\frac{Z_{M}}{|\mathbf{r}_{i}-\mathbf{R}_{M}|}+\frac{1}{2}\sum_{M\neq
N}\frac{Z_{M}Z_{N}}{|\mathbf{R}_{N}-\mathbf{R}_{M}|}.$ (7)

In Eq. (3), the capital indices $M$ and $N$ refer to nuclei of mass and charge
$m_{M}$ and $Z_{N}$, $i$ and $j$ to electrons with mass $m=1$ and charge
$q=-1$. The electronic and nuclear positions are labelled $\mathbf{r}$ and
$\mathbf{R}$ respectively, and $\mathbf{p}$ is the linear momentum of the
particles. The first line of Eq. (3) contains the nonrelativistic kinetic
energy operators for nuclei and electrons, while the second line describes the
Coulomb interaction.

Although it is often assumed that the Hamiltonian in Eq. (3) is sufficient to
model a molecular system, this is not the case. This is apparent if we realise
that the excited states of the system have an infinite lifetime, i.e. they do
not radiate. In Eq. (3), we disregard the electromagnetic degrees of freedom,
retaining only the electrostatic interaction between particles while
neglecting retardation effects. The motion of electric charges generates
radiation fields which in turn influence the particles’ motion: this back-
reaction mechanism is the origin of spontaneous decay. The inclusion of the
electromagnetic degrees of freedom is also fundamental when optical devices
are employed to modify the photon environment. It is then necessary to go
beyond the Coulomb interaction terms and resort to Maxwell’s electrodynamics
to describe self-consistently the interaction between light and matter.162,
163, 102, 164 The particle motion is determined by the Lorenz force
$\mathbf{F}_{L}$

$\mathbf{F}_{L}=m{\mathbf{a}}=q\bigg{(}\mathbf{E}(\mathbf{r},t)+\frac{{\bm{v}}}{c}\times\mathbf{B}(\mathbf{r},t)\bigg{)},$
(8)

where $c$ is the speed of light, $\bm{v}$ and $\mathbf{a}$ are the velocity
and the acceleration of the particle, $q$ is the charge and $m$ is the mass.
The electric $\mathbf{E}$ and magnetic $\mathbf{B}$ fields are determined by
Maxwell’s equations, here expressed in cgs units

$\begin{cases}\nabla\cdot\mathbf{E}=4\pi\rho(\mathbf{r},t)\\\
\nabla\cdot\mathbf{B}=0\\\
\nabla\times\mathbf{E}=-\frac{1}{c}\frac{\partial\mathbf{B}}{\partial t}\\\
\nabla\times\mathbf{B}=\frac{4\pi}{c}\mathbf{j}(\mathbf{r},t)+\frac{1}{c}\frac{\partial\mathbf{E}}{\partial
t}\end{cases},$ (9)

where $\mathbf{r}$ is a vector defining a point in space, $\rho(\mathbf{r},t)$
is the charge density and $\mathbf{j}(\mathbf{r},t)$ is the current density.
The fields are naturally described using potentials $(\phi,\mathbf{A})$, such
that

$\displaystyle\mathbf{B}$ $\displaystyle=\nabla\times\,\mathbf{A}$ (10)
$\displaystyle\mathbf{E}$
$\displaystyle=-\nabla\,\mathbf{\phi}-\frac{1}{c}\frac{\partial\mathbf{A}}{\partial
t}.$ (11)

We notice that, unlike the fields $\mathbf{B}$ and $\mathbf{E}$, the vector
potential $\mathbf{A}$ and the scalar potential $\phi$ are not measurable
quantities and must, therefore, only be seen as mathematical tools. In Eqs.
(10) and (11), the potentials are identified only through their derivatives.
This implies that the same fields, i.e. the same system, can be modelled using
different potentials. The transformation

$\displaystyle\mathbf{A}$ $\displaystyle\to\mathbf{A}+\nabla\,\chi$ (12)
$\displaystyle\mathbf{\phi}$
$\displaystyle\to\mathbf{\phi}-\frac{1}{c}\frac{\partial\chi}{\partial t},$
(13)

known as gauge transformation, leads to unchanged electromagnetic fields for
any well-behaved scalar function $\chi(\mathbf{r},t)$. This gauge freedom can
be exploited to simplify the equations. We will make use of the so-called
Coulomb gauge

$\nabla\cdot\mathbf{A}=0$ (14)

which allows us to account for the Coulomb interactions among the particles
explicitly. The scalar potential is, in this gauge, determined by the Poisson
equation as in the electrostatic case

$\nabla^{2}\phi=-4\pi\rho,$ (15)

while the vector potential is obtained from the following equation:

$\nabla^{2}\mathbf{A}-\frac{1}{c^{2}}\frac{\partial^{2}\mathbf{A}}{\partial
t^{2}}=-\frac{4\pi}{c}\mathbf{j}^{\perp},$ (16)

where $\mathbf{j}^{\perp}$ is the solenoidal (divergence-free) part of the
density current.

The Hamiltonian for the field-particle system is described in terms of the
position of the particles, the vector potential and their conjugate momenta.
The self-consistent interaction between molecules and light is described by
the Pauli-Fierz Hamiltonian, for which Hamilton’s equations can be shown to be
equivalent to Eqs. (8) and (9)162, 163

$\displaystyle H_{PF}$
$\displaystyle=\sum_{M}\frac{1}{2m_{M}}\bigg{(}\mathbf{p}_{M}-\frac{Z_{M}}{c}\mathbf{A}(\mathbf{R}_{M})\bigg{)}^{2}+\sum_{i}\frac{1}{2}\bigg{(}\mathbf{p}_{i}+\frac{1}{c}\mathbf{A}(\mathbf{r}_{i})\bigg{)}^{2}$
$\displaystyle+\frac{1}{2}\sum_{i\neq
j}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}-\sum_{i,M}\frac{Z_{M}}{|\mathbf{r}_{i}-\mathbf{R}_{M}|}+\frac{1}{2}\sum_{M\neq
N}\frac{Z_{M}Z_{N}}{|\mathbf{R}_{N}-\mathbf{R}_{M}|}$
$\displaystyle+\frac{1}{2}\int\left[\mathbf{E}_{\perp}\left(\mathbf{r}\right)^{2}+\mathbf{B}\left(\mathbf{r}\right)^{2}\right]d\mathbf{r}.$
(17)

In the first line of Eq. (3), we have the vector potential $\mathbf{A}$
computed at the positions of nuclei and electrons, and the second line
describes the Coulomb interactions between the particles. The explicit
separation of Coulomb interactions provided by the Coulomb gauge allows for
the separation of Hamiltonian in Eq. (3), providing a natural chemical
description of the system. The last line accounts for the electromagnetic
field energy. Expanding the first terms, we obtain an interaction term between
the particle Hamiltonian in Eq. (3) and the radiation fields mediated by the
vector potential $\mathbf{A}$

$H_{int}=-\sum_{x}\frac{q_{x}}{2m_{x}c}\big{(}\mathbf{p}_{x}\cdot\mathbf{A}(\mathbf{r}_{x})+\mathbf{A}(\mathbf{r}_{x})\cdot\mathbf{p}_{x}\big{)}+\sum_{x}\frac{q^{2}_{x}}{2m_{x}c^{2}}\mathbf{A}^{2}(\mathbf{r}_{x}),$
(18)

where $x$ runs on both electrons and nuclei. Following the QED prescription,
the vector potential in Eq. (3) is promoted to an operator

$\mathbf{A}\longrightarrow\hat{\mathbf{A}}=\sum_{\mathbf{k},\tau}\frac{c\lambda}{\sqrt{2\omega_{k}}}\left(\bm{\epsilon}_{\mathbf{k}\tau}e^{i\mathbf{k}\cdot\mathbf{r}}b_{\mathbf{k}\tau}+\bm{\epsilon}^{*}_{\mathbf{k}\tau}e^{-i\mathbf{k}\cdot\mathbf{r}}b^{\dagger}_{\mathbf{k}\tau}\right),$
(19)

where $b^{\dagger}_{\mathbf{k}\tau}$ and $b_{\mathbf{k}\tau}$ respectively
create and annihilate a photon with wave vector $\mathbf{k}$ and frequency
$\omega_{k}=\frac{|\mathbf{k}|}{c}$. The field polarization
$\bm{\epsilon}_{\mathbf{k}\tau}$ with $\tau=1,2$ spans the 2D plane
perpendicular to $\mathbf{k}$. If we take the polarization
$\bm{\epsilon}_{\mathbf{k}\tau}$ to be real, the fields are described by a
superposition of linearly polarized plane waves. In an equivalent way, we can
take the complex orthogonal unit vectors

$\bm{\epsilon}_{\mathbf{k}\pm}=\frac{1}{\sqrt{2}}(\bm{\epsilon}_{\mathbf{k}1}\pm
i\bm{\epsilon}_{\mathbf{k}2})$ (20)

which describe left and right circularly polarized waves.

#### Remarks

There are now some subtle points to mention. Choosing the Coulomb gauge, the
vector potential only has (two) transverse components
$\mathbf{A}=\mathbf{A}_{\perp}$. At the same time, the longitudinal component
of the electric field (Coulomb field) is described by the scalar potential. In
this gauge, the equations are not manifestly covariant. However, this does not
imply a loss of relativistic invariance since the electromagnetic field
prediction agrees with special relativity. In fact, the retardation of the
interaction between particles is obtained by a cancellation of the
instantaneous interactions from the Coulomb term and the transverse field. For
a different gauge, the Coulomb interactions would also be mediated by the
longitudinal component of the vector potential.162, 163

The canonical momentum $\bm{p}_{i}$, contrary to the standard quantum
mechanical formulation, is not the linear momentum of the particle but here
has a field-dependent component163, 162

$\bm{p}_{i}=m_{i}\bm{v}_{i}+\frac{q_{i}}{c}\mathbf{A}(\bm{r}_{i}).$ (21)

This means that, although in Eq. (3) we can identify the original particle
Hamiltonian in Eq. (3), the physical meaning of these two operators is
different.

The Pauli-Fierz Hamiltonian is self-adjoint and bounded from below165, 166
such that a stable ground state exists. However, the original excited states
of the particle Hamiltonian Eq. (3) become metastable (i.e. with a finite
lifetime) because of their coupling to the continuous photonic spectrum.167,
168, 169, 170

In Eq. (3), we have disregarded the nuclear and electronic spin, which could
be introduced in Eq. (3) with the addition of the following contribution

$H_{spin}=\frac{1}{2}\sum_{i}\bm{\sigma}_{i}\cdot\mathbf{B}(\bm{r}_{i})-\sum_{N}\gamma_{N}\bm{I}_{N}\cdot\mathbf{B}(\bm{R}_{N}),$
(22)

where the first term of Eq. (22) accounts for the interaction of the
electronic spin, described by the Pauli matrices $\bm{\sigma}_{i}$, with the
(internal) magnetic field $\mathbf{B}$, while the second term describes the
energy associated with the nuclear spin $\bm{I}_{N}$ with magnetogyric ratio
$\gamma_{N}$.

Finally, in Eq. (3) particles can in principle interact with infinitely high-
frequency modes. However, the description of the kinetic energy is
nonrelativistic. When the field frequency is comparable with the rest energy
of the particles $\omega\sim mc^{2}$, relativistic effects (e.g. creation of
electron-positron pairs) not included in Eq. 3 become relevant.162 It is then
necessary to introduce a cutoff as the field’s momentum $\mathbf{k}$ increases
(ultraviolet cutoff). Moreover, the masses in Eq. (3) are not the physical
masses usually employed in quantum mechanics, which include the contribution
of the electromagnetic energy created by the particle, but they are their bare
masses.163

##  4 Approximations for ab initio polaritonic chemistry

While the Hamiltonian in Eq. (3) allows for a consistent treatment of light
and matter, approximations are needed to perform computational studies. The
problem of the continuum of photonic modes could be overcome by employing a
fine discretization of the spectrum. Moreover, different approaches can be
developed depending on the treatment of the photonic degrees of freedom. In
the Cavity Born-Oppenheimer approach (CBO),171, 172, 173 the photon
coordinates $q_{\alpha}$ are embedded in the nuclear wave function
$\chi(\bm{R},\bm{q},t)$ and separated from the electronic degrees of freedom,
described by and electronic wave function $\phi(\bm{r};\bm{R},\bm{q})$. The
complete wave function can then be expanded as

$\Psi(\bm{r},\bm{R},\bm{q},t)=\sum_{k}\phi_{k}(\bm{r};\bm{R},\bm{q})\chi_{k}(\bm{R},\bm{q},t).$
(23)

In a polaritonic approach, the electronic and photon coordinates are treated
on the same footing and described by a polaritonic wave function which depends
parametrically on the nuclear coordinates only81, 174, 82, 175, 60, 56, 176

$\Psi(\bm{r},\bm{R},\bm{q},t)=\sum_{k}\phi_{k}(\bm{r},\bm{q};\bm{R})\chi_{k}(\bm{R},t).$
(24)

While these Born-Huang expansions are in principle equivalent, approximations
will lead to different results. In the simplest case, we take into account
only a single term of these expansions as in standard electronic Born-
Oppenheimer approximation. In the discussion of exact response theory in
Section (5) we do not refer to any explicit form of the Hamiltonian. In
Section (6), we focus on the polaritonic approach and describe the QED-CC and
QED-HF response theory. The interaction between the nuclear motion and the
electromagnetic degrees of freedom can lead to novel non-adiabatic coupling
terms that could, in some situations, jeopardize the validity of the BO
approximation.51, 177, 28, 178, 81, 179 Nevertheless, we will not explore such
effects, and we will work in a fixed-nuclei framework. Although these
approximations already reduce the complexity of the problem, further
simplifications are needed to deal with the electromagnetic environment and
treat the molecule and its interaction with the fields with chemical accuracy.

###  4.1 Cavity QED

The Hamiltonian in Eq. (3) allows for simultaneous treatment of light and
matter, where we could, in principle, include any optical device that
contributes to a modification of the electromagnetic environment.
Nevertheless, the molecular complexity and the presence of the optical
apparatus make the problem impractical from a computational point of view.

The explicit modelling of the photonic device is avoided by defining an
effective change in the structure of the electromagnetic modes. To this end,
we impose boundary conditions on the fields (e.g. vanish on the surface of the
optical device) to model the electromagnetic confinement induced by the
cavity. The vector potential is expanded in terms of a complete orthonormal
set of suitable electromagnetic normal mode functions
$\\{\bm{S}_{\alpha}(\bm{r})\\}$.102, 180 For a rectangular parallelepiped with
perfectly conducting walls of sides $L_{x}$, $L_{y}$ and $L_{z}$, the
transverse mode functions are standing waves180

$\displaystyle S_{x}(\bm{r})$
$\displaystyle=\sqrt{\frac{8}{V_{c}}}\cos\big{(}k_{x}x\big{)}\sin\big{(}k_{y}y\big{)}\sin\big{(}k_{z}z\big{)}$
$\displaystyle S_{y}(\bm{r})$
$\displaystyle=\sqrt{\frac{8}{V_{c}}}\sin\big{(}k_{x}x\big{)}\cos\big{(}k_{y}y\big{)}\sin\big{(}k_{z}z\big{)}$
$\displaystyle S_{z}(\bm{r})$
$\displaystyle=\sqrt{\frac{8}{V_{c}}}\sin\big{(}k_{x}x\big{)}\sin\big{(}k_{y}y\big{)}\cos\big{(}k_{z}z\big{)},$

where $V_{c}=L_{x}\times L_{y}\times L_{z}$ is the total volume of the
parallelepiped and $\mathbf{k}$ is the wave vector which is now quantized:

$k_{i}=n\frac{\pi}{L_{i}}\qquad n=1,2,\,\dots$ (25)

In the absence of external charges, if we assume that the walls of the
parallelepiped are grounded, the scalar potential vanishes identically
${\phi}=0$. The quantization is then performed by promoting the expansion
coefficients to quantum operators, and the radiation Hamiltonian for the free
fields inside the perfect optical device is discretized in terms of the cavity
eigenmodes

$H_{rad}=\sum_{\mathbf{\alpha}}\omega_{\mathbf{\alpha}}b^{\dagger}_{\mathbf{\alpha}}b_{\mathbf{\alpha}},$
(26)

where $b^{\dagger}_{\mathbf{\alpha}}$ is the creation operator for the
$\alpha$-mode of frequency $\omega_{\mathbf{\alpha}}$.

In the presence of molecules, which can be considered as a source of
electromagnetic fields, we must also ensure that the scalar potential in Eq.
(15) is consistent with the boundary conditions. We can enforce them by adding
an auxiliary scalar field ${F}$ to the particles’ Coulomb potential. This
potential can be thought of as a potential generated by image charges102, 180
placed outside the volume of quantization (i.e. outside the volume of interest
of the system). This potential can be expressed in terms of a response kernel
$\Theta(\bm{s};\bm{r})$ and the real charge distribution $\rho$ inside the
cavity102, 180

${F}(\bm{r})=\int_{\bar{V}}d^{3}s\;\frac{\int_{V_{c}}d^{3}r^{\prime}\;\rho(\bm{r}^{\prime})\Theta(\bm{s};\bm{r}^{\prime})}{|\bm{r}-\bm{s}|}$
(27)

where $\bm{s}$ is a point outside the cavity volume. The Hamiltonian
describing the molecular system then reads180

$\displaystyle H_{c}=$
$\displaystyle\sum_{i}\frac{1}{2}\mathbf{p}_{i}^{2}+\frac{1}{2}\sum_{i\neq
j}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}-\sum_{i,M}\frac{Z_{M}}{|\mathbf{r}_{i}-\mathbf{R}_{M}|}+\frac{1}{2}\sum_{M\neq
N}\frac{Z_{M}Z_{N}}{|\mathbf{R}_{N}-\mathbf{R}_{M}|}$
$\displaystyle+\sum_{i}\frac{1}{2c}(\mathbf{p}_{i}\cdot\mathbf{A}(\bm{r}_{i})+\mathbf{A}(\bm{r}_{i})\cdot\mathbf{p}_{i})+\sum_{i}\frac{1}{2c^{2}}\mathbf{A}^{2}(\bm{r}_{i})$
$\displaystyle+\sum_{\alpha}\omega_{\alpha}{b}^{\dagger}_{\alpha}{b}_{\alpha}+\frac{1}{2}\sum_{M}Z_{M}{F}(\bm{R}_{M})-\frac{1}{2}\sum_{i}{F}(\bm{r}_{i}).$
(28)

As seen from (27), the last two terms introduced by ${F}$ correspond to a
modified Coulomb interaction kernel due to the boundary conditions, i.e. they
lead to a modified longitudinal interaction among the particles. The factor
$1/2$ appears because it describes an interaction with an image charge
distribution.102 While some authors suggest that this contribution is the
major difference from free space,181, 182, 183 the discussion of these terms
for ab initio QED is often neglected. Nevertheless, we will see that they can
be handled in a more practical way by a suitable unitary transformation.180,
72

#### Remarks

The Hamiltonian in Eq. (4.1) now acts on a Hilbert space $\mathcal{V}$ which
is the direct product of the full-CI Fock space for the electrons
$\mathcal{V}_{FCI}$ and the photon space spanned by the cavity eigenmodes
$\mathcal{V}_{C}$

$\mathcal{V}=\mathcal{V}_{FCI}\otimes\mathcal{V}_{C}.$ (29)

However, each cavity setup corresponds to specific boundary requirements,
leading to different quantum fields: different boundary conditions could
substantially modify both the shape of the eigenmodes
$\bm{S}_{\alpha}(\bm{r})$ and the Coulomb kernel $\Theta(\bm{s};\bm{r})$.

The photonic space $\mathcal{V}_{C}$ is usually truncated, considering only a
few eigenmodes relevant to the system. Therefore, the excited states of the
Hamiltonian in Eq. (4.1) are again true eigenstates with an infinite lifetime.
Notice that, in the limit of an infinitely large cavity, if we retain a
considerable number of modes we obtain again a fine energy discretization of
the photon bath that recovers the decay channels.48

The masses of the particles in (4.1) are dependent on the electromagnetic mode
structure and, in general, are different from their physical and free-space
bare masses.163, 184 However, in QED ab initio calculations, the physical
masses are usually employed, implicitly considering the interaction energy
with the continuum photon bath48. The effects of the modification of the
masses on chemical and physical properties from an ab initio perspective, to
the best of our knowledge, is still to explore.

Finally, because of the states of the optical device $\mathcal{V}_{C}$, the
system’s symmetries now differ from the bare molecule ones due to the cavity’s
modified eigenmodes and potential in Eq. (27). This leads to modifications of
selection rules and other symmetry-related properties.

###  4.2 Dipole approximation

The dipole approximation is commonly employed in ab initio polaritonic
chemistry. It assumes that the (relevant) electromagnetic modes have a
wavelength much larger than the characteristic lengths of the molecules. The
complete shape of the modes $\bm{S}_{\bm{\alpha}}(\bm{r})$ is therefore
irrelevant, and we only need to evaluate the fields at a point internal to the
molecular structure. In Eq. (3), we then set
$\mathbf{A}(\bm{r})\to\mathbf{A}(0)$, and obtain the Pauli-Fierz Hamiltonian
$H_{PF}^{v}$ in dipole approximation and velocity representation

$\displaystyle H_{PF}^{v}$
$\displaystyle=\sum_{i}\frac{1}{2}\mathbf{p}_{i}^{2}+\frac{1}{2}\sum_{i\neq
j}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}-\sum_{i,M}\frac{Z_{M}}{|\mathbf{r}_{i}-\mathbf{R}_{M}|}+\frac{1}{2}\sum_{M\neq
N}\frac{Z_{M}Z_{N}}{|\mathbf{R}_{N}-\mathbf{R}_{M}|}$
$\displaystyle+\sum_{i}\frac{1}{2c}(\mathbf{p}_{i}\cdot\mathbf{A}(0)+\mathbf{A}(0)\cdot\mathbf{p}_{i})+\sum_{i}\frac{1}{2c^{2}}\mathbf{A}^{2}(0)$
$\displaystyle+\sum_{\alpha}\omega_{\alpha}{b}^{\dagger}_{\alpha}{b}_{\alpha}+\frac{1}{2}\sum_{M}Z_{M}\tilde{F}(\bm{R}_{M})-\frac{1}{2}\sum_{i}\tilde{F}(\bm{r}_{i}),$
(30)

where we applied the BO approximation. The dipole approximation allows for the
modelling of the photons only by adequately tuning the frequency
$\omega_{\alpha}$ and the coupling strength ${\lambda}_{\alpha}$ for the
cavity modes.48 The problem is eventually reduced to a single effective cavity
mode. Assuming to work with real field polarization vectors, the final
Hamiltonian for our system is

$\displaystyle H$
$\displaystyle=\sum_{pq}h_{pq}E_{pq}+\frac{1}{2}\sum_{pqrs}g_{pqrs}e_{pqrs}+h_{nuc}$
$\displaystyle+\sqrt{\frac{\omega_{\alpha}}{2}}(\bm{\lambda}_{\alpha}\cdot\bm{d})(b^{\dagger}_{\alpha}+b_{\alpha})+\frac{1}{2}(\bm{\lambda}_{\alpha}\cdot\bm{d})^{2}$
$\displaystyle+\omega_{\alpha}{b}^{\dagger}_{\alpha}{b}_{\alpha}+\frac{1}{2}\sum_{M}Z_{M}\tilde{F}(\bm{R}_{M})-\frac{1}{2}\sum_{i}\tilde{F}(\bm{r}_{i}),$
(31)

where $\bm{\lambda}_{\alpha}=\lambda_{\alpha}\bm{\epsilon}_{\alpha}$ is the
light-matter coupling vector of the $\alpha$ electromagnetic mode, and we used
standard second quantization notation for the electronic Hamiltonian.185 This
formulation is suitable for nonchiral cavities, while for helicity-preserving
devices, we necessarily need to work with complex polarization vectors (see
Eqs (19) and (20)). The Hamiltonian in Eq. (4.2) can be recast in a more
convenient form, called length representation $H_{PF}^{l}$, by employing a
unitary transformation $U$ such that

$U\mathbf{p}U^{\dagger}=\mathbf{p}-\frac{1}{c}\mathbf{A}(0),$ (32)

which will allow us to write the interaction between the molecule and the
fields in a manifestly dipolar fashion. The transformation is

$U=\text{exp}\bigg{[}i\sum_{i}\mathbf{r}_{i}\cdot\frac{1}{c}\mathbf{A}(0)\bigg{]}=\text{exp}\bigg{[}-i\mathbf{d}\cdot\frac{1}{c}\mathbf{A}(0)\bigg{]}$
(33)

where $\mathbf{d}$ is the electronic dipole operator. By using Eq. (33) and a
transformation which changes the phases of the creation/annihilation
operators186, 162, 187, 188

$V=\text{exp}\bigg{[}-i\frac{\pi}{2}\sum_{\alpha}b^{\dagger}_{\alpha}b_{\alpha}\bigg{]},$
(34)

we finally obtain

$\displaystyle H_{PF}^{l}=VUH_{PF}^{v}U^{\dagger}V^{\dagger}$
$\displaystyle=\sum_{i}\frac{1}{2}\mathbf{p}_{i}^{2}+\frac{1}{2}\sum_{i\neq
j}\frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}-\sum_{i,M}\frac{Z_{M}}{|\mathbf{r}_{i}-\mathbf{R}_{M}|}+\frac{1}{2}\sum_{M\neq
N}\frac{Z_{M}Z_{N}}{|\mathbf{R}_{N}-\mathbf{R}_{M}|}$
$\displaystyle+\sum_{\alpha}\sqrt{\frac{\omega_{\alpha}}{2}}(\bm{\lambda}_{\alpha}\cdot\bm{d})(b^{\dagger}_{\alpha}+b_{\alpha})+\frac{1}{2}\sum_{\alpha}(\bm{\lambda}_{\alpha}\cdot\bm{d})^{2}$
$\displaystyle+\sum_{\alpha}\omega_{\alpha}{b}^{\dagger}_{\alpha}{b}_{\alpha}$
(35)

The bilinear light-matter interaction term of this Hamiltonian can be
interpreted as the interaction of the molecular dipole with the displacement
field.162 The quadratic term is called dipole self-energy and ensures the
Hamiltonian leads to well-defined states and properties.186 Moreover,
transformation Eq. (33) cancels the image charge contribution from $\tilde{F}$
so that no reference to the image charges appears in Eq. (4.2).180, 72, 163
Notice that these Hamiltonians are also often referred to as length or
velocity gauge, although in a QED framework this can be misleading, and the
term form or representation is more appropriate. The hierarchy of
approximations that are employed to study quantum light-matter systems is
pictorially summarized in Fig. 2.

#### Remarks

There are several advantages to employing the Pauli-Fierz Hamiltonian in the
length representation. Contrary to the velocity representation Eq. (4.2), Eq.
(4.2) has a light-matter interaction term linear in the field with no
reference to the image charge distribution.

The transformation in Eq. (33) implies that the momentum operator $\mathbf{p}$
now represents the kinetic momentum of particles

$\mathbf{p}=m\bm{v},$ (36)

so that the first line of Eq. (4.2) truly represents the Coulomb and kinetic
energy of the matter subsystem only.

The operator $b^{\dagger}$ is no longer a purely photonic operator since $U$
in Eq. (33) mixes electromagnetic and matter degrees of freedom. The original
photon creation operator in the length representation after the tranformatoin
is186, 188, 162

$VUb^{\dagger}_{\alpha}U^{\dagger}V^{\dagger}=-ib_{\alpha}^{\dagger}-i\frac{1}{\sqrt{2\omega_{\alpha}}}(\bm{\lambda}_{\alpha}\cdot\bm{d}).$
(37)

These operators are now connected to the auxiliary fields of the macroscopic
Maxwell’s equations, where we can recognize the transverse polarization162,
186

$\mathbf{P}_{\perp}=\frac{1}{4\pi}\sum_{\alpha}(\bm{\lambda}_{\alpha}\cdot\bm{d})\bm{\lambda}_{\alpha}.$
(38)

The separation into "matter" and "photon" degrees of freedom is, therefore,
blurred in this representation.162

Figure 2: Graphic summary of the hierarchy of approximations for computational
polaritonic chemistry. The starting point is the open system described by the
nonrelativistic Pauli-Fierz Hamiltonian, which includes the molecular system
and the optical device immersed in the photon continuum. As a first
approximation, the device is assumed to be perfect, and cavity QED formalism
with a limited number of effective photon states is employed. Then, in the
dipole approximation, the field is assumed to be uniform over the molecular
scale. The exact shape of the photon states is irrelevant, and the system is
described through an effective coupling. At the ab initio level, the
polaritonic problem is approximated with a suitable parametrization.60, 56,
64, 58, 65, 59, 62, 61, 63 Alternatively, there are phenomenological models
where the molecular complexity is usually reduced to a few selected reference
states, as, for instance, in the Jaynes-Cummings model.4

###  4.3 Beyond the dipole approximation

Ab initio QED calculations often rely upon the dipole approximation, which
provides a good compromise between accuracy and affordability. Nevertheless,
achieving a better description of the light-matter interaction is an essential
challenge for polaritonic systems, and a few attempts in this direction have
already been reported.189, 190, 191 To go beyond the dipole approximation, we
could replace the vector potential $\mathbf{A}(\bm{r})$ with its first-order
expansion around the molecular origin. For a sinusoidal field, this reduces to
an expansion in terms of the wavevector $\bm{k}$

$e^{i\bm{k}\cdot\bm{r}}\approx 1+i\bm{k}\cdot\bm{r}+\dots$ (39)

This expansion is often employed in the semiclassical description of the
fields to model, for instance, X-ray absorption leading to the well-known
electric quadrupole and magnetic dipole interaction terms.192, 193, 194 The
corresponding QED Hamiltonian can then be obtained in the velocity
representation by the following substitution in Hamiltonian Eq. (4.1)

$\mathbf{A}(\bm{r})\to\mathbf{A}(0)+\bm{r}\cdot\frac{\partial\mathbf{A}}{\partial\bm{r}}\bigg{|}_{\bm{r}=0}+\dots$
(40)

The multipolar Hamiltonian can be obtained similarly to the length
representation of the dipole Hamiltonian, applying the Power-Zienau-Woolley
(PZW) transformation with the multipolar expansion of the vector potential,163
i.e. the transformation of Eq. (33) with the vector potential of Eq. (40).
This transformation also leads to electric (dipole and quadrupole) and
magnetic self-energy terms, which guarantee the boundedness of the
Hamiltonian.163, 186, 190 Overcoming the dipole approximation is particularly
relevant when the molecules interact with chiral fields, which leads to novel
phenomena such as cavity-induced circular dichroism or cavity enantiomeric
discrimination.64, 190 The description of magnetic interactions is then
essential. However, the multipolar expansion of the interaction beyond the
electric dipole presents some difficulties, already in the semiclassical
approximation,192, 193, 194, 195 as the interaction terms are now origin
dependent. The dipole operator also depends on the choice of the origin for
charged systems, but this can be shown to be equivalent to a gauge
transformation on Hamiltonian Eq. (4.2). Manifest origin invariance can then
be recovered by a suitable coherent state transformation, even for a finite
electronic basis set (see Eqs. (201) and (6.1)).60, 61, 64 At the same time,
to the best of our knowledge, a similar solution has yet to be developed for
higher-order Hamiltonians, and therefore using such an expansion could lead to
unphysical results. This issue could be solved by retaining the complete shape
of the field $\bm{S}_{\alpha}(\bm{r})$ instead of relying upon a multipolar
expansion. This is considered computationally and theoretically challenging
also because of the dependence of the field shape on the boundary conditions.
Nevertheless, some progress has been made for a sinusoidal field in the
semiclassical and QED frameworks. 196, 194, 197, 198, 189 Alternatively,
effective reformulations of the QED problem in terms of reduced quantities,
such as in QEDFT, can convey a different perspective for higher-order
approximations.56, 187, 199, 200, 201

##  5 Exact polaritonic response theory

In this section, we highlight the similarities and differences with standard
molecular response theory following the derivation of Jørgensen and Olsen,84
with particular emphasis on the quantities that describe the cross-talk
between light and matter.

From the exact solution of the time-dependent Schrödinger equation

$i\frac{d}{dt}\ket{\Psi}={H}\ket{\Psi}$ (41)

we derive, for an operator $\Omega$, the Ehrenfest theorem202

$\frac{d}{dt}\braket{\Psi}{\Omega}{\Psi}=-i\braket{\Psi}{[\Omega,H]}{\Psi}+\braket{\Psi}{\frac{\partial\Omega}{\partial
t}}{\Psi},$ (42)

where ${H}$ is here a polaritonic Hamiltonian and $\ket{\Psi}$ are polaritonic
states of entangled light-matter character. The operator $\Omega$ can then be
an electronic, photonic or mixed electron-photon operator. The exact
eigenstates of the Hamiltonian are defined by the eigenvalue equation

$H\ket{n}=E_{n}\;\ket{n}\quad n=0,1,2\ldots\;,$ (43)

and are assumed to form a complete basis for the radiation-matter Hilbert
space $\mathcal{V}$. We further assume that at the initial time the system is
in its ground polaritonic state $\ket{0}$. We then apply the following
perturbation operator

$V^{t}=e^{\eta
t}\sum_{n}\bigg{(}V^{\omega_{n}}e^{-i\omega_{n}t}+\big{(}V^{\omega_{n}}\big{)}^{\dagger}e^{i\omega_{n}t}\bigg{)}$
(44)

where $\eta$ is a positive infinitesimal that ensures the perturbation to
vanish for $t\to-\infty$. The time evolution of the initial ground state is
then determined by the time-dependent Hamiltonian

$H+V^{t}$ (45)

and can be modelled by a unitary transformation

$\ket{\tilde{0}}=\text{exp}[iP^{\prime}(t)]\ket{0},$ (46)

where we used a tilde to indicate the time-evolved state and $P^{\prime}(t)$
is a Hermitian operator:

$P^{\prime}(t)=\sum_{n>0}\big{(}P^{\prime}_{n}\ket{n}\bra{0}+P_{n}^{\prime*}\ket{0}\bra{n}\big{)}+(P^{\prime}_{0}+P_{0}^{\prime*})\ket{0}\bra{0}.$
(47)

The ground state phase evolution can be explicitly factorized85

$\ket{\tilde{0}}=e^{i\phi(t)}\;\text{exp}\bigg{(}i\sum_{n>0}\big{(}P_{n}\ket{n}\bra{0}+P_{n}^{*}\ket{0}\bra{n}\big{)}\bigg{)}\ket{0},$
(48)

and if the operator $\Omega$ in Eq. (42) does not involve time
differentiation, it can be disregarded. We will then make use of the phase-
isolated wave function

$\displaystyle\ket{\bar{0}}$
$\displaystyle=\text{exp}\bigg{(}i\sum_{n>0}\big{(}P_{n}\ket{n}\bra{0}+P_{n}^{*}\ket{0}\bra{n}\big{)}\bigg{)}\ket{0}$
(49)

and perform an expansion of the state transfer coefficients in orders of the
perturbation

$P_{n}=P_{n}^{(0)}+P_{n}^{(1)}+P_{n}^{(2)}+P_{n}^{(3)}+\dots$ (50)

which leads to the perturbative expansion of the phase-isolated state

$\displaystyle\ket{\bar{0}}^{(0)}$ $\displaystyle=\ket{0}$ (51)
$\displaystyle\ket{\bar{0}}^{(1)}$
$\displaystyle=i\sum_{n>0}\ket{n}P_{n}^{(1)}$ (52)
$\displaystyle\ket{\bar{0}}^{(2)}$
$\displaystyle=-\frac{1}{2}\ket{0}\sum_{j>0}P_{j}^{(1)}P_{j}^{*(1)}+i\sum_{n>0}\ket{n}P_{n}^{(2)}$
(53) $\displaystyle\ket{\bar{0}}^{(3)}$
$\displaystyle=-\frac{1}{2}\ket{0}\sum_{j>0}\big{(}P_{j}^{(2)}P_{j}^{*(1)}+P_{j}^{(1)}P_{j}^{*(2)}\big{)}+i\sum_{n>0}\ket{n}\big{(}P_{n}^{(3)}-\frac{1}{6}P_{n}^{(1)}\sum_{j>0}P_{j}^{(1)}P_{j}^{*(1)}\big{)}.$
(54)

By defining $P_{-n}\equiv P_{n}^{*}$ and by using in Eq. (42) the state
transfer operators
$\Omega\in\big{\\{}\ket{0}\bra{n}\;;\,\ket{n}\bra{0}\big{\\}}$, we obtain the
following hierarchy of differential equations, for positive and negative
indices84

$\displaystyle i\;\text{\text{sgn}}(k)\dot{P}_{k}^{(1)}-\omega_{k}P_{k}^{(1)}$
$\displaystyle=-iV_{k}^{t}$ (55) $\displaystyle
i\;\text{\text{sgn}}(k)\dot{P}_{k}^{(2)}-\omega_{k}P_{k}^{(2)}$
$\displaystyle=\sum_{n}V^{t}_{k-n}P_{n}^{(1)}$ (56) $\displaystyle
i\;\text{\text{sgn}}(k)\dot{D}_{k}^{(3)}-\omega_{k}D_{k}^{(3)}$
$\displaystyle=\sum_{n}\bigg{[}V^{t}_{k-n}P_{n}^{(2)}+i\theta(kn)V^{t}_{n}P^{(1)}_{-n}P_{k}^{(1)}\bigg{]}+iV^{t}_{k}\sum_{n>0}P^{(1)}_{-n}P_{n}^{(1)},$
(57)

where $\theta(x)$ is the Heaviside step function, while $\omega_{k}$,
$D^{(n)}_{k}$ and $V^{t}_{n}$ are defined as

$\displaystyle\omega_{k}=\omega_{-k}=E_{k}-E_{0}$ (58) $\displaystyle
D_{k}^{(3)}=P^{(3)}_{k}-\frac{2}{3}P_{n}^{(1)}\sum_{j>0}P^{(1)}_{j}P^{(1)}_{-j}$
(59) $\displaystyle V_{k}^{t}=\begin{pmatrix}\braket{k}{V^{t}}{0}\\\
-\braket{0}{V^{t}}{-k}\end{pmatrix};\begin{pmatrix}k>0\\\ k<0\end{pmatrix}$
(60) $\displaystyle
V_{kn}^{t}=\begin{pmatrix}0&\braket{k}{V^{t}}{-n}-\delta_{k,-n}\braket{0}{V^{t}}{0}\\\
\braket{n}{V^{t}}{-k}-\delta_{-k,n}\braket{0}{V^{t}}{0}&0\end{pmatrix};\begin{pmatrix}k,n>0&k>0,n<0\\\
k<0,n>0&k,n<0\end{pmatrix}$ (61)

Notice that the solutions for lower-order coefficients are necessary for
higher-order equations. We can solve these equations, which all have the
structure

$i\;\text{\text{sgn}}(k)\frac{df(t)}{dt}-\omega_{k}f(t)=g(t)$ (62)

with the general solution

$f(t)=-i\;\text{\text{sgn}}(k)e^{-i\omega_{k}\text{\text{sgn}}(k)t}\int_{-\infty}^{t}d\tau
e^{i\omega_{k}\text{\text{sgn}}(k)\tau}g(\tau),$ (63)

which fulfils the vanishing initial conditions for the coefficients.
Expressing the change of the time development average value of an observable
$A$ from its zero-order mean value in orders of the perturbation we obtain

$\displaystyle\delta\braket{A}(t)=A_{V^{t}}(t)-A_{V^{t}}(0)=\braket{\tilde{0}}{A}{\tilde{0}}-\braket{0}{A}{0}$
$\displaystyle=$ $\displaystyle\sum_{n}e^{-i\omega_{n}t+\eta
t}\braket{\braket{A;V^{\omega_{n}}}}_{\omega_{n}+i\eta}$
$\displaystyle+\frac{1}{2}$
$\displaystyle\sum_{mn}e^{-i\omega_{n}t-i\omega_{m}t+2\eta
t}\braket{\braket{A;V^{\omega_{n}},V^{\omega_{m}}}}_{\omega_{n}+i\eta,\omega_{m}+i\eta}+$
$\displaystyle+\frac{1}{6}$
$\displaystyle\sum_{lmn}e^{-i\omega_{l}t-i\omega_{m}t-i\omega_{n}t+3\eta
t}\braket{\braket{A;V^{\omega_{l}},V^{\omega_{m}},V^{\omega_{n}}}}_{\omega_{l}+i\eta,\omega_{m}+i\eta,\omega_{n}+i\eta}+\dots.$
(64)

We then define the linear, quadratic and cubic response functions,
respectively

$\displaystyle\braket{\braket{A;V^{\omega_{n}}}}_{\omega_{n}}$ (65)
$\displaystyle\braket{\braket{A;V^{\omega_{n}},V^{\omega_{m}}}}_{\omega_{n},\omega_{m}}$
(66)
$\displaystyle\braket{\braket{A;V^{\omega_{l}},V^{\omega_{m}},V^{\omega_{n}}}}_{\omega_{l},\omega_{m},\omega_{n}}.$
(67)

These functions express the linear, quadratic and cubic time-variation of an
observable $A$ as resulting from a perturbation in the frequency domain

$\displaystyle\delta\braket{A}^{(1)}(\omega_{1})$
$\displaystyle=\braket{\braket{A;V^{\omega_{1}}}}_{\omega_{1}}$ (68)
$\displaystyle\delta\braket{A}^{(2)}(\omega_{1},\omega_{2})$
$\displaystyle=\frac{1}{2}\braket{\braket{A;V^{\omega_{1}},V^{\omega_{2}}}}_{\omega_{1},\omega_{2}}$
(69) $\displaystyle\delta\braket{A}^{(3)}(\omega_{1},\omega_{2},\omega_{3})$
$\displaystyle=\frac{1}{6}\braket{\braket{A;V^{\omega_{1}},V^{\omega_{2}},V^{\omega_{3}}}}_{\omega_{1},\omega_{2},\omega_{3}},$
(70)

and can therefore be interpreted as molecular (hyper)polarizabilities. For
instance, for an electric dipole perturbation

$V^{\omega}=-\bm{d}\cdot\bm{E}(\omega),$ (71)

the first order variation of the molecular dipole reads

$\delta\braket{\bm{d}}^{(1)}(\omega)=-\braket{\braket{\bm{d};\bm{d}}}_{\omega}\cdot\bm{E}(\omega)$
(72)

and the linear dipole-dipole response function
$\braket{\braket{\bm{d};\bm{d}}}_{\omega}$ can then be interpreted as the
(negative) time-dependent molecular polarizability at frequency $\omega$.
Using the perturbation expansions of the coefficients and the state wave
function Eqs. (50) and (51)-(54), we obtain

$\displaystyle\braket{A}_{V^{t}}=$ $\displaystyle\braket{0}{A}{0}$
$\displaystyle+$ $\displaystyle\braket{0^{(1)}}{A}{0}+\braket{0}{A}{0^{(1)}}$
$\displaystyle+$
$\displaystyle\braket{0^{(1)}}{A}{0^{(1)}}+\braket{0}{A}{0^{(2)}}+\braket{0^{(2)}}{A}{0}$
$\displaystyle+$
$\displaystyle\braket{0^{(3)}}{A}{0}+\braket{0}{A}{0^{(3)}}+\braket{0^{(1)}}{A}{0^{(2)}}+\braket{0^{(2)}}{A}{0^{(1)}}+\dots$
$\displaystyle=$ $\displaystyle\braket{0}{A}{0}$ $\displaystyle-$
$\displaystyle i\sum_{n}A_{-n}P_{n}^{(1)}$ $\displaystyle-$ $\displaystyle
i\sum_{n}A_{-n}P_{n}^{(2)}+\sum_{n,j>0}P_{-j}^{(1)}A_{j-n}P_{n}^{(1)}$
$\displaystyle-i$
$\displaystyle\sum_{n}A_{-n}D^{(3)}_{n}+\sum_{n,k>0}P^{(2)}_{-k}A_{k-n}P_{n}^{(1)}+\sum_{n,k>0}P^{(1)}_{-k}A_{k-n}P_{n}^{(2)}+\dots$
(73)

By comparing Eqs. (5) and (5), using Eq. (63) for the state transfer
parameters, we finally obtain the explicit expressions for the linear,
quadratic and cubic response functions in terms of the eigenstates of the
Hamiltonian84

$\displaystyle\braket{\braket{A;V^{\omega}}}_{\omega}=-\sum_{k}\frac{\text{\text{sgn}}(k)A_{-k}V_{k}^{\omega}}{\omega-\text{\text{sgn}}(k)\omega_{k}}$
(74)
$\displaystyle\braket{\braket{A;V^{\omega_{1}},V^{\omega_{2}}}}_{\omega_{1},\omega_{2}}=\hat{P}(\omega_{1},\omega_{2})\left[-\sum_{k,n}\frac{A_{-k}V^{\omega_{1}}_{k-n}V^{\omega_{2}}_{n}}{(\omega_{1}+\omega_{2}-\text{sgn}(k)\omega_{k})(\omega_{2}-\text{sgn}(n)\omega_{n})}+\sum_{k,n>0}\frac{V^{\omega_{1}}_{-k}A_{k-n}V^{\omega_{2}}_{n}}{(\omega_{1}+\omega_{k})(\omega_{2}-\omega_{n})}\right]$
(75)
$\displaystyle\braket{\braket{A;V^{\omega_{1}},V^{\omega_{2}},V^{\omega_{3}}}}_{\omega_{1},\omega_{2},\omega_{3}}=\hat{P}(\omega_{1},\omega_{2},\omega_{3})\left\\{\sum_{k}\frac{A_{-k}\text{sgn}(k)}{\omega_{1}+\omega_{2}+\omega_{3}-\text{sgn}(k)\omega_{k}}\right.\;\times\left[\sum_{n}\frac{\theta(kn)V^{\omega_{1}}_{n}V^{\omega_{2}}_{-n}V^{\omega_{3}}_{k}}{(\omega_{2}+\text{sgn}(n)\omega_{n})(\omega_{3}-\text{sgn}(k)\omega_{k})}\right.$
$\displaystyle\left.-\sum_{n,m}\frac{V^{\omega_{1}}_{k-n}V^{\omega_{2}}_{n-m}V^{\omega_{3}}_{m}}{(\omega_{2}+\omega_{3}-\text{sgn}(n)\omega_{n})(\omega_{3}-\text{sgn}(m)\omega_{m})}+\sum_{n>0}\frac{V^{\omega_{1}}_{k}V^{\omega_{2}}_{n}V_{-n}^{\omega_{3}}}{(\omega_{2}-\text{sgn}(n)\omega_{n})(\omega_{3}+\text{sgn}(n)\omega_{n})}\right]$
$\displaystyle\left.-\sum_{m,n,k>0}\left(\frac{V^{\omega_{1}}_{-km}V^{\omega_{2}}_{-m}V_{n}^{\omega_{3}}A_{k-n}}{(\omega_{1}+\omega_{2}+\omega_{k})(\omega_{2}+\omega_{m})(\omega_{3}-\omega_{n})}-\frac{V^{\omega_{1}}_{-k}V^{\omega_{2}}_{n-m}V_{m}^{\omega_{3}}A_{k-n}}{(\omega_{2}+\omega_{3}-\omega_{n})(\omega_{1}+\omega_{k})(\omega_{3}-\omega_{m})}\right)\right\\}$
(76)

where $\hat{P}(\omega_{1},\omega_{2})$ and
$\hat{P}(\omega_{1},\omega_{2},\omega_{3})$ sum all the permutations of the
frequencies

$\displaystyle\hat{P}(\omega_{1},\omega_{2})f(\omega_{1},\omega_{2})$
$\displaystyle=f(\omega_{1},\omega_{2})+f(\omega_{2},\omega_{1})$ (77)
$\displaystyle\hat{P}(\omega_{1},\omega_{2},\omega_{3})f(\omega_{1},\omega_{2},\omega_{3})$
$\displaystyle=f(\omega_{1},\omega_{2},\omega_{3})+f(\omega_{1},\omega_{3},\omega_{2})$
$\displaystyle+f(\omega_{2},\omega_{1},\omega_{3})+f(\omega_{2},\omega_{3},\omega_{1})$
$\displaystyle+f(\omega_{3},\omega_{1},\omega_{2})+f(\omega_{3},\omega_{2},\omega_{1}),$
(78)

and we set the switching parameter $\eta$ to zero.

These equations are identical to the expressions obtained in molecular
response theory,84, 85, 86, 87, 88, 89, 90, 91, 92, 93 since they only rely on
the Schrödinger equation. However, the explicit treatment of internal
electromagnetic fields as dynamical variables in our system leads to novel
perspectives. First of all, we must remember that the eigenstates of the
system are mixed matter-photon states, so they belong to a larger Hilbert
space than the usual molecular one. The properties of these field-dressed
states will therefore differ from the bare-molecule ones, leading, for
instance, to different excitation energies and transition moments. Moreover,
additional possibilities for probing the properties of the system arise since
the external perturbation $V^{t}$ can now act on both matter and photon
degrees of freedom. Furthermore, considering the internal electromagnetic
fields allows us to explore additional observables connected to the radiation
fields and their interactions with matter.

The linear response function Eq. (74) describes the first-order variation of
an observable $A$ and it is connected to the linear polarizabilities of the
system. Its poles occur at the polaritonic excitation energies, and the
corresponding residues are connected to the transition moments between the
ground and the excited states

$\displaystyle\lim_{\omega\to\omega_{n}}(\omega-\omega_{n})\braket{\braket{A;B}}_{\omega}$
$\displaystyle=\braket{0}{A}{n}\braket{n}{B}{0}$ (79) $\displaystyle\equiv
A_{0n}B_{n0}$
$\displaystyle\lim_{\omega\to-\omega_{n}}(\omega+\omega_{n})\braket{\braket{A;B}}_{\omega}$
$\displaystyle=-\braket{0}{B}{n}\braket{n}{A}{0}$ (80) $\displaystyle\equiv-
B_{0n}A_{n0}.$

The quadratic response function Eq. (75) is connected to the first
hyperpolarizabilities of the system since it describes the second-order
response to an external perturbation. Its residues are connected to the
transition moments and matrices between the excited states. These can be
obtained once the transition moments between ground and excited states have
been computed

$\displaystyle\lim_{\omega_{2}\to\omega_{m}}(\omega_{2}-\omega_{m})\braket{\braket{A;B,C}}_{\omega_{1},\omega_{2}}=\sum_{k>0}\bigg{[}\frac{\braket{0}{A}{k}\big{(}\braket{k}{B}{m}-\delta_{mk}\braket{0}{B}{0}\big{)}}{(\omega_{1}+\omega_{m}-\omega_{k})}-\frac{\braket{0}{B}{k}\big{(}\braket{k}{A}{m}-\delta_{mk}\braket{0}{A}{0}\big{)}}{\omega_{1}+\omega_{k}}\bigg{]}\braket{m}{C}{0}$
(81)
$\displaystyle\lim_{\omega_{1}\to-\omega_{q}}(\omega_{1}+\omega_{q})\lim_{\omega_{2}\to\omega_{m}}(\omega_{2}-\omega_{m})\braket{\braket{A;B,C}}_{\omega_{1},\omega_{2}}=-\braket{0}{B}{q}\big{(}\braket{q}{A}{m}-\delta_{mq}\braket{0}{A}{0}\big{)}\braket{m}{C}{0}.$
(82)

These expressions also give us access to excited state properties such as the
dipole moment. Finally, the cubic response function describes the third-order
variations of an observable as a consequence of an external perturbation, and
it is therefore connected to the second hyperpolarizabilities of the system.

The response functions exhibit the following symmetry relations:

$\displaystyle\braket{\braket{A;B}}_{\omega}=$
$\displaystyle\braket{\braket{B;A}}_{-\omega}$ $\displaystyle=$
$\displaystyle\braket{\braket{A;B}}_{-\omega}^{*}$ (83)
$\displaystyle\braket{\braket{A;B,C}}_{\omega_{1},\omega_{2}}=$
$\displaystyle\braket{\braket{A;C,B}}_{\omega_{2},\omega_{1}}$
$\displaystyle=$
$\displaystyle\braket{\braket{C;A,B}}_{-\omega_{1}-\omega_{2},\omega_{1}}$
$\displaystyle=$
$\displaystyle\braket{\braket{A;B,C}}_{-\omega_{1},-\omega_{2}}^{*}$ (84)
$\displaystyle\braket{\braket{A;B,C,D}}_{\omega_{1},\omega_{2},\omega_{3}}=$
$\displaystyle\braket{\braket{A;C,B,D}}_{\omega_{2},\omega_{1}\omega_{3}}$
$\displaystyle=$
$\displaystyle\braket{\braket{D;A,B,C}}_{-\omega_{1}-\omega_{2}-\omega_{3},\omega_{1},\omega_{2}}$
$\displaystyle=$
$\displaystyle\braket{\braket{A;B,C,D}}_{-\omega_{1},-\omega_{2},-\omega_{3}}^{*},$
(85)

for Hermitian operators $A$, $B$, $C$, and $D$.

#### Time-independent limit

The static response theory can be obtained from the previous derivation by
setting the external frequencies $\omega$ and the switching parameter $\eta$
to zero. Since the perturbation is now time-independent, it is possible to
relate the static polaritonic properties to the derivatives of the Hamiltonian
eigenvalues. The perturbation $V$ acting on the system is time-independent and
supposed to be smooth and small, and it is usually expressed as an external
field $\bm{F}$ acting on a molecular multipole $\bm{x}$

$V=-\bm{x}\cdot\bm{F}.$ (86)

The perturbed Hamiltonian then reads

$\displaystyle H=H_{0}+V=H_{0}-\bm{x}\cdot\bm{F}.$ (87)

The Hamiltonian in Eq. (87) is endowed with well-defined eigenstates and
eigenvalues, whose small changes compared to the unperturbed eigenvalues and
eigenfunctions are commonly studied via the Rayleigh-Schrödinger (RS)
perturbation theory.203 From the eigenfunctions $\ket{n}$ of the unperturbed
Hamiltonian $H_{0}$ (see Eq. (43)), the ground-state eigenfunction
$\ket{\tilde{0}}$ and the eigenvalue $\tilde{E}$ of the perturbed operator are
expressed as a power series:

$\displaystyle\ket{\tilde{0}}$
$\displaystyle=\ket{0}+\ket{\tilde{0}^{(1)}}+\ldots$ (88)
$\displaystyle\tilde{E}$ $\displaystyle=E_{0}+E^{(1)}+E^{(2)}+\ldots$ (89)

The first- and second-order corrections to the energy read:

$\displaystyle E^{(1)}=\bra{0}V\ket{0}$ (90) $\displaystyle
E^{(2)}=\sum_{k=1}^{\infty}\frac{\bra{0}V\ket{k}\bra{k}V\ket{0}}{E_{0}-E_{k}}=\frac{1}{2}\braket{\braket{V;V}}_{\omega=0}.$
(91)

Eq. (90) is the well-known Helmann-Feynman theorem204, 205 and only requires
the knowledge of the unperturbed wave function, while Eq. (91) also includes a
sum-over-states contribution from each excited state $\ket{k}$ with energy
$E_{k}$. The derivatives of the energy with respect to the external fields
provide us with the ground state molecular properties. From Eqs. (86) and
(90), we then identify the derivative of the energy with respect to the field
$\bm{F}$ as the ground state mean multipole value

$-\frac{d\tilde{E}}{d\bm{F}}\bigg{|}_{\bm{F}=0}=\braket{0}{\bm{x}}{0}$ (92)

From Eq. (91), the second derivative of the energy is then interpreted as the
multipole polarizability to the field $\bm{F}$, and from Eq. (74) we identify
the linear response function for $\omega=0$

$\frac{d^{2}\tilde{E}}{d^{2}\bm{F}}\bigg{|}_{\bm{F}=0}=2\sum_{k=1}^{\infty}\frac{\bra{0}\bm{x}\ket{k}\bra{k}\bm{x}\ket{0}}{E_{0}-E_{k}}=\braket{\braket{\bm{x};\bm{x}}}_{\omega=0}.$
(93)

The second derivatives of the energy are then connected to the molecular
polarizabilities. Higher-order derivatives relate hyperpolarizabilities and
nonlinear response functions. Notice that these expressions are the same as
standard molecular static response theory, and they also hold for exact
polaritonic states since they are based only on the Schrödinger equation.

###  5.1 On the rotational average of computed molecular properties

In this section, we discuss the problem of the orientational average of the
computed polarizabilities, which is straightforwardly accomplished in standard
molecular response theory to connect the single-molecule calculation to the
macroscopic sample response. Since we are now considering states that belong
to the extended light-matter Hilbert space in Eq. (29), we have to consider
the molecule and the electromagnetic environment, which are defined by the
optical device. In general, this introduces anisotropy in the system. The
cavity field can also break the molecular symmetry, possibly allowing for
otherwise symmetry-forbidden transitions. This anisotropy-symmetry breaking
introduces further complications in the computation of molecular properties.

The bare molecular states are independent of the spatial orientation of the
molecule. The energy is invariant for rigid-body rotations of the system,
while the charge density rotates following the molecular structure, and so do
properties such as molecular polarizabilities. However, the strong coupling
between molecular and cavity states depends on the projection of the molecular
dipole operator onto the cavity field, as can be seen in the electron-photon
interaction term in Hamiltonian Eq. (4.2). This also means that there is a
non-trivial energy dependence on the relative molecular orientation in the
cavity field. Excited states will be particularly affected by this orientation
dependence. A bare electronic excitation with a transition moment orthogonal
to the cavity field will have a negligible direct coupling with the photonic
states. On the other hand, if the transition dipole moment is aligned with the
field polarization, there will be mixing, giving rise to intense polaritonic
excitations of hybrid light-matter character. This is illustrated in Fig. (3),
where we plot the upper and lower polaritonic energies of a hydrogen molecule
as a function of the angle between the transition dipole moment and the cavity
field.

Figure 3: Potential energy surfaces for a hydrogen molecule’s upper and lower
polaritons resonantly coupled with a quantum cavity, as a function of the
molecular orientation. The energies are computed by using the time-dependent
quantum electrodynamic Hartree-Fock method with a cc-pVDZ basis set and
coupling constant set to $\lambda=0.01$ a.u. (see section (6)). When the
molecule’s transition dipole moment is aligned with the cavity field (90°),
the Rabi splitting is maximal, while it decreases to zero when the transition
moment is perpendicular to the field’s polarization (180°).

These considerations on the molecular orientation pose some issues for the
computation of molecular properties, as the direction of the cavity field is
fixed in space by the experimental setup. In contrast, inside the cavity, the
molecules are usually randomly or quasi-randomly oriented. Following a
chemical approach, we could compute the sample’s properties using the
frequency-dependent polarizability Eq. (74), which must be averaged over the
different molecular orientations206, 193. The molecular orientation will be
defined by a general set of parameters $\Omega$. For instance, it can be
uniquely defined by the three Euler’s angles $(\phi,\theta,\chi)$. Notice that
a different set of parameter can be convenient depending on the setup. For
instance, for a Fabry-Pérot resonator with two degenerate modes with
wavevector $\mathbf{k}$ and perpendicular polarizations, a rotation of the
molecule along $\mathbf{k}$ does not change the Hamiltonian. Then it is
sufficient to use only two parameters instead of the three Euler’s angles. If
we consider the Boltzmann weight of each different orientation, we have

$\overline{\braket{\braket{A\;;B}}}_{\omega}=\int_{\Omega}e^{-\frac{E_{0}(\Omega)}{k_{B}T}}\;\sum_{k}\bigg{(}\frac{A_{0k}(\Omega)B_{k0}(\Omega)}{\omega-\omega_{k}(\Omega)}-\frac{B_{0k}(\Omega)A_{k0}(\Omega)}{\omega+\omega_{k}(\Omega)}\bigg{)}d\Omega\bigg{/}\int_{\Omega}e^{-\frac{E_{0}(\Omega)}{k_{B}T}}d\Omega,$
(94)

where $\overline{\braket{\braket{A\;;B}}}_{\omega}$ is the averaged response
function, $E_{0}(\Omega)$ is the ground state energy of the polaritonic system
as a function of the molecular orientation, $k_{B}$ is the Boltzmann constant,
and $T$ is the temperature. Even if we disregard the Boltzmann weight and
perform an isotropic average, the result does not simply lead to an isotropic
polarization tensor. The excitation energies in the denominator depend on the
orientation, as well as the numerator. The polarization will be generally
anisotropic because of the introduction of the optical device.

As a concrete example, we can consider the computation of absorption spectra.
The linear absorption spectrum is connected to the dipole-dipole
polarizability, where $A$ and $B$ in Eq. (74) are components of the dipole
operator $d_{i}$84, 85, 206

$\displaystyle\braket{\braket{d_{i}\;;d_{j}}}_{\omega}=\sum_{k>0}\bigg{(}\frac{\braket{0}{d_{i}}{k}\braket{k}{d_{j}}{0}}{\omega-\omega_{k}}-\frac{\braket{0}{d_{j}}{k}\braket{k}{d_{i}}{0}}{\omega+\omega_{k}}\bigg{)}.$

Neglecting the dependence of the polaritons on the molecular orientation, we
obtain the usual expression for absorption by randomly oriented molecules in
terms of the oscillator strength

$f_{k0}=\frac{2}{3}\;\omega_{k}\sum_{i=x,y,z}|\braket{0}{d_{i}}{k}|^{2}.$ (95)

This corresponds to fixing the relative molecule-field orientation. A
fictitious peak-broadening is often applied to the computed spectrum
$S(\omega)$, e.g. a Lorentzian lineshape with bandwidth broadening parameter
$\Delta$

$\Gamma(\omega,\omega_{k})=\frac{1}{\pi}\;\frac{\Delta}{(\omega-\omega_{k})^{2}+\Delta^{2}}$
(96)

as reported in Refs.94, 95 so that the computed spectrum is

$S(\omega)=\sum_{k}f_{k0}\;\Gamma(\omega,\omega_{k}).$ (97)

However, different orientations could play an important role in determining
the properties of the sample, and the choice of a fixed relative orientation
does not give the full picture. For instance, if we are interested in a
specific transition, we could suppose the field polarization to be parallel to
the transition dipole moment. However, this could lead to biased results as
certain effects could be enhanced while others suppressed.

The study of orientational disorder in ab initio polaritonic chemistry is
still to be addressed. A possible approach would be to perform molecular
dynamics simulations that sample the different orientations of the molecule
inside the cavity and perform an average of the spectra obtained for each
snapshot. This corresponds to an incoherent superposition of the spectra of an
ensemble. The appearance of collective effects can introduce further
complications in the computation of molecular properties. If we could perform
an ab initio simulation including a sufficiently large number of molecules
randomly or quasi-randomly oriented, the corresponding computed properties
would include both collective and orientational effects. Alternatively,
collective effects can initially be addressed with a simplified model, such as
the Tavis-Cumming (TC) model,68, 67 which considers an ensemble of identical
two-level non-interacting systems. Moreover, it disregards the dipole self-
energy, the molecular dipoles and the so-called counter-rotating terms of the
Hamiltonian. The Rabi splitting is then predicted to scale with the square
root of the number of molecules and, therefore, it is sufficient to employ a
smaller coupling strength compared to single-molecule calculations. The
generalization of the TC model to an energy-broadened set of two-level
molecules207 predicts a Rabi splitting dependent on the quadratic average of
the coupling strengths and on the energy broadening. The predicted spectrum
would be very different from the one obtained by an average of single-molecule
absorption spectra.

The interplay between single-molecule and collective effects will need a
careful comparison between simplified models4, 2, 68, 67 and more accurate
simulations, and recently some methods to study polaritonic chemistry in an
impurity-like approach by an effective embedding have been proposed.200, 201,
208, 209, 210, 211 We believe both these effects are important. However, which
one dominates will likely depend on the experimental conditions. Both should
be carefully considered when accurately modelling polaritonic properties, and
further studies in this direction are needed.

###  5.2 Equivalent expressions for polaritonic properties

In this section, we derive equations of motion for polaritonic response
functions. Inserting the expansion Eq. (5) into Ehrenfest theorem Eq. (42), we
obtain a set of equations as in standard molecular response theory84

$\displaystyle\omega_{1}\braket{\braket{A;V^{\omega_{1}}}}_{\omega_{1}}$
$\displaystyle=\braket{\braket{[A,H];V^{\omega_{1}}}}_{\omega_{1}}+\braket{0}{[A,V^{\omega_{1}}]}{0}$
(98)
$\displaystyle(\omega_{1}+\omega_{2})\braket{\braket{A;V^{\omega_{1}},V^{\omega_{2}}}}_{\omega_{1},\omega_{2}}$
$\displaystyle=\braket{\braket{[A,H];V^{\omega_{1}},V^{\omega_{2}}}}_{\omega_{1},\omega_{2}}+\hat{P}(\omega_{1},\omega_{2})\braket{\braket{[A,V^{\omega_{1}}];V^{\omega_{2}}}}_{\omega_{2}}$
(99)

Considering Eq. (98) and the position-momentum operator identity

$\bm{p}=-i[\bm{r},H^{l}]$ (100)

which holds for the dipole Hamiltonian in the length representation Eq. (4.2),
we obtain the relation

$\omega\braket{\braket{\bm{r};\bm{r}}}_{\omega}=i\braket{\braket{\bm{p};\bm{r}}}_{\omega},$
(101)

which states that the frequency-dependent dipole polarizability can be
evaluated regarding the position or the conjugate momentum. This ensures the
equivalence of expressions in the dipole and velocity form of transition
dipole moments

$\omega_{n}\braket{0}{\bm{r}}{n}=i\braket{0}{\bm{p}}{n},$ (102)

where $\ket{n}$ is an excited state and $\omega_{n}$ is the excitation energy.
We note that for the velocity representation, in Eq. (4.2), the analogous
relation is

$\bm{p}+\frac{1}{c}\mathbf{A}(\bm{r})=-i[\bm{r},H^{v}],$ (103)

which states, exactly as Eq. (102), an equivalence between the transition
dipole moment and the transition kinetic momentum of the electrons. Indeed,
Eqs. (100) and (103) are connected by the transformation in Eq. (33). We
recall that in second quantization, Eqs. (100) and (103) are fulfilled only in
the limit of a complete basis set.185

An example of a relation introduced by the electromagnetic degrees of freedom
is provided by the commutation relation between the Hamiltonian and the
annihilation operators in the length form

$[b_{\alpha},H^{l}]=\lambda_{\alpha}\sqrt{\frac{\omega_{\alpha}}{2}}(\bm{d}\cdot\bm{\epsilon}_{\alpha})+\omega_{\alpha}b_{\alpha},$
(104)

which through Eq. (98) leads to

$\displaystyle\omega\braket{\braket{b_{\alpha};B}}_{\omega}=$
$\displaystyle\langle{\langle{\lambda_{\alpha}\sqrt{\frac{\omega_{\alpha}}{2}}\bm{d}\cdot\bm{\epsilon}_{\alpha}+\omega_{\alpha}b_{\alpha};B}}\rangle\rangle_{\omega}$
$\displaystyle+\braket{0}{[b_{\alpha},B]}{0}$ (105)

for any operator $B$. We obtain from the residues of Eq. (5.2) the relation

$(\omega_{n}-\omega_{\alpha})\braket{0}{b_{\alpha}}{n}=\lambda_{\alpha}\sqrt{\frac{\omega_{\alpha}}{2}}\braket{0}{\bm{d}\cdot\bm{\epsilon}_{\alpha}}{n}$
(106)

which allows the computation of photonic transition moments in terms of matter
transition moments. This reflects the entanglement between electronic and
photonic degrees of freedom. The analogous relation for the creation operators
is

$(\omega_{n}+\omega_{\alpha})\braket{0}{b^{\dagger}_{\alpha}}{n}=-\lambda_{\alpha}\sqrt{\frac{\omega_{\alpha}}{2}}\braket{0}{\bm{d}\cdot\bm{\epsilon}_{\alpha}}{n}.$
(107)

Moreover, from the relation

$[b_{\alpha}+b_{\alpha}^{\dagger},H^{l}]=\omega_{\alpha}(b_{\alpha}-b^{\dagger}_{\alpha})$
(108)

we obtain

$\displaystyle\omega_{n}\braket{0}{b_{\alpha}+b^{\dagger}_{\alpha}}{n}$
$\displaystyle=\omega_{\alpha}\braket{0}{b_{\alpha}-b_{\alpha}^{\dagger}}{n}$
(109)
$\displaystyle=\lambda_{\alpha}\sqrt{\frac{\omega_{\alpha}}{2}}\braket{0}{\bm{d}\cdot\bm{\epsilon}_{\alpha}}{n}\frac{2\omega_{n}\omega_{\alpha}}{\omega_{n}^{2}-\omega_{\alpha}^{2}}$
(110)

If we introduce the conjugate field displacement $q_{\alpha}$ and momentum
$p_{\alpha}$ coordinates

$\displaystyle p_{\alpha}$
$\displaystyle=i\sqrt{\frac{\omega_{\alpha}}{2}}(b_{\alpha}^{\dagger}-b_{\alpha})$
(111) $\displaystyle q_{\alpha}$
$\displaystyle=\frac{1}{\sqrt{2\omega_{\alpha}}}(b_{\alpha}^{\dagger}+b_{\alpha}),$
(112)

we can write Eq. (109) as

$\omega_{n}\braket{0}{q_{\alpha}}{n}=i\braket{0}{{p}_{\alpha}}{n}$ (113)

which is the equivalent relation of Eq. (102) for the photonic conjugate
momenta. Similarly, by employing the Hamiltonian in the velocity form in Eq.
(4.2) or through the transformation in Eq. (33), we obtain equivalent
relations in the velocity representation. Notice that relation Eq. (109) also
holds for the dipole Hamiltonian in the velocity form, but the physical
meaning of the operators (111) and (112) is changed since they refer to a
different representation. This will be discussed extensively in the following
sections.

Analogous relations can be obtained from the higher-order equations of motion
in Eqs. (99), for instance for the transition dipole moments among excited
states we have

$\displaystyle-i(\omega_{m}-\omega_{n})\big{(}$
$\displaystyle\braket{n}{r_{i}}{m}-\delta_{mn}\braket{0}{r_{i}}{0}\big{)}=$
$\displaystyle\braket{n}{p_{i}}{m}-\delta_{mn}\braket{0}{p_{i}}{0}.$ (114)

Notice that relations similar to Eq. (5.2) also hold for the photonic
conjugate momenta:

$\displaystyle-i(\omega_{m}-\omega_{n})\big{(}$
$\displaystyle\braket{n}{q_{\alpha}}{m}-\delta_{mn}\braket{0}{q_{\alpha}}{0}\big{)}=$
$\displaystyle\braket{n}{p_{\alpha}}{m}-\delta_{mn}\braket{0}{p_{\alpha}}{0}.$
(115)

###  5.3 Response functions in cavity QED

This section provides a discussion of response functions in cavity QED.
Although far from exhaustive, we provide the reader with several examples of
matter-photon and photon-photon response functions. We focus on the
peculiarities introduced by the photon dressing, and we explicitly discuss the
result of using different mathematical representations of the Hamiltonian.

The perturbation is described by a single frequency component of Eq. (44)

$V^{t}=\bigg{(}V^{\omega}\;e^{-i\omega t}+V^{-\omega}e^{i\omega
t}\bigg{)}\;e^{\eta t},$ (116)

where this also includes the case of static perturbations by imposing
$\omega=\eta=0$.

####  5.3.1 External electromagnetic fields

Spectroscopic techniques probe the system by means of an external
electromagnetic field, whose electric and magnetic components are coupled to
the motion of particles. The external probe is treated as classical (non-
quantized) fields, described by its own vector $\mathbf{A}_{e}(\bm{r},t)$ and
scalar $\phi_{e}(\bm{r},t)$ potentials. It is not necessary to describe the
internal and external fields through the same gauge. The Hamiltonian of the
system in the Coulomb gauge and Born-Oppenheimer approximation reads162

$\displaystyle H=\sum_{i}\frac{1}{2}$
$\displaystyle\bigg{(}\bm{p}_{i}-\frac{1}{c}\mathbf{A}(\bm{r}_{i})-\frac{1}{c}\mathbf{A}_{e}(\bm{r}_{i})\bigg{)}^{2}+V_{coul}$
$\displaystyle-\sum_{i}\phi_{e}(\bm{r}_{i})+\sum_{\alpha}\omega_{\alpha}b^{\dagger}_{\alpha}b_{\alpha},$
(117)

where $\alpha$ labels the photon modes and $i$ refers to electrons. Expanding
the first term, we get the polaritonic Hamiltonian in Eq. (3) and the
interaction term $V^{t}$ with the external fields:

$\displaystyle V^{t}=$
$\displaystyle\frac{1}{c^{2}}\sum_{i}\mathbf{A}(\bm{r}_{i})\cdot\mathbf{A}_{e}(\bm{r}_{i})$
$\displaystyle+\frac{1}{2c^{2}}\sum_{i}\mathbf{A}^{2}_{e}(\bm{r}_{i})-\sum_{i}\frac{1}{2c}\big{(}\mathbf{A}_{e}(\bm{r}_{i})\cdot\bm{p}_{i}+\bm{p}_{i}\cdot\mathbf{A}_{e}(\bm{r}_{i})\big{)}$
$\displaystyle-\sum_{i}\phi_{e}(\bm{r}_{i}).$ (118)

While the last three terms of Eq. (5.3.1) are also found in molecular response
theory in the semi-classical approximation,206, 194, 192, 195 the first term
is a purely QED contribution. This term cancels the field contribution from
the $\bm{p}$ terms and ensures the coupling is only to the matter subsystem.
When we perform the dipole approximation for the cavity fields and apply the
transformation in Eq. (33), we obtain the length Hamiltonian $H^{l}$. The
first term in Eq. (5.3.1) is cancelled, leaving only the familiar interaction
terms. By a suitable expansion of the external potential, the interaction term
is obtained in a multipolar fashion as for the standard semiclassical
theory206, 192, 193, 194

$\displaystyle H_{mult}=$
$\displaystyle-\bm{d}\cdot\mathbf{E}_{e}-\bm{m}\cdot\mathbf{B}_{e}$
$\displaystyle-\frac{1}{6}\Theta_{\alpha\beta}\partial_{\alpha}(E_{e})_{\beta}-\frac{1}{2}(B_{e})_{\alpha}(B_{e})_{\beta}\chi_{\alpha\beta}+\dots$
(119)

where we defined respectively the electric dipole $\bm{d}$, magnetic dipole
$\bm{m}$, electric quadrupole $\Theta_{\alpha\beta}$ and magnetic
susceptibility $\chi_{\alpha\beta}$ operators (in atomic units):

$\displaystyle\bm{d}$
$\displaystyle=\sum_{M}Z_{M}\bm{R}_{M}-\sum_{i}\bm{r}_{i}$ (120)
$\displaystyle\Theta_{\alpha\beta}$
$\displaystyle=\sum_{i}(3r_{i\alpha}r_{i\beta}-\delta_{\alpha\beta}r_{i}^{2})$
(121) $\displaystyle\bm{m}$
$\displaystyle=\sum_{i}\frac{1}{2c}\bm{l}_{i}=\sum_{i}\frac{1}{2c}\bm{r}_{i}\times\bm{p}_{i}$
(122) $\displaystyle\chi_{\alpha\beta}$
$\displaystyle=\sum_{i}\frac{1}{4c^{2}}(r_{i\alpha}r_{i\beta}-\delta_{\alpha\beta}r_{i}^{2}).$
(123)

Note that since we employ the dipole approximation in the length form, these
multipolar operators refer exclusively to the matter subsystem and have the
same physical meaning as in standard response theory. On the contrary, in the
velocity form the operators

$\bm{m}=\sum_{i}\frac{1}{2c}\bm{l}_{i}=\sum_{i}\frac{1}{2c}\bm{r}_{i}\times\bm{p}_{i}$
(124)

has a mixed matter-photon character because of the field component in the
conjugate momentum $\bm{p}$ (see Eq. (21)).

If we only retain the dipole interaction term in (5.3.1), the external
electric field $\mathbf{E}_{e}(0)$ computed at the molecular position is a
constant that can be factorized out of the response functions. The dipole-
dipole response function $\braket{\braket{d_{i}\;;d_{j}}}_{\omega}$ can then
be interpreted as the negative molecular polarizability at frequency $\omega$

$\displaystyle\braket{\braket{d_{i}\;;d_{j}}}_{\omega}=\sum_{k>0}\bigg{(}\frac{\braket{0}{d_{i}}{k}\braket{k}{d_{j}}{0}}{\omega-\omega_{k}}-\frac{\braket{0}{d_{j}}{k}\braket{k}{d_{i}}{0}}{\omega+\omega_{k}}\bigg{)}$
(125)

where $i$ and $j$ refer to cartesian components of the dipole operator
$\bm{d}$ and $k$ labels excited polaritonic states. This function describes
the first-order time evolution of the molecular dipole moment subject to an
external electric field:

$\delta\braket{\bm{d}}^{(1)}(t)=-\int d\omega\;e^{-i\omega t+\eta
t}\sum_{j=x,y,z}\braket{\braket{\bm{d}\,;d_{j}}}_{\omega+i\eta}E_{j}(\omega).$
(126) $V^{t}=-\bm{d}\cdot\mathbf{E}(t)=-\bm{d}\cdot\bigg{(}e^{\eta t}\int
d\omega\;e^{-i\omega t}\mathbf{E}(\omega)\bigg{)}$ (127)

The poles occur when the frequency $\omega$ matches the energy difference
between the ground and excited states, and the corresponding residues

$\displaystyle\lim_{\omega\to\omega_{k}}(\omega-\omega_{k})\braket{\braket{d_{i}\;;d_{j}}}_{\omega}$
$\displaystyle=\braket{0}{d_{i}}{k}\braket{k}{d_{j}}{0}$ (128)
$\displaystyle\lim_{\omega\to-\omega_{k}}(\omega+\omega_{k})\braket{\braket{d_{i}\;;d_{j}}}_{\omega}$
$\displaystyle=-\braket{0}{d_{j}}{k}\braket{k}{d_{i}}{0}$ (129)

provide information on the transition dipole moments between the ground and
the excited states. While these expressions are the same as in standard
molecular response theory,84 the states involved are here polaritonic and the
predicted properties will consequently be modified.

If we replace one of the electric dipole operators with the magnetic dipole
$\bm{m}$ we obtain the electric dipole-magnetic dipole response function

$\displaystyle\braket{\braket{d_{i}\;;m_{j}}}_{\omega}=\sum_{k>0}\bigg{(}\frac{\braket{0}{d_{i}}{k}\braket{k}{m_{j}}{0}}{\omega-\omega_{k}}-\frac{\braket{0}{m_{j}}{k}\braket{k}{d_{i}}{0}}{\omega+\omega_{k}}\bigg{)},$
(130)

which describes the electric dipole response when a magnetic field is applied

$\delta\braket{\bm{d}}^{(1)}(t)=-\int d\omega\;e^{-i\omega t+\eta
t}\sum_{j=x,y,z}\braket{\braket{\bm{d}\,;m_{j}}}_{\omega+i\eta}B_{j}(\omega)$
(131) $V^{t}=-\bm{m}\cdot\mathbf{B}(t)=-\bm{m}\cdot\bigg{(}e^{\eta t}\int
d\omega\;e^{-i\omega t}\mathbf{B}(\omega)\bigg{)}.$ (132)

From the residues of the frequency-dependent electric dipole-magnetic dipole
polarizability in Eq. (130), we obtain the rotational strength of optically
active molecules

$\displaystyle\lim_{\omega\to\omega_{k}}(\omega-\omega_{k})\braket{\braket{d_{i}\;;m_{j}}}_{\omega}$
$\displaystyle=\braket{0}{d_{i}}{k}\braket{k}{m_{j}}{0}$ (133)
$\displaystyle\lim_{\omega\to-\omega_{k}}(\omega+\omega_{k})\braket{\braket{d_{i}\;;m_{j}}}_{\omega}$
$\displaystyle=-\braket{0}{m_{j}}{k}\braket{k}{d_{i}}{0}.$ (134)

As previously discussed, the operator $\bm{m}$ can only in the length
representation be interpreted as the molecular magnetic dipole operator, while
it has a different physical meaning in the velocity Hamiltonian. Notice that
in polaritonic systems, the optical activity can appear both from molecular
and field chirality. As discussed in Sec. (2), several groups have recently
developed chiral cavities,147, 148, 149, 150, 115, 151, 152, 153, 154, 155,
156, 157, 158, 159, 160, 161 which can host only one-handedness of light
polarization within themselves. This opens up the possibility of engineering
chiral properties through the chirality of the photon field,212, 36, 37, 213,
214, 147 which is transferred to molecular states via light-matter strong
coupling. For instance, we expect it to induce chirality in a non-optically
active molecule, similar to induced circular dichroism (ICD) when achiral
molecules are placed in a chiral solvent.215, 216, 217, 218 On the other hand,
chiral molecules interacting with a chiral field should give rise to
"polaritonic diastereoisomers", possibly leading to novel approaches to chiral
discrimination.

As the states of the system here belong to the radiation-matter Hilbert space
in Eq. (29), we can study how the perturbation of the matter subsystem leads
to modifications of the photonic properties. If we focus on the perturbation
in Eq. (127), in the velocity representation, the response function

$\displaystyle\braket{\braket{b^{\dagger}_{\alpha}b_{\alpha}\;;d_{j}}}_{\omega}=\sum_{k>0}\bigg{(}\frac{\braket{0}{b^{\dagger}_{\alpha}b_{\alpha}}{k}\braket{k}{d_{j}}{0}}{\omega-\omega_{k}}-\frac{\braket{0}{d_{j}}{k}\braket{k}{b^{\dagger}_{\alpha}b_{\alpha}}{0}}{\omega+\omega_{k}}\bigg{)}$
(135)

describes the time evolution of the photon number in the $\alpha$-mode

$\delta\braket{b^{\dagger}_{\alpha}b_{\alpha}}^{(1)}=-\int
d\omega\;e^{-i\omega t+\eta
t}\sum_{j}\braket{\braket{b^{\dagger}_{\alpha}b_{\alpha}\,;d_{j}}}_{\omega+i\eta}E_{j}(\omega).$
(136)

The residues computed at the polaritonic excitations

$\displaystyle\lim_{\omega\to\omega_{k}}(\omega-\omega_{k})\braket{\braket{b^{\dagger}_{\alpha}b_{\alpha}\;;d_{j}}}_{\omega}$
$\displaystyle=\braket{0}{b^{\dagger}_{\alpha}b_{\alpha}}{k}\braket{k}{d_{j}}{0}$
(137)
$\displaystyle\lim_{\omega\to-\omega_{k}}(\omega+\omega_{k})\braket{\braket{b^{\dagger}_{\alpha}b_{\alpha}\;;d_{j}}}_{\omega}$
$\displaystyle=-\braket{0}{d_{j}}{k}\braket{k}{b^{\dagger}_{\alpha}b_{\alpha}}{0}$
(138)

are now connected to dipole and photon number transition moments between the
ground and the excited state. Note that to correctly model the coupling to an
external field in the velocity representation, we should also include the
first interaction term in Eq. (5.3.1). In the length form, the electric field
mode creation (annihilation) operators are described by Eq. (37). The response
function equivalent to Eq. (135) is then given by

$\displaystyle\langle\langle\big{(}-ib^{\dagger}_{\alpha}-i\sqrt{\frac{\omega_{\alpha}}{2}}\bm{\lambda}_{\alpha}\cdot\bm{d}\big{)}\big{(}ib_{\alpha}+i\sqrt{\frac{\omega_{\alpha}}{2}}\bm{\lambda}_{\alpha}\cdot\bm{d}\big{)};d_{j}\rangle\rangle_{\omega}.$
(139)

The response function Eq. (135) represents the variation of the number of
displacement field modes, in agreement with the transformation in Eq. (33).

In an analogous way, we can consider the time evolution of the field
displacement coordinates Eqs. (111) and (112). Their time development is
obtained by means of the response functions

$\displaystyle\braket{\braket{q_{\alpha}\;;d_{j}}}_{\omega}=\sum_{k>0}\bigg{(}\frac{\braket{0}{q_{\alpha}}{k}\braket{k}{d_{j}}{0}}{\omega-\omega_{k}}-\frac{\braket{0}{d_{j}}{k}\braket{k}{q_{\alpha}}{0}}{\omega+\omega_{k}}\bigg{)}$
(140)
$\displaystyle\braket{\braket{p_{\alpha}\;;d_{j}}}_{\omega}=\sum_{k>0}\bigg{(}\frac{\braket{0}{p_{\alpha}}{k}\braket{k}{d_{j}}{0}}{\omega-\omega_{k}}-\frac{\braket{0}{d_{j}}{k}\braket{k}{p_{\alpha}}{0}}{\omega+\omega_{k}}\bigg{)},$
(141)

with residues

$\displaystyle\lim_{\omega\to\omega_{k}}(\omega-\omega_{k})\braket{\braket{q_{\alpha}\;;d_{j}}}_{\omega}$
$\displaystyle=\braket{0}{q_{\alpha}}{k}\braket{k}{d_{j}}{0}$ (142)
$\displaystyle\lim_{\omega\to-\omega_{k}}(\omega+\omega_{k})\braket{\braket{q_{\alpha}\;;d_{j}}}_{\omega}$
$\displaystyle=-\braket{0}{d_{j}}{k}\braket{k}{q_{\alpha}}{0}$ (143)
$\displaystyle\lim_{\omega\to\omega_{k}}(\omega-\omega_{k})\braket{\braket{p_{\alpha}\;;d_{j}}}_{\omega}$
$\displaystyle=\braket{0}{p_{\alpha}}{k}\braket{k}{d_{j}}{0}$ (144)
$\displaystyle\lim_{\omega\to-\omega_{k}}(\omega+\omega_{k})\braket{\braket{p_{\alpha}\;;d_{j}}}_{\omega}$
$\displaystyle=-\braket{0}{d_{j}}{k}\braket{k}{p_{\alpha}}{0}$ (145)

that provide information on the transition dipole and photon moments.

####  5.3.2 Photonic environment perturbation

Classical charge currents $\bm{J}_{e}(\bm{r},t)$ are a source of
electromagnetic fields, which means that they can directly perturb the cavity
photon field. The interaction term for the Pauli-Fierz Hamiltonian in Eq. (3)
coupled to a classical external current reads162

$H_{\bm{J}_{e}}=-\frac{1}{c}\int\;d^{3}r\,\bm{J}_{e}(\bm{r},t)\cdot\mathbf{A}(\bm{r}).$
(146)

Following Refs. 187, 94, 56, resolving the external currents into modes which
act directly on the $\alpha$ mode of the cavity, in the dipole approximation
and the length form, the interaction Hamiltonian can be written as

$H_{\bm{J}_{e}}=\sum_{\alpha}\frac{j_{\alpha}(t)}{2\omega_{\alpha}}\;p_{\alpha}=\sum_{\alpha}p_{\alpha}\;\frac{1}{2\omega_{\alpha}}\int
d\omega\;\tilde{j}_{\alpha}(\omega)e^{-i\omega t}$ (147)

where $p_{\alpha}$ is the photon conjugate momentum of mode $\alpha$, and
${j}_{\alpha}$ is connected to the $\alpha$-mode of the time-derivative of the
external current $\bm{J}_{e}$.187, 94, 56 The interaction Hamiltonian Eq.
(147) generates mode excitations into the system. Due to the coupling of light
and matter, this perturbation provides an indirect way to probe the matter and
photon-matter correlation properties.

In the length form, the photonic response function

$\displaystyle\braket{\braket{b^{\dagger}_{\alpha}b_{\alpha}\;;p_{\alpha}}}_{\omega}=$
$\displaystyle\sum_{k>0}\bigg{(}\frac{\braket{0}{b^{\dagger}_{\alpha}b_{\alpha}}{k}\braket{k}{p_{\alpha}}{0}}{\omega-\omega_{k}}-\frac{\braket{0}{p_{\alpha}}{k}\braket{k}{b^{\dagger}_{\alpha}b_{\alpha}}{0}}{\omega+\omega_{k}}\bigg{)}$
(148)

describes the time evolution of the displacement field photon number.
Interestingly, the more general response function

$\displaystyle\braket{\braket{b^{\dagger}_{\beta}b_{\beta}\;;p_{\alpha}}}_{\omega}=$
(149)
$\displaystyle\sum_{k>0}\bigg{(}\frac{\braket{0}{b^{\dagger}_{\beta}b_{\beta}}{k}\braket{k}{p_{\alpha}}{0}}{\omega-\omega_{k}}-\frac{\braket{0}{p_{\alpha}}{k}\braket{k}{b^{\dagger}_{\beta}b_{\beta}}{0}}{\omega+\omega_{k}}\bigg{)},$
(150)

which describes the time evolution of the $\beta$-photon number when the
system is coupled to the current $j_{\alpha}(t)$

$\displaystyle\braket{b^{\dagger}_{\beta}b_{\beta}}^{(1)}(t)=\int
d\omega\;e^{-i\omega t+\eta
t}\braket{\braket{b^{\dagger}_{\beta}b_{\beta},p_{\alpha}}}_{\omega+i\eta}\frac{\tilde{j}_{\alpha}(\omega)}{2\omega_{\alpha}}$
(151)

is different from zero even when $\alpha\neq\beta$. This happens because of
the coupling to matter degrees of freedom, which provide an indirect
interaction between different photon modes. This result opens up the
possibility of photon transfer between modes with different frequencies or
polarization. In the limit of zero coupling, the response function Eq. (151)
would differ from zero only when $\alpha=\beta$, as photons do not directly
interact. Analogously, the response function

$\displaystyle\braket{\braket{p_{\alpha}\;;p_{\beta}}}_{\omega}=\sum_{k>0}\bigg{(}\frac{\braket{0}{p_{\alpha}}{k}\braket{k}{p_{\beta}}{0}}{\omega-\omega_{k}}-\frac{\braket{0}{p_{\beta}}{k}\braket{k}{p_{\alpha}}{0}}{\omega+\omega_{k}}\bigg{)},$
(152)

describes the time evolution of the $\beta$-cavity mode momentum when
perturbing the $\alpha$-cavity mode

$\braket{p_{\alpha}}^{(1)}(t)=\int d\omega\;e^{-i\omega t+\eta
t}\braket{\braket{p_{\alpha},p_{\beta}}}_{\omega+i\eta}\frac{\tilde{j}_{\beta}(\omega)}{2\omega_{\beta}}.$
(153)

As discussed in the other sections, the transformation in Eq. (33) links the
response functions for the physical observables in the velocity and length
form. For the creation and annihilation operators of the electric photon
states, we then need to employ the transformation in Eq. (37).

We can also describe the time evolution of matter degrees of freedom when the
system is perturbed through Eq. (147). The response function

$\displaystyle\braket{\braket{d_{j}\;;p_{\alpha}}}_{\omega}=\sum_{k>0}\bigg{(}\frac{\braket{0}{d_{j}}{k}\braket{k}{p_{\alpha}}{0}}{\omega-\omega_{k}}-\frac{\braket{0}{p_{\alpha}}{k}\braket{k}{d_{j}}{0}}{\omega+\omega_{k}}\bigg{)}$
(154)

is analogous to Eq. (141) as it describes the time evolution of the dipole
moment when the photon degrees of freedom are perturbed by ${j}_{\alpha}(t)$

$\braket{\bm{d}}^{(1)}(t)=\int d\omega\;e^{-i\omega t+\eta
t}\braket{\braket{\bm{d},p_{\alpha}}}_{\omega+i\eta}\frac{\tilde{j}_{\alpha}(\omega)}{2\omega_{\alpha}}.$
(155)

The response functions in Eqs. (135), (140), and (154) are examples of
response functions which describe the cross-talk between photon and matter
degrees of freedom. In fact, in the limit of $\lambda=0$, the eigenstates of
the system would be given by the direct product between bare molecular states
and photons states. As a consequence, Eqs. (135), (140), (141), and (154)
would be identically zero. This means that a perturbation acting on the
molecular Hilbert space cannot induce any time evolution in the (decoupled)
cavity-radiation space, and vice versa.

####  5.3.3 Static perturbations and energy derivatives

This section provides several examples of static perturbations on the
polaritonic system and their connection to the energy derivatives.

##### External electric and magnetic fields.

A static external electric field is described by the scalar potential
$\phi_{e}$ that can be expanded in a Taylor series around the origin

$\phi_{e}(\bm{r})=\phi(0)-r_{i}E_{i}^{(0)}-\frac{1}{2}r_{i}(\partial_{i}E_{j}){(0)}r_{j}+\dots$
(156)

where we employed Einstein’s summation convention. If we retain only the
lowest expansion term, the interaction operator in the Hamiltonian is the
familiar electric dipole term

$H_{int}=-\bm{d}\cdot\mathbf{E}_{e}(0).$ (157)

Note that the transformation in Eq. (33) commutes with this operator, which
implies it is the same in both the length and the velocity representation. The
static electric dipole polarizability $\bm{\alpha}_{0}$ then reads

$-\bm{\alpha}_{0}=\frac{d^{2}{E}}{d^{2}\mathbf{E}_{e}}\bigg{|}_{\mathbf{E}_{e}=0}=2\sum_{k>0}\frac{\bra{0}\bm{d}\ket{k}\bra{k}\bm{d}\ket{0}}{E_{0}-E_{k}}=\braket{\braket{\bm{d};\bm{d}}}_{\omega=0}.$
(158)

An analogous interaction term is introduced when the system is immersed in a
static uniform external magnetic field, which can be described by the
following vector potential206

$\mathbf{A}_{e}=\frac{1}{2}\mathbf{B}_{e}\times\bm{r}.$

In the length form, the interaction Hamiltonian reads

$H_{int}=-\bm{m}\cdot\mathbf{B}_{e}-\frac{1}{2}\mathbf{B}_{e}^{T}\bm{\chi}\mathbf{B}_{e}$
(159)

where

$\displaystyle\bm{m}$
$\displaystyle=\sum_{i}\frac{1}{2c}\bm{l}_{i}=\sum_{i}\frac{1}{2c}\bm{r}_{i}\times\bm{p}_{i}$
(160) $\displaystyle\chi_{\alpha\beta}$
$\displaystyle=\sum_{i}\frac{1}{4c^{2}}(r_{i\alpha}r_{i\beta}-\delta_{\alpha\beta}r_{i}^{2})$
(161)

and since we are employing the length form, these quantities refer only to the
matter subsystem. The static molecular magnetizability $\bm{\xi}_{0}$ is
computed as the second derivative of the energy with respect to the magnetic
field

$\bm{\xi}_{0}=-\frac{d^{2}{E}}{d^{2}\mathbf{B}_{e}}\bigg{|}_{\mathbf{B}_{e}=0}=-\braket{\braket{\bm{m};\bm{m}}}_{\omega=0}-\braket{0}{\bm{\chi}}{0}.$
(162)

The first term in Eq. (162) is called paramagnetic contribution, while the
second term is called diamagnetic contribution and usually dominates for
closed-shell molecules. Magnetizabilities have been extensively studied in the
framework of molecular response theory,219, 220, 221, 222, 223, 224, 225 but
the effect of photon-dressing, to the best of our knowledge, is yet to be
explored. The evaluation of magnetic properties for approximate wave functions
is affected by an origin dependence on the choice of the gauge origin of the
external field. In molecular response theory, the origin independence can be
recovered by using gauge invariant atomic orbitals (GIAO).87, 226

##### Spin interactions.

Nuclear magnetic resonance (NMR) and electronic paramagnetic resonance (EPR)
provide important information on the molecular structure and find wide
applications in chemistry. The electron spin is connected to magnetic moment
of the electron

$\displaystyle\bm{m}^{s}_{i}=-g_{e}\mu_{B}\bm{s}_{i},$ (163)

where $g_{e}$ is the electronic $g$ factor, $\mu_{B}$ is the Bohr magneton and
$\bm{s}_{i}$ is the spin operator of electron $i$. In the same way, the
nuclear spin is connected to the nuclear magnetic moment

$\displaystyle\bm{M}^{s}_{N}=\gamma_{N}\bm{I}_{N}.$ (164)

where $\gamma_{N}$ is the magnetogyric ratio and $\bm{I}_{N}$ is the spin
operator of nucleus $N$. The spin degrees of freedom modify the Pauli-Fierz
Hamiltonian in Eq. (4.2), introducing spin-spin, spin-orbit and spin-external
field interaction terms. The cavity environment introduces a twofold
modification on NMR and EPR spectra: it modifies the ground state wave
function, which is now dressed by the photons, and introduces further terms in
the Hamiltonian due to the interaction of the spins with the cavity field. The
effects of these additional QED terms are still to be explored and require
going beyond the dipole approximation.

We now focus on the NMR properties of closed-shell molecules and disregard the
electronic spin contribution to the magnetic field. The total vector potential
$\mathbf{A}^{tot}(\bm{r}_{i})$ has contributions from the cavity field, the
external field, and the nuclear spins

$\displaystyle\begin{split}\mathbf{A}^{tot}(\bm{r})&=\mathbf{A}(\bm{r})+\frac{{\mathbf{B}_{e}}\times\bm{r}}{2}+\frac{1}{c^{2}}\sum_{N}\left(\frac{\bm{M}_{N}^{s}\times(\bm{r}-\bm{R}_{N})}{{|\bm{r}-\bm{R}_{N}}|^{3}}\right),\end{split}$
(165)

and the total magnetic induction is

$\displaystyle\mathbf{B}^{tot}(\bm{r})=\nabla\times\mathbf{A}^{tot}(\bm{r}).$
(166)

Therefore, the following interaction terms need to be included in the
Hamiltonian in Eq. (4.2):

$\displaystyle\begin{split}H_{int}=&-\sum_{i}\bm{m}^{s}_{i}\cdot\mathbf{B}^{tot}(\bm{r}_{i})-\sum_{N}\bm{M}^{s}_{N}\cdot\mathbf{B}^{tot}(\bm{R}_{N})\\\
&+\sum_{M>N}\frac{1}{2c^{2}}\frac{{R}_{MN}^{2}\left({\bm{M}^{s}_{M}}\cdot{\bm{M}^{s}_{N}}\right)-3\left({\bm{M}^{s}_{M}}\cdot\bm{R}_{MN}\right)\left(\bm{R}_{MN}\cdot{\bm{M}^{s}_{N}}\right)}{{R}_{MN}^{5}},\end{split}$
(167)

where the first term mediates the indirect coupling of nuclear spins, the
second term is the nuclear Zeeman interaction, and the last term in Eq. (167)
is the direct couplings between nuclear dipole magnetic moments. Moreover,
using Eq. (165) in the velocity Hamiltonian leads to additional interaction
terms. The NMR properties are then connected to the following energy
derivatives87

$\displaystyle\mathbf{E}^{(1,1)}_{N}=\frac{d^{2}E(\mathbf{B}_{e},\bm{M}^{s})}{d\mathbf{B}_{e}d\bm{M}^{s}_{N}}\bigg{|}_{\mathbf{B}_{e}=0,\bm{M}^{s}=0}$
(168)
$\displaystyle\mathbf{E}^{(0,2)}_{M,N}=\frac{d^{2}E(\mathbf{B}_{e},\bm{M}^{s})}{d\bm{M}^{s}_{M}d\bm{M}^{s}_{N}}\bigg{|}_{\mathbf{B}_{e}=0,\bm{M}^{s}=0},$
(169)

where $\mathbf{E}^{(1,1)}_{N}$ is related to the nuclear shielding tensor
$\bm{\sigma}_{N}$

$\displaystyle\bm{\sigma}_{N}=\bm{1}+\mathbf{E}^{(1,1)}_{N},$ (170)

which describes the electron shielding effects on nucleus $N$. At the same
time, $\mathbf{E}^{(0,2)}_{M,N}$ is connected to the direct $\bm{D}_{M,N}$ and
indirect (electron-mediated) $\bm{K}_{M,N}$ spin-spin coupling tensors

$\displaystyle\mathbf{E}^{(0,2)}_{M,N}=\bm{K}_{M,N}+\bm{D}_{M,N}.$ (171)

As for magnetizabilities, the computation of NMR properties is affected by
origin dependencies, which are usually eliminated using GIAO orbitals. The
modifications of the ground state density induced by the QED environment
change how the nuclei are shielded from the electrons, affecting the NMR
properties. Moreover, the novel interaction terms between spins and the cavity
magnetic field can possibly affect spin-spin coupling, resulting in
modifications to the shape of NMR signals.

##### QED environment perturbations.

A modification of the QED environment is reflected in the coupling strength
and the frequency of the cavity modes. If we consider a single-mode and
perturb the coupling strength

$\bm{\lambda}\to\bm{\lambda}+\Delta\bm{\lambda},$ (172)

the Hamiltonian reads

$\displaystyle H^{\prime}=$
$\displaystyle\sum_{pq}h_{pq}E_{pq}+\frac{1}{2}\sum_{pqrs}g_{pqrs}e_{pqrs}+h_{nuc}$
$\displaystyle+\sqrt{\frac{\omega_{\alpha}}{2}}(\bm{\lambda}\cdot\bm{d})(b^{\dagger}+b)+\frac{1}{2}(\bm{\lambda}\cdot\bm{d})^{2}$
$\displaystyle+\omega_{\alpha}\left\\{{b}^{\dagger}{b}+\frac{1}{2}\right\\}$
$\displaystyle+\sqrt{\frac{\omega_{\alpha}}{2}}(\Delta\bm{\lambda}\cdot\bm{d})(b^{\dagger}+b)+(\bm{\lambda}\cdot\bm{d})(\Delta\bm{\lambda}\cdot\bm{d})$
$\displaystyle+\frac{1}{2}(\Delta\bm{\lambda}\cdot\bm{d})^{2}$ (173)
$\displaystyle=$ $\displaystyle H^{l}+V^{(1)}+V^{(2)},$ (174)

where $H^{l}$ is the unperturbed dipole Hamiltonian in the length
representation Eq. (4.2), and we have defined the first and second-order
perturbators as

$\displaystyle V^{(1)}$
$\displaystyle=\sqrt{\frac{\omega_{\alpha}}{2}}(\Delta\bm{\lambda}\cdot\bm{d})(b^{\dagger}+b)+(\bm{\lambda}\cdot\bm{d})(\Delta\bm{\lambda}\cdot\bm{d})$
(175) $\displaystyle V^{(2)}$
$\displaystyle=\frac{1}{2}(\Delta\bm{\lambda}\cdot\bm{d})^{2}.$ (176)

We therefore obtain expressions for the first and second energy derivatives

$\displaystyle\frac{dE}{d\Delta\bm{\lambda}}\bigg{|}_{\bm{\Delta\lambda}=0}$
$\displaystyle=\braket{0}{\bm{d}\bigg{(}\sqrt{\frac{\omega_{\alpha}}{2}}(b^{\dagger}+b)+(\bm{\lambda}\cdot\bm{d})\bigg{)}}{0}$
(177)
$\displaystyle\frac{d^{2}E}{d\Delta\bm{\lambda}^{2}}\bigg{|}_{\bm{\Delta\lambda}=0}$
$\displaystyle=\braket{0}{\bm{d}\bm{d}}{0}+\frac{1}{2}\braket{\braket{\frac{dV^{(1)}}{d\Delta\lambda};\frac{dV^{(1)}}{d\Delta\lambda}}}_{\omega=0}.$
(178)

These derivatives measure how sensitive the system is to a change in the
coupling strength. In particular, Eq. (178) can be interpreted as the static
ground state polarizability with respect to the cavity field fluctuations. We
note that because of the dipole self-energy in $H^{l}$, at very large coupling
strength, these expressions can be approximated as

$\displaystyle\frac{dE}{d\Delta\bm{\lambda}}\bigg{|}_{\bm{\Delta\lambda}=0}$
$\displaystyle\approx\bm{\lambda}\cdot\braket{0}{\bm{d}\bm{d}}{0}$ (179)
$\displaystyle\frac{d^{2}E}{d\Delta\bm{\lambda}^{2}}\bigg{|}_{\bm{\Delta\lambda}=0}$
$\displaystyle\approx\frac{1}{2}\braket{\braket{(\bm{\lambda}\cdot\bm{d})\bm{d};(\bm{\lambda}\cdot\bm{d})\bm{d}}}_{\omega=0}$
(180)
$\displaystyle=-\sum_{k>0}\frac{|\bm{\lambda}\cdot\bra{0}\bm{d}\bm{d}\ket{k}|^{2}}{E_{k}-E_{0}}$
(181)

As we increase $\bm{\lambda}$, from Eq. (179) the energy increases since
$\braket{0}{\bm{d}\bm{d}}{0}$ is positive definite, but the Hessian in Eq.
(181) is negative. Therefore, the system is endowed with a stable ground
state. Riso et al. showed that for infinitely large coupling strength, the
polaritonic ground state energy will converge to the SC-QED-HF energy.64 On
the other hand, if the dipole self-energy is disregarded, the remaining terms
are

$\displaystyle\frac{dE}{d\Delta\bm{\lambda}}\bigg{|}_{\bm{\Delta\lambda}=0}$
$\displaystyle\approx\braket{0}{\bm{d}\sqrt{\frac{\omega_{\alpha}}{2}}(b^{\dagger}+b)}{0}$
(182)
$\displaystyle\frac{d^{2}E}{d\Delta\bm{\lambda}^{2}}\bigg{|}_{\bm{\Delta\lambda}=0}$
$\displaystyle\approx\frac{\omega_{\alpha}}{4}\braket{\braket{\bm{d}(b^{\dagger}+b);\bm{d}(b^{\dagger}+b)}}_{\omega=0}$
(183)
$\displaystyle=-\frac{\omega_{\alpha}}{2}\sum_{k>1}\frac{|\bra{0}\bm{d}(b^{\dagger}+b)\ket{k}|^{2}}{E_{k}-E_{0}}.$
(184)

While the Hessian is still negative, there is no guarantee from Eq. (182) that
the energy is increasing, and the ground state energy could decrease
indefinitely. The dipole self-energy has been shown to be fundamental for the
stability of the ground state.186
A modification of the cavity environment can effectively change the photon
frequency

$\omega_{\alpha}\to\omega_{\alpha}+\Delta\omega_{\alpha}.$ (185)

From the expansion

$\sqrt{\omega_{\alpha}+\Delta\omega_{\alpha}}\approx\sqrt{\omega_{\alpha}}+\frac{1}{2}\sqrt{\omega_{\alpha}}\frac{\Delta\omega_{\alpha}}{\omega_{\alpha}}-\frac{1}{8}\sqrt{\omega_{\alpha}}\left(\frac{\Delta\omega_{\alpha}}{\omega_{\alpha}}\right)^{2}+\dots$
(186)

for $\Delta\omega_{\alpha}\ll\omega_{\alpha}$, the perturbed Hamiltonian reads

$\displaystyle H^{\prime}$ $\displaystyle=H^{l}$
$\displaystyle+\frac{1}{2}\frac{\Delta\omega_{\alpha}}{\omega_{\alpha}}\sqrt{\frac{\omega_{\alpha}}{2}}(\bm{\lambda}\cdot\bm{d})(b^{\dagger}+b)+\Delta\omega_{\alpha}\;b^{\dagger}b$
$\displaystyle-\frac{1}{8}\left(\frac{\Delta\omega_{\alpha}}{\omega_{\alpha}}\right)^{2}\sqrt{\frac{\omega_{\alpha}}{2}}(\bm{\lambda}\cdot\bm{d})(b^{\dagger}+b)+\dots$
(187)

The first and second-order perturbation operators are therefore

$\displaystyle V^{(1)}$
$\displaystyle=\frac{1}{2}\frac{\Delta\omega_{\alpha}}{\omega_{\alpha}}\sqrt{\frac{\omega_{\alpha}}{2}}(\bm{\lambda}\cdot\bm{d})(b^{\dagger}+b)+\Delta\omega_{\alpha}\;b^{\dagger}b$
(188) $\displaystyle V^{(2)}$
$\displaystyle=-\frac{1}{8}\left(\frac{\Delta\omega_{\alpha}}{\omega_{\alpha}}\right)^{2}\sqrt{\frac{\omega_{\alpha}}{2}}(\bm{\lambda}\cdot\bm{d})(b^{\dagger}+b).$
(189)

The energy first and second energy derivatives read

$\displaystyle\frac{dE}{d\Delta\omega_{\alpha}}\bigg{|}_{\Delta\omega_{\alpha}=0}$
$\displaystyle=\braket{0}{\frac{1}{2\omega_{\alpha}}\sqrt{\frac{\omega_{\alpha}}{2}}(\bm{\lambda}\cdot\bm{d})(b^{\dagger}+b)+\;b^{\dagger}b}{0}$
(190)
$\displaystyle\frac{d^{2}E}{d\Delta\omega_{\alpha}^{2}}\bigg{|}_{\Delta\omega_{\alpha}=0}$
$\displaystyle=-\frac{1}{4}\left(\frac{1}{\omega_{\alpha}}\right)^{2}\sqrt{\frac{\omega_{\alpha}}{2}}\braket{0}{(\bm{\lambda}\cdot\bm{d})(b^{\dagger}+b)}{0}$
$\displaystyle+\frac{1}{2}\braket{\braket{\frac{dV^{(1)}}{d\Delta\omega_{\alpha}};\frac{dV^{(1)}}{d\Delta\omega_{\alpha}}}}_{\omega=0}.$
(191)

These derivatives measure how sensitive the system is to a change in the
photon frequency.

In the following section, we develop response theory for approximate wave
functions, focusing on QED-HF and QED-CC response theory.

##  6 Ab initio approximations

In this section, we derive the linear and quadratic response equations for the
ab initio QED-HF and QED-CC methods.

###  6.1 QED-HF

The QED-HF is the generalization of the HF model to molecules in a QED
environment.60, 61 The system is described through the dipole Hamiltonian

$\displaystyle H=$
$\displaystyle\sum_{pq}h_{pq}E_{pq}+\frac{1}{2}\sum_{pqrs}g_{pqrs}e_{pqrs}+h_{nuc}$
$\displaystyle+\sum_{\alpha}\sqrt{\frac{\omega_{\alpha}}{2}}(\bm{\lambda}_{\alpha}\cdot\bm{d})_{pq}E_{pq}(b^{\dagger}_{\alpha}+b_{\alpha})$
(192)
$\displaystyle+\frac{1}{2}\sum_{pqrs}\sum_{\alpha}(\bm{\lambda}_{\alpha}\cdot\bm{d})_{pq}(\bm{\lambda}_{\alpha}\cdot\bm{d})_{rs}E_{pq}E_{rs}$
$\displaystyle+\sum_{\alpha}\omega_{\alpha}\left({b}^{\dagger}_{\alpha}{b}_{\alpha}+\frac{1}{2}\right)$
(193)

where the dipole self-energy term can be included in the molecular Hamiltonian
by a proper modification of the one and two-electron integrals

$\displaystyle H_{mol}$
$\displaystyle=\sum_{pq}h_{pq}E_{pq}+\frac{1}{2}\sum_{pqrs}g_{pqrs}e_{pqrs}$
(194) $\displaystyle+$
$\displaystyle\frac{1}{2}\sum_{pqrs}\sum_{\alpha}(\bm{\lambda}_{\alpha}\cdot\bm{d})_{pq}(\bm{\lambda}_{\alpha}\cdot\bm{d})_{rs}E_{pq}E_{rs}$
$\displaystyle=\sum_{pq}\bar{h}_{pq}E_{pq}+\frac{1}{2}\sum_{pqrs}\bar{g}_{pqrs}e_{pqrs},$
(195)

where we have defined

$\displaystyle\bar{g}_{pqrs}$
$\displaystyle=g_{pqrs}+\sum_{\alpha}(\bm{\lambda}_{\alpha}\cdot\bm{d})_{pq}(\bm{\lambda}_{\alpha}\cdot\bm{d})_{rs}$
(196) $\displaystyle\bar{h}_{pq}$
$\displaystyle=h_{pq}+\frac{1}{2}\sum_{r}\sum_{\alpha}(\bm{\lambda}_{\alpha}\cdot\bm{d})_{pr}(\bm{\lambda}_{\alpha}\cdot\bm{d})_{rq}.$
(197)

The QED-HF wave function ansatz is a factorized state

$\ket{\mathrm{R}}=\ket{\mathrm{HF}}\otimes\sum_{\bm{n}}\prod_{\alpha}(b^{\dagger}_{\alpha})^{n_{\alpha}}\ket{0}c_{\bm{n}},$
(198)

where $\ket{HF}$ is a single Slater determinant, $\ket{0}$ is the photon
vacuum, and $c_{\bm{n}}$ are expansion coefficients for the photon number
states. The QED-HF state is obtained by minimization of the mean value of the
Hamiltonian

$E_{\text{QED-HF}}=\braket{\mathrm{R}}{H}{\mathrm{R}}$ (199)

with respect to electronic and photonic parameters. Haugland et al.60, 61
showed that the photonic parameters for the ground-state QED-HF wave function
define a coherent state through the following transformation

$\displaystyle\ket{\mathrm{R}}$ $\displaystyle=\ket{\mathrm{HF}}\otimes
U_{\text{QED-HF}}\ket{0}\equiv U_{\text{QED-HF}}\ket{HF,0}$ (200)
$\displaystyle U_{\text{QED-HF}}$
$\displaystyle=\prod_{\alpha}\text{exp}\bigg{(}-\frac{\bm{\lambda}_{\alpha}\cdot\braket{\bm{d}}_{\text{QED-
HF}}}{\sqrt{2\omega_{\alpha}}}(b^{\dagger}_{\alpha}-b_{\alpha})\bigg{)}.$
(201)

By transforming the dipole Hamiltonian using Eq. (201) we get

$\displaystyle U^{\dagger}_{\text{QED-HF}}H\,U_{\text{QED-HF}}$
$\displaystyle=\sum_{pq}h_{pq}E_{pq}+\frac{1}{2}\sum_{pqrs}g_{pqrs}e_{pqrs}+h_{nuc}$
$\displaystyle+\sum_{\alpha}\sqrt{\frac{\omega_{\alpha}}{2}}(\bm{\lambda}_{\alpha}\cdot(\bm{d}-\braket{\bm{d}}_{\text{QED-
HF}}))(b^{\dagger}_{\alpha}+b_{\alpha})$
$\displaystyle+\frac{1}{2}\sum_{\alpha}(\bm{\lambda}_{\alpha}\cdot(\bm{d}-\braket{\bm{d}}_{\text{QED-
HF}}))^{2}$
$\displaystyle+\sum_{\alpha}\omega_{\alpha}{b}^{\dagger}_{\alpha}{b}_{\alpha}.$
(202)

We note that the transformed Hamiltonian Eq. (6.1) is now origin invariant
also for charged systems. Minimizing the photon-averaged coherent state
transformed dipole Hamiltonian in Eq. (6.1), with respect to the orbital
coefficients, we obtain the QED-Fock matrix

$\displaystyle F_{pq}=F_{pq}^{e}$
$\displaystyle+\frac{1}{2}\sum_{\alpha}\bigg{(}\sum_{a}(\bm{\lambda}_{\alpha}\cdot\bm{d}_{pa})(\bm{\lambda}_{\alpha}\cdot\bm{d}_{aq})-\sum_{i}(\bm{\lambda}_{\alpha}\cdot\bm{d}_{pi})(\bm{\lambda}_{\alpha}\cdot\bm{d}_{iq})\bigg{)},$
(203)

where $F^{e}_{pq}$ is the standard electronic Fock matrix and the indices $a$
and $i$ refer respectively to virtual and occupied orbitals. The optimization
condition that defines the QED-HF molecular orbitals is Brillouin’s theorem,
$F_{ia}=0$, as in standard HF theory. However, the eigenvalues of the Fock
matrix $\bm{F}$, usually interpreted as orbital energies, are now origin
dependent for charged systems. Nevertheless, the total QED-HF energy is origin
invariant since the occupied-virtual blocks $F_{ia}$ are origin independent.

####  6.1.1 Time-dependent QED-HF

Once the optimized QED-HF reference state $\ket{\mathrm{R}}$ has been
obtained, the time development due to an external perturbation $V^{t}$ is
parametrized as

$\text{exp}\big{(}-i\Lambda(t)\big{)}{\ket{\mathrm{R}}}=\text{exp}\big{(}-i\Lambda(t)\big{)}\ket{\mathrm{HF}}\otimes
U_{\text{QED-HF}}\ket{0},$ (204)

where $\Lambda(t)$ is the Hermitian operator

$\displaystyle\Lambda(t)$
$\displaystyle=\frac{1}{\sqrt{2}}\sum_{ai}\big{(}\kappa_{ai}E_{ai}+\kappa^{*}_{ai}E_{ia}\big{)}+\sum_{\alpha}\big{(}\gamma_{\alpha}b^{\dagger}_{\alpha}+\gamma^{*}_{\alpha}b_{\alpha}\big{)}$
$\displaystyle=\bm{\kappa}(t)+\sum_{\alpha}\big{(}\gamma_{\alpha}b^{\dagger}_{\alpha}+\gamma^{*}_{\alpha}b_{\alpha}\big{)}.$
(205)

The $\kappa_{ai}$ parameters in Eq. (6.1.1) are orbital rotation parameters,
$\gamma_{\alpha}$ describe the evolution of the QED-HF coherent state, and we
compactly defined the orbital rotation operator as

$\bm{\kappa}(t)=\frac{1}{\sqrt{2}}\sum_{ai}\big{(}\kappa_{ai}E_{ai}+\kappa^{*}_{ai}E_{ia}\big{)}.$
(206)

Notice that electronic and photonic operators commute, implying that the
exponential in Eq. (204) can be split into the product of an orbital rotation
exponential and a time-dependent coherent state for the field

$\displaystyle\text{exp}\big{(}-i\Lambda(t)\big{)}\ket{R}$
$\displaystyle=\text{exp}\big{(}-i\sum_{\alpha}\big{(}\gamma_{\alpha}b^{\dagger}_{\alpha}+\gamma^{*}_{\alpha}b_{\alpha}\big{)}\big{)}\text{exp}\big{(}-i\bm{\kappa}(t)\big{)}\ket{R}$
(207)
$\displaystyle=\text{exp}\big{(}-i\bm{\kappa}(t)\big{)}\text{exp}\big{(}-i\sum_{\alpha}\big{(}\gamma_{\alpha}b^{\dagger}_{\alpha}+\gamma^{*}_{\alpha}b_{\alpha}\big{)}\big{)}\ket{R}$
(208)

Equivalently, we can parameterize the time dependence as

$\ket{\text{QED-
HF}}(t)=\text{exp}\bigg{(}-\sum_{\alpha}\frac{\bm{\lambda}_{\alpha}\cdot\braket{\bm{d}}_{\text{QED-
HF}}}{\sqrt{2\omega}}(b^{\dagger}_{\alpha}-b_{\alpha})\bigg{)}\text{exp}\big{(}-i\sum_{\alpha}\big{(}\gamma_{\alpha}b^{\dagger}_{\alpha}+\gamma^{*}_{\alpha}b_{\alpha}\big{)}\big{)}\text{exp}\big{(}-i\bm{\kappa}(t)\big{)}\ket{\mathrm{HF}}\otimes\ket{0},$
(209)

where we moved the coherent state transformation to the left. Since Eq. (209)
and Eq. (208) differ by an unimportant phase factor, these parameterizations
lead to the same time evolution. By using the Ehrenfest theorem and developing
the equations in orders of the perturbation, we obtain the zero, first and
second-order equations

$\displaystyle\braket{[\Omega^{(0)},H]}_{\mathrm{R}}=0$ (210) $\displaystyle
i\braket{\mathrm{R}}{[\frac{\partial}{\partial
t}\Lambda^{(1)},\Omega^{(0)}]}{\mathrm{R}}=\braket{\mathrm{R}}{[\Lambda^{(1)},[\Omega^{(0)},H]]}{\mathrm{R}}-i\braket{\mathrm{R}}{[\Omega^{(0)},V^{t}]}{\mathrm{R}}-i\braket{\mathrm{R}}{[\Omega^{(1)},H]}{\mathrm{R}}$
(211) $\displaystyle\langle\mathrm{R}|\big{[}\frac{\partial}{\partial
t}\Lambda^{(2)},\Omega^{(0)}\big{]}+\big{[}\frac{\partial}{\partial
t}\Lambda^{(1)},\Omega^{(1)}\big{]}+{i}\big{[}\Lambda^{(1)},\big{[}\frac{\partial}{\partial
t}\Lambda^{(1)},\Omega^{(0)}\big{]}\big{]}|\mathrm{R}\rangle=\braket{\mathrm{R}}{[V^{t},\Omega^{(1)}]}{\mathrm{R}}+\braket{\mathrm{R}}{[H,\Omega^{(2)}]}{\mathrm{R}}$
$\displaystyle+i\langle\mathrm{R}|\big{[}\Lambda^{(2)},\big{[}H,\Omega^{(0)}\big{]}\big{]}|R\rangle+i\langle\mathrm{R}|\big{[}\Lambda^{(1)},\big{[}V^{t},\Omega^{(0)}\big{]}\big{]}|\mathrm{R}\rangle+i\langle\mathrm{R}|\big{[}\Lambda^{(1)},\big{[}H,\Omega^{(1)}\big{]}\big{]}|\mathrm{R}\rangle-\frac{1}{2}\langle\mathrm{R}|\big{[}\Lambda^{(1)},\big{[}\Lambda^{(1)},[H,\Omega^{(0)}]\big{]}\big{]}|\mathrm{R}\rangle,$
(212)

where we also accounted for the possibility that $\Omega$ depends on the
external perturbation $\Omega\equiv\Omega(V^{t})$. The apexes $(0)$-$(2)$
refer to the expansion order and

$\displaystyle\Lambda^{(1)}(t)$
$\displaystyle=\frac{1}{\sqrt{2}}\sum_{ai}\big{(}\kappa_{ai}^{(1)}E_{ai}+\kappa^{*}_{ai}\,{}^{(1)}\;E_{ia}\big{)}+\sum_{\alpha}\big{(}\gamma_{\alpha}^{(1)}b^{\dagger}+\gamma_{\alpha}^{*}\,{}^{(1)}\;b\big{)}$
(213) $\displaystyle\Lambda^{(2)}(t)$
$\displaystyle=\frac{1}{\sqrt{2}}\sum_{ai}\big{(}\kappa_{ai}^{(2)}E_{ai}+\kappa^{*}_{ai}\,{}^{(2)}\;E_{ia}\big{)}+\sum_{\alpha}\big{(}\gamma_{\alpha}^{(2)}b^{\dagger}+\gamma_{\alpha}^{*}\,{}^{(2)}\;b\big{)}.$
(214)

As long as $\Omega$ is a one-electron operator or a purely photonic 1/2
operator, the zero-order condition is satisfied due to the QED-HF optimization
condition. The first-order equation is non-trivial, and its solution gives us
the time evolution of the parameters to first-order in the perturbation.
Following the derivation of the response equations in molecular response
theory,84, 227 we make use of the operators

$\Omega\in\bigg{\\{}b_{\alpha},\,\frac{1}{\sqrt{2}}\,e^{-i\bm{\kappa}}E_{ia}e^{i\bm{\kappa}},\,b^{\dagger}_{\alpha},\,\frac{1}{\sqrt{2}}\,e^{-i\bm{\kappa}}E_{ai}e^{i\bm{\kappa}}\bigg{\\}}.$
(215)

Making this choice, we ensure the response equations are identical to the ones
derived from the time-dependent variational principle.84, 228, 229 It is
convenient to move from the time domain to the frequency domain by performing
a Fourier transformation of the parameters

$\displaystyle\kappa_{ai}(t)$ $\displaystyle=\int d\omega\;e^{-i\omega
t}\kappa_{ai}^{\omega}+\int\int d\omega\,d\omega^{\prime}\;e^{-i\omega
t}e^{-i\omega^{\prime}t}\kappa_{ai}^{\omega,\omega^{\prime}}+\dots$ (216)
$\displaystyle\gamma(t)$ $\displaystyle=\int d\omega\;e^{-i\omega
t}\gamma_{\alpha}^{\omega}+\int\int d\omega\,d\omega^{\prime}\;e^{-i\omega
t}e^{-i\omega^{\prime}t}\gamma_{\alpha}^{\omega,\omega^{\prime}}+\dots$ (217)

such that the operators

$\displaystyle\Lambda^{\omega}$
$\displaystyle=\frac{1}{\sqrt{2}}\sum_{ai}\big{(}\kappa_{ai}^{\omega}E_{ai}+[\kappa_{ai}^{-\omega}]^{*}\;E_{ia}\big{)}$
$\displaystyle+\sum_{\alpha}\big{(}\gamma^{\omega}_{\alpha}b^{\dagger}_{\alpha}+[\gamma^{-\omega}_{\alpha}]^{*}\;b_{\alpha}\big{)}$
(218) $\displaystyle\Lambda^{\omega,\omega^{\prime}}$
$\displaystyle=\frac{1}{\sqrt{2}}\sum_{ai}\big{(}\kappa_{ai}^{\omega,\omega^{\prime}}E_{ai}+[\kappa_{ai}^{-\omega,-\omega^{\prime}}]^{*}\;E_{ia}\big{)}$
$\displaystyle+\sum_{\alpha}\big{(}\gamma^{\omega,\omega^{\prime}}_{\alpha}b^{\dagger}_{\alpha}+[\gamma^{-\omega,-\omega^{\prime}}_{\alpha}]^{*}\;b_{\alpha}\big{)}$
(219)

are the Fourier components of Eqs. (213) and (214)

$\displaystyle\Lambda^{(1)}(t)=\int d\omega\;e^{-i\omega t}\Lambda^{\omega}$
(220) $\displaystyle\Lambda^{(2)}(t)=\int
d\omega\,d\omega^{\prime}\;e^{-i\omega
t}e^{-i\omega^{\prime}t}\Lambda^{\omega,\omega^{\prime}}.$ (221)

##### Linear response equations.

From Eq. (211), performing a Fourier transform and by using Eq. (218) and Eq.
(215), after some algebra, we obtain the standard matrix equation

$\left[\begin{pmatrix}\mathbf{A}&\mathbf{B}\\\
\mathbf{B}^{*}&\mathbf{A}^{*}\end{pmatrix}-\omega\begin{pmatrix}\bm{1}&0\\\
0&-\bm{1}\end{pmatrix}\right]\begin{pmatrix}\bm{X}\\\
\bm{Y}\end{pmatrix}=i\begin{pmatrix}\bm{g}_{1}\\\
\bm{g}_{2}\end{pmatrix}\equiv i\bm{g},$ (222)

where the vectors $\bm{X}$ and $\bm{Y}$ collect the Fourier transformed
parameters

$\bm{X}=\begin{pmatrix}\gamma^{\omega}_{\alpha}\\\
\kappa_{ai}^{\omega}\end{pmatrix}\quad\bm{Y}=\begin{pmatrix}[\gamma^{-\omega}_{\alpha}]^{*}\\\
[\kappa_{ai}^{-\omega}]^{*}\end{pmatrix}.$ (223)

The right hand side of Eq. (222) is the generalized gradient

$\bm{g}_{1}=\begin{pmatrix}\braket{[b_{\alpha},V^{\omega}]}_{\mathrm{R}}\\\
\frac{1}{\sqrt{2}}\braket{[E_{ia},V^{\omega}]}_{\mathrm{R}}\end{pmatrix}\quad\bm{g}_{2}=\begin{pmatrix}\braket{[b^{\dagger}_{\alpha},V^{\omega}]}_{\mathrm{R}}\\\
\frac{1}{\sqrt{2}}\braket{[E_{ai},V^{\omega}]}_{\mathrm{R}}\end{pmatrix}.$
(224)

The matrix

$\mathbf{E}^{[2]}=\begin{pmatrix}\bm{A}&\mathbf{B}\\\
\mathbf{B}^{*}&\bm{A}^{*}\end{pmatrix}$ (225)

is the generalized Hessian matrix while

$\bm{S}^{[2]}=\begin{pmatrix}\bm{1}&0\\\ 0&-\bm{1}\end{pmatrix}$ (226)

is the generalized metric matrix. Therefore, Eq. (222) is a generalization of
the Casida equations of TDHF in molecular response theory.230, 231 Indeed, the
explicit expressions of $\mathbf{A}$ and $\mathbf{B}$ are

$\displaystyle\mathbf{A}$
$\displaystyle=\begin{pmatrix}\omega_{\alpha}\delta_{\alpha\beta}&\sqrt{{\omega_{\alpha}}}(\bm{\lambda}_{\alpha}\cdot\bm{d}_{ib})\\\
\sqrt{{\omega_{\alpha}}}(\bm{\lambda}_{\alpha}\cdot\bm{d}_{bi})&\mathbf{A}_{el}\end{pmatrix}=\begin{pmatrix}\braket{[b_{\alpha},[H,b^{\dagger}_{\alpha}]]}_{\mathrm{R}}&\frac{1}{\sqrt{2}}\braket{[b_{\alpha},[H,E_{ai}]]}_{\mathrm{R}}\\\
\frac{1}{\sqrt{2}}\braket{[b^{\dagger}_{\alpha},[H,E_{ia}]]}_{\mathrm{R}}&\frac{1}{2}\braket{[E_{jb},[H,E_{ai}]]}_{\mathrm{R}}\end{pmatrix}$
(227) $\displaystyle\mathbf{B}$
$\displaystyle=\begin{pmatrix}0&-\sqrt{{\omega_{\alpha}}}(\bm{\lambda}_{\alpha}\cdot\bm{d}_{ib})\\\
-\sqrt{{\omega_{\alpha}}}(\bm{\lambda}_{\alpha}\cdot\bm{d}_{bi})&\mathbf{B}_{el}\end{pmatrix}=\begin{pmatrix}\braket{[b_{\alpha},[H,b_{\alpha}]]}_{\mathrm{R}}&\frac{1}{\sqrt{2}}\braket{[b_{\alpha},[H,E_{ia}]]}_{\mathrm{R}}\\\
\frac{1}{\sqrt{2}}\braket{[b^{\dagger}_{\alpha},[H,E_{ai}]]}_{\mathrm{R}}&\frac{1}{2}\braket{[E_{bj},[H,E_{ai}]]}_{\mathrm{R}}\end{pmatrix}.$
(228)

The electronic blocks

$\displaystyle\big{(}A_{el}\big{)}_{ia,bj}$
$\displaystyle=\frac{1}{2}\braket{\mathrm{R}}{\left[E_{ia},\left[H,E_{bj}\right]\right]}{\mathrm{R}}=\delta_{ij}F_{ab}-\delta_{ab}F_{ij}+2\bar{g}_{aibj}-\bar{g}_{abji}$
(229) $\displaystyle\big{(}B_{el}\big{)}_{ia,jb}$
$\displaystyle=\frac{1}{2}\braket{\mathrm{R}}{\left[E_{ia},\left[H,E_{jb}\right]\right]}{\mathrm{R}}=\bar{g}_{biaj}-2\bar{g}_{aibj},$
(230)

have the same definition as the standard Casida matrices, $\bm{A}$ and
$\mathbf{B}$, in the TDHF theory. Nevertheless, we point out that these blocks
differ from the bare TDHF matrices for two reasons. First, the QED-HF orbitals
and orbital energies differ from the bare molecular ones. Second, the two-
electron integrals now contain the dipole self-energy contributions, as shown
in Eq. (196). This contribution can be explicitly separated95 following Eqs.
(196) and (197)

$\displaystyle\mathbf{A}_{el}$
$\displaystyle=\mathbf{A}_{el}^{\text{HF}}+\bm{\Delta}$ (231)
$\displaystyle\mathbf{B}_{el}$
$\displaystyle=\mathbf{B}_{el}^{\text{HF}}+\bm{\Delta^{\prime}}$ (232)

where the $\bm{\Delta}$ and $\bm{\Delta^{\prime}}$ matrices now include the
dipole self-energy contribution to the Hamiltonian, while
$\mathbf{A}_{el}^{\text{HF}}$ and $\mathbf{B}_{el}^{\text{HF}}$ only include
the standard electronic integrals. The TD-QED-HF matrices have additional
dimensions due to the cavity modes. The coupling between the molecular and the
photonic parameters is given by the projection of the transition dipole
moments onto the coupling strength vector $\bm{\lambda}$ in the off-diagonal
terms of Eqs. (229) and (230). Moreover, the $\mathbf{A}$ matrix has a non-
zero photonic block that contains the frequencies $\omega_{\alpha}$ of the
cavity modes. As expected, in Eq. (227) there is no direct interaction between
the photon modes, which are coupled indirectly through matter degrees of
freedom. We see that these equations are an extension of the familiar TDHF
response equations since in the zero coupling limit we recover TDHF solutions
and the one photon lines of the empty cavity. We can then write the response
function from the time evolution of an operator $\Omega$

$\displaystyle\braket{\Omega}(t)-\braket{\Omega}_{\mathrm{R}}=\braket{\Omega}^{(1)}(t)+\dots=$
$\displaystyle-i\int_{-\infty}^{\infty}e^{-i\omega
t}\left(\sum_{ai}\frac{{\kappa}^{\omega}_{ai}}{\sqrt{2}}\braket{[\Omega,E_{ai}]}_{\mathrm{R}}+\sum_{ai}\frac{[{\kappa}^{-\omega}_{ai}]^{*}}{\sqrt{2}}\braket{[\Omega,E_{ia}]}_{\mathrm{R}}+{\gamma}^{\omega}\braket{[\Omega,b_{\alpha}^{\dagger}]}_{\mathrm{R}}+[\gamma_{\alpha}^{-\omega}]^{*}\braket{[\Omega,b]}_{\mathrm{R}}\right)+\dots$
(233)

Inverting equation Eq. (222), we identify the linear response function

$\displaystyle\braket{\braket{\Omega,V^{\omega}}}_{\omega}=$
$\displaystyle-i\bigg{(}\sum_{ai}\frac{{\kappa}^{\omega}_{ai}}{\sqrt{2}}\braket{[\Omega,E_{ai}]}_{\mathrm{R}}+\sum_{ai}\frac{[{\kappa}^{-\omega}_{ai}]^{*}}{\sqrt{2}}\braket{[\Omega,E_{ia}]}_{\mathrm{R}}+{\gamma}_{\alpha}^{\omega}\braket{[\Omega,b^{\dagger}]}_{\mathrm{R}}+[\gamma^{-\omega}_{\alpha}]^{*}\braket{[\Omega,b]}_{\mathrm{R}}\bigg{)}$
$\displaystyle=$ $\displaystyle
i\bm{g}_{\Omega}^{\dagger}\begin{pmatrix}\bm{X}\\\
\bm{Y}\end{pmatrix}=-\bm{g}_{\Omega}^{\dagger}\left[\begin{pmatrix}\mathbf{A}&\mathbf{B}\\\
\mathbf{B}^{*}&\mathbf{A}^{*}\end{pmatrix}-\omega\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\right]^{-1}\bm{g},$ (234)

where $\bm{g}_{\Omega}$ has the same structure as $\bm{g}$ with the operator
$\Omega$ replacing the perturbation $V^{\omega}$ in Eq. (224). We note that
the response function fulfils the symmetry relation in Eq. (5). From the
generalized eigenvalue equation

$\left[\begin{pmatrix}\mathbf{A}&\mathbf{B}\\\
\mathbf{B}^{*}&\mathbf{A}^{*}\end{pmatrix}-\omega\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\right]\begin{pmatrix}\bm{X}\\\
\bm{Y}\end{pmatrix}=\begin{pmatrix}0\\\ 0\end{pmatrix},$ (235)

we obtain the eigenvectors $(\bm{x}_{i},\bm{y}_{i})^{T}$ and the spectral
decomposition

$\displaystyle\left[\begin{pmatrix}\mathbf{A}&\mathbf{B}\\\
\mathbf{B}^{*}&\mathbf{A}^{*}\end{pmatrix}-\omega\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\right]^{-1}=\sum_{i>0}\left(\frac{1}{\omega_{j}-\omega}\begin{pmatrix}\bm{x}_{i}\\\
\bm{y}_{i}\end{pmatrix}\otimes\begin{pmatrix}\bm{x}_{i}^{\dagger}&\bm{y}_{i}^{\dagger}\end{pmatrix}-\frac{1}{\omega_{j}+\omega}\begin{pmatrix}\bm{y}^{*}_{i}\\\
\bm{x}_{i}^{*}\end{pmatrix}\otimes\begin{pmatrix}\big{(}\bm{y}_{i}^{*}\big{)}^{\dagger}&\big{(}\bm{x}_{i}^{*}\big{)}^{\dagger}\end{pmatrix}\right).$

We identify the excitation energies of the system as the eigenvalues of Eq.
(235). The transition moments are obtained from the residues of the response
function, identified as

$\braket{0}{\Omega}{k}=\begin{pmatrix}\bm{x}_{k}^{\dagger}&\bm{y}_{k}^{\dagger}\end{pmatrix}\cdot\begin{pmatrix}\braket{[b_{\alpha}^{\dagger},\Omega]}_{\mathrm{R}}\\\
\frac{1}{\sqrt{2}}\braket{[E_{ai},\Omega]}_{\mathrm{R}}\\\
\braket{[b_{\alpha},\Omega]}_{\mathrm{R}}\\\
\frac{1}{\sqrt{2}}\braket{[E_{ia},\Omega]}_{\mathrm{R}}\end{pmatrix}.$ (236)

Although the eigenvalue problem in Eq. (235) is non Hermitian, it is possible
to show that if the computed reference state is close to the ground state, Eq.
(235) has solutions with real and positive eigenvalues.84 We note that in Eq.
(222), setting $\omega=0$, we obtain the static coupled-perturbed QED-HF
equations.

Several approximations can be proposed for these response equations. A
hierarchy of approximations has been discussed by Yang et al.95 for their
TDDFT-Pauli-Fierz (TDDFT-PF) model, which defined a eigenvalue problem
analogous to Eq. (222). They define the Tamm-Dancoff approximation (TDA) by
neglecting $\mathbf{B}_{el}$. However, contrary to standard electronic TDA,
this is not equivalent to a CI-singles approach. For this reason, we suggest
it is more natural to define the TDA approximation by neglecting the whole
$\mathbf{B}$ matrix, since this is equivalent to a CI problem with singly
excited determinants with zero photons and the ground state determinant with
one photon. Neglecting $\mathbf{B}$ and the dipole self-energy contribution in
$\bm{\Delta}$, they define the TDDFT-Jaynes-Cummings approximation, similar to
a JC calculation with all the single excited determinants.

##### Quadratic response equations

The quadratic response equations are obtained in a similar fashion as for the
linear response.84, 232, 85 By employing Eq. (212) and the set of operators
(215), we obtain the second-order response equation

$\displaystyle(\omega_{1}+\omega_{2})\braket{\mathrm{R}}{[\Lambda^{\omega_{1},\omega_{2}},\bm{\xi}]}{\mathrm{R}}-\braket{\mathrm{R}}{[\Lambda^{\omega_{1},\omega_{2}},[\bm{\xi},H]]}{\mathrm{R}}=\frac{1}{2}\hat{P}(\omega_{1},\omega_{2})\big{(}\frac{i}{2}\braket{\mathrm{R}}{[\Lambda^{\omega_{1}},[\Lambda^{\omega_{2}},[\bm{\xi},H]]]}{\mathrm{R}}+\braket{\mathrm{R}}{[\Lambda^{\omega_{1}},[\bm{\xi},V^{\omega_{2}}]]}{\mathrm{R}}\big{)},$
(237)

where $\hat{P}(\omega_{1},\omega_{2})$ sums the permutations of $\omega_{1}$
and $\omega_{2}$ and we used the compact notation

$\bm{\xi}=\begin{pmatrix}b_{\alpha}\\\ \frac{1}{\sqrt{2}}E_{ia}\\\
b^{\dagger}_{\alpha}\\\ \frac{1}{\sqrt{2}}E_{ai}\end{pmatrix}.$ (238)

To compute the quadratic response function
$\braket{\braket{A;B,C}}_{\omega_{1},\omega_{2}}$, together with the
$\mathbf{E}^{[2]}$ and $\bm{S}^{[2]}$ matrices defined in Eqs. (225) and (226)
for the linear response, we need the additional matrices

$\displaystyle E^{[3]}_{ijk}$
$\displaystyle=\frac{1}{2}\braket{\mathrm{R}}{[\xi_{i}^{\dagger},[\xi_{j}^{\dagger},[H,\xi_{k}]]]}{\mathrm{R}}$
(239) $\displaystyle X^{[2]}_{ij}$
$\displaystyle=\braket{\mathrm{R}}{[\xi_{i}^{\dagger},[X,{\xi}_{j}]]}{\mathrm{R}},$
(240)

where $X=A,B$ or $C$. The second-order variation of a time-independent
observable $A$ is

$\displaystyle\braket{A}^{(2)}(t)$
$\displaystyle=i\braket{[\Lambda^{(2)},A]}_{\mathrm{R}}-\frac{1}{2}\braket{[\Lambda^{(1)},[\Lambda^{(1)},A]]}_{\mathrm{R}}$
$\displaystyle=\int\int
d\omega_{1}d\omega_{2}\,\frac{1}{2}\hat{P}(\omega_{1},\omega_{2})\big{(}i\braket{[\Lambda^{\omega_{1},\omega_{2}},A]}_{\mathrm{R}}-\frac{1}{2}\braket{[\Lambda^{\omega_{2}},[\Lambda^{\omega_{1}},A]]}_{\mathrm{R}}\big{)}e^{-i(\omega_{1}+\omega_{2})t},$
(241)

where we made the expression explicitly symmetric in the frequencies. From Eq.
(241), we can identify the quadratic response function, which has poles where
the frequencies or their sum match an excitation energy233

$\displaystyle\braket{\braket{A;B,C}}_{\omega_{1},\omega_{2}}=\hat{P}(\omega_{1},\omega_{2})\big{(}i\braket{[\Lambda^{\omega_{1},\omega_{2}},A]}_{\mathrm{R}}-\frac{1}{2}\braket{[\Lambda^{\omega_{2}},[\Lambda^{\omega_{1}},A]]}_{\mathrm{R}}\big{)}.$
(242)

For the residues and the response function, we then need the following
vectors:

$\displaystyle\big{[}\bm{N}^{a}(\omega_{1}+\omega_{2})\big{]}^{\dagger}$
$\displaystyle=\bm{g}_{A}^{\dagger}[\mathbf{E}^{[2]}-(\omega_{1}+\omega_{2})\bm{S}^{[2]}]^{-1}$
$\displaystyle\bm{N}^{b}(\omega_{1})$
$\displaystyle=[\mathbf{E}^{[2]}-\omega_{1}\bm{S}^{[2]}]^{-1}\bm{g}_{B}$
$\displaystyle\bm{N}^{c}(\omega_{2})$
$\displaystyle=[\mathbf{E}^{[2]}-\omega_{2}\bm{S}^{[2]}]^{-1}\bm{g}_{C}.$

As for the linear response equations, these vectors and matrices have
additional dimensions due to the photonic parameters.