Cooperative quantum phenomena in light-matter platforms

Quantum cooperativity is evident in light-matter platforms where quantum emitter ensembles are interfaced with confined optical modes and are coupled via the ubiquitous electromagnetic quantum vacuum. Cooperative effects can find applications, among other areas, in topological quantum optics, in quantum metrology or in quantum information. This tutorial provides a set of theoretical tools to tackle the behavior responsible for the onset of cooperativity by extending open quantum system dynamics methods, such as the master equation and quantum Langevin equations, to electron-photon interactions in strongly coupled and correlated quantum emitter ensembles. The methods are illustrated on a wide range of current research topics such as the design of nanoscale coherent light sources, highly-reflective quantum metasurfaces or low intracavity power superradiant lasers. The analytical approaches are developed for ensembles of identical two-level quantum emitters and then extended to more complex systems where frequency disorder or vibronic couplings are taken into account. The relevance of the approach ranges from atoms in optical lattices to quantum dots or molecular systems in solid-state environments.

Some of the most intriguing phenomena in nature, both in the classical and quantum domains, are a product of cooperative effects, i.e., they cannot be understood by sole consideration of the individual constituents as they arise from the interplay among them. While at the fundamental science level, the understanding of many of those problems in the quantum domain already poses a great intellectual challenge, there is an ever increasing interest to build, control and harness complex cooperative platforms for emerging quantum technologies [1].
Light-matter platforms provide an optimal playground for the observation and exploitation of quantum cooperative effects. Quantum light, either multimode, as naturally arising in the quantum electromagnetic vacuum or single mode, as confined in the small volume of an optical resonator, can induce strong interactions among quantum emitters (QEs). Cooperativity then occurs e.g. in free space under high-density conditions and manifests itself in a strongly modified material response owed to a continuous scattering and rescattering of photons between the matter constituents. Two main aspects brought on by the common coupling of an ensemble of N quantum emitters to an electromagnetic environment are dipole-dipole interactions, stemming from a virtual exchange of photons and collective radiative emission, stemming from the loss of excitation into the infinite number of the electromagnetic vacuum modes. While the former is included as a coherent effect, the latter is an incoherent one, observable either as an increase (superradiance) or a decrease (subradiance) in the collective spontaneous emission rate compared to that of an isolated emitter.
Free-space near-field dipole-dipole interactions are utilized in structured subwavelength arrays of quantum emitters where they allow for the hopping of surface excitations. Such arrays are ideal platforms for achieving strong light-matter interactions and high fidelities for photon storage capabilities [2,3]. Furthermore, these platforms can allow for the propagation of photons or excitations to be protected against disorder and scattering caused by defects, in an approach dubbed topological quantum optics [4][5][6]. An advantage over linear topological photonic systems is the intrinsic nonlinearity of the quantum emitters which could lead to a rich many-body physics dynamics. Individual addressing of a single qubit emitter has been suggested via quantum spin lenses [7] and an approach to quantum networking with composite quantum systems comprised of many atomic arrays has been proposed building on theoretical results showing the production of a Bell entangled superposition quantum state for two distant arrays [8].
Cavity quantum electrodynamics (cavity QED) [19][20][21][22] can additionally mediate and enhance emitter-emitter interactions by interfacing quantum emitters with confined optical resonances. The famous model of Dicke superradiance [23,24], showing the quick burst of spontaneous emission from N indistiguishable quantum emitters can be simulated in the context of cavity QED and finds application in the operation of superradiant lasers [25,26]. The combination of subwavelength arrays with plain mirrors allows for the design of hybrid cavities with one or two frequency-dependent end mirrors, superior in performance to ones made up by plain, frequency-insensitive mirrors [27,28].
While many light-matter platforms are designed at the level of quasi-pure quantum emitters, i.e., considering only electron-photon interactions, a variety of extremely promising directions in quantum engineering are utilizing more complex quantum emitters such as organic molecules, quantum dots or color centers. In particular at the level of organic molecules, the strong confinement of light modes in micro-cavities gives rise to a novel research direction into vacuum-dressed or cavity-dressed materials, i.e., vacuum-hybridized materials with enhanced properties. At the level of mesoscopic systems, changes in charge conductivity [29][30][31][32][33], energy transfer rates [34,35], chemical reactivity [36][37][38][39][40] have been experimentally observed and theoretically studied. At the level of single molecules, the focus is in producing reliable single quantum emitters as single photon sources with applications in entanglement generation or in optical quantum computing [41]. Recent results have shown that a cavity-dressed molecule can act as an almost ideal quantum emitter, exhibiting a closed electronic transition [42]. For two or more closely spaced molecules, cooperativity is naturally manifested in the process of Förster resonance energy transfer based on energy exchange via near-field dipole-dipole interactions followed by quick vibrational relaxation [43,44], a key process in photosynthetic light harvesting.
Many of the aforementioned applications can be understood within the formalism of open quantum system dynamics extended to the more complex problem of correlated matter, such as occurring in a coupled quantum emitter ensemble. To this end, this tutorial utilizes two competing, interconnected approaches, one at the level of the density operator time evolution, i.e., the master equation (ME) and the other following time dynamics of system operators, i.e., the quantum Langevin equations (QLEs) approach. As the two formalisms are standard textbook methods [45][46][47][48][49][50], this tutorial proceeds with more complex aspects of light-matter interactions such as the emergence of a cooperative ME for coupled quantum-emitter systems in free space, introduced in Sec. II A together with the consequential occurrence of subradiance and superradiance, which are tackled in Sec. II B. Particular aspects such as the single excitation subspace and Dicke superradiance [23] are discussed in Sec. II C and Sec. II D, respectively. The master equation allows for the derivation of equations of motion for N mutually coupled emitters, utilized in understanding on how energy dispersion relations and energy band gaps in one-dimensional (1D) arrays can be engineered (in Sec. III A). Also in 1D chain and ring configurations, applications in quantum metrology, quantum information and lasing are presented in Sec. III B. The optical response of 2D subwavelength arrays is derived in Sec. III C to show the perfect reflection of incoming light around certain collective resonances.
Fundamental aspects of cavity QED with correlated quantum emitters are introduced in Sec. IV A followed by the input-output theory for operators detailed in Sec. IV B. The question of frequency disorder, relevant in the case of more complex quantum-emitter ensembles affected, for example, by inhomogeneous broadening, is tackled in Sec. IV C. Applications of the formalism introduced in Sec. IV are described in Sec. V A on 1D arrays towards antiresonance spectroscopy applications and in 2D highly reflective arrays employed in hybrid cavities exhibiting Fano-like narrow resonances. Finally, for three-dimensional mesoscopic ensembles of invertable two-level systems, the theory of Dicke superradiance is applied to the characterization of superradiant lasers inside lossy cavities.
The transition to increase the complexity of the considered quantum emitter systems is undertaken in Sec. VI A by providing a derivation for electron-vibron interactions in the form of a Holstein Hamiltonian [51], and emphasizing the application of QLEs at the level of polarons to characterize molecular emission and absorption [52,53]. The effect of cooperative behavior as mediated by dipole-dipole interactions is tackled in Sec. VI B to analytically characterize the FRET (Förster resonance energy transfer) process as well as for the description of the physics of molecular dimers [54], a standard model used in quantum chemistry.

II. COOPERATIVITY OF LIGHT AND MATTER
We provide in this section a quick review of collective effects in an ensemble of free space standing N quantum emitters with emphasis on the dipole-dipole interactions and collective decay mediated by the electromagnetic vacuum. The emergent properties of subradiance, superradiance and dipole-dipole-induced collective energy shifts are then introduced under the weak excitation assumption (weak driving). Eventually, dynamics in the full extended Hilbert space of the N emitters is treated by introducing the collective Bloch sphere and the resulting Dicke superradiance regime. This formalism is the starting point for the understanding of applications exemplified in Sec. III and Sec. V.
A. Collective radiative emission The standard system we consider in the next sections is that of an ensemble of N identical quantum emitters modeled for simplicity as hydrogenlike atoms with nuclei fixed in random positions R j , where j = 1, . . . , N . The (single) electronic degree of freedom for each emitter is described by quantized momentump j and positionr j (relative to its respective nucleus) operators appearing in the Hamiltonian h j =p 2 j /(2µ) + V (r j ) consisting of the kinetic energy (where the reduced mass is denoted by µ) and an electrostatic potential V (r j ). Out of the infinite set of eigenvectors of h j , we pick the lowest energy one (and set its energy to zero) |g j ⟩ and assume that all the physics described in the following involves transitions to only one excited state |e j ⟩ at energy ω 0 (we set ℏ = 1). This amounts to a two-level system approximation where the unity of the Hilbert space is a sum of only two projectors 1 j = |g j ⟩ ⟨g j | + |e j ⟩ ⟨e j |. To quantify transitions between the two levels, we introduce the standard ladder (Pauli) operators σ j = |g j ⟩ ⟨e j |, σ † j = |e j ⟩ ⟨g j |. Their commutator [σ † j , σ j ] = σ j z is the population difference operator. The free Hamiltonian can then be written as h j = ω 0 σ † j σ j or alternatively as h j = ω 0 σ j z /2 (up to a constant energy shift).
To describe the interaction of the emitters with the electromagnetic vacuum, one introduces a fictitious perfectly reflecting box of volume V = ℓ 3 and follows a standard quantization procedure for the electromagnetic field imposing periodic boundary conditions [45][46][47][48][49][50]55]. This leads to a plane-wave expansion of the electric field operatorÊ where the allowed k vectors are multiples of 2π/ℓ on each Cartesian direction and the index k runs over all possible k vectors and also over the two orthogonal polarizations with unit vectors ϵ k . The bosonic operators follow the commutation [a k , a † k ′ ] = δ kk ′ and their action is to create and destroy excitations in a given field mode. The zero-point electric field amplitude for a given mode with frequency ω k = ck is defined as E k = ω k /(2ϵ 0 V) (where c is the speed of light and ϵ 0 is the vacuum permittivity). The Hamiltonian for the field inside the box can then be written as a sum over an infinite number of bosonic modes H vac = k ω k (a † k a k + 1/2). The coupling of light and matter occurs within the formalism of the minimal coupling Hamiltonian [45][46][47][48][49][50]55] describing the total quantum system of charges and electromagnetic vacuum modes. Notice that the canonical momentum of the electron is shifted top j = µṙ j +eÂ(R j +r j ) whereÂ is the vector potential operator of the electromagnetic field. The observation that the size of the electronic orbital is much smaller than the wavelength λ 0 = 2πc/ω 0 associated with the optical transition between the two electronic orbitals allows for a simplified form of the Hamiltonian In this so-called dipole approximation, the light-matter interaction Hamiltonian H int = N j=1d j ·Ê(R j ) involves only the dipole moment operatord j = −er j and the electric field operator evaluated at the position of the nuclei. Under the two-level assumption, the dipole operator is written asd j = d eg σ j + h.c., where the transition dipole moment is computed between the two states d eg = ⟨g j |d j |e j ⟩ (assuming identical emitters). Moreover, one performs an additional rotating-wave approximation (RWA) where fast oscillating terms σ j a k and σ † j a † k (oscillating under free evolution with ω 0 + ω k ) are neglected. Under these approximations, the light-matter interaction part of the Hamiltonian can be written as FIG. 1: Schematics of the procedure followed to obtain the open-system dynamics formulation for a system of N quantum emitters (denoted by σj and σ j ′ ) subject to the quantum electromagnetic vacuum. The key point is that, for a pair of emitters, coupling to the common reservoir modes (denoted by a k and a k ′ ) gives rise to separation-dependent dipole-dipole interactions at rate Ω jj ′ as well as to collective decay rates γ jj ′ .
describing photon-emitter energy exchanges at rates g k = E k ϵ k · d eg . Notice that this interaction Hamiltonian conserves the number of excitations in the system. The Hamiltonian H governs unitary evolution in an infinite-dimensional Hilbert space (owing to the infinite number of electromagnetic modes). One could then, in principle, use a deterministic Schrödinger equation or equivalently the von Neumann equation for the density operator ρ tot (t) to compute the state of the system at any time. However, this is an extremely complex task requiring an immense computational resource even for small size systems. Instead, in order to drastically reduce complexity, an open-system dynamics approach is usually followed, which consists in reducing the system of interest to the Hilbert space of dimension 2 N of the matter part only (as schematically illustrated in Fig. 1). A main observation allowing such reduction of complexity is that the time-dependent light-matter interaction Hamiltonian can be simply expressed in the interaction picture as describing an excitation exchange between the emitters and the field, with the time-dependent operator acting solely in the Hilbert space of the photon modes. The open-system dynamics procedure is then based on a weak-coupling assumption and implies that the electromagnetic degrees of freedom are traced over to obtain an equation of motion for ρ(t) = Tr em [ρ tot (t)]. The steps for this procedure are summarized in Appendix A1 closely following the standard textbook derivation in Ref. [45] but also widely covered in, among others, Refs. [46][47][48][49][50]. In the final master equation all system properties are derived from two-time correlations of the operators F j (t). For N quantum emitters in the electromagnetic vacuum at zero temperature, the resulting collective quantum master equation readṡ The last term in the equation above describes irreversible loss of excitation into the electromagnetic vacuum and is expressed in the form of a superoperator (an operator acting on density operators) defined as [56] (7) Notice that this is not in standard, diagonal Lindblad form [45][46][47][48][49][50] which, for a collapse operator O and collapse rate γ O , is defined as (8) and describes decay at rate γ O through a single channel with operator O. The decay rates are γ jj ′ = (3γ/2)F (k 0 R jj ′ ) where the relative distance vector is R jj ′ = |R j − R j ′ | and the function of position is defined as with unit vector e d in the direction of d eg . The rate γ = ω 3 0 d 2 eg /(6πc 3 ϵ 0 ) is the spontaneous emission rate of a single emitter (note that the decay rate of the excitedstate population is given by 2γ). The collective radiative decay properties are given by the oscillatory behavior of F (kR) as illustrated in Fig. 2b. Notice that, for emitters separated by much more than λ 0 , the decay becomes purely diagonal as expected for noninteracting, independent emitters.
The coherent term H dd in Eq. (6) instead describes a dipole-dipole interaction characterized by a virtual exchange of excitation via the vacuum modes without loss of photons The dipole-dipole exchange rate Ω jj ′ = −(3γ/2)G(k 0 R jj ′ ) can be obtained from the same function that characterizes the collective decay rates via the definition where P denotes the Cauchy principal value. The explicit functional dependence of F (kR) and G(kR) on distance is detailed in Appendix A2. The behavior illustrated in Fig. 2b shows that the dipole-dipole interaction ceases at large distances, as expected, but diverges at close separations. This is, however, only an artefact of the initial assumptions that the dipole-electric field interaction is valid at any interparticle distance. This is of course not true, as for separations on the order of the size of the orbitals, one ends up with a fundamentally quantum many-body problem where the tunneling of electrons between neighboring emitters has to be taken into account Dipole-dipole coupling Ω12 and collective decay rate γ12 scaling with normalized interatomic distance k0a for parallel (∥) and perpendicular (⊥) orientation of dipoles (with respect to the connection vector R12). (c) Possible experimental realization of strong dipole-dipole couplings for two identical chromophores connected by an insulating bridge of length much smaller than λ0. The green arrows represent the orientation of the transition dipoles. Adapted from Ref. [54].
(leading to molecule formation, hybridization of orbitals, etc). More involved models, based, for example, on a quantum electrodynamics density-functional formalism, can then be employed [57].

