Dynamics and current fluctuations in an ac-driven charge shuttle

The behaviour of a charge shuttle under a pure ac field has been recently considered theoretically and experimentally. If the system presents an asymmetry in the tunnelling amplitudes the device acts as a nanoelectromechanical rectifier, transforming a pure ac voltage field into a direct current. In this paper, we first review the model and the appearance of the rectifying effect for bias voltages below the threshold of self-oscillation. We discuss in some detail the dynamics of the central island that, like the current, presents a strong dependence on the forcing ac field frequency. In the presence of both a constant and a small oscillating bias voltage, we analyse the transition from the static to the self-oscillating solution. We then consider current fluctuations (full counting statistics) for periodic motion of the grain. We explicitly evaluate the current noise numerically and we find that it shows clear signatures of correlated transport at certain locking frequencies. In the adiabatic limit, we obtain a simple expression for the full counting statistics and calculate explicitly the first four moments.


Introduction
The field of nanoelectromechanical (NEMS) devices has been growing very rapidly in recent years both experimentally and theoretically. Important experimental results in this area include the use of a single electron transistor (SET) as displacement sensor [1] or the study of quantum transport through suspended nanotubes [2]- [4], oscillating molecules [5]- [8] and islands [9,10]. On the theoretical side, several works [11]- [23] have highlighted various aspects of the role of Coulomb blockade in NEMS.
An exciting prototype example of a mechanical assisted SET device is the charge shuttle proposed in 1998 by Gorelik et al [11]. As shown in [11] a SET with an oscillating central island can shuttle electrons between the two electrodes with a corresponding low noise level [12,20,21]. Although the experimental realization of the charge shuttle is difficult, promising systems such as C 60 molecules in break junctions [5,8] or silicon structures [9,10] already exist.
Several interesting features in the dynamics of the shuttle emerge when it is driven by an ac bias due to the interplay between the internal frequency of the oscillating island and the external frequency [24]. The most striking is probably a rectification effect when the contact resistances to the reservoirs are asymmetrical. Moreover, since the whole shuttling phenomenon is due to the nonlinear dependence of the electron transition rates on the position of the island, one expects that the mechanical motion can respond also at multiple harmonics both of the external and of the internal resonating frequency.
In the present paper, we expand our results on the ac-driven shuttle presented in [24] and discuss in some detail the behaviour of the average current and its fluctuations together with a detailed analysis of the dynamics of the island. The paper is divided as follows. In section 2, we present the model for the ac forced shuttle. In section 3, we describe the dynamical behaviour of the island and the current dependence on the external frequency. In section 4, the transition from a static to an oscillating solution is studied in the presence of a small ac field. In section 5, we derive an expression for the full counting statistics (FCS) in the adiabatic limit where we also discuss the first four moments of current fluctuations. Section 6 is devoted to the conclusions. force due to the dissipative medium, and an electric force due to the applied bias. The island is connected to the left and right leads through tunnel junctions characterized by resistances Here λ is the tunnelling length of the order of 1 nm, and x is the displacement from the equilibrium position in absence of external drive. The island is coupled to the leads and the gate through two capacitances C g and C. While it is crucial to take into account the (exponential) dependence of the tunnelling resistance on the position of the island, the dependence on x of the capacitances is much weaker and it can be ignored. The system is symmetrically biased at a voltage V(t) (V R = −V L = V/2), the charging energy is E C = e 2 /2C where e is the electron charge, C = C g + 2C, and the gate charge Q g is C g V g .
We will consider the case of low temperatures (k B T E C ), charge degeneracy (Q g = e/2), and voltages |V | < E C /e. In this regime, the grain can accommodate only n = 0 or 1 additional electrons (see right panel of figure 1). We also assume that typical driving frequencies ω are small compared to k B T/h E C /h. This condition on the frequency is rather weak (typically E C /h ∼ 10 THz) but it allows a simple description of the tunnelling in terms of time-dependent rates.
In the simplest approximation, the dynamics of the central island can be derived from the Newton equation [11] together with a stochastic equation governing the charge −en(t) on the island. In equation (1) m is the mass of the grain, ω o is the oscillator eigenfrequency, ηω o is a damping coefficient (η 1), and L is the distance between the two leads. In the regime of incoherent transport the (stochastic) evolution of the charge −en(t) is governed by the following four rates [25]: where equations (2) and (3) refer to n = 1 → 0 transitions while equations (4) and (5) to n = 0 → 1 transitions. Here FL, FR, TL and TR stands for from/to and left/right, indicating the direction for the electron tunnelling associated to the corresponding as shown in the central panel of figure 1.
is the Heaviside function. The current I is then determined by counting the net number of electrons that have passed through the shuttle during a given time t.

