Frequency-Dependent Sternheimer Linear-Response Formalism for Strongly Coupled Light-Matter Systems.

The rapid progress in quantum-optical experiments, especially in the field of cavity quantum electrodynamics and nanoplasmonics, allows one to substantially modify and control chemical and physical properties of atoms, molecules, and solids by strongly coupling to the quantized field. Alongside such experimental advances has been the recent development of ab initio approaches such as quantum electrodynamical density-functional theory (QEDFT), which is capable of describing these strongly coupled systems from first principles. To investigate response properties of relatively large systems coupled to a wide range of photon modes, ab initio methods that scale well with system size become relevant. In light of this, we extend the linear-response Sternheimer approach within the framework of QEDFT to efficiently compute excited-state properties of strongly coupled light-matter systems. Using this method, we capture features of strong light-matter coupling both in the dispersion and absorption properties of a molecular system strongly coupled to the modes of a cavity. We exemplify the efficiency of the Sternheimer approach by coupling the matter system to the continuum of an electromagnetic field. We observe changes in the spectral features of the coupled system as Lorentzian line shapes turn into Fano resonances when the molecule interacts strongly with the continuum of modes. This work provides an alternative approach for computing efficiently excited-state properties of large molecular systems interacting with the quantized electromagnetic field.

Strong interactions between light and matter within enhanced photonic environments such as optical cavities and plasmonic devices have attracted great attention in recent years. The signature of such strong interactions is the formation of new hybrid light–matter states (polaritons), which are manifested by a Rabi splitting in the spectrum of the coupled system. These new states of matter can be used to control, for instance, chemical reactions,[1−3] enhance charge and energy transport,[4−6] and selectively manipulate electronic excited states,[7] to name a few examples. Such coupled light–matter systems have the tendency to exhibit significantly different properties than the uncoupled subsystems even at ambient conditions, which suggests various interesting applications in chemistry and material science.[8−11] These intriguing effects caused by the emergence of polaritons manifest strongly in the excited-state properties of the coupled systems, for example, in the absorption or emission spectra.[7−10]

The occurrence of different effects due to the emergence of polaritons shows the complexity that arises when light and matter strongly mix. Because of this inherent complexity of the coupled light-matter system, the theoretical description of these effects are nontrivial. Quite often, the coupled system is studied with quantum optical models that potentially oversimplify the matter subsystem. One such simplification selects just a few energy levels of an atomic or molecular system and couples it to the photon modes of an optical cavity.[12−14] Another common simplification is on the photon side, where a realistic cavity that is normally open is reduced to just a few modes that cannot account for the finite lifetime of excitations. However, in many cases these phenomenological models are not sufficient to capture important details of the coupled system, for instance, the emergence of bound states in the continuum,[15] how collective strong coupling leads to local modification of chemical properties,[16] and in cavity-modified chemistry where the reaction-rate is reduced under cavity-induced resonant vibrational strong coupling.[2,17] This calls for ab initio methods, which allow one to treat from first-principles complex matter systems interacting with the quantized field[18−20] within nonrelativistic quantum electrodynamics (QED).[21−23] Nonrelativistic QED is the basis of all approaches to theoretically capture the emerging physics of polaritonic chemistry.[22] Yet so far it remains debated which aspects of this highly complex theory are the ones mainly responsible for the observed changes in chemistry.[24] Thus, first-principles approaches provide a mostly unbiased approach to this fundamental question of cavity-modified chemistry. Among the existing first-principles methods for treating strongly coupled light-matter systems, quantum electrodynamical density-functional theory (QEDFT) has become a valuable approach for describing ground- and excited-state properties of complex matter systems coupled to a photonic environment.[25,26] The Casida-like approach[27] common to molecular and quantum chemistry was recently extended to a matter-photon description within the linear-response framework of QEDFT.[26] The feasibility of treating the excited-state properties of a single molecule and an ensemble of molecules coupled to a realistic description of a cavity has been demonstrated.[16,26,28,29] A different approach within QEDFT is to solve the time-dependent Kohn–Sham equations coupled to the Maxwell equations in real time.[3,30−32] On the one hand, one major advantage the real-time approach has is that it scales favorably with the system size, as it involves only occupied Kohn–Sham orbitals, but to obtain a converged response spectrum requires a long time-propagation, which is not favorable for larger systems. On the other hand, the Casida approach requires both occupied and unoccupied orbitals, and it also scales with the number of photon modes considered.

