# Real-Space, Real-Time Approach to Quantum-Electrodynamical Time-Dependent
Density Functional Theory

Justin Malave  Department of Physics and Astronomy, Vanderbilt University,
Nashville, Tennessee, 37235, USA  Alexander Ahrens  Department of Physics and
Astronomy, Vanderbilt University, Nashville, Tennessee, 37235, USA  Daniel
Pitagora  Department of Physics and Astronomy, Vanderbilt University,
Nashville, Tennessee, 37235, USA  Cody Covington  Department of Chemistry,
Austin Peay State University, Clarksville, USA  Kálmán Varga
<EMAIL_ADDRESS>Department of Physics and Astronomy, Vanderbilt
University, Nashville, Tennessee, 37235, USA

###### Abstract

The Quantum-Electrodynamical Time-Dependent Density Functional Theory (QED-
TDDFT) equations are solved by time propagating the wave function on a tensor
product of a Fock-space and real-space grid. Applications for molecules in
cavities show the accuracy of the approach. Examples include the coupling
strength and light frequency dependence of the energies, wave functions,
optical absorption spectra, and Rabi splitting magnitudes in cavities, as well
as a description of high harmonic generation in cavities.

##  I Introduction

The possibility of altering physical and chemical properties by coupling
matter to light has attracted intense experimental Hutchison et al. (2012);
Balili et al. (2007); Schachenmayer et al. (2015); Xiang et al. (2020);
Reserbat-Plantey et al. (2021); Coles et al. (2014); Kasprzak et al. (2006);
Schwartz et al. (2011); Plumhof et al. (2014); Hutchison et al. (2013); Wang
et al. (2021a); Basov et al. (2021) and theoretical interestFeist and Garcia-
Vidal (2015); Vendrell (2018); Schäfer and Johansson (2022); Riso et al.
(2022); Badankó et al. (2022); Mazza and Georges (2019); Di Stefano et al.
(2019); Galego et al. (2015); Herrera and Spano (2016); Galego et al. (2016);
Shalabney et al. (2015); Buchholz et al. (2019); Schäfer et al. (2019);
Ruggenthaler et al. (2018); Flick et al. (2015, 2017); Lacombe et al. (2019);
Mandal et al. (2020a); Flick et al. (2019); Latini et al. (2019); Flick et al.
(2018a); Flick and Narang (2018); Garcia-Vidal et al. (2021); Thomas et al.
(2019); Mordovina et al. (2020); Wang et al. (2021b); DePrince (2021);
Haugland et al. (2021); Hoffmann et al. (2020); Flick and Narang (2020);
Mandal et al. (2020b); Galego et al. (2017); Hoffmann et al. (2020); Flick et
al. (2015, 2018b); Schafer et al. (2018); Sidler et al. (2020); Theophilou et
al. (2020); Sidler et al. (2021); Buchholz et al. (2020); Tokatly (2018);
Rivera et al. (2019); Szidarovszky et al. (2020); Cederbaum (2021); Cho et al.
(2022); Pavosevic and Rubio ; Pavošević and Flick (2021); Pavošević et al.
(2022); Cederbaum and Kuleff (2021). There are several excellent review
articles highlighting the present state of the art of experimental and
theoretical approaches related to light-matter interaction in cavities. These
include reviews about the properties of hybrid light-matter states Ebbesen
(2016); Genet et al. (2021), ab initio calculations Ruggenthaler et al.
(2018); Sidler et al. (2022) and molecular polaritonics Li et al. (2022);
Sánchez-Barquilla et al. (2022); Fregoni et al. (2022).

The theoretical and computational description of the coupled light matter
system is challenging because the already difficult quantum many-body problem
of electrons and nuclei is aggravated with the addition of the photon degrees
of freedom. A plethora of approaches going beyond the simple two-level atom
model Jaynes and Cummings (1962) has been proposed in the last few years. Most
of these approaches are based on successful many-body quantum methods adapted
to include photon interactions . The use of the Pauli-Fierz (PF) non-
relativistic quantum electrodynamics Hamiltonian is found to be the most
useful framework Ruggenthaler et al. (2018); Rokaj et al. (2018); Mandal et
al. (2020a, b); Tokatly (2018) for practical calculations. The PF Hamiltonian
is a sum of electronic and photonic Hamiltonians and a cross term describing
the electron-photon interaction. Due to this cross term, one has to use a
coupled electron-photon wave function,

$\sum_{\vec{n}}\Phi_{\vec{n}}(\mathbf{x})\chi_{\vec{n}}$ (1)

where
$\mathbf{x}=(\mathbf{r}_{1},\mathbf{r}_{2},{\ldots},\mathbf{R}_{1},\mathbf{R}_{2},{\ldots})$
are the spatial coordinates of the electrons and nuclei, and
$\vec{n}=(n_{1},n_{2},{\ldots}N_{p})$ are the quantum numbers of the photon
modes. The occupation number basis,
$\chi_{\vec{n}}=|n_{1},n_{2},{\ldots},N_{p}\rangle$, is used to represent the
bosonic Fock-space of photon modes (see Appendix A for a more detailed
definition).

Wave function based approaches Ahrens et al. (2021); Beutel et al. (2021);
Liebenthal et al. (2022); DePrince (2021); Haugland et al. (2020) typically
use coupled electron-photon wave functions and the product form significantly
increases the dimensionality. In Refs. Liebenthal et al. (2022); DePrince
(2021); Haugland et al. (2020); Mordovina et al. (2020) the coupled cluster
(CC) approach is used by defining a reference wave function as a direct
product of a Slater determinant of Hartree-Fock states and a zero-photon
photon number state. The exponentiated cluster operator acting upon this
product state defines the ground state QED-CC wave function. Systematic
improvability is the greatest benefit of this approach. The stochastic
variational method Ahrens et al. (2021); Beutel et al. (2021) (QED-SVM) also
uses a product of matter and photonic wave functions but in this case, the
matter part is described by explicitly correlated Gaussian basis states Mitroy
et al. (2013). The parameters of the variational ansatz are optimized by a
stochastic selection process leading to highly accurate energies and wave
functions. Both the QED-CC and QED-SVM approaches are limited to small atoms
and molecules.

A pioneering approach to describe the interaction of light and matter in
cavities is the quantum electrodynamical density functional theory (QEDFT)
Haugland et al. (2021); Flick et al. (2018b); Schafer et al. (2018);
Ruggenthaler et al. (2014); Tokatly (2013); Flick and Narang (2020); Yang et
al. (2021). The QEDFT is an exact reformulation of the PF Hamiltonian based
many-body wave theory. In practical applications of QEDFT one has to develop
good approximations of the fields and currents so that the auxiliary non-
interacting system generates the same physical quantities as the interacting
system. To facilitate this need, the development of polaritonic exchange-
correlation functionals is a focus of intense research interest Flick et al.
(2018b); Pellegrini et al. (2015); Schafer et al. (2021); Flick (2021). In
QEDFT the spatial and photon wave functions are factorized allowing the
separation of the electronic and photonic components. One then has to solve
coupled equations for the matter and photon parts. The dimensionality limits
the product basis approaches to a few photon modes, while the QEDFT can be
applied to hundreds of thousands of photon modes Flick et al. (2019). More
details and applications of the QEDFT method and its time-dependent version
can be found in Refs. Schäfer et al. (2019); Ruggenthaler et al. (2018); Flick
et al. (2017); Sidler et al. (2021); Flick et al. (2019).

In this paper, we present a QED-DFT (and QED-TDDFT) approach using a coupled
electron-photon wave function like that shown in Eq. (1). The wave function
will be defined on a tensor product (TP) of a spatial grid and Fock state
representation. To distinguish this approach from the aforementioned QEDFT
approach, we will call the present approach QED-DFT-TP and QED-TDDFT-TP. The
QED-DFT-TP is a particular realization of the QEDFT theory using a different
ansatz – a coupled electron-photon wave function.

The tensor product form increases the dimensionality; however, it preserves
the quantized photon states. In this way, we have direct access to non-
classical observables of the photon field e.g., the photon-number, the purity,
or the Mandel Q number Mandel (1979). The coupled electron-photon wave
function provides a more complex description of the light-matter interaction
by calculating the spatial wave function in each photon sector.

Each orbital of the molecule is coupled with different Fock basis states
(photon states) describing the quantized light. The light-matter coupling part
of the Hamiltonian will describe the interaction between the orbital
components in different photons states. The Fock states are orthogonal,
preserving and even increasing the sparsity of the real-space based DFT
Hamiltonian. This sparsity allows the use of the efficient iterative
diagonalization approaches traditionally used in real-space DFT approaches.

As the coupling part and the product ansatz are simple, the present approach
can be easily implemented in any real-space approaches. Plane-wave or orbital-
based DFT methods can see similarly simple adaptation.

The approach will be tested using time-dependent and time-independent problems
and our results will be compared to QED-SVM, QED-CC, and QEDFT calculations.
One expects that the approach can be used everywhere where regular DFT and
TDDFT have been useful, from bond length and density distribution
calculations, to descriptions of high harmonic generation and optical
absorption.

##  II Formalism

###  II.1 Hamiltonian

We assume that the system is nonrelativistic and the coupling to the light can
be described by the dipole approximation. The Pauli-Fierz non-relativistic QED
Hamiltonian provides a consistent quantum description at this level
Ruggenthaler et al. (2018); Rokaj et al. (2018); Mandal et al. (2020a, b);
Tokatly (2018). The dipole approximation assumes that the spatial variation of
the electric field is negligible across the size of the system – physically
valid if the system size is much smaller than the wavelength of the light. The
PF Hamiltonian in the Coulomb gauge is $H=H_{e}+H_{ep}$ where $H_{e}$ is the
electronic Hamiltonian and

