Light–matter interaction of a molecule in a dissipative cavity from first principles

Recent experiments have explored regimes of light–matter interaction of molecular systems in nanophotonic cavities, where the interaction of a photon with the excitations of a single molecule can be substantially enhanced. When the loss in the cavity, characterized by a loss rate κ, dominates over the light–matter coupling rate g (g < κ), the system is in the weak-coupling regime. In this regime, the cavity mode enhances light–matter interactions, which can result in decay enhancement of molecular excited states via the Purcell effect.¹ This regime has been leveraged to improve spectroscopic detection of molecular species^2–6 without qualitatively modifying the electronic character of the molecular excitations.

A qualitative change in the physico-chemical properties of matter emerges when g > κ, which leads to the hybridization of excitations of matter with cavity photons to form polaritonic states.^7–14 This coherent non-perturbative regime is denoted as strong light–matter coupling (for more detailed discussions of the definition of light–matter coupling regimes, see Refs. 13 and 15–18), resulting in modifications of molecular chemical reactivity via electronic or vibrational strong coupling,^19–26 electrical²⁷ and excitonic²⁸ conductivity, optical properties,^29–35 energy transfer,^36–39 polariton lasing,^40,41 or superconductivity.^42–44 

In the intermediate regime, when g ∼ κ, the distinction between weak and strong light–matter coupling vanishes. This regime has been shown in optical cavities characterized by a small quality factor Q, where strong confinement of light in plasmonic⁴⁵ or conventional⁴⁶ cavities can lead to efficient coupling of light with excitons in many^47–49 or a few organic molecules.⁵⁰ Here, the polaritonic states are no longer well defined long-lived states⁵¹ nor can the light–matter interaction be treated as a weak perturbation to molecular excitations. The intermediate regime has been suggested to influence incoherent dissipative rates governing steady-state populations of molecular electronically excited states and can, for example, protect organic molecules from photo-oxidation processes.⁵² 

Theoretical descriptions of coupling of light with molecular excitations often rely on models based on cavity quantum electrodynamics (cQED) and theories of open-quantum systems.^53–56 Alternatively, parametrically described emitters can be coupled to an electromagnetic environment quantized within the linear-response theory of dielectric lossy media.^57–62 Fully ab initio predictions of the effects of light–matter coupling on the molecular properties require non-perturbative methods combining the power of cQED, theories of open-quantum systems, and many-body electronic structure theory.¹³ While longitudinal electromagnetic fields of dissipative matter excitations, such as plasmons, can be naturally included in time-dependent density functional theories by including their corresponding charge explicitly in the calculation,^63–71 standard theories must be extended to correctly account for transverse electromagnetic fields.

To bridge this critical gap in the field, we build on the formalism of quantum-electrodynamical density functional theory (QEDFT),^{12–14,72–76} briefly described in Sec. II, to incorporate cavity losses. We apply our generalized formalism to study ab initio vacuum light–matter interaction in both the weak- and strong-coupling regimes. Concretely, we study an isolated electronic transition in a benzene molecule in Sec. III A and a toluene molecule in Sec. III B, where many electronic transitions lead to an effective many-level electronic system. We show that light–matter interactions can induce interactions among the electronic states of the molecules that can result in energy transfer between excited states.

The Hamiltonian H for a non-relativistic system of M electronically excited states interacting with the quantized light field of N photon modes, in the absence of an external classical current and under the dipole approximation in the length gauge, is given by^9,14,72

where H_e is the electronic Hamiltonian, the kth quantized photon mode is given by the operators for the photon conjugate momentum $pk=iℏωk2(ak−ak†)$⁠, and the photon displacement coordinate is $qk=ℏ2ωk(ak+ak†)$⁠, where ℏ is the reduced Planck constant. Both photon operators are given in terms of the photon annihilation a_k and creation $ak†$ operators, and ω_k is the frequency of mode k. The electronic system couples to the photon modes through the total position operator $R=∑i=1Mri$ of the electronic system and q_k of the photonic system. This coupling is scaled by the cavity strength λ_k, which is given by