B. Superradiance and subradiance
The superoperator in Eq. (7) describes nontrivial mutual decay characterized by the matrix Γ made up by the elements γ jj ′ . However, this term is not in standard diagonal form (characterized by a single collapse operator) as it is not comprised of N independent decay channels. One can perform a basis transformation with a matrix T (such that T −1 = T ⊤ ) which diagonalizes Γ such that diag (γ 1 , ...,γ N ) = T ⊤ ΓT whereγ k is the kth eigenvalue of the decay matrix. In order to achieve the Lindblad form defined in Eq. (8) one can define a set of collapse operators Π k = j T jk σ j such that now describes N independent decay channels each with an associated collapse operator Π k and associated loss rateγ k . Notice that the preparation of a collective quantum superposition that decays at one of the ratesγ k would require the application of the Hermitian conjugate of the corresponding collapse operator Π † k to the collective ground state |G⟩ = |g⟩ 1 ⊗ |g⟩ 2 ⊗ ... ⊗ |g⟩ N of the system. In practice this is, however, not a straightforward task.
Two coupled emitters -Let us first make the connection between collective radiative effects and the symmetry of quantum superpositions by considering the simplest example of two quantum emitters separated by distance a. Diagonalization of the decay matrix leads to superradiant/subradiant decay channels characterized by γ ± = γ ±γ and at the same time renders the Hamiltonian in diagonal form with eigenenergies ω ± = ω 0 ± Ω. The eigenstates are symmetric/antisymmetric superpositions |±⟩ = (|e⟩ 1 ⊗ |g⟩ 2 ± |g⟩ 1 ⊗ |e⟩ 2 )/ √ 2. This is illustrated in Fig. 2a in the collective basis where the other two states are the fully excited one |E⟩ = |e⟩ 1 ⊗ |e⟩ 2 and the ground state |G⟩ = |g⟩ 1 ⊗ |g⟩ 2 . The splitting between the two levels in the collective basis 2Ω = 2Ω 12 (a) as well as the magnitude and the sign of the mutual decay rateγ = γ 12 (a) depend strongly on the particular choice of dipole orientations as well as on separation (as illustrated in Fig. 2b). However, for distances below half a wavelength a < λ 0 /2, the antisymmetric (symmetric) states always have a subradiant (superradiant) character. One can then understand the connection between state symmetry and radiative properties in terms of an destructive (constructive) interference of radiative paths.
The strong scaling of the near-field dipole-dipole coupling with the interparticle distance renders such a simple system valuable for experimental applications in superresolution imaging. For example, experimentally it has been shown that optical resolution of fluorescent molecules can be achieved at distances as small as 12 nm [58]. In addition, for strongly coupled emitters, chemical or mechanical means can be employed to correct energy shifts and render the closely spaced emitters indistinguishable such that a source of indistinguishable photons can be achieved. Finally, we remark that subwavelength separations (even at the level of less than 10 nm) could be achieved in assembled molecular dimers [54] where two chromophores are coupled via an insulating bridge (illustrated in Fig. 2c). The connection between the dipole-dipole exchange for two-level systems and a vibronic dimer model for chromophores is detailed in Sec. VI B.

C. The single excitation subspace
While the simultaneous diagonalization of the Lindblad term and the Hamiltonian is generally not possible for N > 2, one can still get some intuition in the nature of cooperative decay in the particular case of N equally spaced quantum emitters in a 1D chain configuration. We start by analyzing the eigenstates and eigenvalues of the dipole-dipole Hamiltonian under the nearest-neighbor approximation and only up to a single excitation shared in the whole system (reducing the Hilbert space from 2 N to only contain N + 1 states). We then follow with an analytical derivation of the dissipation rates of such collective states. As a consequence of the fact that the Lindblad term (i.e., the decay matrix Γ) and the Hamiltonian do not commute, dissipation is not diagonal in the collective basis.
Eigenstates of the dipole-dipole Hamiltonian -Let us first inspect the Hamiltonian in the analytically solvable case where a nearest-neighbor (NN) approximation is performed. This is justified as the dipole-dipole interactions scales as R −3 jj ′ for distances smaller than a wavelength and thus the nearest neighbor coupling is almost an order of magnitude larger than that of the next to nearest neighbor one. The resulting Hamiltonian is then in the form of a tridiagonal Toeplitz matrix where Ω is the coupling between two neighbors. The free Hamiltonian has eigenenergies of degeneracy C N n = N !/[(N − n)!n!] for a given level of energy nω 0 where n ranges from 0 for the ground state to N for the highest excited state. We illustrate this in Fig. 3a where we plot the eigenvalues of the Hamiltonian up to the second excitation manifold. The reduced single-excitation manifold on which we focus is characterized by the particle basis |j⟩ = σ † j |G⟩ with j = 1, . . . , N and the collective ground state |G⟩. The resulting eigenenergies of the Hamiltonian are for an index k running from 1 to N (and the trivial energy of the ground state is 0). The corresponding eigenstates of the Hamiltonian are then described in the collective exciton basis Notice that the reverse transformation is straightforward |j⟩ = k f jk| k⟩. The transformation also holds at the level of operatorsσ k = j f jk σ j (in the single excitation manifold) which sees the diagonal form of the Hamiltonian as H = k ϵ kσ † kσ k .
Collective dissipation -To compute the decay rates of the Hamiltonian eigenstates, we use the master equation to derive the equation of motion for the population component ρ kk = ⟨k| ρ |k⟩. We arrive aṫ which shows that the eigenstates of the Hamiltonian couple within the first excitation manifold in addition to the decay to the ground state. The diagonal elements can then be estimated by setting k = k ′ to derivẽ γ k = jj ′ γ jj ′ f jk f j ′ k and more explicitlỹ From this expression, one can derive a scaling of subradiant states with roughly N −3 (see lower panel of Fig. 3b for the scaling of the most subradiant state) which will be utilized for cavity antiresonance spectroscopy applications in Sec. V A. More generally, scalings with N −s up to s = 5 can be reached [12]. Also, superradiance emerges, which for small distances a is characterized by a rate proportional to N but eventually saturates with increasing N (around N = λ 0 /a, as illustrated in Fig. 3b). Notice that generally, for k ̸ = k ′ , the dissipative couplings jj ′ γ jj ′ f jk f j ′ k ′ are not vanishing, which is the direct consequence of the fact that the Lindblad term and Hamiltonian do not commute. Up to here we have only considered the dissipative properties of the eigenstates of the Hamiltonian. Instead, one can directly inspect the Lindblad term by a diagonalization of the decay matrix in Eq. (12). Numerical results in Fig. 3b indicate that in this case subradiant states are present which are characterized by decay rates scaling exponentially down with N . As these states are not eigenstates of the Hamiltonian, they are not easily addressable. The targeting scheme would then involve some excitation scheme that applies the collective operator Π † k to the ground state |G⟩. In Sec. III B we will show how a combination of a magnetic field gradient together with a pulsed-laser excitation could instead be used to perform such an action.
In Sec. III A we provide an alternative treatment in terms of excitations propagating on the 1D support by hopping between neighboring sites via the dipole-dipole exchange. The discrete index k then describes the quasimomentum of an individual collective mode and its location with respect to the light cone distinguishes between superradiant versus subradiant modes.
Entanglement properties -The collective eigenstates of the dipole-dipole Hamiltonian (both super-and subradiant ones) commonly feature non-classical correlations [59,60], rendering them as an interesting resource for quantum information processing; highly subradiant states are of course even more useful due to the increased lifetime of correlations. To this end we analyze the logarithmic negativity [61], which is an entanglement monotone. For a bipartite system consisting of the subsystems A and B, it is defined as where ρ T A denotes the partial transpose with respect to the subsystem A and | · | is the tracenorm. In Fig. 3c, we show the time evolution of the state with k = N . At distinct time points, we compute the logarithmic negativity for each emitter (i.e., we choose our bipartite system to consist of the ith emitter and the rest of the chain). One can see, that the amount of bipartite entanglement is significantly larger in the center of the chain, even in the initial state. Over time, this behavior is retained and correlation is only slowly lost due to excitation loss of the chain. Even at t = 100γ −1 there still is considerable entanglement in the system. Let us now go beyond the single-excitation manifold and consider a famous example introduced by Dicke [23] which shows the generation of a superradiant pulse from an ensemble of indistinguishable quantum emitters. The model assumes an idealized case of N quantum emitters within a very small volume and neglects their dipoledipole interactions. While this is per se an unrealistic assumption (a densely packed ensemble of QEs exhibits large dipole-dipole frequency shifts), we indicate later in Sec. V C how this model can be realized in the context of cavity QED and is relevant for the physics of lasing in what is known as bad cavity superradiant lasers [25].
In the following we will make use of the Bloch-sphere illustration in Fig. 4a, where N two-level systems can be described by a collective spin of length N /2. This follows from 1/2 ⊗ 1/2 ⊗ ... ⊗ 1/2 = 0 ⊕ 1 ⊕ ... ⊕ N /2 for N even and 1/2 ⊗ 1/2 ⊗ ... ⊗ 1/2 = 1/2 ⊕ 3/2 ⊕ ... ⊕ N /2 for N odd (the tensor product space of dimensions 2 N of N spins 1/2 can be written as a direct sum of spaces of different sizes). We can then introduce collective spin operators S z = j σ j z /2 and S = j σ j and the total spin S 2 = S 2 z +(S † S +SS † )/2. As in the standard angular momentum algebra, it is then possible to find a collective basis, known as the Dicke basis of S 2 and S z , indexed by two quantum numbers |s, m s ⟩ In a closed system interacting with the electromagnetic field, the interactions are mediated by S and S † . Thereby the selection rules for optical transitions are given by ∆m s = ±1, ∆s = 0 and the Hilbert space splits into noninteracting subspaces defined by the quantum number s with dimension 2s + 1 (illustrated in Fig. 4a).
Let us now describe the phenomenon of Dicke superradiance in the ideal case where all mutual decay rates are equal and maximal γ jj ′ = γ. A trivial observation is now that the Lindblad decay term in Eq. (7) assumes a very simple form with a single superradiant decay channel and collective collapse operator We will follow the time evolution under such a Lindblad term for an initially fully excited ensemble characterized by the state vector |N /2, N /2⟩. Notice that the action of the collapse operator cannot take the system out of the s = N /2 symmetric manifold containing N + 1 Dicke states. We can derive an equation of motion for the population difference operator average ⟨Ṡ z ⟩ = −2γ ⟨S † S⟩ which simply states that the loss rate is proportional to the emitted intensity. For a single quantum emitter one has ⟨σ † σ⟩ = ⟨σ z ⟩ /2 + 1/2, implying that the spontaneous emission always follows an exponential law. The population difference of a collective Dicke state |s, m s ⟩ decays instead at a state-dependent rate 2γ ⟨S † S⟩ = 2γ(s + m s )(s − m s + 1). For the initially fully excited state, the decay rate is 2γN which is the same as expected for N independently decaying emitters. However, particle-particle correlations start building up during the evolution and by the time the |N /2, 0⟩ state (for even N ) is reached, the emission is superradiant and scales approximately as γN 2 /2 (see Fig. 4b). Notice that the crucial effect of correlations can be distinguished by Using that the sum over populations is given by ⟨ i σ † i σ i ⟩ = s + m s , the dipole-dipole correlation between emitters i and j can be estimated as ⟨σ † i σ j ⟩ = (s 2 − m 2 s )/(N (N − 1)), reaching a maximum of approximately 1/4 for m s = 0 and becoming zero for m s = ±s.

