Stochastic thermodynamics of self-oscillations: the electron shuttle

Self-oscillation is a phenomenon studied across many scientific disciplines, including the engineering of efficient heat engines and electric generators. We investigate the single electron shuttle, a model nano-scale system that exhibits a spontaneous transition towards self-oscillation, from a thermodynamic perspective. We analyze the model at three different levels of description: The fully stochastic level based on Fokker-Planck and Langevin equations, the mean-field 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 mean-field 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.


I. INTRODUCTION
Self-oscillation has been described as "the generation and maintenance of a periodic motion by a source of power that lacks a corre-sponding 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][3][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, selfoscillatory systems provide a useful transduction mechanism for the design of autonomous motors and heat engines. One particularly interesting system exhibiting selfoscillation 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][7][8][9][10][11] to semiclassical [11][12][13] and completely classical descriptions [5,[14][15][16]. The electron shuttle has been experimentally realized by a vibrational fullerene molecule [17], gold grains [18][19][20] as well as nanopillars [21] as molecular junctions between two leads. Also macroscopic electron shuttles, consisting of a pendulum between two capacitor plates, have been investigated [22]. Reviews on the electron shuttle can be found in Refs. [23][24][25][26][27][28].
Classical self-oscillating systems have been analyzed using the tools of deterministic non-linear dynamics [1,29]. 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 [11]. 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, treating the electron shuttle as a paradigmatic isothermal engine that converts direct electric current into periodic mechanical motion. 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 behavior of nanoscale systems. As described in review articles and monographs [30][31][32][33][34][35][36], 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 [37][38][39][40][41][42]. In these models, externally applied time-arXiv:1902.08174v1 [cond-mat.stat-mech] 21 Feb 2019 periodic driving leads to the conversion of thermal fluctuations into work. By contrast, stochastic self-oscillatory systems represent autonomous nano-scale engines, characterized by the absence of an externally imposed cycle. A variety of autonomous engines have recently drawn both theoretical and experimental attention. These include a Brownian gyrator [43,44], a rotor engine [45,46], a heat engine based on Josephson junctions [47] or mechanical resonators [44,48], and an experimental realization of the Feynman's ratchet-and-pawl mechanism [49]. While particular thermodynamic aspects such as nonequilibrium hot electron transport [50], subresonance inelastic electronic transport [51,52] and tip-induced cooling [53] 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 modeled 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 Refs. [5,14]. We find that the abrupt onset of self-oscillatory behavior 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 (Sec. II) followed by mathematical descriptions of the system at the different levels mentioned above (Sec. III): the fully stochastic model in Sec. III A, the mean-field approach in Sec. III B, and the intermediate, perturbative model in Sec. III C. The dynamics at the different levels are discussed and compared in Sec. III D. In Sec. IV the first and second laws of thermodynamics are derived at the different levels of description (Secs. IV A-IV C), followed by a discussion of the thermodynamic behavior of the electron shuttle (Sec. IV D). Finally, in Sec. V, we discuss our findings and point out future applications.

