Molecular optomechanics in the anharmonic regime: from nonclassical mechanical states to mechanical lasing

Cavity optomechanics aims to establish optical control over vibrations of nanoscale mechanical systems, to heat, cool or to drive them toward coherent, or nonclassical states. This field was recently extended to encompass molecular optomechanics: the dynamics of THz molecular vibrations coupled to the optical fields of lossy cavities via Raman transitions. The molecular platform should prove suitable for demonstrating more sophisticated optomechanical effects, including engineering of nonclassical mechanical states, or inducing coherent molecular vibrations. We propose two schemes for implementing these effects, exploiting the strong intrinsic anharmonicities of molecular vibrations. First, to prepare a nonclassical mechanical state, we propose an incoherent analogue of the mechanical blockade, in which the molecular anharmonicity and optical response of hybrid cavities isolate the two lowest-energy vibrational states. Secondly, we show that for a strongly driven optomechanical system, the anharmonicity can suppress the mechanical amplification, shifting and reshaping the onset of coherent mechanical oscillations. Our estimates indicate that both effects should be within reach of existing platforms for Surface Enhanced Raman Scattering.


I. INTRODUCTION
Recently, in response to surprising experimental results [1,2], it has been suggested that Raman scattering of light from molecules in plasmonic cavities can be cast as an optomechanical process [3][4][5][6][7][8], with the molecular vibrations modes playing the role of ultra-high frequency mechanical resonators.This realization brought the vast set of tools developed for canonical cavity optomechanics to the field of Surface-or Tip-Enhanced Raman Scattering (SERS and TERS) research [9][10][11][12].The resulting formalism of molecular optomechanics led to new insights into the correlations of the inelastically scattered Raman light [4,13], control over the quantum-mechanical description of the single-and multi-mode plasmonic cavities [5,6,14,15], or the dynamics of systems with multiple molecules [16].It also enabled the theoretical proposals [17], and experimental demonstrations of new THz detection techniques [18,19].
Simultaneously, molecular optomechanics stretched the landscape of the conventional cavity optomechanics [20] towards the largely unexplored regimes of ultrahigh mechanical frequencies characteristic of molecular vibrations, complex (multimode) optical spectrum, and to systems with hundreds, or thousands of homogeneous, and both directly and indirectly coupled mechanical modes.It also brought optomechanics closer to the elusive limit of the strong single-photon coupling [21,22], by confining molecules in ultra-small volume optical cavities [23].
Unlike in the canonical optomechanical systems, the dynamics of the THz molecular vibrations involves only a * mikolaj.schmidt@mq.edu.au 1.Schematic of the anharmonic molecular optomechanics setup.Molecule, placed in the gap of a plasmonic cavity, experiences off-resonant Raman transitions between the vibrational states of its ground electronic manifold.The anharmonic nature of the potential yields uneven spectrum of the vibrational modes, which can be either used to separate the dynamics of the two lowest-energy states, implementing an acoustic two-level system, or suppress the formation of acoustic lasing transition.
few, lowest-energy mechanical states, thanks to the combination of low thermal population (n th b < 0.1) [9,10,24], large mechanical losses, and small populations of the optical cavities, which render the mechanical amplification mechanism ineffective [3,4,23].Therefore, in the original formulation of molecular optomechanics [3,4], and follow-up contributions, the molecular vibrations was routinely approximated by a harmonic model of the potential (see the schematic representation in Fig. 1).However, recent experiments (see e.g.Ref. [7]) are beginning to explore the more efficient amplification mechanisms, setting up questions about the role of the intrinsic anharmonicity of the mechanical potential [25][26][27].To date, these effects were simply neglected, and their potential experimental observations attributed with other effects, such as the bond breaking [7].At the same time, new types of hybrid plasmonic-dielectric cavities, characterized with small optical mode volumes and sharp spectral features [6,14,28,29] offer realistic designs of systems which could resolve these anharmonicities.
Therefore, in this work we ask if this anharmonicity can be harnessed for new physics, and how do the predictions of cavity optomechanics hold up in a system with a strong mechanical nonlinearity: Can the anharmonic mechanical potential open up a pathway towards vibrational quantum nonlinearity [30], and expand the toolbox of nonlinear mechanical phenomena explored to date in optomechanics [31][32][33][34]?Can we use it to engineer nonclassical mechanical states of vibrations [35] without employing external nonlinear elements [36,37], or can the current experiments induce coherent mechanical lasing [38][39][40] in the presence of the mechanical linearity?
The manuscript is structured as follows: In Section II we introduce the formalism, and corrections to the conventional framework of molecular optomechanics, including the redefined coupling mechanism.We then show how this anharmonicity can be harnessed to prepare the molecular vibrations in a nonclassical state (Section III), and how the anhamornicities modify the mechanical amplification, and the onset of mechanical lasing in molecular systems (Section IV).

