Nonlinear response of a driven vibrating nanobeam in the quantum regime

We analytically investigate the nonlinear response of a damped doubly clamped nanomechanical beam under static longitudinal compression which is excited to transverse vibrations. Starting from a continuous elasticity model for the beam, we consider the dynamics of the beam close to the Euler buckling instability. There, the fundamental transverse mode dominates and a quantum mechanical time-dependent effective single particle Hamiltonian for its amplitude can be derived. In addition, we include the influence of a dissipative Ohmic or super-Ohmic environment. In the rotating frame, a Markovian master equation is derived which includes also the effect of the time-dependent driving in a non-trivial way. The quasienergies of the pure system show multiple avoided level crossings corresponding to multiphonon transitions in the resonator. Around the resonances, the master equation is solved analytically using Van Vleck perturbation theory. Their lineshapes are calculated resulting in simple expressions. We find the general solution for the multiple multiphonon resonances and, most interestingly, a bath-induced transition from a resonant to an antiresonant behavior of the nonlinear response.


Introduction
The experimental realization of nanoscale resonators which show quantum mechanical behavior [1,2,3,4,5] is currently on the schedule of several research groups worldwide and poses a rather non-trivial task. Important key experiments on the way to this goal have already been reported in the literature [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20] and are also reviewed in this Focus Issue. Most techniques to reveal the quantum behavior so far address the linear response in form of the amplitude of the transverse vibrations of the nanobeam around its eigenfrequency. The goal is to excite only a few energy quanta in a resonator held at low temperature. To measure the response, the ultimate goal of the experiments is to increase the resolution of the position measurement to the quantum limit [11,17,18,21,22]. As the response of a damped linear quantum oscillator has the same simple Lorentzian shape as the one of a damped linear classical oscillator [23], a unique identification of the "quantumness" of a nanoresonator in the linear regime can sometimes be difficult.
One possible alternative is to study the nonlinear response of the nanoresonator which has been excited to its nonlinear regime. A macroscopic beam which is clamped at its ends and which is strongly excited to transverse vibrations displays the properties of the Duffing oscillator being a simple damped driven oscillator with a (cubic) nonlinear restoring force [24]. Its nonlinear response displays rich physical properties including a driving induced bistability, hysteresis, harmonic mixing and chaos [24,25,26]. The nonlinear response of (still classical) nanoscale resonators has been measured in recent experiments [9,10,11,16]. In the range of weak excitations, the standard linear response arises while for increasing driving, the characteristic response curve of a classical Duffing oscillator has been identified.
No signatures of a quantum behavior in the nonlinear response of realized nanobeams have been reported up to present. One reason is that a nanomechanical resonator is exposed to a variety of intrinsic as well as extrinsic damping mechanisms depending on the details of the fabrication procedure, the experimental conditions and the used materials [20,27,28,29]. Possible extrinsic mechanisms include clamping losses due to the strain at the connections to the support structure, heating, coupling to higher vibrational modes, friction due to the surrounding gas, nonlinear effects, thermoelastic losses due to propagating acoustic waves, surface roughness, extrinsic noise sources, dislocations, and other material-dependent properties. An important internal mechanism is the interaction with localized crystal defects. Controlling this variety of damping sources is one of the major tasks to be solved to reveal quantum mechanical features. Recent measurement show that in the so far realized devices based on silicon and diamond structures, damping has been rather strong at low frequencies [27,28] indicating even sub-Ohmic-type damping [23] which would make it difficult to observe quantum effects at all. However, using freely suspended carbon nanotubes [30,31] instead could reduce damping at low frequencies due to the more regular structure of the long molecules which can be produced in a very clean manner. Further experimental work is required to clarify this point and to optimize the experimental conditions. Nanoscale nonlinear resonators in the quantum regime have been investigated theoretically starting from microscopic models based on elasticity theory for the beam [32,33,34]. Carr, Lawrence and Wybourne have considered an elastic bar under static longitudinal compression beyond the Euler instability leading to two stable equilibrium positions around which the transverse vibrations of the beam occur. It turned out that quantum tunneling between the two minima is in principle possible in silicon beams and carbon nanotubes. However, the strain has to be controlled with extreme accuracy and the quantum fluctuations in position are of the order of 0.1Å. The detection of such small lengths certainly is challenging. However, a possible method to increase the resolution could be the use of the phenomenon of stochastic resonance [35] for a coherent signal amplification of the nonlinear response of nanomechanical resonators in their bistable regime [36]. Werner and Zwerger [33] have studied a similar setup close to the Euler buckling instability which occurs at a critical strain ǫ c . There, the frequency of the fundamental mode vanishes and quartic terms in the Lagrangian have to be taken into account. An effective Hamiltonian has been derived for the amplitude of the fundamental mode being the dynamical variable which moves in an anharmonic potential. Depending on the strain ǫ being below (ǫ < ǫ c ) or above (ǫ > ǫ c ) the critical value, a monostable or bistable situation can be created. The conditions for macroscopic quantum tunneling to occur have been estimated for the bistable case. In order to measure single-phonon transitions in a nanoresonator, it has been proposed to use its anharmonicity together with a second nanoresonator acting as a transducer for the phonon number in the first one [37]. In this way, the measured signal being the induced current is directly proportional to the position of the read-out oscillator.
In Ref. [34], we have considered a similar setup but restricted to the statically monostable case below the Euler instability, i.e., for ǫ < ǫ c . In addition, we have allowed for a time-dependent periodic driving force F (t) such that an effective monostable quantum Duffing oscillator arises. Possible origin of the driving can be the magnetomotive force when an ac current is applied and the beam is placed in a transverse magnetic field. Moreover, a (weak) influence of the environment has been modeled phenomenologically by a simple Ohmic harmonic bath. The nonlinear response has been determined numerically from solving a Born-Markovian master equation for the reduced density operator of the system after the bath has been traced out. We have identified discrete multiphonon transitions as well as macroscopic quantum tunneling of the fundamental mode amplitude between the two stable states in the driving induced bistability. Moreover, a peculiar multiphonon antiresonant behavior has been found in the numerical results for the damped system [38]. The discrete multiphonon (anti-)resonances are a typical signature of quantum mechanical behavior [34,38] and are absent in the corresponding classical model of the standard Duffing oscillator [24,25,26], also when thermal fluctuations are included [39].
While we have approached the problem in Refs. [34,38] by numerical means, we present in this work a complete analytical investigation of the dynamics of the quantum Duffing oscillator. We intend to elucidate the mechanism behind the reported [38] bath-induced transition from the resonant to the antiresonant nonlinear response of the nanobeam. This is achieved by solving a Born-Markovian master equation for the reduced density operator in the rotating frame. Within the rotating wave approximation (RWA), a simplified system Hamiltonian follows whose eigenstates are the quasienergy states. The corresponding quasienergies show avoided level crossings when the driving frequency is varied. They correspond to multiphonon transitions occurring in the resonator. Moreover, we include the dissipative influence of an environment and find that the dynamics around the avoided quasienergy level crossings is well described by a simplified master equation involving only a few quasienergy states. Around the anticrossings, we find resonant as well as antiresonant nonlinear responses depending on the damping strength. The underlying mechanism is worked out in the perturbative regime of weak nonlinearity, weak driving and weak damping. There, Van Vleck perturbation theory allows to obtain the quasienergies and the quasienergy states analytically. The master equation can then be solved in the stationary limit and subsequently, the line shapes of the resonant as well as the antiresonant nonlinear response can be calculated.
The problem of a driven quantum oscillator with a quartic nonlinearity has been investigated theoretically in earlier works in various contexts. In the context of the radiative excitation of polyatomic molecules [40], Larsen and Bloembergen have calculated the wave-functions for the coherent multiphoton Rabi precession between two discrete levels for a collisionless model. More recently, also Dykman and Fistul [41] have considered the bare nonlinear Hamiltonian under the rotating wave approximation. Drummond and Walls [42] have investigated a similar system occurring for the case of a coherently driven dispersive cavity including a cubic nonlinearity. Photon bunching and antibunching have been predicted upon solving the corresponding Fokker-Planck equation. Vogel and Risken [43] have calculated the tunneling rates for the Drummond-Walls model by use of continued fraction methods. Dmitriev, D'yakonov and Ioffe [44] have calculated the tunneling and thermal transition rates for the case when the associated times are large. Dykman and Smelyanskii [45] have calculated the probability of transitions between the stable states in a quasi-classical approximation in the thermally activated regime. Recently, the role of the detector (in this case, a photon detector) has been studied for the quantum Duffing oscillator in the chaotic regime [46]. The power spectra of the detected photons carry information on the underlying dynamics of the nonlinear oscillator and can be used to distinguish its different modes. However, the line-shape of the multi-phonon resonance which is the central object for studying the nonlinear response remained unaddressed so far. In addition, we start from a microscopic Hamiltonian for the bath and present a fully analytical treatment of system and environment in the deep quantum regime of weak coupling.
Our paper is structured such that we introduce the elasticity model for the doubly clamped nanobeam, derive the effective quartic Hamiltonian and discuss the model for damping in Section 2. Then, we discuss the coherent dynamics of the pure system in terms of the RWA and the Van Vleck perturbation method in Section 3. The dissipative dynamics is studied in Section 4, while the observables are defined in Section 5. The solution for the line shapes are given in Section 6 before the final conclusions are drawn in Section 7.

