Light-Matter Response in Non-Relativistic Quantum Electrodynamics: Quantum Modifications of Maxwell's Equations
Johannes Flick, Davis M. Welakuh, Michael Ruggenthaler, Heiko Appel, Angel Rubio
LLight-Matter Response in Non-Relativistic Quantum Electrodynamics: QuantumModifications of Maxwell’s Equations
Johannes Flick,
1, 2, ∗ Davis M. Welakuh, † Michael Ruggenthaler, ‡ Heiko Appel, § and Angel Rubio
2, 3, ¶ John A. Paulson School of Engineering and Applied Sciences,Harvard University, Cambridge, MA 02138, USA Max Planck Institute for the Structure and Dynamics of Matter andCenter for Free-Electron Laser Science & Department of Physics,Luruper Chaussee 149, 22761 Hamburg, Germany Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA (Dated: June 19, 2019)We derive the full linear-response theory for non-relativistic quantum electrodynamics in the longwavelength limit, show quantum modifications of the well-known Maxwell’s equation in matter andprovide a practical framework to solve the resulting equations by using quantum-electrodynamicaldensity-functional theory. We highlight how the coupling between quantized light and matterchanges the usual response functions and introduces new types of cross-correlated light-matter re-sponse functions. These cross-correlation responses lead to measurable changes in Maxwell’s equa-tions due to the quantum-matter-mediated photon-photon interactions. Key features of treating thecombined matter-photon response are that natural lifetimes of excitations become directly accessi-ble from first principles, changes in the electronic structure due to strong light-matter coupling aretreated fully non-perturbatively, and for the first time self-consistent solutions of the back-reaction ofmatter onto the photon vacuum and vice versa are accounted for. By introducing a straightforwardextension of the random-phase approximation for the coupled matter-photon problem, we calculatethe first ab-initio spectra for a real molecular system that is coupled to the quantized electromag-netic field. Our approach can be solved numerically very efficiently. The presented framework leadsto a shift in paradigm by highlighting how electronically excited states arise as a modification of thephoton field and that experimentally observed effects are always due to a complex interplay betweenlight and matter. At the same time the findings provide a new route to analyze as well as proposeexperiments at the interface between quantum chemistry, nanoplasmonics and quantum optics.
I. INTRODUCTION
Recent years have seen tremendous experimental ad-vances in the nascent field of strongly-coupled light-matter systems [1, 2]. In particular, new experimentaladvances have been demonstrated in polaritonic chem-istry [3–5], solid-state physics [6], biological systems [7],nanoplasmonics [8, 9], two-dimensional materials [10, 11]or optical waveguides [12], among others.In this so-called strong-coupling regime, as a resultof mixing matter and photon degrees-of-freedom [13,14], novel effects emerge such as changes in chemicalpathways [15–17] ground-state electroluminescence [18],cavity-controlled chemistry for molecular ensembles [19,20], or optomechanical coupling in optical cavities [21],new topological phases of matter [22], superradiance [23]or superconductivity [24].Due to the inherent complexity of such coupledfermion-boson problems described in general by quantumelectrodynamics (QED), the theoretical treatment is usu-ally drastically simplified. One common approximation ∗ Electronic address: fl[email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected] ¶ Electronic address: [email protected] is to restrict the description of the system to simplifiedeffective models that heavily rely on input parameters.Current state of the art in the theoretical description ofstrong light-matter coupling very often employs a few-level approximation. This approximation leading to theRabi or Jaynes-Cummings model [25, 26] in the single-emitter case, or the Dicke model [27] in the many-emittercase, is however often not sufficient [28, 29], in particularwhen observables besides the energy are of interest [29],such as in experimental setups involving the modificationof chemical reactivity [1].Alternatively, in linear spectroscopy, the current theo-retical description is built on the semi-classical approxi-mation [30]. Herein, the many-particle electronic systemis treated quantum mechanically and the electromagneticfield appears as an external perturbation. As an externalperturbation, the electromagnetic field probes the quan-tum system, but is not a dynamical variable of the com-plete system (see also supplemental material S1). Sincein the strong-coupling regime light and matter must beon the same level, a semi-classical approximation is notadequate and the feedback between light and matter hasto be considered.It is, however, long known that the radiative life-times are finite. Furthermore, experimentally excited-state properties are usually inferred from (de)excitationsof the photon field, which is in stark contrast to the usualsemi-classical theoretical description based solely on the a r X i v : . [ qu a n t - ph ] J un electronic subsystem.In free-space, this mismatch can be circumvented sinceexcited-state properties such as radiative lifetimes ofatoms and molecules can be calculated perturbativelyusing the theory of Wigner-Weisskopf [31] employingthe Markov approximation. However, this perturbativetreatment of the coupling of light and matter becomesinsufficient in the case that strong light-matter couplingis achieved, e.g., due to many emitters or due to reducingthe mode volume of a cavity. In such cases the Markovapproximation breaks down and the Wigner-Weisskopftheory is not applicable anymore [32]. Additionally it isnot straightforward how to extend the original formula-tion of Wigner-Weisskopf to many electronic levels andhence to an ab-initio treatment of electronic systems.As a consequence, the current literature shows a largegap for situations, where light and matter is stronglycoupled and observables such as excited-state densities,radiative lifetimes, or electron-photon correlated observ-ables of interest. A good example is the control of theradiative lifetimes of single molecules [33, 34] by chang-ing the environment. In such cases the properties of themany-body system are changed, e.g., the excitation en-ergies and lifetimes are strongly modified. This happensbecause certain modes of the photon vacuum field areenhanced which can lead to a strong coupling of lightwith matter. Alternatively, increasing the number of par-ticles leads to an enhancement of the coupling due tothe self-consistent back-reaction of matter onto the pho-ton field and vice versa. It is important to realize thatsuch changes are non-perturbative for the photon field aswell as for the matter subsystem and hence need a self-consistent implementation. This fact is most pronouncedin the appearance of polaritonic states and their influenceon chemical and physical properties of matter [1, 13].In this paper, we close this gap by presenting a prac-tical and general framework that subsumes electronic-structure theory, nanoplasmonics, and quantum optics.We present a new description that challenges our con-ception of light and matter as distinct entities [35] andthat expresses the excited states as modifications of thephoton field. We do so by introducing a linear-responseformalism for coupled matter-photon systems. This for-malism leads naturally to modifications of Maxwell’sequations and the ability to calculate radiative life-times in arbitrary photon environments, including free-space, high-Q optical cavity or nanoplasmonic struc-tures. We make this approach practical by introducing alinear-response framework for quantum-electrodynamicaldensity-functional theory (QEDFT) [13, 14, 36–38]. Thisdevelopment is specifically timely since QEDFT hasnow been successfully applied to real systems in equi-librium [39] – which demonstrates the feasibility of ab-initio strong-coupling calculations – yet an accurate andefficient approach to excited states within QEDFT hasbeen missing. This work therefore furthermore closes agap within the QEDFT framework. II. LIGHT-MATTER INTERACTION IN THELONG WAVELENGTH LIMIT
Our fundamental description of how the charged con-stituents of atoms, molecules and solid-state systems, i.e.,electrons and positively charged nuclei, interact is basedon QED [13, 40–42], thus the interaction is mediatedvia the exchange of photons. Adopting the Coulombgauge for the photon field allows us to single out thelongitudinal interaction among the particles which givesrise to the well-known Coulomb interaction and leavesthe photon field purely transversal. Assuming then thatthe kinetic energies of the nuclei and electrons are rela-tively small, allows us to take the non-relativistic limitfor the matter subsystem of the coupled photon-matterHamiltonian, which gives rise to the so-called Pauli-FierzHamiltonian [13, 37, 42] of non-relativistic QED. In anext step one then usually assumes that the combinedmatter-photon system is in its ground state such that thetransversal charge currents are small and that the cou-pling to the (transversal) photon field is very weak. Be-sides the Coulomb interaction it is then only the physicalmass of the charged constituents (bare plus electromag-netic mass [42]) that is a reminder of the photon fieldin the usual many-body Schr¨odinger Hamiltonian. Inthis work, however, we will not disregard the transver-sal photon field, which makes the presented frameworkmuch more versatile and applicable to situation outsideof standard quantum mechanics (see also appendix B).
A. Novel Spectroscopy from quantum descriptionof light-matter interaction
In the following, we consider cases, in which the semi-classical approximation breaks down, as outlined in theintroduction. From the Pauli-Fierz Hamiltonian, wemake the long-wavelength or dipole approximation inthe length-gauge [43] since the wavelength of the pho-ton modes are usually much larger than the extendof the electronic subsystem which leads (in SI units)to [36, 37, 44] [45]ˆ H ( t ) = ˆ H e + M (cid:88) α =1 (cid:34) ˆ p α + ω α (cid:18) ˆ q α − λ α ω α · R (cid:19) (cid:35) + j α ( t ) ω α ˆ q α , (1)where ˆ H e is the standard many-body electronic Hamil-tonian [46]. We further restrict ourselves to arbitrarilymany but a finite number M of modes α ≡ ( k , s ) with s being the two transversal polarization directions that areperpendicular to the direction of propagation k . The fre-quency ω α and polarization (cid:15) α that enter in λ α = (cid:15) α λ α with λ α = S k ( r ) / √ (cid:15) and mode function S k ( r ) definethese electromagnetic modes. S k ( r ) is normalized, hasthe unit 1 / √ V with the volume V and we choose a ref-erence point r where we have placed the matter subsys-tem to determine the fundamental coupling strength [47].These photon modes couple via the displacement coor-dinate ˆ q α = (cid:113) (cid:126) ω α (ˆ a α + ˆ a † α ), where ˆ q α is given in termsof photon annihilation ˆ a α and creation ˆ a † α operators, tothe total dipole moment R = (cid:80) Ni =1 e r i [48]. The ˆ q α appears in the contribution of mode α to the displace-ment field ˆ D α = (cid:15) ω α λ α ˆ q α [43]. Further, the conju-gate momentum of the displacement coordinate is givenby ˆ p α = − i (cid:113) (cid:126) ω α (ˆ a α − ˆ a † α ). Besides a time-dependentexternal potential v ( r , t ), we also have an external per-turbation j α ( t ) that acts directly on the mode α of thephoton subsystem. Here j α ( t ) is connected to a classicalexternal charge current J ( r , t ) that acts as a source forthe inhomogeneous Maxwell’s equation.Formally, however, due to the length-gauge transfor-mations, the j α ( t ) corresponds to the time-derivativeof this (mode-resolved) classical external charge cur-rent [36, 37] (see also appendix A). Physically the staticpart j α, merely polarizes the vacuum of the photon fieldand leads to a static electric field [38, 49]. The time-dependent part δj α ( t ) then generates real photons in themode α . This term is also known as a source term inquantum field theory [40], where it generates the parti-cles (here the photons) that are studied. From this per-spective it becomes obvious that instead of using δj α ( t )one could equivalently slightly change the initial stateof the fully coupled system by adding incoming pho-tons that then scatter off the coupled light-matter groundstate [42]. B. Linear Response in the Length Gauge
With the Hamiltonian of Eq. (1) in length gaugewe can then in principle solve the corresponding time-dependent Schr¨odinger equation (TDSE) for a giveninitial state of the coupled matter-photon systemΨ ( r σ , ..., r N σ N , q , ..., q M )i (cid:126) ∂∂t Ψ( r σ , ..., t ) = ˆ H ( t )Ψ( r σ , ..., t ) , (2)where σ correspond to the spin degrees-of-freedom. How-ever, instead of trying to solve for the infeasible time-dependent many-body wave function, we restrict our-selves to weak perturbations δv ( r , t ) and δj α ( t ) and as-sume that our system is in the ground state of the cou-pled matter-photon system initial time. In this case,first-order time-dependent perturbation theory can beused to approximate the dynamics of the coupled matter-photon system (for details see supplemental materialS2). This framework gives us access to linear spec-troscopy, e.g., the absorption spectrum of a molecule.Traditionally, if we made a decoupling of light and mat-ter, i.e., we assumed Ψ ( r σ , ..., r N σ N , q , ..., q M ) (cid:39) ψ ( r σ , ..., r N σ N ) ⊗ ϕ ( q , ..., q M ), we would only con-sider the matter subsystem ψ (the photonic part ϕ wouldbe completely disregarded). Physically, we would in- vestigate the classical dipole field that the electrons in-duced due to a classical external perturbation δv ( r , t ).To determine this induced dipole field we would only con-sider the linear response of the density operator ˆ n ( r ) = (cid:80) Ni =1 δ ( r − r i ) which would be given by the usual density-density response function in terms of the electronic wavefunction ψ only [50]In this work however, since we do not assume the de-coupling of light and matter, the full density-density re-sponse is taken with respect to the combined ground-state wave function Ψ and is consequently different tothe traditional density-density response. Further, sincewe can also perturb the photon field in the cavity by δj α ( t ) which will subsequently induce density fluctua-tions, the density response δn gets a further contributionleading to δn ( r t ) = (cid:90) dt (cid:48) (cid:90) d r (cid:48) χ nn ( r t, r (cid:48) t (cid:48) ) δv ( r (cid:48) t (cid:48) ) (3)+ M (cid:88) α =1 (cid:90) dt (cid:48) χ nq α ( r t, t (cid:48) ) δj α ( t (cid:48) ) . Here the response function χ nn ( r t, r (cid:48) t (cid:48) ) corresponds to thedensity-density response but with respect to the coupledlight-matter ground state and χ nq α ( r t, t (cid:48) ) corresponds tothe response induced by changing the photon field. Inthe standard linear-response formulation, due to the de-coupling ansatz, changes in the transversal photon fieldwould not induce any changes in the electronic subsys-tem. Since obviously we now have a cross-talk betweenlight and matter, we accordingly have also a genuinelinear-response of the quantized light field δq α ( t ) = (cid:90) dt (cid:48) (cid:90) d r (cid:48) χ q α n ( t, r (cid:48) t (cid:48) ) δv ( r (cid:48) t (cid:48) ) (4)+ M (cid:88) α (cid:48) =1 (cid:90) dt (cid:48) χ q α q α (cid:48) ( t, t (cid:48) ) δj α (cid:48) ( t (cid:48) ) , where χ q α n ( t, r (cid:48) t (cid:48) ) is the full response of the photons dueto perturbing the electronic degrees, and χ q α q (cid:48) α ( t, t (cid:48) ) is thephoton-photon response function. The response func-tion χ q α n ( t, r t (cid:48) ) is in general not trivially connected to χ nq α ( r t, t (cid:48) ), due to the different time-ordering of t and t (cid:48) .The entire linear-response in non-relativistic QED forthe density and photon coordinate can also be writtenin matrix form [51]. In this form we clearly see thatthe density response of the coupled matter-photon sys-tem depends on whether we use a classical field δv ( r , t ),photons, which are created by δj α ( t ), or combinationsthereof for the perturbation. Furthermore, we can alsodecide to not consider the classical response of the cou-pled matter-photon system due to δn ( r , t ), but ratherdirectly monitor the quantized modes of the photon field δq α ( t ). This response yet again depends on whether wechoose to use a classical field δv ( r , t ) that induces pho-tons in mode α or whether we directly generate thosephotons by an external current δj α ( t ). And we also seethat the different modes are coupled, i.e., that photonsinteract. Similarly as charged particles interact via cou-pling to photons, also photons interact via coupling tothe charged particles. Keeping the coupling to the pho-ton field explicitly therefore, on the one hand, changesthe standard spectroscopic observables, and on the otherhand also allows for many more spectroscopic observablesthan in the standard matter-only theory. C. Maxwell-Kohn-Sham linear-response theory
The problem of this general framework in practice isthat already in the simplified matter-only theory we usu-ally cannot determine the exact response functions of amany-body system. The reason is that the many-bodywave functions, which we use to define the response func-tions, are difficult, if not impossible to determine beyondsimple model systems. So in practice we need a differ-ent approach that avoids the many-body wave functions.Several approaches exist that employ reduced quantitiesinstead of wave functions [52–54]. The workhorse of thesemany-body methods is DFT and its time-dependent for-mulation TDDFT [55–57]. Both theories have been ex-tended to general coupled matter-photon systems withinthe framework of QED [13, 36–38, 58].QEDFT allows us to solve instead of the TDSE equiv-alently a non-linear fluid equation for the charge density n ( r , t ) coupled non-linearly to the mode-resolved inho-mogeneous Maxwell’s equation [36–38, 59]. While theseequations are in principle easy to handle numerically,we do not know the forms of all the different termsexplicitly in terms of the basic variables of QEDFT,i.e. ( n ( r , t ) , q α ( t )). To find accurate approximationsone then employs the Kohn-Sham (KS) scheme, wherewe model the unknown terms by a numerically easyto handle auxiliary system in terms of wave functions.The simplest approach is to use non-interacting fermionsand bosons which lead to a similar set of equations,which are however uncoupled. Enforcing that both givethe same density and displacement field dynamics givesrise to mean-field exchange-correlation (Mxc) potentialsand currents [49, 60, 61]. Formally this Mxc poten-tial and current is defined as the difference of the po-tential/current that generate a prescribed internal pairin the auxiliary non-interacting and uncoupled system( v s ([ n ] , r , t ) , j s α ([ q α ] , t )) and the potential/current thatgenerates the same pair in the physical system defined byEq. (1) which we denote by ( v ([ n, q α ] , r , t ) , j α ([ n, q α ] , t )),i.e., v Mxc ([ n, q α ] , r , t ) = v s ([ n ] , r , t ) − v ([ n, q α ] , r , t ) , (5) j α, M ([ n ] , t ) = j s α ([ q α ] , t ) − j α ([ n, q α ] , t ) (6)= − ω α λ α · R ( t ) . In the time-dependent case we only have a mean-fieldcontribution to the Mxc current [36, 38] where the to-tal dipole moment is written as R ( t ) = (cid:82) d r e r n ( r , t ). Further, we have ignored the so-called initial-state de-pendence because we assume (for notational simplicityand without loss of generality) in the following that wealways start from a ground state [61, 62] of the matter-photon coupled system. In this way we can recastthe coupled Maxwell-quantum-fluid equations in termsof coupled non-linear Maxwell-KS equations for auxil-iary electronic orbitals, which sum to the total density (cid:80) i | ϕ i ( r , t ) | = n ( r , t ), and the displacement fields q α ( t ),i.e., i (cid:126) ∂∂t ϕ i ( r , t ) = (cid:20) − (cid:126) m e ∇ + v KS ([ v, n, q α ] , r , t ) (cid:21) ϕ i ( r , t ) , (7) (cid:18) ∂ ∂t + ω α (cid:19) q α ( t ) = − j α ( t ) ω α + ω α λ α · R ( t ) . (8)Here we use the self-consistent KS potential v KS ([ v, n, q α ] , r , t ) = v ( r , t ) + v Mxc ([ n, q α ] , r , t ) thatneeds to depend on the fixed physical potential v ( r , t ) [61], and instead of the full bosonic KS equa-tion for the modes α we just provide the Heisenbergequation for the displacement field. Although theauxiliary bosonic wave functions might be useful forfurther approximations it is only q α ( t ) that is physicallyrelevant and thus we get away with merely coupledclassical harmonic oscillators, i.e., the mode resolvedinhomogeneous Maxwell’s equation. To highlight theextra self-consistency due to coupling between lightand matter we contrast the traditional electron-onlyKS theory with the Maxwell KS theory in Fig. 1. It isthen useful to divide the Mxc potential into the usualHartree-exchange-correlation (Hxc) potential that weknow from electronic TDDFT and a correction term thatwe call photon-exchange-correlation potential (pxc), i.e., v Mxc ([ n, q α ] , r , t ) = v Hxc ([ n ] , r , t ) + v pxc ([ n, q α ] , r , t ) . Clearly, the correction term v pxc will vanish if we take thecoupling | λ α | to zero and recover the purely electroniccase. Since by construction the Maxwell KS system re-produces the exact dynamics, we also recover the exactlinear-response of the interacting coupled system (see alsosupplemental material S3). We can express this with thehelp of the Mxc kernels defined by the functional deriva-tives of the Mxc quantities f n Mxc ( r t, r (cid:48) t (cid:48) ) = δv Mxc ( r t ) δn ( r (cid:48) t (cid:48) ) , f q α Mxc ( r t, t (cid:48) ) = δv Mxc ( r t ) δq α ( t (cid:48) ) ,g n α M ( t, r (cid:48) t (cid:48) ) = δj α, M ( t ) δn ( r (cid:48) t (cid:48) ) , g q α (cid:48) M ( t, t (cid:48) ) = δj α, M ( t ) δq α (cid:48) ( t (cid:48) ) ≡ . and use the corresponding definitions for the Hxc kernel(that only for the variation with respect to n has a non-zero contribution) and the pxc kernels. We note thatusing Eq. (6) we explicitly find g n α M ( t − t (cid:48) , r ) = − δ ( t − t (cid:48) ) ω α λ α · e r . (9) TDDFT Kohn-Sham: no feedback loop QEDFT Maxwell-Kohn-Sham feedback loop ( 𝜕 𝜕𝑡 + 𝜔 𝛼2 ) 𝐄 𝛼 ( 𝑡 ) = − 𝜕 𝜕𝑡 𝝀 𝜶 ∙ 𝐑(𝒕) ( 𝜕 𝜕𝑡 + 𝜔 𝛼2 ) 𝑞 𝛼 ( 𝑡 )⏞ = 𝒋 𝜶 (𝒕)𝜔 𝛼 + 𝜔 𝛼 𝝀 𝜶 ∙ 𝐑(𝒕) 𝑖ℏ 𝜕𝜕𝑡 𝜑 𝑖 (𝒓, 𝑡) = [− ℏ
2𝑚 𝛁⃗⃗ + 𝑣(𝒓, 𝑡) + 𝑣 Mxc ([𝒏, 𝒒 𝜶 ]; 𝒓, 𝑡)]] 𝜑 𝑖 (𝒓, 𝑡) 𝑖ℏ 𝜕𝜕𝑡 𝜑 𝑖 (𝒓, 𝑡) = [− ℏ
2𝑚 𝛁⃗⃗ + 𝑣(𝒓, 𝑡) + 𝑣 Hxc ([𝒏]; 𝒓, 𝑡)]] 𝜑 𝑖 (𝒓, 𝑡) 𝐄 ⊥𝐞𝐱𝐭 (𝒓, 𝑡) ≈ 𝐄 ⊥𝐞𝐱𝐭 (𝑡) 𝐄 ⊥𝐞𝐱𝐭 (𝒓, 𝑡) ≈ 𝐄 ⊥𝐞𝐱𝐭 (𝑡) FIG. 1. Schematics of the Maxwell KS approach contrasted with schematics of the usual semi-classical KS theory. While inthe semi-classical approach the KS orbitals are used as fixed input into the mode-resolved inhomogeneous Maxwell’s equationin vacuum through the total dipole R ( t ) = (cid:82) d r e r (cid:80) i | ϕ i ( r , t ) | (see also appendix A), in the Maxwell KS framework theinduced field acts back on the orbitals, which leads to an extra self-consistency cycle. and g q α (cid:48) M ( t, t (cid:48) ) vanishes, since j α,M in Eq. (6) has no func-tional dependency on q α . Via these kernels we find with χ nn, s ( r t, r (cid:48) t (cid:48) ) and χ q α q α (cid:48) , s ( t, t (cid:48) ), where χ q α q α (cid:48) , s ( t, t (cid:48) ) ≡ α (cid:54) = α (cid:48) , the uncoupled and non-interacting response func-tions that χ nn ( r t, r (cid:48) t (cid:48) ) = χ nn, s ( r t, r (cid:48) t (cid:48) ) + (cid:90) (cid:90) d x d τ χ nn, s ( r t, x τ ) (cid:18)(cid:90) (cid:90) d τ (cid:48) d y f n Mxc ( x τ, y τ (cid:48) ) χ nn ( y τ (cid:48) , r (cid:48) t (cid:48) ) (10)+ (cid:88) α (cid:90) d τ (cid:48) f q α Mxc ( x τ, τ (cid:48) ) χ q α n ( τ (cid:48) , r (cid:48) t (cid:48) ) (cid:33) ,χ q α q α (cid:48) ( t, t (cid:48) ) = χ q α q α (cid:48) ,s ( t, t (cid:48) ) + (cid:88) β (cid:90) (cid:90) (cid:90) dτ dτ (cid:48) d x χ q α q β,s ( t, τ ) g n β M ( τ, x τ (cid:48) ) χ nq α (cid:48) ( x τ (cid:48) , t (cid:48) ) , (11)and accordingly for the mixed matter-photon response functions χ nq α ( r t, t (cid:48) ) = (cid:90) (cid:90) dτ d x χ nn,s ( r t, x τ ) (cid:32)(cid:90) (cid:90) dτ (cid:48) d y f n Mxc ( x τ, y τ (cid:48) ) χ nq α ( y τ (cid:48) , t (cid:48) ) + (cid:88) α (cid:48) (cid:90) dτ (cid:48) f q α (cid:48) Mxc ( x τ, τ (cid:48) ) χ q α (cid:48) q α ( τ (cid:48) , t (cid:48) ) (cid:33) , (12) χ q α n ( t, r (cid:48) t (cid:48) ) = (cid:88) β (cid:90) (cid:90) (cid:90) dτ dτ (cid:48) d y χ q α q β,s ( t, τ ) g n β M ( τ, y τ (cid:48) ) χ nn ( y τ (cid:48) , r (cid:48) t (cid:48) ) . (13)Here we employed the formal connection between re-sponse functions and functional derivatives χ nn ( r t, r (cid:48) t (cid:48) ) = δn ( r , t ) /δv ( r (cid:48) , t (cid:48) ) as well as χ q α q α (cid:48) ( t, t (cid:48) ) = δq α ( t ) /δj α (cid:48) ( t (cid:48) )and accordingly for the auxiliary system. The Mxc ker-nels correct the unphysical responses of the auxiliarysystem to match the linear response of the interactingand coupled problem. So in practice, instead of the full wave function, what we need are approximationsto the unknown Mxc kernels. Later we will providesuch approximations, show how accurate they performfor a model system and then apply them to real sys-tems. If we decouple light and matter, i.e., Ψ (cid:39) ψ ⊗ ϕ ,and disregard the photon part ϕ (as is usually done inmany-body physics), we recover the response function ofEq. (10) with f q α Mxc ≡
0, and f n Mxc → f n Hxc . The responsefunction, which is calculated with the bare matter ini-tial state ψ , then obeys the usual Dyson-type equationrelating the noninteracting and interacting response inTDDFT [63, 64] with v Mxc ([ n, q α ] , r , t ) → v Hxc ([ n ] , r , t ). D. Excited states as properties of the photon field
Following the above discussion, the usual responsefunctions will change and novel response functions areintroduced if we keep the matter-photon coupling explic-itly. This leads to many exciting consequences. Firstly,we get the completely self-consistent response of the sys-tem including all screening, retardation and other effectsthat become important when either the matter subsys-tem is becoming large [65–68] or when strong-couplingsituations are considered. Since light and matter influ-ence each other non-perturbatively the usual simplifiedapproximations that only treat one part of the systemaccurately become unreliable [28, 29] (see also discus-sion in Sec. IV C). Secondly, due to the matter-mediatedphoton-photon interactions (see appendix A and Fig. 2)the usual Maxwell’s equations are changed. A very inter-esting consequence is that in contrast to a purely classicaltheory we can theoretically distinguish whether a system is perturbed by a free current (that in turn would gen-erate a classical electromagnetic field) or by a free elec-tromagnetic field, e.g., a classical laser pulse. Thirdly,we rectify fundamental failings of standard quantum me-chanics, such as the prediction of infinitely-lived excitedstates. The inclusion of the photon modes introducesthe missing photon bath that leads to finite lifetimes(see appendix B and Sec. IV B). In connection to thisit becomes important that we suddenly have access toa wealth of new observables that describe the photonfield. Most importantly this implies the possibility tocompletely change our perspective of excited states ofatoms and molecules. Indeed, in line with the experi-mental situation where changes in the photon field giveus information on the excited states, we can view excited-state properties as arising from quantum modifications ofthe Maxwell’s equations in matter (cid:18) ∂ ∂t + ω α (cid:19) δq α ( t ) = − δj α ( t ) ω α + ω α λ α · (cid:90) d r e r δn ( r , t ) . The response of the density is then found with help of theresponse functions Eqs. (10)-(13). In the usual case of anexternal classical field δv ( r , t ) and δj α ( t ) = 0 we then findthe induced field by (suppressing detailed dependencieswith (cid:82) d r → (cid:82) and (cid:82) d r (cid:80) α → (cid:80) (cid:82) ) (cid:18) ∂ ∂t + ω α (cid:19) δq α ( t ) = ω α λ α · (cid:90) e r χ nn, s δv + ω α λ α · (cid:90) e r χ nn, s f n Mxc χ nn δv + ω α λ α · (cid:88)(cid:90) e r χ nn, s f q α (cid:48) Mxc δq α (cid:48) . (14)Here the first term on the right-hand side correspondsto the non-interacting matter-response. However, due tothe electron-electron interaction we need to take into ac-count also the self-polarization of interacting matter (sec-ond term). Finally, the third term describes the matter-mediated photon-photon response. The excited statesof the coupled light-matter system are in this descrip-tion changes in the photon field. That this perspectiveis actually quite natural becomes apparent if one consid-ers the nature of the emerging resonances for a real sys-tem (see Fig. 7). These resonances are mainly photonicin nature, as they describe the emission/absorption ofphotons (see appendix B). Let us consider now in moredetail what the terms on the right-hand side of the mod-ified Maxwell’s equations mean physically. First of all, in a matter-only theory the self-consistent solution ofthe Maxwell’s equations together with the response ofthe bare matter-system would correspond approximatelyto the first two terms on the right-hand side (see ap-pendix A). The photon-photon interaction would not becaptured in such an approximate approach. Secondly, tohighlight the physical content of the different terms wecan make the mean-field contributions due to v M ( r t ) = (cid:88) α (cid:18)(cid:90) d r (cid:48) λ α · e r (cid:48) n ( r (cid:48) t ) − ω α q α ( t ) (cid:19) λ α · e r (15)+ (cid:90) d r (cid:48) e n ( r (cid:48) t )4 π(cid:15) | r − r (cid:48) | explicit Maxwell’s equation in matter 𝜕 𝜕𝑡 + 𝜔 𝛼2 𝐄 𝛂𝐭𝐨𝐭 𝑡 = − 𝜕𝜕𝑡 𝜇 𝐣 𝛂 𝒕 + 𝜕𝜕𝑡 𝛌 𝛂 ∙ 𝐑 𝐣 𝛂 , 𝑡 𝜕 𝜕𝑡 + 𝜔 𝛼2 𝐄 𝛂𝐞𝐱𝐭 𝑡 = −𝜇 𝜕𝜕𝑡 𝐣 𝛂 𝒕 𝐑 𝐣 𝛂 , t = e 𝐫 𝜒 n,s𝑛 + 𝜒 n,s𝑛 𝑓 Hxc 𝜒 𝑛𝑛 𝐫 ′ ∙ 𝐄 𝛂𝐞𝐱𝐭 The external field corresponding to 𝐣 𝛂 satisfy the equation The dipole in linear-response TDDFT is 𝜕 𝜕𝑡 + 𝜔 𝛼2 𝑞 𝛼 𝑡 = − 𝐣 𝛂 𝒕 𝜔 𝛼 + 𝜔 𝛼 𝝀 𝜶 ∙ 𝐑 𝐣 𝛂 , 𝑞 𝛼 , 𝑡 𝐑 𝐣 𝛂 , 𝑞 𝛼 , t = e 𝐫 𝜒 n,s𝑛 𝑓 Mxc𝑛 𝜒 q α 𝑛 𝐣 𝛂 𝑡 𝛼 + e 𝐫 𝜒 n,s𝑛 𝑓 Mxcq α′ 𝑞 𝛼 ′ 𝛼 Self-consistent dipole in linear-response QEDFT is
Dependence of R on 𝑞 𝛼 shows self-consistency in Maxwell’s equation 𝛁 ∙ 𝐃 ⊥ 𝐫, 𝑡 = 0; 𝐃 ⊥ 𝐫, 𝑡 = 𝜖 𝐄 ⊥ 𝐫, 𝑡 + 𝐏 ⊥ 𝐫, 𝑡 𝜕 𝜕𝑡 − 𝛁 𝐄 ⊥ 𝐫, 𝑡 = −𝜇 𝜕 𝜕𝑡 𝐏 ⊥ 𝐫, 𝑡 ; 1𝑐 𝜕 𝜕𝑡 − 𝛁 𝐃 ⊥ 𝐫, 𝑡 = −𝛁 𝐏 ⊥ 𝐫, 𝑡 Non-self-consistent Maxwell field Self-consistent Maxwell field In the dipole approximation
FIG. 2. Schematics that contrasts the usual Maxwell’s equation (left) with the fully self-consistent Maxwell’s equation (right).Top: The induced transversal electric field E ⊥ as a consequence of the induced polarization P ⊥ , which can be equivalentlyexpressed in terms of the auxiliary displacement field D ⊥ . Left: mode-resolved non-self-consistent Maxwell’s equation with nobackreaction. The external charge current j α induces the external electric field in E tot α = E α + E ext α which acts as an externalperturbation through the dipole. Since the constituents of ˜ χ nn expressed in TDDFT are purely electronic, the induced fielddoes not couple back to the Maxwell field. Right: self-consistent Maxwell’s equation in which j α induces the internal field q α ( t ) through the electron-photon correlated dipole which has an explicit dependence as seen in the QEDFT form of χ nq α .The self-consistency of the induced field through the dipole introduces nonlinearities in the coupled system thus changes theMaxwell field at the level of linear-response. (cid:18) ∂ ∂t + ω α (cid:19) δq α ( t ) = ω α λ α · (cid:90) e r χ nn, s δv + ω α λ α · (cid:90) e r χ nn, s (cid:34) e π(cid:15) | r (cid:48) − r (cid:48)(cid:48) | + (cid:88) α (cid:48) ( λ α (cid:48) · e r (cid:48)(cid:48) ) λ α (cid:48) · e r (cid:48) (cid:35) χ nn δv − ω α λ α · (cid:88)(cid:90) e r χ nn, s ( ω α (cid:48) λ α (cid:48) · e r (cid:48) ) δq α (cid:48) + ω α λ α · (cid:90) e r χ nn, s f n xc χ nn δv + ω α λ α · (cid:88)(cid:90) e r χ nn, s f q α (cid:48) xc δq α (cid:48) . The second term on the right-hand side then correspondsto the random-phase approximation (RPA) to the instan-taneous matter-matter polarization. Here a new termthat corresponds to the dipole self-energy induced by thecoupling to the photons arises. The third term on theright hand side is the RPA approximation to the dipole-dipole mediated photon interaction. To give these termsfurther physical meaning note that in the usual pertur-bative derivation of the van-der-Waals interaction [41]the first two terms would cancel and leave the photonicdipole-dipole interaction that gives rise to the R − for small distances and the R − for larger distances. The restare exchange-correlation (xc) contributions that arise dueto more complicated interactions among the electronsand photons. The last term effectively describe photon-photon interactions mediated by matter. In addition,we want to highlight that xc contributions are directlyresponsible for multi-photon effects, such as two-photonor three-photon processes (see Fig. 4). If we only keepthe mean-field contributions of the coupled problem, wewill denote the resulting approximation in the followingas photon RPA (pRPA) to distinguish it from the bare FIG. 3. Two-level system (with excitation ω ) coupled to onemode of the radiation field (with frequency ω c ). The mattersubsystem is driven by an external classical field v ( t ) and thephoton mode is driven by an external classical current j ( t )and both subsystems are coupled with a coupling strength λ . RPA of only the Coulomb interaction. We see how theMaxwell’s equations in matter change for bound charges,i.e., fields due to the polarization of matter, only. A newterm, the photon-photon interaction, appears. For freecharges, i.e., due to an external charge current δj α ( t ), wesee similar changes. Clearly, if we would not have a cou-pling to matter, then there would be no induced densitychange and we just find the vacuum Maxwell’s equationscoupled to an external current for the electric field. Inother terms, the displacement field trivially correspondsto the electric field (see appendix A). III. EXAMPLES FOR THE COUPLEDMATTER-PHOTON RESPONSE
In this section, we discuss the new perspective enabledby the linear response formalism of QEDFT in more de-tail for a simple and illustrative model system. We dis-cuss a slight generalization of the Rabi model [69, 70],which is the standard model of quantum optics. TheRabi model describes a single electron on two latticesites/energy levels interacting with a single photon mode.We schematically depict the system in Fig 3 and presentall further details of this system in appendix D.First, let us analyze the optical spectra for such a sys-tem and scrutinize the different approximations to theMxc kernels. We will compare the numerical exact re-sults, with the mean-field (pRPA) and the rotating-waveapproximation (RWA). In Fig. 4 (a), (b) and (c) we seehow the optical spectra of the resonantly coupled system(i.e. δ = ω − ω c = 0) change for an increasing electron-photon coupling strength λ . Already for small coupling,the splitting of the electronic state into an upper andlower polariton becomes apparent. Approximately thesestates are given in terms of the RWA as | + , (cid:105) and |− , (cid:105) .The difference in energy between the lower and upperpolariton is called the Rabi splitting Ω R and is used toindicate the strength of the matter-photon coupling. Inmolecular experiments values of up to Ω R /ω c (cid:39) .
25 have been measured [71, 72]. Up to λ = 0 . λ = 0 . . ω c ) the approxima-tions do not recover the exact results. Increasing furtherleads then to not only a disagreement in transition fre-quencies but also the weights of the transitions becomeincreasingly different.Besides a simple check for the approximations to theMxc kernels, the extended Rabi model also allows us toget some understanding of the novel response functions χ σ x q , χ qσ x and χ qq , where σ x is the expectation-value ofthe corresponding Pauli matrix and describes the den-sity/occupation changes between the two sites/energylevels. This means, we consider mixed spectroscopic ob-servables where we perturb one subsystem and then con-sider the response in the other. We analogously employ χ qσ x ( ω ) and χ σ x q ( ω ), respectively, to determine a “mixedpolarizability” (see supplemental material S5). If weplot this mixed spectrum (see Fig. 4 (c) displayed indotted-red for the numerically exact case), we find thatwe have positive and negative peaks. Indeed, this high-lights that excitations due to external perturbations canbe exchanged between subsystems, i.e., energy absorbedin the electronic subsystem can excite the photonic sub-system and vice versa. The oscillator strength of thephotonic spectrum (based on χ qq ) in Fig. 4 (b) providesus with a measure of how strong the displacement field(and with this also the electric field) reacts to an exter-nal classical charge current with frequency ω . Similarly,the mixed spectrum (based on χ σ x q or χ qσ x ) in Fig. 4 (c)provides us with information of how strong one subsys-tem of the coupled system reacts upon perturbing theother one. The oscillator strength here is not necessar-ily positive. What is absorbed by one subsystem can betransferred to the other.In Fig. 4 (d) and (e), we show specifically the absorp-tion spectra of the Rabi model for ultra-strong coupling,i.e., λ = 0 .
7. In this regime, three new peaks arise forthe exact case accounting for high-lying excited stateswith non-vanishing dipole moments due to the strongelectron-photon coupling. The new absorption peaks inFig. 4 (d), also shown in the inset, describes the resonantcoupling case which the RWA and pRPA fail to capture instrong coupling, since processes beyond one-photon areinvolved. Similarly, Fig. 4 (e) depicts the case were thefield is half-detuned from the electronic resonance indi-cating a two-photon process. Clearly in ultra-strong cou-pling the absorption peaks are merely shifted close to the
FIG. 4. Linear-response spectra for the extended Rabi model (dotted-red) compared to the pRPA (dashed-blue) and RWA(full-orange) approximations and for different coupling strengths λ . (a.) Absorption spectra due to matter-matter response,(b.) spectra due to photon-photon response, (c.) spectra due to matter-photon or photon-matter response. (d.) The case for λ = 0 . ω = 2 and ω c = 1 with strength and energies shifted to frequencies favoring2-photon processes. The insets in (d.) and (e.) zoom into the frequency axis showing many-photon process. bare frequencies of the individual subsystems, but remaindressed by the photon field as new peaks arise due to thecoupling. The pRPA and RWA capture the first of thetwo peaks around ω = 2, which is also the frequency ofthe atom, but fail to capture higher lying non-vanishingcontributions to the spectra. These higher-lying peakscorrespond to multi-photon processes. With more accu-rate approximation for the xc potential results closer tothe exact ones can be obtained. We note at this pointthat the peaks in Fig. 4 are artificially broadened andin reality correspond to sharp transitions due to excitedstates with infinite lifetimes. How to get lifetimes quan-titatively will be discussed in the next section. IV. COUPLED MATTER-PHOTON RESPONSE:REAL SYSTEMS
In this section, we apply the introduced formalism inpRPA approximation to real systems. We make thelinear-response formulation practical by reformulatingthe problem as an eigenvalue equation in the frequency-
FIG. 5. Schematic of absorption spectroscopy in optical cav-ities: Benzene (C H ) molecule and λ α denotes the polariza-tion direction of the photon field. domain. For electron-only problems this formulation isknown as the Casida equation [64]. We refer the reader toappendix C for a derivation of our extension of the Casidaequation, which includes transverse photon fields.For the following discussion, we consider benzene0molecules in an optical cavity. In Fig. 5 we schemati-cally depict the experimental setup for a photoabsorp-tion experiment under strong light-matter coupling fora single molecule. First we study the prototypical cav-ity QED setup where a molecule is strongly coupled toa single cavity mode of a high-Q cavity. In the secondsetup, we lift the restriction of only one mode and insteadcouple the benzene molecule to many modes that sam-ple the electromagnetic vacuum field without enhancingthe coupling to a specific mode by hand. In the thirdsetup, we study the behavior of two molecules in an op-tical cavity, as well as a dissipative situation, where onlya few modes are strongly coupled, embedded in a quasi-continuum of modes. In the last example, we analyzethe strong coupling of a single molecule to a continuum ofmodes. We find a transition from Lorentzian lineshape toa Fano lineshape [73] for increasing electron-photon cou-pling strength. These different setups provide us withthe first ab-initio calculation for the spectrum of a realmolecule in a high-Q cavity, the first ab-initio determina-tion of intrinsic lifetimes and the first ab-initio calcula-tion of the non-perturbative interplay between electronicstructure, lifetime and strong-coupling. The two last sit-uations need a self-consistent treatment of photons andmatter alike and cannot be captured by any availableelectronic-structure or quantum-optical method. All ofthose examples highlight the novel possibilities and per-spectives that the QEDFT framework provides.
A. Strong light-matter coupling
The first results we discuss are a set of calculations,where a benzene molecule is strongly coupled to a singlephoton mode in an optical high-Q cavity. We have imple-mented the linear-response pseudo-eigenvalue equation ofEq. (C16) into the real-space code OCTOPUS [74, 75]and details of the numerical parameters are given in ap-pendix E [76].In the first calculation, we include a single cavitymode in resonance to the Π-Π ∗ transition of the benzenemolecule [74, 77], i.e., ω α = 6 .
88 eV. For the light-mattercoupling strength λ α = | λ α | , we choose five different val-ues, i.e. λ α = (0 , . , . , . , .
09) eV / /nm thatcorrespond to a transition from the weak to the strong-coupling limit and the cavity mode is assumed to be po-larized along the x-direction. In Fig. 6, we show theabsorption spectra for these different values of λ α . Westart by discussing the λ α = 0 case that is shown in black.This spectrum corresponds to a calculation of the ben-zene molecule in free space and the spectrum is withinthe numerical capabilities identical to Ref. [74] [78]. Westress that here the broadening of the peaks is only doneartificially since the photon bath is not included in thecalculation. In the examples of Sec. IV B and IV D weinclude many modes and hence sample the photon bathnon-perturbatively. We tune the electron-photon cou-pling strength λ α in Fig. 6. We find for increasing cou- ¯ hω (eV) s t r e n g t h f un c t i o n ( a r b . u . ) λ α = 0 λ α = 2 . λ α = 5 . λ α = 8 . λ α = 11 . Rabi splittingAbsorption Spectra (Benzene)
FIG. 6. Absorption spectra for the benzene molecule in freespace (black) and under strong light-matter coupling in anoptical cavity to ultra-strong coupling (blue). The value for λ α is given in units of [eV / /nm]. pling strength a Rabi splitting of the Π-Π ∗ peak into twopolaritonic branches. The lower polaritonic branch hashigher intensity, compared to the upper polaritonic peak.Numerical values for the excitation energy E I , the tran-sition dipole moment x I and the oscillator strength f I are given in Tab. I in the appendix. This demonstratesthat ab-initio theory is able to describe excited-stateproperties of strong light-matter coupling situations andcaptures the hybrid character of the combined matter-photon states. Thus predictive theoretical first-principlecalculations for excited-states properties of real systemsstrongly coupled to the quantized electromagnetic fieldare now available. This will allow unprecedented insightsinto coupled light-matter systems, since we have access tomany observables that are not (or not well [29]) capturedby quantum-optical models. B. Lifetimes of excitations from first principles
Next we consider how to obtain lifetimes from QEDFTlinear-response theory. In this example, we explicitlycouple the benzene molecule to a wide range of pho-ton modes similar as in the spontaneous emission cal-culation of Ref. [79]. While in Ref. [79], the system wassimulated with 200 photon modes, we choose here now80.000 photon modes. The energies of the sampled pho-ton modes cover densely a range from 0 .
19 meV, for thesmallest energy up to 30.51 eV for the largest one with aspacing of ∆ ω = 0 .
38 meV. However, we do not samplethe full three-dimensional mode space together with thetwo polarization possibilities per mode but rather con-sider a one-dimensional slice in mode space. This one-dimensional sampling of mode frequencies will changethe actual three-dimensional lifetimes, but for demon-strating the possibilities of obtaining lifetimes this is suf-ficient [80]. The sampling of the photon modes corre-sponds to the modes of a quasi-one dimensional cav-1ity. We choose a cavity of length L x [79] in x -directionwith a finite width in the other two directions that aremuch more confined. Thus we employ ω α = αcπ/L x and λ α = (cid:113) (cid:126) (cid:15) L x L y L z sin( ω α /c x ) e x , where x = L x / x -direction. While wehave a sine mode function in the x -direction, we assume aconstant mode function in the other directions. For thisexample, we choose a cavity of length L x = 3250 µ m in x -direction, L y = 10 . A in y -direction and L z = 2 . A in z -direction.The results of this calculation are shown in Fig. 7. In hω (eV)00.040.080.120.16 s t r e n g t h f un c t i o n ( e V ) (a) Absorption Spectrum (Benzene)6.8 6.85 6.9 6.95¯ hω (eV)00.040.080.120.16 s t r e n g t h f un c t i o n ( e V ) ∆ E (b) Π − Π ∗ transition 0 . . . . . . hω (eV)00.050.10.15 electronicphotonic(c) σ − σ + transition FIG. 7. First principles lifetime calculation of the electronicexcitation spectrum of the benzene molecule in an quasione-dimensional cavity: (a) Full spectrum of the benzenemolecule, (b) zoom to the Π − Π ∗ transition, where the blackarrow indicates the full width at half maximum (FWHM) ∆ E ,(c) zoom to a peak contributing to the σ − σ + transition. Thegray spectrum is obtained by Wigner-Weisskopf theory [31].The dotted spectral data points correspond to many coupledelectron-photon excitation energies which together comprisethe natural lineshape of the excitation. Blue color refers toa more photonic nature of the excitations, vs. red color to amore electronic nature. Fig. 7 (a) we show the full spectrum. The electron-photon absorption function that has been obtainedby coupling the benzene molecule to the quasi one-dimensional cavity with 80.000 cavity modes is plottedin blue. Since we have sampled the photon part densely,we do not need to artificially broaden the peaks anymore.Formulated differently, we can directly plot the oscilla-tor strength and the excitation energies of our resultingeigenvalue equation and do not need anymore to employthe Lorentzian broadening. In Fig. 7 from blue (morephotonic) to red (more electronic) for the electron-photon absorption spectrum we plot the different contributionsof each pole in the response function. These results con-firm our intuition that resonances are mainly photonic innature and that a Maxwell’s perspective of excited statesis quite natural. In (b) we zoom to the Π-Π ∗ transition.Due to quasi one-dimensional nature of the quantizationvolume, we find a broadening of the peak that is largerthan it is for the case of a three-dimensional cavity dueto the sampling of the electromagnetic vacuum. This issimilar to changing the vacuum of the electromagneticfield. Accordingly the lifetimes of the electronic statesare shorter if the electromagnetic field is confined to onedimension and we will discuss this in the next section. C. Connection to standard Wigner-Weisskopftheory
If the coupling between light and matter is very weakand neither subsystem gets appreciably modified dueto the other, in contrast to the previous strong light-matter coupling case, the radiative lifetimes of atomsand molecules can be calculated using the perturbativeWigner-Weisskopf theory [31] in single excitation approx-imation, as well as under the assumption of the Markovapproximation. These approximations are justified in theusual free-space case, where the results of Wigner andWeisskopf reproduce the prior results of Einstein basedon the ad-hoc A and B coefficients. However it does notinclude the treatment of ensembles of molecules that ef-fectively enhance the matter-photon coupling strength,as shown below. Under the assumption of Wigner-Weisskopf theory, the radiative decay rate is given byΓ D = ω | d | π(cid:15) (cid:126) c . (16)For a one-dimensional cavity in x-dimension the resultschange to [32] Γ D = ω | d | L y L z (cid:15) (cid:126) c (17)For comparison, we show in Fig. 7 in grey the peaksthat are predicted by Wigner-Weisskopf theory. Sinceour sampling is very dense, we find for both peaks shownin the bottom a good agreement with Eq. 17.In fact, if we take the continuum limit for the photonmodes, we recover in our framework the lifetimes pre-dicted by Wigner-Weisskopf theory including the diverg-ing energy shifts [81], i.e. the Lamb shift. Due to theLamb shift, our resulting peaks are slightly shifted, dueto the divergencies. These divergencies can be handledby renormalization theory. The lifetimes can now be ob-tained the following way: We measure the full width athalf maximum (FWHM), indicated by the black arrowin (b). In this case, we find ∆ E FWHM = 0 . τ Π − Π ∗ follows by τ Π − Π ∗ = (cid:126) / ∆ E FWHM = 32 .
27 fs. Using the Wigner-Weisskopf for-mula from Eq. 17, and the dipole moments and energies2from the LDA calculation without a photon field, we finda lifetime of 32 .
21 fs. As a side remark, the same tran-sition using Eq. 16 has a free-space lifetime of 0 .
89 ns,roughly in the range of the 2p-1s lifetime of the Hydrogenatom of 1 . σ − σ + transition. We find a narrow ab-initio peak that is notas well sampled as the Π − Π ∗ . We note in passing thatwe find a ionization energy of 9 .
30 eV using ∆-SCF inthe benzene molecule with the LDA exchange-correlationfunctional. We note in passing that we find a ionizationenergy of 9 .
30 eV using ∆-SCF in the benzene moleculewith the LDA exchange-correlation functional. In oursimulation, coupling to peaks higher than the ionizationenergy are broadened by continuum (box) states.
D. Beyond the single molecule limit anddissipation in QEDFT hω (eV)00.040.080.120.16 s t r e n g t h f un c t i o n ( e V ) hω (eV) s t r e n g t h f un c t i o n ( a r b . u . ) (b) Mulit mode - Rabi Splitting (Benzene) FIG. 8. (a) Two molecules of benzene strongly coupled to80.000 cavity modes of an one-dimensional cavity. The furtherapart the molecules are, the closer the peak gets to the singlemolecule peak. Also we notice the doubled peak broadening(shorter lifetime). The gray spectrum is obtained by Wigner-Weisskopf theory [31]. (b) We show the Rabi splitting in asituation of a single strongly coupled mode with 80.000 cavitymodes (green), and three strongly coupled modes with 80.000cavity modes (blue). The red lines correspond to the samesetup as in (a). The dashed lines refer to the frequency ofthe cavity modes. The peaks become broadened due to theinteraction with the continuum.
In contrast to the free-space result, where weak cou-pling as well as the assumption of a dilute gas of molecules are implied, in the case of single-moleculestrong coupling [8] or when nearby molecules or an en-semble of interacting molecules modify the vacuum, theusual perturbative theories break down. Changes inthe electronic and the photonic subsystem become self-consistent and the usual distinction of light and matterbecomes less clear. In such situations the linear-responseformulation of QEDFT as well as the Maxwell’s per-spective of excited-state properties becomes most pow-erful. Consider, for instance, two benzene moleculesweakly coupled to a one-dimensional continuum of pho-ton modes. If the molecules are far apart we just findthe usual Wigner-Weisskopf result. But if we bring themolecules closer (see Fig. 8 (a)), we see that the com-bined resonance shifts and the combined linewidth be-comes broader, implying a shortened lifetime. In Fig. 8(b), we consider the case of single-molecule strong cou-pling, where a few out of the 80.000 modes have an en-hanced coupling strength. In red, we show the spectrumwhere the molecule is coupled to the continuum, as isalso shown in Fig. 7. We then introduce a single stronglycoupled mode at the Π − Π ∗ transition energy and theresulting spectra is shown in green. We note that in thefigure, the cavity frequencies are plotted in dashed lines.The single mode introduces the expected Rabi splittinginto the upper and lower polariton and the peaks of theupper and lower polariton become broadened due to theinteraction with the continuum. Interestingly, we finda different line broadening for the lower and the upperpolaritonic peak, since only the sum of both has to beconserved. The smaller broadening for these two lowerpolaritonic states implies that the radiative lifetime ofthe lower and upper polaritonic state is longer than thelifetime of the excitation in weakly-coupled free-space.In blue, we show the spectra, where we have introducedthree strongly coupled modes in addition to the cavity80.000 modes of the continuum. We tune the two addi-tional cavity modes in resonance to the lower and upperpolariton peak of the green plot. We find additional peaksplitting, but also a shifting of peak positions at 7.8 eV.In the last numerical example, we study the strong cou-pling to the continuum for the case of a single molecule.The results are shown in Fig. 9. Here, we effectively en-hance the light-matter coupling strength by reducing thevolume of the cavity along the y and z direction. Forcomparison, we show in red the setup that is also shownin Fig. 7, where the excitations have Lorentzian line-shape consistent with Wigner-Weisskopf theory as dis-cussed in the previous section. By gradually reducingthe dimensions along the y and z direction, we find dras-tic changes in the lineshape of the excitations. Thesechanges lead to the transition of the lineshape from aLorentzian to a Fano lineshape, as becomes clearly visi-ble for L x L z = 0 . A .As a summary, we have presented in this section, thatlineshapes, as well as lifetimes can be inferred directlyfrom first principle calculations. In case of Lorentzianlineshape, we find that the width of the calculated peaks3 hω (eV) s t r e n g t h f un c t i o n ( a r b . u ) L y L z = 28 . A L y L z = 4 . A L y L z = 1 . A L y L z = 0 . A Changing the continuum (Benzene)
FIG. 9. Ab-intio lifetime calculation of the electronic excita-tion spectrum of the benzene molecule in an one-dimensionalcavity along x -direction with different length in L y and L z direction. The red spectra refer to the same setup as inFig. 7. Effectively the electron-photon strength increaseswith smaller L y and L z length leading to a transition from aLorentzian lineshape to a Fano lineshape. (no need to introduce any artificial broadening as com-monly done) correspond to the lifetimes. These calcula-tions demonstrate that ab-initio theory is able to capturethe true nature of excitations, i.e., resonances with finiteintrinsic lifetimes, without the need of an artificial bathor post-processing. This allows a new perspective ofwell-known results. Furthermore, we find that the ex-citations measured in absorption/emission experimentsare mainly photonic in nature, and it is only the peakposition that is dominated by the matter constituents.This is of course very physical, since what we see is theabsorption/emission of a photon, not of the matter con-stituents. Further, since we describe the photon vacuumon the same theoretical footing as the matter subsys-tem, we have full control over the photon field making itstraightforward to simulate very intricate changes, e.g.,changing the character of a specific mode out of basi-cally arbitrarily many, and investigating its influence onexcited-states properties such as the radiative lifetime.This allows predictive first-principle calculations for in-tricate experimental situations similar to the ones en-countered in Ref. [33, 34]. V. SUMMARY AND OUTLOOK
In this work we have introduced linear-response the-ory for non-relativistic quantum-electrodynamics in thelong wavelength limit. Compared to the conventionalmatter-only response approaches, we have highlightedhow in the coupled matter-photon case the usual re-sponse functions change, how novel photon-photon andmatter-photon response functions are introduced, howthese novel response functions provide a photonic per-spective on excited state properties, how the results lead to changes in the usual Maxwell’s equation in matterand how we can efficiently calculate all these responsefunctions in the framework of QEDFT. By investigatinga simple model system, we have shown how the spec-trum of the matter subsystem is changed upon couplingto the photon field. Further we have demonstrated therange of validity of a simple yet reliable approximation tothe in general unknown mean-field exchange-correlationkernels. Using this approximation we have presented thefirst ab-initio calculations of the spectrum of real systems(benzene molecules) coupled to the modes of the quan-tized electromagnetic field. In one example we have cal-culated the change upon strong coupling to a single modeof a high-Q cavity, which leads to a large Rabi splitting.In the second example we have calculated from first prin-ciples the natural linewidths of benzene coupled to a spe-cific sampling of the vacuum field. In the last examples,we demonstrated the abilities to calculate many-moleculesystems, as well as dissipative strong-coupling situa-tions, as well as strong coupling to the continuum, wherewe find a transition from Lorentzian lineshape to Fanolineshape, where the usual (perturbative) approaches tolight-matter coupling fail. These results demonstrate theversatility and possibilities of QEDFT, where light andmatter are treated on equal quantized footing. In thecontext of strong light-matter coupling, e.g., in polari-tonic chemistry, the presented linear-response formula-tion allows now to determine polaritonically modifiedspectra from first principles. Together with ab-initioground-state calculations [39] QEDFT now provides aworkable first-principle description to analyze and pre-dict photon-dressed chemistry and material sciences. Inparticular, our novel approach provides a unique prac-tical computational scheme to compute photon-dressedexcited-state potential-energy surfaces and non-adiabaticcoupling elements that are required for ab-initio calcu-lations in the emerging field of polaritonic chemistry.Further, in the context of standard ab-initio theory,the linear-response formulation of QEDFT now allowsthe calculation of intrinsic lifetimes and provides ac-cess to quantum-optical observables. Specifically, due tothe non-perturbative nature of the approach, quantum-optical problems where the self-consistent feedback be-tween light and matter has to be taken into account,e.g., that many molecules change the photon vacuumand hence the Markov approximation breaks down, be-come feasible. For optical physics, the presented linear-response framework presents an interesting opportunityto study the modifications of the Maxwell’s equations inmatter from first principles. Finally we want to high-light that although the QEDFT linear-response frame-work is new, its similarity to the usual matter-only linear-response formulation in terms of an pseudo-eigenvalueproblem makes it very easy to include in already exist-ing first-principle codes. This, together with the abovediscussed novel possibilities in different fields of physics,shows that there are many interesting cases that can bestudied with the presented method.4
VI. ACKNOWLEDGEMENTS
We would like to thank Christian Sch¨afer and No-rah Hoffmann for insightful discussions, and SebastianOhlmann for the help with the efficient massive paral-lel implementation. JF acknowledges financial supportfrom the Deutsche Forschungsgemeinschaft (DFG) un-der Contract No. FL 997/1-1 and all of us acknowledgefinancial support from the European Research Council(ERC-2015-AdG-694097).
Appendix A: Modification of the Maxwell’s equation
In this section, we give more details on the modifi-cations of the Maxwell’s equations. The semi-classicaldescription of light-matter interaction is limited as a re-sult of the transverse field being treated as an externalperturbation. This approximation breaks the feedbackloop between light and matter that leads to apparentchanges in the Maxwell’s equation. Let us start from theclassical description and assume that we are interestedin the induced fields due to an external perturbation.If everything is perfectly classical there is no differencewhether we perturb by an external transversal field a ⊥ or an external classical current j ⊥ due to the inhomoge-neous Maxwell’s equation in vacuum (cid:18) c ∂ ∂t − ∇ (cid:19) a ⊥ ( r , t ) = µ c j ⊥ ( r , t ) . (A1)Now, if we have some theory to relate these externalperturbation to the induced current J ⊥ [ a ⊥ ], the inducedfield reads (cid:18) c ∂ ∂t − ∇ (cid:19) A ⊥ ( r , t ) = µ c J ⊥ ([ a ⊥ ] , r , t ) , (A2)from which we can calculate the induced physical fields,e.g., the transversal electric field in Coulomb gauge is E ⊥ ( r , t ) = − c ∂ t A ⊥ ( r , t ) [82]. We can again combinethese two results and look at the total field A tot ⊥ = a ⊥ + A ⊥ , which obeys (cid:18) c ∂ ∂t − ∇ (cid:19) A tot ⊥ ( r , t ) = µ c ( j ⊥ ( r , t ) + J ⊥ ([ j ⊥ ] , r , t )) . (A3)Using the Maxwell relations once more we can equiva-lently find for, e.g., the induced electric field (cid:18) c ∂ ∂t − ∇ (cid:19) E ⊥ ( r , t ) = − µ ∂∂t J ⊥ ([ a ⊥ ] , r , t ) . (A4)We can now make a connection to the Maxwell’s equa-tion in matter, where the j ⊥ is called the free currentand J ⊥ the bound current. Assuming that we can ex-press the transversal induced current locally around the center of charge as J ⊥ ( r , t ) ≈ ∂∂t P ⊥ ( r , t ), where we usethe polarization P ⊥ ( r , t ) = (cid:15) e M (cid:88) α =1 λ α ( r ) (cid:90) d r (cid:48) λ α ( r (cid:48) ) · r (cid:48) n ([ a ⊥ ] , r (cid:48) , t ) , and expand the electric field in the modes λ α ( r ) as E ⊥ ( r , t ) = M (cid:88) α =1 λ α ( r ) E α ( t ) , (A5)we can rewrite the above equation at the center of charge,i.e., λ α ( r ) → λ α , as (cid:18) ∂ ∂t + ω α (cid:19) E α ( t ) = − ∂ ∂t λ α · R ([ a ⊥ ] , t ) . (A6)Using this kind of approach we can connect δn ( r , t ) ofEq. (3) to the induced electric field δ E ⊥ ( r , t ), where weemploy a spatially homogeneous vector potential a ⊥ ( t )that gives rise to the external electric field E ext ⊥ ( t ) = − c ∂∂t a ⊥ ( t ). In a final step, to avoid solving the abovemode-resolved Maxwell’s equations, one often even ig-nores the spatial dependence of the induced field andmerely uses E α ( t ) = − λ α · R ([ a ⊥ ] , t ). If we now de-termine in linear response R ([ a ⊥ ] , t ) we immediately seethat when χ nn is changed due to strong light-matter cou-pling also the induced field is changed. Furthermore, thereformulation of the linear-response kernel in Eq. (10)shows that we get a feedback from the induced pho-ton field onto the matter. Such intrinsic back-reaction(screening) effects are very important for large systems,as is well known from solid-state physics, where the bare(vacuum) electric field as determined by Eq. (A6) doesnot agree with the measured spectrum. One needs to in-clude the self-consistent polarization of the system thatcounter-acts the external perturbing field. This can bedone approximately in linear response by solving self-consistently a Maxwell’s equation with the matter re-sponse as input [65–68]. In the theory of classical elec-trodynamics, a convenient way to do so is to switch tothe Maxwell’s equations in matter. In the above consid-erations this means we introduce the displacement field D ⊥ = (cid:15) E ⊥ + P ⊥ , where now all the knowledge abouthow the system reacts to an external perturbation is en-coded again in P ⊥ such that we find (cid:18) c ∂ ∂t − ∇ (cid:19) D ⊥ ( r , t ) = −∇ P ⊥ ([ a ⊥ ] , r , t ) . (A7)After expanding D ⊥ ( r , t ) = (cid:15) (cid:80) α ω α λ α ( r ) q α ( t ) andthen performing the long wave-length limit we arrive at (cid:18) ∂ ∂t + ω α (cid:19) q α ( t ) = ω α λ α · R ([ a ⊥ ] , t ) , (A8)which is the classical analogue of Eq. 8. In the usual de-coupled light-matter description without self-consistency5we then simply determine R ([ a ⊥ ] , t ) from the electric per-mittivity and ignore any feedback that describes how thematter system affects (screens) the field. Approximateself-consistency is found once the induced field E ⊥ istaken into account to screen the perturbing field E ext ⊥ .But in our case we want to go beyond this simple ap-proximate self-consistency which will break down oncethe coupling between light and matter is strong. Notethat in the macroscopic Maxwell’s equation the electricfield becomes E α ( t ) = ω α q α ( t ) − λ α · R ([ a ⊥ ] , t ), and wesee that if we ignore the spatial dependence in determin-ing E ⊥ we basically assume D ⊥ = E ⊥ .In our description we keep the photon field as a dy-namical variable of the system such that the Maxwellfield couples to the electronic system, leading to a fullyself-consistent description of the light-matter response.Besides the changes in χ nn , which when used as an in-put into Eqs. (A4) or A8, captures the self-consistentresponse of the light-matter system, we can now also di- rectly access the induced electric field by considering theresponse of the displacement field due to χ qn and the useof ˆ E ⊥ = M (cid:88) α =1 λ α ω α (cid:18) ˆ q α − λ α ω α · R (cid:19) . As discussed in Sec. II D, this leads to a complete changeof perspective, since it highlights that the excited statesof the coupled light-matter system can be viewed aschanges in the quantized Maxwell field in accordance tothe usual experimental situation. On the other hand, wecan now also investigate what the quantum descriptionof the coupled light-matter system does to the Maxwell’sequations. We therefore consider the case where the free(time-derivative of the) current δj α ( t ) is non-zero whilethe external classical field is zero, i.e., δv ( r , t ) = 0. Inthis case, we find (cid:18) ∂ ∂t + ω α (cid:19) δq α ( t ) = − δj α ( t ) ω α + ω α λ α · (cid:90) e r χ nn, s f n Mxc χ nq α δj α + ω α λ α · (cid:88)(cid:90) e r χ nn, s f q α (cid:48) Mxc δq α (cid:48) . (A9)If we contrast this to the classical Maxwell’s equation inmatter (cid:18) ∂ ∂t + ω α (cid:19) δq α ( t ) = − δj α ( t ) ω α + ω α λ α · δ R ([ j ⊥ ] , t ) , (A10) where R ([ j ⊥ ] , t ) would be determined from the responseof the matter system due to the corresponding externalfield a ⊥ , we see that besides the self-consistent responseof the matter system (second term on the right hand side)also a genuine new (matter-mediated) photon-photon in-teraction term (third term on the right hand side) ap-pears. Making again the mean-field explicit leads to (cid:18) ∂ ∂t + ω α (cid:19) δq α ( t ) = − δj α ( t ) ω α + ω α λ α · (cid:90) e r χ nn, s (cid:34) e π(cid:15) | r (cid:48) − r (cid:48)(cid:48) | + (cid:88) α (cid:48) ( λ α (cid:48) · e r (cid:48)(cid:48) ) λ α (cid:48) · e r (cid:48) (cid:35) χ nq α δj α − ω α λ α · (cid:88)(cid:90) e r χ nn, s ( ω α (cid:48) λ α (cid:48) · e r (cid:48) ) δq α (cid:48) + ω α λ α · (cid:90) e r χ nn, s f n xc χ nq α δj α + ω α λ α · (cid:88)(cid:90) e r χ nn, s f q α (cid:48) xc δq α (cid:48) . If we ignore the xc contributions to the matter-photonand photon-photon response we get the pRPA approx-imation to the Maxwell’s equation in matter. In thispRPA form we clearly see how the Maxwell’s equa-tion becomes non-linear because of the feedback betweenlight and matter. Such non-linearities of the Maxwell’sequations are investigated in great detail in high-energyphysics in the context of strong-field QED [83]. Inthat case the strong fields lead to particle creation andthus a matter-mediated photon-photon interaction. In our case, we do not need these high energies becausewe consider the photon-photon interaction due to con-densed matter in form of atoms, molecules or solids anduse, e.g., a cavity to enhance the coupling. That thechanges in the Maxwell’s equations are not purely theo-retical concepts but lead to observable effects can be seenin many physical situations. As mentioned before, themost well-known effect are polarization effects in solid-state systems [68], but more strikingly are effects dueto the quantum-matter-mediated photon-photon interac-6tions, see e.g. Ref. [84]. In this context, the presentedab-initio method allows to theoretically investigate thephoton-photon interactions and possibly predict systemswith very strong photon-photon correlations. In suchcases the strong photon correlations could be used to givecomplementary insights into molecular systems or to im-print the photonic correlations on the matter subsystem.Besides these differences we highlight that the quantizedMaxwell’s equation in matter, if we allow for both, a freeexternal current and a free external field, can indeed dis-criminate between these two sorts of perturbations. In apurely classical theory, due to Eq. A1, there can be nodifference. This provides a completely new playground toinvestigate the difference between classical and quantumphysics.
Appendix B: Novel photonic observables andradiative lifetimes
In the presented framework, besides the above high-lighted changes in, e.g., the Maxwell’s equations, novelobservables become accessible. For instance, one canmonitor the response of the matter system due to aperturbation of the photonic subsystem by an externalcurrent. This allows to investigate directly the cross-correlation between the matter and the photon subsys-tem induced by χ nq α . Also note, that this cross-correlationobservable allows to distinguish between the response dueto a purely classical field δv ( r , t ) or due to a quantizedfield, since δj α ( t ) generates photons (which is equivalentto just use a slightly different initial state with an in-coming photon pulse) that then perturb the correlatedmatter-photon system. This makes the presented frame-work applicable to also determine observables due tonovel spectroscopies that use quantum light [85]. Thisarea of spectroscopy is so far not accessible with commonfirst-principle methods. One further important observ-able that can be captured in this approach is the intrin-sic lifetimes of excited states, which is not accessible instandard matter-only quantum mechanics. Let us brieflyexplain what we mean by this. In standard quantummechanics we find besides the ground state also othereigenstates, i.e., excited states. Hereby an eigenstate is asquare-integrable eigenfunction of a self-adjoint, usuallyunbounded Hamiltonian. If we excite a matter systemfrom its ground state into such an excited state, it willremain in this state as long as we do not perturb it. Inquantum mechanics we then also have generalized eigen-states, so-called scattering states, which are not square-integrable and that constitute the continuous spectrum ofsuch a Hamiltonian [86]. The simplest example is the freeelectronic Hamiltonian ˆ T = (cid:80) Ni =1 − (cid:126) m e ∇ i which in in-finite space has a purely continuous spectrum consistingof non-normalizable plane-waves [87]. The physical inter-pretation of such scattering states - as already the nameindicates - is that particles propagate to infinity and donot stay bound anywhere. Thus exciting a matter system from its ground state into such a generalized eigenstatecorresponds to the physical process of ionization. Ioniza-tion, however, is something completely different than theprocess of spontaneous emission. That is, if we put anatom or molecule into an “excited state”, even withouta further perturbation it will relax to the ground stateby emitting radiation. The time the system stays in this“excited state” before emitting a photon is called the life-time. The process of spontaneous emission clearly cannotbe captured by standard quantum mechanics where mat-ter and light are decoupled. Non-relativistic QED, how-ever, does capture this process [42] by coupling the mat-ter system to the quantized electromagnetic field whichconsists of infinitely many harmonic oscillators. In thisway the excited states of the bare matter system turninto resonances and the ground state (usually) remainsthe only eigenstate of the combined matter-photon sys-tem. While formally these resonances are indeed scatter-ing states of the combined matter-photon system, it isonly the photonic part that shows a scattering behavior,i.e., a photon leaves the vicinity of the matter subsys-tem. The matter subsystem just relaxes to the only sta-ble state, its ground state [42]. In linear response such re-laxation processes express themselves as finite linewidthsof excitations, where the linewidth can be associated withthe lifetimes of the different resonances. In our slightlysimplified treatment based on Eq. (1) we only considera finite number of photon modes, and hence we do nothave genuine resonances. However, by including enoughmodes we sample the influence of the vacuum and in-stead of one sharp transition peak (which numerically isusually artificially broadened) we get many that approx-imate the resonance. In this way linear-response theoryfor non-relativistic QED in the long wavelength limit candetermine lifetimes of real systems. We show an exam-ple for such an ab-initio lifetime calculation in Sec. IV B.This provides a further new field of research that thepresented framework makes accessible for ab-initio the-ory. While it is conceptually very interesting to revisitwell-known results for intrinsic radiative lifetimes of gasphase molecules, since we can now study the nature ofresonances in detail (see, e.g., the discussion on the pho-tonic nature of resonances in Sec. IV B), we have nowaccess to even more exciting experimental situations. Bychanging the environment, e.g., putting the molecule in-side a cavity and thus enhance certain modes while sup-pressing others, one can change and control the radiativelifetimes of single molecules [33, 34] (see also Sec. IV D).We can thus theoretically study and predict realistic ex-perimental situations where non-trivial changes in thephotonic vacuum, e.g., due to nearby surfaces or otherphysical entities, directly influence intrinsic lifetimes andproperties of resonances.7 Appendix C: Linear-response theory as apseudo-eigenvalue problem
In this section, we reformulate the linear-response the-ory of coupled electron-photon systems as a pseudo-eigenvalue problem. The entire linear-response in non-relativistic QED for the density and photon coordinatecan be written in matrix form as δnδq δq ... δq M = χ nn χ nq χ nq . . . χ nq M χ q n χ q q χ q q . . . χ q q M χ q n χ q q χ q q . . . χ q q M ... ... ... . . . ... χ q M n χ q M q χ q M q . . . χ q M q M δvδj δj ... δj M (C1)where we imply integration over time and space when ap-propriate. In this form we clearly see that the density re-sponse of the coupled matter-photon system depends onwhether we use a classical field δv ( r , t ), photons, whichare created by δj α ( t ), or combinations thereof for the per-turbation. The explicit coupling between the subsystems(i.e. matter and photons) demonstrates changes in thesubsystems as a result of the back-reaction between mat-ter and photons. The cross-talk between the respectivecoupled subsystems shows up in the cross-correlation re-sponse functions which leads to changes in the respectiveobservables ( n ( r , t ) , q α ( t )). This becomes evident by con-sidering an external perturbation of the coupled systemwith the external potential δv ( r t ) reduces to the coupledset of responses (cid:40) δn ( r t ) = (cid:82)(cid:82) dt (cid:48) d r (cid:48) χ nn ( r t, r (cid:48) t (cid:48) ) δv ( r (cid:48) t (cid:48) ) ,δq α ( t ) = (cid:82)(cid:82) dt (cid:48) d r (cid:48) χ q α n ( t, r (cid:48) t (cid:48) ) δv ( r (cid:48) t (cid:48) ) . (C2)Here, the cross-correlation response function χ q α n ( t, r (cid:48) t (cid:48) )accounts for the action of the matter subsystem on thephoton field which gives rise to a response of the pho-ton field as a result of perturbing the matter. In thesemi-classical approach in which TDDFT is based on,the cross-correlation response function do not show upbut rather just a simplified form (since there the wavefunction describes only the matter subsystem) of the χ nn ( r t, r (cid:48) t (cid:48) ). Similarly, a perturbation of the coupled sys-tem with the external charge current δj α ( t ) results in (cid:40) δn ( r t ) = (cid:80) Mα =1 (cid:82) dt (cid:48) χ nq α ( r t, t (cid:48) ) δj α ( t (cid:48) ) ,δq α ( t ) = (cid:80) Mα (cid:48) =1 (cid:82) dt (cid:48) χ q α q α (cid:48) ( t, t (cid:48) ) δj α (cid:48) ( t (cid:48) ) . (C3)The cross-correlation response function χ nq α ( r t, t (cid:48) ) ac-counts for the action of the photon field on the matterthus specifying the response of the density by perturbing the photon field and χ q α q α (cid:48) ( t, t (cid:48) ) describes how photon in-teracts via matter thus specifying changes in response ofthe photon field.Next, we need to find an efficient way to solve theselinear-response equations in terms of the Maxwell KS sys-tem. First, performing a Fourier transformation fromtime t and t (cid:48) to frequency space ω and using the Hxc andpxc kernels, we write the response functions of Eqs. (10)-(13) in the following compact notation χ nn = χ nn, s + χ nn, s (cid:34)(cid:0) f n pxc + f n Hxc (cid:1) χ nn + (cid:88) α f q α pxc χ q α n (cid:35) , (C4) χ q α n = (cid:88) β χ q α q β, s g n β M χ nn , (C5) χ nq α = χ nn, s (cid:34)(cid:88) α (cid:48) f q α (cid:48) pxc χ q α (cid:48) q α + (cid:0) f n pxc + f n Hxc (cid:1) χ nq α (cid:35) , (C6) χ q α q α (cid:48) = χ q α q α (cid:48) , s + (cid:88) β χ q α q β , s g n β M χ nq α (cid:48) . (C7)Those equations are coupled with respect to the ex-ternal perturbations as seen in Eqs. (C2)-(C3). Theperturbation with respect to the external potential δv ( r t ) results in a coupled set of response functions { χ nn ( r t, r (cid:48) t (cid:48) ) , χ q α n ( t, r (cid:48) t (cid:48) ) } and for the external current δj α ( t ) gives the coupled set (cid:110) χ nq α ( r t, t (cid:48) ) , χ q α q α (cid:48) ( t, t (cid:48) ) (cid:111) .These pairs of coupled response functions have to besolved in a self-consistent way to obtain the exact in-teracting response functions. The response functions ofEqs. (C2) and (C3) can be expressed in frequency spacethrough a Fourier transform that yields δn v ( r , ω ) = (cid:90) d r (cid:48) χ nn ( r , r (cid:48) , ω ) δv ( r (cid:48) , ω ) , (C8) δq α,v ( ω ) = (cid:90) d r (cid:48) χ q α n ( ω, r (cid:48) ) δv ( r (cid:48) , ω ) , (C9) δn j ( r , ω ) = (cid:88) α χ nq α ( r , ω ) δj α ( ω ) , (C10) δq α,j ( ω ) = (cid:88) α (cid:48) χ q α q α (cid:48) ( ω ) δj α (cid:48) ( ω ) . (C11)A perturbation with the external potential δv ( r , ω ) in-duces the responses δn v ( r , ω ) and δq α,v ( ω ). Making asubstitution of Eqs. (C4) and (C5) into the density anddisplacement field response due to an external potential δv ( r , ω ) yields after some algebra the following eigenvalueproblem (for a detailed derivation, we refer the reader toappendix S4)8 L (Ω q ) K (Ω q ) M (Ω q ) M (Ω q ) K ∗ (Ω q ) L (Ω q ) M ∗ (Ω q ) M ∗ (Ω q ) N N ∗ ω α N N ∗ ω α X (Ω q ) Y (Ω q ) A (Ω q ) B (Ω q ) = Ω q − − X (Ω q ) Y (Ω q ) A (Ω q ) B (Ω q ) . (C12)In this equation, X and Y are the contributions tothe full solution in the matter part of the equation,while A and B are the contributions to the solu-tion in the photon part of the equation. Further, Ω q refers to the many-body electron-photon excitation en-ergies. In comparison to the standard linear-responseformulation of TDDFT, new 2 × M -block accounts for the explicit electron-photon interac-tion, the N -block accounts for the dipole coupling of theelectronic system to the photon field and the ω α -blockare the frequencies of the photon field. The quantity L ai,jb (Ω q ) = δ ab δ ij ( (cid:15) a − (cid:15) i )+ K ai,jb (Ω q ) contains the dif-ference of two Kohn-Sham energies (cid:15) a and (cid:15) i , where i refers to occupied orbitals and the index a to unoccupiedorbitals. The coupling-matrix K is given by K ai,jb (Ω q ) = (cid:90) (cid:90) d r d y ϕ i ( r ) ϕ ∗ a ( r ) f nMxc ( r , y , Ω q ) ϕ b ( y ) ϕ ∗ j ( y ) . (C13)The quantity K ai,jb (Ω q ) differs from the electron-onlycase since f nMxc = f nHxc + f npxc . Treating the photonfield only externally reduces this matrix to the standardcoupling matrix in TDDFT linear response with f nMxc = f nHxc . The two new coupling functions appearing, M and N that couple the matter block are given explicitly as M α,ai (Ω q ) = (cid:90) d r ϕ i ( r ) ϕ ∗ a ( r ) f q α Mxc ( r , Ω q ) , (C14) N α,ia = 12 ω α (cid:90) d r ϕ ∗ i ( r ) ϕ a ( r ) g n α M ( r ) . (C15)We emphasize here, that the exact coupling matrix N α,ia has no frequency dependence since the exact kernel forEq. (C15) is equivalent to just the mean-field kernel ofthe photon modes as can be seen from Eq. (9). Giventhe exact kernels, the nonlinear pseudo-eigenvalue prob-lem in Eq. (C12) allows to compute the exact excita-tion energies of the coupled matter-photon system. Ofcourse, in practice, approximations have to be employedfor the matter-photon response kernels as is also requiredin the matter-only response formalism. Since the ex-plicitly known mean-field kernel g n α M ( r ) is already exact,only f q α Mxc ( r , Ω q ) and f nMxc ( r , y , Ω q ) are left to be approx-imated.The above matrix equation of Eq. (C12) can be cast intoa Hermitian eigenvalue form following the same transfor-mations as, e.g., in Ref. [88], where we assume real-valuedorbitals, i.e., K = K ∗ , M = M ∗ and N = N ∗ . Further,we drop the dependency on Ω q for brevity. Then we find the pseudo-eigenvalue equation, reminiscent to the equa-tions found for excitation energies in Hartree-Fock theoryand TDDFT [64]. The eigenvalue problem of Eq. (C12)is now written in a compact Hermitian form as (cid:18) U VW ω α (cid:19) (cid:18) E P (cid:19) = Ω q (cid:18) E P (cid:19) , (C16)where the matrices U , V and W are given by U = ( L − K ) / ( L + K )( L − K ) / , V = 2( L − K ) / M / N / ω / α and W = 2 ω / α N / M / ( L − K ) / . The matrices aregiven explicitly by U qq (cid:48) = δ qq (cid:48) ω q + 2 √ ω q ω q (cid:48) K qq (cid:48) (Ω q ) , (C17) V qα = 2 (cid:113) ω q M αq (Ω q ) N αq ω α , (C18) W αq = 2 (cid:113) ω α N αq M αq (Ω q ) ω q , (C19)where the off-diagonal matrices V qα and W αq are trans-pose of each other, i.e., V qα = W (cid:62) αq . The index q = ( a, i )describes transitions from the electronic occupied ( i ) tounoccupied states ( a ) and thus the difference of Kohn-Sham energies is given by ω q = (cid:15) a − (cid:15) i . With α wedenote the photon modes. The eigenvectors E and P can be used to compute oscillator strengths of the cou-pled matter-photon system (see appendix S5). In thedecoupling limit of light-matter interaction, Eq. (C16)reduces to the well-known Casida equation
Eq.(S4) [64].So far we did not solve anything but have just rewrit-ten the problem in terms of unknown Mxc kernels thatcorrect the uncoupled and non-interacting auxiliary re-sponse functions. To actually solve this problem we needto provide approximations to these unknown quantities.Here it becomes advantageous to have divided the fullMxc kernels in Hxc and pxc terms, such that we can usewell-established approximations from electronic TDDFTfor the Hxc and specifically developed approximationsfor the pxc terms (see Sec. III for more details). In thefollowing, we will employ the above introduced pRPAapproximation, which is a straightforward generalizationof the standard RPA of electronic-structure theory andyields the following kernels f n H ( r , r (cid:48) ) = e π(cid:15) | r − r (cid:48) | , f q α p ( r ) = − ω α λ α · e r ,f n α p ( r , r (cid:48) ) = (cid:88) α ( λ α · e r (cid:48) ) λ α · e r , g n α M ( r ) = − ω α λ α · e r . We note that at the pRPA level the matter-photon cou-pling mediated via g n α M is exact. The influence of the9photon-matter xc contributions f q α xc and f n xc will be high-lighted in the next section. By connecting to the eigen-states E and P , we can assign to each of the individualpoles of the response function, i.e. the excitation ener-gies, the amount of photonic and electronic contributionto that excitation by using σ e = N pairs (cid:88) i =1 | E ,i | , (C20) σ p = M (cid:88) α =1 | P ,α | , (C21)where N pairs corresponds to the number of occupied-unoccupied pairs of KS orbitals, in our case 30 × σ e and σ p is normalized to one, i.e. σ e + σ p = 1.In the pRPA approximation all the frequency depen-dence that we suppressed at times for brevity now gen-uinely vanishes (an adiabatic approximation) which al-lows us to express M and N of Eqs. (C14) and (C15)as M α,ai = − ω α (cid:90) drϕ i ( r ) ϕ ∗ a ( r ) λ α · e r , (C22) N α,ia = − (cid:90) drϕ ∗ i ( r ) ϕ a ( r ) λ α · e r . (C23)Next we want to connect to the standard matter-onlylinear-response framework [57]. In defining the oscilla-tor strength for the density-density response function, wemake use of the relationship between the polarizabilitytensor and susceptibility. The first-order dipole polariz-ability is given by δ R ( t ) = (cid:90) d r e r δn ( r , t ) , (C24)and in frequency space R ( ω ) = ←→ α ( ω ) E ( ω ). The dy-namic polarizability tensor can then be written as ←→ α µν ( ω ) = (cid:90) d r er µ δn ( r , ω ) δE ν ( ω ) , (C25)with µ, ν = (1 , ,
3) denoting all three spatial directions.Connecting to the QEDFT linear-response theory, wefind ←→ α µν ( ω ) = (cid:88) I r † µ S / Z I Z † I S / r ν ω − Ω I , (C26)where r Iµ = (cid:82) d r er µ (cid:80) i,a Φ ia ( r ) is the Kohn-Sham tran-sition dipole matrix element of the many-body transition I . Further, we have used S = ( L − K ) and the transi-tion density is defined as Φ ia ( r ) = ϕ ∗ i ( r ) ϕ a ( r ) in termsof Kohn-Sham orbitals. This then allows to obtain thefull photoabsorption cross section from the trace of thepolarizability tensor through σ ( ω ) = 4 πωc I m Tr ←→ α ( ω ) / . (C27) For the oscillator strength [57, 64], we find f I = 23 (cid:88) µ =1 (cid:12)(cid:12)(cid:12) Z † I S / r Iµ (cid:12)(cid:12)(cid:12) = 23 ω I (cid:88) µ =1 |(cid:104) Ψ | er µ | Ψ I (cid:105)| (C28)and also in the case of QEDFT, the oscillator strengthsatisfy the Thomas-Reiche-Kuhn sum rule (also knownas f -sum rule), i.e. (cid:80) I f I = N , where N is the totalnumber of electrons in the system. At this point, we alsowant to introduce the dipole strength function S ( ω ) [57]that is defined as S ( ω ) = (cid:88) I f I δ ( ω − Ω I ) (C29)and integrates according to the f -sum rule to the totalnumber of electrons. For the non-standard part of ourresponse theory, i.e., matter-photon and photon-photonperturbations, we use similar constructions to display theresults. Their derivations and definitions are given in ap-pendix S5. We will discuss their physical meaning in thenext section where we employ a simple yet illuminatingmodel system. This will not only allow us to explainmany of the so far abstract ideas in a straightforwardmanner, but we can also test the accuracy of the pRPA. Appendix D: Examples for the coupledmatter-photon response: Details on the Rabi Model
In this section, we give more details on the modelsystem that have been employed in Sec. III. The modelHamiltonian we consider is given by (in this section weswitch for simplicity to atomic units)ˆ H R ( t ) = ω σ z + ω c ˆ a † ˆ a + λ ˆ σ x ˆ q + j ( t )ˆ q + v ( t )ˆ σ x , (D1)where ω is the transition frequency between the groundstate | g (cid:105) and excited state | e (cid:105) and ˆ σ x as well as ˆ σ z arethe usual Pauli matrices. We only keep one photon modewith frequency ω c and use the usual photon creation andannihilation operators to represent the harmonic oscilla-tor of this mode. By further compressing the notation,we then describe the coupling between matter and lightby a coupling strength λ and the displacement coordi-nate ˆ q = √ ω c (cid:0) ˆ a + ˆ a † (cid:1) . Finally, we couple the mattersystem to a classical external perturbation v ( t ) and thephoton system to a classical external current j ( t ) (a pic-torial representation of the coupled system is given inFig. 3). We note for consistency with respect to otherworks [37, 44, 49] that in the above Rabi model we canperform a unitary transformation that allows us to ex-change ˆ σ x and ˆ σ z . Both forms of the extended Rabimodel are therefore equivalent. We further note thatwith respect to the full non-relativistic QED problem inthe long wavelength approximation of Eq. (1) the Rabimodel does not include the dipole self-energy term pro-portional to ( λ · R ) . This is because the analogous term0in this model is just a constant energy shift, i.e., it isproportional to σ x = ˆ [43]. For more levels this is nolonger the case [49] and this term has to be taken into ac-count, else the resulting eigenstates do not have a propercontinuum limit [43]. The responses that we want toconsider in the following are those observables that cou-ple to the external perturbations. In our case this is σ x ( t ) = (cid:104) Ψ( t ) | ˆ σ x | Ψ( t ) (cid:105) (in essence the atomic dipole) andthe displacement field q ( t ) = (cid:104) Ψ( t ) | ˆ q | Ψ( t ) (cid:105) .The response of these observables ( δσ x ( t ) , δq ( t )) to per-turbations by the external pair ( δv ( t ) , δj ( t )) can be writ-ten similarly as Eq.(C1) in the collective form (cid:18) δσ x ( t ) δq ( t ) (cid:19) = (cid:90) dt (cid:48) (cid:18) χ σ x σ x ( t, t (cid:48) ) χ σ x q ( t, t (cid:48) ) χ qσ x ( t, t (cid:48) ) χ qq ( t, t (cid:48) ) (cid:19) (cid:18) δv ( t (cid:48) ) δj ( t (cid:48) ) (cid:19) . (D2)Again we find besides the usual matter-matter response χ σ x σ x also matter-photon responses χ σ x q and χ qσ x , respec-tively, as well as a photon-photon response function χ qq .Next, in analogy to Sec. II, we reformulate the coupledmatter-photon problem in form of a Maxwell KS auxil-iary problem. Using by now well-established results ofQEDFT for the extended Rabi model systems [37, 44]we can introduce two effective fields v Mxc ([ σ x , q ]; t ) = v s ([ σ x ]; t ) − v ([ σ x , q ]; t ) , (D3) j M ([ σ x ]; t ) = j s ([ q ]; t ) − j ([ σ x , q ]; t ) , (D4)that force the auxiliary uncoupled, yet non-linearMaxwell KS system to generate the same dynamics ofthe internal pair ( σ x ( t ) , q ( t )) as the corresponding cou-pled reference system. For an uncoupled initial Maxwellstate | Ψ (cid:105) = | ψ (cid:105) ⊗ | ϕ (cid:105) that provides the same initialconditions for the internal pair as the physical initial state[37, 44], we then have to solve self-consistently i ∂∂t | ψ ( t ) (cid:105) = (cid:104) ω σ z + ( v ( t ) + v Mxc ([ σ x , q ]; t )) ˆ σ x (cid:105) | ψ ( t ) (cid:105) , (D5) (cid:18) ∂ ∂t + ω c (cid:19) q ( t ) = − j ( t ) − λσ x ( t ) . (D6)Since the photon subsystem is merely a shifted harmonicoscillator we get away with only solving the classical har- monic oscillator equation coupled to the dipole of thematter subsystem. We can then express the coupled re-sponse functions of Eq. (D2) in analogy to Eqs. (10)-(13)by the uncoupled auxiliary response functions χ σ x σ x ,s and χ qq,s as χ σ x σ x ( t, t (cid:48) ) = χ σ x σ x ,s ( t, t (cid:48) ) (D7)+ (cid:90) (cid:90) dτ dτ (cid:48) χ σ x σ x ,s ( t, τ ) f σ x Mxc ( τ, τ (cid:48) ) χ σ x σ x ( τ (cid:48) , t (cid:48) )+ (cid:90) (cid:90) dτ dτ (cid:48) χ σ x σ x ,s ( t, τ ) f qMxc ( τ, τ (cid:48) ) χ qσ x ( τ (cid:48) , t (cid:48) ) ,χ σ x q ( t, t (cid:48) ) = (cid:90) (cid:90) dτ dτ (cid:48) χ σ x σ x ,s ( t, τ ) f qMxc ( τ, τ (cid:48) ) χ qq ( τ (cid:48) , t (cid:48) )(D8)+ (cid:90) (cid:90) dτ dτ (cid:48) χ σ x σ x ,s ( t, τ ) f σ x Mxc ( τ, τ (cid:48) ) χ σ x q ( τ (cid:48) , t (cid:48) ) ,χ qσ x ( t, t (cid:48) ) = (cid:90) (cid:90) dτ dτ (cid:48) χ qq,s ( t, τ ) g σ x M ( τ, τ (cid:48) ) χ σ x σ x ( τ (cid:48) , t (cid:48) ) , (D9) χ qq ( t, t (cid:48) ) = χ qq,s ( t, t (cid:48) ) (D10)+ (cid:90) (cid:90) dτ dτ (cid:48) χ qq,s ( t, τ ) g σ x M ( τ, τ (cid:48) ) χ σ x q ( τ (cid:48) , t (cid:48) ) . The only real difference is that in the Rabi case we donot have a longitudinal interaction and therefore the Mxccontributions come solely from the matter-photon cou-pling, i.e., f σ x Mxc = f σ x pxc and f qMxc = f qpxc . This allows usto study exclusively the influence of these new terms andhow approximations of them perform.
1. Matter-photon correlation effect in Maxwell’sequations
Let us follow the previous general section II D andbriefly consider the influence of the matter-photon cou-pling on the Maxwell’s equations in this model system,i.e., Eq. (D6). The inhomogeneous Maxwell’s equationhere accounts for the back-reaction of the matter on thefield through the atomic dipole operator σ x ([ v, j ]; t ). Ifwe, for instance, perturb the two-level system directlyvia a δv ( t ), the response of the Maxwell’s equation ex-pressed in terms of the uncoupled problem with the helpof Eq. (D8) becomes (cid:0) ∂ t + ω c (cid:1) δq ( t ) = − λ (cid:90) dt (cid:48) χ σ x σ x ,s ( t, t (cid:48) ) δv ( t (cid:48) ) − λ (cid:90) (cid:90) (cid:90) dt (cid:48) dτ dτ (cid:48) χ σ x σ x ,s ( t, τ ) f σ x Mxc ( τ, τ (cid:48) ) χ σ x σ x ( τ (cid:48) , t (cid:48) ) δv ( t (cid:48) ) − λ (cid:90) (cid:90) dτ dτ (cid:48) χ σ x σ x ,s ( t, τ ) f qMxc ( τ, τ (cid:48) ) δq ( τ (cid:48) ) . (D11)Having no coupling, i.e., the Mxc terms are zero, merelyrecovers the usual inhomogeneous Maxwell’s equation fora classical external current. The matter system evolvesaccording to the perturbation and we can determine its induced Maxwell field without any back-reaction. Thesecond term describes the matter polarization due to theinduced field and leads to an effective self-interaction ofthe two-level system. If there would be more than one1particle this would induce an effective matter-matter in-teraction as well. The third term then accounts for thefield polarization and induces an effective self-interactionin the mode of the light field. That is, the coupling to matter leads to a photon-photon interaction. This canbe made more explicit by separating the mean-field con-tribution v M ( t ) = λq ( t ) and rewriting the above equationas (cid:0) ∂ t + ω c (cid:1) δq ( t ) = − λ (cid:90) dt (cid:48) χ σ x σ x ,s ( t, t (cid:48) ) δv ( t (cid:48) ) − λ (cid:90) (cid:90) (cid:90) dt (cid:48) dτ dτ (cid:48) χ σ x σ x ,s ( t, τ ) f σ x xc ( τ, τ (cid:48) ) χ σ x σ x ( τ (cid:48) , t (cid:48) ) δv ( t (cid:48) ) − λ (cid:90) dτ χ σ x σ x ,s ( t, τ ) δq ( τ ) − λ (cid:90) (cid:90) dτ dτ (cid:48) χ σ x σ x ,s ( t, τ ) f qxc ( τ, τ (cid:48) ) δq ( τ (cid:48) ) . (D12)The third term on the right-hand side is then the pRPAform of photon-photon response. Similar terms also ap- pear for a perturbation induced by an external current δj ( t ) which can be rewritten with the help of Eq. (D9)and the mean-field made explicit as (cid:0) ∂ t + ω c (cid:1) δq ( t ) = − δj ( t ) − λ (cid:90) dτ χ σ x σ x ,s ( t, τ ) δq ( τ ) − λ (cid:90) (cid:90) dτ dτ (cid:48) χ σ x σ x ,s ( t, τ ) f qxc ( τ, τ (cid:48) ) δq ( τ (cid:48) ) − λ (cid:90) (cid:90) (cid:90) dt (cid:48) dτ dτ (cid:48) χ σ x σ x ,s ( t, τ ) f σ x xc ( τ, τ (cid:48) ) χ σ x q ( τ (cid:48) , t (cid:48) ) δj ( t (cid:48) ) . (D13)Here we used that f q M ( τ, τ (cid:48) ) = λδ ( τ − τ (cid:48) ). As is mostobvious in the pRPA limit, both types of perturbationslead to the same resonance conditions, i.e., peaks in theresponses. They are connected to the combined eigen-states of the matter-photon system. However, the de-tailed response can differ strongly. That these resonanceconditions that we get from the pRPA are indeed con-nected to the coupled eigenstates we will show next.
2. Application of the pseudo-eigenvalue problem
As a preparatory step we first rewrite the linear-response problem of the extended Rabi model in termsof the previously introduced pseudo-eigenvalue problemof Eq. (C16). In the two-level one-mode case we considerhere, this reduces to (cid:18) U (Ω q ) V (Ω q ) W (Ω q ) ω c (cid:19) (cid:18) E P (cid:19) = Ω q (cid:18) E P (cid:19) . (D14)Where the matrices in the model system reduceto functions of Ω q as U = ω + 2 ω K (Ω q ), V (Ω q ) = 2 ω / M (Ω q ) / N / ω / c , and W (Ω q ) =2 ω / c N / M (Ω q ) / ω / . The coupling functions aregiven explicitly using Eqs.(C13)-(C15) as K (Ω q ) = f σ x M + f σ x xc (Ω q ) ,M (Ω q ) = f qM + f qxc (Ω q ) ,N = 12 ω c g σ x M , where the Kohn-Sham states is the dipole matrix element ϕ a ϕ ∗ i = (cid:104) g | ˆ σ x | e (cid:105) = 1. The Mxc kernels can be definedusing the inverse of the auxiliary and interacting responsefunctions (see also Eqs. (S7)-(S8) in the appendix) andare given in frequency space by f σ x Mxc ( ω ) = (cid:0) χ σ x σ x ,s ( ω ) (cid:1) − − (cid:0) χ σ x σ x ( ω ) (cid:1) − , (D15) f qMxc ( ω ) = − ( χ qn ( ω )) − . (D16)Here (cid:0) χ σ x σ x ,s ( ω ) (cid:1) − , (cid:0) χ σ x σ x ( ω ) (cid:1) − and ( χ qn ( ω )) − are theinverses of the uncoupled response function of the elec-tronic subsystem, the fully coupled response function ofthe electronic dipole and of the displacement field of theRabi model, respectively. With these quantities we thendetermine spectroscopic observables such as the photoab-sorption cross section. To determine this cross sectionwe first note that the linear polarizability α ( ω ) inducedby the external potential v ( ω ) is related to the “dipole-dipole” response function as α ( ω ) = χ σ x σ x ( ω ). UsingEq. (C27), we can determine the photoabsorption crosssection of the Rabi model (see Fig. 4 (a) displayed indotted-red for the numerically exact case). σ ( ω ) = 4 πωc I m χ σ x σ x ( ω ) . (D17)Here, the mean of the polarizability was not consideredsince the Rabi model is a one-dimensional system. Anal-ogously, we define a linear “field polarizability” β ( ω ) dueto polarizing the photon mode by an external current. Inthe same way, we relate the field polarizability to the re-sponse function of the photon mode as β ( ω ) = χ qq ( ω ) and2then determine a photonic spectrum from (see Fig. 4 (b)displayed in dotted-red for the numerically exact case).˜ σ ( ω ) ≡ πωc I m χ qq ( ω ) . (D18)Finally, we consider mixed spectroscopic observableswhere we perturb one subsystem and then consider theresponse in the other. We analogously employ χ qσ x ( ω )and χ σ x q ( ω ) in Eqs.(D17) and (D18), respectively, to de-termine a “mixed polarizability”. If we plot this mixedspectrum (see Fig. 4 (c) displayed in dotted-red for thenumerically exact case), we find that we have positiveand negative peaks. Indeed, this highlights that exci-tations due to external perturbations can be exchangedbetween subsystems, i.e., energy absorbed in the elec-tronic subsystem can excite the photonic subsystem andvice versa. Next, we want to employ the pRPA ap-proximation to the extended Rabi model and try tosolve it analytically. The pRPA is equivalent to usingthe mean-field approximation in the coupled equations,i.e., approximating the electron-photon coupling termas ˆ σ x ˆ q ≈ (cid:104) ˆ σ x (cid:105) ˆ q + (cid:104) ˆ q (cid:105) ˆ σ x . This corresponds then to acoupled Schr¨odinger-Maxwell treatment of the coupledmatter-photon problem [13]. In the Maxwell KS equa-tions this leads to approximating the full v Mxc by themean-field potential v M = v p = λq . The mean-field cur-rent is known explicitly as j M = λσ x . In the case of thepseudo-eigenvalue problem this amounts to approximat-ing K = f σ x M = 0, M = f q M = λ and N = ω c g σ x M = λ ω c .Consequently we have U = ω , V = W = 2 λ (cid:114) ω , ω α = ω c . The resulting nonlinear eigenvalue equation yields theexcitation frequenciesΩ ( − ) = 12 (cid:0) ω + ω c (cid:1) − (cid:113) ( ω − ω c ) + 8 λ ω , (D19)Ω (+) = 12 (cid:0) ω + ω c (cid:1) + 12 (cid:113) ( ω − ω c ) + 8 λ ω , (D20)and the corresponding normalized eigenvectors can begiven in closed form as E = (cid:18) − sin θ cos θ (cid:19) , and P = (cid:18) cos θ sin θ (cid:19) (D21)The resulting pRPA-approximated spectra are dis-played in Fig. 4 in dashed-blue. We will discuss the re-sults in a little more detail at the end of this section.Before we consider a slightly more advanced approxima-tion based on the rotating-wave approximation (RWA).If we slightly simplify the full Rabi problem by approxi-mating the full coupling as ˆ σ x ˆ q ≈ √ ω c (cid:0) ˆ σ + ˆ a + ˆ σ − ˆ a † (cid:1) weend up with the Jaynes-Cumming Hamiltonian [89] givenas ˆ H JC ( t ) = ω σ z + ω c ˆ a † ˆ a + λ √ ω c (cid:0) ˆ σ + ˆ a + ˆ σ − ˆ a † (cid:1) + j ( t )ˆ q + v ( t )ˆ σ x . (D22) Here we used ˆ σ ± = (ˆ σ x ± i ˆ σ y ) /
2. The above approxima-tion is called the RWA because we ignore quickly oscil-lating terms and thus assume that the excitation of thematter subsystem can only destroy and the de-excitationonly create a photon. This approximation is justified(with respect to the full wave function) if we are in theweak coupling regime, i.e., λ (cid:28) ω c , and near to res-onance, i.e., δ = ω − ω c ≈
0. The ground-state of theJaynes-Cummings model is the uncoupled tensor productof the matter ground-state and the photon ground-statewith ground-state energy of E = − ω /
2. The excitedstates of the Jaynes-Cummings Hamiltonian are knownanalytically and are given by (we only show the lowest ly-ing excited states where a single photon is excited and forwhich the matrix elements are non-zero with the ground-state) |− , (cid:105) = − sin θ | g (cid:105)| (cid:105) + cos θ | e (cid:105)| (cid:105) , (D23) | + , (cid:105) = cos θ | g (cid:105)| (cid:105) + sin θ | e (cid:105)| (cid:105) . (D24)With these eigenstates we find the transition frequen-cies that correspond to the linear response from theground state (due to the approximations involved onlyone photon absorbed or emitted) to beΩ − (0) = 12 ( ω c + ω − Ω ) , (D25)Ω + (0) = 12 ( ω c + ω + Ω ) , (D26)where Ω = √ δ + 4 λ (cid:48) where λ (cid:48) = λ √ ω c . So we al-ready know where the RWA will generate the poles ofthe response function. Since we know analytically theeigenfunctions in the RWA, we can construct the RWAresponse functions analytically. Using the definitions ofthe Mxc kernels of Eq. (D15) we can then analyticallyconstruct the RWA Mxc kernels. These kernels are fre-quency dependent, therefore the resulting Mxc approx-imation is non-adiabatic [57]. Substituting them into M (Ω q ) we recover the known poles Ω q = Ω ± (0) fromEq. (D14). Further, we can then construct the differ-ent spectra associated with the RWA. We show them inFig. 4 in full-orange. Appendix E: Numerical details
We start by discussing the general setup before con-sidering the specialized situations discussed above. Wehave implemented the linear-response pseudo-eigenvalueequation of Eq. (C16) into the real-space code OCTO-PUS [74, 75]. The absorption spectrum of the ben-zene molecule has been very successfully studied withTDDFT calculations [74, 77]. Small organic moleculesand benzene in particular are rewarding systems to bestudied with TDDFT, since the adiabatic approxima-tion in concert with the local-density approximation(LDA) [90, 91] capture the occurring Π-Π ∗ transition ex-ceptionally well [77]. This transition is a characteristic of3carbon conjugate compounds [74] and occurs around 7 eVin the case of a benzene molecule. To calculate the elec-tronic structure of the benzene molecule, we follow closelythe setup of Ref. [74]. Thus, we use a cylindric real spacegrid of 8 ˚A length with the radius of 6 ˚A in the x - y plane,and a spacing of ∆ x = 0 . A . For the benzene nuclearstructure, we use the CC bond length of 1 .
396 ˚A, andCH bond length of 1 .
083 ˚A. We explicitly describe the 30valence electrons, while the core atoms are considered im-plicitly by LDA Troullier-Martins pseudopotentials [92].In the excited state manifold, we include 500 unoccupiedstates in the pseudo-eigenvalue calculation. This numberamounts to 30 ×
500 = 7500 pairs of occupied-unoccupiedstates. Further, to describe the electron-electron interac-tion in the response functions, we apply the adiabaticLDA (ALDA) kernel, i.e. f nMxc → f n Hxc,ALDA + f np .Solving the linear-response pseudo-eigenvalue problem ofEq. (C16) provides us with the transition amplitudes, aswell as the excitation energies of the correlated electron-photon system. These quantities can be used to calcu-late photoabsorption spectra by using e.g. Eq. (C29). Instandard calculations, to obtain such spectra and mimicthe finite lifetime of the excited state usually a peak-broadening is applied. In our case, where necessary, weapply a Lorentzian broadening, i.e. the standard imple-mentation of the OCTOPUS code, that is of the following form Γ( ω, ω I ) = 1 π ∆( ω − ω I ) + ∆ , (E1)where ω I is the excitation frequency, and ∆ the broad-ening parameter. The actual dipole strength function asdefined in Eq. (C29) is then obtained by S ( ω ) = (cid:88) I f I Γ( ω, ω I ) , (E2)where f I denotes the oscillator strength as defined inEq. (C28). We obtain the spectra for systems notimmersed in the photon bath in this paper, where thepeaks have been broadened by applying the broadeningas defined in Eq. (E1) with ∆ = 0 . [1] T. W. Ebbesen, Accounts of Chemical Research , 2403(2016).[2] M. Sukharev and A. Nitzan, Journal of Physics: Con-densed Matter , 443003 (2017).[3] J. George, A. Shalabney, J. A. Hutchison, C. Genet, andT. W. Ebbesen, J. Phys. Chem. Lett. , 1027 (2015),http://dx.doi.org/10.1021/acs.jpclett.5b00204.[4] H. Hiura, A. Shalabney, and J. George, (2018),10.26434/chemrxiv.7234721.v3.[5] A. Thomas, L. Lethuillier-Karl, K. Nagara-jan, R. M. A. Vergauwe, J. George, T. Chervy,A. Shalabney, E. Devaux, C. Genet, J. Moran,and T. W. Ebbesen, Science , 615 (2019),http://science.sciencemag.org/content/363/6427/615.full.pdf.[6] C. Riek, D. V. Seletskiy, A. S. Moskalenko, J. F. Schmidt,P. Krauspe, S. Eckart, S. Eggert, G. Burkard, andA. Leitenstorfer, Science , 420 (2015).[7] D. Coles, L. C. Flatten, T. Sydney, E. Hounslow, S. K.Saikin, A. Aspuru-Guzik, V. Vedral, J. K.-H. Tang, R. A.Taylor, J. M. Smith, and D. G. Lidzey, Small ,1701777 (2017), 1701777.[8] R. Chikkaraddy, B. de Nijs, F. Benz, S. J. Barrow, O. A.Scherman, E. Rosta, A. Demetriadou, P. Fox, O. Hess,and J. J. Baumberg, Nature , 127 (2016).[9] F. Benz, M. K. Schmidt, A. Dreismann, R. Chikkaraddy,Y. Zhang, A. Demetriadou, C. Carnegie, H. Ohadi,B. de Nijs, R. Esteban, J. Aizpurua, and J. J. Baumberg,Science , 726 (2016).[10] M.-E. Kleemann, R. Chikkaraddy, E. Alexeev, D. Kos,C. Carnegie, W. Deacon, A. Casalis De Pury, C. Grosse,B. De Nijs, J. Mertens, A. Tartakovskii, and J. Baum- λ α [eV / /nm] E I [eV] (cid:104) x I (cid:105) [A] f I [a.u.]0 6.88 0.952 0.5462.77 6.69 0.721 0.3042.77 7.03 0.626 0.2415.55 6.49 0.791 0.3555.55 7.18 0.550 0.1908.32 6.28 0.848 0.3958.32 7.30 0.482 0.14911.09 6.06 0.896 0.42611.09 7.41 0.420 0.114TABLE I. Rabi splitting of the Π − Π ∗ transition: electron-photon interaction strength λ α = | λ α | , excitation energy E I ,transition dipole moment x I and the oscillator strength f I .berg, (2017), 10.17863/cam.13084.[11] A. Bisht, J. Cuadra, M. Wers¨all, A. Canales, T. J. An-tosiewicz, and T. Shegai, Nano Letters , 189 (2019),https://doi.org/10.1021/acs.nanolett.8b03639.[12] M. Mirhosseini, E. Kim, X. Zhang, A. Sipahigil, P. B.Dieterle, A. J. Keller, A. Asenjo-Garcia, D. E. Chang,and O. Painter, arXiv e-prints , arXiv:1809.09752 (2018),arXiv:1809.09752 [quant-ph].[13] M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Ap-pel, and A. Rubio, Nature Reviews Chemistry , 0118(2018). [14] J. Flick, N. Rivera, and P. Narang, Nanophotonics ,1479 (2018).[15] J. Galego, F. J. Garcia-Vidal, and J. Feist, Phys. Rev.X , 041022 (2015).[16] J. Galego, F. J. Garcia-Vidal, and J. Feist, Nature Com-munications , 13841 (2016).[17] M. Kowalewski, K. Bennett, and S. Mukamel,J. Phys. Chem. Lett. , 2050 (2016),http://dx.doi.org/10.1021/acs.jpclett.6b00864.[18] M. Cirio, S. De Liberato, N. Lambert, and F. Nori, Phys.Rev. Lett. , 113601 (2016).[19] F. Herrera and F. C. Spano, Phys. Rev. Lett. , 238301(2016).[20] J. Galego, F. J. Garcia-Vidal, and J. Feist, Phys. Rev.Lett. , 136001 (2017).[21] P. Roelli, C. Galland, N. Piro, and T. J. Kippenberg,Nature Nanotechnology , 164 (2015).[22] D. Shin, H. H¨ubener, U. D. Giovannini, H. Jin, A. Ru-bio, and N. Park, Nature Communications (2018),10.1038/s41467-018-02918-5.[23] G. Mazza and A. Georges, Phys. Rev. Lett. , 017401(2019).[24] M. A. Sentef, M. Ruggenthaler, and A. Rubio, ScienceAdvances (2018), 10.1126/sciadv.aau6969.[25] D. Braak, Phys. Rev. Lett. , 100401 (2011).[26] Q. Xie, H. Zhong, M. T. Batchelor, and C. Lee, Journalof Physics A: Mathematical and Theoretical , 113001(2017).[27] B. M. Garraway, Philosophical Transactions of the RoyalSociety of London A: Mathematical, Physical and Engi-neering Sciences , 1137 (2011).[28] D. De Bernardis, P. Pilar, T. Jaako, S. De Liberato, andP. Rabl, Phys. Rev. A , 053819 (2018).[29] C. Sch¨afer, M. Ruggenthaler, and A. Rubio, Phys. Rev.A , 043801 (2018).[30] C. F. Gilbert Grynberg, Alain Aspect, Introduction toQuantum Optics From Semi-classical Approach to Quan-tized Light (Cambridge University Press, New York,2010).[31] V. Weisskopf and E. Wigner, Zeitschrift f¨ur Physik ,54 (1930).[32] V. Buˇzek, G. Drobn´y, M. G. Kim, M. Havukainen, andP. L. Knight, Phys. Rev. A , 582 (1999).[33] R. Lettow, V. Ahtee, R. Pfab, A. Renn, E. Ikonen,S. G¨otzinger, and V. Sandoghdar, Opt. Express ,15842 (2007).[34] D. Wang, H. Kelkar, D. Martin-Cano, T. Utikal,S. G¨otzinger, and V. Sandoghdar, Phys. Rev. X ,021014 (2017).[35] M. Ruggenthaler, Physics (2017),10.1103/physics.10.105.[36] I. V. Tokatly, Phys. Rev. Lett. , 233001 (2013).[37] M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel, I. V.Tokatly, and A. Rubio, Phys. Rev. A , 012508 (2014).[38] M. Ruggenthaler, ArXiv e-prints (2015),arXiv:1509.01417 [quant-ph].[39] J. Flick, C. Sch¨afer, M. Ruggenthaler, H. Ap-pel, and A. Rubio, ACS Photonics , 992 (2018),https://doi.org/10.1021/acsphotonics.7b01279.[40] L. H. Ryder, Quantum field theory (Cambridge universitypress, 1996).[41] D. Craig and T. Thirunamachandran,
Molecular Quan-tum Electrodynamics: An Introduction to Radiation-molecule Interactions , Dover Books on Chemistry Series (Dover Publications, 1998).[42] H. Spohn,
Dynamics of charged particles and their radi-ation field (Cambridge university press, 2004).[43] V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Ru-bio, Journal of Physics B: Atomic, Molecular and OpticalPhysics , 034005 (2018).[44] C. Pellegrini, J. Flick, I. V. Tokatly, H. Appel, andA. Rubio, Phys. Rev. Lett. , 093001 (2015).[45] In principle, QEDFT can be formulated for each levelof theory of QED as presented in Ref. [37]. As a conse-quence, our formalism can be extended to more generalformulations, including full minimal coupling, beyond thedipole approximation.[46] A. Szabo and N. Ostlund, Modern Quantum Chemistry:Introduction to Advanced Electronic Structure Theory ,Dover Books on Chemistry (Dover Publications, 1989).[47] All results presented in this paper are independent of r .[48] Throughout this paper, we use the implicit definition e = −| e | .[49] T. Dimitrov, J. Flick, M. Ruggenthaler, and A. Ru-bio, New Journal of Physics (2017), 10.1088/1367-2630/aa8f09.[50] In the following, we suppress the spin component of thewave function and focus exclusively on the spatial andmode dependence, i.e., Ψ( r , ..., r N , q , ..., q M ; t ).[51] N. Hoffmann, Response formalism in density-functionaltheory for quantum electrodynamics , Master’s thesis,Technical University Berlin (2016).[52] A. L. Fetter and J. D. Walecka,
Quantum theory of many-particle systems (Courier Corporation, 2003).[53] G. Stefanucci and R. van Leeuwen,
NonequilibriumMany-Body Theory of Quantum Systems: A Modern In-troduction (Cambridge University Press, 2013).[54] M. Bonitz,
Quantum kinetic theory (Springer, 1998).[55] R. M. Dreizler and E. K. Gross,
Density functional the-ory: an approach to the quantum many-body problem (Springer Science & Business Media, 2012).[56] E. Engel and R. Dreizler,
Density Functional Theory: AnAdvanced Course , Theoretical and Mathematical Physics(Springer, 2011).[57] C. A. Ullrich,
Time-dependent density-functional theory:concepts and applications (OUP Oxford, 2011).[58] J. Flick and P. Narang, Phys. Rev. Lett. , 113002(2018).[59] M. Ruggenthaler and R. van Leeuwen, EPL (EurophysicsLetters) , 13001 (2011).[60] J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio,Proc. Natl. Acad. Sci. U. S. A. , 203202 (2015).[62] N. T. Maitra, K. Burke, and C. Woodward, Phys. Rev.Lett. , 023002 (2002).[63] M. Petersilka, U. J. Gossmann, and E. K. U. Gross,Phys. Rev. Lett. , 1212 (1996).[64] M. Casida, in Recent Developments and Applications ofModern Density Functional Theory , edited by J. Semi-nario (Elsevier, Amsterdam, 1996).[65] H. Ehrenreich, in
The Optical Properties of Solids (1966)p. 106.[66] W. L. Moch´an and R. G. Barrera, Physical Review B ,4984 (1985).[67] J. J. Maki, M. S. Malcuit, J. Sipe, and R. W. Boyd,Physical review letters , 972 (1991). [68] E. Luppi, H. H¨ubener, and V. V´eniard, Phys. Rev. B , 235201 (2010).[69] I. I. Rabi, Phys. Rev. , 324 (1936).[70] I. I. Rabi, Phys. Rev. , 652 (1937).[71] A. Shalabney, J. George, J. Hutchison, G. Pupillo,C. Genet, and T. W. Ebbesen, Nat. Commun. , 5981(2015).[72] J. George, T. Chervy, A. Shalabney, E. Devaux, H. Hiura,C. Genet, and T. W. Ebbesen, Phys. Rev. Lett. ,153601 (2016).[73] C. Ott, A. Kaldun, P. Raith, K. Meyer,M. Laux, J. Evers, C. H. Keitel, C. H.Greene, and T. Pfeifer, Science , 716 (2013),http://science.sciencemag.org/content/340/6133/716.full.pdf.[74] M. A. Marques, A. Castro, G. F. Bertsch, and A. Rubio,Computer Physics Communications , 60 (2003).[75] X. Andrade, D. Strubbe, U. D. Giovannini, A. H. Larsen,M. J. T. Oliveira, J. Alberdi-Rodriguez, A. Varas,I. Theophilou, N. Helbig, M. J. Verstraete, L. Stella,F. Nogueira, A. Aspuru-Guzik, A. Castro, M. A. L. Mar-ques, and A. Rubio, Phys. Chem. Chem. Phys. , 31371(2015).[76] The routines used to perform all calculations in this workwill be made publicly available. They can be easily trans-ported to any other first principles code that has the mat-ter linear-response equations implemented to make themready to describe the complete QED response, i.e. jointmatter-photon response, as described in this work.[77] K. Yabana and G. F. Bertsch, International Journal ofQuantum Chemistry , 55 (1999).[78] The spectrum in Ref. [74] has been obtained using anexplicit time-propagation with finite time. In the limit ofzero broadening and including all unoccupied states, wewould find identical spectra with very long propagatedspectra.[79] J. Flick, M. Ruggenthaler, H. Appel, andA. Rubio, Proceedings of the NationalAcademy of Sciences , 1 (1976).[82] Some textbooks [41] define the connection of the electricfield to the vector potential without the prefactor c . Weuse the current notation to be consistent with relativisticliterature and Ref. [37].[83] A. Di Piazza, C. M¨uller, K. Hatsagortsyan, and C. Kei-tel, Reviews of Modern Physics , 1177 (2012).[84] O. Firstenberg, T. Peyronel, Q.-Y. Liang, A. V. Gor-shkov, M. D. Lukin, and V. Vuleti´c, Nature , 71(2013).[85] K. E. Dorfman, F. Schlawin, and S. Mukamel, Rev. Mod.Phys. , 045008 (2016).[86] P. Blanchard and E. Br¨uning, Mathematical Methods inPhysics: Distributions, Hilbert Space Operators, Varia-tional Methods, and Applications in Quantum Physics ,Vol. 69 (Birkh¨auser, 2015).[87] G. Teschl,
Mathematical methods in quantum mechanics ,Vol. 157 (American Mathematical Soc., 2014).[88] R. Bauernschmitt and R. Ahlrichs, Chemical Physics Let-ters , 454 (1996).[89] E. Jaynes and F. Cummings, Proceedings of the IEEE , 89 (1963). [90] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965).[91] J. P. Perdew and A. Zunger, Phys. Rev. B , 5048(1981).[92] N. Troullier and J. L. Martins, Phys. Rev. B , 1993(1991).[93] R. Kubo, J. Phys. Soc. Jpn. 12 (1957).[94] M. M. Marques, C. Ullrich, F. Nogueira, A. Rubio,K. Burke, and E. Gross, Time-Dependent Density Func-tional Theory (Springer, Berlin, 2006).[95] E. K. U. Gross, J. F. Dobson, and M. Petersilka,
DensityFunctional Theory , edited by R. F. Nalewajski (Springer,Berlin, 1996).[96] C. Jamorski, M. E. Casida, and D. R. Salahub,The Journal of Chemical Physics , 5134 (1996),http://dx.doi.org/10.1063/1.471140.[97] M. E. Casida and D. R. Salahub, The Jour-nal of Chemical Physics , 8918 (2000),https://doi.org/10.1063/1.1319649.[98] M. Casida and M. Huix-Rotllant, Annual Review ofPhysical Chemistry , 287 (2012), pMID: 22242728,https://doi.org/10.1146/annurev-physchem-032511-143803.[99] Z. Yang and C. A. Ullrich, Phys. Rev. B , 195204(2013). Supplemental Information:Light-Matter Response in Non-Relativistic Quantum Electrodynamics: QuantumModifications of Maxwell’s Equations
Appendix S1: Current state of the art for spectroscopic: semi-classical description
To highlight the many differences of the presented framework to the standard linear-response approach we give herea brief recapitulation of the standard (matter-only) theory. The current theoretical description of linear spectroscopictechniques is built on the semi-classical approximation [30]. Herein, the many-particle electronic system is treatedquantum mechanically while the nuclei are subject to the Born-Oppenheimer approximation and the electromagneticfield appears as an external perturbation. As an external perturbation, the electromagnetic field probes the quantumsystem, but is not a dynamical variable of the complete system. To arrive at the semi-classical description startingfrom the full non-relativistic description of the Pauli-Fierz Hamiltonian [38], several approximations are used tosimplify the problem. In the following, we list these approximations explicitly • The mean-field approximation renders the Pauli-Fierz Hamiltonian as a problem of two coupled equations,i.e. the time-dependent Pauli equation and the inhomogeneous Maxwell’s equations, and is also know as theMaxwell-Pauli equation [13]. • The decoupling of these Maxwell-Pauli equations leads to the inhomogeneous Maxwell’s equation becomingindependent of the electronic system and all field effects are treated as a classical external field that perturbsthe many-electron system. • The dipole approximation, which ensures the uniformity of the external (decoupled) field over the extend of theelectronic system.Based on these approximations the Pauli-Fierz Hamiltonian [13] reduces to the time-dependent semi-classical Hamil-tonian for many-particle systems given asˆ H e ( t ) = N (cid:88) i =1 (cid:18) − (cid:126) m e ∇ i + v ( r i , t ) (cid:19) + e π(cid:15) N (cid:88) i>j | r i − r j | , (S1)including the kinetic energy, time-dependent external potential and the longitudinal Coulomb interaction. The time-dependent external potential has two parts v ( r , t ) = v ( r ) + δv ( r , t ). Here, v ( r ) describes the attractive part ofthe external potential due to the nuclei and δv ( r , t ) = e r · E ⊥ ( t ) with E ⊥ ( t ) being a classical external (transversal)probe field in dipole approximation that couples to the electronic subsystem. In this decoupling limit of light andmatter, the many-particle wavefunction is labeled only by the particle coordinate and spin as Ψ( r σ , ..., r N σ N ). In thedipole approximation we can investigate dipole-related spectroscopic observables such as polarizabiltiy, absorption andemission spectra, etc from linear to all orders in the external perturbation. Consider the particular case of a responseof an electronic system to an external weak probe field. In the dipole limit a key observable in the study of electronicand optical excitations in large many-particle systems is the electron density. Formulated within linear-response, thedensity response to an external perturbation is given as [93]: δn ( r t ) = − i (cid:126) (cid:90) tt dt (cid:48) (cid:90) d r (cid:48) (cid:104) Ψ | [ˆ n I ( r t ) , ˆ n I ( r (cid:48) t (cid:48) )] | Ψ (cid:105) = (cid:90) tt dt (cid:48) (cid:90) d r (cid:48) ˜ χ nn ( r t, r (cid:48) t (cid:48) ) δv ( r (cid:48) t (cid:48) ) . (S2)Here, ˜ χ nn ( r t, r (cid:48) t (cid:48) ) is the density-density function with respect to the ground-state Ψ ( r σ , ..., r N σ N ). Practical calcu-lations for the response of a many-electron system is a considerable challenge due to the large degrees of freedom. Inpractice, time-dependent density functional theory (TDDFT) [94, 95] is one of the most frequently applied theoriesto approach this problem. Knowing the electron density in TDDFT we can in principle calculate all observables ofinterest. Formulated within TDDFT linear-response, the density-density response function of the interacting systemcan be expressed in terms of non-interacting the density-density response function and an exchange-correlation (xc)kernel that has a form of a Dyson-type equation [63]:˜ χ nn ( r t, r (cid:48) t (cid:48) ) = χ nn, s ( r t, r (cid:48) t (cid:48) ) + (cid:90) (cid:90) d x d τ (cid:90) (cid:90) d τ (cid:48) d y χ nn, s ( r t, x τ ) f Hxc ( x τ, y τ (cid:48) ) ˜ χ nn ( y τ (cid:48) , r (cid:48) t (cid:48) ) , (S3)where χ nn, s and f Hxc = (cid:0) χ nn, s (cid:1) − − ( ˜ χ nn ) − . One of the most widely employed approaches to TDDFT linear-responseis the Casida formalism which can be written in a compact matrix form. The Casida equation obtains the exactexcitation energies Ω q of the many-particle system and requires all occupied and unoccupied Kohn-Sham orbitals andenergies including the continuum of states. In practice, the Casida equation is often cast into the following form U E = Ω q E . (S4)The explicit form of the matrix elements is given as (with q = ( i, a )) U qq (cid:48) = δ qq (cid:48) ω q + 2 √ ω q ω q (cid:48) K qq (cid:48) (Ω q ) , (S5) K ai,jb (Ω q ) = (cid:90) (cid:90) d r d r (cid:48) ϕ i ( r ) ϕ ∗ a ( r ) f Hxc ( r , r (cid:48) , Ω q ) ϕ b ( r (cid:48) ) ϕ ∗ j ( r (cid:48) ) . The Casida formalism is well established and has been applied to a variety of systems, see e.g. Refs. [64, 96–99] andreferences therein.The many obvious shortcomings of the approximations that lead to the standard Schr¨odinger equation (S1) arewell-known and discussed to some extend in the main part of the paper (for more details see, e.g., Ref. [13]). Wepoint out that all of the above ubiquitous fundamental equations are modified and the results based on the introducedgeneralized equations can differ strongly, as discussed in Sec. IV of the main article.
Appendix S2: Linear-response in non-relativistic QED
To help the reader with the unfamiliar generalized linear-response framework for coupled light-matter systems, wehere derive the linear-response equations and the ensuing response functions presented in Sec. II. In the non-relativisticsetting of QED, the static and dynamical behavior of the coupled electron-photon systems is given byˆ H ( t ) = ˆ H + ˆ H ext ( t ) . (S1)Where we define the time-independent electron-photon Hamiltonian asˆ H = ˆ T + ˆ W ee + 12 M (cid:88) α =1 (cid:34) ˆ p α + ω α (cid:18) ˆ q α − λ α ω α · R (cid:19) (cid:35) + N (cid:88) i =1 v ( r i ) + M (cid:88) α =1 j α, ω α ˆ q α , (S2)where the kinetic energy operator is ˆ T = − (cid:126) m e (cid:80) Ni =1 ∇ i , the Coulomb potential is ˆ W ee = e π(cid:15) (cid:80) Ni In this section, we present linear-response in QEDFT by employing the maps between interacting and non-interacting system, we express the interacting response functions in terms of two non-interacting response functionsand exchange correlation kernels. The responses due to ( v ( r t ) , j α ( t )) are evaluated at the ground-state ( v ( r ) , j α, )and will not be written explicitly.The non-interacting subsystems moving in an effective potential and current ( v s ( r t ) , j sα ( t )) can be written as atime-dependent problem of the Schr¨odinger i (cid:126) ∂∂t Φ( t ) = ˆ H KS ( t )Φ( t ) . (S1)Here, Φ( t ) is the wave function of the auxiliary non-interacting system and the non-interacting effective Hamiltonianˆ H KS ( t ) = ˆ H (0)KS + ˆ H ( ext )KS ( t ) that is meant to reproduce the exact density and displacement field, is given explicitly asˆ H (0)KS = ˆ T + ˆ H pt + (cid:16) v ( r ) + v (0) Mxc ([ n, q α ]; r ) (cid:17) + (cid:88) α ω α (cid:16) j α, + j (0) α,Mxc [ n, q α ] (cid:17) ˆ q α , and ˆ H ( ext )KS ( t ) = ( v ( r t ) + v Mxc ([ n, q α ]; r t )) + (cid:88) α ω α ( j α ( t ) + j α,Mxc ([ n, q α ]; t )) ˆ q α . Here ˆ H pt = (cid:80) Mα =1 (cid:2) ˆ p α + ω α ˆ q α (cid:3) is the oscillator for the photon mode and the mean-field xc potential and currentare defined as v Mxc ([ n, q α ]; r t ) := v s ([ n ]; r t ) − v ([ n, q α ]; r t ) , (S2) j α,Mxc ([ n, q α ]; t ) := j sα ([ q α ]; t ) − j α ([ n, q α ]; t ) . (S3)In the above definitions of v Mxc ([ n, q α ]; r t ) and j α,Mxc ([ n, q α ]; t ), the initial state dependence of the interacting Ψ and non-interacting Φ system has been dropped. For completeness, the definition of j α,Mxc ([ n, q α ]; t ) accounts for afunctional dependence on q α but this term can be calculated explicitly since it has no xc part as seen in Eq. (8). Thesimplified form of j α,Mxc is shown in Eq. (6).Through similar steps as in Eqs.(S6)-(S8), in first-order the solution of the Schr¨odinger-Kohn-Sham equation readsΦ( t ) (cid:39) ˆ U KS , ( t )Φ − i (cid:126) ˆ U KS , ( t ) (cid:90) tt dt (cid:48) ˆ H ( ext )KS ,I ( t (cid:48) ) ˆ U † KS , ( t )Φ . (S4)where ˆ U KS , = e − i ˆ H (0)KS t/ (cid:126) . Next, the bijective mapping between the interacting and non-interacting system that yieldsthe same density and photon coordinate is given as( v ( r t ) , j α ( t )) ←→ Ψ ( n ( r t ) , q α ( t )) ←→ Φ ( v s ( r t ) , j sα ( t )) , (S5)which can be inverted as ( v s ([ v, j α ]; r (cid:48) t (cid:48) ) , j sα ([ v, j α ]; t (cid:48) )). The response of the electronic subsystem due to the pertur-bations with the external pair ( v ( r t ) , j α ( t )) is δn ( r t ) = − i (cid:126) (cid:90) (cid:90) dτ d x (cid:90) (cid:90) dt (cid:48) d r (cid:48) (cid:104) Φ | [ˆ n I ( r t ) , ˆ n I ( x τ )] | Φ (cid:105) δv s ([ v, j α ]; x τ ) δv ( r (cid:48) t (cid:48) ) δv ( r (cid:48) t (cid:48) ) − i (cid:126) (cid:90) (cid:90) dτ d x (cid:88) α (cid:90) dt (cid:48) (cid:104) Φ | [ˆ n I ( r t ) , ˆ n I ( x τ )] | Φ (cid:105) δv s ([ v, j α ]; x τ ) δj α ( t (cid:48) ) δj α ( t (cid:48) ) . Where (cid:104) Φ | [ˆ n I ( r t ) , ˆ q α,I ( τ )] | Φ (cid:105) = 0 since both, electronic and photonic subsystems, are independent in the non-interacting system. From Eq. (S5), we have ( v s ([ n ]; r t ) , j sα ([ q α ]; t )) such that the above equation becomes δn ( r t ) = (cid:90) (cid:90) dτ d x (cid:90) (cid:90) dt (cid:48) d r (cid:48) (cid:90) (cid:90) dτ (cid:48) d y χ nn,s ( r t, x τ ) δv s ([ n ]; x τ ) δn ( y τ (cid:48) ) δn ([ v, j α ]; y τ (cid:48) ) δv ( r (cid:48) t (cid:48) ) δv ( r (cid:48) t (cid:48) )+ (cid:90) (cid:90) dτ d x (cid:88) α (cid:90) dt (cid:48) (cid:90) (cid:90) dτ (cid:48) d y χ nn,s ( r t, x τ ) δv s ([ n ]; x τ ) δn ( y τ (cid:48) ) δn ([ v, j α ]; y τ (cid:48) ) δj α ( t (cid:48) ) δj α ( t (cid:48) ) , (S6)where χ nn,s ( r t, x τ ) = ( − i/ (cid:126) )Θ( t − τ ) (cid:104) Φ | [ˆ n I ( r t ) , ˆ n I ( x τ )] | Φ (cid:105) is the non-interacting density-density response function.For clarity, the above density response is δn ( r t ) = δn v ( r t ) + δn j ( r t ), where ( δn v ( r t ) , δn j ( r t )) is the density responseto the external pair ( v ( r t ) , j α ( t )), respectively.Using Eqs.(S2) and (S3), we define the mean-field xc kernels as: f nMxc ([ n, q α ]; r t, r (cid:48) t (cid:48) ) = δv s ([ n ]; r t ) δn ( r (cid:48) t (cid:48) ) − δv ([ n, q α ]; r t ) δn ( r (cid:48) t (cid:48) ) , (S7) f q α Mxc ([ n, q α ]; r t, t (cid:48) ) = − δv ([ n, q α ]; r t ) δq α ( t (cid:48) ) , (S8) g nMxc ([ n, q α ]; t, r (cid:48) t (cid:48) ) = − δj α ([ n, q α ]; t ) δn ( r (cid:48) t (cid:48) ) , (S9) g q α (cid:48) Mxc ([ n, q α ]; t, t (cid:48) ) = δj sα ([ q α ]; t ) δq α (cid:48) ( t (cid:48) ) − δj α ([ n, q α ]; t ) δq α (cid:48) ( t (cid:48) ) , (S10)where δv s ([ n ]; r t ) δq α ( t (cid:48) ) = 0 = δj sα ([ q α ]; t ) δn ( r (cid:48) t (cid:48) ) . These kernels are the respective inverse of the interacting and non-interactingresponse functions.From Eq. (S6), density response to δv ( r t ) can be written in terms of the density-density response function given by χ nn ( r t, r (cid:48) t (cid:48) ) = (cid:90) (cid:90) dτ d x χ nn,s ( r t, x τ ) (cid:90) (cid:90) dτ (cid:48) d y f nMxc ([ n, q α ]; x τ, y τ (cid:48) ) δn ([ v, j α ]; y τ (cid:48) ) δv ( r t (cid:48) )+ (cid:90) (cid:90) dτ d x χ nn,s ( r t, x τ ) (cid:90) (cid:90) dτ (cid:48) d y δv ([ n, q α ]; x τ ) δn ( y τ (cid:48) ) δn ([ v, j α ]; y τ (cid:48) ) δv ( r (cid:48) t (cid:48) ) . Making the following substitution in the above equation (cid:90) (cid:90) d y dτ (cid:48) δv ([ n, q α ]; x τ ) δn ( y τ (cid:48) ) δn ([ v, j α ]; y τ (cid:48) ) δv ( r (cid:48) t (cid:48) ) = δ ( x − r (cid:48) ) δ ( τ − t (cid:48) ) − (cid:88) α (cid:90) dτ (cid:48) δv ([ n, q α ]; x τ ) δq α ( τ (cid:48) ) δq α ([ v, j α ]; τ (cid:48) ) δv ( r (cid:48) t (cid:48) ) , where δv ([ n, q α ]; x τ ) /δv ( r (cid:48) t (cid:48) ) = δ ( x − r (cid:48) ) δ ( τ − t (cid:48) ), we obtain the relation χ nn ( r t, r (cid:48) t (cid:48) ) = χ nn,s ( r t, r (cid:48) t (cid:48) ) + (cid:90) (cid:90) (cid:90) (cid:90) dτ d x dτ (cid:48) d y χ nn,s ( r t, x τ ) f nMxc ( x τ, y τ (cid:48) ) χ nn ( y τ (cid:48) , r (cid:48) t (cid:48) )+ (cid:88) α (cid:90) (cid:90) (cid:90) dτ d x dτ (cid:48) χ nn,s ( r t, x τ ) f q α Mxc ( x τ, τ (cid:48) ) χ q α n ( τ (cid:48) , r (cid:48) t (cid:48) ) . (S11)Next, the density response to δj α ( t ) in Eq. (S6) is expressed in terms of the response function as χ nq α ( r t, t (cid:48) ) = (cid:90) (cid:90) dτ d x χ nn,s ( r t, x τ ) (cid:90) (cid:90) dτ (cid:48) d y f nMxc ( x τ, y τ (cid:48) ) δn ([ v, j α ]; y τ (cid:48) ) δj α ( t (cid:48) )+ (cid:90) (cid:90) dτ d x χ nn,s ( r t, x τ ) (cid:90) (cid:90) dτ (cid:48) d y δv ([ n, q α ]; x τ ) δn ( y τ (cid:48) ) δn ([ v, j α ]; y τ (cid:48) ) δj α ( t (cid:48) ) . Using the relation (obtained from δv ([ n, q α ]; x τ ) /δj α ( t (cid:48) )) (cid:90) (cid:90) d y dτ (cid:48) δv ([ n, q α ]; x τ ) δn ( y τ (cid:48) ) δn ([ v, j α ]; y τ (cid:48) ) δj α ( t (cid:48) ) = − (cid:88) α (cid:48) (cid:90) dτ (cid:48) δv ([ n, q α ]; x τ ) δq α (cid:48) ( τ (cid:48) ) δq α (cid:48) ([ v, j α ]; τ (cid:48) ) δj α ( t (cid:48) ) , the response function is given as χ nq α ( r t, t (cid:48) ) = (cid:90) (cid:90) (cid:90) (cid:90) dτ d x dτ (cid:48) d y χ nn,s ( r t, x τ ) f nMxc ( x τ, y τ (cid:48) ) χ nq α ( y τ (cid:48) , t (cid:48) )+ (cid:88) α (cid:48) (cid:90) (cid:90) (cid:90) dτ d x dτ (cid:48) χ nn,s ( r t, x τ ) f q α (cid:48) Mxc ( x τ, τ (cid:48) ) χ q α (cid:48) q α ( τ (cid:48) , t (cid:48) ) . (S12)Similarly, the response to the photonic subsystem to linear perturbations from the external pair ( v ( r t ) , j α ( t )) is δq α ( t ) = − i (cid:126) (cid:88) β (cid:90) tt dτ ω β (cid:104) Φ | [ q α,I ( t ) , q β,I ( τ )] | Φ (cid:105) (cid:90) (cid:90) dt (cid:48) d r (cid:48) δj sβ ([ v, j α ]; τ ) δv ( r (cid:48) t (cid:48) ) δv ( r (cid:48) t (cid:48) ) − i (cid:126) (cid:88) β (cid:90) tt dτ ω β (cid:104) Φ | [ q α,I ( t ) , q β,I ( τ )] | Φ (cid:105) (cid:88) α (cid:48) (cid:90) dt (cid:48) δj sβ ([ v, j α ]; τ ) δj α (cid:48) ( t (cid:48) ) δj α (cid:48) ( t (cid:48) ) , where (cid:104) Φ | [ˆ q α,I ( t ) , ˆ n I ( x τ )] | Φ (cid:105) = 0 in the non-interacting system. By defining the non-interacting photon-photonresponse function as χ q α q β,s ( t, τ ) = ( − i/ (cid:126) )Θ( t − τ )(1 /ω β ) (cid:104) Φ | [ q α,I ( t ) , q β,I ( τ )] | Φ (cid:105) and using Eq. (S5), where we have( v s ([ n ]; r t ) , j sα ([ q α ]; t )), the response can be written as δq α ( t ) = (cid:88) β (cid:90) dτ χ q α q β,s ( t, τ ) (cid:88) β (cid:48) (cid:90) (cid:90) (cid:90) dt (cid:48) d r (cid:48) dτ (cid:48) δj sβ ([ q α ]; τ ) δq β (cid:48) ( τ (cid:48) ) δq β (cid:48) ([ v, j α ]; τ (cid:48) ) δv ( r (cid:48) t (cid:48) ) δv ( r (cid:48) t (cid:48) )+ (cid:88) β (cid:90) dτ χ q α q β,s ( t, τ ) (cid:88) α (cid:48) ,β (cid:48) (cid:90) (cid:90) dt (cid:48) dτ (cid:48) δj sβ ([ q α ]; τ ) δq β (cid:48) ( τ (cid:48) ) δq β (cid:48) ([ v, j α ]; τ (cid:48) ) δj α (cid:48) ( t (cid:48) ) δj α (cid:48) ( t (cid:48) ) . (S13)The above response of the displacement field is δq α ( t ) = δq α,v ( t ) + δq α,j ( t ), where ( δq α,v ( t ) , δq α,j ( t )) is the responseto the external pair ( v ( r t ) , j α ( t )), respectively.From Eq. (S13), the field response to δv ( r t ) can be written in terms of the photon-density response function as χ q α n ( t, r (cid:48) t (cid:48) ) = (cid:88) β (cid:90) dτ χ q α q β,s ( t, τ ) (cid:88) β (cid:48) (cid:90) dτ (cid:48) g q β (cid:48) Mxc ( τ, τ (cid:48) ) χ q β (cid:48) n ( τ (cid:48) , r (cid:48) t (cid:48) )+ (cid:88) β (cid:90) dτ χ q α q β,s ( t, τ ) (cid:88) β (cid:48) (cid:90) dτ (cid:48) δj β ([ n, q α ]; τ ) δq β (cid:48) ( τ (cid:48) ) δq β (cid:48) ([ v, j α ]; τ (cid:48) ) δv ( r (cid:48) t (cid:48) ) . Using the relation (obtained from δj β ([ n, q α ]; τ ) /δv ( r (cid:48) t (cid:48) )) (cid:88) β (cid:48) (cid:90) dτ (cid:48) δj β ([ n, q α ]; τ ) δq β (cid:48) ( τ (cid:48) ) δq β (cid:48) ([ v, j α ]; τ (cid:48) ) δv ( r t (cid:48) ) = − (cid:90) (cid:90) dτ (cid:48) d y δj β ([ n, q α ]; τ ) δn ( y τ (cid:48) ) δn ([ v, j α ]; y τ (cid:48) ) δv ( r (cid:48) t (cid:48) ) , the response function is given as χ q α n ( t, r (cid:48) t (cid:48) ) = (cid:88) β (cid:90) dτ (cid:90) (cid:90) dτ (cid:48) d y χ q α q β,s ( t, τ ) g n β Mxc ( τ, y τ (cid:48) ) χ nn ( y τ (cid:48) , r (cid:48) t (cid:48) ) , (S14)where g n β Mxc = g n β M and g q α Mxc = 0 as determined from the equation of motion for the displacement field. Also, fromEq. (S13), field response to δj α can be written in terms of the photon-photon response function as χ q α q α (cid:48) ( t, t (cid:48) ) = (cid:88) β (cid:90) dτ χ q α q β,s ( t, τ ) (cid:88) β (cid:48) (cid:90) dτ (cid:48) g q β (cid:48) Mxc ( τ, τ (cid:48) ) χ q β (cid:48) q α (cid:48) ( τ (cid:48) , t (cid:48) )+ (cid:88) β (cid:90) dτ χ q α q β,s ( t, τ ) (cid:88) β (cid:48) (cid:90) dτ (cid:48) δj β ([ n, q α ]; τ ) δq β (cid:48) ( τ (cid:48) ) δq β (cid:48) ([ n, q α ]; τ (cid:48) ) δj α (cid:48) ( t (cid:48) ) . Making the following substitution (where δj β ([ n, q α ]; τ ) /δj α (cid:48) ( t (cid:48) ) = δ ( τ − t (cid:48) ) δ β,α (cid:48) ) in the above equation (cid:88) β (cid:48) (cid:90) dτ (cid:48) δj β ([ n, q α ]; τ ) δq β (cid:48) ( τ (cid:48) ) δq β (cid:48) ([ v, j α ]; τ (cid:48) ) δj α (cid:48) ( t (cid:48) ) = δ ( τ − t (cid:48) ) δ β,α (cid:48) − (cid:90) (cid:90) dτ (cid:48) d x δj β ([ n, q α ]; τ ) δn ( x τ (cid:48) ) δn ([ v, j α ]; x τ (cid:48) ) δj α (cid:48) ( t (cid:48) ) , yields the photon-photon response function χ q α q α (cid:48) ( t, t (cid:48) ) = χ q α q α (cid:48) ,s ( t, t (cid:48) ) + (cid:88) β (cid:90) (cid:90) (cid:90) dτ dτ (cid:48) d x χ q α q β,s ( t, τ ) g n β Mxc ( τ, x τ (cid:48) ) χ nq α (cid:48) ( x τ (cid:48) , t (cid:48) ) , (S15)where g q β (cid:48) Mxc = 0 since j α,M in Eq. (6) has no functional dependency on q α . Appendix S4: Matrix formulation of QEDFT response equations In this section we present a matrix formulation of non-relativistic QEDFT response equations which in the no-coupling limit reduces to Casida equation. Through a Fourier transform of Eqs.(S11)-(S12) and Eqs.(S14)-(S15) andmaking a substitution into Eqs.(C8)-(C11), we express the responses in the following form: δn v ( r , ω ) = (cid:88) i,a (cid:104) ϕ a ( r ) ϕ ∗ i ( r ) P (1) ai,v ( ω ) + ϕ i ( r ) ϕ ∗ a ( r ) P (1) ia,v ( ω ) (cid:105) , (S1) δn j ( r , ω ) = (cid:88) i,a (cid:104) ϕ a ( r ) ϕ ∗ i ( r ) P (1) ai,j ( ω ) + ϕ i ( r ) ϕ ∗ a ( r ) P (1) ia,j ( ω ) (cid:105) , (S2) δq α,v ( ω ) = L (1) α,v, − ( ω ) + L (1) α,v, + ( ω ) , (S3) δq α,j ( ω ) = L (1) α,j, − ( ω ) + L (1) α,j, + ( ω ) . (S4)Here, the subscripts ( v, j ) on the first-order responses P (1) ia,v , P (1) ia,j , P (1) ai,v , P (1) ai,j , L (1) α,v, ± and L (1) α,j, ± shows to whatexternal perturbations ( δv ( r , t ) , δj α ( t )) is being considered to induce the coupled responses. In defining Eqs.(S1)-(S4), we used the static KS orbitals in the Lehmann spectral representation of χ nn,s ( r , r (cid:48) , ω ) and photon-photonresponse function χ q α q α,s ( ω ) for a single-photon in Fock number basis are given as χ nn,s ( r , r (cid:48) , ω ) = (cid:88) i,a (cid:18) ψ a ( r ) ψ i ( r (cid:48) ) ψ ∗ i ( r ) ψ ∗ a ( r (cid:48) ) ω − ( (cid:15) a − (cid:15) i ) + iη − ψ i ( r ) ψ a ( r (cid:48) ) ψ ∗ a ( r ) ψ ∗ i ( r (cid:48) ) ω + ( (cid:15) a − (cid:15) i ) + iη (cid:19) ,χ q α q α,s ( ω ) = 12 ω α (cid:18) ω − ω α + iη − ω + ω α + iη (cid:19) . where the summations over occupied and unoccupied Kohn-Sham orbitals are performed according to (cid:80) i = (cid:80) Ni =1 and (cid:80) a = (cid:80) ∞ a = N +1 and from here on lim η → + is implied. The first-order responses P (1) ia,v , P (1) ia,j , P (1) ai,v , P (1) ai,j , L (1) α,v, ± and L (1) α,j, ± are given by [ ω − ω ai ] P (1) ai,v ( ω ) = (cid:90) d r ϕ i ( r ) ϕ ∗ a ( r ) δv (1)KS ,v ( r , ω ) , (S5)[ ω + ω ai ] P (1) ia,v ( ω ) = − (cid:90) d r ϕ a ( r ) ϕ ∗ i ( r ) δv (1)KS ,v ( r , ω ) , (S6)[ ω − ω ai ] P (1) ai,j ( ω ) = (cid:90) d r ϕ i ( r ) ϕ ∗ a ( r ) δv (1)KS ,j ( r , ω ) , (S7)[ ω + ω ai ] P (1) ia,j ( ω ) = − (cid:90) d r ϕ a ( r ) ϕ ∗ i ( r ) δv (1)KS ,j ( r , ω ) , (S8)[ ω − ω α ] L (1) α,v, − ( ω ) = 12 ω α δj (1) α, KS ,v ( ω ) , (S9)[ ω + ω α ] L (1) α,v, + ( ω ) = − ω α δj (1) α, KS ,v ( ω ) , (S10)[ ω − ω α ] L (1) α,j, − ( ω ) = 12 ω α δj (1) α, KS ,j ( ω ) , (S11)[ ω + ω α ] L (1) α,j, + ( ω ) = − ω α δj (1) α, KS ,j ( ω ) , (S12)where ω ai = ( (cid:15) a − (cid:15) i ) and the respective effective potentials and currents ( δv s,ν ( r , ω ) , j sα,ν ( ω )) as δv (1)KS ,v ( r , ω ) = δv ( r , ω ) + (cid:90) d r (cid:48) f nMxc ( r , r (cid:48) , ω ) δn v ( r (cid:48) , ω ) + (cid:88) α f q α Mxc ( r , ω ) δq α,v ( ω ) , (S13) δv (1)KS ,j ( r , ω ) = (cid:90) d r (cid:48) f nMxc ( r , r (cid:48) , ω ) δn j ( r (cid:48) , ω ) + (cid:88) α f q α Mxc ( r , ω ) δq α,j ( ω ) , (S14) δj (1) α, KS ,v ( ω ) = (cid:90) d r g n α M ( r ) δn v ( r , ω ) , (S15) δj (1) α, KS ,j ( ω ) = δj α ( ω ) + (cid:90) d r g n α M ( r ) δn j ( r , ω ) . (S16)The mean-field kernel is given by g n α M ( r ) = − ω α λ α · r . As stated above, the subscripts ( v, j ) on the responses, KSpotentials and currents signifies as to what external perturbations ( δv ( r , t ) , δj α ( t )) is being considered. The Kohn-Sham scheme of QEDFT decouples the interacting system such that the responses are paired as ( δn v ( r , ω ) , δq α,v ( ω ))due to δv ( r , ω ) and ( δn j ( r , ω ) , δq α,j α ( ω )) due to δj α ( ω ). Therefore, substituting Eqs.(S13) and (S15) into Eqs.(S5)-(S6) and Eqs.(S9)-(S10) and after some simplification, we obtain (cid:88) j,b [ δ ab δ ij ( ω ai − ω ) + K ai,jb ( ω )] P (1) bj,v ( ω ) + K ai,bj ( ω ) P (1) jb,v ( ω ) + (cid:88) α δ ab δ ij M α,bj ( ω ) (cid:16) L (1) α,v, − ( ω ) + L (1) α,v, + ( ω ) (cid:17) = − v ai ( ω ) , (S17) (cid:88) j,b [ δ ab δ ij ( ω ai + ω ) + K ia,bj ( ω )] P (1) jb,v ( ω ) + K ia,jb ( ω ) P (1) bj,v ( ω ) + (cid:88) α δ ab δ ij M α,jb ( ω ) (cid:16) L (1) α,v, − ( ω ) + L (1) α,v, + ( ω ) (cid:17) = − v ia ( ω ) , (S18)[ ω α − ω ] L (1) α,v, − ( ω ) + (cid:88) jb (cid:104) N α,jb P (1) bj,v ( ω ) + N α,bj P (1) jb,v ( ω ) (cid:105) = 0 , (S19)[ ω α + ω ] L (1) α,v, + ( ω ) + (cid:88) jb (cid:104) N α,jb P (1) bj,v ( ω ) + N α,bj P (1) jb,v ( ω ) (cid:105) = 0 , (S20)Also, substituting Eqs.(S14) and (S16) into Eqs.(S7)-(S8) and Eqs.(S11)-(S12) and after some simplification, we obtain (cid:88) j,b δ ab δ ij (cid:34) (( ω ai − ω ) + K ai,jb ( ω )) P (1) bj,j ( ω ) + K ai,bj ( ω ) P (1) jb,j ( ω ) + (cid:88) α M α,bj ( ω ) (cid:104) L (1) α,j, − ( ω ) + L (1) α,j, + ( ω ) (cid:105)(cid:35) = 0 , (S21) (cid:88) j,b δ ab δ ij (cid:34) (( ω ai + ω ) + K ia,bj ( ω )) P (1) jb,j ( ω ) + K ia,jb ( ω ) P (1) bj,j ( ω ) + (cid:88) α M α,jb ( ω ) (cid:104) L (1) α,j, − ( ω ) + L (1) α,j, + ( ω ) (cid:105)(cid:35) = 0 , (S22)[ ω α − ω ] L (1) α,j, − ( ω ) + (cid:88) jb (cid:104) N α,jb P (1) bj,j ( ω ) + N α,bj P (1) jb,j ( ω ) (cid:105) = − ω α δj α ( ω ) , (S23)[ ω + ω α ] L (1) α,j, + ( ω ) + (cid:88) jb (cid:104) N α,jb P (1) bj,j ( ω ) + N α,bj P (1) jb,j ( ω ) (cid:105) = − ω α δj α ( ω ) , (S24)where we defined the coupling matrices K ai,jb ( ω ) = (cid:90) (cid:90) d r d y ϕ i ( r ) ϕ ∗ a ( r ) f nMxc ( r , y , ω ) ϕ b ( y ) ϕ ∗ j ( y ) , (S25) M α,ai ( ω ) = (cid:90) d r ϕ i ( r ) ϕ ∗ a ( r ) f q α Mxc ( r , ω ) , (S26) N α,ia = 12 ω α (cid:90) d r ϕ ∗ i ( r ) ϕ a ( r ) g n α M ( r ) , (S27)and v ia ( ω ) = (cid:90) d r ϕ ∗ i ( r ) δv ( r , ω ) ϕ a ( r ) . (S28)The coupling matrix N α,ia has no frequency dependence since this is just the mean-field kernel of the photon modes. Wenow introduce the following abbreviations L ( ω ) = δ ab δ ij ( (cid:15) a − (cid:15) i ) + K ai,jb ( ω ), K ( ω ) = K ai,jb ( ω ), M ( ω ) = M α,bj ( ω ), N = N α,bj , X ( ω ) = P (1) bj,v ( ω ), Y ( ω ) = P (1) jb,v ( ω ), X ( ω ) = P (1) bj,j ( ω ), Y ( ω ) = P (1) jb,j ( ω ), A ( ω ) = L (1) α,v, − ( ω ), B ( ω ) = L (1) α,v, + ( ω ), A ( ω ) = L (1) α,j, − ( ω ), B ( ω ) = L (1) α,j, + ( ω ), V ( ω ) = − v ai ( ω ), J α ( ω ) = − δj α ( ω )2 ω α .Using these notations, we cast Eqs.(S17)-(S20) and Eqs.(S21)-(S24) into two matrix equations given by L ( ω ) K ( ω ) M ( ω ) M ( ω ) K ∗ ( ω ) L ( ω ) M ∗ ( ω ) M ∗ ( ω ) N N ∗ ω α N N ∗ ω α + ω − − X ( ω ) Y ( ω ) A ( ω ) B ( ω ) = V ( ω ) V ∗ ( ω )00 (S29) L ( ω ) K ( ω ) M ( ω ) M ( ω ) K ∗ ( ω ) L ( ω ) M ∗ ( ω ) M ∗ ( ω ) N N ∗ ω α N N ∗ ω α + ω − − X ( ω ) Y ( ω ) A ( ω ) B ( ω ) = J α ( ω ) J α ( ω ) (S30)0Next, we argue that the right hand side of the above matrices remains finite as the frequency ω approaches the exactexcitation frequencies ω → Ω q of the interacting system while the density and displacement field responses on theleft hand side has poles at the true excitation frequencies Ω q . This allows us to cast Eq. (S29) and Eq. (S30) into aneigenvalue problem L (Ω q ) K (Ω q ) M (Ω q ) M (Ω q ) K ∗ (Ω q ) L (Ω q ) M ∗ (Ω q ) M ∗ (Ω q ) N N ∗ ω α N N ∗ ω α X (Ω q ) Y (Ω q ) A (Ω q ) B (Ω q ) = Ω q − − X (Ω q ) Y (Ω q ) A (Ω q ) B (Ω q ) (S31) (Ω q ) K (Ω q ) M (Ω q ) M (Ω q ) K ∗ (Ω q ) L (Ω q ) M ∗ (Ω q ) M ∗ (Ω q ) N N ∗ ω α N N ∗ ω α X (Ω q ) Y (Ω q ) A (Ω q ) B (Ω q ) = Ω q − − X (Ω q ) Y (Ω q ) A (Ω q ) B (Ω q ) (S32)It is convenient to cast Eqs.(S31) and (S32) into a Hermitian eigenvalue problem which is given by (cid:32) U VW ω α (cid:33) (cid:32) E P (cid:33) = Ω q (cid:32) E P (cid:33) , (S33) (cid:32) U VW ω α (cid:33) (cid:32) E P (cid:33) = Ω q (cid:32) E P (cid:33) , (S34)where we assumed real-valued orbitals, i.e., K = K ∗ , M = M ∗ and N = N ∗ , and the matrices are given by U = ( L − K ) / ( L + K )( L − K ) / , V = 2( L − K ) / M / N / ω / α , W = 2 ω / α N / M / ( L − K ) / , and theeigenvectors are E = N / ( L − K ) − / ( X + Y ) and P = M / ω − / α ( A + B ).The pseudo-eigenvalue problem of Eqs.(S33) and (S34) is the final form of QEDFT matrix equation for obtainingexact excitation frequencies and oscillator strengths. Appendix S5: Oscillator Strengths In this section, we derive the oscillator strengths resulting from the eigenvectors of the pseudo-eigenvalue problemof Eqs.(S33) and (S34). Multiplying out Eq. (S29), we write the matrix equation in the form( L + K )( X + Y ) + 2 M ( A + B ) − ω ( X − Y ) = − v , ( L − K )( X − Y ) − ω ( X + Y ) = 0 , N ( X + Y ) + ω α ( A + B ) − ω ( A − B ) = 0 ,ω α ( A − B ) − ω ( A + B ) = 0 . From here on we set S = ( L − K ), the above pair of equations now becomes S ( L + K ) E + 2 SM P − ω E = − S v , ω α N E + ω α P − ω P = 0 . This can be written in matrix form as (cid:34)(cid:32) S ( L + K ) 2 SM ω α N ω α (cid:33) − ω (cid:32) (cid:33)(cid:35) (cid:32) E P (cid:33) = − (cid:32) S v (cid:33) , (S1)where E = X + Y and P = A + B . We perform the same steps as above to make the nonlinear eigenvalueproblem Hermitian and obtain (cid:2) C − ω (cid:3) (cid:32) N / S − / E M / ω − / α P (cid:33) = − (cid:32) N / S / v (cid:33) , (S2)1where C = (cid:32) U VW ω α (cid:33) . We determine the vectors given as E = − S / (cid:2) C − ω (cid:3) − S / v , (S3) P = − ω / α M − / (cid:2) C − ω (cid:3) − N / S / v . (S4)When Z I is normalized, we can use the spectral expansion to get (cid:2) C − ω (cid:3) − = (cid:88) I Z I Z † I Ω I − ω , (S5)where Z I = (cid:32) E I P I (cid:33) . The oscillator strength for the density-density response function which is related to the dynamicpolarizability is given in Eq.(C28). 1. Oscillator strength for the photon-matter response function Next, we substitute the expression of the spectral expansion Eq. (S5) in Eq. (S4) and by substituting P in Eq. (S2)yields δq α,v ( ω ) = (cid:88) I (cid:40) ω / α M − / Z I Z † I N / S / ω − Ω I (cid:41) v ( ω ) . The oscillator strength is given by f pnI,α = 2 ω / α M − / Z I Z † I N / S / . (S6)Also, from Eq.(C9) and using the Lehmann representation of the response function χ q α n ( r (cid:48) , ω ) the response δq α,v ( ω )is given by δq α,v ( ω ) = (cid:90) d r (cid:48) (cid:88) k (cid:20) k (cid:104) Ψ | ˆ q α | Ψ k (cid:105)(cid:104) Ψ k | ˆ n ( r (cid:48) ) | Ψ (cid:105) ω − Ω k (cid:21) δv ( r (cid:48) , ω ) , The oscillator strength of Eq.(S6) can be expressed as matrix elements of the internal pair (ˆ n ( r ) , ˆ q α ) as f α,k ( r (cid:48) ) = 2Ω k (cid:104) Ψ | ˆ q α | Ψ k (cid:105)(cid:104) Ψ k | ˆ n ( r (cid:48) ) | Ψ (cid:105) ≡ f pnI,α . (S7) 2. Oscillator strength for the matter-photon response function Following similar steps as above with Eq. (S30) we obtain E = − S / N − / (cid:2) C − ω (cid:3) − M / ω / α J (cid:48) α , (S8) P = − ω / α (cid:2) C − ω (cid:3) − ω / α J (cid:48) α . (S9)where J (cid:48) α ( ω ) = j α ( ω )2 ω α and J α ( ω ) = − J (cid:48) α ( ω ). By substituting the spectral expansion Eq. (S5) in E and furthersubstituting in Eq. (S3) yields δn j ( r , ω ) = − (cid:88) ia,I Φ ia S / N − / Z I Z † I M / ω / α Φ ai (Ω I − ω ) J (cid:48) α ( ω ) . Following a similar procedure as above, we express the density response to the external charge current as δn j ( r , ω ) = (cid:88) I (cid:40) Φ ia S / N − / Z I Z † I M / ω / α Φ ia ω − Ω I (cid:41) j α ( ω ) ω α , ia ( r ) = ϕ ∗ i ( r ) ϕ a ( r ) and the oscillator strength is given by f npI,α = 1 ω α Φ ia S / N − / Z I Z † I M / ω / α Φ ia . (S10)From Eq.(C10) and using the Lehmann representation of the response function χ nq α ( r , ω ), the response δn j ( r , ω ) isgiven by δn j ( r , ω ) = (cid:88) α,k (cid:20) k (cid:104) Ψ | ˆ n ( r ) | Ψ k (cid:105)(cid:104) Ψ k | ˆ q α | Ψ (cid:105) ω − Ω k (cid:21) δj α ( ω ) ω α , The oscillator strength of Eq.(S10) can be expressed as matrix elements of the internal pair (ˆ n ( r ) , ˆ q α ) as f k,α ( r ) = 2Ω k (cid:104) Ψ | ˆ n ( r ) | Ψ k (cid:105)(cid:104) Ψ k | ˆ q α | Ψ (cid:105) ≡ f npI,α . (S11) 3. Oscillator strength for the photon-photon response function We define a collective photon coordinate for the α modes Q = (cid:80) α q α (in analogy with R = (cid:80) i e r i ). By perturbingthe photon field through the photon coordinate with an external charge current j α ( ω ), we induce a polarization ofthe field of mode α which we denote as Q ( ω ) = (cid:80) α β α ( ω ) j α ( ω ). Where β α ( ω ) is the polarizability of field of the α mode. To first-order, the collective coordinate is given by δQ ( t ) = (cid:88) α δq α ( t ) . (S12)The field polarizability in frequency space can be written as β α ( ω ) = (cid:88) α (cid:48) δq α ( ω ) δj α (cid:48) ( ω ) . (S13)By substituting Eq. (S9) in Eq. (S4) and using the spectral expansion yields δq α,j ( ω ) = − (cid:88) I ω / α Z I Z † I ω / α Ω I − ω J (cid:48) α . By substituting the above relation in Eq. (S13) we obtain β α ( ω ) = − (cid:88) α (cid:48) (cid:88) I ω / α Z I Z † I ω / α Ω I − ω δj α ( ω ) / ω α δj α (cid:48) ( ω ) , which simplifies to β α ( ω ) = − (cid:88) I ω α ω / α Z I Z † I ω / α Ω I − ω . (S14)Eq. (S14) is the field polarizability analogous to the atomic polarizability tensor of Eq. (C26). As in Eq.(C27) in whichthe molecular isotropic polarizability, α ( ω ) is defined as the mean value of three diagonal elements of the polarizabilitytensor, i.e., α ( ω ) = 1 / α xx ( ω ) + α yy ( ω ) + α zz ( ω )), we analogously define an absorption cross section of the fieldgiven by ˜ σ α ( ω ) ≡ πωc I m Tr β α ( ω ) / . (S15)For the oscillator strength, from Eq.(C11) and using the Lehmann representation of the response function χ q α q α (cid:48) ( ω )the response δq α,j ( ω ) is given by δq α,j ( ω ) = (cid:88) α (cid:48) ,k (cid:20) k (cid:104) Ψ | ˆ q α | Ψ k (cid:105)(cid:104) Ψ k | ˆ q α (cid:48) | Ψ (cid:105) ω − Ω k (cid:21) δj α (cid:48) ( ω ) ω α (cid:48) . We find the oscillator strength f ppI,α = 13 ω α (cid:12)(cid:12)(cid:12) Z † I ω / α (cid:12)(cid:12)(cid:12) = 23 Ω I (cid:88) α (cid:48) ω α (cid:48) (cid:104) Ψ | ˆ q α | Ψ I (cid:105)(cid:104) Ψ I | ˆ q α (cid:48) | Ψ (cid:105) ..