Molecular polaritonics in dense mesoscopic disordered ensembles

We study the dependence of the vacuum Rabi splitting (VRS) on frequency disorder, vibrations, near-field effects and density in molecular polaritonics. In the mesoscopic limit, static frequency disorder alone can already introduce a loss mechanism from polaritonic states into a dark state reservoir, which we quantitatively describe, providing an analytical scaling of the VRS with the level of disorder. Disorder additionally can split a molecular ensemble into donor-type and acceptor-type molecules and the combination of vibronic coupling, dipole-dipole interactions and vibrational relaxation induces an incoherent FRET (F\"{o}rster resonance energy transfer) migration of excitations within the collective molecular state. This is equivalent to a dissipative disorder and has the effect of saturating and even reducing the VRS in the mesoscopic, high-density limit. Overall, this analysis allows to quantify the crucial role played by dark states in cavity quantum electrodynamics with mesoscopic, disordered ensembles.


I. INTRODUCTION
The strength of light-matter coherent exchanges is enhanced when confined light modes, such as provided by optical cavities, are utilized. For N ideal two-level quantum emitters, equally coupled to a cavity mode, a collective enhancement proportional to N 1/2 can be obtained [1]. This is evident in the scaling of the collective vacuum Rabi splitting (VRS) in cavity quantum electrodynamics (cQED) [2][3][4]. In the particular case where more complex emitters, such as organic molecules (J-Aggregates, dye molecules, etc.), are collectively coupled to optical or plasmonic resonators, these standard results of cQED have been extensively invoked to describe the collective Rayleigh scattering loss from a cavity [5], the modification of energy transfer and transport [6][7][8][9][10][11], charge transport [12][13][14][15] or chemical reactions in the presence of strong light-matter interactions [16][17][18][19]. However, molecular polaritonics is characterized by emitters with large inhomogenous broadening, coupled to local vibrational baths and with strong near-field interactions, in which case analytical approaches are typically limited to only a few molecules and often with only one vibrational mode [20][21][22][23]. The numerical complexity of treating many electronic and vibrational degrees of freedom renders such problems hard to solve even with extensive simulations [24][25][26].
We propose here a fully analytical approach which allows to quantify the effect of disorder on light-matter interactions in the strong coupling regime. In a first step, we introduce the formalism for the case of pure two-level systems (involving electronic transitions only) with general applicability to cQED with atoms, quantum dots, superconducting qubits, etc. [27][28][29]. In a second step, we exemplify the application of this formalism to more complex systems involving electron-phonon interactions and in particular to molecular polaritonics [18].
The two main conceptual ingredients of our approach consist in the move to a collective basis for N emitters and Figure 1. a) Cavity-enclosed dense, disordered molecular ensemble, with electronic transitions subject to vibronic and near-field couplings is naturally split (owing to frequency disorder) between donor-like (red) and acceptor-like (blue) molecules. b) Schematics of processes leading to the incoherent FRET migration of excitations. The near-field coupling (with distance dependent strength Ω(r)), followed by rapid vibrational relaxation within the vibrational manifold of the excited electronic acceptor state, leads to a unidirectional flow of energy. the occurrence of a natural averaging in the mesoscopic limit (as opposed to averaging over many realizations [30][31][32][33]). As widely acknowledged, polaritons are formed by one bright superposition state hybridized with light while the rest of N − 1 dark states are only indirectly coupled owing to disorder [20,[34][35][36][37][38]. We take an open system dynamics approach to derive an analytical rate for the irreversible loss of energy from polaritonic states into the dark state manifold accompanied by a degradation of the VRS. In the bare basis, we elucidate the reduction of the VRS by showing that particles which are too far detuned or too lossy can fall out of the macroscopic polaritonic superposition.
We then apply our formalism to molecular polaritonics where the interplay between static disorder, near-field couplings and vibrational relaxation leads to a FRET process characterized by incoherent transfer of excitations from energetically higher donor-type to lower frequency, acceptor-type molecules (see Fig. 1). We map this problem into an incoherent dynamics in Lindblad form describing migration of excitation at rates analytically computable and derive the scaling law for the VRS with density applying the open system dynamics previously derived for pure two-level systems.
The paper is structured as follows: we introduce the Tavis-Holstein-Cummings model for N molecules each with two electronic and n vibronic degrees of freedom coupled to a confined cavity mode in Sec. II. We then proceed by analyzing the cavity transmission in the presence of frequency disorder in Sec. III and show how a mesoscopic average leads to a decay of polaritons into the dark state manifold. To characterize the degree of participation of quantum emitters to the collective strong coupling condition, we introduce a measure of macroscopicity of quantum superpositions reaching value N for perfect superpositions and unity for complete mixtures. The effects of dipole-dipole couplings together with vibrational relaxation are taken into account in Sec. IV and the elimination of the dark state reservoir is revisited, this time including the process of incoherent excitation migration within the molecular ensemble. Finally, we conclude and present an outlook in Sec. V.

II. MODEL
We consider N molecules indexed by j = 1, ...N with electronic states |g j and |e j (lowering operator σ j = |g j e| j ) separated by energy splittings ω j ( = 1) inhomogeneously distributed around ω with a distribution function p(δ) normalized to unity ∞ −∞ p(δ)dδ = 1. In particular we choose p(δ) = (1/ √ 2πw 2 )e −δ 2 /(2w 2 ) . We write each molecule frequency splitting as ω +δ j where the average around the central frequency vanishes δ j cl = 0 while the variance is δ 2 j cl = w 2 . The molecules are randomly spatially distributed within a volume V at positions r j . Each molecule exhibits a number n of nuclear coordinates each with frequency ν k (with k = 1...n) with harmonic motion described by the annihilation operators b jk such that b jk , b † jk = 1. The vibronic couplings are modeled as Holstein terms with Huang-Rhys factors λ 2 k stemming from a difference in the equilibrium positions of the ground and excited electronic potential landscapes.
For high densities, the near-field dipole-dipole interactions at rates Ω jj are dependent on the separation (with a standard |r j − r j | −3 dependence) and relative orientation of transition dipoles. The dipole-dipole Hamiltonian is and describes an excitation transfer via a virtual photon exchange. The free Hamiltonian is (see Ref. [23]) and adds to the vibronic coupling Hamiltonian [39] The vibronic coupling is obtained as a harmonic approximation of a Morse potential surface by expanding the electronic potential landscapes around their minima: the difference between the minima in the ground and excited state leads then to the Huang-Rhys factors λ 2 k . Such a model is widely employed [20,21,40,41] especially for molecules in condensed matter environments, as fast vibrational relaxation insures that states with more than one vibrational excitation are never reached.
The cavity mode is described by bosonic operator a at frequency ω c coupled with g j (r j ) ≡ g j to each molecule. The Tavis-Cummings Hamiltonian is then This is a simplification of the Dicke model when neglecting counter-rotating terms such as a † σ † j . While some current experiments operate on the brink of the ultrastrong coupling regime [42][43][44], polariton dynamics is well reproduced within this approximation.
We then proceed by writing the master equation of the system where the dissipative dynamics is included in the Lindblad part. For a collapse operator O with rate γ O the Lindblad term applied to a density operator ρ is All channels of dissipation are then modeled as standard Lindblad superoperators with collapse operators a, σ j , b jk and loss rates κ, γ, Γ k .

III. EFFECTS OF DISORDER
We will first show the effect of frequency disorder as the occurrence of dark state resonances in the (linear) spectral response of the cavity when driven with an external (weak) laser source. The pump is modelled via the following Hamiltonian where the pump frequency is ω and the weak drive amplitude is η. The equations of motion for the averages α = a and β j = σ j then reaḋ In a more compact form we can writė where the vector of amplitudes is v = (β 1 , . . . , β N , α) , driving is included also in vector form as v d = (0, . . . , 0, η) and the drift matrix is explicitly given in Appendix A.

A. Steady state cavity transmission
In steady state, the equations above lead to the normalized cavity amplitude transmission t = κ a /η expressed as valid also for positioning and orientational disorder with randomized couplings g j . The effect of orientational disorder is a trivial renormalization of the collective coupling from g √ N to g N /2 (for a completely random orientation of the molecular dipoles, as discussed in Appendix D).
In the following we restrict the discussion to the case of identical couplings g j = g (for all j). For a given realization of disorder, Fig. 2a shows two polaritonic peaks at ±g √ N obtained by the hybridization of a symmetric collective state to the cavity field. Non-zero disorder introduces couplings to N − 1 orthogonal asymmetric states visible in the cavity transmission as unequal height peaks between the polaritons. In the mesoscopic limit (see Fig. 2b), a natural averaging occurs which leads to a smoothing out of the additional peaks. Also, the polariton's height is decreased while their width is increased suggesting a loss mechanism which we will quantitatively address in the following in a transformed bright-dark basis.

B. Bright-dark state dynamics
We start with w = 0 and 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 generally obtainable by a Gram-Schmidt algorithm that leads to all vectors orthogonal to the bright state one and to each other. However, for the simplest case g j = g a straightforward choice of coefficients is indicated by a discrete Fourier transform We index the dark state manifold for k = 1, . . . , N −1 and note that for k = N we recover the bright stateD N =B. The equations of motion for averages B = B , D = D and α become (in a frame rotating at the central emitter frequency ω) with δ = ω c − ω and the couplings are defined as Fourier transforms of the disorder distribution For k = N , the equations above indicate that the bright state is the only one coupled to the cavity mode with the expected collective rate g N . However, disorder induces couplings to the whole manifold of dark states and within the dark manifold as well.
We can more compactly write the equations above as where the average of the disorder distribution (expected to vanish in the mesoscopic limit) isδ = ∆ kk = N j=1 δ j /N . This matrix can be put in a more convenient form due to the structure of the terms ∆ lN . Considering the relations ∆ N l = ∆ N −lN and ∆ 1N = ∆ N N −1 , we can rewrite the matrix as with p = (0, . . . , 0, g N ) . Notice that here we have separated between the cavity and matter states by introducing the reduced matrix of dimensions N × N (referring to the matter part) which assumes the following form (g) VRS degradation from g √ N to g √ N − 1 as a particle is removed from the polaritonic superposition by increasing its detuning from the cavity resonance.
The eigenvectors of the cyclic matrix M red are given by

C. Elimination of the dark reservoir
The procedure we will employ roughly follows the illustration in Fig. 2c showing first the identification of a bright state and then the elimination of the dark reservoir resulting in an effective unidirectional loss of energy from the polaritonic states. The elimination of the dark state manifold can be done 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 obtain (Appendix C)Ḃ (17) and obtain a memory kernel describing a generally non-Markovian loss process. In the mesoscopic limit, one finds ). Now one can identify a Markovian limit as a particular case of wide frequency distributions w γ. The Markovian result is then simply reproduced by seeing that the kernel f (t−t ) naturally tends to a delta function. In such a case the treatment can be simplified by setting all derivatives to zero in Eqs. (12) to find the dark state amplitudes The matrix M has dimensions (N − 1) × (N − 1) and represents the part of the matrix M red referring to the dark states only Replacing the eliminated variables into the equation of motion for the bright mode we obtain the effective dissipative dynamicṡ where the effect of the reservoir is to induce an effective frequency shift δ dark and loss rate γ dark obtained as the real and imaginary parts, respectively, of the following expression In the mesoscopic limit of large N , one can further simplify the expression of the loss rate to find extremely simple scaling laws of the decay rate induced by the dark state manifold (see Appendix B for details) The analytical results are in excellent agreement with numerical simulations (see Fig. 2d) which also indicate that bothδ and δ dark vanish. From here we can deduce the dependence of the VRS on disorder which can be obtained by diagonalizing the dynamics in the reduced cavity-bright state subspace to lead to This is an important result later generalized to molecules to analytically quantify the effect of time-dependent disorder associated with continuous shifting of electronic resonances via vibronic driving. The VRS is illustrated in Fig. 2e as a function of increasing disorder. The initial increase of the polariton splitting occurs in the case κ > γ as the maximum value of 2g N is reached when γ + γ dark = κ. The complete degradation of the strong coupling condition occurs when the disorder level is of the order of the cavity photon loss. The validity of the Markovian approximation is graphically illustrated in Fig. 2f. In Appendix C, we perform a more in-depth analysis of the non-Markovian regime by means of the quantum Langevin equations approach.