III. SUBWAVELENGTH QUANTUM EMITTER ARRAYS
The interplay between dipole-dipole interactions and collective radiance in quantum-emitter ensembles leads to a multitude of applications of 1D and 2D subwavelength arrays such as in nonlinear quantum optics [18,27], nanooptomechanics [15], the design of quantum metamaterials with magnetic response at optical frequencies [62,63], as platforms for quantum information processing [7,8,64] or as chiral light-matter interfaces [65]. These applications are based on the fact that such structures can support collective surface resonances that can interact in a controllable fashion with impinging fields.
We will start by studying the dispersion relations of the surface modes on 1D platforms by means of a Bloch ansatz, showing the occurrence of band gaps and Dirac points. This direction has recently emerged showing the promise of subwavelength arrays for topological quantum optics implementations [4,5]. We then proceed by providing illustrations on 1D emitter chains and rings aimed at showing the usefulness of subradiance as a resource for: i) improved frequency sensitivity (as shown in Refs. [10,66]), ii) robust quantum memories [59] and iii) the design of nanoscale coherent light sources as recently proposed in Ref. [67]. Finally, we derive the optical response of two dimensional subwavelength arrays around certain confined surface-mode resonances [18,68,69] and describe a regime recently experimentally tackled showing close to unity reflectivity for arrays of optically trapped atoms [13].
A. Band structure and topology of 1D chains Up to this point, we have solved for collective resonances and their associated collective radiation rates by a direct diagonalization of the Hamiltonian with near-field coupling terms. One can however take a solid-state approach instead where the array provides a crystalline structure for the quasiexcitations propagating on its surface. As the simplest example, let us first revisit the one-dimensional equidistant chain of N emitters with identical frequencies as already discussed in Sec. II C described by the Hamiltonian We first assume nearest neighbor (NN) coupling Ω jj ′ = Ωδ j,j ′ ±1 which insures translational symmetry with periodicity a (equivalently stated, the unit cell contains one site only). We impose periodic boundary conditions (PBCs) to simulate the mesoscopic case by asking that at the edge Ω 1N = Ω. We then ask the question: what kind of excitations can propagate in this chain, i.e., what kind of dispersion relations such excitations will exhibit. To this end, we start with the equations of motion for the expectation values β j = ⟨σ j ⟩ which can be obtained from ⟨σ j ⟩ = Tr[σ jρ ] whereρ = −i[H, ρ] (neglecting radiative emission in a first step) We furthermore assume a weak excitation limit ⟨σ j z ⟩ ≈ −1. In a compact matrix formulation, we rewritev = M v where v = (β 1 , · · · , β N ) ⊤ and the drift matrix is expressed as a circulant Toeplitz matrix The direct diagonalization of the matrix above gives the expected N collective modes listed in Sec. II C with a state index from 1 to N . However, we now instead look for a dispersion relationship with a quasimomentum q, to which end we plug the ansatz β j = A q e i(qaj−ωt) into Eq. (23) which straightforwardly leads to the following dispersion relation From the application of the PBCs we have β 1 = β N +1 which leads to qaN = 2πm where m is an integer. This indicates that the allowed quasimomenta are q = 2πm/(N a) and we can fix the first Brillouin zone to m ∈ {−N /2, N /2} corresponding to q from −π/a (left propagating wave) to π/a (right propagating wave). The resulting dispersion relations for two arrangements of dipoles (parallel and perpendicular to the chain axis) are indicated in Fig. 5a.
Beyond NN with dissipation -The results presented above describe a simplified model where only nearestneighbor interactions are employed and spontaneous emission is disregarded. We now proceed by analyzing the full equations of motion including dissipation which are obtained from Eq. (6) and reaḋ These equations take into account the exact behavior of Ω jj ′ , γ jj ′ (as defined in Sec. II A and App. A2) and lead to the following dispersion relation The result is numerically illustrated in Fig. 5b and compared with the approximation analytically obtained for NN coupling in Fig. 5a. The corresponding collective decay is plotted in the lower panel of Fig. 5b. The vertical lines indicate the location of the light cone, where the quasimomentum of the excitation propagating on the surface equals the wave vector of the photon k 0 = 2π/λ 0 , which for the chosen distance a = λ 0 /5 leads to q 0 a/π = 0.4. Notice that the interpretation in terms of the location with respect to the light cone is straightforward as waves with quasimomentum larger (in absolute value) than that of the photon cannot escape the chain and therefore are subradiant (in agreement with results presented in Ref. [3]).
Emerging band gaps -Let us now consider a slightly more complex structure with alternating emitter frequencies ω 1 and ω 2 and identical NN couplings at rate Ω (as illustrated in Fig. 5c). One then immediately notices that this implies a lattice with a double unit cell. We denote the two types of emitters by an upper index such that their amplitudes are (β j ) and the site index j runs from 1 to N /2 (such that we keep the chain length at N a).
The dipole-dipole interaction now couples the two kinds of emitters with each other within the unit cell and also between neighboring unit cells leading to the following set of coupled equationṡ We now ask for propagating waves with the following profile β and plug in this ansatz into the above equations to result in the following eigenvalue problem Two types of eigenvalues ω ± (q) result from diagonalization of the above matrix corresponding to two distinct energy bands The particularity of the system is that the two bands exhibit a band gap |ω 1 − ω 2 | at the edge of the first Brillouin zone (as depicted in Fig. 5c). The PBC now impose that β N /2+1 which leads to the condition q = πm/(N a). The first Brillouin zone is now defined by m ∈ {−N /2, N /2} corresponding to q varying between −π/(2a) to π/(2a).
Dirac points -Further complexity in the band structure can be engineering by assuming alternating coupling strengths (for example from Ω to −Ω for consecutive pairs as in Fig. 5d). Noticing that this situation is again characterized by a double unit cell, we follow the same steps as above to ask for two kind of propagating waves and reach the following eigenvalue problem which leads to the following dispersion relation for the two branches The quasimomentum varies as above between −π/(2a) to π/(2a). The resulting dispersion curve plotted in Fig. 5d shows the emergence of an avoided crossing in the center of the Brillouin zone which exhibits a Dirac point with linear dispersion relation at frequency degeneracy Berry phase -A particularly interesting case occurs for identical frequency emitters (transition frequency ω 0 ) with alternating interaction strengths Ω 1 and Ω 2 . This can be mapped onto the Su-Schrieffer-Heeger model (SSH) [70], which for electrons describes the emergence of bulk insulating phases distinguished by topological invariants and which in nature occurs in polyacetylene molecules. As before we write the equations of motioṅ and testing for the two kinds of propagating waves results in the following eigenvalue problem Here we will focus both on the energy dispersion curves but also on the characteristics of their corresponding eigenstates. The eigenvalues are given by While the band structure is insensitive to the exchange of Ω 1 and Ω 2 , the eigenvectors are not. The normalized eigenstates are analytically expressed in vector form as where the q-dependent phase is given by the following relation The eigenvalue problem can now be written formally as H(q) |φ n (q)⟩ = ω n (q) |φ n (q)⟩ (where n stands for ±). We now proceed by assuming that a path q(t) is taken, with q(t = 0) = −π/(2a) and q(τ ) = π/(2a); this is illustrated in Fig. 6b, particularized to the lower energy band. By moving adiabatically slow, tunneling to the orthogonal eigenstate is not allowed. Generally, we can then write the state of the system at any time t as |ψ(t)⟩ = e −iθ(t) |φ n (q(t))⟩ where the time-dependent phase could be computed directly by the application of the time-ordered evolution operator with a time-dependent Hamiltonian to the initial state. However, a more elegant solution comes from simply writing the Schrödinger equation H(q(t)) |ψ(t)⟩ = i∂ t |ψ(t)⟩ explicitly, to arrive at an equation for the phase We sandwich the equation above with ⟨φ n (q(t))| and integrate to obtain two distinct contributions. The integral of the energy over the whole band vanishes as the final and initial point are equal in energy. The second contribution is path dependent and reads The integration over time can be turned into a path integration such that This describes the path dependent Berry phase where A n (q) = ⟨φ n (q)|∂ q |φ n (q)⟩ can be identified as the Berry potential. The Berry phase is gauge invariant for closed integration loops. In general, the temporal aspect of the adiabatic deformation is a more or less fictitious process but it is helpful to unravel the topological structure. Integrating the Berry phase for our Hamiltonian H(q) over the Brillouin zone (−π/(2a), π/(2a)) for the lower band results in where the matrix element is easily computed from the vector expression of the eigenstates in Eq. (36) to result in A − (q) = (1/2)dϕ(q)/dq.
For differentiation we use the following expression for ϕ(q) = arctan2(Ω 1 + Ω 2 cos(2qa), Ω 2 sin(2qa)) and we note that ν is similar to a winding number. The result of the integral is a topological invariant showing that the two cases Ω 1 > Ω 2 and Ω 1 < Ω 2 with the same band structure outcome are topologically different which is expressed by the Berry phase for a closed parameter path (see Fig. 6c). In order to go from one topological phase ν in the bulk to another here given by smoothly varying Ω 1 and Ω 2 , one needs to cross a point (Ω 1 = Ω 2 ) where the band gap is zero which would violate adiabaticity. This shows that the Hamiltonians for the two insulating phases occurring for Ω 1 > Ω 2 and Ω 1 < Ω 2 are not adiabatically equivalent which is manifested in the difference of the topological invariant. In the case of open boundary conditions one finds that the number of edge states forms a topological invariant as well [71].

B. Applications of quantum-emitter rings and chains
Subradiance can be exploited towards applications in quantum metrology for sensitive frequency detection as well as in quantum information for the engineering of robust quantum memories. Moreover, subradiance in symmetric arrangements, such as rings, can be utilized to design nanoscale light sources acting as thresholdless nanolasers. To describe these effects, we make use here of the collective master equation in the single excitation regime introduced in Sec. II C and of the Bloch sphere representation introduced in Sec. II D.
Quantum metrology -Ramsey interferometry is routinely used in quantum metrology for the most sensitive measurements of optical clock frequencies. The conventional method of separated oscillatory fields [72] assumes an ensemble of atoms initially in the ground state illuminated by a laser of frequency ω ℓ in two consecutive steps separated by time τ . The two pulses are assumed to be instantaneous and tuned as π/2 pulses that can be visualized as π/2 rotations around the y axis on the collective Bloch sphere of radius N /2 (introduced in Sec. II D). After the interrogation time τ the population difference ⟨S z ⟩ (ω ℓ ) is monitored as a function of the laser frequency. Its behavior is sinusoidal where the argument is the total accumulated phase (ω 0 − ω ℓ )τ stemming from the mismatch of the (drifting, variable) laser frequency ω ℓ and the constant frequency separation of the emitter ground and excited states ω 0 . Analysis of the monitored population difference curve indicates a minimal sensitivity: where the minimization is performed with respect to ω = ω 0 − ω ℓ . For weakly decaying emitters (γ ≪ τ −1 ), the sensitivity is simply 1/(τ √ N ) and can obviously be optimized by using longer interrogation times and more emitters. However, longer interrogation times bring decay into play, while operation at higher densities achieved by increasing N indefinitely are marked by the onset of cooperative effects such as dipole-dipole shifts and super/subradiance. For independently decaying emitters, the application of the master equation introduced in Sec. II leads to a simple analytical solution for both ⟨S z ⟩ and ∆S z allowing one to estimate [δω] indep = e γτ /(τ √ N ). Further optimization with respect to the interrogation time gives an optimal τ opt = 1/γ and optimal sensitivity eγ/ √ N , which shows that the main impediment of Ramsey interferometry subject to radiative loss is the limited interrogation times available. This expression hints towards a deterioration of the sensitivity in the case of equally illuminated dense ensembles where symmetric collective states are addressed which are typically characterized by superradiant behavior.
To protect against such detrimental effects, Refs. [10,66] introduced an alternative procedure which uses an additional step in the Ramsey sequence aimed at hiding the collective states into decoherence-free subspaces. To this end, one complements the π/2 pulse with a phase distribution pulse, which for a particular atom j is represented by a rotation around the z-direction with the angle φ . The first generalized Ramsey pulse operator of such an asymmetric Ramsey technique is then The action of this phase distribution pulse is illustrated in Fig. 7a on the single particle Bloch sphere. In the second step, after the free evolution where robustness is now expected owing to the folding of the collective state into the subradiant part of the Bloch sphere, at time τ the phase spread is (instantaneously) reversed and a π/2 pulse follows leading to the second generalized Ramsey pulse Finally, detection takes place as before and the sensitivity is minimized with respect to frequency to obtain results as shown in Fig. 7a. Here, the sensitivity for symmetric illumination is worse than that for independent emitters while the interrogation times and consequently the minimum frequency sensitivity can be considerably increased for phased excitations with m ̸ = 0. Further results and analytical considerations can be found in Refs. [10,66].
Quantum memories -We now restrict the discussion to the single-excitation manifold of a 1D emitter chain. Here, we aim at targeting collective subradiant states as they show both robustness against decoherence while exhibiting multipartite entanglement. As pointed out in Sec. II C, eigenstates of the dipole-dipole Hamiltonian with lower energy typically exhibit subradiance. To access such states, one can proceed by first selecting them with individual addressing by tailoring the laser light amplitude to fit the shape of the collective state one wishes to address. For example, for a collective state| k⟩ introduced in Sec. II C one can provide geometrical matching with the following driving Hamiltonian Moreover, imposing the condition for resonance by setting ω ℓ = ϵ k ensures that states far enough from the desired targeted one are only weakly populated. In Ref. [59] it is shown that enhanced lifetimes much larger than (2γ) −1 can be reached by such tailored excitation. While tailored phase excitation with subwavelength resolution might pose great challenges, alternative methods could be envisioned: for example, symmetric addressing could be combined with the application of a magnetic field gradient. The effect of magnetic field gradient applied along the direction of the emitter chain is to progressively shift the excited state by a quantity ∆ j = 2∆ B (j − 1) from the first emitter with j = 1 to the end of the chain j = N . During the duration τ of an applied laser pulse, this amounts to a rotation around the z axis of the Bloch vector of each emitter by the angle ∆ j τ . For a conveniently chosen duration τ , the effective phase difference between neighboring emitters can be controlled: for example with the choice ∆ B τ = π a completely antisymmetric superposition can be constructed. For example, for two coupled emitters, Fig. 7b shows results for pulsed adiabatic transfer of population from the ground state to the antisymmetric state by a simultaneous off-resonant drive of the symmetric state and the action of the magnetic field gradient. The population then ends up in the protected antisymmetric state which then decays much slower than the timescale defined by (2γ) −1 . More results and analytical calculations for larger systems sizes can be found in Ref. [59].
Nanoscale coherent-light sources -Subwavelength spaced ensembles of quantum emitters in ring-like configurations resemble the structure of certain biological lightharvesting complexes (LHCs), which have been shown to act as an extremely efficient system of antennae in the photosynthetic process [73]. While the complexity of such biological systems is incredibly hard to tackle even at the computational level, a number of theoretical works have proposed simplified phenomenological models where, for example, the combined effects of disorder, vibronic coupling, electron-phonon couplings, are included as timedependent frequency shifts (in the Hamiltonian) or as dephasing (as an additional Lindblad term) [74,75].
Inspired by such naturally occurring systems, a recent theoretical proposal has shown the possibility of designing a thresholdless laser, i.e., a coherent light source [67] in the configuration depicted in Fig. 7d. Here, a central emitter acts as a gain medium and is pumped by incoherent light at some rate γ p . The role of the optical resonator (as present in standard lasing systems) is then taken by the optical modes defined by the geometry of the ring emitters (indexed by 1, . . . , N ). In order to test the properties of the system (emission rate, coherence of light, etc.) one can proceed to solve the master equation under the assumption that the incoherent pumping can be modeled by a Lindblad term with a σ † p collapse operator (that, opposed to the spontaneous emission case, takes population from the ground state to the excited state).
A first observation is that the physics of the system can be reduced to solely the interaction of the symmetric mode |ψ sym ⟩ = (1/ √ N ) j |j⟩ = σ † sym |G⟩ of the ring resonator with the central atom. The reduced Hamiltonian is then (46) where Ω sym = N j=2 Ω 1j is the dipole energy shift of the symmetric state and the decay of the system is governed by . For favorable geometries, it is possible to obtain a subradiant decay rate for the symmetric mode of the ring resonator. Additionally, such configurations can reduce the dissipative coupling of the central atom, strengthening the lossless transport of population from the central atom to the ring. We can now analyze the emitted light properties by using the fact that the far field is proportional to the sum of the dipole operators such that we can define the normalized second-order correlation function at zero time delay Numerical simulations, plotted in Fig. 7c show a close to unity g (2) (0), indicating a coherent state as an output of such a thresholdless laser. For higher pumping rates, antibunching of light is expected. Further results and considerations can be found in Refs. [67].