$\displaystyle H_{ep}$ $\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{N_{p}}{1\over
2}\left[p_{\alpha}^{2}+\omega_{\alpha}^{2}\left(q_{\alpha}-{\boldsymbol{\lambda}_{\alpha}\over\omega_{\alpha}}\mathbf{D}\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\sum_{\alpha=1}^{N_{p}}\left[\omega_{\alpha}\left(\hat{a}_{\alpha}^{+}\hat{a}_{\alpha}+\frac{1}{2}\right)-\omega_{\alpha}q_{\alpha}\boldsymbol{\lambda}_{\alpha}\mathbf{D}+{\frac{1}{2}}\left(\mathbf{\lambda_{\alpha}}\mathbf{D}\right)^{2}\right],$

(atomic units are used in this work). In Eq. (II.1) $\mathbf{D}$ is the dipole
operator, the photon fields are described by quantized oscillators, and
$q_{\alpha}={\frac{1}{\sqrt{2\omega_{\alpha}}}}(\hat{a}^{+}_{\alpha}+\hat{a}_{\alpha})$
is the displacement field. Here $\hat{a}$ and $\hat{a}^{+}$ are, respectively,
the lowering and raising operators of the quantized harmonic oscillator. This
Hamiltonian describes $N_{p}$ photon modes with photon frequency
$\omega_{\alpha}$ and coupling $\boldsymbol{\lambda}_{\alpha}$. The coupling
term is written as Ruggenthaler et al. (2014)
$\boldsymbol{\lambda}_{\alpha}={1/\sqrt{\epsilon_{0}}}\mathbf{u}_{\alpha}(\mathbf{r}_{0})$,
where $\mathbf{u}_{\alpha}(\mathbf{r})$ is the cavity mode function and
$\mathbf{r}_{0}$ is the center of the cavity where the molecule is placed. The
first term in Eq. (II.1) is the Hamiltonian of the photon modes, the second
term couples the photons to the dipole, and the last term is the dipole self-
interaction,

$H_{d}={1\over
2}\sum_{\alpha=1}^{N_{p}}\left(\boldsymbol{\lambda}_{\alpha}\mathbf{D}\right)^{2}.$

For the electronic Hamiltonian we adapt the Kohn-Sham (KS) TDDFT Runge and
Gross (1984) description,

$H_{e}=-{\hbar^{2}\over 2m}\nabla^{2}+V_{KS}(\mathbf{r}),$ (3)

with

$V_{KS}(\mathbf{r})=V_{\mathrm{H}}[\rho](\mathbf{r})+V_{\mathrm{XC}}[\rho](\mathbf{r})+V_{\mathrm{ion}}(\mathbf{r}).$
(4)

Here, $\rho$ is the electron density, $V_{\mathrm{H}}$ is the Hartree
potential, and $V_{\mathrm{XC}}$ is the exchange-correlation potential (Local
density approximation (LDA) is used), and $V_{\mathrm{ion}}$ is the external
potential due to the ions. The potential of the ions can be represented by
employing the norm-conserving pseudopotentials of the form given by Troullier
and Martins Troullier and Martins (1991). We assume that the nuclei are fixed
in their positions.

One should emphasize at this point that the proper way of introducing the
Hamiltonian would be to use mapping theorems following the steps of Ref. Flick
et al. (2015); Pellegrini et al. (2015); Ruggenthaler et al. (2014). See also
Ref. Schafer et al. (2021) for recent developments. Our approach can be
thought to be using the conventional velocity gauge TDDFT Hamiltonian as a
starting point

$H_{V}={1\over 2}\left(-i\hbar\nabla-{e\over
c}\mathbf{A}\right)^{2}+V_{KS}(\mathbf{r}).$ (5)

Then we replace the classical vector potential with the quantized vector
potential of the cavity:

$\hat{\mathbf{A}}(\mathbf{r})=\sum_{\alpha}\left({\hbar\over
2\epsilon_{0}\omega_{\alpha}}\right)^{1/2}\left[\hat{a}_{\alpha}\mathbf{u}_{\alpha}(\mathbf{r})+\hat{a}_{\alpha}^{+}\mathbf{u}_{\alpha}^{*}(\mathbf{r})\right].$
(6)

In the dipole approximation, placing the molecule in the middle of the cavity
at $\mathbf{r}_{0}$ one has

$\hat{\mathbf{A}}(\mathbf{r}_{0})=\sum_{\alpha}\boldsymbol{\lambda}_{\alpha}\hat{a}_{\alpha},$
(7)

and

$H_{V}={1\over 2}\left(-i\hbar\nabla-{e\over
c}\hat{\mathbf{A}}(\mathbf{r}_{0})\right)^{2}+V_{KS}(\mathbf{r}).$ (8)

One can use a unitary transformation to transform the Hamiltonian into length
gauge Rokaj et al. (2018); Taylor et al. (2022)

$H_{L}=U^{\dagger}H_{V}U,$ (9)

with

$U=\exp\left\\{-{i\over\hbar}\hat{\mathbf{A}}(\mathbf{r}_{0})\cdot\mathbf{D}\right\\}.$
(10)

After a straightforward calculation Rokaj et al. (2018) one has

$H_{L}=-{\hbar^{2}\over
2m}\nabla^{2}+V_{KS}(\mathbf{r})+\sum_{\alpha=1}^{N_{p}}{1\over
2}\left[p_{\alpha}^{2}+\omega_{\alpha}^{2}\left(q_{\alpha}-{\boldsymbol{\lambda}_{\alpha}\over\omega_{\alpha}}\mathbf{D}\right)^{2}\right].$
(11)

Either $H_{V}$ or $H_{L}$ can be used in the calculations. The advantage of
$H_{V}$ is that one can use it with periodic boundary conditions. In either
case, the Hamiltonian depends not only on the spatial coordinates, but on the
photon creation and annihilation operators, as well and one has to use a
tensor product of a spatial and number state basis for the orbitals. We will
use $H=H_{L}$ in this work.

###  II.2 Basis functions

The $\omega_{\alpha}q_{\alpha}\boldsymbol{\lambda}_{\alpha}\mathbf{D}$ term
couples the electronic and photonic degrees of freedom. The orbitals of the
coupled electron and photon system at the KS level can be written

$\Phi_{mn}=\phi_{mn}(\mathbf{r})|n\rangle,\ \ \ \ \ (m=1,{\ldots},N_{occ}),\ \
\ \ \ (n=0,{\ldots},N_{F}),$ (12)

where $|n\rangle$ is the Fock space basis for the photons, $N_{F}$ is the
dimension of the Fock space, and $N_{occ}$ is the number of orbitals. In the
following we assume that there is one dominant photon mode – $N_{p}=1$ in Eq.
II.1 – with frequency $\omega$. Extension to $N_{p}>1$ is possible, but it
quickly leads to prohibitively large basis dimensions. We use a real-space
representation Chelikowsky et al. (1994); Yabana and Bertsch (1996); Tancogne-
Dejean et al. (2020); Varga et al. (2004) for the electronic part

$\phi_{mn}(\mathbf{r})=\phi_{mn}(x,y,z),$ (13)

and the orbitals will be defined on a 4 dimensional (4D) product grid of the 3
dimensional real-space grid and the 1 dimensional Fock space. The 4D grid has
$N_{x}\times N_{y}\times N_{z}\times(N_{F}+1)$ grid points, where
$N_{x},N_{y},N_{z}$ are the grid points in real-space and $N_{F}$ is the size
of the Fock basis. The orbitals are written as

$\Phi_{m}(\mathbf{r})=\left(\begin{array}[]{c}\Phi_{m0}(x_{i},y_{j},z_{k})\\\
\Phi_{m1}(x_{i},y_{j},z_{k})\\\ \vdots\\\ \Phi_{mn}(x_{i},y_{j},z_{k})\\\
\vdots\end{array}\right).$ (14)

Due to the orthogonality of the Fock basis states we have

$\left(\Phi_{mn}|\Phi_{m^{\prime}n^{\prime}}\right)=\langle\phi_{mn}|\phi_{m^{\prime}n}\rangle\delta_{nn^{\prime}},$
(15)

where the round bracket stands for the integration over the real and Fock
space and the angle bracket is the integration over the space part,

$\langle\phi_{mn}|\phi_{m^{\prime}n}\rangle=\sum_{ijk}\phi_{mn}(x_{i},y_{j},z_{k})\phi_{m^{\prime}n}(x_{i},y_{j},z_{k}).$
(16)

Using Gram-Schmidt orthogonalization the real-space functions
$\phi_{1n},{\ldots},\phi_{N_{occ}n}$ are orthogonalized for each Fock state
$(n=0,{\ldots},N_{F})$. The new orthogonal set
$\hat{\phi}_{1n},{\ldots},\hat{\phi}_{N_{occ}n}$ is normalized

$\sum_{n=0}^{N_{F}}|\hat{\phi}_{mn}|^{2}=1.$ (17)

Now one can define the electron density as

$\rho(\mathbf{r})=\sum_{m=1}^{N_{occ}}c_{m}\sum_{n=0}^{N_{F}}|\hat{\phi}_{mn}(\mathbf{r})|^{2},$
(18)

where $c_{m}$ is the number of electrons on the $m$th orbital. One can also
define the density belonging to a given Fock state as

$p_{n}(\mathbf{r})=\sum_{m=1}^{N_{occ}}c_{m}|\hat{\phi}_{mn}(\mathbf{r})|^{2},$
(19)

and the photon occupation probability in the Fock space

$P_{n}={1\over N}\int p_{n}(\mathbf{r})d{\mathbf{r}},$ (20)

where $N$ is the number of electrons.

The orbitals will be calculated by iterative minimization – conjugate gradient
in the present case – in the same way as in the conventional real-space
approaches. For this we need to calculate the action of the Hamiltonian on the
wave function. Noting that

$\omega\left(\hat{a}^{+}\hat{a}+{1\over 2}\right)|n\rangle=(n+{1\over
2})\omega|n\rangle,$ (21)

and

$\displaystyle q|n\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\omega}}\left(\hat{a}+\hat{a}^{+}\right)|n\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\omega}}\left(|\sqrt{n}|n-1\rangle+\sqrt{n+1}|n+1\rangle\right),$

