Dynamic of Single Molecules in Collective Light-Matter States from First Principles

The coherent interaction of a large collection of molecules with a common photonic mode results in strong light-matter coupling, a feature that proved highly beneficial for chemistry and termed the research topics polaritonic and QED chemistry. Considering complex microscopic chemical reactions in combination with a macroscopic number of molecules renders existing ab initio approaches inapplicable. In this work, I introduce a simple approach to capture the collective nature while retaining the full ab initio representation of single molecules. By embedding the majority of the molecular ensemble into the dyadic Green tensor, we obtain a computationally cheap and intuitive description of the dynamic of a single molecule in the ensemble - an approach that seems ideal for polaritonic chemistry. The introduced embedding radiation-reaction potential is thoroughly discussed, including prospects, applications and limitations. A first application demonstrates the linear response of single molecules that are part of a larger ensembles of molecules. Then, by virtue of a simple proton-tunneling model, I illustrate that the influence of collective strong coupling on chemical reactions features a nontrivial dependence on the number of emitters. Bridging classical electrodynamics, quantum optical descriptions and the ab initio description of realistic molecules, this work can serve as guiding light for future developments and investigations in the quickly growing fields of QED chemistry and QED material design.

However, most of the so far existing chemically relevant applications couple many molecules collectively to a common cavity mode. On one hand, a predictive theoretical description of chemical processes requires a thorough description of the electronic and vibrational * Electronic address: christian.schaefer.physics@gmail.com structure from first principles. On the other hand, a direct description of large ensembles is prohibited by the quickly increasing computational cost. The vast majority of theoretical descriptions are therefore based on simplified models that posses a limited applicability to complex systems and chemical reactions. Descriptions from first principles have been restricted to small numbers of molecules [47][48][49][50][51][52][53]. Notable representatives are quantum-electrodynamical density-functional theory (QEDFT) [54][55][56][57] and cavity coupled-cluster theory [58,59].
By embedding the majority of the molecular ensemble into a local potential, we will discover a path to open first-principles techniques towards collective strongcoupling with arbitrary numbers of ensemble molecules and species. This approach represents an extension of the radiation-reaction ansatz derived in reference [47] and inherits its computational and conceptual simplicity, allowing an almost effortless implementation into existing time-dependent density-functional theory (TDDFT) libraries.
We start by introducing the radiation-reaction ansatz (section II) that accounts for the fully self-consistent interaction between classical light and quantum matter in a simple and elegant fashion. Based on this construction and the experimental conditions in polaritonic chemistry, the majority of the ensemble will be embedded into the local potential (section II A). Solving the common Kohn-Sham equations [60][61][62] with the addition of this embedding radiation-reaction potential provides then access to the dynamic of single molecules in collective lightmatter states. Its application and physical consequences are demonstrated for electronic and vibrational strong coupling (section III).

