Collective light–matter interactions have been used to control chemistry and energy transfer, yet accessible approaches that combine *ab initio* methodology with large many-body quantum optical systems are missing due to the fast increase in computational cost for explicit simulations. We introduce an accessible *ab initio* quantum embedding concept for many-body quantum optical systems that allows us to treat the collective coupling of molecular many-body systems effectively in the spirit of macroscopic quantum electrodynamics while keeping the rigor of *ab initio* quantum chemistry for the molecular structure. Our approach fully includes the quantum fluctuations of the polaritonic field and yet remains much simpler and more intuitive than complex embedding approaches such as dynamical mean-field theory. We illustrate the underlying assumptions by comparison to the Tavis–Cummings model. The intuitive application of the quantized embedding approach and its transparent limitations offer a practical framework for the field of *ab initio* polaritonic chemistry to describe collective effects in realistic molecular ensembles.

## I. INTRODUCTION

May it be for catalysis, solar energy harvesting, or superconductivity, controlling the state of a material on demand and the atomistic level defines a considerable area of modern science and technology. In addition to the intrinsic electronic and phononic degrees of freedom of the material, their interplay with the electromagnetic environment can play a critical role. Prominent examples include Floquet engineered materials^{1–3} and mode-selective chemistry.^{4} Coherently driving the material directly has two major downsides: (i) it requires energy, and (ii) it results in uncontrolled dissipation, which tends to obscure the path toward the desired target state. Electromagnetic resonators can reshape the continuous optical free-space spectrum into distinct optical modes, thus allowing a material inside to exchange energy only with specific modes. At sufficiently strong interaction, the resonant modes of matter and cavity hybridize and result in polaritonic quasiparticle excitations. Over the recent years, a plethora of changes in energy transfer^{5–16} and chemical reactivity^{17–24} have been observed and established the field of polaritonic chemistry.^{25–27}

The strength of hybridization between light and matter scales approximately as $NE/V$, with the number of collectively coupled emitters *N*_{E} and the effective quantization volume *V* of the confined field. Strong coupling is reached if the hybridization energy exceeds all losses and the two bright polaritonic excitations can be clearly separated from their background. Two prominent strategies exist to reach strong coupling: (i) minimizing *V* by using, e.g., subwavelength confined optical modes in (meta)plasmonic systems,^{28–31} and (ii) increasing the number of emitters *N*_{E} that interact coherently with the confined resonator mode.^{5–7,18,20,32} Polaritonic chemistry has primarily utilized the latter, which imposes, unfortunately, a high burden on the theory as the description of large ensembles becomes computationally overwhelming. As a consequence, the literature is dominated by Tavis–Cummings-like models, which replace the complex electronic structure with two levels and the resonator structure with a single harmonic oscillator. Effective (embedding) extensions of such models have recently gained attention^{33–35} as they promise to compress complexity into intuitive expressions. While simple quantum optical models have demonstrated remarkable success in describing the collective interaction, the goal to simultaneously capture changes in the chemical structure^{24,26,27,36} calls for new approaches that are able to combine sophisticated cavities^{37} with complex molecular restructuring during chemical reactions. Importantly, chemistry is foremost local while collective interaction is intrinsically delocalized, i.e., a local reformulation would greatly ease the theoretical challenge.

In this spirit, quantizing macroscopic Maxwell’s equations for all but a few molecules, or even a single molecule, is an attractive approach. Macroscopic quantum electrodynamics (QED) is capable of describing a complex optical environment seen by a single molecule via the linear response functions of the surrounding molecules and the cavity structure. Screening effects are accounted for via local field corrections, the impact of which on collective strong coupling has not yet been explored. Reference 38 introduced the Embedding radiation-reaction approach (Erra), which offers a simplified workflow compared to macroscopic QED, and demonstrated that, within mean-field theory, the collective interaction can be approximately absorbed into an effective light–matter interaction that is trivially compatible with local *ab initio* methods.

In this paper, we clarify the potential and limitations of quantized embedding approaches based on macroscopic QED for collective strong coupling by connecting them to standard quantum optical models. To do so, we introduce in Sec. II quantized embedding approaches and, specifically, in Sec. II B 2, the Quantized embedding radiation-reaction approach (Qerra). Section III demonstrates how Qerra recovers the Tavis–Cummings Hamiltonian in a suitable limit as well as collective effects (superradiance). For simple cavity structures, we find that Qerra offers a simple and straight-forward workflow (Sec. II C) that provides the quantized description of a macroscopic number of realistic molecules interacting with the cavity in the low-excitation regime. Our results in Sec. IV demonstrate how embedding approaches offer a computationally feasible way to go beyond common approximations employed in standard quantum optics approaches to describe the local dynamic of a single molecule that is collectively strongly coupled to a realistic azopyrrole in chloroform solution. Section V concludes our discussion and provides an outlook over the remaining challenges and possible solution strategies for polaritonic embedding approaches.