Current and mechanical response to the ac field
The problem can be studied numerically by simulating the stochastic process governed by the four rates equations (2)-(5) and by the deterministic evolution equation (1). We used a simple technique that consists of solving analytically Newton equations between two tunnelling events. One divides the trajectory into small intervals where the probability of occurrence of two tunnelling events is much smaller than one, and uses the rates to generate events with the appropriate statistics. In the following, we present results for typically 10 6 tunnelling events for each plotted point. After a transient time, the system reaches a stationary behaviour. In figure 2, the stationary dc current is shown as a function of the frequency of the external bias. The rich frequency landscape structure is generic. We observed a qualitatively identical behaviour in a wide range of parameters for which the dissipation is sufficiently large to prevent self-oscillation at constant voltage. The first conclusion one can draw is the existence of a rectification effect. The effect is of truly NEMS origin and most probably it has already been observed by Scheible and Blick [10]. They considered a nanopilar forced by an ac field, and they observed a rectification effects in the MHz range. The structure of the main resonance strongly resembles the prediction shown in figure 2, with a change of sign of the current at ω o . The theory presented strictly does not apply to the experiment of [10] which was performed at room temperature; nevertheless a qualitative comparison seems to indicate that the main features are present even when the Coulomb blockade is not completely established (E C ≈ 80 K in that experiment). A second important point to emphasize in figure 2 is the presence of a rich structure at low frequency. This is a consequence of the nonlinearity of the problem, since the dips and the peaks appear at rational values of the ratio ω/ω o (cf right panel of figure 2). In order to understand better the whole response, it is useful to exploit an analytical description of the shuttle within the approximation scheme proposed in [11] which we discuss below.

