Quantum theory of collective strong coupling of molecular vibrations with a microcavity mode

We develop a quantum mechanical formalism to treat the strong coupling between an electromagnetic mode and a vibrational excitation of an ensemble of organic molecules. By employing a Bloch-Redfield-Wangsness approach, we show that the influence of dephasing-type interactions, i.e., elastic collisions with a background bath of phonons, critically depends on the nature of the bath modes. In particular, for long-range phonons corresponding to a common bath, the dynamics of the"bright state"(the collective superposition of molecular vibrations coupling to the cavity mode) is effectively decoupled from other system eigenstates. For the case of independent baths (or short-range phonons), incoherent energy transfer occurs between the bright state and the uncoupled dark states. However, these processes are suppressed when the Rabi splitting is larger than the frequency range of the bath modes, as achieved in a recent experiment [Shalabney et al., Nat. Commun. 6, 5981 (2015)]. In both cases, the dynamics can thus be described through a single collective oscillator coupled to a photonic mode, making this system an ideal candidate to explore cavity optomechanics at room temperature.


Introduction
The field of cavity optomechanics explores the interaction between electromagnetic (EM) radiation and the quantized mechanical motion of nano-or micro-oscillators [1][2][3]. Recent developments promise a rich variety of applications such as precision mechanical measurements and coherent control of quantum states (see [4,5] for recent reviews). When the interaction becomes sufficiently large, coherent energy exchange between the optical and mechanical degrees of freedom becomes possible and the system could reach the strong coupling regime. Along this direction, a novel approach consists in coupling a cavity mode to a molecular bond vibration [6,7], which is in the ground state already at room temperature, without requiring any cooling. In order to achieve strong coupling, this has to be done by using an ensemble of molecules. This phenomenon has common ingredients with the collective strong coupling observed when electronic excitations of organic molecules interact with cavity [8][9][10] or plasmonic modes [11], which has been studied extensively during the last years (see [12,13] for recent reviews).
In the strong coupling with an ensemble of vibrational modes, the cavity resonance couples to a collective superposition of the molecular vibrations, the so-called bright state, forming hybrid states called polaritons. In principle, all other superpositions of vibrational excitations (the so-called dark states) remain uncoupled. However, within this setting, a thermal bath of low-frequency rovibrational modes, which normally only introduces dephasing [14,15], interacts with the vibrational excitations and may influence the system dynamics.
It is an open question how the dephasing mechanisms affect the strongly coupled dynamics and what role the dark states play. In particular, it is unknown whether the bright mode in such a system effectively behaves like a single isolated oscillator, which would enable the direct application of protocols developed in cavity optomechanics.
In this paper, we present a fully quantum theory of the phenomenon of collective strong coupling of molecular vibrations with a cavity mode. In contrast to classical transfer matrix calculations, which can be used to fit the experimental spectra [6,7] but only provide very limited information about the incoherent dynamics in the system, our approach naturally incorporates all incoherent processes, in particular those induced by dephasingtype interactions. We describe the system using a quantum-mechanical model of a single photon mode coupled to an ensemble of harmonic oscillators representing molecular vibrations (the C=O bond stretching mode of polyvinyl acetate at an energy of 215 meV in the experiment [6]). In order to incorporate the coupling of the molecular vibrations to low-frequency rovibrational modes, we assume that the oscillators are connected to either a common or to independent phononic baths (see figure 1). We employ Bloch-Redfield-Wangsness (BRW) theory [16,17] to obtain master equations describing the system dynamics under strong coupling. The final master equations only contain a few Lindblad terms, with rates determined by the product of (i) the bath spectral density at the corresponding transition frequency and (ii) algebraic prefactors obtained from transforming the system into its eigenstate basis. Using the specific properties of the system under study, we can evaluate all these prefactors analytically, and are thus able to directly read off the population transfer rates between the system eigenstates. This allows detailed insight into the system dynamics and, specifically, the role of the dark states. Our results demonstrate that for large enough Rabi splitting, the bright state indeed behaves like a single isolated oscillator and the dark states play an almost negligible role. For the case of a common bath, corresponding to long-range bath phonons, this is even true regardless of the Rabi splitting.
2. Theory 2.1. Coherent dynamics We model the system as a set of N molecular vibrational modes coupled to a single electromagnetic (EM) mode in a microcavity, as depicted in figure 1. The mirrors in the experiments [6,7] are actually planar, and the photonic modes form a continuum, with a dispersion relation depending on the in-plane momentum ⃗ ∥ k . However, the assumption of a single EM mode coupling to many molecules is justified when comparing the density of EM modes with the molecule density in the experiment. For the jth transversal mode in a cavity of length L with background refractive index n, the momentum perpendicular to the plane mirrors is fixed to = π ⊥ k j L , so that the dispersion in Inserting the cavity mode dispersion gives a density of states with energy ω ω < max of (here and in the following, we set =  1) where ω π = c n j L j . We assume that modes with j = 1 and ω ω ω There are thus on the order of 10 12 molecular vibrational modes coupled to each photonic mode, and our assumption is clearly justified.
Within the rotating wave approximation, i.e., neglecting processes that create or destroy two excitations, the coherent dynamics of the system is then governed by the Hamiltonian: Figure 1. Sketch of the model system. An ensemble of N molecules interacts collectively with a cavity mode of frequency ω c . The cavity is tuned to a molecular vibrational mode. The remaining rovibrational modes are introduced by either connecting all molecules to a common bath (panel (a)) or by coupling each molecule with an independent bath (panel (b)). In both cases the harmonic bath is characterized by a spectral density ω J ( ).
where a is the annihilation operator for the cavity mode with frequency ω c , and c i is the annihilation operator of the optically active vibrational mode of molecule i, characterized by its frequency ω i and its position ⃗ r i . The cavity-oscillator interaction is given by g i , which depends on the cavity electric field strength and the change of the molecular dipole moment under displacement from the equilibrium position (see [6] for details). For simplicity, in this work we assume a regular configuration in which all the molecules are , as well as zero detuning (ω ω = m c ). Using a direct numerical implementation of BRW theory, which works with arbitrary (but small enough) systems, we have explicitly checked that orientational disorder ( = g g i j ) and inhomogenous broadening (ω ω = i j ) do not significantly affect the results presented below.
We focus on the linear response of the system, so that we can restrict the treatment to the zero-and singleexcitation subspaces 3  . Collective strong coupling emerges when Ω R is larger than the losses of the system. We can distinguish three types of loss mechanisms: cavity losses (rate κ), non-radiative internal losses within the molecules (rate γ nr ), and dephasing-type interactions. Spontaneous radiative decay is very slow (on the scale of milliseconds) due to the low transition frequencies and can be safely neglected 4 .