II. EMBEDDING THE ENSEMBLE INTO THE ELECTROMAGNETIC ENVIRONMENT
In the non-relativistic limit, the Hamiltonian governing the electronic and nuclei dynamic under the influence of classical electromagnetic fields is given byĤ = i 1 2mi (−i ∇ i − q i A(r i t)/c) 2 +Ĥ +ε 0 /2 dr[E ⊥ (rt) 2 + c 2 B(rt) 2 ] with fixed Coulomb gauge ∇ · A = 0. In addition, the self-consistent interaction between light and matter dictates that the electromagnetic fields follow Maxwell's equation. As demonstrated in reference [47], the classical interaction with light can be efficiently incorporated into electronic structure theory via the transverse component of the dyadic Green tensor G ⊥ E r,⊥ (r, ω) = iµ 0 ω V dr G ⊥ (r, r , ω) · (−ej(r , ω)) which is the formal solution of Helmholtz's equation ∇ × −ω 2 µ 0 ε 0 ε r (rω) G(r, r , ω) = δ(r, r ) (1) and characterizes the electromagnetic environment. The microscopic paramagnetic current j serves as driving inhomogeneity, i.e., oscillating charges emit light. In addition, we have the freedom to assume that parts of the system behave as local and potentially isotropic linear medium r , µ r which shape via the Helmholtz equation the electromagnetic environment, represented by G. [63] In this way, G can subsume large fractions of a system while we free computational resources for the microscopic system described by j. The local radiation-reaction potential acts on the microscopic currents via the Schrödinger equation and manifests the self-consistent interaction with transverse fields within the long-wavelength approximation. This allows us to extend for example TDDFT with the self-consistent classical interaction between an electromagnetic environment and matter by simply adding v rr (rt) = er · E r,⊥ (t) to the usual Kohn-Sham equations. Since we anyway have access to the dipolar moment in TDDFT and dr(−ej(rt)) = ∂ t R(t), the radiation-reaction potential is trivial to implement into existing TDDFT libraries.
The longitudinal components of the electric field follow in dipolar approximation the equivalentV (t) = −R· E r, (t), E r, (r, ω) = iµ 0 ω V dr G (r, r , ω) · (−ej(r , ω)) system of equations. The longitudinal components are most relevant only in the near-field, i.e., they can be safely ignored when considering the cavity-mediated action between emitters that are separated by many nanometers. Our focus will rest in the following on the collective strong coupling mediated by transversal cavity fields. Nevertheless, solvation and other near-field effects such as the interaction with plasmonic systems are often critical in chemistry [10,12,64,65] and a wholesome picture should consider longitudinal and transversal fields on equal footing. We will focus in the following on electromagnetic environments that embody Fabry-Pérot-like eigenmodes and a large number of emitters. Clearly, G, and with it the presented embedding radiation-reaction ansatz, is generic and able to account for any kind of complex environment [47].
In the collective strong-coupling regime between many molecules and a cavity, single photon processes dominate and the individual effect of the cavity field acting on a single molecule is comparably weak. The electronic and nuclear excitation structure remains then largely unaffected, the combined light-matter system hybridizes based on those bare excitations. As long as this condition is satisfied it is safe to separate the matter-ensemble into two parts. The first describes a single explicit molecule featuring the paramagnetic current j, which will be obtained as solution to Schrödinger's equation. The second accounts for all the remaining molecules as polarizable material ε r (r, ω) = 1 + χ E (r, ω). Surely this separation is also convenient to account for two different species, e.g., the reactive molecule and its solvent. In this sense, the problem of chemical reactions in collective light-matter coupling can be seen as impurity problem, a single molecule undergoing a reaction is influenced by a large electromagnetic environment.
Knowing the bare Green tensor G 0 (empty cavity), we can obtain the 'ensemble-dressed' dyadic Green tensor with the help of the Dyson equation G(r, r , ω) = G 0 (r, r , ω) We assumed here (and earlier) that the generic susceptibility P(r, ω) = ε 0 dr 3 χ(r, r , ω) · E(r , ω) is local χ(r, r , ω) ≈ χ(r, ω)δ(r−r ). In many situations, the orientation of emitters within a gas or fluid will be random and its average polarizability is consequentially isotropic. [66] The concentration of the ensemble of emitters can vary such that the above equation would require for instance an iterative solution. We can, however, obtain an analytic solution in the special case that the ensemble occupies a small region in space compared to the wavelength's supported by G 0 . [67] The following sections demonstrate the simplicity of the embedding radiationreaction ansatz and illustrate its physical consequences for the collective strong-coupling of electronic and vibrational excitations.
A. The Radiation-Reaction Potential for Collective Strong Coupling For a dilute ensemble satisfying the long wavelength approximation, the response can be approximated as is evaluated at the position of the explicit molecule and V E is the volume occupied by the ensemble of emitters. The remaining task is to identify χ E (ω), in the following exemplified for TDDFT. Density-functional theory is designed to handle longitudinal perturbations, transversal perturbations imply the description of transverse currents which in turn implies the usage of current-density functional-theory [68][69][70][71]. In the long-wavelength approximation however, Maxwell's equations simplify tremendously and the coupling to longitudinal and transversal fields obtains the same formV ⊥, (t) = −R · E ⊥, (t). The (dipolar) polarization can be conveniently obtained from the densitydensity response function [72] as We can find a direct relation between polarization density P(ω) used in Maxwell's equation (χ(r, ω) ≈ χ(ω)) and the dipole R(ω) calculated from TDDFT R(ω) = α(ω) · δE local (ω) This approximation is reasonable if we focus on weak external fields that are consistent with a linear response treatment of the ensemble response. Furthermore, the long-wavelength approximation implies that the local fields acting on the microscopic current j and the ones acting on the ensemble polarizability α(ω) are comparable. In general, the local fields polarizing a single molecule consist of the sum of the external field plus the polarization-field of the surrounding material. Obtaining the self-consistent relation between micro-and macrofields is in general nontrivial, posing for instance notable challenges in solid state systems [73]. The Clausius-Mossotti relation [74] provides a simple and reliable solution for isotropic (α(ω) ≈ α(ω) = 1 3 tr(α(ω))) dielectric gases and fluids comprising of dilute molecules with negligible permanent dipole χ E (ω) = Under the assumption of a dilute molecular environment Obtaining the polarizability α(ω) is a standard task for any time-dependent electronic-structure code. However, it might be not always possible or desirable to obtain the embedding kernel from first principles and instead simplified approaches such as the Drude-Lorentz model present easily accessible alternatives. The following sections will demonstrate the usage of the simple model as well as the full ab initio description.
In combination with the radiation-reaction potential V rr (t) = −R · E r,⊥ (t) = −R · F −1 t (iµ 0 ωG ⊥ (ω)) * dr(−ej(rt)) , the feedback between the ensembledressed electromagnetic environment and the local molecule leads to collective states and adjusted local dynamics. The number of ensemble emitters does not influence the computational cost of this approach, the embedding scheme allows us therefore to reach the alleged macroscopic number of molecules participating in experiments.
The simplicity of the resulting equation becomes apparent from a simple example that will follow us throughout this work.
For this purpose, we chose a simplified one-dimensional Fabry-Pérot cavity The approximate analytic solution for such a transverse Green dyadic is easily obtained as We will add a weak cavity-decay rate η to account for a finite spectral width (see Appendix for numerical details). The poles of the Green tensor are then shifted to the positions of the collective ensemble+cavity hybridized states, routinely called polaritonic states, with hybridization energy Ω ∝ √ N E as detailed in the following. The necessary computational steps for a real-time calculation of the molecular dynamics in collective lightmatter states is illustrated in figure 1 and comprises: (1) Obtain ground-state density (potentially including QEDFT corrections). (2) Calculate or model χ E (ω) (e.g. via α(ω)). (3) Obtain dyadic Green tensor, either by solving for G 0 and dressing it according to equation (3)/(4) or directly solving equation (1) with ε r (r, ω). (4) Perform real-time TDDFT propagation as usual including the embedding radiation-reaction potential equation (2) and potential quantum-corrections from QEDFT as exemplified in App. B. The radiation-reaction potential is trivially combined with existing TDDFT libraries as the integrated current is given via the continuity equation as dr(−ej(rt)) = ∂ t R(t) and therefore naturally obtained during the propagation.
The additional steps 2-3 are performed only once before any real-time propagation, their computational cost Figure 1. Schematic illustration of the embedding radiationreaction approach exemplified for TDDFT. A single molecule described with time-dependent density-functional theory is coupled to an electromagnetic environment dressed by an ensemble of molecules. This paves the way to describe the dynamic of single molecules in large systems that are collectively coupled to electromagnetic environments from first principles.
remains therefore small and is dominated by obtaining α(ω).
Equipped with the necessary framework, we are now able to describe the dynamic of single molecules in collective strong-coupling. We start by demonstrating the emergence of collective states for small artificial but also realistic systems. We then shift our focus to chemistry and illustrate for a simple proton-tunneling model how chemical reactions can be described and modified via vibrational strong-coupling.

