Quantum description of surface-enhanced resonant Raman scattering within a hybrid-optomechanical model

TL;DR

This paper develops a quantum model of surface-enhanced resonant Raman scattering, revealing a regime of coherent interaction between molecular vibrations and electronic states mediated by plasmonic nanocavities, enabling control over vibrational modes.

Contribution

It introduces a quantum hybrid-optomechanical model for SERRS, highlighting the possibility of coherent vibrational-electronic interactions and vibrational mode control.

Findings

Existence of a coherent interaction regime between electronic and vibrational states.

Vibrational pumping can selectively activate specific vibrational modes.

Potential for controlling plasmon-assisted chemical processes.

Abstract

Surface-Enhanced Raman Scattering (SERS) allows for detection and identification of molecular vibrational fingerprints in minute sample quantities. The SERS process can be also exploited for optical manipulation of molecular vibrations. We present a quantum description of Surface-Enhanced Resonant Raman scattering (SERRS), in analogy to hybrid cavity optomechanics, and compare the resonant situation with the off-resonant SERS. Our model predicts the existence of a regime of coherent interaction between electronic and vibrational degrees of freedom of a molecule, mediated by a plasmonic nanocavity. This coherent mechanism can be achieved by parametrically tuning the frequency and intensity of the incident pumping laser and is related to the optomechanical pumping of molecular vibrations. We find that vibrational pumping is able to selectively activate a particular vibrational mode, thus providing a mechanism to control its population and drive plasmon-assisted chemistry.

Figures

Equations

$$H_{\mathrm{om}}=H_{\mathrm{pl}}+H_{\mathrm{vib}}+H_{\mathrm{pl}-\mathrm{vib}}+H_{\mathrm{pump}},$$

$$H_{\mathrm{pl}}$$

$$H_{\mathrm{vib}}$$

$$H_{\mathrm{pl}-\mathrm{vib}}$$

$$H_{\mathrm{pump}}$$

$$H_{\mathrm{bom}}=$$

$$\hslash \Omega (b^{†}+d_{\mathrm{om}}{a}^{†}a)(b+d_{\mathrm{om}}{a}^{†}a),$$

$$H_{\mathrm{bom}}\approx$$

$$H_{\mathrm{om}}^{\mathrm{res}}$$

$$H_{\mathrm{mol}}=$$

$$+\hslash \Omega (b{}^{†}+{\sigma }_{\mathrm{e}}d)(b+{\sigma }_{\mathrm{e}}d),$$

$$H_{\mathrm{pl}-\mathrm{e}}=$$

$${\mathcal{L}}_a[\rho ]$$

$${\mathcal{L}}_b[\rho ]$$

$${\mathcal{L}}_{\sigma }[\rho ]$$

$$n_T^{\mathrm{vib}}=\frac{1}{\exp \left(\frac{\hslash \Omega }{k_{\mathrm{B}}T}\right)-1},$$

$${\mathcal{L}}_{b^{†}}[\rho ]$$

$${\mathcal{L}}_{{\sigma }_z/2}(\rho )=-\frac{{\gamma }_{\varphi }}{4}\left(\{{{\sigma }_z}^{†}{\sigma }_z,\rho \}-2{\sigma }_z\rho {{\sigma }_z}^{†}\right),$$

$$\begin{array}{rl}{s}_{\mathrm{e}}(\omega )& \propto 2\mathrm{Re}{\int }_0^{\infty }⟨⟨a{}^{†}(0)a(\tau )⟩⟩\mathrm{d}\tau .\end{array}$$

$${\Gamma }_{\mathrm{eff}}\approx \frac{g^2{\gamma }_a}{\left[{\left(\frac{{\gamma }_a}{2}\right)}^2+(\delta +\Omega d^2-\Delta {)}^2\right]}\approx \frac{g^2{\gamma }_a}{\left[{\left(\frac{{\gamma }_a}{2}\right)}^2+(\delta -\Delta {)}^2\right]},$$

$$H_{\mathrm{red}}=\hslash \delta \sigma {}^{†}\sigma -\hslash \frac{1}{2}{\mathcal{E}}_{\mathrm{pl}}{\sigma }_x+\hslash \Omega (b{}^{†}+{\sigma }_{\mathrm{e}}d)(b+{\sigma }_{\mathrm{e}}d),$$

$$H_{\mathrm{vib}}=\hslash \Omega b^{†}b+\hslash d\Omega ⟨{\sigma }_{\mathrm{e}}⟩(b^{†}+b),$$

$${\mathcal{L}}_{\mathrm{eff}}^{-}[\rho ]=-\frac{{\Gamma }_{-}}{2}\left(b{}^{†}b\rho +\rho b{}^{†}b-2b\rho b{}^{†}\right)$$

$${\mathcal{L}}_{\mathrm{eff}}^{+}[\rho ]=-\frac{{\Gamma }_{+}}{2}\left(bb{}^{†}\rho +\rho bb{}^{†}-2b{}^{†}\rho b\right)$$

$$⟨b^{†}b{⟩}_{\mathrm{SS},\mathrm{in}}$$

$$\approx \frac{{\Gamma }_{+}}{{\gamma }_b}\propto \mathrm{Re}\{\tilde{S}(-\Omega )\},$$

$$\mathrm{Re}\{\tilde{S}(-\Omega )\}$$

$$\approx |⟨\sigma ⟩{|}^2\mathrm{Re}\left\{{\int }_0^{\infty }⟨⟨\sigma (\tau ){\sigma }^{†}(0)⟩⟩e^{-\mathrm{i}(\Omega -{\omega }_{\mathrm{L}})\tau }\mathrm{d}\tau \right\}$$

$$\approx \underset{{\tilde{S}}_{\mathrm{coh}}^{\mathrm{R}}}{\underbrace{\frac{|{\mathcal{E}}_{\mathrm{pl}}{|}^2}{4[{\delta }^2+({\Gamma }_{\mathrm{tot}}/2{)}^2]}}}\underset{{\tilde{S}}_{\mathrm{in}}^{\mathrm{R}}}{\underbrace{\frac{{\Gamma }_{\mathrm{tot}}/2}{(\delta +\Omega {)}^2+({\Gamma }_{\mathrm{tot}}/2{)}^2}}}.$$

$$H_{\mathrm{TLS}}=\hslash \frac{1}{2}\delta {\sigma }_z-\hslash \frac{1}{2}{\mathcal{E}}_{\mathrm{pl}}{\sigma }_x+\hslash \frac{1}{2}\delta .$$

$$\begin{array}{r}{\sigma }_z=\frac{\delta }{{\lambda }_{\mathrm{TLS}}}{\sigma }_z^{\prime }+\frac{{\mathcal{E}}_{\mathrm{pl}}}{{\lambda }_{\mathrm{TLS}}}{\sigma }_x^{\prime },\\ {\sigma }_x=-\frac{{\mathcal{E}}_{\mathrm{pl}}}{{\lambda }_{\mathrm{TLS}}}{\sigma }_z^{\prime }+\frac{\delta }{{\lambda }_{\mathrm{TLS}}}{\sigma }_x^{\prime }.\end{array}$$

$$\hslash \frac{1}{2}\Omega d{\sigma }_z(b^{†}+b)\to \underset{\text{Rabi term}}{\underbrace{\hslash \frac{1}{2}\Omega d\frac{{\mathcal{E}}_{\mathrm{pl}}}{{\lambda }_{\mathrm{TLS}}}{\sigma }_x^{\prime }(b{}^{†}+b)}}$$

$$+\text{residual polaronic coupling},$$

$$\Omega d\frac{{\mathcal{E}}_{\mathrm{pl}}}{{\lambda }_{\mathrm{TLS}}}⪆|3{\Gamma }_{\mathrm{tot}}/4+{\gamma }_b/2|.$$

$$H_{\mathrm{om}}^{\mathrm{res}}$$

$$H_{\mathrm{pl}}$$

Full text

††preprint: APS/123-QED

I Introduction

Surface plasmon excitations in metallic particles are able to squeeze and enhance electromagnetic fields down to the nanometric scale and thus dramatically enhance the interaction of nearby molecules with the incident light. The plasmonic near-field enhancement has been exploited in plasmon-enhanced spectroscopies, particularly in Surface-Enhanced Raman Spectroscopy (SERS)Kneipp et al. (1996, 1997); Xu et al. (1999); Haslett et al. (2000); Xu et al. (2000); Brolo et al. (2004); Maher et al. (2004); Haes et al. (2005); Moskovits (2005); Le Ru and Etchegoin (2006); Stiles et al. (2008); Tong et al. (2014); Zrimsek et al. (2017), which enables detection of minute quantities of molecular samples. The improved design of plasmonic cavities has allowed for spectroscopic investigation of even single molecules that are placed into ultranarrow plasmonic gaps Xu et al. (1999); Zhang et al. (2013). Current experimental strategies have taken advantage of the properties of plasmonic cavity modes that allowed reaching the plasmon-exciton strong-coupling regime with single molecules Chikkaraddy et al. (2016), as well as intramolecular optical mapping of single-molecule vibrations in SERS Zhang et al. (2013) or in electroluminescence Doppagne et al. (2017); Lee et al. (2019). These results suggest the possibility to push the use of plasmonic modes to further actively control the quantum state of a single molecule and thus influence its chemistry Stipe et al. (1997); Komeda (2005); Crim (1996); Ho (2002); Pascual et al. (2003); Hahna and Ho (2005); Herrera and Spano (2016). Recent theoretical and experimental studies Roelli et al. (2015); Schmidt et al. (2016, 2017); Kamandar Dezfouli and Hughes (2017); Benz et al. (2016) have revealed that off-resonant SERS can be understood as a quantum optomechanical process Aspelmeyer et al. (2014); Roelli et al. (2015); Schmidt et al. (2016) where the single plasmon mode (sustained in a plasmonic cavity) of frequency $\omega_{\rm pl}$ plays the role of the macroscopic optical cavity and the molecular vibration of frequency $\Omega$ plays the role of the macroscopic oscillation of the mirror. The description of such a process requires the development of concepts and methods beyond the standard classical description of SERS Roelli et al. (2015); Schmidt et al. (2016); Benz et al. (2016); Schmidt et al. (2017); Kamandar Dezfouli and Hughes (2017).

In this work we address a quantum mechanical theory of Surface-Enhanced Resonant Raman Scattering (SERRS), where an optical plasmonic mode supported by a metallic nanostructure mediates a coherent laser excitation of a nearby single molecule described as an electronic two-level system (TLS), coupled to a vibrational mode. In SERRS this laser excitation is assumed to be resonant with the electronic transition in the molecule. We discuss the similarities and differences between the SERRS Hamiltonian and the off-resonant quantum optomechanical Hamiltonian, which has been described previously. To that end we adopt a range of optomechanical parameters available in typical resonant situations. We then show that non-trivial phenomena emerge in SERRS under intense laser illumination, when the non-linearities of the molecule can trigger the coherent coupling of molecular electronic and vibrational degrees of freedom Xu et al. (2004); Johansson et al. (2005); del Valle et al. (2009); Delga et al. (2014a); Gu et al. (2010); Trügler and Hohenester (2008); Ridolfo et al. (2010). We further exploit the analogy with quantum optomechanics to propose a mechanism of on-demand frequency-selective pumping of molecular vibrations Kneipp et al. (1996); Haslett et al. (2000); Brolo et al. (2004); Maher et al. (2004); Le Ru and Etchegoin (2006) via the coherent laser illumination. These phenomena may provide a means to drive plasmon-enhanced vibrational spectroscopy to the realm of single-molecule selective chemistry or engineering of single-molecule optomechanical systems involving molecular vibrations on demand.