Incoherent dynamics
We now turn to the description of the incoherent dynamics induced by the different loss mechanisms. Whereas cavity losses and nonradiative internal molecular decay can be seen as pure decay channels and included as constant Lindblad terms when analyzing the system dynamics, dephasing-type interactions must be treated in a more detailed fashion, as we detail below. In the case analyzed in this work, these interactions are due to elastic scattering of low-frequency rovibrational bath modes with the main vibrational modes involved in strong coupling, described by the interaction Hamiltonian The spatial extension or localization of the bath modes determines the character of the coupling. In the following, we focus on two limiting scenarios. In the first scenario, all the vibrational modes are coupled to the same common bath (see figure 1(a)), characterized by delocalized phonons In the second scenario, each molecular vibrational mode is coupled to an independent bath characterized by on-site phonons b ik , with bath Hamiltonian In both scenarios, we assume that all molecules are identical (λ λ = ik k and ω ω = ik k ). The frequently used approach of treating dephasing through a frequency-independent Lindblad superoperator acting only on the vibrational modes is not valid in our case as, within the strong coupling regime, the molecule-cavity coupling (frequency Ω R ) is much faster than the correlation time of the phononic environment (the timescale on which bath correlations decay, i.e., within which the bath 'forgets' about state changes). Instead, the influence of the background modes has to be taken into account in the dressed basis obtained after diagonalizing the strong-coupling interaction. The bath is completely characterized by the spectral density ω in our case). In the following, we assume an Ohmic environment with a quadratic cutoff at frequency ω cut where η is a dimensionless constant that determines the system-bath coupling strength.
If the system-bath coupling is sufficiently weak, BRW theory [16,17] can be safely applied to derive a master equation for the system dynamics [20,21]. This approach relies on the first and second Born approximation, which calculate the effect of the system-bath coupling up to second order perturbation theory and assume that the bath state remains unmodified, i.e., the bath density matrix ρ b is time-independent (and thermally populated in the following), where ρ ρ ρ . Additionally, since the decay of the bath correlations ( ω ∼1 cut , see below) occurs on a much shorter time scale than the dynamics caused by the interaction with the bath, the Markov approximation is used, disregarding all memory effects of the system-bath interaction. In the interaction picture (denoted by a tilde, = , the system density operator ρ then evolves according to where Tr b denotes the trace over the bath degrees of freedom. The bath-dependent part of the system-bath coupling is fully encoded in the bath correlation functions between the modes on molecular sites i and j, ϕ τ , while for independent baths, the off-diagonal terms vanish, . The autocorrelation function ϕ τ ( )can be expressed through the spectral density [22] ∫ ϕ τ ω is the Bose occupation factor and β = k T 1 B , with k B the Boltzmann constant and T the temperature. In the high-temperature limit βω ≫ 1 cut , it can be shown that ϕ τ ω , i.e., the bath correlation time is ω ∼1 cut . We now evaluate equation (5) within the system eigenbasis. Before we proceed, we note that the common approach of using frequency-independent Lindblad terms to describe dephasing is equivalent to neglecting the strong coupling in the incoherent dynamics, i.e., to use only the uncoupled system Hamiltonian † † in the interaction picture on the right-hand side of equation (5). When the molecule-cavity coupling is comparable to or faster than the decay of bath correlations Ω ω ≳ ( ) R cut , this is an invalid approximation, and it is crucial to include the full system Hamiltonian when deriving the master equation to satisfy detailed balance [23]. This point was also stressed in a recent paper in the context of electronic strong coupling [24], based on a phenomenological model of the system.
We thus proceed by expressing the system part of ϕ H (∝c c i i † ) in terms of the dressed eigenbasis (where ) and inserting this expansion in equation (5). This leads to , and the sums over p, q, r and s include all system eigenstates. Furthermore, give the overlaps between system eigenstates and vibrational mode excitations and can be chosen real. Finally, we made the substitution τ ′ = − t t . Inserting ϕ τ ( )from equation (6) leads to integrals of the type where P.V. denotes the Cauchy principal value-we neglect these imaginary parts as they only induce small energy shifts (Lamb shifts) that can be reabsorbed in the coherent dynamics, and arrive to . Note that exists.
In the resulting master equation, the terms in which ω pq differs from ω sr oscillate as a function of time t, which can lead to a violation of the positivity of ρ for long times. If the timescale τ ϕ of the bath-induced system dynamics is much slower than the coherent dynamics, i.e., τ τ ≫ ϕ SC , these terms can be removed by averaging the master equation over a time short compared to τ ϕ , but long compared to τ SC [25,26]. The total bare-molecule width γ = 3.2 meV [6] presents an upper bound for the bare-molecule dephasing rate γ γ ⩽ ϕ , and thus gives a lower bound for τ γ ∼ ϕ ϕ 1 . As a consequence, only the secular terms ω ω = pq sr persist. This secular approximation 5 results in a master equation where populations and coherences are decoupled. We note that the secular approximation aids in the interpretation of the different terms that are obtained, but is not required and indeed only used in some of the results shown in the following.
The secular terms are enumerated in table 1, which lists the states pq rs { , }connected by transitions with frequency ω sr . For simplicity, we define and label rates associated with each transition frequency, which correspond to twice the noise power spectrum evaluated at ω sr , such that γ ω = S 2 ( ) rs sr . These are 'bare' rates in the sense that they do not contain the algebraic prefactors from the basis transformation. The factor of two is included so that for the terms where pq = sr, we obtain a standard Lindblad term with rate γ rs , i.e.
sr rs sr rs is a standard Lindblad superoperator. Positive frequencies ω > 0 sr correspond to phonon emission where the system transitions from a higher-to a lower-energy state, while negative frequencies ω < 0 sr correspond to phonon absorption. As shown in table 1, we obtain secular terms connecting the two polaritons (γ ), terms connecting the polaritons with the dark , and terms connecting states with the same energy (γ = ϕ S 2 (0), equal to the bare-molecule dephasing rate). The latter give pure dephasing for the polaritons | + 〉, | − 〉, but produce coupling between populations and coherences for all dark states.
The final master equations are obtained by using the properties of the basis transformation matrix u ia to evaluate the algebraic prefactors ∑ u u u u i ip iq ir is (independent baths) and ∑ u u u u i j ip iq jr js , (common bath) in equation (9). Specifically, we use that i) the polariton-vibrational mode overlaps are given by , and iii) that dark states are orthogonal to the Note that i and j in the sums are molecule indices (not including the cavity mode), so that , is generally not true. This procedure gives the final secularized master equation for the density operator, given in the Schrödinger picture below. For a common bath with ϕ τ ϕ τ = ( ) ( ) ij , we find that many terms vanish because of the orthogonality relations, giving correspond to incoherent excitation transfer between system eigenstates and are depicted schematically in figure 2(a). A phonon of frequency Ω R may be emitted transferring an excitation from the upper polariton | + 〉 to the lower Label The secular approximation is also sometimes called the rotating wave approximation, but should not be confused with the rotating wave approximation performed in the system Hamiltonian in equation (2). 6 The same γ γ , a e are obtained under strong classical driving of a single two-level system [27]. polariton | − 〉, with rate γ 4 e . Phonon absorption occurs analogously, with characteristic rate γ 4 a . Furthermore, the polaritons undergo pure dephasing with rate γ ϕ 4. We note that the factor of 1 4 in these rates can be easily understood: the bare-molecule dephasing interaction (rate γ ϕ ) is reduced by a factor of two for the polaritons, which consists of an equal superposition of molecules and the cavity mode. Since the interaction in the dressed basis couples each polariton not only with itself, but also with the other polariton, the prefactors are reduced by another factor of two. Taking into account the spectral density at the energy difference between the polaritons modifies the rates for transitions between different polaritons (giving γ a and γ e instead of γ ϕ ), but the prefactor of 1 4 remains. Finally, the last term in equation (12) corresponds to baremolecule dephasing for a common bath, but projected into the degenerate dark-state subspace (using Remarkably, for the case of a common bath, i.e., long-range bath phonons, the dark states are completely decoupled from the polaritons and the bright state behaves identically to a single oscillator interacting with the cavity field. Turning to the case of independent baths, ϕ τ  and Γ 2 e , as shown in figure 2(b). On the other hand, pure dephasing of the polaritons and direct transitions between them through absorption or emission of phonons of frequency Ω R play a negligible role within this dephasing scenario, as their rates scale as N 1 . Additional terms couple between different polariton-dark state coherences (equations (13c), (13d)), without affecting the populations. Finally, the second term in equation (13e) again corresponds to a bare-molecule Lindblad dephasing term that has been restricted to act only within the degenerate darkstate subspace. For → ∞ N , equation (13) simplifies to   (panel (b)). The arrows indicate decay processes while wavy lines represent elastic events that produce dephasing. The fuzzy halo around the dark states indicates excitation transfer and dephasing acting within this manifold.
As is clear from the expressions for the different decay rates, Ω R is the parameter that controls the interaction between the dressed excitation and the phonon bath. Hence the importance of these decoherence mechanisms is dictated by the spectral density evaluated at Ω R or Ω 2 R . This is in clear contrast to what we would obtain by treating the interaction ϕ H in the uncoupled basis. In that case, the Markov approximation would have resulted (independent baths). In both cases, the dephasing interaction would be totally controlled by just the zero-frequency limit of the spectral density, resulting in an overestimation of the rates γ a e , and Γ a e , . This fact emphasizes the key importance of deriving Lindblad terms in the strongly coupled basis when considering non-flat reservoirs [23].
Apart from the decay and dephasing mechanisms induced by the dephasing-type interactions, which conserve the number of excitations, the excited states may decay through nonradiative molecular decay and cavity losses into the ground state | 〉 0 . These decay process, which we have neglected in the theory up to now, could be described within the same framework by dissipative coupling with thermal baths (e.g., the photonic modes outside the cavity). However, as the energy shifts induced by the strong coupling are small compared to the transition frequency (Ω ω ≪ R m ), we can include them through Lindblad terms for the bare cavity and molecular vibrational modes. After removing nonsecular terms in the eigenstate basis, this gives new Lindblad terms with decay rates of γ κ + ( 2 2) nr for the two polaritons and γ nr for the dark states. We can now obtain the total decay rates by collecting terms that transfer excitations out of each state (i.e., excluding pure dephasing terms), giving where Γ + , Γ − , and  Γ are the decay rates of the upper polariton, lower polariton, and dark states, respectively. The intrinsic lifetimes for each state in the two bath scenarios are given by τ While the decay of the polaritons is dominated by cavity losses (κ 2) for realistic parameters, the upper polariton does have a slightly shorter lifetime than the lower polariton in both scenarios, as Γ Γ > e a and γ γ > e a . The exact amount of asymmetry depends on the temperature T and Rabi splitting Ω R . We note that in contrast to the conceptionally similar model for electronic strong coupling of Canaguier-Durand et al [24], our model does not predict a polariton lifetime that is much longer than the cavity mode lifetime.

Results
In the following, we apply our theoretical framework to the experimental results of Shalabney et al [6]. For the vibrational mode of the bare molecules, they report a linewidth of γ = 3.2 meV, with negligible inhomogeneous broadening. This linewidth has contributions from nonradiative decay and dephasing, γ γ γ = + ϕ nr , which can not be distinguished in the absorption spectrum. Although direct information about the weights of nonradiative and dephasing channels is thus not available, dephasing is in general expected to provide a significant contribution for vibrational transitions [14,15]. Therefore, in our calculations we will use a factor ] to measure the relative importance of the two channels. In this way, the factor η in equation (4), which quantifies the strength of the system-bath coupling, is simply given by η K. The cut-off frequency for the thermal bath of low-frequency rovibrational excitations is chosen as ω = 6 cut meV, corresponding to the range of low-frequency phonon modes in the system [6]. Finally, we use a cavity loss rate of κ = 17 meV as estimated in [6] through fitting of the transmission spectrum in the strong coupling regime.
We calculate the absorption spectra by introducing a weak driving term, , which coherently pumps the cavity mode. The density matrix in the steady state, ρ ss , can be calculated in the frame rotating with the driving frequency ω. In this frame, H d is time-independent, but the system frequencies are shifted and the density matrix ρ ss depends on ω. The absorption spectrum is then obtained ss For the numerical implementation, we employ the open-source QuTiP package [28]. This also implies that for both bath scenarios and within the experimental conditions reported in [6,7], the ensemble of vibrational modes behaves effectively as just a single collective molecular oscillator, the bright state, coupled to the cavity mode. Consequently, the dark modes are effectively decoupled from the system dynamics under external driving. We note here that the fit used to extract the cavity width κ from the experimental data in [6] is performed under strong coupling (using a transfer matrix method). Thus, the change of linewidths predicted by our model would already be present in the observed spectra, and is consequently absorbed in the extracted cavity linewidth. This unfortunately precludes a direct test of our model based on the experimental absorption data, for which the bare cavity linewidth would need to be available.
Furthermore, it is interesting to note that in the case of individual baths, the dephasing contributions to the polariton modes are completely suppressed, analogous to the well-known suppression of inhomogeneous broadening under strong coupling [29]. This can be understood by the fact that dephasing and inhomogeneous broadening are closely related, corresponding respectively to temporal fluctuations or to a static distribution of the transition energies.
The dark modes could play a bigger role in a situation with smaller Rabi splitting, for which, in order to still achieve strong coupling, the cavity losses also need to be reduced as compared to the experiments (e.g., by using thicker mirrors). The resulting absorption spectra are shown in figure 3(b), for the same parameters as in figure 3(a) but now with κ = 1 meV and Ω = 6.5 R meV. In both bath scenarios, a slight asymmetry between upper and lower polariton is now noticeable, as the emission of phonons from the upper polariton γ Γ ω ∝ + n ( , ( ) 1) e e is more likely than the absorption of phonons in the lower polariton γ Γ ω ∝ n ( , a a . However, as the involved transition frequencies are smaller than the thermal energy k ( B = = T T K 25.9 meV for 300 ), the thermal occupation ω n ( )is significant and the rates for phonon emission and absorption are comparable. Furthermore, although the phenomenology in the observed absorption spectra is quite similar in both dephasing scenarios, with a reduction of the polariton linewidths for increasing dephasing f, the underlying physics are now quite distinct. This is demonstrated clearly when inspecting the population dynamics, shown in figure 4 for the case of the upper polariton being initially excited. In both cases, the dynamics are dominated initially by the fast decay of the upper polariton through cavity losses (rate κ 2, cf equation (15)). However, the more interesting aspect is the dynamics within the excited subspace: for the common bath, shown in figure 4(a), the dark states are completely decoupled from the dynamics and population is only transferred between the polaritons. In contrast, in the case of independent baths, shown in figure 4(b), the excited population is quite efficiently transferred from the upper polariton to the dark states, resulting in a fast rise of the dark state population. The dark states subsequently decay through nonradiative molecular loss (γ nr , cf equation (15)), which is smaller than the polariton loss rate. We also note that the small but visible population transfer to the lower polariton only shows up because we show results for N = 100, not in the limit → ∞ N (cf equation (14)).

Conclusions
In summary, by using a fully quantum framework we have studied in detail the phenomenon of collective strong coupling when an ensemble of molecular vibrational modes interacts with a cavity EM mode, as realized experimentally in two recent papers [6,7]. We have demonstrated that dephasing-type interactions with a thermal bath of background modes in such a system have to be treated beyond the usual Lindblad approximation in order to represent the effects of the spectral density of bath modes correctly. We have investigated two 'extreme' scenarios for the bath, with either a common bath for all molecules, or independent baths for each molecule. For the experimentally relevant parameters, we find that the dark modes are almost totally decoupled from the polaritons in both scenarios, and the bright state behaves almost like a single isolated oscillator. Our findings thus suggest that this type of system is an ideal and simple platform to explore the exciting possibilities of cavity optomechanics at room temperature.