Semi-classical dynamics of nano-electromechanical systems
aa r X i v : . [ c ond - m a t . m e s - h a ll ] J un Semi–classical dynamics of nano–electromechanical systems
R. Hussein, ∗ A. Metelmann, P. Zedler, and T. Brandes
Institut f¨ur Theoretische Physik, TU Berlin, Hardenbergstr. 36, Berlin D-10623, Germany (Dated: November 19, 2018)We investigate the dynamics of a single phonon (oscillator) mode linearly coupled to an electronicfew–level system in contact with external particle reservoirs (leads). A stationary electronic currentthrough the system generates non–trivial dynamical behaviour of the oscillator. Using Feynman–Vernon influence functional theory, we derive a Langevin equation for the oscillator trajectory thatis non–perturbative in the system–leads coupling and from which we extract effective oscillatorpotentials and friction coefficients. For the two simplest cases of a single and two coupled electroniclevels, we discuss various regimes of the oscillator dynamics.
PACS numbers: 71.38.-k, 73.21.La, 85.85.+j
I. INTRODUCTION
Nano–electromechanical systems (NEMS) are test–beds for the observation of fundamental quantum be-haviour of objects which are huge on the scale of indi-vidual atoms. For example, recent experiments have al-lowed a detailed study and control of single phonons bycooling a macroscopic resonator mode close to its groundstate and coupling it to single electronic degrees of free-dom.One fascinating aspect of NEMS is their concep-tual simplicity that nevertheless can give rise to highlycomplex physics, and the links that can be estab-lished to other fields such as molecular electronics oroptomechanics . One of the challenges are the de-tails of the oscillator–electron coupling, i.e. in the lan-guage of measurement theory to understand, utilize andcontrol the ‘back–action’ effects of the detector (e.g.,superconducting single–electron transistors ) onto theoscillator.In many cases, even if no further approximations insimple theoretical models are made, the coupling toexternal electronic reservoirs that link the detector tothe outer world is treated perturbatively, i.e., in theframework of (quantum) Master equations. This hasturned out to be a highly successful approach, in par-ticular to describe such various systems as NEMS cou-pled to single electron transistors (SETS) , Franck–Condon blockades in transport through moleculeswith strong electron–phonon coupling, or quantumshuttles . From the theory of electronic transportthrough nanostructures , however, it is known thatsuch approximations usually are reliable only in thelimit of high external voltage bias, where non–Markovianeffects due to quantum coherences between the exter-nal reservoirs and the electronic system (SET, quantumdot etc.) can be neglected. It is therefore desirable to de-velop tools that allow a description of NEMS beyond theMaster equation regime (weak electron–leads coupling)and at the same time are not merely perturbative in thecoupling of the oscillator to the electronic environment .In the past, the coupling of electrons to a single bosonicmode has been solved exactly for the case where only one single electron is present , i.e. in an empty band ap-proximation. The inclusion of Fermi sea reservoirs at dif-ferent chemical potential transforms this into a difficultmany–body problem out of equilibrium, and approxima-tions are necessary .Our approach in this paper is to combine exact solu-tions of the electronic system with a semi–classical ex-pansion, together with an adiabatic approximation forthe oscillator dynamics within the Feynman–Vernon in-fluence functional (double path integral) theory . Werevise this method, which has first been used for simpleNEMS models by Mozyrsky and co–workers , and ex-tend it to allow for the description of a relatively largeclass of non–equilibrium electronic environments. Thekey idea is a systematic expansion around the classicalpath in order to obtain a Langevin equation for the oscil-lator. Already at the simplest level of this approximation(neglecting quadratic fluctuations around the diagonalpath in the reduced density matrix of the oscillator), thecoupling to the electronic non–equilibrium environmentgives rise to non–trivial effects such as effective oscilla-tor potentials and non–linear friction coefficients lead-ing to both positive and negative damping . Gaussianfluctuations around the classical path are built into thetheoretical description, but they have to be evaluated bynumerical solutions of the underlying Langevin equationswhich is not done in this paper.We compare two non–interacting electronic ‘quantumdot’ models with one and two levels between source anddrain reservoirs: a single dot, and two dots in series. Theoscillator couples linearly to the dot occupation (singledot) or to the occupation difference (double dot). Oneparticular feature of the double dot case (where quan-tum superpositions of the electrons become important)is the occurence of limit cycles in phase space caused bya negative damping.The paper is organized as follows: after introducingthe path integral formalism with a generic model in Sec.II, we present the single dot case in Sec. III and thedouble dot case in Sec. IV. Detailed derivations of theimportant formulae can be found in the appendices. II. GENERIC MODEL
A large class of NEMS can be described as a compo-sition of an electronic system H e , a mechanical system H osc and a linear coupling between the two. Thus we setup a generic Hamiltonian H gen by H gen = H e + H osc − ˆ F ˆ q, (1) H osc = 12 m ˆ p + V osc (ˆ q ) . (2)Here, H osc describes a single oscillator with ˆ p (momen-tum) and ˆ q (position) operators. The oscillator modeis confined in a potential V osc ( q ), ˆ F denotes an elec-tronic force operator, and m labels the oscillator mass.Whithin this paper the reduced Planck constant is set toone ( ~ = 1).Our generic model does not include an additional os-cillator damping mechanism . In the usual Master equa-tion treatment of NEMS, Lindblad-form damping due toexternal degrees of freedoms is included phenomenologi-cally. In the path integral formalism used here, such de-grees of freedom can be easily included at least for linearor weak coupling to the oscillator. In order to elucidatethe effect of the electronic environment that we treat inall orders in the coupling to external electronic reservoirs(contained in H e ), we choose not to include additionaldamping terms in our model here. A. Stochastic equation of motion
We describe the oscillator dynamics by the reduceddensity matrix of the oscillator in position representation ρ osc ( q, q ′ , t ) = h q | ρ osc ( t ) | q ′ i , for which we derive a semi-classical equation of motion for the oscillator position byusing Feynman–Vernon influence functional theory sim-ilar to Mozyrsky and co-workers . Assuming that thetotal density matrix χ ( t ) factorises at the initial time t into a system and a bath part χ ( t ) = ρ osc ( t ) ⊗ ρ B , thepropagation of the reduced oscillator density h q | ρ osc ( t ) | q ′ i at time t is given by a double path integral , cf. ap-pendix A. A transformation to center–of–mass and rela-tive coordinates x t = q t + q ′ t , y t = q t − q ′ t (3)has the notion to detach the classical trajectory x t fromthe quantum mechanical deviations. Within a Born–Oppenheimer approximation, change of variables allowsus to study a slow oscillator by an adiabatic approxima-tion of the classical trajectory x t ≈ x + t ˙ x , (4)cf. appendix A. In this approach, the typical timescaleof the oscillator movement is slow compared to the elec-tronic transition rates. In the subsequent, when having introduced the angular oscillator frequency by ω andelectron transition rates by Γ L , Γ R we have to satisfy thecondition ω ≪ Γ L , Γ R .In the next step, we derive a stochastic equation ofmotion for the classical trajectory, taking into accountthe propagation of the initial reduced density matrix,cf. appendix A. The key step here is a cluster expansionto quadratic order in the off-diagonal path y t that de-scribes the Gaussian fluctuations around the classical os-cillator trajectory, where the fluctuations are determinedby the properties of the non-equilibrum environment. Toachieve a self–consistent equation of motion, we then re-insert the full time–dependence of the fixed classical tra-jectory in accordance with the adiabatic approximationand end up with the Langevin equation m ¨ x t + V ′ osc ( x t ) − h ˜ F [ x ]( t ) i + ˙ x t A [ x ]( t ) = ξ t . (5)Here the interaction picture of the electronic operatoris given by ˜ F [ x ]( t ) = exp[ i ( H e − ˆ F x ) t ] ˆ F exp[ − i ( H e − ˆ F x ) t ], and ξ t is a stochastic force with zero mean andthe correlation function h ξ t ξ t ′ i = 2 Re h δ ˜ F [ x ]( t ) δ ˜ F [ x ]( t ′ ) i . (6)The fluctuation of the electronic operator is defined by δ ˜ F [ x ]( t ) = ˜ F [ x ]( t ) − h ˜ F [ x ]( t ) i . The friction A [ x ] is givenby A [ x ]( t ) = 2 Z tt dt ′ t ′ Im h δ ˜ F [ x ]( t ) δ ˜ F [ x ]( t ′ ) i . (7)In the following we choose t = −∞ as initial time. Forthe specific cases of single and double dots, we checkedthat the upper integration boundary can be extended toinfinity. III. ANDERSON–HOLSTEIN MODEL (AHM)
The AHM combines a single bosonic mode with a sim-ple electronic transport system. We describe the bosonicpart in first quantisation as H osc = 12 m ˆ p + 12 mω ˆ x (8)with the bosonic position and momentum operators ˆ x and ˆ p , the phonon frequency ω and the phonon mass m .The oscillator length is defined by l ≡ [ mω ] / . Theelectronic part is a single dot level confined between twoleads: H e = ε d ˆ d † ˆ d + X kα ε kα ˆ c † kα ˆ c kα + X kα (cid:0) V kα ˆ c † kα ˆ d + H. c. (cid:1) . The dot level has energy ε d and creation/annihilationoperators d † / d . The operators ˆ c † kα /ˆ c kα create/anihilate occupation zeroaverage occupationoccupation unity x = 0 x = g Γ L Γ l x = gl µ L µ R V Bias F Hooke F e FIG. 1: Sketch of the AHM model. Two leads (grey bars) withchemical potentials µ L , µ R = µ L − eV Bias embed the dot levelwhich couples to an oscillator. Here V Bias is the bias voltageand e the electron charge. The level energy ε d is shifted bythe oscillator position x , the shifted level energy ε x = ε d − λx is stationary, if the oscillator force F Hooke = mω x and theelectron force F e = λ h ˆ n [ x ] i are in balance. Here m and ω describe the oscillator mass and angular frequency, and λ (g)denotes the coupling constant, cf. eq.(10). For sufficientlysmall tunneling rates Γ L , Γ R the stable stationary level en-ergies (dashed lines) are located at x = 0 , gl Γ L / Γ , gl . Theinstable points of ε x (dotted lines) are located around thechemical potentials. electrons with momentum k and energy ε ka in a free elec-tron gas in the left ( α = L) or right ( α = R) lead. Elec-tronic transitions are possible with amplitude V ka be-tween the dot and a state in the lead. Between the twosubsystems there is a simple linear coupling with couplingconstant λ , such that the total Hamiltonian reads H AHM = H e + H osc − λ ˆ d † ˆ d ˆ x. (9)For convenience we introduce the dimensionless couplingconstant: g = λmω l . (10)Here, we regard spinless electrons. Note that the general-ization to a model including both spin directions requiresan onsite interaction term like U ˆ n ↑ ˆ n ↓ . In the following,we only consider non–interacting electrons. A physicalrealization of the model discussed here would correspondto either spin polarized electrons or a coupling to theoscillator that effects only electrons of one certain spindirection. A. Langevin equation
When applying our generic form (1) to the AHM, weobtain the Langevin equation m ¨ x t − F eff ( x t ) + ˙ x t A [ x ]( t ) = ξ t (11) PSfrag replacements Position h x i /l B i a s v o l t ag e V B i a s / ω F Hookes F e U eff × − − g Γ L / Γ g g Γ L / Γ g Γ L / Γ g Γ L / Γ ggg FIG. 2: left ) Density plot of the effective oscillator poten-tial U eff as a function of oscillator position h x i and V Bias in units of ω and l with the parameters Γ L , R /ω = 0 . ε d /ω = 2 . βω = 10 and g = 2 .
4. For increasing biasvoltage, the effective potential shows two, three, two and oneminima. This different regions are separated by white linesat V Bias /ω = 2 . / . / . Right ) U eff and the two con-tributions F Hooke , F e to the force F eff at bias voltages (sym-metric choice) V Bias /ω = 0 . / . / . U eff exhibits extremawhere the two force contributions (grey and dashed line) arein balance. The width of the center plateau in the electronicforce contribution F e ∝ h n i grows with increasing bias volt-age V Bias , cf. Fig. (1) for the intermediate occupation h n i . Forsufficiently high bias voltage only one minimum remains. with the effective force F eff and the friction A [ x ]: F eff ( x t ) = − mω x t + λ h ˜ n ( t ) i = − ω l (cid:20) xl − g h ˜ n ( t ) i (cid:21) ,A [ x ]( t ) = 2 λ Z tt dt ′ t ′ Im h δ ˜ n ( t ) δ ˜ n ( t ′ ) i . (12)The effective force has two contributions; a term propor-tional to the elongation (Hooke’s law) and the electronforce which is proportional to the dot occupation; whichonly contributes if the dot is occupied. The friction term x t A [ x ] results from stochastic electron jumps between theleads and the dot.For finite bias voltage, figure 1 shows the positions of thedot energy level and the points of instable balance, whichresult from the balance of both forces.The effective force and therewith the oscillator poten-tial are determined by the dot occupation n ( t ), whereasthe friction and the stochastic force correlation h ξ t ξ t ′ i = λ h δ ˜ n ( t ) δ ˜ n ( t ′ ) i depend on the imaginary/real partof the dot correlation function. The dot correlation func-tion in terms of the lesser and greater Green’s function(which are derived in appendix B) reads h δ ˜ n ( t ) δ ˜ n (0) i = G < ( − t ) G > ( t ) . (13) B. Effective potential
The occupation of the dot is calculated in an adia-batic approach with the help of the lesser Green’s func-tion G < ( ω ), cf. appendix B, where we assume constanttunneling rates Γ α . For finite temperatures, the dot oc-cupation reads h ˜ n ( t ) i = − i π Z dω G < ( ω )= 12 − π X α Γ α Γ Im Ψ (cid:18)
12 + β Γ4 π + i β ( ε x − µ α )2 π (cid:19) , (14)whereby β stands for the inverse temperature, ε x is ashort notation for ε d − λ ˆ x , µ α denotes the chemical pote-tials and Γ ≡ Γ L + Γ R . Ψ designates the Digamma func-tion. By integration we obtain the effective potential U eff ( x ) = − Z x dx ′ F eff ( x ′ ) = (15)2 β X α Γ α Γ Re (cid:20) ln Γ (cid:18) ξ + i β ( ε d − λx − µ α )2 π (cid:19) − ln Γ (cid:18) ξ + i β ( ε d − µ α )2 π (cid:19)(cid:21)(cid:12)(cid:12)(cid:12)(cid:12) ξ = + β Γ4 π + 12 xl (cid:20) xl − g (cid:21) ω , (16)with Γ( · ) denoting the Gamma function. In the absenceof friction the dynamics of the oscillator is determinedby the effective potential U eff ( x ). The high temperaturecase is of minor interest, because the temperature washesout the structures of U eff ( x ).In Fig. 2 we present features of the effective potential U eff ( x ) at zero temperature. U eff ( x ) has minima whenthe effective force F eff ( x ) is zero. This is the case whenthe oscillator force and the electron force are in balance.We plot the two contributions to the effective force: theforce F Hooke is proportional to the displacement h x i , thecontribution F e results from the dot occupation (scaledwith g ) and has two steps at x α l = 1 g (cid:20) ε d ω − µ α ω (cid:21) , α ∈ { L , R } . (17)In the upper part of Fig. 2 we easily recognize that forlarge bias voltage there will be only one minimum at x/l = g Γ L / Γ. By changing the bias voltage or the cou-pling strength, we can reach situations where two minimaat x/l = 0 and at x/l = g are added, like in the middlepart of Fig. 2 and also situations, where only the twominima at the side remain and the one in the middlevanishes, like in the bottom part of the 2 (this agreeswith Mozyrsky et al. ). Increasing bias voltage shiftsthe steps in F e apart, whereas increasing the couplingconstant minimises the distance. The positions of theminima/steps are exact in the zero rate limit (Γ L → R →
0) where the averaged occupation is step–like. Forfinite rates, the steps smoothes out. The form of the oscillator potential already hints to-wards the phase space spanned by position h x i and mo-mentum m h ˙ x i . In the absence of friction, the minima ofthe potential correspond to the fixed points of the oscil-lator motion. C. Phase space portrait
The classical trajectory h x t i is obtained by neglectingfluctuations due to the stochastic force ξ t in (11); there-fore the equation of motion reads m h ¨ x t i + h ˙ x t i A [ h x i ]( t ) − F eff ( h x t i ) = 0 . (18)In Fig. 3 we show the effective potential and the fric-tion, as well as the solutions of the classical equation ofmotion in phase space without and with friction. Therows correspond to three different voltages leading to thethree typical cases with one, three and two minima ofthe potential (the latter case is investigated in ). Theinitial conditions are chosen for each phase diagram sep-arately in order to make the characteristic shapes visible.The phase diagrams without friction follow directly fromthe shape of the effective potential. The trajectories in-cluding the friction follow the ones without friction fora while until a position with a peak in the friction isreached that produces a kink-like damping feature. Thefriction causing the kinks is maximal, when the shiftedenergy level ε x is in resonance with the chemical poten-tials µ L , R = ± V Bias / V Bias /ω = 0 . PSfrag replacements
Position h x i /l F r i c t i o n A l / ω P o t e n t i a l U e ff / ω M o m e n t u m h ˙ x i / [ ω l ] M o m e n t u m h ˙ x i / [ ω l ] A) B) C) D) − − − − − − − FIG. 3: Effective oscillator potentials (A), frictions A ( x ) (B) and the phase space portraits without (C) and with (D) frictionin units of ω and l with the parameters Γ L , R /ω = 0 . ε d /ω = 3 . g = 2 .
45 at zero temperature. From the bottom tothe top the values of the bias voltage (symmetric choice) read V Bias /ω = 0 . / . / . Right ). The peaks of the friction are located at x = [ ε d ∓ V Bias / /λ where the shifted level energies ε x are in resonancewith the chemical potentials µ L , R , cf. Fig. 1. IV. DOUBLE QUANTUM DOT SYSTEM (DQD)
The DQD that we treat in this section consits of twosingle dot levels coupled by a tunnel barrier. Again weassume a coupling to a single bosonic mode. The totalHamiltonian is composed of the oscillator part H osc , cf.eq. (8), the electronic part H e and an interaction partwhich describes the coupling between the oscillator andthe two dots. In contrast to the AHM, here the oscillatorcouples to the difference of the occupation numbers withthe coupling strength λ . The total Hamiltonian thereforereads H DQD = H e + H osc − λx ( ˆ d † L ˆ d L − ˆ d † R ˆ d R ) , (19)containing the electronic part H e = X kα ε kα ˆ c † kα ˆ c kα + X α ν α ˆ d † α ˆ d α + T c ˆ d † L ˆ d R + T ∗ c ˆ d † R ˆ d L + X kα (cid:2) V kα ˆ c † kα ˆ d α + V ∗ kα ˆ d † α ˆ c kα (cid:3) . (20) In this model, ˆ d † α / ˆ d α denotes the creation/annihilationoperator of the α th ( α ∈ { L , R } ) dot. ν α designates thecorresponding level energy and T c describes the tunnelcoupling matrix element between the two dots. Note thatwe include no Coulomb interaction terms here.The stationary average of the population difference h σ z i corresponds to the electronic force operator ˆ F ofthe generic model, eq. (1). The occupation of the α thdot h ˆ n α i = − i π Z dω G <α,α ( ω ) (21)can be calculated by using Keldysh’s equation (C5), seesection C 1. Therewith the occupation difference of the m o m en t u m < p > / [ l ] position
4; 0 .
43; 0 .
49; 1 . ω . We used equal ratesΓ L = Γ R = 1 . ω . For the chemical potentials we assumed µ L = 1 ~ ω and µ R = − ~ ω . The dimensionless coupling constantis chosen as g = 2 .
5, and the internal bias voltage as V int = 5 ~ ω , wheras ν L = − ν R = eV int / left/right dot follows from h σ z i = h n L i − h n R i = Γ4 π " µ L Z −∞ dω ( ω − e ν R ) + (cid:0) Γ4 (cid:1) − | T C | ω + 2 Aω + B − µ R Z −∞ dω ( ω − e ν L ) + (cid:0) Γ4 (cid:1) − | T C | ω + 2 Aω + B , (22)with the abbreviations ˜ ν L , R = ν L , R ∓ λx and A = − (cid:2) | T c | + ˜ ν − (Γ / (cid:3) , B = (cid:2) | T c | + ˜ ν + (Γ / (cid:3) , (23)whereas we assumed e ν R = − e ν L . Calculating the integrals leads to h σ z i = Γ4 π X α ∈ L , R " e ν L I ( µ α , A, B )+ sgn( ν α ν L ) (cid:26) I ( µ α , A, B ) + (cid:16)e ν + (Γ / − | T C | (cid:17) I ( µ α , A, B ) (cid:27) . (24)They are expressed in terms of the auxiliary functions I j defined in appendix D. A. Langevin equation
When applying the generic model, cf. eq. (5), to theDQD Hamiltonian H DQD , we obtain the Langevin equa-tion m ¨ x t − F eff ( x t ) + ˙ x t A [ x ]( t ) = ξ t (25)with the effective force F eff and the friction A [ x ]( t ). Incontrast to the AHM, the effective force is affected bythe population difference and not only by the occupationnumber. Explicitly, F eff ( x t ) = − mω x t + λ h ˜ σ z ( t ) i , (26) A [ x ]( t ) = 2 λ Z tt dt ′ t ′ Im h δ ˜ σ z ( t ) δ ˜ σ z ( t ′ ) i . (27)In appendix C 2 we derive the explicit expression for thefriction. In the case of infinite bias (IB) and e ν R = − e ν L we obtain A IB [ x ]( t ) = − | T c | λ Γ e ν L e ν + | T c | + 5 Γ (cid:16)e ν + | T c | + Γ (cid:17) . (28)In contrast to the AHM the friction does not disappearfor infinite bias. The second difference that we can recog-nise by regarding the prefactor e ν L is that we obtain re-gions where the friction is negative . The latter also holdsfor the finite bias case.The real part of the correlation function according tothe population difference σ z = d † L d L − d † R d R determinesthe correlation function of the stochastic force h ξ t ξ t ′ i = λ h δ ˜ σ z [ x ]( t ) δ ˜ σ z [ x ]( t ′ ) i , (29)with fluctuation δ ˜ σ z [ x ]( t ) = ˜ σ z [ x ]( t ) − h ˜ σ z [ x ]( t ) i . B. Fixed point analysis
The effective potential determines the behavior of theoscillator trajectories in the phase space ( h p i – h x i –plane),as seen in section III B for the Anderson Holstein model.In the following we examine the differential equation ofthe system by studying its fixed points . Some fur-ther theoretical details of this kind of investigation areexplained in appendix C 3.For the double dot we obtain the dynamical system h ˙ x t i = 1 m h p t ih ˙ p t i = ω l (cid:20) − h x t i l + g h σ z ( t ) i − h p t i ω l l A [ h x i ]( t ) (cid:21) . (30)Fixed points occur under the condition h ˙ p t i = h ˙ x t i = 0,i.e. h p t i = 0 and following from that, the fixed pointsposition coordinates are equal to the roots of the effectiveforce F eff ( h x t i ) = − h x t i l + g h σ z ( t ) i = 0 . The Jacobian matrix is obtained from J ∗ = ω l h − l + g ∂∂ h x t i h σ z ( t ) i (cid:12)(cid:12) h x ∗ i i − A [ h x ∗ i ] ! evaluated at the fixed point h x ∗ i , whereby h p ∗ i = 0.Determinant ∆ and trace τ become∆ = ω l (cid:20) l − g ∂∂ h x t i h σ z ( t ) i (cid:12)(cid:12) h x ∗ i (cid:21) , τ = − A [ h x ∗ i ] . (31)The trace decides about the stability of a fixed pointand is equal to the negative friction here. For the casewithout friction the trace is equal to zero, therefore onlycenters occur in the phase plane. In the case with frictionthe trace can be either positive or negative leading toboth, stable and unstable fixed points. For comparisonsee Figure 4, where row A depicts the results of a fixedpoint analysis.The effective force F eff is plotted in the upper partof each plot in row A together with trace and determi-nant. Therewith the characteristics of the fixed pointsare determined and it is possible to predict the shape ofthe phase space portrait. In the lower parts of the plotsin row A these predictions are illustrated, fixed pointsare marked by black lines. For small tunnel coupling T c (diagram A1) we obtain seven fixed points. These canbe characterised as three stable and one unstable spirals,each seperated by one of the three saddle points. Sta-ble spirals correspond to the different rest positions forthe oscillator and the saddle points to the points of in-stable balance. Increasing the tunnel coupling leads tofive fixed points. In graph A2 we obtain a stable spiralin the middle enclosed by two saddle points and followedby a stable spiral on each side. In the diagrams A3 andA4 the mean point changes to an unstable spiral. Withfurther increased tunnel coupling (A3) there remain onlythree fixed points and for even higher values of T c (A4)the number of zeros in the effective force reduces to one.For the oscillator the latter means that it is shifted to anew rest position independent from its inital position. C. Phase space portraits
In Figure 4 the rows B and C depict phase space por-traits for the double dot system without and with fric-tion. The four columns correspond to four different val-ues of the tunnel coupling T c . As initial condition, themomentum was set to zero and the positions were chosenin order to show the various shapes of the trajectories.In the case without friction we recognize periodic cycleswhich are stable and run around one or more fixed points.These centers were also expected from the analysis in sec-tion IV B. The fixed points correspond to certain statesof the double dot system. If the left dot is occupiedthe rest position of the oscillator is shifted to the rightand correspondingly to the left for an occupied right dot.These points turn to stable spirals when the friction isturned on. The states when both dots are occupied orempty correspond to the fixed point in the middle.The lowermost row shows what happens when we in-clude the friction in our calculations. In contrast to thesingle dot system, here the friction has positive as wellas negative values depending on the position of the os-cillator. This means that the oscillator is either deceler-ated or accelerated. Both can be interpreted as inelasticjumps of the electrons, where energy is transfered be-tween the electrons and the oscillator in both directionslike it has been observed in . There, the authors considera resonator coupled to a superconducting single electrontransistor (SSET). As a result of the interplay of positiveand negative damping in a certain parameter range theyobserved limit cycles and bistability in the phase plane,in our work we obtain a likewise behaviour for the os-cillator. In contrast to our semiclassical approximationthey investigate the Wigner function of the system withnumerical master equations. They compare these resultsto a mean field evaluation of the expectation value ofthe oscillator position . For weak coupling the meanfield approach gives quantitatively correct results, andfor higher coupling it still describes the dynamics qual-itatively correct. These results suggest that our use ofaverage oscillator positions and momenta is qualitativelycorrect for a description of the oscillator instabilities.Consider again row C in Figure 4, where the resultsfor the dynamics of the oscillator with friction are plot-ted. The outer left and right stable spiral do not changeby increasing the tunnel coupling T c . By contrast theoscillator’s behaviour between these stable rest positionschanges a lot. In graph C1 we observe a stable and anunstable spiral, like we expect from the fixed point anal-ysis. In the neighbourhood of the unstable fixed pointthe friction is negative, so the oscillator trajectory is re-pelled and ends up in the right stable spiral. In diagramC2 we recognize that the latter path becomes stable. Thelimit cycle appears when the unstable spiral has disap-peared and exists as long as the tunneling coupling isin the range of 0 . ≤ | T c | /ω ≤ .
5. There we ob-serve a bistability: as the initial position gets closer tothe fixed point the limit cycle turns into a stable spi-ral. By further increasing T c , the middle spiral becomesunstable and a second limit cycle appears (C3). Thislimit cycle with a smaller radius exists in the range of0 . ≤ | T c | /ω ≤ .
54. For | T c | /ω ≃ .
49 the systemundergoes a Hopf bifurcation , which happens when apair of complex eigenvalues from the dynamical system,which determine the evolution in the phase plane, seesection C 3, cross the imaginary axis from the left to theright half-plane. In other words, the trace τ changesits sign and at the bifurcation point the eigenvalues arepurely imaginary λ , = ∓ i √ ∆, see equation (C16). Inthe last graph of row C both limit cycles have disappearedand the oscillators path ends up in the left or right stablepoint.In our calculations we choose Γ = 3 ω , standing insome contrast to our adiabatic approach, which impliesa slow oscillator (Γ ≫ ω ). If Γ ∼ ω the interactionbetween the current and the oscillator is strongest , be-cause both act on the same timescale. Interesting effectsstill appear with a slightly enlarged Γ like in our plots,but for Γ ≫ ω there remains only one stable spiral. This means, that the oscillator rest position is shifted from itsground position caused by the stochastic processes initi-ated by the current. We presume that our approach isuseful also for a comparative fast oscillator and we willaccomplish further investigations with a non-adiabaticapproximation to reconsider our results. V. SUMMARY
We have derived a stochastic equation of motion thatdescribes the dynamics of a single oscillator coupled toan electronic environment out of equilibrium. We studiedtwo cases, namely the single dot level and double dot(two–level) electronic system. For both cases we haveexplained the features of effective potential and frictionfor the ensemble averaged oscillator motion. The effectswe recognize fit well together with former works. In theDQD model limit cycles and bistabilities appear.Until now the master equation has been used for mostinvestigations of the oscillator behaviour in NEMS. Wehave used a method that gives us access to regions wherethe master equation has problems: we naturally includefinite bias, and arbitrary electron coupling to externalreservoirs.We had to stay in a regime with a relatively fast oscil-lator in order not to miss the interesting physical effects.The validity of our method in this regime could still beimproved with a non-adiabatic calculation.
VI. ACKNOWLEDGEMENTS
This work was supported by project DFG BR 1528/5-2, the WE Heraeus foundation and the Rosa Luxemburgfoundation. We acknowledge discussions with C. Emary.
Appendix A: Path integral representation of thereduced density matrix
The reduced density matrix ρ ( t ) = tr B χ ( t ), obtainedfrom the total density matrix χ ( t ) by tracing out thebath degrees of freedom, describes the mechanic sub-system H osc , cf. (1). By using a factorising initialcondition χ ( t ) = ρ ( t ) ⊗ ρ B the elements of the re-duced density matrix can be expressed in a path integralrepresentation h q | ρ osc ( t ) | q ′ i = Z dq dq ′ h q | ρ osc ( t ) | q ′ i× Z q ( t ) q ( t ) D q ( τ ) Z q ′ ( t ) q ′ ( t ) D ∗ q ′ ( τ ) e i ( S q − S q ′ ) F [ q, q ′ ]( τ ) . (A1)Hereby S q = R tt dt ′ (cid:2) m ˙ q t ′ − V osc ( q t ′ ) (cid:3) denotes the clas-sical action in the path integral for eq. (A1) and theFeynman–Vernon influence functional is defined by F [ q, q ′ ]( t , t ) = tr B (cid:8) ˆ U † [ q ′ ]( t, t ) ˆ U [ q ]( t, t ) ρ B (cid:9) . (A2)The time–evolution operators are defined byˆ U [ q ]( t, t ) = T ← exp (cid:2) − i Z tt dt ′ H res [ q ]( t ′ ) (cid:3) (A3)corresponding to the reservoirs H res [ q ]( t ) = H e − ˆ F q t . Inthe next step we establish an effective interaction pictureby ˜ U [ q ]( t, t ) = e i H ( t − t ) ˆ U [ q ]( t, t )= T ← exp (cid:2) − i Z tt dt ′ ˜ V [ q ]( t ′ ) (cid:3) , ˜ V [ q ]( t ) = e i H ( t − t ) V [ q ]( t ) e − i H ( t − t ) . (A4)Hereby the reservoir is decomposed in an unperturbedpart H = H e − ˆ F x and the perturbation V [ q ]( t ) = − ˆ F ( t ˙ x + y t ) and V [ q ′ ]( t ) = − ˆ F ( t ˙ x − y t ), where wewrote x t ≈ x + t ˙ x using the adiabatic approximation(4).By using eq. (A4), the influence functional can be ex-panded to second order in perturbation theory. Denoting h · i = tr B { ρ B · } , ˜ V ′ := ˜ V [ q ′ ]( t ) and ˜ V := ˜ V [ q ]( t ), theFeynman–Vernon influence functional reads F [ q, q ′ ]( t , t ) = h ˜ U [ q ′ ]( t, t ) ˜ U [ q ]( t, t ) i = h (cid:2) T → e i R tt dt ˜ V ′ (cid:3)(cid:2) T ← e − i R tt dt ˜ V (cid:3) i = 1 + i Z tt dt h ˜ V ′ − ˜ V i + Z tt dt Z t t dt h ( ˜ V ′ − ˜ V ) ˜ V − ˜ V ′ ( ˜ V ′ − ˜ V ) i . (A5)When plugging in the definition of the perturbation (in-teraction picture) and taking into account that the forceoperator ˆ F is hermitian, we obtain F [ x + y/ , x − y/ t , t ) = 1 + i Z tt dt y t × (cid:20) h ˜ F ( t ) i − x Z t t dt t Im h ˜ F ( t ) ˜ F ( t ) i (cid:21) − Z tt dt Z t t dt y t Re h ˜ F ( t ) ˜ F ( t ) i y t . (A6)Performing a cluster expansion , the influence func-tional can be expressed in terms of the influence phaseΦ[ x, y ]( t , t ) ≡− ln F [ x + y/ , x − y/ t , t ) = − i Z tt dt y t × (cid:20) h ˜ F ( t ) i − x Z t t dt t Im h δ ˜ F ( t ) δ ˜ F ( t ) i (cid:21) + Z tt dt Z t t dt y t Re h δ ˜ F ( t ) δ ˜ F ( t ) i y t , (A7) whereby δ ˜ F ( t ) = ˜ F ( t ) − h ˜ F ( t ) i denotes the fluctuationaround h ˜ F ( t ) i . This equation can be easily verified by ex-panding the exponential exp[ − Φ[ x, y ]( t , t )] to second or-der in perturbation theory and comparing with eq. (A6).Together with the classical action, we define an effec-tive action functional by A [ x, y ]( t , t ) := S x + y/ − S x − y/ + i Φ[ x, y ]( t , t )= Z tt dt (cid:2) m ˙ x ˙ y − V osc ( x + y V osc ( x − y (cid:3) + i Φ[ x, y ]( t , t ) , (A8)whereby the reduced density matrix expressed in termsof the new variables, eq. (3), reads h x + y | ρ ( t ) | x − y i = Z dx dy h x + y | ρ ( t ) | x − y i× Z x ( t ) x ( t ) D x ( τ ) Z y ( t ) y ( t ) D ∗ y ( τ ) e i A [ x,y ]( t ,t ) . (A9)In a semiclassical approach we assume small deviationsof the off–diagonal trajectories y from the diagonal ones.Thus the potential difference leads in second order in y to V osc ( x + y − V osc ( x − y V ′ osc ( x ) + O ( y ) . (A10)This approximation is exact for quadratic potentials, aswe treat in this paper. Furthermore, with the boundaryconditions y ( t ) = y ( t ) = 0 and integration by parts theeffective action functional is quadratic in y , A [ x, y ]( t , t ) = − Z tt dt y (cid:20) m ¨ x + V ′ osc ( x ) − h ˜ F ( t ) i + 2 ˙ x Z t t dt t Im h δ ˜ F ( t ) δ ˜ F ( t ) i (cid:21) + i Z tt dt Z t t dt y Re h δ ˜ F ( t ) δ ˜ F ( t ) i y ] ≡ − Z tt dt y K [ x ] + i Z tt dt Z t t dt y L , [ x ] y . (A11)Completing the square of the integral kernel R D ∗ y exp { i A [ x, y ] } , the resulting path integral de-scribes a stochastic process with Langevin equation K t [ x ] = m ¨ x t + V ′ osc ( x t ) − h ˜ F ( t ) i + 2 ˙ x Z tt dt ′ t ′ Im h δ ˜ F ( t ) δ ˜ F ( t ′ ) i = ξ t (A12)with ξ t a Gaussian stochastic force. To obtain a selfcon-sistent equation of motion, we finally replace x by x t .This substitution concerns also the unperturbed Hamil-tonian which leads to H → H e − ˆ F x t . The Langevin0equation then reads m ¨ x t + V ′ osc ( x t ) − h ˜ F ( t ) i + 2 ˙ x t Z tt dt ′ t ′ Im h δ ˜ F ( t ) δ ˜ F ( t ′ ) i = ξ t . (A13)The term quadratic in y in eq. (A11) determines the cor-relation function of the stochastic force, h ξ t ξ t ′ i = 2 Re h δ ˜ F ( t ) δ ˜ F ( t ′ ) i . (A14) Appendix B: Green’s functions of the single dot
The single dot lesser Green’s function without couplingin energy space is G < ( ω ) = i Γ L f L ( ω ) + Γ R f R ( ω )( ω − ε x ) + Γ / . (B1)Here we have used the adiabatic approximation for thecenter of mass coordinate, thus the coupling to the os-cillator simply shifts the level ε d by λ h x i , so we have toreplace ε d by the shifted energy ε x = ε d − λ h x i . Thegreater Green’s function is correspondingly G > ( ω ) = − i Γ L [1 − f L ( ω )] + Γ R [1 − f R ( ω )]( ω − ˜ ε ) + Γ / G < ( ω ) − i Γ( ω − ˜ ε ) + Γ / . (B2)Hereby f α ( ω ) = f ( ω − µ α ) = 1 / [exp { β ( ω − µ α ) } + 1] de-note Fermi functions with lead index α , inverse temper-ature β and chemical potential µ α . The time-dependentGreen’s functions are obtained by Fourier transformationas G ≶ ( t ) = R dω π e − iωt G ≶ ( ω ). In the zero–temperaturelimit we have to replace f α ( ω ) = Θ( µ α − ω ). Then for t = 0 an integration leads to G ≶ ( t ) = e − i ˜ εt π X α Γ α Γ (cid:20) − e +Γ / | t | E (cid:8) [ i Ω α sgn( t ) + Γ / | t | (cid:9) sgn( t )+ e − Γ / | t | E (cid:8) [ i Ω α sgn( t ) − Γ / | t | (cid:9) sgn( t ) ± πie − Γ / | t | Θ( ± Ω α ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) Ω α = µ α − ε x , (B3)where the first exponential integral is defined byE ( x ) = Z ∞ dt e − xt t . (B4)For the case t = 0 the lesser Greens function reads − iG < ( t = 0) = 12 π Z ∞−∞ dω X α Γ α Θ( µ α − ω )( ω − ε ) + Γ /
4= 12 − π X α Γ α Γ arctan (cid:2)
2Γ (˜ ε − µ α ) (cid:3) . (B5) For finite temperature one has to regard Fermi func-tions instead of the Heaviside theta function. By virtueof the residue theorem one finds − iG < ( t = 0)= 12 − π X α Γ α Γ Im Ψ (cid:18)
12 + β Γ4 π + i β (˜ ε − µ α )2 π (cid:19) , (B6)with the digamma function Ψ. Appendix C: Characteristics and calculations for thedouble dot1. Green’s function of the DQD
The Green’s functions in the frequency domain accord-ing to the Hamiltonian in equation (19) is derived via theequation of motion method. There the Green’s function G is defined as resolvent of the Hamiltonian H via( ω − H ) G ( ω ) = | φ λ i = | λ i the previousequation yields h λ | ( ω − H ) G ( ω ) | λ ′ i = δ λ,λ ′ (C2)Taking the matrix elements and inserting the Hamilto-nian leads to a set of equations, from whom the Green’sfunctions are derived.Here, the electron-phonon coupling is described adiabat-ically, so the interaction part H SB came in by shiftingthe dot level energies ν L , R → e ν L , R with e ν L , R = ν L , R ∓ λx .Finally the dot Green’s function is given through G D ( ω ) = (cid:18) G LL ( ω ) G LR ( ω ) G RL ( ω ) G RR ( ω ) (cid:19) with the elements G LL ( ω ) = ω − e ν R − Σ R ( ω )[ ω − e ν L − Σ L ( ω )] [ ω − e ν R − Σ R ( ω )] − | T c | G RR ( ω ) = ω − e ν L − Σ L ( ω )[ ω − e ν L − Σ L ( ω )] [ ω − e ν R − Σ R ( ω )] − | T c | G LR ( ω ) = T c ω − e ν L − Σ L ( ω ) G RR ( ω ) G RL ( ω ) = T ∗ c ω − e ν R − Σ R ( ω ) G LL ( ω ) (C3)Hereby we introduced the self energy Σ α , ( α ∈ R , L) cor-responding to the left or the right dot withΣ α ( ω ) = X k | V kα | g kα,kα ( ω ) . (C4) g kα,kα is the undisturbed Green’s function for the leads.We derive the associated retarded and the advancedGreen’s function by using the continuation rules G ( ω ± i + ) → G R,A ( ω )1The derivation of the lesser/greater Green’s function istaking usage of the Keldysh equation: G ≶ αβ = X γ G Rα,γ ( ω ) Σ ≶ γ ( ω ) G Aγ,β ( ω ) (C5)whereas we assume both dots initially unoccupied. Thelesser/greater self energy follows fromΣ ≶ γ ( ω ) = X k | V kγ | g ≶ kγ,kγ ( ω ) , ⇒ Σ <γ ( ω ) = i Γ γ ( ω ) f γ ( ω ) , Σ >γ ( ω ) = − i Γ γ ( ω ) [1 − f γ ( ω )] . (C6)
2. Calculation of the friction for the DQD
The friction is determined by the imaginary part of thecorrelation function of the double dot. From equation(27), we start with expressing the correlation functionthrough Green’s functions and use their Fourier trans-forms: A [ x ]( t )= 2 λ X α,β [2 δ α,β − Z dt ′ t ′ i ∗ h G <βα ( t ′ − t ) G >αβ ( t − t ′ ) − G <βα ( t − t ′ ) G >αβ ( t ′ − t ) i = λ π Im X α,β [2 δ α,β − Z dω Z dω G <βα ( ω ) ∗ G >αβ ( ω ) e i ( ω − ω ) t ( − i ) ∂∂ω Z dt ′ π e − i ( ω − ω ) t ′ = λ π X α,β [2 δ α,β − Z dωG <βα ( ω ) ∂∂ω G >αβ ( ω ) . (C7) In the last step a term was identified as the derivative ofDirac’s delta, this result agrees with the solution in thework by Mozyrsky et al. .The lesser/greater Green’s function Keldysh equation fol-low from the Keldysh equation (C5) and we consider thetunneling rates to be frequency independent, thereforethe self energies are ( T = 0):Σ <γ ( ω ) = i Γ γ Θ( µ γ − ω ) , Σ >γ ( ω ) = − i Γ γ Θ( ω − µ γ ) . (C8)Because of the lengthy calculation, we just outline thecalculation for the LL-term T LL . The other terms( T RL , T LR and T RR ) can be derived in a similar way. Here,the product of retarded and advanced Green’s functionsequates the squared modulus of G α,β , so we can write G R LL ( ω ) G A LL ( ω ) = | G R LL ( ω ) | = | G A LL ( ω ) | ,G R LR ( ω ) G A RL ( ω ) = | G R LR ( ω ) | = | G A RL ( ω ) | . (C9)For simplicity in the following we omit the superscript ofthe retarded Greens function and abbreviate the moduliby | G LL ( ω ) | ≡ | G R LL ( ω ) | , | G LR ( ω ) | ≡ | G R LR ( ω ) | . We obtain for the first term of the friction T LL = Z dω (cid:20) G < LL ( ω ) ∂∂ω G > LL ( ω ) (cid:21) = Γ Z dω (cid:20) | G LL ( ω ) | Θ( µ L − ω ) + | G LR ( ω ) | Θ( µ R − ω ) (cid:21) ∗ " (cid:20) ∂∂ω | G LL ( ω ) | (cid:21) Θ( ω − µ L ) + | G LL ( ω ) | δ ( ω − µ L ) + (cid:20) ∂∂ω | G LR ( ω ) | (cid:21) Θ( ω − µ R ) | G LR ( ω ) | δ ( ω − µ R ) . (C10)Some parts of the integration terms can directly be eval-uated with the help of the delta and Heavyside functions.By choosing the condition µ L > µ R , we get rid of a casedistinction, which would be necessary for two productterms which included a delta-function. After performing also an integration by parts, we arrive at T LL = Γ " | G LL ( µ L ) | + | G LR ( µ R ) | − | G LL ( µ L ) | | G LR ( µ L ) | + 4 | G LL ( µ R ) | | G LR ( µ R ) | +4 µ L Z µ R dω | G LL ( ω ) | (cid:20) ∂∂ω | G LR ( ω ) | (cid:21) . (C11)2In the following the Green’s functions, derived in sec-tion C 1, were inserted, whereas we assume equal tunnel-ing rates for the left and the right side, Γ L = Γ R = Γ.Then a number of integrations by parts is performed todispose of the derivation in the integral term. So we ob-tain a closed expression, whereas N ( ω ) abbreviates thedenominator of the Green’s function T LL = Γ " h ( µ L − e ν R ) + Γ i N ( µ L ) + | T c | N ( µ R ) +2 | T c | h ( µ R − e ν R ) + Γ i N ( µ R ) − | T c | µ L Z µ R dω ( ω − e ν R ) N ( ω ) . (C12) In an analogue way the other terms in (C7) can be de-rived, and finally the solution for the friction is A [ x ]( t ) = λ π Γ " h ( µ L − e ν R ) + Γ i N ( µ L ) + h ( µ R − e ν L ) + Γ i N ( µ R ) + | T c | N ( µ L ) + | T c | N ( µ R ) + 2 | T c | h ( µ R − e ν R ) + Γ i N ( µ R ) − | T c | h ( µ L − e ν R ) + Γ i N ( µ L ) − | T c | h ( µ R − e ν R )( µ R − e ν L ) − Γ i N ( µ R ) − | T c | ( e ν L − e ν R ) µ L Z µ R dω N ( ω ) . (C13)With the Assumption e ν R = − e ν L we get: N ( ω ) = ω + 2 Aω + B , with A = − ( e ν + | T c | − Γ / ,B = e ν + | T c | + Γ / , (C14)and the integral in equation (C13) is given trough I ( A, B ) in section D.
3. Fixed point analysis for the double dot
The fixed points of a nonlinear two dimensional systemcan be investigated with standard methods for linear dy-namical system .The general solution for a two dimensional linear system˙ x = A x is x ( t ) = c e λ t v + c e λ t v (C15)and so determined by the eigenvalues λ , of the matrixA. The constants c , depend on the initial conditionsand v , are eigenvectors.The eigenvalues can be obtained from λ , = 12 ( τ ± p τ − , (C16) thereby τ corresponds to the trace and ∆ to the deter-minant of A. These two qualities determine the evolutionof the trajectories in the phase plane.For a fixed point x ∗ the condition ˙ x = 0 must be ful-filled. Various classes of fixed points exist, whereas thedeterminant ∆ decides which kind of point appears. Incase of saddle points the determinant is negative and itis positive for spirals or nodes. The difference between aspiral and a node is that for the second one the eigenval-ues have no imaginary part.The trace τ defines the stability of nodes and spirals,this is caused by the fact that τ determines the sign ofthe eigenvalue’s real part, for instance with negative realpart decaying oscillations occur and the fixed point is sta-ble, see equation (C15). There also exist some borderlinecases, whereas the centers are the most significant ones,they occur when the trace is equal to zero.This analysis can be assigned to a two dimensional non-linear system ˙ x = f ( x ). By assuming a small disturbance u = x − x ∗ from a fixed point, we can invest if this distur-bance grows or decays by performing a Taylor expansion˙ u = f ( x ∗ , x ∗ ) + u ∂f ∂x (cid:12)(cid:12) x ∗ + u ∂f ∂x (cid:12)(cid:12) x ∗ + h . t . ˙ u = f ( x ∗ , x ∗ ) + u ∂f ∂x (cid:12)(cid:12) x ∗ + u ∂f ∂x (cid:12)(cid:12) x ∗ + h . t . (C17)3The first term is zero and higher terms (h . t . ) can be ne-glected because the disturbance is small. So we get alinearised system ˙ u = J ∗ u , containing the Jacobi matrix J ∗ evaluated at the fixed point coordinates. The above explained analysis can be performed for this system. Thisis valid as long no borderline cases occur, then the higherterms may be more important. Appendix D: Auxiliary integrals (DQD) I ( µ, A, B ) = Z µ −∞ dω ω + 2 Aω + B = 12 √ A − B (cid:20) + 1 p A − √ A − B arctan (cid:0) µ p A − √ A − B (cid:1) − p A + √ A − B arctan (cid:0) µ p A + √ A − B (cid:1) + π (cid:18)s A − √ A − B − s A + √ A − B (cid:19)(cid:21) , I ( µ, A, B ) = Z µ −∞ dω ωω + 2 Aω + B = 12 √ A − B (cid:20) + 12 ln (cid:0) µ + A − p A − B (cid:1) −
12 ln (cid:0) µ + A + p A − B (cid:1)(cid:21) , I ( µ, A, B ) = Z µ −∞ dω ω ω + 2 Aω + B = 12 √ A − B (cid:20) − q A − p A − B arctan (cid:0) µ p A − √ A − B (cid:1) + q A + p A − B arctan (cid:0) µ p A + √ A − B (cid:1) − π (cid:18) q A −√ A − B − q A + √ A − B (cid:19)(cid:21) , I ( A, B ) = Z dω ω + 2 Aω + B ] = 18 B ( B − A ) (cid:20) ω ( B − A − Aω ) B + 2 Aω + ω − ( A − B + A p A − B ) arctan ω √ A −√ A − B √ A − B p A − √ A − B + ( A − B − A p A − B ) arctan ω √ A + √ A − B √ A − B p A + √ A − B (cid:21) . (D1)4 ∗ [email protected] A. D. O’Connell, M. Hofheinz, M. Ansmann, R. C. Bial-czak, M. Lenander, E. Lucero, M. Neeley, D. Sank, H.Wang, M. Weides, J. Wenner, J. M. Martinis and A. N.Cleland, Nature , 697 (2010). S. Groblacher, J. B. Hertzberg, M. R. Vanner, G. D. Cole,S. Gigan, K. C. Schwab, and M. Aspelmeyer, Nat. Phys. , 485 (2009). F. Marquardt and S. M. Girvin, Physics , 40 (2009). A. Naik, O. Buu, M. D. LaHaye, A. D. Armour, A. A.Clerk, M. P. Blencowe, and K. C. Schwab, Nature ,193 (2006). A. A. Clerk, F. Marquardt, and K. Jacobs, New J. Phys. , 095010 (2008). A. D. Armour and M. P. Blencowe, New J. Phys. ,095004 (2008). M. D. LaHaye, J. Suh, P. M. Echternach, K. C. Schwab,and M. L. Roukes, Nature , 960 (2009). M. Blencowe, Phys. Rep. , 159 (2004). D. A. Rodrigues, J. Imbers, and A. D. Armour, Phys. Rev.Lett. , 067204 (2007). D. A. Rodrigues and A. D. Armour, New. J. Phys. , 251(2005). A. D. Armour, M. P. Blencowe, and Y. Zhang, Phys. Rev.B , 125313 (2004). J. Koch and F. von Oppen, Phys. Rev. Lett. , 206804(2005); J. Koch, M. E. Raikh, and F. von Oppen,Phys. Rev. Lett. , 056801 (2005); J. Koch, F. von Op-pen, and A. V. Andreev, Phys. Rev. B , 205438 (2006). D. Fedorets, L. Y. Gorelik, R. I. Shekhter, and M. Jonson,Phys. Rev. Lett. , 166801 (2004). A. D. Armour and A. MacKinnon, Phys. Rev. B , 035333(2002). T. Novotn´y, A. Donarini, and A.-P. Jauho,Phys. Rev. Lett. , 256801 (2003). T. Novotn´y, A. Donarini, C. Flindt, and A.-P. Jauho,Phys. Rev. Lett. , 248302 (2004). T. Brandes, Phys. Rep. , 315 (2005). A. Braggio, J. K¨onig, and R. Fazio, Phys. Rev. Lett. , 026805 (2006); C. Flindt, T. Novotn´y, A. Braggio,M. Sassetti, and A.-P. Jauho, Phys. Rev. Lett. , 150601(2008). A. Mitra, I. Aleiner, and A. J. Millis, Phys. Rev. B ,245302 (2004). L. Glazman and R. Shehkter, Sov. Phys. JETP , 163(1988). N. S. Wingreen, K. W. Jacobsen, and J. W. Wilkins,Phys. Rev. B , 11834 (1989). K. Flensberg, Phys. Rev. B , 205323 (2003). A. A. Clerk, Phys. Rev. B , 245306 (2004). D. Mozyrsky, M. B. Hastings, and I. Martin, Phys. Rev. B , 035104 (2006). R. P. Feynman and F. L. Vernon, Annals of Physics ,118 (1963). D. Mozyrsky, I. Martin, and M. B. Hastings, Phys. Rev.Lett. , 018303 (2004). S. D. Bennett and A. A. Clerk, Phys. Rev. B , 201301(2006). L. S. Schulman,
Techniques and Applications of Path In-tegration (Dover Publications, Inc., 2005). U. Weiss,
Quantum Dissipative Systems , vol. 13 (WorldScientific Publishing, 2008). M. A. Armen and H. Mabuchi, Phys. Rev. A , 063801(2006). D. A. Rodrigues, J. Imbers, T. Harvey, and A. D. Armour,New. J. Phys. , 84 (2007). S. Strogatz,
Nonlinear Dynamics and Chaos: With Ap-plications to Physics, Biology, Chemistry and Engineering (Westview Press, 2000). N. G. Van Kampen,
Stochastic processes in physics andchemistry (Elsevier, 2008), 3rd ed. H. J. Haug and A.-P. Jauho,