C. Optical response of 2D subwavelength mirrors
Let us now consider the situation depicted in Figs. 8a,b where a 2D periodic array of quantum emitters is positioned in the x-y plane. To derive the optical response of such a structure, we consider excitation in the form of a plane wave with wavenumber k = ω ℓ /c = 2π/λ impinging at normal incidence, along the z axis.
The electric field (source field) at some position R emitted by a collection of emitters, each located at R j , takes the following expression in the far-field limit [17,69] (only positive frequency component and in the direction of the emitter dipole moment) The total amplitude is given by the sum of dipole-radiated fields, which, in the linear regime, are proportional to the expectation value of the individual particle dipole operators β j (in a frame rotating at the laser frequency). We aim at describing regimes where plane waves propagating through the 2D structure either get reflected or transmitted only in the z direction while scattering in other directions is inhibited. To this end, we assume a constant phase illumination (denoting this regime as a symmetric driving case) where the electric field amplitude does not depend on x and y. The total field can then be written as the sum of the incident field and the field radiated by the emitters Let us first look at the case of just a single emitter where the radiated field is a dipole pattern with the z-direction amplitude falling off with increasing distance z. We have introduced the single-atom polarizability α p = −d eg β/E in which can be expressed in terms of the resonant wavelength, linewidth and ω ℓ as The extinction cross section of a quantum emitter σ ext (i.e., the effective area seen by an impinging photon) can be related to the polarizability via σ ext (ω ℓ ) = kℑ(α p )/ϵ 0 which in the case of resonant illumination simply becomes σ ext (ω ℓ = ω 0 ) = 3λ 2 0 /(2π). Notice the interesting aspect that the cross section is much larger (proportional to λ 2 0 ) than the square of the actual size of the electronic orbital which can be a milion times smaller.
For many emitters, the driving for the symmetric illumination case can be included as Transforming into a frame rotating at the laser frequency ω ℓ yields the equations of motion for the dipole amplitudeṡ Here, the collective effects coming from the coherent and incoherent emitter-emitter interactions are contained in the off-diagonal elements of the matrix M jj ′ = (iΩ jj ′ + γ jj ′ ) (while the diagonal elements M jj = γ jj give the independent decay dynamics). The assumption of constant drive η over the whole array combined with the assumption that the array is quasi-infinite leads to a trivial solution where β j = β for any j = 1, ...N can be easily estimated from the equation above as . This in practice means that only a symmetric superposition of all emitters is excited (only the symmetric surface mode is activated). Notice that the sum N j ′ =1 M jj ′ does not depend on the index j which is why we set it to 1 above and the result can be cast into the form i The term Ω eff leads to a shift of the collective resonance of the array (as also observed in 1D configurations [76]), while γ eff describes the effective decay rate of the array. Also notice that, as we will also detail more in the end of the subsection, the solution β j = β holds solely un-der constant illumination conditions. Generally, for a laser drive of the form E(z)f (x, y), i.e., showing phase or spatial imprinting in the transverse direction, different surface modes of the array would be driven (depending on the overlap integral between the function f (x, y) and the transverse profile of the surface mode.
The sum over the spatial distribution along the zdirection can be estimated by a plane-wave expansion as (for derivation see Appendix A3) where k mn = k 2 − q 2 m − q 2 n with q m = (2πm/a) (q n analogously). Here, q m and q n represent the quasimomenta of the surface modes on the array propagating along the x-and y-direction (Appendix A3). The wave number k mn is real if (a/λ) 2 ≥ m 2 + n 2 . For a subwavelength lattice, this can then obviously only be fulfilled with the choice m, n = 0: this means that only a mode with a k vector equal in amplitude to the impinging laser can propagate in the z-direction. This is exactly the symmetric mode which is also the only one that constant illumination can activate. Notice that, even if the illumination phase would slightly vary over the array, any surface mode which is accidentally excited would radiate in directions other than z. One can finally express the source field radiated by the dipoles in the far-field limit as with a renormalized effective polarizability summed over the whole ensemble response as α eff p = −d eg j β j /(N E in ) and expressed as We illustrate in Fig. 5c the behavior of the collective rates Ω eff , γ eff as a function of the lattice constant a/λ 0 . We note that similarly to the two-particle interactions discussed in Sec. II A, the real and imaginary parts are not independent but can be connected by a Kramers-Kronig relation. Regions with γ eff < γ (negative values in Fig. 8c) correspond to collective subradiant behavior and are of particular interest. One can see that a special operation point occurs e.g. at a/λ 0 ≈ 0.8 where the collective frequency shift vanishes while a pronounced subradiant behavior remains, facilitating the experimental realization of subradiant optical mirrors with cold atoms in an optical lattice [13]. Notice that the effective cross section per emitter can be considerably increased as, on resonance, one has σ eff ext (ω ℓ ) = σ ext (ω ℓ )γ/γ eff . For a < λ 0 one can furthermore approximate the effective decay rate as γ eff = 3γ(λ 0 /a) 2 /(4π) [see Fig. 8c] [69], showing a decrease by a factor of roughly 2.68 around the optimal operation point at a/λ 0 ≈ 0.8. This can also be connected to an increase in the overall reflectivity of the array. To derive this, we write E dip (z) = r(ω ℓ )E in e ik|z| where the complex reflectivity amplitude reads (considering small detunings around the resonance) while the transmission amplitude is obtained as t = 1 + r.
In Fig. 8d we plot the absolute square of these quantities as a function of the lattice constant for resonant illumination ω ℓ = ω 0 . One can see that, while generally the reflectivity of the array is high over a broad range of separations, for certain values of a/λ 0 ≈ (0.2, 0.8) the atomic dipoles can even act as a perfect mirror and reflect the entire input field with unit efficiency [69]. Moreover, for the whole region where the approximation γ eff = 3γ(λ 0 /a) 2 /(4π) holds, the mirror shows no losses (quantified by the scattering outside the z axis mode), i.e., |r| 2 + |t| 2 = 1. This is of course only valid in the absence of any other channels of nonradiative decay (at rate γ nr ), in which case the denominator of r acquires an extra term i(γ eff + γ nr ). The resulting intensity distribution of the total electric field is shown in Fig. 8(e), revealing that the emitters can indeed shut off the transmission of an incident plane wave for z > 0. Let us stress that the expressions above are only valid under symmetric illumination conditions, which is the relevant experimental situation as tackled in Ref. [13]. As optical lattices have interatomic distances at the level of λ 0 /2 or larger, the simplest operation is around a point where Ω eff vanishes (at a = 0.8λ 0 ) and unit reflectivity is reached as soon as the laser is resonant to the emitter transition ω ℓ = ω 0 . However, the subradiance of such a symmetrically excited collective state is only a factor of around 4π × 0.8 2 /3 ≈ 2.68 (at a = 0.8λ 0 ) smaller than γ. To fully exploit the scalability of subradiance with N , one could instead assume antisymmetric, or phased illumination conditions. For example, Fig. 8f shows that extremely narrow resonances scaling as N −5.5 can be reached for very dense arrays at a = 0.1λ 0 with antisymmetric phases in a checkerboard pattern where f nm = (−1) n+m (for n = 1, ..., √ N and m = 1, ..., √ N ) (deviations from the scaling are due to imperfect addressing of subradiant states). Let us only sketch how the derivation above will change in such a case. First, Eq. (49) is changed to include the x and y spatial dependence of the incoming field in some function f (x, y) with the periodicity equal to the lattice constant a. Then, the combination of surface modes which are excited by such an illumination is computed from the steady state of Eqs. (51) where each β nm is driven by an amplitude ηf nm . Finally, replacing the new solution for β nm into Eq. (48), the sum in Eq. (52) needs to be recomputed. Finally, the specific modes which can propagate into the far field in the z direction will depend on the specific phased illumination pattern chosen.

D. Further remarks
One-and two-dimensional ensembles of coupled quantum emitters arranged in regular patterns provide an ideal platform for achieving strong light-matter interactions and high fidelities for photon storage capabilities [2,3]. Their subradiance properties are an important resource for applications ranging from quantum-information processing [9] and metrology [2,10,11] to excitation transfer [77][78][79] and it is envisioned that one can build quantum matter in a bottom-up approach from nanoscopic lattices of atoms and photons [9]. A single subradiant array has been shown to act as an ideal quantum memory with efficient storage and retrieval [11] and it has been suggested [13] that this geometry could lead to vast improvements in the error bound of quantum memories. Regular or honeycomb lattices of three-level systems in a V -configuration, where time-reversal symmetry is broken by the application of a magnetic field, have been proposed as platforms for studying topological phenomena with strongly interacting photons [4,5] and further improvements have been proposed in the form of interfacing the arrays with two-dimensional photonic crystals [6]. A crucial aspect of these proposals that distinguish them from linear topological photonic systems is the intrinsic nonlinearity of the quantum emitters which could lead to a rich many-body physics dynamics on such subradiant lattices. Furthermore, composite quantum systems comprised of many atomic arrays could find applications in quantum networking: at the level of two distant layers, nonlocal entangled Bell superposition states have been shown to exist [8]. While quantum information processing at the level of atomic layers require qubit encoding on delocalized spin states over the whole array, quantum spin lenses have been recently introduced, where incoming flying qubit photons can be mapped and stored in single atoms [7].

IV. COOPERATIVITY IN CAVITY QED
Light-matter interactions can be greatly enhanced by using optical elements which confine electromagnetic fields in very small volumes. This is obvious from the scaling of the zero-point electric field amplitude E k = ω k /(2ϵ 0 V) which indicates that stronger field amplitudes per photon mode are achieved for smaller mode volumes. This led to the development of cavity QED as a subfield of quantum optics specializing in the description of coherent light-matter interactions in cavities [19][20][21][22]. This section provides fundamental concepts and tools of cavity QED with coupled quantum emitter systems, which are then utilized to describe applications in Sec. V.
We start by introducing a master equation approach to intracavity light-matter interactions encompassing both collective emitter loss and cavity photon loss. We then move on to a quantum Langevin equations approach sup-plemented with an input-ouput formalism allowing access to the correlations of cavity output field operators. Finally, to address inhomogeneous emitter ensembles as often present in experimental setups, we tackle the question of strong frequency disorder effects on light-matter interactions.
A. Cavity QED with coupled quantum emitters The simplest example of an optical cavity is the coplanar design known as a Fabry-Pérot cavity comprised of two highly reflective parallel mirrors. For a distance ℓ between the mirrors and assuming at first perfect reflectivity, such a setup defines resonances conditioned by ℓ = nλ/2. The fundamental mode has, therefore, a wavelength of 2ℓ with frequency πc/ℓ and higher harmonics are multiples of this mode. The quantity ω FSR = πc/ℓ is known as the free spectral range and gives the frequency difference between consecutive optical resonances of the structure. Assuming two mirrors positioned at z = 0 and z = ℓ and a transverse area S of the supported optical resonances, one can proceed in quantizing the field inside the optical resonator by associating bosonic operators a n to all resonances and writing the total Hamiltonian as H c = n ω n a † n a n . The electric field operator can then be decomposed aŝ where the sum runs over all allowed wave vectors k n = nπ/ℓ and the zero-point electric field amplitude E n = ω n /(ϵ 0 ℓS) and all polarizations ϵ n are assumed orthogonal to the z-direction. For all situations we will consider in the following, a single cavity mode suffices. We denote the mode by the operator a at frequency ω c such that the free Hamiltonian becomes H c = ω c a † a.
ME for a driven, lossy cavity mode -Perfect mirrors define infinitely sharp resonances; in order to allow the in-coupling of light into the optical resonator side mirrors with slightly less than unit reflectivity are used. For a double-sided optical resonator, the localized resonances can then couple to the infinite number of modes to the left and right of the cavity, leading to a loss of the intracavity photons. This loss can be described within the master equation formalism by following a phenomenological model where photons from mode a can tunnel to the left b(ω) and right c(ω) continuum of free radiation modes via an excitation exchange Hamiltonian. The standard open-system dynamics approach detailed in Appendix A1 leads to a Lindblad form of the cavity photon decay where κ = κ R + κ L encompasses total losses via both the right-side and left-side mirrors.
To allow for the driving of the cavity mode let us now assume a continuous-wave laser with power P entering the cavity from the left side. This can be included in the following Hamiltonian where the drive amplitude is η = 2κ L P/ω ℓ . The equation of motion for the expectation value of the cavity operator can be derived from the master equation (with Hamiltonian H c + H ℓ and Lindblad term L κ [ρ]) leading to The steady state of the equation above (reached by requiring that ⟨ȧ⟩ = 0) describes the expected Lorentzian response ⟨a⟩ ss = η/(κ + i(ω c − ω ℓ )) of the cavity field exhibiting a linewidth κ and resonance frequency ω c .
Tavis-Cummings Hamiltonian -Let us consider now N nonidentical quantum emitters with transition frequencies ω j placed inside the optical cavity around a single optical mode of interest. Assuming the field to be varying only along the cavity axis, for emitter j positioned at z j within the cavity volume, the dipolar interaction is characterized by a position-dependent coupling strength g j = Ed eg sin(kz j ). The light-matter Hamiltonian describing excitation exchange between the emitters and the cavity mode can then be cast in the standard Tavis-Cummings [80] form This interaction is a particular form of the free-space case in Eq. (3) where the coupling to the cavity mode a is enhanced while the coupling to all other modes is inhibited. At the single emitter level, the interaction is known as the Jaynes-Cummings Hamiltonian [81]. Notice that such Hamiltonian conserves the excitation number and can be exactly solved [80]. The excitation nonconserving case (known also as counter-rotating terms) for N = 1 is the exactly integrable [82], Rabi Hamiltonian [83]. The extension of the Rabi model to N > 1 for identical emitters ω j = ω 0 and identical couplings g j = g is a simplified case describing exchange between a single collective operator j σ j and the cavity mode. This is known as the Dicke model, generally nonintegrable [84,85] except for when ω 0 = 0, and in the thermodynamic limit with N → ∞ [86][87][88] where phase transitions have been identified [89].
Linear optical response -For dense, interacting ensembles, the dipole-dipole interactions H dd defined in Eq. (10) and the collective decay terms as defined in Eq. (7)  From here, one can derive equations of motion for single operator amplitudes. The emerging set of coupled equations is not closed, i.e., expectation values of single operators are coupled to two-operator correlations which again couple to three or more operator correlations and so on. However, in a first approximation, assuming weak driving conditions (η ≪ γ) some factorization of operators under the approximation ⟨σ j z ⟩ ≈ −1 can occur. This approximation is sufficient to describe the linear optical response as used throughout Sec. V A. In Sec. V C instead, we will go beyond this approximation to describe population-inverted systems exhibiting lasing behavior.
Let us denote the expectation values by α = ⟨a⟩ and β j = ⟨σ j ⟩ and derivė In Sec. IV B these equations will be extended to the quantum regime within the quantum Langevin equations formalism.
Regimes of cavity QED -Let us quickly review a few standard concepts of cavity QED with non-interacting ensembles by setting Ω jj ′ = 0 and γ jj ′ = γδ jj ′ . In steady state, the normalized transmission of the cavity t = κα/η can be easily computed from Eqs. (61) to lead to The cavity transmission can then show different regimes either characterized by a single Lorentzian peak (empty cavity or weak light-emitter coupling), a peak with a narrow dip (the Purcell, or antiresonant regime for intermediate couplings) or a polaritonic regime, where two cavity resonances are resolved in transmission. These regimes are dependent on the magnitude of the loss rates relative to the coupling strength. Mathematically, they can be easily deduced from the eigenproblem defined by the evolution matrix of Eqs. (61). For N = 1, diagonalization of Eqs. (61) leads to normal mode splitting when g > |κ − γ|/2 with new frequencies ω ± = ω 0 ± g 2 + (κ − γ) 2 /4 (on resonance ω 0 = ω c ). The difference ω + − ω − is known as the vacuum Rabi splitting (VRS). The newly formed quantum states in such a strong coupling regime are hybrid light-matter ones, i.e., polaritons (|g1⟩±|e0⟩)/ √ 2 and can be read in the cavity transmission t as two well-separated peaks. For N > 1, the collective strong coupling regime can be defined when g N > |κ − γ|/2 with a collective coupling rate g N = j |g j | 2 . The loss rates of the polaritons are equal in this regime (κ + γ)/2. One can also identify a weak coupling regime for which γ < g N < κ but with a strong cooperativity C N = g 2 N /(κγ) ≫ 1. For many emitters equally coupled to the cavity mode, the cooperativity shows a collective enhancement proportional to N . For a single emitter, this regime (known as the Purcell regime) leads to a renormalization of the intrinsic spontaneous emission rate to γ(1 + C) understood as the addition of a new decay channel via loss through the cavity. A scan of t as a function of laser frequency shows in such a case a dip in the cavity transmission around ω ℓ = ω 0 known as a cavity antiresonance. The case of N coupled emitters will be presented in Sec. V A showing that the addressing of collective subradiant resonances can lead to a superlinear ∝ N 4 scaling of C N .

B. Input-output formalism for operators
The master equation formalism previously utilized is based on discarding the state of the environment and focusing on the dynamics of the reduced density operator of the much smaller system of interest. An alternative approach are quantum Langevin equations (also called Heisenberg-Langevin equations) [48,90] which follow the evolution at the level of system operators instead of system states, i.e., its density operator. In essence, the formalism consists in supplementing the Heisenberg equations of motion with the proper dissipation and fluctuation terms.
Mapping the ME onto QLEs -For open-system dynamics described by loss in Lindblad form, a direct transformation between the master equation and the QLE formulation exists that reads [48] The equation above for a generic system operator O is applied to each individual Lindblad collapse operator c acting at rate γ c and with associated input noise c in . The input noise is a stochastic operator which is zero averaged

all other correlations vanish) and satisfies the following
. For linear evolution, such an equation can be formally integrated and expectation values of correlations containing any number of operators can be obtained.
We will proceed by first exemplifying the derivation and solution to a QLE applied to a lossy bosonic field corresponding to a decaying and externally driven optical cavity mode. Then we introduce the input-output relations allowing the readout of intracavity light-matter interactions via the cavity output field. Finally, we introduce a set of coupled linear QLEs for an ensemble of coupled, collectively decaying quantum emitters within the same optical cavity mode and show how quantum correlations can be computed. We will then utilize QLEs in Sec. V A for the description of the optical response of subradiant arrays in optical cavities and in Sec. V B for the description of hybrid cavities. An extension of the QLEs to include electron-vibron coupling will be introduced in Sec. VI to analytically tackle molecule-photon dynamics.
QLEs and input-output relations. -For a driven, lossy, single-sided optical cavity, Eq. (63) yields the following QLĖ Notice that an expectation value of the equation above simply reproduces Eq. (59) which has been alternatively derived directly from the master equation. The more general Langevin equation also allows for an exact solution obtainable by formal integration. We separate this solution into a transient and steady state part a(t) = a tr (t)+a ss (t).
While the transient part a tr (t) decays away on a timescale defined by the cavity loss κ −1 , the steady state part a ss (t) dominates in the long time (steady state) limit It is important to notice that the inclusion of the noise terms in such equations is crucial: while the commutator [a tr (t), a † tr (t)] = e −2κt decays away, the commutator [a ss (t), a † ss (t)] = 1 − e −2κt maintains the standard bosonic commutation relation [a(t), a † (t)] = 1 at every time during evolution. In addition to the equation of motion for the cavity field amplitude operators, the full description of the optical cavity includes the input-output relation which allows for the derivation of the properties of the cavity output mode. For example, the average shows that the cavity transmission can be simply written as: QLEs for cooperative cavity QED -We now assume a cavity mode coupled to N interacting emitters, as a system described by the Hamiltonian H TC + H dd + H ℓ and loss rates incorporated in the terms L κ [ρ] + L e [ρ] (as introduced in Sec. IV A). In a first step, we notice that the decay terms are not in Lindblad form, thus we perform a diagonalizing transformation to bring it into the form of N independent decay channels as described by Eq. (12) (such that Eq. (63) can then be directly applied). We also only consider the low excitation limit ⟨σ j z ⟩ ≈ −1. Notice that this approximation discards the nonlinear response of the dipoles, which could, for example, lead to a collective Kerr nonlinear effect [17]. One then obtains a set of coupled linear differential equationṡ with the emitter and cavity detunings defined as ∆ = ω ℓ − ω 0 and ∆ c = ω ℓ − ω c respectively. The input noise terms σ j in (stemming from the coupling of the emitters to the electromagnetic modes outside the solid-state angle covered by the cavity field) are considered here as zeroaveraged and delta-correlated in time ⟨σ j in (t)σ j ′ , † in (t ′ )⟩ = (γ jj ′ /γ)δ(t − t ′ ) as a consequence of the low excitation assumption (for the more complex case of high excitation we refer to Refs. [17,48]). We can write the equations above in a compact form in terms of the matrix M(∆) = −i∆1 + iΩ + Γ and the coupling vector G = (g 1 , g 2 , ..., g N ) ⊤ . Additionally defining the vectors σ = (σ 1 , . . . , σ N ) ⊤ and σ in = (σ 1 in , . . . , σ N in ) ⊤ , Eqs. (67) express aṡ This is the starting point for the derivations in Sec. V A where we compute the classical and quantum response of a weakly excited transverse 1D or 2D array of coupled quantum emitters to a driven cavity mode.
Quantum correlations -From the equations of motion for operators Eqs. (68), one can go a step further and fully analyze the quantum properties (second-order correlation function for light g (2) , bipartite entanglement, squeezing properties, etc.) of both the photon and matter counterparts. To this end, let us first consider only fluctuation operators only by expansion around expectation values a = α + δa, σ j = β j + δσ j . We can then cast Eqs. (68) in convenient vector form asv = Av + Nv in for the vectors v = (δa, δa † , δσ, δσ † ) ⊤ , v in = (a in , a † in , σ in , σ † in ) ⊤ with the drift matrix expressed in compact form as (for details also see Ref. [17]) We have defined the vector 0 containing N zeros and 0 is a N × N matrix with only zeros. The matrix multiplying the input noise operators is (0) FIG. 9: Second-order optical correlation function g (2) (0) of a driven emitter-cavity system as a function of the cavityemitter detuning ωc − ω0. Parameters: γ = 0.1κ, η = 5 · 10 −3 κ, g = 0.3κ.