## II. QUANTIZED EMBEDDING APPROACHES FOR COLLECTIVE STRONG COUPLING

We first give an overview of field quantization in general absorbing and dispersive media via macroscopic quantum electrodynamics (Sec. II A). In Sec. II B, we then demonstrate how this framework can be employed to examine the collective strong coupling of a macroscopic number of complex emitters through embedding approaches. Section II C and Fig. 1 present a summary of the underlying workflow of these embedding approaches.

### A. Field quantization in absorbing environments

To find the electromagnetic field in the presence of linear macroscopic media,^{39,40} a fruitful approach is offered by macroscopic electrodynamics. Here, the macroscopic number of matter degrees of freedom is encompassed into a dielectric function, namely, the permittivity *ɛ*(**r**, *ω*).^{41}^{,} *ɛ* relates the linear polarization field generated by the media to the external electric field. Canonically quantizing classical macroscopic electrodynamics^{42} leads to the theoretical framework of macroscopic quantum electrodynamics,^{39,40} which is capable of finding the quantized electromagnetic field in general absorbing and dispersive optical environments given by their permittivity.

^{39,40}

*c*is the speed of light in vacuum,

*ɛ*

_{0}is the vacuum permittivity, and

**f**

^{(†)}(

**r**′,

*ω*) are bosonic creation and annihilation operators satisfying

**r**−

**r**′| →

*∞*. The magnetic field operator can be found from Eq. (1) via

**B**= i∇ ×

**E**/

*ω*.

**r**

_{i}and with dipole operators

**d**

_{i}, the corresponding light–matter Hamiltonian in the multipolar coupling scheme and the long-wavelength approximation reads

^{40,43}

*H*

_{mat}, including the dipole self-energy term,

^{44}and

*H*

_{F}is the Hamiltonian of the medium-assisted field,

*ɛ*(

**r**,

*ω*) for the given setup and, subsequently, determine the Green tensor $G$ by solving Eq. (4). $G$ and

*ɛ*then fully determine the quantized electromagnetic field. For example, the spectral density of the field reads

**r**.

^{45}Beyond a single emitter, the dynamics of a few emitters (possibly strongly) coupled to the medium-assisted field can analogously be found using macroscopic QED. However, finding the dynamics of macroscopically many emitters coupled to the four-dimensional continuum of field operators

**f**

^{†}(

**r**,

*ω*) poses great computational difficulties in general. Embedding approaches, as outlined in Sec. II B, suggest a suitable path to compress the complexity of many-body systems to an effective single- or few-particle problem.

### B. *Ab initio* embedding approaches

*N*

_{E}≫ 1 emitters (e.g., molecules) is placed inside a cavity. The basic idea of the embedding approaches is now to separate the system of

*N*

_{E}emitters and the cavity into microscopic components—which are solved directly—and macroscopic components—which are described effectively. The microscopic part is the part of interest, and in the following is considered to consist of just a single molecule (or a few) that can be described in detail by

*ab initio*methods. The macroscopic part consists of the remaining

*N*=

*N*

_{E}− 1 molecules and the cavity. Both the cavity mirrors and the remaining

*N*molecules are assumed to dress the electromagnetic environment of the microscopic part via linear response theory, i.e., via an effective linear susceptibility

**. Here,**

*χ***=**

*χ*

*χ*_{mir}+

*χ*_{mol}consists of a part describing the cavity mirror (

*χ*_{mir}) and the molecules (

*χ*_{mol}); compare Fig. 1.

*χ*_{mol}can be obtained by calculating the polarizability

**of a single molecule via**

*α**ab initio*methods and then employing a Clausius–Mossotti relation to find

*χ*_{mol}, e.g., in the dilute gas limit one finds

^{46}

*V*is the volume in which the emitters are located, and

*θ*

_{V}is the Heaviside theta function, i.e.,

*θ*

_{V}(

**r**) = 1 if

**r**∈

*V*and

*θ*

_{V}(

**r**) = 0 if

**r**∉

*V*.

