Ab Initio Calculations of Quantum Light–Matter Interactions in General Electromagnetic Environments

The emerging field of strongly coupled light–matter systems has drawn significant attention in recent years because of the prospect of altering both the physical and chemical properties of molecules and materials. Because this emerging field draws on ideas from both condensed-matter physics and quantum optics, it has attracted the attention of theoreticians from both fields. While the former often employ accurate descriptions of the electronic structure of the matter, the description of the electromagnetic environment is often oversimplified. In contrast, the latter often employs sophisticated descriptions of the electromagnetic environment while using oversimplified few-level approximations of the electronic structure. Both approaches are problematic because the oversimplified descriptions of the electronic system are incapable of describing effects such as light-induced structural changes in the electronic system, while the oversimplified descriptions of the electromagnetic environments can lead to unphysical predictions because the light–matter interactions strengths are misrepresented. In this work, we overcome these shortcomings and present the first method which can quantitatively describe both the electronic system and general electromagnetic environments from first principles. We realize this by combining macroscopic QED (MQED) with Quantum Electrodynamical Density-Functional Theory. To exemplify this approach, we consider the example of an absorbing spherical cavity and study the impact of different parameters of both the environment and the electronic system on the transition from weak-to-strong coupling for different aromatic molecules. As part of this work, we also provide an easy-to-use tool to calculate the cavity coupling strengths for simple cavity setups. Our work is a significant step toward parameter-free ab initio calculations for strongly coupled quantum light–matter systems and will help bridge the gap between theoretical methods and experiments in the field.


Introduction
When quantum light and matter interact strongly, hybrid light-matter polariton states emerge.[4][5][6][7][8] Polariton formation can be realized in resonant electromagnetic environments, the paradigmatic example of which is a Fabry-Perot cavity.In some setups, it is possible to enhance the vacuum fluctuations of the electromagnetic field sufficiently to drive the formation of polaritons even in the absence of actual light in the cavity.To realize such large lightmatter coupling strengths it is, however, necessary to go beyond the paradigmatic cavity setup, consisting of two parallel mirrors, and often requires intricate meta-and nano-optical setups.Examples include plasmonic nanocavities, 9 optical cavities using four-wave mixing schemes , 10 metasurface systems, 11 self-assembled Casimir microcavities, 12 deep-strong coupling in plasmonic nanoparticles, 13 and many more. 1,14These cavity setups will generally feature complicated mode structures, especially with losses present, and therefore require descriptions beyond simple single mode descriptions.
Because this emerging field lives at the interface between condensed matter physics and quantum optics, a proper theoretical treatment simultaneously requires a quantitative description of the both electronic system and the quantum electromagnetic environment.However, to date, no method can properly account for both sides of the problem.Existing methods either apply few-level approximations to the electronic system, or oversimplified descriptions of the electromagnetic environment essentially treating the light-matter coupling strengths as free parameters.This is problematic because the light-matter interaction is very sensitive to both phase matching and energetic alignment of field and matter, as well as the spatial overlap of the electronic wave functions and the modes of the electromagnetic field. 15Furthermore, simplified descriptions of the electronic system are incapable of describing effects such as light-induced structural changes in the electronic system and are generally not applicable in the ultra-and deep strong coupling regimes.Finally, the oversimplified descriptions of the electromagnetic environments can lead to unphysical predictions because the light-matter interaction strengths are misrepresented. 15 this paper we present the first method that accounts quantitatively for both the electro-magnetic environment and the full electronic structure of matter by combining Macroscopic Quantum Electrodynamics (MQED) 16 with Quantum Electrodynamical Density-functional Theory (QEDFT). 17We have previously shown how standard DFT can be combined with MQED to provide a quantitative, first principles description of the quantum light-matter interactions for real cavity setups in the weak and intermediate coupling regimes beyond the dipole approximation. 15However, the previous work relied on a wave-function ansatz which only considered a subset of the electronic structure.This is expected to become problematic in the ultra-and deep strong coupling regimes.Furthermore, our previous method would also fail capture light-induced structural changes in the electronic system.In this work, we overcome these limitations and present a general method that is applicable in all regimes of light-matter coupling.Our new methodology allows us to study the interaction of the full electronic structure of electronic systems with realistic electromagnetic environments.We exemplify our approach on a spherical micro-cavity and highlight that the intricate interplay between the cavity geometry and material composition, and the electronic structure of the molecules has a profound impact on the light-matter coupling and the transition from weak to strong coupling.These results highlight the need for a quantitative description of both the electromagnetic environment and molecular system.This work is a significant step towards parameter free ab initio calculations for strongly coupled quantum light-matter systems and will help bridge the gap between theoretical methods and experiments in the field.

