|  Theory and modeling of light-matter interactions in chemistry: current and
future†
---|---
|  Braden M. Weighta,b, Xinyang Li,a, and Yu Zhang∗‡a
|  Light-matter interaction not only plays an instrumental role in
characterizing the material’s properties via various spectroscopic techniques
but also provides a general strategy to manipulate material properties via the
design of novel nanostructures. This perspective summarizes recent theoretical
advances in modeling light-matter interactions in chemistry, mainly focusing
on plasmon chemistry and polariton chemistry. The former utilizes the highly
localized photon, plasmonic hot electrons, and local heat to drive chemical
reactions. In contrast, polariton chemistry tune chemistry by modifying the
potential energy curvatures via polariton formation. The perspective starts
with the basic background of light-matter interactions, molecular quantum
electrodynamics theory, and the challenges of modeling light-matter
interactions in chemistry. Then, the recent advances in modeling plasmon and
polariton chemistry are described, and future directions toward multiscale
simulations of light-matter interaction-mediated chemistry are discussed.

##

††footnotetext: a Theoretical Division, Los Alamos National Laboratory, Los
Alamos, NM, 87545, USA. Tel: +1 505 606 2149; E-mail:
<EMAIL_ADDRESS>b Department of Physics and Astronomy, University
of Rochester, Rochester. NY, 14627, U.S.A.

##  1 Introduction

In 2007, a report compiled for the Office of Science of the Department of
Energy (DOE) identified five grand challenges in basic energy science 1,
including 1) controlling material processes at the level of electrons, 2)
designing and perfecting atom- and energy-efficient syntheses of new forms of
matter with tailored properties, 3) understanding and controlling the
remarkable properties of matter that emerge from complex correlations of
atomic and electronic constituents, 4) mastering energy and information on the
nanoscale to create new technologies with capabilities rivaling those of
living things, and 5) characterizing and controlling matter away—especially
far away—from equilibrium. One of the emerging techniques to address (some of)
these challenges is light-matter interaction, which can be used to monitor,
manipulate, and design materials’ properties through multiple interactions
between electrons, photons, and phonons 2, 3.

In conventional chemistry, chemists modify the functionalities of molecules
based on different functional groups. For centuries, such a vision has been
widely used in chemical synthesis to modify the chemical and physical
properties of molecules or molecular materials. However, more is different 4,
and chemistry in a complex electromagnetic (EM) environment has emergent
properties due to multiple couplings that can be used to control chemistry.
Light-matter interactions have been instrumental in many branches of physics,
chemistry, materials, and energy science 5, 3, 6. In most of the previous
applications, light-matter interaction is within the weak coupling scenario
and can usually be treated at the lowest order in quantum electrodynamics via
many-body perturbation theory 7. Such treatments are widely used in different
types of spectroscopy techniques 8, 9, optoelectronics 10, 11, 12, 13, 14,
quantum sensing 15, quantum information 16, light harvesting 17, 18, and
beyond.

Alternatively, light-matter interactions open multiple new avenues for
manipulating matter through novel emerging elementary excitations 19,
including plasmons and polaritons. Plasmons (either surface plasmons or
localized surface plasmons) are the collective oscillations of conduction
electrons in nanostructures, which can be excited when the frequency of
external light matches the plasmon resonant energies. Plasmon excitation
results in significantly amplified absorption and scattering cross-sections.
Moreover, plasmon excitations can overcome the diffraction limit and
concentrate the incident light into a highly localized volume, enhancing the
electromagnetic (EM) field (or photons) by several orders of magnitude in the
near field. The subsequent plasmon decay process generates hot electrons or
heat through electron-electron and electron-phonon scatterings at different
length and time scales 20. Nevertheless, the resulting locally enhanced EM
field, hot electrons, and heat can stimulate chemical reactions through
various mechanisms 21, 22. In fact, plasmon-mediated chemical reactions
(namely plasmon chemistry) have become a promising strategy to drive chemical
processes over the past decade 21.

On the other hand, when molecules are collectively and resonantly coupled with
plasmon excitation, a new quasiparticle (namely polaritons) can be formed in
the strong coupling regime. The term "strong" coupling is relative, meaning
that the light-matter coupling is large enough to compete with or overcome
dissipation or dephasing (i.e., when the coherent energy exchange between a
confined light mode and quantum matter is faster than the decay and
decoherence time scales of each part). Such strong coupling can be achieved
when either a plasmonic mode is coupled with a few molecules 23 or many
molecules are collectively coupled to a single cavity mode 24. In the strong
coupling regime, photons and electronic/excitonic excitations in matter become
equally important and are strongly coupled on an equal quantized footing. As a
consequence, individual "free" particles no longer exist. Instead, the
fundamental excitations of the light-matter interacting system are polaritons,
which are hybrid light-matter excitations (superpositions of quantized light
and matter) (Fig 1a)3 and possess both light and matter
characteristics/topologies. Experiments have shown that matter properties can
be modified with the formed polaritons, resulting in different photophysics
and photochemistry25. Since photon energies and light-matter coupling
strengths are relatively tunable through cavity control, light-matter
interaction in the strong coupling regime provides a fundamentally new way to
manipulate matter properties for various desired applications, including
lasing 26, 27, electronics 28, long-range energy transfer 29, 30, 31, 32,
Bose-Einstein condensates 33, 34, 35, and chemical reactions 36, 37, 37, 38,
39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56.
Nevertheless, in polariton chemistry, light and matter cannot be treated as
separate entities, as the strong coupling between them dresses each of them.
Consequently, previous quantum chemistry methods for treating the electronic
structure of matter become invalid. Brand new theoretical and modeling
capabilities for describing polariton chemistry are required.

Fig. 1:  Light-matter interaction is a fundamentally multiscale/multiphysics
problem. It involves multiple interactions among electrons, photons, and
phonons at different time and length scales. As a result, multiple
quasiparticles may be excited due to the light-matter interactions, including
exciton (due to electron-hole interaction), (excitonic/vibrational) polariton
(due to coupling between photon and exciton/phonon), polaron (due to electron-
phonon coupling), and polaron-polaritons. In addition, all quantum systems are
fundamentally open systems, leading to dissipation/dephasing due to the bath.

Despite the attractive applications of strong light-matter interaction, there
are many open questions and fundamental problems about the mechanisms of
physics and chemistry mediated by the strong light-matter interaction. The
experimental progress should be complemented by theoretical advances. Progress
in understanding the coupling between photons and elementary quasiparticles
(plasmons, phonons, and excitons) in materials requires a generalized
treatment of photons as one of the core degrees of freedom (DOFs) in light-
matter interaction. Unfortunately, the strong light-matter interaction is
fundamentally a multiscale and multiphysics problem that involves multiple
interactions between many DOFs and their interplay with environments across
different spatial and time scales. The traditional perturbation methods in the
weak coupling regime are not applicable to the strong coupling regime, making
it urgent to develop new theoretical and modeling techniques, especially
multiscale methods, to understand and ultimately predict strong light-matter
interaction-mediated physics and chemistry.

This perspective reviews recent theoretical efforts toward understanding the
underlying mechanisms of plasmon and polariton chemistry due to the complex
light-matter interactions and discusses our thoughts on future work. The
article is structured as follows: Section 2 introduces the framework and basic
mathematical structure of molecular quantum electrodynamics theory (QED),
Section 3 explores the theoretical development in simulating plasmonic
cavities and their strong interactions with molecules, Section 4 surveys
recent progress in polariton chemistry in Fabry-Pérot-like cavities, and
Section 5 concludes the discussion and provides a final perspective on future
work in all realms of strong light-matter interaction.

##  2 Brief introduction to light-matter interactions and molecular QED theory

This section briefly introduces the quantum theories of light-matter
interactions in the nonrelativistic limit 57. All light-matter interactions
arise from the interplay of matter DOFs (degrees of freedom) (electrons or
spins, nuclei) and an EM environment. Hence, a full ab initio theory for
light-matter interactions fundamentally requires the electronic structure
theory of matter and principles of electrodynamics. Quantum electrodynamics
(QED) is the indispensable and most precise theory for describing the
interactions of charged particles and the dynamics of the EM field in mutual
interaction.

###  2.1 Molecular quantum electrodynamics theory

The minimally coupled Coulomb Hamiltonian governs the nonrelativistic dynamics
of matter in an EM environment,

$\hat{H}=\sum^{N_{e}+N_{n}}_{i}\frac{1}{2m_{i}}\left[\hat{P}_{i}-z_{i}\mathbf{A}(\mathbf{r}_{i}t)/c\right]^{2}+\hat{V}+\hat{H}_{EM},$
(1)

with the Coulomb gauge $\nabla\cdot\mathbf{A}=0$.
$\hat{P}_{i}=-i\hbar\nabla_{i}$ and $z_{i}$ are the momentum operator and
charge of particle $i$, respectively. $N_{e}$ and $N_{n}$ are the number of
electrons and nuclei in the systems.
$\hat{V}=\hat{V}_{ee}+\hat{V}_{nn}+\hat{V}_{en}+V_{ext}$ contains all of the
Coulomb interactions
($\frac{1}{8\pi\epsilon_{0}}\frac{q_{i}q_{j}}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}$)
between the electronic and nuclear DOFs and another external potential. The
Hamiltonian that describes the EM fields is

$\hat{H}_{EM}=\frac{\epsilon_{0}}{2}\int
d\mathbf{r}\left[\hat{\mathbf{E}}^{2}(\mathbf{r})+c^{2}\hat{\mathbf{B}}^{2}(\mathbf{r})\right].$
(2)

It obeys Maxwell’s equations of motion and couples with the Schrödinger
equation self-consistent through the vector potential
$\mathbf{A}(\mathbf{r}_{i})$. With multiple modes available, each field
(including the vector potential $\hat{\mathbf{A}}$) can be expressed as a sum
over all possible radiation modes 58, 59

$\displaystyle\hat{\mathbf{A}}(\mathbf{r},t)=$
$\displaystyle\sum_{\alpha}\mathbf{u}_{\alpha}\Big{(}\hat{A}_{\alpha}e^{i(\mathbf{K}\cdot\mathbf{r}-\omega_{\alpha}t)}+\hat{A}^{\dagger}_{\alpha}e^{i(\mathbf{K}\cdot\mathbf{r}+\omega_{\alpha}t)}\Big{)}$
(3) $\displaystyle\hat{\mathbf{E}}(\mathbf{r},t)=$ $\displaystyle
i\sum_{\alpha}\omega_{\alpha}\mathbf{u}_{\alpha}\Big{(}\hat{A}_{\alpha}e^{i(\mathbf{K}\cdot\mathbf{r}-\omega_{\alpha}t)}+\hat{A}^{\dagger}_{\alpha}e^{i(\mathbf{K}\cdot\mathbf{r}+\omega_{\alpha}t)}\Big{)}$
(4) $\displaystyle\hat{\mathbf{B}}(\mathbf{r},t)=$ $\displaystyle
i\sum_{\alpha}(\nabla\times\mathbf{u}_{\alpha})\Big{(}\hat{A}_{\alpha}e^{i(\mathbf{K}\cdot\mathbf{r}-\omega_{\alpha}t)}+\hat{A}^{\dagger}_{\alpha}e^{i(\mathbf{K}\cdot\mathbf{r}+\omega_{\alpha}t)}\Big{)}$
(5)

where $\nabla\times\mathbf{u}=\mathbf{K}\times\mathbf{u}$,
$\alpha\equiv\lambda\mathbf{K}$ denotes the photon mode with momentum
$\mathbf{K}$ and polarization $\lambda\in\\{-1,1\\}$, $\mathbf{u}$ is the unit
vector denoting the direction of the vector potential, and $\hat{A}_{\alpha}$
($\hat{A}_{\alpha}^{\dagger}$) is the mode decomposition coefficients that
create (annihilate) the $\alpha^{\mathrm{th}}$ radiation mode. Note that
$\hat{\mathbf{E}}=-\frac{1}{c}\partial_{t}\hat{\mathbf{A}}(\mathbf{r})$ and
$\hat{\mathbf{B}}=\frac{1}{c}\nabla\times\hat{\mathbf{A}}(\mathbf{r})$ in the
Coulomb gauge (i.e., $\nabla\cdot\hat{\mathbf{A}}=0$). The photon modes in a
given nanophotonic/nanoplasmonic structure can be readily calculated via
standard mode decomposition techniques 60.

In the second quantization, the photonic Hamiltonian can then be written as,

$\hat{H}_{p}=\sum_{\alpha}\frac{\hbar\omega}{2}\left[\hat{a}_{\alpha}\hat{a}^{\dagger}_{\alpha}+h.c.\right]=\sum_{\alpha}\hbar\omega(\hat{a}^{\dagger}_{\alpha}\hat{a}_{\alpha}+1/2).$

with the equivalence,

$\displaystyle A_{\alpha}\equiv\sqrt{\frac{\hbar}{2\epsilon_{0}\omega
V}}\hat{a}_{\alpha},\ \ \ \ \ \
A^{*}_{\alpha}\equiv\sqrt{\frac{\hbar}{2\epsilon_{0}\omega
V}}\hat{a}^{\dagger}_{\alpha},$ (6)

Thus, the fields are quantized by the isomorphic association of a quantum
mechanical harmonic oscillator to each radiation mode $\alpha$. Introducing
the canonical position and momentum operators as,
$\hat{q}_{\alpha}=\sqrt{\frac{\hbar}{2\omega_{\alpha}}}(\hat{a}^{\dagger}_{\alpha}+\hat{a}_{\alpha})$,
$\hat{p}_{\alpha}=i\sqrt{\frac{\hbar\omega_{\alpha}}{2}}(\hat{a}^{\dagger}_{\alpha}-\hat{a}_{\alpha})$.
The photonic Hamiltonian can now be rewritten as,

$\hat{H}_{p}=\frac{1}{2}\sum_{\alpha}(\hat{p}^{2}_{\alpha}+\omega^{2}\hat{q}^{2}_{\alpha}).$
(7)

Then the full light-matter Hamiltonian in Coulomb gauge can now be fully
constructed and written as,

$\hat{H}_{c}=\sum^{N_{e}+N_{n}}_{i}\frac{1}{2m_{i}}\left(\hat{P}_{i}-z_{i}\hat{\mathbf{A}}(\mathbf{r}_{i})\right)^{2}+\hat{V}+\hat{H}_{p},$
(8)

which is often referred to as the minimal coupling Hamiltonian.

We now aim to apply a unitary transformation on Eq. 8 to achieve an expression
where the momenta of the molecular DOFs ($\hat{p}_{i}$) are decoupled from the
vector potential ($\hat{\mathbf{A}}(\mathbf{r}_{i})$). In other words, we will
shift the light-matter coupling from momentum fluctuations into displacement
fluctuations . This transformation, referred to as the PZW transformation, can
be written as
$\hat{U}_{\mathrm{PZW}}=exp[-\frac{i}{\hbar}\hat{\mathbf{D}}\cdot\hat{\mathbf{A}}]$
where
$\hat{\mathbf{D}}=\sum_{i}^{N_{n}}z_{i}\hat{\mathbf{R}}_{i}-\sum_{i}^{N_{e}}e\hat{\mathbf{r}}_{i}$
is the molecular dipole moment. The PZW transformation
$\hat{U}(\mathbf{d})\hat{H}_{c}\hat{U}^{\dagger}(\mathbf{d})$ is nothing but a
reduction in matter momentum such that
$\hat{p}_{i}-z_{i}\mathbf{A}\rightarrow\hat{p}_{i}$ and a boost in photonic
momentum by
$\hat{p}_{\alpha}\rightarrow\hat{p}_{\alpha}+\sqrt{2\omega/\hbar}\mathbf{d}\cdot\mathbf{A}_{0}$.
Applying PZW transformation results in the QED Hamiltonian, referred to as the
dipole gauge Hamiltonian,

$\hat{H}=\hat{T}+\hat{V}+\hat{H}_{ep}$ (9)

where

$\displaystyle\hat{H}_{ep}=\hat{U}_{\mathrm{PZW}}^{\dagger}\hat{H}_{p}\hat{U}_{\mathrm{PZW}}$
$\displaystyle=\frac{1}{2}\sum_{\alpha}\left[(\hat{p}_{\alpha}+\bm{\lambda}_{\alpha}\cdot\mathbf{D})^{2}+\omega^{2}_{\alpha}(\hat{q}_{\alpha})^{2}\right]$
(10)
$\displaystyle=\frac{1}{2}\sum_{\alpha}\left[\hat{p}^{2}_{\alpha}+\omega^{2}_{\alpha}(\hat{q}_{\alpha}-\frac{\bm{\lambda}_{\alpha}}{\omega_{\alpha}}\cdot\mathbf{D})^{2}\right].$

The second line is reached via canonical transformation between coordinate and
momentum operators, $\hat{p}\rightarrow-\omega\hat{q},\hat{q}\rightarrow
1/\omega\hat{p}$. Here, $\bm{\lambda}_{\alpha}=\sqrt{\frac{1}{\epsilon
V}}\mathbf{u}_{\alpha}$ and
$g_{\alpha}\equiv\bm{\lambda}_{\alpha}\cdot\mathbf{D}$ defines the light-
matter coupling strength, which depends on the volume of quantized photon $V$
in radiation mode $\alpha$. The total light-matter Hamiltonian can now be re-
partitioned after the unitary transformation as,

$\displaystyle\hat{H}_{\mathrm{PF}}$
$\displaystyle=\hat{H}_{\mathrm{M}}+\hat{H}_{\mathrm{p}}+\hat{H}_{\mathrm{ep}}+\hat{H}_{\mathrm{DSE}}$
(11)
$\displaystyle=\hat{H}_{\mathrm{M}}+\sum_{\alpha}\Big{[}\omega_{\alpha}(\hat{a}^{\dagger}_{\alpha}\hat{a}_{\alpha}+\frac{1}{2})+\sqrt{\frac{\omega_{\alpha}}{2}}\bm{\lambda}_{\alpha}\cdot\hat{\mathbf{D}}(\hat{a}^{\dagger}_{\alpha}+\hat{a}_{\alpha})+\frac{1}{2}(\bm{\lambda}_{\alpha}\cdot\hat{\mathbf{D}})^{2}\Big{]}.$
(12)

Here, $\hat{H}_{\mathrm{M}}=\hat{T}_{\mathrm{n}}+\hat{T}_{\mathrm{e}}+\hat{V}$
is the bare molecular Hamiltonian which includes all Coulomb interactions
$\hat{V}$ between electrons and nuclei as well as the kinetic energy operators
of both, $\hat{T}_{\mathrm{e}}$ and $\hat{T}_{\mathrm{n}}$, respectively. This
Hamiltonian is often referred to as the Pauli-Fierz (PF) Hamiltonian.

