Theory and modeling of light-matter interactions in chemistry: current and future

Abstract

Light–matter interaction not only plays an instrumental role in characterizing materials' 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 and polariton chemistry. The former utilizes the highly localized photon, plasmonic hot electrons, and local heat to drive chemical reactions. In contrast, polariton chemistry modifies the potential energy curvatures of bare electronic systems, and hence their chemistry, via forming light-matter hybrid states, so-called polaritons. 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.

---

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,¹ 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 interactions, which can be used to monitor, manipulate, and design materials' properties through the complex interplay 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 molecules or molecular materials' chemical and physical properties. However, in the current era of technology,⁴ previously “standard” chemistry, now in a complex electromagnetic (EM) environment, has emergent properties due to multiple new types of couplings between the light and matter degrees of freedom that can be used to control chemistry. Light-matter interactions have been instrumental in many branches of physics, chemistry, materials, and energy science.^3,5,6 In most previous applications, the magnitude and scale of the light–matter interaction fall within the so-called weak coupling regime and can usually be treated at the lowest order in quantum electrodynamics via many-body perturbation theory.⁷ Such treatments are widely used in different types of spectroscopy techniques,^8,9 optoelectronics,^10–14 quantum sensing,¹⁵ quantum information,¹⁶ light harvesting,^17,18 and beyond.

Alternatively, light-matter interactions open multiple new avenues for manipulating matter through novel emerging elementary excitations,¹⁹ 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.²⁰ 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.²¹

On the other hand, when molecules are collectively and resonantly coupled with plasmon excitations, 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 the 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²³ or many molecules are collectively coupled to a single cavity mode.²⁴ In the strong coupling regime, photons and electronic/excitonic excitations in matter become equally important and are strongly coupled on an equally 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. 1)³ 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 photochemistry.²⁵ 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,²⁸ long-range energy transfer,^29–32 Bose-Einstein condensates,^33–35 chemical reactions,^36–57 internal conversion,^58–60 singlet fission,^61–63 manipulation of conical intersection,^64–74etc. Nevertheless, in polariton chemistry, light and matter cannot be treated as separate entities, as the strong coupling between them dresses each of them. Consequently, previously developed quantum chemistry methods for treating the electronic structure of matter alone become invalid. Brand new theoretical and modeling capabilities for describing polariton chemistry are required.

Fig. 1 Light-matter interactions are 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 degrees of freedom (DOFs).

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 that can offer atomistic insights into the problems that are not accessible by experiments alone. 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 DOFs in light-matter interactions. Unfortunately, the strong light–matter interaction is fundamentally a multiscale and multiphysics problem that involves multiple interactions between many DOFs and their complex interplay with environments across different spatial and time scales (Fig. 1). 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 across different coupling regimes in the nonrelativistic limit.⁷⁵ In principle, 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. Hence, we start with a brief introduction to the universal molecular QED theory, followed by its derivation of ab initio methods for plasmon and polariton chemistry (in different coupling regimes), respectively.

2.1 Molecular quantum electrodynamics theory

The minimally coupled Coulomb Hamiltonian governs the nonrelativistic dynamics of matter in an EM environment,with the Coulomb gauge ∇·A = 0. P̂_i = −iħ∇_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. V̂ = V̂_ee + V̂_nn + V̂_en + V_ext contains all of the Coulomb interactions between the electronic and nuclear DOFs and another external potential. The Hamiltonian that describes the EM fields isIt obeys Maxwell's equations of motion and couples with the Schrödinger equation self-consistent through the vector potential A(r_i).

The general quantization of electromagnetic field involves three steps: (1) classical mode description of the electromagnetic field, (2) discretization of classical modes via standard mode decomposition techniques, and (3) quantization via the correspondence principle.⁷⁶ After the quantization of the electromagnetic field, each field (including the vector potential Â) can be expressed as a sum over all possible modes^76,77

where ∇ × u = K × u, α ≡ λK denotes the photon mode with momentum K and polarization λ ∈ {−1,1}, u is the unit vector denoting the direction of the vector potential, and Â_α (Â^†_α) is the mode decomposition coefficients that create (annihilate) the αth radiation mode. Note that and in the Coulomb gauge (i.e., ∇·Â = 0). The photon modes in a given nanophotonic/nanoplasmonic structure can be readily calculated via standard mode decomposition techniques.⁷⁸

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

with the equivalence,Thus, the fields are quantized by the isomorphic association of a quantum mechanical harmonic oscillator to each radiation mode α. Introducing the canonical position and momentum operators as, . The photonic Hamiltonian can now be rewritten as,Then the full light-matter Hamiltonian in Coulomb gauge can now be fully constructed and written as,which is often referred to as the minimal coupling Hamiltonian.

We now aim to apply a unitary transformation on eqn (8) to achieve an expression where the momenta of the molecular DOFs (p̂_i) are decoupled from the vector potential (Â(r_i)). In other words, we will shift the light-matter coupling from momentum fluctuations into displacement fluctuations. This transformation is written as where is the molecular dipole moment.⁷⁹ This transformation Û(D)Ĥ_cÛ^†(D) is nothing but a reduction in matter momentum such that p̂_i − z_iA → p̂_i and a boost in photonic momentum by . Applying the transformation results in the QED Hamiltonian, referred to as the dipole gauge Hamiltonian,

Ĥ = T̂ + V̂ + Ĥ_ep (9)

whereThe second line is reached via canonical transformation between coordinate and momentum operators, p̂ → −ωq̂,q̂ → 1/ωp̂. Here, and g_α ≡ λ_α·D defines the light-matter coupling strength, which depends on the volume of quantized photon V in radiation mode α. The total light-matter Hamiltonian can now be re-partitioned after the unitary transformation as,

Ĥ_PF = Ĥ_M + Ĥ_p + Ĥ_ep + Ĥ_DSE (11)

Here, Ĥ_M = T̂_n + T̂_e + V̂ is the bare molecular Hamiltonian which includes all Coulomb interactions V̂ between electrons and nuclei as well as the kinetic energy operators of both, T̂_e and T̂_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 nm³ or even ų (namely picocavities^23,80,81).

The derivation of eqn (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.^82,83 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.⁸²

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 is particularly strong when the molecule is placed next to a metallic surface or nanostructure.⁸⁴ The plasmon can tune photophysics and photochemistry (i.e., plasmon chemistry) via various possible pathways in such a regime, leading to a promising approach for accelerating and manipulating chemical reactions, which will be the topics of Section 3. In this regime, the light–matter interaction can be treated semiclassically via the classical treatment of the electromagnetic field,where E(r) and B(r) are classical electromagnetic fields. Afterward, the conventional quantum chemistry method can be utilized and extended to describe the plasmon-induced chemical processes. But it should be noted that, even in the weak coupling regime, there is growing interest in exploring the quantum properties of plasmon for applications.^76,85–87 More details will be discussed in Section 3.

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 because the coupling is strong enough to compete with the dissipation. In this regime, polariton formation 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). The theoretical advances towards understanding polariton chemistry and our perspective on future multiscale modeling of polariton chemistry are discussed in Section 4.

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⁸⁸ has led to a wide range of applications, such as plasmon-enhanced molecular spectroscopy,^8,89–92 photovoltaics,^93,94 nanophotonic lasers and amplifiers,^95–98 quantum information,¹⁶ and many others.^18,99,100

Following the ultrafast plasmon excitation, nonradiative plasmon decay leads to the formation of energetic electron–hole pairs (namely hot-carriers),^101–105 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,^101,102 reaching a quasi-thermal equilibrium but with an effective high 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 pathways (more details in Section 3.4). 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.^{20–22,92,106–116} One extraordinary feature of Plasmon chemistry in the realm of chemical reaction mechanisms is its blend of facets from thermochemistry, photochemistry, and photocatalysis.⁹² Unlike traditional reactions that typically focus on singular mechanisms, plasmon chemistry showcases the interaction between different mechanisms in complex electronic and optical environments. At the same time, the mixture of different mechanisms makes it very challenging to comprehend the underlying principles of plasmon chemistry, making it considerably more complicated than transition chemistry. This unique interplay often results in a diverse distribution of reactive zones on substrates.⁹² To fully understand plasmon chemistry's nuances, an interdisciplinary approach is essential, which should consider multi-scale processes, the current state of the field, and the need for advanced experimental methods.

3.1 Experimental endeavours in plasmon-mediated chemical reactions: the need for theoretical insights

Plasmon chemistry stands at the intersection of innovation and discovery for boosting chemical transformation. However, despite the spreading interest, inherent complexities make it still very challenging for efficiency optimization. By leveraging knowledge from established domains such as plasmon-enhanced spectroscopy, many experimental strategies to bolster plasmon chemistry efficiency can be extrapolated. In this section, we discuss several typical experimental investigations into plasmon chemistry across diverse applications and illustrate the essential role of theoretical/numerical modeling in understanding the underlying mechanisms for further improvement. For a more expansive disquisition on the plasmon-mediated chemical reactions for various applications, enthusiasts are redirected to several recent reviews.^{112,117–119}

