Dynamical instabilities of a resonator driven by a superconducting single-electron transistor

We investigate the dynamical instabilities of a resonator coupled to a superconducting single-electron transistor (SSET) tuned to the Josephson quasiparticle (JQP) resonance. Starting from the quantum master equation of the system, we use a standard semiclassical approximation to derive a closed set of mean field equations which describe the average dynamics of the resonator and SSET charge. Using amplitude and phase coordinates for the resonator and assuming that the amplitude changes much more slowly than the phase, we explore the instabilities which arise in the resonator dynamics as a function of coupling to the SSET, detuning from the JQP resonance and the resonator frequency. We find that the locations (in parameter space) and sizes of the limit cycle states predicted by the mean field equations agree well with numerical solutions of the full master equation for sufficiently weak SSET–resonator coupling. The mean field equations also give a good qualitative description of the set of dynamical transitions in the resonator state that occur as the coupling is progressively increased.


Introduction
Nanoelectromechanical systems in which a high frequency mechanical resonator is coupled to a mesoscopic conductor [1] have been predicted to display a wide variety of different dynamical behaviours depending on the nature of the conductor. When a mechanical resonator is linearly coupled to the transport electrons in either a quantum point contact [2]- [4], or a normal state single-electron transistor [5]- [7], the electrons act on the resonator like an effective thermal bath [2]- [6], [8]. Under such circumstances, the resonator is damped and reaches an almost Gaussian steady state whose width is set by the fluctuations in the motion of the charges through the conductor. In contrast, where the electro-mechanical coupling is nonlinear [9,10], or the conductor is close to a transport resonance [11]- [13], the mechanical resonator can be driven by the electrons into states of self-sustaining oscillation.
In this paper, we analyse the instabilities that arise in the dynamics of a mechanical resonator coupled to a superconducting single-electron transistor (SSET) [11], [14]- [16] operated in the vicinity of a transport resonance. In the SSET, transport resonances occur when states of the SSET island differing by one Cooper pair are degenerate so that coherent Cooper pair tunnelling between the island and one of the leads is possible. The simplest such resonance, which we concentrate on here, is called the Josephson quasiparticle (JQP) resonance [11,14,17,18] and involves a cycle of processes in which current flows via a combination of coherent Josephson tunnelling between the SSET island and one of the leads, followed by incoherent quasiparticle decays into the other junction.
When a SSET tuned to the vicinity of the JQP resonance is coupled to a mechanical resonator, the resonator dynamics is very sensitive to the precise choice of bias point for the SSET [11,14]. For operating points on one side of the JQP resonance, the SSET damps the resonator and can be regarded as an effective thermal bath, behaviour which was confirmed in recent experiments [16].

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
In contrast, for operating points on the other side of the resonance the electrical degrees of freedom pump the resonator leading to the possibility of states of self-sustaining oscillation. In this regime the resonator dynamics can be investigated either by numerical solution of the quantum master equation [12], or provided the resonator is much slower than the electrical degrees of freedom (as was the case in recent experiments [7,16]), an effective Fokker-Planck equation for the resonator can be derived [11,13].
Numerical solution of the SSET-resonator master equation [12] revealed interesting similarities with a quantum optical device known as the micromaser. In a micromaser [19,20], a cavity resonator is driven by the passage of a steady stream of excited two-level atoms. As the atom-cavity interaction is increased, the resonator undergoes a sequence of dynamical transitions leading in some cases to non-classical steady states. A corresponding set of dynamical transitions was found to occur in the SSET-resonator system together with regions of non-classicality.
Here we use an alternative approach, namely a mean field description to analyse the instabilities in the dynamics of the SSET-resonator system. This kind of approach has been used extensively in the analysis of nonlinear quantum optical systems where its usefulness has been well established [21,22]. The mean field equations we derive provide a relatively compact description of the system and their stability properties are readily analysed using the techniques of classical dynamical systems theory. Although the mean field description neglects some of the correlations in the system, comparison with numerical results from the full master equation show that for sufficiently weak coupling the mean field theory is close to being quantitatively correct. Although quantitative agreement is poor at stronger couplings, the mean field equations still give a good qualitative description of the system's dynamics displaying a similar sequence of transitions to that found numerically [12].
Although we will use language appropriate to a nanomechanical resonator throughout this paper, we anticipate that much of our analysis will also be relevant to the case of a superconducting resonator coupled to a SSET. Photon assisted satellites of JQP peaks have been observed in experiments in which a SSET was coupled to microstrip transmission line [23]. Furthermore, recent experiments have demonstrated coherent coupling between a superconducting stripline resonator and a Cooper-pair box with a coupling Hamiltonian between the resonator and the charges on the box very similar to the mechanical case [24]. Such systems would only need to be modified to allow quasiparticle tunnelling off the Cooper-pair box [25] to become analogous to the device considered here and hence it seems likely that a range of dynamical instabilities similar to that which we find for the mechanical case could also be produced in a superconducting resonator.
The outline of this paper is as follows. In section 2, we introduce our model of the SSETresonator system and give the appropriate quantum mechanical master equation. The mean field equations are derived in section 3. In section 4, we transform the mean field equations into plane polar coordinates and exploit the fact that the amplitude of the resonator oscillations is slowly changing to derive an effective amplitude dependent damping of the resonator due to the SSET. We then show that this quantity can be used to predict the presence of limit cycles in the resonator dynamics. We compare the limit cycle solutions predicted by the mean field theory with numerical calculations using the full master equation in section 5. Then in section 6, we briefly consider the implications of our theoretical calculations for experiments on nanomechanical-SSET systems. Finally, in section 7, we draw our conclusions. In the appendices, we give additional details about the derivation of the master equation for the SSET-resonator system and the stability analysis of the mean field equations.