Now it is clear that the light-matter coupling strength is determined by the
two quantities: the molecular dipole strength and the quantized cavity volume.
Hence, there are two general strategies to enter the strong coupling regime:
(I) reduce the cavity volume and (II) increase the number of molecules
(increasing the total dipole moment). Therefore, there are two major
experimental nanocavity designs for strong light-matter coupling. The first is
the nanophotonic cavity that leverages a large number of molecules to achieve
strong coupling. The other is the nanoplasmonic cavity that leverages a
locally enhanced electric field confined in a small volume to enhance the
coupling. This local field is effectively confined to the nearby surrounding
of a spherical nanoparticle residing on a lattice of nanoparticles on the
scale of nm3 or even Å3 (namely picocavities 61, 23, 62).

The derivaiton of Eq. 10 assumes dipole approximation. However, the dipole
approximation may fail in the ultra-confined nanoplasmonic cavities where the
EM fields are confined within nanometric volumes 63, 64. Consequently, the
size of molecules becomes comparable to the cavity volumes, and the widely
used dipole approximation breaks down. Position-dependent coupling strength
that requires the spatial distribution of excitonic and photonic quantum
states is found to be a key aspect in determining the dynamics in ultrasmall
cavities both in the weak and strong coupling regimes 63.

###  2.2 Weak and strong coupling

The light-matter interactions in nanoplasmonic environments can be split into
weak and strong coupling regimes. The weak-coupling regime is associated with
the Purcell enhancement of spontaneous emission. This effect has been found to
be particularly strong when the molecule is placed next to a metallic surface
or nanostructure 65. In such a regime, the plasmon has been found to be able
to tune photophysics and photochemistry (i.e., plasmon chemistry) via various
possible pathways.

The other regime is the strong coupling regime, where the light-matter
interaction cannot be treated perturbatively. The strong coupling is
characterized by a reversible (coherent) exchange of energy (known as Rabi
oscillation) between the matter and the cavity photon, as the coupling is
strong enough to compete with the dissipation. In this regime, the formation
of polariton requires theoretical methods that treat the matter and photonic
DOFs on equal footing and describe the multiple interactions between photons,
electrons, and nuclei on different length and time scales (polariton
chemistry).

##  3 Plasmon chemistry

Although confined plasmonic modes result in localized and enhanced fields,
plasmons usually suffer from strong dissipation, making it hard to enter the
strong-coupling regime in plasmonic systems. Nevertheless, even without strong
coupling, plasmonics provides a unique setting for manipulating light via the
confinement of EM (below the diffraction limit). Such extreme concentration of
EM field 66 has led to a wide range of applications, such as plasmon-enhanced
molecular spectroscopy 8, 67, 68, 69, 70, photovoltaics 71, 72, nanophotonic
lasers and amplifiers 73, 74, 75, 76, quantum information 16, and many others
77, 18, 78.

Following the ultrafast plasmon excitation, nonradiative plasmon decay leads
to the formation of energetic electron–hole pairs (namely hot-carriers) 79,
80, 81, 82, 83, which are highly nonthermal and can have considerably higher
energies than those rising from thermal equilibrium. The hot electrons (HEs)
(and their corresponding holes) redistribute their energies quickly as a
result of electron–electron scattering 79, 80, reaching a quasi-thermal
equilibrium but with a high effective temperature. Further cooling of the hot
electrons takes place via energy dissipation into the phonon modes of the
nanoparticle, and the energy is ultimately dissipated to the surroundings via
thermal conduction. Nevertheless, investigations in the last two decades have
found that chemical reactions can be stimulated by localized EM fields (or
photons), and electronic and/or thermal energies (that result from plasmon
decay) via various different pathways (more details in Sec. 3.3). This leads
to an emerging field of plasmon chemistry that designs nanostructure-based
surface plasmons as mediators to redistribute and convert photon energy in
various time, space, and energy scales to drive chemical reactions 70, 21, 84,
85, 20, 86, 87, 88, 89, 22, 90, 91, 92, 93, 94, 95.

In most scenarios, the nanoplasmonic properties are obtained by solving
Maxwell’s equations with proper dielectric function and boundary conditions.
However, Maxwell’s equations fail to describe the quantum effects (such as the
tunneling in charge transfer mode and quantum confinement 96), and quantum
theory for plasmon (quantum plasmonics) is required in nanoscale.

###  3.1 Ab initio method for quantum plamsonics

From a quantum mechanical point of view, a plasmon is nothing but a specific
collective excitation. Hence, the quantum chemistry method for excited states
can be used to predict the optical properties of plasmonic nanostructures 97.
The many-body perturbation theory, time-dependent density functional theory
(TDDFT), or time-dependent density functional tight-binding (TDDFTB) are
methods of choice to compute the plasmon excitations 98, 99, 100, 101, 102,
103, 104, 105, 106. Without losing generality, here we use TDDFT theory as an
example to demonstrate the computation of plasmon excitations.

TDDFT is the formal extension of the Hohenberg-Kohn-Sham density functional
theory (DFT). Within the DFT, The Hohenberg-Kohn (HK) theorem 107 states that
the ground-state electron density unambiguously defines the many-electron
ground state for an $N-$electron system under the influence of an external
potential. And there exists an energy functional that guarantees the ground
state energy can be reached by the variational principle though the exact
functional is an unsolved problem,

$E[\rho]=T+V_{H}(\mathbf{r})+E_{xc}[\rho]+\int
V(\mathbf{r})\rho(\mathbf{r})d\mathbf{r}.$ (13)

Where
$V_{H}[\rho]=\int\frac{\rho(\mathbf{r}^{\prime})}{|\mathbf{r}-\mathbf{r}^{\prime}|}d\mathbf{r}^{\prime}$
is the Hartree potential, $E_{xc}$ is the exchange-correlation functional,
which contributes $<10\%$ to the total energy, but describes the most critical
correlation effects 108. But DFT is a ground-state theory. Light-matter
interactions usually result in many elementary excitations within the matter.
Among various many-body methods, TDDFT (the formal extension of DFT theory) is
the preferred tool to evaluate the excited states and optical properties in
extended molecular or condensed matter systems.

In the linear response regime, the induced electron density due to the light-
matter interaction $\delta V(\mathbf{r})$ is given by

$\delta\rho(\mathbf{r})=\int\chi(\mathbf{r},\mathbf{r}^{\prime})\delta
V_{ext}(\mathbf{r}^{\prime})d\mathbf{r}^{\prime}.$ (14)

Where $\chi(\mathbf{r},\mathbf{r}^{\prime})$ is the polarizability describing
the density response of the many-body ground state with respect to an external
perturbation. In addition, according to Runge-Gross theorem 109,
$\delta\rho(\mathbf{r})$ is also the induced density of the KS system, but due
to a perturbation

$\delta\rho(\mathbf{r})=\int\chi_{KS}(\mathbf{r},\mathbf{r}^{\prime})\delta
V_{KS}(\mathbf{r}^{\prime})d\mathbf{r}^{\prime}.$ (15)

$\chi_{KS}$ here is the linear response of the KS electrons, which can be
trivially evaluated from the KS orbitals and energies.

$\delta V_{KS}(\mathbf{r})=\delta
V_{ext}(\mathbf{r})+V_{H}[\delta\rho(\mathbf{r})]+\delta V_{xc}(\mathbf{r}).$
(16)

The induced XC potential is given by $\delta V_{xc}(\mathbf{r})=\int
f_{xc}(\mathbf{r},\mathbf{r}^{\prime})\delta\rho(\mathbf{r}^{\prime})d\mathbf{r}^{\prime}$,
where $f_{xc}$ is the dynamical XC kernel $\frac{\delta
V_{xc}(\mathbf{r})}{\delta n(\mathbf{r}^{\prime})}$ 110,

$f_{xc}(\mathbf{r},\mathbf{r}^{\prime})=\frac{\delta
V_{xc}(\mathbf{r})}{\delta\rho(\mathbf{r}^{\prime})},$ (17)

Hence, a linear equation for the induced density can be derived from Equations
14-17,

$[1-\chi_{KS}(\omega)f_{hxc}(\omega)]\chi(\omega)\delta
V_{ext}(\omega)=\chi_{KS}(\omega)\delta V_{ext}(\omega).$ (18)

where
$f_{hxc}(\mathbf{r},\mathbf{r}^{\prime})=\frac{1}{|\mathbf{r}-\mathbf{r}^{\prime}|}+f_{XC}(\mathbf{r},\mathbf{r}^{\prime})$.
Casting Equation 18 into the matrix of coupled KS single excitations, it is
possible to calculate the excitation energies $\omega_{s}$ of the system (the
poles of the response function) and transition densities (and oscillator
strengths). Consequently, the response function $\chi(\omega)$ can be
rewritten as

$\chi(\mathbf{r},\mathbf{r}^{\prime},\omega)=2\sum_{s}\rho_{s}(\mathbf{r})\rho_{s}(\mathbf{r}^{\prime})\zeta_{s}(\omega),$
(19)

where
$\zeta_{s}(\omega)=\frac{1}{\omega-\omega_{s}+i\eta}-\frac{1}{\omega+\omega_{s}-i\eta}$
and $\eta$ represents a positive inifinitesimal. The factor of 2 comes from
summation over spin indices. The excitation energies in the system are denoted
as $\omega_{s}$, which are calculated from the Casida method 110.
$\rho_{s}(\mathbf{r})$ is the transition density, which can be expanded in
terms of electronic transitions between occupied state $i$ to unoccupied
states $a$ 111,

$\rho_{s}(\mathbf{r})=\sum_{ia}X^{s}_{ia}\psi_{i}(\mathbf{r})\psi_{a}(\mathbf{r})\left(\frac{\epsilon_{a}-\epsilon_{i}}{\omega_{s}}\right)^{1/2},$
(20)

where $\epsilon_{i}$ are the KS eigenvalues and the corresponding molecular
orbitals are $\psi_{i}$. And $X^{s}_{ia}$ are the Casida transition
coefficient from the occupied $i$ state to the unoccupied $a$ state of the
$s^{th}$ excitation. By examining the nature of $X^{s}_{ia}$, it is possible
to distinguish between normal excitation and plasmonic excitation. Plasmonic
excitations generally are characterized by collective transitions from
occupied to virtual orbitals 112, 113.

###  3.2 Hot electron generation and relaxation

Once the plasmon excitation is obtained, the electron-plasmon interaction can
be described by the Hamiltonian

$\hat{H}_{int}=\frac{e}{2m_{e}}\int
d\mathbf{r}\hat{\Psi}^{\dagger}V_{\textit{eff}}(\mathbf{r})\hat{\Psi},$ (21)

where $V_{\textit{eff}}(\mathbf{r})$ is the effective potential induced by the
excitation, which can be calculated from the transition density
$\rho_{s}(\mathbf{r})$ by following Lundquvist’s approach 112, 111. The
polarizability of the NP upon external excitation is given by Equation 19. And
the induced potential is given by 111,

$V_{\textit{eff}}(\mathbf{r})=\epsilon^{-1}(\mathbf{r},\mathbf{r}^{\prime})\delta
V_{ext}(\mathbf{r}^{\prime}).$ (22)

Where the dielectric function $\epsilon^{-1}(\mathbf{r},\mathbf{r}^{\prime})$
is given by

$\epsilon^{-1}(\mathbf{r},\mathbf{r}^{\prime})=\delta(\mathbf{r},\mathbf{r}^{\prime})+\int
d\mathbf{r}_{1}f_{hxc}(\mathbf{r},\mathbf{r}_{1})\chi(\mathbf{r}_{1},\mathbf{r}^{\prime},\omega).$
(23)

Thus, within the second quantization, electron-excitation Hamiltonian
$\hat{H}_{eps}$ reads

$\hat{H}_{eps}=\sum_{ij}\left[\mathcal{M}^{s}_{ij}\hat{c}^{\dagger}_{i}\hat{c}_{j}\hat{b}_{s}+\text{h.c.}\right],$
(24)

where $\mathcal{M}^{s}_{ij}$ is the electron-excitation coupling strength,
describing the scattering of quasiparticles from state $i$ into state $j$ via
the emission or absorption of excitation in the state $s$. These elements are
given by (see SI for details)

$\mathcal{M}^{s}_{ij}=\zeta_{s}(\omega)\mathbf{E}\cdot\mathbf{f}^{s}\langle\psi_{i}(\mathbf{r}_{1})|V^{s}_{H}(\mathbf{r}_{1})|\psi_{j}(\mathbf{r}_{1})\rangle$
(25)

where
$\mathbf{f}^{s}=\sqrt{\frac{2m_{e}\omega_{s}}{2\hbar^{2}}}\sum_{ia}X^{s}_{ia}\left(\frac{\epsilon_{a}-\epsilon_{i}}{\omega_{s}}\right)^{1/2}\mathbf{d}_{ia}$
is the oscillator vector and
$|\mathbf{f}^{s}|^{2}=\frac{2m_{e}}{3\hbar^{2}}\omega_{s}|\langle\Psi_{0}|\mathbf{r}|\Psi_{s}\rangle|^{2}$
is the oscillator strength. $V^{s}_{H}(\mathbf{r})=\int
d\mathbf{r}^{\prime}\frac{\rho_{s}(\mathbf{r}^{\prime})}{|\mathbf{r}-\mathbf{r}^{\prime}|}$
is the Hartree potential induced by the transition density. Detailed
derivation can be found in the SI. Above equation actually includes both
plasmon excitation ($\zeta_{s}(\omega)\mathbf{E}\cdot\mathbf{f}^{s}$),
screening effect and hot electron-hole pairs generation
($\langle\psi_{i}(\mathbf{r}_{1})|V^{s}_{H}(\mathbf{r}_{1})|\psi_{j}(\mathbf{r}_{1})\rangle$).
Neglecting the screening effect would result in a much larger electron-plasmon
coupling matrix and shorter HC lifetime distribution 114. In contrast, our
model connects the widely used semiempirical model and the recently developed
quantum model with the plasmon excitation and HC generation treated on equal
footing. Equation 24 is the general formalism that describes the coupling
between electrons and photoexcitation. If the external field matches the
plasmon energy, plasmon resonance will be excited, and Equation 24 reduces to
the electron-plasmon coupling.

HC generation. The electron-plasmon coupling describes the HC generation
following plasmon decay. After the electron-plasmon coupling matrix is
obtained, the HC generation can be readily calculated from the Fermi golden
rule 83, 79,

$\Gamma^{ex}_{i\rightarrow
a}=\frac{4}{\hbar}\sum_{s}|\mathcal{M}^{s}_{ia}|^{2}\frac{\gamma_{ex}}{(\epsilon_{i}-\epsilon_{a}+\omega_{s})^{2}+\gamma^{2}_{ex}}\delta(\omega-\omega_{s}).$
(26)

$\gamma_{ex}$ refers to the linewidth of the excitation. $\omega$ is the
excitation energy. $\delta(\omega-\omega_{s})$ describes the generalized
photoexcitation. Hence, Equation 26 describes the photoexcitation and HC
generation. When the energy of the external field $\omega$ matches the plasmon
energy, Equation 26 describes the HC generation from plasmon decay. Otherwise,
it describes the HC generation from regular excitation.

HC relaxation. After generation from plasmon decay, the HCs will undergo a
relaxation process mainly due to the electron-electron and electron-phonon
scatterings. The dynamics of the HCs are investigated by propagating the
density matrix $\rho(t)$. The equation of motion (EOM) of $\rho$ is subject to
the following quantum Liouville Von-Neumann equation 115,

$i\partial_{t}\rho(t)=i\partial_{t}\rho(t)|_{\text{coh}}+\mathcal{L}_{ee}[t]+\mathcal{L}_{ph}[t]+\mathcal{L}_{s}[t]$
(27)

The first part on the right-hand side (RHS) of the above equation describes
the coherent evolution of the density matrix between different states. The
second and third parts on the RHS of Equation 27 represent the dissipations
induced by the electron-electron and electron-phonon scattering effects,
respectively. Finally, the last term on the RHS of Equation 27 described the
extraction of HCs. In general, in the presence of electron-electron and
electron-phonon interaction, the dissipation can be described by a Liouville
equation by employing the many-body perturbation theory (MBPT) 116,
$\mathcal{L}(t)=\mathcal{Q}(t)-\mathcal{Q}^{\dagger}(t)$. Where the
dissipation matrix $\mathcal{Q}_{t}$ can be written in terms of Green’s
functions $G$ and self-energies $\Sigma$ 116, 79, i.e.,

$\mathcal{Q}(t)=i\int^{t}dt^{\prime}_{-\infty}\left[\Sigma^{<}(t,t^{\prime})G^{>}(t^{\prime},t)-\Sigma^{>}(t,t)G^{<}(t^{\prime},t)\right].$
(28)

Hence, the key point is to develop approximations to the self-energies
$\Sigma^{<,>}(t,t^{\prime})$ and efficient numerical methods to compute the
$\mathcal{Q}(t)$. The detailed formalisms of $\mathcal{Q}(t)$ for electron-
electron and electron-phonon scatterings can be found in Ref. 79. But it
should be noted that the approximation used for electron-electron scattering
should conserve both particles and energy.

###  3.3 Plasmon-mediated chemical reactivities and theoretical challenges

Fig. 2: (a) Proposed mechanisms of plasmon-mediated reactivities. 1) Enhanced
intramolecular excitation, 2) direct charge transfer, 3) indirect charge
transfer, and 4) local heating.

