Theory for the stationary polariton response in the presence of vibrations

TL;DR

This paper develops a quantum optical model for surface plasmon polaritons coupled with molecules, incorporating molecular vibrations to explain asymmetric emission and decoherence in hybrid systems, relevant for applications like coherent light harvesting.

Contribution

It introduces a novel theoretical framework that includes molecular vibrations via $P(E)$ theory, enhancing understanding of polariton response and emission asymmetries in hybrid systems.

Findings

Accounts for the effect of molecular vibrations on polariton response

Explains asymmetric emission from upper and lower polariton modes

Provides a link between hybrid mode response and molecular fluorescence

Abstract

We construct a model describing the response of a hybrid system where the electromagnetic field - in particular, surface plasmon polaritons - couples strongly with electronic excitations of atoms or molecules. Our approach is based on the input-output theory of quantum optics, and in particular it takes into account the thermal and quantum vibrations of the molecules. The latter is described within the $P(E)$ theory analogous to that used in the theory of dynamical Coulomb blockade. As a result, we are able to include the effect of the molecular Stokes shift on the strongly coupled response of the system. Our model then accounts for the asymmetric emission from upper and lower polariton modes. It also allows for an accurate description of the partial decoherence of the light emission from the strongly coupled system. Our results can be readily used to connect the response of the hybrid modes to the emission and fluorescence properties of the individual molecules, and thus are relevant in understanding any utilization of such systems, like coherent light harvesting.

Equations

$$H_{s-c}={\omega }_c{c}^{†}c+\sum _{j=1}^N\text{quantity}({\omega }_m{\sigma }_j^{†}{\sigma }_j+g_j{\sigma }_j^{†}c+g_j^{\ast }{c}^{†}{\sigma }_j).$$

$$H_v=\sum _{j=1}^N{\omega }_v{b}_j^{†}{b}_j+\sum _{j=1}^N{\omega }_v\sqrt{S}{\sigma }_j^{†}{\sigma }_j\text{quantity}(b_j^{†}+b_j).$$

$$\dot{c}$$

$${\dot{\sigma }}_j^S$$

$${\Sigma }_{\mathrm{out},s/p}^{T/R}=\text{quantity}({\delta }_R^{T/R}{c}_{\mathrm{in}}+\sqrt{{\kappa }_o^{T/R}}c){\delta }_p^{s/p}+\sum _j{\eta }_{j,s/p}^{T/R}{\sigma }_j.$$

$$S_{s/p}^{T/R}(\omega ;{\omega }_d)=\frac{1}{2\pi }\int \text{differential}t{e}^{i\omega t}\text{expectationvalue}{\Sigma }_{\mathrm{out},s/p}^{T/R†}(0){\Sigma }_{\mathrm{out},s/p}^{T/R}(t),$$

$$P(E)=\frac{1}{2\pi }\int \text{differential}t{e}^{iEt}P(t).$$

$$P(t)=e^{\mathcal{T}⟨(i\sqrt{S}[p(t)-p(0)]{)}^2⟩/2}=e^{S⟨[p(t)-p(0)]p(0)⟩},$$

$$\begin{array}{rl}\dot{x}(t)& ={\omega }_{v}p(t)\\ \dot{p}(t)& =-{\omega }_{v}x(t)-\gamma p(t)+\xi (t),\end{array}$$

$$⟨\xi (t)\xi (t^{\prime })⟩=\int \text{differential}\omega \exp [-i\omega (t-t^{\prime })]S_{\xi }(\omega ),$$

$$S_{\xi }(\omega )=\frac{\gamma \omega }{\pi {\omega }_v}\left[\coth \left(\frac{\omega }{2k_{B}T}\right)+1\right]$$

$$\left(\begin{array}{c}x(\omega )\\ p(\omega )\end{array}\right)=\frac{1}{{\omega }^2-{\omega }_v^2+i\omega \gamma }\left(\begin{array}{c}-{\omega }_v\\ i\omega \end{array}\right)\xi (\omega ).$$

$$J(t)=S\text{expectationvalue}p(t)p(0)=\frac{S\gamma }{\pi {\omega }_v}\int \text{differential}\omega e^{-i\omega t}\frac{{\omega }^3}{({\omega }^2-{\omega }_v^2{)}^2+{\omega }^2{\gamma }^2}\left[\coth \left(\frac{\omega }{2k_{B}T}\right)+1\right].$$

$$J(t)=S(n_{\mathrm{th}}+1)e^{-i{\omega }_{v}t-\frac{\gamma }{2}\text{absolutevalue}t}+S{n}_{\mathrm{th}}{e}^{i{\omega }_{v}t-\frac{\gamma }{2}\text{absolutevalue}t}.$$

$$\mathbb{E}(E^n)=i^n\text{evaluated}{\text{derivative}[n]P(t)t}_{t=0}=i^n\text{evaluated}{\text{derivative}[n]e^{J(t)}t}_{t=0}.$$

$$P(-E)=\exp \left(-\frac{E}{k_{B}T}\right)P(E).$$

$$L(t_1,t_2,t_3)=\text{expectationvalue}{Q}^{†}(t_1)Q(t_2)Q^{†}(0)Q(t_3).$$

$$L(t_1,t_2,t_3)=e^{-S\mathcal{T}⟨[p(t_1)-p(t_2)+p(0)-p(t_3){]}^2⟩/2},$$

$$L(t_1,t_2,t_3)=\frac{P(t_1-t_2)P(t_1-t_3)P(t_2)P(-t_3)}{P(t_1)P(t_2-t_3)}.$$

$$\tilde{P}(\omega )=\frac{1}{\pi }\frac{\gamma {\omega }^2}{({\omega }_v^2-{\omega }^2{)}^2+{\omega }^2{\gamma }^2}$$

$$J(t)=S\text{quantity}(n_{\mathrm{th}}+1)e^{-i{\omega }_{v}t}+S{n}_{\mathrm{th}}{e}^{i{\omega }_{v}t}.$$

$$P(E)$$

$$p_k(S)$$

$$\mathbb{E}(E)=S{\omega }_v\text{mbox}and\text{var}(E)=(2n_{\mathrm{th}}+1)S{\omega }_v^2.$$

$$1/P(t)$$

$$L({\omega }_1,{\omega }_2,{\omega }_3)$$

$$=\sum _{k_1,k_2,k_3,k_4,k_5,k_6}(-1{)}^{k_1+k_2}{p}_{k_1}{p}_{k_2}{p}_{k_3}{p}_{k_4}{p}_{k_5}{p}_{k_6}{e}^{4S(2n_{\mathrm{th}}+1)}\delta ({\omega }_1-\text{quantity}[k_1+k_3+k_4]{\omega }_v)$$

$$\times \delta ({\omega }_2-\text{quantity}[k_2-k_3+k_5]{\omega }_v)\delta ({\omega }_3+\text{quantity}[k_2+k_4+k_6]{\omega }_v).$$

$$\text{expectationvalue}{Q}_j^{†}({\omega }_1)Q_j({\omega }_2)=P({\omega }_1)\delta ({\omega }_1+{\omega }_2)$$