A. Illustration of the Single-Molecule Response
The moment we embed a macroscopic susceptibility into the dyadic Green tensor G ⊥ , the spectral response of the now dressed environment will change. If the excitation energies of χ E (ω) and G 0 overlap and their interaction is non-zero, the dressed G ⊥ will feature polaritonic excitations that are separated by the hybridization energy. Our individual molecule that hybridizes with those new quasi-particle excitations will be forced to oscillate at the excitation energies of the collective states. The larger the number of ensemble emitters N ensemble , the smaller its participation in the bright upper and lower polaritonic states. The spectral weight of the single molecules moves into an increasing number of dark states which are located at the original excitation energy. This bright/dark state behavior represents the hallmark of collective strong coupling and is at the heart of polaritonic chemistry [1,75]. By embedding the ensembledressed cavity into the electronic structure calculations, we obtain simultaneous access to the complex electronic dynamics and its interplay with macroscopic collective states.
A more complex bright-state dynamic appears if cavity, ensemble and explicit emitter are energetically detuned (see figure 3).
Cavity and ensemble will create the usual (large symmetric) splitting to which the detuned hydrogen contributes weakly. At larger N ensemble , the bright ensemble states start to hybridize with the bare excitation of hydrogen, leading to upper, middle and lower polaritonic bright states. The hybridization strength between middle and lower polariton is dominated by the fundamental coupling strength between hydrogen and cavity. Close to the avoided crossing, the contribution of hydrogen to the hybridized state can be substantially larger than the contribution of an average ensemble emitter -at this point hydrogen and the bright ensemble state contribute equally to the full polaritonic states. This becomes apparent when comparing figure 2 to figure 3. Similar effects have been also observed with quantum-electrodynamical density-functional theory [76]. Appendix A illustrates how the embedding radiation-reaction approach and the here discussed observations can be understood from the perspective of simple Photoabsorption cross-section σ(ω) of one-dimensional hydrogen coupled coupled with strength g0/ ωc = 0.0539 to a single cavity mode which is detuned at ωc = 11.7 eV with loss η = 10 −2 ωeg and dressed by N ensemble ensemble emitters. The ensemble response was represented by a Drude-Lorentz model quantum optical models.
An important limitation of this simplified approach rests in the self-consistency of the embedding procedure and the strength of any external perturbation. We made two assumptions, first, the ensemble can be described as linear medium, i.e., external fields are weak. While we could apply the radiation-reaction potential without issue for the full ensemble at once to treat e.g. highharmonic generation [47], i.e., ε r = 1, j = j i and G = G 0 , the resulting computational cost would be prohibitive for larger systems. The second, the cavity field exerts only weak forces on the ensemble. This limit is fulfilled in the weak and strong coupling regimes for light-matter interaction, i.e., cavity and ensemble can hybridize but the underlying matter-response experiences only minor changes. When reaching the ultra-strong coupling regime between light and matter, the underlying electronic or nuclear excitation structure will be modified considerably [7,42,56,59,77]. We can account for those change by e.g. using the QEDFT localdensity approximation proposed in [56] which is fully compatible with the here proposed radiation-reaction potential. We demonstrate such a combination in App. B and show that, for the here selected parameter regime, the QEDFT corrections are negligible. Furthermore, we can self-consistently update all components if extensions into stronger coupling regimes are envisioned. Imagine we couple a small set of molecules very strongly to the cavity but embed all but one into the Green tensor, substantial light-matter interaction will change the response of each molecule, which in turn affects the dressed Green tensor -a self-consistent cycle emerges in which χ, j affect G and the latter modifies χ, j. Such a generalized embedding scheme, which appears particularly appealing in the linear response regime, and its combination with the local-density approximation will be the subject of future work. The real-time propagation employed in the present work provides natural access out-of-equilibrium molecular dynamics and conveniently integrates cavity and ensemble lifetimes into the microscopic description.
Collective strong coupling experiments are safely within the strong coupling regime and typically performed in the 'dark', i.e., without or with only weak external drive [2,3,78]. The conditions of the illustrated approach are therefore satisfied. Most importantly, both limitations are only mildly limiting the dynamic of the explicitly treated molecule. It is entirely within the limits of the approximation if our explicit molecule undergoes for instance a charge reorganization or a chemical reaction -both considered to be of major interest in QED chemistry and QED material design.