In addition to these methods, there is another successful scheme that can combine the strengths of the previously mentioned methods known as the Sternheimer approach.[33] The Sternheimer approach has been used for a long time in the context of density-functional perturbation theory, for instance, for calculating phonon spectra.[34] Recent applications of the Sternheimer equation have also been used to compute the frequency-dependent electronic response.[35−38] The Sternheimer equation has been formulated within the framework of time-dependent density-functional theory (TDDFT), which allows one to study the dynamic response of much larger complex systems, as it includes only occupied orbitals.[38,39] One of a few advantages the Sternheimer approach has over real-time TDDFT is that it is formulated in the frequency space, and the responses at different frequencies can be computed independently of each other allowing for the use of parallelization schemes that speed up the computation. Another advantage is that, since the responses at different frequencies can be treated independently, we can compute any part of the spectrum without necessarily starting from the zero frequency. In this work, we extend the frequency-dependent Sternheimer approach of TDDFT to the framework of QEDFT. An advantage the electron-photon Sternheimer approach has over the Casida approach is that it scales well with the system size, since only occupied orbitals are treated explicitly, and the arbitrarily many but finite photon modes that can be included do not add to this scaling. This approach becomes useful for investigations in polaritonic chemistry or materials in nanoplasmonic cavities that usually consider a large number of particles interacting with the electromagnetic field. We start by showing the applicability of the method in capturing not only the absorption properties of strongly coupled light-matter system but also the modified dispersion properties of the coupled system for the case of an azulene molecule. In addition, we show the spectra of the photon field that capture similar features of strong light-matter coupled systems indicating how the hybrid characteristics can be viewed from either of the subsystems at the same time highlighting the cross-talk between the subsystems. To show the efficiency of the Sternheimer approach, we study the absorption spectrum of a lithium hydride (LiH) molecule coupled to a continuum of photon modes. For the case of coupling the molecule weakly to half a million photon modes, we recover the spectrum of the free space case. By effectively enhancing the light-matter coupling strength of the bath modes to the molecule, we observe changes in the spectrum as the Lorentzian line shape turns into Fano resonances.

This article is structured as follows. First, we present the physical setting of a many-electron system coupled to photons in nonrelativistic QED and subsequently present the linear-response setting of this framework. Second, we present in Section Three a derivation of the frequency-dependent Sternheimer approach for electron–photon coupled systems in the linear-response regime formulated within the framework of QEDFT and discuss numerical details of the Sternheimer scheme. In the next section, we investigate the complex polarizability of a molecular system coupled to a high-Q optical cavity and highlight how the absorption and dispersion properties get modified due to strong light–matter coupling. Also, we show for the same molecular system the polaritonic features that arise in the spectra of the photon field. Furthermore, we demonstrate the efficiency and low computational cost of the Sternheimer approach by coupling a LiH molecule to a (discretized) continuum of states of photon modes and show the physical effects that arise. Finally we present a conclusion and an outlook.

From Microscopic Fields to the Quantum Description of Light–Matter Interaction

We are interested in the dynamics of matter interacting with the electromagnetic field within the setting of nonrelativistic QED where both constituents of the coupled system are treated on an equal quantized footing. While this setting of slowly moving charged particles can be deduced from QED, concepts from classical electrodynamics are equally instructive to arrive at this description. In this regard, we lay emphasis on the full description of the electromagnetic field that couples to a matter system. Our starting point is the inhomogeneous microscopic Maxwell equations for the transverse part of the electromagnetic field[40]where E(r, t) and B(r, t) are the classical electric and magnetic fields, respectively. The transverse charge current j(r, t) represents both the free and bound current. If we consider j(r, t) to represent only the bound current, then it is related to the polarization P(r, t) of the matter as . The Maxwell’s equations take into account the back-reaction of the matter on the electromagnetic field. For a quantum mechanical description, the field variables are promoted to field operators in the Heisenberg picture. In this representation, the energy of the transverse electromagnetic field can be expressed as[41]where we introduced the displacement field . Equation will differ from other forms only in the choice of canonical variables, and here we impose a commutation relation between and to be where ϵ is the Levi-Civita symbol. For any photonic environment of varying geometry, the fields in eq can be expanded in the modes[41]The expansion coefficients and are, respectively, the field amplitudes of the electric displacement and the magnetic field, Sα(r) are the mode functions, and is the total electronic dipole operator. In eq , we employed the dipole approximation when the electromagnetic field interacts with the matter system via the electronic dipole. We will later (see end of Section 3) briefly discuss how to go beyond this common simplification. Making a substitution of eqs –(6) into (3) results in the following expression of the electromagnetic energy.[41]The displacement coordinate and momentum operator are related to the amplitudes as and where they satisfy the commutation relation . Equation tells us that what would normally be the purely photonic Hamiltonian is now a mixture of matter and photon degrees. The term λα in eq is the light–matter coupling strength given aswhere the mode function is evaluated at the center of charge.[42] In deriving eq we assumed a finite photonic environment with appropriate boundary conditions. For instance, we can assume a planar cavity in the z-direction, while in the x- and y-directions we have the usual free-space or periodic boundary conditions (see also Figure ). In the z-direction we would then have where n is a unit vector normal to the cavity surfaces. For real systems we, however, have usually a continuum of modes; that is, the cavity geometry is open to free infinite space. We can approximate this situation by extending the quantization volume of the electromagnetic field beyond the photonic environment and thus work with a discretized continuum. By making this discretization finer and finer, that is, by taking the quantization volume to infinity, we can approximate the open-cavity situation arbitrarily well. The discrete continuum description of the photon field has the advantage that it accounts for the emission or absorption of a photon in real space[43] and allows for modeling an open photonic environment.[44] Together with the Hamiltonian representing the bound charged particles, that is, the kinetic energy, binding, and interaction potentials, eq constitutes the so-called Pauli-Fierz Hamiltonian in the length gauge.[17,45] In the case where we include time-dependent external perturbations, the length gauge Hamiltonian is given byHere the N electrons are described by the electronic coordinates and the momentum operator , which satisfy the commutation relation . The interaction due to the longitudinal part of the photon field can be written as a mode-expansion in Coulomb gauge, which for the free-space case results in the standard Coulomb interactionwhere k = 2πn/L are the allowed wave vectors of the photon field for an arbitrarily large but finite box of length L.[46] For the transverse field, we consider an arbitrarily large but finite number of photon modes M. It is important to note that, when we sample a large number of modes to describe the photon continuum, we might need to use the bare mass of the electrons instead of the renormalized physical mass.[21,47] In Section 4.2 of this work, we make the common assumption that only the sampled continuum due to a cavity or photonic nanostructure is changed with respect to the free space case. The rest of the continuum of modes not affected by the cavity is subsumed in the already renormalized physical mass of the electrons. The coupled light-matter system can be perturbed externally using the time-dependent external potential and current in eq , which can be split intowhere v(r) describes the attractive potentials of the nuclei, and δv(r, t) indicates a classical external probe field that couples to the electronic subsystem. For the external perturbing charge current, the static part jα merely polarizes the vacuum, and the time-dependent part δjα(t) then generates photons in the mode α. The physical implication of an external current that acts on the photon field can be best understood from eqs and (2). An external current will generate photons, which constitute a magnetic and electric field. In contrast to the classical external scalar potential v(r, t), these induced fields are fully quantized. Either of these perturbations can be used to probe the coupled light–matter system.