Master equation
The SSET-resonator system we consider is shown schematically in figure 1(a). The SSET consists of an island linked to leads by tunnel junctions with resistances R J and capacitances C J which are taken to be equal for simplicity. The mechanical resonator is treated as a single-mode harmonic oscillator, with frequency ω, which forms a movable gate capacitor. The equilibrium position of the resonator is a distance d from the SSET island and we assume that the displacement of the resonator with respect to this equilibrium position, x, is always small in comparison (i.e. |x| d). Under these circumstances the gate capacitance can be expanded up to just linear order and we have C g (x) = C g (1 − x/d), where C g is the capacitance at x = 0.
The JQP resonance [11,14,17,18] which we are interested in here occurs when two conditions are met. Firstly, at the centre of the resonance there is no change in the energy of the system when a Cooper pair tunnels between one of the leads and the SSET island. Secondly, the bias voltage, V ds , must be large enough to allow quasiparticle decay processes between the island and the other lead to occur. The charge processes involved in the JQP resonance we consider are summarized in figure 1(b). Josephson tunnelling across the left junction leads to coherent oscillations between the SSET island states |0 and |2 which differ by a single Cooper pair. These oscillations are interrupted by the tunnelling of a quasiparticle from the island into the right lead which takes the island into charge state |1 . A further quasiparticle then tunnels from the island to the right lead, returning the island to state |0 and the cycle begins again. The large electrostatic charging energy of the small SSET island and the carefully controlled bias voltage ensure that at low temperatures all other charge processes are strongly suppressed.
The master equation for the SSET island charge and resonator is obtained from the full Hamiltonian of the system by tracing out the quasiparticle degrees of freedom [17]. The steps involved in this derivation (together with details of the approximations and simplifications involved) are sketched in appendix A. The final result is a master equation of the form [12]

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
The evolution consists of a coherent part, described by the Hamiltonian, H co , together with two dissipative terms L leads and L damping which describe quasiparticle decay from the island and the surroundings of the resonator respectively. The effective Hamiltonian is given by [12] where E is the energy difference between states |2 and |0 , E J is the Josephson energy of the superconductor, x and p are the canonical position and momentum operators of the resonator and m is the effective resonator mass. The resonator-SSET coupling is described by the length-scale x s which measures the shift in the equilibrium position of the resonator brought about by adding a single electronic charge to the SSET island [14]. The tunnelling of quasiparticles from the island is described by where is the quasiparticle decay rate. Note that this is a simplified expression which uses a single rate for the two decay processes and neglects both the position dependence of the quasiparticle rates and their dependence on bias point (these approximations are discussed in more detail in appendix A). Simplified in this way, the master equation is also essentially equivalent to that used to describe a resonator coupled to a double quantum dot [26]. The term which describes the effect of the resonator's surroundings is where γ ext represents the external damping and n the equilibrium average resonator occupation that would occur in the absence of the SSET. We use this form of the oscillator damping kernel as opposed to the quantum optical form which we have used elsewhere [12] as although it is less convenient for numerical calculations, it leads to the correct (translationally invariant) classical limit, which is essential if we wish to derive appropriate mean field equations [27,28]. The behaviour of the resonator depends very sensitively on the sign of E. For E > 0, Cooper pairs tend to absorb energy from the resonator, damping its motion [11,14]. 2 In contrast, for E < 0 and when the resonator is slow (i.e. ω ) the resonator tends to absorb energy from the Cooper pairs and it is in this regime that instabilities can occur. For faster resonator speeds [12] (ω/ 1), absorption of energy by the resonator from the SSET, and hence the location of instabilities, becomes concentrated around points where E −jhω with j an integer.