II SE(R)RS as an optomechanical process

Here we extend the analogy between optomechanics and SERS and describe SERRS in the framework of cavity quantum electrodynamics (QED) as a hybrid-optomechanical process Ramos et al. (2013); Akram et al. (2015). For further convenience we now briefly describe and compare the resonant and off-resonant scenarios.

II.1 Off-resonant SERS

In off-resonant SERS all electronic transitions in the molecule are far detuned from the frequency of the probing laser. The optomechanical Hamiltonian describing off-resonant SERS reads Roelli et al. (2015); Schmidt et al. (2016, 2017); Kamandar Dezfouli and Hughes (2017)

[TABLE]

where

[TABLE]

Here operators $a$ ($a^{\dagger}$) and $b$ ($b^{\dagger}$) are the annihilation (creation) operators for plasmons and vibrations, respectively, $\mathcal{E}$ is the pumping amplitude (with $\lvert\mathcal{E}\rvert^{2}$ proportional to the laser intensity) that characterizes the interaction Hamiltonian between the plasmon mode and the classical laser illumination of frequency $\omega_{\rm L}$, and $g_{\rm om}$ is the optomechanical coupling constant that can be connected to the Raman tensor of the molecule Roelli et al. (2015); Schmidt et al. (2017, 2016) and to the plasmonic near field. The Raman tensor contains the influence of the off-resonant electronic transitions in the molecule that are not considered explicitly in the Hamiltonian. The pumping term $H_{\rm pump}$ is considered in the rotating-wave approximation (RWA) assuming that $\mathcal{E}$ does not reach more than $\mathcal{E}\approx 0.1\omega_{\rm pl}$.

In what follows it will be convenient to interpret the optomechanical interaction as a displacement of the vibrational mode by a dimensionless value $n_{\rm pl}d_{\rm om}$, dependent on the number of excitations (plasmons) in the cavity, $n_{\rm pl}$. This becomes apparent after rearranging the bare optomechanical Hamiltonian $H_{\rm bom}=H_{\rm pl}+H_{\rm vib}+H_{\rm pl-vib}$ into the form

[TABLE]

with the dimensionless displacement $d_{\rm om}=-g_{\rm om}/\Omega$. The first line of Eq. (3) is the non-linear Hamiltonian of the cavity excitations (plasmons). In the limit of a weakly populated cavity and $d_{\rm om}\ll 1$ we can neglect the small nonlinear term $-\hbar\Omega d_{\rm om}^{2}a^{\dagger}a^{\dagger}aa$, and we re-define the plasmon frequency as $[\hbar\omega_{\rm pl}-\Omega d_{\rm om}^{2}]\equiv\hbar\omega^{\prime}_{\rm pl}$ and recover the linear plasmon Hamiltonian $\hbar\omega^{\prime}_{\rm pl}a^{\dagger}a$, so that:

[TABLE]

The second line in Eq. (3) has the sought form of a vibrational mode displaced by an amount that depends on the number of plasmonic excitations in the cavity $a^{\dagger}a$.

The level structure of the bare optomechanical Hamiltonian $H_{\rm bom}$ [Eq. (4)] is visualized in Fig. 1(a). The large grey dashed parabola illustrates an effective potential supporting the plasmonic mode. The vibrational potential is represented by the small parabolas that are displaced along the dimensionless normal coordinate $q$ by a magnitude $n_{\rm pl}d_{\rm om}$ proportional to the plasmonic number state $|n_{\rm pl}\rangle$Nunnenkamp et al. (2011, 2012). The energies of the plasmon Hamiltonian form an equidistant ladder, schematically drawn for the three lowest plasmon number states $|n_{\rm pl}\rangle$, and contain a fine-structure of molecular vibrational sub-levels.

II.2 SERRS

We complete the previous picture of SERS by addressing the scenario where the frequency of the incident laser approaches the molecular electronic resonance. To describe SERRS we can consider the molecule as a TLS that interacts with the vibrational modes via a polaronic coupling term Li et al. (1994); Herrera and Spano (2017); Galperin et al. (2009); Ramos et al. (2013); Ćwik et al. (2016); Zeb et al. (2018) and with the plasmonic cavity via the Jaynes-Cummings coupling term [see Fig. 1(b)]. The molecular vibrations are modeled as bosons within the Born-Oppenheimer approximation, where the effective harmonic vibrational potential is given by the ground state [$E_{\rm g}(q)$] and the excited state [$E_{\rm e}(q)$] potential energy surfaces (PESs) along a normal vibrational coordinate $q$, respectively Klessinger and Michl (1995); May and Kühn (2008); Galego et al. (2015, 2016). We consider that the vibrational energies of the molecule, $\hbar\Omega$, are the same for the ground and for the excited state.

The Hamiltonian of the SERRS system $H_{\rm om}^{\rm res}$, using the rotating wave approximation, can be expressed as Betzholz et al. (2014); Ramos et al. (2013); Jaehne et al. (2008); del Pino et al. (2015):

[TABLE]

with

[TABLE]

and $H_{\rm pl}$, $H_{\rm pump}$ as defined in Eq. (2) above. Here the operator $\sigma$ ($\sigma^{\dagger}$) is the lowering (raising) operator of the TLS, with $\sigma_{\rm e}=\sigma^{\dagger}\sigma$ the TLS number operator. The parameter $d$ is the dimensionless displacement between the minima of the ground- and excited-state PESs, which is related to the Huang-Rhys factor Klessinger and Michl (1995); May and Kühn (2008); Ćwik et al. (2016); Herrera and Spano (2017); del Pino et al. (2015), $S$, as $S=d^{2}$, and is a measure of the coupling between the molecular vibration and the excitonic transition. The interaction of the localized plasmon excitation and the molecular electronic levels is mediated by the plasmon-exciton coupling constant $g$. More details about the model Hamiltonian are given in Appendix A.

The level diagram describing the SERRS Hamiltonian $H_{\rm om}^{\rm res}$ [Eq. (5)] is sketched in Fig. 1 (b). Strikingly, both the off-resonant Hamiltonian, $H_{\rm bom}$, and the molecular Hamiltonian in SERRS, $H_{\rm mol}$, can be represented as a series of mutually displaced harmonic vibrational PESs. The electronic states in SERRS thus play the role of the plasmon number states in the off-resonant case, an analogy which can be identified from the comparison of the Hamiltonians in Eq. (6) and Eq. (4).

In the limit of single-photon optomechanics Nunnenkamp et al. (2011, 2012); Akram et al. (2010), where the plasmon Hilbert space is limited to the vacuum state and the singly excited state, the molecular Hamitonian $H_{\rm mol}$ and $H_{\rm bom}$ become formally identical. However, as we detail later, if the incident laser is strong, the non-linear character of the excitonic TLS Hamiltonian [Eq. (6)] will lead to novel physical phenomena that cannot be achieved in the off-resonant SERS situation [Eq. (4)]. This perspective makes it particularly attractive to study SERRS.

II.3 Dissipative processes in SE(R)RS

In realistic systems, excitations undergo decay, pumping and dephasing processes. Losses and thermal pumping are considered in the dynamics of the system by solving the master equation for the density matrix, $\rho$, with incoherent damping introduced via the Lindblad-Kossakowski terms Breuer and Petruccione (2003) for the plasmon and the vibration

[TABLE]

and in case of SERRS for the electronic TLS

[TABLE]

where $\gamma_{b}$ is the vibrational, $\gamma_{a}$ the plasmonic and $\gamma_{\sigma}$ the electronic decay rate, respectively, and $n^{\rm vib}_{T}$ is the equilibrium thermal population determined by the reservoir according to the Bose-Einstein distribution:

[TABLE]

with the Boltzmann constant $k_{\rm B}$ and thermodynamic temperature $T$. Furthermore, for finite temperatures we consider the thermal pumping of the vibrations via:

[TABLE]

In most of the paper we consider the low-temperature limit ($T=0$ K) where the thermal populations can be neglected. Finite temperatures are considered in Section VIII where we study the effects of finite temperature on the steady-state vibrational populations. We further include into the model the pure dephasing of the molecular electronic excitations in the form of the Lindblad term Breuer and Petruccione (2003); del Valle et al. (2009); Esteban et al. (2014).

[TABLE]

with $\sigma_{z}=\sigma^{\dagger}\sigma-\sigma\sigma^{\dagger}$. This description of pure dephasing is valid in the plasmon-exciton weak-coupling regime considered below Neuman and Aizpurua (2018).

The decay of the state-of-the-art plasmonic cavities $\gamma_{a}$ is ultimately limited by the material properties of the metal Wang and Shen (2006); Delga et al. (2014b); Kristensen and Hughes (2013); Sauvan et al. (2013); Koenderink (2010), reaching quality factors $Q$ ($Q=\omega_{\rm pl}/\gamma_{a}$) of up to $Q\approx 50$ ($\hbar\gamma_{a}\approx 40$ meV). However, the values of $Q$ commonly achieved in plasmonic systems are usually smaller ($Q\approx 1-20$). For example, the leaky gap mode formed between a tip of a scanning tunneling microscope and a metallic substrate often used for single-molecule spectroscopy Zhang et al. (2013, 2016); Doppagne et al. (2017); Imada et al. (2017); Zhang et al. (2017); Lee et al. (2019) can be strongly damped, $\hbar\gamma_{a}\sim 10^{2}$ meV, and can thus be regarded as a low-$Q$ plasmonic cavity.

On the other hand, a typical decay rate of molecular excitations decoupled from the plasmonic cavity is much smaller than that of plasmons (as small as $\hbar\gamma_{\sigma}\sim 10^{-2}$  meV $\ll\gamma_{a}$). The linewidth of the molecular resonance is thus mostly limited by the pure dephasing $\gamma_{\phi}$ which strongly scales with temperature and is highly dependent on the environment surrounding the molecule. It is thus possible to engineer conditions (low-temperature vacuum experiment) under which the pure dephasing becomes small and the linewidth of the molecular electronic excitation decreases below $<10$ meV and may even be limited only by the spontaneous decay.

II.4 Setting the SERRS regime

Let us consider a value of $\gamma_{a}$ representing a leaky gap plasmonic mode formed between a tip of a scanning tunneling microscope and a metallic substrate as typically used for single-molecule spectroscopy Doppagne et al. (2017); Imada et al. (2017); Zhang et al. (2017, 2013, 2016), for which $\hbar\gamma_{a}\sim 10^{2}$ meV. We consider the bad-cavity limit (weak-coupling regime) where the plasmon-exciton coupling $g$ is small compared to the plasmonic losses, but large with respect to the intrinsic decay $\gamma_{\sigma}$ ($\gamma_{\sigma}\ll g\ll\gamma_{a}$). In the bad-cavity limit, the molecular levels are only weakly perturbed by the presence of the plasmon that boosts the decay of the molecular electronic excitation via the Purcell effect, and focuses the incident light on the resonant molecule via the plasmonic near-field enhancement. The limit of the plasmon-exciton strong-coupling regime, where the coupling between the plasmon and the TLS dominates over the plasmonic losses, will be detailed elsewhere.