one has

$H\Phi_{mn}(\mathbf{r})=-{1\over
2}\nabla^{2}\Phi_{mn}(\mathbf{r})+\left(V_{KS}(\mathbf{r})+\mu(\boldsymbol{\lambda}\mathbf{r})+\omega\left(n+{1\over
2}\right)\right)\Phi_{mn}(\mathbf{r})-\sqrt{\omega\over
2}(\boldsymbol{\lambda}\mathbf{r})\left(\sqrt{n}\Phi_{mn-1}(\mathbf{r})+\sqrt{n+1}\Phi_{mn+1}(\mathbf{r})\right),$
(23)

where

$\mu=\int\boldsymbol{\lambda}\mathbf{r}\rho(\mathbf{r})d\mathbf{r}.$ (24)

In the first part of this equation, the Hamiltonian acts on $\phi_{mn}$ in
each Fock space in the same way as in conventional real-space approaches. The
kinetic energy operator is represented by nine-point finite differencing, and
the nonlocal part of the pseudopotential is calculated by a summation in the
pseudopotential core radius around the atomic position, leading to a very
sparse Hamiltonian matrix. In the second part of the equation the photon
spaces $n$ and $n\pm 1$ are connected and the coupling Hamiltonian matrix is
diagonal. This large sparse system is ideal for iterative approaches.

###  II.3 Ground state calculation

The ground state calculation is similar to the conventional DFT calculations:

  1. 1.

Initialization of the orbitals. $\phi_{mn}^{(0)}$ are approximated e.g. by
atomic orbitals. The components belonging to the lowest ($n=0,1$) Fock spaces
are expected to be dominant so a weight factor is used to enhance those
components, $\phi^{(0)}_{mn}\rightarrow w_{n}\phi^{(0)}_{mn}$. Orthogonalize
the orbitals to generate the starting wave function

$\hat{\Phi}_{m}^{(0)}(\mathbf{r})=\left(\begin{array}[]{c}\hat{\Phi}^{(0)}_{m0}(\mathbf{r})\\\
\hat{\Phi}^{(0)}_{m1}(\mathbf{r})\\\ \vdots\\\
\hat{\Phi}^{(0)}_{mn}(\mathbf{r})\\\ \vdots\end{array}\right),$ (25)

where
$\hat{\Phi}^{(0)}_{mn}(\mathbf{r})=\hat{\phi}^{(0)}_{mn}(\mathbf{r})|n\rangle$.

  2. 2.

Iterative minimization step. Use

${\Phi}_{m}^{(k+1)}(\mathbf{r})=\sum_{j}a_{j}H^{j}\hat{\Phi}_{m}^{(k)}(\mathbf{r}),$
(26)

where the $a_{j}$ coefficients define the iterative procedure, e.g. steepest
descent, conjugate gradient, imaginary time-propagation, etc. Eq. 23 is used
to calculate $H\hat{\Phi}^{(k)}_{mn}$.

  3. 3.

Calculate $\hat{\Phi}_{m}^{(k+1)}(\mathbf{r})$ by Gram-Schmidt
orthogonalization ${\Phi}_{m}^{(k+1)}(\mathbf{r})$.

  4. 4.

Calculate $\rho$ and update $V_{KS}$.

  5. 5.

k=k+1, go to step 2 until convergence is reached.

###  II.4 Time propagation

The ground state orbitals will be used to initialize the time propagation

$\hat{\Phi}_{m}(\mathbf{r},t=0)=\hat{\Phi}_{m}(\mathbf{r}).$ (27)

Any time propagation method typically used for TDDFT can be used here, and in
this work we use Taylor time propagation Yabana and Bertsch (1996)

$\hat{\Phi}_{m}(\mathbf{r},t+\Delta t)={\rm e}^{-iH\Delta
t}\hat{\Phi}_{m}(\mathbf{r},t)=\sum_{j=0}^{4}{(-i\Delta t)^{j}\over
j!}H^{j}\hat{\Phi}_{m}(\mathbf{r},t),$ (28)

where $\Delta t$ is chosen to be sufficiently small to conserve the norm of
the orbitals during propagation.

##  III Applications

Atomic units are used in the calculations, and we drop the “a.u.” notation
after every number. The coupling will be defined as
$\boldsymbol{\lambda}=\lambda\mathbf{u}$ where $\mathbf{u}$ is a unit vector,
e.g. (1,0,0). $\lambda$ can be related to the effective cavity volume
$\lambda=1/\sqrt{\epsilon_{0}V_{eff}}$ ($\epsilon_{0}$ is the permittivity of
vacuum) Galego et al. (2019). Sub nm3 volumes have been reached in
picocavities Benz et al. (2016); Carnegie et al. (2018) corresponding to
$\lambda<0.05$, which will be the highest value that we consider in this work,
except for an example for benzene in which we calculate the wave function in
higher photon spaces.

###  III.1 Hydrogen molecule

Figure 1: Comparison of the results of the QED-DFT-TP, QED-CC and QED-SVM
calculations for energy dependence on bond length (a), $\omega$ (b) and
$\lambda$(c) . The QED-CC results are from the calculation of Ref. Mallory and
au2 (2022) using cc-pVQZ basis. $\omega=0.08$ is used for a and c,
$\lambda$=0.05 is used for a and b. $N_{x}=N_{y}=N_{z}=81$ and $h=0.5$ grid
spacing used in the calculations.

7

In the first example, we calculate the energy of a hydrogen molecule in a
cavity with fixed nuclei. Two-photon spaces, $|0\rangle$ and $|1\rangle$, are
used. Higher photon spaces have negligible contribution for this coupling
strength region. The results are compared to an SVM and a QED-CC calculation.
As mentioned in the introduction, both the SVM and the QED-CC are very
accurate wave function based approaches and are in complete agreement with
each other in these calculations. Fig. 1 shows the total energy as a function
of bond length (denoted as $d$) for $\lambda=0.05$
($\boldsymbol{\lambda}=(\lambda,0,0))$ and $\omega=0.081$. The molecular bond
lay parallel to the photon polarization vector. The parameters are chosen to
make a comparison to the QED-CC calculations of Ref. Mallory and au2 (2022).
The calculated QED-DFT-TP energy is shifted with respect to QED-CC and SVM
energies, due to the LDA and pseudopotential approximations, but the overall
trend is very similar. The only noticeable difference is in the 2 $\leq$ $d$
region where the QED-DFT-TP curve has a larger slope. This is mainly due to
the LDA approximation which performs well near the equilibrium separation, but
inaccurate when the bond is stretched as it is discussed e.g. in Ref. Becke
(2014).

At this coupling strength, the effect of the coupling on the ground state is
mainly an energy shift. It is more interesting to study the $\omega$ and
$\lambda$ dependence for a given bond length. This is shown in Figs 1 b and c.
The SVM and QED-CC results are in excellent agreement, and the QED-DFT-TP
results show the same behavior except for the energy shift.

Fig. 2 compares photon occupation number $P_{n}$ as a function of bond length
$\lambda$ and $\omega$. $N_{F}=1$ is used, so $P_{0}+P_{1}=1$ and only $P_{0}$
is shown in the figures. In the case of bond length, SVM and QED-DFT agree
well up to $d=1.7$. For larger distances, the two probabilities are different,
probably due to the aforementioned LDA issue. The SVM and QED-DFT-TP
description of the $\lambda$ and $\omega$ dependence agrees well.

Figure 2: Photon occupation number as a function of bond length (a), $\omega$
(b) and $\lambda$(c) . The QED-CC results are from Mallory and au2 (2022).
$\omega=0.08$ is used for a and c, $\lambda$=0.05 is used for a and b.
$N_{x}=N_{y}=N_{z}=81$ and $h=0.5$ used in the calculations.

###  III.2 The LiH and the HF molecules

Figs. 3a and b shows the $\lambda$ and $\omega$ dependence of the energy of
the LiH molecule. The $\omega$ dependence calculated by QED-SVM and QED-DFT-TP
are very similar. Due to the pseudopotential approach the total energies
calculated by QED-DFT-TP and QED-SVM are different and the QED-SVM energies
are shifted up by a constant to fit them in the same figure. We chose a larger
$\omega$ ($\omega=0.5)$ for the calculation of the $\lambda$ dependence to
show the importance of the coupling of the $|0\rangle$ and $|1\rangle$ photon
spaces. Fig. 3b shows that the energy in the $|0\rangle$+$|1\rangle$ case is
significantly lower than in the $|0\rangle$ case for higher $\lambda$ values
(in the $|0\rangle$ space only the dipole self interaction is added to the
Hamiltonian). The $\lambda$ dependence calculated by QED-SVM and the QED-DFT-
TP are in good agreement.

Figs. 3c and d shows the $\lambda$ and $\omega$ dependence for a HF molecule.
The $\omega$ dependence is very similar in both the QED-DFT-TP and QED-CC
cases. The energy shift is due to the pseudopotential treatment of the core
electrons in the QED-DFT-TP case. In this case the weight of the $|0\rangle$
space is 95% and the weight of the $|1\rangle$ space is 5% for the largest
$\lambda$, so in overall the $|0\rangle$ space dominates but the $|1\rangle$
space also contributes to the energy as Fig. 3d shows.