II. PHENOMENOLOGY
In this section we explain the basic mechanism of the electron shuttle (see also Fig. 1) before introducing the mathematical descriptions in Sec. III. The shuttle is composed of a metallic grain [18][19][20] or molecular cluster [21,54] and a nanomechanical oscillator (e.g. a cantilever FIG. 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 center. For weak friction the shuttle may pass the center 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). [55] or an oscillating molecule [17,25]), 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 Fig. 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 the oscillator restoring force the shuttle returns (see Fig. 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.

III. MODELLING
In this section we discuss different levels of description of the electron shuttle, i.e., a fully stochastic description in Sec. III A, a mean-field approximate description in Sec. III B and a perturbative description based on time scale separation in Sec. III C. We will then compare and discuss the dynamics of the system at the different levels in Sec. III D.
In the literature there exist proposals to describe the electron shuttle fully quantum mechanically [6][7][8][9][10][11], semiclassically [11][12][13] or fully classically [5,[14][15][16]. 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 [24]. The latter is the case for an experimental realization of the shuttle with a gold grain [18] and is sometimes referred to as classical shuttling of particles [24]. The underlying mechanism (tunneling of electrons via Fermi's golden rule), is nevertheless intrinsically quantum.

A. 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 quantum dot (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 [56][57][58]. We use this simplified model of a single electron shuttle for illustrational and numerical purposes, but the phenomenology discussed in Sec. II 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 ∈ R and velocity v ∈ R. 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 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 F el = αV q acts on the shuttle when it is charged (q = 1), pushing it towards the reservoir with lower chemical potential (see Fig. 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,17,24,25]. 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 [17]. The restoring force acting on the shuttle is then given by F harm = −kx with spring constant k, and the shuttle is damped by F damp = −γ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 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 = γ (β osc m 2 ). The first line of Eq. (1) describes the underdamped evolution of the oscillator in the potential U q (x) = 1 2 kx 2 − αV qx, at fixed electron charge q. The second line 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 of electrons between the QD and the fermionic lead ν ∈ {L, R}. The transition rates R ν qq ′ (x) are given by and R ν qq (x) = − ∑ q ′ ≠q R ν q ′ 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 center of mass of the shuttle [5,12,[59][60][61], 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 ν (ω) ≡ [exp(β ν (ω − µ ν )) + 1] −1 with inverse temperature β ν . Note that the quantity ε − αV x enters the Fermi functions in Eq. (2), as the energy of the shuttle depends on both the QD energy ε and the electrostatic potential −αV x (see also Sec. IV A).
Eq. (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, β osc = β ν = β. Nonequilibrium conditions then arise solely due to the applied bias voltage, i.e., µ L ≠ µ R .
Since the space of dynamical variables defined by the triple (x, v, q) is large, solving Eq. (1) numerically is ex-pensive. We therefore turn to the trajectory representation of the single electron shuttle. The coupled stochastic differential equations produce the FPE (1) at the ensemble level, as we show in Appendix A by explicitly looking at the evolution of averages. In Eq. (4) the thermal fluctuations are taken into account by a Wiener process dB(t) with zero mean E [dB(t)] = 0 and variance E (dB(t)) 2 = dt. Here, denotes an average of the stochastic process. Eqs. (3) and (4) represent for a fixed q the Langevin equation of an underdamped particle moving in the shifted harmonic potential U q (x). Eq. (5) describes changes in the charge state q due to the stochastic tunnelling of electrons. The independent Poisson increments dN ν q ′ q (x, t) ∈ {0, 1} obey the statistics: 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 R ν q ′ q (x). The second line in Eq. (6) enforces that only one tunnelling event per time interval can occur, i.e., either all dN ν q ′ q (x, t) are zero or dN ν q ′ 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 Eqs. (3) and (4) describe a simple underdamped harmonic oscillator. However, the tunnelling of electrons still depends on x [see Eq. (2)] and therefore the QD remains coupled to the oscillator. Second, a complete decoupling of the QD and the oscillator is achieved in the limit λ → ∞ and α → 0. In that case the QD coupled to the fermionic leads describes the well known single electron transistor (SET) [62][63][64].

B. Mean-field approximation
In order to understand the nonlinear dynamics of the compound system of QD and oscillator, we first look at the mean-field equations derived from the full stochastic evolution. From the FPE, Eq. (1), we obtain for the ensemble averaged position ⟨x⟩ and velocity ⟨v⟩: Here and throughout the paper, denotes an ensemble average, and are the probabilities for the QD to be empty and occupied, respectively. Eqs. (7) are exact, but in order for them to form a closed set we need an expression for dp 1 dt. Integrating Eq. (1) over x and v we obtain (11) Due to the nonlinearity of the tunnelling rates with respect to x [see Eq. (2)] we approximate and we refer to this as the mean-field (MF) approximation. Note that if the tunnelling rates R ν qq ′ (x) were linear in x, Eq. (12) would be an equality and Eq. (7) would be closed without the MF approximation.
The MF approximation is thus described by the nonlinear differential equations wherep ≡ (p 0 ,p 1 ) ⊺ andq ≡p 1 . The entries of the rate matrix R ν (x) are definded by Eq.
. 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 Eqs. (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 2π 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 Eq. (14). Eqs. (13)-(15) are identical to those found in the original proposal of Gorelik et al. [5] for the case of one excess electron.

C. Multiple scale perturbation theory
To improve on the MF approximation, which only captures the average dynamics of the electron shuttle, we can perturbatively solve the FPE, Eq. (1), by assuming a separation of time scales between the short dwell-time of electrons and the slow movement of the oscillator, and applying multiple scale (MS) perturbation theory [65][66][67]. Specifically, we assume that during one oscillation there are many electron tunnelling events: Γ ≫ k m. We provide details of the MS calculation in Appendix B, and summarize the result here.
Working to first order in the perturbation we obtain where the vector (π 0 (x), π 1 (x)) T denotes the stationary state of the rate matrix R(x) = ∑ ν R ν (x), i.e. it is the right eigenvector corresponding to the zero eigenvalue, normalized to unity: ∑ q π q (x) = 1. The probability density to find the oscillator at position x with velocity v at time t is given byp(x, v, t), which obeys the FPE: Here, q eq (x) ≡ π 1 (x) is the instantaneous stationary charge of the QD given by (18) As seen in Eq. (17), the effects of the electronic degrees of freedom on the evolution of the oscillator are incorporated into an effective potential U eff (x), along with position dependent friction and diffusion coefficients: where The zeroth order perturbation (see Appendix B) corresponds to an adiabatic approximation, i.e., infinite time scale separation, Γ → ∞. In that limit we haveγ(x) → γ andD(x) → D, and Eq. (17) describes underdamped Brownian motion in an effective potential U eff (x), at inverse temperature γ (Dm 2 ) = β 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 Eq. (17), the x-dependence ofγ (Dm 2 ) breaks detailed balance, giving rise to non-equilibrium behavior and allowing for the possibility of self-oscillations.
To solve Eq. (17) approximately, we parametrize x and v by the energy E and the oscillation phase θ, such that 1 2 kx 2 + 1 2 mv 2 = E, and we assume that the probability density does not depend on the phase [68]: p(E, θ, t) ≈p(E, t). With the transformations of Eq. (22) and the latter assumption we find a FPE for the energy distribution by averaging over the angle θ: wherep(E, t) is the transformed probability distributioñ p(x, v, t). The effective friction and diffusion parameters take after the transformation the form Solving for the steady state of Eq. (23), i.e. ∂p ∂t = 0, we get [69]p where N is a normalization constant.

D. 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 C.

Stochastic dynamics
We start by looking at the fully stochastic model, given by Eq. (1). Rather than solving the FPE directly, we generated stochastic trajectories evolving under the Langevin Eqs. Also, the tunnelling events in Fig. 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 Sec. II.
The directed shuttling of electrons coincides with selfoscillation, as illustrated in Fig. 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 Eqs. , the probability density is concentrated around a circular orbit, revealing selfoscillatory harmonic motion with some amplitude and phase noise. These behaviours are consistent with the trajectories shown in Fig. 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 Eq. (2)].

Mean-field dynamics
We now turn to the mean-field (MF) dynamics. The white dot and circle in Fig. 2 c) and d) correspond to the solutions of the MF model given by Eqs. (13)- (15). As we can see the MF solutions coincide very well with the stochastic phase-space distribution. As shown in previous extensive studies [11], when the parameter V crosses a critical valueV 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 βV cr = 15.0 (see Appendix D). Three dynamical regimes can be characterized, as we discuss below.
Single electron transistor (SET) regime (I): The point (x fix ,v fix ) = (αVq k, 0) is a fixed point of Eqs. (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 Fig. 2 c)]. In   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 Fig. 2 b) and Fig. 3 c)]. When the shuttle is occupied by an electron, the electrostatic force is sufficient to overcome friction, leading to 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 βV = 55.0 there are still more tunnelling events than from shuttling electrons one by one, which in Fig. 3 c) can be seen from the fact that both currentsĪ ν M are finite when x ≈ 0. We also see this in the stochastic case very clearly [see Fig. 2 Crossover regime (II): When V ≈V cr the system exhibits both SET and shuttle behaviour. Above the critical bias voltageV 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 positionx and velocityv -as in the shuttling regime -as well as chargeq and matter currentsĪ ν M [see Fig. 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.

Perturbative dynamics
To gain further insight into the transition to selfoscillation, we solve the full FPE, Eq. (1), perturbatively by imposing a time scale separation between the rapid tunnelling events of electrons and the slow movement of the oscillator (see Sec. III C). In Fig. 4 a) we plot the steady state probability densityp ss (E) obtained from this calculation [see Eq.

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) =Ā cos(ωt +φ 0 ) +x fix , and we choose A MF =Ā 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 Fig. 2 with fit parameters 1 σ 2 , c, and x 0 . We then define the self-oscillation amplitude A FPE in terms of the value(s) x at which the function f (x) = g(x + c) has a maximum: when x 0 ≤ σ, f (x) has a unique maximum at x = A FPE = 0, and when x 0 > σ, f (x) has distinct maxima at x = ±A FPE . Note that this definition is not sensitive to small oscillations of the stochastic system, as it gives A FPE = 0 when x 0 ≤ σ, even though the shuttle may be self-oscillating. Finally, for the MS perturbative solution, we define the self-oscillation amplitude as A MS = 2E max k, where E max is the value of E at which the functionp ss (E) is maximized [see Eqs. (22) and (25)]. Similarly to A FPE , and for the same reason, A MS is not sensitive to small oscillations.
In Fig. 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 onset of oscillation at a critical value of the voltage. The specific values are given by βV cr = 15.0, βV * FPE = 13.2 and βV * 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 A FPE and A MS 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 Fig. 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 Sec. III D 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, A FPE is quite noisy close to the onset; see inset of Fig. 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.

IV. 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 behavior of the electron shuttle at the different levels of description, focusing on the thermodynamic signatures of the onset of spontaneous self-oscillation.

A. Stochastic thermodynamics
We start by looking at the full stochastic model. The total energy of the coupled system is given by 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 [31] where ○ denotes Stratonovich-type calculus. In the last term, dq = ∑ ν dq ν and dq ν = ∑ q ′ (q ′ − q)dN ν q ′ q (x, t) denotes an electron jump with respect to reservoir ν. The second term involving the velocity can be re-expressed by multiplying Eq. (4) with v. This yields mv ○ dv = (−kxv − γv 2 + αV qv)dt + 2γ β osc v ○ dB(t) (28) Inserting Eq. (28) into Eq. (27), we get Here, we have introduced the chemical work δW chem = ∑ ν µ ν dq ν and the heat flow to the oscillator from its thermal reservoir due to friction and thermal noise δQ osc = −γv 2 dt + 2γ β osc v ○ dB(t) [31]. The remaining terms in Eq. (29) are identified as heat exchanged with the reservoir ν, defined as δQ ν = (ε − αV x − µ ν ) ○ dq ν . 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 Eq. (29) over many realizations, equivalently by averaging with respect to the probability density (see Sec. III A and Appendix A): We note that derivatives with respect to time ( d dt ) denote exact (or complete) differentials whereas a dot (⋅) denotes inexact ones. The average heat absorbed from reservoir ν is given by Here, ⟨I ν M ⟩ ≡ E [dq ν dt] is the matter current from reservoir ν, and represents the position-current correlation. The average chemical work is given by and the heat current entering from the reservoir of the oscillator is The latter equation is formally equivalent to the definition of heat flow for underdamped Langevin dynamics [31]. Similarly, the definition of the chemical work flow [see Eq. (34)] is consistent with the corresponding definition for the SET (see, e.g., Refs. [70] and [35] and references therein). However, the definition of heat with respect to left and right leads [see Eq. (31)] differs by the additional contribution of −αV ⟨xI ν M ⟩, 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, Eq. (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 where p(x, v, q, t) is the solution of the FPE, Eq. (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, lim x→±∞ xp = lim v→±∞ vp = 0) we obtain where is a probability current. LettingṠ 1 (t) andṠ 2 (t) denote the two terms on the right side of Eq. (37), we integrate by parts to rewrite the first term as follows: From Eq. (35) we obtain Summing Eqs. (39) and (40) we arrive aṫ Next, we rewrite the second term on the right side of Eq. (37) as follows: From the property of (local) detailed balance obeyed by the electron tunnelling rates [see Eq. (2)], i.e.
we derive the identity where the first term on the right relates to heat exchange with the fermionic leads [see Eqs. (31) - (33)]. Summing Eqs. (43) and (45) and rearranging terms, we obtaiṅ wherė (47) Here, non-negativity follows from the log-sum inequality.
Adding Eqs. (41) and (46), we find that the rate of change of the system's Shannon entropy is given by Here, the entropy flow rateṠ e is the rate at which the total entropy of the reservoirs decreases due to heat exchange with the system [71]. The quantitẏ is the total entropy production rate, which can be expressed as the sum of two independently non-negative contributions [Eq. 50)], from the continuous [Eq. (42)] and discrete [Eq. (47)] degrees of freedom. The nonnegativity ofΣ, Eq. (50), shows that the second law holds in our system. We note that the two separate parts of the total entropy production rate [see Eqs. (42) and (47)] are formally equivalent to the definitions derived for an independent underdamped harmonic oscillator and an independent SET [35,70]. However, as is apparent especially at steady states, where d dt S(t) = 0, the entropy production rate involves the heat flows [see Eq. (51)], for which the definitions for the independent systems differ from the definition for the electron shuttle [see Eqs. (31) and (35)].

B. Mean-field thermodynamics
At the mean field level, the total energy of the system (denoted by an overbar) corresponding to Eqs. (13)- (15) is given byĒ and its rate of change is whereĪ Note thatĪ ν M > 0 denotes the flow of matter from reservoir ν into the system, and vice versa forĪ ν M < 0. Using Eqs. (13)- (15), the first law of thermodynamics takes the form HereQ is the heat flow into the QD from reservoir ν, is the rate at which chemical work is performed on the system, andQ is the heat flow into the oscillator from its bath. Note thatQ osc is always negative, in contrast to the stochastic case [see Eq. (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 Eqs. (13)- (14), while the quantum dot is represented by a probability distributionp = (p 0 ,p 1 ) T evolving under a master equation, Eq. (15). We therefore define the entropy of the system to be the Shannon entropy of the QD probability distribution: As in Sec. IV A, the total entropy production rate is the sum of the rates of change of the entropies of the system and the reservoirs: We now analyze the three terms on the right side of this equation.
The rate of change of the system entropy is given by By Eq. (58), the second term on the right side of Eq. (60) is equal to γβ oscv2 . To analyze the third term we use Eqs. (44), (54) and (56) to write Combining results, we obtaiṅ (63) in agreement with the second law.

C. Perturbative thermodynamics based on multiple scales
Multiple scale perturbation theory gives the following result for the shuttle probability density (see Sec. III C): where π q (x) denotes the instantaneous equilibrium distribution of the QD andp(x, v, t) is the solution to Eq. (17). Eq. (64) represents an approximate solution of the FPE, Eq. (1). Therefore, to compute thermodynamic quantities such as heat flows and chemical work at this level of approximation, we use Eq. (64) to evaluate the relevant averages introduced in Sec. IV A. Note that an alternative attempt, which we will not follow here, could be to derive the laws of thermodynamics based on the effective FPE, Eq. (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 [72].
To evaluate the heat flow into the oscillator, Eq. (35), we first write where ⟨•⟩ MS denotes an average taken with the densityp(x, v, t). Using the transformation to energy E and oscillation phase θ given by Eq. (22), and assuminĝ p(E, θ, t) ≈p(E, t) (see Sec. III C), we get after averaging over θ. Here,p(E, t) is the solution to Eq. (23), which at steady state is given by Eq. (25). For the chemical work and the heat exchanges with fermionic reservoirs, we need the matter currents ⟨I ν M ⟩ and position-current correlations, ⟨xI ν M ⟩. Substituting Eq. (64) into Eq. (32) we get Transforming to (E, θ)-space and averaging over θ gives For the position-current correlations, we similarly get Combining results with Eqs. (31), (34) and (35) we obtain As in Sec. IV A, the first law is expressed by Eq. (30), but the heat flows and chemical work are now given by Eqs. (72)-(74). Using Eq. (64), the system entropy [Eq. (36)] becomes a sum of distinct contributions from the harmonic oscillator and the quantum dot: where S QD (x) = − ∑ q π q (x) ln π q (x). The decompositionṠ =Ṡ e +Σ [see Eq. (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 π q (x)p(x, v, t) in Eq. (42), we finḋ Transforming to (E, θ)-space and averaging over θ then results iṅ Applying a similar analysis for the discrete degrees of freedom [see Eq. (47)] we geṫ which after transforming to (E, θ)-space and averaging over θ becomeṡ where Σ disc,MS is non-negative since σ(E) ≥ 0 by the log-sum inequality, andp(E, t) ≥ 0 is a probability density. Combining results, we geṫ again in agreement with the second law.

D. 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 C.

Entropy production rate
We first look at the steady state entropy production rates as a function of the applied bias voltage V . Fig. 5 shows the entropy production rate for the stochastic case (orange solid), the MF case (red dash-dotted) and the MS description (blue dotted). Also shown is the single electron transistor entropy production rateΣ SET (green thin), i.e. Eq. (47) for the case where the position x is fixed at the origin [35,70] or, alternatively, in the completely decoupled limit α → 0 , λ → ∞ (see Sec. III A).
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 =  (51), and with our choice of setting all reservoir temperatures to be equal, β ν = β osc = β, we obtaiṅ The last equality follows from the preservation of electron number, ⟨I R M ⟩ + ⟨I L M ⟩ = 0. In the MF case, we arrive at the analogous result using Eqs. (55), (57) and (60).
At the MF level, we see in Fig. 5 that below βV cr = 15.0 the entropy production rateΣ is essentially the same as for the SET. This is understandable: below the bifurcation, the system evolves to a stable fixed point, with the quantum dot at rest near the origin (see Sec. III D 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Σ fromΣ SET for V >V cr .
Interestingly, we see in Fig. 5 that self-oscillation lowers the rate of entropy production, relative to the value it would have taken had the quantum dot remained at rest; that is,Σ <Σ SET for V >V cr . In effect, above the critical voltage, when faced with a "choice" between two modes of behavior -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, Eq. (82). For βV ≫ 1 the SET current approaches I SET M = Γ 2, as our choice of chemical potentials produces a SET steady state in which p 0 = p 1 = 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Ī M = ω 2π, where ω is the oscillation frequency. Our parameter choices give I SET M = 0.5 andĪ 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 Fig. 5. By the same token, if the parameters were chosen such thatĪ M > I SET M (roughly, if ω = k m > πΓ), then we would getΣ >Σ SET .
In the stochastic case the oscillator undergoes thermal motions in the steady state, the matter current through the quantum dot is lower than the corresponding SET current, and as a resultΣ <Σ SET . Note that the onset of oscillations is not as clearly marked as in the MF case, rather the entropy production rate transitions smoothly from one regime to the other. From Fig. 5 we see that for both small and large bias voltages, the MF and stochastic entropy production rates agree quite well: both approach the SET value when V ≪V cr , and when V ≫V cr the MF value only slightly underestimates the entropy production rate. However, around the bifurcation atV cr the stochastic entropy production rate deviates substantially from the MF prediction.
We finally note that the stochastic entropy production rate is well approximated by the multi-scale (MS) results, up to V ≈V cr . As argued in Sec. III D, for V >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 Fig. 5.

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⟩ = dĒ dt = ⟨dE dt⟩ MS = 0, as we have verified numerically. In Fig. 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 µ ν = ε ± V 2, as can be seen from Eq. (31): ⟨Q ν ⟩ = (ε − µ ν ) ⟨I ν M ⟩ − αV ⟨xI ν M ⟩. By the conservation of the matter, the first term on the right is the same for the left and the right reservoir and differences arise only from the correlation of x and I ν M . Since α = 0.06 ≪ 1, the difference is barely noticeable in Fig. 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 (Fig. 5), at the MF level the onset of oscillations at βV cr is clearly reflected in the heat flowsQ ν (Fig. 6) andQ osc (Fig. 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 ν osc ⟩ ≈ ⟨Q ν 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.
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 voltageV 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.

V. 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 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 single electron transis-tor 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 single electron transistor.
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. It can be used 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 ). 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 Ref. [73]. Similarly, one can then use the mechanical motion in order to pump electrons and generate a matter current against an externally applied electric field. The thermodynamic analysis of such a device is subject of our further research in this direction. With the present work, we hope to stimulate a discussion about the thermodynamic possibilities of such devices.

ACKNOWLEDGMENTS
The authors thank T. Brandes for initiating this project. C. W. acknowledges fruitful discussions with S. Restrepo  Here, E [•] denotes averages of the stochastic process. We now evaluate the sum in Eq. (A2): Since all powers of the Poisson increment are of order dt, we have to evaluate the sum exactly. Since dN ν q ′ q (x, t) k = dN ν q ′ q (x, t) for all k ∈ N we can rewrite the expectation value as follows (omitting any dependencies on x and v): where the last equality follows from a general identity of point/Poisson processes (see, e.g., Eq. (B.54) in [74]). Hence, up to order O(dt) we write Since Eq. (A1) is equivalent to Eq. (A4) 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, Eq. (1), and the stochstic differential Eqs. (3)-(5) describe the same process.

Appendix B: Multiple scale perturbation theory
In this section we derive Eq. (17) from the full FPE, Eq. (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.
Since the fast time scale describes the dynamics of the two state system, we treat the L-part as perturbation and introduce the bookkeeping parameter ≪ 1 such that we can rewrite Eq. (B8) as The idea of MS perturbation theory is now to introduce two time scales, a fast one (t) and a slow one τ = t such that the probability density is a function of both times scales, i.e.p(x, v,t, τ ) = p(x, v, t). The temporal derivative transforms to a sum provided by the chain rule: Then, Eq. (B9) is given by From this point on, we will refer tot as t and tõ p(x, v, t, τ ) as p(t, τ ). Assuming that we can express p as a series of orders of , 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, τ ) ≈ p (0) (t, τ ) + p (1) (t, τ ). (B13) We start with the governing equation for O( 0 ): The simplest solution of the ordinary differential Eq. (B14) is given by assuming that the left hand side of Eq. (B14) 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, τ ) = p (0) (τ ). This means p (0) (τ ) must be the eigenvector of R with eivenvalue 0, i.e. p (0) (τ ) = π(x)p (0) (τ ).
The specific form of p (0) (τ ) will be determined by the first order of the perturbation hierarchy. Note that, if we stop the perturbation theory here, Eq. (B15) implies an infinite time scale separation, which is equivalent to an adiabatic approximation.
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