Formally integrating this system of linearly coupled equations gives the solution
If all eigenvalues of the drift matrix are negative, the system is stable and will go towards a steady state (i.e., the transient term e At vanishes). From the steady state solution only, one can then define a correlation matrix V = ⟨v(t)v ⊤ (t)⟩ which can be expressed as where we used that ⟨v in (t ′ )v ⊤ in (t ′′ )⟩ = Cδ(t ′ − t ′′ ) and defined the diffusion matrix as D = NCN ⊤ . The matrix containing the correlations can be readily obtained as For the covariance matrix, one can then derive the Lyapunov equation using integration by parts Alternatively, one can also perform a Fourier analysis to turn the system of differential equations to an algebraic set of coupled equations iωv(ω) = Av(ω) + Nv in (ω), which allows to relate the intracavity operators to the input noise via Furthermore, one can relate the input fields directly to the output fields by using the fact that in the time domain v out (t) = N ⊤ v(t) − v in (t) and therefore one can express v out (ω) = F(ω)v in (ω), where In the frequency domain, the δ-correlations are preserved for the output fields ⟨v out (ω)v ⊤ out (ω ′ )⟩ = S out (ω)δ(ω + ω ′ ), where the two-operator correlations are completely encoded in the spectrum matrix given by While this suffices to characterize the two-operator correlations of the system, also higher-order correlations (e.g. four-operator correlations) can be of interest. We note that all higher order correlations can be expressed as a sum over products of two-point correlations via the Isserlis' theorem. An example of such a correlation is the the g (2) -function already considered in Sec. III C, characterizing the photon statistics of the system. It is defined in steady state as Nonclassical light sources show sub-Poissonian statistics g (2) (0) < 1 (antibunching), implying that it becomes unlikely for photons to be detected in pairs. In Fig. 9 this is illustrated for the coupled cavity-emitter system where the mechanism originates from the anharmonicity of the Jaynes-Cummings ladder, allowing only single photons inside the cavity.

C. Cavity QED with disordered ensembles
In the most general case of Eqs. (61) the frequencies of the emitters are not identical. This is often encountered in experiments, as the coupling of electronic systems to embedding matrices can lead to strong inhomogeneous broadening, such as it is the case for quantum dots or molecular systems. To quantify frequency disorder, we assume that the frequencies ω j are distributed around ω 0 according to the distribution function p(δ) normalized to unity ∞ −∞ p(δ)dδ = 1. In particular we choose a Gaussian distribution of frequencies p(δ) = (1/ √ 2πw 2 )e −δ 2 /(2w 2 ) . We write the transition frequencies as ω j = ω 0 + δ j with vanishing classical average ⟨δ j ⟩ cl = 0 and variance ⟨δ 2 j ⟩ cl = w 2 . For simplicity, we restrict the discussion to the case of identical couplings g j = g (for all j) while more generally we refer the reader to Ref. [91].
For w = 0 we note that the cavity couples only to a symmetric superpositionB = j σ j / √ N , i.e., a bright state, with a collective coupling strength g N = √ N g. The other N −1 combinations define dark states which are obtainable by a Gram-Schmidt algorithm that leads to all vectors orthogonal to the bright state and to each other. However, a simple choice of coefficients is indicated by a discrete Fourier transformD k = 1/ √ N N j=1 e −i2πjk/N σ j . We index the dark-state manifold for k = 1, . . . , N − 1 and note that for k = N we haveD N =B. The equations of motion for all variables D k = ⟨D k ⟩ and α = ⟨a⟩ become (in a frame rotating at ω 0 ) where the couplings between collective states are defined as ∆ kk ′ = 1/N N j=1 δ j e −i2πj(k−k ′ )/N . For k = N , the equations above indicate that the bright state couples to the cavity mode with the standard g N rate, while also being coupled to all dark states.
Non-Markovian regime -We proceed (see Fig. 10) to eliminate the dark-state manifold in an exact way without making a Markovian approximation (which would imply that the dark state reservoir has no memory and therefore it would allow to set all derivatives of D k to zero). Instead, we formally integrate the equations for D k to obtaiṅ where the memory kernel generally describes a non-Markovian loss process. In the mesoscopic limit N → ∞, one finds This now allows to identify the Markovian regime discussed in the next paragraph, where the condition that the disorder w ≫ γ implies that the kernel f (t − t ′ ) tends towards a delta function.
Markovian limit -In the Markovian limit, the elimination of the dark-state manifold is straightforward as all the derivativesḊ k for k ̸ = N can be set to zero. This , containing an effective frequency shift δ dark and a loss rate γ dark derived from δ dark + iγ dark = In the mesoscopic limit, one can further simplify the expression of the loss rate and derive scaling laws where γ dark = w 2 /γ for γ ≫ w and γ dark = πw/4 for w ≫ γ. Moreover, both δ and δ dark vanish in this regime. The dependence of the VRS on disorder can then be obtained by diagonalizing the dynamics in the reduced cavity-bright state subspace The expression above shows that a large degree of disorder (on the order of κ) can lead to a strong reduction of the VRS and consequently pull the system out of the collective strong coupling regime. In Ref.
[91] a more detailed derivation is provided that shows that particles which have a large detuning with respect to the cavity frequency, or a very lossy behavior are effectively pulled out of the macroscopic superposition participating in the strong coupling condition. This in turn leads to the degradation of the VRS described by the equation above.

V. APPLICATIONS IN CAVITY QED
We now discuss a few applications of cavity QED with coupled quantum-emitter ensembles. For 1D and 2D arrays of emitters, we now have systems exhibiting cooperative surface resonances, as already analyzed for free space applications in Sec. III, but with an extra degree of freedom introduced by a cavity-confined field mode. In the limit of low reflectivity, such systems have been shown to act as quick frequency switchers [16] and to give rise to enhanced optical nonlinearities [17]. In the high reflectivity limit, such as theoretically predicted in Refs. [68,69] and experimentally proven in Ref. [13], we introduce the necessary steps, as laid out in Ref. [28], to derive the correct input-output formalism for optical resonators made up of subwavelength arrays as end-mirrors.
In the weak reflectivity limit, the array reacts as a strongly dispersive optical element within a short frequency window but it does not, at the same time, considerably change the spatial profile of the cavity mode (owing to its weak reflectivity). The treatment is then perturbative, within the standard Tavis-Cummings formalism and nontrivial effects occur such as enhanced optical nonlinearities and the reach of an enhanced collective Purcell effect with a cooperativity scaling up to ∝ N 4 (for a 1D arrangement).
In the opposite case, where the emitter array acts as a near-unity reflectivity end-mirror, a hybrid cavity design emerges where asymmetric transmission profiles can be achieved, potentially much narrower than those obtained with frequency-independent mirrors of comparable reflectivity. For such designs, the standard input-output theory and the master equation for the photon mode losses (from Sec. IV A) loses validity, as the tunneling rate of photons strongly depends on their frequency. We present a general roadmap to derive the correct input-output relations for optical cavities comprised of either one or two of such end-mirrors. We remark that this is not limited to subradiant emitter arrays but also extendable to patterned subwavelength gratings or photonic crystal structures [92][93][94][95] and semiconducting monolayers [96][97][98].
Finally, for ensembles of cavity-embedded quantum emitters, externally and incoherently pumped, we provide a minimal theory of lasing to illustrate how the laser threshold condition can be easily derived using the formalism introduced in Sec. IV A and by going beyond the linear regime. In particular, we analyze the physics of superradiant lasers, where coherence is being stored in the gain medium instead of in the cavity photon field. Finally, we connect the dynamics of such a laser to the Dicke superradiance model introduced in Sec. II D.
A. Antiresonance spectroscopy with 1D or 2D arrays Let us first consider a weakly reflecting array placed transversely to the axis of a standard optical cavity. In the presence of near-field couplings, the system undergoes dynamics describable at the level of averages by Eqs. (61) and at the fully quantum level by Eqs. (67). We will assume a vector G containing all the cavity-emitter couplings and later particularize to symmetric coupling (where all couplings are equal to g) and antisymmetricallyphased coupling (with alternating (−1) j g couplings). Assuming a single cavity drive, one can then easily deduce the cavity transmission as where the effective ∆-dependent collective energy shifts and linewidths are derived from the matrix M(∆) = −i∆1 + iΩ + Γ as real and imaginary parts These expressions are the equivalents of the decay rates and shifts for free space arrays given in Sec. III C where Enhanced cooperativity -Notice that one can now proceed by introducing a modified N emitter effective cooperativity As mentioned in Sec. IV A, for a single emitter the cooperativity is independent on the emitter properties (dipole moment, decay rate) and simply depends on the cavity finesse. For N uncoupled emitters the same holds true and only a linear increase with N will be obtained. For coupled emitters, an interesting decoupling of the dipole moment in the direction of the cavity (quantified by the term G ⊤ G) from the collective radiative properties γ eff (∆) can be achieved. As γ eff is not a natural constant of the ensemble, but strongly dependent on the relative positioning and phase of individual emitters, one can reach subradiant states with γ eff ≪ γ. By proper design of the cavity transverse field amplitude profile, the numerator can at the same time be maximized, resulting in a scaling up of C eff well above the independent emitter case N g 2 /(κγ). This is illustrated in Fig. 11a where a scaling with N 4 for a 1D chain is shown possible. In Figs. 11b and c, we consider scans of the cavity transmission for a system of four coupled emitters. The system will have three subradiant and one superradiant state with resonances given by Eq. (14) (the superradiant state is located at ω 0 + 2Ω cos(π/5)). We adress the system either antisymmetrically with G = (g, −g, g, −g) ⊤ [ Fig. 11b] or symmetrically with G = (g, g, g, g) ⊤ [Fig. 11c]. While antisymmetric excitation overlaps with two subradiant states, symmetric excitation overlaps with one subradiant and one superradiant state. In Ref. [16] it is shown that phased excitation could be realized by using the higher-order TEM-modes of the optical cavity. Alternatively, one could use symmetric addressing around the optimal points utilized experimentally in Ref. [13] or symmetric addressing combined with a magnetic field gradient as in Ref. [59].