Having absorbed most of the emitters into the environmental susceptibility *χ*_{mol}, we can use macroscopic QED as outlined in Sec. II A to quantize the effective electromagnetic environment $G$ for the given ** χ**, which now represents the cavity dressed by N emitters. What hampers the direct application of such an approach is the fact that the spectral density in Eq. (7) diverges in the coincidence limit (for

**r**=

**r**′). Physically, this feature is caused by the coarse graining of the medium, which neglects that the field seen locally by the emitter differs from the macroscopic one due to the discreteness of the individual emitters.

^{47}To account for these effects, local field corrections have to be employed, e.g., via the so-called real cavity model

^{48}as outlined in Sec. II B 1.

#### 1. Full macroscopic QED embedding

*R*

_{C}centered around the position of the emitter inside the surrounding media; see Fig. 1 (left). The Green tensor of the new geometry, including the small cavity surrounding the emitter, can be determined in the limit that the radius

*R*

_{C}satisfies

*R*

_{C}

*k*≪ 1, where

*k*=

*ɛ*(

*ω*)

*ω*/

*c*is the wave-vector of the field. In this limit, one finds that the Green tensor inside the small cavity reads

^{49}

^{,}

^{48–50}medium-assisted van der Waals interactions,

^{51}and strong-coupling of a single emitter in a cavity field with a background medium.

^{52}However, it is not clear how to connect the full macroscopic QED embedding approach based on the real cavity model outlined here to standard quantum optics frameworks, such as the Tavis–Cummings Hamiltonian. Moreover, in the full macroscopic QED embedding approach, the scattering Green tensor of the cavity and the emitters $G(1)$ must be calculated for each type of emitter and each emitter density separately, as each new configuration of emitters inside the cavity results in slightly different mode-structures. Both limitations can be lifted by introducing the streamlined embedding approach of Qerra.

#### 2. Qerra: Good cavity, dilute emitter, and small volume approximation

*V*

_{mic}that is much smaller than the resonant wavelength of the cavity. We further consider the good cavity and dilute emitter limits in which one can safely ignore the diverging bulk contribution to the Green tensor as well as local field corrections and approximate the bare Green tensor of the empty cavity $G\u0304\u2248G\u0304(1)$ by the scattering contribution $G\u0304(1)$. This neglects all direct interactions between the emitters that are not mediated by the cavity. In this limit, we can follow Ref. 38 to obtain the Green tensor inside the small volume

*V*

_{mic}accounting for the surrounding

*N*emitters and the cavity via the following Lippmann–Schwinger equation:

^{53,54}

**r**,

**r**′ ∈

*V*

_{mic}, and $G\u0304$ is the Green tensor of the empty cavity that we approximate by its scattering contribution, i.e., $G\u0304\u2248G\u0304(1)$. Assuming that $G\u0304$ is approximately constant over the spatial extension of the cloud of molecules

*V*

_{mic}, this equation can be solved via the dipole approximation, i.e., assuming that

*χ*_{mol}(

**r**″,

*ω*) =

*χ*_{mol}(

*ω*)

*δ*(

**r″**), such that

*ab initio*QED to macroscopically large values of

*N*. This scheme can further be easily extended to not just a single type of molecule but, e.g., also the scenario where a molecule of interest is placed in a solvent. To this end, the polarizability would consist of a sum of the polarizabilities of the molecules of interest and the solvent molecule weighted by their relative density. Note that Qerra provides solutions for any arbitrary number and mixture of emitters inside the resonator structure once the Green tensor of the bare resonator structure $G\u0304$ and the emitter polarizability

**have been determined. In contrast, in the full embedding approach accounting for local field corrections outlined in Sec. II B 1, the Green tensor of the cavity and the emitters must be calculated for every new concentration of the emitters.**

*α*So far, Qerra’s classical equivalent has been used to find the classical electromagnetic environment of a molecule inside a cloud of *N* ≫ 1 molecules that is placed inside a cavity, given by a bare Green tensor consisting of just a single Lorentzian,^{38} which physically corresponds to a perfectly single mode cavity. In the following, we use macroscopic QED to quantize the field in the embedding ansatz and show how this connects Qerra to standard quantum optics approaches in Sec. III. Before doing so, we summarize the workflow of the macroscopic QED embedding approaches for collective strong coupling in Sec. II C.

### C. Workflow of quantized embedding approaches

