Polariton condensation with saturable molecules dressed by vibrational modes

Polaritons, mixed light-matter quasiparticles, undergo a transition to a condensed, macroscopically coherent state at low temperatures or high densities. Recent experiments show that coupling light to organic molecules inside a microcavity allows condensation at room temperature. The molecules act as saturable absorbers with transitions dressed by molecular vibrational modes. Motivated by this we calculate the phase diagram and spectrum of a modified Tavis-Cummings model, describing vibrationally dressed two-level systems, coupled to a cavity mode. Coupling to vibrational modes can induce re-entrance, i.e. a normal-condensed-normal sequence with decreasing temperature and can drive the transition first order.

Microcavity polaritons (superpositions of photons and excitons) are light mass bosonic quasiparticles that have been observed to form a Bose-Einstein condensate, and the non-thermalized cousin, a polariton laser [1,2]. The onset of coherence ranges from around 20K in CdTe [3] to room temperature for materials such as GaN and ZnO with large exciton dipole moments [4][5][6], as well as organics such as anthracene [7,8] or J-aggregates [9]. Organic molecules have localized (Frenkel) excitons that naturally have strong optical coupling, but additionally have strong coupling to localized molecular vibrational modes, which introduces a new set of behaviors and possibilities not normally present in inorganic systems. These vibronic couplings are mostly studied only at the single exciton level, but for example have recently been shown to play a complex and crucial role in energy transfer in light harvesting complexes [10][11][12]. Organic polariton condensation provides a natural venue to explore many-particle effects arising from simultaneous strong exciton-phonon and exciton-photon coupling. In addition, since polariton condensates and polariton lasers can be seen as an exotic variety of ultra-low threshold laser, understanding polariton lasing with organic molecules may provide a pathway to improve the performance of solid state organic lasers [13].
The theory of excitons in molecular crystals has been extensively developed by Davydov [14] and more recently Agranovich [15]. However, modelling organic molecules with strong coupling to light is less developed. Work to date has focused on using rate equations to calculate the luminescence spectrum [16] and relaxation processes [17], as well as the effects of disorder on the spectrum [18,19]. Focusing on low densities, these works have generally used a theory of a weakly interacting gas of polaritons, in some cases explicitly derived [20] from models of molecules as saturable absorbers. Such theories well describe the very low density regime, but are less accurate at higher densities. Instead, one may deal with a model of molecules as saturable absorbers (two-level systems) coupled to radiation in an optical cavity, by extending the Dicke [21] or Tavis-Cummings [22] model [23,24], which allows one to work at higher densities.
In this letter, we consider the effects of strong molecular exciton-phonon coupling on collective behavior. We investigate how the phase transition to a superradiant state [25][26][27] at finite density [23,24] is modified by the coupling to phonons. We find that the phase transition can become first order for sufficiently strong phonon coupling. We also investigate an apparent paradox raised by the existence of phonon replicas in the luminescence spectrum: Bose-Einstein condensation occurs when the chemical potential reaches the bottom of the excitation spectrum, at which point a mode becomes macroscopically occupied. However, at finite temperature phonon replicas exist both above and below the bare polariton frequency, and as these are bosonic modes, it appears that condensation of these replicas should occur. We find that instead, condensation is avoided as the photon spectral weight vansihes when the chemical potential crosses such modes. Nonetheless, at large exciton-phonon coupling, the polariton mode which does condense contains a strong admixture of different phonon replicas.
The model we study generalizes the Dicke model [21], which describes two-level systems (molecules) coupled to a single photon mode in the microcavity. To this we add a coupling between two-level systems and vibrational modes of the molecules. We thus have: Hereψ † is the creation operator for cavity photon, where ω c is the cavity frequency. Twolevel systems are described by the Pauli matrices σ i n , ǫ = ǫ − µ where ǫ is the optical transition frequency of the molecules, the molecule-photon coupling strength is denoted g. The chemical potential µ controls the total excitation density (excited molecules and photons), which is conserved [24]. Localized vibrational modes of the molecules are created by the operatorsâ † n , these have a frequency Ω and their coupling to the optical transition is given by the Huang-Rhys parameter S. We assume an homogeneous system of N two-level systems, and note that the characteristic polariton splitting is given by the combination g √ N which is intensive since g scales as the inverse root of the cavity mode volume. This model neglects finite momentum photon modes, as well as Forster transfer of excitations between molecules. The latter is unlikely to affect the thermodynamics, as the bandwidth due to Forster transfer is small compared to the scale of the photon dispersion. Ignoring finite momentum photon modes is valid for the high density regime, where mean-field approaches (which neglect finite momentum fluctuations) become valid [28].
In the following we consider the case where the cavity frequency is detuned above the molecular transition, by a value of the order of g We focus on this regime for a number of reasons. Firstly, the clear observation of equilibrated Bose-Einstein condensation of polaritons in inorganic materials [3] required such detuning, so that the lower polariton has a strong excitonic fraction, increasing the polaritonpolariton scattering leading to thermalization [29,30]. Secondly, the Dicke model can show multiple normalsuperradiant phase transitions [24] when ∆ > 0; this is linked to the existence of "Mott-lobes" in the Jaynes-Cummings-Hubbard model [31]. In the Dicke model only two such Mott lobes exist. Observation of such behavior is not attainable for Wannier excitons, as one cannot easily saturate the excitonic occupancy [32], but may well be observable in molecular systems; hence we aim to see whether such effects survive coupling to phonons.
Regarding the vibrational mode frequency we consider two cases, Ω ∼ g √ N /2 and Ω ∼ g √ N /20. The first value is relevant to anthracene [7] for which g √ N ∼ 300meV and the Ω ∼ 180meV. This is a relatively large vibrational frequency, due to the stiffness of the fully saturated anthracene molecule, and so we also present illustrations of the behavior at lower phonon frequency. The value of the Huang-Rhys parameter varies significantly between different materials; typical values are around S = 1. In order to show more clearly the behaviour arising from strong coupling, we show below also results that would arise with even larger S.
To calculate the phase diagram of Eq. (1) we use a semiclassical treatment of the photon field, known to be exact [24] in the thermodynamic limit N → ∞, g √ N → const., i.e. we replace ψ → λ √ N , subject to the self-consistency conditionω c λ = −g √ N σ − . Here σ − is found by diagonalizing the on-site problem h = 1 2 ǫ + Ω √ S(â +â † ) σ z + g √ N (λσ + + H.c.) + Ωâ †â numerically, truncating the maximum number of vibrational excitations at n max ≫ S, and thermally populat-ing the resulting eigenstates. Anticipating possible first order transitions, one must also compare the free energies of the normal (λ = 0) and condensed (λ = 0) solutions to determine the global minimum free energy.  Figure 1 shows phase diagrams, i.e. critical temperature or coupling strength as a function of chemical potential. As for the regular Dicke model, the behavior at zero temperature is reminiscent of normal state Mott lobes, separated by a superradiant (condensed) state [24] near to µ = ǫ. At µ = ǫ, T = 0, the critical coupling strength vanishes. This is analogous to the existence of superfluidity at the boundary between two Mott lobes [33]. Coupling to phonons adds several additional features not present in the regular Dicke model. Visible in Fig. 1(b) is re-entrant behaviour as a function of temperature, discussed below. At larger values of S (see Fig. 2 and later discussion) there can be first order transitions into the superradiant state. For yet larger S (not shown) a first order jump can also occur within the condensed phase. (First order transitions have also recently been noted within other variants of the Dicke model [34]).
We next discuss these various features in more detail. The existence of the re-entrant behavior in Fig. 1(b) can be understood from the effect of phonon sidebands of the optical transition meaning that a lower chemical potential is required to reach the condensed state, i.e. the optical transition frequency is reduced due to coupling to phonons. For this to happen we require that the transition starts from a multiphonon ground state and ends in an excited state with fewer vibrations. Hence, such a feature is only visible for k B T > Ω. Since the superradiant bubble vanishes for k B T ≫ g √ N , this feature is only visible in cases where Ω ≪ g √ N , Fig. 1(b). The existence of a first order normal to superradiant transition can be understood by considering the extent of polaron formation at large S, i.e. entanglement of the phonon and electronic states of the molecule. Such entanglement reduces the overlap between ground and excited states, and so reduces the coupling to light, but also lowers the phonon energy. Such a state is clearly favored in the normal state. In the condensed state, it is instead favorable to increase the optical polarization by suppressing the entanglement, and having similar phonon configurations for both electronic states. First order transition arises due to a sudden switch between these states. At small T the strength of the first order jump is largest at points near to, but not exactly at, µ = ǫ.
In order to gain a clearer understanding of the first order transition, we introduce a variational ansatz which captures the behavior at large S. This ansatz can be framed as a mean-field treatment of the phonons within a coherent state path integral representation. To allow for possibile entanglement between phonons and the electronic state, we make a variational polaron transform [35,36] . This transformation is exact, but if followed by a mean-field approximationâ → α,then η, α, and λ become variational parameters. This phonon mean-field ansatz is valid if the vibrational states are approximately coherent states, i.e. when the typical phonon occupations, which are ∼ S, are ≫ 1.
The free energy of this phonon mean-field theory can be written as: where ζ = δ 2 + (gλ) 2 is written in terms of the effective molecular transition frequency δ and phonon-suppressed optical couplingg, given by: Minimizing this free energy with respect to the variational parameters λ, α and η we obtain the gap equationω c λ =g 2 λtanh(βζ)/2ζ. Defining κ = g 2 λ 2 tanh(βζ)/[ζ 2 − (δ tanh(βζ)) 2 ] we may write the equations for α, η as: The equation for η, describing the extent of polaron state formation is instructive. At small λ or large Ω, η → 1, and one has fully developed polarons (vibrational state fully entangled with the electronic state, as in the absence of any applied field). If, on the other hand, λ is large then η → 0 and the polaron formation is suppressed; i.e. strong driving suppresses polaron formation [35,36]. However, an additional level of self-consistency appears in the current problem, not present for normal variational polaron approaches. The photon field λ depends on the polarization of the molecules, and the effective coupling strengthg. Therefore, when the bare coupling g is small, the photon field is small, and polarons are well developed, Comparison between critical g √ N from exact on-site diagonalization (left) and phonon mean-field ansatz (right) for S = 6, Ω = 1, ∆ = 4, temperatures as indicated. The color scale shows λ = ψ / √ N at the phase boundary. Light (yellow) colors imply second order transition whereas dark (blue) first order transition. NB: around resonance the transition is always first order, though sometimes weakly further suppressing the effective couplingg. At larger couplings, polaron formation is suppressed, producing a stronger coupling. At zero temperature, this leads to a jump within the condensed phase, between a weakly and strongly polarized phase. Within the variational approach, such a jump occurs near ǫ = µ if S > 27/8. At non-zero temperature the same effect leads to the first order normal to condensed transitions, as a competition between states with different values of η.
In the limit S = 0, the phonons are irrelevant, and the both phonon mean-field theory and the numerical diagonalization reduce to the Dicke model. Figure 2 shows how the two approaches match in the large S limit. In this limit, one expects the phonon state can be approximately described as a coherent state. However, for smaller S, the phonon mean-field theory predicts more strongly first order transitions than the numerical approach. Numerics shows that the first order jump reduces and the transition becomes second order as S → 0.
When the phase transition is second order, the mechanism of the phase transition is as for Bose-Einstein condensation: macroscopic occupation of a bosonic mode occurs when the chemical potential crosses this mode. However, this presents an apparent paradox when for a system with phonon replicas. The luminescence spectrum shows replicas both above the bare polariton mode, associated with creating additional phonons, and at T > 0, also below the bare polariton, associated with destroying existing thermally populated phonons. Thus at any non-zero temperature, an infinite number of (increasingly weak) bosonic modes exist below the bare polariton, and these appear susceptible to become macroscopically occupied at low chemical potential. This clearly does not occur in the phase diagrams presented so far.
The resolution of this apparent paradox lies in the changing nature of the luminescence spectrum as one varies the chemical potential. The luminescence spec- where ρ is the total excitation density, (b) phonon spectral weight vs phonon change for the critical mode for various T trum arises from the hybridization of the photon with the various phonon replicas of the molecular transition, and can be found from analytic continuation of the photon thermal Green's function [37], D(ω n ), given by where Σ +− is the photon self-energy, Z is the partition function, and α pq = | p|σ − |q | 2 is the overlap between the polaron states corresponding to p excitations in the electronic ground state and q excitations in the electronic excited state. The poles of the Green's function correspond to the excitation spectrum visible in luminescence, thermally occupied; i.e. the luminescence spectrum P (ω) = −n B (ω)ℑ [πD(iω n = ω − µ + i0)].
The luminescence spectrum calculated from these expressions is shown in Fig. 3(a), along with the chemical potential. It is clear there are many places where the chemical potential appears to cross the phonon replicas, but crucially, all these crossings are avoided because the photon spectral weight of the replica vanishes at these points. Even in the regular Dicke model such behavior occurs: one can reach an inverted state of the two-level systems while avoiding condensation if the chemical potential crosses the exciton energy at a point where the normal modes are pure exciton and pure photon [24]. The same behavior occurs repeatedly in the case with phonons. Mathematically this behavior arises from collision between the zeros and poles of the Green's functions. The zeros ω * of the Green's function are given by ω * −µ =ǫ+(q−p)Ω. Ifǫ+(q−p)Ω = 0 then one may note that the numerator in the summand in Eq. (6) disappears -at such points poles and zeros cancel, and ω * = µ, thus whenever a pole and zero coincide, the chemical potential crosses a pole. Note however the reverse is not guaranteed; the chemical potential may reach a pole without an associated zero, at which point condensation occursthis happens at the right edge of Fig. 3(a).
In a similar manner, one may also find phonon composition of the various modes, and in particular, the mode which eventually attains a macroscopic occupation. To do this, we may rewrite the above Green's function as the inverse of a matrix in the space of photon and exciton-phonon states, so that one may determine the residue (spectral weight) for a given state, and thus find the composition of a given normal mode. Introducing A p,q = g 2 N e βǫ/2−βpΩ − e −βǫ/2−βqΩ α pq /Z this matrix has the form: The previous photon Green's function corresponds to D 00 , while the various elements D (pq),(pq) correspond to the weights of transitions involving particular initial (p) and final (q) phonon number states. Figure 3(b) shows the composition as a function of q − p, i.e. the number of absorbed phonons of the mode which acquires a macroscopic occupation, exactly at the critical point, for various temperatures. At low temperature, as expected no phonon emission is possible (no phonons are thermally populated in the ground state), and the distribution depends purely on the overlaps α 0,p , controlled by the Huang-Rhys factor. For the parameters used in Fig. 3(b), if g √ N > k B T > Ω, then re-entrant behavior as shown in Fig. 1(b) occurs. This corresponds to the transition at k B T > Ω being associated with macroscopic occupation of a mode involving phonon emission, as clearly illustrated in the figure.
In conclusion, we have shown that strong excitonphonon coupling modifies the phase diagram and the spectrum of a polariton condensate, in a parameter regime that should be accessible for experiments in organic microcavities.
JAC acknowledges support from EPSRC. SR acknowledges support from the Cambridge Commonwealth trust. SR, PBL and JK acknowledge discussions with S. Yarlagadda. JK acknowledges discussions with B. Lovett and R. Gomez-Bombarelli and financial support from EP-SRC program "TOPNES" (EP/I031014/1) and EPSRC (EP/G004714/2). Argonne National Laboratory's work supported by the U.S. Department of Energy, Office of Basic Energy Sciences under contract no. DE-AC02-06CH11357. JAC performed the numerical calculations, and SR performed the variational calculations. All authors contributed to the conception of the project, and the preparation of the manuscript.