Figure 1

Schematic setup of an azulene molecule confined within a high-Q optical cavity. The cavity field is polarized along the x-axis with mode coupling λα, and the photon propagation vector is along the cavity axis of length L in the z-direction. The frequency of the photon mode is ωα.

Linear Response Formulation in the Length Gauge

To characterize the properties of a system, one can investigate the system’s response to an external perturbation. In the case of a weak external perturbation, we have access to linear response properties of the system such as its polarizability, which gives access to its excitation energies and oscillator strengths. Usually this requires knowledge of the linear density response, and in the case of a coupled matter-photon system, we have access to the displacement field.[26] In the length gauge, the linear response of the electron density n(r, t) to the external potential δv(r, t) and charge current δjα(t) yields the response equation[26]Because of the coupling between light and matter, we can equally compute the linear response of the photon displacement coordinate qα(t) due to the external potential δv(r, t) and current δjα(t) that results in the response equation[26]

The response functions , , , and are intrinsic properties of the electron–photon coupled system, which can be computed to obtain excited-state properties of the system. However, computing these response equations or the response functions directly is usually very challenging even for the electron-only system. One possible way to do this efficiently is to reformulate the response equations using the Maxwell-Kohn–Sham system of QEDFT that reproduces the same response of the density and photon coordinate.[18,19,26] In such a setting, the response functions can be computed approximately giving access to, for instance, excitation energies and oscillator strengths. We recently extended the Casida equation[27] within the framework of QEDFT to treat electron-photon coupled systems.[26] This approach computes the excitation energies and oscillator strengths of either of the coupled response functions and by diagonalizing a pseudoeigenvalue equation.[26] The Casida approach, which requires both occupied and unoccupied Kohn–Sham orbitals in addition to the sampled photon modes, is efficient for small coupled systems.[16] However, for larger electronic systems coupled to many photon modes, this leads to a rapid increase in computational effort in the Casida approach, as the Casida matrix equation increases in size.

An alternative approach, which rather computes the response equations instead of the response function, is the frequency-dependent Sternheimer equation.[33] Formulated within the framework of TDDFT, this method computes the linear density response due to an external weak perturbation[38,39,48] as well as nonlinear responses.[38,49] This approach has several advantages, the main one being that it relies only on the occupied Kohn–Sham orbitals, thereby relieving the computation complexity for very large systems. In the following we extend this approach to treat an electron-photon coupled system within the framework of QEDFT.

The Sternheimer Approach for Electron–Photon Coupled Systems