Model for the driven nanoresonator
We consider a freely suspended nanomechanical beam of total length L and mass density σ = m/L which is clamped at both ends (doubly clamped boundary conditions) and which is characterized by its bending rigidity µ = Y I being the product of Young's elasticity modulus Y and the moment of inertia I. In addition, we allow for a mechanical force F 0 > 0 which compresses the beam in longitudinal direction. Moreover, the beam is excited to transverse vibrations by a time-dependent driving field F (t) =f cos(ω ex t). In a classical description, the transverse deflection φ(s, t) characterizes the beam completely, where 0 ≤ s ≤ L. Then, the Lagrangian of the vibrating beam follows from elasticity theory as [33] Before we study the dynamics of the driven beam, we consider first the undriven system with F (t) ≡ 0. For the case of small deflections |φ ′ (s)| ≪ 1, the Lagrangian can be linearized and the time-dependent Euler-Lagrange equations can be solved by the eigenfunctions φ(s, t) = n φ n (s, t) = n A n (t)g n (s), where g n (s) are the normal modes which follow as the solution of the characteristic equation. For the doubly clamped nano-beam, we have φ(0) = φ(L) = 0 and φ ′ (0) = φ ′ (L) = 0. However, it turns out that this situation is closely related to the simpler case that the nano-beam is also fixed at both ends but its ends can move such that the bending moments at the ends vanish, i.e., φ(0) = φ(L) = 0 and φ ′′ (0) = φ ′′ (L) = 0 (free boundary conditions). For the case of free boundary conditions, the characteristic equation yields the normal modes g free n (s) = sin(nπs/L) and the corresponding frequency of the n−th mode follows as ω free At the critical force F c = µ(π/L) 2 , the fundamental frequency ω free 1 (F 0 → F c ) vanishes as √ ǫ, where ǫ = (F c − F 0 )/F c is the distance to the critical force, and the well-known Euler instability occurs.
For the case of doubly clamped boundary conditions, the characteristic equation yields a transcendental equation for the normal modes which cannot be solved analytically. However, close to the Euler instability F 0 → F c , the situation simplifies again. After expanding, one finds for the fundamental frequency Approaching the Euler instability, the frequencies of the higher modes ω n≥2 remain finite, while the fundamental frequency Hence, the dynamics at low energies close to the Euler instability will be dominated by the fundamental mode alone which simplifies the treatment of the nonlinear case, see below. The fundamental mode g 1 (s) can also be expanded close to the Euler instability and one obtains in zero-th order in ǫ g 1 (s) ≃ sin 2 πs L . (3)