Mean field equations
The mean field equations for the SSET-resonator system consist of the set of equations of motion for the expectation values of all the relevant SSET and resonator operators. They are derived by 6 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT multiplying the master equation by each operator in turn and taking the tracė where the SSET charge operators are defined as p ij = |i j| with i, j = 0, 1, 2 and the expectation values are defined by · · · = Tr[ρ . . . ]. Note that, as discussed in appendix A, the p 12 , p 10 components of the density operator decouple from the rest and can safely be neglected from the mean field equations (and also the master equation itself) as they do not affect the resonator. In contrast to the simpler case of a resonator coupled to a normal state SET [6,29], the equations of motion of the first moments do not form a closed set. In other words, the dynamics of the first moments depend on the behaviour of higher moments leading to an infinite hierarchy of equations of motion for progressively higher order moments of the system operators. In order to derive a simple set of dynamical equations we need to make a semiclassical approximation whereby the expectation value of a product of two operators is replaced by the corresponding product of the expectation values of the individual operators [21,22]. Thus in this case we substitute x p 02 for xp 02 . This approximation cannot be justified rigorously. Indeed, dropping the correlations between x and ρ 02 means that we lose the ability to describe the noise in the system (which is determined by the behaviour of the higher moments). Nevertheless, the semiclassical approximation is well-known as a way of deriving a set of dynamical equations which typically capture many of the important elements of the dynamics of the corresponding quantum system [22]. In this case, we find that the much simpler mean field equations which result from the semiclassical approximation provide a very useful qualitative and, for sufficiently weak SSETresonator coupling, quantitative description of the different dynamical transitions which the resonator can undergo.
Thus, having made the semiclassical approximation, and using the conservation of probability ( p 00 + p 11 + p 22 = 1) to eliminate one of the charge variables, we obtain the following closed set of mean field equationṡ where x 2 q =h/(2mω) and we have dropped the angled brackets for convenience. The quantities α and β are defined as the real and imaginary parts of p 02 respectively.
Despite the decoupling of the second moment in the above equations, the mean field equations reproduce many of the features found in the dynamics of the full master equation. The time evolution of these equations reveals the resonator relaxing to a fixed point for some values of the parameters. For other values, the oscillations grow at first, before settling into a limit cycle, i.e. an oscillation at fixed amplitude. The mean field equations also show clear evidence of bistability for certain parameter values: in these cases the long-time behaviour of the resonator depends on the choice of initial conditions.
The fixed point solution of the mean field equations is obtained by setting the time derivatives to zero in equations (11)- (16). As we discuss in appendix B, the stability of the fixed point can be established using standard techniques [22].

Analysis in radial coordinates
Although a straightforward stability analysis of the mean field equations can tell us quite a lot about the dynamics of the system, we find that transforming the mean field equations into planepolar coordinates and making one further simplifying assumption allows us to proceed much further. The assumption we make is that the SSET-resonator coupling and external damping are sufficiently weak that the resonator's energy changes much more slowly than either its phase or the charge state of the SSET, conditions which are readily met for most practical implementations of the SSET-resonator system. This type of approximation has already proved useful in describing a variety of nanoelectromechanical and optomechanical systems [10,14,30]. In terms of our analysis of the mean field equations, the assumption of a wide separation of time scales between the evolution of the amplitude and phase of the resonator allows us to derive an effective equation of motion for the resonator amplitude from which we can determine the number and location of stable limit cycle solutions [30].
We proceed by rewriting the mean field equations (11)- (16) in plane polar co-ordinates (A, φ) which describe the amplitude and phase of the resonator, defined through the relations x − x fp = A cos φ and v = ωA sin φ, where x fp is the fixed point resonator position. In terms of 8 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT (A, φ) the mean field equations take the form, Note that x fp is equal to −x s (p 11 + 2p 22 ) when p 11 and p 22 have their fixed point values.
We now assume that the evolution of the resonator amplitude, A, is sufficiently slow that during a single period of oscillation it can be treated as a constant and that the phase evolution over the same time can be approximated by φ = ωt. The dynamics of the electronic degrees of freedom can then be obtained for this constant amplitude and steadily evolving phase. The resulting forced dynamics of the SSET charge variables are then averaged over the resonator period to calculate their effect on the resonator amplitude A which we characterize by an amplitude dependent damping term [10,11,30], γ SSET (A). This slow-A approximation will certainly be appropriate in the vicinity of a limit cycle solution of the mean field equations and will be valid more generally so long as the free (uncoupled) evolution of both the resonator and charge degrees of freedom is much faster than the rate of change of the resonator amplitude in the coupled system, i.e. for ω, γ ext and sufficiently weak SSET-resonator coupling. Whilst there is no simple way of evaluating the strength of the SSET resonator coupling at which this approximation breaks down, from equations (17) and (18), we expect that this approach will be valid for a given amplitude, A, if x s A (note that if the populations p 11 and p 22 remain much less than unity then the conditions on x s will be less restrictive). Furthermore, it is clear that in order to be consistent with these assumptions we should obtain an effective damping γ SSET (A) whose magnitude is much less than ω and .