where E_k is the amplitude of the electric field at the center of the electronic charge density and e is the elementary charge. In a lossless cavity, λ_k represents a set of discrete modes whose excitation energies are well separated, and thus, usually only one or a few modes need to be considered in the description of light–matter coupling. Note that, in principle, molecular vibrational modes may be included in Eq. (1)^12,77 in future studies.

To calculate λ_k, as a representative example, we couple the molecular electronic degrees of freedom with a lossy transverse optical mode. Following the reasoning of Refs. 78 and 79, we replace the cavity strength of the single mode of a lossless 3D cavity λ_c at energy ℏω_c with many photon modes (the cavity mode coupled with the “modes of the universe”) whose spectral density follows a Lorentzian profile L(Δω, κ, ω_kc),

where

In principle, the appropriate realistic spectral density can be obtained using any method of quantizing lossy electromagnetic environments,^{56,57,62,78,80} although within QEDFT, longitudinal fields, such as plasmonic fields, should be included in the matter many-body Hamiltonian. The Lorentzian profile distributes the total cavity strength λ_c over a range of frequency values ω_k with a Lorentzian line shape around the central mode frequency (ω_kc = ω_k − ω_c), where the degree of broadening is controlled by the loss rate κ, which can be related, e.g., to the reflectivity of cavity mirrors. Here, the modes are uniformly spaced, that is, the frequency spacing Δω = ω_k+1 − ω_k for all k is constant. In general, transforming to or from another energy spacing Δω(ω_k) may offer computational benefits for an arbitrary electromagnetic environment; the detailed mathematical formalism is given in  Appendix B. Equation (3) also implies that the cavity strength obeys the sum rule,

where Δω(ω) → D(ω)dω in the continuum limit, D(ω) is the density of photon modes, and dω → 0.

In the following, we also define the coupling rate g_i,k that relates to the cavity strength λ_k and is a measure of the coupling strength between the cavity mode k and a particular electronic excitation i connecting the electronic ground state |g⟩ with an electronically excited state |e_i⟩ of the molecule,

Here, d_i = e⟨g|R|e_i⟩ are transition dipole moments that can be, for example, obtained together with their corresponding excitation energies ℏω_i from a standard linear-response time-dependent density-functional theory (TDDFT).⁸¹ These quantities can then be used as input parameters for cQED models describing light–matter interactions. This effective interaction rate naturally emerges in parametric cQED models, such as the Jaynes–Cummings model.⁸² We adapt these approaches by including the Lorentzian profile of the cavity strength to formulate a Fano-type cQED model that describes the linear-response of a discrete many-level electronic system coupled to a continuous photonic reservoir, permitting analysis in both the time and frequency domains. The relevant Hamiltonian H_cQED, in the rotating-wave approximation, is given by

where $σi†$ (σ_i) and $ak†$ (a_k) are creation (annihilation) operators for the ith of M excited electronic states and kth of N photon modes, respectively, H.c. is the Hermitian conjugate, $ℏωiel$ and ℏω_k are the mode energies, and g_i,k is the coupling defined as in Eq. (6).

The state of the system can be described in the single excitation subspace as

where the coefficients c₀, $ciel$⁠, and $ckph$ are time dependent. We plug this ansatz into the Schrödinger equation to obtain the following system of linear differential equations:

Note that the coefficient c₀ drops from the system of equations as the Hamiltonian conserves the number of excitations due to the rotating-wave approximation. This system of linear differential equations can be solved to obtain the eigenvalues, ℏω_l (i.e., energies), and eigenvectors, |v_l⟩ (i.e., polaritonic states), of the coupled system. With these parameters, we can calculate the electronic absorption spectrum and weight of the original, unmixed electronic and photonic states. In this cQED model, we describe the light–matter coupling via a many-mode Jaynes–Cummings term (adopting the rotating-wave approximation) and neglect the term ∝ R² appearing in the full light–matter coupling Hamiltonian given in Eq. (1). Several studies^9,83–86 have analyzed the effect of this term, which becomes, in particular, relevant in the ultrastrong coupling limit. In  Appendix A, we compare our full QEDFT calculations with R² contributions included to the cQED model without R² contributions and comment on their influence.