$$\text{expectationvalue}{Q}_j^{†}({\omega }_1)Q_j({\omega }_2)Q_j^{†}({\omega }_3)Q_j({\omega }_4)$$

$$=L({\omega }_1,{\omega }_2,{\omega }_4)\delta ({\omega }_1+{\omega }_2+{\omega }_3+{\omega }_4),$$

$$A(\Delta )=\int \text{differential}EP(E)\chi (\Delta -E)$$

$$F=\int \text{differential}{\omega }_1\text{differential}{\omega }_{2}L({\omega }_1,{\omega }_d-\omega -{\omega }_1,{\omega }_2)\chi ({\omega }_1-\Delta )\chi ({\omega }_2+\Delta ).$$

Full text

I Introduction

Photonic structures, such as optical cavities or surface plasmon polaritons can modify electromagnetic vacuum field by confining the light to smaller volumes and restricting the number of available photonic modes. Any electronic excitation inside such modified vacuum can interact much stronger with the confined light mode. This interaction can become strong enough for the coupling energy to show up in absorption and emission spectra of such systems, suggesting formation of hybrid light-matter states, called polaritons. Common examples studied in this strong coupling limit are single atomsRempe et al. (1987), excitons in semiconductorsWeisbuch et al. (1992), and photoactive molecules Törmä and Barnes (2014); Yu et al. (2019).

More recently strong coupling of molecules with confined light modes has been in the focus of interest, because the hybridization between light and matter into polaritons not only delocalizes the excitation over many molecules but also changes their potential energy surface, and thus provides a new way to control chemistryHertzog et al. (2019). Experiments on strongly coupled molecules have already shown (i) suppression of photo-oxidation of TBDC J-aggregates coupled to plasmonic nano-prismsMunkhbat et al. (2018), and of photo-isomerization of Spiropyran inside an optical cavity Hutchison et al. (2012); (ii) enhanced electronic conductivity in organic semiconductorsOrgiu et al. (2015); (iii) inter-molecular excitation energy transfer over large distances inside optical cavitiesZhong et al. (2016, 2017); and (iv) enhanced decay of triplet states in Erythrosine B moleculesStranius et al. (2018). Since polaritons are like interacting dressed photons with mass they can undergo Bose-Einstein condensation even at room temperatureDeng et al. (2010), which further enables very efficient and thresholdless polariton lasingDeng et al. (2003); Christopoulos et al. (2007).

Strong coupling between a single molecule and electromagnetic field is very hard to achieve Chikkaraddy et al. (2016); Wang et al. (2019). The common way to circumvent this problem is to couple multiple molecules to the same photonic mode. Often these systems are still described within effective two-state models accounting only the two polaritonic states Schwartz et al. (2013). However, such a description disregards the fact that the visible polariton modes are now superpositions of several molecular excitations and the photonic mode, and they are not the only eigenmodes of the system. The response of the whole system also depends on the presence of ”dark modes”, i.e. superpositions having no photonic component. These dark modes become relevant especially when dissipation processes within the molecules, such as those linked to vibrations, are included. In this case they can dramatically affect the predicted guiding of the chemistry, and even the validity of the whole concept. They have been taken into account in some multiscale simulations coupling the investigated molecules to thermal environmentsGroenhof et al. (2019). However, those simulations often consider the transient response, whereas majority of experiments on light-matter coupling concern stationarily driven setups.

Here we construct a detailed description of the stationarily driven response of the strongly coupled system, taking into account the effect of inhomogeneous broadening of the molecular response due to quantum and thermal vibrations of the molecules. We take the vibrations into account via the $P(E)$ theory analogous to that used in Coulomb blockade Ingold and Nazarov (1992); Heikkilä (2013). This theory describes the probability of absorbing (for $E>0$) or emitting ($E<0$) the energy $E$ to/from the vibrations. For the specific models of harmonic vibrations, such a $P(E)$ function can be calculated exactly. In general, we find how this $P(E)$ is related to the absorption and emission spectra of individual molecules. Therefore, an alternative approach is to deduce an effective $P(E)$ for the measured spectra of individual molecules. The resulting fluorescence spectrum is similar to that found via quantum many-body theory Mahan (2013). However, in these approaches only the transient response is considered which requires assumptions on the initial state of the system. In our work, the focus is on the stationary response, in which case different conservation laws and their explicit breaking come into the focus. What is more, we connect this fluorescence spectrum directly to the absorption/emission spectrum of the strongly coupled system. In a certain limit of parameters, the resulting inhomogeneous broadening of the molecular absorption/emission then determines the linewidth of the polariton modes. In particular, our model explains the asymmetric emission spectra of upper and lower polaritons seen in many experiments Baieva et al. (2017); Bellessa et al. (2004); Hakala et al. (2009); Baieva et al. (2012); Koponen et al. (2013); Guebrou et al. (2012), as well as the varying polarization of that emission depending on the quantum coherence of the system as shown recently Baieva et al. (2017).

Besides the detailed description of the vibrations we include polarization of the confined light field and the positions of the molecules. The position of a molecule with respect to the light field only amounts to a phase factor to the light-matter coupling but collectively it leads to experimentally observable effects. Perhaps the most striking effect is the superradiance due to coherent emission described by Dicke in the 1950’sDicke (1954), but conversely it is seen in the usual experiments where many molecules are distributed over a region larger than the light wavelength. On the other hand, the polarization of the light and the transition dipole moment of a molecule determine whether there is coupling at all: if these directions are perpendicular, the coupling vanishes. This provides another way to control light-matter interaction which could be used in applications Lodahl et al. (2017). In this article we describe both the incoherent and coherent limits of polaritonics.

Although many aspects of our theory can be generalized to any confined light mode, like resonances of Fabry-Pérot cavities, here we focus in particular on surface plasmon polaritons (below, plasmons) driven by an external light field. Plasmons are evanescent like electromagnetic modes propagating along metal-dielectric surface with a 2D momentum ${\vec{k}}$ along the surface. In general they have a non-linear dispersion $\omega({\vec{k}})$, and due to the evanescent nature their electromagnetic field is highly confined to the surface. Because of this confinement, the dipolar coupling to molecular excitations residing at the surface can be made strong Bellessa et al. (2004); Törmä and Barnes (2014) leading to the observed avoided crossing between the two systems and thus offering possibility to control photochemical reactions. A typical way to launch plasmons is via Kretschmann configuration, i.e. coupling an external electromagnetic field to the surface modes via a prism Kretschmann and Raether (1968). In this setup the angle with which the light enters the prism determines a specific plasmon ${\vec{k}}$ vector. Hence, in this work we concentrate on a single plasmon mode with defined ${\vec{k}}$ and a generic frequency $\omega_{c}$.