The workflow of the embedding ansatz is summarized in Fig. 1: In the first step, the polarizability of the individual molecules and the solvent is calculated in free-space from first principles using, e.g., time-dependent density-functional theory. The susceptibility *χ*_{mol} of the ensemble of *N* molecules is then found, for example, via a Clausius–Mossotti type relation [compare Eq. (8)]. The Green tensor $G(rNE,rNE,\omega )\u2261G(\omega )$ of the optical environment of the singled-out emitter at position $rNE$ can be obtained either via the full macroscopic QED embedding approach (see Sec. II B 1) or via Qerra (see Sec. II B 2). Once $G(\omega )$ and the susceptibility of the environment are determined, the resulting quantized field is given by macroscopic QED via Eq. (1). We have thus reduced the problem of a molecule interacting with a cavity of arbitrary geometry and material as well as a macroscopic number of other molecules to the interaction of a single molecule with a dressed multimode quantized field.

The resulting multimode macroscopic QED field can be utilized directly in *ab initio* calculations.^{55,56} Alternatively, one can further reduce the complexity of the four-dimensional continuum of field modes given by **f**(**r**, *ω*) to only a few dominant modes of the field using a suitable few-mode model.^{57–62}

The embedding approaches outlined here allow one to include realistic cavity models (described through a general Green tensor), go beyond the rotating-wave approximation, and account for the full molecular structure provided from first principles, yet the computational complexity remains constant for any number of molecules *N*. The approximations of these embedding approaches are summarized and compared to those used to derive the Tavis–Cummings Hamiltonian (see Sec. III) in Table I.

Framework . | Approximations . |
---|---|

Full MQED embedding | Linear response for embedded molecules, local field corrections (real cavity model) |

Qerra | Linear response for embedded molecules, good cavity, dilute emitter, small volume approximation |

Tavis–Cummings | Two-level approximations, rotating-wave approximation, single-mode approximation |

Framework . | Approximations . |
---|---|

Full MQED embedding | Linear response for embedded molecules, local field corrections (real cavity model) |

Qerra | Linear response for embedded molecules, good cavity, dilute emitter, small volume approximation |

Tavis–Cummings | Two-level approximations, rotating-wave approximation, single-mode approximation |

## III. COMPARISON TO THE TAVIS–CUMMINGS MODEL

Collective strong coupling is routinely studied via the Tavis–Cummings model and extensions thereof, such as the Holstein–Tavis–Cummings model.^{63,64} In such approaches, the electronic degrees of freedom of the emitters are approximated by two-level systems that are all located at the same position inside the cavity, and the cavity field is given by a lossless single harmonic mode.

*N*two-level systems resonantly coupled to a single-mode cavity and assume that the rotating-wave approximation applies. We consider only one polarization direction, such that we include only one diagonal element $G$ of the Green tensor. For a single mode field, it reads

^{37,65}

*ω*

_{c}is the frequency of the cavity,

*γ*

_{c}is its decay rate, and the constant

*f*

_{1}determines the resulting light–matter coupling strength; see below.

We can now follow different approaches: (i) start with a single mode approximation of the electric field resulting from Eq. (14) to find the dynamics via the resulting Tavis–Cummings Hamiltonian; (ii) use the Qerra embedding workflow from above, i.e., find the new, polaritonic multimode field for the cavity coupled to *N* molecules; and (iii) we can use the Tavis–Cummings Hamiltonian as in (i), but additionally apply the Holstein–Primakoff approximation, which is valid in the low-excitation limit. These three approaches are outlined and contrasted in the subsequent three sections.

### A. Exact diagonalization of the Tavis–Cummings model

*N*

_{E}two-level systems with transition frequency

*ω*

_{A}=

*ω*

_{c}to the field given by the Green tensor in Eq. (14) leads to the following Tavis–Cummings model:

^{54,57}

*d*being the dipole moment of one of the two-level systems. Here, we have neglected the dissipation of the cavity field by considering the limit

*γ*

_{c}→ 0. The Hamiltonian in Eq. (15) preserves the total number of excitations, such that

*H*

_{TC}can be decomposed into subspaces with fixed numbers of excitations. In the following, we focus on the single-excitation subspace of

*H*

_{TC}given by $HTC(1)$. $HTC(1)$ consists of the single excited state |1⟩ of the cavity; the bright state |

*B*⟩ ≡

*N*

^{−1/2}

*∑*

_{i}|

*e*

_{i}⟩, where |

*e*

_{i}⟩ is the excited state of the

*i*th emitter; the dark state manifold |

*D*⟩ that spans the subspace of the Hilbert space of the

*N*emitters with one excitation orthogonal to |

*B*⟩; and the excited state $|eNE\u232a$ of the singled-out emitter. As the dark states |

*D*⟩ are not coupled to the other degrees of freedom, we neglect them in the following. As outlined in the supplemental material of Ref. 38, the resulting single-excitation Hamiltonian can be expressed in the basis ${(|1\u232a+|B\u232a)/2,(|1\u232a\u2212|B\u232a)/2,|eNE\u232a}$ as