Within QEDFT, on the other hand, the electronic and photonic excitations, in the presence of nuclei, and their interactions are treated on the same quantized footing.^72,73,75,76 In particular, in this article, we adopt a linear-response time-dependent QEDFT that generalizes the Casida equation of TDDFT, first introduced in Refs. 14 and 76, where the electron–electron interactions included in TDDFT and the electron–photon interactions are solved simultaneously in an iterative solver, whereas the cQED approach involves adding electron–photon coupling to the TDDFT output before solving. We use the QEDFT implementation in a publicly available modified version of the pseudopotential, real-space DFT code OCTOPUS.^87–89 

In Sec. III A, we first demonstrate the present method on a model system by coupling an isolated electronic excitation of benzene to a cavity with varying cavity strength and losses. We note the qualitative similarities between the results calculated from first principles and those expected from the well-known, parametric Jaynes–Cummings model.⁸² In Sec. III B, we study a more complex system that consists of many electronically states of toluene coupled to a cavity with varying cavity strength and loss, respectively. Importantly, we observe interactions between electronically excited states mediated by the cavity that can result in population transfer between excited states.

To demonstrate our generalized QEDFT approach, we consider an example of a single benzene molecule placed in a lossy 1D optical cavity with a linearly polarized electric field along the x axis. The orientation of the molecule is shown in Figs. 1(a) and 1(b). In the first principles QEDFT calculation, we consider the full electronic structure with all single-particle excitations. Nevertheless, since the target excitation |e₁, 0⟩ with energy $ℏω1el=6.93$ eV and strong x-polarized transition dipole moment d_1,x = 0.96 eÅ tuned in resonance with the cavity mode is well separated in energy from other optically active excited states, benzene can be effectively understood as a two-level electronic system composed of the ground state |g, 0⟩ and an excited state |e₁, 0⟩ for the range of parameters we consider. We verify this simplification by calculating an electronic x-polarized absorption spectrum A_x of the molecule in the framework of linear-response QEDFT,^14,76

where $C=2me/(3ℏ2)$ (with m_e being the electron mass) is a frequency-independent pre-factor, ℏω_l is the mode energy, d_i,x is the x-component of the transition dipole moment of an electronic excitation i, and $Cilel$ (⁠$Cklph$⁠) is the projection of an original, unmixed electronic (photonic) state |e_i, 0⟩ (|g, 1_k⟩) to a resulting polaritonic state |v_l⟩. M excited electronic states of benzene are coupled with N photon modes to produce M + N hybrid electron–photon states called polaritons. For presentation purposes, we broaden the delta function δ(ω − ω_l) with a Lorentzian,

The spectral broadening ℏΓ can be interpreted as, for instance, the resolution of the spectrometer. We note that this artificial broadening is not included in our ab initio calculations and does not represent properties of the studied system. In  Appendix D, we show the influence of ℏΓ on the visibility of fine features.

The absorption spectrum with no light–matter coupling is shown in Fig. 1(c) and features a single dominant absorption peak corresponding to an electronic transition to |e₁, 0⟩ from |g, 0⟩, agreeing well with computational predictions in Ref. 76 and with experiments presented in Ref. 90. We verify that the excited eigenstate |e₁, 0⟩ is separated in energy at least 3 eV from the nearest eigenstates with a substantial transition dipole moment in the x-direction. For the range of coupling rates considered in this paper, |g, 0⟩ and |e₁, 0⟩ thus effectively form a two-level electronic system for the considered range of parameters.