B. Hybrid cavities with 2D subwavelength mirrors
We now aim to provide a theoretical framework for hybrid cavities where subwavelength arrays are to be used as end mirrors as suggested in Refs. [13,28,69]. The standard approach for cavity loss and input-output relations taken in Sec. IV is based on the assumption that the tunneling rate of photons between a cavity-confined mode a and the outside continuum of modes (b(ω) to the left and c(ω) to the right) is flat around the cavity resonance. This assumption is strongly modified in the case of subwavelength reflective arrays and the presence of a surface-confined resonance (mode d) has to explicitly be taken into account (see illustration in Fig. 12a). The roadmap to construct an input-output quantum theory for such hybrid cavities follows the steps described in Ref. [28] where it has been applied to photonic crystal mirrors. The procedure consists of three steps: i) the derivation of equations of motion for cavity field operators based on a coupled modes model where the photon-exchange processes between modes a and d are included, ii) the derivation of the cavity amplitude transmission from a classical transfer-matrix approach based on the expression of the reflectivity of the array r(ω) derived in Sec. III C and iii) the extraction of the phenomenologically introduced parameters by matching the predictions of the two theories. We apply this procedure here to a hybrid cavity made of one flat mirror and one subwavelength array and point out expected improvements in the cavity finesse where the hybrid design makes use of two emitter array mirrors.
Coupled modes theory -We consider the situation depicted in Fig. 12a where the cavity mode a is coupled to the surface-confined mode d at some complex rate G and the photon tunneling rates U, V, W ∈ R are frequency independent. Tunneling will then give rise to loss rates κ L = πV 2 , κ R = πW 2 , γ d = πU 2 (see Appendix A1). The dynamics is then described by the Langevin equations showing that the mode a is subject to two types of noises: b in , entering through the left and c in , entering through the right mirror. In a first approximation, we will consider that the surface-confined mode reacts to the same noise term b in consisting of modes on the left of the cavity but extra independent noise terms stemming for example from nonradiative losses could be added. The total decay rate is κ = κ L + κ R and the input fields have the usual nonvanishing correlation functions ⟨b in (t)b † in (t ′ )⟩ = δ(t − t ′ ) [and similarly for c in (t)]. The associated output fields follow the input-output and c out = c in − √ 2κ R a. We are mainly interested here in deriving the free parameters in the coupled modes model, i.e., G, ω a and κ L . For more in-depth discussions highlighting the non-Markovianity of such hybrid cavities we refer to Ref. [28] where a full analysis of photonic crystal mirror cavities is provided.
The transmission coefficient t(ω) = ⟨c out (ω)⟩ / ⟨b in (ω)⟩ is easily expressed as t(ω) = − √ 2κ R ⟨a(ω)⟩ / ⟨b in (ω)⟩ in the case of sole driving through the left cavity port and can be analytically derived from Eqs. (86) in the frequency domain where we introduce the susceptibilities ϵ −1 a (ω) = κ + i(ω a − ω) and ϵ −1 . This is the main result of the coupled modes theory, which at this point still retains complete generality; in the following we will apply it to the particular scenario of a subwavelength reflective array.
Transfer matrix results -The transfer matrix approach, on the other hand, consists in solving the classical onedimensional wave propagation in a one-dimensional setup with two mirrors parametrized by the polarizabilities ζ 0 (right mirror) and ζ L (ω) (left, subwavelength array mirror). In linear response theory, one can find the transmission function of the setup for any incoming plane wave at a given frequency ω. We use the parametriza- for the subwavelength array [cf. Eq. (55)] and choose a fixed ζ 0 ≫ 1 (corresponding to a close to unity reflectivity) for the right mirror. The classical transmission coefficient (see Appendix A4) then reads (88) where θ = ωℓ/c. For ζ L (ω) and ζ 0 infinite, the expression above reaches unit absolute value at resonances ω m = m × ω FSR for any positive integer m. The resonances are separated by the free spectral range ω FSR = cπ/ℓ. Let us first fix a given resonance number m such that ω m lies in the neighborhood of the mirror resonance and see that for finite ζ 0 and assuming ζ L (ω) is flat in frequency and equal to ζ 0 , the transmission can reach unity at a shifted ω ′ m = ω m + ω FSR /(πζ 0 ). The linewidth of such a resonance is then κ flat = ω FSR /(2πζ 2 0 ). When the two mirrors have unequal reflectivities, the transmission of the optical resonator is always less than unity. Therefore, for variable susceptibility ζ L we ask now that a resonanceω should be reached for which |ζ L (ω)| = ζ 0 such that the transmission is unity. There are two solutions for this equation and we pick the one where the mirror resonance sits to the left of the cavity resonance. The reason is that the reflectivity at ω d is exactly unity which means the cavity transmission reaches zero. With the choiceω = ω d + γ d /ζ 0 the zero of the cavity sits on the left of the cavity resonance. We also fixω = ω m , which in practice means that one first determines γ d , ζ 0 and then adjusts the cavity length.
For varying γ d , Fig. 12b shows a scan of the cavity resonance. For γ d ≫ κ flat , the expected Lorentzian response is obtained (gray dashed line) while for decreasing γ d at the level of κ flat , a very narrow (around γ d /ζ 0 ) asymmetric Fano profile is obtained. The zero of the hybrid cavity transmission is at ω d while unity is reached at ω m .
Notice that the double hybrid cavity has an even narrower linewidth as illustrated in Fig. 12c for the same regimes as provided for Fig. 12b. In addition, the transmission profile is symmetric and has the advantage of presenting two Fano resonances, situated symmetrically with respect to the cavity resonance.
Extraction of parameters -Comparison of the two expressions in Eq. (88) and Eq. (87) can allow for the identification of G, ω a and κ L . In a first step, we ask that the zero of the transmission at the mirror resonance vanishes t(ω d ) = 0 corresponding to the point at which the mirror reflectivity is unity. This leads immediately to the identification G = √ κ L γ d .
The derived coupling between the surface resonance and the cavity-confined mode is purely dissipative, in stark contrast with the situation treated in Ref. [28], where an additional real photon exchange process between modes a and d occurs. In the next step, we ask that unit transmission |t(ω)| 2 = 1 is reached at the resonance ω m = ω d + γ d /ζ 0 . From Eq. (87), after a few calculations, we obtain a solution for the left mirror loss as As the loss rate is defined in the real domain, we ask for the imaginary part to vanish which gives ω a = ω m + ζ 0 κ R . This also leads to κ L = κ R (1 + ζ 2 0 ). Further fitting of the transmission profiles also identifies κ R = (κ flat /2)γ d ζ 0 /(ω FSR + γ d ζ 0 ).
With the proper definition of the parameters appearing in Eqs. (86), one can proceed in solving various problems involving hybrid cavities with quantum-emitter ensembles (eventually in the direction of Fano cavity lasing as in Ref. [99]) or with movable mirrors. In the direction of quantum optomechanics, Ref. [28] has shown that the Fano profile of photonic crystal mirror cavities can be constructively utilized to lead to cooling in the absence of heating. The mechanism is based on sideband resolution by fitting a Stokes, heating sideband inside the Fano resonance. For subwavelength emitter arrays the same can be realized by tailoring the relationship between the cavity length and the array resonance. In addition, the double-sided hybrid cavity provides a scenario in which both sidebands can be inhibited as it presents symmetrically placed Fano dips.

C. Superradiant lasers
The Tavis-Cummings Hamiltonian introduced in Sec. IV A together with the collective spin algebra from Sec. II D are sufficient for a basic understanding of the main advantages presented by superradiant lasers (as treated in Refs. [25,26,100]). First, we will assume a minimal model for a laser comprised of a cavity containing an incoherently pumped gain medium. Under certain conditions (implying the existence of a pumping threshold) a non-zero intracavity field amplitude can be obtained, which signifies a coherent-light output, i.e., lasing. The distinction between the good and bad cavity regimes (behavior quantified by the cavity losses) will then be seen to characterize the difference between standard and superradiant lasers. Finally, we will make the connection between superradiant lasing and superradiance as described by Eq. (21) in Sec. II D by deriving an effective master equation for the gain medium after eliminating the cavity field.
The minimal model assumes N identical emitters undergoing (independent) decay with the standard Lindblad term (as defined in Eq. (8)) L γ (each emitter has a decay rate γ with corresponding collapse operator σ i ). Additionally, the incoherent pump assumes a Lindblad term L γp (each emitter is pumped at a rate γ p with corresponding collapse operator σ † i which brings population to the excited state). The pump can for example be realized as illustrated in Fig. 13a by excitation into an intermediate level |i⟩ followed by fast decay to the lasing level |e⟩. For γ p > γ population inversion can be produced, however without exciting the emitter dipoles ⟨σ j ⟩ = 0 as all processes are incoherent. To produce a non-zero average dipole moment, which would correspond to the generation of nonthermal light, the coupling to the optical resonator is crucial. From the Hamiltonian in Eq. (60) (with g j = g and ω 0 = ω c ) supplemented with the two Lindblad terms (and in terms of collective operators defined in Sec. II D), we can derive the equations of motion: The above set of equations is generally not solvable. However, an expansion around averages a = α+δa, S = s+δS and S z = s z + δS z where all δa, δS and δS z are zeroaveraged quantum noise terms, can lead to a simplification. The observation (which one can eventually infer from numerical simulations) is that ⟨δaδS⟩ ≪ αs (and similar for all other two-operator correlations) in the limit that N is large. In such a case, a much simpler set of equations can be obtained for averageṡ In a standard laser, stimulated emission produces a large intracavity field amplitude exhibiting coherence past threshold. In a superradiant laser, the intracavity field amplitude is instead kept small and coherence is stored in the collective atomic dipole thus protecting the laser linewidth from thermal effects. (b) Behavior of intracavity photon number |α| 2 and coherence |s| 2 in the standard, good cavity regime (κ = 0.1g) and (c) in the superradiant, bad cavity regime (κ = 40g). Other parameters are γ = 10 −2 g and N = 200. Notice that the upper threshold γ + p ≈ N g 2 /κ is very different in both regimes.
Let us now assume that a steady state can be reached (by setting all derivatives to zero) in which case the solution of the above equations yields a population difference s z = 1/(2C p ), with C p = g 2 /[κ(γ+γ p )], the cavity-field intensity (photon number) and the dipole coherence |s|/|α| = κ/g. The expression above is positive only above a lasing threshold (signifying the existence of a lasing steady state), which can be derived by finding the required values for the pumping strength γ p . Requiring that in the lasing regime one has |α| 2 > 0, yields a quadratic equation with upper and lower thresholds given by If N g 2 /κ is large compared to γ, the thresholds can be approximated as γ − p ≈ γ, γ + p ≈ N g 2 /κ. Note that the lower threshold corresponds to achieving population inversion while the upper threshold implies that a too large incoherent pumping rate γ p prevents the formation of coherence in the system. While the analysis above is valid in any regime, the great advantage brought on by superradiant lasers is that the coherence is stored in the gain medium which means that the laser linewidth is unaffected by common problems in optical resonators operating at large intracavity powers (such as thermal vibrations of the end mirrors). This can be seen in the ratio |s|/|α| = κ/g which indicates that for κ ≫ g the intracavity field can be negligible while the gain medium coherence is very large [see Figs. 13b,c].
Let us now provide an alternative understanding of the superradiant lasing problem, by connecting to Sec. II D which we do by focusing on sketching the derivation of an effective master equation for the N emitters in the case of a lossy cavity. We assume at time t = 0 the emitters to be fully inverted, while the cavity mode is in the |vac⟩ vacuum state. We neglect spontaneous emission into free space as the cavity will provide a dominant channel of energy loss by imposing that its decay rate κ is much larger than g √ N . The evolution of the system is then described bẏ We notice first that the collapse operator a only couples to the total spin ladder operators S and S † which insures that the dynamics only takes place within the symmetric manifold with maximal spin length N /2. The second observation is that, owing to the quick loss rate κ which is much larger than the coherent light-matter exchange rate g √ N , the cavity field is permanently in a very low occupancy state. This allows us to apply the standard master equation procedure (see Appendix A1) where we assume the crude factorization ρ tot (t) = ρ(t) ⊗ |vac⟩ ⟨vac| at all times. The next step is in evaluating two-time correlation terms with the new time-dependent field operator F(t) = ga(t) (equivalent to the expression in Eq. (5)). The situation is now much simpler than in the free-space spontaneous emission case, as one can use the two-time correlations ⟨a(t ′ )a † (t ′′ )⟩ = e −κ(t ′ −t ′′ ) to find an effective damping rate g 2 /κ. The dynamics can then be reduced to the atomic subspace and incorporated in a master equationρ with H 0 = ω 0 [S z + (N /2)1] and the cavity-induced collective spontaneous emission The resulting dynamics is of course identical to the one described in Sec. II D showing the emergence of quick bursts of light as superradiant pulses but with a rate established by the cavity g 2 /κ. The effect is the spontaneous synchronization of dipoles by the cavity to give rise to the superradiant lasing regime. Here, the lossy cavity acts mainly as a communication bus to drive the synchronization between the atoms and the average photon number inside the cavity is typically well below one. For example, Ref. [25] reports the realization of a superradiant laser with an average photon number of less than 0.2 and a single-atom cooperativity C = 7.7 · 10 −3 , while the lasing regime achieved in Ref [26] uses a MHz linewidth optical clock transition as its gain medium and is envisioned to act as an active atomic clock insensitive to fluctuations in reference cavity length (a fundamental limitation in conventional lasers).

D. Further remarks
Composite systems made up of flat, standard mirrors and two-dimensional single emitter-thick regular arrays can find a multitude of applications in the direction of quantum technologies. When used in a standard cavity QED scenario, they can act as quick phase switchers with applications in precision spectroscopy of quantum network characterization [16] or in hybrid quantum optomechanical setups with enhanced photon-phonon couplings [101]. In the strong reflectivity limit, subradiant arrays have been interfaced with two-dimensional semiconductor monolayers to show enhanced quantum nonlinear optical properties [27]. The identification of dark collective states between the two mirrors has been proposed to allow for the preparation of Bell superpositions states between the two subradiant layers [8], in a double-mirror hybrid setup as discussed above. As such subradiant mirrors are of extremely small mass, their zero-point motion is much larger than that of traditional dielectric mirrors used in standard optomechanics: this opens new opportunities for quantum optomechanics at the single photon-phonon level [14,15].

VI. QUANTUM OPTICS WITH MOLECULES
A relatively recent novel direction of research with great promise in the direction of quantum technologies is the engineering of cavity-dressed materials. In the field of molecular polaritonics, embedding organic semiconductors in optical cavities allows for the design of novel materials with enhanced properties such as exciton and charge transport [29][30][31][32][33], superconductive behavior [102,103] or modified chemical reactivity [36-38, 40, 104] due to light-modified energy potential surfaces. This approach is also extended to study, control and design phase transitions in quantum materials by quantum light [105][106][107][108][109] etc.
Quantum emitters widely utilized in such experiments are of a much more complex nature than the two-level system approximation. We exemplify here how the previously introduced methods (master equation, quantum Langevin equations) and concepts can be extended to additionally include couplings between electronic transitions and vibrations or phonon modes. While this approach is in principle amenable to a wide range of solid-state emitters where phononic couplings play a crucial role, such as quantum dots or vacancy centers in diamond, we focus here on the specific case of molecules. In particular, organic molecules have emerged as tunable and efficient light-matter interfaces due to their relatively large dipole moments, wide range of transition frequencies and narrow linewidths at cryogenic temperatures [110]. As such, single-molecule impurities embedded in solid-state host matrices hold promise as single-photon sources [111][112][113][114] or nonlinear optical elements [115]. There also is interest in exploiting the strong inherent coupling between electronic transitions and the molecular vibrational degrees for the realization of quantum optomechanical effects at the single molecule level [116][117][118].
In this section, we will mainly follow the approach introduced in Ref. [52] by considering the time evolution of a polaron operator (i.e. an "effective", vibrationallydressed electronic dipole operator). We proceed with a first-principle derivation of the Holstein Hamiltonian for electron-vibron interactions and derive the standard Franck-Condon physics and absorption/emission properties of molecules. We then illustrate how one can circumvent a major downfall of molecular systems, i.e., they do not have closed transitions, by employing the Purcell effect: a cavity-dressed molecule can then behave as an ideal quantum emitter. Next we analyze molecular dimers as good candidates for the observation and exploration of controllable near-field couplings. Finally we introduce and analyze cooperative processes occurring between near-field coupled disordered molecules leading to what is known as the Förster resonance energy transfer.

