Langevin Approach to Quantum Optics with Molecules

TL;DR

This paper develops a quantum Langevin equations framework to analytically study light-molecule interactions, revealing effects on emission, cavity transmission, and energy transfer in quantum optical systems with vibrational modes.

Contribution

It introduces an analytical approach to model molecular quantum optics, including emission profiles, cavity effects, and energy transfer, advancing understanding of light-molecule interactions.

Findings

Analytical expressions for emission branching ratios in the Purcell regime.

Characterization of cavity transmission asymmetries in strong coupling.

Insights into F"{o}rster energy transfer mediated by vacuum or cavity modes.

Abstract

We investigate the interaction between light and molecular systems modeled as quantum emitters coupled to a multitude of vibrational modes via a Holstein-type interaction. We follow a quantum Langevin equations approach that allows for analytical derivations of absorption and fluorescence profiles of molecules driven by classical fields or coupled to quantized optical modes. We retrieve analytical expressions for the modification of the radiative emission branching ratio in the Purcell regime and for the asymmetric cavity transmission associated with dissipative cross-talk between upper and lower polaritons in the strong coupling regime. We also characterize the F\"{o}rster resonance energy transfer process between donor-acceptor molecules mediated by the vacuum or by a cavity mode.

Figures

Equations

$$h_0^{(j)}=\left({\omega }_{\text{e}}+\sum _{\text{k}=1}^n{\lambda }_{\text{k}}^2{\nu }_{\text{k}}\right){\sigma }_j^{†}{\sigma }_j+\sum _{\text{k}=1}^n{\nu }_{\text{k}}{b}_{j\text{k}}^{†}{b}_{j\text{k}}.$$

$$h_H^{(j)}=-\sum _{\text{k}=1}^n{\lambda }_{\text{k}}{\nu }_{\text{k}}\left(b_{j\text{k}}^{†}+b_{j\text{k}}\right){\sigma }_j^{†}{\sigma }_j,$$

$$H_c={\omega }_{\text{c}}{a}^{†}a+i\sum _{j=1}^N\left(g_j{a}^{†}{\sigma }_j-g_j^{\ast }{\sigma }_j^{†}a\right).$$

$$H=\sum _{j=1}^N\left(h_0^{(j)}+h_H^{(j)}\right)+\sum _{i\ne j}^N{h}_{dd}^{(ij)}+H_c+H_{\ell }.$$

$$\frac{d{b}_{\text{k}}(t)}{dt}$$

$$\frac{d\tilde{\sigma }(t)}{dt}$$

$$\sigma (t)$$

$$+\sqrt{2\gamma }{\int }_0^{t}d{t}^{\prime }{e}^{-\left[\gamma -i({\omega }_{\ell }-{\omega }_{\text{e}})\right](t-t^{\prime })}{\Sigma }_{\text{in}}(t^{\prime })\mathcal{D}(t){\mathcal{D}}^{†}(t^{\prime }).$$

$$P_{\text{e}}^{\text{ss}}={\eta }_{\ell }^2\sum _{\left\{m_{\text{k}}\right\}=0}^{\infty }\frac{{\text{s}}_{m_1}^{{\lambda }_1}...{\text{s}}_{m_n}^{{\lambda }_n}\left[\gamma +{\Gamma }_{\left\{m_{\text{k}}\right\}}\right]/\gamma }{{\left[\gamma +{\Gamma }_{\left\{m_{\text{k}}\right\}}\right]}^2+{\left[{\omega }_{\ell }-{\omega }_{\text{e}}-{\nu }_{\left\{m_{\text{k}}\right\}}\right]}^2}.$$

$$S_{\text{em}}(\omega )=\sum _{\left\{m_{\text{k}}\right\}=0}^{\infty }\frac{2{\text{s}}_{m_1}^{{\lambda }_1}...{\text{s}}_{m_n}^{{\lambda }_n}\left[\gamma +{\Gamma }_{\left\{m_{\text{k}}\right\}}\right]}{{\left[\gamma +{\Gamma }_{\left\{m_{\text{k}}\right\}}\right]}^2+{\left[\omega -{\omega }_{\text{e}}+{\nu }_{\left\{m_{\text{k}}\right\}}\right]}^2}.$$

$$\text{M}=\left(\begin{array}{cc}-\left[\gamma -i({\omega }_{\ell }-{\omega }_{\text{e}})\right]& g\\ -g& -\left[\kappa -i({\omega }_{\ell }-{\bar{\omega }}_{\text{c}})\right]\end{array}\right),$$

$$v(t)={\int }_{-\infty }^{t}d{t}^{\prime }\text{T}{e}^{{\Lambda }_M(t-t^{\prime })}{\text{T}}^{-1}\mathcal{D}(t){\tilde{v}}_{\text{in}}(t^{\prime }).$$

$${\text{t}}_{\text{c}}$$

$$-\frac{{\text{T}}_{12}{\text{T}}_{21}}{\det \text{T}}\sum _{\left\{m_k\right\}=0}^{\infty }\frac{\kappa {\text{s}}_{m_1}^{{\lambda }_1}...{\text{s}}_{m_n}^{{\lambda }_n}}{{\gamma }_{+}+{\Gamma }_{\left\{m_{\text{k}}\right\}}-i\left({\omega }_{\ell }-{\omega }_{+}-{\nu }_{\left\{m_{\text{k}}\right\}}\right)}.$$

$${\kappa }_{\text{UL}}=\frac{1}{2}\frac{{\lambda }^2{\nu }^2\Gamma }{{\Gamma }^2+{\left({\omega }_{+}-{\omega }_{-}-\nu \right)}^2}.$$

$${\alpha }_{\text{cav}}=\frac{(1+C_{\text{00}})e^{-{\sum }_{\text{k}=1}^n{\lambda }_{\text{k}}^2}}{1+C_{\text{00}}{e}^{-{\sum }_{\text{k}=1}^n{\lambda }_{\text{k}}^2}},$$

$${\kappa }_{\text{ET}}=\sum _{n_{\text{A}},n_{\text{D}}}\frac{2{\text{s}}_{n_{\text{D}}}^{{\lambda }_{\text{D}}}{\text{s}}_{n_{\text{A}}}^{{\lambda }_{\text{A}}}(n_{\text{A}}{\Gamma }_{\text{A}}+n_{\text{D}}{\Gamma }_{\text{D}})\Omega (r{)}^2}{(n_{\text{A}}{\Gamma }_{\text{A}}+n_{\text{D}}{\Gamma }_{\text{D}}{)}^2+(\Delta -n_{\text{A}}{\nu }_{\text{A}}-n_{\text{D}}{\nu }_{\text{D}}{)}^2}.$$

$$H=\left[{\omega }_e+\frac{p^2}{2\mu }+\frac{1}{2}\mu {\nu }^2{\left(R-R_e\right)}^2\right]{\sigma }^{†}\sigma +\left[\frac{p^2}{2\mu }+\frac{1}{2}\mu {\nu }^2{\left(R-R_g\right)}^2\right]\sigma {\sigma }^{†}.$$