Effective single-particle Hamiltonian
Since the fundamental mode vanishes when F 0 → F c , one has to include the contributions beyond the quadratic terms ∝ φ ′2 , φ ′′2 of the transverse deflections in the Lagrangian. The next higher order is quartic and yields terms ∝ φ ′4 , φ ′2 φ ′′2 . Inserting again the normal mode expansion in the Lagrangian generates self-coupled modes k A 4 k as well as couplings terms k,l A 2 k A 2 l between the modes. This interacting field-theoretic problem cannot be solved any longer. However, since the normal mode dominates the dynamics at low energies closed to the Euler instability, one can neglect the higher modes in this regime. Hence, we choose the ansatz φ(s, t) = A 1 (t)g 1 (s) in the regime F 0 → F c and restrict the discussion in the rest of our work to this regime. The so-far classical field theory can be quantized by introducing the canonically conjugate momentum P ≡ −i ∂/∂A 1 and the time-dependent driving force can straightforwardly be included. Note that when the driving frequency is close to the fundamental frequency of the beam, the fundamental mode will dominate also in absence of a static longitudinal compression force. However, a compression force helps to enhance the nonlinear effects which are in the focus of this work. After all, an effective quantum mechanical time-dependent Hamiltonian results which describes the dynamics of a single quantum particle with "coordinate" X ≡ A 1 in a time-dependent anharmonic potential. It reads with the effective mass m eff = 3σL/8 and the nonlinearity parameter α = (π/L) 4 F c L(1+ 3ǫ). The classical analogous system is the Duffing oscillator [24] (when, in addition, damping is included, see below). It shows a rich variety of features including regular and chaotic motion. In this work, we focus on the parameter regime where only regular motion occurs. For weak driving strengths, the response as a function of the driving frequency ω ex has the well-known form of the harmonic oscillator with the maximum at ω ex = ω 1 . For increasing driving strength, the resonance grows and bends away from the ω ex = ω 1 -axis towards larger frequencies (since α > 0). The locus of the maximal amplitudes is often called the backbone curve [24]. The corresponding nonlinear response of the quantum system shows clear signatures of sharp multi-phonon resonances whose line shapes will be calculated below.