To be specific, we consider the plasmon–molecule system in the strong-coupling regime. We describe the plasmon by a single bosonic mode $c$ of frequency $\omega_{c}$ and a given polarization $\hat{u}_{pl}$ with respect to its wavevector $\vec{k}$. A concrete example of such a plasmon is the surface plasmon polariton traveling along an interface in the $xy$-plane in the $y$-direction with $\hat{u}_{pl}=\quantity(0,\sin{\beta},\cos{\beta})$ as in Fig. 1a. The plasmon interacts with $N$ identical molecules scn which we approximate as two-level systems with transition frequency $\omega_{m}$. We denote the rising (lowering) operator of a molecule with $\sigma_{j}^{\dagger}$ ($\sigma_{j}$). As in typical experiments, we assume that the electric dipole moments of the molecules point in uniformly random directions $\hat{n}_{j}$. Following the standard approach of quantum optics Walls and Milburn (2007); González-Tudela et al. (2013), the Hamiltonian of the strong-coupled system is in the rotating wave approximation ($\hbar=1$)

[TABLE]

The position $\vec{r}_{j}$ of a molecule affects the coupling $g_{j}$ in two ways: it contains a complex phase factor due to the phase of the plasmon, and the coupling strength depends on the distance to the interface. If this distance is independent of the polarization the latter effect may be disregarded and the average value used. Also, the coupling strength depends on the angle between the plasmon polarization and dipole moment of a molecule. Thus, we write $g_{j}=ge^{i\vec{k}\cdot\vec{r}_{j}}\quantity(\hat{n}_{j}\cdot\hat{u}_{pl})$.

In addition to the strong-coupled system we include the vibrational modes of the molecules. We assume a single vibration mode $b_{j}$ per molecule with eigenfrequency $\omega_{v}$ but the generalization to multiple modes is straightforward (Appendix A). These vibrations and their interactions are described by

[TABLE]

The coupling between electronic and vibrational modes is quantified with a dimensionless parameter $\sqrt{S}$, the Huang–Rhys factor Huang and Rhys (1950), which is related to the Stokes shift measured in fluorescent emission.

We seek an approach to find the response of the strongly-coupled plasmon–molecule system in the presence of vibrations. To this end, we employ the input-output formalism of quantum optics Gardiner and Collett (1985); met . We assume that there are separate bosonic baths for each molecule, vibration and the plasmon to which the coupling is linear in $\sigma_{j},b_{j}$ and $c$, respectively. In the Markov approximation these couplings are described by the dissipation rates $\tilde{\kappa}_{j},\gamma_{j}$, and $\kappa$ of the molecules, vibrations, and plasmon. In the following, we suppose identical molecules and vibrations so that $\gamma_{j}=\gamma$ and $\tilde{\kappa}_{j}=\tilde{\kappa}$. We neglect the thermal fluctuations of plasmons and molecules here as $\hbar\omega_{m},\hbar\omega_{c}\gg k_{B}T$ even at room temperature. We simplify the molecule–vibration Hamiltonian by introducing a new polaron operator $\sigma^{S}_{j}=e^{\sqrt{S}\quantity(b_{j}^{\dagger}-b_{j})}\sigma_{j}\equiv Q_{j}\sigma_{j}$. Lastly, we assume a low driving power which corresponds to the single-excitation limit $\sigma^{\dagger}\sigma\approx 0$. We find that when $\frac{\hbar\omega_{v}}{k_{B}T}>\frac{\gamma S}{2\kappa_{m}}$, where $\kappa_{m}=\tilde{\kappa}+\gamma S$ is the total effective damping rate of the individual molecules, the dynamics of the vibrational modes $b_{j}$ are approximately uncoupled from the plasmon–molecule system as shown in Appendix B. This allows us to use the Caldeira–Leggett model Caldeira and Leggett (1981) for the vibrational dynamics. The plasmon and molecular equation are in this case

[TABLE]

where $\tilde{\omega}_{m}=\omega_{m}-S\omega_{v}$ is the renormalized molecular frequency while $\kappa_{\text{ext}}$ and $\kappa_{m}^{\rm ext}$ are the couplings to external driving fields. For the plasmon–molecule system we assume that only the plasmon is driven so that $c_{\rm in}=\alpha e^{-i\omega_{d}t}$ and $\sigma_{{\rm in},j}=0$.

We model a measurement on the plasmon–molecule system so that the incoming light $c_{\rm in}$ produces a reflected $\Sigma^{R}_{\rm{out}}$ and transmitted field $\Sigma^{T}_{\rm{out}}$. These fields contain both the plasmon and the molecular emission but not the emission of phonons from the vibrations, because they are usually not measured. Phonon emission hence allows for a loss of energy in the process, so that the power in the output fields can be lower than the one in the input. We also separately include coupling to s- and p-polarized light represented by $\hat{u}_{p}=\hat{u}_{y}$ and $\hat{u}_{s}=\hat{u}_{x}$ (Fig. 1a). Since the propagating plasmon cannot emit s-polarized light to the direction perpendicular to the interface but the molecules have no directional preference, we consider s- and p-polarized output fields separately. The output fields obey a general expression

[TABLE]

In this equation $\delta_{R}^{T}=0$ and $\delta_{R}^{R}=1$ meaning that only the reflected field interferes with the input field. The $\delta_{p}^{s/p}$ is defined similarly because the plasmon couples only to p-polarized modes. The constants $\eta_{j,s/p}^{T/R}$ describe the coupling of the molecule electronic states to the environmental s- and p-polarized free space modes, and thus $\eta_{j,s/p}^{T/R}=\sqrt{\kappa_{m}^{T/R}}\quantity(\hat{n}_{j}\cdot\hat{u}_{s/p})$. These fields and couplings to the system are represented schematically in Fig. 1b. The output spectral density is obtained from

[TABLE]

where $\omega$ is the frequency of the output field and $\omega_{d}$ the driving frequency.

We note that the Markov approximation leading to Eqs.  (3) disregards the heating of the various baths of the plasmons, molecules and the vibrations. These heating effects can be disregarded when the heat conductance from those baths to other degrees of freedom exceeds that due to the losses described by $\kappa,\kappa_{m}$, and $\gamma$.

II $P(E)$ theory

The presence of the vibrations makes the input-output equations (3) non-linear as they contain products of different dynamical fields. This non-linearity leads to an inelastic (fluorescent) response of the molecules to the light field, where the emitted light from the molecules takes place at lower frequencies than the absorption. This is often referred as the Stokes shift. In order to take this non-linearity into account in the output spectra we introduce $P(E)$ theory similar to the one in dynamical Coulomb blockade Ingold and Nazarov (1992). Recently, the same problem has been discussed in Ref. Reitz et al., 2019 using similar methods, but only in a specific limit of vibrations (see below). The identification of $P(E)$ allows for a more general approach enabling also resolving of polariton emission, which is lacking from Ref. Reitz et al., 2019. Let us first define $P(t)=\expectationvalue{Q_{j}^{\dagger}(t)Q_{j}(0)}$ and its Fourier transform

[TABLE]