However, the computations of plasmon excitation and distribution of hot
electrons (from plasmon decay) are insufficient to give insights into plasmon-
mediated chemistry because plasmon-induced photocatalysis is a complex
dynamical process that involves multiple interactions and mechanisms. To date,
four major microscopic mechanisms have been proposed to explain how plasmonics
facilitates various chemical reactions through the concentration of light and
hot-carrier dynamics (Figure. 2). The first mechanism involves the extreme
concentration of light, which significantly promotes the excitation of
electrons within the adsorbed molecule. Such an excitation takes place within
the adsorbate only but is significantly enhanced by plasmonics, which is
referred to as the enhanced intramolecular excitation mechanism 92, 90, 117.
The second mechanisms involve the transfer of plasmonic hot carriers. Due to
the hybridization of the adsorbate–nanoparticle system, the generated hot
electrons (or holes) can transfer to the absorbed molecules. Such transfers
result in the placement of an adsorbate–nanoparticle system onto a manifold of
excited potential energy surfaces (PESs), where the adsorbed molecule
experiences strong forces that activate its internal vibrational excitations
and, ultimately, chemical transformation. This mechanism is attributed to an
indirect HE transfer from metal nanoparticles to the adsorbed molecule 91, 22.
Third, due to the hybridization, direct excitation of electrons from metal
states near the Fermi level to the unoccupied molecular orbital (LUMO) of the
adsorbed molecule can occur when the plasmon frequency is resonant with the
excitation energy between the metallic occupied orbitals and hydrated metal-
molecular orbitals. Such a reaction pathway is referred to as the direct
charge transfer (CT) mechanism. Such excitation circumvents the thermalization
of HEs 94, 118, 119, 120, but requires matching between the plasmon energy and
the energy gap between the metal states and the molecular LUMO. Finally, local
heating resulting from hot-carrier relaxation can thermally activate a
reaction 121, 93, 122.

These mechanisms share some common features that render their clear
distinction very challenging. Indeed, different mechanisms have been proposed
even for the same chemical reaction, leaving a very confusing situation.
Hence, in order to obtain a comprehensive understanding of plasmon-mediated
chemistry, it’s essential for the theoretical methods to fulfill the following
requirements:

  1. 1.

Efficient computation of a dense manifold of excited states because plasmonic
excitations in metallic nanostructures are usually not low-lying states. Even
for a small nanoparticle (2 nm diameter, for example), thousands of excited
states are to be computed to reach the plasmon excitations 87.

  2. 2.

The ability to capture the nonadiabatic transition between excited states
since the indirect transfer mechanism involves the transition between the hot
electron state (excitations localized within the plasmonic system only) and
the charge transfer state 89.

  3. 3.

Inclusion of electron-vibrational couplings that lead to the hot electron
relaxation and local heating.

  4. 4.

Inclusion of electron-electron scattering that leads to the redistribution of
hot electrons 123, 82.

In particular, recent debates on the thermal impact underpin the necessity of
atomistic and dynamical insights via nonadiabatic simulations of the HE
generation, transfer, and relaxation processes on equal footing.

###  3.4 Nonadiabatic simulation of Plasmon-mediated chemical reactivities

Ehrenfest dynamics and Surface Hopping. With the Born-Oppenhermier
approximation, the electronic wave function $\Theta(\mathbf{r},\mathbf{R})$ is
expanded in the basis of adiabatic BO states, which depend on the electronic
coordinates $\mathbf{r}$ and the nuclear coordinates $\mathbf{R}(t)$ according
to

$\Theta(\textbf{r},\textbf{R})=\sum_{n=1}^{N_{st}}c_{n}(t)\ket{\phi_{n}(\textbf{r},\textbf{R}(t))}.$
(29)

Here $N_{st}$ is the total number of adiabatic electronic states,
$\ket{\phi_{n}(\textbf{r},\textbf{R}(t))}$ is the adiabatic electronic
wavefunction of state n and $c_{n}(t)$ are the time-dependent complex
expansion coefficients.
$\textbf{R}(t)=\left\\{\textbf{R}_{A}(t)\right\\}_{A=1}^{A=N_{A}}$ ($N_{A}$ is
the total number of atoms in the system) are the nuclear trajectories which
are obtained by solving the classical Newton’s equations of motion (EOMs) with
the mixed quantum-classical framework 124,

$M_{A}\frac{d^{2}\textbf{R}_{A}}{dt^{2}}=-\nabla_{\textbf{R}_{A}}E(\textbf{R}),$
(30)

where $M_{A}$ is the mass of $A^{th}$ atom. $E(\textbf{R})$ is the potential
energy surface (PES), which can be an averaged one
$E(\mathbf{R})=\sum_{n}C_{n}(t)E_{n}(\mathbf{R})$ within the mean-field
Ehrenfest dynamics or a single PES of a certain state $E_{k}(\mathbf{R})$
within the Surface hopping (SH) framework 125, 126, 124.

By substituting Eq. 29 into the time-dependent Schrödinger equation and
keeping only the first-order nonadiabatic coupling terms, a set of EOMs for
the coefficients $c_{n}(t)$ along a given classical trajectory can be obtained
124, 127

$i\hbar\frac{\partial c_{n}(t)}{\partial
t}=c_{n}(t)E_{n}(\textbf{R})-i\hbar\sum_{m}c_{m}(t)\dot{\textbf{R}}\cdot\textbf{d}_{nm}.$
(31)

Here the orthogonal condition of adiabatic states
$\langle\Psi_{n}|\Psi_{m}\rangle=\delta_{nm}$ is used.
$\mathbf{d}_{nm}=\langle\Psi_{n}|\nabla_{\textbf{R}}|\Psi_{m}\rangle$ is the
nonadiabatic derivative coupling term (or nonadiabatic coupling vector, NACR).
A key variable in Eq. 31 is the time-derivative nonadiabatic coupling scalar
(NACT) between two adiabatic states

$\dot{\textbf{R}}\cdot\textbf{d}_{nm}=\langle\Psi_{n}|\frac{\partial}{\partial
t}|\Psi_{m}\rangle,$ (32)

which is responsible for the nonadiabatic transitions between different
adiabatic states and can be easily calculated with many ab initio methods for
excited states.

With the SH algorithm, the state that governs the dynamics of nuclei is
determined by a stochastic process. The time-dependent coefficients $c_{n}(t)$
obtained in Eq. 31 are used to calculate the hopping probabilities between
different electronic excited states within the framework of the FSSH
algorithm. The hopping probabilities between excited states $n$ and $m$ are
given by128

$g_{n\rightarrow m}(\textbf{R},t)=\frac{\int_{t}^{t+N_{q}\delta t}dt\
b_{mn}(\textbf{R},t)}{a_{nn}(t)},$ (33)

where $N_{q}=\frac{\Delta t}{\delta t}$, with $\Delta t$ and $\delta t$
correspond to the time steps for evolving motions of nuclei and electrons in
Eq. 30 and Eq. 31, respectively. The chosen value of $\delta t$ must be small
enough to resolve strongly localized peaks in NACT in order to avoid
underestimation of transition probabilities. This is particularly important
when crossings between adiabatic states are encountered in the trajectory, in
which $\delta t$ will be further refined. $a_{nn}(t)=c_{n}(t)c_{n}^{*}(t)$
defines the time-dependent density matrix elements, and
$b_{mn}(\textbf{R},t)=-2Re(a_{nm}^{*}\ \dot{\textbf{R}}\cdot\textbf{d}_{nm})$.
Note that $g_{n\rightarrow m}=-g_{m\rightarrow n}$ and $g_{n\rightarrow n}=0$
since $\textbf{d}_{nm}$ are antisymmetric. Hopping between adiabatic states is
determined stochastically by comparing $g_{n\rightarrow m}$ to a random number
$\xi(\xi\in(0,1))$. A hop from state $n$ to state $m$ is performed if

$\sum_{l=1}^{m-1}g_{n\rightarrow l}<\xi\leq\sum_{l=1}^{m}g_{n\rightarrow l},$
(34)

where states are assumed to be ordered with increasing transition energy. On
the other hand, the system remains in state n when
$\sum_{l=1}^{N_{st}}g_{n\rightarrow l}<\xi<1$. If $g_{n\rightarrow m}<0$, the
hop is unphysical, and the probability is set to zero. Finally, if a hop to a
higher energy state is predicted, there must be sufficient nuclear kinetic
energy along the direction of NACR. Otherwise, the hop is rejected. After a
successful hop, the total electron-nuclear energy is conserved by rescaling
the nuclear velocity in the direction of the NACR according to the procedure
described in reference 126, 129. In addition, during the dynamics, we monitor
the relative phase of the ground to excited state transitions and maintain the
same phase (sign) to avoid a sudden sign change in the NACT. This is done by
enforcing the sign of the largest component of the Casida eigenvectors to the
same along the trajectory.

Despite broad popularity in the community, either Ehrenfest dynamics or the
surface hopping (SH) approach have well-known limitations such as generating
artificial electronic coherence or giving incorrect long-time population 126.
To address these challenges and provide accurate dynamics, a
multiconfigurational Ehrenfest (MCE) dynamics approach 130 and ab initio
multiple cloning (AIMC) 131, 132, 133, 134 are developed accordingly.

Multiconfigurational Ehrenfest (MCE). MCE generalizes EHR formalism by
representing the wave function as a linear combination of Ehrenfest
configurations. Each MCE configuration moves along its own Ehrenfest (mean-
field) trajectory. Within the MCE formalism, the molecular wavefunction
$\ket{\Psi}$ is expressed in the trajectory-guided Gaussian basis functions
(TBF) representation ($\ket{\psi_{n}}$),

$\ket{\Psi(t)}=\sum c_{n}\ket{\psi_{n}(t)}.$ (35)

And each configuration (or TBF) is described by the product of nuclear and
electronic parts,

$\ket{\psi_{n}(t)}=\ket{\chi_{n}(\mathbf{R},t)}\sum_{I}a^{n}_{I}\ket{\phi^{n}_{I}(\mathbf{r},\mathbf{R}(t))}.$
(36)

Where $\ket{\phi^{n}_{I}(\mathbf{r},\mathbf{R}(t))}$ is the adiabatic state of
configuraiton $n$. $\ket{\chi_{n}}$ are Gaussian nuclear basis functions

$\ket{\chi_{n}}=\left(\frac{2\alpha}{\pi}\right)^{N_{d}/4}e^{\left\\{-\alpha(R-\bar{R})+\frac{i}{\hbar}P(R-\bar{R})+\frac{i}{\hbar}\gamma_{n}(t)\right\\}}.$
(37)

The couplings between TBFs in the MCE approach are described by the EOM of
$c_{n}(t)$, which can be readily obtained by substituting Eq. 35 into the Schr
"odinger equation:

$i\hbar\sum_{n}S_{mn}\dot{c}_{n}=\sum_{n}\left[H_{mn}-i\hbar\bra{\psi_{m}}\frac{d\psi_{m}}{dt}\rangle\right]c_{n}.$
(38)

where

$H_{mn}=\sum_{I,J}(a^{m}_{I})^{*}a^{n}_{J}\bra{\chi_{m}\phi^{m}_{I}}T+V\ket{\chi_{n}\phi^{n}_{J}},$
(39)

and the overlap $S_{mn}$ is

$S_{mn}=\bra{\psi_{m}}\psi_{n}\rangle=\langle{\chi_{m}}\ket{\chi_{n}}\sum_{I,J}(a^{m}_{I})^{*}a^{n}_{J}\langle{\phi^{m}_{I}}\ket{\phi^{n}_{J}}.$
(40)

The nuclear part of Eq. 39 can be obtained analytically,

$\bra{\chi_{m}\phi^{m}_{I}}V\ket{\chi_{n}\phi^{n}_{J}}=\bra{\chi_{m}}-\frac{\hbar^{2}}{2}\nabla_{\mathbf{R}}M^{-1}\nabla_{\mathbf{R}}\ket{\chi_{n}}\langle\phi^{m}_{I}\ket{\phi^{n}_{J}}.$
(41)

While the electronic part (or the potential energy matrix elements) are
approximated by 134

$\displaystyle\bra{\chi_{m}\phi^{m}_{I}}V\ket{\chi_{n}\phi^{n}_{J}}=\frac{1}{2}\bra{\phi^{m}_{I}}\phi^{n}_{J}\rangle\bra{\chi_{m}}\chi_{n}\rangle\times$
$\displaystyle\Big{\\{}(V^{m}_{I}+V^{n}_{J})+\frac{i}{4\alpha\hbar}(\mathbf{P}_{n}-\hat{P}_{m})\cdot(\nabla_{\mathbf{R}}V^{m}_{I}+\nabla_{\mathbf{R}}V^{n}_{J})$
$\displaystyle-\frac{1}{2}(\mathbf{R}_{m}-\mathbf{R}_{n})\cdot(\nabla_{\mathbf{R}}V^{m}_{I}-\nabla_{\mathbf{R}}V^{n}_{J})\Big{\\}}.$
(42)

As shown above, within the MCE schemes, the electronic states are now
different for different configurations. The overlaps between TBFs have to be
calculated and taken into account.

Ab initio multiple cloning (AIMC). The AIMC method combines the best features
of ab initio Multiple Spawning (AIMS) 135 and Multiconfigurational Ehrenfest
(MCE) methods. Similar to the MCE method, the individual trajectory basis
functions (TBFs) of AIMC follow Ehrenfest equations of motion. However, the
basis set is expanded in a similar manner to AIMS when these TBFs become
sufficiently mixed. Consequently, AIMC avoids prolonged evolution on the mean-
field potential energy surface (PES).

Within the MCE formalism, the Ehrenfest basis set is guided by an average
potential, which is accurate for dynamical processes where the coupling
between states persists in time between nearly parallel PES. But the mean-
field treatment can be unphysical when the PES of two or more populated
electronic states become different in shape, which leads to wave packet
branching after leaving the nonadiabatic coupling region. To deal with these
cases, the AIMC algorithm is applied to expand the original basis set of TBFs
by "cloning" one TBF into two copies in a way that does not alter the original
wave function 133, 131, 132 This is done by creating one of the clones
$\ket{\psi_{n_{1}}}$ in a pure state and the other clone $\ket{\psi_{n_{2}}}$,
which includes contributions from all other electronic states. The
corresponding MCE amplitudes $\\{c_{n_{1}},c_{n_{2}}\\}$ are adjusted to
conserve the original wavefnction 134, 133.

It should be noted the MQC method can be complemented with any electronic
structure solvers as long as the gradients and NACs are available. Practical
applications need to balance accuracy and numerical efficiency.

Fig. 3: (a) Schematic diagram of H2 molecule adsorbed on the NP surface. b)-c)
Spatial distribution of the bonding and antibonding orbitals of H2 on Jellium
NP. d) Local DOS of bonding and antibonding orbitals of H2 modules as a
function of H-H bond length. e) Depending on the symmetry, the H2 dissociation
can be suppressed or restored in the plasmonic dimer. Figures are adapted with
permission from Ref. 88. Copyright 2018 American Chemical Society.

###  3.5 Simulation of plasmon-mediated phenomena

####  3.5.1 Jellium model-based NAMD simulation

Compared to the costly atomistic ab initio calculations for plasmonic
nanoparticles, the optical absorption of simple sp-metal nanoparticles
described by the Jellium model can be easily calculated and analyzed 136, 137,
97, 138, 83, 88. In general, good agreement with experimental optical
properties can be achieved by tuning the Jellium radius. Using the Jellium
model, we investigated the generation and relaxation of plasmonic hot carriers
79, as well as the atomic-scale mechanism of plasmonic hot-carrier-mediated
chemical processes such as H2 dissociation 88. Our numerical simulations
showed that after photoexcitation, hot carriers transfer to the antibonding
state of the H2 molecule from the nanoparticle, leading to a repulsive-
potential-energy surface and H2 dissociation (Figure LABEL:fig:_Jellium(b-d)).
This process occurs when the molecule is close to a single nanoparticle.
However, in a plasmonic dimer, dissociation can be inhibited due to sequential
charge transfer that effectively reduces the occupation of the antibonding
state, as shown in Figure 3(e). When the molecule is asymmetrically positioned
in the gap, the symmetry is broken, and dissociation is restored by
significantly suppressing additional charge transfer. Thus, these models
illustrate the potential for structurally adjustable photochemistry through
plasmonic hot carriers.

####  3.5.2 TDDFT calculations.

However, the Jellium model oversimplifies the electronic structures of the
plasmonic nanostructures, and it lacks atomistic details. Insights from PESs
built on the ab initio atomistic model are essential for a more in-depth
understanding of plasmon chemistry. Thus, we employed the linear response
TDDFT (LR-TDDFT) calculations within the Casida formalism to compute the
adiabatic PESs of the H2 molecule adsorbed on an Au6 cluster (H2@Au6) 89, in
order to explore key pathways in LSPR-promoted chemical reactions. Despite the
model system is too small to support plasmonic mode and thus cannot describe
the dephasing of plasmons that produce hot electrons, the key point of using
the model system is to capture key aspects of the later stages of plasmon-
facilitated photocatalysis, thus providing mechanistic insights.

Our findings based on DFT calculations indicate that in the ground state, the
Au6 cluster supports the adsorption of the H2 molecule at the tip site. We
further use LR-TDDFT calculations within the Casida formalism to determine the
adiabatic excited states and corresponding oscillator strengths. We calculate
two-dimensional PESs in the desorption and dissociation reaction coordinates
and observe that the adiabatic excited PESs bear similarity to the ground
state PES in the Franck-Condon region. Therefore, the H2 adsorbate is
comparatively stable, as corroborated by the relatively low dissociation and
desorption probabilities from our quantum dynamics simulations. By developing
the orbital wave function overlapping (OWO) diabaization scheme, we are able
to divide the dense manifold of the excited state into two groups: 1) one
group is dominated by electronic excitations confined to the Au6 cluster,
which can be likened to HE states in metal nanoclusters; 2) the other group of
excited states has the antibonding $\sigma^{*}$ characteristics of the H2
adsorbate due to the hybridization between H2 and Au6 cluster, which are
denoted as CT states. The crossings among the HE and CT states provide
pathways leading from the excited HE states to CT states via nonadiabatic
transitions. Quantum dynamics simulations on the diabatic PESs demonstrate
that the CT diabatic states are able to drive H2 dissociation efficiently and
thus are responsible for the experimentally observed HD formation on Au
nanoparticles. Our results nevertheless give a clear physical picture of
photoinduced H2 dissociation on Au clusters.

Fig. 4:  Schematic diagram of competition among different pathways in plasmon-
mediated chemical reactions on a dense manifold of excited states. The excited
states can be divided into HE and CT states. The CT states are responsible for
the chemical reaction, which can be triggered by the nonadiabatic transition
between HE and CT states. Such nonadiabatic transitions have to compete with
hot carrier relaxations that lead to local heating.

The presence of HE-CT crossings is not unique to plasmonic catalysis. Such
features have also been found in other nonplasmonic catalysis, such as those
on the surface of semiconductors 139. However, plasmonic materials are unique
in generating a high concentration of HE due to their larger absorption cross
sections. The combination of a high concentration of HE and HE-CT crossings is
thus the distinctive feature of plasmonic catalysis.

####  3.5.3 TDDFTB-based NAMD simulations