Macroscopic Quantum Electrodynamics
Macroscopic Quantum Electrodynamics (MQED) is a framework for quantizing the electromagnetic field in the presence of arbitrary absorbing or dispersing environments. 15,16,18The central object in MQED is the classical Green's function that solves the Helmholtz equation for a point source, the so-called dyadic Green's function (DGF), 19 Here c is the speed of light, ω is the angular frequency, ϵ(r, ω) and κ(r, ω) = µ −1 (r, ω) are the spatially dependent dielectric function and inverse magnetic permability respectively, and I is the unit dyad.The DGF is of central importance to the quantized theory of electromagnetic fields in lossy environments because it simultaneously carries the information about the electromagnetic boundary conditions, and serves as a projector from the coupled light-matter system onto the electromagnetic degrees of freedom. 16r a spatially local magnetoelectric medium in the nonrelativistic limit, the MQED expansion of the electric field in the Power-Zienau-Woolley (PZW) frame 20,21 (multipolar gauge) can be written as, 16,18 Here fλ (r ′ , ω) are the spatially resolved polaritonic field operators of MQED which fulfill the commutation relations of the quantum harmonic oscillator, and G e (r, r ′ , ω) and G m (r, r ′ , ω) are the electric and magnetic components of the DGF respectively, In the following we neglect magnetic interactions and consider the coupling between light and matter within the dipole approximation.Therefore, if we consider a set of emitters i with positions (centers of charge) r i , the interaction only samples the electromagnetic field at these positions.In this sense, the full electric field Ê(r) in Eq. 2 contains more information than strictly necessary to describe the light-matter interaction completely.As discussed in Refs, 18,22,23 it is therefore possible to arrive at a significantly more compact expression by alternatively expanding the electric field in terms of a set of N explicitly orthogonalized bright modes at each frequency, Here Ĉ( †) i (ω) destroys (creates) a photon in the ith bright mode.E i (r, ω) describes the spatial mode function of the electric field associated with mode i.The normalisation factor G j (ω) is the square root of the dipole spectral density, 24 Finally, the matrix V ij (ω) is the transformation matrix which obeys V (ω)S(ω)V † (ω) = I where and it is a result of the mode-orthogonalization inherent to the emitter-centered representation. 18The number of emitter centered modes per frequency, N , is equal to the number of emitter positions times the number of dipole orientations considered.In this work, we consider a single emitter position r 0 and the full three-dimensional space of dipole orientations.
This results in three emitter-centered modes per frequency.Importantly, while the bright mode representation of the MQED field is able to account for the spatial dependence of the inter-emitter interaction, the local coupling of each emitter to the field is described within the dipole approximation.Even though the electric field in Eq. 6 retains spatial dependence it is therefore not a fully beyond dipole approximation representation of the MQED field.We note that the description of the local coupling within the dipole approximation is in principle problematic for the high frequency modes which are inherently included in Eq. 6.Because no real electronic system is truly point-like, and the field expansion in principle includes modes of arbitrarily large frequencies, the coupling of the emitters to some of these modes would inherently require a beyond dipole approximation.However, as discussed further below, it is necessary in practice to truncate the mode expansion at some upper frequency, and for the modes included in the actual calculations the dipole approximation is justified.6][27] Therefore we write (see Supplementary Note A), where H Mat is the standard Coulomb gauge matter Hamiltonian describing the electronic system, and λ i (ω) is the cavity field strength of the i'th bright mode at frequency ω, where r 0 denotes the electronic center of charge and E i (r 0 , ω) is defined in terms of the dyadic Greens function in Eq. 7. The cavity field strengths of an arbitrary electromagnetic environment are thus fully determined by the DGF.We have further introduced the photon field quantities pi and qi that are connected to the magnetic-and electric field in their corresponding mode and are given explicitly by qi (ω) = h 2ω 1/2 ( Ĉi (ω) + Ĉ † i (ω)) and pi (ω) = hω 2 1/2 ( Ĉi (ω) − Ĉ † i (ω)).In terms of these new quantities, the electric field expansion at the center of charge r 0 reads, The light-matter interaction further contains the dipole moment operator for N e electrons with position ri , R = Ne i=1 ri , which gives rise to an explicit electron-photon interaction and the dipole-self energy term ( ).While the dipole self-energy term is often neglected or absorbed into the matter Hamiltonian in the context of MQED, 16,18,28 recent works have shown that its inclusion is critical to ensure gauge invariance, and the stability of the coupled light-matter system in the strong-and ultra strong coupling regimes. 29,30When included explicitly, the dipole self-energy term is most commonly expanded in terms of the modes of the electromagnetic field. 17,25As discussed in Supplementary Note A, we follow the same procedure in this work to arrive at the dipole self-energy term in Eq. 10.In practice, most numerical implementations apply some truncation of the photonic Hilbert space.This is also true for our framework as discussed in the following sections.It has been shown that any truncation of the photonic Hilbert space must be accompanied by a consistent truncation of the dipole self-energy term to avoid unphysical predictions of both ground and excited state properties. 31In this work we therefore include the same modes in the expansion of the light-matter interaction term and the dipole self-energy term.
The interaction of the electronic system with the electric field within the dipole approximation can thus be expressed as the interaction of the electronic system and a continuous set of quantum harmonic oscillator modes via the dipole moment of the electronic system.We thus arrive at a Hamiltonian which is able to describe the full electronic structure of matter in the presence of a realistic electromagnetic environment.While the electronic system is described fully ab initio, it is important to note that within MQED, the environment is de-scribed via its spatially dependent dielectric properties.These can for example be calculated themselves using ab initio methodology for the materials making up the cavity structure, [32][33][34][35] or be described using simpler models of dielectric response such as e.g. the Drude model or the Lorentz Oscillator model. 24,36Note that our formulation addresses both the problem of how to formulate the length gauge Hamiltonian in the presence of optical losses, and allows for the explicit calculation of the cavity field strengths in terms of the boundary conditions set by the cavity via the DGF of the electromagnetic field.