*N*

_{E}th emitter with coupling strength $g/2$. This is an exact solution of the model constrained to the single excitation manifold.

### B. Embedding approach

^{40}

*χ*

_{mol}from the polarization in Eq. (18). Inserting

*χ*

_{mol}and the single-mode Green tensor in Eq. (14) into Eq. (13), we find

^{66}

*ω*

_{c}=

*ω*

_{A}, this expression can be rewritten as

*ω*

_{±}. The coupling strengths

*g*

_{±}of these two field modes to the singled-out emitter can be obtained by comparing Eq. (20) with Eq. (14) as well as using Eq. (16),

*N*

_{E}coupled to the embedding environment, consisting of the cavity field mode and the

*N*other emitters, are governed by the effective Hamiltonian,

*N*of the emitters, the dynamics restricted to the single-excitation manifold obtained via the embedding approach and the exact diagonalization of the Tavis–Cummings Hamiltonian are exactly the same, as long as one does not enter the ultra-strong coupling regime (Ω

_{R}/

*ω*

_{c}≪ 1). Note that the difference here is the renormalization of the transition, which is present in the embedding approach but missing in the Tavis–Cummings Hamiltonian.

### C. Low-excitation approximation of the Tavis–Cummings model

We further simplify the Tavis–Cummings Hamiltonian in Eq. (15) using the low-excitation approximation. To this end, we employ the Holstein–Primakoff transformation.^{67}

*N*two-level systems, i.e., $Jz=\u2211j=1N\sigma z(j)$ and $J\xb1=\u2211j=1N\sigma \xb1(j)$ such that

*J*operators in terms of bosonic operators

*b*via

*J*

_{z}=

*b*

^{†}

*b*−

*N*/2 →

*b*

^{†}

*b*(in the last step we dropped a constant energy shift), $J+=b\u2020N\u2212b\u2020b$, and $J\u2212=N\u2212b\u2020bb$ we find

*b*

^{†}

*b*≪

*N*such that $N\u2212b\u2020b\u2248N$. We find

*N*two-level systems and the bare cavity mode by introducing new polaritonic, bosonic annihilation operators $f+=(a+b)/2$ and $f\u2212=(a\u2212b)/2$. We find

*N*if one is restricted to the single-excitation manifold considered in Sec. III A, i.e.,

*N*and any number of excitations, as long as Ω

_{R}/

*ω*

_{c}≪ 1. This correspondence shows that the embedding ansatz is valid whenever the low-excitation approximation is valid, i.e., for

*b*

^{†}

*b*≪

*N*. Most experiments in polaritonic chemistry are performed in this low-excitation limit, many of them even in the “dark” for

*N*≫ 10

^{6}, such that Qerra promises reasonably accurate predictions for as long as direct intermolecular interactions are of minor relevance.

Collectively coupled polaritonic systems at low pumping, i.e., in the low-excitation regime, represent, therefore, harmonic, i.e., classical, environments. Qerra extends this statement from simple models to the realms of *ab initio* QED by allowing for realistic response functions of the emitters and the cavity mirrors with only a little additional computational cost. The harmonic nature of the environment also validates the usage of macroscopic QED, as the latter is based on the linear optical response.

### D. Collectivity

We have used Qerra to single out one molecule acting as an “impurity,” which is embedded in the environment dressed by all other molecules via linear response theory. What is unclear so far is how much collectivity remains by following this procedure.

*N*

_{E}two-level systems, superradiance occurs when we prepare the two-level systems in a coherent superposition state $|+\u232a=(1/NE)\u2211i=1NE\sigma +i|{0}\u232a=(1/NE)(J++\sigma +(NE))|{0}\u232a$. When calculating the transition rate from | +, 0⟩ → |0, 1⟩ and from |1

_{j}, 0⟩ → |0, 1⟩, where |1

_{j}⟩ is the state where atom

*j*is in the excited state and all other atoms are in the ground state and |0⟩ and |1⟩ correspond to the zero and one-photon states with respect to

*a*, we find the transition matrix elements,

*H*

_{HP}in Eq. (31) that | + ⟩ couples with coupling strength $gN$ to the single field mode. Since the Holstein–Primakoff Hamiltonian is equivalent to the one from our embedding approach, it immediately follows that collective effects such as superradiance are also included in the latter approach. To show this explicitly, we have to express | + ⟩ in terms of the new eigenmodes defined by $fi(\u2020)$and $\sigma \xb1(NE)$. We find in the low-excitation limit,