Solving the electronic dynamics
With a fixed amplitude, A, and a harmonically oscillating phase, φ = ωt, the electronic degrees of freedom form a set of four coupled differential equations with time dependent coefficients. The dynamics of even this simplified system is still non-trivial, but we can make progress by making use of our assumption that the electronic degrees of freedom relax rapidly compared to 9 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT the resonator amplitude. We assume that transients in the charge dynamics can be neglected and hence that the effect on the amplitude of the resonator is dominated by the periodic response of the charge degrees of freedom to the harmonic drive.
To extract the relevant periodic solutions for the charges we rewrite the equations of motion in terms of a Fourier series consisting of harmonics of the resonator frequency, defined by, e.g.
The resulting equations for the Fourier coefficients of the electronic variables are, By solving for p n 11 , p n 22 , α n in terms of β n , we can rewrite these equations so that we have an equation for β n in terms of β n+1 , β n−1 , β n+2 , β n−2 . This is equivalent to a matrix equation involving a band-diagonal matrix with five nonzero diagonals. The matrix equation we must solve is where β is the vector of the coefficients β n (with n running from −∞ to +∞), M represents the matrix of terms derived from equations (23)- (26), and d is a vector representing the Kronecker delta, δ n,0 , i.e. d(0) = 1 and d(n = 0) = 0. Solving the Fourier series is equivalent to inverting the matrix M and the coefficients β n are given by, The components decay rapidly towards zero for |n| 1 and so the matrix can be truncated and solved numerically with little error. In practice we found that calculating the coefficients up to n = ± 80 proved more than adequate for the parameters we considered.
Examples of the oscillations in the SSET charge driven by the resonator are shown in figure 2 for the case where the resonator frequency matches the quasiparticle decay rate. We use a single simplified set of SSET parameters in all numerical calculations to illustrate the behaviour of our model [12,14]: = V ds /eR J , R J = h/e 2 and E J = hV ds /(16eR J ). The SSET-resonator coupling strength is parameterized by the dimensionless quantity κ = mω 2 x 2 s /eV ds . The charge oscillations in figure 2 resemble the behaviour of a periodically kicked damped oscillator, but more specifically, at larger amplitudes they are very similar to an oscillator whose frequency is time-dependent. Analogous oscillations have been seen in the amplitude of radiation within an optomechanical cavity [30,31]. We can establish the unforced behaviour of the charge oscillations for a given amplitude using equations (19)- (22) and treating the phase as constant (φ = 0). The unforced SSET charge oscillations have a frequency which increases linearly with the amplitude, A. For sufficiently small E J , we can approximate the frequency by and the decay rate is /2. Thus as the amplitude of the resonator increases, the number of oscillations in the SSET charge degrees of freedom during each resonator period increases accordingly.

Effective damping
Having calculated the response of the SSET charge to the periodic driving provided by the resonator, we can now calculate how the resonator, in turn, responds to the dynamics of the SSET. Our assumption is that the change in the amplitude of the resonator oscillations is much slower than the oscillations themselves. We can therefore average the effect of the SSET charge dynamics on the resonator over a single resonator period, over which time the amplitude of where the bar on the amplitude,Ā, indicates that the equation is only valid on time scales longer than 2π/ω. In terms of the Fourier series, the integration eliminates all of the terms except those with n = ±1, hence the equation of motion forĀ can be written as The second term can be interpreted as an amplitude dependent damping arising from the interaction with the SSET [10,11], γ SSET (Ā)Ā = 2ωx s Im[p 1 11 (Ā) + 2p 1 22 (Ā)] as it appears on the same footing as the external damping rate, γ ext . However, when γ SSET (Ā) < 0 it means that the charges transfer energy to the resonator. The resonator has a constant amplitude (but not phase) whenever dĀ/dt = 0, and hence supports limit cycle solutions whenever this condition is met forĀ = 0. Figure 3 shows 3 γ SSET (Ā)Ā calculated numerically as a function ofĀ, along with various values of γ ext . Limit cycle solutions exist for the values ofĀ where the curve of −γ SSET (Ā)Ā is crossed by a given line representing γ extĀ . The stability of the limit cycle solutions depends on the gradient of dĀ/dt with respect toĀ in the usual way, hence a limit cycle is stable where this gradient is negative [32]. Thus we can see by inspection from figure 3 that for the case where γ ext / = 0.001, there are three possible limit cycle solutions two of which are stable (the largest and smallest amplitude ones).
It is clear from figure 3 that it is the oscillations in γ SSET (Ā) as a function ofĀ that lead to the existence of more than one limit cycle solution 4 for sufficiently weak γ ext . The oscillations in γ SSET (Ā) can in turn be understood as arising from changes in the commensurability of the SSET charge oscillations with the resonator frequency as the amplitude is increased. This effect can be seen clearly by comparing the oscillations in figure 2 with those in figure 3. At the first peak in −γ SSETĀ , the charge undergoes one oscillation during the resonator period, a number which increases by two at each subsequent peak in −γ SSETĀ .