Mean-field approximation
If the electric force is much smaller than the mechanical one, = eV o /(ω 2 o mλL) 1, one can average the force generated by the stochastic variable n(t) over many periods of oscillation. One should in fact realize that what matters are the hopping events. If the charge on the island does not vary, the electric force is constant and thus no energy can be pumped on the system. It is thus the charge fluctuations that induce the changes in the dynamics. In order to have a sizable effect at 1, we need to take into account many hopping events. This allows us to neglect the fluctuations with respect to the average, even if the timescales of the hopping and of the oscillations are of the same order. Formally this consists in substituting into equation (1) its statistical mean: n(t) = n nP n (t), where P n (t) is the probability that the occupation number of the island is n (in our case only n = 0 and 1 are allowed values).
The charge dynamics in the central island is then determined by the the coupled set of equation which include the dynamical equation for the position of the island together with the master equation [25] for the probability vector |P = {P 0 , P 1 }.
The matrixˆ (t) is given bŷ The probability is conserved, i.e. P 0 + P 1 = 1. The rates depend on time through the extenal voltage V(t) and through the position of the island x(t). The instantaneous (average) current through the structure is 6 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT where the subscript a in the current indicates that the fluctuations of the force acting on the shuttle, due to the discrete nature of the charge tunnelling, are neglected. As shown in [11] the shuttle instability, at constant bias, is controlled by the ratio /η. It is thus possible to assume that both and η are small, with their ratio arbitrary.
We are now ready to discuss the dependence of the current on the external frequency. The most prominent structure, also observed in the experiments of [10], is present at ω ≈ ω o , and it corresponds to the main mechanical resonance. The current changes sign across the resonance. The shape of the resonance can be explained as an effect of the phase relation between the driving voltage and the displacement of the grain. A simple way to prove this fact is to solve equation (7) for a given function x(t) = A sin(ωt − φ) therefore neglecting higher harmonics. The crucial parameter is then the phase φ and its frequency dependence. We found that if we substitute the usual relation for a damped harmonic oscillator: we could reproduce qualitatively the behaviour of figure 2 near the resonance, and in particular the change of sign. We verified from the numerical simulations that φ(ω) around the main resonance agrees very well with the simple expression of the damped harmonic oscillator given above.
More insight can be gained by the solution of the coupled set of equations (6) and (7) with periodic boundary conditions in order to obtain the stationary behaviour. The results for three cases are shown in figure 3. The landscape near ω o of the rectified current can be understood by means of the following argument. The total resistance of the structure, has a minimum for positive x due to the asymmetry of the barriers. Then, even if the transport is not given by Ohm's law, the fact that R T (x) has a minimum shifted from the symmetric point allows larger currents to flow when the shuttle is near the right lead. For φ ≈ 0, the voltage is positive when the shuttle is on the right, while the opposite is true for φ ≈ π. The phase shift crosses sharply from 0 to π when ω is swept through the resonance. This explains while on the left of the resonance the rectified current is positive and on the right it is negative. The maximum is due to a combination of the fast frequency dependence of the phase and of amplitude that display the usual Lorenzian maximum at resonance. Between these two values, the current has to vanish for a value of the frequency that is very near ω o .
From the right panel of figure 2, it is clear that additional structures appear for ω ≈ ω o /ν with ν = 2, 3, . . . . Our numerical simulations indicate that, except for the fundamental frequency, even ν are favourite with respect to odd ones. As we already anticipated, the motion of the shuttle and the oscillating source become synchronized at commensurate frequencies whose ratio is µ/ν (with µ positive integer). This ratio, also known as the winding number, can be defined as where θ(t) is the accumulated angle of rotation of the representative point (ẋ,ẍ) in the phase space. In practice θ(t)/2π gives the number of oscillations performed by the shuttle during the time t. It is however important to choose as a representative point the velocityẋ and the accelerationẍ instead of the position x and the velocityẋ. In this way, the angle really gives the number of 'oscillations' performed by the shuttle. At large value of w the structure can be quite complex as can be seen from figure 4 and the above definition succeeds in giving the right number of closed cycles in each case. Again these plots have been obtained in the meanfield approximation. At the main resonances for ω/ω o ≈ 1/ν, the shuttle performs ν nearly harmonic oscillations per period of the forcing field. Since different w implies different topology in the phase space, it is interesting to follow the transition from one winding number to another.

I(t)
Left ω/ω o = 1.02 In particular the transition from w = 2 to w = 4 is shown in figure 4. This is performed without passing through the w = 3 state, that appears actually as a small plateau in the middle of the w = 2 plateau.
When the system is frequency locked at a winding number w, it is possible to define the phase shift φ(t) = θ(t) − w ωt. For an accurate analysis, the value used for w in the definition of φ(t) must be the best µ/ν approximation to w, since even corrections of the order of 10 −3 to the exact µ/ν results would accumulate in the calculation of φ, and give an arbitrary shift of the order of π. After a transient time for perfect locking, if the oscillation is perfectly harmonic, φ should no longer depend on t. Thus an additional important quantity to analyse is the phase shift variance ( φ) 2 = φ 2 − φ 2 . This is calculated by sampling 20 points per cycle over 10 3 cycles. For a given w, the smaller is the value of φ, the better the system locks to that external frequency. A zero value of the φ indicates that the oscillation is also harmonic. A periodic, but non-harmonic, oscillation will induce a fluctuation of the phase in time, but the fluctuations are correlated and φ will be smaller than the completely uncorrelated result. The numerical results for w and φ are shown in figure 5. It shows the dependence of the winding number as  a function of the external frequency (left panel) together with the analysis of φ (right panel) calculated from the stochastic simulation, green points, and from the average approximation, black line. The locking at rational winding numbers is confirmed by the presence of plateaux of decreasing width. As expected the most stable plateau is at w = 1: the system locks very well at this frequency. We find that this holds up to very high frequency in the average approximation, where w = 1 seems the only possible winding number. The stochastic simulation would indicate instead that for ω/ω o > 1.3 the locking with w = 1, as defined from equation (10), is no more established. We verified that this upper limit to the w = 1 plateaux depends essentially linearly on the asymmetry R R /R L . Quite generally plateau size increases by increasing the asymmetry. However, by studying φ for w = 1 in the whole frequency range, one actually obtains that a correlation is always present (i.e. φ < π/ √ 3 ≈ 1.81) even when equation (10) gives a value of w different from one. The stochastic fluctuations thus unlock the shuttle globally, but not locally. Looking at the presence of local locking at other winding numbers we find that for instance w = 2 is clearly present for ω > 1 and reversely w = 1 is present around ω = 1/2 (see red dashed lines in the right panel of figure 5). For global locking, only one phase variance is minimal. It corresponds to the 'dominant' winding number. The presence of correlations of other winding numbers may indicate partial locking at these winding numbers (as in the region 0.6 ω 0.9 for w = 1 and 2) or the contribution of higher harmonics of x(t) (as in the region ω > 1).