$$H=\frac{P^2}{2\mu }+\frac{1}{2}\mu {\nu }^2{Q}^2+{\omega }_e{\sigma }^{†}\sigma -\mu {\nu }^{2}Q{R}_{\text{ge}}{\sigma }^{†}\sigma +\frac{1}{2}\mu {\nu }^2{R}_{\text{ge}}^2{\sigma }^{†}\sigma .$$

$$H=\nu b^{†}b+({\omega }_{\text{e}}+{\lambda }^2\nu ){\sigma }^{†}\sigma -\lambda \nu (b^{†}+b){\sigma }^{†}\sigma .$$

$$\dot{\mathcal{O}}=i\left[H,\mathcal{O}\right]-\left[\mathcal{O},c^{†}\right]\left\{{\gamma }_{c}c-\sqrt{2{\gamma }_c}{c}_{\mathrm{in}}\right\}+\left\{{\gamma }_c{c}^{†}-\sqrt{2{\gamma }_c}{c}_{\mathrm{in}}^{†}\right\}\left[\mathcal{O},c\right],$$

$$\frac{db}{dt}$$

$$\frac{d\tilde{\sigma }}{dt}$$

$$\frac{d{\sigma }^{†}\sigma }{dt}$$

$$\frac{db}{dt}$$

$$\frac{d\tilde{\sigma }}{dt}$$

$$\frac{d{\sigma }^{†}\sigma }{dt}$$

$$\lambda \Gamma ⟨(b^{†}-b)\tilde{\sigma }⟩-\sqrt{2\Gamma }\lambda ⟨(b_{\text{in}}^{†}-b_{\text{in}})\tilde{\sigma }⟩=-{\lambda }^2\Gamma \sigma e^{-{\lambda }^2/2},$$

$$b(t)=\sqrt{2\Gamma }{\int }_{-\infty }^{t}d\tau e^{-(\Gamma +i\nu )(t-\tau )}{b}_{\text{in}}(\tau ).$$

$$b^{†}(t)=\sqrt{2\Gamma }{\int }_{-\infty }^{t}d\tau e^{-(\Gamma -i\nu )(t-\tau )}{b}_{\text{in}}^{†}(\tau ),$$

$$\displaystyle[b(t),b^{\dagger}(t^{\prime})]=\left\{\begin{array}[]{cc}e^{-(\Gamma+i\nu)(t-t^{\prime})}&t\geq t^{\prime}\\ e^{-(\Gamma-i\nu)(t^{\prime}-t)}&t^{\prime}\geq t\end{array}\right..$$

$$F(\omega )=2\Re \left\{{\int }_0^{\infty }d\tau e^{i\omega \tau }⟨F(\tau )F(0)⟩\right\}.$$

$$⟨b(t)b^{†}(t^{\prime })⟩=2\Gamma {\int }_{-\infty }^{t}d\tau {\int }_{-\infty }^{t^{\prime }}d{\tau }^{\prime }{e}^{-(\Gamma +i\nu )(t-\tau )}{e}^{-(\Gamma -i\nu )(t^{\prime }-{\tau }^{\prime })}\delta (\tau -{\tau }^{\prime })=e^{-(\Gamma +i\nu )(t-t^{\prime })},$$

$$F(\omega )=2{\lambda }^2{\nu }^2\Re \left\{\frac{1}{\Gamma +i(\omega -\nu )}\right\}=\frac{2{\lambda }^2{\nu }^2\Gamma }{{\Gamma }^2+(\omega -\nu {)}^2}.$$

Full text

Recent experimental progress Wang et al. (2019) has shown that the Purcell enhancement of the zero-phonon line of a single molecule can strongly alter the branching ratio of spontaneous emission between the line of interest and additional Stokes lines thus turning the molecule into an ideal quantum emitter. At the mesoscopic level, experiments in the collective strong coupling regime of organic molecules with cavities or delocalized plasmonic modes have shown important light-induced modifications of material properties. Experimental and theoretical endeavors go into the direction of charge and energy transport enhancement Orgiu et al. (2015); Schachenmayer et al. (2015); Feist and Garcia-Vidal (2015); Hagenmüller et al. (2017, 2018), Förster resonance energy transfer (FRET) enhancement Zhong et al. (2016, 2017); Garcia-Vidal and Feist (2017); Du et al. (2018); Reitz et al. (2018), modified chemical reactivity Hutchison et al. (2012); Galego et al. (2017); Herrera and Spano (2016); Martínez-Martínez et al. (2018); Kampschulte and Denschlag (2018), polariton dynamics Schwartz et al. ; Canaguier-Durand et al. (2015) etc. Oftentimes however, experiments rely on theoretical models developed for standard cavity quantum electrodynamics (cavity QED) Haroche and Kleppner (1989); Berman (1994); Walther et al. (2006) with two-level systems where one distinguishes between i) the Purcell regime, characterized by modifications of the spontaneous emission rates and ii) the strong coupling regime leading to the occurrence of hybrid light-matter states referred to as polaritons. Recent theoretical efforts aim at covering this gap by solving a generalized light-electronic-vibrations problem modeled as a Holstein-Tavis-Cummings Hamiltonian. Investigations aim at providing an understanding of the vibrationally induced cavity polariton asymmetry Canaguier-Durand et al. (2015); Neuman and Aizpurua (2018), vibrationally dressed polaritons Zeb et al. (2018), dark vibronic polaritons  Herrera and Spano (2017, 2018), developing a cavity Born-Oppenheimer theory Flick et al. (2017a, b) or deriving relevant simplified models for large scale numerics in the mesoscopic limit del Pino et al. (2018).

We provide here an alternative path based on solving the Holstein-Tavis-Cummings dynamics at the level of operators rather than states. The basic model considers a molecular box (see Fig. 1(b)) comprised of an internal electronic transition coupled to any number of vibrational modes. Radiative decay and vibrational relaxation are included as stochastic source terms in a set of coupled standard quantum Langevin equations Gardiner and Zoller (2004); Gardiner and Collett (1985); Collett and Gardiner (1984); Breuer and Petruccione (2002); Walls and Milburn (2006) for vibrational and polaron operators (similarly applied in optomechanical systems Genes et al. (2008); Rabl (2011); Nunnenkamp et al. (2011)). The method is numerically tested showing correct predictions for molecular emission/absorption lines. In the strong coupling regime of cavity QED, polaron-polariton dynamics provides insight in the upper-lower polariton dissipative coupling. In the Purcell regime we analytically derive the modified branching ratio of radiative transitions between the zero-phonon line and all other Stokes lines Wang et al. (2019). For different molecular species, we derive the donor–acceptor FRET rate occurring either via direct dipole-dipole virtual photon exchanges or in the presence of a cavity. We find regimes where the direct transfer can be increased by a Purcell-enhancement of the acceptor’s absorption lines.