Phenomenological model for damping
In our approach, we do not intend to focus on the role of the microscopic damping mechanisms as this depends on the details of the experimental device. Instead, we introduce damping phenomenologically in the standard way [23] by coupling the resonator Hamiltonian Eq. (4) to a bath of harmonic oscillators described by the standard Hamiltonian with the spectral density with damping constant γ s and cut-off frequency ω c . Our results discussed below are valid for an Ohmic (s = 1, γ 1 ≡ γ) as well as for super-Ohmic (s > 1) baths. Sub-Ohmic baths will not be considered here since the weak-coupling assumption which allows the Markov approximation does not hold any longer. Formally, the coefficients in the master equation would diverge in the sub-Ohmic case, see Eq.
To proceed, we scale H tot (t) to dimensionless quantities such that the energies are in units of ω 1 while the lengths are scaled in units of x 0 ≡ m eff ω 1 . Put differently, we formally set m eff = = ω 1 = 1. The nonlinearity parameter α is scaled in units of α 0 ≡ ω 1 /X 4 0 , while the driving amplitudes are given in units of f 0 ≡ ω 1 /x 0 . Moreover, we scale temperature in units of T 0 ≡ ω 1 /k B while the damping strengths are measured with respect to ω 1 .

Coherent dynamics and rotating wave approximation (RWA)
Let us first consider the resonator dynamics without coupling to the bath. For convenience, we switch to a representation in terms of creation and annihilation operators a and a † , such that X = x 0 (a + a † )/ √ 2. Moreover, it is convenient to switch to the rotating frame by formally performing the canonical transformation R = exp [−iω ex a † at]. We are interested in the nonlinear response of the resonator around its fundamental frequency, i.e., for ω ex ≈ ω 1 , and will not consider the response at higher harmonics. We further assume that the driving amplitude f is not too large such that the nonlinear effects are small enough in order not to enter the chaotic regime. This suggests to use a rotating wave approximation (RWA) of the full system Hamiltonian H(t) in Eq. (4) as the fast oscillating terms will be negligible around the fundamental frequency for weak enough driving. By eliminating all the fast oscillating terms from the transformed Hamiltonian, one obtains the Schrödinger equation in the rotating framẽ H|φ α = ε α |φ α with the Hamiltonian in the RWÃ Here, we have introduced the detuningω = ω 1 − ω ex , the nonlinearity parameter ν = 3 α/(4ω 2 1 ), f =f (8 ω 1 ) −1/2 andn = a † a. In the static frame, an orthonormal basis (at equal times) follows as The Hamiltonian (7) has been studied in Refs. [40,41]. The quasi-energy levels for a given number N of phonons are pairwise degenerate, ε N −n = ε n for n ≤ N, vanishing f → 0 andω = −ν(N + 1)/2. For a finite driving strength f > 0, the exact crossings turn then into avoided crossings which is a signature of multiphonon transitions [34,41]. A typical quasienergy spectrum is shown in Fig. 1 for the parameters ν = 10 −3 and f = 10 −4 . The dashed vertical lines indicate the multiple avoided level crossings which occur all for the same driving frequency. For |ε| = |2f /[ν(N + 1)]| ≪ 1, each pair of degenerate levels interacts only weakly with the other levels, and act effectively like a two-level Rabi system [40]. The Rabi frequency is related to the minimal splitting of the levels and is calculated perturbatively with ε as a small parameter in the following Section.

Van Vleck perturbation theory
Let us therefore consider the multiphonon resonance atω = −ν(N + 1)/2. In addition, we are interested in the response around the resonance and therefore introduce the small deviation ∆. We formally rewriteH as Let us then first discuss the dynamics at resonance (∆ = 0). We divide it in the unperturbed part H 0 and the perturbation εV according to respectively. The unperturbed Hamiltonian is diagonal and near the resonance its spectrum is divided in well separated groups of nearly degenerate quasienergy eigenvalues. An appropriate perturbative method to diagonalize this type of Hamiltonian is the Van Vleck perturbation theory [49,50,51]. It defines a unitary transformation yielding the HamiltonianH in an effective block diagonal form. The effective Hamiltonian has the same eigenvalues as the original one, with the quasidegenerate eigenvalues in a common block. The effective Hamiltonian can be written asH In our case, each block is a two by two matrix corresponding to a subspace formed by a couple of quasienergy states forming an anticrossing. Let us consider the effective Hamiltonian H ′ n corresponding to the involved levels |n and |N − n , being eigenstates of the harmonic oscillator. The degeneracy in the corresponding block is lifted at order N − 2n in Van Vleck perturbation theory. The block Hamiltonian then reads where This is the lowest order of the perturbed Hamiltonian which allows to calculate the corresponding zero-th order eigenstates. By diagonalizingH ′ n in Eq. (12), one finds the minimal splitting for the N−phonon transition as For the case away from resonance, we consider a detuning ∆ = ε N δ. Within the Van Vleck technique, only the zero-th block is influenced according tõ the other blocks given in Eq. (12) for n = 0 are not influenced by this higher-order correction. The eigenvectors for the HamiltonianH at zero-th order are obtained by diagonalizingH ′ n in Eq. (12). One finds |φ n = |n for n ≥ N + 1 or |φ n = (|n +|N −n )/ √ 2 and |φ N −n = (|n −|N −n )/ √ 2 for 0 < n < N/2 and |φ N/2 = |N/2 if N is even. Moreover, where we have introduced the angle θ via tan θ = −2Ω N,0 /[ν(N + 1)N∆].