Quantum Electrodynamical Density Functional Theory
QEDFT is a generalization of density functional theory (DFT) for electronic systems interacting strongly with quantized modes of the electromagnetic field. 17The method of QEDFT can describe the full electronic structure of the matter as well as the interaction of electrons with photons, and represents a good compromise between accuracy and computational cost.
QEDFT has been successfully applied to describe both the ground-state, 37,38 and excitedstates 26,39 of single (few) molecules strongly coupled to quantized modes of light, as well as for applications in polaritonic chemistry. 40The existence of the QEDFT formulation can be proven under very general conditions 17 and it in principle allows for the treatment of coupled light-matter systems with many electrons and many photonic modes under very general conditions.However, currently, most practical implementations of QEDFT are based on the dipole approximation for the light-matter coupling and a discrete mode expansion of the electromagnetic field.The latter implies that the material or molecular system of interest is embedded in a lossless electromagnetic medium.There have been previous studies applying the QEDFT formulation to lossy optical cavities to describe e.g.photon losses through cavity mirrors. 26,27,41Here, different models of the optical cavity were used, but no general connection between the cavity field parameters and the optical environment for absorbing and dispersing magnetoelectric bodies has yet been established.As a result, the electron-photon coupling parameters, while in principle connected to the physical quantity of the vacuum electric field at the center of charge of the system, are then in practice often treated as free parameters.This highlights the other current practical limitation of QEDFT: Treating the cavity coupling parameters as essentially free parameters make quantitative calculations and direct comparison with experiments hard.Here we overcome these practical limitations by representing the electromagnetic environment in terms of the MQED field expansion which allows us to accurately account for both optical losses and the quantitative magnitude of the cavity field strengths.The result is a parameter-free description of arbitrary electronic systems coupled to general, lossy electromagnetic environments.
To calculate the excited-state properties of the coupled electron-photon problem defined by the Hamiltonian in Eq. 10, we employ the linear-response formulation of QEDFT, which has been described in detail in Refs. 26,39,42Importantly, this method considers the full electronic structure of the matter system and goes beyond the rotating wave approximation.
We specifically employ the generalized Casida formulation of linear response QEDFT to calculate the oscillator strengths of the many-body excitations of the coupled light-matter (polaritonic) system.The generalized Casida equation reads, Here ω α is a diagonal matrix with the frequency of photon modes in the diagonal, U accounts for the coupling between the electrons and V the coupling between electrons and photons.
The eigenvalues Ω 2 S of the Casida equation are the square of the electron-photon excitation energies ω I .F and P describe the Casida eigenvectors of the matter and photon block, respectively.The square norm of the two gives the electronic-and photonic fraction of the excitation respectively.Introducing the pair orbital index S = (ia), corresponding to the pair orbital Φ S (r) = ϕ i (r)ϕ * a (r) with energy ϵ S = ϵ a − ϵ i , U can be expanded as, If N pair pair-orbitals are included in the calculation, U is therefore an N pair × N pair matrix.
V αS is the matrix accounting for the coupling between the electrons and the photon modes, The size of these matrices will be N modes × N pair .We note that these equations include the exchange correlation kernels f n M xc , f q M xc and g nα M which have to be approximated in practice.
So far, the available exchange-correlation functionals for QEDFT in the time-domain are still limited, 37,38 but recent developments for efficient density functionals based on the photon-free formulation of QEDFT 43 or the QEDFT fluctuation-dissipation theorem 44 are promising.
Finally, note that when neglecting the quantized light, only the top left block in Eq. 13 survives and one recovers the standard Casida formulation of linear response TDDFT for finite systems. 45 general, Eq. 13 requires a self-consistent solution as U and V depend on the eigenvalues Ω S .In this work, we neglect any exchange-correlation contribution in the photonic exchangecorrelation kernels and apply the mean-field photonic random phase approximation.We refer to Ref. 26 for a thorough discussion of this approximation.In this case, the connection between the photonic exchange-correlation kernels and the cavity field strengths in mode α We emphasize here the connection of the linearresponse exchange-correlation kernels of QEDFT with the dyadic Greens function of the MQED framework via the cavity field strength λ α defined in Eqs.7 and 11.Importantly, the Casida formulation works for a discrete set of cavity field strengths.In practice, it is therefore necessary to employ a dense, discretized sampling of the continuous frequency expressions for the coupling strengths.As discussed in Supplementary Note C, we employ uniform sampling in this work.7][48] While it is not immediately clear that these approaches can be integrated with the framework presented here because the quasi-normal modes are made interacting in these schemes, it would be an interesting avenue of future research to investigate whether such sampling methods can be used in our formalism.
As discussed in Ref. 26 it is possible to calculate the oscillator strengths of the coupled system in terms of the eigenvectors and eigenvalues of Eq. 13 and these fully describe the linear response of the coupled light-matter system.In this work, we characterize the response of the system using the oscillator strengths of the electronic polarizability, f I .
For the transition between the many body ground state Ψ 0 and the excited state Ψ I the oscillator strength is formally given in terms of the transition dipole matrix element as Importantly, we note that even in the presence of quantum light-matter interaction, these oscillator strengths still fulfill the f -sum rule I f I = N e , 26 where N e denotes the number of electrons.The excitation spectrum is then characterized by the strength function, 26 where the sum index I runs over the many-body excitations of the coupled light-matter (polaritonic) system.We further emphasize that in this work we do not apply any broadening to the spectra and that the width of the peaks is determined completely by the electromagnetic environment.This highlights a further advantage of using the QEDFT-MQED combination, namely that this provides a natural description of the transition linewidths as they relate to decay induced by the electromagnetic environment.This is relevant beyond cavity QED settings and paves the way towards TDDFT without artificial linewidths.
3 Results and discussion