Model — We consider $N$ molecules ($j$ running from $1$ to $N$) with ground $\ket{g_{j}}$ and excited $\ket{e_{j}}$ electronic states and Pauli lowering operators $\sigma_{j}=\ket{g_{j}}\bra{e_{j}}$. Each molecule’s excited state orbital is perturbed by its coupling to $n$ vibrational modes (frequencies $\nu_{\text{k}}$), corresponding to stretching along different nuclear coordinates. These are described by creation/annihilation operators $b^{\dagger}_{j\text{k}}$ and $b_{j\text{k}}$ ($k$ running from $1$ to $n$). The vibronic couplings (see Supplemental Material (SM) SM ) are modeled via the standard Holstein Hamiltonian Holstein (1959) (with $\hbar=1$) with coupling parameters $\lambda_{\text{k}}=\mu_{\text{k}}\nu_{\text{k}}R^{(k)}_{\text{ge}}r^{(k)}_{\text{zpm}}$ (the square root of the Huang-Rhys factor) where the spread of the zero point motion wavepacket is $r^{(k)}_{\text{zpm}}=\sqrt{1/(2\mu_{\text{k}}\nu_{\text{k}})}$ (with $\mu_{\text{k}}$ – the effective mass for the given k-vibration). The free dynamics of a single molecule is then governed by

[TABLE]

Holstein contributions from all vibrational modes add as

[TABLE]

and lead to the standard Franck-Condon physics illustrated in Fig. 1(a). Notice that for each molecule one could define a single collective vibrational bright mode that couples to the electronic degrees of freedom and diagonalize the free vibrational part in terms of dark modes orthogonal to it and bright-dark couplings. In some cases such chain transformation (as used in Ref. del Pino et al. (2018)) can provide a great simplification to the problem especially for numerical simulations. This is not necessary in our treatment as we can include all vibrational modes on equal footing at the analytical level. The molecules can be embedded in a cavity (see Fig. 1(a)) with a single mode of interest - annihilation operator $a$ with frequency $\omega_{\text{c}}$. Each molecule is coupled with $g_{j}$ to the cavity field; the interaction is the standard Tavis-Cummings Hamiltonian

[TABLE]

At high densities, molecules sit in the near-field of their neighbors which leads to vacuum mediated dipole-dipole exchanges modeled as $h^{(ij)}_{dd}=\Omega_{ij}\sigma_{i}\sigma^{\dagger}_{j}$ (with a strongly distance dependent exchange rate $\Omega_{ij}$). In addition, either the molecules or the cavity can be driven by a laser at frequency $\omega_{\ell}$ with pump amplitude $\eta_{\ell}$ (or $\eta_{\text{c}}$ - for cavity pumping). Adding all the parts together, the system coherent dynamics is governed by

[TABLE]

Laser driving is included as $H_{\ell}=i\eta_{\ell}\sum_{j=1}^{N}\left(\sigma^{\dagger}_{j}e^{-i\omega_{\ell}t}-\sigma_{j}e^{i\omega_{\ell}t}\right)$ (for molecule pumping with amplitude $\eta_{\ell}$) or as $H_{\ell}=i\eta_{\text{c}}\left(a^{\dagger}e^{-i\omega_{\ell}t}-ae^{i\omega_{\ell}t}\right)$ (with $\eta_{\text{c}}$ for cavity pumping). Dissipation is added as Lindblad terms: radiative decay (at rates $\gamma_{j}$ with collapse operators $\sigma_{j}$), cavity field decay (rate $\kappa$, collapse operator $a$) and internal vibration relaxation (rates $\Gamma_{jk}$, collapse operators $b_{jk}$). For a collapse operator $\cal{O}$ with rate $\gamma_{\cal{O}}$ the Lindblad term applied to a density operator $\rho$ is: ${\cal{L}}_{\cal{O}}[\rho]=\gamma_{\cal{O}}\left\{2\cal{O}\rho\cal{O}^{\dagger}-\rho\cal{O}^{\dagger}\cal{O}-\cal{O}^{\dagger}\cal{O}\rho\right\}$. For the optical cavity field the Lindblad term is exact even in the zero temperature limit. For vibrations, this is an approximation as the physics of bulk-vibrations interactions could be better included as a Brownian noise dissipation model Hu et al. (1992). We compare the two dissipation models in the SM SM showing agreement for spectral line frequencies and identifying terms leading to disagreement in linewidth predictions.

Single molecule spectroscopy— Setting $g=0$ and assuming a single molecule $N=1$ coupled to many vibrational modes $n>1$, the molecular Hamiltonian can be diagonalized by a standard (level dependent) polaron transformation $U^{\dagger}=({\cal{D}^{\dagger}})^{\sigma^{\dagger}\sigma}$. This is a collective transformation as ${\cal{D}^{\dagger}}={\cal{D}}^{\dagger}_{1}...{\cal{D}}^{\dagger}_{n}$ where each individual displacement is defined as ${\cal{D}}_{\text{k}}=\exp(-i\sqrt{2}\lambda_{\text{k}}p_{\text{k}})=\exp[\lambda_{\text{k}}(b^{\dagger}_{\text{k}}-b_{\text{k}})]$ with the momentum quadrature $p_{\text{k}}$ as a generator. The driving part is transformed accordingly $\tilde{H}_{\ell}=i\eta_{\ell}\left(\sigma{\cal{D}^{\dagger}}{e^{i\omega_{\ell}t}}-{\cal{D}}\sigma^{\dagger}e^{-i\omega_{\ell}t}\right)$ and describes the generation of a superposition of coherent wavepackets with amplitudes $\lambda_{\text{k}}$ – accompanying stimulated absorption and inverse displacement (with $-\lambda_{\text{k}}$) for stimulated emission. However, we only indirectly make use of this transformation by shifting to a Langevin description where the time dynamics of a polaron operator $\tilde{\sigma}(t)=\sigma(t){\cal{D}^{\dagger}}(t)$ is followed. Under weak driving conditions $\eta_{\ell}\ll\gamma$ the exact Langevin equations (see SM SM ) can be simplified to:

[TABLE]

where we condensed the noise affecting the electronic transition into a term of non-vanishing average $\Sigma_{\text{{in}}}(t)=\eta_{\ell}/\sqrt{2\gamma}+\sigma_{\text{in}}(t)$. For zero temperature $T=0$, the only non-zero noise correlations are $\braket{\sigma_{\text{in}}(t)\sigma_{\text{in}}^{\dagger}(t^{\prime})}=\delta(t-t^{\prime})$. The reservoir model is consistent with models used elsewhere Neuman and Aizpurua (2018): the power spectrum of the vibrational force $F=\lambda\nu(b+b^{\dagger})$ is given by a Lorentzian centered around the vibrational frequency $F(\omega)=2\lambda^{2}\nu^{2}\Gamma/\left[\Gamma^{2}+(\omega-\nu)^{2}\right]$. For simplicity of presentation we will proceed with $\text{T}=0$; however the SM SM presents a general approach to include thermal effects by including a finite T via $\braket{b_{\text{in}}^{\dagger}(t)b_{\text{in}}(t^{\prime})}=\Gamma\bar{n}\delta(t-t^{\prime})$ and $\braket{b_{\text{in}}(t)b_{\text{in}}^{\dagger}(t^{\prime})}=\Gamma(\bar{n}+1)\delta(t-t^{\prime})$ where the thermal occupancy is given by $\bar{n}=[\exp({\hbar\nu/(k_{\text{B}}\text{T})})-1]^{-1}$. The first equation can be integrated to give a solution for ${\cal{D}}(t)$ in terms of vibrational input noises $b_{\text{in}}^{\text{k}}(t)$ and $b_{\text{in}}^{\text{k},\dagger}(t)$. We can then derive the solution for $\tilde{\sigma}(t)=\sigma(t){\cal{D}^{\dagger}}(t)$ and subsequently for the bare Pauli operator