Next, we place the molecule in a cavity whose central mode energy ℏω_c is in resonance with this transition (⁠$ωc=ω1el$⁠). In the single excitation subspace, this system thus corresponds closely to the Jaynes–Cummings model, which describes the interaction of a two-level system with a single resonant bosonic mode. In this scenario, schematically shown in Figs. 1(a) and 1(b), the two-level system can be represented by a ground state |g, 0⟩ and an excited electronic state |e₁, 0⟩ of a molecule whose transition dipole moment d₁ allows for efficient coupling with light and whose frequency is well separated from other electronic excitations. The bosonic mode is then represented by a single electromagnetic mode of an optical cavity (of quantized electric-field amplitude E_c and experiences losses with decay rate κ) that interacts with the electronic states of the molecule via the Jaynes–Cummings coupling rate $g=−1ℏd1⋅Ec$⁠. If the coupling rate is small compared to the loss rate of the cavity, when g < κ, the system is in the weak-coupling regime in which the cavity enhances the decay rate of the electronically excited state, i.e., the Purcell effect. This enhancement manifests as broadening of the electronic absorption spectrum of the molecule, as illustrated in Fig. 1(a). In contrast, when g > κ, the system is in the strong-coupling regime in which the electronically excited state and the cavity photon hybridize and form new hybrid light–matter polariton states denoted as |U⟩, “upper,” and |L⟩, “lower,” polaritons in Fig. 1(b). The polariton states are manifested in the molecular absorption spectrum as a doublet of peaks whose splitting is ∝g.

By tuning the cavity strength λ_c (which we assume to be polarized in the x-direction) and loss ℏκ, and therefore changing its spectral profile as shown in Fig. 2(a), we are able to control the regime of light–matter coupling in accordance with the predictions of the Jaynes–Cummings model. We first calculate a series of absorption spectra of the molecule coupled with the cavity field: In the calculations shown in Fig. 1(d), N = 5000 cavity modes whose cavity strength is distributed according to Eq. (3) are used to represent the single leaky cavity mode coupled with a single dominant electronically excited state among the M excited states of benzene, and the absorption spectrum is calculated with spectral broadening ℏΓ = 0.001 eV.

We first assume that the cavity is characterized by a constant cavity loss ℏκ = 0.001 eV and vary the value of cavity strength λ_c. The calculated absorption spectra are shown in Figs. 1(d–i) for λ_c = (0.001, 0.002, 0.008) e$V12$/nm or g/κ = (0.18, 0.36, 1.43), respectively. The largest ℏg = 0.23 meV $≪0.1ℏω1el$⁠, so the coupling strength remains far below the ultrastrong regime.^91,92 As we increase the coupling strength, we observe the transition from the weak-coupling regime in Fig. 1(d–i) where the molecular absorption spectrum features a single broadened Lorentzian peak to the strong-coupling regime in Fig. 1(d–iii) featuring the polaritonic peak doublet.

As the regime of light–matter coupling is a function of the ratio g/κ, we show in Figs. 1(d–iii) that the strong light–matter coupling regime can transition back into the weak-coupling regime when cavity losses are increased for a given cavity strength λ_c = 0.008 e$V12$/nm. Concretely, we increase the cavity loss rate from ℏκ = 0.001 eV in (iii) to ℏκ = 0.004 eV in (iv) to ℏκ = 0.008 eV in (v). We observe that with increasing loss rate, the doublet merges back into a single broadened Lorentzian peak as we reenter the weak-coupling regime. The peak in Fig. 1(d–v) is broader than the peak in Fig. 1(d–i) as the Purcell decay ∝g²/κ is larger in the former case.

To develop further physical intuition on how the N photonic states hybridize with the M electronic states to form M + N polaritonic states, we calculate the weight $Wilel=|Cilel|2$ (⁠$Wklph=|Cklph|2$⁠) of an original, unmixed electronic (photonic) state |e_i, 0⟩ (|g, 1_k⟩) in a resulting lth polaritonic state |v_l⟩. We plot the total weight from the electronic states $wlel=∑i=1MWilel$ and the total weight from the photonic states $wlph=∑k=1NWklph=1−wlel$ as the overlaid color of the curve in Fig. 1(d). Notably, we observe concentrations of the electronic character near the absorption peaks, even when separated in energy as in Fig. 1(d–iii), which proves the hybrid polaritonic character of the states composing the peaks.