Dissipative dynamics in presence of the bath
Having discussed the coherent dynamics, we include now the influence of the harmonic bath coupled to the driven system. We therefore assume that the coupling is weak enough such that the standard Markovian master equation for the reduced density operator ρ(t) can be applied. The influence of the bath enters in the superoperator with the correlators and The kernels are given by where T is the environment temperature. Moreover, U(t, t ′ ) = T exp(i t t ′ H(t)dt) is the propagator with the time order operator T . Next, we project the density matrix on the orthonormal set |ϕ α (t) = exp [−iω ex a † at]|φ α , such that the matrix elements read Performing the derivative one obtainṡ For the dissipative term, we need to compute with X αβ,+1 = X * βα,−1 = x 0 φ α |a|φ β / √ 2 being the matrix element of the destruction operator in the rotating frame. We need, moreover, In an analogous way, we have Here, we have defined the time derivative ± and Here, N αβ,±1 are defined as in terms of the bath density of states J(|ε|), the bosonic thermal occupation number and the Heaviside function θ(x). Eq. (29) illustrates why we have to restrict to (super-)Ohmic baths, since N(ε) would diverge for s < 1 at low energies. During the calculation, the τ -integration in the double integrals in Eqs. (26) and (25) has been evaluated by using the representation ∞ 0 dτ exp (iωτ ) = πδ(ω) + iP p (1/ω), where P p denotes the principal part. The contributions of the principal part result in quasienergy shifts of the order of γ s which are the so-called Lamb shifts. As usual, these have also been neglected here.
The ingredients can now be put together to obtain the Markovian master equation in the static frame aṡ Next, we perform a 'moderate rotating wave approximation' consisting in averaging the time-dependent terms in the bath part over the driving period T ωex = 2π/ω ex . This is consistent with the assumption of weak coupling which assumes that dissipative effects on the dynamics are noticeable only on a time scale much larger than T ωex . Under this approximation, the master equation becomeṡ with the dissipative transition rates It is instructive to compare this master equation to the one in Ref. [52], given in terms of the full Floquet quasienergy states. The key difference here is that the density matrix is projected onto the approximate eigenvectors exp (−iω ex a † at)|φ α rather than onto the exact Floquet solutions. As a consequence of the RWA, the sums in Eq. (33) only include the n = ±1 terms indicating that only one-step transitions are possible where n = −1 refers to emission and n = +1 to absorption. Being consistent with the RWA, we can assume that |ν|, |f |, |ω ex −ω 1 | ≪ ω 1 which yields to |ε α −ε β | ≪ ω ex . Hence, N αβ,+1 is the product of the bath density of states and the bosonic occupation number at temperature T . This corresponds to the thermally activated absorption of a phonon from the bath. On the other hand, N αβ,−1 given in Eq. (29) contains the temperature-independent term .. + J(ω ex ) describing spontaneous emission.