Instability in the presence of the ac-drive
In the previous section, we saw the response of a shuttle to a pure ac field for subcritical values of /η, i.e. when the shuttle instability is not present for static voltage. In this section, we describe the transition from the stable to the unstable state in presence of a small forcing ac field. When the instability develops continuously from the static solution (soft excitation), one can study analytically the effect of an ac field by expanding around the static solution. We begin by solving the evolution equation in Fourier series The electric force for a static voltage can also be written as a Fourier series as follows: where C ν = (ω/π) 2π/ω 0 sin(νωt) n P n (t)n and D ν has an identical expression with the cosine substituting the sine. Plugging these expressions into equation (1), one readily obtains the linear Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT relation between the strength of the forcing harmonics (C and D) and the response of the system (A and B) for ν > 1 The presence of Lorenzian resonances at ω ν = ω o /ν is clearly visible in this form. As soon as a forcing exists at a certain ν, the response will be stronger near the resonance frequency ω o /ν. This is the reason why weak nonlinearity induces structures in the rectified current at frequencies ω o /ν, as shown in figure 1. For a static voltage there is no reference phase, and thus we can set, for instance, B 1 = 0. The coefficients C 1 and D 1 depend through the master equation (7) strongly on A 1 and more weakly on the other coefficients for ν > 1. Equation (13) may then be solved by non-vanishing values of A 1 only for a precise value of the frequency ω = ω R . This is very near the bare oscillation value ω o : while the self-consistent equation for A 1 becomes (near the resonance) The physical interpretation of this last equation is very simple, it represents the energy balance between the dissipated energy (left-hand side) and the pumped energy (right-hand side) during one oscillation [11]. For small value of a = A 1 /λ, it is possible to calculate analytically C 1 and D 1 , it suffices to solve the master equation by expanding in the amplitude of oscillation a. In our case, defining L (0) = R (0) =˜ ω o , at order a 3 we obtain and The energy balance equation (15) gives the condition > c = η(1 + 4˜ 2 )/˜ for the value a = 0 to become locally unstable. The a 3 term in equation (17) determines the nature of the instability. If it is negative then a increases continuously from zero when crosses the value c . Note also that once the problem has been solved for the principal harmonics, the amplitude of the higher harmonics can be obtained by using equation (13) the solution of the first harmonics to calculate C ν and D ν for ν > 1.
In the presence of a small ac additional forcing field the equations (12) and (13) remain valid with the replacement → [1 + α sin(ωt + δ)] and P n (t) → P n (t)[1 + α sin(ωt + δ)], with α = V ac /V dc 1. We can again set B 1 = 0, but we need to take into account the phase shift δ of the forcing field with respect to the oscillation. By keeping only the principal harmonics we find that the equations for C 1 and D 1 are modified by the appearance of the terms α cos δ/2 and α sin δ/2 on the right-hand side of equations (16) and (17), respectively. The two unknown quantities a and δ are determined by equations (14) and (15). For c , one finds the usual result for a damped harmonic oscillator: a has a Lorenzian shape and the maximum value is ≈ α/[2η( − c )] for δ = π/2. According to this expression the amplitude diverges for approaching c , but its range of validity stops when the cubic term into equation (17) becomes important. At that point the sin δ term changes sign and the topological structure of the resonance changes drastically. In figure 6, the evolution of ω(a) as a function of is shown as obtained from a numerical solution of the system of equations (14) and (15) where D 1 and C 1 include the linear term in α. Figure 6 shows the transition from a standard harmonic oscillator response (top-left panel), to a self-oscillating system where periodic solutions occur only for a small range of the amplitude a and the frequency ω (bottom-right panel). The evolution is characterized by a change of the topology of the solution in the a-ω plane. Increasing the static electric field ( ), the nonlinearity becomes more important and at some point (slightly larger than c ) the Lorenzian shape is so deformed that an independent island develops (evident in the last two plots). In order to fully understand this evolution, we need to study also the stability of the solution. This can be done by keeping a slow time dependence of the phase δ and the amplitude a in deriving the equation of motion. Keeping only the first derivatives of a and δ one obtains Written in this way it is clear that the solution of equations (14) and (15) gives stationary states (the right vector is zero) for the time evolution of the phase and the amplitude. Studying the real Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT part of the eigenvalues of the linear expansion of equation (18) around these solutions allows the determination of their stability. Stable solutions are indicated in red in figure 6. Therefore by increasing the first part of the solutions to become unstable are those for ω far from the resonance. We cannot predict the behaviour of the system in this case, we only can say that no periodic solution at the forcing frequency are possible. The stable regions restricts more and more around the resonance, and the shape also of the resonance changes drastically, as discussed above. In the end, the only stable part of the diagram is a short segment (whose size is controlled by α) around the solution ω R , a R that would be present if no ac field was applied.