In the bad-cavity limit, the parameters determining the regime of the off-resonant optomechanical coupling, $d_{\rm om}$, as well as that defining the exciton-vibration coupling, $d$, in the resonant model, describe formally the same physical phenomena under weak-illumination conditions. This follows from the formal similarity between $H_{\rm bom}$ [Eq. (4)] and $H_{\rm mol}$ [Eq. (6)] established in Subsection II.2. In off-resonant SERS, the condition $|d_{\rm om}\Omega|>\gamma_{a}/2$ sets the so-called optomechanical strong coupling. In such situation the optomechanical non-linearity $-\Omega d_{\rm om}^{2}a^{\dagger}a^{\dagger}aa$ becomes important and the system becomes interesting for quantum applications. It has been estimated that the optomechanical coupling can reach up to $d_{\rm om}\sim 10^{-1}$ for some molecular species Roelli et al. (2015).

On the other hand in SERRS, for relevant dye molecules with electronic excitations in the visible, $d$ ranges from $d\sim 0.1$ for rigid molecules (such as porphyrines Huang et al. (2005)) up to values of $d\sim 1$ for soft organic molecules Herrera and Spano (2017); Sánchez-Carrera et al. (2010). SERRS might thus offer relatively high optomechanical coupling strengths even for a single organic molecule. Moreover, under the conditions of small molecular dephasing the width of the excitonic resonance becomes much smaller than that of the plasmon and SERRS may open the possibility to achieve large optomechanical coupling compared to the relevant linewidth $|d\Omega|>\Gamma_{\rm tot}/2+\gamma_{\phi}$ (strong optomechanical coupling), with $\Gamma_{\rm tot}=\gamma_{\sigma}+\Gamma_{\rm eff}$ and $\Gamma_{\rm eff}$ the width of the TLS due to the Purcell effect induced by the plasmons that we discuss later. In off-resonant SERS it is the cavity width $\gamma_{a}$ which determines the regime of optomechanical coupling and achieving the condition $|d_{\rm om}\Omega|>\gamma_{a}/2$ is more challenging.

In this work we consider a range of relatively large coupling strength between the plasmon and the molecular TLS ($\hbar g\approx$9 meV to $\hbar g=$50 meV), although still small enough to be in the bad-cavity limit Chikkaraddy et al. (2016); Trügler and Hohenester (2008). These selected values allow us to explore different regimes of plasmon-assisted interaction between molecular excitons and vibrations, as detailed below. We also adopt $\hbar\gamma_{a}=500$ meV, $\hbar\gamma_{\sigma}=0.02$ meV, $\hbar\omega_{\rm pl}=2$ eV, and $\hbar\omega_{\rm eg}=2\,$eV, as typical representative values of realistic molecules and plasmonic systems Trügler and Hohenester (2008); del Pino et al. (2015). We assume vibrational frequencies ranging between $\hbar\Omega=10$ meV and $\hbar\Omega=50$ meV, describing low-energy vibrations of typical organic molecules. In particular, when discussing phenomena emerging under weak laser illumination (as discussed below) we adopt the value of $\hbar\Omega=50$ meV and $\hbar\gamma_{b}=2$ meV to describe the molecular vibration, whereas to address the regime of strong laser illumination we choose $\hbar\Omega=10$ meV and $\hbar\gamma_{b}=1$ meV. This lower value of the vibrational frequency allows us to access the non-linear response of the system even for realistic laser intensities that comply with the rotating-wave approximation assumed in our model. With the exception of Section VIII, throughout the manuscript we consider zero ambient temperature $T=0$ K.

In the following we describe the inelastic emission spectra and the vibrational pumping in SERRS for (i) the linear response regime (relatively weak laser illumination in Secs. III and IV) and (ii) strong laser illumination where the molecular levels are dressed by the intense laser field and form a qualitatively new set of light-matter states (Secs. V and VI).

III Photon emission spectra in the linear regime

We first discuss the spectral response and the physics of hybrid-optomechanical vibrational pumping in SERRS systems in the limit of weak incident laser intensities, for which the system can be treated within the linear-response theory (further denoted as weak illumination). In this regime we limit the description to $T=0$ K as the thermal effects would mask the optomechanical pumping and damping processes that are weak for the low laser intensities discussed here. We relax this assumption later when we address the regime where the system is pumped by an intense laser.

For convenience, we define the detuning parameters $\Delta=\omega_{\rm pl}-\omega_{\rm L}$ and $\delta=\omega_{\rm eg}-\omega_{\rm L}$ (with the exciton frequency $\omega_{\rm eg}=E_{\rm e}(d)-E_{\rm g}(0)$) and define the coherent amplitude of the plasmon annihilation operator induced by the incident monochromatic illumination $\alpha_{\rm S}=\frac{-\mathcal{E}}{\Delta-{\rm i}\gamma_{a}/2}$. The solution of the dynamics of the hybrid optomechanical Hamiltonian and the respective Lindblad terms with the parameters described above, allow for calculating the steady-state emission spectrum $s_{\rm e}(\omega)$ from the plasmonic cavity using the quantum regression theorem (QRT):

[TABLE]

In this spectrum we remove the elastic scattering contribution and use the notation $\langle\langle O_{1}O_{2}\rangle\rangle=\langle O_{1}O_{2}\rangle-\langle O_{1}\rangle\langle O_{2}\rangle$, where $O_{1}$ and $O_{2}$ are operators. More details about the specific implementation of this expression can be found in Appendix B.

To illustrate the emission properties of typical molecules we plot in Fig. 2 (a) the inelastic spectra in a SERRS system [normalized to $s_{\rm e}(\hbar\omega_{\rm eg})$ and vertically shifted], calculated for weak illumination, $\hbar\mathcal{E}=1\times 10^{-2}$ meV (roughly corresponding to a laser power density of $W\approx 1\times 10^{-4}$ $\mathrm{\SIUnitSymbolMicro}\mathrm{W}\mathrm{/}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}$ ), from a monochromatic laser of frequency $\hbar\omega_{\rm L}=1.975$ eV (green dashed line) and for exciton-plasmon coupling $\hbar g=13$ meV. The excitonic energy is $\hbar\omega_{\rm eg}=2$ eV. We calculate the spectra for two large values of $d=1,\,0.5$ (top spectra) representing soft organic molecules and for a small value of $d=0.1$ of a rigid molecule (two bottom spectra). We choose $\hbar\gamma_{\phi}=20$ meV for the three top spectra to demonstrate the effect of pure dephasing on the emission of molecules interacting with a decoherence-inducing environment. In the bottom spectrum no dephasing is considered, $\hbar\gamma_{\phi}=0$ eV.

The bottom spectrum in Fig. 2 (a), calculated for weak exciton-vibration coupling, $d=0.1$, and considering no pure dephasing ($\hbar\gamma_{\phi}=0$ meV), features two sharp emission peaks. The fluorescence peak appears at frequency $\omega=\omega_{\rm eg}$ regardless of the incident laser frequency. The second peak, appearing at $\omega=\omega_{\rm L}-\Omega$, is the Raman-Stokes emission line. The anti-Stokes line is not visible because the vibrations are not thermally populated for $T=0$ K. The Raman (SERRS) line always appears at a constant detuning from the laser frequency which facilitates its identification in the spectrum. When the pure dephasing is increased, the fluorescence line starts to broaden and also increases in absolute intensity (the latter is not manifested in Fig. 2 (a) due to normalization). The SERRS emission becomes hardly distinguishable on top of the strong fluorescence background for $d=0.1$. As $d$ increases, the fluorescence background becomes asymmetrical and broadens towards lower energies due to radiative transitions allowed by the simultaneous exchange of energy between electronic and vibrational states (hot luminescence). This so-called vibrational progression of the luminescence spectrum thus consists of a series of broad peaks, each peak positioned at frequency $\omega_{\rm eg}-n\Omega$ (with $n$ a positive integer), with its amplitude determined by the overlap of the vibrational wave functions in the electronic ground and excited states, respectively (Franck-Condon factors) Nunnenkamp et al. (2011, 2012); Akram et al. (2010); Klessinger and Michl (1995); May and Kühn (2008). The Raman-Stokes lines appear on top of the fluorescence peaks at frequencies $\omega_{\rm L}-n\Omega$. The strength of the Raman lines is determined from a combination of the Franck-Condon overlaps, as in the case of the hot luminescence, and from the further enhancement due to the proximity of the molecular electronic resonance that (i) enhances the interaction of the incident laser with the molecular transition and (ii) boosts the efficiency of the Raman-Stokes emission. The Raman peaks are also notably narrower than the fluorescence peaks when dephasing is large ($\hbar\gamma_{\phi}=20$ meV), which facilitates their identification on top of the broad and intense fluorescence background.

The difference between the physical origin of the SERRS lines and the fluorescence lines becomes clearer if we plot the emission spectra as a function of the laser detuning from the exciton frequency, $\delta=\omega_{\rm eg}-\omega_{\rm L}$. The emission spectra are shown in Fig. 2 (b) for a rigid molecule ($d=0.1$) and for no dephasing, similar to the bottom spectrum in Fig. 2 (a). The molecule is pumped by an incident laser of amplitude $\hbar\mathcal{E}=1\times 10^{-2}$ meV and we consider an exciton-plasmon coupling, $\hbar g=13$ meV. The incident laser frequency is marked in the spectra as a white dashed line diagonally crossing the color plot. The color plot with the spectra shows both the Raman-Stokes peaks that appear red-detuned from the incident-laser frequency $\omega_{\rm L}$ by $n\Omega$, and the fluorescence peaks emerging at the energy of the excitonic transition regardless of $\delta$.

Furthermore, the emission shown in Fig. 2(b) is enhanced when the exciton is resonantly pumped, $\hbar\delta=0$ eV, or when when the frequency of the first-order Raman-Stokes line, $\omega=\omega_{\rm L}-\Omega$, corresponds to the bare excitonic resonance, $\hbar\delta=-\hbar\Omega=-50$ meV. When the laser frequency is tuned to the molecular exciton ($\hbar\delta=0$ eV), the incident laser coherently (coherent population $n_{\sigma}^{\rm coh}=|\langle\sigma\rangle|^{2}$) and incoherently (incoherent population $n_{\sigma}^{\rm incoh}=\langle\sigma^{\dagger}\sigma\rangle-|\langle\sigma\rangle|^{2}$) populates the electronic excited state. Thereafter, the molecule efficiently emits both the Raman-Stokes ($\propto n_{\sigma}^{\rm coh}$) and the hot luminescence ($\propto n_{\sigma}^{\rm incoh}$) photons at $\omega=\omega_{\rm eg}-\Omega$. On the other hand, when $\hbar\delta=-\hbar\Omega=-50$ meV, the spectral position of the first-order Raman-Stokes line coincides with the resonance frequency of the molecular exciton. In this case, the molecular fluorescence peaks, now appearing at the spectral positions of the Raman-Stokes lines, are suppressed since the off-resonance illumination does not efficiently populate the excited electronic state and the emission peaks appear mainly due to the SERRS mechanism. Both of these mechanisms of SERRS enhancement are closely related to the process of optomechanical vibrational pumping, described for the off-resonant case Roelli et al. (2015); Schmidt et al. (2016, 2017); Kamandar Dezfouli and Hughes (2017).

Finally we remark that the spectral map in Fig. 2(b) also features lines appearing due to higher-order Raman scattering and hot luminescence. These lines show much lower intensity than the lines of lower orders, but they exhibit the same mechanism of emission enhancement, as is apparent from the spectral map.