[TABLE]

This expression allows for the derivation of any system property both in the transient regime as well as in steady state. For the absorption (extinction) profile, we evaluate $\braket{\sigma^{\dagger}(t)\sigma(t)}$ in steady state under continuous driving conditions (see SM SM ). As the input noises commute at any time, one is left with the task of evaluating $\braket{{\cal{D}}_{\text{k}}(t){\cal{D}}^{\dagger}_{\text{k}}(t^{\prime})}=\exp(-\lambda^{2}_{\text{k}})\exp\left(\lambda_{\text{k}}^{2}[\exp(-\left(\Gamma_{\text{k}}+i\nu_{\text{k}})(t-t^{\prime})\right)]\right)$. After the subsequent integration the excited state population $P^{\text{ss}}_{\text{e}}=\braket{\sigma^{\dagger}\sigma(t)}$ is given by

[TABLE]

The sum above runs over all the indices $\left\{m_{\text{k}}\right\}\equiv{m_{1},...,m_{n}}$ and the numerator contains a product of Poissonian coefficients $\text{s}_{m_{\text{k}}}^{\lambda_{\text{k}}}=\exp({-\lambda_{\text{k}}^{2}})\lambda_{\text{k}}^{2m_{\text{k}}}/{m_{\text{k}}}!$. The denominator indicates a series of resonances corresponding to frequency up-shifted absorption lines with vibrationally modified linewidths $\Gamma_{\left\{m_{\text{k}}\right\}}=\sum_{\text{k}=1}^{n}m_{\text{k}}\Gamma_{\text{k}}$ and corresponding vibrationally-induced energy shifts $\nu_{\left\{m_{\text{k}}\right\}}=\sum_{\text{k}=1}^{n}m_{\text{k}}\nu_{\text{k}}$ (see Fig. 1(c)). Notice that the absorption profile in the pure two-level system limit is immediately recovered by setting all $\lambda_{\text{k}}=0$, which implies that only terms with all $m_{\text{k}}=0$ contribute such that the expected simple Lorentzian is obtained $\braket{\sigma^{\dagger}\sigma(t)}=\eta_{\ell}^{2}/[\gamma^{2}+(\omega_{\ell}-\omega_{\text{e}})^{2}]$.

To evaluate the fluorescence spectrum (see SM SM ), we assume initial preparation of the molecule in a mixed state between ground and excited electronic states with no vibrations and compute the two-time correlations $\braket{\sigma^{\dagger}(\tau)\sigma(0)}$. The spectrum is defined as $S_{\text{em}}(\omega)=2\Re{\int_{0}^{\infty}d\tau\braket{\sigma^{\dagger}(\tau)\sigma(0)}e^{-i\omega\tau}}$ Auffèves et al. (2008) and is expressed as

[TABLE]

The Stokes lines are present in the resonances of the denominator requiring $\omega=\omega_{\text{e}}-\nu_{\left\{m_{\text{k}}\right\}}$ (see Fig. 1(c) - plotted in a frame rotating with $\omega_{\ell}$). Setting all $m_{\text{k}}=0$ reveals the zero-phonon line strength in the emission proportional to $\text{s}_{0}^{\lambda_{1}}...\text{s}_{0}^{\lambda_{n}}=e^{-\left(\lambda_{1}^{2}+...+\lambda_{n}^{2}\right)}$. Notice that for a molecule with many vibrations, even if all Huang-Rhys factors are small, the collective vibrationally induced reduction of the oscillator strength on the zero-phonon line can be considerably large. Also, the branching ratio of the molecules defined as $\alpha=\gamma_{\text{00}}/(\gamma_{\text{00}}+\gamma^{\prime})$ ($\gamma_{\text{00}}$ is the radiative rate of the zero-phonon line and $\gamma^{\prime}$ represents the sum of rates to all other states in the ground state manifold) can be immediately estimated to $\alpha=e^{-\left(\lambda_{1}^{2}+...+\lambda_{n}^{2}\right)}$.

Cavity-molecule spectroscopy – We set $N=1$ and consider laser driving of the cavity mode instead. The vibrational mode follows the same Langevin equation as in Eq. (5) while the cavity-molecule coupling is described in terms of dressed polaron operators $\tilde{\sigma}(t)=\sigma(t){\cal{D}^{\dagger}}(t)$ and $\tilde{a}(t)=a(t){\cal{D}^{\dagger}}(t)$. We can formulate the dynamics in matrix form (see SM SM ) $\partial_{t}\mathbb{\tilde{v}}(t)=\textbf{M}\mathbb{\tilde{v}}(t)+\mathbb{\tilde{v}}_{\text{in}}(t)$ for the vector of operators $\mathbb{\tilde{v}}(t)=(\tilde{\sigma}(t),\tilde{a}(t))^{\intercal}$. The dynamics is described by the drift matrix (with renormalized cavity frequency $\bar{\omega}_{\text{c}}=\omega_{\text{c}}-\sum_{\text{k}=1}^{n}\lambda_{\text{k}}^{2}\nu_{\text{k}}$)

[TABLE]

and the vector of input noises contains the contributions from the displaced electronic and optical driving noises and is given by $\mathbb{\tilde{v}}_{\text{in}}(t)=(\sqrt{2\gamma}{\cal{D}^{\dagger}}(t)\sigma_{\text{in}}(t),\sqrt{2\kappa}{\cal{D}^{\dagger}}(t)A_{\text{in}}(t))^{\intercal}$ where we condensed here the input of the cavity in a single term $A_{\text{in}}=\eta_{c}/\sqrt{2\kappa}+a_{\text{in}}$. Diagonalization of the M matrix implies that $\textbf{M}=\textbf{T}\Lambda_{M}\textbf{T}^{-1}$ where $\Lambda_{M}$ is a diagonal matrix with $\epsilon_{\pm}=-\left(\gamma_{\pm}+i\omega_{\pm}\right)$ on the diagonal. The T matrix is made of the eigenvector components. A formal solution in steady state for the bare operator vector $\mathbb{v}(t)=\left(\sigma(t),a(t)\right)^{\intercal}$ can be obtained

[TABLE]

In general one is interested in deriving the cavity transmission function as the laser detuning is swept around the cavity resonance which we define as $\text{t}_{\text{c}}=\kappa\braket{a(t)}/\eta_{\text{c}}$. To this end we evaluate the second element of the operator vector and perform a classical averaging which leads to

[TABLE]

where ${\det\textbf{T}}=\text{T}_{11}\text{T}_{22}-\text{T}_{12}\text{T}_{21}$. The poles of the denominator show that the maximum transmission occurs at $\omega_{\pm}$ with linewidths $\gamma_{\pm}$. For $\lambda=0$ (pure quantum emitter case) the strong coupling regime with $g\gg(\kappa,\gamma)$ predicts two eigenfrequencies roughly positioned at $\omega_{\text{e}}\pm g$ with linewidths $(\kappa+\gamma)/2$. For $\lambda>0$, dissipatively-induced polariton cross-talk leads to an asymmetric transmission profile as illustrated in Fig. 1(d). A first order perturbative approach (as derived in SM SM ) including upper-lower polariton exchanges via single vibrations leads to an approximate dissipative rate:

[TABLE]

For $g=\nu/2$, the upper polariton lifetime can be simply casted as $2\kappa^{\prime}=2\kappa[1+\lambda^{2}g^{2}/(\kappa\Gamma)]$. The upper polariton is a Lorentzian with an effectively reduced maximum corresponding to an increased linewidth $(\kappa+\gamma)/2+\kappa_{\text{UL}}/2$.

Another interesting aspect is the modification of the radiative branching ratio in the Purcell regime (experimentally employed to turn a molecule into an efficient quantum emitter Wang et al. (2019)). We analyze the fluorescence spectrum of a partially excited molecule for $\kappa\gg g_{\text{00}}\gg\gamma$ but large cooperativity $C_{\text{00}}=g_{\text{00}}^{2}/(\kappa\gamma)\gg 1$. We denote by $g_{\text{00}}=ge^{-\sum_{\text{k}=1}^{n}\lambda_{\text{k}}^{2}/2}$ the coupling rate between the cavity mode and the zero-phonon line reduced from the total rate $g$ by the reduction in oscillator strength. Assuming very good resolution of the zero-phonon line with $\omega_{\text{c}}=\omega_{\text{e}}$ and $\kappa\ll\nu_{\text{k}}$ the emergence of the Purcell regime is accompanied by $\omega_{+}=\omega_{-}=\omega_{\text{e}}$, $\gamma_{+}\approx\kappa$ but an enlarged $\gamma_{-}=\gamma^{\prime}+\gamma_{\text{00}}(1+C_{\text{00}})\gg\gamma$. While emission into levels other than the zero-phonon line is unperturbed and sums up to $\gamma^{\prime}=\gamma-\gamma_{\text{00}}$, the zero-phonon line decay is strongly enhanced and the cavity-modifed branching ratio can be expressed as

[TABLE]

going to unity in the limit that $C_{\text{00}}\rightarrow\infty$ and reproducing the bare molecule result for $C_{\text{00}}=0$.

Förster resonance energy transfer – We consider two molecules (donor-D and acceptor-A) with frequency mismatch $\omega_{\text{D}}-\omega_{\text{A}}=\Delta\gg\nu_{\text{D,A}}$ [see Fig. 2(a)]. As Fig. 2(b) illustrates, their emission-absorption spectra contain a small overlap which allows an energy transfer from the donor to the acceptor via near field strongly distance dependent dipole-dipole exchanges $\Omega(r)\propto 1/r^{3}$. For $\Omega(r)>\gamma_{D}$, an initially excited donor will lose energy both via radiative emission at rate $2\gamma_{\text{D}}$ as well as via a resonant exchange with the acceptor’s higher excited state vibrational manifold. In a standard regime where $\Omega(r)\ll\Gamma_{\text{A}}$ a uni-directional process emerges as the acceptor quickly relaxes to the zero vibrational state and radiates around its peak emission frequency. We analyze the process in a transient pump-probe scenario where the donor is quickly excited and its natural decay rate is modified by a rate $\kappa_{\text{ET}}$ at which it spills energy into the acceptor via the emission-absorption overlap. We proceed by deriving effective rate equations: $\partial_{t}P_{\text{D}}=-(2\gamma_{\text{D}}+\kappa_{\text{ET}})P_{\text{D}}$ and $\partial_{t}P_{\text{A}}=-2\gamma_{\text{A}}P_{\text{A}}+\kappa_{\text{ET}}P_{\text{D}}$ for acceptor/donor excited state populations $P_{\text{A/D}}$ under the perturbation condition $\Omega(r)\ll\Gamma_{\text{A}}$. From the Langevin equations of motion (see SM SM ) we identify the term responsible for energy migration as $-2\Omega(r)\Im\braket{\sigma^{\dagger}_{\text{A}}\sigma_{\text{D}}}\approx-\kappa_{\text{ET}}P_{\text{D}}(t)$ and estimate (for simplicity we set $\gamma_{\text{D}}=\gamma_{\text{A}}$) the Förster resonance energy transfer as

[TABLE]

The scaling with distance is straightforward as $\Omega(r)^{2}\propto r^{-6}$. Optimization of the transfer rate requires maximizing the Poissonian coefficients while bringing the Lorentzian into resonance. This is none other than the usual overlap integral widely discussed in literature but expressed in terms of discrete emission/absorption lines. The resonance transfer occurs when emission meets absorption i.e., fluorescence lines of the donor $\omega_{\text{D}}-n_{\text{D}}\nu_{\text{D}}$ overlap absorption lines of the acceptor $\omega_{\text{A}}+n_{\text{A}}\nu_{\text{A}}$. The formalism allows for extensions to many vibrational modes, temperature effects as well as for the inclusion of dissipative collective effects (subradiance/superradiance) which have been shown to play a role in enhancing the energy transfer rate between quantum emitters Reitz et al. (2018). The pump-prope dynamics is illustrated in Fig. 2(c) where the donor evolution sees an accelerated rate in the presence of a nearby placed acceptor (line tangent show the convergence of the early time evolution with the analytically estimated transfer rate).

In addition to the direct dipole-dipole interaction cavity-mediated processes via couplings $g_{\text{D,A}}$ between donor/acceptor and mode $a$ participate in the energy transfer. We focus on deriving scaling laws in the bad-cavity regime (large cooperativity) via the same perturbative method applied in the case of the direct Förster energy transfer (see SM SM ): this leads to an extra term in the acceptor population rate equation $\partial_{t}P_{\text{A}}=-2\gamma_{\text{A}}P_{\text{A}}+\kappa_{\text{ET}}(r)P_{\text{D}}+\kappa_{\text{ET}}^{\text{cav}}P_{\text{D}}$. The first term contains a cavity-modified dipole-dipole energy transfer rate springing from $2\Omega(r)\Im\braket{\sigma^{\dagger}_{\text{A}}\sigma_{\text{D}}}$. The second term includes purely cavity mediated processes coming from $2g_{A}\Im\braket{\sigma^{\dagger}_{\text{A}}a}$. While the analytical expressions are cumbersome, an immediate result emerges which implies an acceleration of the transfer rate when the cavity is matched to the peak absorption line of the acceptor roughly at $\omega_{c}=\omega_{\text{A}}+n_{\text{A}}\nu_{\text{A}}$ with $n_{\text{A}}$ around $\lambda_{\text{A}}^{2}\nu_{\text{A}}$. Further optimization requires the particular details of the donor/acceptor emission/absorption lines and their relevant vibrational levels involved in the transfer.

Conclusions - A treatment based on quantum Langevin equations allows for analytical insight into the role of vibrations in cavity QED with molecules. Its strength lies in the inclusion of an arbitrary number of vibrations in thermal equilibrium at any T. We have derived spectroscopic quantities and presented general expressions for energy transfer rates between donor-acceptor molecules outside and inside cavities as well as a perturbative expression for dissipative polariton coupling. The approach can be extended to include the dynamics of bulk phonons and their coupling to vibrations leading to vibrational relaxation and phonon wings. Future endeavors will also aim at treating mesoscopic systems in the weak excitation limit with the input-output formalism of Ref. Plankensteiner et al. (2019) providing a quantum description of the cavity output field. This could allow for efficient numerical simulations of the collective strong coupling regime of cavity QED with molecules, identifying the role of dipole-dipole interactions, super/subradiant and collective dark states.