II. FORMALISM
In the elementary formulation of molecular optomechanics, we consider a single quantized optical mode, coupled through nonlinear interaction to a single mechanical mode [3,4].The dynamics of this setup is described by the sum of the optical, mechanical, and interaction Hamiltonians: The optical mode has resonant frequency ω a , and is coherently driven with amplitude Ω and frequency ω l , so that in the frame rotating with ω l we have Ĥopt = (ω a − ω l )â † â + Ω(â + â † ).To explicitly consider the harmonic or anharmonic characteristics of molecular vibrations, here we explicitly write the mechanical Hamiltonian as with V M representing the Morse potential The eigenfrequency of the kth, out of N bound states of ĤM , is with δω b = ω 2 b /(4D e ), and N = (ω b /δω b − 1)/1.The optomechanical interaction between molecular system (characterized by the position operator x) and the optical cavity mode (with electric field operator Ê) is mediated by the Raman dipole, induced in the molecule with Raman polarizability tensor R by the optical cavity field as pR = Rx Ê. ( The explicit connection to the molecular vibrations is given by representing x in the basis of the eigenstates |φ k of the mechanical Hamiltonian: x k,j σk,j where σkj = |φ k φ j | denotes the transition operator.The analytical expressions for the matrix elements x kj = φ k | x |φ j are given in Appendix A. Note that in contrast to the harmonic model of optomechanics, the anharmonic potential introduces diagonal components to the x operator, represented by x(0) .The interaction Hamiltonian of the system takes the form [5,15,41]: where Ê(r) = E 0 (r)â † + h.c., and we carried out the rotating wave approximation to remove the optical modesqueezing terms (∝ â2 , (â † ) 2 ).
We note that this framework explicitly assumes an off-resonant nature of the Raman scattering from a single molecule -see Refs.[8,38,42] for the models of the Surface-Enhanced Resonant Raman Scattering (SERRS), and Ref. [15] for the extension of the off-resonant molecular optomechanics towards multiple molecules.
Since the cavity is coherently driven, and only weakly coupled to the mechanical system, the optical mode fluctuates around a coherent state with amplitude â where κ is the optical decay rate.We can then linearize the interaction Hamiltonian by separating the coherent part of the optical mode from the fluctuations â = δâ + â ≈ δâ + α as: (8) where G 0 = E 0 (r d )RE * 0 (r d ), with r d denoting the position of the molecule.
From here, the optomechanical linearization neglects the terms ∝ δâ † δâ to write Ĥom ≈ Ĥ(0 where we assumed, without the loss of generality, that α can be made real.The diagonal contributions to the displacement operator x(0) define The second term in om introduces a shift of the frequency of each vibrational level |φ i , mentioned also in Ref. [43], proportional α 2 G 0 x k,k .We thus dress the Morse potential eigenfrequencies given in Eq. ( 4) as For the Morse potential we can find a good approximation for the diagonal matrix elements x k,k (see Appendix A).For reference, we note that the difference between the neighboring eigenfrequencies is and can be approximated using the analytical expressions for the diagonal matrix elements x k,k .In particular, as we show in Appendix A, the shift due to the diagonal term is nearly constant for all k.

A. Master equation for the anharmonic molecular vibrations
From here, we can formulate the effective description of the mechanical state by embracing the quantum noise approach [44,45], treating the driven optical mode as a structured reservoir, and following the evolution of the reduced density matrix ρ of the mechanical system.The coupling to the mechanical mode is then determined by the interaction Hamiltonian [8,17], including both Ĥ(±) om and the first term in Ĥ(0) om in Eq. (10).In the secular approximation, interaction Ĥ(±) om dictates that the mechanical excitation and decay rates will be given by the spectra of two-time correlators [δâ(τ )] † δâ(0) , calculated at the dressed transition frequencies ωk − ωj : erator, and g = αG 0 x zpf is the effective optomechanical coupling rate, and x zpf = /(2mω b ) is the zero-point fluctuation of the harmonic oscillator with frequency ω b .For the more direct comparison with the canonical cavity optomechanics, we normalize the jump operators by the corresponding matrix elements x k,j /x zpf , and separate the first two terms which describe the Stokes (mechanical excitation), and anti-Stokes (mechanical relaxation) processes, respectively.The remaining two terms describe the effects of the coupling to a thermal bath.For simplicity, we assume that those would be completely described by transitions between neighboring eigenstates, through the constant mechanical decay rate γ, and the thermal bath population approximated by the Bose-Einstein population n th b at frequency ω b .
Finally, we note that the mechanical system will exhibit an entirely incoherent dynamics, and thus omit the effect of the interaction term −G 0 αx (0) (δâ † + δâ) in the interaction Hamiltonian Ĥ(0) om (Eq.( 10)), which does not change the mechanical state, and yields an irrelevant, dephasing-like term.

Multiple optical modes
Realistic optical systems used in SERS, like the plasmonic nano-and pico-cavities [7,18,19,23], or hybrid metallic-dielectric systems [6,14], support multiple overlapping and interacting optical modes, which significantly influence the optomechanical dynamics.In particular, the driving, Stokes and anti-Stokes processes (identified in canonical optomechanics with phonon heating and cooling) can all be mediated via different optical modes.
An extension of the single-mode model to a multi-mode one presents several difficulties: for example, the incident laser can couple to more than one cavity mode, complicating the definition of coherent amplitude α; similarly, the optomechanical coupling parameters g 0 needs to be redefined to explicitly account for coupling with modes of the cavity with different field distributions.These difficulties are typically addressed by generalizing the master equations presented above, following the prescriptions from [5,15,43], which introduce the explicit Stokes Γ − rates, calculated using the Green's function of the system, which completely account for the complex position-and frequency-dependent field distributions of the electric field in the cavity.(see Appendix B for the definitions and discussion).Adapting this approach to the anharmonic systems, and assuming that the transitions are limited to the neighboring mechanical states, we can rewrite Eq. ( 13) as In particular, for the system with a single optical mode, we can formally write where S opt,single (ω) = dτ exp(iωτ ) [δâ(τ )] † δâ(0) , simplifying Eq. ( 14) to Eq. ( 13).To maintain this intuitive picture and the connection to the canonical optomechanics in the multi-mode case, throughout this work we will assume that the laser selectively couples to only one particular mode of the cavity, and define |α| 2 as the cavity population.Furthermore, we will fold the entire frequency dependence of the Stokes and anti-Stokes rates into the generalized optical spectrum S opt (ω), while keeping g 0 constant.This is a fairly informal step, but it will allow us to develop an intuitive picture of the anharmonic molecular dynamics in multi-mode optical systems.The master equation in Eq. ( 14) sets up the dynamics of the system in terms of the diagonal elements of the mechanical density matrix, or the populations of the vibrational states ρ k,k , and the total transition rates from the kth state, which include the contributions from the unstructured thermal bath 16) We list the dynamical equations for these population in Appendix D and, in the following sections, investigate their solutions in two cases: weakly pumped, strongly anharmonic system, and strongly pumped system with weaker anharmonicity.

III. NONCLASSICAL MECHANICAL STATES
To prepare nonclassical states of molecular vibrations in an incoherently driven anharmonic system, we need to suppress its the excitation beyond the two lowestorder states {|φ 0 , |φ 1 } [30] -that is, form a phonon blockade.Similar schemes have been explored in other branches of physics, most famously in circuit QED, where the Kerr nonlinearity enables the generation of nonclassical states of superconducting circuits [46].However, that functionality is enabled by the presence of the coherent microwave drive, which induces transitions between specific levels of the anharmonic ladder.In the molecular optomechanics, as well as the canonical cavity optomechanics, all transitions are due to incoherent processes, and so we turn to engineering the rates of these incoherent processes Γ (k) ± defined above, by structuring the optical spectrum of the system.
An example of a system with the desirable spectrum is depicted in Fig. 2(c), where we show a (not to scale) schematic of a hybrid metallic-dielectric cavity, explored in several recent studies [6,14,28,29].Here, a dimer of gold nanoparticles supports a lossy (Q ∼ 10) plasmonic mode with dramatically reduced effective mode volumes ∼ 10 −6 λ 3 [23], and coupled to a high-Q dielectric microresonator in the form of toroidal [6], or nanobeam cavity [14].The optical response S opt of this system, defined in the previous section, is shown in Fig. 2(a), and exhibits a strong Fano feature due to the off-resonant interaction between the high-and low-Q optical modes.Expressions used to model S opt in the presence, and absence of coupling between the modes (depicted with dashed line in Fig. 2(a)), and parameters used in this work, are listed in Appendix E.
In an idealized scheme, we realize the incoherent blockade by tuning the first Stokes frequency ω l − (ω 1 − ω0 ) to matche the peak of the Fano feature, so that state |φ 1 can be efficiently populated from |φ 0 ; additionally, if the second Stokes frequency ω l −(ω 2 − ω1 ) matches the dip in the optical spectrum, transitions to the |φ 2 state become suppressed.This incoherent blockade should be therefore governed by the contrast between the optical spectra calculated at S opt (ω l − ω1 + ω0 ) and S opt (ω l − ω2 + ω1 ).A similar scheme, involving multiple high-Q optical modes for controlling the state of a MHz mechanical oscillator with Kerr nonlinearity, was proposed by Rips et al. [35].
To characterize the population of the mechanical states, and its non-classical statistics, in a manner most resembling the usual second quantization language of populations and statistics of the harmonic system, in Appendix C we employ the position operator x as the key observable, used to define the steady-state mechanical populations and intensity correlations We note that as the system becomes anharmonic, and we turn away from the second quantization framework, n x losses its exact definition as the phonon population.Nevertheless, we embrace this language for simplicity, and a direct mapping to the harmonic setup.Furthermore, if we measured the direct IR emission from the transitions between the mechanical states, these magnitudes would characterize the intensity, and statistics of the emitted Γ - Γ + Γ + (0) IR light.Further details about the calculations of n x and g (2) x (τ ) are listed in Appendix D. In Fig. 2(d,e) we plot with solid lines the populations n x and intensity correlations g x (0), as a function of the laser frequency ω l .The former are visibly suppressed when laser is chosen to have its first Stokes frequency ω l − (ω 1 − ω0 ) match the dip in the optical cavity spectrum, around 498 GHz.Conversely, when the second Stokes transition at ω l − (ω 2 + ω1 ) matches the dip in S opt for the laser around 495 THz, we approach the blockade condition described above, the system demonstrates sub-Poissonian statistics g (2) x (0) < 1. Away from these features, the system acquires the thermal statistics g (2) x (0) ∼ 2. In Fig. 2(f,g) we show the same magnitudes, calculated as a function of nonlinearity δω b , and coherent cavity population |α| 2 .Here again the statistics diverges from the thermal one towards sub-Poissonian (denoted as a blue region), when the second Stokes transition frequency is tuned to the dip of the Fano feature in the optical spectrum S opt .Since that anti-Stokes frequency explicitly depends on |α| 2 due to the dressing by the coherent field (see Eq. ( 12)), the region of sub-Poissonian statistics is largely diagonal.In Appendix F we discuss how this region of sub-Poissonian statistics changes with the laser frequency.
To gain analytical insight into these effects, we consider the analytical solution to the coupled equations for the populations ρ k,k of the three lowest-energy states (see derivation in Appendix D 1), approximate the numerator and denominator in the definition of g x (0) (Eq.( 18)) by 2ρ 2,2 and ρ 1,1 , respectively, to find g (2)  x (0) ≈ 2 ≈ 2 Γ(1) We plot the function given in the first line in Fig. 2(e) with the dashed line, finding a qualitative agreement of the spectral range corresponding to the sub-Poissonian statistics with the full calculations (solid line).We have verified that the discrepancy is due to the truncation to the three states.
The second line represents a far more crude approximation, where we assume that the anti-Stokes transitions rates Γ(k) − are largely independent of k, and far larger than the second anti-Stokes rate Γ(1) + .The first fraction in this expression directly characterizes the contrast in the anti-Stokes due to the Fano feature of the optical cavity.
As we discuss in more detail in Appendix F, we can further suppress the intensity correlations by increasing the relative role of the optomechanical feedback over the thermal pumping.This can be achieved by employing a larger intensity of the optical driving, although one needs to account for the |α| 2 dependence of the dressed mechanical frequencies (Eq.( 12)), to ensure that the transitions |φ 0 → |φ 1 and |φ 1 → |φ 2 match the peak, and the trough of the optical spectrum.One could also explore hybrid resonators featuring larger contrasts of the peak and troughs of the Fano resonance.
In numerical modelling, we assumed a single molecule exhibiting optomechanical coupling g 0 /2π = 2 THz, consistent with the values reported for picocavities in Ref. [23], and below the estimates for coupling with a small ensemble of about 100 molecules in nanocavities in Ref. [7].The non-classical state of vibrations in Fig. 2(g) was reported for aharmonicities as small as δω b /2π = 0.1 THz, which yields the anharmonicity parameter ξ = δω b /ω b ≈ 5 × 10 −3 -arguably large, but reportedly accessible with molecular systems investigated in the context of SERS [43].We also choose to plot the results against the population of the driven optical (plasmonic) mode, rather than the input laser intensity.The low populations explored here are typical for a strongly driven, lossy plasmonic cavities [7,23], and should offer good approximation to the characteristics of hybrid systems under the assumption that the laser predominantly drives the plasmonic mode.
Finally, we note that since our scheme is based around harnessing the changes to the optical spectrum between the two Stokes transitions, it does require the mechanical system to exhibit a strong nonlinearity δω b to resolve the features of the optical spectrum.However, it does not require the mechanical system to operate in the conventional phonon blockade regime δω b > γ.

IV. MECHANICAL AMPLIFICATION AND LASING
Anharmonicity of molecular vibrations should also have a strong effect on the opposite regime of molecular optomechanics, where the system is strongly optically driven, in the effort to amplify the mechanical mode, and boost the intensity of Stokes emission [10,12].In this scenario, the amplification is likely to be suppressed until it saturates the ladder of bound states of the vibrations, or the non-resonant model of Raman scattering breaks down.We explore these effects in Section IV A.
Moreover, beyond the amplification regime, canonical optomechanical systems exhibit a transition to mechanical lasing or dynamical instability, where the mechanical component exhibits coherent oscillations [38,39].In Section IV B we explore similar effects in the context of molecular optomechanics, analyzing the impact of anharmonicity on the threshold, amplitude, and the trajectory of these oscillations.

A. Amplification
Using the formulation developed in Section III, we can readily explore the mechanical amplification by analyzing the steady-state populations mechanical n x as a function of the input pump power, or the cavity population.To this end, we consider a simpler optical setup, with a single optical mode (S opt = S opt,single ) being driven with a blue-detuned laser shown schematically in Fig. 3(a), to promote the Stokes emission, and mechanical amplification.To describe the anharmonicity, we approximate the rates of the non-thermal Stokes and anti-Stokes processes by expanding the cavity spectrum S opt (ω) as and where we used Eq. ( 12), and the explicit form of the diagonal matrix elements (Eq.( A8)).Parameters η ± are defined as the derivatives of the optical spectra near the frequencies given as arguments of the optical spectra in second lines of the above equations.Per the definition of the Raman transition rates (see Eq. ( 15)), we can approximate them using the above expansion as: where ].The exemplary optical spectrum with exaggerated anharmonic shifts, and the corresponding rates, is schematically depicted in Fig. 3(a).Using this formulation, in Appendix D 2 we show that the rate equation for the mechanical population n x is thus modified from its conventional, linear form, into In a single-mode optical cavity, the mechanical lasing setup would require driving on the blue side of the optical resonance (see Fig. 3(a)), in which case the slope of spectrum S opt at ω l + ω b should be negative (η − < 0).Thus, we can interpret the additional terms in the rate equation as damping with the linear and quadratic dependence on the mechanical population n x , and a constant term; all are also dependent on the driving intensity, through the coherent cavity population |α| 2 in g 2 , and the frequencies at which derivatives η ± are calculated.− ), between every two neighboring mechanical levels depend on the spectra of the singlemode optical (plasmonic) cavity.The shifts are exaggerated for clarity.(b) Phonon populations nx calculated using the full model (solid lines), analytical solution to the Eq. ( 27) (dashed lines) and standard deviations of the mechanical oscillations σx (filled areas) as a function of the optical cavity population.Black line denotes results for the harmonic system, while the red, and teal lines correspond to increasing anharmonicity δω b /2π = 0.1, and 0.2 THz.
In Fig. 3(c) we show the steady-state mechanical population n x as a function of the cavity population for several values of δω b .Dashed lines denote the n x derived from Eq. ( 27), by setting its LHS to 0, and solving the resulting quadratic equation for n x .Solid lines are calculated by solving the complete set of rate equation, as described in Appendix D. For reference, we include the result obtained with the harmonic oscillator (dashed black line), which clearly demonstrates divergence, signalling the conventional threshold of the mechanical lasing.
From Fig. 3, we can derive several observations: (i) Even for the smallest anharmonicity (red lines, δω b /2π = 0.1 THz) the mechanical amplification is significantly suppressed, compared to the harmonic case; (ii) n x exhibits a limited superlinear dependence on the cavity population and thus driving intensity; for the stronger nonlinearity (teal lines, δω b /2π = 0.2 THz), n x is either linear or supralinear with |α| 2 , in a significant deviation from the harmonic optomechanical models, (iii) for neither anharmonic system does the population exhibit a clear mechanical lasing threshold.
Despite the anharmonicity suppressing the heating of the vibrations, they appear to build up a substantial population > 10.While this is still far below the limit of N bound states, the Raman transitions between the higher-energy states are likely to couple to the electronic excited states of the molecule, in a resonant Raman fashion.This effect is naturally absent in the canonical cavity optomechanics.