IV Vibrational pumping in the linear regime

Inasmuch as the emission of Raman-Stokes photons is accompanied by the creation of a vibrational quantum, the enhanced Raman-Stokes emission is reflected in the incoherent steady-state vibrational population $\langle b^{\dagger}b\rangle_{\rm SS,in}$. To elucidate the role of Raman-Stokes scattering in the process of vibrational pumping in SERRS, we plot in Fig. 2 (c) the vibrational population as a function of incident laser detuning for two values of plasmon-exciton coupling. The upper panel corresponds to $\hbar g=\hbar\sqrt{\gamma_{b}\gamma_{a}/6}\approx 13$ meV, for which the broadening of the electronic resonance due to the plasmonic Purcell effect (see also Appendix C),

[TABLE]

becomes comparable to the broadening of the vibrational line $\gamma_{b}$, and in the lower panel we use $\hbar g=50$ meV, ensuring that $\Gamma_{\rm eff}>\gamma_{b}$ (with $\hbar\Gamma_{\rm eff}\approx 20$ meV). In the calculations we set $\hbar\mathcal{E}=1\times 10^{-2}$ meV to make sure that we stay in the linear regime. The numerically calculated values of $\langle b^{\dagger}b\rangle_{\rm SS,in}$ (black lines) are qualitatively similar in both cases. A set of peaks are clearly observed which correspond to the enhancement of the Raman-Stokes emission for detunings of $\hbar\delta=0$ eV, $\hbar\delta=-\hbar\Omega=-50$ meV, and higher orders ($\delta=-n\Omega$, $n>1$), respectively. The effect of the larger plasmon-exciton coupling $g$ is to broaden the peaks and to smear off the population maxima associated with enhancement of the higher-order Raman-Stokes emission ($\delta=-n\Omega$, $n>1$).

To shed light on the mechanism of vibrational pumping in SERRS, we derive the effective vibrational dynamics which results from the elimination of the plasmon and the TLS dynamics, following standard methods from the theory of open quantum systems Breuer and Petruccione (2003), in close analogy with the procedure developed in hybrid quantum optomechanics Jaehne et al. (2008) (see description of the procedure in Appendix C and Appendix D). Upon elimination of the plasmon, the effective reduced TLS-vibrational Hamiltonian, $H_{\rm red}$, becomes

[TABLE]

where $\mathcal{E}_{\rm pl}=-2g\alpha_{\rm S}$ is the coherent pumping of the molecule mediated by the plasmon, and $\sigma_{x}$ is the Pauli $x$ operator. Moreover, in the bad-cavity limit, the molecular excitonic TLS is effectively broadened due to the plasmon via the Purcell effect. The total TLS decay rate thus becomes $\gamma_{\sigma}\rightarrow\Gamma_{\rm tot}=\Gamma_{\rm eff}+\gamma_{\sigma}$.

By further eliminating the TLS from the vibrational dynamics, assuming that the broadening of the electronic levels $\Gamma_{\rm tot}$ is larger than $|d\Omega|$ so that the Markovian approach applies, we obtain an effective vibrational Hamiltonian that includes the coherent pumping due to the TLS excited-state population

[TABLE]

which is accompanied by the effective incoherent damping $\Gamma_{-}$ and pumping $\Gamma_{+}$ rates (see explicit expressions in Appendix D), which need to be added to the intrinsic vibrational dissipation rate [described by the original Lindblad term in Eq. (8)] via the following (new) Lindblad terms:

[TABLE]

for the effective damping, and

[TABLE]

for the effective pumping. These rates are defined as $\Gamma_{-}=2\left(\Omega d\right)^{2}\mathrm{Re}\{\tilde{S}(\Omega)\}$ and $\Gamma_{+}=2\left(\Omega d\right)^{2}\mathrm{Re}\{\tilde{S}(-\Omega)\}$, respectively, where $\mathrm{Re}$ indicates the real part, and $\tilde{S}(s)=\int_{0}^{\infty}\langle\langle\sigma_{\rm e}(\tau)\sigma_{\rm e}(0)\rangle\rangle e^{{\rm i}s\tau}{\rm d}\tau$ is the spectral function corresponding to the one-sided Fourier transform of the correlation function of the TLS for a generic frequency $s$. The latter is calculated for the TLS decoupled from the vibrations, but coupled with the plasmon, which effectively broadens the TLS (details about the analytical calculation of $\tilde{S}(s)$, and thus $\Gamma_{+}$ and $\Gamma_{-}$, are provided in Appendix D). Finally, the incoherent steady-state vibrational population induced by the effective pumping of the vibrations via the TLS in this approximation becomes

[TABLE]

where the last approximation originates from the fact that under weak pumping $\gamma_{b}\ll\Gamma_{-}-\Gamma_{+}$. We note that in the linear regime $\langle b^{\dagger}b\rangle\approx\langle b^{\dagger}b\rangle_{\rm SS,in}$ which allows for direct comparison of the numerically calculated and the analytically obtained populations. From Eq. (19) it follows that the behavior of the spectral function $\tilde{S}(-\Omega;\delta)$ as a function of the incident laser frequency (i.e. $\delta$) determines the conditions for which the vibrational pumping occurs. In the linear regime we can simplify the expression for $\mathrm{Re}\{\tilde{S}(-\Omega)\}$, in analogy with the description of the off-resonant model Roelli et al. (2015); Schmidt et al. (2016, 2017); Kamandar Dezfouli and Hughes (2017) as:

[TABLE]

The two terms, $\tilde{S}^{\rm R}_{\rm coh}\approx|\langle\sigma\rangle|^{2}$ and $\tilde{S}^{\rm R}_{\rm in}\approx\mathrm{Re}\{\int_{0}^{\infty}\langle\langle\sigma(\tau)\sigma^{\dagger}(0)\rangle\rangle e^{-{\rm i}(\Omega-\omega_{\rm L})\tau}{\rm d}\tau\}$, can then be interpreted as the efficiency of the coherent driving ($\tilde{S}^{\rm R}_{\rm coh}$ resonant at $\delta=0$) and the efficiency of the spontaneous Stokes-Raman emission ($\tilde{S}^{\rm R}_{\rm in}$ resonant at $\delta=-\Omega$), respectively.

The effective vibrational dynamics and steady-state values derived in this section are an accurate approximation to the exact problem only if the decay rate, $\Gamma_{\rm tot}$, of the dressed TLS is significantly larger than the intrinsic vibrational decay rate $\gamma_{b}$, and the exciton-vibration coupling is weak (moderate values of $d$). Realistic situations in molecular spectroscopy might not satisfy these conditions and in those situations it would be necessary to adopt a numerical treatment to obtain accurate results. Nevertheless, the properties of the TLS spectral function reveal the origin of the vibrational pumping even beyond the limits of validity of the analytical-model.

We plot the analytical result for the evolution of the vibrational populations as a function of detuning $\delta$ with a red dotted line in Fig. 2 (c). The analytical vibrational populations share with their numerically calculated counterparts (black lines) the same dominant peaks appearing for zero laser detuning from the exciton frequency, $\delta=0$, and for detuning $\delta=-\Omega$, when the frequency of the first-order Raman-Stokes line coincides with the excitonic frequency. These two values of the laser detuning also lead to an enhancement of the Raman-Stokes emission [Fig. 2 (b)], which is here the driving mechanism of the optomechanical vibrational pumping.

Although the analytical model nicely describes the main features of the fully numerically calculated vibrational populations, it cannot explain the presence of the weaker higher-order peaks. This is due to the Markov approximation leading to Eq. (16), Eq. (17) and Eq. (18) which treats the exciton-vibration interaction perturbatively. In the full model, the vibrational pumping mechanism is also present for the vibrational transitions responsible for higher-order Raman scattering and hot luminescence. Moreover, in the case that $\hbar g=13$ meV, the analytical model overestimates the vibrational populations induced by the optomechanical amplification for $\hbar\delta=-\hbar\Omega=-50$ meV, since the effective broadening of the TLS, $\Gamma_{\rm tot}=\Gamma_{\rm eff}+\gamma_{\sigma}$, is similar to the vibrational broadening $\gamma_{b}$, and the Markov approximation becomes less accurate. For $\hbar g=50$ meV the effective broadening $\Gamma_{\rm tot}>\gamma_{b}$ and the analytical model describes the low-order features of the vibrational populations accurately.

V Photon emission spectra for strong laser intensities

Let us explore in the following the regime where the system is illuminated by a strong-power incident laser, which induces the non-linear response of the molecule, and thus requires a treatment beyond the standard optomechanical description applicable to off-resonant SERS. To observe such non-linear effects already for realistic values of intensity of the resonant laser illumination (and retaining the validity of the RWA) in this section we further consider values of vibrational frequencies $\hbar\Omega=10$ meV at the lower end of a typical vibrational spectrum of an organic molecule.

V.1 Influence of laser intensity

The influence of the incident laser amplitude $\mathcal{E}$ is shown in Fig. 3 (a,b), where color maps of the emission spectra are displayed as a function of $\mathcal{E}$. The incident laser frequency is tuned to the TLS electronic transition, and we consider the results for $d=0.2$, which show the well-known Mollow triplet, a spectral structure resulting from the resonance fluorescence (RF) of the dressed TLS Mollow (1969); Gu et al. (2010); Ridolfo et al. (2010); Ge et al. (2013); Wrigge et al. (2008); Muñoz et al. (2018). The Mollow triplet consists of a strong emission line centered at the incident laser frequency and two side spectral peaks of similar spectral width that shift away as the laser intensity is increased [Figs. 3 (a) and (b)]. At a specific pumping amplitude $\mathcal{E}$ [white dashed lines in Fig. 3 (a) and (b)], the detuning of the Mollow triplet side peaks matches the vibrational frequency $\pm\Omega$ of the molecule and thus coherent effects emerge due to the interaction between the electronic RF and the vibrational Raman scattering. The visibility and nature of these effects depends on the width of the lines, which is dominated by the Purcell effect and hence on the coupling $g$. We thus consider again the two representative situations of interaction analysed in this work: (i) $\hbar g=20$ meV, where the electronic peak is spectrally broader than the vibrational line, as the Purcell effect strongly broadens the former [Figs. 3 (a), (c)], and (ii) $\hbar g=9$ meV, a situation where the broadening of the electronic peaks is approximately equal to the vibrational broadening [Figs. 3 (b), (d)].

For clarity, the emission spectra for the selected values of $\mathcal{E}$ that provide the matching ($\hbar\mathcal{E}\approx 65$ meV for $\hbar g=20$ meV, and $\hbar\mathcal{E}\approx 130$ meV for $\hbar g\approx 9$ meV; roughly corresponding to a pumping power density of the order of $W\approx 1$ $\mathrm{m}\mathrm{W}\mathrm{/}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}$ and 10 $\mathrm{m}\mathrm{W}\mathrm{/}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}$, respectively) are shown in Fig. 3 (c) and (d), respectively. When the RF peak is much broader than the width of the Raman line [Fig. 3 (c)], the interference results in small but sharp features that might be detectable in experimental spectra and remind of Fano-resonancesMiroshnichenko et al. (2010). On the other hand, when the linewidth of the RF is similar to the linewidth of the Raman lines, the two spectral lines exhibit a clear anticrossing [inset in Fig. 3 (b)] that results in a splitting of the spectral features of each branch of the Mollow triplet [Fig. 3 (d)]. This splitting occurs as a result of the strong coupling between the molecule’s electronic (TLS) and its vibrational degrees of freedom Ramos et al. (2013), as predicted in the context of light emission from semiconductor quantum dots Kabuss et al. (2011).