Acknowledgments - We acknowledge financial support from the Max Planck Society and from the German Federal Ministry of Education and Research, co-funded by the European Commission (project RouTe), project number 13N14839 within the research program ”Photonik Forschung Deutschland” (C. S. and C. G.). We acknowledge fruitful discussions with T. Ebbesen, C. Genet, V. Sandoghdar, S. Götzinger, D. Martin-Cano and B. Gürlek.

Appendix A Ad-hoc derivation of the Holstein model for a single nuclear coordinate

Let us justify the form of the Holstein Hamiltonian Holstein (1959) by following a first-principle derivation for a single nuclear coordinate $R$ of effective mass $\mu$. We assume that, along the nuclear coordinate, the equilibria for ground (coordinate $R_{\text{g}}$, state vector $\ket{g}$) and excited (coordinate $R_{\text{e}}$ and state vector $\ket{e}$) electronic orbitals are different. Notice that, for the simplest case of a homonuclear diatomic molecule made of nuclei each with mass $m$, a single vibrational mode exists corresponding to the relative motion mode with effective mass $\mu=m/2$; the equilibria coordinates $R_{\text{g}}$ and $R_{\text{e}}$ correspond to the bond length in ground/excited states. Assuming that both the ground and the excited state have the same parabolic shape around the minima (with second derivative leading to a vibrational frequency $\nu$) one can write the total Hamiltonian:

[TABLE]

Notice that the Hamiltonian is written in a Hilbert space spanning both electronic dynamics (via the Pauli operators) and mechanical dynamics. Introducing small oscillations around the equilibria $Q=R-R_{\text{g}}$ and subsequently $R-R_{\text{e}}=Q+R_{\text{g}}-R_{\text{e}}=:Q-R_{\text{ge}}$ we obtain

[TABLE]

The last term is a renormalization of the bare electronic transition frequency energy (absent in some theoretical treatments) which will naturally go away when diagonalizing the Hamiltonian via the polaron transformation resulting in $\omega_{\text{e}}$ as the natural electronic transition. We can now rewrite the momentum and position operators in terms of bosonic operators $Q=r_{\text{zpm}}(b^{\dagger}+b)$, $P=ip_{\text{zpm}}(b^{\dagger}-b)$ by introducing the zero-point-motion $r_{\text{zpm}}=1/\sqrt{2\mu\nu}$ and $p_{\text{zpm}}=\sqrt{\mu\nu/2}$. Reexpressing the terms above via the Huang-Rhys factor $\lambda=\mu\nu R_{\text{ge}}r_{\text{zpm}}$ yields the Holstein-Hamiltonian

[TABLE]

One can bring this Hamiltonian into diagonal form $\nu b^{\dagger}b+\omega_{e}\sigma^{\dagger}\sigma$ via the polaron transformation $U^{\dagger}=({\cal{D}^{\dagger}})^{\sigma^{\dagger}\sigma}$ where the displacement is defined as ${\cal{D}}=\exp(-i\sqrt{2}\lambda p)=\exp[\lambda(b^{\dagger}-b)]$ with the dimensionless momentum quadrature $p=i(b^{\dagger}-b)/\sqrt{2}$ as generator (we also define the dimensionless position quadrature $q=(b^{\dagger}+b)/\sqrt{2}$). Notice that the polaron transformation also removes the vibronic shift of the excited state.

Appendix B Quantum Langevin equations

We will list the exact and approximated quantum Langevin equations for a single vibration with the observation that the many vibrations case is a straightforward generalization. The derivation is based on the property that an exact transformation exists linking any master equation in Lindblad form to an equivalent equation of motion that properly includes noise terms associated with the Lindblad collapse operators Gardiner and Zoller (2004). Explicitly, for any operator $\cal{O}$

[TABLE]

where an input noise occurs for any collapse operator $c$ at rate $\gamma_{c}$. In our case we have collapse operators $\sigma$ and $b$ with corresponding rates $\gamma$ and $\Gamma$. After working out the commutation relations, one obtains:

[TABLE]

Under the assumption of weak driving $\eta_{\ell}\ll\gamma$ the transition dipole scales at most as $\eta_{\ell}/\gamma$ while the population scales as $\eta_{\ell}^{2}/\gamma^{2}$. The back-action term leading to a change in the vibrational dynamics is then small compared to $\Gamma$ or $\nu$ which allows us to work in perturbation theory up to second order in $\eta_{\ell}$. We also can neglect terms $\sigma^{\dagger}\sigma$ with respect to unity. This leads to the following perturbative Langevin equations (where now we transform into a frame rotating at the laser frequency) which will be the basis of our analytical treatment:

[TABLE]

Here, in (20b) we have neglected the term $\lambda\Gamma(b^{\dagger}-b)\tilde{\sigma}$ as well as the corresponding noise term $-\sqrt{2\Gamma}\lambda(b_{\text{in}}^{\dagger}-b_{\text{in}})\tilde{\sigma}$. We can however estimate the contribution of those terms by averaging (tracing) over the vibrational ground state $\ket{0}$ (assuming fast vibrational relaxation):

[TABLE]

which has to be compared with the bare electronic contribution $-\gamma\braket{\tilde{\sigma}}=-\gamma\sigma e^{-\lambda^{2}/2}$. We can see that the terms effectively lead to an additional dephasing of the dipole operator with $\gamma+\lambda^{2}\Gamma$. This would imply disagreement with experimental observations of lifetime limited transitions for the zero-phonon line of molecules as vibrational relaxation rates are much larger than radiative rates. However, we show in (SM D) that a modelling of non-unitary dynamics based on a Brownian noise model removes this inconsistency at the level of linewidths while showing perfect agreement in spectroscopic line positions. However, for simplicity of derivations and comparison with numerical simulations we rely on the simple Lindblad form master equation in the following.

The procedure we will follow to derive molecular spectra is to first integrate the vibrational dynamics $b(t)$ and subsequently replace the time evolved displacement operator ${\cal{D}}(t)$ in the equations for $\tilde{\sigma}(t)$ and $\sigma^{\dagger}\sigma(t)$. Most of the physics is contained in the two-time correlations $\braket{{\cal{D}}(t){\cal{D}^{\dagger}}(t^{\prime})}$ which we evaluate in Appendix C.

Appendix C Dynamics of the vibrational mode. Properties of the displacement operators

We integrate here the equations of motion for the vibrational operators and compute the two-time correlation functions for displacement operators as relevant to the further computation of any electronic degree of freedom operator.

Simplified model: creation and annihilation operators.

First, let us analyze the properties of the vibrational mode operators. We will follow the dynamics on timescales defined by the slow radiative evolution with a timescale $\gamma^{-1}\gg\Gamma^{-1}$. This allows us to consider that the vibrations are practically always in steady state and operators at any time $t>0$ are given by a formal integration with initial conditions extended to $-\infty$ such that

[TABLE]

The corresponding Hermitian conjugate operator is given by

[TABLE]

and the commutator of such operators results in

[TABLE]