B. Dynamical instability
The linearized coupling theory of optomechanics predicts a natural limit to the mechanical amplification, when the system reaches a threshold of the dynamical instability (or phonon lasing), at which point the incoherent phonon population should diverge.In reality near that threshold the linearized formulation fails, the phonon population remains finite, and the mechanical mode begins to exhibit coherent oscillations.These oscillations grow quickly until the system reaches a steady state, with their amplitude dependent on the optical driving strength and detuning from the optical cavity, and the optomechanical coupling [39,40,47].
Here, we ask if an optomechanical system with an anharmonic, Morse potential, would similarly exhibit steady-state mechanical solutions near the mechanical lasing threshold.While the complete mapping of this effects -possible in the harmonic case -goes beyond the scope of this work, we explore it in the range of parameters relevant for molecular optomechanics.

Harmonic potential
The classical trajectories for the optomechanical system can be identified from the dynamical equations for the optical amplitude α, and the position and momentum of the mechanical oscillator x and p: In the mechanical lasing regime, the mechanical trajectory settles into oscillations x(t) = x 0 + A cos(ω b t), and thus the nonvanishing amplitude A is typically used to characterize both the onset, and magnitude of the oscillations [20,39,48].We reproduce that result numerically, by rewriting the above coupled equations into a form more suitable for numerical simulations (see Appendix H and Ref. [48]), and solve them with initial conditions α(0) = x(0) = p(0) = 0, until the mechanical mode settles into steady-state oscillations.We then characterize these oscillations by their standard deviation σ x = (x(t) − x t ) 2 t , with the temporal average ... t calculated in the steady-state (and nominally equivalent to A).Since these oscillations represent coherent dynamics, we introduce the equivalent coherent mechanical population defined as n x,coh = (σ x /x zpf ) 2 /4.The normalized deviation σ x /x zpf , and the corresponding n x,coh for the harmonic system are denoted in Fig. 3(c) as a filled gray area.