Finally, we remark that in Fig. 5(a) in  Appendix A, we show that the results obtained from QEDFT in this parameter regime can be reproduced well using the cQED model, as expected in a model system.

First principles-based techniques enable quantitative predictions a priori of complex systems. Here, we systematically describe the effect of weak-to-strong light–matter coupling to an effective many-level electronic system, toluene. Although we again include all M electronically excited states from the ab initio calculation, here we discuss three x-polarized optically active excited states clustered in energy: |e₁, 0⟩ with energy $ℏω1el=6.64$ eV and x-polarized transition dipole moment d_1,x = 0.76 eÅ, |e₂, 0⟩ with $ℏω2el=6.71$ eV and d_2,x = 0.11 eÅ, and |e₃, 0⟩ with $ℏω3el=6.78$ eV and d_3,x = 0.08 eÅ. The x-polarized absorption spectrum of toluene in free space is plotted in Fig. 3(c) where we observe absorption peaks that scale as $∝di,x2$ corresponding to |e₁, 0⟩, |e₂, 0⟩, and |e₃, 0⟩. There are additional eigenstates that are further off-resonant or have smaller transition dipole moments. We neglect them in our discussion in the main text because they do not interact strongly with the cavity. Nonetheless, we still include their effects in our calculation. As we discuss in  Appendix D, for small enough spectral broadening ℏΓ, these excitations appear in the spectra as additional Fano resonances.

When the molecule is placed in a 1D lossy cavity with central mode frequency $ωc=ω1el$ and M = 36 000 photon modes, we observe the transition from the weak-to-strong coupling regime in the electronic absorption spectra in Fig. 3(d) with spectral broadening ℏΓ = 0.001 eV upon varying cavity strength λ_c and loss ℏκ, as shown in Fig. 2(b). At suitable cavity strength λ_c and loss ℏκ, we also observe Fano resonances resulting from photon-mediated interactions between electronic states.

We first assume that the cavity is characterized by cavity strength λ_c = 0.01e$V12$/nm and loss ℏκ = 0.01 eV or g/κ = 0.14. As we increase the coupling strength, we observe the transition from the weak-coupling regime in Fig. 3(d–i) where the molecular absorption spectrum features a broadened Lorentzian peak at $ℏω1el$ to the strong-coupling regime in Fig. 3(d–iii) where λ_c = 0.10 e$V12$/nm or g/κ = 1.4 features the polaritonic peak doublet. This doublet corresponds to the lower |L⟩ and upper |U⟩ polaritons in Fig. 3(a). The relative electronic character w^el of the doublet is lower than those of the peaks at $ℏω2el$ and $ℏω3el$⁠, which do not mix substantially with the photon modes.

In a many-level electronic system, electronic states may interact through the photon modes. We demonstrate the effects of this interaction by increasing the cavity strength λ_c to 0.43 e$V12$/nm and maintaining loss ℏκ = 0.01 eV (g/κ = 6.0) in Fig. 3(d–iii) from 0.10 e$V12$/nm in Fig. 3(d–ii). At this cavity strength, $ωp≈ω12el=ω2el−ω1el$⁠, so the upper polariton becomes nearly degenerate |e₂, 0⟩ with energy $ℏω2el$⁠, resulting in the upper polariton mixing with |e₂, 0⟩ to form an upper–lower polariton |UL⟩ and an upper–upper polariton |UU⟩. These conditions are schematically illustrated in Fig. 3(b).

Increasing loss ℏκ for the toluene-lossy cavity system such that |e₁, 0⟩ and |e₂, 0⟩ interact through the photon modes results in Fano resonances in the absorption spectrum. To capture this effect, from the strong-coupling regime in Fig. 3(d–iii) where the cavity strength λ_c = 0.43 e$V12$/nm and loss ℏκ = 0.01 eV, we increase loss ℏκ to 0.08 eV and maintain λ_c = 0.43 e$V12$/nm, or g/κ = 0.74, as shown in Fig. 3(d–iv). Further increasing loss ℏκ enables resonant electronic eigenstates to interact with other electronic eigenstates even further in energy: in Fig. 3(d–v), where λ_c = 0.43 e$V12$/nm and loss ℏκ is increased to 0.32 eV, the peak at $ℏω3el$ broadens. In addition, the individual polaritonic peaks merge into a single unresolved one, a characteristic of the weak-coupling regime.