A. Optical response of a single molecule
To include the electron-vibration (vibronic) coupling into our formalism, we start with a first-principle derivation for a single electron coupled to a single nuclear coordinate of reduced mass µ. Notice that, for the simplest case of a homonuclear diatomic molecule made of nuclei each with mass m, a single vibrational mode exists corresponding to relative motion with effective mass µ = m/2; the equilibria coordinates then correspond to the bond length in the ground and excited states. Under the Born-Oppenheimer approximation, the ground and excited state potentials are obtained by solving the Schrödinger equation for each given fixed internuclear distance. We assume that, along the nuclear coordinate, the equilibrium positions for ground (coordinate R g , state vector |g⟩) and excited (coordinate R e , state vector |e⟩) electronic orbitals are slightly shifted with respect to each other (as illustrated in Fig. 14). We can then write the total Hamiltonian H mol = H el ⊗ H vib as where the potential surfaces are described by the quantized position and momentum operator satisfying [R,P ] = i. The dynamics takes place now in an infinite-dimensional Hilbert space: in the electronic subspace the algebra is that of a spin 1/2 and is described by projectors σ † σ (into the excited state) and σσ † (into the ground state) while in the motional subspace one can for example describe the dynamics in terms of plane waves. However, a great simplification can be obtained by performing a harmonic expansion of the potential surfaces around the minima Finally, we can introduce small oscillations around the equilibriaQ =R − R g and subsequentlyR − R e =Q + R g − R e =:Q − R ge leading to The last term is a renormalization of the bare electronic transition frequency energy which will naturally go away when diagonalizing the Hamiltonian via the polaron transformation resulting in ω 0 as the natural electronic transition frequency. We can now rewrite the momentum and position operators in terms of bosonic operatorsQ = r zpm (b † + b),P = ip zpm (b † − b) by introducing the zero-point motion r zpm = 1/ √ 2µν and p zpm = µν/2. Reexpressing the terms above via the factor √ S = µνR ge r zpm (S is the Huang-Rhys factor) yields the Holstein-Hamiltonian [51] A molecular box illustration is shown in Fig. 15a as a minimal model for a molecule. One can bring this Hamiltonian into diagonal form via the polaron transformation U † pol = (D † ) σ † σ = |g⟩ ⟨g| + D † |e⟩ ⟨e| where the displacement is defined as 2 as generator. The position quadrature is defined analogously as The displacement operator creates a coherent state with amplitude √ S when applied to the vibrational ground state D |0⟩ vib = | √ S⟩ vib . Note that the vibrational creation and annihilation operators transform in the polaron picture as and furthermore the following commutation relations are fulfilled [p, The Hamiltonian in the diagonal basis then becomesH mol = ω 0 σ † σ + νb † b and has eigenvectors {|g, n⟩ , |e, n⟩} while the eigenvectors in the original basis can be obtained by an inverse polaron transform and are given by {|g, n⟩, D |e, n⟩}.
Franck-Condon physics -Let us now ask what is the optical response of the molecule to a weak laser excitation modelled via H ℓ = iη(σ † e −iω ℓ t −σe iω ℓ t ). In an interaction picture, reached after performing a transformation of the Hamiltonian with e iω ℓ σ † σt , one has Starting with the molecule in the absolute ground state |g, 0⟩, the result of the drive is to excite the electron to state |e⟩ while exciting motion on the nuclear coordinate to a coherent motional state. One can immediately compute the probability of absorption showing a Poissonian distribution f n S = e −S S n /n! of the occupancies of number states |n⟩ with coefficients known as Franck-Condon factors. In particular, the so-called zero-phonon line (transition corresponding to n = 0) is reduced by a factor e −S at zero temperature. The Huang-Rhys factor S should be interpreted as the average number of vibrational quanta created upon excitation or deexciation of the molecule. The probability of the emission process is similarly computed between |e, 0⟩ and final . . .   FIG. 15: (a) The dynamics of a molecule with a single nuclear coordinate can be seen as the coupling of an electronic transition operator σ to a single vibrational mode b with a strength given by the Huang-Rhys factor S. (b) For many vibrational modes, a molecule can be constructed by coupling an electronic transition operator σ to n vibrational modes with strengths S k , frequencies ν k and linewidths Γ k described by the (Lorentzian) spectral density function J(ω). (c) Energy eigenstates of the electron-vibron system denoted by |g, n⟩ and |e, n⟩. The zero phonon line is the |g, 0⟩ to |e, 0⟩ transition. For small S < 1 the width of the arrows indicate the strength of the optical transitions. For a molecule with many vibrational levels, spontaneous emission outside the zero phonon line at summed rate γ ′ can be larger than into it at rate γ00. (d) The enhancement of spontaneous emission via resonant coupling to an optical cavity can lead to a great modification of the strength of the zero phonon line thus turning a molecule into a closed quantum system. states |g, n⟩ and yields the same Poissonian distribution in n. This is a known result in molecular spectroscopy where the absorption and emission spectra are mirror images of each other. For small S < 1 the physics of absorption and emission in a molecule is illustrated in Fig. 15c where the width of the arrows indicate the strength of the optical transitions between eigenstates of the free Hamiltonian (without the Holstein interaction).
The Holstein Hamiltonian in Eq. (99) can be generalized to a more complex molecular box (is illustrated in Fig. 15b) involving n v independent nuclear coordinates with frequencies ν k and coupling strengths √ S k . The total Hamiltonian is then expressed as where the accumulated frequency shift isω 0 = ω 0 + nv k=1 S k ν k . Diagonalization of this Hamiltonian is straightforward as it involves a collective displacement for all nuclear coordinates The expectation value of the displacement operator ⟨D k ⟩ = ⟨D † k ⟩ = e −S k ⟨p 2 k ⟩ (assuming the vibrational mode to be in a thermal state) can be expressed as e −S k (n k +1/2) (expanding the displacement operator and making use of the fact that Gaussian states are characterized by their second-order moments and odd moments are vanishing) where the average thermal occupancy is given by 2n k = coth(βν k /2) − 1 with β = 1/(k B T ).
Branching ratio manipulation -The branching ratio quantifies the rate of emission into the zero-phonon line versus all other lines. For ideal quantum emitters this is automatically unity in the absence of coupling to any other channels of de-excitation. For an isolated molecule this is instead given by α = ⟨D⟩ 2 = k ⟨D k ⟩ 2 = e − k S k (at zero temperature) stemming from loss of excitation into other vibrational levels. Complex molecules have many vibrations: even for small Huang-Rhys factors the sum in the exponential gives a considerable reduction of the oscillator strength of the zero-phonon line. The Purcell effect (introduced in Sec. IV A), stemming from the resonant interaction of the zero-phonon line transition with a lossy confined optical resonance has been experimentally proven to turn molecular emitters into almost perfect two-level quantum emitters [42]. Assuming that the cavity linewidth is narrow compared to the vibrational frequency κ ≪ ν and only couples to the zero-phonon line transition, we can define by g 00 = g ⟨D⟩ the effective reduced coupling between the cavity and the zero-phonon line with a corresponding cooperativity C 00 = g 2 00 /(κγ). This will lead to an enlarged zero-phonon transition linewidthγ 00 = γ 00 (1 + C 00 ) as illustrated in Fig. 15 while other transitions are unperturbed and sum up to γ ′ = γ − γ 00 . The cavity-modified branching ratio can then be expressed as which reproduces the bare molecule result for C 00 = 0 and goes to unity in the limit that C 00 → ∞. Experimentally, values of C 00 close to 100 have been reached [42], showing that almost unity branching ratios can be achieved.
Quantum Langevin equations for polarons -We will now proceed by deriving a quantum Langevin equation for a polaron operatorσ = σD † (i.e., a vibrationally dressed electronic dipole operator). This will allow us to include all coherent and incoherent effects at the level of the equations of motion. The collective displacement operator D † can refer to a large set of molecular vibrational modes D † = nv k=1 D † k . The method can also be extended to include coupling to an (eventually infinite) set of phonon modes of the host environment [53,119]. We will assume a Brownian noise dissipation model (see Appendix A5) for the molecular vibrational modes. Starting with the Holstein Hamiltonian one can derive the following equations of motion for vibrationṡ In addition to the environment induced damping at rates 2Γ k and associated thermal input noise ξ k (see Appendix A5 for noise properties), vibrations are driven by their coupling to the electronic degree of freedom. For the electronic transition, the polaron operator can be computed fromσ =σD † + σḊ † where the full Langevin equation for σ is derived from the Holstein Hamiltonian and the Linbdlad form radiative emission with the transformation indicated by the relation in Eq. (63). We obtain thenσ The equation is derived under the assumption of weak driving η ≪ γ and therefore small occupancy of the electronic excited state ⟨σ † σ⟩ ≪ 1 and with special attention devoted to the rules for taking the time derivative of any exponential operator e A(t) with [Ȧ(t), A(t)] ̸ = 0. Moreover, some of the terms cancel (see Ref. [52] for full details) and an extremely simple equation of motion for the polaron operator is obtaineḋ Here we have combined the input fields affecting the electronic transition into the expression Σ in = η/ √ 2γ + σ in and the displaced input field (drive plus noise) is Σ in = Σ in D † . A solution of the bare electronic dipole operators can be obtained by a formal integration as a sum σ(t) = σ tr (t) + σ ss (t) between transient and steady state solutions To derive the optical response of the electronic degree of freedom in the presence of its coupling to vibrations one then has to evaluate two-time correlation functions of the displacement operator. This is factorizable as we consider that all vibrational modes are independent from each other, i.e., ⟨D(t)D † (t ′ )⟩ = k ⟨D k (t)D † k (t ′ )⟩. By expanding the displacement operators and applying the Isserlis' theorem (or Wick's probability theorem for multivariate normal distributions) one obtains where we made use of the fact that the momentum variance is stationary, i.e., ⟨p 2 k ⟩ = ⟨p k (t) 2 ⟩ = 1/2 +n k . This can be easily computed by using the properties for the Brownian noise model (see Appendix A5) with τ = t − t ′ . While expressions can be derived for any temperature of the bath, in the following we will, for the sake of simplicity, focus on the zero temperature case for all vibrational modes. This is typically a very good approximation as molecular vibrations typically lie in the range of THz frequencies and are therefore barely occupied as the condition ℏν k ≪ k B T is typically fulfilled. The displacement correlation function in this case simplifies to which we will make use of in the next subsection to derive analytical expressions for absorption and emission lines and linewidths.
Absorption and emission -The absorption and emission spectra of a driven molecule can be readily obtained from Eq. (108). Considering large times t ≫ 1/γ, the first term in Eq. (108) goes towards zero and the average dipole moment can be expressed as (taking the average over both electronic and vibrational degrees of freedom) with the laser detuning ∆ = ω ℓ − ω 0 . The population of the electronic excited state evolves according to ∂ t ⟨σ † σ⟩ = −2γ ⟨σ † σ⟩ + 2ηℜ ⟨σ(t)⟩ meaning that the steady state absorption profile can be directly computed from the expectation value of the coherence ⟨σ(t)⟩. We define S abs (∆) = ⟨σ † σ⟩ as the absorption profile for a varying laser frequency with expression To evaluate this integral, we expand the displacement correlation function from Eq. (111) in a power series The absorption spectrum can then be readily expressed as a sum over all indices {n k } = n 1 , . . . , n nv weighted by the Franck-Condon factors f S k where the denominator indicates a series of blue-shifted absorption sidebands with increased linewidths Γ {n k } = For the definition of the emission spectrum one can consider the transient dynamics of a molecule which is partially excited at time t = 0: ⟨σ † σ(0)⟩ = p 0 . The emission spectrum is then defined as the Fourier transform S em (ω) = 2ℜ ∞ 0 dτ ⟨T σ † (τ )σ(0) ⟩ e −iωτ (where T denotes time ordering) and can be computed as (normalized by p 0 ) The denominator indicates a series of red-shifted Stokes lines at frequencies ω = ω 0 − ν {n k } .

B. Near-field coupled molecules
Let us now consider two electronically coupled molecules. We start by discussing the vibronic dimer, a well-known model to describe the interplay between electronic and vibrational interactions in molecular aggregates which has been studied both theoretically [120][121][122] and experimentally [54,123]. As opposed to the pure electronic dimer discussed in Sec. II B, this model gives rise to a more complex energy landscape with states possessing both electronic and vibrational character. Effects showing quantum coherence in electronically excited dimers have been experimentally investigated in a combination of single molecule spectroscopy and quantum chemistry techniques [54]. We then move on to describe the dipole-dipole mediated unidirectional energy transfer process (FRET) between two molecules.
Vibronic dimer model -The vibronic dimer model consists of two molecules (indexed here by 1 and 2), each with a single electronic ground and excited state and each coupled to a single harmonic nuclear coordinate q 1/2 with Huang-Rhys factor S. The molecules are coupled to each other electronically (e.g. via dipole-dipole interaction) with strength Ω. Assuming that the excited state energies of the two molecules are identical and using the single excitation manifold with states |e⟩ 1 ⊗ |g⟩ 2 and |g⟩ 1 ⊗ |e⟩ 2 , the Hamiltonian of the dimer can be represented by where 1 is the 2×2 unity matrix. Introducing symmetrized coordinates Q ± = (q 1 ± q 2 )/ √ 2 and electronic operators σ ± = (σ 1 ± σ 2 )/ √ 2, the Hamiltonian transforms tõ from which one can see that only the coordinate Q − couples between electronic states in the collective basis. The coordinate Q + instead leads to a constant shift and can eventually be eliminated by a proper choice of the origin. The Hamiltonian in Eq. (118) can be easily diagonalized leading to the eigenvalues for the excited state potential energy surfaces The potential energy surfaces are plotted in Fig. 16 for various coupling strengths Ω. One can see that the coupling lifts the degeneracy of the potential surfaces at Q − = 0 and leads to a splitting of 2Ω. Depending on Ω, the lower energy surface either has two minima or just a single minimum. An important quantity is the reorganization energy ∆ϵ = 2Sν defined as the vertical energy separation between the potential curves at each minimum for zero coupling Ω = 0 [as illustrated in Fig. 16]. One finds that a double minimum of the lower energy surface located at Q ± − = ± S[1 − Ω 2 /(S 2 ν 2 )] exists when the level splitting 2Ω is smaller than the reorganisation energy |2Ω| < ∆ϵ.
Förster resonance energy transfer -Förster resonance energy transfer (FRET) is the main excitation transfer mechanism in photosynthetic light harvesting complexes, where energy is transfered with high efficiency from a light-harvesting antenna through a network of molecules (chromophores) to a reaction centre [43,44]. Although a lot of progress has been made in understanding and unravelling the reasons for the high efficiency of the photosynthetic energy transfer occuring in nature over the last years, there is still ongoing debate e.g. about the role of quantum coherence or entanglement in the energy transfer process [124,125]. In the sense of this tutorial, FRET can be seen as a truly cooperative process since it involves a complex interplay between electronic and nuclear degrees of freedom as well as both coherent and incoherent couplings. Generally speaking, the requirement for FRET is a spectral overlap between the emission profile of the donor molecule and the absorption profile of the acceptor, giving rise to a resonant transfer process typically accompanied by quick vibrational relaxation [126].
Here we exemplify the application of the QLEs formalism to derive a perturbative expression for the FRET rate between two near-field coupled molecules D and A with energy mismatch ∆ = ω D − ω A , first by considering a simple configuration where each molecule has a single vibrational mode b D and b A (with corresponding vibrational frequencies ν D and ν A and Huang-Rhys factors S D and S A ). We will assume an initially excited donor molecule which can undergo a resonant exchange with the acceptor's excited state vibrational manifold. In a standard scenario, a unidirectional process emerges as the acceptor quickly relaxes to its vibrational ground state Γ D ≫ Ω. We start with the equations of motion for the polaron-transformed dipole operators which shows coherent coupling induced by the dipole-dipole near-field exchange of virtual photonṡ To derive a perturbative expression for the energy transfer rate, we consider a scenario with an initial full excitation of D and no excitation of A. Then from above one can derive an equation of motion for the acceptor's populatioṅ which shows how the two-particle correlations ⟨σ † A σ D ⟩ appear as source terms for the acceptor population. The procedure to compute 2Ωℑ ⟨σ † A σ D ⟩ involves then formal integration of the equation of motion for the acceptor polaron operator (as detailed in the Appendix A6) under the assumption that Γ D , Γ A ≫ γ D , γ A . Finally, one arrives at an expression for the energy transfer rate showing propor-tionality to the donor population 2Ω ⟨σ † A σ D ⟩ = κ ET P D (t) where, for simplicity, we have set γ D = γ A . The denominator of the above expression asks that the resonance condition ω D −n D ν D = ω A +n A ν A is fulfilled, i.e., fluorescence lines of the donor |e D , 0⟩ → |g D , n D ⟩ have to overlap with absorption lines of the acceptor |g A , 0⟩ → |e A , n A ⟩, both of which are weighted by the respective Franck-Condon factors.
In the case of many vibrational modes n for donor and acceptor, we can generalize the result by writing general displacements D A = nv k=1 D k A and D D = nv k=1 D k D for all vibrational modes. The rate is computed by assuming multiple paths of energy transfer between the two molecules involving all vibrational modes. We assume an initially electronically excited state with no vibrations present |e D ; 0 1 , 0 2 ...0 nv ⟩ of molecule D and ground state without vibrations |g A ; 0 1 , 0 2 ...0 nv ⟩ for molecule A. The emission of molecule D leads it into state |g D ; m 1 , ...m k , ...m nv ⟩ and resonant interactions can occur with state |e A ; l 1 , ...l k ′ , ...l nv ⟩ of molecule A. Summing over all these processes leads to an analytical cumbersome expression for κ ET (listed in Appendix A6) seen simply as a generalization of the energy transfer rate from above. Such an analytical result is the discrete version of the well established integral formulation [126] describing the overlap between the emission spectrum of molecule D and absorption spectrum of molecule A as illustrated in Fig. 17a for a donor-acceptor pair with 8 vibrational modes and spectral density J(ω) = k 2S k ν 2 k Γ k /[Γ 2 k + (ω − ν k ) 2 ] (assumed to be identical for both molecules) stemming from the coupling of the molecular vibrations to some external phonon bath allowing for vibrational relaxation at rates Γ k . The process is unidirectional as shown in Fig. 17b as ω D − ω A dictates the direction of the energy flow.