The $\lambda$ dependence only agrees in the tendency that by increasing
$\lambda$ the energy is increased, but the $\lambda$ dependence is much larger
in the QED-DFT-TP case.

###  III.3 Na2 cluster

The next example is a Na2 cluster with a bond length $d$=5.8. After the ground
state calculation the initial wave function is excited with a delta kick (see
Appendix D for details) and then time propagated. The photon coordinate $q(t)$
and the dipole moment $D_{x}(t)$ are shown in Fig. 4. We have calculated these
quantities by the QED-TDDFT-TP and with the QEDFT approaches (see Appendix B
for details). The time dependence of $q(t)$ calculated by the two methods
differ greatly. This is understandable: in the present calculation $q(t)$ is
calculated by solving the time-dependent Schrödinger equation and calculating
the expectation value from the time-dependent wave function comprising
different photon spaces. In QEDFT, the expectation value of $q(t)$ is
calculated using Eq. (39). The dipole moments calculated by the two approaches
show some similarity. But these quantities are also calculated in different
ways: one is solving Eq. (38) the other is using Eq. (23) in the coupled
orbital and photon spaces.

What we are really interested in is the excitation spectra shown Fig. 4c. The
single peak of the Na2 spectrum when no cavity is present splits into two
polariton peaks in the cavity in both QEDFT and the present calculation. The
QEDFT predicts a somewhat bigger splitting but both approaches show a larger
upper polariton peak. The position of the disodium peak while not in a cavity
seems to be at the position of the lower polariton peak of the present
calculation. This is just a coincidence: the position of the single peak
shifts to higher energies with nonzero $\lambda$, and the present calculation
also predicts upper and lower polaritons with respect to the single peak. We
will show this next.

Fig. 4d shows the evolution of the Rabi splitting as a function of $\lambda$.
Using the $|0\rangle$ space only, the dipole self-interaction changes the
energy of the system and changes the position of the peak in the absorption
cross section (Fig. 4d) compared to the ”no cavity” case (no dipole self-
interaction). By coupling the $|0\rangle$ and $|1\rangle$ photon spaces, the
Rabi splitting appears in the form of two peaks: one above and one below the
single peak of the $|0\rangle$ photon space case. Increasing $\lambda$ the
splitting is increasing as it is expected. The figure also shows the
contribution from the $|0\rangle$ component in the $|0\rangle$+$|1\rangle$
coupled case. For small $\lambda$ the $|0\rangle$ component is dominant (but
there is no splitting without the coupling), but for larger $\lambda$ (see
Fig. 4c for $\lambda=0.05$) the $|1\rangle$ space significantly contributes;
in this case to the lower peak. Note that the width of the peaks is related to
the length of the simulation time. The present model excludes physical
mechanisms of line broadening. More time steps at a fixed temporal resolution
result in narrower, Dirac-delta like peaks.

Figure 3: Dependence of the energy of LiH and HF molecules on $\omega$ and
$\lambda$. The Li-H distance is 3.015 and the H-F distance is 1.7229,
$\lambda=0.05$ was used to calculate the $\omega$ dependence for both
molecules. For LiH $\omega=0.5$, for HF $\omega=0.0813$ was used to calculate
the $\lambda$ dependence. $N_{x}=N_{y}=N_{z}=51$ grid points is used with
$h=0.5$ grid spacing. The QED-CC results are from Mallory and au2 (2022). In
the case of the LiH molecule the QED-SVM results are shifted up by 7.3.

Figure 4: QED-TDDFT-TP calculation for the Na2 dimer. (a) Time dependence of
the photon displacement coordinate. (b) Time dependence of the dipole moment.
(c) Rabi splitting calculated by the QEDFT and by the present approach.
$|0\rangle$ is the contribution to the absorption cross section from the 0
photon space in a calculation using $|0\rangle+|1\rangle$ space,
$|0\rangle+|1\rangle$ is the total absorption cross section.
$\boldsymbol{\lambda}=(\lambda,0,0)$, $\lambda=0.05$. (d) Dependence of the
Rabi splitting on $\lambda$. DSI is the dipole self-interaction contribution.
In this case the molecule is in a cavity but coupled only to the $|0\rangle$
photon space. The Na2 dimer is oriented along the $x$ direction and the two
atoms are d=5.8 distance apart. $N_{x}=N_{y}=N_{z}=51$, h=0.5, dt=0.07,
$N_{t}$=10000 d=5.8. $\omega=0.07$ is used in the calculations, $N_{t}$=30000.

###  III.4 Benzene and p-nitroaniline molecules

In this section, we show how the electron density changes in different photon
spaces in cavities Riso et al. (2022). The first example is a benzene molecule
(see Fig. 5). In this case we try a larger $\lambda$ to couple higher photon
spaces in the wave function. The $|0\rangle$, $|1\rangle$, and $|2\rangle$
spaces are coupled and the weights of these components are 0.94, 0.056, and
0.004. The electron density of the benzene in the $|0\rangle$ photon space
(Fig. 5a) is very similar to the ground state electron density. The density in
the $|1\rangle$ photon space (Fig. 5b) is somewhat more compact than that in
$|0\rangle$ and a density hole appears in the middle of the ring. The
$|2\rangle$ density (Fig. 5c) is split perpendicular to the polarization
vector $\boldsymbol{\lambda}=(\lambda,\lambda,0)$. This split is the
consequence of the symmetry breaking introduced by the $\boldsymbol{\lambda}$
dependent parts of the $H_{ep}$ term in Eq. (II.1).

Figs 5d-i show the HOMO-2 and the HOMO orbitals. The orbitals in the
$|0\rangle$ photon space remain very similar to the ground state orbitals (the
slight modification due to the dipole self-interaction is not visible in the
figures). The $|1\rangle$ and the $|2\rangle$ components of the HOMO-2 and
HOMO orbitals (Figs. 5e, f, and 5g, i) are significantly different from the
$|0\rangle$ component. In the case of HOMO-2 only the charge distribution
changes and the shape of the orbital components are similar, but in the case
of the HOMO the lobes are also different. Other orbitals (not shown) have
similar shape changes.

The next example is a p-nitroaniline molecule which was studied by QED-CC in
Ref. Haugland et al. (2020) and we have adapted the parameters from that work.
This is an example of a molecule with a low lying charge transfer excitation.
To measure the charge transfer one can define Haugland et al. (2020)

$\Delta
q(x)=\int_{-\infty}^{x}dx\int_{-\infty}^{\infty}dy\int_{-\infty}^{\infty}dz\Delta\rho(x,y,z),$
(29) $\Delta\rho(x,y,z)=\rho^{cavity}(x,y,z)-\rho^{nocavity}(x,y,z).$ (30)

$\Delta q(x)$ measures the amount of charge moved in the $x$ direction – the
principle axis of the molecule. The density difference induced in the molecule
by the cavity is shown in Fig. 6. The charge redistribution is very similar to
what is shown in Ref. Haugland et al. (2020). Electrons transfer from the
nitro (NO2) side to the aniline (NH2) side. The $\Delta q(x)$ function shown
in Fig. 6 is also very similar to QED-CC calculation in Ref. Haugland et al.
(2020).

Figure 5:  Electron density ISO surfaces of benzene in the (a) $|0\rangle$,
(b) $|1\rangle$, (c) $|2\rangle$ photon spaces. The three shades of the ISO
surfaces have values of 0.01, 0.001, and 0.0001 HOMO-2 orbital isosurfaces of
benzene in the (d) $|0\rangle$, (e) $|1\rangle$, (f) $|2\rangle$ photon
spaces. HOMO orbital isosurfaces of benzene in the (g) $|0\rangle$, (h)
$|1\rangle$, (i) $|2\rangle$ photon spaces. $N_{x}=N_{y}=N_{z}=51$ grid points
are used with $h=0.5$ grid spacing, $\boldsymbol{\lambda}=(0.1,0.1,0)$,
$\omega=0.5$. The darker (magenta) color shows the positive and the lighter
(cyan) color shows the negative values. The two shades of the isosurfaces have
values of 0.0005 and 0.0002.

Figure 6:  (top)The cavity induced density difference between of a
p-nitroaniline molecule in cavity shown as isosurfaces. The two shades of the
isosurfaces have values of 0.0005 and 0.0002. (bottom) Charge distribution of
a p-nitroaniline molecule in cavity. $\boldsymbol{\lambda}=(0.05,0,0)$,
$\omega=0.178$ values were used in the calculation. $N_{x}=N_{y}=N_{z}=55$
grid points is used with $h=0.5$ grid spacing. The lighter (cyan) and darker
(magenta) regions represent charge depletions and accumulations,respectively.
Figure 7: Rabi splitting for benzene for different $\lambda$ values.
$N_{x}=N_{y}=N_{z}=55$ grid points is used with $h=0.5$ grid spacing. The time
step is 0.05 and the total propagation time is 1000.
$\boldsymbol{\lambda}=(0,\lambda,0)$ and $\omega=0.25$ was used in the
calculations.

We have also calculated the optical absorption spectrum of the benzene
molecule in a cavity by time propagating the dipole moment after an initial
delta kick perturbation. The Rabi splittings for benzene is shown in Fig. 7.
The calculated splittings are 0.22 (0.33), 0.63 (0.69), 0.81 (1.02), 1.03
(1.35) for $\lambda=0.01,0.02,0.03,0.04$, respectively. The numbers in the
parenthesis are from Ref. Flick et al. (2019). Just as in the case of the Na2
cluster our QED-TDDFT-TP approach gives somewhat smaller Rabi splitting than
the QEDFT approach, but the overall trend, a direct proportionality between
the splitting and $\lambda$, is the same.