D. Reduction of VRS owing to far detuned particles
Let us now provide further clarifications of the mechanism of VRS degradation owed to inhomogeneous broadening by reverting our analysis to the alternative bare basis approach. To this end we start with the standard scenario of N identical particles coupled equally to the cavity fieldβ and notice immediately that the equations can be described in terms of a single bright state such thaṫ This indicates that the problem is simply described by a single collective mode strongly coupled at g √ N coupling while the dark manifold is completely decoupled and does not play any role in the dynamics. Now instead we assume N − 1 particles identically and resonantly coupled to the cavity mode while an additional particle is detuned by δ. We can rewrite the equations above now in terms of the N − 1 bright statė The bright state is similarly defined as: A close inspection of the above equations shows that when δ is detuned from the cavity resonance, the corresponding VRS shows a drop from g √ N for δ = 0 to g √ N − 1 for δ κ (see Fig. 2g). For large w, this behavior indicates that, when disorder is strong more particles are likely to have frequencies very far from the cavity resonance, which finally leads to the degradation of the strong coupling condition as clearly illustrated in Fig. 2e.

E. Macroscopicity of mesoscopic quantum superposition states
As the VRS can be derived from a Hamiltonian formulation restricted to the single excitation subspace it is interesting to also investigate the connection between VRS, as observed for example in the cavity transmission and the properties of the quantum superposition state. To this end, we propose a simple measure of quantum macroscopicity that aims at describing the number of quantum emitters actively and identically participating in an extended quantum state. The measure is restricted to the single excitation subspace spanned by states where one individual emitter is excited by |j = |g, . . . , e j , . . . , g, 0 for j ∈ {1, . . . , N }. The cavity mode excitation is represented by |N + 1 = |g, . . . , g, 1 and the ground state is |0 = |g, . . . , g, 0 . The averaging in the reduced subspace is performed with the reduced density operator The projectors P 0 = 1 N +2 − |0 0|, P N +1 = 1 N +2 − |N + 1 N + 1| simply eliminate the parts of the density matrix containing the ground and the photonic state. Notice that can be computed simply as 1 − ρ 00 − ρ N +1,N +1 . This allows us to rewrite the macroscopicity with respect to