When the molecules are identical, $P(E)$ does not depend on the molecule index $j$; this assumption is easily lifted if needed. The $P(E)$ function normalizes to unity and is real for stationary vibrations, i.e. $\expectationvalue{Q_{j}^{\dagger}(t+\tau)Q_{j}(\tau)}=\expectationvalue{Q_{j}^{\dagger}(t)Q_{j}(0)}$ for any time $\tau$. We can thus interpret the $P(E)$ function as a probability distribution of transforming energy $E$ to the vibrations ($E>0$) or vice versa ($E<0$). This $P(E)$ function is characterized by four parameters: vibration eigenfrequency $\omega_{v}$, their linewidth $\gamma$, Huang–Rhys factor $S$, and temperature $T$ of their bath. It constitutes a full description of the response of the vibrations.

We present a general derivation of the $P(E)$ function and a related $L$ function assuming Gaussian thermal fluctuations. Then, we derive $P(E)$ analytically in the limit in which $\gamma$ vanishes. In this regime $P(E)$ is related to the absorption function defined by Huang and Rhys Huang and Rhys (1950). However, our analytic results for the response apply also in the case of general $\gamma$ and can be used for different models of vibrations.

II.1 Derivation of $P(E)$

We now derive the $P(E)$ function analytically in a similar manner as in the context of dynamical Coulomb blockade Ingold and Nazarov (1992). To establish notation, we omit the molecular index $j$ here and denote $x=b^{\dagger}+b$ and $p=i\quantity(b^{\dagger}-b)$, the dimensionless position and momentum operator, respectively. Then, we may write $Q^{\dagger}(t)=e^{i\sqrt{S}p(t)}$ in the correlator $P(t)$, which is the inverse Fourier transform of $P(E)$. This correlator can be evaluated for thermal vibrations. If the vibrations are described by a harmonic oscillator Hamiltonian, the fluctuations are Gaussian and the weak version of the Wick’s theorem (see e.g. Ref. Ingold and Nazarov, 1992 and an example of a non-Gaussian $P(E)$ in Ref. Heikkilä et al., 2004) applies. We identify $P(t)$ as the characteristic function of fluctuations of the stochastic quantity $p(t)-p(0)$, where everywhere in the calculations $p(t)$ should be ordered to the left of $p(0)$. We assume the thermal vibrations to be stationary, and therefore the expectation value of $p(t)-p(0)$ vanishes (as $\expectationvalue{p(t)}=\expectationvalue{p(0)}$ for stationary vibrations). Consequently for Gaussian fluctuations we can write the characteristic function in terms of the variance alone. In that case Heikkilä (2013),

[TABLE]

where the latter equality uses the fact that $\langle p(t)^{2}\rangle=\langle p(0)^{2}\rangle$. The operator ${\cal T}$ takes care of ordering $p(t)$ before $p(0)$, but that operator is no longer needed in the second equality because there $p(t)$ always precedes $p(0)$ in operator products.

Now, $p(t)$ can be obtained by solving the quantum Langevin equations without rotating wave approximation (also known as the Caldeira–Leggett model Caldeira and Leggett (1981))

[TABLE]

where $\gamma$ is the linewidth of vibrations, and $\xi$ is a Langevin force describing the thermal fluctuations. It has the correlator

[TABLE]

where the noise correlator is given by

[TABLE]

for thermal noise Giovannetti and Vitali (2001).

The Langevin equations (8) can be solved via Fourier transform. The result is

[TABLE]

After some Fourier analysis with the help of Eqs. (9), (10) and (11) we find $P(t)=e^{J(t)-J(0)}$ according to Eq. (7) with

[TABLE]

The resulting $P(E)$ is thus governed by three dimensionless parameters: the Huang–Rhys factor $S$, the quality factor of vibrations $\omega_{v}/\gamma$, and the relative temperature $k_{B}T/\omega_{v}$. Note that in the Caldeira–Leggett model $J(t)$ is related to the vibrational spectral density $J_{v}(t)=S\omega_{v}^{2}\expectationvalue{x(t)x(0)}$ via their respective Fourier transforms by $J_{v}(\omega)=J(\omega)/\omega^{2}$.

A simpler expression for $J(t)$ is obtained if instead of the Caldeira–Leggett model one uses the usual quantum optical equation $\dot{b}=-(i\omega_{v}+\frac{\gamma}{2})b+\sqrt{\gamma}b_{\rm in}$ as in Reitz et al. (2019) for example. Then for white noise $\expectationvalue{b_{\rm in}(t)b_{\rm in}^{\dagger}(t^{\prime})}=\quantity(n_{\rm th}+1)\delta(t-t^{\prime})$ we find

[TABLE]

Here $n_{\rm th}=\quantity(e^{\omega_{v}/(k_{B}T)}-1)^{-1}$ is the Bose factor, i.e. the mean number of thermal phonons at the vibrational frequency $\omega_{v}$. We arrive at the same solution from the Caldeira–Leggett model by using the method of residues to calculate the integral (12) and then approximating $\gamma\ll\omega_{v}$. This is hence the limit where Ref. Reitz et al., 2019 is valid. However, typical multiscale quantum chemistry calculations assume the opposite limit of a large $\gamma\gtrsim\tilde{\kappa}$, where molecular vibrations decay before the photon excitation.

Lastly, there are two general properties of the $P(E)$ function worth noting. First, since $P(t)$ may be regarded as a characteristic function of the probability distribution $P(E)$, the raw moments of the energy can be expressed as

[TABLE]

With the help of this formula the mean and variance of $P(E)$ can be found. Second, the Kubo–Martin–Schwinger (KMS) relation for thermal fluctuations at temperature $T$ leads to the detailed balance condition (or emission-absorption asymmetry) for $P(E)$

[TABLE]

This asymmetry in $P(E)$ is relevant for the anti-Stokes part of the spectrum. Some approximations, such as the white noise approximation, break this balance condition.

II.2 $L$ function

Another quantity we encounter that is relevant for the emission spectrum of a molecule is the Fourier transform of the four-point correlator

[TABLE]

This function is clearly related to $P(t)$ as for certain time arguments it coincides with the definition of $P(t)$, e.g. $L(t,0,0)=P(t)$. Using the same assumptions as in the derivation of $P(E)$ we may write

[TABLE]

where $\cal T$ orders operator products so that they are in the same order as in Eq. (16). Now, since the vibrations are stationary, we may write $L$ in terms of $P(t)$’s

[TABLE]

Even if we can fully calculate $J(t)$, the Fourier transform of $L$ is not straightforward to evaluate numerically in the general case.

II.3 $\gamma=0$ limit

Next, we consider the limit in which the dissipation rate of vibrations vanishes and derive expressions for both $P(E)$ and $L$. We note that the definition of $J(t)$, Eq. (12), contains a nascent delta function

[TABLE]

which in the limit $\gamma\rightarrow 0$ reduces to $\tilde{P}(\omega)=\frac{1}{2}\quantity(\delta(\omega-\omega_{v})+\delta(\omega+\omega_{v}))$. Therefore

[TABLE]

