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 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) that 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.


I. INTRODUCTION
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][2][3], enhance charge and energy transport [4][5][6], 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][9][10][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][8][9][10].
The occurrence of different effects due to the emergence of polaritons shows the complexity that arises when was recently extended to a matter-photon description within the linear-response framework of QEDFT [25]. 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,25,27,28]. A different approach within QEDFT is to solve the time-dependent Kohn-Sham equations coupled to the Maxwell equations in real time [3,[29][30][31]. 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 [32]. The Sternheimer approach has been used for a long time in the context of density-functional perturbation theory, for instance, for calculating phonon spectra [33]. Recent applications of the Sternheimer equation have also been used to compute the frequencydependent electronic response [34][35][36][37]. The Sternheimer equation has been formulated within the framework of time-dependent density-functional theory (TDDFT) that allows to study dynamic response of much larger complex systems as it includes only occupied orbitals [37,38]. 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 does not add to this scaling. This approach becomes useful for investigations in polaritonic chemistry or materials in nano-plasmonic 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 which 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 coupling of the bath modes to the molecule, we observe changes in the spectrum as the Lorentzian line shape turn into Fano resonances.
This article is structured as follows. First, we present the physical setting of a many-electron system coupled to photons in non-relativistic QED and subsequently present the linear-response setting of this framework. Secondly, we present in Sec. IV a derivation of the frequency-dependent Sternheimer approach for electronphoton 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 efficient 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 arises. Finally we present a conclusion and an outlook.

II. 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 non-relativistic 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 emphasizes 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 [39]: where E(r, t) and B(r, t) are the classical electric and magnetic fields, respectively. The transverse charge current j(r, t) represent 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 j(r, t) = ∂ ∂t P(r, t). The Maxwell's equations takes 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 [40]: where we introduced the displacement fieldD(r) = 0Ê (r) +P(r). Equation (1) will differ from other forms only in the choice of canonical variables and here we impose a commutation relation betweenB andD to be The expansion coefficientsd α andb α are respectively the field amplitudes of the electric displacement and the magnetic field, S α (r) are the mode functions andR = i er i is the total electronic dipole operator. In Eq. (4), we employed the dipole approximation when the electromagnetic field interacts with the matter system via the electronic dipole. Making a substitution of Eqs. (2)-(4) into (1) results to the following expression of the electromagnetic energy [40]: The displacement coordinateq α and momentum operator p α are related to the amplitudes asd α = √ 0 ω αqα and b α = 1/ 0pα where they satisfy the commutation relation [q α ,p α ] = i δ α,α . Equation (5) tells us that what would normally be the purely photonic Hamiltonian is now a mixture of matter and photon degrees. The term λ α in Eq. (5) is the light-matter coupling strength given as where the mode function is evaluated at the center of charge [41]. In deriving Eq. (5) we have assumed a finite photonic environment with appropriate boundary conditions. For real systems we, however, have usually a continuum of modes, i.e., 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, i.e. take 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 [42] and allows for modeling an open photonic environment [43]. Together with the Hamiltonian representing the bound charged particles, i.e. the kinetic energy, binding and interaction potentials, equation (5) constitutes the so-called Pauli-Fierz Hamiltonian in the length gauge [17,44]. In the case where we include time-dependent external perturbations, the length gauge Hamiltonian is given bŷ Here the N electrons are described by the electronic coordinatesr i and the momentum operatorp i , which satisfy the commutation relation r i ,p j = i δ ij . The interaction due to the longitudinal part of the photon field w(|r i −r j |) can be written as a mode-expansion in Coulomb gauge which for the free-space case results to the standard Coulomb interaction w(|r −r |) = n e 2 k n e ikn(r−r ) where k n = 2πn/L are the allowed wave vectors of the photon field for an arbitrarily large but finite box of length L [45]. 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,46]. In Sec. V B 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 time-dependent external potential and current in Eq. (7) can be split into v ext (r, t) = v(r) + δv(r, t), where v(r) describes the attractive potentials of the nuclei and δv(r, t) a classical external probe field that couples to the electronic subsystem. Also, the static part j α merely polarizes the vacuum and the time-dependent part δj α (t) then generates photons in the mode α. Either of these perturbations can be used to probe the coupled light-matter system.

III. LINEAR RESPONSE FORMULATION IN THE LENGTH GAUGE
To characterize the properties of a system, one can investigate its 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 [25]. 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 [25] δn(r, t) = dt d 3 r χ n n (r, t; r , t )δv(r , t ) Due to 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 to the response equation [25] δq α (t) = dt d 3 r χ qα n (t; r , t )δv(r , t ) The response functions χ n n , χ n qα , χ qα n and χ qα q α 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,25]. 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 [26] within the framework of QEDFT to treat electron-photon coupled systems [25]. This approach computes the excitation energies and oscillator strengths of either of the coupled response functions {χ n n , χ qα n } and χ n qα , χ qα q α by diagonalizing a pseudo-eigenvalue equation [25]. The Casida approach which requires both occupied and unoccupied Kohn-Sham orbitals in addition with the sampled photon modes scales favorably 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 [32]. Formulated within the framework of TDDFT, this method computes the linear density response due to an external weak perturbation [37,38,47] as well as non-linear responses [37,48]. This approach has several advantages the main one being that it relies only on the occupied Kohn-Sham orbitals thereby restricting the computation complexity for very large systems. In the following we extend this approach to treat electron-photon coupled system within the framework of QEDFT.