Figure 8: High harmonics spectrum of the HF molecule (black lines). The red
dashed lines show the contribution from the $|1\rangle$ photon space. The blue
dotted line shows the high harmonics of the HF molecule without the cavity.
The solid vertical lines show the position of $n\omega_{L}$, the dashed
vertical lines are at $(n+1/2)\omega_{L}$ (where $n$ is an integer). The axis
of the molecule molecule is in the $x$ direction, the H-F distance is 1.7229,
$N_{x}=N_{y}=N_{z}=51$, $h=0.5$, $\boldsymbol{\lambda}=(0.05,0,0)$, $\Delta
t=0.05$ and $\omega=0.081$ is used. The laser is polarized in the $x$
direction and the laser parameters are $\omega_{L}$=0.057 (800 nm wavelength),
$T_{L}=2/\omega_{L}$ and $E_{x}$=0.005.

###  III.5 High harmonic generation

As a final example we present the high harmonic spectrum of the HF molecule
with and without a cavity. The molecule is excited with a continuous laser
pulse,

$E(\mathbf{r},t)=(E_{x},0,0)\sin(\pi t/(6T_{L}))^{2}\sin(\omega_{L}t)$ (31)

and after 30 laser cycles the high harmonic spectrum is calculated using the
dipole acceleration:

$I(\omega)=\left|\int_{0}^{T}{\partial^{2}\mathbf{D}(t)\over\partial
t^{2}}{\rm e}^{-i\omega t}dt\right|^{2}.$ (32)

Fig. 8 shows the calculated high harmonics. We are mostly interested in the
first few harmonics and use a relatively small computational box to speed up
the calculations. To resolve more harmonics one needs a larger box, which
increases the computational cost.

If no cavity is present, the harmonic spectrum of the HF molecule oriented in
the direction of the polarization of the laser consists of peaks at
$n\omega_{L}$ (see Fig 8 for the $\omega=0$ case). Since HF does not have
inversion symmetry, both the even and odd harmonics are present. Placing the
molecule in a cavity using the $|0\rangle$ photon space only, the high
harmonic spectrum of the HF molecule does not change (the solid black and
dotted blue lines are barely distinguishable in Fig 8).

Using the $|0\rangle+|1\rangle$ photon spaces new peaks appear in the harmonic
spectrum at $\omega+n\omega_{L}$ positions. For example if
$\omega=0.5\omega_{L}$ then the new peaks are at $0.5\omega_{L}+n\omega_{L}$
(Fig. 8). These peaks have significant intensity, and the peaks are present at
weaker and stronger couplings as well. If $\omega=n\omega_{L}$ then the new
(cavity induced) peaks coincide with peaks of the laser and the spectrum
remains very similar to the harmonic spectrum without a cavity. Fig. 8 also
shows the contribution from the $|1\rangle$ photons spaces. The dipole moment
is the sum of the components from the $|0\rangle$ and $|1\rangle$ photon
spaces, and the harmonic spectrum is a result of an interference of these
terms. Fig 8 shows that in the case of $\omega=(k+1/2)\omega_{L}$ (k=0,1,2…)
the $|1\rangle$ component has pronounced peaks at $\omega$,$\omega+\omega_{L}$
and $\omega+2\omega_{L}$. The $|1\rangle$ photon space not only induces new
peaks, but significantly contributes to their intensity as well.

##  IV Summary

We have developed and implemented a QED-TDDFT-TP approach that uses a direct
product of a Fock space for light and a real-space grid for electrons. We have
tested the approach by comparing it to accurate QED-SVM and QED-CC
calculations and the results are very promising. The coupling strength and
light frequency dependence of energies calculated by the present approach is
very similar to the QED-SVM and QED-CC results.

One of the main advantages of the approach is that it preserves the quantized
photon spaces giving direct access to non-classical photon observables. A
further important advantage is that as the electron and the photon coordinates
are separated one can asses the contribution to physical quantities from
different photon spaces.

The light-matter coupling strongly modifies the molecular orbitals, and the
orbitals have very different shapes in different photon spaces. This also
leads to charge transfer and redistribution, and our results are similar to
the QED-CC calculation Haugland et al. (2020).

We have also studied the Rabi splitting in molecules in optical cavities and
the results show that the splitting is caused by the presence of light states:
there is no splitting without the presence of $|0\rangle$ and $|1\rangle$
photon spaces.

The approach can also be used to study high harmonic generation in strongly
coupled light-matter coupled systems. There are several recent experimental
and theoretical works studying the influence of polaritons on nonlinear
optical processes Kockum et al. (2017); Chervy et al. (2016); Barachati et al.
(2018); Liu et al. (2019); Welakuh and Narang (2022). We plan to extend the
present study to describe the high harmonic generation for more complex
molecules. We also plan to calculate the nonlinear susceptibilities of
molecules in cavities using combining the present approach with the method
described in Ref. Goncharov and Varga (2012).

The tensor product of the real-space and Fock-space increases by the number of
Fock spaces. For experimentally realizable strong coupling one needs $N_{F}$=2
or $N_{F}$=3 Fock space dimensions for a single photon frequency. This
increases the cost of the calculation in time and memory requirements by a
factor of $N_{F}$. If there are more than one photon frequency is important
then the calculation cost increases by $N_{F}^{N_{p}}$ where $N_{p}$ is the
number of relevant photon modes. For larger $N_{p}$ this will be prohibitively
expensive and one needs alternatives, such as the QEDFT approach, which can
handle any number of photons.

The increased dimensionality is not the only complicating factor. The
Hamiltonian (23) (except from a trivial $n\omega$ shift) acts the same way on
every real-space component $\phi_{mn}(\mathbf{r},t)$ for a given photon
number. If there would be no coupling in the Hamiltonian then one would have
the same set of $\phi_{mn}(\mathbf{r},t)$, ($m=1,{\ldots},N_{occ}$) orbitals
in each $|n\rangle$ photon space. In the self-consistent iterative
diagonalization of the nonlinear Hamiltonian, one has to carefully mix the
density so that the photon components are gradually optimized. The process is
very sensitive to the coupling (the off-diagonal part of the Hamiltonian),
because the coupling makes the orbitals different in different photon spaces.
Sudden changes in the density leads to slow damped oscillatory convergence.
Future work is necessary on optimizing the self-consistent solution for this
tensor product approach to have efficient calculations.

The development of accurate exchange-correlation functionals are also
important for future calculations. In the present work we have compared our
calculations to QED-SVM and QED-CC results, and these approaches may help in
developing appropriate functionals for QED-TDDFT-TP.

###### Acknowledgements.

This work has been supported by the National Science Foundation (NSF) under
Grant No. IRES 1826917. J.M., A.A. and D.P. equally contributed to this work.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflict of interest to disclose.

##  Appendix A Photon space

In this Appendix we give a brief overview of the properties of the Fock space
and occupation number basis. The photon Hamiltonian (from Eq. (II.1))

$H_{p}=\sum_{\alpha=1}^{N_{p}}{1\over
2}\left(p_{\alpha}^{2}+\omega_{\alpha}^{2}q_{\alpha}^{2}\right)=\sum_{\alpha=1}^{N_{p}}\omega_{\alpha}\left(\hat{a}_{\alpha}^{+}\hat{a}_{\alpha}+\frac{1}{2}\right),$
(33)

where

$q_{\alpha}={\frac{1}{\sqrt{2\omega_{\alpha}}}}(\hat{a}_{\alpha}+\hat{a}^{+}_{\alpha})\
\ \ \
p_{\alpha}={\frac{1}{i\sqrt{2\omega_{\alpha}}}}(\hat{a}_{\alpha}-\hat{a}^{+}_{\alpha}).$
(34)

Denoting the vacuum state as $|0\rangle$, any eigenstate of
$\hat{a}_{\alpha}^{+}\hat{a}_{\alpha}$ can be calculated by multiple
applications of the creation operators

$|n_{\alpha}\rangle={1\over\sqrt{n!}}(\hat{a}_{\alpha}^{+})^{n}|0\rangle.$
(35)

This is called the occupation number representation, $n_{\alpha}$ defines the
excited state of the $\omega_{\alpha}$ photon mode. The bosonic operators
satisfy the commutation relations

$[\hat{a}_{\alpha},\hat{a}_{\alpha^{\prime}}]=[\hat{a}_{\alpha}^{+},\hat{a}_{\alpha^{\prime}}^{+}]=0,\
\ \ \ \
[\hat{a}_{\alpha},\hat{a}_{\alpha^{\prime}}^{+}]=\delta_{\alpha,\alpha^{\prime}}\hat{1}$
(36)

where $[a,b]=ab-ba$.

Now one can define an occupation number basis

$\chi_{\vec{n}}=|n_{1},n_{2},{\ldots},N_{p}\rangle={1\over\sqrt{n_{1}!{\ldots}n_{N_{p}}!}}(\hat{a}_{1}^{+})^{n_{1}}{\ldots}(\hat{a}_{N_{p}}^{+})^{n_{N_{p}}}|0\rangle$
(37)

where $n=n_{1}+{\ldots}N_{p}$. In this abstract space representation the
symmetry requirement for bosons is satisfied by allowing that mode $i$ is
occupied by $n_{i}$ photons. For the product state in Eq. (37) the restriction
to the symmetric subspace is tacitly assumed, that is one has to symmetrize
the states for identical harmonic oscillators.

##  Appendix B QEDFT