C(t)
W-state mixed state |1 Figure 3. Time dynamics of macroscopicity. Time dynamics of macroscopicity C(t) for N = 10 particles for three different initial collective states of the quantum emitter ensemble: perfect superposition |W , completely mixed state and single particle excitation (e.g., |1 ). a) Undriven cavity shows conservation of macroscopicity. b) Driven cavity without disorder and η = γ shows production of maximal macroscopicity in all cases. c) Driven cavity with disorder at the level w = 30γ. Other parameters are: g = 40γ, κ = 20γ. Partial reduction of macroscopicity is obtained as an effect of disorder.
components of the original density matrix ρ as Notice that the measure is tuned such that it equals N for a perfect W-state |W = j |j / √ N and it drops to unity for completely mixed states such as described by a density operator ρ = j |j j| /N . We illustrate the time dynamics of the introduced measure of macroscopicity in three distinct cases for three different initial states. We distinguish between an initial state with maximal macroscopicity (the W-state) and two states with single particle participation (mixed state versus single excitation state). The important result illustrated in Fig. 3b shows that by coherently driving the cavity mode one can create maximal macroscopicity independently on the initial state. Also, in the presence of disorder, the macroscopicity is diminished for all initial states by the same amount. According to the interpretation obtained in the VRS case, one can see that for this given realization of disorder two far-detuned particles fall out of the macroscopic superposition thus diminishing C by 2 [ Fig. 3c]. In conclusion, we find that cavity driving can create macroscopicity while disorder and strong donor behavior can destroy it. It will therefore be interesting to extend such a measure beyond the single excitation subspace and to pursue in the future an analysis of the connection between collective strong coupling and the macroscopicity of quantum superposition states.