C. Further remarks
We have covered here the electron-vibron interaction and described the optical response to a weak classical optical drive. We have however hinted that coupling to confined optical modes can lead to strong modifications on the matter side, such as the described model of a molecule turned into an efficient quantum emitter [42]. Indeed, many experimental and theoretical efforts have been recently targeted into the direction of cavity charge and energy transport enhancement [29][30][31][32][33], Förster resonance energy transfer (FRET) enhancement [34,35,[127][128][129], modified chemical reactivity [36-38, 40, 104], polariton dynamics [130][131][132], etc. On the theory side a large effort is aimed at translating standard cavity QED concepts (a) Overlap between acceptor emission spectrum and donor absorption spectrum leads to the FRET process. The results assume a given vibrational spectral energy J(ω) (assumed here to be identical for both molecules) shown in the inset, for 8 vibrational modes considered. (b) Numerical analysis of the energy transfer rate (for Ω = 30γ) shows a unidirectional flow from the higher energy molecule (donor) to the lower energy one (acceptor) optimized around a detuning corresponding to the overlap of the donor electronic state to a high density of vibrational states in the acceptor excited state manifold.
(as introduced in Sec. IV A) with two-level systems into the realm of molecules. To this end, theoretical efforts consider a generalized light-electronic-vibrations problem modeled as a Holstein-Tavis-Cummings Hamiltonian. Investigations aim at providing an understanding of the vibrationally induced cavity polariton asymmetry [132,133], vibrationally dressed polaritons [134], dark vibronic polaritons [135,136], developing a cavity Born-Oppenheimer theory [57,137] or deriving relevant simplified models for large scale numerics in the mesoscopic limit [138].
We also remark that molecules are natural quantum mechanical platforms interfacing electronic and vibrational degrees of freedom via an interaction that resembles the radiation pressure Hamiltonian (via a boson-spin replacement) utilized in quantum optomechanics; as vibration at terahertz frequencies are practically in the ground state even at room temperature and the strength of the coherent coupling can be comparable to the vibrational frequency, molecules could be a good platform for studying quantum state transfer between light and motion [118]. Moreover, molecules are good single photon emitters as they are characterized by large emission rates and they allow for highly efficient collection schemes, thus utilizable in applications for photon-photon entanglement generation and optical quantum computing [41].

VII. CONCLUSIONS
Quantum systems, comprised of many particles, inherently coupled via common reservoirs (the electromagnetic vacuum, optical resonators, optical waveguides, etc.) can exhibit cooperative behavior with some specific properties scaling more favorably than what would be expected from a collection of uncoupled particles. In order to analytically and numerically understand the cooperativity of light-matter platforms, we have detailed two alternative approaches based on either the time evolution of the density operator (master equation approach) or of system operators (stochastic quantum Langevin equations). Based on this toolbox, this tutorial introduced a set of equations successfully applied to chains, rings and arrays of quantum emitters or ensembles of molecules within the confined mode of optical resonators.
A first set of applications that we have tackled are based on effects stemming from the dipole-dipole interaction occurring in dense quantum emitter ensembles. This interaction allows for the hopping of excitations in one or two-dimensional regular arrays and for the engineering of energy bands with special properties such as band gaps or Dirac points. Thus, such systems are ideal platforms for the investigation of topological quantum photonics with built-in nonlinearities. Their subradiant properties allow for applications in quantum metrology, in the design of high-fidelity photon-storage platforms and in the generation of entangled Bell states of two or more quantum emitter arrays (hinting towards quantum networks of spatially distant subwavelength arrays). Their enhanced reflectivity properties render them useful as extremely light mechanical resonators for quantum nano-optomechanical applications or as metasurfaces with enhanced optical nonlinearities at the level of a single photon. Also, subradiant rings illuminated by incoherent light can provide a natural gain medium via supported waveguide resonances leading to thresholdless nanolasers.
A second set of applications that we have discussed are based on the enhancement of the photon-emitter coupling in cavity QED. We have demonstrated that the light-matter cooperativity can be greatly enhanced by addressing collective subradiant states of one or two-dimensional arrays which are characterized by narrow antiresonances in transmission. A hybrid cavity approach, where flat mirrors are replaced with frequency-dependent mirrors (such as subwavelength arrays or photonic crystal mirrors) can lead to narrow cavity resonances useful for example in resolved sideband optomechanics. We provided an input-output theory, valid outside Markovian regimes with wide applicability. In the same context of cavity QED, we provided a simple theory of superradiant lasers with potential applications in enhancing the stability of atomic clocks.
As molecular systems, quantum dots, vacancy centers in solid-state hosts are increasingly important for quantum technology applications, it is crucial to understand and to tailor the electron-vibration interaction. We provided here a first-principle derivation of the Holstein Hamiltonian for electron-vibron interactions in simple molecular systems and extended the QLEs approach to polaron physics. This theory leads to very simple analytical insight in the understanding of molecular spectral lines, their absorption and emission profiles and processes such as FRET migration of energy in donor-acceptor configurations.
We also made the analytical connection between the molecular dimer model, widely used in quantum chemistry for studying vibronic quantum coherence and the dipole-dipole interaction Hamiltonian.

Spontaneous emission of a single atom
Let us consider an interaction of the form H int = σ † F(t) + F † (t)σ. For spontaneous emission of a single atom, the term driving the emitter expresses as F(t) = k g k a k e i(ω0−ω k )t , where the coupling is given by Note that in our notation, the sum over k vectors implicitly contains a sum over the two orthogonal polarizations ϵ (1) k and ϵ (2) k . Inserting the interaction Hamiltonian into Eq. (A6), using the cyclic property of the trace (e.g., Tr B [F(t)ρ B (0)F † (t − s)] = Tr B [F † (t − s)F(t)ρ B (0)] = ⟨F † (t − s)F(t)⟩) and assuming a bath ρ B = |vac⟩ ⟨vac|, the only remaining terms areρ In consistence with the main text of the tutorial, we now denote the system density matrix simply by ρ. One is left with evaluating the term ⟨F(t)F † (t − s)⟩ = k |g k | 2 e i(ω0−ω k )s since ⟨a k a † k ′ ⟩ = δ kk ′ . Identifying the minimal volume in k-space which is occupied by two modes of orthogonal polarization and same k vector as (2π/ℓ) 3 , we can replace the sum by an integral and go into spherical coordinates Additionally, one still has to account for the sum over the two orthogonal polarizations which can be performed as For the integration over s, one can make use of the Sokhotski-Plemelj theorem where P denotes the Cauchy principal value. While the real part of the above expression gives rise to decay, the imaginary part corresponds to the Lamb shift which will eventually give rise to a small frequency renormalization of the atom. Neglecting this small correction and only considering the decay, one can finally arrive at a master equation in Lindblad formρ where γ = (ω 3 0 d 2 eg )/(6πc 3 ϵ 0 ) is the spontaneous emission rate.

The master equation for collective dynamics of closely spaced quantum-emitter ensembles
Let us now consider the case of N identical emitters located at positions R j mutually coupled to the electromagnetic vacuum described by the interaction H int = j (σ † j F j (t) + F † j (t)σ j ) with F j (t) = k g k a k e ikRj e −i(ω k −ω0)t , where we assume equal orientation of all transition dipoles d j eg = d eg . In this case, the master equation takes the forṁ where we denote the correlations as and used furthermore that ⟨F j (t)F † j ′ (t − s)⟩ = ⟨F j ′ (t − s)F † j (t)⟩ * . After summing over the two orthogonal polarizations we now again perform the continuum limit as in the previous paragraph and obtain for the correlation function Eq. (A13) where the function F (kR) is defined as the solid angle integral F (kR) = 1/(4π) dΩ k (1 − (d eg · k) 2 )/(d 2 eg k 2 ) · e ikR . This integral can be calculated as F (kR) = 1 4π 1 + (e d · ∇ R ) 2 k 2 π 0 dθ k e ikR cos θ k 2π 0 dϕ k = 1 + (e d · ∇ R ) 2 k 2 sin(kR) kR . (A15) Finally, resolving the integral over s by means of Eq. (A10) results in where All together, this yields the master equation Eq. (6) of the main text.

The master equation for cavity-induced Dicke superradiance
For N identical emitters placed within the volume of an optical resonator, the interaction Hamiltonian is given by H int = S † F(t) + F † (t)S with F(t) = ga(t). Using that the correlation of the field operator is given by (assuming the cavity field to be in the vacuum field due to its fast decay) the term appearing in the master equation iṡ

Damped cavity mode
The tunneling of photons between a confined cavity mode a to the outside continuum (modes on the left b(ω) and modes on the right c(ω)) can be phenomenologically described via the interaction Hamiltonian H int = a † F(t) + F † (t)a where F(t) = dω[V (ω)b(ω) + W (ω)c(ω)]e −i(ω−ωc)t with frequency-dependent rates V (ω) and W (ω). Evaluating the correlation (approximating the outside field as ρ B = ρ and subsequently the integral over s, gives rise to loss rates via left and right mirrors κ L = πV (ω c ) 2 , κ R = πW (ω c ) 2 and a master equation in Lindblad formρ where κ = κ L + κ R sums the loss rates via the left and right mirror.
A2. Functional dependence of F (kR), G(kR) The expression for F (kR) in Eq. (A15) can be brought into a more convenient form widely used in the literature by expressing the nabla operator in spherical coordinates Assuming the dipoles to be aligned along the z direction e d = e z with e z = cos θe R − sin θe θ , we can express Calculating the derivatives, one obtains F (kR) = (1 − cos 2 θ) sin(kR) kR + (1 − 3 cos 2 θ) cos(kR) (kR) 2 − sin(kR) (kR) 3 .

A3. Lattice sum of spherical waves
The sum over spherical waves of a 2D array with atoms periodically arranged at positions R j in the x-y plane can be expressed along the z direction by going into the continuum as where we denote the in-plane components by R ∥ = (x, y). The spherical wave can be expanded by means of a Weyl expansion as [140] S with in-plane momentum k ∥ = (k x , k y ) and k z = k 2 − |k ∥ | 2 . Poisson's summation formula allows one to turn the sum over the real-space lattice into a sum over the reciprocal lattice propagating), as illustrated in Fig. 18. We assume the scatterer to have a complex reflectivity r and transmittivity t. The two are connected as t = 1 + r and in the absence of absorption one has |t| 2 + |r| 2 = 1. The outgoing fields can be related to the incoming fields as A = rB + tC, (A30a) D = tB + rC. (A30b) Note that in the case where one has no incoming wave from the right (C = 0), this simply gives D = tB and A = rB.
The fields on the left can be connected to the fields on the right via For the two-mirror arrangement considered in Sec. V B, each mirror can be characterized by its polarizability ζ j , j = R, L, where the polarizability generally depends on frequency ζ j = ζ j (ω). The transmission and reflection coefficients for each mirror are given by t j = 1/(1 − iζ j ) and r j = iζ j /(1 − iζ j ), respectively and the transfer matrix of each mirror reads while the matrix for free-space propagation is given by T f = diag(e iθ , e −iθ ) where θ = ωℓ/c. The transfer matrix of the whole system is then obtained by multiplying the three individual matrices from which the cavity transmission coefficient can be obtained as t = 1/T 22 .

A5. Brownian noise model for vibrational relaxation
A particular case of dissipation arises in the case of a quantum harmonic oscillator with free energy H = ν(p 2 + q 2 )/2 interacting with a heat bath of (mutually independent) harmonic oscillators. In the Brownian motion model the decay 2Γ stems from the correlations of the stochastic input noise ξ which is only affecting the momentum quadrature. This model is a consequence of displacement-displacement interactions between the oscillator and an infinite surrounding bath H int = −q k α k Q k with coupling coefficients α k [46]. The correlations of the input noise ξ at a given bath temperature in the time domain read [141] where we defined the thermal noise spectrum S th (ω) While generally, the noise is not delta-correlated, one can approximate Eq. (A35) in the limit of large temperatures as (using that in this limitn + 1/2 ≈ k B T /(ℏν) with the average thermal occupancyn = [exp(ℏν/k B T ) − 1] −1 and coth(x) ≈ 1/x for x ≪ 1) where δ ′ (t − t ′ ) is the time derivative of the delta function. From the spectrum evaluated at ±ν one can obtain the cooling and heating contributions: S th (ν) = 4Γ(n + 1) and S th (−ν) = 4Γn. Generally, from these properties one can always estimate the damping rate as 2Γ = (S th (ν) − S th (−ν))/2 and the equilibrium thermal occupancy as n = S th (−ν)/[S th (ν) + S th (−ν)]. In the case of small temperatures k B T /(ℏν) ≪ 1, the thermal spectrum becomes very asymmetric S th (ω) = [4ωΓ/ν] θ(ω) as it vanishes for negative frequencies (we introduced the Heaviside function θ(ω)). Generally, one can proceed with a Fourier domain analysis of the steady state of the equations above which leads to p(ω) = ϵ(ω)ξ(ω) and q(ω) = iνp(ω)/ω with the mechanical susceptibility describing the response of the phonon bath One can also make use of the fact that the noise is always delta-correlated in Fourier domain ⟨ξ(ω)ξω ′ )⟩ = S th (ω)δ(ω+ω ′ ). One can then proceed by computing the time-domain correlations (using that ϵ(−ω) = ϵ * (ω)) Generally, these expressions can be solved by a contour integral in the complex plane. We note that the integral of the momentum correlation function diverges and generally requires the introduction of a cutoff frequency Λ > ν to keep it finite. However, for Γ ≪ ν, the expression for |ϵ(ω)| 2 S th (ω) is sharply peaked around ±ν. In this case, one can expand around the poles ω = ν + δ and ω = ν − δ with |δ| ≪ |ν| and only keep leading order terms in δ. The integrals in Eq. (A39) can then be approximated as the Fourier transforms of Lorentzian lineshapes which is much easier to solve and results in (for t > t ′ ) with τ = t − t ′ .

A6. Derivation of the FRET rate
Let us evaluate the term 2Ωℑ ⟨σ † A σ D ⟩ which is the one that give rise to energy transfer between the two molecules. Formal integration of the equation of motion for the acceptor gives Due to the quick vibrational relaxation, we can make a great simplification by neglecting the backaction of the acceptor's dipole moment onto the donor and assume free evolution for the donoṙ Integrating this, one can calculate the term 2Ωℑ ⟨σ † A σ D ⟩ which requires the correlation function (assuming t ≥ t ′ ≥ 0 and t ≫ 1/Γ D ) where the four-operator correlation function can be reduced to a two-operator correlation function by commuting D † D (0) with D D (t) and D † D (t ′ ) and using that D D (0)D † D (0) = 1. Under the assumption of Γ D , Γ A ≫ γ D , γ A , one can finally arrive at an expression for the energy transfer rate showing proportionality to the donor population 2Ω ⟨σ † A σ D ⟩ = κ ET P D (t) (assuming γ D = γ A ):