Spectrum of the Holstein force. Non-zero temperature bath.

Let us estimate the power spectrum of the force acting on the excited electronic orbital by defining (in analogy with Ref. Neuman and Aizpurua (2018)) an effective hermitian force $F=\lambda\nu(b+b^{\dagger})$. The power spectrum is defined as:

[TABLE]

which is none other than the spectrum of the displacement variable in equilibrium with the environment. Notice that the only non-vanishing correlations at different times are:

[TABLE]

which leads to a spectrum of the effective force acting on the excited electronic orbital of the molecule

[TABLE]

One can also include a $T>0$ bath in the correlations of the noise operators $\braket{b_{\text{in}}^{\dagger}(t)b_{\text{in}}(t^{\prime})}=\Gamma\bar{n}\delta(t-t^{\prime})$ and $\braket{b_{\text{in}}(t)b_{\text{in}}^{\dagger}(t^{\prime})}=\Gamma(\bar{n}+1)\delta(t-t^{\prime})$ where the thermal occupancy $\bar{n}=[\exp({\hbar\nu/(k_{\text{B}}\text{T})})-1]^{-1}$. This would lead to a temperature dependent spectrum

[TABLE]

Two-time correlations of displacement operators.

Now we can estimate properties of the displacement operator. The trace over the displacement operators in the Heisenberg picture is done with the initial vacuum density operator $\ket{0}\bra{0}$ for the vibrational mode (as well as for all the outside modes encompassed in the input noise). We can then deduce the two-time correlations for a given time ordering $t>t^{\prime}$

[TABLE]

This leads to the main expressions that we will be using in the following especially in the sum form expressed below:

[TABLE]

and

[TABLE]

where $s^{\lambda}_{m}=e^{-\lambda^{2}}\lambda^{2m}/m!$.

Appendix D Brownian noise model.

The dynamics of a mechanical oscillator subject to thermal noise (modeled as Brownian motion dynamics) is encompassed in the following equations of motion for the dimensionless quadratures $q=(b^{\dagger}+b)/\sqrt{2}$ and $p=i(b^{\dagger}-b)/\sqrt{2}$

[TABLE]

combined with the correlations of the input noise at a given bath temperature in the time domain Genes et al. (2008)

[TABLE]

For an arbitrary temperature, the average thermal occupancy is $\bar{n}=1/(\exp(\hbar\omega/k_{\text{B}}T)-1)$. The spectrum at the $\pm\nu$ then gives the cooling and heating contributions: $S_{\text{th}}(\nu)=4\Gamma(\bar{n}+1)$ and $S_{\text{th}}(-\nu)=4\Gamma\bar{n}$. Generally, from these properties one can always estimate the damping rate $2\Gamma=(S_{\text{th}}(\nu)-S_{\text{th}}(-\nu))/2$ and equilibrium occupancy as $\bar{n}=S_{\text{th}}(-\nu)/[S_{\text{th}}(\nu)+S_{\text{th}}(-\nu)]$. Close to zero temperature $T=0$, the thermal spectrum becomes very asymmetric $S_{\text{th}}(\omega)=\left[4\omega\Gamma/\nu\right]\theta(\omega)$ as it vanishes for negative frequencies (we introduced the Heaviside function $\theta(\omega)$). Generally, one can proceed with a Fourier domain analysis of the steady state of the equations above which lead to:

[TABLE]

The zero point energy can be estimated by using $\int_{0}^{\infty}d\omega|\epsilon(\omega)|^{2}=\pi/(4\Gamma)$ which leads to $\braket{q^{2}}=1/2$ as expected. However, for the momentum the integral $\int_{0}^{\infty}d\omega|\omega\epsilon(\omega)|^{2}$ diverges and requires the introduction of a cutoff frequency $\Lambda>\nu$ to keep it finite. Under this assumption and noticing that the susceptibility is sharply peaked at $\nu$ one can evaluate the integral by fixing $\omega=\nu$ in the numerator and obtaining the expected equipartition rule. Furthermore, one can obtain the following two properties

[TABLE]

which is identical to the Langevin case (when replacing $p=i(b^{\dagger}-b)/\sqrt{2}$, $q=(b^{\dagger}+b)/\sqrt{2}$). This also implies that the correlation function for the displacement operator $\braket{{\cal{D}}(t){\cal{D}^{\dagger}}(t^{\prime})}$ is identical in both models. Calculating the equation of motion for $\dot{\tilde{\sigma}}=\dot{\sigma}\mathcal{D}^{\dagger}+\sigma\dot{\mathcal{D}}^{\dagger}$ one finds

[TABLE]

We can again estimate the contribution of the off-diagonal terms which gave rise to a dephasing term in the Langevin model. Here we obtain

[TABLE]

where we used that $\braket{p\mathcal{D}^{\dagger}}=\frac{1}{i\sqrt{2}}\frac{d}{d\lambda}\braket{\mathcal{D}^{\dagger}}$ and $\braket{\zeta\mathcal{D}^{\dagger}}=\sum_{m=0}^{\infty}\frac{(i\sqrt{2}\lambda)^{m}}{m!}\braket{\zeta(t)p(t)^{m}}=\sqrt{2}i\Gamma\lambda e^{-\lambda^{2}/2}$ (obtained by an application of the Isserlis’ theorem and with relation $\braket{\zeta(t)p(t)}=\Gamma$).

Appendix E Absorption profile (steady state analysis)

We will analyze the absorption profile of the system under continuous (weak) driving in steady state. This amounts to computing the expression of the excited state population $\braket{\sigma^{\dagger}\sigma}$ as a function of the incoming laser frequency.

Single vibrational mode.

From Eqs. (20) we find:

[TABLE]

To proceed we will need the expression of the Pauli operator in steady state which we obtain as:

[TABLE]

To find the average molecular dipole one then performs a trace over the right hand side of the operator expression:

[TABLE]

Now we can first see that the last term is vanishing while the general expression can be evaluated as

[TABLE]

Many vibrational modes.

Extending to many modes we have the same equation as before

[TABLE]

except that now the displacement is produced by all vibrations

[TABLE]

For the calculation of absorption and emission profiles we have to evaluate terms like

[TABLE]

The expectation value can be factorized as modes are independent from each other at any given time and we use for each mode

[TABLE]

so we have to compute:

[TABLE]

We can now perform the integral to find

[TABLE]

where $\Gamma_{\{m_{k}\}}=\sum_{k=1}^{n}m_{k}\Gamma_{k}$, $\nu_{\{m_{k}\}}=\sum_{k=1}^{n}m_{k}\nu_{k}$ and $\left\{m_{k}\right\}\equiv{m_{1},...m_{n}}$ is a list of all indices. This leads to the corresponding population expression

[TABLE]

Appendix F Emission spectrum

Transient regime.

A simple way to obtain the emission spectrum Auffèves et al. (2008) is to consider the transient dynamics of a molecule initially partially excited with $\braket{\sigma^{\dagger}(0)\sigma(0)}=p_{0}$. From Eq. (6) we can compute two-time correlations and eliminate all the vanishing ones (between transient and steady state contributions) to arrive at

[TABLE]

Performing the Fourier transform gives

[TABLE]

Steady state regime.