Practical ab initio methods for computing optical excitation spectra are usually achieved by applying many-body methods that solve the correlated problem exactly or in an approximate way. A few of the most popular ab initio methods to determine the electronic structure of atoms or molecules in quantum chemistry are Hartree–Fock theory, configuration interaction (CI), coupled cluster (CC), or (TD)DFT.[50] In terms of accuracy, CI and CC[51] are both favorable over (TD)DFT. Because of the improved accuracy of CC, this has led to its extension to quantum electrodynamics coupled cluster theory (QED-CC)[20,52] to treat strongly coupled light-matter systems. QED-CC is, however, limited to small matter systems and only a few photon modes. To overcome this limitation in the matter system size and photon modes, we need to employ other electronic structure methods that scale favorably with system size. One such many-body methods is TDDFT, which is considered a very promising methodology, since it provides a good balance between accuracy and computational cost. Within the context of TDDFT there exist different formalisms for computing optical excitation spectra.[49] The Sternheimer formalism is a standard method in electronic structure theory for computing the spectra of many-body systems.[34−39] The frequency-dependent Sternheimer method formulated within TDDFT is a perturbative approach on the Kohn–Sham orbitals that computes the density response without relying explicitly on unoccupied Kohn–Sham orbitals.[38,39] On the basis of this advantage, an extension of this approach to the setting of QEDFT to treat complex atomic and molecular systems coupled to an arbitrary large but finite number of photon modes is an important alternative method to existing QEDFT methods.[26,30] The derivation presented here is solely in the frequency space following that of ref (48). In an electron-only description, the Sternheimer approach obtains only electronic observables such as the electron density response. However, when this method is formulated within the QEDFT framework we have, in addition to the density response, the response of the photon displacement coordinate (field). The mode-resolved response of the field gives access to physical processes such as the absorption or emission process. Starting with the reformulation of the density and photon displacement coordinate responses in the QEDFT framework, the coupled responses due to a weak external potential δv(r, ω) are[26]where the first-order Kohn–Sham potential and currents in eqs and (13) are given in terms of the interacting density, photon coordinate responses, and kernels.

The mean-field plus exchange-correlation kernels and are defined to be the variation of the mean-field plus exchange-correlation potential (i.e., vMxc = vM + vxc) with respect to the density and photon coordinate, respectively, while the mean-field kernel is the variation of the current with respect to the electron density.[26] These kernels account for the correlations in the Kohn–Sham setting of QEDFT in linear response. Given the exact Mxc kernels we recover the exact response of the coupled light–matter system. In practice we will need to approximate the xc part of the kernels. We note that, for , the exact xc conctribution is zero.[26] The noninteracting response functions of the decoupled electronic and photonic subsystems of eqs and (13) are given explicitly as[26]Here, ϵ and φ(r) are the ground-state energies and orbitals of the Kohn–Sham system, and ωα is the frequency of the α mode. The parameters η and η′ shift the poles (excitation energies) of eqs and (17) to the lower half of the complex plane and are, in general, not equal in both uncoupled systems.

Since the Sternheimer method is a perturbative approach to the Kohn–Sham orbitals, we start by describing the unperturbed equilibrium setting of the coupled electron-photon system, as this corresponds to the zeroth-order of a perturbation expansion, for example, of the density. For this case, we start by describing the static Kohn–Sham system of ground-state QEDFT,[53] where we must solve the coupled Kohn–Sham equationswhere is the ground-state Kohn–Sham Hamiltonian. The ground-state density can be obtained from the Kohn–Sham orbitals as n(r) = ∑|φ(r)|2 and the photon coordinate from eq . The mean-field plus exchange-correlation potential vMxc(r) represents the longitudinal interactions between the electrons as well as all the transversal interactions of the electrons with the photon field.

To solve for the linear density response and photon displacement coordinate of eqs and (13), we first start by substituting eq into the density response n(r, ω) of eq . The density response can now be written in a form that includes a sum over only occupied orbitals aswhere the first-order response of the Kohn–Sham orbitals in eq are given byHere, solving for the Kohn–Sham orbital responses is highly involved, since we need to first determine infinitely many Kohn–Sham orbitals and evaluate an infinite sum over all these orbitals. However, this can be circumvented by acting with and on eqs and (22), which results in the following equations.Using the static Kohn–Sham eq in the above two equations simplifies the right-hand sides to

In the next step, we take advantage of the completeness of the infinite set of ground-state Kohn–Sham orbitals, that is, in eqs and (24), which simplifies to the frequency-dependent Sternheimer equations of the following formwhere the first-order Kohn–Sham potential δvKS(r, ω) is given by

The response of the photon coordinate δqα(ω) in eq to the external potential δv(r, ω) can be expressed in the following formwhere we substituted eq into (12). The first-order responses of the photon coordinates and are given explicitly as