Note that this coincides with the limit $\gamma\rightarrow 0$ in the white noise model Eq. (13). The corresponding characteristic function $P(t)$ is known in probability theory to be that of the Skellam distribution Skellam (1946). It is a distribution that describes the difference of two independent Poisson processes. In our case, these processes are the emission and absorption of phonons. $P(E)$ then describes the total number of phonons transferred from/to vibrations to/from their environment. The resulting $P(E)$ function is

[TABLE]

where $I_{k}(x)$ is the modified Bessel function of the first kind. In the zero-temperature limit $p_{k}(S)=e^{-S}\frac{S^{k}}{k!}$ for $k\geq 0$ and $p_{k}(S)=0$ for $k<0$, i.e. the probability to emit phonons becomes Poissonian and the absorption probability vanishes.

We find the average and variance of the $\gamma=0$ distribution by using Eqs. (14) and (20)

[TABLE]

The variance depends on the temperature so that for high temperatures $k_{B}T\gg\omega_{v}$ the variance is directly proportional to the temperature; $\text{var}(E)\approx 2S\omega_{v}k_{B}T$. It should be noted that both the variance and the average are proportional to $S$ which also holds for a Poissonian quantity. The physical picture is that the mean of $E$ describes the Stokes shift in the molecules whereas the variance (or standard deviation) is connected with the inhomogeneous broadening of the molecular linewidth due to vibrations.

Finally, we derive $L$ in the $\gamma\rightarrow 0$ limit using Eq. (18). It is necessary to simplify $1/P(t)$ in order to find the Fourier transform of $L$. Since $P(t)=e^{J(t)-J(0)}$ and $J(t)\propto S$ we may find $1/P(t)$ by changing $S\rightarrow-S$ in Eq. (21). Using the parity of the modified Bessel function of the first kind $I_{k}(-x)=(-1)^{k}I_{k}(x)$ we can express the inverse as

[TABLE]

Below, we omit the $S$-dependence and denote $p_{k}(S)=p_{k}$.

The Fourier transform of $L$ is straightforward with the help of Eq. (23). We obtain

[TABLE]

This result can be used to obtain the fluorescence spectrum of a molecule. The expression is slightly cumbersome to use because the six sums obtain values from $-\infty$ (or from [math] when $T=0$) to $\infty$. This problem is alleviated by the rapid decrease of $p_{k}$ as a function of $k$. Consequently, Eq. (24) is straigthforward to compute numerically.

III Stokes shift

Before solving the full plasmon–molecule problem we illustrate how the $P(E)$ theory is used to model a measurement of the Stokes shift in a molecule–vibration system. This is achieved by removing the plasmon term from Eq. (3b) and driving the molecules, i.e., adding the term $\sigma_{{\rm in},j}=\frac{\alpha e^{i\theta_{j}}}{\sqrt{N}}\delta(\omega-\omega_{d})$ where $\theta_{j}$ represents the phase of the driving field for molecule $j$. The driving is scaled so that the total input power spectral density is given by $I_{\rm in}=\absolutevalue{\alpha}^{2}\delta(\omega-\omega_{d})$. Then, we solve Eq. (3b) with Fourier transform and convolution theorem. The spectra $S^{T/R}$ are found from Eq. (5) when the output fields are changed to $\Sigma^{T/R}_{\rm out}=\sum_{j}(\sqrt{\kappa_{m}^{T/R}}\sigma_{j}+\delta^{T/R}_{R}\sigma_{{\rm in},j})$. Here, the ’reflected’ field should not be understood literally but rather as the field that contains the driving field. The ’transmitted’ field is fully from the molecular fluorescence. Since $\sigma_{j}=Q_{j}^{\dagger}\sigma_{j}^{S}$ and the solution $\sigma_{j}^{S}$ of Eq. (3b) depends on $Q_{j}$ we encounter a four-point correlator $\expectationvalue{Q^{\dagger}_{j}(\omega_{1})Q_{j}(\omega_{2})Q_{k}^{\dagger}(\omega_{3})Q_{k}(\omega_{4})}$ in the calculation of $S^{T/R}$. Here, $Q^{\dagger}_{j}(\omega)$ refers to the Fourier transform of $Q_{j}^{\dagger}(t)$. Assuming that the vibration modes are independent and identical in different molecules, the correlator factorizes into two-point correlators when $j\neq k$. These resulting two- and four-point correlators are related to $P(E)$ by

[TABLE]

where $L(\omega_{1},\omega_{2},\omega_{4})$ is the Fourier transform of Eq. (18) in the general case.

When discussing the response of molecules to driving it is useful to introduce a frequency $\Delta=\omega_{d}-\tilde{\omega}_{m}$ which is the detuning between the driving and renormalized molecular frequency. Without vibrations the molecular response is characterized by $\chi(\Delta)=(i\Delta-\frac{\kappa_{m}}{2})^{-1}$ which describes Lorentzian absorption and emission spectra. However, in the presence of vibrations, the information about the spectral properties is contained in

[TABLE]

and

[TABLE]

These functions are associated with absorption and fluorescence of molecules, respectively, and play an important role in the plasmon–molecule problem.

III.1 Incoherent limit

Let us assume that molecules are randomly arranged, so that the phase $\theta_{j}$ is random. Averaging over them, the resulting spectra are

[TABLE]

The emission spectrum $S^{T}$ is determined by $F$ which describes inelastic scattering (output field frequency $\omega$ different from driving frequency $\omega_{d}$). In the ’reflected’ field $S^{R}$ we find also the input power spectral density $I_{\rm in}$ and a term proportional to $A$ representing absorption.

Both $F$ and $A$ are straightforward to determine from Eqs. (21) and (24), i.e., when the vibrational linewidth $\gamma\rightarrow 0$. Then, $F$ is also delta-peaked at frequencies $\omega-\omega_{d}=m\omega_{v}$ with an integer $m$. The absorption spectrum is obtained from power conservation $S^{A}(\omega_{d})=I_{\rm in}-S^{T}-S^{R}$ evaluated at the driving frequency $\omega=\omega_{d}$ and it is mostly determined by $A$. In Fig. 2b we have plotted the emission spectrum $S^{T}$ along with the absorption spectrum $S^{A}$. The absorption maximum is at the bare molecular frequency $\omega_{m}$ while the emission maximum is at approximately $\omega_{m}-2S\omega_{v}$. The difference is the Stokes shift. The spectra correspond to the results describing the transient response obtained with Green functions Mahan (2013). However, in our stationary model the absorption is not a mirror image of the emission because the emission may happen also from the excited vibrational states.

III.2 Coherent limit

Besides the experimentally more typical incoherent situation we look at the coherent limit. Then the phase $e^{i\theta_{j}}$ is fixed. This happens for instance when the distance between the molecules is much smaller than the wavelength of the driving field or the molecules are in a suitably chosen lattice. We renormalize the input in this case to be $\sigma_{{\rm in},j}=\frac{\alpha}{N}e^{-i\omega_{d}t}$ so that again the total input power is distributed evenly and is independent of the number of molecules $N$. Then, the calculation can be repeated to give

[TABLE]