The onset of strong coupling between the electronic and vibrational degrees of freedom, and thus the clear line splitting, can be understood with the help of a simplified Hamiltonian of the system (see Appendix C). This Hamiltonian is a result of an elimination of the plasmon cavity from the original Hamiltonian. Disregarding for the moment the vibrational part in Eq. (15), the simplified Hamiltonian ${H}_{\rm TLS}$ becomes

[TABLE]

Here $\mathcal{E}_{\rm pl}=-2g\alpha_{\rm S}\propto\mathcal{E}$ again corresponds to the amplitude of the plasmon-enhanced electric field. The Hamiltonian in Eq. (21) can be diagonalized by a unitary transformation that rotates the TLS Pauli matrices in the $x-z$ plane, (${\sigma}_{x},\,{\sigma}_{z}\rightarrow{\sigma}^{\prime}_{x},\,{\sigma}^{\prime}_{z}$):

[TABLE]

Under those operations, the simplified Hamiltonian, $H_{\rm TLS}=\hbar\lambda_{\rm TLS}\sigma^{\prime}_{z}+\hbar\frac{1}{2}{\delta}$, describes the dynamics of an effective electronic TLS dressed by the incident coherent illumination with an effective frequency $\lambda_{\rm TLS}=(\mathcal{E}_{\rm pl}^{2}+{\delta}^{2})^{1/2}$. According to this simplified Hamiltonian, the effective electronic frequency $\lambda_{\rm TLS}$ can be tuned by either changing the intensity of the incident laser (i.e. $\mathcal{E}_{\rm pl}\propto\mathcal{E}$) or by detuning the incident laser frequency ${\delta}$. This dressed TLS interacts with the molecular vibrations via the resonant Rabi interaction term Ramos et al. (2013)

[TABLE]

which becomes resonant if $\lambda_{\rm TLS}\approx\Omega$, the condition for the Mollow side peaks to coincide with the spectral position of the SERRS lines.

On top of the effect of dressing the molecular levels, the plasmonic cavity increases the effective damping rate $\Gamma_{\rm eff}$ of the TLS by means of the Purcell effect (causing the broad peaks of the Mollow triplet). The condition to reach strong coupling in this situation can be derived by relating the decay rate of the molecular vibration and that of the dressed electronic transition with the exciton-vibration coupling strength:

[TABLE]

This condition is at the origin of the strong coupling observed in the peaks of Fig. 3 (d) ($\hbar g=9$ meV), but it is not reached in the case presented in Fig. 3 (c) ($\hbar g=20$ meV) where Fano-type features appear as a sign of weak coupling. When the strong coupling between the vibrational Raman scattering and the RF pathways is reached, the peak-splitting in the emission spectra in Fig. 3 (d) can be also interpreted using the dressed-atom picture originally introduced by Cohen-Tannoudji Cohen-Tannoudji and Reynaud (1977a); Agarwal and Jha (1979). The dressed-atom picture allows for interpreting the splitting of the Raman and resonance-fluorescence peaks in terms of the coherent interaction among molecular vibronic states, induced by the incident coherent laser illumination. This approach shows that the final emission peaks emerge from a coherent combination of both the RF-type transitions and the Raman-type transitions, making the two mechanisms inseparably connected. We elaborate on the dressed-molecule picture assuming small $d$ in Appendix E, and after that provide a more general result allowing large $d$ in the next section.

V.2 Influence of large vibrational displacement $d$

We have so far used a moderate value of the displacement, $d=0.1$, however, in realistic molecules significant exciton-vibration coupling can lead to larger values of $d$. In Fig. 3 (e,f) we show the evolution of the spectrum with $d$ for the same two values of the plasmon-TLS coupling, i.e. $\hbar g=20$ meV [Fig. 3(e)] and $\hbar g\approx 9$ meV [Fig. 3(f)]. For all the values of $d$ considered, the laser intensity is chosen such that the resonance-fluorescence (RF) lines match the position of the Raman lines. For small values of $d\approx 0.1$, the RF profile follows the behaviour described in Figs. 3(c,d). As $d$ gradually increases, the spectra start to exhibit additional features due to the increasing importance of higher-order vibronic transitions. When the RF line is significantly broader than the Raman lines [Fig. 3(e)], an increase of the coupling $d$ gradually increases the spectral dip located at the frequency of the Stokes and anti-Stokes emission. Additional weak spectral features appear as $d$ is increased at larger laser detuning. When the linewidth of the Mollow side peaks is similar to the width of the Raman lines [Fig. 3(f)], the splitting of the strongly coupled hybrid lines becomes larger as $d$ increases. For large values of $d$, all of the spectra in Figs. 3(e,f) acquire a more complex structure due to the generally complicated coherent interaction between the molecular vibrational and the electronic degrees of freedom, with the emergence of the additional (weak) peaks originating from higher-order Raman and resonance-fluorescence transitions.

VI Vibrational pumping for strong laser intensities

In this section, we extend the treatment of linear-response SERRS introduced previously towards the case of strong incident illumination, where non-linear effects become important.

To that end we invoke the effective vibrational Hamiltonian introduced in Eq. (16) together with the incoherent damping $\Gamma_{-}=2\left(\Omega d\right)^{2}\mathrm{Re}\{\tilde{S}(\Omega)\}$ and pumping $\Gamma_{+}=2\left(\Omega d\right)^{2}\mathrm{Re}\{\tilde{S}(-\Omega)\}$ rates in Eq. (17) and Eq. (18), respectively. As above, the spectral function $\tilde{S}(s)$ is obtained from the effective dynamics of the TLS, which is effectively broadened by the plasmon via the Purcell effect. The hierarchy of approximations considered in this section is schematically depicted in Fig. 4 (a). These effective rates are dependent on the spectral function $\tilde{S}(s)$ of the reservoir evaluated at frequencies $\Omega$ and $-\Omega$, respectively. Note that the analytical model is limited to cases where the electron-vibration coupling $\Omega d$ is smaller than the effective broadening $\Gamma_{\rm eff}$ of the electronic resonance. We thus perform full numerical calculations to obtain the results (i.e., populations) spanning the full range of model parameters, however we use the analytical model for qualitative discussion.

The value of the spectral function $\tilde{S}(s)$ at frequencies $\pm\Omega$ determines the strength of the effective vibrational pumping ($\Gamma_{+}$) or damping ($\Gamma_{-}$), which modify the vibrational populations as can be seen from the first expression in Eq. (19). It is therefore possible to achieve different regimes of interaction with the vibrations which range from pumping to damping by simply modifying the illumination conditions (laser intensity and frequency detuning) that provoke a variation of the shape of the spectral function. When the laser intensity is large, the reservoir function $\tilde{S}(s)$ reflects the structure of the TLS dressed by the incident laser, and therefore it becomes qualitatively different from the weak-illumination case.

In Fig. 4(b) the spectral function, $\mathrm{Re}\{\tilde{S}(s)\}$, for $d=0.2$, $\hbar\mathcal{E}=100$ meV, and $\hbar g=20$ meV, is shown to peak around the effective frequencies of the dressed TLS ($s=\pm\lambda_{\rm TLS}$). When the incident laser is detuned from the TLS transition ($\tilde{\delta}=\delta+d^{2}\Omega\neq 0$) the spectral function changes symmetry. For $\tilde{\delta}>0$ (red detuning marked with a red line) a regime of vibrational damping can be reached ($\mathrm{Re}\{\tilde{S}(\Omega)\}>\mathrm{Re}\{\tilde{S}(-\Omega)\}$), whereas for $\tilde{\delta}<0$ (blue detuning marked with a blue line) a regime of vibrational pumping ($\mathrm{Re}\{\tilde{S}(\Omega)\}<\mathrm{Re}\{\tilde{S}(-\Omega)\}$) is achieved. This effect is more pronounced for a situation where $\tilde{S}(\pm\Omega)$ corresponds to the maxima of $\tilde{S}(s)$.

To illustrate the possibility to achieve a controlled excitation of molecular vibrations on demand, we numerically solve the full Hamiltonian of the system [Eq. (5)], and show in Figs. 4(c,d) the steady-state vibrational population $\langle b^{\dagger}b\rangle$, for an electron-plasmon coupling of $\hbar g\approx 9$ meV and two different values of the dimensionless displacement ($d=0.2$ and $d=0.5$). The vibrational pumping is non-trivially influenced by both the detuning ${\delta}$ of the incident laser frequency from the TLS transition frequency and by the incident laser intensity ($\propto\mathcal{E}^{2}$), so that the optimal laser intensity depends on the laser detuning ${\delta}$. When the electron-vibration coupling is large ($d=0.5$), the population reaches multiple intense local maxima. In this case, by adequately tuning the laser frequency one can efficiently excite Franck-Condon transitions involving a change of more than one vibrational transition (higher-order processes). As expected, the population maxima are found when the spectral position of the side peaks of the electronic spectral function matches the frequency of the higher-order vibrational transitions ($\lambda_{\rm TLS}\approx n\Omega$, with $n$ an integer), a condition traced by the blue lines in Fig. 4 (d), and displayed only for negative detuning ${\delta}$.

VII Selective vibrational pumping

The potential to control the activation of molecular vibrations can be exploited in the selective excitation of different vibrational modes. Let us consider the coupling of a plasmonic system with a molecule supporting two vibrations at frequencies $\hbar\Omega_{1}=10$ meV and $\hbar\Omega_{2}=17.5$ meV, both coupled to independent reservoir modes (baths) with $\hbar\gamma_{\rm vib,\,1}=\hbar\gamma_{\rm vib,\,2}=1$ meV, as schematically depicted in Fig. 5 (a). We simplify the description of the system and use the effective Hamiltonian where the plasmonic degrees of freedom are eliminated (see Appendix F). We assume that the vibrational modes are coupled to the TLS via a polaronic coupling term ($d_{1}=d_{2}=0.2$), and do not consider the direct coupling between the two vibrational modes. However, this model Hamiltonian naturally couples the two vibrational modes indirectly via the electronic TLS of the molecule. Our model thus partially accounts for thermalization effects, without considering the effect of the surrounding environment that may further incoherently couple the vibrational modes.

The resulting vibrational populations, $\langle b_{i}^{\dagger}b_{i}\rangle$, are shown in Fig. 5 (b) as a function of the intensity and detuning of the incoming laser. The color map depicting the population of the vibrational mode at frequency $\Omega_{1}$ is displayed together with a dashed contour plot that shows the corresponding results for the mode at $\Omega_{2}$. Each mode presents a clear maximum for suitable illumination conditions. Noticeably, the maxima are shifted with respect to each other both in frequency and amplitude, so that changing the illumination conditions serves to pump more efficiently one mode or another. To highlight the selectivity of the vibrational pumping mechanism, we extract line cuts of Fig. 5 (b) for constant laser pumping, $\hbar\mathcal{E}=100$ meV, [Fig. 5 (c)], and for constant laser detuning $\hbar\delta=-9$ meV [Fig. 5 (d)]. As observed in Fig. 5 (c,d) the conditions of intensity and detuning for maximum population of one mode give a much weaker population of the other mode (solid versus dashed lines). This scheme of interactions makes it possible to achieve selective vibrational pumping by either tuning the laser frequency for a given illumination intensity or modifying the laser intensity for a fixed illumination frequency.