IV. THE STERNHEIMER APPROACH FOR ELECTRON-PHOTON COUPLED SYSTEMS
Practical ab-initio methods for computing optical excitation spectra are usually achieved by applying manybody methods that solve the correlated problem in an approximate way. One of 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 formalism for computing optical excitation spectra [48]. 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 [37,38]. Based on this advantage, an extension of this approach to the setting of QEDFT to treat complex atomic and molecular systems coupled to arbitrary large but finite number of photon modes is an important alternative method to existing QEDFT methods [25,29]. The derivation presented here is solely in the frequency space following that of Ref. [47]. 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 moderesolved 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 response due to a weak external po-tential δv(r, ω) are [25]: where the first-order Kohn-Sham potential and currents in Eqs. (10) and (11) are given in terms of the interacting density and photon coordinate responses including meanfield exchange-correlation kernels: The mean-field exchange-correlation kernels f n Mxc and f qα Mxc are defined to be the variation of the mean-field exchange-correlation potential with respect to the density and photon coordinate, respectively while the meanfield kernel g nα M is the variation of the mean-field current with respect to the electron density [25]. These kernels account for the correlations in the Kohn-Sham setting of QEDFT in linear response. The non-interacting response functions of the decoupled electronic and photonic subsystems of Eq. (10) and (11) are given explicitly as [25]: Here, k and ϕ k (r) are the ground-state energies and orbitals of the Kohn-Sham system and ω α is the frequency of the α mode. The parameters η and η shifts the poles (excitation energies) of Eqs. (14) and (15) 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 on the Kohn-Sham orbitals, we start by describing the unperturbed equilibrium setting of the coupled electron-photon system as this corresponds to the zerothorder 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 [49] where we have to solve the coupled Kohn-Sham equationŝ whereĥ =ĥ ([v, n, q α ]) is the ground-state Kohn-Sham Hamiltonian. The ground-state density can be obtained from the Kohn-Sham orbitals as n(r) = k |ϕ k (r)| 2 and the photon coordinate from Eq. (17). The meanfield exchange-correlation potential v Mxc (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 Eq. (10) and (11), we first start by substituting Eq. (14) into the density response n(r, ω) of Eq. (10). The density response can now be written in a form that includes a sum over only occupied orbitals as where the first-order response of the Kohn-Sham orbitals ϕ (±) k (r, ω) in Eq. (18) are given by Here, 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 ω −ĥ + k + iη and ω +ĥ − k + iη on Eqs. (19) and (20), respectively and using the static Kohn-Sham equation of Eq. (16) in concert with the completeness of the infinite set of ground-state Kohn-Sham orbitals, i.e., ∞ j=1 ϕ j (r )ϕ * j (r) = δ(r − r). After some algebra [30] this leads to the frequency-dependent Sternheimer equations of the following form where the first-order Kohn-Sham potential δv KS (r, ω) is given by δv KS (r, ω) = δv(r, ω) + d 3 r f n Mxc (r, r , ω)δn(r , ω) The response of the photon coordinate δq α (ω) in Eq. (10) to the external potential δv(r, ω) can be expressed in the following form where we substituted Eq. (15) into (10). The firstorder responses of the photon coordinates δq (+) α (ω) and δq (−) α (ω) are given explicitly as To obtain the response of the density and photon coordinate of Eqs. (18) and (24), we have to solve Eqs. (21)-(26) self-consistently. The self-consistency in solving these equations becomes evident by noting that the right-hand side of the Sternheimer equations (21) and (22) depends on the solution through δv KS (r, ω) which in turn depends on δn(r, ω) and δq α (ω). These two quantities depends on the first-order perturbed Kohn-Sham orbitals ϕ (±) k (r, ω) and photon responses δq (±) α (ω). It is important to note that the first-order response of the Kohn-Sham orbitals must satisfy the orthogonality condition with the groundstate Kohn-Sham orbitals [37,38]: 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 density and is related to the photo-absorption cross-section as σ(ω) = (4πω/3c) Trᾱ µµ [47]. Since the solutions ϕ ± k (r, ω) of Eqs. (26) and (25) are complex-valued, the density response of Eq. (18) become complex as well. This gives rise to the polarizability α µν (ω) having real and imaginary parts. The imaginary part of 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 [50].
In the decoupling limit of light and matter when |λ α | → 0, the Sternheimer equations (21) and (22) still retain the same form, however, the potential δv KS (r, ω) simplifies to that of an electron-only interacting system as f n Mxc → f n Hxc and f qα Mxc → 0. Also, the ground-state Kohn-Sham Hamiltonian in Eqs. (21) and (22) reduces toĥ =ĥ([v, n]) as v Mxc ([n, q α ]; r) → v Hxc ([n]; r), thus, decoupling the photon contribution of Eq. (17). 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 [30].
Details about the numerical treatment of Eqs. (21) and (22) have been discussed in the TDDFT framework of the frequency-dependent Sternheimer method [37,38]. 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 linewidth 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. (25) and (26). 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 N v and N c denote the number of occupied and unoccupied Kohn-Sham orbitals, respectively [25] and M describes the number of photon modes. The dimensionality of the matrix increases with N c 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 [25,30]. 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 [34,37,38]. This is evident since we can substitute Eqs. (24)-(26) into (23) such that the complexity rests in solving the Sternheimer equations (21) and (22). We have implemented the linear-response frequency-dependent Sternheimer equations (18) and (21)-(26) into the real-space code OCTO-PUS [37,51].

V. APPLICATIONS OF THE FREQUENCY-DEPENDENT STERNHEIMER APPROACH
In this section, we now apply the introduced electronphoton 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 [29,30]. This makes the frequency-dependent Sternheimer method of QEDFT a valid alternative for studying excited states properties of strongly coupled light-matter systems.
In the following, we first investigate the physical 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.

A. 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 Fig. (2)) and the photon displacement field (in Eq. (4)). 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 electronphoton system (for the photon displacement field).
The molecular system considered here is an azulene (C 10 H 8 ) molecule which is a bicyclic, nonbenzenoid aromatic hydrocarbon studied in Ref. [24]. We describe in detail how we compute the electronic structure of azulene in App. A. 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,25].
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  Fig. (1). The optical absorption spectra of the azulene molecule has been computed with TDDFT which captures the π − π * transition occurring at around 4.825 eV [52,53]. In Fig. (2 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. Based on the Kramers-Kronig relations, an absorption usually occurs simultaneously with an anomalous dispersion [50]. The anomalous dispersion describes a sudden change in the material's dispersion spectrum in the vicinity of a resonant absorption. We also find in real part of the polarizability an anomalous dispersion around the π − π * excitation which shows how its dispersive properties decreases 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 [54]. We now couple a single cavity mode in resonance to the π − π * excitation and tune the coupling λ = |λ| as in Fig. (2 b-d). For the Im {α xx (ω)}, an increasing coupling strength results to 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-Cumming model [28]. This splitting which is a characteristic of strong light-matter coupling shows how excitedstate properties of matter get modified when strongly coupled to a cavity mode. For the Re {α xx (ω)}, 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 Fig. (3  a) 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 {α xx (ω)} highlights that the dispersion properties of the matter system becomes modified due to strong light-matter coupling. The modifi- cation of dispersion properties for strongly coupled lightmatter 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 U = − 1 2ε0c Re {α xx (ω)} I, where I is the field intensity [55]. 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 [55]. 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 Fig. (2) where the Re {α xx (ω)} is modified under strong coupling and highlights a new perspective with potential applications in engineering optical dipole traps for neutral atoms or molecules.
Next, we study the spectral properties of the photon field when we probe the matter subsystem. This observ-able δ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 χ qα n is a commutator between photonic and electronic quantities [25]. 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 Fig. (4), 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 Fig. (2 a). 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 the hybrid light-matter features occur not only in the matter subsystem but also in the photon subsystem due to the selfconsistent interaction. For the Re {q α (ω)}, 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 {α xx (ω)} where the anomalous dispersion get modified due to strong light-matter coupling, the same holds true for the anomalous region in the spectrum of Re {q α (ω)} as shown in Fig. (3 b). Due to the self-consistent back-reaction between subsystems, we expect that the Re {q α (ω)} 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 App. A). The results presented here demonstrates that the electron-photon Sternheimer formalism is able to describe excited-state properties of strong light-matter coupled system.