Observable for the nonlinear response
Assuming an ergodic dynamics of the full system, or equivalently that there is just one eigenvector ̺ ∞ of the superoperator S in Eq. (32), corresponding to a vanishing eigenvalue, and that all the other eigenvalues have negative real part, the asymptotic solution of Eq. (32) is To simplify notation, we omit in the following the superscript ∞ but only refer to the stationary state ̺ αβ ≡ ̺ ∞ αβ . We are interested in the mean value of the position operator in the stationary state according to Using Eq. (23) yields X = A cos (ω ex t + ϕ), with the oscillation amplitude and the phase shift with θ being the Heaviside function.
6. Analytical solution for the lineshape of the multiphonon resonance in the perturbative regime When the driving frequency ω ex is varied, the amplitude A shows characteristic multiphonon resonances at those values for which the quasienergy levels form avoided level crossings [34]. While in Ref. [34] these resonances have been studied numerically, it is the central result of this work to calculate their line shape analytically by solving the corresponding master equation in the Van Vleck perturbative regime. Within the limit of validity of the RWA, i.e., |ν|, |f |, |ω ex − ω 1 | ≪ ω 1 , we have |ε α − ε β | ≪ ω ex . In the regime of low temperature k B T ≪ ω ex , it follows from Eq. (29) that N αβ,−1 ≃ J(ω ex ) and N αβ,1 ≃ 0 entering in the transition rates in Eq. (33). This approximation corresponds to consider spontaneous emission only and yields the dissipative transition rates Here, we have defined A αβ ≡ φ α |a|φ β . Note that it is consistent with the previous approximation to set ω ex /ω 1 ≈ 1. Hence, all the following results are valid for Ohmic as well as super-Ohmic baths. In the following we will use this simplified transition rates to solve the master equation near the multiple multiphonon resonances. The transition between the groundstate and the N-phonon state is the narrowest. Hence, it will be affected first when a finite coupling to the bath is considered. In particular, it is interesting to consider the case when the damping constant γ s is larger than the minimal splitting Ω N 0 between the two quasienergy states but smaller than all the minimal splittings of the other, i.e., Ω N 0 < γ s ≪ Ω N n for n ≥ 1. In this case, we can assume a partial secular approximation: We set all the off-diagonal elements to zero except for ̺ 0N and ̺ N 0 = ̺ * 0N . In this regime the stationary solutions are determined by the conditions 0 = β L αα,ββ ̺ ββ + L αα,0N 2 Re(̺ 0N ) , For very weak damping, i.e., when γ s is smaller than all minimal splittings (γ s ≪ Ω N n ), the off-diagonal elements of the density matrix are negligibly small and can be set to zero. Within this approximation, the stationary solution for the density matrix is determined by the simple kinetic equation In this regime, a very simple physical picture arises. The bath causes transitions between different quasienergy states, but here, the transition rates are independent from the quasienergies. It is instructive to express the quasienergy solutions in terms of the harmonic oscillator (HO) solutions as |φ α = n c αn |n with some coefficients c αn . The transition rates between two quasienergy states then read This formula illustrates simple selection rules in this low-temperature regime: Only those components of the two different quasienergy states contribute to the transition rate whose energy differs by one energy quantum (n ↔ n + 1).

One-phonon resonance vs. antiresonance
Before we consider the general multiphonon case, we first elaborate on the one-phonon resonance. This, in particular, allows to make the connection to the standard linear response of a driven damped harmonic oscillator which is resonant at the frequency ω 1 + ν. We will illustrate the mechanism how this resonant behavior is turned into an antiresonant behavior when the damping is reduced (and the driving amplitude f is kept fixed). The corresponding effective HamiltonianH ′ 0 follows from Eq. (15) and is readily diagonalized by the quasienergy states |φ 0 and |φ 1 which are of zero-th order in ε and which are given in Eq. (16). The master equation (39) can be straightforwardly solved in terms of the rates L αβ,α ′ β ′ for which one needs the ingredients A 00 = −A 11 = sin(θ/2) cos(θ/2), A 01 = cos 2 (θ/2) and A 10 = − sin 2 (θ/2). The general solution follows as , Imρ 01 = Ω(∆) L 01,01 − L 01, 10 Reρ 01 , where Ω(∆) = ε 0 − ε 1 . In the following, we calculate the amplitude A according to Eq. (36) to zero-th order in ε. In Fig. 2, we show the nonlinear response for the parameter set (in dimensionless units) f = 10 −5 and ν = 10 −3 . Moreover, the one-phonon resonance condition reads ω ex = ω 1 + ν. The transition from the resonant to antiresonant behavior depends on the ratio γ/Ω 10 = γ/(2f ). For the case of stronger damping γ/(2f ) = 10, we find that the response shows a resonant behavior with a Lorentzian form similar to the response of a damped linear oscillator. In fact, the corresponding standard classical result is also shown in Fig. 2 (black dashed line). The only effect of the nonlinearity to lowest order perturbation theory is to shift the resonance frequency by the nonlinearity parameter ν. The resonant behavior turns into an antiresonant one if the damping constant is decreased to smaller values. A cusp-like line profile arises in the limit of very weak damping when the damping strength is smaller than the minimal splitting, i.e., γ/(2f ) ≪ 1. Then, the response follows from the master equation (40) as This antiresonance lineshape is also shown in Fig. 2 (see dotted-dashed line). At resonance ∆ = 0, we have an equal population of the quasienergy states: ρ 00 = ρ 11 = 1/2 and both add up to a vanishing oscillation amplitude A since A 00 = −A 11 . Note that we show also the solution from the exact master equation containing all orders in ǫ, for the case γ/(2f ) = 0.5 and s = 1 (blue dashed line in Fig. 2), in order to verify the validity of our perturbative treatment.