In the approach of Refs. Schäfer et al. (2019); Ruggenthaler et al. (2018);
Flick et al. (2017); Sidler et al. (2021); Flick et al. (2019) two coupled
equation are solved

$i\hbar\frac{\partial}{\partial t}\varphi_{i}(\mathbf{r},t)=(-{1\over
2}\nabla^{2}+V_{KS}(\mathbf{r},t)+V_{P}(\mathbf{r},q,t))\varphi_{i}(\mathbf{r},t)$
(38) $\left(\frac{\partial^{2}}{\partial
t^{2}}+\omega_{\alpha}^{2}\right)q_{\alpha}(t)=-\frac{j_{\alpha}(t)}{\omega_{\alpha}}+\omega_{\alpha}\boldsymbol{\lambda}_{\alpha}\cdot\mathbf{R}(t)$
(39)

where $\mathbf{R}(t)=\int d\mathbf{r}e\mathbf{r}\rho(\mathbf{r},t)$. The
photon exchange potential is

$V_{\mathrm{P}}(\mathbf{r},q,t)=\sum_{\alpha=1}^{M}\left(\int
d\mathbf{r}^{\prime}\boldsymbol{\lambda}_{\alpha}\cdot\mathbf{r}^{\prime}\rho\left(\mathbf{r}^{\prime},t\right)-\omega_{\alpha}q_{\alpha}(t)\right)\boldsymbol{\lambda}_{\alpha}\mathbf{r}$
(40)

The photon is propagated as

$\left({\partial^{2}\over\partial
t^{2}}+\omega_{\alpha}^{2}\right)q_{\alpha}(t)=\omega_{\alpha}\boldsymbol{\lambda}_{\alpha}\mathbf{R}(t)$
(41)

##  Appendix C Stochastic variational approach

In this case the full Coulombic Hamiltonian is used, the Hamiltonian of an $N$
electron system interacting with a Coulomb interaction and confined in an
external potential $V_{c}$ is

$H_{e}=-\sum_{i=1}^{N}{{\nabla^{2}_{i}}\over
2m_{i}}+\sum_{i<j}^{N}{q_{i}q_{j}\over|\mathbf{r}_{i}-\mathbf{r}_{j}|}+\sum_{i=1}^{N}V_{c}(\mathbf{r}_{i}),$
(42)

where $\mathbf{r}_{i},q_{i}$, and $m_{i}$ are the coordinate, charge, and mass
of the $i$th particle ($m_{i}=1$ for electrons in atomic units).

Introducing the shorthand notations
$\vec{r}=(\mathbf{r}_{1},...,\mathbf{r}_{N})$, the variational trial wave
function is written as a linear combination of products of spatial and photon
space basis functions

$\Psi(\vec{r})=\sum_{n}\sum_{k=1}^{K_{n}}c_{k}^{n}\psi_{k}^{n}(\vec{r})|n\rangle.$
(43)

The spatial part of the wave function is expanded into ECGs for each photon
state $|n\rangle$ as

$\psi_{k}^{n}(\vec{r})={\cal A}\\{{\rm e}^{-{1\over
2}\sum_{i<j}^{N}\alpha_{ij}^{k}(\mathbf{r}_{i}-\mathbf{r}_{j})^{2}-{1\over
2}\sum_{i=1}^{N}\beta_{i}^{k}(\mathbf{\mathbf{r}}_{i}-\mathbf{s}_{i}^{k})^{2}}\chi_{S}\\}$
(44)

where ${\cal A}$ is an antisymmetrizer, $\chi_{S}$ is the $N$ electron spin
function (coupling the spin to $S$), and $\alpha_{ij}^{k}$,$\beta_{i}^{k}$ and
$\mathbf{s}_{i}^{k}$ are nonlinear parameters. The basis functions are
optimized as described in Ref. Ahrens et al. (2021).

##  Appendix D Absorption cross section

The polarizability tensor $\boldsymbol{\alpha}(\boldsymbol{\omega})$ can be
calculated by time propagation of the electron orbitals states. The initial
state is first perturbed by a delta kick potential

$V_{ext}=-e\mathbf{r}\cdot\mathbf{k}_{i}\delta\left(t_{0}\right),$ (45)

where $\mathbf{k}_{i}$ is the electric field in the $i$ direction ($i=x,y,z$.
The magnitude of the electric field should be sufficiently small to initiate a
linear-response.

The polarizability tensor in frequency space is defined as the Fourier
Transform of the time-dependent dipole moment

$\alpha_{ij}(\omega)=\frac{1}{k_{i}}\int_{0}^{\infty}dt\left[D_{j}(t)-D_{j}\left(t_{0}\right)\right]e^{-i\omega
t}.$ (46)

The time-dependent dipole moments can be calculated from the time-dependent
electronic density $(\rho(\boldsymbol{r},t)))$ :

$\mathbf{D}(t)=\int\mathbf{r}\cdot\rho(\boldsymbol{r},t)d\mathbf{r}$ (47)

The imaginary part of the dynamic polarizability tensor can be related to the
oscillator strength of each transition Ullrich (2012)

$\displaystyle\frac{1}{3}\operatorname{Tr}(\Im m\boldsymbol{\alpha}(\omega))$
$\displaystyle=$
$\displaystyle\frac{\pi}{3}\sum_{n=1}^{\infty}\sum_{\nu=1}^{3}\left|\left\langle\Psi_{0}\left|\hat{r}_{\nu}\right|\Psi_{n}\right\rangle\right|^{2}\delta\left(\omega-\Omega_{n}\right)$
(48) $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{\pi}{2\Omega_{n}}f_{n}\delta\left(\omega-\Omega_{n}\right),$

where the $\left|\Psi_{n}\right\rangle$ is the wave function of the $n$-th
excited state, $\Omega_{n}$ is the corresponding excitation energy, and
$f_{n}$ is the oscillator strength belonging to the transition probability for
the excitation. The dynamic polarization tensor
$(\boldsymbol{\alpha}(\omega))$ and the transition dipole moment are now
connected as
$\left(\mathbf{D}_{n0}=\left\langle\Psi_{0}|\hat{\boldsymbol{r}}|\Psi_{n}\right\rangle\right).$

From the imaginary part of the dynamic polarizability tensor one can extract
the photo-absorption cross-section:

$\sigma(\omega)=\frac{4\pi\omega}{3c}\operatorname{Tr}\left[\Im
m(\boldsymbol{\alpha}(\omega)\right]=\frac{2\pi^{2}}{c}S(\omega)$ (49)

where $c$ is the speed of light and
$S(\omega)=\sum_{n=1}^{\infty}f_{n}\delta\left(\omega-\Omega_{n}\right)$ is
the dipole spectral function.

A damping function is added to the Fourier transform of the photo-absorption
spectrum for finite propagation times. The damping function introduces an
artificial decay of the excited population to smooth spectrum. We use a
Gaussian damping function $\left(e^{-\eta^{2}t^{2}}\right)$ and the parameter
$\eta$ is chosen to make the damping very small at the end of propagation.

## References

  * Hutchison et al. (2012) J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux, and T. W. Ebbesen, Angewandte Chemie International Edition 51, 1592 (2012), eprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/anie.201107033, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/anie.201107033.
  * Balili et al. (2007) R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, and K. West, Science 316, 1007 (2007), ISSN 0036-8075, eprint https://science.sciencemag.org/content/316/5827/1007.full.pdf, URL https://science.sciencemag.org/content/316/5827/1007.
  * Schachenmayer et al. (2015) J. Schachenmayer, C. Genes, E. Tignone, and G. Pupillo, Phys. Rev. Lett. 114, 196403 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.196403.
  * Xiang et al. (2020) B. Xiang, R. F. Ribeiro, M. Du, L. Chen, Z. Yang, J. Wang, J. Yuen-Zhou, and W. Xiong, Science 368, 665 (2020), ISSN 0036-8075, eprint https://science.sciencemag.org/content/368/6491/665.full.pdf, URL https://science.sciencemag.org/content/368/6491/665.
  * Reserbat-Plantey et al. (2021) A. Reserbat-Plantey, I. Epstein, I. Torre, A. T. Costa, P. A. D. Gonçalves, N. A. Mortensen, M. Polini, J. C. W. Song, N. M. R. Peres, and F. H. L. Koppens, ACS Photonics 8, 85 (2021), URL https://doi.org/10.1021/acsphotonics.0c01224.
  * Coles et al. (2014) D. M. Coles, N. Somaschi, P. Michetti, C. Clark, P. G. Lagoudakis, P. G. Savvidis, and D. G. Lidzey, Nature Materials 13, 712 (2014), ISSN 1476-4660, URL https://doi.org/10.1038/nmat3950.
  * Kasprzak et al. (2006) J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, et al., Nature 443, 409 (2006), ISSN 1476-4687, URL https://doi.org/10.1038/nature05131.
  * Schwartz et al. (2011) T. Schwartz, J. A. Hutchison, C. Genet, and T. W. Ebbesen, Phys. Rev. Lett. 106, 196405 (2011), URL https://link.aps.org/doi/10.1103/PhysRevLett.106.196405.
  * Plumhof et al. (2014) J. D. Plumhof, T. Stöferle, L. Mai, U. Scherf, and R. F. Mahrt, Nature Materials 13, 247 (2014), ISSN 1476-4660, URL https://doi.org/10.1038/nmat3825.
  * Hutchison et al. (2013) J. A. Hutchison, A. Liscio, T. Schwartz, A. Canaguier-Durand, C. Genet, V. Palermo, P. Samorì, and T. W. Ebbesen, Advanced Materials 25, 2481 (2013), eprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.201203682, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/adma.201203682.
  * Wang et al. (2021a) K. Wang, M. Seidel, K. Nagarajan, T. Chervy, C. Genet, and T. Ebbesen, Nature Communications 12, 1486 (2021a), ISSN 2041-1723, URL https://doi.org/10.1038/s41467-021-21739-7.
  * Basov et al. (2021) D. N. Basov, A. Asenjo-Garcia, P. J. Schuck, X. Zhu, and A. Rubio, Nanophotonics 10, 549 (2021), URL https://doi.org/10.1515/nanoph-2020-0449.
  * Feist and Garcia-Vidal (2015) J. Feist and F. J. Garcia-Vidal, Phys. Rev. Lett. 114, 196402 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.196402.
  * Vendrell (2018) O. Vendrell, Phys. Rev. Lett. 121, 253001 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.253001.
  * Schäfer and Johansson (2022) C. Schäfer and G. Johansson, Phys. Rev. Lett. 128, 156402 (2022), URL https://link.aps.org/doi/10.1103/PhysRevLett.128.156402.
  * Riso et al. (2022) R. R. Riso, T. S. Haugland, E. Ronca, and H. Koch, Nature Communications 13, 1368 (2022), ISSN 2041-1723, URL https://doi.org/10.1038/s41467-022-29003-2.
  * Badankó et al. (2022) P. Badankó, O. Umarov, C. Fb́ri, G. J. Halász, and A. Vibók, International Journal of Quantum Chemistry 122, e26750 (2022), eprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.26750, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/qua.26750.
  * Mazza and Georges (2019) G. Mazza and A. Georges, Phys. Rev. Lett. 122, 017401 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.122.017401.
  * Di Stefano et al. (2019) O. Di Stefano, A. Settineri, V. Macrì, L. Garziano, R. Stassi, S. Savasta, and F. Nori, Nature Physics 15, 803 (2019), ISSN 1745-2481, URL https://doi.org/10.1038/s41567-019-0534-4.
  * Galego et al. (2015) J. Galego, F. J. Garcia-Vidal, and J. Feist, Phys. Rev. X 5, 041022 (2015), URL https://link.aps.org/doi/10.1103/PhysRevX.5.041022.
  * Herrera and Spano (2016) F. Herrera and F. C. Spano, Phys. Rev. Lett. 116, 238301 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.238301.
  * Galego et al. (2016) J. Galego, F. J. Garcia-Vidal, and J. Feist, Nature Communications 7, 13841 (2016), ISSN 2041-1723, URL https://doi.org/10.1038/ncomms13841.
  * Shalabney et al. (2015) A. Shalabney, J. George, J. Hutchison, G. Pupillo, C. Genet, and T. W. Ebbesen, Nature Communications 6, 5981 (2015), ISSN 2041-1723, URL https://doi.org/10.1038/ncomms6981.
  * Buchholz et al. (2019) F. Buchholz, I. Theophilou, S. E. B. Nielsen, M. Ruggenthaler, and A. Rubio, ACS Photonics 6, 2694 (2019).
  * Schäfer et al. (2019) C. Schäfer, M. Ruggenthaler, H. Appel, and A. Rubio, Proceedings of the National Academy of Sciences 116, 4883 (2019), ISSN 0027-8424, eprint https://www.pnas.org/content/116/11/4883.full.pdf, URL https://www.pnas.org/content/116/11/4883.
  * Ruggenthaler et al. (2018) M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel, and A. Rubio, Nature Reviews Chemistry 2, 0118 (2018), ISSN 2397-3358, URL https://doi.org/10.1038/s41570-018-0118.
  * Flick et al. (2015) J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Proceedings of the National Academy of Sciences 112, 15285 (2015), ISSN 0027-8424, eprint https://www.pnas.org/content/112/50/15285.full.pdf, URL https://www.pnas.org/content/112/50/15285.
  * Flick et al. (2017) J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Proceedings of the National Academy of Sciences 114, 3026 (2017), ISSN 0027-8424, eprint https://www.pnas.org/content/114/12/3026.full.pdf, URL https://www.pnas.org/content/114/12/3026.
  * Lacombe et al. (2019) L. Lacombe, N. M. Hoffmann, and N. T. Maitra, Phys. Rev. Lett. 123, 083201 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.083201.
  * Mandal et al. (2020a) A. Mandal, S. Montillo Vega, and P. Huo, The Journal of Physical Chemistry Letters 11, 9215 (2020a), pMID: 32991814.
  * Flick et al. (2019) J. Flick, D. M. Welakuh, M. Ruggenthaler, H. Appel, and A. Rubio, ACS Photonics 6, 2757 (2019).
  * Latini et al. (2019) S. Latini, E. Ronca, U. De Giovannini, H. Hübener, and A. Rubio, Nano Letters 19, 3473 (2019), pMID: 31046291.
  * Flick et al. (2018a) J. Flick, N. Rivera, and P. Narang, Nanophotonics 7, 1479 (2018a), URL https://doi.org/10.1515/nanoph-2018-0067.
  * Flick and Narang (2018) J. Flick and P. Narang, Phys. Rev. Lett. 121, 113002 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.113002.
  * Garcia-Vidal et al. (2021) F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, Science 373 (2021), ISSN 0036-8075, eprint https://science.sciencemag.org/content/373/6551/eabd0336.full.pdf, URL https://science.sciencemag.org/content/373/6551/eabd0336.
  * Thomas et al. (2019) A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. A. Vergauwe, J. George, T. Chervy, A. Shalabney, E. Devaux, C. Genet, J. Moran, et al., Science 363, 615 (2019), ISSN 0036-8075, eprint https://science.sciencemag.org/content/363/6427/615.full.pdf, URL https://science.sciencemag.org/content/363/6427/615.
  * Mordovina et al. (2020) U. Mordovina, C. Bungey, H. Appel, P. J. Knowles, A. Rubio, and F. R. Manby, Phys. Rev. Research 2, 023262 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.023262.
  * Wang et al. (2021b) D. S. Wang, T. Neuman, J. Flick, and P. Narang, The Journal of Chemical Physics 154, 104109 (2021b).
  * DePrince (2021) A. E. DePrince, The Journal of Chemical Physics 154, 094112 (2021), URL https://doi.org/10.1063/5.0038748.
  * Haugland et al. (2021) T. S. Haugland, C. Schäfer, E. Ronca, A. Rubio, and H. Koch, The Journal of Chemical Physics 154, 094113 (2021), URL https://doi.org/10.1063/5.0039256.
  * Hoffmann et al. (2020) N. M. Hoffmann, L. Lacombe, A. Rubio, and N. T. Maitra, The Journal of Chemical Physics 153, 104103 (2020).
  * Flick and Narang (2020) J. Flick and P. Narang, The Journal of Chemical Physics 153, 094116 (2020), URL https://doi.org/10.1063/5.0021033.
  * Mandal et al. (2020b) A. Mandal, T. D. Krauss, and P. Huo, The Journal of Physical Chemistry B 124, 6321 (2020b), pMID: 32589846.
  * Galego et al. (2017) J. Galego, F. J. Garcia-Vidal, and J. Feist, Phys. Rev. Lett. 119, 136001 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.136001.
  * Flick et al. (2018b) J. Flick, C. Schäfer, M. Ruggenthaler, H. Appel, and A. Rubio, ACS Photonics 5, 992 (2018b).
  * Schafer et al. (2018) C. Schafer, M. Ruggenthaler, and A. Rubio, Phys. Rev. A 98, 043801 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.98.043801.
  * Sidler et al. (2020) D. Sidler, M. Ruggenthaler, H. Appel, and A. Rubio, The Journal of Physical Chemistry Letters 11, 7525 (2020), pMID: 32805122, URL https://doi.org/10.1021/acs.jpclett.0c01556.
  * Theophilou et al. (2020) I. Theophilou, M. Penz, M. Ruggenthaler, and A. Rubio, Journal of Chemical Theory and Computation 16, 6236 (2020), pMID: 32816479, URL https://doi.org/10.1021/acs.jctc.0c00618.
  * Sidler et al. (2021) D. Sidler, C. Schäfer, M. Ruggenthaler, and A. Rubio, The Journal of Physical Chemistry Letters 12, 508 (2021), pMID: 33373238, URL https://doi.org/10.1021/acs.jpclett.0c03436.
  * Buchholz et al. (2020) F. Buchholz, I. Theophilou, K. J. H. Giesbertz, M. Ruggenthaler, and A. Rubio, Journal of Chemical Theory and Computation 16, 5601 (2020), pMID: 32692551, URL https://doi.org/10.1021/acs.jctc.0c00469.
  * Tokatly (2018) I. V. Tokatly, Phys. Rev. B 98, 235123 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.235123.
  * Rivera et al. (2019) N. Rivera, J. Flick, and P. Narang, Phys. Rev. Lett. 122, 193603 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.122.193603.
  * Szidarovszky et al. (2020) T. Szidarovszky, G. J. Halász, and Á. Vibók, New Journal of Physics 22, 053001 (2020), URL https://doi.org/10.1088/1367-2630/ab8264.
  * Cederbaum (2021) L. S. Cederbaum, The Journal of Physical Chemistry Letters 12, 6056 (2021), pMID: 34165990, URL https://doi.org/10.1021/acs.jpclett.1c01570.
  * Cho et al. (2022) D. Cho, B. Gu, and S. Mukamel, Journal of the American Chemical Society 144, 7758 (2022), pMID: 35404593, URL https://doi.org/10.1021/jacs.2c00921.
  * Pavosevic and Rubio (0) F. Pavosevic and A. Rubio, The Journal of Chemical Physics 0, null (0), URL https://doi.org/10.1063/5.0095552.
  * Pavošević and Flick (2021) F. Pavošević and J. Flick, The Journal of Physical Chemistry Letters 12, 9100 (2021), pMID: 34520211, URL https://doi.org/10.1021/acs.jpclett.1c02659.
  * Pavošević et al. (2022) F. Pavošević, S. Hammes-Schiffer, A. Rubio, and J. Flick, Journal of the American Chemical Society 144, 4995 (2022), pMID: 35271261, URL https://doi.org/10.1021/jacs.1c13201.
  * Cederbaum and Kuleff (2021) L. S. Cederbaum and A. I. Kuleff, Nature Communications 12, 4083 (2021), ISSN 2041-1723, URL https://doi.org/10.1038/s41467-021-24221-6.
  * Ebbesen (2016) T. W. Ebbesen, Accounts of Chemical Research 49, 2403 (2016), URL https://doi.org/10.1021/acs.accounts.6b00295.
  * Genet et al. (2021) C. Genet, J. Faist, and T. W. Ebbesen, Physics Today 74, 42 (2021), URL https://doi.org/10.1063/PT.3.4749.
  * Sidler et al. (2022) D. Sidler, M. Ruggenthaler, C. Schäfer, E. Ronca, and A. Rubio, The Journal of Chemical Physics 156, 230901 (2022), URL https://doi.org/10.1063/5.0094956.
  * Li et al. (2022) T. E. Li, B. Cui, J. E. Subotnik, and A. Nitzan, Annual Review of Physical Chemistry 73, 43 (2022), URL https://doi.org/10.1146/annurev-physchem-090519-042621.
  * Sánchez-Barquilla et al. (2022) M. Sánchez-Barquilla, A. I. Fernández-Domínguez, J. Feist, and F. J. García-Vidal, ACS Photonics 9, 1830 (2022), URL https://doi.org/10.1021/acsphotonics.2c00048.
  * Fregoni et al. (2022) J. Fregoni, F. J. Garcia-Vidal, and J. Feist, ACS Photonics 9, 1096 (2022), URL https://doi.org/10.1021/acsphotonics.1c01749.
  * Jaynes and Cummings (1962) E. Jaynes and F. W. Cummings (1962).
  * Rokaj et al. (2018) V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 034005 (2018), URL https://doi.org/10.1088/1361-6455/aa9c99.
  * Ahrens et al. (2021) A. Ahrens, C. Huang, M. Beutel, C. Covington, and K. Varga, Phys. Rev. Lett. 127, 273601 (2021), URL https://link.aps.org/doi/10.1103/PhysRevLett.127.273601.
  * Beutel et al. (2021) M. Beutel, A. Ahrens, C. Huang, Y. Suzuki, and K. Varga, The Journal of Chemical Physics 155, 214103 (2021), URL https://doi.org/10.1063/5.0066427.
  * Liebenthal et al. (2022) M. D. Liebenthal, N. Vu, and A. E. DePrince, The Journal of Chemical Physics 156, 054105 (2022), URL https://doi.org/10.1063/5.0078795.
  * Haugland et al. (2020) T. S. Haugland, E. Ronca, E. F. Kjønstad, A. Rubio, and H. Koch, Phys. Rev. X 10, 041043 (2020), URL https://link.aps.org/doi/10.1103/PhysRevX.10.041043.
  * Mitroy et al. (2013) J. Mitroy, S. Bubin, W. Horiuchi, Y. Suzuki, L. Adamowicz, W. Cencek, K. Szalewicz, J. Komasa, D. Blume, and K. Varga, Rev. Mod. Phys. 85, 693 (2013), URL https://link.aps.org/doi/10.1103/RevModPhys.85.693.
  * Ruggenthaler et al. (2014) M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel, I. V. Tokatly, and A. Rubio, Phys. Rev. A 90, 012508 (2014), URL https://link.aps.org/doi/10.1103/PhysRevA.90.012508.
  * Tokatly (2013) I. V. Tokatly, Phys. Rev. Lett. 110, 233001 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.110.233001.
  * Yang et al. (2021) J. Yang, Q. Ou, Z. Pei, H. Wang, B. Weng, Z. Shuai, K. Mullen, and Y. Shao, The Journal of Chemical Physics 155, 064107 (2021), URL https://doi.org/10.1063/5.0057542.
  * Pellegrini et al. (2015) C. Pellegrini, J. Flick, I. V. Tokatly, H. Appel, and A. Rubio, Phys. Rev. Lett. 115, 093001 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.115.093001.
  * Schafer et al. (2021) C. Schafer, F. Buchholz, M. Penz, M. Ruggenthaler, and A. Rubio, Proceedings of the National Academy of Sciences 118, e2110464118 (2021), eprint https://www.pnas.org/doi/pdf/10.1073/pnas.2110464118, URL https://www.pnas.org/doi/abs/10.1073/pnas.2110464118.
  * Flick (2021) J. Flick, _Simple exchange-correlation energy functionals for strongly coupled light-matter systems based on the fluctuation-dissipation theorem_ (2021), URL https://arxiv.org/abs/2104.06980.
  * Mandel (1979) L. Mandel, Opt. Lett. 4, 205 (1979), URL http://opg.optica.org/ol/abstract.cfm?URI=ol-4-7-205.
  * Runge and Gross (1984) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
  * Troullier and Martins (1991) N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
  * Taylor et al. (2022) M. A. D. Taylor, A. Mandal, and P. Huo, Opt. Lett. 47, 1446 (2022), URL http://opg.optica.org/ol/abstract.cfm?URI=ol-47-6-1446.
  * Chelikowsky et al. (1994) J. R. Chelikowsky, N. Troullier, and Y. Saad, Phys. Rev. Lett. 72, 1240 (1994), URL https://link.aps.org/doi/10.1103/PhysRevLett.72.1240.
  * Yabana and Bertsch (1996) K. Yabana and G. F. Bertsch, Phys. Rev. B 54, 4484 (1996), URL https://link.aps.org/doi/10.1103/PhysRevB.54.4484.
  * Tancogne-Dejean et al. (2020) N. Tancogne-Dejean, M. J. T. Oliveira, X. Andrade, H. Appel, C. H. Borca, G. Le Breton, F. Buchholz, A. Castro, S. Corni, A. A. Correa, et al., The Journal of Chemical Physics 152, 124119 (2020), URL https://doi.org/10.1063/1.5142502.
  * Varga et al. (2004) K. Varga, Z. Zhang, and S. T. Pantelides, Phys. Rev. Lett. 93, 176403 (2004), URL https://link.aps.org/doi/10.1103/PhysRevLett.93.176403.
  * Galego et al. (2019) J. Galego, C. Climent, F. J. Garcia-Vidal, and J. Feist, Phys. Rev. X 9, 021057 (2019), URL https://link.aps.org/doi/10.1103/PhysRevX.9.021057.
  * Benz et al. (2016) F. Benz, M. K. Schmidt, A. Dreismann, R. Chikkaraddy, Y. Zhang, A. Demetriadou, C. Carnegie, H. Ohadi, B. de Nijs, R. Esteban, et al., Science 354, 726 (2016), eprint https://www.science.org/doi/pdf/10.1126/science.aah5243, URL https://www.science.org/doi/abs/10.1126/science.aah5243.
  * Carnegie et al. (2018) C. Carnegie, J. Griffiths, B. de Nijs, C. Readman, R. Chikkaraddy, W. M. Deacon, Y. Zhang, I. Szabó, E. Rosta, J. Aizpurua, et al., The Journal of Physical Chemistry Letters 9, 7146 (2018), URL https://doi.org/10.1021/acs.jpclett.8b03466.
  * Mallory and au2 (2022) J. D. Mallory and A. E. D. I. au2, _Reduced-density-matrix-based ab initio cavity quantum electrodynamics_ (2022), eprint 2204.00725.
  * Becke (2014) A. D. Becke, The Journal of Chemical Physics 140, 18A301 (2014), URL https://doi.org/10.1063/1.4869598.
  * Kockum et al. (2017) A. F. Kockum, V. Macrì, L. Garziano, S. Savasta, and F. Nori, Scientific Reports 7, 5313 (2017), ISSN 2045-2322, URL https://doi.org/10.1038/s41598-017-04225-3.
  * Chervy et al. (2016) T. Chervy, J. Xu, Y. Duan, C. Wang, L. Mager, M. Frerejean, J. A. W. Münninghoff, P. Tinnemans, J. A. Hutchison, C. Genet, et al., Nano Letters 16, 7352 (2016), pMID: 27960510, URL https://doi.org/10.1021/acs.nanolett.6b02567.
  * Barachati et al. (2018) F. Barachati, J. Simon, Y. A. Getmanenko, S. Barlow, S. R. Marder, and S. Kéna-Cohen, ACS Photonics 5, 119 (2018), URL https://doi.org/10.1021/acsphotonics.7b00305.
  * Liu et al. (2019) B. Liu, M. Crescimanno, R. J. Twieg, and K. D. Singer, Advanced Optical Materials 7, 1801682 (2019), eprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/adom.201801682, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/adom.201801682.
  * Welakuh and Narang (2022) D. M. Welakuh and P. Narang, _Tunable nonlinearity and efficient harmonic generation from a strongly coupled light-matter system_ (2022), URL https://arxiv.org/abs/2203.00691.
  * Goncharov and Varga (2012) V. A. Goncharov and K. Varga, The Journal of Chemical Physics 137, 094111 (2012), URL https://doi.org/10.1063/1.4749793.
  * Ullrich (2012) C. A. Ullrich, _Time-Dependent Density-Functional Theory: Concepts and Applications_ (Oxford University Press, USA, 2012).