Spherical microcavity
To exemplify the developed approach, we consider the spherically layered microcavity also considered in Refs. 49,50As shown in Fig. 1a, the spherical cavity consists of three concentric spherical layers, each characterized by a frequency dependent dielectric function, ϵ n (ω).
Due to spherical symmetry of the cavity setup, the DGF is most efficiently expanded onto a set of vector spherical harmonics. 19For a general emitter position inside the cavity, it is necessary to carefully converge the number of vector spherical harmonics used in the calculation of the DGF.However, as we discuss in Supplementary Note B, if the emitter is placed in the center of the cavity, the situation simplifies significantly.In this case only the lowest order transverse magnetic mode of the cavity contributes and we can write, 50 where r TM n=1 (ω) is the reflection coefficient for the lowest order transverse magnetic (TM) mode at the interface between the cavity region and the metal shell.r TM n=1 (ω) can be calculated by invoking the standard electromagnetic boundary conditions at the interfaces between the different regions of the cavity geometry.Notice that n i • ImG(0, 0, ω) • n j ∝ δ ij which means that the mode orthogonalization is trivial in this case, and the cavity field strengths can be derived directly using Eq.11,

Drude metal shell
We now consider the case where the inner and outer regions consist of vacuum.For the middle region we consider a simple but realistic model of a metallic mirror, namely a Drude metal with dielectric function, where ω p is the metal plasma frequency and γ is the Drude dampening rate.As a concrete example of a metal, we use the Drude parameters for gold taken from Ref. 51 Within the Drude mode, the plasma frequency of gold is around 8.5 eV1 .This results in the dielectric function shown in Figure .1b. Below the plasma frequency, the real part of the dielectric function will be negative and the material surface will consequently be highly reflective.Above the plasma frequency, the real part of the dielectric function becomes positive and the material will loose its metallic characteristics resulting in a significant loss of surface reflectivity.We note that because the V matrices in Eq. 13 grows with the number of photon modes, it is necessary in practice to truncate the electromagnetic environment at some frequency.For the spherical cavity, ω p provides a natural cutoff and we therefore include photon modes up to ω p in the calculations.Further details of the QEDFT calculations are given in Supplementary Note D.
In Figure 1c the mode structure of the cavity setup is shown for two different cavity radii, R = 140nm and R = 450nm.It is clearly observed that the number of modes, as well as their spectral position, is directly linked to the radius of the microcavity.Furthermore, we clearly observe that above the gold plasma frequency, the mirrors lose their reflectivity which results in the loss of the sharp mode structure which is replaced by a continuum.This highlights that the formalism we present is able to directly link the cavity field strengths to the real cavity setup made of real materials.Figure 1d zooms in on the mode around 7.1 eV in the cavity with R = 140nm and shows the effect of changing the Drude dampening parameter.We clearly observe that the width of the cavity mode increases with increasing dampening in the metal which highlights the connection between the width of the cavity modes and the losses in the gold.
This example highlights how the use of the emitter centered representation of MQED allows us to directly and uniquely relate the light-matter coupling strength to a real electromagnetic environment and connect it to the QEDFT formalism.While this is a relatively simple example, the approach is general and works analogously for an arbitrary electromagnetic environment provided that the DGF can be determined.We next add a benzene molecule to the cavity.Benzene is chosen mainly because of it prevalence as a test system in the existing TDDFT and QEDFT literature on strong coupling, 26,27,38,54 but we emphasize that the method can treat arbitrary finite electronic systems.We focus on finding a cavity configurations with a mode resonant with the Π → Π * transition of the benzene molecule.The first step is to determine the spectral position of this transition.Using the Casida linear response QEDFT framework without photons, we find that the transition occurs at 6.808 eV (182 nm) in free space and that the transition dipole moment in-plane is |d| = 0.096 e nm.This is consistent with previous TDDFT calculations for benzene. 26,27 shown in Fig. 2a, it is possible to find different radii of the gold microcavity for which there is a cavity mode resonant with the benzene Π → Π * electronic transition.As expected, the cavity field strength increases as the cavity is made smaller.All but the smallest cavity are optical cavities in the sense that the characteristic dimension of the cavity, the radius, is larger than half the wavelength of the transition.For the smallest cavity of radius 16 nm, a significant increase in the coupling strength relative to the other sizes is observed.
This happens exactly because this cavity is sub-wavelength sized and therefore significant near-field coupling to the surface plasmon mode of the gold starts to occur.We mention in passing that while we solve the coupled system using QEDFT, the approach presented here for the calculation of cavity field strengths is applicable regardless of the method used to solve the coupled system.For completeness, we therefore note that in terms of the lightmatter coupling strengths, g, more commonly used in quantum optics, the three cavities shown in Fig. 2a respectively correspond to hg = 3.81 eV, hg = 0.62 eV and hg = 0.21 eV for the Π → Π * transition of the benzene molecule.Importantly, these light-matter coupling strengths are a property of the coupled system and they are therefore specific to both the cavity and molecular transition under consideration.One should therefore be careful with general conclusions on the magnitude of light-matter interactions based on these numbers.
Fig. 2b shows the linear absorption spectra of the coupled emitter-cavity system calculated for the different cavity radii using the linear response QEDFT method.We emphasize that all linewidths in the figure are true linewidths in the sense that they are not related to any broadening parameters in the QEDFT calculation and only reflect the density of states in the optical environment.A radius-dependent Purcell enhancement with decreasing radius is clearly observed, reflecting the reduction in radiative lifetime resulting from the altered optical environment.However, it is not possible to achieve strong coupling with a single benzene molecule using the gold-shell cavity.We attribute this to the fact that as the coupling strength gets larger with decreasing radius, the optical losses also increase, resulting in a broader cavity resonance.We note in passing that the Purcell enhancement for the smallest R = 16 nm cavity is around 3500 which means that the local field enhancement at the center of the spherical microcavity is comparable to what is found in experiments with plasmonic microcavities. 55e reason that it was impossible to reach strong coupling with benzene in the gold cavity was the losses of the cavity mode.In an attempt to reach the SC regime, we therefore next seek to reduce the losses in the cavity.As already discussed above, the width of the cavity resonance is reduced for smaller Drude dampening parameters, γ.For this reason, Fig. 3a shows the absorption spectrum for the case with the true gold dampening, as well as 25%, 10% and 5% of the dampening respectively.We mention in passing that in practice one could imagine realizing these lower losses using for example metals specifically engineered to show weaker losses. 56Furthermore, one could also imagine exploring different cavity setups potentially leveraging the lower losses in dielectric nano-optical setups. 57Considering the wide range of available materials this design space becomes enormous. 35,58As shown in the inset of Fig. 3a, we find that reducing the losses results in a narrower cavity mode without a significant reduction in the overall cavity field strength.At around 25% of the Another way to engineer the coupling strength is by changing the emitter.There are two ways to do this, either by changing the number of emitters or by changing the emitter itself.To investigate the first option, we take the cavity with 25% of the true gold losses, and compute the absorption spectrum for one, two and three benzene molecules, all of which we place in the center of the cavity.0][61][62] To perform this analysis, one needs to solve a coupled many electron, many photon problem and it again highlights the strength of the method.
We first use the Casida method without photons (standard TDDFT) to characterize the spectral properties of the aromatic molecules.fails to describe the linewidth of the transitions.We observe that with increasing number of aromatic rings N , the transition energy redshifts and the transition dipole moment increases linearly (Fig. 4b).The linearly increasing transition dipole moment would suggest that the light-matter coupling strength can be monotonically increased by simply using a larger aromatic molecule.However, because the transition energy also redshifts with increasing molecule length, the cavity has to be re-optimized to be resonant with the transition for each molecule, as shown in Fig. 4c.Specifically, focusing on modifications of the R = 16 nm cavity, we find that this re-optimization of the cavity means that the cavity radius needs to be increased.This increase in radius leads to a reduced field concentration via an increased effective mode volume.This behavior highlights the important point that the light-matter coupling strength is a joint property of both the electronic system and the electromagnetic environment.A proper treatment of both is therefore essential for quantitative predictions.
We can characterize the intricate interplay between the electromagnetic environment and the electronic structure by looking at an effective coupling strength for the cavity modes which we define as, 27 Here the total cavity field strength parameter is defined as the coupling strength averaged across the cavity peak p, λ c = p dω|λ d (ω)| 2 , and ω c is the center frequency of the cavity mode.The d subscript indicates that we take the field strength parameter for dipole orientation d/|d|.g eff would thus be the true light-matter coupling if the total spectral weight was concentrated in a single mode.As shown in Fig. 4d, we observe that the increase in transition dipole moment is counteracted by the reduced field concentration for the larger molecules, effectively resulting in a weaker light-matter coupling strength.This is in stark contrast to the intuitive argument based solely on the increased dipole moment of the longer molecules.It should be noted that g eff is not a perfect measure of the light-matter coupling strength and it only gives a rough idea of the cavity re-optimization's effect.This is because the widths of the modes are not taken into account which, as we have seen in Fig. 3a, is very important for the nature of the light-matter coupling.
Importantly, the coupling strength does not reach a significant fraction of the transitions frequencies for any of the cases considered above.This highlights that reaching the ultra-and deep strong coupling with single or few molecules is in general hard.For the coupling strengths we find in this paper, the results from the combined QEDFT-MQED method would therefore agree with those found using the method presented in Ref. 15 With the presented framework it would be possible to perform further engineering of the electromagnetic environment to increase the coupling.The application of the framework to general electromagnetic environments is discussed further below.