*H*

_{Qerra}in Eq. (24) [which can be connected directly to the Holstein–Primakoff Hamiltonian

*H*

_{HP}in its form (32)], we get the same enhancement of $NE$: The single photon state can be expressed as $|1,0\u232a=(f+\u2020+f\u2212\u2020)/2|0,0\u232a$ and we find

*ω*

_{±}=

*ω*

_{c}±Ω

_{R}with $\Omega R=gN$, compare Eq. (21). This conclusion is similar to Ref. 33 and linear dispersion theory,

^{68}which has been shown to cover strong coupling

^{69}as well as collective effects.

## IV. APPLICATIONS

Having established a connection between the Qerra embedding approach and the Tavis–Cummings Hamiltonian, we next illustrate how Qerra can be used to go beyond the limitations of the Tavis–Cummings model with only a little extra computational effort. First, in Sec. IV A, we show how Qerra can be used to account for detuning of the two-level emitters and emitter losses. Then, we use Qerra to study the collective strong coupling of a more realistic molecular system in Sec. IV B. In Sec. IV C, we apply the full macroscopic QED embedding approach (Sec. II B 1) to a setup in which the emitters are not localized in the center of the cavity but are evenly distributed between the mirrors of a planar Fabry–Pérot resonator. We account for dispersion and absorption in the gold mirrors of the cavity and further include local field corrections. For all applications, we only consider a single polarization direction of the field and focus on analyzing the spectral density of the optical environment that is given in terms of the Green tensor via Eq. (7).

### A. Qerra: Two-level emitters

*ω*

_{c}and decay rate

*γ*

_{c}, and the

*N*two-level emitters with frequency

*ω*

_{A}, decay

*γ*

_{A}, and polarizability

*γ*

_{A}(Purcell effect) and the frequency

*ω*

_{A}(environment-induced frequency shift) are affected by the cavity and

*N*

_{E}emitters via the embedding approach. We see in Fig. 2(a) that for zero detuning Δ ≡

*ω*

_{c}−

*ω*

_{A}= 0, the single Lorentzian of the cavity (dashed black line) splits into two polaritonic peaks separated by the Rabi frequency Ω

_{R}. The slight difference in the height of the upper and lower polariton indicates that in general, the coupling strength of the singled-out emitter to the upper and lower polariton is not the same, i.e.,

*g*

_{+}≠

*g*

_{−}; compare Eq. (27). If the two-level systems are red or blue shifted from the cavity frequency, we find that the resonance frequencies of the polaritons are also red or blue shifted, respectively. Furthermore, the relative height of the two polaritonic peaks changes, meaning that a singled-out emitter in a blue shifted (red shifted) ensemble will couple more strongly to the lower polariton (upper polariton) [see blue and green line in Fig. 2(a)].

When considering the impact of the decay rate *γ*_{A} of the two-level systems on the spectral density, we find that with a decreasing decay rate of the emitter, the two polaritonic resonances in the spectral density become narrower. This shows that with decreasing decay rates of the emitters, also the decay rate of the polaritonic modes decreases, as expected.

### B. Qerra: Realistic molecular example using Qerra

Having established Qerra as a tool to study the collective strong coupling of two-level systems to a single-mode cavity, we next show how it can be used to analyze collective strong coupling of realistic molecular ensembles. To this end, we have calculated the polarizability of *ortho*-azopyrrole in *trans*- and *cis*-configuration (15% and 5% of the total number of coupled molecules) in chloroform solution (80% of coupled molecules) using time-dependent density-functional theory in combination with an implicit polarization model for the chloroform solution; see Subsection 2 of the Appendix for details. The isotropic average Im[*α*_{ave}(*ω*)] is shown in Fig. 2(e).

The cavity again consists of a single linearly polarized cavity mode that is tuned in resonance with one of the resonances of the azopyrrole solution at 227.8 nm; compare the black dashed line in Fig. 2(e). The resulting spectral density of the Qerra workflow is shown as the black line in the same figure. We can clearly see that the single resonance of the cavity splits into a lower and upper polariton. As before, the lower polariton has a lower height, resulting in a lower coupling to the singled-out emitter. The line shape of these two polaritonic modes differs from simple Lorentzians, indicating the necessity of multiple modes to fully describe them in an effective few-mode model.^{57} On top of these two polaritonic resonances, we find that the dressed spectral density also contains additional features corresponding to the other resonances of the azopyrrole solution in the vicinity of the cavity frequency. As all of these features only consist of a single peak per resonance, we find that none of the other resonances couple collectively strongly to the cavity mode (the decay rate is bigger than the collective coupling strength, resulting in no visible Rabi splitting).