As shown in Figure 4, plasmon-mediated chemistry involves a dense manifold of
excited states, nonadiabatic transitions between HE and CT states, and their
competition with various dissipation channels. Such complicated processes
require nonadiabatic simulations that treat the HE generation, relaxation, and
HE-CT transitions treated on equal footing. To this end, larger clusters
should be used in modeling plasmonic catalysis, and more efficient
semiempirical methods, such as those based on time-dependent density
functional tight binding (TDDFTB) 140, might be needed for the efficient
electronic structure calculations. To this aim, We recently developed an
efficient NAMD method by combing the TSH algorithm and LR-TDDFTB method. 141,
128, 142, 143 LR-TDDFTB is the tight-binding extension of TDDFT, 144, 141
which has been successfully applied to investigate the optical properties of
plasmonic NPs. 145, 146 It enables calculations of plasmonic excitation along
with several hundred excitation states per time step. The relaxation processes
of plasmon-induced by electron-phonon interactions are treated by the TSH
algorithm. With the NAMD-DFTB method, we demonstrate the plasmon relaxation of
an Au20 cluster. Our simulations show that Au20 can support plasmon-like
excitation, which includes the superposition of multiple single-particle
excitation components.

The numerically efficient LR-TDDFTB method allows us to address a dense
manifold of excited states to ensure the inclusion of plasmon excitation.
Starting from the photoexcited plasmon states in Au20 cluster, we find that
the time constant for relaxation from plasmon excited states to the lowest
excited states is about 2.7 ps, mainly resulting from a step-wise decay
process caused by low-frequency phonons of the Au20 cluster. Furthermore, our
simulations show that the lifetime of the phonon-induced plasmon dephasing is
$\sim$10.4 fs, and such a swift process can be attributed to the strong
nonadiabatic effect in small clusters. Our simulations demonstrate a detailed
description of the dynamic processes in nanoclusters, including plasmon
excitation, hot carrier generation from the plasmon excitation dephasing, and
the subsequent phonon-induced relaxation process.

The NEXMD-DFTB method was employed to investigate the plasmon-induced bond
activation of CO adsorbed on Au${20}$ cluster 87. The simulations provide a
comprehensive and accurate depiction of the multiple dynamic processes
involved in plasmon-mediated chemistry, including plasmon excitation, HEs
relaxation, direct/indirect HE transfer, and the activation of CO vibration
mode induced by HE transfers. These simulations reveal the critical role of
charge transfer (CT) states in plasmon-induced CO activation. The CT states
excite both direct and indirect HE transfer, which leads to the activation of
CO stretching mode, an essential component of plasmon energy relaxation. Our
simulations demonstrate the high efficiency of CO vibrational mode activation,
achieving a success rate of approximately 40%. Notably, the HE transfer occurs
at a faster rate than the conventional scattering process of the Au${20}$,
completing within approximately 100 fs, while the energy relaxation takes
place over a timescale of approximately 1 ps. Furthermore, our direct
atomistic simulations provide detailed insights into the potential energy
evolution during the plasmon-mediated chemical transformation, thereby
enabling a comprehensive understanding of the energy relaxation and HE
transfers during the reaction, as well as elucidating the reaction pathways.

##  4 Polariton chemistry

###  4.1 Introduction and background

When the coupling of the light and molecular DOFs is strong compared to the
dissipation of the excited state population (e.g., cavity leakage or
spontaneous/stimulated electronic emission), the system is said to be strongly
coupled. Generally, the strong-coupling regime is characterized by the
coherent energy exchange between the photon field and the electronic emitter
(i.e. Rabi oscillations). In this regime, the electronic and photon subsystems
can no longer be treated separately, making the accurate simulation of
polariton dynamics very challenging. In this section, we explore various
theoretical approaches in modeling the light-matter interactions present in
optical cavities, such as Fabry–Pérot, where the coherent exchange in energy
between all DOFs directly affects the resulting chemistry. Here, the
electronic and photonic DOFs must be treated quantum mechanically on equal
footing.

Recent experiments have shown a propensity to change chemistry via the
coupling of quantized radiation and various molecular DOFs, namely electronic
and vibrational strong coupling. The Fabry–Pérot-like cavities offer an
extremely tunable cavity frequency (via the effective length between the
cavity mirrors) while exhibiting widely varying coupling strengths highly
dependent on the experimental setup 147, 147, 53, 148, 149, 150, 151, 152,
153, 154, 155. There are many open questions regarding these experiments, such
as the collective effects (i.e., many-molecule or many-mode effects), which
are present in many of the recent works due to the complexity of performing
single-molecule experiments.

###  4.2 Theoretical and Computational Challenges

There are conceptual and technical challenges in performing simulations of
polaritonic systems. The first theoretical hurdle is that of the Hamiltonian
itself. For the Pauli-Fierz Hamiltonian (Eq. LABEL:EQ:_H_PF), the light-matter
interaction requires explicit knowledge of the molecular dipole operator
$\hat{D}$ in the working basis. It turns out that a significant simplification
can be made if one neglects all the contributions from the dipole matrix
except for the ground-to-excited transition dipole matrix elements. There are
two primary reasons for making this approximation: (I) The entirety of this
quantity is not usually printed by default when executing standard electronic
structure software for electronically excited states; however, the usual
information given from these calculations is the ground-to-excited transition
dipole moment, $D_{0J}$, which, for example in linear response calculations
(e.g., TD-HF or TD-DFT), is a trivial quantity to achieve since the result of
such schemes is the ground-to-excited transition density.156 (II)
Historically, the field of quantum optics was narrowly focused on the light-
matter interaction between single atoms and the quantized cavity field. In
such cases, the diagonal elements of the dipole operator are zero by
construction, and the electronic energy differences between excited states
were large, leading to neglecting of the excited-to-excited dipole coupling.
With these approximations, one would find that a reduced Hamiltonian can be
achieved, commonly referred to as the Jaynes-Cumming (JC) Hamiltonian.

###  4.3 Common Approximations Toward Historic Quantum Optics Models

The PF Hamiltonian (Eq. 11) is derived by applying the dipole approximation,
which assumes the wavelength of the cavity fields is substantially larger than
the matter system so that the spatial dependence of the transverse fields is
neglected. In addition, many of the recent works on the simulation of ab
initio polaritons have relied on approximate versions of the PF Hamiltonian
that stem from historical applications in the quantum optics community.157,
158, 159, 160, 161, 162, 163 In these heavily approximated Hamiltonians,
usually a large truncation of the electronic and photonic subspaces is also
performed such that only the ground and a single electronic excited state,
$|g\rangle$ and $|e\rangle$ are included while only including the vacuum and
singly excited Fock states, $|0\rangle$ and $|1\rangle$ in a single-mode
cavity $\alpha=0$. The total basis for this simple model is then confined to
$\\{|g,0\rangle,|g,1\rangle,|e,0\rangle,|e,1\rangle\\}$. Note here that
another approximation is that there is no permanent dipole in the ground
$|g\rangle$ or excited $|e\rangle$ electronic states, i.e.,
$\mathbf{D}_{gg}=\mathbf{D}_{ee}=0$. The most commonly used Hamiltonian for
modeling ab initio polaritons is the Jaynes-Cummings Hamiltonian
$\hat{H}_{\mathrm{JC}}$, which can be written as,

$\hat{H}_{\mathrm{JC}}=\hat{H}_{m}+\hat{H}_{p}+\sqrt{\frac{\omega_{\mathrm{c}}}{2}}\bm{\lambda}\cdot\mathbf{D}_{ge}(\hat{\sigma}\hat{a}^{\dagger}+\hat{\sigma}^{\dagger}\hat{a}),$
(43)

where $\hat{\sigma}^{\dagger}$ ($\hat{\sigma}$) is the creation (annihilation)
operator for the molecule excitation between the ground $g$ and excited $e$
states. Here, two approximations have been made: (I) the rotating wave
approximation (RWA) – which is to say, neglecting the highly oscillatory
$\hat{\sigma}^{\dagger}\hat{a}^{\dagger}_{\alpha}$ and
$\hat{\sigma}\hat{a}_{\alpha}$ terms in the light-matter interaction – and
(II) neglecting the DSE $\hat{H}_{\mathrm{DSE}}$. This Hamiltonian is valid at
ultra-low coupling strengths where the splitting between the one-photon-
dressed ground state $|g,1\rangle$ and the excited state with zero photons
$|e,0\rangle$ exhibit linear splitting (e.g. Rabi splitting) with an increase
in the coupling strength $\bm{\lambda}$. The other two basis states
$|g,0\rangle$ and $|e,1\rangle$ are completely decoupled from the interaction.

Two other common approximations to the PF Hamiltonian are the Rabi
$\hat{H}_{\mathrm{Rabi}}$ and the rotating wave approximation
$\hat{H}_{\mathrm{RWA}}$. These two approximate Hamiltonians can be written
as,

$\displaystyle\hat{H}_{\mathrm{Rabi}}$
$\displaystyle=\hat{H}_{\mathrm{M}}+\hat{H}_{\mathrm{p}}+\sqrt{\frac{\omega}{2}}\bm{\lambda}\cdot\hat{\mathbf{D}}(\hat{\sigma}+\hat{\sigma}^{\dagger})(\hat{a}^{\dagger}+\hat{a}),$
(44) $\displaystyle\hat{H}_{\mathrm{RWA}}$
$\displaystyle=\hat{H}_{\mathrm{M}}+\hat{H}_{\mathrm{p}}+\sqrt{\frac{\omega}{2}}\bm{\lambda}\cdot\mathbf{D}_{ge}(\hat{\sigma}_{ge}\hat{a}^{\dagger}+\hat{\sigma}_{ge}^{\dagger}\hat{a})+\hat{H}_{\mathrm{DSE}},$
(45)

where the Rabi Hamiltonian
$\hat{H}_{\mathrm{Rabi}}=\hat{H}_{\mathrm{PF}}-\hat{H}_{\mathrm{DSE}}$ and the
RWA Hamiltonian applies the RWA to the PF Hamiltonian. Here,
$\hat{\sigma}_{ge}=|g\rangle\langle e|$ is the annihilation operator of the
two-level electronic system. A more in-depth discussion on these Hamiltonians
and the effects of the dipole self-energy term can be found in Refs. 164, 157,
165, and 166.

Note that in a two-state electronic system with only $g$ and $e$ states with
no permanent dipole – e.g., historically, a single-atom system – then the DSE
term only provides a uniform energy shift to the overall energy, similar to
the zero-point energy of the photon Hamiltonian. In this sense, neglecting
this term historically was well-motivated, but in systems with permanent
dipoles in the ground or excited states (e.g., Li-F) or a system with many
electronic states (e.g., molecules or solid-state materials), then the DSE
term provides non-trivial physics and must be included. More recently, a many-
level, many-mode generalization of these approximated Hamiltonians has been
used, which retains the main approximations in each but is now not identical
to the historical definition of each.157

Since the correct PF Hamiltonian in the dipole gauge has been known, the
question remains of why the community returns to the approximated Hamiltonians
for ab initio as well as model calculations. There are many subtleties to
using the full PF Hamiltonian. From the electronic structure perspective, the
many-level dipole operator needs to be computed as well as its square. The DSE
term provides a very complicated description of the system for large coupling
strengths since the dipole matrix, $\hat{\mathbf{D}}$, in realistic molecules
is far from sparse with its square leading to further complications.164, 167
This allows for strong coupling between arbitrary states that is not trivial
to know a priori based on chemical or physical intuition. Often, these
Hamiltonians are parameterized based on cavity-free electronic structure
calculations to obtain the energies and dipoles of the electronic adiabatic
states (see more details in Sec. 4.4.2). In this case, the number of included
electronic (as well as photonic) basis states should be treated as a
convergence parameter. In this case, the DSE causes mixing between far-
separated-in-energy electronic states, which leads to slow convergence in the
basis set size for the matter DOFs,157 and in general poor results compared to
benchmarks, especially at large light-matter coupling strengths when using a
small electronic basis.168, 169, 170 This is precisely why many authors are
developing self-consistent formulations to construct the low-energy
polaritonic states without the need to calculate all the high-energy
electronic states (see Sec. 4.4.2).168, 171, 170, 172, 167, 173, 174, 175.
Further, when considering the quantum dynamical propagation of polaritons in
ab initio systems, additional nuclear gradients are required compared to
cavity-free simulations where only the gradients of the adiabatic state
energies are required to propagation the quantum dynamics, namely the nuclear
gradients on the adiabatic dipole matrix and its square (i.e.,
$\nabla_{R}\hat{\mathbf{D}}$ and $\nabla_{R}\hat{\mathbf{D}}^{2}$).158

###  4.4 Theoretical modeling of polariton chemistry

####  4.4.1 Cavity Born-Oppenheimer approximation

Within the Born-Oppenheimer (BO) approximation,176 the total electronic-
photonic-nuclear can be factorized as,

$\Phi(\mathbf{r},\mathbf{R},{\bf
q}_{\alpha})=\chi(\mathbf{R})\Psi(\mathbf{r},{\bf q}_{\alpha};\mathbf{R}),$
(46)

where $\chi(\mathbf{R})$ and $\Psi(\mathbf{r},{\bf q}_{\alpha};\mathbf{R})$
are the nuclear and polaritonic (i.e., electronic and photonic) wavefunctions,
respectively. Note here that the polaritonic wavefunction is parameterized by
the nuclear positions, exactly like the case without photonic DOFs outside the
cavity. Further, one can invoke the usual Born-Huang-like expansion over the
Born-Oppenheimer factorization as,

$\Phi(\mathbf{r},\mathbf{R},{\bf
q}_{\alpha})=\sum_{\mu}\chi_{a}(\mathbf{R})\psi_{\mu}(\mathbf{r},{\bf
q}_{\alpha};\mathbf{R}),$ (47)

where $\psi_{\mu}(\mathbf{r},{\bf q}_{\alpha};\mathbf{R})$ are the BO
wavefunctions analogous to those as outputted in standard electronic structure
packages for the ground and excited adiabatic states. In this basis, which we
will call the adiabatic polaritonic basis to draw a direct connection to the
bare electronic case, we will discuss various ways to calculate such
polaritonic wavefunctions $\psi_{\mu}(\mathbf{r},{\bf q}_{\alpha};\mathbf{R})$
from ab initio calculations in a variety of approaches and levels of
approximation.

In the following three sections, we will explore ways to obtain rigorous
nuclear-position-parameterized wavefunctions for the entangled adiabatic
electron-photon states $\psi_{\mu}(\mathbf{r},{\bf q}_{\alpha};\mathbf{R})$.
First a brief description of a direct diagonalization approach (Sec. 4.4.2)
with Hamiltonians parameterized with information from standard electronic
structure, while the following sections (Secs. 4.4.3 and 4.4.4) will focus on
the self-consistent approach toward re-developing the standard many-body
schemes in electronic structure theory (e.g., HF, DFT, CCSD, etc.) for the QED
Hamiltonian (e.g., QED-HF, QED-DFT, etc.).

In principle, one can keep a less restrictive expression for the total
wavefunction which is a single exact product of the two distributions in three
possible factorization schemes,

$\displaystyle\Phi(\mathbf{r},\mathbf{R},{\bf q}_{\alpha};t)$
$\displaystyle=\chi(\mathbf{R};t)\Psi(\mathbf{r},{\bf q};\mathbf{R},t)$
$\displaystyle=\chi(\mathbf{r};t)\Psi(\mathbf{R},{\bf q};\mathbf{r},t)$ (48)
$\displaystyle=\chi({\bf q};t)\Psi(\mathbf{r},\mathbf{R};{\bf q},t)$

where $\chi$ is the marginal probability distribution chosen to be the exact
distribution for either of the three unique DOFs. The other factor $\Psi$
describes the combined wavefunction of the leftover DOFs that, for example,
when $\Psi=\Psi(\mathbf{r},{\bf q};\mathbf{R},t)$, satisfies,

$1=\prod_{\alpha}\int d{\bf q}_{\alpha}d\mathbf{r}\Psi(\mathbf{r},{\bf
q}_{\alpha};\mathbf{R},t).$ (49)

In this formalism, the exact time-dependent dynamics is solved via a set of
coupled equation that involves the time-dependent potential energy surface
(TDPES)
$\epsilon(\mathbf{R},t)=\langle\Psi(\mathbf{r};\mathbf{R},t)|\hat{H}_{\mathrm{el}}(\mathbf{r};\mathbf{R},t)-i\partial_{t}|\Psi(\mathbf{r};\mathbf{R},t)\rangle$
dictated by the instantaneous $\Psi$ at $\mathbf{R}$ and the time-dependent
vector potential which takes the form
$A(\mathbf{R};t)=\langle\Psi(\mathbf{r};\mathbf{R},t)|-i\nabla|\Psi(\mathbf{r};\mathbf{R},t)\rangle$.
More details can be found in the original works.177, 178, 179, 180, 181, 182,
183, 184, 185

####  4.4.2 Direct Diagonalization

Historically, the application of QED theory was focused on atomic physics.
Here, single atoms were subjected to a cavity photon yielding Jaynes-Cummings-
like linear splitting (see Eq. 43) since the atomic dipole matrix was
approximated as a single matrix element between the ground and excited state.
In this approximation, the QED Hamiltonian was “parameterized” by the energies
of two electronic adiabatic states, $E_{g},E_{e}$, and transition dipole
element $\bm{\mu}_{eg}$ connecting the ground and excited state. For ab initio
systems, this single dipole matrix element can be approximated by a Taylor-
Series expansion around the Frank-Condon point.159, 163, 161 In this way, one
can directly diagonalize this 4x4 Hamiltonian in the
$\\{|g,0\rangle,|g,1\rangle,|e,0\rangle,|e,1\rangle\\}$ using standard, exact
diagonalization techniques. The left basis $\\{|g\rangle,|e\rangle\\}$ states
are the electronic adiabatic states that are eigenfunctions of the electronic
Hamiltonian
$\hat{H}_{\mathrm{el}}=\hat{T}_{e}+\hat{V}=\hat{H}_{\mathrm{M}}-\hat{T}_{R}$,
where $\hat{V}$ includes all Coulomb interactions between nuclei and
electrons. The right basis states $\\{|0\rangle,|1\rangle\\}$ are the vacuum
and single-photon Fock states, which are eigenstates of the photonic
Hamiltonian
$\hat{H}_{\mathrm{p}}=\hbar\omega_{\alpha}(\hat{a}_{\alpha}^{\dagger}\hat{a}_{\alpha}+\frac{1}{2})$.
In this minimal basis, the JC Hamiltonian can be written as,

