Theory for the stationary polariton response in the presence of vibrations

We construct a model describing the response of a hybrid system where the electromagnetic ﬁeld—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 ﬂuorescence properties of the individual molecules, and thus are relevant in understanding any utilization of such systems, such as coherent light harvesting.


I. INTRODUCTION
Photonic structures, such as optical cavities or surface plasmon polaritons, can modify an 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 more strongly with the confined light mode. This interaction can become strong enough for the coupling energy to show up in the absorption and emission spectra of such systems, suggesting the formation of hybrid light-matter states, called polaritons. Common examples studied in this strong-coupling limit are single atoms [1], excitons in semiconductors [2], and photoactive molecules [3,4].
More recently, strong coupling of molecules with confined light modes has been 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 chemistry [5]. Experiments on strongly coupled molecules have already shown (i) suppression of photo-oxidation of TBDC J-aggregates coupled to plasmonic nanoprisms [6], and of photoisomerization of Spiropyran inside an optical cavity [7]; (ii) enhanced electronic conductivity in organic semiconductors [8]; (iii) intermolecular excitation energy transfer over large distances inside optical cavities [9,10]; and (iv) enhanced decay of triplet states in erythrosine B molecules [11]. Since polaritons are like interacting dressed photons with mass, they can undergo Bose-Einstein condensation even at room temperature [12], which further enables very efficient and thresholdless polariton lasing [13,14].
Strong coupling between a single molecule and electromagnetic field is very hard to achieve [15,16]. 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 for the two polaritonic states [17]. 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 environments [18]. However, those simulations often consider the transient response, whereas the 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 [19,20]. 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 [21]. 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 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 [22][23][24][25][26][27], as well as the varying polarization of that emission depending on the quantum coherence of the system, as shown recently [22].
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 lightmatter 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 1950s [28], 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 [29]. In this paper, we describe both the incoherent and coherent limits of polaritonics.
Although many aspects of our theory can be generalized to any confined light mode, such as 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 a metal-dielectric surface with a two-dimensional momentum k along the surface. In general, they have a nonlinear dispersion ω( 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 [3,23], leading to the observed avoided crossing between the two systems and thus offering the possibility to control photochemical reactions. A typical way to launch plasmons is via the Kretschmann configuration, i.e., coupling an external electromagnetic field to the surface modes via a prism [30]. In this setup, the angle with which the light enters the prism determines a specific plasmon k vector. Hence, in this work, we concentrate on a single plasmon mode with defined k and a generic frequency ω 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 ω c and a given polarizationû pl with respect to its wave vector 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 witĥ u pl = (0, sin β, cos β ) as in Fig. 1(a). The plasmon interacts with N identical molecules [31,32], which we approximate as two-level systems with transition frequency ω m . We denote the rising (lowering) operator of a molecule with σ † j (σ j ). As in typical experiments, we assume that the electric dipole moments of the molecules point in uniformly random directionsn j . Following the standard approach of quantum optics [33,34], the Hamiltonian of the strong-coupled system is in the rotating wave approximation (h = 1), The position 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 k· r j (n j ·û 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 ω v , but the generalization to multiple modes is straightforward (Appendix A). These vibrations and their interactions are described by The coupling between electronic and vibrational modes is quantified with a dimensionless parameter √ S, the Huang-Rhys factor [35], 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 [36,37]. We assume that there are separate bosonic baths for each molecule, vibration, and plasmon to which the coupling is linear in σ j , b j , and c, respectively. In the Markov approximation, these couplings are described by the dissipation ratesκ j , γ j , and κ of the molecules, vibrations, and plasmon. In the following, we suppose identical molecules and vibrations so that γ j = γ andκ j =κ. We neglect the thermal fluctuations of plasmons and molecules here ashω m ,hω c k B T even at room temperature. We simplify the molecule-vibration Hamiltonian by introducing a new polaron operator we assume a low driving power which corresponds to the single-excitation limit σ † σ ≈ 0. We find that when¯h ω v k B T > γ S 2κ m , where κ m =κ + γ 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 [38] for the vibrational dynamics. The plasmon and molecular equation are, in this case, whereω m = ω m − Sω v is the renormalized molecular frequency, while κ ext and κ ext m are the couplings to external driving fields. For the plasmon-molecule system, we assume that only the plasmon is driven so that c in = αe −iω d t and σ in, j = 0.
We model a measurement on the plasmon-molecule system so that the incoming light c in produces a reflected R out and transmitted T out field. 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û p =û y andû s = u x [ Fig. 1(a)]. 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, In this equation, δ T R = 0 and δ R R = 1, meaning that only the reflected field interferes with the input field. The δ s/p p is defined similarly because the plasmon couples only to ppolarized modes. The constants η T /R j,s/p describe the coupling of the molecule electronic states to the environmental s-and p-polarized free-space modes, and thus η T /R j,s/p = κ T /R m (n j · u s/p ). These fields and couplings to the system are represented schematically in Fig. 1(b). The output spectral density is obtained from where ω is the frequency of the output field and ω d is 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 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 κ, κ m , and γ .

II. P(E ) THEORY
The presence of the vibrations makes the input-output equations (3) nonlinear as they contain products of different dynamical fields. This nonlinearity 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 nonlinearity into account in the output spectra, we introduce P(E ) theory similar to the one in a dynamical Coulomb blockade [19]. Recently, the same problem has been discussed in Ref. [39] using similar methods, but only in a specific limit of vibrations (see below). The identification of P(E ) allows for a more general approach, also enabling the resolution of polariton emission, which is lacking from Ref. [39]. Let us first define P(t ) = Q † j (t )Q j (0) and its Fourier transform, 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 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 ω v , their linewidth γ , 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 γ vanishes. In this regime, P(E ) is related to the absorption function defined by Huang and Rhys [35]. However, our analytic results for the response also apply in the case of general γ and can be used for different models of vibrations.

A. Derivation of P(E )
We now derive the P(E ) function analytically in a similar manner as in the context of dynamical Coulomb blockade [19]. To establish notation, we omit the molecular index j here and denote , the dimensionless position and momentum operator, respectively. Then, we may write Q † (t ) = e i √ Sp(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. [19] and an example of a non-Gaussian P(E ) in Ref. [40]) 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 p(t ) = p(0) for stationary vibrations]. Consequently, for Gaussian fluctuations, we can write the characteristic function in terms of the variance alone. In that 245426-3 case [20], where the latter equality uses the fact that p(t ) 2 = p(0) 2 . The operator 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 the rotating wave approximation (also known as the Caldeira-Leggett model [38]), where γ is the linewidth of vibrations, and ξ is a Langevin force describing the thermal fluctuations. It has the correlator where the noise correlator is given by for thermal noise [41]. The Langevin equations (8) can be solved via Fourier transform. The result is After some Fourier analysis with the help of Eqs. (9)-(11), we find P(t ) = e J (t )−J (0) according to Eq. (7), with The resulting P(E ) is thus governed by three dimensionless parameters: the Huang-Rhys factor S, the quality factor of vibrations ω v /γ , and the relative temperature [39] for example.
Here, n th = (e ω v /(k B T ) − 1) −1 is the Bose factor, i.e., the mean number of thermal phonons at the vibrational frequency ω 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 γ ω v . This is hence the limit where Ref. [39] is valid. However, typical multiscale quantum chemistry calculations assume the opposite limit of a large γ κ, 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 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 ), 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.