B. Embedding Radiation-Reaction for Realistic Systems
Our aspiration is to describe the dynamic of realistic molecules in collective light-matter states. For this purpose, the embedding radiation-reaction ansatz has been implemented into GPAW [79] -a simple task as the local potential is easily combined with the TDDFT framework. This allows us to validate the previous conceptions for a realistic system, here chosen to be a chain of sodium dimers in H-aggregate configuration. The dimers as well as the cavity are oriented in z-direction, meaning that T c · α(ω) · c = α zz (ω). The ensemble dressed dyadic for this example reads then We follow now the in section II A illustrated flow. First, we obtain the bare response of a single sodium dimer using the radiation-reaction potential [47] with the cross sectional area A = 10 2Å . Next, we obtain the polarizability α zz (ω) which defines G ⊥ (ω) and with it the corresponding embedding radiation-reaction potential. Finally, a linear-response calculation of a single sodium dimer that is affected by the ensemble dressed cavity is performed. Figure 4 sets this embedding approach (red) in relation with the direct response of N sodium dimers (black) coupled to the cavity mode. Even when varying the number of sodium dimers and the cavity volume (g ∝ 1/ √ V where V = Aπc/ω c ), the embedding approach accurately recovers the direct calculations. The most notable difference is of course that the embedding approach describes the response of a single dimer, i.e., most of the spectral density is located in the dark state which is not visible in the direct N-dimers+cavity calculations. The polaritonic excitations in the embedding calculations show however the expected 1/N trend, e.g., (2 N a 2 )/(4 N a 2 ) ≈ 2, as well as the correct asymmetry between upper and lower polariton. Weak Coulomb mediated dipole-dipole interactions result in a slight blueshift, as expected from H-aggregates [47]. Increasing the distance between the dimers to 20Å alleviates this effect and we find excellent agreement between direct and embedding approach. As described earlier, Coulomb mediated couplings can be included in the longitudinal component of the embedding radiation-reaction ansatz.
Using the embedding radiation-reaction ansatz is substantially cheaper than calculating the full dimer chain explicitly. While the cost of the latter grows quickly beyond the manageable, our embedding approach retains the cost of a single-dimer calculation. This computational advantage opens up the possibility to investigate far more complex systems, incorporating the experimentally relevant realizations of collective strong coupling which will be the subject of future work. While we have limited ourselves thus far on the spectral properties of collective strong coupling, we will demonstrate in the following how chemistry itself is altered.