$(\hat{H}_{\mathrm{JC}})_{IJ,nm}\dot{=}\begin{bmatrix}E_{g}&0&0&0\\\
0&E_{e}&\sqrt{\frac{\omega_{\mathrm{c}}}{2}}\bm{\lambda}\cdot\mathbf{D}_{eg}&0\\\
0&\sqrt{\frac{\omega_{\mathrm{c}}}{2}}\bm{\lambda}\cdot\mathbf{D}_{eg}&E_{g}+\hbar\omega_{0}&0\\\
0&0&0&E_{e}+\hbar\omega_{0}\end{bmatrix},$ (50)

where the lowest and highest-energy basis states
$\\{|g,0\rangle,|e,1\rangle\\}$ are completely decoupled from the other basis
states as well as each other. In this sense, the polaritonic ground state
$|\Psi_{\mu=0}\rangle=|g,0\rangle$ and the fourth polaritonic state
$|\Psi_{\mu=4}\rangle=|e,1\rangle$, which is valid for ultra-low light-matter
coupling strengths where the linear splitting is valid at the resonance
condition. For non-linear effects, such as those including permanent dipoles,
off-resonant excitations, or simply the exploration of the ground polaritonic
state. To make this leap, one needs to make two major changes: (I) invoke the
more rigorous Pauli-Fierz (PF) Hamiltonian $\hat{H}_{\mathrm{PF}}$ (see Eq.
11) and (II) treat the number of basis electronic and photonic states
convergence parameters.

The matrix elements of the single-mode PF Hamiltonian in this same adiabatic-
Fock basis – but now extended to include, in principle, an infinite number of
electronic and photonic states – can be written as,

$\displaystyle(\hat{H}_{\mathrm{PF}})_{IJ,nm}$
$\displaystyle=\bigg{[}E_{I}+\omega_{\alpha}(n+\frac{1}{2})\bigg{]}\delta_{IJ}\delta_{nm}$
(51)
$\displaystyle+\sqrt{\frac{\omega_{\alpha}}{2}}\bm{\lambda}\cdot\mathbf{D}_{IJ}(\sqrt{m+1}\delta_{n,m+1}+\sqrt{m}\delta_{n,m-1})$
$\displaystyle+\frac{1}{2}\sum_{K}^{\mathcal{N}_{\mathrm{el}}}(\bm{\lambda}\cdot\mathbf{D}_{IK})(\bm{\lambda}\cdot\mathbf{D}_{KJ})\delta_{nm},$

where $\mathcal{N}_{\mathrm{el}}$ is the number of adiabatic electronic states
included in the basis. It is important to note that in this basis, the PF
Hamiltonian is extremely sparse since the coupling elements only connect
adjacent Fock states via the molecular dipole matrix since the matrix elements
of the photonic coordinate $\hat{q}$ are that of the harmonic oscillator
(i.e., only have nonzero super- and sub-diagonal elements). Further, the DSE
contributions, in general, connect all of the electronic states of the system
(with the same photon number) and is the most non-trivial aspect of this
Hamiltonian and will vary strongly between molecular systems.

In this section, we have chosen to expand on the JC and PF Hamiltonians,
confining ourselves to the adiabatic-Fock basis; however, other options exist,
such as the generalized coherent states186, 172, 187 or polarized Fock
states.188 In both of these bases, the photonic states are chosen such that
the molecular dipole parameterizes the photonic state, thereby, in principle,
reducing the convergence of the photonic basis. Further, this “direct
diagonalization” approach to solving the PF Hamiltonian is dependent on this
basis convergence, and the electronic basis is much more rigid since the
adiabatic basis is ubiquitously used for its convenience. However, the
convergence of this basis has yet to be thoroughly tested for a wide range of
systems in solving the PF Hamiltonian, but it is expected to converge slowly
due to the contributions from the DSE term.167, 157 This evidences the need to
move to a more rigorous self-consistent solution for the polaritonic adiabatic
states as is done for the electronic adiabatic states themselves.

####  4.4.3 Self-consistent Polaritonic Single-particle Approaches

Ground State scQED Hartree-Fock
The mean-field approach to the QED Hamiltonian can be first cast in an
identical way as the standard Hartree-Fock procedure in a larger Hilbert
space, including the photonic states. A useful basis, referred to as
generalized coherent states (GCS), can be performed such that the light-matter
interaction part of the PF Hamiltonian can be shifted away. In the following
discussion, we use the notation of Ref. 172; although many works use a similar
scheme.187, 189, 172 Assuming that the ground state HF wavefunction is
$|\mathrm{HF}\rangle$ without the influence of the cavity, then the PF
Hamiltonian can be partially diagonalized in the electronic sub-space as,

$\displaystyle\langle\mathrm{HF}$
$\displaystyle|\hat{H}_{\mathrm{PF}}|\mathrm{HF}\rangle=E_{\mathrm{HF}}+\omega_{\alpha}\big{(}\hat{a}^{\dagger}_{\alpha}\hat{a}_{\alpha}+\frac{1}{2}\big{)}$
(52)
$\displaystyle+\sqrt{\frac{\omega_{\alpha}}{2}}\langle\bm{\lambda}_{\alpha}\cdot\hat{\mathbf{D}}\rangle_{\mathrm{HF}}\cdot(\hat{a}^{\dagger}_{\alpha}+\hat{a}_{\alpha})+\frac{1}{2}\langle(\bm{\lambda}_{\alpha}\cdot\hat{\mathbf{D}})^{2}\rangle_{\mathrm{HF}},$

where
$\langle\cdots\rangle_{\mathrm{HF}}=\langle\mathrm{HF}|\cdots|\mathrm{HF}\rangle$
is the HF ground state expectation value of the electronic subsystem.
Introducing the GSC transformation as,

$\hat{U}({z_{\alpha}})=e^{z_{\alpha}\hat{a}_{\alpha}^{\dagger}-z_{\alpha}^{*}\hat{a}_{\alpha}},$
(53)

where ${\bf z}=(z_{0},z_{1},...)$ is a vector of arbitrary complex numbers
specific to each cavity mode, we can completely remove the light-matter
interaction from the Hamiltonian with a specific choice of $z_{\alpha}$ as,

${\bf
z}\rightarrow-\frac{\langle\bm{\lambda}_{\alpha}\cdot\hat{\mathbf{D}}\rangle_{\mathrm{HF}}}{\sqrt{2\omega_{\alpha}}}.$
(54)

This rotation results in a shift of the photonic operators by ${\bf z}$.
Transforming $\hat{H}_{\mathrm{PF}}$ as $\hat{U}({\bf
z})\langle\mathrm{HF}|\hat{H}_{\mathrm{PF}}|\mathrm{HF}\rangle\hat{U}^{\dagger}({\bf
z})=\hat{H}_{\mathrm{pl}}({\bf z})$, which can be written as,

$\displaystyle\hat{H}_{\mathrm{pl}}({\bf z})$
$\displaystyle=\langle\mathrm{HF}|\hat{U}({\bf
z})\hat{H}_{\mathrm{PF}}\hat{U}^{\dagger}({\bf z})|\mathrm{HF}\rangle$ (55)
$\displaystyle=E_{\mathrm{HF}}+\omega_{\alpha}\big{(}\hat{a}_{\alpha}^{\dagger}\hat{a}_{\alpha}+\frac{1}{2}\big{)}+\frac{1}{2}\langle\big{(}\bm{\lambda}_{\alpha}\cdot\Delta\hat{\mathbf{D}}\big{)}^{2}\rangle_{\mathrm{HF}},$

where the light-matter coupling is inherently removed in this basis since
$\langle\Delta\hat{\mathbf{D}}\rangle_{\mathrm{HF}}=\langle\hat{\mathbf{D}}-\langle\hat{\mathbf{D}}\rangle_{\mathrm{HF}}\rangle_{\mathrm{HF}}=0$.
Note that this removal only exists for mean-field-like methods, which do not
include direct interactions between the ground and excited states. Now, the
light-matter interaction only appears through the DSE-like term, which now
depends not on the dipole $\mathbf{D}$ but rather on the variance in the
ground state permanent dipole,

$\langle(\Delta\hat{\mathbf{D}})^{2}\rangle_{\mathrm{HF}}=\langle\hat{\mathbf{D}}^{2}\rangle_{\mathrm{HF}}-\langle\hat{\mathbf{D}}\rangle_{\mathrm{HF}}^{2},$
(56)

which can be understood as dipole fluctuations inducing light-matter
interactions. In this basis, the Hamiltonian can now be diagonalized in the
photonic sub-space with photon number states. The ground polaritonic state is
then defined as the product $|\Psi_{0}({\bf r},{\bf
q_{\alpha}};R)\rangle=|\mathrm{HF}\rangle\otimes|0\rangle$ The HF energy can
be written as,

$\displaystyle E_{\mathrm{QED-HF}}$
$\displaystyle=E_{\mathrm{HF}}+\sum_{\alpha}\frac{1}{2}\langle(\bm{\lambda}_{\alpha}\cdot\Delta\hat{\mathbf{D}})^{2}\rangle_{\mathrm{HF}}$
(57)
$\displaystyle=E_{\mathrm{HF}}+\sum_{i\nu}\sum_{\alpha}(\bm{\lambda}_{\alpha}\cdot\mathbf{D}_{ia})^{2}$

where $i$ and $a$ are the occupied and virtual single-particle state labels.
This energy is variationally minimized to achieve convergence in the basis of
canonical spin-orbitals.

QED Density Functional Theory

Similarly to the scQED-HF approach, the self-consistent QED density functional
theory (scQED-DFT) is composed in a similar way, where the main elements of
DFT remain, such as the exchange-correlation functional of the density.
However, there are additional photonic DOF that influence the electronic
subsystem. As for the scQED-HF case, the generalized coherent state basis
removes the light-matter coupling term while leaving the DSE term. In this
sense, the energy of the Kohn-Sham energy can be written as,

$E_{\mathrm{QED-
DFT}}=E_{\mathrm{DFT}}+\frac{1}{2}\sum_{\alpha}\langle(\bm{\lambda}_{\alpha}\cdot\hat{\mathbf{D}})^{2}\rangle,$
(58)

which is then minimized with respect to the Kohn-Sham spin-orbitals. There are
many ways to set up the scQED-DFT problem, such as employing a novel exchange-
correlation functional to account for light-matter correlation effects190,
191, 192, 193 or working with electron-only exchange-correlation functionals
using the coherent state basis for the photonic DOFs.171, 157, 172 In any
case, the resulting ground state is uncorrelated by nature of the DFT
formalism.

Excited States scQED-TD-(HF,DFT)

The time-dependent analogues to the aforementioned single-particle approaches
are powerful tools to probe non-equilibrium densities that give rise to
electronic excited states. One of the most popular approaches is one of linear
response (LR), resulting in LR-TD-HF and LR-TD-DFT for bare molecular systems
in the random phase approximation (RPA). Although, it should be noted that the
real-time propagation of the single-particle density matrix – leading to the
real-time TD-HF and real-time TD-DFT approaches – is, in principle, a more
robust approach but one that is usually more costly than that of linear
response. Such schemes have already been developed for the simulation of
molecular polaritons using classical photon DOFs.194, 195 Here, we will focus
our attention on the LR formalism, specifically using a Casida-like approach
to writing the random phase approximation (RPA), originally formulated by
Flick and co-workers168 and rewritten in the language of Casida by the groups
of Shao157 and DePrince.171 It should be noted that other formulations of CIS-
like excited states can be found in the community, such as the non-Hermitian
CIS aimed at simulating cavity loss via a complex photon frequency170.

As per usual, and following the notation of Ref. 171, the LR-TD-HF and LR-TD-
HF eigenvalue equations using the Casida formalism can be written to include
the QED components that satisfy the Pauli-Fierz Hamiltonian (Eq. 11) as,

$\displaystyle\begin{bmatrix}{\bf A}+\bm{\Delta}&{\bf
B}+\bm{\Delta}^{\prime}&\mathbf{g}^{\dagger}&\mathbf{g}^{\dagger}\\\ {\bf
B}+\bm{\Delta}^{\prime}&{\bf
A}+\bm{\Delta}&\mathbf{g}^{\dagger}&\mathbf{g}^{\dagger}\\\
\mathbf{g}&\mathbf{g}&{\bf\omega}&{\bf 0}\\\ \mathbf{g}&\mathbf{g}&{\bf
0}&{\bf\omega}\\\ \end{bmatrix}$ $\displaystyle\begin{bmatrix}{\bf X}\\\ {\bf
Y}\\\ {\bf N}\\\ {\bf M}\end{bmatrix}$
$\displaystyle=\Omega\begin{bmatrix}{\bf 1}&{\bf 0}&{\bf 0}&{\bf 0}\\\ {\bf
0}&{\bf-1}&{\bf 0}&{\bf 0}\\\ {\bf 0}&{\bf 0}&{\bf 1}&{\bf 0}\\\ {\bf 0}&{\bf
0}&{\bf 0}&{\bf-1}\end{bmatrix}\begin{bmatrix}{\bf X}\\\ {\bf Y}\\\ {\bf N}\\\
{\bf M}\end{bmatrix}$ (59)

where A and B are the usual excitation and de-excitation block matrices
written using the TD-HF and TD-DFT approximations,

$\displaystyle A^{\mathrm{TD-HF}}_{ia,jb}$
$\displaystyle=\delta_{ij}\delta_{ab}(\epsilon_{i}-\epsilon_{a})+\langle
ja||bi\rangle,\ \ \ \ B^{\mathrm{TD-HF}}_{ia,jb}=\langle ij||ab\rangle$ (60)
$\displaystyle A^{\mathrm{TD-DFT}}_{ia,jb}$
$\displaystyle=\delta_{ij}\delta_{ab}(\epsilon_{j}-\epsilon_{a})+f^{\mathrm{h/xc}}_{ja,bi},\
\ \ \ B^{\mathrm{TD-DFT}}_{ia,jb}=f^{\mathrm{h/xc}}_{ij,ab}$ (61)

with $\langle ja||bi\rangle=\langle ja|bi\rangle-\langle ja|ib\rangle$ as the
standard two-body interaction terms, $f^{\mathrm{h/xc}}_{ia,jb}$ as the
combination of the Hartree (h) and exchange-correlation (xc) terms from a
standard DFT formalism, and $\epsilon_{a}$ is the energy of single-particle
orbital $a$. Recall that $i,j$ denote occupied orbitals while $a,b$ denote
unoccupied orbitals. The QED-specific terms can be written similarly as,

$\displaystyle g_{ia}$
$\displaystyle=\sum_{\alpha}\sqrt{\frac{\omega_{\alpha}}{2}}D^{\bm{\lambda}_{\alpha}}_{ia}$
(62) $\displaystyle\Delta_{ia,jb}$
$\displaystyle=\sum_{\alpha}D^{\bm{\lambda}_{\alpha}}_{ia}D^{\bm{\lambda}_{\alpha}}_{bj}-D^{\bm{\lambda}_{\alpha}}_{ij}D^{\bm{\lambda}_{\alpha}}_{ab}$
(63) $\displaystyle\Delta_{ia,jb}^{\prime}$
$\displaystyle=\sum_{\alpha}D^{\hat{u}_{\alpha}}_{ia}D^{\bm{\lambda}_{\alpha}}_{jb}-D^{\bm{\lambda}_{\alpha}}_{ib}D^{\bm{\lambda}_{\alpha}}_{aj}$
(64)

Here,
$D^{\bm{\lambda}_{\alpha}}_{ia}=\bm{\lambda}_{\alpha}\cdot\mathbf{D}_{ia}$
with $\hat{{\bf u}}_{\alpha}$ as the polarization of the
$\alpha_{\mathrm{th}}$ cavity mode. Further, $\Delta=\Delta^{\prime}$ if the
molecular orbitals are real-valued.171

An important distinction between various implementations of the QED-TD-DFT
approaches in the community is whether the single-particle orbitals used in
the formulation are “relaxed” in the presence of the cavity or are simply the
bare electronic single-particle states. For example, in Ref. 157, the orbitals
are not relaxed while in the Refs. 168 and 171 the orbitals are relaxed. While
it is clear that using a relaxed reference state for the basis of the RPA
equations would provide a more rigorous result, it is not clear whether
identical results can be obtained in the infinite basis limit of both
approaches, i.e., including more single-particle states in the CIS-like
expansion in excited Slater determinants. Since the RPA equations are
iteratively solved, the expansion coefficients of the excited Slater
determinants may result in the same excited state observables in the infinite
basis limit, while for a finite basis, it may not.

####  4.4.4 Self-consistent Polaritonic Coupled Cluster Approaches (CC,EOM-CC)

Despite computational efficiency, DFT or mean-field HF methods usually
underestimate the correlations. In particular, the mean-field method cannot
describe the electron-photon correlation, and the exchange-correlation
function for electron-photon interaction is unknown. Consequently, the QED
counterpart of coupled-cluster theory (QED-CC) is proposed.174, 173, 175, 196,
36, 37, 197 Similar to the conventional CC theory, QED-CC employs an
exponential wavefunction Ansatz to derive the ground state,

$\ket{\Psi_{CC}}=e^{\hat{T}}\ket{\Phi_{0}},$ (65)

where $\ket{\Phi_{0}}$ is the reference wave function which is usually chosen
to be the tensor product of Hartree-Fock (HF) determinant and photon vacuum
state, i.e., $\ket{\Phi_{0}}=\ket{\varphi_{0}}\otimes\ket{0}$. $\hat{T}$ is
defined in excitation configurations. Within QED-CC theory, the generalized
excitation operator includes three components,
$\hat{T}=\hat{T}_{e}+\hat{T}_{p}+\hat{T}_{ep}$, including electronic
($\hat{T}_{e}$), photonic ($\hat{T}_{p}$), and coupled electronic-photonic
($\hat{T}_{ep}$) excitations,