Artificial photosynthesis. One of the most enticing applications of plasmon chemistry echoes the wonders of photosynthesis.^18,120–122 The goal is to convert solar energy into useful chemicals, a task intricate due to its kinetic and thermodynamic demands. However, in advancements in the realm of solar conversion, a groundbreaking study has showcased the potential of autonomous plasmonic catalysts capable of water splitting under visible light,¹²⁰ as shown in Fig. 2a, present promising trajectories for renewable energy solutions. This paradigm shift focuses on deriving charge carriers from these plasmonic structures. The research highlights the introduction of a solar water-splitting device engineered from a gold nanorod array. Unlike conventional methods, this device leverages the power of “hot electrons” initiated by stimulating surface plasmons within the gold nanostructure. A standout feature is each nanorod's autonomous operation, eliminating the need for external wiring and its impressive capability to produce H₂ molecules under sunlight. The device's resilience and long-term stability further cement its potential. While the initial solar-to-hydrogen efficiency aligns with that of early semiconductor-based water splitters, it remains on the lower side for real-world applications. Nevertheless, the research suggests prospective structural enhancements, emphasizing the modification of the nanorod configuration, that could potentially escalate its efficiency. An important note for consideration, despite the current efficiency metrics, is the device's remarkable longevity, which overtakes some of the top-performing semiconductor-based alternatives. This study marks a pivotal step towards the innovative utilization of plasmonic devices in solar conversion.

Fig. 2 Experimental examples: (a) structure and mechanism of operation of the autonomous plasmonic solar water splitter. Adapted with permission from ref. 120. Copyright 2013 Springer Nature Limited. (b) Light-driven CO₂ reduction with single atomic site antenna-reactor plasmonic photocatalysts. Adapted with permission from ref. 123. Copyright 2020 the author(s), under exclusive license to Springer Nature Limited. (c) Schematic diagram of the proposed mechanism for the oxidation of PTAP to DMAB. Hot electrons generated from the plasmon transfer to the adsorbed O₂ molecules, generating ³O₂ that participated in the PATP oxidation to DMAB. Adapted with permission from ref. 124. Copyright 2015 John Wiley and Sons. (d) Schematic of the three-electrode electrochemical system for distinguishing thermal from hot electron effects by utilizing photoelectrochemical characterization under the chopped light. Adapted from ref. 125 under Creative Commons CC BY.

CO₂ reduction. Given the increasing global CO₂ levels, plasmon chemistry offers a promising solution. Plasmonic structures have shown potential in boosting the photo-reduction of CO₂,^{119,123,126,127} which is crucial for environmental cleanup and advancing green chemistry. In a new development, researchers have presented a plasmonic photocatalyst that blends a copper (Cu) nanoparticle “antenna” with individual atomic ruthenium (Ru) sites (Fig. 2b). This design allows methane dry reforming to happen at room temperature, using the energy from light. Unlike traditional thermocatalytic methods, this photocatalyst works effectively with light at normal conditions and has notable stability and selectivity in its actions. This difference from standard thermally-driven reactions is due to the creation of hot carriers, which improve the rate of carbon-hydrogen activation on Ru sites and speed up hydrogen release. It should also be noted that, different from conventional plasmonic structure, the design of antenna—reactor systems utilize plasmonic metal (the antenna) to collect and concentrate visible light energy and transfer that energy to a catalytic metal (that is, the reactor) to drive a chemical reaction.^128–130 Because the reactors were not in direct contact with the plasmonic nanoparticles, it was argued that the energy transfer resulting in increased reaction rates could only take place via a field effect where the field from the plasmonic metal was felt by the catalytic reactors resulting in the excitation of charge carriers in reactors.

Organic transformation. Beyond the catalysis of inorganic reactions, as previously addressed, plasmonic NPs have showcased commendable performance in catalyzing organic transformations. These hold potential applications, spanning from the synthesis of basic commodity molecules to intricate pharmaceutical compounds. This approach boasts benefits such as superior selectivity, accelerated reaction rates, and comparatively lenient reaction conditions. An illustrative instance is the oxidation of p-aminothiophenol (PATP) or p-nitrothiophenol (PNTP) to dimercaptoazobenzene (DMAB). This reaction has been the focal point of myriad studies in recent years.^{118,124,131–133} The hypothesized microscopic mechanism posits that plasmon-excited electrons transfer from plasmonic nanostructures to O₂ molecules, resulting in the emergence of activated oxygen species, which, in turn, catalyze the formation of DMAB, as shown in Fig. 2(c). Such a mechanism found experimental validation through experiments toggling specific conditions. The obvious evidence of PATP's transformation into DMAB, when air permeates the PATP and plasmonic nanostructured system, is corroborated by the ascendant Raman modes of DMAB over time. In contrast, in the absence of O₂ or electron transfer (attained by coating the plasmonic NPs with an insulator like SiO₂), merely the Raman signals of PATP are.¹³³ In addition, it was found that the plasmon-induced hot holes can trigger the oxidation of PATP, even in the absence of oxygen.¹³⁴ Notably, despite the absence of a tangible PATP transformation on gold nanoparticles due to inefficient charge transfer from the plasmonic nanostructure directly to the adsorbates, PATP oxidation has been witnessed either in the presence or absence of O₂, specifically with Au nanostructures fortified with a TiO₂ shell. These reactions outperform their counterparts in sheer gold nanoparticle systems in terms of efficiency.¹¹⁸ Fascinatingly, it was unveiled that achieving the selective oxidation of PATP to either DMAB or PNTP is possible, depending on whether the plasmonic NPs are reinforced by TiO₂ and the presence of UV-illumination.¹²⁴ Under UV exposure on the bare plasmonic NPs, PATP's oxidation to DMAB is observed. In contrast, there is no catalytic effect with the TiO₂ support. However, after applying UV light exposure, a transition to PNTP is detected with the TiO₂ support. For an expansive disquisition on the plasmon-mediated organic transformation, enthusiasts are redirected to a contemporary review in ref. 117.

Experimental limitations. Though there are many demonstrated applications of plasmon chemistry as a promising approach to accelerate or manipulate chemical processes, the underlying mechanisms are still unclear. The advances in experimental techniques have played a pivotal role in elucidating plasmon chemistry mechanisms. Techniques, such as transient absorption spectroscopy, offer insights into the domain of plasmonic hot-carriers.¹³⁵ Nonetheless, the vast expanse of these reactions remains to be traversed for a holistic comprehension. Traversing the landscape of plasmon chemistry requires a comprehensive understanding of multiple coherent and dissipative processes, including the enhancement of electromagnetic near-fields, local heating effects, and charge-carrier excitation/transfer. These processes occur across varied scales, from femtoseconds to nanoseconds. Since all these effects can drive chemical processes simultaneously, the individual contributions of these elements to overarching reactions are challenging to discern. On the other hand, decoupling individual contributions is vital for designing better plasmon catalysts that synchronize the orchestration of these effects to improve efficiency. In particular, distinguishing the thermal effect from the hot electron effect has become an important topic and is the subject of many debates. Despite the complexity of distinguishing these effects, there are growing efforts to solve these key questions via various designs.^{125,136–140} One of these efforts utilized the unique photoelectrochemical behavior of a plasmonic Au nanoelectrode array as shown in Fig. 2d.¹²⁵ The plasmonic photocurrent can be intensified at negative and positive potentials, with its direction determined by the selected potential. This allows the electrode to favor either reduction or oxidation reactions. The photocurrent is composed of a fast photoelectronic component and a slower photothermal component. The photoelectronic current aligns with the plasmon absorption spectrum and intensifies with light intensity, whereas the photothermal current displays a linear relationship only within specific light intensities. Consequently, this proposed design can differentiate between the photoelectronic and photothermal effects. However, the photoelectronic response measured in ref. 123 happens at a time scale several orders of magnitude larger than the hot-electron lifetime and thermal conduction time scales. More advanced time-resolved techniques are required to provide more comprehensive evidence.

Overall, the distinctiveness of plasmon chemistry is rooted in the interplay between molecules, incident photons, and plasmonic nanostructures. Though many recent breakthroughs have been made to demonstrate plasmon chemistry's application in various areas, further improvements require a deep microscopic understanding of the synergy of multiple effects at the electron level. Such understanding will play an indispensable role in the strategic design and fabrication of plasmonic nanostructures, optimization of surface or interface intermediaries, and their harmonized integration. To elucidate the intricate mechanism that governs plasmon chemistry, comprehensive investigations encompassing various time, spatial, and energy scales are essential, which cannot be achieved by experiments alone. It's crucial that advanced experimental methods, boasting high space-time resolutions, are synergized with microscopic theoretical modeling. Specifically, theoretical models capable of addressing all potential mechanisms on equal footing provide invaluable perspectives into quandaries that remain elusive to experimental endeavors alone. Intuitively, the convenience of modeling allows for the toggling of individual effects to scrutinize their ramifications on holistic reactions. The real theoretical challenge, however, lies in considering these multifaceted effects equitably, spanning varied time and length scales. In the following sections, we will review our theoretical efforts toward simulating all conceivable mechanisms on an equal footing, with an aspiration to understand plasmon chemistry at the electron level.

3.2 Ab initio method for quantum plasmonics

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⁸⁷), and quantum theory for plasmon (quantum plasmonics) is required in nanoscale.

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.⁸⁵ 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.^141–149 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¹⁵⁰ 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 variational principle can reach the ground state energy, though the exact functional is an unsolved problem,

Where 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.¹⁵¹ 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 δV(r) is given by

where χ(r,r′) 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,¹⁵² δρ(r) is also the induced density of the KS system, but due to a perturbationχ_KS here is the linear response of the KS electrons, which can be trivially evaluated from the KS orbitals and energies.

δV_KS(r) = δV_ext(r) + V_H[δρ(r)] + δV_xc(r). (18)

The induced XC potential is given by where f_xc is the dynamical XC kernel ,¹⁵³Hence, a linear equation for the induced density can be derived from eqn (16)–(19),

[1 − χ_KS(ω)f_hxc(ω)]χ(ω)δV_ext(ω) = χ_KS(ω)δV_ext(ω). (20)