To obtain the response of the density and photon coordinate of eqs and (28), we must solve eqs –(30) self-consistently. The self-consistency in solving these equations becomes evident by noting that the right-hand side of the Sternheimer eqs and (26) depends on the solution through δvKS(r, ω), which in turn depends on δn(r, ω) and δqα(ω). These two quantities depend on the first-order perturbed Kohn–Sham orbitals and photon responses . It is important to note that the first-order response of the Kohn–Sham orbitals must satisfy the orthogonality condition with the ground-state Kohn–Sham orbitals.[38,39]

From solving the self-consistent Sternheimer equations we can compute the dynamic polarizability of the coupled system, which is given in terms of the variation of the densityand is related to the photoabsorption cross-section as .[48] Since the solutions of eqs and (29) are complex-valued, the density response of eq becomes complex as well. This gives rise to the polarizability α(ω) having real and imaginary parts. The imaginary part of the polarizability describes the absorption of radiation, and the real part defines the refraction properties of the matter system due to a perturbation from an external electromagnetic field.[54]

In the decoupling limit of light and matter when |λα| → 0, the Sternheimer eqs and (26) still retain the same form; however, the potential δvKS(r, ω) simplifies to that of an electron-only interacting system as and . Also, the ground-state Kohn–Sham Hamiltonian in eqs and (26) reduces to as vMxc([n, qα]; r) → vHxc([n]; r), thus decoupling the photon contribution of eq . The derivation of the Sternheimer scheme for the electron density and photon displacement coordinate responses in the QEDFT framework due to a weak external charge current δjα(ω) follows the same steps as above.[31]

Details about the numerical treatment of eqs and (26) have been discussed in the TDDFT framework of the frequency-dependent Sternheimer method.[38,39] Therefore, we only summarize features in the numerical application of these equations. First, the positive infinitesimal parameter η is required for numerical stability for the solution of the Sternheimer equations close to the resonance frequencies, as it removes the divergences. It is also necessary to obtain the imaginary part of the polarizability. In addition, this parameter accounts for the artificial line width that represents the finite lifetimes of the excitations. Our extension of the Sternheimer method to treat electron-photon coupled systems introduced the small positive infinitesimal parameter η′ that enters the self-consistent Sternheimer equations as in eqs and (30). This parameter is necessary to ensure that the poles at ωα are finite. In our simulations we found that ℏη′ = 0.001 eV is the ideal value to obtain converged results, and we used ℏη = 0.1 eV.

For the electron-photon Casida approach, the resulting dimension of the coupled but truncated matrix is where Nv and Nc denote the number of occupied and unoccupied Kohn–Sham orbitals, respectively,[26] and M describes the number of photon modes. The dimensionality of the matrix increases with Nc and M-photon modes. We have been so far able to treat a finite matter system coupled to 150 000 modes with an efficient massive parallel implementation of the Casida equation.[26,31] In terms of scaling with system size, the electron-photon Sternheimer approach is better when compared to the Casida approach, since it still scales the same as the electron-only Sternheimer case.[35,38,39] This is evident since we can substitute eqs –(30) into (27) such that the complexity rests in solving the Sternheimer eqs and (26). We implemented the linear-response frequency-dependent Sternheimer eqs and (25)–(30) into the real-space code OCTOPUS.[38,55] Let us finally comment on the restriction to dipole light-matter coupling. The full Pauli-Fierz Hamiltonian of nonrelativistic QED uses the full minimal-coupling prescription and hence includes all multipole interactions.[17,21,45] And also for the full theory, QEDFT[32,53] has been formulated and applied.[32] This shows that a linear-response formulation of QEDFT with full minimal-coupling is possible. A detailed derivation and implementation of linear-response QEDFT for minimal coupling is, however, currently still missing.

Applications of the Frequency-Dependent Sternheimer Approach

In this section, we now apply the introduced electron–photon frequency-dependent Sternheimer approach for studying excited-state properties of molecular systems coupled to a photon mode or a continuum of modes. This approach has been validated by comparing the optical absorption spectrum of a single benzene ring coupled to photons to that obtained using the electron-photon Casida and time-propagation methods of QEDFT.[30,31] This makes the frequency-dependent Sternheimer method of QEDFT a valid alternative for studying excited-state properties of strongly coupled light-matter systems.

In the following, we first investigate a cavity QED setup in which a single molecule is strongly coupled to a photon mode of a high-Q cavity where we expect to capture the hallmark of strong light–matter coupling (Rabi splitting). In the next setup, we include a large but finite number of photon modes that simulates the electromagnetic vacuum and investigate situations where a molecular system couples weakly and strongly to the continuum.

Single-Molecule Strong Coupling

The first example studies intrinsic properties of a strongly coupled light-matter system that is commonly not considered, for instance, the real part of the polarizability (in Figure ) and the photon displacement field (in eq ). These quantities are particularly interesting, as they give insight into the dispersive properties of the coupled system (for the real part of the polarizability) and how energy is exchanged between the electron-photon system (for the photon displacement field).

Figure 2