B. 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 sample densely the electromagnetic vacuum. Since the Sternheimer approach for electron-only system is known to scale favorably with the system size [37,38], 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 λ α = 2 0LxLy Lz sin(ω α /c x 0 )e x , where x 0 = L x /2 is the position of the molecule in xdirection and ω α = αcπ/L x are the frequencies of the modes [25]. The volume V = L x L y L z with L x = 3250 µm, L y = 10.58Å and L z = 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 cut-off of 190.74 eV and a spacing between modes of 0.38 meV. Sampling the continuum of modes serve to constitute the line width of the excitations and also represent dissipation channels in the coupled system [25,43]. The one-dimensional sampling of mode frequencies that couple weakly to the matter subsystem will not capture the actual three-dimensional lifetimes. In the matteronly (uncoupled) case, we use a broadening η = 0.1 eV (as in Sec. V A) 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 Fig. (5) 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 X 1 Σ + → A 1 Σ + around 3.2 eV that corresponds to an electronic transition from the bonding to the antibonding σ-orbital [56,57]. 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. Now, we effectively enhance the coupling strength |λ| by reducing the cavity volume along the y-and zdirection. The results are shown in Fig. 6 where the blue line is the result shown in Fig. 5 that has a Lorentzian profile. We find that upon reducing the area L y L z , 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 [58,59]. 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 [60]. 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 [61]. Also, the broadening of the spectra (see Fig. (5 b)) and decrease in amplitude is a consequence of the interference [61]. 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.