### C. Full macroscopic QED embedding with realistic cavity and molecular system

*ɛ*

_{Au}of a Fabry–Pérot cavity. The singled-out emitter is assumed to be located in the center between the two mirrors. To treat this problem, we make use of the full macroscopic QED embedding approach, including local field corrections; see Sec. II B 1. We show the resulting spectral density for the singled-out emitter in Figs. 2(g) and 2(f); see Subsection 1 of the Appendix for numerical details. Local field corrections dominate the spectral density in this system, suggesting that Coulombic polarization-corrections routinely used in quantum chemistry for the description of solute–solvent dynamics (also used in Sec. IV B) provide a major contribution to the dynamic of the impurity molecule. In Fig. 2(g), we further distinguish the scattering part of the spectral density, given by

^{70}

^{,}

*J*

_{0}(

**r**,

*ω*) =

*J*(

**r**,

*ω*) −

*J*

_{scat}; compare Eqs. (9) and (7). Note that the bulk part

*J*

_{0}gives the spectral density in the absence of the cavity, such that

*J*

_{scat}represents all changes in the optical environment of the singled-out molecule induced by the cavity. We see that

*J*

_{scat}≪

*J*

_{0}, indicating that the single emitter “feels” the cavity only to a small extent, and its dynamics are instead dominated by its direct chemical and free-space environment. This is no surprise, as for Tavis–Cummings-like models of collectively coupled ensembles, each individual molecule contributes only to a minor extent and is also affected only to a minor extend. In such simple models, the effective coupling strengths acting on a single molecule have an upper bound defined by the bare cavity coupling, as is apparent from Eq. (25) and following (see also Ref. 38). Under which conditions the cavity is able to induce local changes remains a matter of active debate.

^{38,71–74}

Still, as shown in Fig. 2(h), the local-field corrected scattering part of the spectral density *J*_{scat} (38) shows the splitting of the cavity resonance at 227.8 nm into a lower and upper polariton. Furthermore, as the cavity resonances of the Fabry–Pérot cavity are very broad, the influence of other emitter resonances onto *J*_{scat} is clearly visible.

## V. CONCLUSIONS

A comprehensive description of collective strong coupling from first principles is a defining challenge in the theoretical description of polaritonics. Established models require a strong simplification of the molecular structure and even then become prohibitively expensive to solve numerically for increasing numbers of emitters.

Here, we discussed quantized embedding approaches for collective strong coupling that are based on macroscopic QED and in which the ensemble of emitters located in the cavity are split into an “impurity” and “environment” subsystem. All environmental degrees of freedom, including their (in)homogeneous broadening, are then used to dress the electromagnetic mode-structure encoded in the dyadic Green tensor $G[\chi ]$. Usage of the entire dyadic results in an embedding version of full macroscopic QED but necessitates local field corrections and the consistent treatment of free-space, direct longitudinal interactions between the emitters, and cavity induced interactions, which complicates its direct application. If the scattering contribution to $G$ dominates the mode-structure, i.e., working with high-Q cavities, and direct (mostly longitudinal) interactions can be ignored in the electromagnetic treatment (dilute limit), we can derive the simplified description Qerra. Qerra is easy to combine with *ab initio* quantum chemistry and offers, therefore, an intuitive approach to transition into the realm of *ab initio* QED. In particular, many of the recently developed quantum chemistry^{75–78} and non-adiabatic dynamic^{12,79–82} methodologies could be pushed closer to experimental reality by (partially) invoking quantum embedding approaches.

Following macroscopic QED allows then to quantize the dressed environment and couple it to the remaining impurity subsystem that is solved in full microscopic rigor. The impurity is coupled effectively to a set of quantized effective oscillator modes that represent the polaritonic eigenstates of the environment. We further illustrated this by demonstrating that Qerra is identical to the Tavis–Cummings model in the low-excitation limit and that collective (superradiant) effects are included in the single-excitation space. However, quantum embedding approaches for collective strong coupling can be readily extended to more complex environments, such as a cavity filled with a solute–solvent mixture that overlaps spectrally. We illustrated how such a complex environment affects the spectral density and how lossy cavities and emitters can be described with the full macroscopic QED embedding approach. Future work should focus on exploring possibilities to account for longitudinal inter-molecular interactions, thus moving Qerra closer to full macroscopic QED embedding, anharmonic corrections to the dynamic of the embedded ensemble, and full self-consistency between impurity and environment.