Interestingly, we obtain an explicit dependence on the number $N$ of molecules for two terms. One of those terms is the inelastic emission term $F$ which means that for large $N$ the spectra are mostly elastic. However, if the vibrations are absent, i.e. $S=0$ which leads to $P(E)=\delta(E)$, $F=\absolutevalue{A(\Delta)}^{2}\delta(\omega-\omega_{d})$ and the $1/N$ dependent terms cancel. Therefore, this coherent effect is not related to sub- or superradiance of molecules described by Dicke Dicke (1954). Rather, it is related to vibrations and their enhanced non-radiative emission which shows up as a diminishing fluorescence as the number of molecules increases.

IV Plasmon–molecule system

With the tools developed in Sec. II and III we can return to the problem of strongly coupled plasmon–molecule system and find the polarized spectra of the system using Eqs. (3) (with only the plasmon being driven, i.e. $\sigma_{{\rm in},j}=0$ in this case). We integrate Eq. (3b) from an initial time $t_{i}\rightarrow-\infty$ to $t_{f}=t$ and neglect the initial condition $\sigma_{j}^{S}(t_{i})$ which has no role in a stationary situation. We substitute this into Eq. (3a) which leads to

[TABLE]

At this point we average the equation over the fluctuating vibrations and use a mean-field approximation. This leads to the $P(E)$ function since $\expectationvalue{Q_{j}^{\dagger}(t)Q_{j}(t^{\prime})}c(t^{\prime})=P(t-t^{\prime})c(t^{\prime})$. Consequently, the elastic response of the plasmon is given by $c(t)=\alpha r(\omega_{d})e^{-i\omega_{d}t}$ where

[TABLE]

Vibrations provide a channel of relaxation broadening the response which is associated with the real part of $A$. The imaginary part contains information about the frequencies of the polariton modes. When the vibrations are absent, $S=0$ and $P(E)=\delta(E)$, the usual strong-coupling response is obtained as $A\rightarrow\chi$ with Rabi splitting proportional to $\sqrt{\sum_{j}\absolutevalue{g_{j}}^{2}}$ at $\omega_{c}=\omega_{m}$. When the vibrations are present, especially the upper polariton branch is perturbed as Fig. 3 shows.

Finally, $\sigma_{j}^{S}$ can be solved from Eq. (3b) in terms of $Q_{j}$ by Fourier transformation using the convolution theorem and $c(\omega)=\alpha r(\omega_{d})\delta(\omega-\omega_{d})$. Then we have all we need to evaluate the output spectra with Eqs. (4) and (5).

IV.1 Incoherent polaritonic response

Let us consider a large number $N$ of identical molecules with random dipole moment directions $\hat{n}_{j}$. In this case we can replace the sums over the molecule index with an integral over a surface of a sphere $\sum_{j}\rightarrow\frac{N}{4\pi}\int\differential{\Omega}$. Then, because $g_{j}=ge^{i\vec{k}\cdot\vec{r}_{j}}\quantity(\hat{n}_{j}\cdot\hat{u}_{pl})$ where $\hat{u}_{pl}=\quantity(0,\sin{\beta},\cos{\beta})$ is the plasmon polarization vector, the square of the Rabi splitting in Eq. (31) is $\sum_{j}\absolutevalue{g_{j}}^{2}=Ng^{2}/3\equiv g_{N}^{2}$. We assume that the positions of the $N$ molecules are random over a region large compared to the wavelength of the plasmon so that we may replace $e^{i\vec{k}\cdot(\vec{r}_{j}-\vec{r}_{k})}\rightarrow\delta_{jk}$ for an ensemble average. Using these assumptions the polarization dependence shows up in the spectra as the coefficients

[TABLE]

where the upper/lower line is for s/p. Above, only the terms where $j=k$ contribute in the sum which results in four-point correlators as in Eq. (25b).

The s- and p-polarized emission spectra are

[TABLE]

The main difference between the s- and p-polarized spectra is because the plasmon emits only p-polarized light. The molecular fluorescence is also slightly enhanced in this polarization for $\beta\neq 0$.

Both the s- and p-polarized emission spectra are now represented with $F$ and $A$ found in molecular fluorescence Eq. (28). Therefore, the emission of strongly coupled plasmon–molecule system is related to the properties of the plasmon and the molecules separately with a few parameters describing the plasmon–molecule coupling strength and their intrinsic decay rates. Note also that these results hold for any $P(E)$, i.e., these results are independent of the vibrational model.

Since only the terms with $F$ contribute to the inelastic emission, for a given driving frequency $\omega_{d}$ the ratio of p- and s-polarized emission at $\omega\neq\omega_{d}$ is $S^{T}_{p}/S^{T}_{s}=2-\cos(2\beta)$. For example, for an interface of vacuum and silver in the Drude model Yang et al. (2015); González-Tudela et al. (2013) $\beta\approx\frac{\pi}{12}$ and $S^{T}_{p}/S^{T}_{s}\approx 1.13$ at $\omega_{c}=2\,\rm{eV}$. This ratio is otherwise independent of the system.

In Fig. 4 we have plotted the elastic emission spectra from Eqs. (33) using the $P(E)$ function (21). The s- and p-polarized emission are similar for small $S$ but for larger values the competition between the plasmon and molecule emission becomes more noticeable. The ratio between elastically emitted p- and s-polarized power is controlled by the ratio between $\kappa_{o}^{T}$ and $C_{s/p}$ and the detuning $\omega_{c}-\omega_{m}$.

We find that the upper and lower polariton modes emit asymmetrically as in other approaches Neuman and Aizpurua (2018); del Pino et al. (2018); Herrera and Spano (2018, 2017a, 2017b); Zeb et al. (2017); Reitz et al. (2019); Michetti and La Rocca (2008) and experiments Baieva et al. (2017); Bellessa et al. (2004); Hakala et al. (2009); Baieva et al. (2012); Chovan et al. (2008). This is caused by the asymmetry of the effective dissipation rate $\Gamma=\kappa/2-g^{2}_{N}\Re{A(\Delta)}$ which is related to the molecule’s absorption spectrum. The number of molecules affects the dissipation rate only via the size of the Rabi splitting. Analytical insight can be obtained when $k_{B}T\ll\hbar\omega_{v}$ and $S\ll 1$. Then, we may consider only single-phonon processes. When $g_{N}>\omega_{v}$, we find the polariton frequencies from the response function (31) for zero detuning $\omega_{c}=\omega_{m}$ to be approximately

[TABLE]

Thus, the vibrations affect both the position of the polariton peaks as well as the size of the Rabi splitting. At these frequencies the dissipation rate in the first order of $\kappa_{m}$ is given by

[TABLE]

Due to vibrations the dissipation rate of the upper polariton ($\Gamma_{+}$) is larger than the dissipation rate of the lower polariton ($\Gamma_{-}$) which suppresses the upper polariton emission compared to the lower polariton.