Multiphonon resonance vs. antiresonance
In this subsection we want to investigate the multiple multiphonon resonances N > 1.
In order to illustrate the physics, we start with the simplest case at resonance and within the secular approximation.
The transition rates between states belonging to the same pair are zero with the exception L (N −1)/2(N −1)/2,(N +1)/2(N +1)/2 = γ s (N + 1)/8. The dynamics can be illustrated with a simple analogy to a double-well potential. Each partner of the pair |φ n and |φ N −n of the quasienergy states consists of a superposition of two harmonic oscillator states |n and |N −n which are the approximate eigenstates of the static anharmonic potential in the regime of weak nonlinearity. In our simple picture, |n and |N − n should be identified with two localized states in the two wells of the quasienergy potential, see Fig. 3 for illustration. Note that a quasipotential can be obtained by writing the RWA Hamiltonian in terms of the two canonically conjugated variables X and P [41]. The right/left well should be identified with the internal/external part of the quasieneryg surface shown in Ref. [41].
In the figure, we have chosen N = 8. Within our analogy, the states |0 , |1 , ..., |N/2−1 are localized in one (here, the left) well, while |N , |N −1 , ..., |N/2+ 1 are localized in the other well (here, the right). The fact that the true quasienergy states are superpositions of the two localized states is illustrated by a horizontal arrow representing tunneling.
From Eq. (44) follows that a bath-induced transition is only possible between states belonging to two different neighboring pairs. As discussed after Eq. (41), the only contribution to the transition rates come from nearby HO states. In our case, we consider only spontaneous emission which corresponds to intrawell transitions induced by the bath. This is shown schematically in Fig. 3 by the vertical arrows with their thickness being proportional to the transition rates. We emphasize that the bathinduced transitions occur towards lower lying HO states. Consequently, in our picture, spontaneous decay happens in the left well downwards but in the right well upwards.
The driving field excites the transition from |0 to |N while the bath generates transitions between HO states towards lower energies according to |N → |N − 1 → ... → |0 when only spontaneous emission is considered.
As a consequence, the ratio of the occupation numbers of two states belonging to two neighboring pairs is simply given by the ratio of the corresponding transition rates according to Hence, the unpaired state |φ N/2 (for N even) or the states |φ (N −1)/2 and |φ (N +1)/2 (for N odd) are the states with the largest occupation probability. By iteration, one finds 6.2.2. Density matrix around the resonance So far, we have discussed the dynamics exactly at resonance. Next, we consider the situation around the resonance and for an increased coupling to the bath. Therefore, we compute the stationary solution using the conditions in Eq. (39) and the general leading order solution for the quasienergy states given in Eq. (16). The expressions for the rates which are modified compared to before follow straightforwardly and are given in the Appendix. Similarly, the only three equations which change compared to the previous situation are also presented there. These equations can be straightforwardly solved by Away from the resonance (|θ| ≪ 1), the density matrix follows as In the limit of strong coupling (γ s ≫ Ω(∆)), one finds for any θ. This nicely illustrates that when the coupling to the bath is strong enough, the possibility of resonant tunneling between |0 and |N is destroyed and a trivial asymptotic state results. This is true even if tunneling transitions between the other states are possible. Moreover, this also shows that moving away from resonance also suppresses multiphonon tunneling transitions. In other words, the only requirement for the multiphonon transition to occur in the stationary limit is the possibility of the tunneling transition |0 → |N .
6.2.3. Lineshape around the resonance Within our partial secular approximation, the lineshape of the oscillator nonlinear response given in Eq. (36) reduces to The leading order is given by the zero-th order expression for ̺ and the first-order expressions for A αα , A N 0 and A 0N . In order to compute these matrix elements, we determine the first order eigenvectors using Van Vleck perturbation theory according to where S 1 is the first order component in the expansion of S with respect to ε given in Eq. (11). The matrix elements of its off-diagonal blocks are given by Here, E α are the eigenenergies of the unperturbed Hamiltonian H 0 given in Eq. (10). This yields for N = 2 The corresponding result for the nonlinear response for N = 2 is shown in Fig. 4 for the case ν = 10 −3 and f = 10 −4 for different values of γ s /Ω 20 . For strong damping γ s /Ω 20 = 5, the resonance is washed out almost completely. Decreasing damping, a resonant lineshape appears whose maximum is shifted compared to the resonance condition ω ex = ω 1 + 3ν/2. Note that the dashed line refers to the result which includes all orders in ε and which follows from the numerical solution of the master equation for an Ohmic bath at temperature T = 0.1T 0 . The picture which arises for the behavior is the following: For weak damping (γ s ≪ Ω 20 ), the equilibrium state is a statistical mixture of quasienergy states. At resonance, the most populated state is |φ 1 which oscillates with a phase difference of −π in comparison with the driving. This is due to the negative sign of A 11 in Eq. (54). Hence, at resonance the overall oscillation of the observable occurs with a phase difference of ϕ = −π. Far away from resonance, the most populated state is |φ 0 , see Eq. (49), which oscillates in phase with the driving. Thus, the overall oscillation occurs in phase, i.e., ϕ = 0. If no off-diagonal element of the density matrix is populated (which is the case for weak damping), the overall phase is either ϕ = 0 or ϕ = −π. Hence, increasing the distance from resonance, the amplitude A has to go through zero yielding a cusp-like line-shape. This implies the existence of a maximum in the response. For slightly larger damping, the finite population of the off-diagonal elements leads to a smearing of the cusp. For larger damping, the resonance is washed out completely, as has been already discussed, see Eq. (50). In this regime, the oscillation is in phase with the driving. By decreasing the damping, the population of the out-of-phase state starts to increase near the resonance resulting in a reduction of the in-phase-phase oscillation and thus producing a minimum of the response. This mechanism is effective for a broad range of parameters including larger N, larger ε, and larger temperature T , as also shown numerically in Refs. [34,38]. Note that a calculation up to first order in ε for the density matrix is required for N odd, since the matrix elements A (N ±1)/2(N ±1)/2 have a zero-th order term, in order that the overall result for A is again of first order in ε. Since one obtains more complicated expressions than before, we omit to present them in their full lengths. In Fig. 5, we show the behavior for N = 3 for various damping constants γ s /Ω 30 for the case f = 0.5 × 10 −4 and ν = 10 −3 . For a large value for γ s /Ω 30 , the resonance is washed out completely. When the damping is decreased, a dip appears which corresponds to an antiresonance.
Decreasing the damping further, the antiresonance turns into a clear resonance. This behavior is opposite to the case N = 1 as discussed above, but similar to the case N = 2.