Approximate analytic solution
Although the SSET charge is bounded within a narrow range of values, its oscillations (as shown in figure 2) nevertheless resemble those of a driven harmonic system. Essentially this is because the parameters chosen are such that the Josephson energy is not sufficient to allow the charge to saturate over the time scale of the oscillations. This similarity suggests that it should be possible to find an approximate solution by reducing equations (19)- (22) to an equation for an appropriately driven harmonic oscillator. We can then follow the approach used in [30] where the harmonic dynamics of an optical cavity coupled to an oscillating mirror was analysed.
For an uncoupled SSET [11,14], it is straightforward to show that in the limit where E J /h 1 the populations of the upper charge states, p 11 and p 22 always remain much less than unity. Furthermore, examining the numerical evolution of the mean field equations we find that this remains true in the coupled system. This suggests that we can simplify the analysis in the limit where E J /h 1 by neglecting the p 11 and p 22 terms in equation (22). Making this approximation, and again assuming that the resonator damping is much slower than the other time scales, we find that equations (21) and (22) reduce to the equations of motion for a damped harmonic oscillator with a time dependent frequency term. Using p 02 = α + iβ, we obtain a single complex equation of motion, which can be solved by use of a Fourier series, p 02 = e iθ e iωntpn 02 where the tilde indicates that the Fourier series includes a global phase shift θ(t) = z sin ωt and z = Ax s /x 2 q . The value ofp n 02 can be written in terms of Bessel functions of the first kind,p n 02 =ψ n J n (−z), where the parameterψ n is given by,

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
As before, the damping is calculated from the Fourier coefficients of p 11 and p 22 , which in turn can be calculated from β n using equations (23) and (24), In order to make use of this Bessel function solution, we must convert between Fourier series with and without global phase shifts. We find, This approximate solution can be used to calculate p 11 (t), p 22 (t) (figure 2) and hence the effective damping ( figure 3). It is noticeable in figure 3 that even for our relatively large choice of E J /h 0.4, the analytical approximation agrees well with the numerics at large amplitudes (i.e. for A/x s 1). This is because as the amplitude of the resonator increases, the driven oscillations which develop in the SSET charge degrees of freedom become progressively faster and the populations p 11 (t) and p 22 (t) become ever smaller as can be seen in figure 2. Using smaller values of E J leads rapidly to better agreement at small amplitudes.

Comparison with numerical results
Much of the usefulness of our analysis of the mean field equations rests on the degree to which this simplified description of the SSET-resonator system faithfully reproduces the behaviour seen in the full numerical solution of the master equation [12]. The mean field equations allow us to calculate the amplitude of limit cycles in the resonator dynamics as a function of all the various parameters of the system, to determine which of them are stable, and calculate the associated SSET current. Of course the mean field description does not describe the noise in the SSETresonator system and hence it cannot tell us the degree to which a particular limit cycle may be occupied in regions of parameter space where there is more than one stable limit cycle (or a stable fixed point solution coexists with one stable limit cycle).