In Fig. 5(b) in  Appendix A, we discuss how the results obtained from QEDFT can be interpreted using the cQED model and compare the absorption spectra obtained from both models. We show that up to small discrepancies arising from the rotating-wave approximation applied in cQED and the mean-field approximation applied in QEDFT, as well as the exclusion of the R² term arising from expansion of Eq. (1) in the cQED model, the two methods are in good quantitative agreement.

To elucidate the physical implications of the absorption spectra, we plot the weights $Wilel$ from QEDFT calculations in Figs. 4(a–i). The conditions correspond to the same cavity conditions as in Figs. 3(d–i), respectively, except that the input cavity mode frequency range is restricted from 6.35 eV to 6.95 eV in Figs. 4(a) and 4(b) for computational tractability as opposed to 4.85 eV–8.45 eV in Fig. 3(d). Absorption spectra are nearly identical for both energy ranges. Since the absorption spectra calculated from first principles in Fig. 3(d) agree with those of the cQED model in Fig. 5 parameterized with first principles calculations of the electronic structure and cavity mode profiles, we can use the cQED model to calculate the explicit time dependence of population transfer in Fig. 4(b) for a toluene-lossy cavity system initially prepared in |e₁, 0⟩.

In the weak-coupling regime in Fig. 4(a–i), we see that all electronic characters from |e₁, 0⟩, |e₂, 0⟩, and |e₃, 0⟩ are localized in one peak each. To demonstrate weak-coupling behavior in the time domain, the time decay of the toluene-lossy cavity system initialized in |e₁, 0⟩ is plotted in Fig. 4(b–i). As expected for an electronic system weakly coupled to photon modes, we observe exponential decay of the population of |e₁, 0⟩.

In the strong-coupling regime where the upper |U⟩ and lower |L⟩ polaritons are well-resolved in Fig. 3(d–ii), in Fig. 4(a–ii), the electronic character of |e₁, 0⟩ is concentrated within the two peaks corresponding to |U⟩ and |L⟩. We calculate the time decay of the toluene-lossy cavity system initialized in |e₁, 0⟩ in Fig. 4(b–ii) for the same cavity parameters as in Fig. 3(d–ii) and observe a signature of strong-coupling behavior: vacuum Rabi oscillations. In this regime, the difference in energy ℏω_p between the polaritons and |e₁, 0⟩ is much lower than the difference in energy $ℏω12el$ between |e₁, 0⟩ and |e₂, 0⟩, so there is negligible interaction between |e₁, 0⟩ and |e₂, 0⟩. This lack of interaction is apparent in both Figs. 4(a–ii) and 4(b–ii). In Fig. 4(a–ii), the two polaritonic peaks have no contribution from |e₂, 0⟩, and in Fig. 4(b–ii), we observe no population transfer to |e₂, 0⟩.

In Fig. 3(d–iii), λ_c is further increased such that electronic states |e₁, 0⟩ and |e₂, 0⟩ interact. In Fig. 4(a–iii), the peaks corresponding to |UU⟩ and |UL⟩ have contributions from both |e₁, 0⟩ and |e₂, 0⟩. The consequence in the time domain is population transfer between |e₁, 0⟩ and |e₂, 0⟩, as well as Rabi oscillations with the photon modes, suggesting that this toluene-lossy cavity system experiences both strong light–matter and strong matter–matter coupling via the cavity modes. As loss ℏκ increases to 0.10 eV in Fig. 3(d–iv), although population still transfers between |e₁, 0⟩ and |e₂, 0⟩, the excitation lifetime commensurately decreases compared to when ℏκ = 0.01 eV. Finally, as κ is increased further to 0.32 eV in Fig. 3(d–v), where |e₃, 0⟩ interacts non-negligibly with the cavity modes, we observe slight population transfer to |e₃, 0⟩ in Fig. 4(b–v).