Spectrum of an azulene molecule in free space (i.e., λ = 0) and coupled to a high-Q optical cavity (i.e., λ > 0) showing the line shapes characteristic of the real and imaginary parts of the polarizability near the π–π* resonance at 4.825 eV. (a) Region near the resonance where the Re is asymmetric about the resonance while the Im is symmetric about the resonance. Coupling the cavity mode resonantly to the π–π* transition and increasing the coupling strength continuously as in (b–d) results in a Rabi splitting into lower and upper polariton branches, each of which has an asymmetric line shape for the different Re .

The molecular system considered here is an azulene (C10H8) molecule, which is a bicyclic, nonbenzenoid aromatic hydrocarbon studied in ref (25). We describe in detail how we compute the electronic structure of azulene in the Supporting Information. Before looking at how these observables get modified due to strong light-matter coupling, we will first present the absorption spectra (obtained from the imaginary part of the polarizability) of the molecular system strongly coupled to photons that captures the Rabi splitting between polaritonic peaks.[8,26]

To study the spectral properties of the coupled system we now confine the azulene molecule inside an optical high-Q cavity that couples to a photon mode with increased strength. The cavity field is polarized along the x-direction with a coupling strength λα as shown in Figure . The optical absorption spectra of the azulene molecule has been computed with TDDFT, which captures the π–π* transition occurring at 4.825 eV.[56,57] In Figure a, we show the x-component of the polarizability of the uncoupled azulene molecule. The imaginary part of the polarizability captures a sharp peak occurring at 4.825 eV due to the π–π* excitation. On the basis of the Kramers–Kronig relations, an absorption usually occurs simultaneously with an anomalous dispersion.[54] The anomalous dispersion describes a sudden change in the material’s dispersion spectrum in the vicinity of a resonant absorption. We also find in the real part of the polarizability an anomalous dispersion around the π–π* excitation, which shows how its dispersive properties decrease when the excitation energy increases. This is characterized by the asymmetric line shape about this resonance, while the imaginary part is symmetric as usually observed.[58] We now place the molecule at the center of the high-Q optical cavity and make the common assumption to describe the cavity by one effective mode. The coupling defined in eq in this particular case of a planar cavity is where L is the length of the cavity, and A is the surface corresponding to the mode volume. The values for λ that are normally used are for cavity volumes on the order of 103 μm3.[59] With this mode volume, the strong coupling regime is achieved by collectively coupling an ensemble of emitters to the photon mode.[1] In the single-molecule limit, recent experiments in picocavity setups have demonstrated effective volumes less than 1 nm3 for achieving strong light-matter coupling.[60,61] On the theory side, investigations into the effective volumes for enhancement of optical fields have been explored[62] with suggestions for nanoplasmonic structures with volumes as small as 0.15 nm3.[63] To explore the strong light-matter coupling regime in this setup, we choose values of λ = 0.01, 0.03, 0.05 au, which correspond to effective volumes LA = 17.5, 2.1, 0.74 nm3, respectively. For the Im , an increasing coupling strength results in an increased Rabi splitting of the π–π* peak into lower and upper polaritonic branches, where the lower branch has more intensity, compared to the upper polaritonic peak as measured in experiments[8] and not captured by common phenomenological models such as the Jaynes-Cummings model.[29] This splitting, which is a characteristic of strong light-matter coupling, shows how excited-state properties of matter get modified when strongly coupled to a cavity mode. For the Re , we find for each of the lower and upper polariton peaks for different λ, asymmetric line shapes about their respective excitation energies indicating anomalous dispersion usually occurs simultaneously with absorption even for strongly coupled systems. In addition, the anomalous dispersion can be controlled for strongly coupled systems by varying the coupling strength. This is clearly shown in Figure where the anomalous dispersion (in particular, for the lower polariton) is smaller for the coupled case when compared to the uncoupled result. The emergence of polaritonic features in the Re highlights that the dispersion properties of the matter system become modified due to strong light-matter coupling. The modification of dispersion properties for strongly coupled light-matter systems has potential in controlling optical dipole traps. This can be made clear by considering the interaction potential of the induced dipole moment normally expressed as , where I is the field intensity.[64] The standard approach for realizing optical dipole traps is by laser detuning from a specific resonance of the bare matter system, for instance, laser detuning from an atomic resonance such that the dipole potential minima occur at regions with maximum intensity for red-detuned traps.[64] For polaritonic resonances that emerge in strongly coupled light-matter systems, the optical dipole traps that can be realized by detuning the external field from these polaritonic resonances can be controlled by strongly coupling to the photon field. This is evident in Figure where the Re is modified under strong coupling and highlights a new perspective with potential applications in engineering optical dipole traps for neutral atoms or molecules.

Figure 3

(a) The real part of the polarizability of azulene showing the change in the anomalous dispersion in free space (λα = 0) and when coupled to a cavity mode (λα > 0). (b) The analogous anomalous dispersion in the photon spectrum occurs only when both subsystems are coupled. (a, b) This feature can be controlled by coupling to a cavity mode.