Anharmonic oscillator
For the anharmonic oscillator, we can write down similar equations for the amplitude of the optical mode, and coordinates of the mechanical oscillator; the only difference will be the change to the term in Eq. (28c) which describes the force F = −V (x) derived from the harmonic potential (−x ω b /(2x 2 zpf )) to the corresponding force associated with the Morse potential (F = −V M (x) = −2D e a[1 − exp(−ax)] exp(−ax)): Unlike for the harmonic oscillators, we do not have an analytical solution for the amplitude of the oscillations in the anharmonic system.Nevertheless, we can again characterize them by the standard deviation σ x /x zpf calculated in the steady state, and the corresponding n x,coh .We depict these results in Fig. 3(c) as the filled areas corresponding to the anharmonic systems with δω b /2π = 0.1 (red area), and 0.2 THz (teal area).These results indicate that mechanical lasing should be possible even within the very limited space of the very few bound states of the mechanical system.For example, for δω b /2π = 0.2 THz (teal area), the Morse potential supports about N = 50 states, and the coherent oscillation have the equivalent coherent population n x,coh 5, comparable to the incoherent populations n x .Much like in Raman or Brillouin lasers, this should lead to observable narrowing of the Raman lines [49,50], should such measurements be possible.In Appendix H we show how the thresholds, and amplitudes of these coherent oscillations change with anharmonicity, finding that the anharmonicity does not introduce a qualitative change to the amplitude of mechanical oscillations, except for a more pronounced and shifted threshold behaviour, and lowered amplitudes.