VIII Effects of temperature on vibrational pumping

Practical SE(R)RS experiments are usually performed at a finite temperature, a situation where the thermal populations of the molecular vibrations can become considerable. We thus study this effect in this section. We calculate the steady-state vibrational populations for a temperature of $T=77$ K (liquid nitrogen temperature) and of $T=300$ K (room temperature) for different detuning and intensity of the incident laser for the same range of parameters as in Fig. 4 (c,d). We assume that the electronic states of the molecule interact with a single vibrational mode. The results of the vibrational population as a function of laser detuning and intensity are shown in Fig. 6 as two-dimensional color maps.

For the liquid-nitrogen temperature ($T=77$ K) the maps of the vibrational population shown in Fig. 6 (a,b) for $d=0.2$ and $d=0.5$, respectively, are qualitatively similar to the results obtained when considering $T=0$ K. However, we observe that the total vibrational populations reach larger maximal values due to the thermal population, which for a vibrational excitation of energy $10$ meV reaches at $T=77$ K a value of $n_{77{\rm\,K}}\approx 0.3$. Interestingly, for positive values of detuning $\delta$ we observe vibrational populations smaller than $n_{77{\rm\,K}}$ which is a result of the optomechanical cooling process that is not observable for $T=0$ K. For room temperature ($T=300$ K) the thermal population of the vibrational mode of energy $10$ meV reaches $n_{300{\rm\,K}}\approx 2$, which is larger than the optomechanically induced steady-state population reached when thermal effects are disregarded. The map of vibrational populations shown in Fig. 6 (c) for small $d=0.2$ and $T=300$ K shows, as in the case of $T=77$ K, that the optomechanical process can lead either to an increase or decrease of the steady-state vibrational population with respect to the thermal equilibrium state. When the detuning fulfils the condition of optomechanical vibrational pumping (i.e. $\delta\approx-\lambda_{\rm TLS}$) the population can reach values larger than the thermal population, as in this case the optomechanical pumping process enhances the effect of thermal population pumping. Conversely, when the optomechanical damping process is enhanced ($\delta\approx\lambda_{\rm TLS}$) the steady-state vibrational population features cooling. These effects are more pronounced in Fig. 6 (d) where we consider $T=300$ K and a stronger electron-vibration coupling of $d=0.5$ [notice the different scale of the color map in Figs. 6 (a-d)]. As pointed out, the steady-state vibrational population is pumped above the value determined by the thermal equilibrium for $\delta<0$. This effect is resonantly enhanced whenever $\delta\approx-n\lambda_{\rm TLS}$ ($n$ integer). On the other hand, the system is cooled for the opposite detuning ($\delta\approx n\lambda_{\rm TLS}$).

We thus conclude that thermal effects become important if the thermal vibrational population becomes comparable to the steady-state population induced by the optomechanical pumping or damping processes. Furthermore, when a finite temperature is considered, it is possible to observe optomechanical vibrational cooling.

IX Conclusion

In conclusion, we have used the formalism of quantum cavity electrodynamics to discuss SERRS as a quantum optomechanical process. First, we have shown the similarities and fundamental differences between the resonant and off-resonant SERS processes. In the linear response regime, the electronic transition appearing in the SERRS plays a role analogous to that of the plasmonic cavity in the off-resonant SERS, however, the former offers generally larger values of optomechanical coupling and smaller values of the molecular-exciton damping $\gamma_{\sigma}+\Gamma_{\rm eff}$, as compared to the plasmon broadening $\gamma_{a}$. We have further described vibrational pumping in SERRS and identified its two principal mechanisms: the efficient pumping of the molecule by laser resonant with the molecular exciton (optimized when $\hbar\delta\approx 0$ eV) and the resonant enhancement of the Raman-Stokes emission ($\delta\approx-\Omega$).

For strong laser intensities, the non-linear character of the molecular electronic TLS introduces non-trivial effects associated with the dressing of the molecular levels by the intense laser. We have shown that a strong exciton-vibration interaction appears when a side peak of the RF (Mollow triplet) of a single molecule Wrigge et al. (2008) is tuned to match the frequency of vibrational Raman lines. The fingerprint of this interaction is accessible in optical emission spectra through the presence of interference features or peak splitting. Finally, the regime of vibrational pumping achieved in SERRS can exploit the spectral features of the molecular electronic resonance, which are often narrower than the plasmonic resonances used in off-resonant SERS. This opens the possibility for selective pumping of different vibrational modes of the molecule. Selective vibrational pumping can offer new ways to control chemical reactivity of molecules Crim (1996) or to engineer vibrational quantum states of experimental interest.

Acknowledgements

The authors acknowledge projects FIS2016-80174-P from Spanish MINECO, ELKARTEK KK-2018/00001, H2020-FETOPEN project ”THOR” Nr. 829067 of the European Commission, PI-2017-30 and PI-2016-41 of the Departamento de Educación, Política Lingüística y Cultura of the Basque government.

Appendix A System Hamiltonian

In this appendix we describe the Hamiltonian used in the main text to describe the resonant surface-enhanced Raman scattering (SERS) process [Eq. (6)], where we consider the coupling between a plasmonic mode and the electronic and vibrational levels of a molecule. When a monochromatic laser excites the plasmonic mode coupled to the molecule, the system Hamiltonian can be written in the form:

[TABLE]

with

[TABLE]

Here $\sigma$ ($\sigma^{\dagger}$) are the lowering (raising) operators of the two level system (TLS) representing the electronic structure of the molecule (with $\sigma^{\dagger}\sigma=\sigma_{\rm e}$), $b$ ($b^{\dagger}$) are the annihilation (creation) operators of the vibrational mode and $a$ ($a^{\dagger}$) are the annihilation (creation) operators of the single bosonic mode representing the plasmonic cavity. The constants appearing in the Hamiltonian have the following meaning: $g$ represents the coupling between the plasmon and the electronic TLS, $\omega_{\rm pl}$ is the plasmon frequency, $E_{\rm g}$ and $E_{\rm e}$ are the energies of the ground and excited electronic state of the molecule, respectively, $\Omega$ is the vibrational frequency and $d$ introduces the electron-phonon coupling. The monochromatic coherent laser illumination of frequency $\omega_{\rm L}$ is coupled to the plasmonic mode via the constant $\mathcal{E}$ (that we take as real), which is proportional to the electric field amplitude of the incident light.

The system Hamiltonian is formally split into the Hamiltonian describing the plasmonic mode as a bosonic oscillator $H_{\rm pl}$, the term describing the level structure of the bare molecule $H_{\rm mol}$, the coupling between the plasmon and the molecule $H_{\rm pl-e}$ and the pumping of the plasmon by the classical laser field $H_{\rm pump}$ in the rotating wave approximation (RWA), assuming that the pumping amplitude $\mathcal{E}$ is sufficiently small compared to the plasmon frequency ($\mathcal{E}\lesssim 0.1\omega_{\rm pl}$).

The molecular Hamiltonian $H_{\rm mol}$ contains the energy splitting of the molecule’s electronic levels, $\hbar[E_{\rm e}-E_{\rm g}]\sigma{{}^{\dagger}}\sigma$, and the vibrational term which depends on the electronic state, $\hbar\Omega(b+\sigma_{\rm e}d)(b^{\dagger}+\sigma_{\rm e}d)$. This coupling is obtained from the Born-Oppenheimer approximation, that takes into account that the vibrational states of the ground states have different equilibrium position than the vibrations defined on the excited state potential energy surface. Due to this displacement, $d$, the vibrational eigenstates in the ground state are not orthogonal to the ones in the excited state. Therefore, when the molecule is excited from the ground electronic state to the excited electronic state, it also simultaneously changes the vibrational state (according to the Franck-Condon principle).

The coupling between the molecule and the plasmon is expressed in the RWA as $H_{\rm pl-e}=\hbar ga\sigma{{}^{\dagger}}+\hbar g^{\ast}a{{}^{\dagger}}\sigma$. The RWA is justified in situations where $g<<\omega_{\rm pl}\approx E_{\rm e}-E_{\rm g}$. Importantly, the RWA allows for further simplifying transformations of the Hamiltonian.

To introduce incoherent effects we employ the approach based on the solution of the quantum master equation for the density matrix $\rho$:

[TABLE]

where the bracket symbolizes the commutator and the Lindblad terms $\mathcal{L}_{c}[\rho]$ introduce the incoherent damping. In particular we use the following Lindblad superoperators:

[TABLE]

where $\gamma_{\sigma}$ is the electronic decay rate, $\gamma_{b}$ the vibrational decay rate, and $\gamma_{a}$ the plasmonic decay rate.

To facilitate the numerical calculation, we apply a unitary transformation $\tilde{H}_{\rm om}^{\rm res}=U_{\omega_{\rm L}}H_{\rm om}^{\rm res}U_{\omega_{\rm L}}^{\dagger}-{\rm i}\hbar U_{\omega_{\rm L}}\dot{U}_{\omega_{\rm L}}^{\dagger}$ with $U_{\omega_{\rm L}}=\exp({\rm i}\sigma_{\rm e}\omega_{\rm L}t+{\rm i}a^{\dagger}a\omega_{\rm L}t)$ that transforms the Hamiltonian in Eq. (24) into the rotating frame (interaction picture). As a consequence, the Lindblad terms remain unchanged, but the Hamiltonian in Eq. (24) is modified by simply replacing (for brevity we keep the same notation for the transformed operators as for the original operators throughout the text):

[TABLE]

The resulting Hamiltonian is time independent which facilitates the numerical solution. However, it includes the direct pumping of the plasmon mode. In the numerical implementation, where we represent the plasmonic states by the number states of the plasmon Hamiltonian $H_{\rm pl}$, we would need a large number of plasmon states to correctly describe the excitation by a strong incident laser. We avoid this problem by redefining the plasmon creation and annihilation operators:

[TABLE]

where $\alpha_{\rm S}=\frac{-\mathcal{E}}{\Delta-{\rm i}\frac{\gamma_{a}}{2}}$. This particular choice of $\alpha_{\rm S}$ allows us to define a new Hamiltonian $H_{\rm om}^{{\rm res},{\alpha}}$ [see also Eq. (1) of the main text]:

[TABLE]

with

[TABLE]

The Lindblad terms appear unchanged provided that the plasmon operators are expressed in the shifted basis. The final form of the Master equation is thus:

[TABLE]

where we have to keep in mind that we are working in the interaction picture and in the displaced basis when we evaluate the correlation functions and the operator mean values.

Appendix B Calculation of emission spectra

We calculate the emission spectra numerically from the quantum regression theorem (QRT) Breuer and Petruccione (2003) using the expression

[TABLE]

In order to solve the dynamics of the damped system we need to solve the quantum master equation [Eq. (34)] for the density matrix $\rho$ given by the Hamiltonian and by the Lindblad terms. To that end, we rewrite the quantum master equation in a form where the density matrix $\rho$ appears as a column vector (see e.g. Ref. Am-Shallem et al. (2015)):

[TABLE]

In the equations, there often appear expressions where the operators act on the density matrix from the right or from the left (e.g. the term $a\rho a^{\dagger}$). In the technical implementation, the expressions are transformed as:

[TABLE]