In steady state the calculation of the emission spectrum is more cumbersome as it involves the evaluation of four-point correlation functions. Using the solutions

[TABLE]

in the steady state regime the autocorrelation function is given by

[TABLE]

The four-point correlation function of the displacement operators can be easily evaluated to be

[TABLE]

Here we have used the conditions $t\geq\tau$ and $t^{\prime}\geq\tau^{\prime}$. The final derivations involve time ordering. For an ordering given by $(1)$ $t>\tau>t^{\prime}>\tau^{\prime}$ we obtain the expression

[TABLE]

Similarly, we obtain

[TABLE]

for $(2)$ $t>t^{\prime}>\tau>\tau^{\prime}$ and

[TABLE]

for $(3)$ $t>t^{\prime}>\tau^{\prime}>\tau$. Injecting these sequences into the autocorrelation function we obtain the general result

[TABLE]

For $t-t^{\prime}=0$ we obtain the result in Eq. 56 for the evolution of the population $\braket{\sigma^{\dagger}(t)\sigma(t)}$. For $t-t^{\prime}\gg 1/\Gamma$ the saturation value

[TABLE]

is reached.

To calculate the emission spectrum we need to evaluate the Fourier transform of the autocorrelation function

[TABLE]

where

[TABLE]

Appendix G Molecule-cavity spectroscopy

Langevin equations.

The master equation to Langevin equation procedure is now applied to the dynamics including the Jaynes-Cummings interaction and decay at rate $\kappa$ of the cavity mode. Notice that the natural basis is that of polaron operators both for the field mode $\tilde{a}=a{\cal{D}}^{\dagger}$ and molecule $\tilde{\sigma}=\sigma{\cal{D}}^{\dagger}$:

[TABLE]

With the simplifications assumed before (weak excitation) the following set of solvable equations emerges:

[TABLE]

Cavity transmission.

While the first equation is independent, the last two are coupled and we write them in matrix form:

[TABLE]

where the drift matrix is

[TABLE]

accompanied by a transformation into a frame rotating with the laser frequency. The input vector is

[TABLE]

Diagonalization of the M matrix implies that $\textbf{M}=\textbf{T}\bm{\Lambda}_{M}\textbf{T}^{-1}$ where $\bm{\Lambda}_{M}$ is a diagonal matrix with $\epsilon_{\pm}$ on the diagonal. The T matrix is made of the eigenvector components. We can now write

[TABLE]

and solve for steady state

[TABLE]

One can write the field component (now transforming back from the polaron operators to the bare basis) in a compact form

[TABLE]

If we want to estimate expectation values, remember that everything is zero except for the driving contribution in the input terms so that (generalizing to many vibrational modes)

[TABLE]

Decay of the upper polariton into lower polariton.

In the single excitation subspace the Jaynes-Cummings part of the Hamiltonian can be diagonalized by defining annihilation operators for upper and lower polariton (assuming resonance $\omega_{c}=\omega_{e}$):

[TABLE]

We start with the equations of motion for the bare polaritons

[TABLE]

as well as the equation of motion for the population of the lower polariton

[TABLE]

from which we can see that the term $\Im\braket{L^{\dagger}U\left(b^{\dagger}+b\right)}$ is the one responsible for population transfer between the polaritons. Making the notation $\epsilon_{+}=(\gamma_{+}+i\omega_{+})$ and $\epsilon_{-}=(\gamma_{-}+i\omega_{-})$ we will analyze the decay dynamics of the upper polariton assumming initial conditions $U(0)$ and $L(0)$. The first order solution for upper and lower polaritons are (neglecting the noise terms as they do not contribute to the process):

[TABLE]

with this we can calculate $\braket{L^{\dagger}U(b^{\dagger}+b)}$ (using that $\braket{b(t^{\prime})b^{\dagger}(t)}=e^{-\left(\Gamma-i\nu\right)(t-t^{\prime})}$):

[TABLE]

Assuming $\Gamma\gg\gamma_{+},\gamma_{-}$ and $\gamma_{+}=\gamma_{-}$ we can write the population transfer as a rate proportional to the population of the upper polariton:

[TABLE]

Appendix H Energy transfer without cavity

We can write the Langevin equations for dipole-dipole coupled donor/acceptor molecules in the transient regime, assuming some initial population of donor $p_{D}(0)$ and zero initial population of the acceptor:

[TABLE]

Notice that we have chosen $\Delta=\omega_{D}-\omega_{A}>0$. We will perform a perturbative treatment where the evolution of $\tilde{\sigma}_{D}(t)$ is in a first approximation given by

[TABLE]

Let us write the hermitian conjugate as well on the two components

[TABLE]

The evolution of $\sigma_{A}(t)$ is then similarly obtained as:

[TABLE]

We will also split the acceptor operator in a two parts $\sigma_{A}(t)=\sigma^{t}_{A}(t)+\sigma^{s}_{A}(t)$. Again the hermitian conjugate on components

[TABLE]

From the Eq. (20) we see that

[TABLE]

The last term containing correlations between acceptor input noise and the evolved acceptor operator does not contribute. The energy transfer is mediated by the $\Im{\braket{\sigma^{\dagger}_{A}\sigma_{D}}}$ term. We will need to evaluate four terms:

[TABLE]

We can now write

[TABLE]

and estimate the integral as:

[TABLE]

Due to the rapid decay of terms harboring $\Gamma_{A}$, $\Gamma_{D}$ we can rewrite this term and obtain the expression

[TABLE]

Ignoring the $n_{1}=n_{4}=0$ components, the real part of the integral is given by

[TABLE]

and therefore we obtain

[TABLE]

Appendix I Energy transfer with cavity

The Langevin equations for a pair of cavity and dipole-dipole coupled donor/acceptor molecules in the transient regime is given by

[TABLE]

where we assume that we have some initial population $P_{D}(0)$ of the donor and zero initial population in the cavity and of the acceptor. Again we define $\Delta=\omega_{D}-\omega_{A}>0$ and $\Delta_{c}=\omega_{D}-\omega_{c}$. Following this perturbative ansatz we obtain

[TABLE]

Using these result in the differential equation of the acceptor population

[TABLE]

we are left with deriving the terms $\Im{\braket{\sigma^{\dagger}_{A}\sigma_{D}}}$, $\Im{\braket{\sigma^{\dagger}_{A}a}}$. The dipole-dipole term equates to

[TABLE]

while the cavity term is given by

[TABLE]

where the term with $a^{\dagger}(0)a(0)$ goes to zero. Evaluating the integrals for the dipole-dipole expression results in

[TABLE]

which simplifies to

[TABLE]

if we take into consideration that terms that contain $\Gamma_{D}$ or $\Gamma_{A}$ decay rapidly. By focussing on the terms that decay with $e^{-2\gamma_{D}t}$ and ignoring all oscillatory terms we finally obtain

[TABLE]

Evaluating the integrals for the cavity induced term results in

[TABLE]

with

[TABLE]

Again by using the fact that $\Gamma_{D}$ or $\Gamma_{A}$ decay rapidly, we obtain

[TABLE]

By focussing on the terms that decay with $e^{-2\gamma_{D}t}$ or $e^{-2\kappa t}$ and ignoring all oscillatory terms we obtain

[TABLE]