Conclusions
We have studied the nonlinear response of a vibrating nanomechanical beam to a timedependent periodic driving. Thereby, a static longitudinal compression force is included and the system is investigated close to the Euler buckling instability. There, the fundamental transverse mode dominates the dynamics whose amplitude can be described by an effective single-particle Hamiltonian with a periodically driven anharmonic potential with a quartic nonlinearity. Damping is modeled phenomenologically by a bath of harmonic oscillators. We allow for an Ohmic as well as for a super-Ohmic spectral density and have considered the regime of weak system-bath coupling. In this regime, the dynamics is captured by a Born-Markovian master equation formulated in the frame which rotates with the driving frequency. The pure driven Hamiltonian shows avoided level crossings of the quasienergies which correspond to multiphonon transitions in the resonator. In fact, a transition between a resonant and an antiresonant behavior at the avoided level crossings has been found which depends on the coupling to the bath. Concentrating to driving frequencies around the avoided level crossings, the dynamics can be simplified considerably by restricting to a few quasienergy levels. In order to illustrate the basic principles governing the resonance-antiresonance transition, we investigate the perturbative regime of weak nonlinearity and weak driving strength. Then, Van Vleck perturbation theory allows to calculate the quasienergies and the quasienergy states and an analytic solution of the master equation becomes possible yielding directly the nonlinear response. For the one-phonon case, we find a simple expression for the nonlinear response which displays a Lorentzian resonant behavior for strong damping. Reducing the damping strength, an antiresonance arises. For the multiphonon transitions, first an antiresonance arises when damping is reduced. For even smaller values of the damping constant, the antiresonance turns into a resonant peak. This is due to a subtle interplay of varying populations of quasienergy states which is affected by the bath.
Finally, a comment on the observability of this effect is in order. The amplitude A measuring the nonlinear response is of the order of the oscillator length scale x 0 . This makes it challenging to measure the effect directly since the deterministic vibrations are on the same length scale like the quantum fluctuations. In turn, more subtle detection strategies have to be worked out, for instance, the capacitive coupling of the resonator to a Cooper pair box [54,55] or single electron transistors [56,57], the use of squeezed states in this setup [55], the use of a second nanoresonator as a transducer for the phonon number in the first one [37], or the coherent signal amplification by stochastic resonance [36]. In any case, the experimental confirmation of the theoretically predicted effects remains to be provided.

Acknowledgments
This work has been supported by the DFG-SFB/TR 12.
Appendix: Density matrix around the multiphonon resonance For completeness, we present in this Appendix the calculation of the density matrix around the multiphonon resonance which is required for Section 6.2.2. The expressions for the rates which are modified compared to before are readily calculated to be