Where . Casting eqn (20) into the matrix of coupled KS single excitations, it is possible to calculate the excitation energies ω_s of the system (the poles of the response function) and transition densities (and oscillator strengths). Consequently, the response function χ(ω) can be rewritten aswhere and η represents a positive infinitesimal. The factor of 2 comes from summation over spin indices. The excitation energies in the system are denoted as ω_s, which are calculated from the Casida method.¹⁵³ρ_s(r) is the transition density, which can be expanded in terms of electronic transitions between occupied state i to unoccupied states a,¹⁵⁴where ε_i are the KS eigenvalues, and the corresponding molecular orbitals are ψ_i. And X^s_ia are the Casida transition coefficient from the occupied i state to the unoccupied a state of the sth 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.^155,156

3.3 Hot electron generation and relaxation

Once the plasmon excitation is obtained, the electron–plasmon interaction can be described by the Hamiltonianwhere V_eff(r) is the effective potential induced by the excitation, which can be calculated from the transition density ρ_s(r) by following Lundquvist's approach.^154,155 The polarizability of the NP upon external excitation is given by eqn (21). And the induced potential is given by,¹⁵⁴

V_eff(r) = ε⁻¹(r,r′)δV_ext(r′), (24)

where the dielectric function ε⁻¹(r,r′) is given byThus, within the second quantization, electron-excitation Hamiltonian Ĥ_eps readswhere 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 ESI for details)where is the oscillator vector and is the oscillator strength. is the Hartree potential induced by the transition density. Detailed derivation can be found in the ESI. Above equation actually includes both plasmon excitation (ζ_s(ω)E·f^s), screening effect and hot electron–hole pairs generation (〈ψ_i(r₁)|V^s_H(r₁)|ψ_j(r₁)〉). Neglecting the screening effect would result in a much larger electron–plasmon coupling matrix and shorter HC lifetime distribution.¹⁵⁷ 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. Eqn (26) 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 eqn (26) 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,^101,105

γ_ex refers to the linewidth of the excitation. ω is the excitation energy. δ(ω − ω_s) describes the generalized photoexcitation. Hence, eqn (28) describes the photoexcitation and HC generation. When the energy of the external field ω matches the plasmon energy, eqn (28) 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 ρ(t). The equation of motion (EOM) of ρ is subject to the following quantum Liouville Von-Neumann equation,¹⁵⁸

The above equation's first part on the right-hand side (RHS) describes the coherent evolution of the density matrix between different states. The second and third parts on the RHS of eqn (29) represent the dissipations induced by the electron–electron and electron–phonon scattering effects, respectively. Finally, the last term on the RHS of eqn (29) describes 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),¹⁵⁹. Where the dissipation matrix can be written in terms of Green's functions G and self-energies Σ,^101,159i.e.,Hence, the key point is to develop approximations to the self-energies and efficient numerical methods to compute the . The detailed formalisms of for electron–electron and electron–phonon scatterings can be found in ref. 101. But it should be noted that the approximation used for electron–electron scattering should conserve both particles and energy.

3.4 Plasmon-mediated chemical reactivities and theoretical challenges

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 (Fig. 3). 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.^111,113,160 The second mechanism involves 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.^22,112 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,^{115,161–163} 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.^114,164,165

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

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. For example, The chemical decomposition of methylene blue molecules on the surface of plasmonic NPs was reported by different research groups,^{115,160,161,166,167} direct hot electron transfer mechanism was proposed in ref. 115 and 161 and indirect hot electron transfer is proposed in ref. 167. The same situation also applies to the plasmon-mediated reduction of PNTP, where both the nonthermal and local heating effects were proposed.^168,169 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. 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.¹⁰⁸

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.¹¹⁰

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

4. Inclusion of electron–electron scattering that leads to the redistribution of hot electrons.^104,170

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.5 Nonadiabatic simulation of plasmon-mediated chemical reactivities

Ehrenfest dynamics and Surface Hopping. With the Born–Oppenhermier approximation, the electronic wave function Θ(r,R) is expanded on the basis of adiabatic BO states, which depend on the electronic coordinates r and the nuclear coordinates R(t) according toHere N_st is the total number of adiabatic electronic states, |ϕ_n(r,R(t))〉 is the adiabatic electronic wavefunction of state n and c_n(t) are the time-dependent complex expansion coefficients. (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,¹⁷¹where M_A is the mass of Ath atom. E(R) is the potential energy surface (PES), which can be an averaged one within the mean-field Ehrenfest dynamics or a single PES of a certain state E_k(R) within the Surface hopping (SH) framework.^171–173

By substituting eqn (31) 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^171,174

Here the orthogonal condition of adiabatic states 〈Ψ_n|Ψ_m〉 = δ_nm is used. d_nm = 〈Ψ_n|∇_R|Ψ_m〉 is the nonadiabatic derivative coupling term (or nonadiabatic coupling vector, NACR). A key variable in eqn (33) is the time-derivative nonadiabatic coupling scalar (NACT) between two adiabatic stateswhich 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 eqn (33) 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 by¹⁷⁵

where , with Δt and δt correspond to the time steps for evolving motions of nuclei and electrons in eqn (32) and (33), respectively. The chosen value of δ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 δt will be further refined. defines the time-dependent density matrix elements, and . Note that g_n→m = −g_m→n and g_n→n = 0 since d_nm are antisymmetric. Hopping between adiabatic states is determined stochastically by comparing g_n→m to a random number ξ(ξ ∈ (0,1)). A hop from state n to state m is performed ifwhere states are assumed to be ordered with increasing transition energy. On the other hand, the system remains in state n when . If g_n→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 ref. 173 and 176. 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.¹⁷³ To address these challenges and provide accurate dynamics, a multiconfigurational Ehrenfest (MCE) dynamics approach¹⁷⁷ and ab initio multiple cloning (AIMC)^178–181 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 |Ψ〉 is expressed in the trajectory-guided Gaussian basis functions (TBF) representation (|Ψ_n〉),

And each configuration (or TBF) is described by the product of nuclear and electronic parts,Where |ϕⁿ_I(r,R(t))〉 is the adiabatic state of configuration n. |χ_n〉 are Gaussian nuclear basis functions

The couplings between TBFs in the MCE approach are described by the EOM of c_n(t), which can be readily obtained by substituting eqn (37) into the Schrödinger equation:

whereand the overlap S_mn isThe nuclear part of eqn (41) can be obtained analytically,While the electronic part (or the potential energy matrix elements) are approximated by¹⁸¹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)¹⁸² 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^178–180 This is done by creating one of the clones |ψ_n1〉 in a pure state and the other clone |ψ_n2〉, which includes contributions from all other electronic states. The corresponding MCE amplitudes {c_n1,c_n2} are adjusted to conserve the original wavefnction.^180,181

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.

3.6 Simulation of plasmon-mediated phenomena