Overall, we demonstrate that tuning the cavity strength λ_c and loss ℏκ can generate polaritonic states with contributions from several excited electronic states and can control population transfer to higher-lying states in energy and excitation lifetimes.

To summarize, we study ab initio correlated optical interactions in matter ranging from the weak-coupling to strong-coupling regime. As an example, we calculate excited-state dynamics and spectral responses of benzene and toluene as effective two-level and many-level electronic systems, respectively, under variable light–matter coupling controlled by the cavity strength λ_c and loss ℏκ. By tuning the cavity parameters, we notice transitions between the weak-coupling and the strong-coupling regimes where polaritonic states can be resolved. In the many-level electronic system, we observe Fano resonances in the electronic absorption spectrum resulting from interactions between electronically excited states mediated by the cavity. These interactions enable cavity-mediated population transfer between electronically excited states where the lifetimes and the degree of population transfer are controlled by the cavity parameters. We reproduce the first principles results using a cQED model parameterized with the QEDFT data and generally note excellent agreement.

This generalized QEDFT formalism is especially useful for predictions where interactions of single molecules with arbitrary electromagnetic environments dominate over vibrationally mediated losses. Looking forward, in electronic structure theory, we anticipate that the development of improved TDDFT methods and integration of light–matter interactions of molecular vibrations from first principles¹² will improve prediction accuracy for higher-lying excited states and enable studies of correlated cavity-electron–nuclei interactions,⁹³ respectively. In cQED, further understanding the effects of the correct inclusion of the R² term in the QEDFT formalism vs a conventional cQED model invoked in the present study will enable the extension of this method to the ultrastrong coupling regime. Finally, a natural extension of this work is to study the interactions of many molecules in lossy cavities from first principles, as formation of collective dark states leads to modifications of the polariton dynamics.^{19,25,30,39,41,52,94–96}

This work demonstrates a predictive technique for single-molecule experiments involving the weak- and strong-coupling light–matter regimes, including modifications of the excited-state potential energy surfaces in the field of polaritonic chemistry. This extension to QEDFT moves toward closing the loop between first principles calculations in electronic structure theory and parametric models of the quantum optics community.

This work was supported by the DOE “Photonics at Thermodynamic Limits” Energy Frontier Research Center under Grant No. DE-SC0019140. D.S.W. is an NSF Graduate Research Fellow. The Flatiron Institute is a division of the Simons Foundation. P.N. is a Moore Inventor Fellow.

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

Normalized x-polarized absorption spectra calculated with the cQED model for benzene and toluene with the same cavity conditions in Figs. 2(a) and 2(b) are plotted in Figs. 5(a) and 5(b), respectively. We note that the QEDFT and cQED agree quantitatively well for both benzene and toluene. There are more substantial differences between the spectra, specifically a blue shift of the cQED spectra, for the toluene study, which used higher coupling strengths than the benzene study. This trend agrees with a comparison between QEDFT and the Rabi model in Ref. 76. To understand the effect of the lack of the R² term in the cQED model that arises from expanding the last term in Eq. (1), in Fig. 5(b), we also plot the absorption spectra calculated with the present QEDFT formalism without the R² term. The effect of the R² is small for the coupling strengths here, so the differences in spectra can be more directly attributed to differences in the rotating-wave approximation made in the cQED model and the mean-field approximation made in QEDFT. Finally, Eqs. (8) and (9) can also be explicitly propagated in time given a set of initial conditions to plot population transfer, as shown in Fig. 4(b).

Here, we demonstrate how to transform the frequency spacing of a photon spectral profile from an exact calculation of the electromagnetic environment of a cavity to an arbitrary mode frequency spacing. Such an equivalency has clear numerical advantages whereby frequencies of interest can be more densely sampled at the cost of lowered frequency density at less relevant ranges for the same computational cost.