Next, we study the spectral properties of the photon field when we probe the matter subsystem. This observable δqα(ω) is now accessible, since we treat the photon field as a dynamical part of the coupled light-matter system. We note that the displacement field in this case represents a mixed (matter and photon) spectroscopic observable, since its response function is a commutator between photonic and electronic quantities.[26] The observable δqα(ω) indicates how the photon field reacts in a standard absorption or emission measurement when the system is probed by an external field represented by the potential δv(r, ω). In Figure , we show the spectrum of the photon displacement coordinate in free space (when λ = 0) and coupled to a cavity mode (when λ > 0). As expected the free space case has no response, since light and matter decouple, and we have access only to the observables in Figure . However, coupling to the photon mode and increasing the coupling strength λ > 0 we observe in the imaginary part of δqα(ω) a Rabi splitting into lower and upper polaritons peaks. The polaritonic peaks are asymmetric about the π–π* excitation energy to which the mode was initially coupled to, and the lower polariton peaks are negative with more intensity compared to the upper polarition. Physically, this result highlights that excitations due to an external perturbation from δv(r, ω) can be exchanged between the coupled subsystems and that the hybrid light-matter features occur not only in the matter subsystem but also in the photon subsystem due to the self-consistent interaction. For the Re , we also find for each of the lower and upper polariton branches an asymmetric line shape about the energies of the respective polariton peaks with varying strengths for different λ. In analogy to the Re where the anomalous dispersion gets modified due to strong light-matter coupling, the same holds true for the anomalous region in the spectrum of Re as shown in Figure . Because of the self-consistent back-reaction between subsystems, we expect that the Re can be made to influence the optical dipole potential thereby controlling how the matter subsystem is trapped in the field. It is important to note that, for the responses of the subsystems, the excitation energies of the strongly coupled system are the same but with differing oscillator strengths (see the Supporting Information). The results presented here demonstrate that the electron-photon Sternheimer formalism is able to describe excited-state properties of strong light-matter coupled system.

Figure 4

Spectrum of the photon displacement coordinate of an azulene molecule in free space (i.e., λ = 0) and coupled to a high-Q optical cavity (i.e., λ > 0). (a) No response, as the photons are decoupled. Coupling the cavity mode resonantly to the π–π* transition and increasing the coupling strength lead to a splitting into lower and upper polaritonic branches in the photonic spectrum as shown in (b–d) for Im . The Re for these cases show an antisymmetric line shape opposite to Re in Figure .

Changes in the Matter Spectral Features

In this section, we consider the case where a molecular system is coupled explicitly to a wide range of photon modes and show how spectral features of the system change when we effectively increase its coupling to the continuum of the electromagnetic field. This computation will at the same time show the advantages the Sternheimer approach has over the Casida approach in terms of scaling with the number of photon modes.

We now consider as matter system a lithium hydride (LiH) molecule coupled to a wide range of photon modes that densely sample the electromagnetic vacuum. Since the Sternheimer approach for an electron-only system is known to scale favorably with the system size,[38,39] the focus here will be to demonstrate that the photon modes do not add to this scaling. Here we sample modes of a quasi one-dimensional mode space by employing the coupling , where x0 = L/2 is the position of the molecule in the x-direction, and ωα = αcπ/L are the frequencies of the modes.[26] The volume V = LLL with L = 3250 μm, L = 10.58 Å, and L = 2.65 Å are chosen such that the sampled modes couple weakly to the molecular system, and we assume a constant mode function in the y- and z-directions.

In this first example, we couple the molecule to 500 000 photon modes of a one-dimensional mode space with an energy cutoff of 190.74 eV and a spacing between modes of 0.38 meV. Sampling the continuum of modes serves to constitute the line width of the excitations and also represents dissipation channels in the coupled system.[26,44] The one-dimensional sampling of mode frequencies that couple weakly to the matter subsystem will not capture the actual three-dimensional lifetimes. In the matter-only (uncoupled) case, we use a broadening ℏη = 0.1 eV (as in Section 4.1) to account for the finite lifetime of the excited states. When the molecule is coupled weakly to the photon continuum, we obtain the same spectral broadening as the uncoupled case. The results of this calculation is shown in Figure , where we compare the photoabsorption cross-section of the uncoupled LiH molecule and the case when it is weakly coupled to 500 000 photon modes. We find that the two results are qualitatively the same, which is evident for the lowest electronic transition X1Σ+ → A1Σ+ around 3.2 eV that corresponds to an electronic transition from the bonding to the antibonding σ-orbital.[65,66] This result shows that the weak coupling of the molecule to the continuum of modes reproduces the results of the matter-only case. We note that obtaining this result using the electron-photon Casida approach will increase the computational cost drastically even for the case of coupling to 100 000 photon modes. Computationally, this result demonstrates that the electron-photon Sternheimer method scales favorably not only with system size but also with the number of photon modes.