Comment on Fabry Perot cavities
A common example of a cavity in the literature for both theory and experiments is the layered Fabry Perot cavity (FPC).As discussed in Supplementary Note F, the FPC is a layered system and consequently its DGF is expanded in terms of in-plane plane waves, augmented by a function accounting for the reflection at the interfaces between the layers. 19Because the FPC only constrains the electromagnetic modes in one direction, it retains significant dispersion of the modes in plane.Consequently, the concentration of electromagnetic density of states is significantly less efficient than in the case of e.g. the spherical microcavity.This means that the resulting coupling strength is weaker and the FPC will therefore generally not be suited for single-or few emitter strong coupling. 63For this reason, we do not perform explicit QEDFT calculations for this cavity setup.However, the FPC can still be suited for collective strong coupling 64 and coupling to extended systems 65 where the extended modes of the electromagnetic environment can be sampled more effectively.

Tool for cavity field parameters in simple cavities
As a part of this work, we are making the code to generate the cavity field strengths available for everyone to use as part of the new PhotonPilot tool.This tool currently allows the user to calculate cavity field strengths for the spherical and layered cavity setups, and we plan to expand its capabilities in the future.

General electromagnetic environments
We emphasize that the method we have presented here is general and applicable to any electromagnetic environment as long as the DGF can be determined.However, in general electromagnetic environments with lower symmetry it is not possible to write down an analytical expression for the DGF.In such cases, the DGF must be constructed numerically from e.g. a mode expansion based on finite element simulations. 19,24We note in passing that in the general setting the Helmholtz equation is not a Hermitian operator.Special care is therefore needed when constructing the spectral representation of the DGF from the modes.
One solution is to use the biorthonormal construction discussed in Ref. 66 We envision that the method presented in this paper will eventually be integrated with the existing Maxwell solver in the Octopus code. 67,68Such an integration would allow for the treatment of general electromagnetic environments completely within Octopus.