V. CONCLUSION
In this work, we show that by the framework of molecular optomechanics offers pathways towards genuine quantum engineering of molecular vibrations, far beyond the typical applications to mechanical amplification or Surface-Enhanced Raman Scattering.
In particular, by embracing the anharmonicity of the molecular vibrations, and engineering the optical spectrum of hybrid cavities, we show how mechanical systems can be driven into the weakly populated, nonclassical states.We present a complete theoretical framework for designing and describing such systems, and formal connection that can serve to characterize the nonclassical statistics of the emitted THz photons.
In the opposite regime of strong driving, we show that the current experimental setups should be capable of inducing coherent oscillations, or mechanical lasing, of the molecular vibrations.This effect should lead to an observable narrowing of the Raman scattering lines, and further the mapping between molecular, and canonical optomechanics.Furthermore, we show that even weak anharmonicities can dramatically change the optical response of the system, reshaping the dependence of the Stokes intensity on driving power.
These findings should significantly expand the toolbox and the impact of the formalism of molecular optomechanics, towards applications in quantum sensing and microscopy, and call for revisiting of the reported experimental results, to analyze the possible impact of the anharmonicities.
Schrödinger equation with Morse potential: has eigenstates denoted as |φ n and corresponding eigenfrequencies In the Dunham expansion we have Y 1,0 = ω b , and Y 2,0 = −( ω b ) 2 /(4D e ), and δω b = ω 2 b /(4D e ).The matrix elements for the eigenstates of the Morse potential, given in Ref. [51], and normalized by the zeropoint fluctuation, are for n > m, with N = (ω b /δω b − 1)/2, and using the Gamma function Γ.We note that these elements are symmetric x n,m = x m,n .
In general, we find that these elements do not vanish for non-neighboring states (n − m = ±1), but drop-off quickly with |n − m|.For the neighboring states, we can simplify the above formulation by using the recursive property of the Γ function Γ(z+1) = zΓ(z).In particular, for m = n + 1, we find For n = m we find with the digamma function Ψ.In the limit of weak anharmonicity and for low-order states N n, we can use its asymptotic property Ψ(z) ≈ ln(z) − 1/(2z) to approximate the diagonal terms x n,n as 1. Shifts of Raman spectra due to the diagonal terms Dropping the small, n-independent constant term in Eq.A9, we can use it to approximate the Stokes transition frequency (see Eq. ( 12 by the position operator x.We use it as an observable which carries information about the excitations, and nonclassical nature of the mechanical state, much like the electric field operator used to define the spectrum and statistics of emission from the system.We thus define and the intensity correlations as g (2)  x (τ ) = G x (τ ) n 2 x .

(C3)
For the mechanical system in the mixed state ρ = i p i |φ i φ i |, we can express the two magnitudes as  (D2) We will now construct an algebraic formulation of the problem of finding the steady-state solution to the finite set of equations for {p 0 , p 1 , ..., p K }.In the equation for p 0 we drop the two terms that describe emission from, and excitation to |φ 0 : Finally, we can turn these coupled equations into an inhomogeneous system by replacing the equation for ṗK with the condition k p k = 1.This system of ODEs can be expressed as Equations ( 28) can be rewritten into a form more suitable for numerical solution by normalizing coordinates x = x/x zpf , p = px zpf / , and introducing time parameter τ = tω b : where we introduced ã = ax zpf .Notably, as ã approaches 0, we recover the harmonic restoring force.

Mechanical lasing amplitudes
In Fig. 6 we revisit the results included in Fig. 3(c), to show how the amplitudes of mechanical lasing change with small anharmonicity.We note that the threshold clearly shifts towards larger cavity populations (included in the definition of the effective optomechanical coupling g), and the oscillations becomes initially more steep.
FIG.1.Schematic of the anharmonic molecular optomechanics setup.Molecule, placed in the gap of a plasmonic cavity, experiences off-resonant Raman transitions between the vibrational states of its ground electronic manifold.The anharmonic nature of the potential yields uneven spectrum of the vibrational modes, which can be either used to separate the dynamics of the two lowest-energy states, implementing an acoustic two-level system, or suppress the formation of acoustic lasing transition.
Appendix D: Solving the rate equationsLet us recall the definitions of rates given in Eq. (16), and the final form of the master equation (Eq.14), to write down the rate equations for the populations p k which define the mixed state density matrix ρ = k p k |φ k φ k |: for the harmonic system, wherex k+1,k /x zpf = √ k + 1 and x k,k+1 /x zpf = √ k + 1,and Γ(k) − are independent of k, the above simplifies to d dt p k =(k + 1) Γ− p k+1 − k Γ− + (k + 1) Γ+ p k + k Γ+ p k−1 .

FIG. 6 .
FIG.6.Dependence of the amplitude of the coherent mechanical oscillations on the anharmonicity δω b , and populations of the optical cavity near the threshold.As in Fig.3(c), we quantify the oscillations through the equivalent coherent population n x,coh , and the normalized standard deviation σx/x zpf .