$\displaystyle T_{e}=$
$\displaystyle\sum_{ia}t^{a}_{i}\hat{a}^{\dagger}_{a}\hat{a}_{i}+\sum_{ijab}t^{ab}_{ij}\hat{a}^{\dagger}_{a}\hat{a}^{\dagger}_{b}\hat{a}_{j}\hat{a}_{i}+\cdots\equiv\sum^{N_{e}}_{\mu}t_{\mu}\hat{\tau}_{\mu},$
(66) $\displaystyle\quad\hat{T}_{p}=$
$\displaystyle\sum_{\alpha}\gamma_{\alpha}\hat{b}^{\dagger}_{\alpha}+\frac{1}{2}\sum_{\alpha\beta}\gamma_{\alpha\beta}\hat{b}^{\dagger}_{\alpha}\hat{b}^{\dagger}_{\beta}+\cdots\equiv\sum^{N_{p}}_{n}\gamma_{n}\hat{B}_{n},$
(67) $\displaystyle\quad\hat{T}_{ep}=$
$\displaystyle\sum_{ia,\alpha}t^{a}_{i}\hat{a}^{\dagger}_{a}\hat{a}_{i}\hat{b}^{\dagger}_{\alpha}+\frac{1}{2}\sum_{ijab,\alpha\beta}t^{ab}_{ij}\hat{a}^{\dagger}_{a}\hat{a}^{\dagger}_{b}\hat{a}_{j}\hat{a}_{i}\hat{b}^{\dagger}_{\alpha}\hat{b}^{\dagger}_{\beta}+\cdots\equiv\sum_{\mu,n}\chi_{\mu,n}\hat{\tau}_{\mu}\hat{B}_{n}.$
(68)

where $N_{e}$ and $N_{f}$ are the numbers of electrons and photon Fock states,
respectively. $\hat{\tau}_{\mu}=\prod^{\mu}_{k}\hat{E}^{a_{k}}_{i_{k}}$ and
$\hat{E}^{a}_{i}=\hat{a}^{\dagger}_{a}\hat{a}_{i}+h.c.$ are the $\mu$-body and
single-body excitation operators, respectively.
$\hat{B}_{n}=\prod^{n}_{\alpha}b^{\dagger}_{\alpha}$ is the $n$-body photonic
excitation. The parameters ($t_{\mu},\gamma_{n},\chi_{\mu,n}$) are the cluster
amplitudes, which are determined through a set of CC equations. The cluster
operator
$\mu\in\\{\hat{\tau}_{\mu},\hat{B}_{n},\hat{\tau}_{\mu}\hat{B}_{n}\\}$
generates a set of orthogonal excited configurations ($\\{\ket{\mu}\\}$),

$\ket{\mu}=\hat{\mu}\ket{\Phi_{0}}.$ (69)

The combined set ($\\{\ket{\Phi_{0}},\ket{\mu}\\}$) defines the Hilbert
subspace for both ground and excited states. Using the CC ansatz in Eq. 65,
the time-dependent Schrödinger equation leads to

$\hat{H}_{PF}\ket{\Psi_{CC}}=E_{cc}\ket{\Psi_{CC}},$ (70)

A set of Slater determinants ($\\{\ket{HF},\ket{\mu}\\}$) is used as the basis
set to solve the Eq. 70, which leads to

$E_{CC}=\bra{HF}\bar{H}_{PF}\ket{HF}.$ (71)

Where $\bar{H}_{PF}=e^{-\hat{T}}\hat{H}_{PF}e^{\hat{T}}$ is the dressed
Hamiltonian. The cluster amplitudes are determined from the projected
equation,

$\Omega_{\mu}\equiv\bra{\mu}e^{-\hat{T}}\hat{H}e^{\hat{T}}\ket{\Phi_{0}}=0$
(72)

The scQED-CC theory leads to the FCI solution if no truncation on the
excitation operator $\hat{T}$ is applied. However, in order to trade-off
between accuracy and computational efficiency, the excitation operator is
usually truncated at the doubles (CCSD) level. Some tests have been performed
on the level of truncation in the photonic excitations (up to 10) in the CC
operator for model systems.196 However, the scaling of scQED-CCSD is
$\mathcal{O}(N_{el}^{6}N^{M_{p}}_{p})$ in general, where $M_{p}$ is the number
of photon modes. Such a scaling makes it usually a bit expensive to rigorously
test in extended molecular systems that include many electrons and/or when
many photon modes and Fock states are used in the calculations. To further
reduce the computational cost, the generalized coherence state 198 may be used
to reduce the photon basis in the scQED-CC calculations.

Fig. 5: Polariton states of LiH as a function of Li-H bond length. (a) CCSD
(black circles) and scQED-CCSD (solid red curve) ground states at zero light-
matter coupling strength $\lambda=0$ a.u. (b) EOM-CCSD (black circles) triplet
$|T_{1}\rangle$ and singlet $|S_{1}\rangle$ states as well as the low-lying
states ($|T_{1}\rangle$, $|\omega\rangle$, and $|S_{1}|rangle$) calculated at
scQED-EOM-CCSD with light-matter coupling strength $\lambda=0$ a.u. The plot
clearly shows the photon line ($|\omega\rangle$) and perfect overlap of the
$S_{1}$ states. (c) Upper ($|P_{U}\rangle$) and lower ($|P_{L}\rangle$)
polariton states with $\lambda=0.088$ a.u. with the triplet state
$|T_{1}\rangle$ uncoupled.

It should be noted that CC methods exhibit size extensivity even in its
truncated forms, meaning that the sum of the energies of non-interacting
subsystems is equal to the total energy. This is in contrast to CI approaches,
where errors in extensivity can become increasingly large with an increasing
number of subsystems. In addition, CC theory has the advantage of size
extensivity in that the excitation energies of each subsystem do not vary with
the size of the total system when the subsystems are non-interacting.

In addition, the excited polariton states can be computed with the
corresponding EOM formalism, which parameterizes a neutral or charged
excitation by applying an excitation operator to the CC ground state 199

$\ket{\hat{R}}=\hat{R}\ket{\Phi_{cc}}=\hat{R}e^{\hat{T}}\ket{\Phi_{0},0}.$
(73)

Where $\hat{R}$ also include the electronic ($\hat{R}_{e}$), photonic
($\hat{R}_{p}$), and coupled electronic-photonic ($\hat{R}_{ep}$) excitations.
Because the excitation operator, $\hat{R}$, commutes with the excitation
operators in $\hat{T}$, solving this eigenvalue problem is equivalent to
finding a right eigenvector of the similarity transformed Hamiltonian,

$\bra{\mu}\bar{H}\hat{R}^{n}\ket{\Phi_{0},0}=E_{n}\hat{R}^{n}_{\mu}.$ (74)

Here, $E_{n}$ is the energy of the nth excited state, and $\mu$ indexes an
element of the excitation operator $\hat{R}$. The excitation operator,
$\hat{R}$, can be chosen to access charged or neutral excitations.

In practice, the eigenvalue problem is solved by iterative diagonalization.
The necessary equations for the sigma vectors are derived and implemented
efficiently.

The sigma equations of polaritonic EOM-CCSD are defined as 199

$\displaystyle\sigma_{m\mu}$
$\displaystyle=\braket{\mu,m}{\bar{H}\hat{R}}{\Phi_{0},0}$
$\displaystyle=\braket{\Phi_{0},0}{\tau^{\dagger}_{\mu}(\hat{b})^{m}\bar{H}\hat{R}}{\Phi_{0},0}.$
(75)

(note, $\ket{\mu,m}=\tau_{\mu}(\hat{b}^{\dagger})^{m}\ket{\Phi_{0},0}$ and
$\tau_{\mu}$ is single or double excitation operator. While this approach
provides the correct equations, the resulting equations contain hundreds of
terms, even for simple EOM theories, and are, therefore, difficult to
optimize. Hence, Garnet’s paper suggested separately generating the different
sectors of the similarity transformed Hamiltonian,

$\bar{H}_{m\mu,n\nu}=\bra{\Phi_{0},0}\tau^{\dagger}_{\mu,m}\bar{H}\tau_{\nu,n}\ket{\Phi_{0},0}.$
(76)

for all relevant excitations. Then generate equations for the different
sectors of the configuration-interaction-like sigma vector as

$\sigma_{m\mu}=\sum_{n\nu}\bar{H}_{m\mu,n\nu}R_{n\nu}.$ (77)

We recently implemented such QED-CCSD/EOM-CCSD method to compute the polariton
states. With the diagrammatic technique, an auto code generator and optimizer
are developed to generate the QED-CCSD/EOM-CCSD equations to arbitrary photon
order. Fig. 5 is an example of the Li-H bond dissociation curves in a single-
mode optical cavity. Fig. 5a shows the ground state Born-Oppenheimer potential
energy surface for the bare electronic system (black circles) and for the
polaritonic system (red curve) at zero light-matter coupling strength
($\lambda=0$). The excited states inside (red curves) and outside (black
circles) of the cavity are shown at zero coupling strength ($\lambda=0$),
where two multiplicities (singlet $|S_{1}\rangle\equiv|S_{1},0\rangle$ and
triplet $|T_{1}\rangle\equiv|T_{1},0\rangle$) of electronic states are shown
as well as a single cavity state ($|\omega\rangle\equiv|S_{0},1\rangle$).
Here, we have used the notation
$|S_{J}\rangle\otimes|n\rangle=|S_{J},n\rangle$ where the left ket in the
product signifies the electronic DOFs while the right signifies the photonic
Fock state label (i.e., the number of photons in the photon-dressed electronic
state). At finite light-matter coupling of $\lambda=0.088$ a.u., the cavity
state with one photon $|\omega\rangle$ couples strongly with the singlet
excited electronic state with zero photons $|S_{1}\rangle$ while the excited
electronic triplet state $|T_{1}\rangle$ is negligibly affected. The so-called
Rabi splitting appears at the degenerate point between the singlet electronic
state and the cavity state, $R_{\mathrm{Li-H}}\approx 3.75$, forming the upper
($|P_{U}\rangle$) and lower ($|P_{L}\rangle$) polaritonic states with mixed
electronic and photonic character.

The computational efficiency of many-body electronic (or polaritonic)
structure codes is nearly as important as the method itself. We have
implemented various backend options for our code (presented in Fig. 5). The
efficiency of Numpy einsum (CPU-accessible) and Torch einsum (GPU-accessible)
functions are shown in Fig. 6 on a vertical log scale as a function of the
number of orbitals (i.e., electrons) included in the calculation. The GPU
hardware allows the QED-CC code to consistently operate with an order of
magnitude less wall time than the CPU version. These results indicate that the
conversion toward GPU hardware is required for exploring the frontier science
of molecular polaritons.

Fig. 6: Example of running calculations with our modular QED-CC code on
different architectures (CPU vs GPU on a desktop) by simply setting the
backend (see the inset code). $10\times-20\times$ speed up is observed on
Torch GPU.

##  5 Perspective and summary

###  5.1 Multiscale method for light-matter interaction

Though there have been extensive developments in QED methods, most of the
current implementations are based on the dipole approximation. However, the
dipole approximation can fail in many cases. The nanoplasmonic cavities allow
for light localization into deeply subwavelength dimensions, leading to
effective mode volumes as small as a few nm3 or even Å3 (microcavities) 23.
Consequently, the size of molecules becomes comparable to the cavity volumes,
and the widely used dipole approximation breaks down. Moreover, the cavity
confines the photon field in a certain direction when coupling molecules with
many modes inside a cavity (either nanophotonic or nanoplasmonic ones). For a
given probing angle $\theta$ (relative to the normal vector of the cavity),
the photon energy has a certain dispersion function. Under this situation, the
dipole approximation along the in-plane direction of the cavity no longer
holds as well 159, 200. Hence, a general method that can compute the multiple
spatial-dependent cavity eigenmodes $\lambda_{n\mathbf{k}}(\mathbf{r})$ and
light-matter couplings (from Maxwell’s equations) is required.

As argued above, the EM field within the cavity is, in principle, not
homogeneous, and spatial dependence matters in light-matter interactions. The
spatial distribution of the EM field can be solved via standard computational
electromagnetic methods, such as Finite-difference time-domain (FDTD) 201, for
Maxwell’s equations. Consequently, the light-matter coupling strength can stem
from a delicate balance between the spatial dependence of the electronic
wavefunctions and the photonic fields. Therefore, only a quantum model fully
incorporates the inhomogeneities of the exciton transition charge density can
quantitatively describe this interplay. Using a fully first-principles
methodology to describe the quantum chemistry of molecules placed inside the
cavity, we can reveal the limitations of the point-dipole approach to address
the exciton dynamics in weak and strong coupling regimes. Besides, current
methods usually don’t consider the feedback of molecular systems on the
cavities. In many situations, modeling of molecular systems only is
insufficient as the molecular response can significantly affect the EM
distributions, especially in nanoplasmonic cavities. Even though the ab initio
QED methods we previously developed take into account the interplay between
the electronic and photonic DOFs, the molecular response on the EM environment
is not considered.

In principle, we could solve the Schrödinger and Maxwell equations
simultaneously in order to obtain access to the radiated fields and, with it
to the self-consistent evolution of light and matter, i.e., the cavity photon
modes can be affected by the modular dipoles due to the Ampere’s Law
(Helmholtz’s equation):

$\displaystyle\left[\nabla\times\frac{1}{\mu_{r}(\mathbf{r}\omega)}\nabla\times-\omega^{2}\mu_{0}\epsilon_{0}\epsilon(\mathbf{r}\omega)\right]\mathbf{E}(\mathbf{r})=-e\mathbf{J}(\mathbf{r}),\text{
or }$ (78)
$\displaystyle\left(\nabla^{2}-\frac{1}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}\right)\mathbf{A}=-\mu_{0}\mathbf{J}.$ (79)

Where $\mathbf{J}$ is the paramagnetic current density due to the molecular
dipoles. And according to Maxwell’s equations, each current induces an
electromagnetic field for which its precise spatial and polarization structure
depends on the electromagnetic environment—oscillating charges emit light.
Hence, a fully self-consistent QED method should consider the feedback of
molecular dipoles on the cavity properties. The flowchart is demonstrated in
Figure 7.

Fig. 7: Self-consistent loop between two sets of equations. The current
density from the QED solvers is used to update EM fields. The self-consistency
between Maxwell and Schrödinger equation can be bypassed (in the linear-
response regime) by using the pre-calculated Dyadic Green’s functions.

Alternatively, the generated field can be expressed with the help of the
dyadic Green’s tensor,

$\mathbf{E}(\mathbf{r},\omega)=i\mu_{0}\omega\int
d\mathbf{r}^{\prime}G(\mathbf{r},\mathbf{r}^{\prime},\omega)[-e\mathbf{J}(\mathbf{r}^{\prime},\omega)],$
(80)

where Green’s function is the formal solution of Helmholtz’s equation

$\left[\nabla\times\frac{1}{\mu_{r}(\mathbf{r}\omega)}\nabla\times-\omega^{2}\mu_{0}\epsilon_{0}\epsilon(\mathbf{r}\omega)\right]G(\mathbf{r},\mathbf{r}^{\prime},\omega)=\delta(\mathbf{r},\mathbf{r}^{\prime}),$
(81)

with linear media $\epsilon(\mathbf{r})$. The Green’s function in Eq. 81 can
be solved from the Frequency domain FDTD methods.

Utilizing Dyadic Green’s functions offers a powerful approach to simplify the
description of certain components of a matter system while focusing on the
computation of electronic structure only. Such a concept has been proposed in
the past in a semiclassical treatment of light-matter interaction 202, 203.
Introducing Green’s function is particularly advantageous in multicomponent
systems that span different length scales, such as a microscopic molecule
(which is described by the current density $\mathbf{J}$) and macroscopic
solvent (represented by a parameterized dielectric function
$\epsilon(\mathbf{r})$). The EM environment embedded via Dyadic Green’s
function can be obtained numerically using many standard Maxwell solvers with
certain boundary conditions, such as FDTD 201 and method of moments 204,
providing the photon mode structures that are coupled with QED electronic
structure solvers. By computing $\mathbf{G}$ beforehand, the self-consistency
is embedded solely via the current density, eliminating the need to treat QED
electronic structures and Maxwell’s equation simultaneously. The bypassing of
the self-consistency between Maxwell-Schrödinger equations is a significant
strength of Dyadic Green’s function approach. This makes the approach more
accessible and easier to implement in various electronic structure solvers.
But, it should be noted that introducing linear media restricts the evaluation
in the linear response regime.

Moreover, the multiscale methods described above can have a significant
influence other than polariton chemistry. It could affect the study of light-
matter interactions in chemistry, physics, materials science, and energy
science at large, such as plasmon chemistry 21, quantum plasmonics 205,
quantum information transduction 206, and photonic-microelectronics
integration and polaritonic devices 207. The integration with open quantum
system method, particularly time-dependent quantum transport theories 208,
209, 210, will also open the door to simulating light-driven quantum transport
phenomena 211, 212 in either weak or strong coupling regimes.

###  5.2 Summary

Controlling chemistry or molecular properties has been a long-standing Holy
Grail for many decades. Only recently, new possibilities have emerged in the
context of light-matter interactions, via either plasmon or polariton-mediated
chemistry, which holds the promise of providing fundamental new strategies to
control chemical reactions that are completely distinct from traditional ones
(such as electrochemistry, thermochemistry, and photochemistry). Consequently,
the application of light-matter interactions in manipulating chemistry has
attracted increasing experimental and theoretical attention. However, the
inherent multiscale nature of light-matter interaction problems poses
significant challenges for both experimental and theoretical investigations of
the underlying processes. Besides, multiscale processes at different length
and time scales make it difficult to optimize the performance of light-matter
interaction-mediated chemistry, and precise modeling of these processes is
essential for obtaining a comprehensive understanding of the fundamental
mechanisms or current experiments. Hence, the development of multiscale
theoretical and computational models that can precisely describe the dynamical
processes in light-matter interaction-mediated chemistry across different time
and length scales, in conjunction with massively parallel algorithms for
large-scale simulations, is essential for obtaining comprehensive insights
into the fundamental mechanisms of current experiments. These simulations
allow for the modeling of complex, multiscale processes that would be
otherwise difficult or impossible to observe experimentally alone. Such
development will ultimately facilitate the discovery of new reaction pathways
enhanced by light-matter interactions.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

The authors acknowledge support from the Laboratory Directed Research and
Development (LDRD) program of Los Alamos National Laboratory.(Grant No.
20220527ECR). XL and YZ acknowledge support from the US DOE, Office of
Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and
Biosciences Division under Triad National Security, LLC (“Triad”) contract
Grant 89233218CNA000001 (FWP: LANLE3F2). LANL is operated by Triad National
Security, LLC, for the National Nuclear Security Administration of the U.S.
Department of Energy (contract no. 89233218CNA000001).