C. Modifying chemistry with collective strong coupling
We concluded previously that the participation of a single molecule to the collective state will usually decrease with increasing number of molecules, as expected from quantum-optical models. Nevertheless, a single molecules that exhibits a detuned excitation will contribute substantially stronger at the avoided crossing compared to a 'regular' molecule of the resonant ensemble (e.g. figure 3). During a chemical reaction, usually triggered by thermal fluctuations, structural changes of the molecule will alter its vibrational spectrum. In between reactant and product state, the molecule will spend time in intermediate configurations which will likely feature different vibrational energies and oscillator strengths. The closer this detuned vibration is to the collective polaritonic states, the stronger its contribution will out-weight the individual contribution of an ensemble molecule. While the reactant configuration will therefore play a negligible role in the polaritonic state, during the reaction our individual molecule might play a dominant role. Even if the number of dark-states might be large, the bright states can play a dominant role in this scenario. Quantitative estimates depend certainly on the specific characteristics of the molecule and its reaction.   Figure 4. Photoabsorption of N a2 chains coupled to a single cavity mode with ωc = 2.98 eV and ηc = 10 −2 ωc. The response is calculated either directly using TDDFT+radiationreaction for all N dimers (black) or using the embedding radiation-reaction (red) where only a single dimer is calculated explicitly. Different combinations of light-matter coupling strengths (g ∝ 1/ √ A) and particle numbers N are illustrated. Vertical gray-dashed lines serve as guide to the eye for specific polaritonic resonances. Note that the embedding approach describes the contribution of a single dimer to the collective state, i.e., the majority of the spectral weight is located in the dark states which is not visible in the direct calculations (as it is dark). Excitation energies and spectral weight are in excellent agreement otherwise. The direct calculations exhibit a weak blue-shift that can be removed by increasing the distance between dimers from 10Å to 20Å (magenta). Its origin lies in the Coulomb-mediated dipole-dipole interactions between the chain-elements. The excitations have been artificially broadened, numerical details in appendix D.
As we will see in the following, especially relevant is its vibrational structure and shifts thereof in intermediate configurations as well as reaction energetics, Stark-shifts and kinetics in relation to light-matter exchange times.
Let us illustrate this concept for a simple protontunneling model that is designed to represent groundstate reactions. A tilted double well potential v(x, t) = x·10 −3 −x 2 ·1.25·10 −3 +x 4 ·(1+  overcomes the barrier will be considered as an indication for the reactivity in our model. We start by setting cavity and ensemble excitation in resonance. Figure 6 shows the influence of collective light-matter coupling on the accumulated proton-tunneling. The influence of the ensemble on the reaction depends largely on the single-particle coupling strength to our impurity molecule and even for N ensemble = 0 is the reactivity already noticeably increased. Smaller coupling values (e.g. g 0 / ω c = 0.0027) show the tendency to decrease this effect with increasing ensemble size. However, even this comparably stringent trend features oscillations which can result in an inhibition of the reaction for larger N ensemble . Increasing the fundamental light-matter coupling does not only shift the resulting curve to smaller N ensemble values but furthermore intensifies the overall trend. In the extreme case of g 0 / ω c = 0.0135, we observe initially an increase in catalysis followed by an inhibiting effect for large N ensemble . Both, catalyzing [80] and inhibiting effect have been observed in experiments using either vibrational [19,20] or electronic strong coupling [9,12]. Figure 7 (top) illustrates that the influence of collective light-matter coupling is sensitive to the energetic alignment between ensemble, cavity mode and reacting system. The associated charge-transfer dynamic is shown in figure 7 (bottom). Such a complex behavior for our simple model demonstrates that the influence of collective strong coupling on chemical reactivity can be highly situational -a strong argument for the importance of ab initio approaches.
The dynamic that gives rise to this behavior emerges from the interplay of single-particle and many-body excitations in the correlated system and can be best understood by investigating the Green tensor of the ensembledressed cavity in figure 8 (here g 0 / ω c = 0.0135). Increasing Rabi-splitting (∝ √ N ensemble ) results in a shortening of the beating period observed in F −1 t iωG(ω). As The dependence of the cavity influence on the number of ensemble molecules N ensemble is non-trivial and oscillates between facilitating and inhibit the reaction. The cavity-influence on the chemical reactivity is far from a simple 1/N ensemble behavior. Especially intriguing, for small fundamental coupling strength, the inhibiting effect of collective strong coupling increases with increasing ensemble size. Up to N ensemble ≈ 10 4 , the absolute influence on the reaction remains comparably stable, followed by a slow decay when reaching the thermodynamic limit. Higher fundamental coupling strengths (g0/ ωc = 0.0135, 0.0108) feature even in the thermodynamic limit (N ensemble 10 8 ) clearly visible resonant effects. It should be noted, however, that the collective hybridization strength reaches substantial values in this case (see e.g. figure 12) and the underlying approximations are expected to limit the domain for reliable predictions to a combination of smaller fundamental couplings and (N ensemble 10 6 ). We use the in figure 2 introduced Drude-Lorentz model with ωp = 6.387 · 10 −4 eV . The first excitations are aligned at t=0 ωc = 0.046721 eV = ωE = ω0−1, the cavity decay rate is fixed to η = 10 −4 ωc.
we limited the reaction-time to approximately 1.2 ps, the qualitative influence of the ensemble will depend on how strong F −1 t iωG(ω) is noticeably altered on this timescale. This defines a minimal hybridization strength dxn(x, t))] − CR(g0 = 0) for increasing number of ensemble emitters for different energetic alignments between cavity/ensemble frequency and bare excitation frequency of our impurity molecule. Bottom, protontunneling The complex excitation-structure leads to a non-trivial influence on the proton-tunneling. In contrast to the other parameters, the smallest frequency tends to inhibit the proton-tunneling process. We fixed the cavity volume V in all calculations to the value of the resonant alignment. It is apparent from the behavior with increasing ensemble number N ensemble , that the specific effect of the cavity on the chemical reactivity is sensitive to the energetic alignment between cavity, ensemble and reacting molecule as well as the hybridization between those systems.
of G(ω) for a given reaction-speed 1/ω Rabi ∼ T Rabi < τ reaction -a quick reaction demands a large hybridization. How a chemical reaction will be influenced by collective strong coupling will therefore depend not only on the energetic alignment between reactant, product, solvent and cavity but also its hybridization strength in relation to the reaction time (see also reference [22]).
Quicker deformation of the potential-energy surface injects more energy into the system, leading to a larger reactivity and an enhanced cavity influence. Moreover, quicker ensemble decoherence leads to broad spectral response which tends to flatten the influence of collective strong coupling on the proton tunneling. Both effects are shown in appendix C.
An important realization is now that the mere necessity of the existence of collective polaritonic states defines a resonant condition. In other words, for the cavity to  have an effect dependent on N ensemble , collective strong coupling has to be present and an avoided crossing and energy-redistribution between shifted intermediate configuration and bare molecules has to be possible. Utilizing collective strong-coupling for chemistry along the here presented mechanism would therefore be limited to a subset of reactions with the correct intermediate-state dynamics. While experimental efforts strive to extend the realm of QED chemistry, among many experimental attempts only few reactions (all with similar molecular characteristics) have been successfully identified. Such a limited applicability of polaritonic chemistry would agree with the here suggested complex pre-requisites and calls for continued effort to identify further reactions in order to allow a conclusive understanding. This effort should consider the question to which extend for instance solvation can influence the number of effectively coupled molecules in collective states.
Identifying the specific influence of collective strong coupling on a chemical reaction remains an open and theoretically challenging problem. While recent work [22][23][24] showed promising progress, much remains to be investigated and a conclusive understanding will likely demand a collaborative effort involving experiment, ab initio theory and simplified models.