B. 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, 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 where 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, Even if we can fully calculate J (t ), the Fourier transform of L is not straightforward to evaluate numerically in the general case.
C. γ = 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 ), given by Eq. (12), contains a nascent δ function, which in the limit γ → 0 reduces toP(ω) Note that this coincides with the limit γ → 0 in the whitenoise model given by Eq. (13). The corresponding characteristic function P(t ) is known in probability theory to be that of the Skellam distribution [42]. 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 where I k (x) is the modified Bessel function of the first kind. In the zero-temperature limit, p k (S) = e −S S k k! for k 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 γ = 0 distribution by using Eqs. (14) and (20), The variance depends on the temperature so that for high temperatures k B T ω v , the variance is directly proportional to the temperature: var(E ) ≈ 2Sω 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 γ → 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 ) ∝ S, we may find 1/P(t ) by changing S → −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 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 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 −∞ (or from 0 when T = 0) to ∞. This problem is alleviated by the rapid decrease of p k as a function of k. Consequently, Eq. (24) is straightforward 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 σ in, j = αe iθ j √ N δ(ω − ω d ) where θ 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 in = |α| 2 δ(ω − ω 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 T /R out = j ( κ T /R m σ j + δ T /R R σ 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 σ j = Q † j σ S j and the solution σ S j of Eq. (3b) depends on Q j , we encounter a four-point in the calculation of S T /R . Here, Q † j (ω) refers to the Fourier transform of Q † j (t ). Assuming that the vibration modes are independent and identical in different molecules, the correlator factorizes into two-point correlators when j = k. These resulting two-and four-point correlators are related to P(E ) by where L(ω 1 , ω 2 , ω 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 = ω d −ω m , which is the detuning between the driving and renormalized molecular frequency. Without vibrations, the molecular response is characterized by χ ( ) = (i − κ 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 and These functions are associated with absorption and fluorescence of molecules, respectively, and play an important role in the plasmon-molecule problem.

A. Incoherent limit
Let us assume that molecules are randomly arranged, so that the phase θ j is random. Averaging over them, the resulting spectra are The emission spectrum S T is determined by F , which describes inelastic scattering (output field frequency ω different from driving frequency ω d ). In the "reflected" field S R , we also find the input power spectral density I 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 γ → 0. Then, F is also δ peaked at frequencies ω − ω d = mω v with an integer m. The absorption spectrum is obtained from power conservation S A (ω d ) = I in − S T − S R evaluated at the driving frequency ω = ω d and it is mostly determined by A. In Fig. 2(b), we have plotted the emission spectrum S T along with the absorption spectrum S A . The absorption maximum is at the bare molecular frequency ω m , while the emission maximum is at approximately ω m − 2Sω v . The difference is the Stokes shift. The spectra correspond to the results describing the transient response obtained with Green functions [21]. However, in our stationary model, the absorption is not a mirror image of the emission because the emission may also happen from the excited vibrational states.

B. Coherent limit
Besides the experimentally more typical incoherent situation, we look at the coherent limit. Then the phase e iθ 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 σ in, j = α N e −iω 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 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 ) = δ(E ), F = |A( )| 2 δ(ω − ω d ) and the 1/N-dependent terms cancel. Therefore, this coherent effect is not related to the sub-or superradiance of molecules described by Dicke [28]. Rather, it is related to vibrations and their enhanced nonradiative emission, which shows up as a diminishing fluorescence as the number of molecules increases.

IV. PLASMON-MOLECULE SYSTEM
With the tools developed in Secs. II and III, we can return to the problem of a strongly coupled plasmon-molecule system and find the polarized spectra of the system using Eqs. (3) (with only the plasmon being driven, i.e., σ in, j = 0 in this case). We integrate Eq. (3b) from an initial time t i → −∞ to t f = t and neglect the initial condition σ S j (t i ) which has no role in a stationary situation. We substitute this into Eq. (3a), which leads tȯ 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 Q † j (t )Q j (t ) c(t ) = P(t − t )c(t ). Consequently, the elastic response of the plasmon is given by Vibrations provide a channel of relaxation, broadening the response which is associated with the real part of A. Finally, σ S j can be solved from Eq. (3b) in terms of Q j by Fourier transformation using the convolution theorem and c(ω) = αr(ω d )δ(ω − ω d ). Then we have all we need to evaluate the output spectra with Eqs. (4) and (5).

A. Incoherent polaritonic response
Let us consider a large number N of identical molecules with random dipole moment directionsn j . In this case, we can replace the sums over the molecule index with an integral over a surface of a sphere, j → N 4π d . Then, because g j = ge i k· r j (n j ·û pl ) whereû pl = (0, sin β, cos β ) is the plasmon polarization vector, the square of the Rabi splitting in Eq. (31) is j |g j | 2 = Ng 2 /3 ≡ g 2 N . We assume that the positions of the N molecules are random over a region that is large compared to the wavelength of the plasmon so that we may replace e i k·( r j − r k ) → δ jk for an ensemble average. Using these assumptions, the polarization dependence shows up in the spectra as the coefficients 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 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 β = 0.
Both the s-and p-polarized emission spectra are now represented with F and A found in molecular fluorescence given by Eq. (28). Therefore, the emission of a 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 FIG. 4. Elastic emission spectra for s and p polarization (red and blue curves, respectively) from a plasmon-molecule system with ω c = ω m . The different curves are offset and scaled for clarity. The parameters are the same as in Fig. 3, except for 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 ω d , the ratio of p-and s-polarized emission at ω = ω d is S T p /S T s = 2 − cos(2β ). For example, for an interface of vacuum and silver in the Drude model [34,43], β ≈ π 12 and S T p /S T s ≈ 1.13 at ω c = 2 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 κ T o and C s/p and the detuning ω c − ω m .
We find that the upper and lower polariton modes emit asymmetrically as in other approaches [39,[44][45][46][47][48][49][50] and experiments [22][23][24][25]51]. This is caused by the asymmetry of the effective dissipation rate = κ/2 − g 2 N Re[A( )], 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 hω v and S 1. Then, we may consider only singlephonon processes. When g N > ω v , we find the polariton frequencies from the response function (31) for zero detuning ω c = ω m to be approximately 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 κ m is given by Due to vibrations, the dissipation rate of the upper polariton ( + ) is larger than the dissipation rate of the lower polariton ( − ), which suppresses the upper polariton emission compared to the lower polariton. 245426-7 FIG. 5. Comparison of experimental polarization ratio data and corresponding theoretical fits for two different molecules. For TDBC, ω m ≈ 2.10 eV and 2Sω v ≈ 5 meV, and for R6G, ω m ≈ 2.27 eV and 2Sω v ≈ 97 meV. We estimate the plasmon linewidth to be κ = 250 meV. The Rabi splitting for TDBC is 167 meV and, for R6G, it is 337 meV. Our data on fluorescence of TDBC are limited below 2.6 eV so we cannot produce an estimate for the polarization ratio above 2.6 eV. In the legend, LP (UP) refers to lower (upper) 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 line shape of absorption is related to Re(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 ω c and its linewidth κ 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 κ T o are needed to evaluate the spectra. However, if we are only interested in the relative magnitudes, it is enough to fix the ratios κ T o /C s and C p /C s . The latter ratio is given by the polarization angle β, which can be evaluated with the dielectric functions of the materials at the interface where the plasmon is excited. The former ratio κ T o /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 the TDBC and R6G molecules from Ref. [22] to our model. Here, the polarization ratio P T p /P T s of the lower (upper) polariton is defined as the ratio of the p-and s-polarized emission peak intensity of the 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ω 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 κ T o /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 peaks, 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 ω d ≈ ω c in Fig. 3. While these trends can be seen in Fig. 5, except for the R6G upper polariton branch, some features differ from the theoretical description. From an 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 line shape far from the absorption maximum is also important.
The fitted ratio κ T o /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 κ T o /C s is larger for R6G than for TDBC. This implies that if the plasmonic emission rate remains the same for the TDBC and R6G samples, there is more molecular emission for TDBC than for R6G.

B. Coherent polaritonic response
If we assume that the plasmon couples strongly to the molecules and drives them coherently, the phase factor e i k· r j is fixed to a constant. This leads to the introduction of two different sums over the molecular indices, where again the upper/lower line is for s/p. The sum iñ C s vanishes because the plasmon polarization vectorû pl is orthogonal to the s-polarization vectorû s , making the product antisymmetric under reflection through the plane orthogonal toû 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 to the destructive interference of emitted light [22], 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 C 2 p and √ N fromC p . Increasing N leads to mostly elastic p-polarized emission, as the inelastic terms and spolarized 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 [28]. 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 κ T o ,C p , and C p .

V. CONCLUSIONS
To summarize, we have constructed a model that allows one to describe 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 also be extended to the case of molecule-cavity systems [4,17,[51][52][53], plasmonic lattices [54], and/or higher-order correlation functions of the emitted light [16]. Our quantum Langevin equation approach allows one to describe the stationary driven system, and hence it complements the often-used computational methods usually concentrating on transient response [18,45].

ACKNOWLEDGMENTS
We acknowledge the support from the Academy of Finland Center of Excellence program (Project No. 284594) and Projects No. 289947, No. 290677, and No. 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 λ i jkl b † i j b kl + H.c., where b i j corresponds to the jth vibrational mode of the ith 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 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 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σ withω m = ω m − j S j ω v, j . The equation for b j again decouples from the dynamics of σ 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 ) = 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 Q † (ω 1 )Q(ω 2 ) = P tot (ω 1 )δ(ω 1 + ω 2 ), (A3a) Q † (ω 1 )Q(ω 2 )Q † (ω 3 )Q(ω 4 ) = L tot (ω 1 , ω 2 , ω 4 ) 245426-9 where P tot (E ) = [P 1 * P 2 * · · · * P M ](E ) is a convolution over M different modes, and similarly for L tot . For example, 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 σ S j . The full dynamical equation of σ S j for the strongly coupled plasmonmolecule system, using the approach of [36], is given bẏ In the main text, we assumed the single-excitation limit in which σ z, j ≈ −1. In addition, we neglect the thermal fluctuations and set σ in, j = 0. Next we discuss when we can neglect the two last terms that are generated by the coupling of σ S j to the vibrational baths. This approximation effectively uncouples the vibrational dynamics from the dynamics of the polaron operator σ S j . As this approximation is related to the molecule-vibration system, the coupling to the plasmon may also be neglected, i.e. g j = 0. For notational brevity, we omit the molecular index j. Let us consider an expansion σ S = σ S 0 +σ S 1 , where σ S 0 is the solution oḟ on Eq. (B1) with the above simplifications. Consequently, the dynamics ofσ S 1 is given bẏ We may now construct the next order of the expansion by settingσ S 1 = σ S 1 +σ S 2 and fixing σ S 1 to be the solution oḟ This equation can be solved with the solution of σ S 0 . The dynamics ofσ S 2 is then determined by an equation similar to Eq. (B3), whereσ S 1 is replaced byσ S 2 and σ S 0 by σ S 1 . Continuing this process gives then the expansion of σ S = ∞ j=0 σ S j . However, we focus only on the first order of the expansion.
Consider now that the molecule is driven coherently, σ in = αe −iω d t , so that the solution of Eq. (B2) is Consequently, we obtain, from Eq. (B4), 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 σ S 0 and σ S 1 . When the input operators of the vibrations represent thermal noise, the b in terms do not contribute to the average. The expectation value of σ S 0 can be expressed as For a thermal ensemble, Q(0) = exp[−S(n th + 1 2 )]. In the calculation of the average of σ S 1 , we need the generalized Wick theorem to write S p n (t )p(0) = nJ (t ) p n−1 (t ) , where J (t ) is the function defined in Eq. (12) and we define its Fourier transform by J (t ) = dωe −iωt J (ω). We obtain Now we have a necessary condition for the consistency of the simplification: The parameter 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 γ = 0 calculation. Then (denoting = ω d −ω m ), In the last approximation, we have written the renormalized linewidth κ m in terms of the bare linewidth of the moleculeκ and neglected the smaller term χ (−2ω 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, ω v k B T > ln(2) ≈ 0.69. It should be remembered that this is only a crude estimate and larger values of γ can diminish the value of 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 θ of incoming light with respect to the normal of the interface then determines the plasmon eigenfrequency ω c . In our model,  which is based on a coherent single-frequency driving at ω d , we can introduce a distribution ρ(ω d ) for light intensity. Then we can relate our theoretical model to the observed spectra by an integral relation S T /R s/p,obs (ω) = dω d ρ(ω d )S T /R s/p (ω; ω d ). (C1) The distribution ρ(ω 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 parts by S T /R s/p (ω; ω d ) = S el δ(ω − ω d ) + S inel (ω; ω d ).
Now, if we assume that the distribution of light is uniform and its bandwidth large compared to the plasmonic linewidth κ, we have S T /R s/p,obs (ω) = S el (ω) + dω d S inel (ω; ω d ). (C2) The inelastic contribution is due to vibrations and molecular fluorescence [i.e., the function F (ω = ω d )]. In the fit for Fig. 5, we disregard the contribution from S inel because its main body is clearly separated from the elastic emission coming around the polariton frequencies.