Size of limit cycles
We begin by considering sufficiently weak couplings (for a given value of γ ext ) that the resonator is limited to at most a single limit cycle state [12] and examine when such states develop and their sizes as a function of the detuning from resonance, E, and the relative speed of the resonator, ω/ . In figure 4 the size of the stable limit cycles calculated using the mean field equations (i.e. using equation (31)) is compared with the numerical solution of the full master equation as a function of E for various resonator frequencies.
The numerical solution of the master equation gives us the full steady-state density matrix of the system, ρ ss , in terms of which the probability of finding the resonator in a given Fock state, |n , is simply P(n) = Tr[|n n|ρ ss ]. Because of the ensemble average and the presence of phase noise, limit cycles do not appear as periodic features in the dynamics of the density matrix, but they can be identified as peaks in the distribution P(n) above the ground state energy (i.e. n > 0). (These peaks typically correspond to individual rings in the associated Wigner function representation of the density matrix [12].) We have used two quantities from the P(n) distribution to compare with the mean field results: the average number of resonator quanta, n = n nP(n), and the n-value of the peak in the distribution, which we define as n lc . In order to compare these values with the stable limit cycle amplitudes calculated in the mean field theory, we express the latter in terms of resonator quanta, n mf = A 2 (mω/2h).
It is clear from figure 4 that the mean field equations prove to be a rather good predictor of the locations and sizes of the limit cycles in the full system dynamics as can be seen by the relatively good agreement between n mf and n lc . It is interesting to note that a comparison of the sizes of the limit cycles with n works well in regions where the resonator energy distribution P(n) is relatively concentrated. Thus n mf is quite close to n when there is a single, large, limit cycle state present, but not in the vicinity of the transitions where the limit cycles begin to form or in regions where a limit cycle state coexists with a stable fixed point (bistable regions).
Despite the generally good agreement between n mp and n lc , figure 4 reveals that there are small differences in the locations of the values of E at which the limit cycles appear (or disappear) predicted by the mean field equations and the master equation. This can be partly explained by the fact that our analysis is most appropriate when the limit cycles are large on a scale set by x s (as discussed in section 4). However, the differences in onset points for the limit cycles exist not just in cases where the limit cycles grow continuously from zero (around E = 0), but also when they emerge with a relatively large radius. An interesting possible explanation for these observations is that the exact location of the dynamical transitions may in fact depend quite sensitively on the noise in the system-an element which is of course missing in our mean field analysis [22,13].
We now turn to compare the mean field analysis with numerical predictions for a given detuning, E, as the SSET-resonator coupling is increased. For resonator frequencies of the order of the quasiparticle decay rate, numerical solution of the master equation showed that increasing the coupling could lead to a sequence of transitions marked by the appearance of increasing numbers of peaks in the steady-state distribution P(n) [12].
In figure 5 we plot the size of the stable limit cycles calculated both numerically using the Fourier series solution of the mean field equations and using the Bessel function expression (equation (34)) as a function of κ 1/2 for ω/ = 1. The predictions of the mean field equations are compared with the locations of the peaks in the numerically calculated resonator distribution, P(n). The mean field calculation shows good qualitative agreement with the full numerics, showing a series of bifurcations and multiple limit cycles as the coupling is increased. For weaker coupling, κ 1/2 0.1, the mean field calculation accurately predicts the size of the limit cycle. Notice, however, that the Bessel function approximation for the limit cycle differs appreciably from the full mean field for the first limit cycle, but otherwise matches the mean field result closely. As before, reducing the value of E J improves the accuracy of the Bessel function approximation.
The appearance of successive stable limit cycle solutions as the coupling is increased in figure 5 is readily understood in terms of the analysis of the charge oscillations given in section 4.1 and the associated oscillations in γ SSET (A) (illustrated in figure 3). From equation (29), we see that increasing the SSET-resonator coupling (i.e. increasing x s ) increases the frequency of the charge oscillations thus changing their commensurability with the mechanical period. This effectively compresses the oscillations of γ SSET (A) as a function of A (i.e. they occur with a progressively smaller period measured in terms of amplitude). Thus for fixed γ ext , increasing the coupling means that more and more stable limit cycle solutions occur and those already present move to smaller and smaller sizes.

SSET current
It is also possible to calculate the current through the SSET using the mean field equations. The current is generated by quasiparticle tunnelling out of the states |1 and |2 , which leads to a time dependent tunnel current across the right junction, I(t) = (p 11 (t) + p 22 (t)). Therefore, when the resonator is in a limit cycle state of a particular amplitude, the corresponding oscillations in p 11 (t) and p 22 (t) (see figure 2) will be passed on to the tunnel current. The average current (defined as either an average over one period of mechanical oscillation, or over an ensemble of systems) of course will not reflect these oscillations, but it will depend on the amplitude of the resonator's motion. For a given resonator amplitude, the corresponding current is given by the Fourier coefficients, I /e = (p 0 11 (A) + p 0 22 (A)). The average current calculated using the mean field equations only agrees well with that calculated using the full master equation when the steady state of the latter predicts a narrow width to the distribution P(n). However, we expect that the current calculated within the mean field picture will be useful in understanding the current noise in the SSET [33]. For example, we might expect signatures of the different frequencies of charge oscillations corresponding to different resonator states to appear in the current noise spectrum.