Figure 5

(a) Photoabsorption cross-section of a LiH molecule coupled to 500 000 photon modes (blue dashed) in a quasi one-dimensional cavity and its comparison to the uncoupled case (red solid). (b) Enlarged view of the X1Σ+ → A1Σ+ transition around 3.2 eV where we observe a slight deviation in the peak amplitude between the uncoupled and the case coupled to the continuum.

Now, we effectively enhance the coupling strength |λα| by reducing the cavity volume along the y- and z-directions. For this purpose, we choose four different areas LL = 28, 0.35, 0.23, and 0.12 Å2, and the length L is fixed as given above with the same number of modes. We chose very small areas to be able to obtain the desired transition between spectral line shapes for the single-molecule case studied here. This will not be the case in collective coupling, since the coupling strength scales as the square-root of the number of identical particles. The results are shown in Figure , where the blue line is the result shown in Figure that has a Lorentzian profile. We find that, when we reduce the area LL, this effectively enhances the coupling to the photon continuum such that the symmetric Lorentz line shapes turn into asymmetric Fano line shapes. Fano resonances occur due to the interference of discrete quantum states with a continuum of states.[67,68] The asymmetry is characterized as a ratio of the transition amplitude to a given discrete state and that of a transition to a continuum state.[69] As this ratio becomes finite due to strong coupling to the continuum, this indicates the onset of a competition between constructive and destructive interference that gives rise to the asymmetric line shape.[70] Also, the broadening of the spectra (see Figure ) and decrease in amplitude are consequences of the interference.[70] These results show the changes in the spectral features of excited states of a matter system strongly coupled to the electromagnetic continuum. Thus, the electron–photon Sternheimer approach is a valid alternative method for studying excited-state properties of real systems strongly interacting with the quantized electromagnetic field.

Figure 6

(a) Photoabsorption spectrum of a LiH molecule coupled to a continuum of photon modes where the coupling is effectively enhanced by changing the cavity volume via the area LL. The Lorentzian line shapes turn into Fano line shapes for increasing effective coupling strength. (b) Magnified view of the X1Σ+ → A1Σ+ transition around 3.2 eV where we observe clearly the asymmetry of the Fano resonances when compared to Lorentzian line shape (blue solid line).

Conclusion and Outlook

In this work we presented a linear-response method that solves the response equations of nonrelativistic QED in the length gauge setting. The approach is based on the Sternheimer equation formulated within the framework of QEDFT that is capable of computing excited-state properties of strongly coupled light-matter systems. This approach serves as an alternative linear-response method for studying response properties of large systems coupled to the quantized electromagnetic field, since it scales favorably with the system size, as it utilizes only the occupied Kohn–Sham orbitals, and it also scales favorably with the number of photon modes. Using the Sternheimer approach we computed different observables of strongly coupled systems. These observables showed how both the dispersion and absorption properties of the matter system changes with potential applications in modifying and controlling optical dipole traps. Also, we showed examples where we lift the restriction to one cavity mode in the dipole approximation and sampled densely the electromagnetic continuum. In one case we showed that, when a LiH molecule is weakly coupled to the photon continuum, we reproduce the free space absorption spectrum of the molecule. When the coupling strength between the light and matter is effectively enhanced, we find changes in the absorption spectrum as symmetric Lorentzian line shapes turn into asymmetric Fano line shapes.

Our investigations in this work employed the adiabatic local-density approximation (ALDA) to treat the Hartree exchange-correlation kernel that accounts for the correlation between electrons. The reason for this choice was to show that, even with the simplest functional (ALDA), the extended electron-photon Sternheimer approach still captures the hallmark of strong light-matter coupling (Rabi splitting) and other features as shown in Figure . It is, however, important to investigate how ALDA performs in comparison to hybrid functionals such as B3LYP or PBE0 in describing the peak position of excitation energies, oscillator strengths, and lifetimes of the polaritonic resonances. This is particularly important, as it will, on the one hand, provide information on how electron correlation affects properties of the Rabi splitting and, on the other hand, scrutinize the reliability of ALDA in describing correlations in strongly coupled electron-photon systems. The electron-photon Sternheimer method presented here is a suitable approach for studying excited-state properties of large systems coupled to a single mode or to the electromagnetic continuum. In the fast-growing field of polaritonic chemistry, where there is an ongoing debate about the mesoscopic scale of quantum-collectively of coupled molecules,[16,24] ab initio methods such as the electron-photon Sternheimer method become desirable to capture intricate details of the complex interactions between the coupled subsystems. Another important property of the Sternheimer approach is that it can be generalized to higher orders to obtain higher-order polarizabilities by solving a hierarchy of Sternheimer equations.[38] For the coupled electron-photon system, this will give access to higher-order polarizabilities with signatures of strong light–matter coupling.