Here $\otimes$ represents the Kronecker product, T denotes transposition and $O_{1}$, $O_{2}$ are generic operators. In practise, if the dimension of the truncated Hilbert space is set to $N=N_{\rm vib}\times N_{\rm pl}\times 2$ (with $N_{\rm vib}$ [$N_{\rm pl}$] the maximal number of vibrational [plasmon] number states considered in the calculation), the vectorised density matrix will have a length of $N^{2}$ and the matrix $\left(O_{2}^{\rm T}\otimes O_{1}\right)$ will be of dimension $N^{2}\times N^{2}$.

Equation (34) becomes in this representation:

[TABLE]

where the $N^{2}\times N^{2}$ square matrix $\mathscr{L}(t)$ represents the Liouville superoperator. In general $\mathscr{L}(t)$ can depend on time, but in our model it is time independent. The steady-state density matrix $\rho_{\rm ss}$ is then obtained from the quantum master equation [Eq. 38] as the eigenvector belonging to the zero eigenvalue of the matrix $\mathscr{L}$.

The two time correlation function [Eq. (35)] that defines the spectrum is calculated using the quantum regression theorem. Utilizing Laplace transform techniques and with the substitution of the Laplace parameter $s\rightarrow-{\rm i}\omega$:

[TABLE]

The direct implementation of Eq. (39) requires to invert the Liouvillian matrix $\mathscr{L}$ [time independent in the frame rotating with $\exp(-{\rm i}\omega_{\rm L}t)$] for each frequency of interest. However, this procedure becomes inefficient for calculations of emission spectra. In such a case, and assuming that the Liouvillian is represented by a diagonalizable matrix (which we verify numerically), it is often more convenient to expand the time-dependent solution into an exponential series with the exponents being the eigenvalues of the Liouvillian superoperator Tan (1999).

The time-dependent solution of the density matrix [or of the vector $\{(a{{}^{\dagger}})^{\rm T}\otimes I\}\vec{\rho}$ since according to the QRT they obey the same differential equation, Eq. (34)] is formally given by the exponential of the Liouvillian as:

[TABLE]

Equivalently, the exponential can be expressed using the eigenvalue decomposition of the Liouvillian $\mathscr{L}={S}{D}{S}^{-1}$ as follows:

[TABLE]

where the operator $D$ is represented by a diagonal matrix so that its exponentiation simply indicates exponentiation of the matrix diagonal elements one by one. We can then write the correlation function $\langle a{{}^{\dagger}}(0)a(\tau)\rangle$ as

[TABLE]

where, for convenience, we have defined a full matrix $A$, a diagonal matrix $B$ containing the exponentiated eigenvalues, and a vector $v$. Note that the explicit dependence on the laser frequency, $\omega_{\rm L}$, appears because we express the operators in the interaction picture. The trace operator is defined in the original representation where the density operator has the form of a square matrix (the left hand side in Eq. (36)). The matrix $B$ can be seen element-wise as $B_{ii}=\delta_{ij}\exp(d_{i}t)$, where $d_{i}$ are the diagonal elements of the matrix $D$. The vector $c=ABv$ is therefore equal to the following exponential series:

[TABLE]

The trace can be represented as a scalar product of the vectorised operator with vector ${\nu}_{tr}$ where $(\nu_{tr})_{j}=1$ for $j=1,\,n+2,\,2n+3,\,...,\,n^{2}$ ,and 0 otherwise. The trace in Eq. (42) above becomes:

[TABLE]

with

[TABLE]

Last, inserting this expression into Eq. (35) we get:

[TABLE]

Appendix C Effective two-level-system Hamiltonian under coherent laser illumination

We discuss in the following how to eliminate the plasmon cavity to obtain a new effective Hamiltonian of the molecule. We assume that the plasmonic cavity (after transforming out the coherent displacement $\alpha_{\rm S}$) acts as a fluctuating reservoir that effectively damps the molecule via the Purcell effect. To describe the effects of the reservoir, we use the standard methods of the quantum noise approach Breuer and Petruccione (2003) and eliminate the plasmon under the assumption that the broadening of the plasmon excitation is considerably larger than the coupling strength between the plasmon and the electronic TLS. We obtain the electronic decay of the molecule into the plasmonic mode (the Purcell effect), assuming that the plasmon cavity is unpopulated (after removing the coherent contributions as described in Appendix A), expressed by the following Lindblad term:

[TABLE]

with

[TABLE]

and where

[TABLE]

which, assuming that the plasmon obeys a dynamics unperturbed by the presence of the molecule, yields

[TABLE]

We neglect the slight frequency shift of the TLS due to the action of the cavity which is formally given by the imaginary part of $S_{a}(s)$. The Hamiltonian of the reduced system thus has the form [Eq. (15)]:

[TABLE]

where we defined $\mathcal{E}_{\rm pl}=-2g\alpha_{\rm S}$ and $\sigma_{x}$ is the Pauli $x$ operator.

The Hamiltonian in Eq. (48) [Eq. (15)], from which the plasmon has been eliminated, can be recast into a form where an effective electronic TLS dressed by the incident coherent illumination interacts with the molecular vibrations via the Rabi interaction term Ramos et al. (2013). We first start by grouping the terms that correspond to the TLS Hamiltonian under strong laser illumination. In particular, we consider as the TLS Hamiltonian the terms that are free of the vibrational operators (do not contain the electron-vibration coupling, i.e. $d=0$):

[TABLE]

where we used the fact that the operator $\sigma_{\rm e}=\sigma^{\dagger}\sigma$ can be rewritten with the help of the standard Pauli $z$ operator as $\sigma_{\rm e}=\frac{1}{2}(\sigma_{z}+I)$.

We first apply the following rotation in the space of the standard Pauli operators $\sigma_{x}$, $\sigma_{y}$ and $\sigma_{z}$ (with the primed operators being the new Pauli operators, and $\lambda_{\rm TLS}=\sqrt{{\delta}^{2}+{\cal E}_{\rm pl}^{2}}$) that diagonalizes Eq. (49):

[TABLE]

The resulting Hamiltonian for the molecule after this transformation is:

[TABLE]

where the TLS Hamiltonian from Eq. (49) corresponds to the diagonal term $\hbar\frac{1}{2}\lambda_{\rm TLS}\sigma^{\prime}_{z}$.

For completeness we note that $\sigma_{y}$ has not been changed by the transformation and from Eq. (50) it follows that the original lowering operator [$\sigma=\frac{1}{2}(\sigma_{x}-{\rm i}\sigma_{y})$] transforms into the following form:

[TABLE]

Appendix D Vibrational pumping in the quantum-noise approximation

The Hamiltonian in Eq. (48) describes the dynamics of the molecule after the plasmon is eliminated. We describe here how it is also possible to eliminate the TLS degrees of freedom to focus on the dynamics of the vibrations, under the assumption that the plasmon-enhanced decay rate of the TLS ($\Gamma_{\rm eff}+\gamma_{\sigma}$) is much larger than the decay rate of the vibrations $\gamma_{b}$ and that the electron-phonon coupling is weak (for example, $d=0.1$ and $\hbar g=50$ meV). We therefore divide the Hamiltonian in Eq. (21) into the part representing the system, the reservoir and the system-reservoir interaction as follows:

[TABLE]

where $\langle\sigma_{e}\rangle$ is the steady-state average of the TLS excited-state population calculated for the TLS decoupled from the vibrations, and $\tilde{\delta}=\delta+d^{2}\Omega$. The elimination of the reservoir can be performed using the quantum noise approach Breuer and Petruccione (2003); Jaehne et al. (2008); Rabl (2010) to Eq. (48), giving the Hamiltonian:

[TABLE]

The effective damping of the vibrations (appearing aside of the intrinsic damping $\gamma_{b}$) can be expressed via the Lindblad superoperator [Eqs. (17) and (18)]

[TABLE]

and the effective vibrational pumping via the superoperator

[TABLE]

The damping and pumping rates that appear in the Lindblad superoperators are given by

[TABLE]

and

[TABLE]

respectively. Here $\tilde{S}(s)=\int_{0}^{\infty}\langle\langle\sigma_{\rm e}(\tau)\sigma_{\rm e}(0)\rangle\rangle e^{{\rm i}s\tau}{\rm d}\tau$ (with $\sigma^{\dagger}\sigma=\sigma_{\rm e}$) is the spectral function that comprises the properties of the electronic TLS ”bath” uncoupled from the vibrations and broadened by the plasmon.

We show below that the analytical expression for the spectral function $\mathrm{Re}\{\tilde{S}(s)\}$ can be expressed in the form:

[TABLE]

Here $\Gamma_{\rm tot}=\Gamma_{\rm eff}+\gamma_{\sigma}$ [see Eq. (47) for definition of $\Gamma_{\rm eff}$] and $\langle\langle\sigma_{\rm e}(\tau)\sigma_{\rm e}(0)\rangle\rangle=\langle\sigma_{\rm e}(\tau)\sigma_{\rm e}(0)\rangle-\langle\sigma_{\rm e}\rangle\langle\sigma_{\rm e}\rangle$ is the part of the correlation function that corresponds to the fluctuations of the operators around the steady-state value. We generally define the fluctuating correlation function as $\langle\langle O_{1}(\tau)O_{2}(0)\rangle\rangle=\langle O_{1}(\tau)O_{2}(0)\rangle-\langle O_{1}\rangle\langle O_{2}\rangle$.

We obtain the spectral function $\mathrm{Re}\{\tilde{S}(s)\}$ in Eq. (55) following the procedure described in references Jaehne et al. (2008); Rabl (2010). We start with the Liouville equation describing the TLS dynamics after effective elimination of the plasmonic cavity. In the case of no electron-phonon coupling, it contains the Hamiltonian:

[TABLE]

together with the Lindblad term

[TABLE]

We obtain the following equations of motion for the mean values of the operators:

[TABLE]

According to the quantum regression theorem Breuer and Petruccione (2003), the correlation functions $\langle\langle\sigma(\tau)\sigma_{\rm e}(0)\rangle\rangle$, $\langle\langle\sigma^{\dagger}(\tau)\sigma_{\rm e}(0)\rangle\rangle$ and $\langle\langle\sigma_{\rm e}(\tau)\sigma_{\rm e}(0)\rangle\rangle$ obey the same time evolution as $\langle\sigma(\tau)\rangle$, $\langle\sigma^{\dagger}(\tau)\rangle$ and $\langle\sigma_{\rm e}(\tau)\rangle$, respectively:

[TABLE]

with the initial values:

[TABLE]

Here all the mean values are evaluated in the steady state. The inhomogeneous part in Eq. (78) disappears because of the conveniently chosen value of the correlation functions when $\tau\rightarrow\infty$.

The direct solution of Eq. (78) yields the result in Eq. (55) for $\mathrm{Re}\{\tilde{S}(s)\}$, which in turn provides the analytical expressions for the effective vibrational pumping $\Gamma_{+}$ and damping $\Gamma_{-}$. We can then solve the effective vibrational dynamics:

[TABLE]

from which we obtain the steady-state values:

[TABLE]

Eq. (82) shows two different sources of vibrational pumping that are involved in the process. The first one, represented by the term $|\langle b\rangle|^{2}$ is the coherent pumping of the vibrations due to the TLS. The second term, $\langle b^{\dagger}b\rangle_{\rm SS,in}=\frac{\Gamma_{+}}{\gamma_{b}+\Gamma_{-}-\Gamma_{+}}$, is due to the ”incoherent” vibrational pumping induced by the fluctuating part of the TLS quantum dynamics.