VI. CONCLUSION AND OUTLOOK
In this work we presented a linear-response method that solves the response equations of non-relativistic 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 linearresponse 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 en-hanced, we find changes in the absorption spectrum as symmetric Lorentzian line shapes turn into asymmetric Fano line shapes.
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,62], 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 [37]. For the coupled electron-photon system, this will give access to higher order polarizabilities with signatures of strong light-matter coupling.

VII. ACKNOWLEDGMENTS
We acknowledge financial support from the European Research Council (ERC-2015-AdG-694097) and the SFB925 "Light induced dynamics and control of correlated quantum systems". This work was supported by the Excellence Cluster "CUI: Advanced Imaging of Matter" of the Deutsche Forschungsgemeinschaft (DFG), EXC 2056, project ID 390715994. The Flatiron Institute is a division of the Simons Foundation. To describe the azulene molecule studied in Sec. V, we compute the electronic structure of the azulene molecule using the same setup as that in Ref. [25] for the benzene molecule. That is, a cylindrical real space grid of 8Å in length with the radius of 6Å in the x-y plane and a spacing ∆x = ∆y = ∆z = 0.22Å is used. We treat the core electrons of the carbon and hydrogen atoms using the local-density approximation Troullier-Martins pseudopotentials [63]. Since the Sternheimer method considers only occupied orbitals, the 48 valence electrons that amounts to 24 doubly occupied orbitals are used in the computation. When we use the Casida approach, we include 500 unoccupied states which amounts to N v * N c = 24 * 500 = 12000 pairs of occupiedunoccupied states. To account for the electron-electron and electron-photon interaction in the linear-response QEDFT framework [25], we apply the adiabatic LDA (ALDA) to the Hartree exchange-correlation kernel (i.e. f n Hxc → f n Hxc,ALDA ) and the photon random-phase approximation to the photon exchange-correlation kernels (i.e. f n pxc , f qα pxc → f n p , f qα p ), respectively. We compute the photo-absorption cross-section of the molecule in free FIG. 7. Comparison of the photo-absorption cross-section of an azulene molecule using the electron-photon Sternheimer and Casida methods. (a) Absorption cross-section in free space (i.e. λα = 0). Coupling the cavity mode resonantly to the π − π * transition and increasing the coupling strength continuously as in (b) to (d) leads to a Rabi splitting into lower and upper polariton branches. Both approaches are in good agreement for the computed spectra.
space and inside the high-Q optical cavity and the results are shown in Fig. (7. This result is related to the imaginary part of Fig. (2) and it shows more excitations at higher energies.
In the linear-response setting of non-relativistic QED in the length gauge [25], the response quantities that arise have the same excitation energies but different oscillator strengths as shown in Fig. (8) for two observables of interest. This is because their respective response functions have the same excitation energies and different oscillator strengths [25].