IV. REDUCTION OF THE VRS IN MOLECULAR POLARITONICS
Light-matter interactions in molecular ensembles are strongly modified in the presence of electron-vibron coupling as well as by the incoherent dynamics of molecular vibrations. In addition, in standard experimental situations densities are very high meaning that near field effects such as dipole-dipole couplings can play an important role. We will provide here a semi-analytical approach incorporating the competition between static disorder, vibronic and dipole-dipole couplings together with vibrational relaxation which can lead to a migration of excitation from a higher energy molecule (donor) to a lower energy one (acceptor) (as illustrated in Fig. 1b). As such ensembles are typically subject to strong inhomogeneous broadening, an automatic separation into donor-like and acceptor like molecules will then take place. To quantify the emergent incoherent FRET migration behavior we provide a phenomenological model which allows analytical and numerical insight into the scaling of the VRS of molecular ensembles with density.

A. FRET migration of excitation
For two near-field coupled adjacent molecules j and j , each with a single vibrational mode b j and b j we perform a polaron transformation which leads to the following transformed operatorsσ j = Q † j σ j andσ j = Q † j σ j (with displacement operators defined as Q j = e λj (b † j −bj ) ). Under the assumption of low population of the excited electronic levels, the dipole-dipole interaction couples the quantum Langevin dynamics of the two moleculeṡ This coupled dynamics can be solved in perturbation theory, assuming that Ω jj is small compared to the vibrational relaxation rates. The solution indicates an effective, largely unidirectional, energy transfer at rate κ jj ET from molecule j to molecule j . 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 j ; 0 1 , 0 2 ...0 n of molecule j and ground state without vibrations |g j ; 0 1 , 0 2 ...0 n for molecule j . The emission of molecule j leads it into state which is the discrete version of the well established integral formulation [45] describing the overlap between the emission spectrum of molecule j and absorption spectrum of molecule j . This is illustrated in Fig. 4a for a donor-acceptor pair with 8 vibrational modes and spectral density J(ω) = k 2λ 2 k ν 2 k Γ k /[Γ 2 k + (ω − ν k ) 2 ] (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. 4b as δ j −δ j dictates the direction of the energy flow.

B. VRS scaling at high densities
The expression of the energy transfer rate in Eq. (32) greatly simplifies numerical simulations as it allows one to introduce an effective model of loss where for each pair of molecules j and j a collapse operator σ j σ † j with corresponding rate κ jj ET is introduced. For molecule j, summing over all the paths of incoherent migration modifies the inherent radiative rate γ to an increased one γ + N j =j κ jj ET . Owing to the random spatial positioning within the ensemble, the FRET migration leads to an effective disorder in dissipation rates. This is directly incorporated in Eqs. (8) by amending the diagonal elements of the evolution matrix Results are then possible for large systems (where a direct simulation of the evolution of the master equation is untractable) by a simple diagonalization of this matrix. The obtained scaling presented in Fig. 4c shows that the FRET mechanism can lead to a strong deviation from the standard one ubiquitous in cavity QED with √ N . Beyond numerical estimates, a fully analytical approach based on the formalism introduced in the previous section is possible allowing one to compute the effective polariton loss rateκ ET +γ dark stemming from the competition between static disorder and FRET migration. The average FRET dissipation ratē is performed over the whole ensemble. In addition, the previously definedγ dark is derived from the matrix of Fourier transformed detunings and FRET rates (Appendix F) and it suffers modifications from the purely frequency disordered case. The new expression, seen as an extension of Eq. (23) (see Appendix F for more details of the derivation) is An averaging over disorder and position is possible assuming homogeneous media showing a scaling ofκ ET with N 2 . This allows for a fit of the result in Fig. 4c showing that saturation stems from the strong increase ofκ ET with density squared.

C. Reduction of VRS owing to strong donors
A very simple analysis in terms of bright and dark states can then shed insight into this scaling by assuming the effect of large decay onto the VRS. Assuming N − 1 molecules with identical decay rates γ and a single lossier molecule with γ similar conclusions as in section III referring to disorder are obtained. We can again cast the equations aṡ When γ = γ the system's bright state leads to polaritons at roughly ±g √ N while with increasing γ the polariton frequencies decrease to ±g √ N − 1 (as illustrated in Fig. 4d). For the many strong donors case we illustrate in Fig. 4e the gradual degradation of the VRS in time, from g √ N to zero, when successively particles are turned from weak to strong donors. Following this interpretation, we added histograms in the inset of Fig. 4c showing the distribution of dissipation rates within the ensemble. This provides a qualitative means to count out the number of lossy donors that fall out of the collective strong coupling condition.

V. CONCLUSIONS AND OUTLOOK
We proposed here an analytical approach which allows to quantify the effect of disorder on light-matter interactions in the strong coupling regime and which is extendable to molecular ensembles characterized by vibronic effects as well as near field effects owing to the electromagnetic vacuum. Employing a phenomenological model that incorporates all these aspects in an effective incoherent FRET migration of energy, scalings of the VRS with increasing density have been obtained, showing a strong divergence from the standard √ N scaling in the absence of particle-particle interactions.
The models used throughout the paper, albeit of limited validity, are standard and widely employed. For example, while the Holstein model for electron-vibron coupling is limited to molecules with large vibrational relaxation (such that anharmonicity is not reached) it provides a proper description to light-molecule [21,23,46] and molecule-molecule interactions [23,47]. The Tavis-Cummings model has also been universally used to predict and explain effects such as cavity mediated energy transfer [8][9][10]23], energy and charge transport [6,7,[12][13][14][15], cavity chemistry [16][17][18][19] etc. However, as some experiments are showing effects brought on by the onset of the ultrastrong coupling regime (such as for example the work function of a material [48]), recent theoretical works suggest that such regime is challenging and interesting [49,50] and approaches are greatly interdisciplinary mixing aspects of quantum optics with quantum chemistry methods [51]. In the future, we will extend our formalism based on linear quantum Langevin equations to include counter-rotating terms in the light-matter interaction Hamiltonian.
Appendix B: Elimination of the dark reservoir in the Markovian limit The expression for the frequency shift and decay rate induced by the dark state reservoir is: Some more insight can be obtained by evaluating the norm of the coupling vector for the dark states m = (∆ 1N , . . . , ∆ N −1N ) . Here, we obtain which is the variance of the frequency distribution. In the case of a Gaussian distribution p(δ) = (1/ √ 2πw 2 )e −δ 2 /(2w 2 ) for the disorder we obtain Var(δ) = w 2 . We can use this result to rewrite m = wm, wherem is the normalized vector of m. This allows us to further evaluate the expression in Eq. B1 to where M = TDT † with the diagonal matrixD andĉ = T †m . SinceD −1 is diagonal, we can find the expression which is the weighted average over the eigenvaluesλ −1 j of the matrix M −1 times the variance of the distribution. Since the eigenvalues of M follow the form λ 1,...,N −1 =δ − iγ −λ 1,...,N −1 whereλ 1,...,N −1 ∈ R we obtain For large N whereδ → 0 we finally obtain the expression In the case that γ λ 1,...,N −1 which is given if γ w, we can obtain a solution for γ dark which is obtained from since N −1 j=1 |ĉ j | 2 = 1. Since the eigenvaluesλ 1,...,N −1 are equally distributed around zero as shown in Fig. 3a and the weights can be roughly approximated by |ĉ j | 2 ≈ 1/(N − 1) (due to the fact that |ĉ j | 2 = 1/(N − 1)), the frequency shift becomes δ dark → 0 for large N and γ dark = w 2 /γ. A rough approximation can be performed in the caseλ j γ for most j ∈ {1, . . . , N − 1}. Here we set again |ĉ j | 2 ≈ 1/(N − 1) and we obtain where we have assumed that the eigenvaluesλ j are linearly distributed from −qw to qw where q is an adjustment or fitting parameter. This is a rather rough approximation of the real eigenvalue distribution depicted in Fig. 5a. In reality the sorted eigenvalue distribution followsλ j = √ 2πw erf −1 (2j/N − 1) for j ∈ {1, . . . , N − 1} when N → ∞, which is the quantile function of the Gaussian distribution for the energy disorder. The best approximation is given for q = 2 resulting in γ dark = (π/4)w as shown in Fig. 5b while δ dark = 0.
where we have definedv = (D 1 , . . . ,D N −1 ) we obtain by injecting the steady state solution for the dark modeŝ into Eq. (C1b) the reduced equation of motion for the bright modė where we have defined the noise term emerging from the dark reservoir by ξ in (s). Following the same steps introduced in the previous subsection where we approximate |c j | 2 ≈ w 2 /(N − 1) and assumingλ j being linearly distributed between −2w and 2w we obtain for the convolution kernel which results iṅ Before taking the limit for large width w γ, we evaluate the noise correlation term for the dark states Here, we have used the definition∆ = M + iγ1. For our approximation, this equates to For very large width w γ we can make the approximation sin(2w( which describes delta correlated noise (Markovian noise) in this limit. With the correlation relation 2γ B in (t)B in (t ) = 2γδ(t − t ) for the noise of the bright mode and simplifying the term in Eq. (C5), we obtain for the equation of motion in the Markovian limit the expressioṅ The Markovian limit can be obtained much more directly and irrespective of the given shape of the disorder distribution in the case where γ w. By using the relation (γ/2) exp(−γ|t − s|) ≈ δ(t − s) for large γ we can rewrite while simultaneously we obtain for the noise correlation term This results in γ dark = Var(δ)/γ for Eq. (C10). In this regime in particular we have identified In the case that we have molecules with randomly oriented dipole moments the cavity coupling strength g j varies from molecule to molecule. The effect of this disorder manifest itself in the definition of the bright and dark modes where the bright mode is now given byB = 1/ N j=1 |g j | 2 N j=1 g * j σ j which results in the equation of motion for the cavityα whereg N = N j=1 |g j | 2 and δ = ω c − ω l . Taking the limit for large numbers of molecules N where we describe the random orientation with respect to the electric field of the cavity mode a by g j = g cos(θ j ) where θ j ∈ [0, π] we obtaiñ The dark modesD k = 1/ N j=1 |g j | 2 N j=1 d * kj σ j can be obtained by employing the Gram-Schmidt orthogonalization procedure starting with the vector 1/ N j=1 |g j | 2 (g * 1 , . . . , g * N ) containing the coefficients of the bright mode. This results inḊ where ∆ From the equations of motion we obtain the matrix for the dark states where δ D k = ∆ (g) kk . Assuming steady state for the dark manifold we obtain the reduced equations of motioṅ For simplicity, for the numerical simulations, we will assume parallel orientation of all dipoles. The dipole-dipole interaction then expresses as Ω jj = 3 2 γ/(k|r j − r j |) 3 with wavenumber k = 2π/λ 0 (λ 0 is the wavelength of the electronic transition) [52]. We will consider initial excitation of molecule j and assume molecule j to be in the electronic ground state initially. The equation of motion for the acceptor's population reads: We therefore have to evaluate the term 2Ω jj σ † j σ j which signals the energy transfer. Formal integration of the equation of motion for the acceptor gives The correlation σ † j (t)σ j (t) can then be expressed as (assuming free evolution of j) where we defined P j (0) = σ † j σ j (0) . We therefore have to evaluate the two correlation functions Q j (t )Q † j (t) and Q j (0)Q † j (t )Q j (t)Q † j (0) (we will assume identical Huang-Rhys factors for both molecules λ := λ j = λ j ): where to obtain the expression for the latter one, we commuted Q † j (0) with Q j (t) and Q † j (t ), respectively and assumed large times t 1/Γ. We can now evaluate the integral I(t): Due to the fast decay of terms containing Γ (we assume Γ γ), we can approximate the energy transfer rate as (n 2 = n 3 = 0) where we introduced the Poissonian coefficients s λ n = e −λ 2 λ 2n /n!.
In the case of many vibrational modes n for donor and acceptor, we can generalize the result by writing general displacements Q j = n k=1 Q k j and Q j = n k=1 Q k j for all vibrational modes. The equations of motion can then be expressed in the same form as in Eqs. (E1) We further more assume that different vibrational modes are independent of each other, i.e., we assume factorizability of all correlation functions, e.g.
where a in , σ in and ξ in,ij are the noise operators for the collapse operators a, σ and σ † i σ j , respectively. For low population, we can linearize the equations of motion for the molecule-cavity systeṁ which results in the cavity transmission Additionally, we obtain for the equations of motion in the bright dark basis the expressionṡ where∆ kk = (1/N ) N j=1 δ j − i j =j κ jj ET e −i2πj(k−k )/N . By performing the same steps as introduced in previous chapters where we have injected the solution for the dark modes at steady state into the bright mode, we finally obtaiṅ Considering a homogeneous distribution of molecules, in the term for the dipole-dipole interaction Ω jj = (3/2)γ/(k|r j − r j |) 3 we can set r j to zero, since each molecule witnesses the same surrounding environment. Additionally, by using the probability distribution p(r, δ) to find a molecule at a distance r with detuning δ we exchange the sum over j with the integration where R is the radius of a spherical volume V and r min is the minimal radius that follows from the volume V /N that each molecule occupies individually. Assuming that r and δ are independent we can rewrite p(r, δ) = p(r)p(δ), where p(r) = 4πr 2 /V and p(δ) = (1/ √ 2πw 2 )e −(δ/ √ 2w 2 ) 2 . This allows us to obtain jj κ jj ET ≈ ∞ n1 ∞ n2 s λ n1 s λ n2 (n 1 + n 2 )ΓN R rmin drp(r) 3γ 2(kr) 3 2 ∞ −∞ dδ p(δ) [(n 1 + n 2 )Γ] 2 + [δ j − δ j − (n 1 + n 2 )ν] 2 .