PAPER • OPEN ACCESS Recent citations Stochastic thermodynamics of self-oscillations: the - Leaking elastic capacitor as model for electron shuttle active matter Robert Alicki et al - Autonomous Implementation of To cite this article: Christopher W Wächtler et al 2019 New J. Phys. 21 073009 Thermodynamic Cycles at the Nanoscale Philipp Strasberg et al - Stochastic thermodynamic cycles of a mesoscopic thermoelectric engine R. David Mayrhofer et al View the article online for updates and enhancements. This content was downloaded from IP address 84.78.139.232 on 09/07/2021 at 16:06 New J. Phys. 21 (2019) 073009 https://doi.org/10.1088/1367-2630/ab2727 PAPER Stochastic thermodynamics of self-oscillations: the electron shuttle OPEN ACCESS Christopher W Wächtler1 , Philipp Strasberg2,3, Sabine H L Klapp1, Gernot Schaller1 and RECEIVED 18 March 2019 Christopher Jarzynski4 1 Institute of Theoretical Physics, Secr. EW 7-1, Technical University Berlin, Hardenbergstr. 36, D-10623, Berlin, Germany REVISED 2 10 May 2019 Complex Systems and Statistical Mechanics, Physics and Materials Science, University of Luxembourg, L-1511, Luxembourg, Luxembourg ACCEPTED FOR PUBLICATION 3 5 June 2019 Física Teòrica: Informació i Fenòmens Quàntics, Departament de Física, Universitat Autònoma de Barcelona, E-08193, Bellaterra (Barcelona), Spain PUBLISHED 4 Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742, United States of America 1 July 2019 E-mail: christopher.w.waechtler@campus.tu-berlin.de Original content from this Keywords: stochastic thermodynamics, electron shuttle, self-oscillation, nanoscale, nonlinear systems, electron transport, autonomous work may be used under the terms of the Creative engines Commons Attribution 3.0 licence. Any further distribution of this work must maintain Abstract attribution to the Self-oscillation is a phenomenon studied across many scientific disciplines, including the engineering of author(s) and the title of the work, journal citation efficient heat engines and electric generators. We investigate the single electron shuttle, a model nano- and DOI. scale system that exhibits a spontaneous transition towards self-oscillation, from a thermodynamic perspective. We analyse the model at three different levels of description: The fully stochastic level based on Fokker–Planck and Langevin equations, the mean-field (MF) level, and a perturbative solution to the Fokker–Planck equation that works particularly well for small oscillation amplitudes. We provide consistent derivations of the laws of thermodynamics for this model system at each of these levels. At the MF level, an abrupt transition to self-oscillation arises from a Hopf bifurcation of the deterministic equations of motion. At the stochastic level, this transition is smeared out by noise, but vestiges of the bifurcation remain visible in the stationary probability density. At all levels of description, the transition towards self-oscillation is reflected in thermodynamic quantities such as heat flow, work and entropy production rate. Our analysis provides a comprehensive picture of a nano-scale self-oscillating system, with stochastic and deterministic models linked by a unifying thermodynamic perspective. 1. Introduction Self-oscillation has been described as ‘the generation and maintenance of a periodic motion by a source of power that lacks a corresponding periodicity’ [1]. As opposed to resonant systems, in which the driving source is modulated externally, the energy required to sustain self-oscillations is supplied by a constant source. The phenomenon is familiar from everyday life, e.g. the human voice and the sound of a violin string. Autonomous oscillations appear in a wide range of biological systems and chemical and biochemical processes [2–4] controlling, e.g. the beating of the heart, circadian cycles in body temperature or the Belousov–Zhabotinsky reaction. By converting direct current into stable oscillations, self-oscillatory systems provide a useful transduction mechanism for the design of autonomous motors and heat engines. One particularly interesting system exhibiting self-oscillation is the electron shuttle, first proposed by Gorelik et al [5], where the mechanical oscillation of a metallic grain is achieved by sequential electron tunnelling between the grain and two connecting leads. This coupled system of mechanical and electronic degrees of freedom has drawn considerable theoretical and experimental attention since its original proposal. Theoretical descriptions of the system range from full quantum mechanical models of the coherent dynamics [6–12] to semiclassical [12–14] and completely classical descriptions [5, 15–17]. The electron shuttle has been experimentally realized by a vibrational fullerene molecule [18], gold grains [19–21] as well as nanopillars [22] as molecular junctions between two leads. Also macroscopic electron shuttles, consisting of a pendulum between two capacitor plates, have been investigated [23]. Reviews on the electron shuttle can be found in [24–29]. © 2019 The Author(s). Published by IOP Publishing Ltd on behalf of the Institute of Physics and Deutsche Physikalische Gesellschaft New J. Phys. 21 (2019) 073009 C W Wächtler et al Classical self-oscillating systems have been analysed using the tools of deterministic nonlinear dynamics [1, 30]. For the electron shuttle in particular, a sharp transition from stationarity to self-oscillation arises due to a Hopf bifurcation as the voltage difference between the leads crosses a threshold value [12]. While the dynamical description is well understood, the electron shuttle has not yet been thoroughly investigated from a thermodynamic perspective. Our aim in this paper is to provide such a perspective and lay out the groundwork for further thermodynamic analysis of the electron shuttle as a paradigmatic isothermal engine that converts direct electric current into periodic mechanical motion [31]. Because the electron shuttle is a nanoscale device, fluctuations play a central role in our analysis, in contrast with deterministic classical models. In our analysis we will apply the tools of stochastic thermodynamics, a framework that formulates the laws of thermodynamics at the single-trajectory level and is particularly useful for investigating the thermodynamic behaviour of nanoscale systems. As described in review articles and monographs [32–38], stochastic thermodynamics has been applied to a wide range of topics, including far-from-equilibrium fluctuation theorems, the operation of biomolecular machines, feedback control of nanoscale systems, the thermodynamic arrow of time, and the thermodynamic implications of information processing. In recent years a number of models of non-autonomous stochastic heat engines have been proposed and investigated within the stochastic thermodynamic framework [39–44]. In these models, externally applied time- periodic driving leads to the conversion of thermal fluctuations into work. By contrast, autonomous nano-scale engines are characterized by the absence of an externally imposed cycle. A variety of such autonomous engines have recently drawn both theoretical and experimental attention. First, thermoelectric devices, which use the interplay of thermal and chemical gradients to perform useful tasks, were proposed [45–50] and experimentally realized using quantum dot (QD) structures [51–53]. Second, stochastic self-oscillatory engines were analysed, including a Brownian gyrator [54, 55], a rotor engine [56, 57], a heat engine based on Josephson junctions [58], solar cells [59–61] or mechanical resonators [55, 62], and an experimental realization of the Feynman’s ratchet- and-pawl mechanism [63]. While particular thermodynamic aspects such as nonequilibrium hot electron transport [64], subresonance inelastic electronic transport [65, 66] and tip-induced cooling [67] have been investigated, a systematic thermodynamic description of a nano-scale self-oscillating system, such as the one we provide for the electron shuttle, is still missing. We investigate the electron shuttle at three different levels of description: the fully stochastic level modelled by a Fokker–Planck equation (FPE) and the equivalent Langevin equation, a mean-field (MF) model described as a deterministic dynamical system, and an intermediate perturbative model based on multiple scale (MS) perturbation theory, containing both deterministic and stochastic elements. We study the dynamics and obtain statements of the first and second laws of thermodynamics at all three levels of description. In doing so, we draw a direct line between our stochastic thermodynamic model of the electron shuttle and the nonlinear dynamic model of [5, 15]. We find that the abrupt onset of self-oscillatory behaviour observed at the deterministic level, appears at the stochastic level as a smoothed but nevertheless discernible transition from stationarity to self- oscillation. At all three levels of description, this transition is reflected in thermodynamic quantities such as the rates of heat flow and entropy production. Outline: The article starts with a short review of the basic idea of an electron shuttle (section 2) followed by mathematical descriptions of the system at the different levels mentioned above (section 3): the fully stochastic model in section 3.1, the MF approach in section 3.2, and the intermediate, perturbative model in section 3.3. The dynamics at the different levels are discussed and compared in section 3.4. In section 4 the first and second laws of thermodynamics are derived at the different levels of description (sections 4.1–4.3), followed by a discussion of the thermodynamic behaviour of the electron shuttle (section 4.4). Finally, in section 5, we discuss our findings and point out future applications. 2. Phenomenology In this section we explain the basic mechanism of the electron shuttle (see also figure 1) before introducing the mathematical descriptions in section 3. The shuttle is composed of a metallic grain [19–21] or molecular cluster [22, 68] and a nanomechanical oscillator (e.g. a cantilever [69] or an oscillating molecule [18, 26]), which hosts the grain or cluster and can oscillate. Furthermore, the shuttle is tunnel-coupled to two leads, such that electrons can jump between the leads and the grain. Here, the rate of tunnelling depends on the position of the shuttle— the closer the shuttle is to the lead, the larger is the rate of tunnelling. A bias voltage applied to the two leads then generates an electric field. The shuttle mechanism works as follows: when the shuttle is close to the reservoir with higher chemical potential, electrons are loaded onto the grain. The electrostatic force due to the electric field between the leads pushes the negatively charged shuttle towards the reservoir with lower chemical potential similar to a charged particle in a capacitor (see figure 1 left). As the shuttle approaches the positively biased reservoir with lower chemical potential, the electrons are unloaded from the grain, leaving it uncharged. Due to 2 New J. Phys. 21 (2019) 073009 C W Wächtler et al Figure 1. Illustration of the model and shuttling mechanism: a single-level dot is coupled to two electronic reservoirs with chemical potentials μL and μR and (inverse) temperatures βL/R. If an electron tunnels into the dot (q = 1) the electrostatic force generated by the bias voltage between the leads pushes the oscillator towards the right lead (left figure). If the dot is unoccupied (q = 0) only the oscillator restoring force acts on the system, pushing the oscillator back towards the centre. For weak friction the shuttle may pass the centre and approach the left lead, closing the cycle (right figure). The tunnelling rates into and out of the left and right reservoirs depend exponentially on the position x of the shuttle, such that electrons tunnel more likely between the QD and the closer lead (see main text), as indicated by the thickness of the arrows. Additionally, the shuttle is subject to thermal noise (not shown). the oscillator restoring force the shuttle returns (see figure 1 right) and the cycle starts again. Above a critical value of the applied bias voltage the damping due to friction is overturned by the electrostatic force. As a result, oscillations of the shuttle are sustained and in each cycle a number of electrons are transported from one lead to the other. 3. Modelling In this section we discuss different levels of description of the electron shuttle, i.e. a fully stochastic description in section 3.1, a MF approximate description in section 3.2 and a perturbative description based on time scale separation in section 3.3. We will then compare and discuss the dynamics of the system at the different levels in section 3.4. In the literature there exist proposals to describe the electron shuttle fully quantum mechanically [6–10, 12], semiclassically [12–14] or fully classically [5, 15–17]. In this work we describe the system classically, which is justified if the intra-grain electronic relaxation time is much shorter than the tunnelling charge relaxation time [25]. The latter is the case for an experimental realization of the shuttle with a gold grain [19] and is sometimes referred to as classical shuttling of particles [25]. The underlying mechanism (tunnelling of electrons via Fermi’s golden rule), is nevertheless intrinsically quantum. 3.1. Fully stochastic description We here introduce the specific model of an electron shuttle considered in this work. In contrast to the original proposal [5] we idealize the QD by assuming Coulomb blockade. That is, we assume the QD can accept no more than a single excess electron, due to Coulomb repulsion. Hence the QD charge state can take the two values q=0 (empty) and q=1 (occupied). In this scenario electrons can be transferred one by one between the two reservoirs [70–72]. We use this simplified model of a single electron shuttle for illustrational and numerical purposes, but the phenomenology discussed in section 2 does not change if multiple electrons are allowed on the QD. The QD with on-site energy ε is hosted by a nanomechanical oscillator. In the following we will refer to this combined system of QD and oscillator as a ‘shuttle’. We describe the movement of the oscillator in one dimension with position x Î  and velocity v Î  . The charge q of the shuttle (setting the electron charge e≡1) can change due to electron tunnelling with one of the two electronic leads, left or right, with chemical potentials μL=ε+V/2 and μR = ε − V/2 for the left and right reservoir, respectively. The bias voltage between the two fermionic reservoirs is then given by V = mL - mR . The QD charge state and the motion of the shuttle are coupled by the electric field that is generated by the bias voltage and assumed to be homogeneous between the leads [5]. Thus an electrostatic force Fel=αVq acts on the shuttle when it is charged (q=1), pushing it towards the reservoir with lower chemical potential (see figure 1 left). Here α is an effective inverse distance between the leads. When there is no excess electron on the shuttle (q = 0) this electrostatic force is absent. The mechanical vibrations of the QD are modelled as a harmonic oscillator with an effective mass m [5, 18, 25, 26]. From a classical point of view, this restoring force can be explained through interactions between the shuttle, its anchor and the leads, which can be approximated by a harmonic potential [18]. The restoring force acting on the shuttle is then given by Fharm = -kx with spring constant k, and the shuttle is damped by 3 New J. Phys. 21 (2019) 073009 C W Wächtler et al Fdamp=−γv with friction coefficient γ. We assume underdamped motion to enable the possibility of oscillatory shuttling. Additionally, we connect the oscillator to its own heat bath at inverse temperature βosc stemming from a dissipative medium in equilibrium. The state of the shuttle is described by the triple (x, v, q). Combining the electron jumps with the underdamped oscillations and thermal fluctuations, we describe the dynamics of the shuttle by a generalized FPE ¶p ⎡ ¶ ¶ ⎛k g aV ⎞ ¶2 ⎤ ¶t = ⎢-v ⎣ ¶x + ⎜ ¶v ⎝ m x+ v- m q⎟ + D 2 ⎥ p + m ⎠ ¶v ⎦ å Rqqn ¢ (x ) p (x, v, q¢, t ). (1) q¢n Here, p≡p(x, v, q, t) denotes the joint probability density to find the shuttle at position x and velocity v with q Î {0, 1} electrons at time t, and we have introduced a velocity diffusion coefficient D = g (b oscm2). The first term of equation (1) describes the underdamped evolution of the oscillator in the potential Uq (x ) = 1 2 kx 2 - aVqx , at fixed electron charge q. The second term couples the mechanical variables (x, v) to the charge state (q), through a rate equation describing transitions from state q¢ to state q, corresponding to the tunnelling n of electrons between the QD and the fermionic lead n Î {L, R}. The transition rates Rqq ¢ (x ) are given by L R10 (x ) = Ge-x lf L (e - aVx ) , R R10 (x ) = Ge+x lf R (e - aVx ) , L R 01 (x ) = Ge-x l [1 - f L (e - aVx )] , R R 01 (x ) = Ge+x l [1 - f R (e - aVx )] , (2) and n Rqq (x )= -åq ¢¹ qRqn¢ q (x ), which guarantees the conservation of probability. Here, Γ denotes the bare transition rate, which for simplicity we take to be equal for the two fermionic reservoirs. The probability for quantum mechanical tunnelling is exponentially sensitive to the tunnelling distance, such that the tunnelling amplitudes are modulated by the dimensionless displacement x/λ of the centre of mass of the shuttle [5, 13, 73–75], where λ is a characteristic tunnelling length. Furthermore, the rates depend on the probability of an electron (hole) with a matching energy in the reservoir, i.e. on the Fermi distribution f n (w ) º [exp (b n (w - mn )) + 1]-1 with inverse temperature β ν. Note that the quantity e - aVx enters the Fermi functions in equation (2), as the energy of the shuttle depends on both the QD energy ε and the electrostatic potential −αVx (see also section 4.1). Equation (1) describes a system connected to three reservoirs at generally different temperatures: a thermal reservoir of the oscillator at inverse temperature βosc and two fermionic reservoirs with inverse temperature β ν and chemical potential μ ν. While the derivations in this work are general, when solving the dynamics numerically we will focus on the case of equal temperatures, b osc = b n = b . Nonequilibrium conditions then arise solely due to the applied bias voltage, i.e. mL ¹ mR . Since the space of dynamical variables defined by the triple (x, v, q) is large, solving equation (1) numerically is expensive. We therefore turn to the trajectory representation of the single electron shuttle. The coupled stochastic differential equations dx = v dt , (3) mdv = ( - kx - gv + aVq) dt + 2Dm 2 dB (t ) , (4) dq = å dqn = å (q ¢ - q) dNqn¢q (x, t ) (5) n nq ¢ produce the FPE(1) at the ensemble level, as we show in appendix A by explicitly looking at the evolution of averages. In equation (4) the thermal fluctuations are taken into account by a Wiener process dB(t) with zero mean  [dB (t )] = 0 and variance  [(dB (t ))2] = dt . Here, [•] denotes an average of the stochastic process. Equations (3) and (4) represent for a fixed q the Langevin equation of an underdamped particle moving in the shifted harmonic potential Uq(x). Equation (5) describes changes in the charge state q due to the stochastic tunnelling of electrons. The independent Poisson increments dNqn¢ q (x , t ) Î {0, 1} obey the statistics:  [dNqn¢ q (x , t )] = R qn¢ q (x ) dt , dNqn¢ q (x , t ) dNqq n˜ n ˜ (x , t ) = dq ¢ q˜dnn˜ dNq ¢ q (x , t ). (6) The first equation specifies that the average number of jumps into state q¢ from a state q in a time interval dt is given by the tunnelling rate Rqn¢ q (x ). The second line in equation (6) enforces that only one tunnelling event per time interval can occur, i.e. either all dNqn¢ q (x , t ) are zero or dNqn¢ q (x , t ) = 1 for precisely one set of indices q, q¢ and ν. Well-known models emerge as a simple limit of our description. First, for α → 0 the motion of the oscillator becomes independent of the charge state q, and equations (3) and (4) describe a simple underdamped harmonic oscillator. However, the tunnelling of electrons still depends on x (see equation (2)) and therefore the QD remains coupled to the oscillator. Second, a complete decoupling of the QD and the oscillator is achieved in the 4 New J. Phys. 21 (2019) 073009 C W Wächtler et al limit l  ¥ and α→0. In that case the QD coupled to the fermionic leads describes the well known single electron transistor (SET) [76–78]. 3.2. MF approximation In order to understand the nonlinear dynamics of the compound system of QD and oscillator, we first look at the MF equations derived from the full stochastic evolution. From the FPE, equation (1), we obtain for the ensemble averaged position á xñ and velocity á v ñ: d á xñ = á vñ , dt d m á v ñ = - k á x ñ - g á v ñ + aVp1. (7) dt Here and throughout the paper, á •ñ = ò dx dv åq • p (x, v, q, t ) (8) denotes an ensemble average, and p0 = ò dx dv p (x, v, 0, t ) (9) p1 = ò dx dv p (x, v, 1, t ) (10) are the probabilities for the QD to be empty and occupied, respectively. Equation (7) are exact, but in order for them to form a closed set we need an expression for dp1/dt. Integrating equation (1) over x and v we obtain ¶ ò dx dv å n n p = R qq ¢ (x ) p (x , v , q ¢, t ) = å á R qq ¢ (x )ñ. (11) ¶t q q¢n n Due to the nonlinearity of the tunnelling rates with respect to x (see equation (2)) we approximate n n á R qq ¢ (x )ñ » R qq ¢ (á x ñ) , (12) n and we refer to this as the MF approximation. Note that if the tunnelling rates Rqq ¢ (x ) were linear in x, equation (12) would be an equality and equation (7) would be closed without the MF approximation. The MF approximation is thus described by the nonlinear differential equations d x¯ = v¯, (13) dt d m v¯ = - kx¯ - gv¯ + aVq¯ , (14) dt d dt p¯ = å Rn (x¯) p¯ , (15) n where p¯ º ( p¯0 , p¯1 ) and q¯ º p¯1. The entries of the rate matrix Rn (x¯) are definded by equation (2), i.e. n [Rn (x¯)]qq ¢ º Rqq ¢ (x¯). The overbars denote that the quantities are governed by MF equations. Within the MF description, the QD is still described by a probability and therefore behaves stochastically, whereas the oscillator is fully deterministic. The temperature of the oscillator bath, βosc, does not appear in equations (13)–(15). In effect the MF approximation describes a macroscopic system for which thermal fluctuations are negligible, as would be expected in the limit of large oscillator mass. In this limit the oscillation period 2p m k becomes much longer than the time scale associated with changes in the charge state q, hence the charge state can be replaced by its local-in-time average, as reflected in equation (14). Equations (13)–(15) are identical to those found in the original proposal of Gorelik et al [5] for the case of one excess electron. 3.3. MS perturbation theory To improve on the MF approximation, which only captures the average dynamics of the electron shuttle, we can perturbatively solve the FPE, equation (1), by assuming a separation of time scales between the short dwell-time of electrons and the slow movement of the oscillator, and applying MS perturbation theory [79–81]. Specifically, we assume that during one oscillation there are many electron tunnelling events: G  k m . We provide details of the MS calculation in appendix B, and summarize the result here. 5 New J. Phys. 21 (2019) 073009 C W Wächtler et al Working to first order in the perturbation we obtain p (x , v , q , t ) » pq (x ) p˜ (x , v , t ) , (16) where the vector (p0 (x ), p1 (x ))T denotes the stationary state of the rate matrix R (x ) = ån Rn (x ), i.e. it is the right eigenvector corresponding to the zero eigenvalue, normalized to unity: åq pq (x ) = 1. The probability density to find the oscillator at position x with velocity v at time t is given by p˜ (x , v , t ), which obeys the FPE: ¶p˜ ⎡ ¶ ¶ ⎛k g˜ (x ) aV ⎞⎤ ¶ 2p˜ = ⎢-v + ⎜ x+ v- qeq (x ) ⎟ ⎥ p˜ + D˜ (x ) 2 . (17) ¶t ⎣ ¶x ¶v ⎝ m m m ⎠⎦ ¶v Here, qeq (x ) º p1 (x ) is the instantaneous stationary charge of the QD given by f L (e - aVx ) - f R (e - aVx ) qeq (x ) = + f R (e - aVx ). (18) 1 + e 2x l As seen in equation (17), the effects of the electronic degrees of freedom on the evolution of the oscillator are incorporated into an effective potential Ueff(x), along with position dependent friction and diffusion coefficients: Ueff (x ) = kx - aVqeq (x ) , (19) aV ¶qeq (x ) g˜ (x ) = g - , (20) c (x ) ¶x a 2V 2qeq (x ) D˜ (x ) = D - [1 - qeq (x )] , (21) m 2c (x ) where c (x ) = -2G cosh (x l ) is the non-zero eigenvalue of R(x). The zeroth order perturbation (see appendix B) corresponds to an adiabatic approximation, i.e. infinite time scale separation, G  ¥. In that limit we have g˜ (x )  g and D˜ (x )  D , and equation (17) describes underdamped Brownian motion in an effective potential Ueff(x), at inverse temperature g (Dm2) = b osc . In this situation, detailed balance is satisfied and the oscillator relaxes to an effective equilibrium state, with no self- sustained oscillations. By contrast, in the first order perturbation represented by equation (17), the x- dependence of g̃ (Dm ˜ 2) breaks detailed balance, giving rise to non-equilibrium behaviour and allowing for the possibility of self-oscillations. To solve equation (17) approximately, we parametrize x and v by the energy  and the oscillation phase θ, 2 2 x= sin q , v = cos q , (22) k m such that 12 kx 2 + 12 mv 2 =  , and we assume that the probability density does not depend on the phase [82]: pˆ ( , q, t ) » pˆ ( , t ). With the transformations of equation (22) and the latter assumption we find a FPE for the energy distribution by averaging over the angle θ (see appendix C): ¶ ¶ ⎡ ⎛ gˆ ¶ ⎞ ⎤ pˆ ( , t ) = ⎢2 ⎜ + mDˆ ⎟pˆ ( , t ) ⎥ , (23) ¶t ¶ ⎣ ⎝ m ¶ ⎠ ⎦ where pˆ (E , t ) is the transformed probability distribution p˜ (x , v , t ). The effective friction and diffusion parameters take after the transformation the form 1 2p gˆ ( ) = 2p ò0 dq g˜ (x ) cos2 q , 1 2p Dˆ ( ) = 2p ò0 dq D˜ (x ) cos2 q. (24) Solving for the steady state of equation (23), i.e. ¶pˆ ¶t = 0, we get [83] ⎛  gˆ (¢) ⎞ pˆss ( ) =  exp ⎜ - ⎝ ò0 d¢⎟ , m Dˆ (¢) ⎠ 2 (25) where  is a normalization constant. 3.4. Dynamics on the different levels of description In this section we discuss and compare the dynamical behaviour of the single electron shuttle on the different levels of description introduced before. All numerical results in this section are obtained using the parameter values specified in appendix D. 6 New J. Phys. 21 (2019) 073009 C W Wächtler et al Figure 2. Top: exemplary trajectories of the position (orange) together with electron jumps between the QD and the left (red) and right (blue) reservoir for (a) bV = 1.0 and (b) β V=40.0 showing clearly the shuttling for a large bias voltage. Bottom: probability density of the oscillator in phase-space for (c) βV=1.0 and (d) βV=40.0 simulated from equations (3)–(5) (see also appendix D). The circular orbit indicates self-oscillations. The white dot and circle correspond to MF solutions (see equations (13)–(15)). 3.4.1. Stochastic dynamics We start by looking at the fully stochastic model, given by equation (1). Rather than solving the FPE directly, we generated stochastic trajectories evolving under the Langevin equations (3)–(5). Figure 2 shows trajectory segments x(t) for two different values of the applied bias voltage: βV=1.0 (figure 2(a)) and βV=40.0 (figure 2(b)). The two figures show quite different behaviour of the stochastic position (orange) as well as the tunnelling of electrons schematically indicated by red (left lead) and blue (right lead) bars. Here, a negative value n of IM = dqn dt denotes the jump of an electron from the QD into the reservoir ν (dqn = -1) whereas a positive n value indicates the reverse process (dqn = 1). Note that the stochastic current IM along a trajectory shows up as delta-peaks. In figures 2(a) and (b) we plot dqn for clearness of the figures. For a small bias voltage (panel (a)) tunnelling events are frequent and the position x(t) oscillates irregularly around the equilibrium position. In contrast, figure 2(b) shows regular oscillations of the position and fewer tunnelling events. Also, the tunnelling events in figure 2(b) are synchronized with the shuttling: during each period of oscillation, the empty shuttle picks up one electron from the left reservoir (red bar) as it moves past the origin in a leftward direction, dx/dt<0, and it releases that electron to the right reservoir on its way back (blue bar), as it moves past the origin in a rightward direction, dx/dt > 0. (Typically, immediately after releasing the electron the shuttle picks up another electron from the left reservoir and quickly delivers it to the right reservoir.) This behaviour reflects the mechanism of single electron shuttling discussed in section 2. The directed shuttling of electrons coincides with self-oscillation, as illustrated in figures 2(c) and (d), which show the stationary probability density of the oscillator, p (x , v ) = åq p (x , v , q ), obtained by simulating a long trajectory evolving under equations (3)–(5) and assuming ergodicity (see appendix D). For small bias voltage (panel (c)), the probability density is peaked close to the origin, showing no sign of regular oscillations. When the applied voltage is 7 New J. Phys. 21 (2019) 073009 C W Wächtler et al Figure 3. MF position x̄ (orange solid) as well as right (blue dotted) and left (red dashed) matter current, I¯M R L and I¯M , during one period in the asymptotic limit. From (a)–(c) the bias voltage V is increased. Below bV¯cr = 15.0 the system is equivalent to a SET and x¯ = aVq¯ k is constant as indicated also by the illustration in (a) (regime (I)). Above the critical voltage the system oscillates and after a crossover regime (II) (panel (b)) the system acts as an electron shuttle (III) transporting one electron per cycle. The illustrations in (c) indicate the position of the oscillator during the cycle. larger (panel d), the probability density is concentrated around a circular orbit, revealing self-oscillatory harmonic motion with some amplitude and phase noise. These behaviours are consistent with the trajectories shown in figures 2(a) and (b), respectively, and they suggest that there exists a value of the applied bias voltage V above which the shuttle oscillates, as we will discuss below. Similar oscillator distributions in phase space have been observed for Wigner functions in semiclassical descriptions of the single electron shuttle [9, 10]. Note that V enters the equations governing the dynamics via the coupling to the electrostatic field and via the chemical potentials (see equation (2)). 3.4.2. MF dynamics We now turn to the MF dynamics. The white dot and circle in figures 2(c) and (d) correspond to the solutions of the MF model given by equations (13)–(15). As we can see the MF solutions coincide very well with the stochastic phase-space distribution. As shown in previous extensive studies [12], when the parameter V crosses a critical value V¯cr , the MF system undergoes a Hopf bifurcation from a stable fixed point to a stable limit cycle. For our choice of parameters the bifurcation takes place at bV¯cr = 15.0 (see appendix E). Three dynamical regimes can be characterized, as we discuss below. SET regime (I): the point (x¯ fix , v¯fix ) = (aVq¯ k , 0) is a fixed point of equations (13)–(14), and below the critical value of the applied voltage this fixed point is stable: from any initial conditions the oscillator spirals into this point, hence at steady state the MF system does not oscillate (see figure 2(c)). In figure 3(a) we find the steady state solution for the MF position at x¯ = x¯ fix = 0.006l and electron currents IML = -I¯MR = 0.122G, describing a fixed oscillator and a constant matter current from left to right lead. The electrostatic force cannot overcome friction, n and the transition rates Rqq ¢ (x¯) (see equation (2)) are constant since x̄ is constant at steady state. The dynamics of the QD are then described by a simple rate equation equivalent to the classical master equation of the SET [37], leading to a net electron current I¯M = I¯ML = -I¯MR = G 2 [ f L (e - aVx¯ fix ) - f R (e - aVx¯ fix )] sech (x¯ fix l ). Note that at the stochastic level (see figure 2(a)) the oscillator is not fixed—only the average position and velocity are equal to the fixed point values. Shuttling regime (III): for a bias voltage V  V¯cr the system is self-oscillating and therefore acts as a shuttle transporting one electron from one lead to the other during each cycle (see figures 2(b) and 3(c)). When the shuttle is occupied by an electron, the electrostatic force is sufficient to overcome friction, leading to 8 New J. Phys. 21 (2019) 073009 C W Wächtler et al Figure 4. Panel (a): numerical solutions of pˆss (E ) (MS, dotted) together with the equivalent plots of the stochastic solution p(E) (FPE, solid) as a function of the energy E of the oscillator. A maximum larger than zero indicates that the system is oscillating. For large values of V the MS solution deviates from the stochastic solution due to the break down of the time scale separation. Panel (b): ‘order parameter’ A of the oscillations as a function of V for the FPE (orange solid), the MF (red dashed) and the MS description (blue dotted): all three descriptions predict an onset of oscillation and agree quite well for the chosen set of parameters. The inset shows a zoom into the onset region. self-sustained oscillations. Note that perfect shuttling, i.e. transport of one electron per oscillation, only occurs at very large bias voltages. For a bias of bV¯ = 55.0 there are still more tunnelling events than from shuttling n electrons one by one, which in figure 3(c) can be seen from the fact that both currents I¯M are finite when x≈0. We also see this in the stochastic case very clearly (see figure 2(b)). Crossover regime (II): when V » V¯cr the system exhibits both SET and shuttle behaviour. Above the critical bias voltage V¯cr the MF fixed point is unstable and a small perturbation to the system causes variations in the charge of the QD. The electrostatic force acting on these charge variations provides positive feedback on the oscillator and compensates for losses due to friction. The asymptotic MF state is characterized by periodic oscillations of the position x̄ and velocity v̄ —as in the shuttling regime—as well as charge q̄ and matter currents n I¯M (see figure 3(b)). However, throughout the entire period of oscillation the QD is able to exchange electrons with both leads, as in the SET regime. As the bias voltage is increased, the amplitude of oscillations increases, and the time during which the shuttle exchanges electrons with the reservoirs decreases, and finally only one electron is transferred per cycle, which corresponds to the pure shuttling regime. 3.4.3. Perturbative dynamics To gain further insight into the transition to self-oscillation, we solve the full FPE, equation (1), perturbatively by imposing a time scale separation between the rapid tunnelling events of electrons and the slow movement of the oscillator (see section 3.3). In figure 4(a) we plot the steady state probability density p˜ss ( ) obtained from this calculation (see equation (25)) for different applied bias voltages (dotted). We also plot the corresponding energy distributions pss ( ) = ò d ( - kx 2 2 - mv 2 2) pss (x , v ) dx dv determined from numerical simulations of the full stochastic evolution (solid). The two sets of distributions show similar behaviour: for small voltages the maximum occurs at  = 0 but for larger values of V the distributions are peaked at non-zero values of the energy, corresponding to self-oscillation as discussed earlier. For larger values of the voltage the maximum of the probability density occurs at smaller values of the energy for the stochastic case, when compared with the MS results. This deviation can be understood in terms of the underlying assumption of time scale separation for the MS perturbation theory: as the bias voltage is increased the system transitions from the SET regime (with clear time scale separation) to the shuttling regime (where time scales are comparable). 3.4.4. Comparison Finally, we compare all three levels of description in terms of an order parameter A that quantifies the magnitude of self-oscillation. In the MF case the position at long times performs oscillations of the form x¯ (t ) = A¯ cos (w¯ t + j¯ 0) + x¯ fix , and we choose AMF = A¯ as our order parameter. In the stochastic case we consider the probability density at v=0, i.e. p(x, v=0), which in the case of large self-oscillations resembles a pair of well-separated peaks (see figure 2(d)). We fit p(x, v = 0) to a normalized sum of Gaussians, g (x ) = N {exp [-(x - c - x 0 )2 2s 2] + exp [-(x - c + x 0 )2 2s 2]} with fit parameters5 σ2, c, and x0. We 5 We have to include a shift c because the the orbit is not exactly centred around the origin. This shift is the counterpart of xfix discussed in the MF context in section 3.4.2. 9 New J. Phys. 21 (2019) 073009 C W Wächtler et al then define the self-oscillation amplitude AFPE in terms of the value(s) x at which the function f (x ) = g (x + c ) has a maximum: when x 0  s , f (x) has a unique maximum at x=AFPE=0, and when x0>σ, f (x) has distinct maxima at x = AFPE . Note that this definition is not sensitive to small oscillations of the stochastic system, as it gives AFPE = 0 when x 0  s , even though the shuttle may be self-oscillating. Finally, for the MS perturbative solution, we define the self-oscillation amplitude as AMS = 2 max k , where  max is the value of  at which the function pˆss ( ) is maximized (see equations (22) and (25)). Similarly to AFPE, and for the same reason, AMS is not sensitive to small oscillations. In figure 4(b) we show the amplitudes of oscillation A for the different levels of description as a function of the applied voltage: FPE (orange solid), MF (red dashed) and MS (blue dotted). All three descriptions show an * = 13.2 onset of oscillation at a critical value of the voltage. The specific values are given by bV¯cr = 15.0, bVFPE * and bV MS = 13.6. These values are surprisingly close to each other, in particular the MS analysis accurately reflects the onset seen in the full FPE simulations. Recall, however, that since AFPE and AMS are not sensitive to small oscillations, the onset to self-oscillation in the FPE and MS cases may not be as abrupt as suggested by the data in figure 4(b). As V is increased, the perturbative solution deviates from the stochastic amplitude of oscillation due to the lack of time scale separation, as discussed earlier (see section 3.4.3). On the other hand, for a large bias voltage the MF and stochastic description coincide quite well, as the deterministic component of the dynamics becomes dominant and the fluctuations become less important. Note that the onset of oscillations in the stochastic case may vary somewhat according to the choice of fitting function g(x). Also, due to fitting of the probability density, AFPE is quite noisy close to the onset; see inset of figure 4(b). Despite these caveats, we see that a transition towards self-oscillation can be identified at all three levels of description. Next, we investigate whether this transition is reflected in thermodynamic quantities such as chemical work rate, heat flow, and entropy production rate. 4. Thermodynamics In this section we formulate the first and second law of thermodynamics at the different levels of description introduced above—stochastic, mean field and multiscale perturbative. In each case we introduce precise definitions of essential thermodynamic quantities, namely heat, work, and entropy production. With these definitions, we compare the thermodynamic behaviour of the electron shuttle at the different levels of description, focusing on the thermodynamic signatures of the onset of spontaneous self-oscillation. 4.1. Stochastic thermodynamics We start by looking at the full stochastic model. The total energy of the coupled system is given by mv 2 kx 2 E= + + eq - aVxq , (26) 2 2 where the first two terms correspond to the kinetic and potential energy of the harmonic oscillator and the third term is the energy of the QD. The last term describes the interaction energy of the oscillator with the QD, which is given by an electrostatic energy analogous to that of a charged particle in a capacitor with a constant electrostatic field of strength αV. By the first law of thermodynamics, a change in the total energy of the system is due either to exchange of heat or to work performed by (on) the system. The change of total energy in the electron shuttle is expressed as [33] dE = kx ◦ dx + mv ◦ dv - aVq ◦ dx + (e - aVx ) ◦ dq , (27) where ◦ denotes Stratonovich-type calculus. In the last term, dq = ån dqn and dqn = åq ¢(q¢ - q ) dNqn¢ q (x , t) denotes an electron jump with respect to reservoir ν. The second term involving the velocity can be re-expressed by multiplying equation (4) with v. This yields mv ◦ dv = ( - kxv - gv 2 + aVqv ) dt + 2g b osc v ◦ dB (t ) = - kx ◦ dx + aVq◦dx - gv 2dt + 2g b osc v ◦ dB (t ). (28) Inserting equations (28) into (27), we get dE = (e - aVx ) ◦ dq - gv 2dt + 2g b osc v ◦ dB (t ) = dQ L + dQ R + dW chem + dQ osc. (29) 10 New J. Phys. 21 (2019) 073009 C W Wächtler et al Here, we have introduced the chemical work dW chem = ån mn dqn and the heat flow to the oscillator from its thermal reservoir due to friction and thermal noise dQ osc = -gv 2dt + 2g b osc v ◦ dB (t ) [33]. The remaining terms in equation (29) are identified as heat exchanged with the reservoir ν, defined as dQ n = (e - aVx - mn ) ◦ dqn . We use the convention that work performed on the system is positive as is heat transferred from a reservoir into the system. With these definitions of heat and work we can derive a consistent second law as we will show later in this section. The average change in total energy is given by averaging equation (29) over many realizations, equivalently by averaging with respect to the probability density (see section 3.1 and appendix A): dE = á Q˙ Lñ + á Q˙ Rñ + á W˙ chemñ + á Q˙ oscñ. (30) dt () d We note that derivatives with respect to time dt denote exact (or complete) differentials whereas a dot (·) denotes inexact ones. The average heat absorbed from reservoir ν is given by á Q˙ n ñ = (e - mn ) á IM n n ñ - aV á xIM ñ. (31) n Here, á IM ñ º  [dqn dt ] is the matter current from reservoir ν, n á IM ñ= ò dx dv [R10n (x ) p (x, v, 0, t ) - R01n (x ) p (x, v, 1, t )], (32) and ò dx dvx [R10 (x ) p (x, v, 0, t ) - R01(x ) p (x, v, 1, t )] n n n á xIM ñ= (33) represents the position-current correlation. The average chemical work is given by á W˙ chemñ = å mn á IMn ñ , (34) n and the heat current entering from the reservoir of the oscillator is6 ⎛ 1 ⎞ á Q˙ oscñ = - g ⎜á v 2ñ - ⎟. (35) ⎝ mb osc ⎠ The latter equation is formally equivalent to the definition of heat flow for underdamped Langevin dynamics [33]. Similarly, the definition of the chemical work flow (see equation (34)) is consistent with the corresponding definition for the SET (see, e.g. [84] and [37] and references therein). However, the definition of heat with respect n to left and right leads (see equation (31)) differs by the additional contribution of -aV á xIM ñ, which stems from the interaction of QD and oscillator. Note that the above definitions of average chemical work flow and average heat flows can also be derived by use of the FPE, equation (1). To establish that the second law holds, i.e. that the average total entropy production rate is non-negative, we consider the evolution of the Shannon entropy ò S (t ) = - dx dv å p (x , v , q , t ) ln p (x , v , q , t ) , q (36) where p(x, v, q, t) is the solution of the FPE, equation (1). Taking the time derivative of S(t), introducing the shorthand notation p(q) ≡ p(x, v, q, t), p (q¢) º p (x , v , q¢, t ), and ⨋ º ò dx ò dv åq , and using the conservation of probability as well as partial integration (assuming vanishing boundary contributions, limx ¥ xp = lim v ¥ vp = 0 ) we obtain d n S (t ) = ⨋[¶v J (x , v , q , t )] lnp (q) - ⨋å R qq ¢ (x ) p (q ¢) lnp (q ) , (37) dt q¢n where g J (x , v , q , t ) = - vp (q) - D¶v p (q) (38) m is a probability current. Letting S˙1 (t ) and S˙2 (t ) denote the two terms on the right side of equation (37), we integrate by parts to rewrite the first term as follows: ⎧g [¶v p (q)]2 ⎫ S˙1 (t ) = ⨋ ⎨ v ¶v p (q) + D ⎬. (39) ⎩m p (q ) ⎭ 6 This can be seen by the connection v ◦ dB = (v + dv 2) · dB , where · refers to It-type calculus. 11 New J. Phys. 21 (2019) 073009 C W Wächtler et al From equation (35) we obtain ⎡ g ⎤ 0 = b osc á Q˙ oscñ + ⨋ ⎢b oscgv 2p (q) + v ¶v p (q)⎥. (40) ⎣ m ⎦ Summing equations (39) and (40) we arrive at S˙1 (t ) = b osc á Q˙ oscñ + S ˙ cont , (41) where ˙ cont = ⨋ [gvp (q) + Dm¶v p (q)]  0. 2 S (42) Dm 2p (q) Next, we rewrite the second term on the right side of equation (37) as follows: 1 n n S˙2 (t ) = - ⨋å [R qq ¢ p (q ¢) ln p (q ) + R q ¢ q p (q ) ln p (q ¢)]. (43) 2 q¢n From the property of (local) detailed balance obeyed by the electron tunnelling rates (see equation (2)), i.e. n R 01 n n n = e b (e - aVx - m ), (44) R10 we derive the identity n 1 Rq¢q 0= å b n á Q˙ nñ - 2 q¢n n ⨋å [R qq n ¢ p (q ¢) - R q ¢ q p (q )] ln n , R qq (45) n ¢ where the first term on the right relates to heat exchange with the fermionic leads (see equations (31)–(33)). Summing equations (43) and (45) and rearranging terms, we obtain S˙2 (t ) = å b n á Q˙ nñ + S˙ disc , (46) n where n R qq ¢ p (q ¢) ˙ disc = 1 ⨋å [R qq S n n ¢ p (q ¢) - R q ¢ q p (q )] ln  0. (47) n 2 q¢n R q ¢ q p (q ) Here, non-negativity follows from the log-sum inequality. Adding equations (41) and (46), we find that the rate of change of the system’s Shannon entropy is given by d S (t ) = S˙e + S ˙, (48) dt with S˙e = b osc á Q˙ oscñ + å b n á Q˙ nñ (49) n S ˙ =S ˙ disc  0. ˙ cont + S (50) Here, the entropy flow rate S˙e is the rate at which the total entropy of the reservoirs decreases due to heat exchange with the system [85]. The quantity ˙ = d S - b osc á Q˙ oscñ - S å b n á Q˙ nñ (51) dt n is the total entropy production rate, which can be expressed as the sum of two independently non-negative contributions (equation (50)), from the continuous (equation (42)) and discrete (equation (47)) degrees of freedom. The non-negativity of Ṡ, equation (50), shows that the second law holds in our system. We note that the two separate parts of the total entropy production rate (see equations (42) and (47)) are formally equivalent to the definitions derived for an independent underdamped harmonic oscillator and an d independent SET [37, 84]. However, as is apparent especially at steady states, where dt S (t ) = 0, the entropy production rate involves the heat flows (see equation (51)), for which the definitions for the independent systems differ from the definition for the electron shuttle (see equations (31) and (35)). 12 New J. Phys. 21 (2019) 073009 C W Wächtler et al 4.2. MF thermodynamics At the mean field level, the total energy of the system (denoted by an overbar) corresponding to equations (13)– (15) is given by mv¯ 2 kx¯ 2 E¯ = + + eq¯ - aVxq ¯ ¯, (52) 2 2 and its rate of change is dE¯ dv¯ = mv¯ + kxv ¯ ¯ + eI¯M - aVxI ¯ ¯M - aVqv ¯ ¯, (53) dt dt where dq¯ I¯M = dt = å (R10n p¯0 - R01n p¯1) º å I¯Mn . (54) n n n n Note that I¯M > 0 denotes the flow of matter from reservoir ν into the system, and vice versa for I¯M < 0. Using equations (13)–(15), the first law of thermodynamics takes the form dE¯ = å (e - aVx¯ - mn ) I¯M n n + mn I¯M - gv¯ 2 dt n L R osc = Q¯˙ + Q¯˙ + W¯˙chem + Q¯˙ . (55) Here n n Q¯˙ = (e - aVx¯ - mn ) I¯M (56) is the heat flow into the QD from reservoir ν, å mnI¯Mn chem W¯˙ = (57) n is the rate at which chemical work is performed on the system, and osc Q¯˙ = - gv¯ 2 (58) ˙ osc is the heat flow into the oscillator from its bath. Note that Q¯ is always negative, in contrast to the stochastic case (see equation (35)). To establish the second law within the MF approximation, we must first define the system entropy at this level of description. In the MF equations of motion, the state of the harmonic oscillator (x¯ , v¯) evolves deterministically under equations (13)–(14), while the QD is represented by a probability distribution p¯ = ( p¯0 , p¯1 )T evolving under a master equation, equation (15). We therefore define the entropy of the system to be the Shannon entropy of the QD probability distribution: S¯ (t ) = - å p¯q (t ) ln p¯q (t ). (59) q As in section 4.1, the total entropy production rate is the sum of the rates of change of the entropies of the system and the reservoirs: ˙¯ = d S¯ - b oscQ¯˙ osc - n S dt å b nQ¯˙ . (60) n We now analyse the three terms on the right side of this equation. The rate of change of the system entropy is given by d ¯ n p¯ S = - å å R qq ¯q ¢ lnp¯q = ¢p å (R10n p¯0 - R01n p¯1) ln p¯0 . (61) dt n q, q ¢ n 1 By equation (58), the second term on the right side of equation (60) is equal to gb oscv¯ 2. To analyse the third term we use equations (44), (54) and (56) to write n n n n n R 01 b nQ¯˙ = b n (e - aVx¯ - mn ) I¯M = (R10 p¯0 - R 01 p¯1) ln n . (62) R10 Combining results, we obtain R n p¯ ¯˙ = S å (R10n p¯0 - R01n p¯1) ln R10n p¯0 + gb oscv¯ 2  0, (63) n 01 1 in agreement with the second law. 13 New J. Phys. 21 (2019) 073009 C W Wächtler et al 4.3. Perturbative thermodynamics based on MFs MF perturbation theory gives the following result for the shuttle probability density (see section 3.3): p (x , v , q , t ) » pq (x ) p˜ (x , v , t ) , (64) where πq(x) denotes the instantaneous equilibrium distribution of the QD and p˜ (x , v , t ) is the solution to equation (17). Equation (64) represents an approximate solution of the FPE, equation (1). Therefore, to compute thermodynamic quantities such as heat flows and chemical work at this level of approximation, we use equation (64) to evaluate the relevant averages introduced in section 4.1. Note that an alternative attempt, which we will not follow here, could be to derive the laws of thermodynamics based on the effective FPE, equation (17). This effective FPE goes, however, beyond a simple adiabatic approximation and hence, its associated entropy production rate does not match the original entropy production rate evaluated with the approximated solution [86]. To evaluate the heat flow into the oscillator, equation (35), we first write á v 2ñ = ò dx dv åq v 2p (x, v, q, t ) » ò dx dv v 2p˜ (x, v, t ) = áv 2ñMS , (65) where á •ñMS denotes an average taken with the density p˜ (x , v , t ). Using the transformation to energy  and oscillation phase θ given by equation (22), and assuming pˆ ( , q, t ) » pˆ ( , t ) (see section 3.3), we get 1 1 ò á v 2ñMS = d  pˆ ( , t ) = á ñMS , m m (66) after averaging over θ. Here, pˆ ( , t ) is the solution to equation (23), which at steady state is given by equation (25). n For the chemical work and the heat exchanges with fermionic reservoirs, we need the matter currents á IM ñ n and position-current correlations, á xIMñ. Substituting equations (64) into (32) we get n n n á IM ñ » á [R10 (x ) p0 (x ) - R 01 (x ) p1 (x )]ñMS . (67) Transforming to ( , q )-space and averaging over θ gives n n ò dRˆ ( ) pˆ (, t ) = áRˆ ( )ñMS n á IM ñMS = (68) with 1 2p n Rˆ ( ) = 2p ò0 n [R10 n (x ) p0 (x ) - R 01 (x ) p1 (x )] dq. (69) For the position-current correlations, we similarly get n á xIM n ( )ñMS ñ » á xR (70) with 2p  n ( ) = 1 xR 2p ò0 n x [R10 n (x ) p0 (x ) - R 01 (x ) p1 (x )] dq. (71) Combining results with equations (31), (34) and (35) we obtain n n ñMS , á Q˙ n ñMS = (e - mn ) á Rˆ ñMS - aV á xR (72) n á W˙ chemñMS = å mn á Rˆ ñMS , (73) n g⎛ 1 ⎞ á Q˙ oscñMS = - ⎜á ñMS - osc ⎟. (74) m⎝ b ⎠ As in section 4.1, the first law is expressed by equation (30), but the heat flows and chemical work are now given by equations (72)–(74). Using equation (64), the system entropy (equation (36)) becomes a sum of distinct contributions from the harmonic oscillator and the QD: S (t ) » á- ln p˜ + S QD ñMS , (75) where SQD (x ) = -åq pq (x ) ln pq (x ). The decomposition S˙ = S˙e + S ˙ (see equation (48)) remains valid in the MS approximation. To show that the entropy production rate Ṡ is non-negative at this level of approximation, we first look at the contribution from the continuous degrees of freedom. Replacing p(x, v, q, t) by pq (x ) p˜ (x , v , t ) in equation (42), we find 14 New J. Phys. 21 (2019) 073009 C W Wächtler et al Figure 5. Entropy production rate for the full stochastic description (equation (50), orange solid), the MF approximation (equation (63), red dashed–dotted), the MS analysis (equation (81), blue dotted), and for the SET (thin green, equation (47) with a  0, l  ¥), as a function of bias voltage. The MF entropy production shows a clear signature of the underlying bifurcation at the critical value bV¯cr = 15.0 . The MS entropy production rate deviates from (specifically, underestimates) the full entropy production rate as V is increased, due to the breakdown of time scale separation. [gvpq p˜ + Dm¶v (pq p˜ )]2 [gvp˜ + Dm¶v p˜ ]2 S ˙ cont » ò dx dv åq Dm 2pq p˜ = ò dx dv Dm 2p˜  0. (76) Transforming to ( , q )-space and averaging over θ then results in [g  pˆ + D  m 2¶ pˆ ]2 S ˙ cont,MS = ò d Dm3pˆ  0. (77) Applying a similar analysis for the discrete degrees of freedom (see equation (47)) we get n R qq ¢ pq ¢ ˙ disc » 1 ⨋å (R qq S n p ¢ q¢ - R n p q¢q q ) ln p˜ , (78) n 2 q¢n R q ¢ q pq which after transforming to ( , q )-space and averaging over θ becomes ˙ disc,MS = 1 d s ( ) pˆ ( , t )  0, S 2 ò (79) where n 1 2p R qq ¢ pq ¢ ò0 n n s ( ) = dq å (R qq ¢ pq ¢ - R q ¢ q pq ) ln . (80) 2p qq ¢ n R qn¢ q pq Ṡdisc,MS is non-negative since s ( )  0 by the log-sum inequality, and pˆ ( , t )  0 is a probability density. Combining results, we get S ˙ MS = S ˙ disc,MS  0, ˙ cont,MS + S (81) again in agreement with the second law. 4.4. Discussion Having derived the general laws of thermodynamics at the different levels, we now discuss the thermodynamic properties of the model at hand and compare the stochastic, MF and MS solutions, focusing on the transition towards self-oscillation. All numerical results are obtained using the parameters specified in appendix D. 4.4.1. Entropy production rate We first look at the steady state entropy production rates as a function of the applied bias voltage V. Figure 5 shows the entropy production rate for the stochastic case (orange solid), the MF case (red dashed–dotted) and the MS description (blue dotted). Also shown is the SET entropy production rate ṠSET (green thin), i.e. equation (47) for the case where the position x is fixed at the origin [37, 84] or, alternatively, in the completely decoupled limit a  0, l  ¥ (see section 3.1). The steady state entropy production rate at the MF and the stochastic level is equal to the chemical work, aside from a factor β. To see this in the stochastic case, note that in the steady state the system’s Shannon entropy and average energy are constant: dS dt = 0 = ádE dt ñ. Combining this observation with equations (30), (34) and (51), and with our choice of setting all reservoir temperatures to be equal, b n = b osc = b , we obtain 15 New J. Phys. 21 (2019) 073009 C W Wächtler et al S ˙ = b á W˙ chemñ = bV á IM L ñ = - bV á IM R ñ. (82) The last equality follows from the preservation of electron number, á IMRñ + á IMLñ = 0 . In the MF case, we arrive at the analogous result using equations (55), (57) and (60). At the MF level, we see in figure 5 that below bV¯cr = 15.0 the entropy production rate S̄ ˙ is essentially the same as for the SET. This is understandable: below the bifurcation, the system evolves to a stable fixed point, with the QD at rest near the origin (see section 3.4.2). Above the bifurcation the shuttle oscillates, hence the resulting entropy production deviates from that of the SET. The sharp transition to self-oscillation is clearly reflected in the deviation of S̄˙ from Ṡ for V > V¯ . SET cr Interestingly, we see in figure 5 that self-oscillation lowers the rate of entropy production, relative to the value it would have taken had the QD remained at rest; that is, S ˙¯ < S˙ SET for V > V¯cr . In effect, above the critical voltage, when faced with a ‘choice’ between two modes of behaviour—oscillatory or at rest—the shuttle adopts the one that generates entropy more slowly. To understand this point quantitatively, note that the entropy production rate is determined by the matter current flowing from left to right through the device, equation (82). For βV?1 the SET current approaches IMSET = G 2, as our choice of chemical potentials produces a SET steady state in which p0=p1=1/2. For the MF case, recall that for V  V¯cr our system approaches the perfect shuttling regime in which one electron is transferred per oscillation period, which implies I¯M = w 2p , where ω is the oscillation frequency. Our parameter choices give IMSET = 0.5 and I¯M » 0.1, hence the SET generates entropy at a rate about five times that of the MF shuttle, in the high-bias limit. These results are consistent with the asymptotic slopes of the SET and MF entropy production rates shown in figure 5. By the same token, if the parameters were chosen such that I¯M > IMSET (roughly, if w = k m > p G), then we would get S ˙¯ > S ˙ SET . In the stochastic case the oscillator undergoes thermal motions in the steady state, the matter current through the QD is lower than the corresponding SET current, and as a result S ˙ V¯cr the key assumption of time scale separation is no longer valid and the perturbative solution breaks down. This is seen in the large deviation of ṠMS from Ṡ in figure 5. 4.4.2. Heat flows Next we look at the heat flows between the shuttle and the three reservoirs at steady state. At steady state the average energy of the system is constant at all levels, i.e. á dE dt ñ = dE¯ dt = á dE dt ñMS = 0, as we have verified numerically. In figure 6 we show the left and right steady state heat flows. At all three levels of description the heat flows between the QD and the left and right reservoirs are nearly indistinguishable. This is a consequence of our parameter choices of small α and mn = e  V 2, as can be seen from equation (31): á Q˙ n ñ = (e - mn ) á IM n ñ - aV á xIM n ñ. By the conservation of the matter, the first term on the right is the same for n the left and the right reservoir and differences arise only from the correlation of x and IM . Since α=0.06=1, the difference is barely noticeable in figure 6. Note that all heat flows are negative while the chemical work is positive, indicating that chemical work is performed on the system, and energy is transferred as heat into all three reservoirs. As with the case of entropy production (figure 5), at the MF level the onset of oscillations at bV¯cr is clearly n osc reflected in the heat flows Q¯˙ (figure 6) and Q¯˙ (figure 7). The latter vanishes below the bifurcation (where the shuttle is at rest), but becomes negative above the bifurcation, where the shuttle’s oscillatory motion gives rise to dissipation due to friction. Below the bifurcation, the MF heat flows agree with the corresponding SET values, as expected. In the stochastic case the transition to self-oscillation is gradual rather than sharp, as seen in the behaviours of á Q˙ Lñ, á Q˙ Rñ and á Q˙ oscñ. Away from the transition—that is, for small and large values of V—these heat flows are well approximated by the MF values, as was the case for the entropy production rate. For the MS solution we see that á Q˙ n oscñ » á Q˙ n oscñMS at small values of V. At larger values of the applied bias voltage the MS heat flows deviate from the full stochastic heat flows, due to the breakdown of time scale separation as previously discussed. 16 New J. Phys. 21 (2019) 073009 C W Wächtler et al Figure 6. Heat flows between the fermionic reservoirs and the QD as a function of the bias voltage V: Solid lines correspond to heat flow of the left reservoir, dotted lines to heat flow of the right reservoir. We plot the full stochastic heat flows (orange, equation (31)), the MF heat flows (red, equation (56)) as well as the heat flows derived by MS perturbation theory (blue, equation (72)). For comparison we also plot the heat flows of the SET (thin green). The MF heat flows show a clear signature of the underlying bifurcation whereas the full stochastic and the MS heat flows transition smoothly between the two regimes of operation. Figure 7. Heat flow between the oscillator and its respective heat bath as a function of the bias voltage V: the solid (green) line corresponds to fully stochastic heat flow see equation (35), the (red) dashed–dotted line shows the MF heat flow (see equation (58)) and the dotted (blue) line shows the effective heat flow by use of MS (see equation (74)). Again the MF heat flow shows a clear signature of the underlying bifurcation whereas the full stochastic and the MS heat flow transition smoothly between the two regimes of operation. To summarize this section, the thermodynamic quantities we have studied—the entropy production rate and the heat exchanges with reservoirs—all bear the signature of the onset of self-oscillation. In particular, all of these quantities deviate substantially from the corresponding SET values above the critical voltage V¯cr , as the system approaches the pure shuttling regime and its oscillatory motion influences the exchange of energy and matter. The transition is abrupt in the mean field approximation, but gradual in the fully stochastic case. 5. Conclusion and future applications We have provided a classical stochastic description of the single electron shuttle based on coupled Langevin equations including thermal and Poissonian noise terms The average dynamics can be well approximated by MF equations away from the onset of self-oscillations for our specific choice of parameters. However, we expect deviations between the stochastic and MF description if fluctuations are strong, i.e. for a system with a smaller mass. Within the MF approximation the system undergoes a Hopf bifurcation from a stable fixed point to a stable limit cycle by changing the applied bias voltage, where a limit cycle corresponds to self-oscillation of the shuttle. By introducing a time scale separation between the frequent tunnelling events and the slow dynamics of the oscillator we were able to perturbatively solve the FPE corresponding to the coupled Langevin equations. This MS solution approximates the full stochastic solution very well below and around the critical bias voltage of the 17 New J. Phys. 21 (2019) 073009 C W Wächtler et al MF description. However, as the applied bias voltage is increased and the system starts shuttling, the underlying assumption of many jumps per cycle is not valid anymore and the perturbation approach breaks down. An order parameter defined in terms of a mean amplitude of oscillation shows clear signatures of the Hopf bifurcation found in the deterministic MF description. This order parameter is not sensitive to small oscillations of the shuttle in the stochastic and MS approaches, hence the system transitions smoothly towards self-oscillation in those cases, and the abrupt onset is realized only in the MF case. In the classical description chosen here, we identify three different regimes: below the bifurcation, i.e. for V < V¯cr the shuttle acts as a SET with additional noise. For very large bias voltages V  V¯cr the system oscillates and transports one electron from the reservoir with higher chemical potential to the reservoir with lower chemical potential per period. In this regime the system truly serves as a shuttle for single electron transport. Between the two limits there is a crossover regime. Additionally, we performed a thermodynamic analysis of the shuttling mechanism. Using stochastic thermodynamics we derived the first and second law at the stochastic as well as the MF level. At the perturbative level, we used the solution from MS perturbation theory to perform ensemble averages in order to define effective thermodynamic quantities. The thermodynamic quantities as entropy production rate, heat flows and chemical work rate show clear signatures of the underlying bifurcation within the MF approximation. The corresponding stochastic and MS quantities lack such an abrupt transition from noisy movement to self- oscillation, but show a smoothed transition, which suggest that the abrupt transition seen in the dynamical description is an artefact of the chosen order parameter. However, the thermodynamic quantities do reflect the transition towards pure shuttling if compared to the SET. The thermodynamic analysis of the electron shuttle provides an exemplary discussion of the thermodynamics of self-sustained oscillations, especially for highly fluctuating systems at the nano-scale, and which is also experimentally realizable. The combined system of QD and oscillator has rich dynamics and is capable of realizing various thermodynamic objectives. Our thermodynamic analysis from different perspectives may be used as the starting point to analyse the performance of the shuttle as a heat pump or refrigerator, with applied chemical work used to transfer heat from cold fermionic reservoirs to the hot reservoir of the oscillator (β ν>βosc). A quantum heat engine using the electron shuttle has recently been discussed in [88]. Alternatively, the shuttle can be transformed into an engine, which uses the chemical gradient in order to perform mechanical work. Such an engine might be constructed as a nanoscale rotor driven by single electron tunnelling, as proposed in [89] and [31]. Similarly, one can then use the mechanical motion in order to pump electrons and generate a matter current against an externally applied electric field [31]. The thermodynamic analysis of such a device is the subject of our further research in this direction. With the present work, we hope to pave the way for proper thermodynamic analyses for realistic engines based on the concept of the electron shuttle and stimulate further discussion about the thermodynamic possibilities of such devices. Acknowledgments The authors thank T Brandes for initiating this project. CW acknowledges fruitful discussions with S Restrepo and J Cerrillo. The authors are thankful for stimulating discussions with H Engel. This work has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Projektnummer 163436311— SFB 910 and the Graduate Research Training Group RTG 1558. PS acknowledges financial support by the European Research Council project Nano Thermo (ERC-2015-CoG Agreement No. 681456) and from the DFG through project STR 1505/2-1. C J was supported by the US National Science Foundation under grant DMR- 1506969. Appendix A. Equivalence of FPE and stochastic differential equations In this section we will show that the FPE, equation (1), and the stochastic differential equations (3)–(5) describe the same system, i.e. both representations reproduce the same expectation values up to order (dt ). Taking an arbitrary differentiable scalar function f one finds for the expectation value of f (q, x, v) by employing the FPE, equation (1), (assuming that boundary contributions vanish, i.e. f (q, x, v)p(x, v, q, t) → 0 as x  ¥ or v  ¥) 18 New J. Phys. 21 (2019) 073009 C W Wächtler et al ¶á f (q , x , v )ñ ⎡ ¶ ¶ ⎛k g aV ⎞ ¶2 ⎤ = ⨋f (q , x , v ) ⎢ - v + ⎜ x+ v- q⎟ + D 2 ⎥ p ( x , v , q , t ) ¶t ⎣ ¶x ¶v ⎝ m m m ⎠ ¶v ⎦ n n + ⨋å [ f (q , x , v ) R qq ¢ (x ) p (x , v , q ¢ , t ) - f (q , x , v ) R q ¢ q (x ) p (x , v , q , t )] q¢n ⎡ ¶f (q , x , v ) ¶f (q , x , v ) ⎛ k g aV ⎞ ¶ 2f (q , x , v ) ⎤ = ⨋⎢ v- ⎜ x+ v- q⎟ + D ⎥p ⎣ ¶x ¶v ⎝m m m ⎠ ¶v 2 ⎦ + ⨋å [ f (q¢ , x , v ) - f (q , x , v )] R qn¢ q (x ) p (x , v , q , t ) q¢n ¶f ¶f ⎛ k g aV ⎞ ¶ 2f = ¶x v- ⎜ ¶v ⎝ m x+ v- m q⎟ + m ⎠ ¶v 2 D + å [ f (q¢) - f (q)] Rqn¢q , (A.1) q¢n where ⨋ º ò dx ò dv åq . Our aim is to show that equation (A.1) is also obtained by the use of the stochastic differential equations (3)–(5). In terms of the stochastic process defined by equations (3)–(5) we can write for the expectation value of the increment of the arbitrary scalar function f the following (using Itô’s Lemma): ⎡ ¶f ⎤ ⎡ ¶f ⎛ k g aV ⎞ ⎤  [ f (q + dq , x + dx , v + dv ) - f (q , x , v )] =  ⎢ v ⎥ dt -  ⎢ ⎜ x + v - q ⎟ ⎥ dt ⎣ ¶x ⎦ ⎣ ¶v ⎝ m m m ⎠⎦ ⎡ ¶ 2f ⎤ ⎡¥ 1 ⎤ +  ⎢ 2 D⎥ dt +  ⎢ å f (k) (dq)k ⎥ +  (dt 2) , (A.2) ⎣ ¶v ⎦ ⎣k = 1 k! ⎦ where we have used the statistical properties of the Wiener increment, i.e.  [dB (t )] = 0 and  [(dB (t ))2] = dt . Note that mixed terms of dB(t) and dNqn¢ q (x , t ) exceed the leading order of dt because  [dNqn¢ q (x , t )] µ dt . Here, [•] denotes averages of the stochastic process. We now evaluate the sum in equation (A.2): since all powers of the Poisson increment are of order dt, we have to evaluate the sum exactly. Since [dNqn¢ q (x , t )]k = dNqn¢ q (x , t ) for all k Î  we can rewrite the expectation value as follows (omitting any dependencies on x and v): ⎡ ¥ f (k) (q) ⎤ ⎡ ¥ f (k) (q) ⎤ ⎡ ⎤ ⎢ å ( dq ) k ⎥ =  ⎢ å å (q¢ - q)k dNqn¢ q ⎥ =  ⎢å { f (q¢) - f (q)} dNqn¢ q ⎥ ⎣k = 1 k! ⎦ ⎢⎣k = 1 k! q ¢ n ⎥⎦ ⎢⎣ q ¢ n ⎥⎦ ⎡ ⎤ =  ⎢å { f (q¢) - f (q)} R qn¢ q ⎥ dt , (A.3) ⎢⎣ q ¢ n ⎥⎦ where the last equality follows from a general identity of point/Poisson processes (see, e.g., equation (B.54) in [90]). Hence, up to order (dt ) we write ¶ ⎡ ¶f ⎤ ⎡ ¶f ⎛ k g aV ⎞ ⎤ ⎡ ¶ 2f ⎤  [ f (q , x , v )] =  ⎢ v ⎥ -  ⎢ ⎜ x + v - q⎟ ⎥ +  ⎢ 2 D ⎥ ¶t ⎣ ¶x ⎦ ⎣ ¶v ⎝ m m m ⎠⎦ ⎣ ¶v ⎦ ⎡ ⎤ +  ⎢å [ f (q¢) - f (q)] R qn¢ q ⎥. (A.4) ⎢⎣ q ¢ n ⎥⎦ Since equation (A.1) is equivalent to equation (A.4) we can conclude that, for the same initial conditions, expectation values of an arbitrary function f with respect to the probability density and with respect to realizations of the stochastic process evolve equally. Hence, the FPE, equation (1), and the stochstic differential equations (3)–(5) describe the same process. Appendix B. MF perturbation theory In this section we derive equation (17) from the full FPE, equation (1). The idea of MS perturbation theory is to impose a time scale separation of the frequent electron tunnelling events and the slow evolution of the oscillator and, furthermore, to demand that those terms of the approximated solution that grow with time, vanish. By imposing the latter condition we ensure that the MS solution of the full probability density will be valid on the long time scale. First we will state some useful properties of the matrix R (x ) = ån Rn (x ) (see equation (2)), which we will use throughout the derivation. Since R(x) is a 2×2 rate matrix it has two eigenvalues: 0 and χ<0. Accordingly, there are two (right) eigenvectors, p (x ) and c , for which R (x ) p (x ) = 0 and R (x ) c = cc holds, respectively. Here, 19 New J. Phys. 21 (2019) 073009 C W Wächtler et al ⎛ p (x ) ⎞ p (x ) = ⎜ 0 ⎟ (B.1) ⎝ p1 (x ) ⎠ is the (instantaneous) stationary solution of R(x). Note that since R=R(x), the stationary state p (x ) is also a function of the position x and we impose the normalization condition π0(x)+π1(x)=1. Furthermore, c= (-11) (B.2) and ⎛x⎞ c = - 2G cosh ⎜ ⎟. (B.3) ⎝l ⎠ There are additionally left eigenvectors p† = (1, 1) (B.4) and c† (x ) = (p1 (x ) , - p0 (x )) , (B.5) satisfying p†R (x ) = 0 and c† (x ) R (x ) = cc† (x ). We start the derivation of equation (17) by considering equation (1) in matrix representation, i.e. ¶p ⎡ ¶ ¶ ⎛k g ⎞⎤ aV 0 0 ¶p ¶ 2p ¶t = ⎢-v ⎣ ¶x + ⎜ ¶v ⎝ m x + v⎟ ⎥ p + m ⎠⎦ m 0 - 1 ¶v ¶v ( ) + D 2 + R (x ) p , (B.6) where p = ( p0 (x , v , t ), p1 (x , v , t )). Introducing the differential operator  ⎡ ¶ ¶ ⎛k g ⎞ ⎤ aV 0 0 ¶ ¶2 L = ⎢-v ⎣ ¶x + ⎜ ¶v ⎝ m x + v⎟ ⎥ + m ⎠⎦ m 0 - 1 ¶v + D 2. ¶v ( ) (B.7) Equation (B.6) takes the compact form ¶p (x , v , t ) = p (x , v , t ) + R p (x , v , t ). (B.8) ¶t Since the fast time scale describes the dynamics of the two state system, we treat the  -part as perturbation and introduce the bookkeeping parameter ò = 1 such that we can rewrite equation (B.8) as ¶p =  p + R p. (B.9) ¶t The idea of MS perturbation theory is now to introduce two time scales, a fast one (t˜ ) and a slow one t =  t˜ such that the probability density is a function of both times scales, i.e. p˜ (x , v , t˜, t ) = p (x , v , t ). The temporal derivative transforms to a sum provided by the chain rule: ¶ ¶ ¶ = + . (B.10) ¶t ¶t˜ ¶t Then, equation (B.9) is given by ¶p˜ ¶p˜ + =  p˜ + R p. ˜ (B.11) ¶t˜ ¶t From this point on, we will refer to t˜ as t and to p˜ (x , v , t , t ) as p (t , t ). Assuming that we can express p as a series of orders of ò, p (t , t ) = p(0) (t , t ) +  p(1) (t , t ) +  2 p(2) (t , t ) +  ( 3) , (B.12) we find a hierarchy of equations for the different orders of ò. The goal of the MS perturbation theory is now to find an approximate solution, such that, after setting ò to 1, it holds p (t , t ) » p(0) (t , t ) + p(1) (t , t ). (B.13) We start with the governing equation for ( 0): ¶p(0) (t , t ) = R p(0) (t , t ). (B.14) ¶t The simplest solution of the ordinary differential equation (B.14) is given by assuming that the left hand side of equation (B.14) is equal to 0, i.e. assuming that the probability density at zeroth order is independent of the fast time scale t: p(0) (t , t ) = p(0) (t ). This means p(0) (t ) must be the eigenvector of R with eivenvalue 0, i.e. 20 New J. Phys. 21 (2019) 073009 C W Wächtler et al p(0) (t ) = p (x ) p(0) (t ). (B.15) Here, p(0) (t ) is a scalar function which represents the probability density of the oscillator alone, i.e. tracing out the charge state q of p(0) (x , v , t ) results in p0 (x ) p(0) (x , v , t ) + p1 (x ) p(0) (x , v , t ) = p(0) (x , v , t ) , (B.16) which is the probability density to find the oscillator at position x with velocity v at time τ (at zeroth order). The specific form of p(0) (t ) will be determined by the first order of the perturbation hierarchy. Note that, if we stop the perturbation theory here, equation (B.15) implies an infinite time scale separation, which is equivalent to an adiabatic approximation. The equation of motion of p (t , t ) at ( ) is given by ¶p(1) ¶p(0) (t ) (t , t ) = - + p(0) (t ) + R p(1) (t , t ). (B.17) ¶t ¶t Since p(0) (x , v , t ) is independent of t (see equation (B.15)), the solution of the latter equation is formally given by t ⎡ ¶p(0) (t ) ⎤ p(1) (t , t ) = e Rt p˜ (1) (t ) + e Rt ò0 e-Rsds ⎢ - ⎣ ¶t + p(0) (t )⎥ , ⎦ (B.18) where p˜ (1) (t ) is a probability vector which does not depend on the fast time scale t but is unspecified at this moment. Next we look at the integral of equation (B.18): we can expand the exponential of the rate matrix by use of the eigenvalues and eigenvectors of R (see equations (B.1)–(B.5)), i.e. t t ò0 e-Rsds = ò0 p (x ) p† + e-cscc†ds. (B.19) Evaluating the integral gives t e-ct - 1 † ò0 e-Rsds = p (x ) p†t - c cc . (B.20) Inserting equations (B.20) into (B.18) we find that there are terms in the solution which grow linearly with t for long times, i.e. ⎡ ¶p(0) ⎤ e Rt p (x ) p†t ⎢ - + p(0)⎥. (B.21) ⎣ ¶t ⎦ Those terms, which will subsequently be referred to as secular terms, prohibit a steady state solution of the perturbation hierarchy. We therefore demand the secular terms to vanish such that we find a stable solution. At the first order perturbation the latter condition is satisfied if ⎡ ¶p(0) ⎤ p† ⎢ - + p(0)⎥ = 0. (B.22) ⎣ ¶t ⎦ Inserting p(0) (see equation (B.15)) yields ⎡ ¶p(0) ⎤ p† ⎢ - p (x ) + p (x ) p(0) ⎥ = 0. (B.23) ⎣ ¶t ⎦ It holds that p†p (x ) = 1 (see equations (B.1) and (B.4)). Therefore ¶p(0) - + p†p (x ) p(0) = 0. (B.24) ¶t Evaluating the action of the differential operator  on p (x ) results in ¶ ¶ ⎛k g ⎞ aV ¶ ¶2 p†p (x ) = - v + ⎜ x + v⎟ - qeq (x ) + D 2, (B.25) ¶x ¶v ⎝ m m ⎠ m ¶v ¶v where qeq (x ) = p1 (x ). Furthermore, we have used that (see equations (B.1) and (B.4)) ¶ ¶ † ¶ p† p (x ) = p p (x ) = 1 = 0. (B.26) ¶x ¶x ¶x The condition for secular terms to vanish in the first order perturbation (see equation (B.22)) is now given by ¶p(0) ¶ ¶ ⎡k g aV ⎤ ¶2 = - v p(0) + ⎢⎣ x + v - qeq (x ) ⎥ p(0) + D 2 p(0) . (B.27) ¶t ¶x ¶v m m m ⎦ ¶v The latter equation is a FPE for p(0) (x , v , t ) describing the underdamped evolution in an effective potential Ueff(x) with ¶x Ueff = kx - aVqeq (x ). This corresponds to our ansatz for the 0th order, where we have assumed 21 New J. Phys. 21 (2019) 073009 C W Wächtler et al that the QD is in its instantaneous equilibrium state at all times. The harmonic potential is therefore altered and the effective FPE describing the oscillator is simple diffusion within the effective potential. We now return to the perturbation of ( ) (equation (B.18)): after removing the secular terms, the first order solution is given by 1 - e-ct ⎡ ¶p(0) ⎤ p(1) = e Rt p˜ (1) + e Rt cc† ⎢- + p(0)⎥. (B.28) c ⎣ ¶t ⎦ The exponential of the rate matrix R can again be expressed in terms of the eigenvectors, e Rt = p (x ) p† + e ct cc†. (B.29) Inserting equations (B.29) into (B.28) results in 1 - e-ct ⎡ ¶p(0) ⎤ p(1) = p (x ) p† p˜ (1) + e ct cc† p˜ (1) + p (x ) p†cc† ⎢- + p(0)⎥ c ⎣ ¶t ⎦ e ct - 1 ⎡ ¶p(0) ⎤ + cc† ⎢- + p(0)⎥ , (B.30) c ⎣ ¶t ⎦ where we have used that c†c = 1 (see equations (B.2) and (B.5)). In the long-time limit terms proportional to e χt will approach zero, since χ<0. Using the latter as well as the fact that p†c = 0, equation (B.30) simplifies to 1 ⎡ ¶p(0) ⎤ p(1) = p (x ) p(1) - cc† ⎢ - + p(0)⎥ , (B.31) c ⎣ ¶t ⎦ where p(1) º p† p˜ (1). Substituting now p(0) (t ) by p (x ) p(0) (t ) (see equation (B.15)) and the differential operator  by its definition (see equation (B.7)) we can approximate the first order solution (equation (B.28)) by 1 ⎧ ¶p (x ) aV ⎜⎛ 0 ⎟⎞ ¶ ⎫ (0) p(1) = p (x ) p(1) - cc† ⎨ - v + ⎬p , (B.32) c ⎩ ¶x m ⎝- p1 (x )⎠ ¶v ⎭ using c†p (x ) = 0 (see equations (B.1) and (B.5)). We now define a new vector ⎧ ¶p (x ) aV ⎜⎛ 0 ⎟⎞ ¶ ⎫ (0) k º ⎨-v + ⎬p (B.33) ⎩ ¶x m ⎝- p1 (x )⎠ ¶v ⎭ such that equation (B.32) becomes 1 p(1) = p (x ) p(1) - cc†k. (B.34) c Note that c†k ¹ 0 in general. Similar to the procedure for the zeroth order perturbation we now look at the second order of p and demand secular terms to vanish. By this condition we will find a differential equation for p(1), which describes the effective evolution of the oscillator at a first order perturbation level without taking the electronic degrees of freedom specifically into account. The governing equation at ( 2) is similar to equation (B.17) and reads ¶p(2) ¶p(1) =- + p(1) + R p(2). (B.35) ¶t ¶t Again, the general solution can be written as t ⎡ ¶p(1) ⎤ p(2) = e Rt p˜ (2) + e Rt ò0 e-Rsds ⎢ - ⎣ ¶t + p(1)⎥. ⎦ (B.36) With the same calculation as above we find that in order for the secular terms in the second order to vanish, the following condition must hold: ⎡ ¶p(1) ⎤ p† ⎢ - + p(1)⎥ = 0 ⎣ ¶t ⎦ ¶p(1) 1 - + p†p (x ) p(1) - p†cc†k = 0. (B.37) ¶t c The term p†p (x ) appears again and is given by equation (B.25). We now evaluate the third term on the left hand side of equation (B.37): first we note that p†c = 0 (see equations (B.2) and (B.4)). Therefore, we only have to evaluate the action of  on c in equation (B.37), which is 22 New J. Phys. 21 (2019) 073009 C W Wächtler et al aV 0 ¶ c = m 1 ¶v , () (B.38) because c is a constant vector (see equation (B.2)). It then holds (see equation (B.4) and (B.5)) 1 † aV † 0 ¶ † c p cc†k = cm p 1 ¶v (c k ) () aV † ¶ = c k. (B.39) cm ¶v Lastly, we look at the derivative of k with respect to v, that is ¶ ⎛¶ ⎞ ¶ aV 0 ¶ 2p(0) ¶v k = - ⎜ p (x ) ⎟ (vp(0)) + ⎝ ¶x ⎠ ¶v m - p1 ¶v 2 . ( ) (B.40) Using p1 = 1 - p0 = qeq , one can show that ¶p1 ¶p ¶qeq - p1 0 = p0 . (B.41) ¶x ¶x ¶x Furthermore it holds that p0 p1 = qeq (x ) - qeq2 (x ). With the latter two simplifications we can rewrite (see equation (B.5)): ⎛¶ ⎞ ¶qeq ¶ aV ¶ 2p(0) c† ⎜ k⎟ = (vp(0)) + qeq (1 - qeq ) . (B.42) ⎝ ¶v ⎠ ¶x ¶v m ¶v 2 Putting everything together, the condition for the secular terms of the solution at second order in the perturbation (see equation (B.37)) is given by ¶p(1) aV ⎡ ¶qeq ¶(vp(0)) aVqeq ¶ 2p(0) ⎤ =  0 p(1) - ⎢ + (1 - qeq ) ⎥, (B.43) ¶t cm ⎣ ¶x ¶v m ¶v 2 ⎦ where ¶ ¶ ⎛k g aV ⎞ ¶2  0 º -v + ⎜ x+ v- qeq⎟ + D 2 . (B.44) ¶x ¶v ⎝ m m m ⎠ ¶v Putting zeroth and first order together (see equations (B.27) and (B.43)), i.e. p˜ (x , v , t ) = v , t ) +  p(1) (x , v , t ) and setting ò to 1 (τ=òt → t), we find that the full probability density p(x, v, q, t) p(0) (x , can be approximated by p (x , v , q , t ) » pq (x ) p˜ (x , v , t ) , (B.45) where p˜ (x , v , t ) is the probability density of solely the oscillator, obtained by tracing out the fast electronic degrees of freedom. The dynamics of p(x, v, t) is governed by a FPE: ¶p˜ ⎡ ¶ ¶ ⎛k g˜ (x ) ⎞ aVqeq (x ) ¶ ⎤ ¶2 = ⎢-v + ⎜ x+ v⎟ - ⎥ p˜ + D˜ (x ) 2 p˜ , (B.46) ¶t ⎣ ¶x ¶v ⎝ m m ⎠ m ¶v ⎦ ¶v where the effective friction and diffusion coefficients are now position dependent and are given by aV ¶qeq (x ) g˜ (x ) = g - , c ¶x a 2V 2qeq (x ) D˜ (x ) = D - (1 - qeq (x )). (B.47) m 2c (x ) The final equation (B.46) is equivalent to equation (17) of the main text. Appendix C. Transformation to energy space In this section we derive the transformed FPE, equations (23), from (17). We first rewrite equation (17) as follows: ¶p˜ ¶p˜ g˜ (x ) ⎛k g˜ (x ) aV ⎞ ¶p˜ ¶ 2p˜ = -v + p˜ + ⎜ x + v- qeq (x ) ⎟ + D˜ (x ) 2 . (C.1) ¶t ¶x m ⎝m m m ⎠ ¶v ¶v Using the transformation of equation (22) as well as the assumption pˆ ( , q, t ) » pˆ ( , t ), we find that derivatives with respect to x and v transform as follows: 23 New J. Phys. 21 (2019) 073009 C W Wächtler et al ¶p˜ ¶pˆ  2k sin q , ¶x ¶ ¶p˜ ¶pˆ  2m cos q , ¶v ¶ ¶ 2p˜ ¶pˆ 2q¶ p . 2ˆ  m + 2 m cos (C.2) ¶v 2 ¶ ¶ 2 In energy space, equation (C.1) then takes the form ¶ ⎡ g˜ (x ) ⎛ g˜ (x ) 2 ⎞ ¶ pˆ ( , t ) = ⎢ +⎜ 2 cos2 q - aVqeq (x ) cos q + D˜ (x ) m⎟ ¶t ⎢⎣ m ⎝ m m ⎠ ¶ ¶2 ⎤ + D˜ (x ) 2m cos2 q ⎥ pˆ ( , t ) , (C.3) ¶ 2 ⎦ where x = 2 k sin q . Upon averaging over θ, the term proportional to cos q will vanish: as qeq(x) is an analytic function of x it can be Taylor expanded in a power series and the individual contributions of all terms 2p vanish due to ò sinn (q ) cos (q ) dq = 0 for all n Î  . By inspection of equation (24) and partial integration one 0 can show that the following relations hold: 1 2p ¶gˆ ( ) 2p ò0 dqg˜ (x ) = 2gˆ ( ) + 2 ¶ , 1 2p ¶Dˆ ( ) 2p ò0 dqD˜ (x ) = 2Dˆ ( ) + 2 ¶ . (C.4) Then, averaging equation (C.3) results in equation (23). Appendix D. Computational methods For all numerical investigations we set b osc = b n = b º 1 as well as Γ≡1 and λ≡1. The other parameters used in this work are given in units of the latter three: αλ=0.06, mλ2Γ2β=12.0, kλ2β=5.0 and γλ2Γβ=0.2. Since the probability space of the coupled system of harmonic oscillator and QD is very large, we assume that the system is ergodic, such that we can sample the steady state probability density of the system by a single long trajectory. Additionally this means that an ensemble average of an arbitrary quantity A in the steady state is calculated by 1 T á Añ = T ò0 A (t ) , (D.1) which is exact for ergodic systems in the limit of T  ¥. We simulate the trajectories after a relaxation time of Γt=1000 until ΓT=5000 000, where we have also checked that further relaxation time or simulation time does not change the probability density or averaged quantities. Note that we have also investigated different initial conditions and have not seen any dependency of the outcome on the initial conditions (after the relaxation time). Finally we note that the time step used in the simulations is ΓΔt=0.0001. Appendix E. Hopf bifurcation of the MF model In order to determine the critical value V¯cr , for which the electron shuttle bifurcates from a stable fixed point into a stable limit cycle, we perform a linear stability analysis around the fixed point of the MF equations (13)–(15): x¯˙ = v¯, k g aV v¯˙ = - x¯ - v¯ + p¯ , m m m 1 p¯˙1 = å R10 n (x¯)(1 - p¯1) - R 01n (x¯) p¯1 , (E.1) n where we have eliminated one equation compared to equations (13)–(15) by use of probability conservation, i.e. p¯0 + p¯1 = 1. The stability of the fixed point is determined by the eigenvalues of the Jacobian J0 of the right hand side of the latter equation evaluated at steady state (x¯˙ = 0, v¯˙ = 0, p¯˙1 = 0): If one or more eigenvalues have positive real part the fixed point is unstable. Figure E1 shows the real part (left) and imaginary part (right) of the three eigenvalues as a function of the applied bias voltage. There exists a pair of conjugate eigenvalues (solid red and 24 New J. Phys. 21 (2019) 073009 C W Wächtler et al Figure E1. Eigenvalues of J0 evaluated at the steady state as a function of the applied bias voltage. As V is increased the two complex conjugated eigenvalues (solid red and dotted blue) become purely imaginary at a critical value of bV¯cr = 15, which denotes the bifurcation point and the stability of the fixed point changes. Further increase of V results in a stable limit cycle and an unstable fixed point. dotted blue) which become purely imaginary at a critical value bV¯cr = 15.0. At this point the Hopf bifurcation sets in and the MF system undergoes the transition from a stable fixed point to a stable limit cycle. ORCID iDs Christopher W Wächtler https://orcid.org/0000-0002-0441-138X Christopher Jarzynski https://orcid.org/0000-0002-3464-2920 References [1] Jenkins A 2013 Phys. Rep. 525 167–222 [2] Van Der Pol B and Van Der Mark J 1928 Phil. Mag. 7 6 763–75 [3] Novák B and Tyson J J 2008 Nat. Rev. Mol. Cell Biol. 9 981 [4] Field R J and Noyes R M 1974 J. Chem. Phys. 60 1877–84 [5] Gorelik L, Isacsson A, Voinova M, Kasemo B, Shekhter R and Jonson M 1998 Phys. Rev. Lett. 80 4526 [6] Boese D and Schoeller H 2001 Europhys. Lett. 54 668 [7] Armour A and MacKinnon A 2002 Phys. Rev. B 66 035333 [8] McCarthy K D, Prokof’ev N and Tuominen M T 2003 Phys. Rev. B 67 245415 [9] Novotny` T, Donarini A and Jauho A P 2003 Phys. Rev. Lett. 90 256801 [10] Novotny` T, Donarini A, Flindt C and Jauho A P 2004 Phys. Rev. Lett. 92 248302 [11] Donarini A, Novotny` T and Jauho A P 2005 New J. Phys. 7 237 [12] Utami D W, Goan H S, Holmes C and Milburn G 2006 Phys. Rev. B 74 014303 [13] Fedorets D, Gorelik L, Shekhter R and Jonson M 2002 Europhys. Lett. 58 99 [14] Nocera A, Perroni C, Ramaglia V M and Cataudella V 2011 Phys. Rev. B 83 115420 [15] Isacsson A, Gorelik L, Voinova M, Kasemo B, Shekhter R and Jonson M 1998 Physica B 255 150–63 [16] Weiss C and Zwerger W 1999 Europhys. Lett. 47 97 [17] Nord T, Gorelik L, Shekhter R and Jonson M 2002 Phys. Rev. B 65 165312 [18] Park H, Park J, Lim A K, Anderson E H, Alivisatos A P and McEuen P L 2000 Nature 407 57 [19] Moskalenko A V, Gordeev S N, Koentjoro O F, Raithby P R, French R W, Marken F and Savel’ev S E 2009 Phys. Rev. B 79 241403 [20] Moskalenko A, Gordeev S, Koentjoro O, Raithby P, French R, Marken F and Savel’ev S 2009 Nanotechnology 20 485202 [21] König D R and Weig E M 2012 Appl. Phys. Lett. 101 213111 [22] Scheible D V and Blick R H 2004 Appl. Phys. Lett. 84 4632–4 [23] Kim C, Prada M, Qin H, Kim H S and Blick R H 2015 Appl. Phys. Lett. 106 061909 [24] Joachim C, Gimzewski J K and Aviram A 2000 Nature 408 541 [25] Shekhter R, Galperin Y, Gorelik L Y, Isacsson A and Jonson M 2003 J. Phys. Condens. Matter 15 R441 [26] Galperin M, Ratner M A and Nitzan A 2007 J. Phys. Condens. Matter 19 103201 [27] Galperin M, Ratner M A, Nitzan A and Troisi A 2008 Science 319 1056–60 [28] Shekhter R I, Gorelik L Y, Krive I V, Kiselev M, Parafilo A and Jonson M 2013 Nanoelectromech. Syst. 1 1–25 [29] Lai W, Zhang C and Ma Z 2015 Front. Phys. 10 59–86 [30] Strogatz S H 2018 Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (Boca Raton, FL: CRC Press) [31] Wächtler C, Strasberg P and Schaller G 2019 arXiv:1903.07500 [32] Esposito M, Harbola U and Mukamel S 2009 Rev. Mod. Phys. 81 1665 [33] Sekimoto K 2010 Stochastic energetics Lecture Notes in Physics (Berlin: Springer) [34] Campisi M, Hänggi P and Talkner P 2011 Rev. Mod. Phys. 83 771 [35] Jarzynski C 2011 Annu. Rev. Condens. Matter Phys. 2 329–51 25 New J. Phys. 21 (2019) 073009 C W Wächtler et al [36] Seifert U 2012 Rep. Prog. Phys. 75 126001 [37] Schaller G 2014 Open quantum systems far from equilibrium Lecture Notes in Physics (Cham: Springer) [38] Van den Broeck C and Esposito M 2015 Physica A 418 6–16 [39] Schmiedl T and Seifert U 2007 Europhys. Lett. 81 20003 [40] Blickle V and Bechinger C 2011 Nat. Phys. 8 143 [41] Rana S, Pal P S, Saha A and Jayannavar A M 2014 Phys. Rev. E 90 042146 [42] Holubec V 2014 J. Stat. Mech. P05022 [43] Martínez I A, Roldán É, Dinis L, Petrov D, Parrondo J M and Rica R A 2016 Nat. Phys. 12 67 [44] Restrepo S, Cerrillo J, Strasberg P and Schaller G 2018 New J. Phys. 20 053063 [45] Segal D 2008 Phys. Rev. Lett. 100 105901 [46] Esposito M, Lindenberg K and Van den Broeck C 2009 Europhys. Lett. 85 60010 [47] Sánchez R and Büttiker M 2011 Phys. Rev. B 83 085428 [48] Strasberg P, Schaller G, Brandes T and Esposito M 2013 Phys. Rev. Lett. 110 040601 [49] Sothmann B, Sánchez R and Jordan A N 2014 Nanotechnology 26 032001 [50] Benenti G, Casati G, Saito K and Whitney R S 2017 Phys. Rep. 694 1–124 [51] Feshchenko A V, Koski J V and Pekola J P 2014 Phys. Rev. B 90 201407(R) [52] Hartmann F, Pfeffer P, Höfling S, Kamp M and Worschech L 2015 Phys. Rev. Lett. 114 146805 [53] Thierschmann H, Sánchez R, Sothmann B, Arnold F, Heyn C, Hansen W, Buhmann H and Molenkamp L W 2015 Nat. Nanotechnol. 10 854–8 [54] Filliger R and Reimann P 2007 Phys. Rev. Lett. 99 230602 [55] Chiang K H, Lee C L, Lai P Y and Chen Y F 2017 Phys. Rev. E 96 032123 [56] Roulet A, Nimmrichter S, Arrazola J M, Seah S and Scarani V 2017 Phys. Rev. E 95 062131 [57] Fogedby H C and Imparato A 2018 Europhys. Lett. 122 10006 [58] Marchegiani G, Virtanen P, Giazotto F and Campisi M 2016 Phys. Rev. Appl. 6 054014 [59] Alicki R, Gelbwaser-Klimovsky D and Szczygielski K 2015 J. Phys. A 49 015002 [60] Alicki R, Gelbwaser-Klimovsky D and Jenkins A 2017 Ann. Phys. 378 71–87 [61] Alicki R 2016 Entropy 18 210 [62] Serra-Garcia M, Foehr A, Molerón M, Lydon J, Chong C and Daraio C 2016 Phys. Rev. Lett. 117 010602 [63] Bang J, Pan R, Hoang T M, Ahn J, Jarzynski C, Quan H T and Li T 2018 New J. Phys. 20 103032 [64] Koch J, von Oppen F, Oreg Y and Sela E 2004 Phys. Rev. B 70 195107 [65] Galperin M, Saito K, Balatsky A V and Nitzan A 2009 Phys. Rev. B 80 115427 [66] Romano G, Gagliardi A, Pecchia A and Di Carlo A 2010 Phys. Rev. B 81 115438 [67] Schulze G et al 2008 Phys. Rev. Lett. 100 136801 [68] Scorrano A and Carcaterra A 2013 Mech. Syst. Signal Pr. 39 489–514 [69] Isacsson A 2001 Phys. Rev. B 64 035326 [70] Giaever I and Zeller H 1968 Phys. Rev. Lett. 20 1504 [71] Kulik I and Shekhter R 1975 Zh. Eksp. Teor. Fiz. 68 623–40 [72] Averin D and Likharev K 1986 J. Low Temp. Phys. 62 345–73 [73] Lai W, Cao Y and Ma Z 2012 J. Phys. Condens. Matter 24 175301 [74] Lai W, Xing Y and Ma Z 2013 J. Phys. Condens. Matter 25 205304 [75] Fedorets D, Gorelik L Y, Shekhter R I and Jonson M 2004 Phys. Rev. Lett. 92 166801 [76] Bonet E, Deshmukh M M and Ralph D 2002 Phys. Rev. B 65 045317 [77] Bagrets D and Nazarov Y V 2003 Phys. Rev. B 67 085316 [78] Harbola U, Esposito M and Mukamel S 2006 Phys. Rev. B 74 235309 [79] Kevorkian J and Cole J 1996 Multiscale and Singular Perturbation Methods Applied Mathematical Sciences (New York: Springer) [80] Nayfeh A H 2008 Perturbation Methods (New York: Wiley) [81] Bender C M and Orszag S A 2013 Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory (New York: Springer) [82] Blanter Y M, Usmani O and Nazarov Y V 2004 Phys. Rev. Lett. 93 136802 [83] Risken H 1996 The Fokker–Planck Equation (Berlin: Springer) [84] Esposito M, Lindenberg K and Van den Broeck C 2009 Europhys. Lett. 85 60010 [85] Van den Broeck C et al 2013 Phys. Complex Colloids 184 155–93 [86] Esposito M 2012 Phys. Rev. E 85 041125 [87] Herpich T, Thingna J and Esposito M 2018 Phys. Rev. X 8 031056 [88] Tonekaboni B, Lovett B W and Stace T M 2018 arXiv:1809.04251 [89] Croy A and Eisfeld A 2012 Europhys. Lett. 98 68004 [90] Wiseman H M and Milburn G J 2010 Quantum Measurement and Control (Cambridge: Cambridge University Press) 26