Conclusion
In this paper, we have presented a methodology combining macroscopic quantum electrodynamics with quantum-electrodynamical density-functional theory which provides a fully ab initio description of coupled quantum light-matter systems.Importantly, while we describe the electronic system at the DFT level, it is also possible to employ our approach to quantitatively describe the electromagnetic environment within standard few-level models of strongly coupled light-matter systems such as e.g. the Jaynes-Cummings model, the Rabi model, or the Travis-Cummings model.
To exemplify this approach, we have considered a benzene molecule strongly coupled to a metallic spherical cavity and investigated the impact of the both cavity radius and cavity loss on the nature of the light-matter coupling.We have further investigated the effect of adding more molecules and exchanging benzene with larger aromatic molecules.Together, these results highlight the intricate interplay between the electronic structure of the emitter and the environment in determining the nature of the light-matter coupling.Our work therefore illustrates the importance of having a proper description of both the electronic system and the electromagnetic environment for a proper description of quantum light-matter interactions.This work sets out the direction for more quantitative calculations in the future and also opens the possibility for the proper treatment of real experimental setups.We emphasize that the connection between the optical environment and the DGF is not limited to setups similar to cavities, but instead provides a general way to determine the electromagnetic spectral density of an arbitrary environment.In addition to the QED setup, our method therefore also provides a way to perform time-dependent density-functional theory in a lossy optical environment without the need for artificial spectral broadening.
Finally, we have provided an easy to use tool that everyone can use to generate cavity parameters for simple cavities such as the spherical microcavity or a layered cavity.