For specificity, we study the Hamiltonian H_cQED,cont in the cQED model, Eq. (7), expressed in terms of a continuously, constantly spaced frequency variable ω for photons and discrete frequencies $ωiel$ for electronically excited states,

We desire to change the integration variable ω to a new one Ω such that ω(Ω) has a more favorable spacing once the Hamiltonian is necessarily discretized for numerical calculations. Assuming a general functional dependence connecting the two frequency variables, note that the transformation of variables leads to modification of all integrals,

where we have defined the density of states D(Ω). This transformation is accompanied by re-definition of creation and annihilation operators of the continuum,

which follows from the requirement that the commutation relation $[ã(Ω),ã†(Ω′)]=δ(Ω−Ω′)$ holds.

We can rewrite H_cQED,cont as

where we have defined a new coupling constant $G̃i(Ω)$⁠. Equation (B4) can be discretized by assuming a constant spacing in variable Ω,

where ΔΩ = Ω_i+1 − Ω_i is the step of the discretization. The annihilation and creation operators of the continuum must be transformed back to the operators of the discrete set of modes as

which follows from the requirement that $[Ãi,Ãj†]=δij$⁠. We write the discrete version of H_cQED,cont, $H̃cQED$⁠, with new mode frequency spacing,

We use a publicly available modified version of the pseudopotential, real-space DFT code OCTOPUS.^87–89 The calculation of electronically excited states via the Casida TDDFT method first requires a converged ground state electronic structure on optimized molecular geometries. Molecular geometries are optimized with norm-conserving, plane-wave Hamann, Schlüter, Chiang, and Vanderbilt (HSCV) pseudopotentials⁹⁷ and a local density approximation (LDA) functional.⁹⁸ The real-space simulation box is parameterized with a mesh spacing of 0.16 Å and consists of spheres with 6 Å radius around each atom.

We calculate the converged ground state electron densities with a standard LDA functional for benzene and a long-range adapted LDA functional for toluene;⁹⁹ the long-range LDA (LRLDA) functional requires the ionization potential from a ΔSCF calculation. Ground state self-consistent field (SCF) energies converged within ∼1 meV/atom with 0.08 Å and 0.14 Å real-space mesh spacing for benzene and toluene, respectively, and a simulation box consisting of a 6 Å radius around each atom.

The energies of the excited states of interest converged within ∼1 meV/atom with 500 extra states in the Casida method. We further confirm the reliability of the excited electronic eigenstates that we couple to the cavity modes by noting that all excited states studied lie below the ionization potentials determined with a ΔSCF calculation.

Molecules can have excitations densely spaced in energy that all interact through the cavity photons. However, because excitations with low transition dipole moments or that are sufficiently off-resonant from the line width of the cavity modes do not interact strongly with the cavity, the effects of their interactions are small. Nonetheless, it is possible to observe their effects by increasing the measurement quality.

For instance, for toluene in the energy range shown in Fig. 3(c), there are several excitations beyond |e₁, 0⟩, |e₂, 0⟩, and |e₃, 0⟩ that do not have an appreciable effect on the absorption spectra when coupled to photons, as shown in Fig. 3(d). The only observable impact is a slight dip on the side of the lower polariton in Figs. 3(d–iii) caused by the eigenstate at 6.58 eV with transition dipole moment d_x = 0.01 eÅ. We study the impact of varying the spectral broadening ℏΓ on the absorption spectra for the conditions in Fig. 3(c), as demonstrated in Fig. 6. Physically, the spectral broadening can be lowered by decreasing the temperature. As the spectral broadening ℏΓ decreases from ℏΓ = (0.005, 0.002, 0.001, 0.0005, 0.0002) eV in Figs. 6(a)–6(e), respectively, we observe the appearance of an additional Fano resonance at the eigenenergy at 6.58 eV, a signature of the presence of an electronic eigenstate interacting with others through the photon modes. Given sufficient measurement resolution, we expect that further decreasing the spectral broadening ℏΓ would enable the observation of the remaining eigenstates in Fig. 3(c).