IV. CONCLUSION
By embedding an ensemble of molecules into the recently proposed radiation-reaction potential [47], we have been able to tackle the conundrum of a first-principles description of the material structure combined with the collective light-matter interaction involving a large number of emitters. The embedding ansatz is trivial to implement into existing time-dependent density-functional theory libraries as it adds a simple local potential. This allowed us to treat collective strong coupling featuring arbitrarily large number of emitters with marginal additional computational cost compared to ordinary realtime TDDFT calculations. We discussed a hierarchy of possible approximations and their associated limitations. Their utilization was first demonstrated for the linear response of a single molecule embedded in the electromagnetic environment of a cavity dressed by an ensemble of emitters. Integrating the embedding radiation-reaction potential into GPAW allowed us to tackle ensembles of realistic molecules from first principles. Lastly, by tuning the cavity in the infrared regime, we demonstrated that the modification of proton-tunneling reactions via collective strong coupling can posses a non-trivial dependence on the number of ensemble molecules.
This work paves the way to tackle new regimes of theoretical chemistry, accounting for the self-consistent interaction between light with macroscopic systems while focusing on the microscopic dynamic of single subsystems -entirely from first principles. Possible extensions include a consistent treatment of near-field effects that are relevant for solvation, plasmonic environments and J/H-aggregates. Equipped with this tool-set, describing consistently the seminal experimental work in polaritonic chemistry [1][2][3] moves within reach. This allows us to relax previous simplifications to only single or few molecules with artificially increased coupling strength [22,51,53,76]. Already our here illustrated simple proton-tunneling model exhibited a nontrivial behavior with the number of emitters. We can expect that complex multi-step chemical reactions will further complicate this trend. The seamless combination of embedding radiation-reaction and quantum-electrodynamical density-functional theory, especially the recently developed local-density approximation [56], ensures that quantum-corrections of the light-matter interaction can be properly accounted for if necessary. Collective strong coupling for chemistry and material design represent timely theoretical challenges to which the here presented embedding radiation-reaction potential provides an important and novel pathway. We illustrate here how the embedding procedure can be understood from the perspective of the quantumoptical Hopfield model for hybridizing modes. This widely used and strongly simplified model assumes the dipolar approximation, the rotating-wave approximation and emitters are approximated as bare 2-level systems. Expanding the corresponding HamiltonianĤ sH = ω câ †â into the excitation basis leads to a matrix equation that is commonly fitted in order to analyze experimental data.
The interaction between an ensemble of collectively interacting modes, describing a large set of ensemble molecules, a selected single molecule, and a cavity mode, will be given under the above approximations by the simplified Hamiltonian-matrix for which the N ensemble molecules posses the identical excitation energy ω E,1 = ω E,N and light-matter coupling g. In this case, arrow-head structure allows us to remove all those degenerate eigenvalues from the matrix by separating the problem into symmetric (also known as bright) and antisymmetric (also known as dark) sub-blocks. For example, the simplest case of this blockdiagonalization is given for 2 emitters coupled to one mode. The Hamiltonian matrix where the cavity couples now only to the symmetric/bright state with amplified strength √ 2g while the antisymmetric/dark state is decoupled from the cavity. The same can be achieved by noticing that for identical coupling and excitation energy the sum of excitation operators can be simplified to collective excitation operators S ± = N i σ ± i , also called Dicke operators. We reduce ourselves now to the bright-state submanifold described by with the excitation structure illustrated in figure 9. At the moment, the ensemble and the single molecule couple still to the bare mode. Let us assume ω c = ω E = 1 for simplicity. The upper-left 2-by-2 block possesses the eigenvalues ω ± = 1 ± g √ N and eigenvectors v T and the single molecule couples now to the collective bright states (upper and lower polariton ω ± ). The embedding procedure via the ensemble-dressed Green tensor is conceptually similar. The dressed G has poles at the bright ensemble+environment eigenvalues with adjusted strength, dark states are not present as they do not provide any spectral intensity from the perspective of the cavity. The microscopic material described via the Schrödinger equation is then affected by radiationreaction forces which are maximal at the polaritonic resonances. The here introduced embedding scheme can be seen as a generalization of the simplified Hopfield model to realistic systems. This generalization allows us to investigate time-dependent local phenomena for realistic systems, i.e., it provides access to chemical reactions, charge-transfer processes and many other features while retaining the full collective interaction. Clearly, the underlying assumptions that the majority of the ensemble is perfectly identical and responds only weak/linear are strong limitations to the embedding approach. More complex systems will demand a more sophisticated separation of the different embedded species, the specific path chosen here is designed to address the conditions for collective strong-coupling in polaritonic chemistry.
Appendix B: Ultra-strong coupling effectscombining radiation-reaction with QEDFT-LDA Every light-matter coupling strength greater than zero will alter the electronic and nuclear structure, reorganizing charge and breaking symmetries [7]. Typically, those effects are referred to as ultra-strong coupling features and their second-order scaling O(g 2 0 ), compared to the O(g 1 0 ) scaling of polariton hybridization, renders them of minor relevance for the vast majority of experimental realizations. Exceptions are located in the domain of plasmonics [35] and circuit-QED [81], the latter allows even coupling strengths on the order of the excitation energy g/ ω ≈ 1. Quantumelectrodynamical density-functional theory (QEDFT) allows to consider such quantum-corrections also for realistic systems. We employ here the recently proposed photon-exchange local-density approximation (pxLDA) [56] in combination with the radiation-reaction potential, an easy task as the latter accounts for the classical Maxwell-components and the pxLDA potential provides access to quantum-fluctuations. Figure 10 (top) exemplifies the renormalization of the electronic ground-state density for the in figure 2 performed calculation with N ensemble = 0. The effect is minute for the chosen couplings which transfers into a negligible spectral shift (figure 10 (bottom), red dashed) based on this change of the initial-state. Usage of the adiabatic pxLDA potential during the time-propagation (magenta dashed-dotted) provides a very weak blue-shift. Consequently, the pxLDA is overall of minor relevance in our calculations and we will omit quantum-effects in the rest of this work. Note however, that such a conclusion is system and parameter specific and adding pxLDA or comparable corrections should be considered for each realization individually.
Appendix C: Influence of loss and perturbation-speed Figure 11 shows the with figure 6 consistent calculations compared to calculations with quicker deformation (5 fs instead of 0.1 ps) or larger ensemble-loss (bottom). Quicker deformation of the potential deposits more energy in the system, resulting in overall larger reactivity and larger absolute cavity influence. Larger loss-rates smooth-out or flatten the cavity-influence but remain RR-potential RR-potential, n pxLDA (x, t = 0) RR-potential + TD-pxLDA Figure 10.
(Top) Ground-state electronic density renormalization ∆n(x) = n g 0 / ω=0.0563 (x) − n g 0 =0 (x) using the pxLDA potential [56] for the system investigated in figure 2 with N ensemble = 0. (Bottom) Photoabsorption cross-section σ(ω) excluding the adiabatic pxLDA functional (black solid), accounting only for the quantum-effects on the initial-state density (red dashed), and accounting for quantum-corrections also during the time-propagation with the help of the adiabatic pxLDA approximation (magenta dashed-dotted).
qualitatively similar to the values chosen in the main text. Figure 12 presents | G(ω)| consistent with figure 8.