Full counting statistics
As we have seen the asymmetry of the resistances allows us to unveil the complex behaviour of the shuttle by studying the current dependence on the forcing frequency. The correlated charge transfer in shuttles appears to be much more visible in the correlations of the current fluctuations.
For the shuttle, at a constant bias voltage, the noise has been considered for the first time in [12] neglecting the correlation among different cycles. Within a similar model, but taking into account the correlation between different cycles, the noise and the FCS of charge transfer was then calculated by one of the author [20]. The noise taking into account the quantum nature of the oscillator was obtained by Novotný et al [21] by resorting to a numerical solution of generalized Bloch equations. An analytical approach was proposed at the same time in [26]. The method to calculate the FCS in NEMS systems proposed in [20] was extended to the quantum description of the oscillation in [27]. The FCS for a superconducting shuttle was also considered [28]. Knowledge of the FCS gives a deeper insight on the dynamical evolution of the nanomechanical systems. For instance in [20] the distribution of transmitted particles appeared to be strongly asymmetrical, due to the fact that different parameters control the transfer of one or two electrons per cycle.
In the case of the ac-driven shuttle, an example of the behaviour of the noise as a function of the external frequency is given in figure 7. Additional interesting behaviour is expected in higher moments. The same method discussed in [20] can be applied to study the statistics of current fluctuations also in presence of an ac field. Current fluctuations can be described by the probability P t (n) of transferring n elementary charges during a measurement time t. After the seminal work of Levitov et al [29] the theory to obtain the FCS has been developed and applied to many different systems (for a review see for instance [30]). It is easier to calculate the generating function S that is the logarithm of the Fourier transform of P t : e −S t (χ) = n e inχ P t (n, χ).
Using the definition (19) it is easy to show that the the nth derivative with respect to iχ of −S gives the nth cumulant of transmitted particles. Explicitly, the first two are the average particles n and the its variance (n −n) 2 .
For the static SET, the problem has been solved by Bagrets and Nazarov [31]. They developed a technique to calculate in a very simple and efficient way the FCS for an arbitrary network of contacts in the Coulomb blockade regime (the so-called 'orthodox' regime, which is the one we are considering in the present paper). In [20] this technique was extended to take into account the motion of the grain, which reflects into a time dependence of the rates. It is possible to show that the generating function can be written as follows: here |P(0) is the vector introduced in equation (7) to describe the state of the system, and q| = {1, 1, . . . , 1}. The operatorˆ χ (t ) is the rate matrix defined by equation (8) where some elements have been modified. Explicitly, if we want to count particles crossing a certain tunnel junction, we need to attach a e iχ factor to the off-diagonal elements of (8) that induces a transfer of charges from the leads to the grain, and a factor e −iχ for the elements involving transfer of particles in the opposite direction. In our case if we count particles at the tunnelling contact with the left reservoir we obtain The problem is now reduced to solving the master equation (7) in the presence of the counting field and coupled to the Newton equation for x(t). When the solution is periodic the expression for the generating function can be further simplified. In this case, it is sufficient to integrate the master equation during one period 14 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT The generating function for the transfer of N electrons is determined by the Nth power of A. For large N the eigenvalue λ M of A with largest absolute value will dominate, and the initial conditions can be neglected. In this way, the generating function becomes simply The problem is thus reduced in finding this eigenvalue.
In the following, we discuss the adiabatic limit of ω ω o , ηω o and L (0) + R (0). At leading order in ω equation (20) becomes where λ 1 (χ, t) is the (right) eigenvalue of χ (t) at time t with minimum real part. The conservation of the number of particles implies that the real part is non-negative. Using equation (24) we calculate as an explicit example the first four cumulants of the current fluctuation for the ac-forced shuttle. For this purpose, we need the eigenvalue with minimal real part as a function of time forˆ defined by equation (8). The variable x(t) is given by the solution of Newton equation neglecting the time derivatives: where n(t) is given by the stationary solution of equation (7) n(t) = P 1 (t) = FL (t) + FR (t) We consider the case 1, then x/λ = v(t)n(t) 1, where v(t) = V(t)/V o = sin(ωt). The exponentials inˆ can be expanded. It is convenient to separate the evolution into two parts, the range 0 ωt < π with positive v, and the range π ωt < 2π, with negative v. The position satisfies at order : The eigenvalue has the form: λ = λ 0 + λ 1 . λ 0 is the eigenvalue with minimal real part ofˆ 0 . λ 1 can be obtained by first-order perturbation theory: λ 1 = f |ˆ 1 |e , where |e and |f are the right and left eigenvectors of 0 corresponding to λ 0 . The calculation is straightforward; we give the result for the first four cumulants C 1 is the average number of electrons transmitted per cycle. It is clear from the expression that C 1 and all odd cumulants vanish for symmetric structure β = 0. This can be verified in general since we find that λ(χ) for the first half cycle is equal to λ(−χ) for the second half.It can be interesting also to note that odd cumulants start at order . This shows that odd cumulants probe a truly nanomechanical effect, since parameterizes the response of the island to the electric field. The explicit form of the first moment, the average current, is The corresponding value is shown in figure 2 as a dashed line. The (small) difference with the numerical result comes from the expansion, ( = 0.5 in the simulation) as is clear with a comparison with the black triangle obtained by solving numerically the adiabatic equations for the current.

Conclusions
In this work, we analysed the properties of the charge shuttle under time-dependent driving showing that the response of the shuttle is rather rich due to the nonlinear dynamics of the grain. Rectification, synchronization, response to multiple of the external frequency appear in the electrical properties of the shuttle. All the results obtained here can be tested experimentally. Indeed, a recent paper by Scheible and Blick [10] already reports on some of the properties which we saw in the average current. Part of these results was presented elsewhere [20,24]. Here, we provided a new description on the dynamical behaviour of the system, the transition from the stationary to the self-oscillating solution, and the FCS in the presence of an ac-driving.