When polarons meet polaritons: Exciton-vibration interactions in organic molecules strongly coupled to confined light fields
aa r X i v : . [ c ond - m a t . m e s - h a ll ] N ov When polarons meet polaritons: Exciton-vibration interactions in organic moleculesstrongly coupled to confined light fields
Ning Wu, ∗ Johannes Feist, † and Francisco J. Garcia-Vidal
1, 2, ‡ Departamento de F´ısica Te´orica de la Materia Condensada and Condensed Matter Physics Center (IFIMAC),Universidad Aut´onoma de Madrid, E-28049 Madrid, Spain Donostia International Physics Center (DIPC), E-20018 Donostia/San Sebasti´an, Spain
We present a microscopic semi-analytical theory for the description of organic molecules interact-ing strongly with a cavity mode. Exciton-vibration coupling within the molecule and exciton-cavityinteraction are treated on an equal footing by employing a temperature-dependent variational ap-proach. The interplay between strong exciton-vibration coupling and strong exciton-cavity couplinggives rise to a hybrid ground state, which we refer to as the lower polaron polariton. Explicit ex-pressions for the ground-state wave function, the zero-temperature quasiparticle weight of the lowerpolaron polariton, the photoluminescence line strength, and the mean number of vibrational quantaare obtained in terms of the optimal variational parameters. The dependence of these quantitiesupon the exciton-cavity coupling strength reveals that strong cavity coupling leads to an enhancedvibrational dressing of the cavity mode, and at the same time a vibrational decoupling of the darkexcitons, which in turn results in a lower polaron polariton resembling a single-mode dressed barelower polariton in the strong-coupling regime. Thermal effects on several observables are brieflydiscussed.
PACS numbers: 71.38.-k, 71.35.-y, 71.36.+c, 81.05.Fb
I. INTRODUCTION
Over the last couple of decades, there has beena renewed interest in organic molecular materials be-cause of their high relevance to organic light-emittingdiodes [1, 2], organic lasers [3], organic solar cells [4],organic field-effect transistors [5], and natural/artificiallight-harvesting systems [6, 7]. Organic materials arealso ideal systems to achieve strong coupling with con-fined light fields due to their large dipole moments andpossible high molecular densities. The strong-couplingregime is entered when the coherent energy exchange be-tween emitters and light modes becomes faster than de-cay and decoherence processes in either constituent. Thisleads to the formation of two polariton modes, i.e., hy-brid eigenstates that have mixed light-matter character,separated by the Rabi splitting. Strong coupling of or-ganic molecules has been studied in a wide variety of pho-tonic systems, among them dielectric microcavities [8–13], metallic microcavities [14, 15], plasmonic modes onflat [16–18] and holey surfaces [19, 20], and nanoparti-cle arrays supporting surface lattice resonances [21, 22].Additionally, the strong field confinement in plasmonicsystems also allows strong coupling with localized sur-face plasmon resonances [23–25], even down to the single-molecule level [26]. The Rabi splitting in these vastlydifferent systems is all quite similar, with typical val-ues of hundreds of meV and reaching up to more than ∗ [email protected]; Present Address : Center for QuantumTechnology Research, School of Physics, Beijing Institute ofTechnology, Bejing 100081, China † [email protected] ‡ [email protected] II. MODEL AND METHODOLOGYA. Hamiltonian
A typical organic microcavity setup consists of layer-structured organic materials sandwiched between two di-electric mirrors that form the microcavity [46]. Mostrecently, strongly coupled organic microcavities withsingle/few molecules have been realized experimentally,where the volume of the microcavity can be scaledto less than 40 cubic nanometers by employing ananoparticle-on-mirror geometry [26]. For simplicity, weconsider an organic microcavity composed of a single one-dimensional organic molecule located in a single-modecavity. The single organic molecule is assumed to consistof N chromophores. Such a system is described by theHamiltonian H = H mat + H c + H e − c ,H mat = H e + H v + H e − v ,H e = X n ε n a † n a n + X n = m J nm a † n a m , J nm = J mn ,H v = X n ω n b † n b n , H c = ω c c † c,H e − v = X n λ n ω n a † n a n ( b n + b † n ) ,H e − c = g X n ( a † n c + c † a n ) . (1)The material part H mat of H is the Holstein Hamilto-nian that describes the organic molecule with intramolec-ular vibrations. In principle, the molecule also interactswith the continuous phonon modes from its surroundingenvironment. Usually, such molecule-phonon coupling isweak, and we henceforth neglect the continuous phononbath for the sake of simplicity. Here, a † n creates an exci-ton state | n i on site n with on-site energy ε n , and J mn isthe hopping matrix element between two distinct sites m and n . The intra-molecular vibrational mode on site n with frequency ω n is created by the boson creation opera-tor b † n . H e − v is the linear exciton-vibration coupling withstrength measured by the Huang-Rhys factor λ n . Theradiation part is described by H c with photon creationoperator c † and cavity frequency ω c . The last term in1 represents the uniform exciton-cavity interaction withinteraction strength g . Here, we have employed the ro-tating wave approximation (RWA) such that no counter-rotating terms are present and H e − c conserves the totalnumber of excitations. This approximation is valid pro-vided the ultrastrong-coupling regime is not reached, i.e.,as long as the Rabi splitting is significantly smaller thanthe excitation energies ε n , ω c (see Ref. [32] for a discus-sion of possible effects caused by the breakdown of theRWA).In this work, we will consider one-dimensionalmolecules with uniform on-site energies and nearest-neighbor electronic couplings, i.e., we set ε n = ε and J mn = Jδ m,n ± . Important examples include linear J-aggregates [35] and the light-harvesting complex II witha ring-like structure [36, 37]. We assume periodic bound-ary conditions in the former case. We have checked thattypical amounts of static disorder and inhomogeneousbroadening do not significantly affect the results pre-sented here. For simplicity, the vibrational modes aremodeled by Einstein oscillators with a single frequency ω n = ω and uniform exciton-vibration coupling λ n = λ .In the following, we will restrict ourselves to the single-excitation subspace with P n a † n a n + c † c = 1, such that(within the RWA) we can truncate the number of cav-ity photons to be, at most, one. In turn, we can write a † n = | n ih vac | and c † = | c ih vac | , where | vac i is the com-mon vacuum of all the annihilation operators appearingin 1, hence an eigenstate of H with vanishing energy.Note that in the absence of the vibrations and phonons,the excitonic and cavity part of H resembles an interact-ing central spin model with spins 1 / a n = 1 √ N X k e ikn a k , b n = 1 √ N X q e iqn b q . (2)We see that only the exciton state with zero momentum, | k = 0 i = a † | vac i , couples to the cavity field, so that thetotal crystal momentum P tot = X k ka † k a k + X q qb † q b q (3)is still a good quantum number. B. The generalized Merrifield transformation
In this work, in order to treat the exciton-vibrationcoupling and exciton-cavity coupling at finite tempera-tures on an equal footing, we employ an extended varia-tional Merrifield transformation [38] determined by min-imizing the Bogoliubov upper bound for the free energy.As demonstrated for the Holstein model in Ref. [40],and more recently in Ref. [48], these kinds of variationalcanonical transformation methods could offer an accu-rate description of both static properties (e.g., the groundstate, the optical spectra, etc.) and dynamical properties(e.g., the exciton transport mechanisms) from intermedi-ate to strong exciton-vibration coupling regimes.When the cavity field is introduced, it couples onlyto the bright exciton, and there is no direct interactionbetween the cavity mode and the vibrations (though ex-plicit cavity-vibration coupling has been considered inRefs. [49–52]). However, as we will show below, in theframework of the canonical transformation, the interplayof the light-matter interaction with the exciton-vibrationcoupling will induce an effective cavity-vibration coupling in the residue interaction in the Merrifield frame. It isstraightforward to extend the present method to nonuni-form or disordered systems.To obtain an optimal zero-order representation of theHamiltonian (1) for a wide range of parameters, we pro-pose the following generalized Merrifield transformation˜ H = e S He −S , S = − X n a † n a n B n − c † cB c , (4)with vibrational operators B n = X l f l ( b n + l − b † n + l ) , B c = h X l ( b l − b † l ) . (5)The variational parameters { f l } and h are chosen to bereal and are determined self-consistently by minimizingthe free energy of the transformed system using Bogoli-ubov’s inequality [53], F ≤ F + hH i H (6)for a generic Hamiltonian H = H + H , where F and F are the free energies of H and H , respectively, and h ... i H represents the thermal average over the canonicalensemble defined by H .Physically, the coefficient f l quantifies the degree ofdressing of an exciton at site n by the vibrational modeat site n + l , while h measures the degree of dressingof the cavity photon by the vibrational mode on eachexcitonic site, though the cavity is not directly coupled tothe vibrations. The usual small polaron transformationfor the Holstein model corresponds to the case of f l = δ l λ and h = 0. By introducing the Fourier transformsof { f l } , ˜ f q = X n e iqn f n , (7)the extended Merrifield generator can be written in themomentum space as S = − √ N X kq a † k + q a k ˜ f q ( b q − b †− q ) − c † c √ N h ( b − b † ) , (8)which clearly converses the total crystal momentum P tot of the transformed states.Besides the circular symmetry, the exciton-vibrationsystem also holds inversion symmetry, which reduces thenumber of independent variational parameters from N to N + 1 ( N +12 ), i.e., { f , f = f N − , ..., f N − = f N +1 , f N } ( { f , f = f N − , ..., f N − = f N +12 } ) for even (odd) N [40]. Using the tilde to indicate the Merrifield frame,the transformed Hamiltonian can be separated in a con-ventional way as ˜ H = ˜ H S + ˜ V + H v , (9)where the system part reads˜ H S = X k =0 E k a † k a k + [ ε a † a + ˜ g √ N ( a † c + c † a ) + ˜ ω c c † c ] . (10)Here, E k = ε + ω ( X m f m − λf ) + 2 ˜ J cos k (11)is the vibrationally renormalized exciton dispersion andwe have introduced the renormalized parameters˜ J = J Θ , ˜ g = g Θ , ˜ ω c = ω c + N h ω , (12)with Θ = h e B c − B n i v = e − coth βω P l ( f l − h ) , Θ | n − n ′ | = h e B n − B n ′ i v = e − coth βω P l ( f l − n − f l − n ′ ) . (13)where h ... i v = Z v Tr v { e − β H v ... } is the thermal averagewith respect to the vibrational modes, with β = 1 /k B T the inverse temperature and Z v = Tr v e − β H v the vibra-tional partition function. It is clear from Eqs. (12) and(13) that the interaction with vibrations will decrease both the effective hopping integral J and the cavity cou-pling g , but will increase the effective cavity frequency ω c .The residue interaction part is of the form˜ V = X k ,k ( P k k + T k − k ) a † k a k + X k ( a † k cR k + c † a k R † k )+ ω √ N hc † c ( b + b † ) , (14)where P k k = JN X nm δ | m − n | , ( e B m − B n − Θ ) e − ik n + ik m ,R k ≡ g √ N X n e − ikn ( e B c − B n − Θ) ,T q ≡ ω √ N ( λ − ˜ f q )( b †− q + b q ) , (15)are operators of vibrational degrees of freedom and sat-isfy P k k = P † k k and T q = T †− q . We see that theextended Merrifield transformation leads to an effectivecavity-vibration interaction ω √ N hc † c ( b + b † ) with thezero-momentum vibrational mode.By examining the zero-order system Hamiltonian ˜ H S ,it is clear that only the bright state, the single-excitonstate with zero momentum | k = 0 i = a † | vac i , couples tothe cavity photon. The N − | k i = a † k | vac i ( k = 0) are themselves eigenstates of ˜ H S . The interactionbetween the bright exciton and the cavity mode resultsin two eigenmodes which bring ˜ H S into a diagonal form,˜ H S = X k =0 E k a † k a k + E U a † U a U + E D a † D a D , (16) where a † U = Ca † − Sc † and a † D = Sa † + Cc † are the cre-ation operators of two new quasiparticles, and the corre-sponding eigenenergies are E U / D = E + ˜ ω c ± s N ˜ g + (cid:18) E − ˜ ω c (cid:19) . (17)Here, the mixing coefficients C = cos θ and S = sin θ aredetermined by tan θ = 2˜ g √ N / (˜ ω c − E ) . (18)Although the two branches of eigenmodes resemble thelower and upper exciton polaritons [45], we have to keepin mind that these structures appear in the Merrifieldframe, and hence do not correspond to physical quasiex-citations. Actually, transforming back to the originalframe from the Merrifield frame will yield physical quasi-particles which are mixtures of excitonic, photonic, andvibrational degrees of freedom. In the following, we willrefer to the U and D quasiparticles as Merrifield polari-tons .To obtain the optimal zero-order Hamiltonian ˜ H =˜ H S + H v , we proceed by minimizing the Bogoliubov upperbound for the free energy of ˜ H at inverse temperature β , F B = − β ln Tre − β ˜H + h ˜V i , (19)where h ... i = Tr { ... e − β ˜H } / Tr { e − β ˜H } . By construction, h ˜ V i = 0, so that the second term in Eq. (19) vanishes.In the single-excitation subspace, the Bogoliubov boundcan be expressed in terms of single-particle eigenenergiesof ˜ H as F B = − β ln Z S + F v , (20)where Z S = X η = { k ( =0) , U , D } e − βE η , (21)is the partition function for ˜ H S , and F v = − β ln Z v is thefree energy of the free vibrational modes. As F v is notdependent on the variational parameters, we only needto minimize the first term of Eq. (20). To this end, thesaddle-point conditions { ∂F B /∂f n = 0 } and ∂F B /∂h =0 should be solved self-consistently. Two forms of theresultant saddle-point equations are listed in AppendixA. We emphasize that the such obtained F B gives anupper bound for the intrinsic free energy of the system.It is convenient to write the residue interaction ˜ V inthe basis {| η i = a † η | vac i} ( η = k ( = 0) , U , D),˜ V = X η η | η ih η | ˜ V η η , ˜ V η η = x η x η [ P k ( η ) ,k ( η ) + T k ( η ) − k ( η ) ]+ x η y η R k ( η ) + y η x η R † k ( η ) + y η y η ω h √ N ( b + b † ) , (22)where { x η } = { , ..., , C, S } , { y η } = { , ..., , − S, C } , (23)and { k ( η ) } = { k ( = 0) , , } . III. GROUND STATE: THE LOWERPOLARON-POLARITON
For the Holstein model without the cavity, one canintroduce the adiabaticity ratio γ = ω / | J | and thedimensionless exciton-vibration coupling strength α = γλ . Then, γ < >
1) defines the adiabatic (an-tiadiabatic) regime, and α > <
1) defines thestrong (weak) exciton-vibration coupling regime [54]. Be-sides the method employed in the present work, theground state of the Holstein model has been widelystudied by various analytical/numerical methods, includ-ing numerical diagonalization based on the two-particleapproximation [55–58], quantum Monte Carlo simula-tion [59], density matrix renormalization-group tech-nique [60], exact-diagonalization method [61, 62], andvariational ansatz [63].In the presence of the cavity and in the zero-temperature limit, the Bogoliubov bound to be mini-mized becomes the zero-order ground-state energy E D ,i.e., the eigenenergy of the lower Merrifield polariton, | D i = a † D | vac i . (24)In this case, it can be shown that the optimal varia-tional parameters are given by (see Appendix A) λN h = 1 S − ω ˜ g √ N CS , (25)˜ f N h = 1 − ω ˜ g √ N CS , (26)and λ ˜ f q = 1 − ˜ g √ Nω CS − − cos q ) ˜ Jω , (27)for q = 0.In the absence of the vibrational modes, the groundstate of the exciton-cavity system is simply obtained bydiagonalizing the Hamiltonian H e + H c + H e − c , yieldingthe bare upper polariton and lower polariton, | φ UP i = C | i − S | c i , | φ LP i = S | i + C | c i , (28)where C = cos θ , S = sin θ with tan θ =2 g √ N / ( ω c − ε − J ). When the vibrational bath ispresent, the corresponding physical ground state can beobtained by transforming | D i back to the original frame, | ψ LPP i = e −S | D i = SN X n X k e − ikn e − √ N P q ( ˜ f q e − iqn b † q − ˜ f − q e iqn b q ) | k i + Ce − h √ N ( b † − b ) | c i . (29) We see that | ψ LPP i has a similar structure to the freeLP state given by Eq. (28), but includes the vibration-induced effects through the vibrational coherent states.The first term on the right-hand side of Eq. (29) mimicsa polaron state with amplitude S , while the second termcorresponds to the vibrational dressing of the cavity statewith amplitude C . Furthermore, unlike the bare LP statewhich only has a component of the bright exciton | i , theexciton-vibration coupling also mixes the excitonic darkstates | k i ( k = 0) into | ψ LPP i . For this reason, we referto the quasiparticle corresponding to | ψ LPP i as the lowerpolaron polariton (LPP).The above form of the LPP state can be comparedwith the following generalized Toyozawa ansatz [54, 64]: | ψ TA i = 1 √ N X n X k e − ikn Φ k e P q ( ξ q e − iqn b † q − ξ ∗ q e iqn b q ) | k i + 1 √ N Φ c e ξ c b † − ξ ∗ c b | c i , (30)which recovers | ψ LPP i in the case ofΦ k = S √ N , ξ q = − ˜ f q √ N , Φ c = C √ N , ξ c = − h √ N . (31)The Toyozawa ansatz (TA) is believed to provide ac-curate results for the ground-state wave function andground-state energy of the Holstein model [54, 63, 65].Since there are more variational parameters in | ψ TA i thanthose in | ψ LPP i , the ground-state energy E TA obtainedby TA is slightly lower than E D . However, we emphasizethat the Merrifield transformation based on minimizingthe Bogoliubov free energy also applies to finite temper-atures .From the variational principle, E D provides an upperbound for the true ground-state energy of the system.A conventional procedure for obtaining a lower approx-imated ground-state energy is to calculate the second-order energy correction in terms of the residue interac-tion ˜ V [40]. For small systems, we have checked numeri-cally that the second-order correction to E D gives a moreaccurate approximation for g √ N /ω ≪
1, but underesti-mates the true ground-state energy in the strong exciton-cavity coupling regime with relatively large g √ N /ω .This can be illustrated by studying a molecular dimerwith N = 2 chromophores, for which the ground-stateenergy can be obtained exactly by numerically diagonal-izing the Hamiltonian in a truncated vibrational spacewith P i =1 , b † i b i = M max vibrations.Figure 1 shows the calculated ground-state energy ofa molecule dimer by the exact numerical diagonalizationwith up to M max = 20 vibrations ( E exact ), the varia-tional Merrifield transformation without and with thesecond-order energy correction ( E D and E corr ; see Ap-pendix B), as well as the Toyozawa ansatz ( E TA ). Weset λ = 1, ω = 1 eV , and ε = ω c = 0 eV , namely,a cavity frequency resonant with the on-site excitonic (a) g √ /ω -2.2-2-1.8-1.6-1.4 G r ound s t a t e e n e r g y ( e V ) E exact E D E corr E TA0.06 0.08 0.1 g √ /ω -1.29-1.28-1.27 1.3 1.35 1.4 g √ /ω -2.1-2.05-2-1.95 J/ω = − / (i) (ii) (b) g √ /ω -4.95-4.9-4.85-4.8-4.75-4.7-4.65-4.6-4.55 G r ound s t a t e e n e r g y ( e V ) E exact E D E corr E TA0.06 0.08 0.1 g √ /ω -4.56-4.558-4.556 1.3 1.35 1.4 g √ /ω -4.94-4.92-4.9-4.88 (i) (ii) J/ω = − FIG. 1. The ground-state energy of a molecular dimer inter-acting with a single cavity mode as a function of the collectiveexciton-cavity coupling g √ N /ω for two different sets of exci-tonic coupling: (a) J/ω = − / J/ω = −
2. Resultsfrom the exact numerical diagonalization (red line), the varia-tional Merrifield transformation without and with the second-order energy correction (solid black and dotted black lines),and the Toyozawa ansatz (blue line) are presented. The in-sets in each figure show the magnifications in the (i) weak and(ii) strong exciton-cavity regimes. Other parameters: λ = 1, ω = 1 eV , ε = ω c = 0 eV , and M max = 20. transition. The results for two sets of nearest-neighborinteractions J/ω = − / J/ω = − E D from the Merrifield transforma-tion overestimates the ground state energy in the weakexciton-cavity regime g √ /ω ≪
1, while the second-order corrected energy E corr gives a more accurate one.However, E corr begins to show large deviation from theexact value E exact and underestimates the true ground- state energy when one enters the strong cavity couplingregime g √ /ω ∼
1. In contrast, both the zeroth-orderenergy E D and the Toyozawa variational energy E TA be-come closer to E exact in this regime. Since the first-orderenergy correction vanishes by construction, while thesecond-order energy correction is always negative, it isexpected that higher order corrections are needed to geta more accurate ground-state energy for large g √ N /ω .As we are mainly interested in the strong exciton-cavitycoupling regime, we will henceforth take E LPP ≈ E D (32)as an approximation of the ground-state energy. Cor-respondingly, we take the LPP wave function given byEq. (29) as an approximated ground state in order toobtain simple and intuitive analytical expressions for ob-servables discussed in the following.Figure 2(a) shows the ground-state energy E LPP of theLPP as a function of the dimensionless exciton-cavitycoupling strength g √ N/ω for two different excitoniccouplings J/ω = − / J/ω = −
2. Other molec-ular parameters are taken as N = 100, ω = 0 . eV , ε = 2 eV , and λ = 1. Note that J/ω = − / small-polaron limit, i.e., the strongexciton-vibration coupling antiadiabatic limit [54]. Thecavity frequency is set to be resonant with the groundstate energy in the absence of the cavity [which is cal-culated by using the Toyozawa ansatz given by Eq. (30)in the g → N = 100], i.e., ω c = E TA ( N =100 , g = 0) = 1 . eV and 1 . eV for J/ω = − / J/ω = −
2, respectively. As expected, the couplingbetween the bright exciton and the cavity mode leads tothe formation of the LPP state which lies below the purepolaron state of the molecule.In Fig. 2(b), we plot the evolution of ˜ f = P n f n andthe cavity dressing parameter N h with the exciton-cavitycoupling g √ N /ω . The decreasing of ˜ f with increasing g clearly indicates a reduced vibrational dressing of ex-citons in the strong-coupling regime. We note that thedecoupling of vibrational degrees of freedom from the ex-citons by strong cavity coupling has also been reportedby Spano and co-workers [29, 30] by using numerical di-agonalization of the Holstein Hamiltonian. It is intrigu-ing to note that for this resonant case, the cavity dress-ing parameter N h increases monotonically as g increases,which means that the cavity mode, even though it is not coupled to the vibrations directly , might become moredressed by the vibrations as the exciton-cavity couplingincreases.However, as can be seen from Eq. (29), the degree ofthe vibrational dressing of the cavity is measured by boththe amplitude C and the dressing parameter N h . Tothis end, we plot in Fig. 3 the absolute values of theamplitudes C and S ( C and S ) in the LPP | ψ LPP i (thebare LP | φ LP i ) as functions of g √ N /ω for both resonantand nonresonant cases. The behavior in the weak cavitycoupling region g √ N /ω ≪ (a) g √ N/ω E L PP ( e V ) J/ω = − / J/ω = − (b) g √ N/ω ˜ f ,J/ω = − / Nh,J/ω = − / f ,J/ω = − Nh,J/ω = − FIG. 2. (a) The ground-state energy E LPP of the LPP asa function of the collective exciton-cavity coupling g √ N /ω for two different excitonic couplings: J/ω = − / J/ω = −
2. (b) The evolution of the collective vibrationaldressing parameter ˜ f of the excitons and the cavity dressingparameter Nh with g √ N/ω . Other parameters: N = 100, ω = 0 . eV , ε = 2 eV , and λ = 1. In both cases, the cavityfrequency is set to be resonant with the ground-state energyof the molecule in the absence of the cavity. investigating the saddle-point equations (25) and (26).As g √ N /ω → + , we have N h ≈ λ θ (˜ ω c − E )1 + ω / [2(˜ ω c − E )] , (33)and ˜ f ≈ λ E − ˜ ω c ) θ ( E − ˜ ω c ) /ω , (34)where θ ( x ) is the Heaviside step function.For the resonant case, the polaron part and the dressedcavity part of | ψ LPP i are roughly equally weighted forweak exciton-cavity coupling [Fig. 3(a)] with | C | ≈ | S | ≈ / √
2. As g increases, the cavity dressing parameter N h increases monotonically, while the amplitude | C | has nodramatic change even up to the strong-coupling regime. g √ N /ω | C || S || C || S | Z LP P
N h ˜ f ω c =1.7864 eV ω c =2 eV ω c =1 eV (a)(b)(c) FIG. 3. Absolute values of the mixing coefficients C and S ( C and S ) in the LPP state | ψ LPP i (the bare LP state | φ LP i ),and the quasiparticle weight Z LPP as functions of g √ N/ω for (a) a resonant cavity mode with ω c = 1 . eV , (b) ω c =2 eV , (c) ω c = 1 eV . Also shown are the dressing parameters Nh and ˜ f . Other parameters: N = 100, ω = 0 . eV , ε = 2 eV , J/ω = − /
4, and λ = 1. For the non-resonant case with the cavity frequency ω c relatively large, so that the condition ˜ ω c − E > N h withincreasing g . However, the increase of the amplitude | C | from 0 + to its saturated value in the strong-coupling re-gion might still indicate an enhanced dressing of the cav-ity field. For the nonresonant case with a red-detunedcavity frequency, both the collective molecular dressingparameter ˜ f and the cavity dressing parameter increaseas g √ N /ω increases, so it is expected that the exciton-cavity coupling can enhance the dressing of both the ex-citon and the cavity. As we will see in the next section,a better measure for the degree of vibrational dressing isthe mean vibration number on a specific exciton/cavitystate, which involves both the amplitude and the dressingparameter.In order to see the crossover with the cavity detuningin the weak cavity coupling region more clearly, we plot inthe left panel of Fig. 4 the ground-state energy as a func- ω c (e V) E L PP ( e V ) ω c (eV) | C || S | Nh ˜ f FIG. 4. Left panel: The ground-state energy E LPP as a func-tion of the cavity frequency ω c . Right panel: Absolute valuesof the amplitudes | C | and | S | , as well as the parameters ˜ f and Nh as functions of ω c . The exciton-cavity coupling is setto be in the weak-coupling regime as g √ N /ω = 0 .
2. Otherparameters: N = 100, ω = 0 . eV , ε = 2 eV , J/ω = − / λ = 1. tion of the cavity frequency ω c for fixed exciton-cavitycoupling strength g √ N/ω = 0 .
2. A crossover can beseen around the resonant frequency 1 . eV . The rightpanel of Fig. 4 shows the behavior of the amplitudes | C | and | S | , and the dressing parameters N h and ˜ f . For ared-detuned cavity with frequency ω c < . eV , theLPP state | ψ LPP i is dominated by the cavity componentwith a small vibrational dressing. For a blue-detunedcavity, the ground state behaves more like a polaron.However, the cavity becomes more dressed though itsamplitude | C | decreases with increasing ω c . Note thatthe profiles of N h and ˜ f approach the limiting formsgiven by Eqs. (33) and (34) if g √ N /ω is lowered downfurther.Intriguing enough, it can be seen from Fig. 3 that C and S tend to be consistent with C and S as g √ N /ω increases for all the three cases considered. This behaviorcan be better understood by introducing the followingquasiparticle weight for the LPP: Z LPP = |h φ LP | ψ LPP i| , (35)which measures how similar the LPP wavefunction | ψ LPP i is to the vibration-free LP wave function | φ LP i .It is easy to show that Z LPP = | S S Θ + C Ce − h N | , (36)where Θ = Θ( h = 0). As can be seen in Fig. 3, Z LPP ap-proaches nearly unity monotonically as the exciton-cavitycoupling increases, which means that the LPP behaveslike a vibration-free LP in the strong-coupling regime.Actually, by investigating Eqs. (25)–(27) in theultrastrong-coupling limit g √ N /ω → ∞ , we have C ≈− S ≈ / √ N h ≈ ˜ f ≈ λS ≈ λ/
2, and ˜ f q ≈ q = 0), which gives the asymptotic form of | ψ LPP i , | ψ LPP i ≈ √ | c i − | i ) e − λ √ N ( b † − b ) | vac v i≈ | φ LP i e − λ √ N ( b † − b ) | vac v i , (37)where | vac v i denotes the vibrational vacuum state.Equation (37) indicates that in the ultrastrong-coupling limit, the LPP state tends to be a separablestate , which is consistent with the bare LP state dressedby the zero-momentum vibrational mode. Furthermore,the vibrational dressing part becomes negligible for largeaggregates with N ≫
1. Correspondingly, the quasipar-ticle weight approaches Z LPP ≈ e − λ N , g √ N /ω → ∞ . (38) IV. THE PHOTOLUMINESCENCE LINESTRENGTH AND THE MEAN NUMBER OFVIBRATIONSA. Zero temperature
In order to better understand the influence of strongexciton-cavity coupling on the molecular system, it isinstructive to study the variation of several observableswith the cavity coupling strength. In this section, we cal-culate the photoluminescence line strength and the meannumber of vibrations using the results obtained in the lastsection.After photoexcitation, a J-aggregate loses its excess en-ergy and reaches the bottom of the exciton band quickly,so that the emission process originates mainly near theband bottom. When the cavity mode is present, the LPPstate | ψ LPP i takes the role of such a band bottom exciton.The 0 − ξ photoluminescence line strength I − ξ arisingfrom transitions between | ψ LPP i and the excitonic groundstate with ξ vibrations is defined as [29] I − ξ = 1 µ X P q n q = ξ |h{ n q }| ˆ µ | ψ LPP i| , (39)where the transition dipole moment operator is given byˆ µ = µ X n ( | n ih vac | + | vac ih n | ) . (40)By inserting Eqs. (29) and (40) into Eq. (39), we obtain(see Appendix C) I − ξ = ( S Θ ) ξ ! N ξ X n ( G n ) ξ , (41)where G n = G ∗ n = X q e iqn ˜ f q ˜ f − q . (42) g √ N/ω I − g √ N/ω I − g √ N/ω I − g √ N/ω N I − /I − ω c =2 eV ω c =1.7864 eV ω c =1 eV FIG. 5. The zero-temperature photoluminescence linestrength I − , I − , I − and the strength ratio NI − /I − calculated using Eq. (41). The results for three sets of cav-ity frequencies ω c = 2 eV (solid line), ω c = 1 . eV (dashedline), and ω c = 1 eV (dotted line) are shown. Other param-eters: N = 100, ω = 0 . eV , ε = 2 eV , J/ω = − /
4, and λ = 1. The first three cases for ξ = 0 ,
1, and 2 can be calculatedas I − = N ( S Θ ) , (43) I − = ( S Θ ˜ f ) , (44) I − = ( S Θ ) N X q ( ˜ f q ˜ f − q ) . (45)The zero-temperature 0 − ξ photoluminescence linestrength I − ξ for ξ = 0 ,
1, and 2 is shown in the first threepanels of Fig. 5, where results for both the resonant andnonresonant cases are presented. For the blue-detunedcavity frequency ω c = 2 eV , we observe a nonmonotonicbehavior in the 0 − I − with increasing g . For the resonant and red-detuned cavity, an amplifi-cation of I − is observed. The last panel of Fig. 5 showsthe line strength ratio, N I − I − = ˜ f , (46)which is proportional to an effective Huang-Rhys fac-tor [29] and decreases with increasing exciton-cavity cou-pling for resonant and blue-detuned cavity frequencies.However, the increase of this ratio for a red-detuned cav-ity shows that the cavity coupling can actually increasethe effective exciton-vibration coupling from weak to in-termediate cavity coupling regions.The above results can be further demonstrated bystudying the mean number of vibrations in the LPP state, (a) g √ N/ω N v −3 g √ N/ω ω c =1.7864 eV ω c =2 eV N (c)v ( ω c =1.7864 eV)N (site)v ( ω c =1.7864 eV)N (c)v ( ω c =2 eV)N (site)v ( ω c =2 eV) (b) g √ N/ω N v −3 g √ N/ω ω c =1 eV ω c =1.5 eV N (c)v ( ω c =1 eV)N (site)v ( ω c =1 eV)N (c)v ( ω c =1.5 eV)N (site)v ( ω c =1.5 eV) FIG. 6. The total mean number of vibrations (left panels)and mean vibration numbers on the cavity state and the localexciton state (right panels) for (a) ω c = 1 . eV , (b) ω c = 1 and 1 . eV . Other parameters: N = 100, ω = 0 . eV , ε = 2 eV , J/ω = − /
4, and λ = 1. which is an important measure of vibrational dressing ofspecific exciton-cavity states [66], N v = h ψ LPP | X q b † q b q | ψ LPP i . (47)Straightforward calculation gives (see Appendix C) N v = ( S Θ ) G N e G N + N C h = N N (site)v + N (c)v (48)where N (c)v = N C h , (49)is the mean number of vibrations projected onto the cav-ity state | c i , and N (site)v = ( S Θ ) G N e G N (50)is the mean number of vibrations in the cloud surround-ing a local exciton, which is identical for all sites due tothe circular symmetry of the molecule.The left panels of Fig. 6(a) and 6(b) show the totalmean number of vibrations N v in | ψ LPP i for ω c = 1 . eV , and ω c = 1 and 1 . eV , respectively. In con-trast to the monotonic decreasing of N v for the resonantand blue-detuned cases, a red-detuned cavity induces an0increase of N v from the weak- to intermediate- couplingregion. As the exciton-cavity coupling is increased fur-ther up to the strong coupling regime, N v will decreasefrom its maximal value at a crossover coupling strength.The decrease of N v with increasing g indicates that theoverall vibrational dressing in | ψ LPP i tends to fade awayin the strong-coupling regime.In the right panels of Fig. 6, we plot the correspond-ing results for the mean vibration numbers on the cavitystate and on the local exciton state. For all cases consid-ered, the mean number of vibrations N (c)v on the cavitystate increases monotonically with increasing g √ N /ω ,indicating an enhancement of vibrational dressing of thecavity mode. This in turn leads to a drop of N (site)v in thestrong-coupling regime, implying an exciton-cavity cou-pling induced vibrational decoupling of excitons. It isinteresting to note that the total mean vibration number N v shows a similar trend as N (site)v for both the blue-detuned and red-detuned cases. B. Thermal effects
As mentioned in Sec. I, the temperature-dependentvariational Merrifield transformation method also allowsus to study static properties of the system at thermalequilibrium. To calculate the thermal average of observ-able ˆ O at inverse temperature β = 1 /k B T , we turn tothe Merrifield frame where the zero-order density matrix˜ ρ ( β ) is both separable and diagonal,˜ ρ ( β ) ≈ ˜ ρ ( β ) = 1 Z S Z v ( X η e − βE η | η ih η | ) e − βH v , (51)where Z S and Z v = 1 / (1 − e − βω ) N are the partitionfunctions for the exciton-cavity system and the vibra-tional bath, respectively. The thermal average of ˆ O canthus be calculated by O ( β ) = Tr S Tr v [˜ ρ ( β )e S ˆOe −S ] . (52)After a straightforward calculation, we arrive at the fol-lowing expressions for the finite-temperature 0 − − I − = Θ Z S Z v X n X η e − βE η x η e ik ( η ) n e − GnN e − βω , (53)and I − = Θ Z S Z v X n X η e − βE η x η e ik ( η ) n e − GnN e − βω e − βω (cid:20) N + 2 N ( G n coth βω − G ) (cid:21) . (54)The finite-temperature extension of the mean number ofvibrations on the cavity state and on the local excitonstate has simple forms, N (c)v = Z y Z S N (¯ n + h ) , (55) and N (site)v = 1 N Z x Z S ( N ¯ n + G N ) , (56)where Z x and Z y are defined in Eq. (A10), and ¯ n =1 / ( e βω −
1) is the mean occupation number of the freevibrational bath. In the zero-temperature limit, we have Z x /Z S → S , Z y /Z S → C , and ¯ n →
0, and hencethe zero-temperature results given by Eqs. (48)–(50) arerecovered. Actually, the factor P η e − βE η x η e ik ( η ) n /Z S in I − and I − , and the ratios Z x / y /Z S in N (c)v and N (site)v , are determined by the energy gap between theLPP state and the lowest excitonic dark state with wavevector | k = 2 π/N i , namely, ∆ E = E π/N − E D , which ismuch higher than the thermal energy k B T for relativelystrong exciton-cavity coupling g √ N /ω ≥
1. We thusexpect that the variational parameters { ˜ f q } and h areclose to those in the zero-temperature limit and almosttemperature independent.Fig. 7(a) shows the line-strength ratio N I − /I − asa function of temperature T for fixed exciton-cavity cou-pling g √ N /ω = 4. At low temperatures, the ratio ap-proaches the zero-temperature result, ˜ f . As the temper-ature increases, this ratio increases due to the decreaseof I − and increase of I − . This temperature depen-dence of I − and I − originates from the reductionof the LPP population and thermal excitation of vibra-tions at high temperatures. Figure 7(b) shows the tem-perature dependence of the mean number of vibrations.We see that N (c)v , N (site)v , and N v all increase with in-creasing temperature. At low temperatures, the meanoccupation number ¯ n is much smaller than h , so that N (c)v ≈ N C h . As the temperature is increased acrossa turning point at which ¯ n is comparable with h [theinset of Fig. 7(b)], the thermal occupation of vibrationsdominates and N (c)v increases rapidly to N (c)v ≈ N C ¯ n in the high-temperature limit.Thus, due to the large energy gap formed in the strongexciton-cavity coupling region, the static properties ofthe system in thermal equilibrium are almost tempera-ture independent below a crossover temperature which isrelated to the degree of vibrational dressing of the cavityand excitons. However, thermal excitation of vibrationsdominates the behavior above the crossover temperature. V. CONCLUSIONS
In this work, we developed a microscopic theory fordescribing organic molecules coupled to a single cavitymode. The molecule is modeled by the Holstein Hamil-tonian that explicitly includes the intramolecular vibra-tions. By employing a temperature-dependent varia-tional approach combining a generalized Merrifield trans-formation with the Bogoliubov inequality, we could treatthe exciton-vibration coupling and exciton-cavity cou-pling on an equal footing. The generalized canonical1 (a)
50 100 150 200 250 30000.511.522.533.54
T (K) NI − /I − ˜ f (b)
50 100 150 200 250 30000.010.020.030.040.050.060.070.080.09
T (K) N (c)v N N (site)v N v
100 200 30000.511.5 x 10 −3 T (K) ¯ nh FIG. 7. (a) The line-strength ratio NI − /I − . (b) Thetotal mean number of vibrations N v , N (c)v , and NN (site)v asfunctions of temperature T for fixed exciton-cavity coupling g √ N /ω = 4. The dotted curve in (a) represents the valueof NI − /I − in the zero-temperature limit, ˜ f . The insetin (b) shows the temperature dependence of the vibrationoccupation number ¯ n and h . Other parameters: N = 50, ω = 0 . eV , ε = 2 eV , J/ω = − / λ = 1, and ω c =1 . eV . transformation we proposed takes the vibrational dress-ing of both the excitons and the cavity into account. Theground state of the system (within the single-excitationsubspace), which we refer to as a lower polaron polari-ton, is shown to be a hybrid state of excitonic, photonic,and vibrational degrees of freedom, and contains bothpolaronlike and polaritonlike structures.Using the above results, explicit expressions forthe quasiparticle weight, the photoluminescence linestrength, and the mean number of vibrations are ob-tained in terms of the optimal variational parameters.The dependence of these quantities upon the exciton-cavity coupling strength shows that the cavity state gainsa profound vibrational dressing in the strong cavity cou-pling regime, while the excitons tend to decouple fromthe vibrations. Finally, we study the temperature depen-dence of the photoluminescence line strength and meannumber of vibrations and show that these quantities arenot affected by the temperature at relatively low temper-atures, but mainly controlled by the thermal excitationof vibrations at high temperatures. Acknowledgements:
This work has been funded by theEuropean Research Council (ERC-2011-AdG, ProposalNo. 290981), by the European Union Seventh FrameworkProgramme under Grant Agreement No. FP7-PEOPLE-2013-CIG-618229, by the Spanish MINECO under Con-tract No. MAT2014-53432-C5-5-R, and by the Marade Maeztu programme for Units of Excellence in R&D(Grant No. MDM-2014-0377). [1] A. J. Heeger, Solid State Commun. , 673 (1998).[2] J. Shinar,
Organic Light Emitting Devices: A Survey (American Institute of Physics, New York, 2002).[3] F. Garnier, G. Horowitz, P. Valat, F. Kouki, and V.Wintgens, Appl. Phys. Lett. , 2087 (1998).[4] M. Granstr¨om, K. Petritsch, A. C. Arias, A. Lux, M. R.Andersson, and R. H. Friend, Nature (London) , 257(1998).[5] A. Dodabalapur, L. Torsi, and H. E. Katz, Science ,270 (1995).[6] H. van Amerongen, L. Valkunas, and R. van Gron-delle, Photosynthetic Excitons (World Scientific, Singa-pore, 2000).[7] G. D. Scholes, G. R. Fleming, A. Olaya-Castro, and R.van Grondelle, Nat. Chem. , 763 (2011).[8] V. M. Agranovich, H. Benisty, and C. Weisbuch, SolidState Commun. , 631 (1997).[9] D. G. Lidzey, D. D. C. Bradley, M. S. Skolnick, T. Virgili,S. Walker, and D. M. Whittaker, Nature (London) , 53 (1998).[10] D. G. Lidzey, D. D. C. Bradley, T. Virgili, A. Armitage,M. S. Skolnick, and S. Walker, Phys. Rev. Lett. , 3316(1999).[11] R. J. Holmes and S. R. Forrest, Phys. Rev. Lett. ,186404 (2004).[12] S. K´ena-Cohen, M. Davan¸co, and S. R. Forrest, Phys.Rev. Lett. , 116401 (2008).[13] P. Michetti, L. Mazza, and G. C. La Rocca, in OrganicNanophotonics (Springer, Berlin Heidelberg, 2015), Vol.39.[14] T. Schwartz, J. A. Hutchison, C. Genet, and T. W. Ebbe-sen, Phys. Rev. Lett. , 196405 (2011)[15] J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux, andT. W. Ebbesen, Angew. Chemie , 1624 (2012).[16] J. Bellessa, C. Bonnand, J. C. Plenet, and J. Mugnier,Phys. Rev. Lett. , 036404 (2004).[17] T. K. Hakala, J. J. Toppari, A. Kuzyk, M. Pettersson,H. Tikkanen, H. Kunttu, and P. T¨orm¨a, Phys. Rev. Lett. , 053602 (2009).[18] A. Berrier, R. Cools, C. Arnold, P. Offermans, M. Crego-Calama, S. H. Brongersma, and J. G´omez-Rivas, ACSNano , 6226 (2011).[19] J. Dintinger, S. Klein, F. Bustos, W. L. Barnes, and T.W. Ebbesen, Phys. Rev. B , 035424 (2005).[20] P. Vasa, R. Pomraenke, S. Schwieger, Y. I. Mazur, V.Kunets, P. Srinivasan, E. Johnson, J. E. Kihm, D. S.Kim, E. Runge, G. Salamo, and C. Lienau, Phys. Rev.Lett. , 116801 (2008).[21] S. R. K. Rodriguez, J. Feist, M. A. Verschuuren, F. J.Garc´ı, and J. G´omez Rivas, Phys. Rev. Lett. , 166802(2013).[22] A. I. V¨akev¨ainen, R. J. Moerland, H. T. Rekola, A.-P.Eskelinen, J.-P. Martikainen, D.-H. Kim, and P. T¨orm¨a,Nano Lett. , 1721 (2014).[23] G. A. Wurtz, P. R. Evans, W. Hendren, R. Atkinson, W.Dickson, R. J. Pollard, A. V. Zayats, W. Harrison, andC. Bower, Nano Lett. , 1297 (2007).[24] G. Zengin, M. Wers¨all, S. Nilsson, T. J. Antosiewicz, M.K¨all, and T. Shegai, Phys. Rev. Lett. , 157401 (2015).[25] E. Eizner, O. Avayu, R. Ditcovski, and T. Ellenbogen,Nano Lett. , 6215 (2015).[26] R. Chikkaraddy, B. de Nijs, F. Benz, S. J. Barrow, O. A.Scherman, E. Rosta, A. Demetriadou, P. Fox, O. Hess,and J. J. Baumberg, Nature (London) , 127 (2016).[27] S. K´ena-Cohen, S. A. Maier, and D. D. C. Bradley, Adv.Opt. Mater. , 827 (2013).[28] J. A. ´Cwik, S. Reja, P. B. Littlewood, and J. Keeling,Europhys. Lett. , 47009 (2014).[29] F. C. Spano, J. Chem. Phys. , 184707 (2015).[30] F. Herrera and F. C. Spano, Phys. Rev. Lett. , 238301(2016).[31] J. Galego, F. J. Garcia-Vidal, and J. Feist, Phys. Rev. X , 041022 (2015).[32] J. A. ´Cwik, P. Kirton, S. De Liberato, and J. Keeling,Phys. Rev. A , 033840 (2016).[33] M. A. Zeb, P. G. Kirton, and J. Keeling,arXiv:1608.08929.[34] T. Holstein, Ann. Phys. (NY) , 325 (1959); , 343(1959).[35] T. Kobayashi, J-Aggregates (World Scientific, Singapore,1996).[36] K. Mukai, S. Abe, and H. Sumi, J. Phys. Chem. B ,6096 (1999).[37] A. Damjanovi´c, I. Kosztin, U. Kleinekath¨ofer, and K.Schulten, Phys. Rev. E , 031919 (2002).[38] R. E. Merrifield, J. Chem. Phys. , 445 (1964).[39] D. R. Yarkony and R. J. Silbey, J. Chem. Phys. , 5818(1977).[40] Y. C. Cheng and R. J. Silbey, J. Chem. Phys. ,114713 (2008).[41] D. P. S. McCutcheon and A. Nazir, New J. Phys. ,113042 (2010).[42] D. P. S. McCutcheon, N. S. Dattani, E. M. Gauger, B.W. Lovett, and A. Nazir, Phys. Rev. B , 081305(R)(2011).[43] C. Roy and S. Hughes, Phys. Rev. X , 021009 (2011).[44] P. Kaer, T. R. Nielsen, P. Lodahl, A. -P. Jauho, and J.Mork, Phys. Rev. B , 085302 (2012).[45] C. Weisbuch, M. Nishioka, A. Ishikawa, and Y. Arakawa,Phys. Rev. Lett. , 3314 (1992).[46] D. M. Coles, P. Michetti, C. Clark, W. C. Tsoi, A. M. Adawi, J-S. Kim, and D. G. Lidzey, Adv. Funct. Mater., , 3691 (2011).[47] N. Wu, A. Nanduri, H. Rabitz, Phys. Rev. A , 062105(2014).[48] E. A. Bloemsma, M. H. Silvis, A. Stradomska,and J. Knoester, Chem. Phys. in press (2016),doi:10.1016/j.chemphys.2016.06.018.[49] A. Shalabney, J. George, J. Hutchison, G. Pupillo, C.Genet, and T. W. Ebbesen, Nat. Commun. , 5981(2015).[50] J. del Pino, J. Feist, and F. J. Garcia-Vidal, New J. Phys. , 081001 (2015).[51] J. del Pino, J. Feist, and F. J. Garcia-Vidal, J. Phys.Chem. C , 29132 (2015).[52] A. Strashko and J. Keeling, Phys. Rev. A , 023843(2016).[53] R. P. Feynman, Statistical Mechanics: A Set Of Lectures ,Advanced Book Classics (Westview, Boulder, CO, 1998).[54] L. -C Ku, S. A. Trugman, and J. Bonˇca, Phys. Rev. B , 174306 (2002).[55] M. R. Philpott, J. Chem. Phys. , 2039 (1971).[56] M. Hoffmann and Z. G. Soos, Phys. Rev. B , 024305(2002).[57] F. C. Spano, J. Chem. Phys. , 5877 (2002).[58] H. Yamagata and F. C. Spano, J. Phys. Chem. Lett. ,622 (2014).[59] H. De Raedt and A. Lagendijk, Phys. Rev. B , 1671(1984).[60] E. Jeckelmann and S. R. White, Phys. Rev. B , 6376(1998).[61] A. S. Alexandrov, V. V. Kabanov, and D. K. Ray, Phys.Rev. B , 9915 (1994).[62] G. Wellein and H. Fehske, Phys. Rev. B , 4513 (1997).[63] D. W. Brown, K. Lindenberg, and Y. Zhao, J. Chem.Phys. , 3179 (1997).[64] Y. Toyozawa, Prog. Theor. Phys. , 29 (1961).[65] V. M. Stojanovi´c, T. Shi, C. Bruder, J. I. Cirac, Phys.Rev. Lett. , 250501 (2012).[66] D. Yarkony and R. Silbey, J. Chem. Phys. , 1042(1976). Appendix A: The saddle point equations
We assume N is even; the analysis for odd N is sim-ilar. The minimized Bogoliubov free energy given byEq. (19) can be reached at the saddle point that is de-termined by the stationary conditions ∂F B /∂f α = 0, for α = 0 , , ..., N , which result in A α Z S = 0 , (A1)3where A = 2 X k =0 e − βE k (cid:20) ( f − λ ) ω − f − f ) ˜ J cos k coth βω (cid:21) + e − βE U { (1 + cos θ ) (cid:20) ( f − λ ) ω − f − f ) ˜ J coth βω (cid:21) + sin θ ( f − h )˜ g √ N coth βω } + e − βE D { (1 − cos θ ) (cid:20) ( f − λ ) ω − f − f ) ˜ J coth βω (cid:21) − sin θ ( f − h )˜ g √ N coth βω } , (A2) A N = 2 X k =0 e − βE k (cid:20) f N ω − f N − f N − ) ˜ J cos k coth βω (cid:21) + e − βE U { (1 + cos θ ) (cid:20) f N ω − f N − f N − ) ˜ J coth βω (cid:21) + sin θ ( f N − h )˜ g √ N coth βω } + e − βE D { (1 − cos θ ) (cid:20) f N ω − f N − f N − ) ˜ J coth βω (cid:21) − sin θ ( f N − h )˜ g √ N coth βω } , (A3)and A n = 2 X k =0 e − βE k (cid:20) f n ω − (2 f n − f n − − f n +1 ) ˜ J cos k coth βω (cid:21) + e − βE U { sin θ ( f n − h )˜ g √ N coth βω θ ) (cid:20) f n ω − (2 f n − f n − − f n +1 ) ˜ J coth βω (cid:21) } + e − βE D {− sin θ ( f n − h )˜ g √ N coth βω − cos θ ) (cid:20) f n ω − (2 f n − f n − − f n +1 ) ˜ J coth βω (cid:21) } , (A4)for n = 1 , , ..., N − ∂F B /∂h gives A h Z S = 0 , (A5)where A h = e − βE U [(1 − cos θ ) N hω − sin θ ( ˜ f − N h )˜ g √ N coth βω e − βE D [(1 + cos θ ) N hω + sin θ ( ˜ f − N h )˜ g √ N coth βω . (A6)An equivalent alternative form of the saddle-pointequations can be obtained through linearly combining these equations, yielding λN h = Z S Z x − ω ˜ g √ N coth βω Z y Z xy , (A7)˜ f N h = 1 − ω ˜ g √ N coth βω Z y Z xy , (A8)and λ ˜ f q = 1 − ˜ g √ N coth βω ω Z xy Z x − − cos q ) ˜ J coth βω ω − Z x X k =0 e − βE k (1 − cos k ) , (A9)for q = 0. Here, Z x = X η x η e − βE η ,Z y = X η y η e − βE η ,Z xy = X η x η y η e − βE η . (A10)Note that ˜ f q = ˜ f − q is real, and Eqs. (A7)–(A9) shouldalso be solved self-consistently. In the zero-temperaturelimit, only terms related to E D in Eqs. (A7)–(A9) survive;we then obtain Eqs. (25)–(27) in the main text. Appendix B: Second-order energy correction to E D From second-order time-independent perturbation the-ory, the second-order corrected energy E corr of the LPPstate is E corr = E D + ∞ X ν =1 ∆ ν , (B1)with the ν -vibration contribution∆ ν = X η X P q n q = ν |h η ; { n q }| ˜ V | D i| E D − ( E η + νω ) , (B2)where the dressed state | η ; { n q }i = Q q ( b † q ) nq √ n q ! | η i has n q vibrations in mode q . After a tedious but straightforwardcalculation, we obtain∆ = − C ω N h (cid:18) − S ω / ( E U − E D ) (cid:19) , (B3)for ν = 1, and∆ ν = 1 ν ! N ν X η E D − ( E η + νω ) { Sx η ˜ J ) [1 + ( − ν cos( k ( η ))] F ν,k ( η ) +2 Sx η ˜ J ˜ g √ N [ Cx η + ( − ν Sy η ] G ′ ν,k ( η ) +(˜ g √ N ) [ Cx η + ( − ν Sy η ] K ′′ ν,k ( η ) } , (B4)4for ν ≥
2, where we have defined F ν,p = 1 N X n e ipn (2 G n − G n − − G n +1 ) ν , G ′ ν,p = 1 N X n e ipn [( G ′ n − G ′ n − ) ν + ( G ′ n − G ′ n +1 ) ν ] , K ′′ ν,p = 1 N X n e ipn G ′′ νn , (B5)Here, G n is given by Eq. (42) and G ′ n = X q e iqn ˜ f q ˜ f ′− q , G ′′ n = X q e iqn ˜ f ′ q ˜ f ′− q , (B6)with ˜ f ′ q = ˜ f q − δ q hN. (B7) Appendix C: Derivation of the finite-temperaturephotoluminescence line strength I − ξ and the meannumber of vibrations At finite temperatures, we will adopt the zero-orderthermal equilibrium density matrix for the calculation ofthermal averages of observable ˆ O , which is assumed tobe in a separable form,ˆ O = ˆ O S ˆ O v , (C1)where ˆ O S and ˆ O v are operators of exciton/cavity and vi-brational degrees of freedom, respectively. Important ex-amples of observables in the above form include the 0 − ξ photoluminescence line strength I − ξ (with ˆ O S = | ih | and ˆ O v = P P q n q = ξ |{ n q }ih{ n q }| ) and the mean vibra-tion number projected onto state | η i (with ˆ O S = | η ih η | and ˆ O v = P q b † q b q ). It turns out to be convenient towork in the Merrifield frame where the zero-order den-sity matrix is diagonal and separable,˜ ρ ( β ) ≈ ˜ ρ ( β ) = ˜ ρ S ˜ ρ v , ˜ ρ S = 1 Z S X η e − βE η | η ih η | , ˜ ρ v = 1 Z v e − βH v . (C2)where β = 1 /k B T is the inverse temperature, and Z S and Z v = 1 / (1 − e − βω ) N are the partition functionsfor the exciton-cavity system and the vibrational bath,respectively. The representation of ˆ O in the Merrifieldframe is˜ˆ O = | c ih c | ˆO Scc e − B c ˆO v e B c + X mn | m ih n | ˆO Smn e − B m ˆO v e B n + X n ( | c ih n | ˆ O Scn e − B c ˆ O v e B n + H . c . ) (C3)where ˆ O Sxy = h x | ˆO S | y i for x, y = c or { n } , and H.c.stands for the Hermitian conjugate. The thermal average O ( β ) at inverse temperature β then can be calculated inthe Merrifield frame as O ( β ) = Tr S Tr v [˜ ρ ( β ) ˜ˆO]= ˆ O Scc ˜ ρ S , cc h e − B c ˆ O v e B c i v + X mn ˆ O Smn ˜ ρ S , nm h e − B m ˆ O v e B n i v +2 ℜ X n ˆ O Scn ˜ ρ S , nc h e − B c ˆ O v e B n i v . (C4)Let us first calculate the finite temperature 0 − ξ pho-toluminescence line strength I − ξ = 1 Z S X P q n q = ξ X mn X η e − βE η x η N e ik ( η )( n − m ) h{ n q }| e B n ˜ ρ v e − B m |{ n q }i . (C5)To calculate the matrix element in the second line of theabove equation, we invoke the following two identities e αb † − α ∗ b e − σb † b e − ( γb † − γ ∗ b ) = e − ( | α | + | γ | ) e − α ∗ γe − σ e ( α − γe − σ ) b † e − βb † b e ( γ ∗ − α ∗ e − σ ) b (C6)with b a bosonic annihilation operator, and h{ n q }| e P q α q b † q e − β P q b † q b q e P q γ q b q |{ n q }i = ξ X χ =0 e − β ( ξ − χ ) X P m q = χ Y q ( α q γ q ) m q m q ! n q ! m q !( n q − m q )! . (C7)After some algebra, we arrive at I − ξ = 1 Z S Z v X n X η e − βE η x η e ik ( η ) n e − P q | fq | N (1+ e iqn − βω ) ξ X χ =0 e − βω ( ξ − χ ) X P q n q = ξ X P m q = χ Y q | ˜ f q | N ! m q ( − e − βω + e iqn e − βω + e − iqn ) m q m q ! n q ! m q !( n q − m q )! . (C8)At zero temperature T = 0, only the term with χ = ξ survives in the summation over χ , and we hence recoverEq. (41) in the main text. For large ξ , it is difficult toobtain I − ξ in closed form. However, I − and I − canbe easily calculated and we thus obtain Eqs. (53) and(54) in the main text.We next turn to the calculation of the mean numberof vibrations N ( a )v in an arbitrary exciton/cavity state | a i . By inserting ˆ O S = | a ih a | and ˆ O v = P q b † q b q intoEq. (C4), and using the following identity:1 Z b Tr b [e − σ b † b e α b † − α ∗ b b † be − ( γ b † − γ ∗ b) ]= e − ( | α | + | γ | )(1+2 n b ) e α ∗ γ (1+ n b )+ αγ ∗ n b n b [1 + α ∗ γ (1 + n b ) e γ + αγ ∗ n b − (1 + n b )( | α | + | γ | )] , (C9)5with n b = 1 / ( e σ − N ( a )v = N |h c | a i| Z y Z S (¯ n + h )+ 1 N Z S X mn h m | a ih a | n i X η e − βE η x η e ik ( η )( n − m ) e − P q | ˜ fq | N [1 − cos( m − n ) q ](1+2¯ n ) X q e i | fq | N sin( m − n ) q ( ¯ n + | ˜ f q | N h e iq ( m − n ) (1 + ¯ n ) + e − iq ( m − n ) ¯ n − n (1 + ¯ n ) i) + 1 √ N Z S ℜ X n h c | a ih a | n i X η e − βE η x η y η e ik ( η ) n e − P q | ˜ fq | N (1+2¯ n ) e − h ( Nh − f )(1+2¯ n ) { ¯ n " N − (1 + ¯ n ) X q | f q | N + h [ f (1 + ¯ n ) + f ¯ n − ¯ n (1 + ¯ n ) N h ] } . (C10) where ¯ n = 1 / ( e βω −
1) is the mean occupation numberof the free vibrational bath.By choosing a = c and a = nn