Figure 1 :
Figure 1: Spherical microcavity setup: (a) Illustration of the spherical microcavity setup.(b) The Drude model dielectric function of gold.(c) The modes of the spherical microcavity plotted for two different radii, R = 140nm in blue and R = 450nm in orange.The mode structure is shown with a sampling density of 10 points/meV.(d) The impact from the Drude dampening parameter on the cavity resonances, illustrated by focusing on the mode around 7.1 eV in the cavity with R = 140nm.

Figure 2 :
Figure 2: (a) Different radii of the spherical microcavity which host a resonance that aligns with the Π → Π * of the benzene molecule.The total light-matter coupling strength of the modes is also shown in the legend, and it can be clearly observed that the cavity-coupling strength grows with decreasing cavity radius.The mode structure is shown with a sampling density of 10 points/meV.(b) The linear absorption spectrum of benzene as a function of cavity radius showing a clear radius dependent Purcell enhancement.

Figure 3 :
Figure 3: (a) The linear absorption spectrum of benzene in the gold cavity from Fig. 2a with an inner radius R = 16nm as a function of the Drude dampening in the metallic mirror region.(b) The linear absorption spectrum in a cavity with with an inner radius R = 16nm and Drude dampening γ = γ gold /4 as a function of the number of benzene molecules.

Figure 4 :
Figure 4: (a) The Casida spectrum without photon modes for different aromatic compounds.The spectra here are shown with an artificial broadening of 0.1361 eV because the photonfree calculations fail to naturally describe the linewidth of the transitions.(b) Transition dipole moment of the Π → Π * -transition as a number of aromatic rings.(c) Mode structure of the smallest cavities which host resonances aligned with the Π → Π * -transition of the different molecules.The mode structure is shown with a sampling density of 10 points/meV.(d) Effective light-matter coupling strength as a function of the number of benzene rings in the aromatic molecule.

Fig. 4 (
a,b) shows respectively the spectrum and transition dipole moment of the Π → Π * -transition (HOMO-LUMO) as a function of the number of aromatic rings N .Note that in Fig. 4a the spectra are shown with the Octopus default artificial broadening of 0.1361 eV.This artificial broadening is necessary because