An alternative method of using the equations for polarized emission spectra (33) is to use experimental molecular absorption and fluorescence data. Then, to a good accuracy the lineshape of absorption is related to $\real(A)$ and fluorescence to $F$, as seen in Eq. (28). From the real part of $A$, the imaginary part may be found numerically by Hilbert transform (due to Kramers–Kronig relations). The response function is then determined from the plasmon eigenfrequency $\omega_{c}$ and its linewidth $\kappa$ together with the strong coupling constant $g_{N}$. Although $g_{N}$ (and the magnitude of $A$) is unknown, it can be fixed so that it corresponds to a given Rabi splitting. Lastly, the coupling coefficients $C_{s/p}$ and $\kappa_{o}^{T}$ are needed to evaluate the spectra. However, if we are only interested in the relative magnitudes, it is enough to fix the ratios $\kappa_{o}^{T}/C_{s}$ and $C_{p}/C_{s}$. The latter ratio is given by the polarization angle $\beta$ which can be evaluated with the dielectric functions of the materials at the interface where the plasmon is excited. The former ratio $\kappa_{o}^{T}/C_{s}$ is difficult to determine directly from the experiments but it can be found by fitting to experimental data.

In Fig. 5 we employ the above method to compare the experimental results for the polarization ratio of TDBC and R6G molecules from Ref. Baieva et al., 2017 to our model. Here, the polarization ratio $P^{T}_{p}/P^{T}_{s}$ of lower (upper) polariton is defined as the ratio of the p- and s-polarized emission peak intensity of lower (upper) polariton. In our numerical analysis, we approximate the fluorescence data by a mirror image of the absorption data over the zero-phonon frequency $\tilde{\omega}_{m}$. We assume that the fluorescence is effectively independent of the driving so that we can calculate and consider only the elastic emission (see also Appendix C). Adding the inelastic emission can only diminish the polarization ratio. Then we calculate the polarization ratio using Eq. (33) and fit the coupling rate ratio $\kappa_{o}^{T}/C_{s}$ to the experimental data for each branch separately.

For lower polariton peaks we find reasonable agreement with the fitted theoretical curves and the experimental data. For upper polariton the correspondence is very limited. Theoretically, we would assume that the polarization ratio increases for the upper polariton for positive detunings while for the lower polariton the ratio increases for negative detunings. This is caused by the polaritonic state becoming more plasmonic, which is seen from the response function being peaked at $\omega_{d}\approx\omega_{c}$ in Fig. 3. While these trends can be seen in Fig. 5, except for R6G upper polariton branch, some features differ from the theoretical description. From experimental point of view the polarization ratio might be affected by any external noise especially in regions of small s-polarized emission. In our modeling we neglect the possible dependence of the plasmon linewidth and coupling rate on the plasmon eigenfrequency. Also, since we use the experimental absorption data to determine the spectra the lineshape far from the absorption maximum is also important.

The fitted ratio $\kappa_{o}^{T}/C_{s}$ controls, generally speaking, the magnitude of the polarization ratio. The effect of Stokes shift seems to be important only at plasmon eigenfrequencies below the fluorescence frequency of the molecules while the external coupling rates control the polarization ratio at higher frequencies. From the fit to the experimental data we find that the ratio $\kappa_{o}^{T}/C_{s}$ is larger for R6G than for TDBC. This implies that if the plasmonic emission rate remains the same for TDBC and R6G samples, there is more molecular emission for TDBC than for R6G.

IV.2 Coherent polaritonic response

If we assume that the plasmon couples strongly to the molecules and drives them coherently, the phase factor $e^{i\vec{k}\cdot\vec{r}_{j}}$ is fixed to a constant. This leads to the introduction of two different sums over the molecular indices

[TABLE]

where again the upper/lower line is for s/p. The sum in $\tilde{C}_{s}$ vanishes because the plasmon polarization vector $\hat{u}_{pl}$ is orthogonal to the s-polarization vector $\hat{u}_{s}$ making the product antisymmetric under reflection through the plane orthogonal to $\hat{u}_{s}$. This is what breaks the symmetry between the polarization directions emitted from the strongly coupled mode.

Then, we find the polarized emission spectra to be

[TABLE]

The s-polarized emission vanishes now fully when the vibrations are absent. This result shows that coherence is crucial to the destructive interference of emitted light Baieva et al. (2017) while the vibrations still provide a mechanism for a partial loss of coherence. On the other hand, the p-polarized spectrum contains the interference terms between the plasmon and molecular output fields. Similar to the Stokes shift case, there are terms with different powers of the number $N$ of molecules. Considering the Rabi splitting (or $g_{N}$) to be fixed, there is one term with an extra $N$ factor from ${\tilde{C}_{p}}^{2}$ and $\sqrt{N}$ from $\tilde{C}_{p}$. Increasing $N$ leads to mostly elastic p-polarized emission, as the inelastic terms and s-polarized emission are independent of $N$. In contrast to the Stokes shift case, the absence of vibrations does not remove the $N$ dependence. Therefore, the result corresponds to superradiance Dicke (1954).

Figure 6 shows that the coherence and the resulting interference between plasmonic and molecular emission have a qualitative effect on the elastic spectra. The difference between s- and p-polarized spectra becomes more evident. While the s-polarized emission is likely to occur on the lower polariton frequency, for the p-polarized emission the upper polariton frequency may be favored depending on the relative magnitudes of $\kappa_{o}^{T},\tilde{C}_{p}$ and $C_{p}$.

V Conclusions

To summarize, we have constructed a model that allows describing the effect of vibrations on the strongly coupled stationary response of driven coupled light-matter modes. Depending on the case, one can either find the $P(E)$ function describing the absorption and emission of vibrations in a given model system, or relate the measured absorption and fluorescence of uncoupled molecules to $P(E)$. With small modifications, this approach can be extended also to the case of molecule-cavity systems Yu et al. (2019); Canaguier-Durand et al. (2013); Schwartz et al. (2013); Chovan et al. (2008); Virgili et al. (2011), plasmonic lattices Hakala et al. (2018) and/or higher-order correlation functions of the emitted light Wang et al. (2019). Our quantum Langevin equation approach allows describing the stationary driven system, and hence it complements the often-used computational methods usually concentrating on transient response del Pino et al. (2018); Groenhof et al. (2019).

Acknowledgements.

We acknowledge the support from the Academy of Finland Center of Excellence program (project no. 284594) and project numbers 289947, 290677, and 317118.

Appendix A Generalization to many nonidentical and interacting vibrational modes

The $P(E)$ theory is straightforward to generalize to multiple vibrational modes when the modes couple linearly to the molecule. A general interaction term in the Hamiltonian is then $\lambda_{ijkl}b_{ij}^{\dagger}b_{kl}+\text{h.c.}$ where $b_{ij}$ corresponds to the $j$th vibrational mode of the $i$th molecule. The molecule-vibrational Hamiltonian is then diagonalized by first diagonalizing the vibrational Hamiltonian and then using the polaron transformation. For simplicity, let us now discuss the case of a single molecule. After the diagonalization of the vibrational part, we may write the interaction Hamiltonian in terms of the new diagonal vibrational modes $b_{j}$ as

[TABLE]