Practical implications
The mean field analysis can be used to provide simple estimates of the kinds of dynamical instabilities that may be seen in a particular experiment. Recent experiments on a 22 MHz nanomechanical resonator coupled to a SSET were carried out [16] in the regime where ω/ 1. Previous calculations [11,13,14] showed that it should be possible to see the transition from the fixed point to a limit cycle state in such systems and a similar conclusion is reached using the mean field equations. Indeed some evidence for instabilities in the resonator motion was found in these experiments, although the primary focus of the work was on the stable regimes [16,34].
It is natural to use the mean field approach to estimate whether the regime where the resonator dynamics involves more than one limit cycle solution is accessible experimentally. For the nanomechanical resonator and SSET parameters in the experiments of Naik et al [16] we find that the wide separation of time-scales, ω/ 1, means that the values of the SSET-resonator coupling and γ ext necessary to reach this regime are well beyond currently achievable values. However, the interaction between the SSET and resonator is much stronger when ω/ ∼ 1, which suggests that the observation of more than one limit cycle may be possible when the resonator and quasiparticle decay timescales are more closely matched.
Although nanomechanical resonators with frequencies ∼1 GHz have been produced, [35] 5 integration with SET electronics has only been achieved for resonators up to ∼100 MHz [7]. However, it is possible to reduce the quasiparticle tunnelling rate, , substantially by increasing the resistance of the relevant tunnel junction and rates 2 × 10 8 s −1 have been demonstrated [25]. Therefore, as an example, we consider a nanomechanical resonator with fundamental frequency ω/2π = 100 MHz coupled to a SSET in which all the charge processes are rather slower than usual. Quasiparticle tunnelling in the SSET is assumed to occur across a junction with very high resistance, 5 M , and a Josephson energy at the other junction which is tuned (e.g. using the method employed by Nakamura et al [25]) to be E J = 1 × 10 −3 eV ds . For the other parameters of the system we choose values which are within the general range explored in recent experiments [7,15,16]. The other SSET parameters we assume are E C = 175 µeV, eV ds = 700 µeV, and = 200 µeV. For the resonator we assume a mass m = 6.8 × 10 −16 kg, with a SSET-resonator separation d = 100 nm, a SSET-resonator capacitance C g = 34 aF and coupling voltage V g = 1 V. We note that this choice of parameters takes us outside the limit E J /h 1 and we must rely on the numerical solution to the full Fourier series (rather than use the Bessel function approximation).
In figure 6, we plot the effective SSET damping as a function of resonator amplitude for a detuning E = −1 × 10 −3 eV ds . From figure 6, we see that we first get two stable limit cycles around γ ext / 2.5 × 10 −4 , corresponding to a resonator quality factor of 2.9 × 10 4 , which should be accessible experimentally. The observation of multiple limit cycles of course requires that the fluctuations in resonator amplitude are not so large as to wash out the difference in amplitudes between the cycles. Although a full calculation of the noise in the system is beyond the mean field theory as presented, we can make an estimate of the length scale of the fluctuations due to external thermal noise which, for the parameters in figure 6 and a temperature of 30 mK we find to be about δA 136x s i.e. about 10% of the separation in amplitude of the limit cycles. This suggests that multistability will occur within a region accessible by experiment though the conditions required to see it are quite demanding.

Conclusions
The mean field analysis presented here provides a simplified description of the SSET-resonator system. In particular, the mean field approach provides a good description of the dynamical instabilities which the resonator can undergo. For relatively weak SSET-resonator couplings the mean field equations describe the onset and size of the first limit cycle state quite accurately for a wide range of resonator frequencies. As the coupling is increased, the mean field equations predict the emergence of a succession of limit cycle states of different sizes. Numerical calculations based on the full master equation reveal that although the mean field equations become progressively less accurate as the coupling is increased they still give a good qualitative description of the dynamical behaviour. Furthermore, the mean field analysis allows us to relate the appearance of further limit cycles to changes in the commensurability of oscillations in the SSET charge with the period of the resonator.
is simply that of a harmonic oscillator, For the bias configuration we consider, Josephson tunnelling is only important between the left lead and the island as tunnelling between the right lead and the island involves energy differences ∼ 2eV ds E J because of the bias voltage applied (we neglect the possibility of the resonator having large enough energy to assist in overcoming this barrier). Thus we include just the coherent tunnelling between the island the left lead (i.e. coherent tunnelling changes N but not n), The quasiparticle Hamiltonian is given by [17] H qp = α=L,R,I k,σ kα c † kασ c kασ , (A.5) where the sum α runs over the three pieces of superconductor (left lead L, right lead R, and island I) and the operator c † kσL creates a quasiparticle in the left lead with energy kL , momentum k and spin state σ. Quasiparticle tunnelling between the island and the leads is described by the Hamiltonian, creates a quasiparticle in state k in lead j and destroys one in state q in the island with tunnelling amplitude T kq . The operators e ±iφ j /2 describe the associated change in the macroscopic charge variables, in terms of which they are written The master equation for the macroscopic island charge N and resonator is obtained by tracing out the quasiparticle degrees of freedom and the count variable n. The procedure is very similar to that used to derive the analogous master equation for a normal-state SET [29]. The derivation uses the standard Born-Markov approach [21] which involves the assumptions that the tunnelling Hamiltonian H T is a weak perturbation on the quasiparticle Hamiltonian and that the quasiparticles relax back to their unperturbed distributions after tunnelling faster than any other time-scale in the problem.
Further simplification is achieved by limiting the analysis to include only those states which are accessible to the system in the zero temperature limit. The charging energy of the SSET island, E c is typically much larger than E J , hence the Josephson tunnelling can be neglected for all states except the two charge states which are selected at the JQP resonance (by tuning 7 N g , see equation A.2) to be almost degenerate. For simplicity we consider states |0 and |2 , which differ by one Cooper pair, corresponding to resonance (for a fixed gate) at N g = 1. For quasiparticle decay to occur, the energy gained when a particle tunnels from the island to a lead must be 2 , where is the superconducting gap [14]. At the JQP resonance the voltage V ds is chosen so that only two processes are allowed: tunnelling from the island into the right lead between states |2 and |1 , and between |1 and |0 . The displacement of the resonator produces changes in the electrostatic energy differences involved in quasiparticle tunnelling which are assumed to be small with respect to the values for a fixed gate. This is the reason for our choice of a dimensionless of coupling, κ = mω 2 x 2 s /eV ds 1, which measures the energy associated with coupling to the resonator in terms of the typical energy scale associated with the unperturbed quasiparticle tunnelling rates. The question of when other charge states become accessible (i.e. for very strong SSET-resonator couplings or high enough temperatures), and what effect they have on the dynamics is an interesting one, but we do not consider it here.
Taking all these factors into account, we obtain a master equation confined to the space of the three charge states |2 , |1 and |0 , where (E) is the quasiparticle tunnel rate for energy E (this is the energy gained by the quasiparticle when it tunnels from island to lead) with the relevant energies for quasiparticle tunnelling from |2 to |1 and from |1 to |0 given by [14], respectively. We have assumed that the position dependent changes in the energy differences are small with respect to the typical energy scale eV ds , so that we can expand the tunnelling rates to first order [14], The master equation we have obtained (equation (A.10)) can be analysed numerically quite easily. However, it is useful to make further simplifications which make it easier to identify the essential physics of the system. The off-diagonal elements of the density matrix 0|ρ|1 , 1|ρ|2 (together with their complex conjugates) decouple from the dynamics of the other parts of the density matrix and can be dropped, together with the E 1 term in H co , as they play no role in the resonator dynamics. We also assume that the difference between the two quasiparticle decay rates and their variations with bias point can be neglected. While the difference between these rates (and their dependence on the bias point) may be important in the analysis of a given experimental system, it is not essential to our theoretical analysis which seeks to describe the basic features of the system in the simplest possible way.
Finally, we also drop the position dependent parts of the quasiparticle tunnelling terms. All the important dynamics arises from the position dependence in the coherent part of the master equation which modulates the detuning from resonance E (which can be arbitrarily small). In contrast, the position dependent coupling which appears in the quasiparticle decay terms is expected in practice to be much weaker than the coherent position coupling and provide only a small modulation of the quasiparticle rates [11,13,14]. Although numerically (within the master equation formalism) we find that including the position dependence of the quasiparticle rates can eventually lead to quite large changes in the steady state distribution function for sufficiently large κ, the general pattern and range of resonator behaviours (including the existence of regions where the resonator state is number-squeezed) remains essentially the same. Furthermore, corrections arising from the position dependence of the quasiparticle rates become progressively smaller as R J /(h/e 2 ) is increased [11,13]. Note, however, that the parameters used in the main text (and in [12]) are chosen to best illustrate the behaviour of our simplified model and we have not attempted to in addition minimize the corrections that would arise if the position dependence of the quasiparticle rates were included.