## Notes and references

  * Fleming and Ratner 2008 G. R. Fleming and M. A. Ratner, _Phys. Today_ , 2008, 61, 28–33.
  * Frisk Kockum _et al._ 2019 A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta and F. Nori, _Nat. Rev. Phys._ , 2019, 1, 19–40.
  * Rivera and Kaminer 2020 N. Rivera and I. Kaminer, _Nat. Rev. Phys._ , 2020, 2, 538–561.
  * Anderson 1972 P. W. Anderson, _Science_ , 1972, 177, 393–396.
  * Novotny and Hecht 2006 L. Novotny and B. Hecht, _Annual Review of Physical Chemistry_ , 2006, 57, 303–331.
  * Weiner and Nunes 2012 J. Weiner and F. Nunes, _Light-Matter Interaction: Physics and Engineering at the Nanoscale_ , Oxford University Press, 2012.
  * Mahan 2000 G. Mahan, _Many-Particle Physics_ , Springer, 2000.
  * Stiles _et al._ 2008 P. L. Stiles, J. A. Dieringer, N. C. Shah and R. P. Van Duyne, _Annu. Rev. Anal. Chem._ , 2008, 1, 601–626.
  * Mukamel 2000 S. Mukamel, _Annu. Rev. Phys. Chem._ , 2000, 51, 691–729.
  * Zhang _et al._ 2014 Y. Zhang, L. Meng, C. Yam and G. Chen, _J. Phys. Chem. Lett._ , 2014, 5, 1272–1277.
  * Yam _et al._ 2015 C. Yam, L. Meng, Y. Zhang and G. Chen, _Chem. Soc. Rev._ , 2015, 44, 1763–1776.
  * Meng _et al._ 2015 L. Meng, C. Yam, Y. Zhang, R. Wang and G. Chen, _J. Phys. Chem. Lett._ , 2015, 6, 4410–4416.
  * Meng _et al._ 2017 L. Meng, Y. Zhang and C. Yam, _J. Phys. Chem. Lett._ , 2017, 8, 571–575.
  * Wang _et al._ 2015 R. Wang, Y. Zhang, G. H. Chen and C. Y. Yam, _Prog. Electromagn. Res._ , 2015, 154, 163–170.
  * Degen _et al._ 2017 C. L. Degen, F. Reinhard and P. Cappellaro, _Rev. Mod. Phys._ , 2017, 89, 035002.
  * Flamini _et al._ 2018 F. Flamini, N. Spagnolo and F. Sciarrino, _Rep. Prog. Phys._ , 2018, 82, 016001.
  * Zhou _et al._ 2018 L. Zhou, Y. Zhang, Z. Zhuo, A. J. Neukirch and S. Tretiak, _J. Phys. Chem. Lett._ , 2018, 9, 6915–6920.
  * Mirkovic _et al._ 2017 T. Mirkovic, E. E. Ostroumov, J. M. Anna, R. van Grondelle, Govindjee and G. D. Scholes, _Chem. Rev._ , 2017, 117, 249–293.
  * Forn-Díaz _et al._ 2019 P. Forn-Díaz, L. Lamata, E. Rico, J. Kono and E. Solano, _Rev. Mod. Phys._ , 2019, 91, 025005.
  * Brongersma _et al._ 2015 M. L. Brongersma, N. J. Halas and P. Nordlander, _Nat. Nanotechnol._ , 2015, 10, 25–34.
  * Zhan _et al._ 2023 C. Zhan, J. Yi, S. Hu, X.-G. Zhang, D.-Y. Wu and Z.-Q. Tian, _Nat. Rev. Methods Primers_ , 2023, 3, 12.
  * Kazuma and Kim 2019 E. Kazuma and Y. Kim, _Angew. Chem. Int. Ed._ , 2019, 58, 4800–4808.
  * Baumberg 2022 J. J. Baumberg, _Nano Lett._ , 2022, 22, 5859–5865.
  * Ribeiro _et al._ 2018 R. F. Ribeiro, L. A. Martínez-Martínez, M. Du, J. C. Gonzalez-Angulo and J. Yuen-Zhou, _Chem. Sci._ , 2018, 9, 6325–6339.
  * Ebbesen 2016 T. W. Ebbesen, _Acc. Chem. Res._ , 2016, 49, 2403–2412.
  * Kéna-Cohen and Forrest 2010 S. Kéna-Cohen and S. R. Forrest, _Nat. Photon._ , 2010, 4, 371–375.
  * Kang _et al._ 2019 J.-W. Kang, B. Song, W. Liu, S.-J. Park, R. Agarwal and C.-H. Cho, _Sci. Adv._ , 2019, 5, eaau9338.
  * Orgiu _et al._ 2015 E. Orgiu, J. George, J. A. Hutchison, E. Devaux, J. F. Dayen, B. Doudin, F. Stellacci, C. Genet, J. Schachenmayer, C. Genes, G. Pupillo, P. Samorì and T. W. Ebbesen, _Nat. Mater._ , 2015, 14, 1123–1129.
  * Zhong _et al._ 2017 X. Zhong, T. Chervy, L. Zhang, A. Thomas, J. George, C. Genet, J. A. Hutchison and T. W. Ebbesen, _Angew. Chem. Int. Ed._ , 2017, 56, 9034–9038.
  * Georgiou _et al._ 2021 K. Georgiou, R. Jayaprakash, A. Othonos and D. G. Lidzey, _Angew. Chem. Int. Ed._ , 2021, 60, 16661–16667.
  * Wang _et al._ 2021 M. Wang, M. Hertzog and K. Börjesson, _Nat. Commun._ , 2021, 12, 1874.
  * Coles _et al._ 2014 D. M. Coles, N. Somaschi, P. Michetti, C. Clark, P. G. Lagoudakis, P. G. Savvidis and D. G. Lidzey, _Nat. Mater._ , 2014, 13, 712–719.
  * Dusel _et al._ 2020 M. Dusel, S. Betzold, O. A. Egorov, S. Klembt, J. Ohmer, U. Fischer, S. Höfling and C. Schneider, _Nat. Commun._ , 2020, 11, 2863\.
  * Zasedatelev _et al._ 2019 A. V. Zasedatelev, A. V. Baranikov, D. Urbonas, F. Scafirimuto, U. Scherf, T. Stöferle, R. F. Mahrt and P. G. Lagoudakis, _Nat. Photon._ , 2019, 13, 378–383.
  * Kavokin _et al._ 2022 A. Kavokin, T. C. H. Liew, C. Schneider, P. G. Lagoudakis, S. Klembt and S. Hoefling, _Nat. Rev. Phys._ , 2022, 4, 435–451.
  * Pavošević _et al._ 2022 F. Pavošević, R. L. Smith and A. Rubio, _Catalysis in Click Chemistry Reactions Controlled by Cavity Quantum Vacuum Fluctuations: The Case of endo/exo Diels-Alder Reaction_ , 2022, http://arxiv.org/abs/2208.06925, arXiv:2208.06925 [physics].
  * Pavošević _et al._ 2022 F. Pavošević, S. Hammes-Schiffer, A. Rubio and J. Flick, _J. Am. Chem. Soc._ , 2022.
  * Schäfer _et al._ 2022 C. Schäfer, J. Flick, E. Ronca, P. Narang and A. Rubio, _Nat. Commun._ , 2022, 13, 7817.
  * Schäfer 2022 C. Schäfer, _J. Phys. Chem. Lett._ , 2022, 13, 6905–6911.
  * Cave and Newton 1997 R. J. Cave and M. D. Newton, _J. Chem. Phys._ , 1997, 106, 9213–9226.
  * Martínez-Martínez _et al._ 2017 L. A. Martínez-Martínez, R. F. Ribeiro, J. C. González-Angulo and J. Yuen-Zhou, _ACS Photonics_ , 2017, 5, 167–176.
  * Yang and Cao 2021 P.-Y. Yang and J. Cao, _J. Phys. Chem. Lett._ , 2021, 12, 9531–9538.
  * Climent and Feist 2020 C. Climent and J. Feist, _Phys. Chem. Chem. Phys._ , 2020, 22, 23545–23552.
  * Wang _et al._ 2022 D. S. Wang, T. Neuman, S. F. Yelin and J. Flick, _J. Phys. Chem. Lett._ , 2022, 13, 3317–3324.
  * Imperatore _et al._ 2021 M. V. Imperatore, J. B. Asbury and N. C. Giebink, _J. Chem. Phys._ , 2021, 154, 191103.
  * Galego _et al._ 2016 J. Galego, F. J. Garcia-Vidal and J. Feist, _Nat. Commun._ , 2016, 7, 13841.
  * Campos-Gonzalez-Angulo _et al._ 2019 J. A. Campos-Gonzalez-Angulo, R. F. Ribeiro and J. Yuen-Zhou, _Nat. Commun._ , 2019, 10, 4685.
  * Philbin _et al._ 2022 J. P. Philbin, Y. Wang, P. Narang and W. Dou, _J. Phys. Chem. C_ , 2022, 126, 14908–14913.
  * Efrima and Bixon 1974 S. Efrima and M. Bixon, _Chem. Phys. Lett._ , 1974, 25, 34–37.
  * Phuc _et al._ 2020 N. T. Phuc, P. Q. Trung and A. Ishizaki, _Sci. Rep._ , 2020, 10, 7318\.
  * Davidsson and Kowalewski 2020 E. Davidsson and M. Kowalewski, _J. Chem. Phys._ , 2020, 153, 234304\.
  * Galego _et al._ 2017 J. Galego, F. J. Garcia-Vidal and J. Feist, _Phys. Rev. Lett._ , 2017, 119, 136001.
  * Mauro _et al._ 2021 L. Mauro, K. Caicedo, G. Jonusauskas and R. Avriller, _Phys. Rev. B_ , 2021, 103, 165412.
  * Vurgaftman _et al._ 2020 I. Vurgaftman, B. S. Simpkins, A. D. Dunkelberger and J. C. Owrutsky, _J. Phys. Chem. Lett._ , 2020, 11, 3557–3562.
  * Hiura and Shalabney 2021 H. Hiura and A. Shalabney, _ChemRxiv_ , 2021.
  * Hiura _et al._ 2019 H. Hiura, A. Shalabney and J. George, _ChemRxiv_ , 2019.
  * Craig and Thirunamachandran 1998 D. Craig and T. Thirunamachandran, _Molecular Quantum Electrodynamics: An Introduction to Radiation-molecule Interactions_ , Dover Publications, 1998\.
  * Archambault _et al._ 2010 A. Archambault, F. m. c. Marquier, J.-J. Greffet and C. Arnold, _Phys. Rev. B_ , 2010, 82, 035411.
  * Andrews _et al._ 1989 D. L. Andrews, D. P. Craig and T. Thirunamachandran, _Int. Rev. Phys. Chem._ , 1989, 8, 339–383.
  * Snyder and Love 1983 A. W. Snyder and J. Love, _Optical Waveguide Theory_ , Springer, 1st edn, 1983\.
  * Chikkaraddy _et al._ 2016 R. Chikkaraddy, B. de Nijs, F. Benz, S. J. Barrow, O. A. Scherman, E. Rosta, A. Demetriadou, P. Fox, O. Hess and J. J. Baumberg, _Nature_ , 2016, 535, 127–30.
  * Benz _et al._ 2016 F. Benz, M. K. Schmidt, A. Dreismann, R. Chikkaraddy, Y. Zhang, A. Demetriadou, C. Carnegie, H. Ohadi, B. de Nijs, R. Esteban, J. Aizpurua and J. J. Baumberg, _Science_ , 2016, 354, 726–729.
  * Neuman _et al._ 2018 T. Neuman, R. Esteban, D. Casanova, F. J. García-Vidal and J. Aizpurua, _Nano Lett._ , 2018, 18, 2358–2364.
  * Wang _et al._ 2021 D. S. Wang, T. Neuman and P. Narang, _J. Phys. Chem. C_ , 2021, 125, 6222–6228.
  * Kelly _et al._ 2003 K. L. Kelly, E. Coronado, L. L. Zhao and G. C. Schatz, _J. Phys. Chem. B_ , 2003, 107, 668–677.
  * Schuller _et al._ 2010 J. A. Schuller, E. S. Barnard, W. Cai, Y. C. Jun, J. S. White and M. L. Brongersma, _Nat. Mater._ , 2010, 9, 193–204.
  * Jiang _et al._ 2003 Jiang, K. Bosnick, M. Maillard and L. Brus, _J. Phys. Chem. B_ , 2003, 107, 9964–9972.
  * Brus 2008 L. Brus, _Acc. Chem. Res._ , 2008, 41, 1742–1749.
  * Nie and Emory 1997 S. Nie and S. R. Emory, _Science_ , 1997, 275, 1102–1106.
  * Zhan _et al._ 2018 C. Zhan, X.-J. Chen, J. Yi, J.-F. Li, D.-Y. Wu and Z.-Q. Tian, _Nat. Rev. Chem._ , 2018, 2, 216–230.
  * Atwater and Polman 2010 H. A. Atwater and A. Polman, _Nat. Mater._ , 2010, 9, 205–213.
  * Zhang _et al._ 2016 Y. Zhang, C. Yam and G. C. Schatz, _J. Phys. Chem. Lett._ , 2016, 7, 1852–1858.
  * Ma _et al._ 2013 R.-M. Ma, R. F. Oulton, V. J. Sorger and X. Zhang, _Laser Photonics Rev._ , 2013, 7, 1–21.
  * Berini and De Leon 2012 P. Berini and I. De Leon, _Nat. Photon._ , 2012, 6, 16–24.
  * Guan _et al._ 2021 J. Guan, R. Li, X. G. Juarez, A. D. Sample, Y. Wang, G. C. Schatz and T. W. Odom, _Adv. Mater._ , 2021, n/a, 2103262.
  * Guan _et al._ 2022 J. Guan, J.-E. Park, S. Deng, M. J. H. Tan, J. Hu and T. W. Odom, _Chem. Rev._ , 2022, 122, 15177–15203.
  * Cushing and Wu 2016 S. K. Cushing and N. Wu, _J. Phys. Chem. Lett._ , 2016, 7, 666–675.
  * Li _et al._ 2015 J. Li, S. K. Cushing, F. Meng, T. R. Senty, A. D. Bristow and N. Wu, _Nat. Photon._ , 2015, 9, 601–607.
  * Zhang 2021 Y. Zhang, _J. Phys. Chem. A_ , 2021, 125, 9201–9208.
  * Clavero 2014 C. Clavero, _Nat. Photon._ , 2014, 8, 95–103.
  * Bernardi _et al._ 2015 M. Bernardi, J. Mustafa, J. B. Neaton and S. G. Louie, _Nat. Commun._ , 2015, 6, 7044.
  * Brown _et al._ 2016 A. M. Brown, R. Sundararaman, P. Narang, I. William A. Goddard and H. A. Atwater, _ACS Nano_ , 2016, 10, 957–966.
  * Manjavacas _et al._ 2014 A. Manjavacas, J. G. Liu, V. Kulkarni and P. Nordlander, _ACS Nano_ , 2014, 8, 7630–7638.
  * Yuan _et al._ 2022 L. Yuan, J. Zhou, M. Zhang, X. Wen, J. M. P. Martirez, H. Robatjazi, L. Zhou, E. A. Carter, P. Nordlander and N. J. Halas, _ACS Nano_ , 2022, 16, 17365–17375.
  * Zhou _et al._ 2020 L. Zhou, J. M. P. Martirez, J. Finzel, C. Zhang, D. F. Swearer, S. Tian, H. Robatjazi, M. Lou, L. Dong, L. Henderson, P. Christopher, E. A. Carter, P. Nordlander and N. J. Halas, _Nat. Energy_ , 2020, 5, 61–70.
  * Zhou _et al._ 2018 L. Zhou, D. F. Swearer, C. Zhang, H. Robatjazi, H. Zhao, L. Henderson, L. Dong, P. Christopher, E. A. Carter, P. Nordlander and N. J. Halas, _Science_ , 2018, 362, 69–72.
  * Wu _et al._ 2023 X. Wu, T. v. d. Heide, T. Frauenheim, S. Tretiak, C. Yam and Y. Zhang, _arXiv Preprint_ , 2023, arXiv:2211.12625, 0.
  * Zhang _et al._ 2018 Y. Zhang, T. Nelson, S. Tretiak, H. Guo and G. C. Schatz, _Acs Nano_ , 2018, 12, 8415–8422.
  * Wu _et al._ 2020 Q. Wu, L. Zhou, G. C. Schatz, Y. Zhang and H. Guo, _J. Am. Chem. Soc._ , 2020, 142, 13090–13101.
  * Kazuma _et al._ 2018 E. Kazuma, J. Jung, H. Ueba, M. Trenary and Y. Kim, _Science_ , 2018, 360, 521–526.
  * Zhang _et al._ 2018 Y. Zhang, S. He, W. Guo, Y. Hu, J. Huang, J. R. Mulcahy and W. D. Wei, _Chem. Rev._ , 2018, 118, 2927–2954.
  * Tesema _et al._ 2019 T. E. Tesema, B. Kafle and T. G. Habteyes, _J. Phys. Chem. C_ , 2019, 123, 8469–8483.
  * Adleman _et al._ 2009 J. R. Adleman, D. A. Boyd, D. G. Goodwin and D. Psaltis, _Nano Lett._ , 2009, 9, 4417–4423.
  * Boerigter _et al._ 2016 C. Boerigter, R. Campana, M. Morabito and S. Linic, _Nat. Commun._ , 2016, 7, 10545–10545.
  * Zhang _et al._ 2019 Y. Zhang, T. Nelson and S. Tretiak, _Computational Photocatalysis: Modeling of Photophysics and Photochemistry at Interfaces_ , ACS Publications, 2019, pp. 239–256.
  * Marinica _et al._ 2012 D. Marinica, A. Kazansky, P. Nordlander, J. Aizpurua and A. G. Borisov, _Nano Lett._ , 2012, 12, 1333–1339.
  * Varas _et al._ 2016 A. Varas, P. García-González, J. Feist, F. García-Vidal and A. Rubio, _Nanophotonics_ , 2016, 5, 409–426.
  * Perera _et al._ 2020 T. Perera, S. D. Gunapala, M. I. Stockman and M. Premaratne, _J. Phys. Chem. C_ , 2020, 124, 27694–27708.
  * Rossi _et al._ 2017 T. P. Rossi, M. Kuisma, M. J. Puska, R. M. Nieminen and P. Erhart, _J. Chem. Theory Comput._ , 2017, 13, 4779–4790.
  * Asadi-Aghbolaghi _et al._ 2020 N. Asadi-Aghbolaghi, R. Rüger, Z. Jamshidi and L. Visscher, _J. Phys. Chem. C_ , 2020, 124, 7946–7955.
  * Pandeya and Aikens 2021 P. Pandeya and C. M. Aikens, _J. Phys. Chem. A_ , 2021, 125, 4847–4860.
  * Alkan and Aikens 2018 F. Alkan and C. M. Aikens, _J. Phys. Chem. C_ , 2018, 122, 23639–23650.
  * Alkan and Aikens 2021 F. Alkan and C. M. Aikens, _J. Phys. Chem. C_ , 2021, 125, 12198–12206.
  * Della Sala 2022 F. Della Sala, _J. Chem. Phys._ , 2022, 157, 104101.
  * D’Agostino _et al._ 2018 S. D’Agostino, R. Rinaldi, G. Cuniberti and F. Della Sala, _J. Phys. Chem. C_ , 2018, 122, 19756–19766.
  * You and Panoiu 2019 J. W. You and N. C. Panoiu, _IEEE J. Multiscale and Multiphys. Comput. Techn._ , 2019, 4, 111–118.
  * Hohenberg and Kohn 1964 P. Hohenberg and W. Kohn, _Phys. Rev._ , 1964, 136, B864–B871.
  * Parr and Yang 1995 R. G. Parr and W. Yang, _Annu. Rev. Phys. Chem._ , 1995, 46, 701–728.
  * Runge and Gross 1984 E. Runge and E. K. U. Gross, _Phys. Rev. Lett._ , 1984, 52, 997–1000.
  * Casida and Huix-Rotllant 2012 M. E. Casida and M. Huix-Rotllant, _Annu. Rev. Phys. Chem._ , 2012, 63, 287–323.
  * Tiago and Chelikowsky 2006 M. L. Tiago and J. R. Chelikowsky, _Phys. Rev. B_ , 2006, 73, 205334\.
  * Román Castellanos _et al._ 2019 L. Román Castellanos, O. Hess and J. Lischner, _Commun. Phys._ , 2019, 2, 47.
  * Wu _et al._ 2022 X. Wu, B. Liu, T. Frauenheim, S. Tretiak, C. Yam and Y. Zhang, _J. Chem. Phys._ , 2022, 157, 214201.
  * Liu _et al._ 2018 J. G. Liu, H. Zhang, S. Link and P. Nordlander, _ACS Photonics_ , 2018, 5, 2584–2595.
  * Breuer _et al._ 2002 H. Breuer, F. Petruccione and S. Petruccione, _The Theory of Open Quantum Systems_ , Oxford University Press, 2002.
  * Haug and Jauho 1998 H. Haug and A. Jauho, _Quantum Kinetics in Transport and Optics of Semiconductors_ , Springer Berlin Heidelberg, 1998.
  * Tesema _et al._ 2017 T. E. Tesema, B. Kafle, M. G. Tadesse and T. G. Habteyes, _J. Phys. Chem. C_ , 2017, 121, 7421–7428.
  * Boerigter _et al._ 2016 C. Boerigter, U. Aslam and S. Linic, _ACS Nano_ , 2016, 10, 6108–6115.
  * Long and Prezhdo 2014 R. Long and O. V. Prezhdo, _J. Am. Chem. Soc._ , 2014, 136, 4343–4354.
  * Long _et al._ 2012 R. Long, N. J. English and O. V. Prezhdo, _J. Am. Chem. Soc._ , 2012, 134, 14238–14248.
  * Bora _et al._ 2016 T. Bora, D. Zoepfl and J. Dutta, _Sci. Rep._ , 2016, 6, 26913\.
  * Sivan _et al._ 2019 Y. Sivan, I. W. Un and Y. Dubi, _Faraday Discuss._ , 2019, 214, 215–233.
  * Brown _et al._ 2017 A. M. Brown, R. Sundararaman, P. Narang, A. M. Schwartzberg, W. A. Goddard and H. A. Atwater, _Phys. Rev. Lett._ , 2017, 118, 087401.
  * Nelson _et al._ 2020 T. R. Nelson, A. J. White, J. A. Bjorgaard, A. E. Sifain, Y. Zhang, B. Nebgen, S. Fernandez-Alberti, D. Mozyrsky, A. E. Roitberg and S. Tretiak, _Chem. Rev._ , 2020, 120, 2215–2287.
  * Tully 1990 J. C. Tully, _J. Chem. Phys._ , 1990, 93, 1061–1071.
  * Tully 2012 J. C. Tully, _J. Chem. Phys._ , 2012, 137, 22A301.
  * Curchod _et al._ 2013 B. F. Curchod, U. Rothlisberger and I. Tavernelli, _Chem. Phys. Chem._ , 2013, 14, 1314–1340.
  * Malone _et al._ 2020 W. Malone, B. Nebgen, A. White, Y. Zhang, H. Song, J. A. Bjorgaard, A. E. Sifain, B. Rodriguez-Hernandez, V. M. Freixas, S. Fernandez-Alberti _et al._ , _J. Chem. Theory Comput._ , 2020, 16, 5771–5783.
  * Fabiano _et al._ 2008 E. Fabiano, T. W. Keal and W. Thiel, _Chem. Phys._ , 2008, 349, 334–347.
  * Makhov _et al._ 2017 D. V. Makhov, C. Symonds, S. Fernandez-Alberti and D. V. Shalashilin, _Chem. Phys._ , 2017, 493, 200–218.
  * Makhov _et al._ 2014 D. V. Makhov, W. J. Glover, T. J. Martinez and D. V. Shalashilin, _J. Chem. Phys._ , 2014, 141, 054110.
  * Makhov _et al._ 2015 D. V. Makhov, K. Saita, T. J. Martinez and D. V. Shalashilin, _Phys. Chem. Chem. Phys._ , 2015, 17, 3316–3325.
  * V. Makhov _et al._ 2016 D. V. Makhov, T. J. Martinez and D. V. Shalashilin, _Faraday Discuss._ , 2016, 194, 81–94.
  * Song _et al._ 2021 H. Song, V. M. Freixas, S. Fernandez-Alberti, A. J. White, Y. Zhang, S. Mukamel, N. Govind and S. Tretiak, _J. Chem. Theory Comput._ , 2021, 17, 3629–3643.
  * Ben-Nun _et al._ 2000 M. Ben-Nun, J. Quenneville and T. J. Martínez, _J. Phys. Chem. A_ , 2000, 104, 5161–5175.
  * Besteiro _et al._ 2017 L. V. Besteiro, X.-T. Kong, Z. Wang, G. Hartland and A. O. Govorov, _ACS Photonics_ , 2017, 4, 2759–2781.
  * Chang _et al._ 2019 L. Chang, L. V. Besteiro, J. Sun, E. Y. Santiago, S. K. Gray, Z. Wang and A. O. Govorov, _ACS Energy Lett._ , 2019, 4, 2552–2568.
  * Yan _et al._ 2016 L. Yan, F. Wang and S. Meng, _ACS Nano_ , 2016, 10, 5452–5458.
  * Akimov _et al._ 2013 A. V. Akimov, A. J. Neukirch and O. V. Prezhdo, _Chem. Rev._ , 2013, 113, 4496–4565.
  * Trani _et al._ 2011 F. Trani, G. Scalmani, G. Zheng, I. Carnimeo, M. J. Frisch and V. Barone, _J. Chem. Theory Comput._ , 2011, 7, 3304–3313.
  * Niehaus _et al._ 2001 T. A. Niehaus, S. Suhai, F. Della Sala, P. Lugli, M. Elstner, G. Seifert and T. Frauenheim, _Phys. Rev. B_ , 2001, 63, 085108.
  * Song _et al._ 2020 H. Song, S. A. Fischer, Y. Zhang, C. J. Cramer, S. Mukamel, N. Govind and S. Tretiak, _J. Chem. Theory Comput._ , 2020, 16, 6418–6427.
  * Wu _et al._ 2022 X. Wu, S. Wen, H. Song, T. Frauenheim, S. Tretiak, C. Yam and Y. Zhang, _J. Chem. Phys._ , 2022, 157, 084114.
  * Hourahine _et al._ 2020 B. Hourahine, B. Aradi, V. Blum, F. Bonafé, A. Buccheri, C. Camacho, C. Cevallos, M. Deshaye, T. Dumitrică, A. Dominguez _et al._ , _J. Chem. Phys._ , 2020, 152, 124101.
  * Douglas-Gallardo _et al._ 2019 O. A. Douglas-Gallardo, M. Berdakin, T. Frauenheim and C. G. Sánchez, _Nanoscale_ , 2019, 11, 8604–8615.
  * Berdakin _et al._ 2019 M. Berdakin, O. A. Douglas-Gallardo and C. G. Sanchez, _J. Phys. Chem. C_ , 2019, 124, 1631–1639.
  * Hutchison _et al._ 2012 J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux and T. W. Ebbesen, _Angew. Chem. Int. Ed._ , 2012, 51, 1592–1596.
  * George _et al._ 2016 J. George, T. Chervy, A. Shalabney, E. Devaux, H. Hiura, C. Genet and T. W. Ebbesen, _Phys. Rev. Lett._ , 2016, 117, 153601.
  * Vergauwe _et al._ 2019 R. M. A. Vergauwe, A. Thomas, K. Nagarajan, A. Shalabney, J. George, T. Chervy, M. Seidel, E. Devaux, V. Torbeev and T. W. Ebbesen, _Angew. Chem. Int. Ed._ , 2019, 58, 15324–15328.
  * Lather _et al._ 2019 J. Lather, P. Bhatt, A. Thomas, T. W. Ebbesen and J. George, _Angew. Chem. Int. Ed._ , 2019, 58, 10635–10638.
  * Sau _et al._ 2021 A. Sau, K. Nagarajan, B. Patrahau, L. Lethuillier-Karl, R. M. A. Vergauwe, A. Thomas, J. Moran, C. Genet and T. W. Ebbesen, _Angew. Chem. Int. Ed._ , 2021, 60, 5712–5717.
  * Hirai _et al._ 2021 K. Hirai, H. Ishikawa, T. Chervy, J. A. Hutchison and H. Uji-i, _Chem. Sci._ , 2021, 12, 11986–11994.
  * Takele _et al._ 2021 W. M. Takele, F. Wackenhut, Q. Liu, L. Piatkowski, J. Waluk and A. J. Meixner, _J. Phys. Chem. C_ , 2021, 125, 14932–14939.
  * 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 and T. W. Ebbesen, _Science_ , 2019, 363, 615–619.
  * Lather _et al._ 2022 J. Lather, A. N. K. Thabassum, J. Singh and J. George, _Chem. Sci._ , 2022, 13, 195–202.
  * Tretiak and Mukamel 2002 S. Tretiak and S. Mukamel, _Chem. Rev._ , 2002, 102, 3171–3212.
  * Yang _et al._ 2021 J. Yang, Q. Ou, Z. Pei, H. Wang, B. Weng, Z. Shuai, K. Mullen and Y. Shao, _J. Chem. Phys._ , 2021, 155, 064107.
  * Zhang _et al._ 2019 Y. Zhang, T. Nelson and S. Tretiak, _J. Chem. Phys._ , 2019, 151, 154109\.
  * Tichauer _et al._ 2021 R. H. Tichauer, J. Feist and G. Groenhof, _J. Chem. Phys._ , 2021, 154, 104112.
  * Groenhof and Toppari 2018 G. Groenhof and J. J. Toppari, _J. Phys. Chem. Lett._ , 2018, 9, 4848–4851.
  * Luk _et al._ 2017 H. L. Luk, J. Feist, J. J. Toppari and G. Groenhof, _J. Chem. Theory Comput._ , 2017, 13, 4324–4335.
  * Groenhof _et al._ 2019 G. Groenhof, C. Climent, J. Feist, D. Morozov and J. J. Toppari, _J. Phys. Chem. Lett._ , 2019, 10, 5476–5483.
  * Tichauer _et al._ 2022 R. H. Tichauer, D. Morozov, I. Sokolovskii, J. J. Toppari and G. Groenhof, _J. Phys. Chem. Lett._ , 2022, 13, 6259–6267.
  * Mandal _et al._ 2022 A. Mandal, M. Taylor, B. Weight, E. Koessler, X. Li and P. Huo, _Theoretical Advances in Polariton Chemistry and Molecular Cavity Quantum Electrodynamics_ , 2022, https://chemrxiv.org/engage/chemrxiv/article-details/63903eee92f084c9612a9086.
  * Rokaj _et al._ 2018 V. Rokaj, D. M. Welakuh, M. Ruggenthaler and A. Rubio, _J. Phys. B: At. Mol. Opt. Phys._ , 2018, 51, 034005.
  * Mandal _et al._ 2020 A. Mandal, T. D. Krauss and P. Huo, _J. Phys. Chem. B_ , 2020, 124, 6321–6340.
  * Weight _et al._ 2023 B. Weight, T. Krauss and P. Huo, _Investigating Molecular Exciton-Polaritons using Many-body Electronic Structure Theory with Cavity Quantum Electrodynamics_ , 2023, https://chemrxiv.org/engage/chemrxiv/article-details/63ba35271699ca2dd7e74772.
  * Flick and Narang 2020 J. Flick and P. Narang, _J. Chem. Phys._ , 2020, 153, 094116.
  * Wang _et al._ 2021 D. S. Wang, T. Neuman, J. Flick and P. Narang, _J. Chem. Phys._ , 2021, 154, 104109.
  * McTague and Foley 2022 J. McTague and J. J. Foley, _J. Chem. Phys._ , 2022, 156, 154103\.
  * Vu _et al._ 2022 N. Vu, G. M. McLeod, K. Hanson and A. E. DePrince, _J. Phys. Chem. A_ , 2022, 126, 9303–9312.
  * Haugland _et al._ 2020 T. S. Haugland, E. Ronca, E. F. Kjønstad, A. Rubio and H. Koch, _Phys. Rev. X_ , 2020, 10, 041043.
  * Haugland _et al._ 2021 T. S. Haugland, C. Schäfer, E. Ronca, A. Rubio and H. Koch, _J. Chem. Phys._ , 2021, 154, 094113.
  * DePrince 2022 A. E. DePrince, _J. Chem. Phys._ , 2022, 154, 094112.
  * Liebenthal _et al._ 2022 M. D. Liebenthal, N. Vu and A. E. DePrince, _J. Chem. Phys._ , 2022, 156, 054105.
  * Flick _et al._ 2017 J. Flick, H. Appel, M. Ruggenthaler and A. Rubio, _J. Chem. Theory Comput._ , 2017, 13, 1616–1625.
  * Vindel-Zandbergen _et al._ 2021 P. Vindel-Zandbergen, L. M. Ibele, J.-K. Ha, S. K. Min, B. F. E. Curchod and N. T. Maitra, _J. Chem. Theory Comput._ , 2021.
  * Lee _et al._ 2021 I. S. Lee, J.-K. Ha, D. Han, T. I. Kim, S. W. Moon and S. K. Min, _J. Comput. Chem._ , 2021, 42, 1755–1766.
  * Abedi _et al._ 2010 A. Abedi, N. T. Maitra and E. K. U. Gross, _Phys. Rev. Lett._ , 2010, 105, 123002.
  * Min _et al._ 2017 S. K. Min, F. Agostini, I. Tavernelli and E. K. U. Gross, _J. Phys. Chem. Lett._ , 2017, 8, 3048–3055.
  * Agostini _et al._ 2016 F. Agostini, S. K. Min, A. Abedi and E. K. U. Gross, _J. Chem. Theory Comput._ , 2016, 12, 2127–2143.
  * Min _et al._ 2015 S. K. Min, F. Agostini and E. K. U. Gross, _Phys. Rev. Lett._ , 2015, 115, 073001.
  * Ha and Min 2022 J.-K. Ha and S. K. Min, _J. Chem. Phys._ , 2022, 156, 174109.
  * Vindel-Zandbergen _et al._ 2022 P. Vindel-Zandbergen, S. Matsika and N. T. Maitra, _J. Phys. Chem. Lett._ , 2022, 13, 1785–1790.
  * Hoffmann _et al._ 2018 N. M. Hoffmann, H. Appel, A. Rubio and N. T. Maitra, _Eur. Phys. J. B_ , 2018, 91, 180.
  * Philbin 2014 T. G. Philbin, _Amer. J. Phys._ , 2014, 82, 742–748.
  * Riso _et al._ 2022 R. R. Riso, T. S. Haugland, E. Ronca and H. Koch, _Nat. Commun._ , 2022, 13, 1368.
  * Mandal and Huo 2019 A. Mandal and P. Huo, _J. Phys. Chem. Lett._ , 2019, 10, 5519–5529.
  * Riso _et al._ 2022 R. R. Riso, T. S. Haugland, E. Ronca and H. Koch, _J. Chem. Phys._ , 2022, 156, 234103.
  * Pellegrini _et al._ 2015 C. Pellegrini, J. Flick, I. V. Tokatly, H. Appel and A. Rubio, _Phys. Rev. Lett._ , 2015, 115, 093001.
  * Ruggenthaler _et al._ 2014 M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel, I. V. Tokatly and A. Rubio, _Phys. Rev. A_ , 2014, 90, 012508.
  * Ruggenthaler _et al._ 2018 M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel and A. Rubio, _Nat. Rev. Chem._ , 2018, 2, 1–16.
  * Flick _et al._ 2015 J. Flick, M. Ruggenthaler, H. Appel and A. Rubio, _Proc. Nat. Acad. Sci._ , 2015, 112, 15285–15290.
  * Li _et al._ 2022 T. E. Li, Z. Tao and S. Hammes-Schiffer, _J. Chem. Theory Comput._ , 2022, 18, 2774–2784.
  * Li and Hammes-Schiffer 2023 T. E. Li and S. Hammes-Schiffer, _Electronic Born-Oppenheimer Approximation in Nuclear-Electronic Orbital Dynamics_ , 2023, http://arxiv.org/abs/2301.04076, arXiv:2301.04076 [physics].
  * Mordovina _et al._ 2020 U. Mordovina, C. Bungey, H. Appel, P. J. Knowles, A. Rubio and F. R. Manby, _Phys. Rev. Research_ , 2020, 2, 023262.
  * Fregoni _et al._ 2021 J. Fregoni, T. S. Haugland, S. Pipolo, T. Giovannini, H. Koch and S. Corni, _Nano Lett._ , 2021, 21, 6664–6670.
  * Philbin 2014 T. G. Philbin, _American Journal of Physics_ , 2014, 82, 742–748.
  * White _et al._ 2020 A. F. White, Y. Gao, A. J. Minnich and G. K.-L. Chan, _J. Chem. Phys._ , 2020, 153, 224112.