where $x_{j}$ and $p_{j}$ are the position and momentum operator of vibrations. The term $u_{j}p_{j}$ follows from the fact that the molecule couples to the bare vibrational modes. Because in the Caldeira–Leggett model the position operator $x_{j}$ couples to the position operator of an environmental (harmonic) mode, the diagonalization is incommensurate with this model unless $u_{j}=0$. In the single-excitation limit we may then introduce the operator $\sigma^{S}=Q\sigma\equiv\prod_{j}Q_{j}\sigma$ where $Q_{j}=e^{\sqrt{S}\quantity(b_{j}^{\dagger}-b_{j})}$. Introducing many molecules into the situation only adds one external index to each operator. When there is no coupling to the plasmon, by following the same approximations as in the main text, we find that the dynamics is given by the input-output equation

[TABLE]

with $\tilde{\omega}_{m}=\omega_{m}-\sum_{j}S_{j}\omega_{v,j}$. The equation for $b_{j}$ again decouples from the dynamics of $\sigma^{S}$ in the single-excitation limit. Similarly to the case of a single vibrational mode we find the two- and four-point correlators of $Q$ in the calculation of the spectra. However, since $Q(t)=\prod_{j}Q_{j}(t)$, the Fourier transform of $Q$ is always a convolution. After diagonalization we may treat the modes as independent so this structure shows up as convolutions of $P(E)$’s and $L$’s defined in the earlier sections. Thus, we have

[TABLE]

where $P_{\rm tot}(E)=[P_{1}*P_{2}*\dots*P_{M}](E)$ is a convolution over $M$ different modes and similarly for $L_{\rm tot}$. For example,

[TABLE]

where $L_{1}$ and $L_{2}$ are defined for a single mode as the Fourier transform of Eq. (16).

Appendix B Approximation to the polaron equation

We discuss the consistency of the approximation that allows us to simplify the input-output equation of $\sigma^{S}_{j}$. The full dynamical equation of $\sigma^{S}_{j}$ for the strongly coupled plasmon–molecule system, using the approach of Gardiner and Collett (1985), is given by

[TABLE]

In the main text, we assumed the single-excitation limit in which $\sigma_{z,j}\approx-1$. In addition, we neglect the thermal fluctuations and set $\sigma_{{\rm in},j}=0$. Next we discuss when we can neglect the two last terms that are generated by the coupling of $\sigma^{S}_{j}$ to the the vibrational baths. This approximation effectively uncouples the vibrational dynamics from the dynamics of the polaron operator $\sigma^{S}_{j}$. As this approximation is related to the molecule–vibration system, also the coupling to the plasmon may be neglected, i.e. $g_{j}=0$. For notational brevity, we omit the molecular index $j$. Let us consider an expansion $\sigma^{S}=\sigma^{S}_{0}+\tilde{\sigma}^{S}_{1}$, where $\sigma^{S}_{0}$ is the solution of

[TABLE]

on Eq. (42) with the above simplifications. Consequently, the dynamics of $\tilde{\sigma}^{S}_{1}$ is given by

[TABLE]

We may now construct the next order of the expansion by setting $\tilde{\sigma}^{S}_{1}=\sigma^{S}_{1}+\tilde{\sigma}^{S}_{2}$ and fixing $\sigma^{S}_{1}$ to be the solution of

[TABLE]

This equation can be solved with the solution of $\sigma_{0}^{S}$. The dynamics of $\tilde{\sigma}^{S}_{2}$ is then determined by an equation similar to Eq. (44) where $\tilde{\sigma}^{S}_{1}$ is replaced by $\tilde{\sigma}^{S}_{2}$ and $\sigma^{S}_{0}$ by $\sigma^{S}_{1}$. Continuing this process gives then the expansion of $\sigma^{S}=\sum_{j=0}^{\infty}\sigma^{S}_{j}$. However, we focus only on the first order of the expansion.

Consider now that the molecule is driven coherently, $\sigma_{\rm in}=\alpha e^{-i\omega_{d}t}$ so that the solution of Eq. (43) is

[TABLE]

Consequently, we obtain from Eq. (45)

[TABLE]

where $*$ denotes a convolution in the Fourier space.

We are now interested in the consistency of the expansion but it is not straightforward to see the effect of the convolution and the underlying dynamics of the vibrations. For this reason, we compare the mean values of $\sigma^{S}_{0}$ and $\sigma^{S}_{1}$. When the input operators of the vibrations represent thermal noise, the $b_{\rm in}$ terms do not contribute to the average. The expectation value of $\sigma^{S}_{0}$ can be expressed as

[TABLE]

For a thermal ensemble $\expectationvalue{Q(0)}=\exp(-S(n_{\rm th}+\frac{1}{2}))$. In the calculation of the average of $\sigma^{S}_{1}$ we need the generalized Wick theorem to write

[TABLE]

where $J(t)$ is the function defined in Eq. (12) and we define its Fourier transform by $J(t)=\int\differential{\omega}e^{-i\omega t}J(\omega)$. We obtain

[TABLE]

Now we have a necessary condition for the consistency of the simplification: The parameter $\mathcal{C}$ should be small compared to unity for the expansion to be sensible. It can be estimated by using the same approximation in Eq. (19) as in the $\gamma=0$ calculation. Then (denoting $\Delta=\omega_{d}-\tilde{\omega}_{m}$)

[TABLE]

In the last approximation we have written the renormalized linewidth $\kappa_{m}$ in terms of the bare linewidth of the molecule $\tilde{\kappa}$ and neglected the smaller term $\chi(-2\omega_{v})$ for clarity. Now, it is clear that the consistency of the approximation is related to the temperature and the linewidths. This condition is always fulfilled when $n_{th}<1$ or alternatively $\frac{\omega_{v}}{k_{B}T}>\ln(2)\approx 0.69.$ It should be remembered that this is only a crude estimate and larger values of $\gamma$ can diminish the value of $\mathcal{C}$.

Appendix C Correspondence to experiments

In the experimental Kretschmann setup a prism is used together with white light to excite the plasmons. The angle $\theta$ of incoming light with respect to the normal of the interface then determines the plasmon eigenfrequency $\omega_{c}$. In our model, which is based on a coherent single frequency driving at $\omega_{d}$, we can introduce a distribution $\rho(\omega_{d})$ for light intensity. Then we can relate our theoretical model to the observed spectra by an integral relation

[TABLE]

The distribution $\rho(\omega_{d})$ can include features of the driving light as well as the prism that couples the light to the interface. We can generally divide the theoretical spectrum into elastic and inelastic part by $S^{T/R}_{s/p}(\omega;\omega_{d})=S_{\mathrm{el}}\delta(\omega-\omega_{d})+S_{\mathrm{inel}}(\omega;\omega_{d})$. Now, if we assume that the distribution of light is uniform and its bandwidth large compared to the plasmonic linewidth $\kappa$ we have

[TABLE]

The inelastic contribution is due to vibrations and molecular fluorescence (i.e. the function $F(\omega\neq\omega_{d})$). In the fit for Fig. 5, we disregard the contribution from $S_{\mathrm{inel}}$, because its main body is clearly separated from the elastic emission coming around the polariton frequencies.