Appendix B. Stability analysis of the fixed point
In this appendix we describe the analysis of the fixed point solutions of the mean field equations (equations (11)- (16) It is obvious from equation (B.1) that there is only a single solution for β in the limit κ → 0. We have checked numerically that there remains only one physically acceptable solution over the range of parameters studied here. In order to establish the stability of the fixed point solution we need to calculate the Jacobian matrix of the mean field equations at the fixed point [32]. The stability of the system is then determined by the eigenvalues of this matrix, λ, which are solutions of the characteristic equation which can be written as λ 6 + a 5 λ 5 + a 4 λ 4 + a 3 λ 3 + a 2 λ 2 + a 1 λ 1 + a 0 = 0. (B. 3) The coefficients a 0 , a 1 , etc. are functions of the various parameters of the system determined by taking the determinant of the appropriate Jacobian. Whilst it is possible to determine the stability of the system by simply calculating all the eigenvalues numerically it turns out that we can establish the stability of the fixed point solution using a somewhat simpler approach. Assuming we start from a stable region and slowly change one of the parameters of the system, a limit cycle of the resonator begins to form when a pair of complex eigenvalues cross the real axis (Hopf bifurcation) [22,32,36]. The locations of these transitions can be determined by using the fact that at these points the characteristic equation must support a solution of the form λ = iν (with ν a positive real number). This requirement can be re-expressed in terms of the condition on the coefficients [22,36] f = a 5 a 2 1 + a 5 (a 0 a 3 − a 1 a 2 ) a 5 (a 0 a 5 − a 1 a 4 ) + a 1 a 3 2 − a 3 a 2 1 + a 5 (a 0 a 3 − a 1 a 2 ) a 5 (a 0 a 5 − a 1 a 4 ) + a 1 a 3 + a 1 = 0. (B.4) The regions of stability marked by arrows in figure 4 were determined by evaluating f for each set of parameters. In stable regions f < 0, changing sign at the Hopf bifurcations. Notice that our analysis is centred on the fixed point solution; it tells us nothing about the possible coexistence of stable fixed point and limit cycle solutions and hence does not describe the regions of bistability in figure 4.