Appendix D: Numerical details
The radiation-reaction potential and the according convolution have been implemented into a small python3 library and the open-source code GPAW. The timeevolution can be efficiently calculated by finite-difference approximations. The implementation has been checked for consistency against pre-existing Maxwell-QEDFT calculations [56].
Close to lossless electromagnetic environments feature sharp peaks in the spectral response, e.g. a onedimensional Fabry-Pérot cavity will be described by a set of harmonic oscillators. In order to allow a reliable numerical description of Fabry-Pérot cavities, we calcu- dxn(x, t))] − CR(g0 = 0) for increasing number of ensemble emitters. Top, CR for varying lightmatter coupling-strength between single perturbed molecule and cavity (the first excitations are aligned at t=0 ωc = 0.046721 eV = ωE = ω0−1, cavity decay rate fixed to η = 10 −4 ωc). The dependence of the cavity influence on the number of ensemble molecules N ensemble is non-trivial, i.e., does not simply decrease with 1/N as often believed. A critical feature is the appearance of beating signatures in G(t) (see figure 8). Shorter perturbation times (dashed, v(x, t) = x·10 −3 −x 2 ·1.25·10 −3 +x 4 ·(1+0.4/[1+e −(t−5f s.)/f s ])·10 −4 H) inject more energy into the system, leading to a larger reactivity and an enhanced cavity influence. Bottom, CR for different ensemble emitter-linewidths (solid, dashed-dotted, dotted) as well as for a shorter deformation period (dashed). Quicker ensemble decoherence leads to broad spectral response which tends to flatten the influence of collective strong coupling on the proton tunneling. lated the Green tensor via g 0 (ω) = 2 c 2 V F t e −ηt k ω −1 k sin(ω k t) .
An apparent but critical feature emerges when obtaining g(ω). The frequency resolution has to be sufficient to properly resolve g 0 and χ, both can become sharp for small losses. It is therefore important to converge the spectral representation of g(ω), for instance by artificially extending the time put into the construction of g(ω). A good indication of under-convergence for the Fabry-Pérot example is when the imaginary part of g 0 (ω) (Lorentzian shape) changes sign around a resonance. In this work, we increased the frequency resolution, compared to the natural resolution given by the propagation-time, by a factor of 2000 for all reactivity calculations. The Green tensor G(ω) and G(t) in figure 8 have been obtained with a 150 times higher resolution as they are visually identical from this point on and have only been used for visualization. Note, that this resolution is not sufficient to obtain an accurate description of the rather small and thus highly sensitive proton-tunneling process. Linear response calculations have been performed with a 10 times higher resolution as they feature already from the start a higher resolution. All GPAW calculations used a 1000 times higher resolution. The time-stepping remained unchanged such that the F −1 t G(ω) and the TDDFT propagation are consistent on the same time-interval. The additional computational cost can become noticeable if the resolution is pushed to extreme values as the FFT becomes then time-consuming. For most applications the cost will remain small or even negligible as FFTs are well optimized and G is constructed only once before the actual TDDFT calculation.
The subsequent real-time object F −1 t iωG(ω) can be contaminated with overtones close to t = 0 + , the explicit variant F −1 t iωG(ω) = −∂ t F −1 t G(ω) is preferred. Attention should be devoted to the different possible definitions of the Fourier transformation and their discrete counterparts. The macroscopic χ E is represented either from first-principles or via a Drude-Lorentz model with parameters as indicated in the captions. The presented cross sections in figure 2 and 3 feature artificial sign-flips towards lower and higher energies (seen from the resonances) as consequence of the finite time-interval. The artificial features were ignored here as they appear at small cross section value and do not affect the conclusions. Stronger damping or longer propagation would remove them if deemed necessary. Figure 2 and 3 have been calculated on a grid with 301 grid-points, 0.1 a 0 spacing and 4th-order finite-difference. The linear-response spectrum has been obtained by perturbing the electronic system with a Lorentzian-shaped delta-kick v(x, t) = −10 −4 1 π 10 −2 (t−1) 2 +10 −4 x (in a.u.) and time-propagation for 800001 time-steps with ∆t = 0.01 a.u. and 4th-order Runge-Kutta. Figure 4 was obtained in multiple steps as described in the main text. The polarizability of a single ensemble sodium dimer was calculated by real-time propagation via α zz (ω) = R z (ω)/(2πK) (K = 10 −5 being the kick-strength) using the LCAO routine of GPAW [79,82] including the radiation-reaction potential [47] with a cross-sectional area A = 10 2Å2 . The subsequently obtained spectra using the embedding radiationreaction calculations used the default real-time to spectrum tool of GPAW with a Lorentzian broadening of 0.05 eV (0.08 eV for 8N a 2 with 20Å distance). All calculations used a spacing of 0.3Å, a box of size 24+chain-length×24×24Å 3 , the atomic p-valence dz basis and the LDA functional. The sodium dimers have a fixed bonddistance of 1.104Å.
All reactivity calculations have been performed on a grid with 301 grid-points, 0.04 a 0 spacing and 4th-order finite-difference. The particle-mass was adjusted to the proton-mass. The real-time dynamics involved the deformation of the potential energy surface as specified in the main text with 100001 time-steps with ∆t = 0.5 a.u. and 4th-order Runge-Kutta.