Under the considered conditions, the incoherent pumping term dominates over the weak coherent pumping term so that the shape of the correlation function $\mathrm{Re}\{\tilde{S}(s)\}$ governs the population of the vibrations as indicated by Eqs. (53) and (54). As we have demonstrated in Section VIII, cooling is also possible Jaehne et al. (2008); Rabl (2010).

Appendix E Small vibrational displacement $d$: dressed-molecule picture

The regime where the linewidth of the RF peaks is comparable to the width of the Raman peaks is a limiting case of Raman scattering in intense fields that has been studied in the context of atomic physics Cohen-Tannoudji and Reynaud (1994, 1977b); Agarwal and Jha (1979); Cohen-Tannoudji and Reynaud (1977a). To understand the splitting of the lines that appear when the Mollow triplet side peaks have the frequency of the Raman lines, it is useful to rewrite the Hamiltonian into a form where the coupling among vibrational states is explicitly present.

This Hamiltonian can be derived from the reduced Hamiltonian ${H}_{\rm red}$ appearing in Eq. (15) by applying the so called small polaron transformation ${H}_{\rm red}\rightarrow{{H}_{\rm red}^{{}^{\prime}}}={U}_{\sigma_{\rm e}}{H}_{\rm red}{U}_{\sigma_{\rm e}}^{\dagger}$, which is represented by the unitary matrix in the form of a displacement operator ${U}_{\sigma_{\rm e}}=\exp\left[d{\sigma}_{\rm e}({b}^{\dagger}-{b})\right]$. This transformation has two effects on the Hamiltonian. First, the vibrational term in ${H}_{\rm red}$ transforms as:

[TABLE]

and, second, the pumping term of the TLS acquires an additional factor that includes the vibrational operators (which yield the well known Franck-Condon factors):

[TABLE]

In this approach we assume that the influence of the plasmon (given by ${H}_{\rm pl}$ and ${H}_{\rm pl-e}$ in Eq. (5)) is effectively included as an enhancement of the incident laser field $\mathcal{E}_{\rm pl}$. Furthermore we consider only the case where the splitting of the fluorescence and Raman spectral peaks [see e.g. Fig. 3 (b,d,f)] is larger than their broadening (the so called secular limit). In such a case, the incoherent broadening does not influence the peaks positions and, therefore, in the following we consider that the system can be described only by the the simplified Hermitian Hamiltonian (we do not consider the Lindblad terms as we are mainly interested in the nature of the transitions). We further assume weak electron-phonon coupling in the molecule and expand the exponential terms containing the vibrational operators to the first order: $\exp\left[\pm d({b}^{\dagger}-{b})\right]\approx I\pm d({b}^{\dagger}-{b})$. Last, we reduce the system comprising the vibrations and the TLS into an effective four-level system that consists of the ground and excited electronic states considering zero or one vibrational excitation for each electronic state.

The diagram of the resulting effective system is drawn in Fig. 7 (a) and (b). By diagonalizing this $4\times 4$ Hamiltonian we achieve a new level structure of the system that, in the dressed-molecule picture, provides the positions of the emission peaks (for detailed discussion of the dressed-molecule (atom) picture see e.g. chapter 10 of reference Scully and Zubairy (1997)). Below we briefly describe how the emission spectra can be understood using this dressed-molecule picture.

In the dressed-molecule picture we consider the simplified Hamiltonian, which can be formally defined in the basis of states $[|{\rm g},\mathcal{N},0\rangle,|{\rm e},\mathcal{N}-1,0\rangle,|{\rm g},\mathcal{N},1\rangle,|{\rm e},\mathcal{N}-1,1\rangle]$, with e (g) labelling the electronic excited (ground) state, $\mathcal{N}$ labelling the ”photon” number state of the exciting field and 0 (1) labelling the number of vibrational excitations, respectively. The exciting field is not quantized explicitly in the original Hamiltonian [Eq. (5)] where it is represented by the plasmon coherent-state amplitude $\alpha_{\rm S}$. We therefore assume that the exciting field is a highly populated bosonic field which peaks sharply around a (mean) occupation number $\tilde{\mathcal{N}}$ yielding $\alpha_{\rm S}=\sqrt{\tilde{\mathcal{N}}}g_{\rm PL-L}$, with $g_{\rm PL-L}$ formally defined as a small coupling constant (such that $\tilde{\mathcal{N}}\gg 1$) between the equivalent exciting field and the molecule. Note that the formal definition of the exciting field is not important for the following discussion as by introducing the quantized exciting field we only aim at mimicking the action of the semi-classical pumping term. However, the number states $|\mathcal{N}\rangle$ of the exciting field are convenient to discuss the dressing of the molecular excited states in terms of the hybridization of the quantum-mechanical states.

We further define the total number of excitations as $n=\mathcal{N}+\delta_{i{\rm e}}$, with $i={\rm e,\,g}$ and $\delta_{ij}$ the Kronecker delta. We consider that the electronic levels, carrying the fine vibrational structure, are dressed by the strong laser illumination. The Hamiltonian can be expressed in the interaction picture of the incident laser field which is exactly tuned to the electronic transition, $\hbar\delta=0$ eV. In the basis $[|{\rm g},\mathcal{N},0\rangle,|{\rm e},\mathcal{N}-1,0\rangle,|{\rm g},\mathcal{N},1\rangle,|{\rm e},\mathcal{N}-1,1\rangle]$ the Hamiltonian, ${H}_{\rm red}^{{}^{\prime}}$, can be represented by the matrix:

[TABLE]

For vanishing electron-phonon coupling, $d=0$, the Hamiltonian in Eq. (89) reduces to the form describing a pair of TLSs dressed by the incident laser illumination. The process of dressing (i.e. diagonalization of the above Hamiltonian with $d=0$) can be viewed as a mixing of the electronic states with the high number states of the exciting laser field, giving rise to the basis of hybridized states $[|n_{-},0\rangle,|n_{+},0\rangle,|n_{-},1\rangle,|n_{+},1\rangle]$ where the first quantum number, $n$, labels the total number of electronic plus laser excitations, and the second quantum number, $m$, belongs to the vibrational states [see Fig. 8 (a) for schematics of the corresponding energy levels]. The hybrid states are defined as $|n_{\pm},m\rangle\equiv(|{\rm g},\mathcal{N},m\rangle\pm|{\rm e},\mathcal{N}-1,m\rangle)/\sqrt{2}$ with $+$ labelling the state with higher energy. In the new basis of such dressed states, we can represent the Hamiltonian as:

[TABLE]

For $d=0$, the splitting of the states $|n_{\pm}\rangle$ for the TLS in each vibrational Fock state is $|2g\alpha_{\rm S}|=|\mathcal{E}_{\rm pl}|$. The two dressed TLS defined for each vibrational Fock state are mutually shifted by the vibrational frequency $\Omega$ along the energy axis. In the absence of electron-phonon coupling $d$, the RF emission (dominating in this case the inelastic emission) is given purely by the transitions conserving the vibrational number state and changing the total number of excitations, $n$, by one. In particular, the central Mollow peak is given by transitions between $|(n+1)_{+},0(1)\rangle\rightarrow|n_{+},0(1)\rangle$ and $|(n+1)_{-},0(1)\rangle\rightarrow|n_{-},0(1)\rangle$, while the side peaks contain transitions $|(n+1)_{-},0(1)\rangle\rightarrow|n_{+},0(1)\rangle$ (red detuned) and $|(n+1)_{+},0(1)\rangle\rightarrow|n_{-},0(1)\rangle$ (blue detuned), respectively. The respective transitions and their corresponding emission peaks (the Mollow triplet) are schematically marked in Fig. 8 (a), where the colouring of the spectral emission peaks (bottom) corresponds to the colour of the respective arrows marking the transitions (top).

If we switch on the electron-phonon interaction $d$, a mixing between the levels belonging to the two vibrational Fock states is introduced, simultaneously allowing additional transitions yielding the Raman emission (i.e. changing the vibrational Fock state). The details of the level mixing and the subsequent emission spectra depend on the particular choice of pumping strength, $\mathcal{E}_{\rm pl}$, in combination with the value of the electron-phonon coupling, $d$. In the following we consider a particular case where the Mollow triplet side peaks overlap with the Raman lines with the laser frequency exactly tuned to the TLS energy splitting ($\hbar\delta=0$ eV and $|\mathcal{E}_{\rm pl}|=\Omega$). Upon diagonalization, the Hamiltonian in Eq. (89) [Eq. (99)] yields the following spectrum of energy levels:

[TABLE]

where we used the assumption that $d\ll 1$ to perform the Taylor expansion of the square root up to the first order. The states having energy $\lambda_{\pm}$ are a coherent admixture of states containing zero and one vibrational excitation $|n_{-},1\rangle$ and $|n_{+},0\rangle$, as discussed above, and the states of energy $\lambda_{3,4}$ can be identified (up to small $o(d)$ admixtures of other states) with $|n_{-},0\rangle$ and $|n_{+},1\rangle$ whose energy is renormalized due to the off-resonant electron-phonon coupling. This level structure of the molecule does not explicitly contain the quantized electromagnetic field of the incident laser. However, in the dressed-molecule picture the molecular level structure [Eq. (100)] is periodically repeated for each manifold represented by a specific number of excitations, $n$, and thus appears repeated along the energy axis displaced by integer values of the laser frequency $\omega_{\rm L}$, as schematically illustrated in Fig. 8 (b). In this picture, the emission events are represented by transitions between manifolds that differ by one excitation quantum of the electronic and effective photonic states, i.e. transitions between the manifolds containing $n$ and $n+1$ excitation quanta [represented by colored lines in Fig. 8 (b)].

The dressed-molecule picture above nicely allows us to identify the spectral peaks which appear in the complex photon emission spectra obtained from the numerical calculation of the complete Hamiltonian in Eq. (5) with its corresponding Lindblad terms. Figure 8 (c) shows such a situation where two vibrational levels corresponding to the ground vibrational state and the first excited vibrational state are considered. Nine main frequencies [coloured lines in Fig. 8 (c)] are identified in the spectrum which nicely coincide with the nine transitions marked in the energy diagram of Fig. 8 (b) [vertical lines marking $\lambda_{i}-\lambda_{j}$, where $i,j\,\in\{+,-,3,4\}$, and $\lambda_{i}$ are defined in Eq. (100)]. For comparison, we show in Fig. 8 (d) the results obtained using a sufficiently large number of vibrations to achieve results converged with respect to the size of the vibrational subspace. In this case, more spectral features appear [we observe higher order transitions and further peak splitting when compared with Fig. 8 (c)]. Nonetheless, the simple model introduced in this section still explains very satisfactorily the spectral positions of the strongest peaks.

Appendix F Description of the system for two vibrational modes

We consider next a molecule with two vibrational modes. To that end, we consider the plasmon as a bath (as described in Appendix C), which allows writing the following reduced Hamiltonian describing the dynamics of the molecule:

[TABLE]

where $b_{1(2)}[b^{\dagger}_{1(2)}]$ is the annihilation (creation) operator of the vibrational mode 1 (2), $\Omega_{1(2)}$ is the vibrational frequency of the respective mode and $d_{1(2)}$ is the displacement of the respective excited PESs with respect to the ground state ones. According to Eq. (25) we define two Lindblad superoperators that effectively account for the Markovian damping of the two vibrational modes:

[TABLE]

with $\gamma_{b_{1}}$ ($\gamma_{b_{2}}$) being the respective decay rates of the vibrational modes.