3.6.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.^{85,105,109,183–185} 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,¹⁰¹ as well as the atomic-scale mechanism of plasmonic hot-carrier-mediated chemical processes such as H₂ dissociation.¹⁰⁹ Our numerical simulations showed that after photoexcitation, hot carriers transfer to the antibonding state of the H₂ molecule from the nanoparticle, leading to a repulsive-potential-energy surface and H₂ dissociation (Fig. 4(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 Fig. 4(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.

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

3.6.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 H₂ molecule adsorbed on an Au6 cluster (H2@Au₆),¹¹⁰ in order to explore key pathways in LSPR-promoted chemical reactions. Despite the model system being 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 Au₆ cluster supports the adsorption of the H₂ 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 H₂ 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 Au₆ cluster, which can be likened to HE states in metal nanoclusters; (2) the other group of excited states has the antibonding σ* characteristics of the H₂ adsorbate due to the hybridization between H₂ and Au₆ 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 H₂ 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 H₂ dissociation on Au clusters.

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.¹⁸⁶ 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.6.3 TDDFTB-based NAMD simulations. As shown in Fig. 5, 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),¹⁸⁷ might be needed for the efficient electronic structure calculations. To this aim, We recently developed an efficient NAMD method by combining the TSH algorithm and LR-TDDFTB method.^{175,188–190} LR-TDDFTB is the tight-binding extension of TDDFT,^188,191 which has been successfully applied to investigate the optical properties of plasmonic NPs.^192,193 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 Au₂₀ cluster. Our simulations show that Au₂₀ can support plasmon-like excitation, which includes the superposition of multiple single-particle excitation components.

Fig. 5 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 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 Au₂₀ 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 Au₂₀ cluster. Furthermore, our simulations show that the lifetime of the phonon-induced plasmon dephasing is ∼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₂₀ cluster.¹⁰⁸ 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₂₀, 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

Since there are several recent experimental and theoretical reviews on polariton chemistry,^25,194–198 only a brief introduction to the recent theoretical development of polariton chemistry will be reviewed in this perspective. Instead, we will focus on our perspective on the theoretical challenges and ongoing and future development of quantum many-body and multiscale methods toward solving the mysteries of polariton chemistry and beyond.

4.1 Introduction and background

As introduced in Section 2, in the strong light-matter coupling limit (strong enough to compete with individual DOF's dissipation rate), the formed quasiparticles (polaritons) can alter the potential energy landscape and consequently manipulate the chemical processes. 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.^53,199–207 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. Like the situation in plasmon chemistry, the challenges in understanding polariton chemistry still lie in the multiple coherent and dissipative processes across various time/length scales.

1. The first theoretical hurdle is that of the Hamiltonian itself. For the Pauli–Fierz Hamiltonian (eqn (12)), the light–matter interaction requires explicit knowledge of the molecular dipole operator 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.²⁰⁸ (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 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.

2. Generally, the strong coupling relies on a larger number of molecules collectively coupled to a cavity, introducing significant challenges to computing the hybridized polaritonic states of many (10⁶ or even larger) molecules in the cavity.

3. While plasmonic cavities are able to achieve strong coupling in a few (or even single) molecule limits, these cavities are highly heterogeneous and dispersive. Consequently, multiple models, dispersion of the cavities, coupling beyond the dipole approximation, and the inherent strong dissipation should be taken into account, introducing different computational complexity. Besides, the molecules inside the nanoplasmonic cavities can also interact with the ones outside the cavity. For example, a recent 2DES investigation of J-aggregates in the plasmonic cavity reveals rich photophysical dynamics between the polaritonic states inside the cavity and uncoupled exciton states outside the cavity,²⁰⁹ which indicates that the molecules outside the cavity cannot be trivially ignored.

4. The damping of the polariton state is usually ignored or introduced by an empirical lifetime of cavity photon (cavity leakage). In principle, cavity leakage depends on the photonic density of states, which can be computed from the first principles. Besides, the light-matter coupling could affect the photonic DOF and, in return, the photon lifetime. Hence, one open question is to what extent the strong coupling changes the photon lifetime, which requires a multiscale method to address the feedback of molecular systems on the cavities and first-principles calculations of cavity leakage.

5. Polariton-mediated phenomena are fundamentally kinetic procedures involving the dynamic interplay between electronic, nuclear, and photonic DOFs at various time and length scales. Indeed, the accurate simulation of the quantum dynamics of a molecule itself is already a challenging task.^171,210,211 The inclusion of photonic DOFs in polariton dynamics adds further complexity. To accurately simulate polariton dynamics, efficient and precise calculations of the gradients and derivative coupling of polariton states are essential. Despite the recent advancements in accurately computing gradients of polariton states,¹⁹⁴ more efforts are needed for large-scale ab initio simulations of polariton dynamics of many molecules.

6. Polariton-mediated phenomena are fundamentally kinetic procedures involving the dynamic interplay between electronic, nuclear, and photonic DOFs at various time and length scales. Indeed, the accurate simulation of the quantum dynamics of a molecule itself is already a challenging task. The inclusion of photonic DOFs in polariton dynamics adds further complexity. To accurately simulate polariton dynamics, efficient and precise calculations of polariton states' gradients and derivative coupling are essential. Despite the recent advancements in accurately computing gradients of polariton states,¹⁹⁴ more efforts are needed for large-scale ab initio simulations of polariton dynamics in multiple molecules.

Hence, simulating polariton chemistry is non-trivial, and state-of-the-art theoretical models and efficient numerical methods are crucial for understanding the polariton dynamics of many molecules.

4.3 Common approximations toward historic quantum optics models

The PF Hamiltonian (as written in eqn (12)) 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.^212–218 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〉 and |e〉 are included while only including the vacuum and singly excited Fock states, |0〉 and |1〉 in a single-mode cavity α = 0. The total basis for this simple model is then confined to {|g,0〉, |g,1〉, |e,0〉, |e,1〉}. Note here that another approximation is that there is no permanent dipole in the ground |g〉 or excited |e〉 electronic states, i.e., D_gg = D_ee = 0. The most commonly used Hamiltonian for modeling ab initio polaritons is the Jaynes–Cummings (JC) Hamiltonian Ĥ_JC, which can be written as,where [small sigma, Greek, circumflex]^† ([small sigma, Greek, circumflex]) is the creation (annihilation) operator for the molecule excitation between the ground g and excited e states. And the Tavis–Cummings (TC)²¹⁹ model is the many molecules generalization of the JC model. (4) Besides, the feedback of molecular systems on the cavity modes (via the current density) is usually ignored. Here, two approximations have been made to derive JC/TC models: (I) the rotating wave approximation (RWA) – which is to say, neglecting the highly oscillatory [small sigma, Greek, circumflex]^†â^†_α and [small sigma, Greek, circumflex]â_α terms in the light–matter interaction – and (II) neglecting the DSE Ĥ_DSE. This Hamiltonian is valid at ultra-low coupling strengths where the splitting between the one-photon-dressed ground state |g,1〉 and the excited state with zero photons |e,0〉 exhibit linear splitting (e.g. Rabi splitting) with an increase in the coupling strength λ. The other two basis states |g,0〉 and |e,1〉 are completely decoupled from the interaction.

Two other common approximations to the PF Hamiltonian are the Rabi model Ĥ_Rabi after explicitly dropping the dipole self-energy (DSE) term and the RWA Ĥ_RWA. These two approximate Hamiltonians can be written as,

where the Rabi Hamiltonian Ĥ_Rabi = Ĥ_PF − Ĥ_DSE and the RWA Hamiltonian applies the RWA to the PF Hamiltonian. Here, [small sigma, Greek, circumflex]_ge = |g〉〈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 ref. 194, 212, 220 and 221.

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, D̂, in realistic molecules is far from sparse with its square leading to further complications.^57,194 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 Section 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,²¹² and in general poor results compared to benchmarks, especially at large light-matter coupling strengths when using a small electronic basis.^222–224 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 Section 4.4.2).^{57,222,224–229} 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., ∇_RD̂ and ∇_RD̂²).²¹³

4.4 Theoretical modeling of polariton chemistry

4.4.1 Cavity Born–Oppenheimer approximation. Within the Born–Oppenheimer (BO) approximation,²³⁰ the total electronic–photonic–nuclear can be factorized as,

Φ(r,R,q_α) = χ(R)Ψ(r,q_α;R), (48)

where χ(R) and Ψ(r,q_α;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,where ψ_μ(r,q_α;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 ψ_μ(r,q_α;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 ψ_μ(r,q_α;R). First, a brief description of a direct diagonalization approach (Section 4.4.2) with Hamiltonians parameterized with information from standard electronic structure, while the following sections (Sections 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 for the QED Hamiltonian (e.g., QED-HF, QED-DFT, etc.).

4.4.2 Direct diagonalization. The polaritonic state can be readily computed via exact diagonalization of polariton Hamiltonian (either JC, TC, Rabi, RWA, or rigorous PF) in a certain basis set. For instance, the matrix elements of the single-mode PF Hamiltonian in the widely used adiabatic-Fock basis can be written as,where 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 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.

Except for the adiabatic-Fock basis, other options for expanding the JC and PF Hamiltonians exist, such as the coherent states^226,231,232 or polarized Fock states.²³³ 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.^57,212 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.
scQED Hartree–Fock and density functional theories. 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.^{224,228,229,232,234} A useful basis, referred to as coherent states (CS), can be performed such that the light–matter interaction part of the PF Hamiltonian can be shifted away.

For a given bare electronic ground state wavefunction (|HF〉), the photonic Hamiltonian can be obtained by integrating out the electronic DOF,

where 〈⋯〉_HF = 〈HF|⋯|HF〉 is the HF ground state expectation value of the electronic subsystem. The photonic Hamiltonian can be trivially diagonalized by introducing the CS transformation Û(z)Ĥ_pÛ^†(z),where z = {z_α} is a vector of complex numbers specific to each cavity mode,Thus, considering a QED-HF ansatz as , the QED-HF energy and corresponding Fock matrix can be derived. Compared to the HF theory for bare electrons, the QED-HF method introduces DSE-mediated one-body and two-body integrals.¹⁹⁵

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. 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 effects^235–238 or working with electron-only exchange–correlation functionals using the coherent state basis for the photonic DOFs.^212,225,226 In any case, the resulting ground state is uncorrelated by the 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.^239,240 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-workers²²² using the QEDFT (or scQED-DFT in the notation of this0 work) method in the language of Casida and further used by the groups of Shao²¹² and DePrince.²²⁵ 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 frequency.²²⁴

As per usual, and following the notation of ref. 225, 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 (eqn (12)).

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. 212, the orbitals are not relaxed while in the ref. 222 and 225 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 and photon-mediated electron–electron correlations, and the exchange–correlation function for electron–photon interaction is unknown. Consequently, the QED counterpart of coupled-cluster theory (QED-CC) is proposed.^{36,37,227–229,241,242} Similar to the conventional CC theory, QED-CC employs an exponential wavefunction Ansatz to derive the ground state,

|Ψ_CC〉 = e^T̂|Φ₀〉, (54)

where |Φ₀〉 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., |Φ₀〉 = |φ₀〉 ⊗ |0〉. T̂ is defined in excitation configurations. Within QED-CC theory, the generalized excitation operator includes three components, T̂ = T̂_e + T̂_p + T̂_ep, including electronic (T̂_e), photonic (T̂_p), and coupled electronic–photonic (T̂_ep) excitations,where N_e and N_f are the numbers of electrons and photon Fock states, respectively. and Ê^a_i = â^†_aâ_i + h.c. are the μ-body and single-body excitation operators, respectively. is the n-body photonic excitation. The parameters (t_μ,γ_n,χ_μ,n) are the cluster amplitudes. By projecting CC wavefunction into a set of orthogonal excited configurations ({|μ〉}),

|μ〉 = [small mu, Greek, circumflex]|Φ₀〉, (58)

where μ ∈ {[small tau, Greek, circumflex]_μ,B̂_n,[small tau, Greek, circumflex]_μB̂_n} is the cluster operator, The cluster amplitudes can be determined from the projected equation by following the standard CC procedure,

Ω_μ ≡ 〈μ|e^−T̂Ĥe^T̂|Φ₀〉 = 0 (59)

Though the scQED-CC theory leads to the FCI solution if no truncation on the excitation operator T̂ is applied, the excitation operator is usually truncated at the doubles (CCSD) level in order to trade-off between accuracy and computational efficiency. Some tests have been performed on the level of truncation in the photonic excitations (up to 10) in the CC operator for model systems.²⁴¹ However, the scaling of scQED-CCSD is in general, where M_p is the number of photon modes. Such a scaling makes it 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 coherence state²⁴³ may be used to reduce the photon basis in the scQED-CC calculations.

It should be noted that CC methods exhibit size extensivity even in their 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, making it an ideal candidate for computing the polaritonic eigenstates of many (non-interacting) molecules.

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²⁴⁴

|R̂〉 = R̂|Φ_cc〉 = R̂e^T̂|Φ₀,0〉. (60)

where R̂ also include the electronic (R̂_e), photonic (R̂_p), and coupled electronic–photonic (R̂_ep) excitations. Because the excitation operator, R̂, commutes with the excitation operators in T̂, solving this eigenvalue problem is equivalent to finding a right eigenvector of the similarity transformed Hamiltonian,

〈μ|H̄R̂ⁿ|Φ₀,0〉 = E_nR̂ⁿ_μ. (61)

Here, E_n is the energy of the nth excited state, and μ indexes an element of the excitation operator R̂. The excitation operator, R̂, can be chosen to access charged or neutral excitations.

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. 6 is an example of the Li–H bond dissociation curves in a single-mode optical cavity. Fig. 6a 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 (λ = 0). The excited states inside (red curves) and outside (black circles) of the cavity are shown at zero coupling strength (λ = 0), where two multiplicities (singlet |S₁〉 ≡ |S₁,0〉 and triplet |T₁〉 ≡ |T₁,0〉) of electronic states are shown as well as a single cavity state (|ω〉 ≡ |S₀,1〉). Here, we have used the notation |S_J〉 ⊗ |n〉 = |S_J,n〉 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 λ = 0.088 a.u., the cavity state with one photon |ω〉 couples strongly with the singlet excited electronic state with zero photons |S₁〉 while the excited electronic triplet state |T₁〉 is negligibly affected. The so-called Rabi splitting appears at the degenerate point between the singlet electronic state and the cavity state, R_Li–H ≈ 3.75, forming the upper (|P_U〉) and lower (|P_L〉) polaritonic states with mixed electronic and photonic character.

Fig. 6 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 λ = 0 a.u. (b) EOM-CCSD (black circles) triplet |T₁〉 and singlet |S₁〉 states as well as the low-lying states (|T₁〉, |ω〉, and |S₁〉) calculated at scQED-EOM-CCSD with light-matter coupling strength λ = 0 a.u. The plot clearly shows the photon line (|ω〉) and perfect overlap of the S₁ states. (c) Upper (|P_U〉) and lower (|P_L〉) polariton states with λ = 0.088 a.u. with the triplet state |T₁〉 uncoupled.

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. 6). The efficiency of Numpy einsum (CPU-accessible) and Torch einsum (GPU-accessible) functions are shown in Fig. 7 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. 7 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 × −20 × 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 nm³ or even ų (picocavities).²³ 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 θ (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.^214,245 Hence, a general method that can compute the multiple spatial-dependent cavity eigenmodes λ_nk(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),²⁴⁶ 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 that 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):

Where J is the paramagnetic current density due to the molecular dipoles. 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 Fig. 8.

Fig. 8 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,

where Green's function is the formal solution of Helmholtz's equationwith linear media ε(r). The Green's function in eqn (65) 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.^247,248 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 J) and macroscopic solvent (represented by a parameterized dielectric function ε(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²⁴⁶ and method of moments,²⁴⁹ providing the photon mode structures that are coupled with QED electronic structure solvers. By computing 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,²¹ quantum plasmonics,⁸⁶ quantum information transduction,²⁵⁰ and photonic-microelectronics integration and polaritonic devices.²⁵¹ The integration with open quantum system method, particularly time-dependent quantum transport theories,^252–254 will also open the door to simulating light-driven quantum transport phenomena^255,256 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 fundamentally 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) and the 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: LANLECF7). 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

1. G. R. Fleming and M. A. Ratner, Phys. Today, 2008, 61, 28–33.
2. A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta and F. Nori, Nat. Rev. Phys., 2019, 1, 19–40.
3. N. Rivera and I. Kaminer, Nat. Rev. Phys., 2020, 2, 538–561.
4. P. W. Anderson, Science, 1972, 177, 393–396.
5. L. Novotny and B. Hecht, Annu. Rev. Phys. Chem., 2006, 57, 303–331.
6. J. Weiner and F. Nunes, Light-Matter Interaction: Physics and Engineering at the Nanoscale, Oxford University Press, 2012.
7. G. Mahan, Many-Particle Physics, Springer, 2000.
8. P. L. Stiles, J. A. Dieringer, N. C. Shah and R. P. Van Duyne, Annu. Rev. Anal. Chem., 2008, 1, 601–626.
9. S. Mukamel, Annu. Rev. Phys. Chem., 2000, 51, 691–729.
10. Y. Zhang, L. Meng, C. Yam and G. Chen, J. Phys. Chem. Lett., 2014, 5, 1272–1277.
11. C. Yam, L. Meng, Y. Zhang and G. Chen, Chem. Soc. Rev., 2015, 44, 1763–1776.
12. L. Meng, C. Yam, Y. Zhang, R. Wang and G. Chen, J. Phys. Chem. Lett., 2015, 6, 4410–4416.
13. L. Meng, Y. Zhang and C. Yam, J. Phys. Chem. Lett., 2017, 8, 571–575.
14. R. Wang, Y. Zhang, G. H. Chen and C. Y. Yam, Prog. Electromagn. Res., 2015, 154, 163–170.
15. C. L. Degen, F. Reinhard and P. Cappellaro, Rev. Mod. Phys., 2017, 89, 035002.
16. F. Flamini, N. Spagnolo and F. Sciarrino, Rep. Prog. Phys., 2018, 82, 016001.
17. L. Zhou, Y. Zhang, Z. Zhuo, A. J. Neukirch and S. Tretiak, J. Phys. Chem. Lett., 2018, 9, 6915–6920.
18. T. Mirkovic, E. E. Ostroumov, J. M. Anna, R. van Grondelle, Govindjee and G. D. Scholes, Chem. Rev., 2017, 117, 249–293.
19. P. Forn-Díaz, L. Lamata, E. Rico, J. Kono and E. Solano, Rev. Mod. Phys., 2019, 91, 025005.
20. M. L. Brongersma, N. J. Halas and P. Nordlander, Nat. Nanotechnol., 2015, 10, 25–34.
21. C. Zhan, J. Yi, S. Hu, X.-G. Zhang, D.-Y. Wu and Z.-Q. Tian, Nat. Rev. Methods Primers, 2023, 3, 12.
22. E. Kazuma and Y. Kim, Angew. Chem., Int. Ed., 2019, 58, 4800–4808.
23. J. J. Baumberg, Nano Lett., 2022, 22, 5859–5865.
24. 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.
25. T. W. Ebbesen, Acc. Chem. Res., 2016, 49, 2403–2412.
26. S. Kéna-Cohen and S. R. Forrest, Nat. Photonics, 2010, 4, 371–375.
27. J.-W. Kang, B. Song, W. Liu, S.-J. Park, R. Agarwal and C.-H. Cho, Sci. Adv., 2019, 5, eaau9338.
28. 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.
29. 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.
30. K. Georgiou, R. Jayaprakash, A. Othonos and D. G. Lidzey, Angew. Chem., Int. Ed., 2021, 60, 16661–16667.
31. M. Wang, M. Hertzog and K. Börjesson, Nat. Commun., 2021, 12, 1874.
32. 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.
33. M. Dusel, S. Betzold, O. A. Egorov, S. Klembt, J. Ohmer, U. Fischer, S. Höfling and C. Schneider, Nat. Commun., 2020, 11, 2863.
34. A. V. Zasedatelev, A. V. Baranikov, D. Urbonas, F. Scafirimuto, U. Scherf, T. Stöferle, R. F. Mahrt and P. G. Lagoudakis, Nat. Photonics, 2019, 13, 378–383.
35. A. Kavokin, T. C. H. Liew, C. Schneider, P. G. Lagoudakis, S. Klembt and S. Hoefling, Nat. Rev. Phys., 2022, 4, 435–451.
36. 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, arXiv, 2022, preprint, arXiv:2208.06925 [physics] DOI:10.48550/arXiv.2208.06925.
37. F. Pavošević, S. Hammes-Schiffer, A. Rubio and J. Flick, J. Am. Chem. Soc., 2022, 144, 4995–5002.
38. C. Schäfer, J. Flick, E. Ronca, P. Narang and A. Rubio, Nat. Commun., 2022, 13, 7817.
39. C. Schäfer, J. Phys. Chem. Lett., 2022, 13, 6905–6911.
40. R. J. Cave and M. D. Newton, J. Chem. Phys., 1997, 106, 9213–9226.
41. L. A. Martínez-Martínez, R. F. Ribeiro, J. C. González-Angulo and J. Yuen-Zhou, ACS Photonics, 2017, 5, 167–176.
42. P.-Y. Yang and J. Cao, J. Phys. Chem. Lett., 2021, 12, 9531–9538.
43. C. Climent and J. Feist, Phys. Chem. Chem. Phys., 2020, 22, 23545–23552.
44. D. S. Wang, T. Neuman, S. F. Yelin and J. Flick, J. Phys. Chem. Lett., 2022, 13, 3317–3324.
45. M. V. Imperatore, J. B. Asbury and N. C. Giebink, J. Chem. Phys., 2021, 154, 191103.
46. J. Galego, F. J. Garcia-Vidal and J. Feist, Nat. Commun., 2016, 7, 13841.
47. J. A. Campos-Gonzalez-Angulo, R. F. Ribeiro and J. Yuen-Zhou, Nat. Commun., 2019, 10, 4685.
48. J. P. Philbin, Y. Wang, P. Narang and W. Dou, J. Phys. Chem. C, 2022, 126, 14908–14913.
49. S. Efrima and M. Bixon, Chem. Phys. Lett., 1974, 25, 34–37.
50. N. T. Phuc, P. Q. Trung and A. Ishizaki, Sci. Rep., 2020, 10, 7318.
51. E. Davidsson and M. Kowalewski, J. Chem. Phys., 2020, 153, 234304.
52. J. Galego, F. J. Garcia-Vidal and J. Feist, Phys. Rev. Lett., 2017, 119, 136001.
53. L. Mauro, K. Caicedo, G. Jonusauskas and R. Avriller, Phys. Rev. B, 2021, 103, 165412.
54. I. Vurgaftman, B. S. Simpkins, A. D. Dunkelberger and J. C. Owrutsky, J. Phys. Chem. Lett., 2020, 11, 3557–3562.
55. H. Hiura and A. Shalabney, ChemRxiv, 2021, preprint DOI:10.26434/chemrxiv.7234721.v5.
56. H. Hiura, A. Shalabney and J. George, ChemRxiv, 2019, preprint DOI:10.26434/chemrxiv.7234721.v4.
57. B. M. Weight, T. Krauss and P. Huo, J. Phys. Chem. Lett., 2023, 14, 5901–5913.
58. B. F. E. Curchod, C. Rauer, P. Marquetand, L. González and T. J. Martínez, J. Chem. Phys., 2016, 144, 101102.
59. W. Wu, A. E. Sifain, C. A. Delpo and G. D. Scholes, J. Chem. Phys., 2022, 157, 161102.
60. A. G. Avramenko and A. S. Rury, J. Phys. Chem. Lett., 2020, 11, 1013–1021.
61. C. Climent, D. Casanova, J. Feist and F. J. Garcia-Vidal, Cell Rep. Phys. Sci., 2022, 3, 100841.
62. B. Gu and S. Mukamel, J. Phys. Chem. Lett., 2021, 12, 2052–2056.
63. L. A. Martnez-Martnez, M. Du, R. F. Ribeiro, S. Kéna-Cohen and J. Yuen-Zhou, J. Phys. Chem. Lett., 2018, 9, 1951–1957.
64. D. Cho, B. Gu and S. Mukamel, J. Am. Chem. Soc., 2022, 144, 7758–7767.
65. A. Csehi, O. Vendrell, G. J. Halász and Á. Vibók, New J. Phys., 2022, 24, 073022.
66. B. Gu and S. Mukamel, J. Phys. Chem. Lett., 2020, 11, 5555–5562.
67. B. Gu and S. Mukamel, Chem. Sci., 2020, 11, 1290–1298.
68. T. Szidarovszky, G. J. Halász, A. G. Császár, L. S. Cederbaum and Á. Vibók, J. Phys. Chem. Lett., 2018, 9, 6215–6223.
69. C. Fábri, G. J. Halász and Á. Vibók, J. Phys. Chem. Lett., 2022, 13, 1172–1179.
70. A. Natan, M. R. Ware, V. S. Prabhudesai, U. Lev, B. D. Bruner, O. Heber and P. H. Bucksbaum, Phys. Rev. Lett., 2016, 116, 143004.
71. A. Hofmann and R. de Vivie-Riedle, Chem. Phys. Lett., 2001, 346, 299–304.
72. M. H. Farag, A. Mandal and P. Huo, Phys. Chem. Chem. Phys., 2021, 23, 16868–16879.
73. C. Arnold, O. Vendrell, R. Welsch and R. Santra, Phys. Rev. Lett., 2018, 120, 123001.
74. I. S. Ulusoy, J. A. Gomez and O. Vendrell, J. Phys. Chem. A, 2019, 123, 8832–8844.
75. D. Craig and T. Thirunamachandran, Molecular Quantum Electrodynamics: An Introduction to Radiation-molecule Interactions, Dover Publications, 1998.
76. A. Archambault, F. M. C. Marquier, J.-J. Greffet and C. Arnold, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 035411.
77. D. L. Andrews, D. P. Craig and T. Thirunamachandran, Int. Rev. Phys. Chem., 1989, 8, 339–383.
78. A. W. Snyder and J. Love, Optical Waveguide Theory, Springer, 1st edn, 1983.
79. R. G. Woolley, Phys. Rev. Res., 2020, 2, 013206.
80. 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–130.
81. 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.
82. T. Neuman, R. Esteban, D. Casanova, F. J. Garca-Vidal and J. Aizpurua, Nano Lett., 2018, 18, 2358–2364.
83. D. S. Wang, T. Neuman and P. Narang, J. Phys. Chem. C, 2021, 125, 6222–6228.
84. K. L. Kelly, E. Coronado, L. L. Zhao and G. C. Schatz, J. Phys. Chem. B, 2003, 107, 668–677.
85. A. Varas, P. García-González, J. Feist, F. García-Vidal and A. Rubio, Nanophotonics, 2016, 5, 409–426.
86. M. S. Tame, K. R. McEnery, S. K. Ozdemir, J. Lee, S. A. Maier and M. S. Kim, Nat. Phys., 2013, 9, 329–340.
87. D. Marinica, A. Kazansky, P. Nordlander, J. Aizpurua and A. G. Borisov, Nano Lett., 2012, 12, 1333–1339.
88. J. A. Schuller, E. S. Barnard, W. Cai, Y. C. Jun, J. S. White and M. L. Brongersma, Nat. Mater., 2010, 9, 193–204.
89. Jiang, K. Bosnick, M. Maillard and L. Brus, J. Phys. Chem. B, 2003, 107, 9964–9972.
90. L. Brus, Acc. Chem. Res., 2008, 41, 1742–1749.
91. S. Nie and S. R. Emory, Science, 1997, 275, 1102–1106.
92. C. Zhan, X.-J. Chen, J. Yi, J.-F. Li, D.-Y. Wu and Z.-Q. Tian, Nat. Rev. Chem., 2018, 2, 216–230.
93. H. A. Atwater and A. Polman, Nat. Mater., 2010, 9, 205–213.
94. Y. Zhang, C. Yam and G. C. Schatz, J. Phys. Chem. Lett., 2016, 7, 1852–1858.
95. R.-M. Ma, R. F. Oulton, V. J. Sorger and X. Zhang, Laser Photonics Rev., 2013, 7, 1–21.
96. P. Berini and I. De Leon, Nat. Photonics, 2012, 6, 16–24.
97. J. Guan, R. Li, X. G. Juarez, A. D. Sample, Y. Wang, G. C. Schatz and T. W. Odom, Adv. Mater., 2021, 2103262.
98. J. Guan, J.-E. Park, S. Deng, M. J. H. Tan, J. Hu and T. W. Odom, Chem. Rev., 2022, 122, 15177–15203.
99. S. K. Cushing and N. Wu, J. Phys. Chem. Lett., 2016, 7, 666–675.
100. J. Li, S. K. Cushing, F. Meng, T. R. Senty, A. D. Bristow and N. Wu, Nat. Photonics, 2015, 9, 601–607.
101. Y. Zhang, J. Phys. Chem. A, 2021, 125, 9201–9208.
102. C. Clavero, Nat. Photonics, 2014, 8, 95–103.
103. M. Bernardi, J. Mustafa, J. B. Neaton and S. G. Louie, Nat. Commun., 2015, 6, 7044.
104. A. M. Brown, R. Sundararaman, P. Narang, I. William, A. Goddard and H. A. Atwater, ACS Nano, 2016, 10, 957–966.
105. A. Manjavacas, J. G. Liu, V. Kulkarni and P. Nordlander, ACS Nano, 2014, 8, 7630–7638.
106. 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.
107. 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.
108. X. Wu, T. van der Heide, S. Wen, T. Frauenheim, S. Tretiak, C. Yam and Y. Zhang, Chem. Sci., 2023, 14, 4714–4723.
109. Y. Zhang, T. Nelson, S. Tretiak, H. Guo and G. C. Schatz, ACS Nano, 2018, 12, 8415–8422.
110. Q. Wu, L. Zhou, G. C. Schatz, Y. Zhang and H. Guo, J. Am. Chem. Soc., 2020, 142, 13090–13101.
111. E. Kazuma, J. Jung, H. Ueba, M. Trenary and Y. Kim, Science, 2018, 360, 521–526.
112. Y. Zhang, S. He, W. Guo, Y. Hu, J. Huang, J. R. Mulcahy and W. D. Wei, Chem. Rev., 2018, 118, 2927–2954.
113. T. E. Tesema, B. Kafle and T. G. Habteyes, J. Phys. Chem. C, 2019, 123, 8469–8483.
114. J. R. Adleman, D. A. Boyd, D. G. Goodwin and D. Psaltis, Nano Lett., 2009, 9, 4417–4423.
115. C. Boerigter, R. Campana, M. Morabito and S. Linic, Nat. Commun., 2016, 7, 10545.
116. Y. Zhang, T. Nelson and S. Tretiak, Computational Photocatalysis: Modeling of Photophysics and Photochemistry at Interfaces, ACS Publications, 2019, pp. 239–256.
117. A. Gellé, T. Jin, L. de la Garza, G. D. Price, L. V. Besteiro and A. Moores, Chem. Rev., 2020, 120, 986–1041.
118. C. Zhan, Z.-Y. Wang, X.-G. Zhang, X.-J. Chen, Y.-F. Huang, S. Hu, J.-F. Li, D.-Y. Wu, M. Moskovits and Z.-Q. Tian, J. Am. Chem. Soc., 2019, 141, 8053–8057.
119. Y. Li, Y. Lei, D. Li, A. Liu, Z. Zheng, H. Liu, J. Guo, S. Liu, C. Hao and D. He, ACS Catal., 2023, 13, 10177–10204.
120. S. Mubeen, J. Lee, N. Singh, S. Krämer, G. D. Stucky and M. Moskovits, Nat. Nano., 2013, 8, 247–251.
121. C. Xiao, H. Hu, X. Zhang and D. R. MacFarlane, ACS Sustainable Chem. Eng., 2017, 5, 10858–10863.
122. S. Yu and P. K. Jain, ACS Energy Lett., 2019, 4, 2295–2300.
123. 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.
124. J. Wang, R. A. Ando and P. H. C. Camargo, Angew. Chem., Int. Ed., 2015, 54, 6909–6912.
125. C. Zhan, B.-W. Liu, Y.-F. Huang, S. Hu, B. Ren, M. Moskovits and Z.-Q. Tian, Nat. Commun., 2019, 10, 2671.
126. S. Singh, R. Verma, N. Kaul, J. Sa, A. Punjal, S. Prabhu and V. Polshettiwar, Nat. Commun., 2023, 14, 2551.
127. M. Dhiman, A. Maity, A. Das, R. Belgamwar, B. Chalke, Y. Lee, K. Sim, J.-M. Nam and V. Polshettiwar, Chem. Sci., 2019, 10, 6594–6603.
128. U. Aslam, S. Chavez and S. Linic, Nat. Nanotechnol., 2017, 12, 1000–1005.
129. C. Zhang, H. Zhao, L. Zhou, A. E. Schlather, L. Dong, M. J. McClain, D. F. Swearer, P. Nordlander and N. J. Halas, Nano Lett., 2016, 16, 6677–6682.
130. D. F. Swearer, H. Zhao, L. Zhou, C. Zhang, H. Robatjazi, J. M. P. Martirez, C. M. Krauter, S. Yazdi, M. J. McClain, E. Ringe, E. A. Carter, P. Nordlander and N. J. Halas, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 8916–8920.
131. A. G. M. da Silva, T. S. Rodrigues, V. G. Correia, T. V. Alves, R. S. Alves, R. A. Ando, F. R. Ornellas, J. Wang, L. H. Andrade and P. H. C. Camargo, Angew. Chem., Int. Ed., 2016, 55, 7111–7115.
132. Y.-F. Huang, H.-P. Zhu, G.-K. Liu, D.-Y. Wu, B. Ren and Z.-Q. Tian, J. Am. Chem. Soc., 2010, 132, 9244–9246.
133. Y.-F. Huang, M. Zhang, L.-B. Zhao, J.-M. Feng, D.-Y. Wu, B. Ren and Z.-Q. Tian, Angew. Chem., Int. Ed., 2014, 53, 2353–2357.
134. L.-B. Zhao, M. Zhang, Y.-F. Huang, C. T. Williams, D.-Y. Wu, B. Ren and Z.-Q. Tian, J. Phys. Chem. Lett., 2014, 5, 1259–1266.
135. H. Reddy, K. Wang, Z. Kudyshev, L. Zhu, S. Yan, A. Vezzoli, S. J. Higgins, V. Gavini, A. Boltasseva, P. Reddy, V. M. Shalaev and E. Meyhofer, Science, 2020, 369, 423–426.
136. Y. Yu, V. Sundaresan and K. A. Willets, J. Phys. Chem. C, 2018, 122, 5040–5048.
137. 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.
138. X. Zhang, X. Li, M. E. Reish, D. Zhang, N. Q. Su, Y. Gutiérrez, F. Moreno, W. Yang, H. O. Everitt and J. Liu, Nano Lett., 2018, 18, 1714–1723.
139. W. Ou, B. Zhou, J. Shen, T. W. Lo, D. Lei, S. Li, J. Zhong, Y. Y. Li and J. Lu, Angew. Chem., Int. Ed., 2020, 59, 6790–6793.
140. G. Baffou, I. Bordacchini, A. Baldi and R. Quidant, Light: Sci. Appl., 2020, 9, 108.
141. T. Perera, S. D. Gunapala, M. I. Stockman and M. Premaratne, J. Phys. Chem. C, 2020, 124, 27694–27708.
142. T. P. Rossi, M. Kuisma, M. J. Puska, R. M. Nieminen and P. Erhart, J. Chem. Theory Comput., 2017, 13, 4779–4790.
143. N. Asadi-Aghbolaghi, R. Rüger, Z. Jamshidi and L. Visscher, J. Phys. Chem. C, 2020, 124, 7946–7955.
144. P. Pandeya and C. M. Aikens, J. Phys. Chem. A, 2021, 125, 4847–4860.
145. F. Alkan and C. M. Aikens, J. Phys. Chem. C, 2018, 122, 23639–23650.
146. F. Alkan and C. M. Aikens, J. Phys. Chem. C, 2021, 125, 12198–12206.
147. F. Della Sala, J. Chem. Phys., 2022, 157, 104101.
148. S. D'Agostino, R. Rinaldi, G. Cuniberti and F. Della Sala, J. Phys. Chem. C, 2018, 122, 19756–19766.
149. J. W. You and N. C. Panoiu, IEEE J. Multiscale and Multiphys. Comput. Techn., 2019, 4, 111–118.
150. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864–B871.
151. R. G. Parr and W. Yang, Annu. Rev. Phys. Chem., 1995, 46, 701–728.
152. E. Runge and E. K. U. Gross, Phys. Rev. Lett., 1984, 52, 997–1000.
153. M. E. Casida and M. Huix-Rotllant, Annu. Rev. Phys. Chem., 2012, 63, 287–323.
154. M. L. Tiago and J. R. Chelikowsky, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 205334.
155. L. Román Castellanos, O. Hess and J. Lischner, Commun. Phys., 2019, 2, 47.
156. X. Wu, B. Liu, T. Frauenheim, S. Tretiak, C. Yam and Y. Zhang, J. Chem. Phys., 2022, 157, 214201.
157. J. G. Liu, H. Zhang, S. Link and P. Nordlander, ACS Photonics, 2018, 5, 2584–2595.
158. H. Breuer, F. Petruccione and S. Petruccione, The Theory of Open Quantum Systems, Oxford University Press, 2002.
159. H. Haug and A. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors, Springer Berlin Heidelberg, 1998.
160. T. E. Tesema, B. Kafle, M. G. Tadesse and T. G. Habteyes, J. Phys. Chem. C, 2017, 121, 7421–7428.
161. C. Boerigter, U. Aslam and S. Linic, ACS Nano, 2016, 10, 6108–6115.
162. R. Long and O. V. Prezhdo, J. Am. Chem. Soc., 2014, 136, 4343–4354.
163. R. Long, N. J. English and O. V. Prezhdo, J. Am. Chem. Soc., 2012, 134, 14238–14248.
164. T. Bora, D. Zoepfl and J. Dutta, Sci. Rep., 2016, 6, 26913.
165. Y. Sivan, I. W. Un and Y. Dubi, Faraday Discuss., 2019, 214, 215–233.
166. Z. Fusco, K. Catchpole and F. J. Beck, J. Mater. Chem. C, 2022, 10, 7511–7524.
167. K.-H. Chen, Y.-C. Pu, K.-D. Chang, Y.-F. Liang, C.-M. Liu, J.-W. Yeh, H.-C. Shih and Y.-J. Hsu, J. Phys. Chem. C, 2012, 116, 19039–19045.
168. A. A. Golubev, B. N. Khlebtsov, R. D. Rodriguez, Y. Chen and D. R. T. Zahn, J. Phys. Chem. C, 2018, 122, 5657–5663.
169. E. L. Keller and R. R. Frontiera, ACS Nano, 2018, 12, 5848–5855.
170. A. M. Brown, R. Sundararaman, P. Narang, A. M. Schwartzberg, W. A. Goddard and H. A. Atwater, Phys. Rev. Lett., 2017, 118, 087401.
171. 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.
172. J. C. Tully, J. Chem. Phys., 1990, 93, 1061–1071.
173. J. C. Tully, J. Chem. Phys., 2012, 137, 22A301.
174. B. F. Curchod, U. Rothlisberger and I. Tavernelli, ChemPhysChem, 2013, 14, 1314–1340.
175. W. Malone, B. Nebgen, A. White, Y. Zhang, H. Song, J. A. Bjorgaard, A. E. Sifain, B. Rodriguez-Hernandez, V. M. Freixas and S. Fernandez-Alberti, et al. , J. Chem. Theory Comput., 2020, 16, 5771–5783.
176. E. Fabiano, T. W. Keal and W. Thiel, Chem. Phys., 2008, 349, 334–347.
177. D. V. Makhov, C. Symonds, S. Fernandez-Alberti and D. V. Shalashilin, Chem. Phys., 2017, 493, 200–218.
178. D. V. Makhov, W. J. Glover, T. J. Martinez and D. V. Shalashilin, J. Chem. Phys., 2014, 141, 054110.
179. D. V. Makhov, K. Saita, T. J. Martinez and D. V. Shalashilin, Phys. Chem. Chem. Phys., 2015, 17, 3316–3325.
180. D. V. Makhov, T. J. Martinez and D. V. Shalashilin, Faraday Discuss., 2016, 194, 81–94.
181. 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.
182. M. Ben-Nun, J. Quenneville and T. J. Martnez, J. Phys. Chem. A, 2000, 104, 5161–5175.
183. L. V. Besteiro, X.-T. Kong, Z. Wang, G. Hartland and A. O. Govorov, ACS Photonics, 2017, 4, 2759–2781.
184. 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.
185. L. Yan, F. Wang and S. Meng, ACS Nano, 2016, 10, 5452–5458.
186. A. V. Akimov, A. J. Neukirch and O. V. Prezhdo, Chem. Rev., 2013, 113, 4496–4565.
187. F. Trani, G. Scalmani, G. Zheng, I. Carnimeo, M. J. Frisch and V. Barone, J. Chem. Theory Comput., 2011, 7, 3304–3313.
188. T. A. Niehaus, S. Suhai, F. Della Sala, P. Lugli, M. Elstner, G. Seifert and T. Frauenheim, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 085108.
189. 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.
190. X. Wu, S. Wen, H. Song, T. Frauenheim, S. Tretiak, C. Yam and Y. Zhang, J. Chem. Phys., 2022, 157, 084114.
191. B. Hourahine, B. Aradi, V. Blum, F. Bonafé, A. Buccheri, C. Camacho, C. Cevallos, M. Deshaye, T. Dumitrica and A. Dominguez, et al. , J. Chem. Phys., 2020, 152, 124101.
192. O. A. Douglas-Gallardo, M. Berdakin, T. Frauenheim and C. G. Sánchez, Nanoscale, 2019, 11, 8604–8615.
193. M. Berdakin, O. A. Douglas-Gallardo and C. G. Sanchez, J. Phys. Chem. C, 2019, 124, 1631–1639.
194. A. Mandal, M. A. Taylor, B. M. Weight, E. R. Koessler, X. Li and P. Huo, Chem. Rev., 2023, 123, 9786–9879.
195. J. J. Foley IV, J. F. McTague and A. E. DePrince III, Ab initio methods for polariton chemistry, arXiv, 2023, preprint, arXiv:2307.04881 [physics] DOI:10.48550/arXiv.2307.04881.
196. F. Herrera and J. Owrutsky, J. Chem. Phys., 2020, 152, 100902.
197. B. S. Simpkins, A. D. Dunkelberger and I. Vurgaftman, Chem. Rev., 2023, 123, 5020–5048.
198. A. D. Dunkelberger, B. S. Simpkins, I. Vurgaftman and J. C. Owrutsky, Annu. Rev. Phys. Chem., 2022, 73, 429–451.
199. J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux and T. W. Ebbesen, Angew. Chem., Int. Ed., 2012, 51, 1592–1596.
200. J. George, T. Chervy, A. Shalabney, E. Devaux, H. Hiura, C. Genet and T. W. Ebbesen, Phys. Rev. Lett., 2016, 117, 153601.
201. 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.
202. J. Lather, P. Bhatt, A. Thomas, T. W. Ebbesen and J. George, Angew. Chem., Int. Ed., 2019, 58, 10635–10638.
203. 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.
204. K. Hirai, H. Ishikawa, T. Chervy, J. A. Hutchison and H. Uji-i, Chem. Sci., 2021, 12, 11986–11994.
205. W. M. Takele, F. Wackenhut, Q. Liu, L. Piatkowski, J. Waluk and A. J. Meixner, J. Phys. Chem. C, 2021, 125, 14932–14939.
206. 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.
207. J. Lather, A. N. K. Thabassum, J. Singh and J. George, Chem. Sci., 2022, 13, 195–202.
208. S. Tretiak and S. Mukamel, Chem. Rev., 2002, 102, 3171–3212.
209. D. Timmer, M. Gittinger, T. Quenzel, S. Stephan, Y. Zhang, M. F. Schumacher, A. Lützen, M. Silies, S. Tretiak, J.-H. Zhong, A. De Sio and C. Lienau, arXiv, 2023, preprint, arXiv:2307.14708 DOI:10.48550/arXiv.2307.14708.
210. B. F. E. Curchod and T. J. Martínez, Chem. Rev., 2018, 118, 3305–3336.
211. R. Crespo-Otero and M. Barbatti, Chem. Rev., 2018, 118, 7026–7068.
212. J. Yang, Q. Ou, Z. Pei, H. Wang, B. Weng, Z. Shuai, K. Mullen and Y. Shao, J. Chem. Phys., 2021, 155, 064107.
213. Y. Zhang, T. Nelson and S. Tretiak, J. Chem. Phys., 2019, 151, 154109.
214. R. H. Tichauer, J. Feist and G. Groenhof, J. Chem. Phys., 2021, 154, 104112.
215. G. Groenhof and J. J. Toppari, J. Phys. Chem. Lett., 2018, 9, 4848–4851.
216. H. L. Luk, J. Feist, J. J. Toppari and G. Groenhof, J. Chem. Theory Comput., 2017, 13, 4324–4335.
217. G. Groenhof, C. Climent, J. Feist, D. Morozov and J. J. Toppari, J. Phys. Chem. Lett., 2019, 10, 5476–5483.
218. R. H. Tichauer, D. Morozov, I. Sokolovskii, J. J. Toppari and G. Groenhof, J. Phys. Chem. Lett., 2022, 13, 6259–6267.
219. M. Tavis and F. W. Cummings, Phys. Rev., 1968, 170, 379–384.
220. V. Rokaj, D. M. Welakuh, M. Ruggenthaler and A. Rubio, J. Phys. B: At., Mol. Opt. Phys., 2018, 51, 034005.
221. A. Mandal, T. D. Krauss and P. Huo, J. Phys. Chem. B, 2020, 124, 6321–6340.
222. J. Flick and P. Narang, J. Chem. Phys., 2020, 153, 094116.
223. D. S. Wang, T. Neuman, J. Flick and P. Narang, J. Chem. Phys., 2021, 154, 104109.
224. J. McTague and J. J. Foley, J. Chem. Phys., 2022, 156, 154103.
225. N. Vu, G. M. McLeod, K. Hanson and A. E. DePrince, J. Phys. Chem. A, 2022, 126, 9303–9312.
226. T. S. Haugland, E. Ronca, E. F. Kjønstad, A. Rubio and H. Koch, Phys. Rev. X, 2020, 10, 041043.
227. T. S. Haugland, C. Schäfer, E. Ronca, A. Rubio and H. Koch, J. Chem. Phys., 2021, 154, 094113.
228. A. E. DePrince, J. Chem. Phys., 2022, 154, 094112.
229. M. D. Liebenthal, N. Vu and A. E. DePrince, J. Chem. Phys., 2022, 156, 054105.
230. J. Flick, H. Appel, M. Ruggenthaler and A. Rubio, J. Chem. Theory Comput., 2017, 13, 1616–1625.
231. T. G. Philbin, Am. J. Phys., 2014, 82, 742–748.
232. R. R. Riso, T. S. Haugland, E. Ronca and H. Koch, Nat. Commun., 2022, 13, 1368.
233. A. Mandal and P. Huo, J. Phys. Chem. Lett., 2019, 10, 5519–5529.
234. R. R. Riso, T. S. Haugland, E. Ronca and H. Koch, J. Chem. Phys., 2022, 156, 234103.
235. C. Pellegrini, J. Flick, I. V. Tokatly, H. Appel and A. Rubio, Phys. Rev. Lett., 2015, 115, 093001.
236. M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel, I. V. Tokatly and A. Rubio, Phys. Rev. A, 2014, 90, 012508.
237. M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel and A. Rubio, Nat. Rev. Chem., 2018, 2, 1–16.
238. J. Flick, M. Ruggenthaler, H. Appel and A. Rubio, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 15285–15290.
239. T. E. Li, Z. Tao and S. Hammes-Schiffer, J. Chem. Theory Comput., 2022, 18, 2774–2784.
240. T. E. Li and S. Hammes-Schiffer, J. Chem. Phys., 2023, 158, 114118.
241. U. Mordovina, C. Bungey, H. Appel, P. J. Knowles, A. Rubio and F. R. Manby, Phys. Rev. Res., 2020, 2, 023262.
242. J. Fregoni, T. S. Haugland, S. Pipolo, T. Giovannini, H. Koch and S. Corni, Nano Lett., 2021, 21, 6664–6670.
243. T. G. Philbin, Am. J. Phys., 2014, 82, 742–748.
244. A. F. White, Y. Gao, A. J. Minnich and G. K.-L. Chan, J. Chem. Phys., 2020, 153, 224112.
245. T. E. Li, B. Cui, J. E. Subotnik and A. Nitzan, Annu. Rev. Phys. Chem., 2021, 73, 43–71.
246. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, 3rd edn, 2005.
247. C. Schäfer and G. Johansson, Phys. Rev. Lett., 2022, 128, 156402.
248. C. Schäfer, F. Buchholz, M. Penz, M. Ruggenthaler and A. Rubio, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2110464118.
249. A. F. Peterson, S. L. Ray and R. Mittra, IEEE Trans. Antennas Propag., 1998, 46, 276–294.
250. A. Gaita-Ariño, F. Luis, S. Hill and E. Coronado, Nat. Chem., 2019, 11, 301–309.
251. D. Sanvitto and S. Kéna-Cohen, Nat. Mater., 2016, 15, 1061–1073.
252. X. Zheng, F. Wang, C. Y. Yam, Y. Mo and G. Chen, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 195127.
253. Y. Zhang, S. Chen and G. Chen, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 085110.
254. Y. Zhang, C. Y. Yam and G. Chen, J. Chem. Phys., 2013, 138, 164121.
255. L. Chen, Y. Zhang, G. Chen and I. Franco, Nat. Commun., 2018, 9, 1–12.
256. T. Boolakee, C. Heide, A. Garzón-Ramírez, H. B. Weber, I. Franco and P. Hommelhoff, Nature, 2022, 605, 251–255.

---