Quantized embedding approaches promise a theoretically rigorous description from first principles for large ensembles of molecules collectively coupled to electromagnetic resonators—paving a way to interweave the intuitive description in quantum optical models with the complexity of *ab initio* QED in a local framework suitable for QED chemistry.

## ACKNOWLEDGMENTS

We thank Andreas Buchleitner and Edoardo Carnio for their insightful discussions. This work was funded by the Spanish Ministry for Science and Innovation—Agencia Estatal de Investigación (AEI) through Grant No. EUR2023-143478. D.L. gratefully acknowledges the financial support from the Georg H. Endress Foundation and the DFG funded Research Training Group “Dynamics of Controlled Atomic and Molecular Systems” (Grant No. RTG 2717). C.S. acknowledges the funding from the Horizon Europe research and innovation program of the European Union under the Marie Skłodowska-Curie Grant Agreement No. 101065117.

This work was partially funded by the European Union. Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or REA. Neither the European Union nor the granting authority can be held responsible for them.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Frieder Lindel**: Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). **Dominik Lentrodt**: Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). **Stefan Yoshi Buhmann**: Supervision (equal); Writing – review & editing (equal). **Christian Schäfer**: Conceptualization (lead); Investigation (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: NUMERICAL DETAILS

#### 1. Green tensor of a gold Fabry–Pérot cavity filled with azopyrrole solution

*ɛ*

_{as}= 1 +

*χ*

_{mol}; compare Eq. (9). We here only consider one polarization component $Gxx(1)=G(1)$ of the field that is in-plane with the cavity mirrors. $G(1)(r,r\u2032,\omega )$in the coincidence limit in the center of the cavity is given by

^{40}

*L*is the distance between the two mirrors, $k=\epsilon as\omega /c$, $k\u22a5=k2\u2212k\Vert 2$ with Im[

*k*

_{⊥}] > 0, $D\sigma =1\u2212r\sigma 2e2ik\u22a5L$ account for multiple reflections between the mirrors, and

*r*

_{p}and

*r*

_{s}are the Fresnel reflection coefficients that read

*ɛ*

_{Au}is given by a Drude model,

*ω*

_{p}= 2.067 × 2

*π*PHz and

*γ*

_{Au}= 4.4491 × 2

*π*THz.

^{83}We further introduce polar coordinates

*k*

_{x}=

*k*

_{‖}sin(

*ϕ*) and

*k*

_{y}=

*k*

_{‖}cos(

*ϕ*) and perform the resulting

*ϕ*integral analytically. The remaining

*k*

_{‖}integral has been evaluated numerically.

#### 2. Polarizability for azopyrrole solution

All polarizabilities have been calculated using Casida time-dependent density-functional theory with the ORCA5.0 code.^{84} Structures for *trans*-azopyrrole, *cis*-azopyrrole, and chloroform have been relaxed using the CAM-B3LYP functional^{85} in a def2-TZVPD basis, including the implicit CPCM solvation model for chloroform.^{86}

^{87}

**d**

_{k0}and transition energy

*ℏω*

_{k0}between ground state and S-state

*k*. We include an artificial broadening of

*η*= 5 ⋅ 10

^{−3}Hartree.

## REFERENCES

*Ab initio*nonrelativistic quantum electrodynamics: Bridging quantum chemistry and quantum optics from weak to strong coupling

*Dispersion forces I: Macroscopic Quantum Electrodynamics and Ground-State Casimir, Casimir–Polder and Van Der Waals Forces*

For simplicity, we assume throughout the manuscript that the permeability *μ*(**r**, *ω*) of the media is given by *μ*(**r**, *ω*) ≈ 1. The embedding approaches discussed here, however, can also be used to study magnetically responding media.

*The Theory of Open Quantum Systems*

*Dispersion Forces II: Many-Body Effects, Excited Atoms, Finite Temperature and Quantum Friction*

*Ab initio*calculations of quantum light–matter interactions in general electromagnetic environments

*Ab initio*few-mode theory for quantum potential scattering problems

*ω*≫

*ω*

_{c}. This has a simple interpretation: Far from any resonances the two-level systems lead to a constant refractive index in the cavity, which shifts its resonances.

We only consider a single polarization direction in-plane with the cavity mirrors for simplicity.

*ab initio*modeling of polaritonic chemistry: The role of non-equilibrium effects and quantum collectivity

*ab initio*QED functional(s): Nonperturbative and photon-free effective frameworks for strong light–matter coupling

*ab initio*cavity-Born-Oppenheimer potential energy surfaces

*Time-dependent Density-Functional Theory: Concepts and Applications*