A Numerical Method

Current noise of a superconducting single electron transistor coupled to a resonator

Abstract

We analyze the current and zero-frequency current noise properties of a superconducting single electron transistor (SSET) coupled to a resonator, focusing on the regime where the SSET is operated in the vicinity of the Josephson quasiparticle resonance. We consider a range of coupling strengths and resonator frequencies to reflect the fact that in practice the system can be tuned to quite a high degree with the resonator formed either by a nanomechanical oscillator or a superconducting stripline fabricated in close proximity to the SSET. For very weak couplings the SSET acts on the resonator like an effective thermal bath. In this regime the current characteristics of the SSET are only weakly modified by the resonator. Using a mean field approach, we show that the current noise is nevertheless very sensitive to the correlations between the resonator and the SSET charge. For stronger couplings, the SSET can drive the resonator into limit cycle states where self-sustained oscillation occurs and we find that regions of well-defined bistability exist. Dynamical transitions into and out of the limit cycle state are marked by strong fluctuations in the resonator energy, but these fluctuations are suppressed within the limit cycle state. We find that the current noise of the SSET is strongly influenced by the fluctuations in the resonator energy and hence should provide a useful indicator of the resonator’s dynamics.

pacs:
85.85.+j, 85.25.Hv, 73.50.Td

I Introduction

The dynamics of a resonator coupled to a superconducting single-electron transistor (SSET) is rather rich with a range of different behaviors expected to occur. Recent experiments in which the resonator was formed by a nanomechanical beam demonstrated ultra-sensitive displacement detection and cooling of the mechanical motion LaHaye et al. (2004); Naik et al. (2006). In another experiment Astafiev et al. (2007) the resonator was formed by a superconducting stripline and the SSET was used to drive it into a laser-like state Rodrigues et al. (2007a); Bennett and Clerk (2006); Hauss et al. (2008) of self-sustained oscillation.

A SSET consists of a superconducting island which is connected to two superconducting leads via tunnel junctions and capacitively coupled to a gate electrode, as shown schematically in Fig. 1. The gate electrode forms part of the resonator, either as a metal wire deposited on top of a nanomechanical beamNaik et al. (2006) or by forming the central electrode of a superconducting striplineAstafiev et al. (2007). Depending on the gate and bias voltages applied, the SSET can support a wide variety of current carrying processes Fitzgerald et al. (1998). Here we focus on a particular current resonance, the Josephson quasiparticle (JQP) cycle Averin and Aleshkin (1989); Fulton et al. (1989); Aleshkin and Averin (1990). At the JQP resonance current flows via a combination of the coherent tunneling of a Cooper pair at one of the junctions followed by the successive tunneling of two quasiparticles across the other junction. The center of the resonance occurs when the electrostatic energy of the two states linked by Cooper pair tunneling is the same. When the SSET is biased away from the center of the resonance the charges flowing through the SSET can either absorb energy from the resonator Clerk and Bennett (2005); Blencowe et al. (2005); Naik et al. (2006) or emit energy into it Lu et al. (2002); Clerk and Bennett (2005); Rodrigues et al. (2007b); Bennett and Clerk (2006); Astafiev et al. (2007). Interestingly, the charge dynamics in the JQP cycle is closely related to that of another mesoscopic conductor namely a double quantum dot in the Coulomb blockade regime and so-called wide bias limitBrandes and Lambert (2003); Aguado and Brandes (2004a).

When the SSET is detuned from resonance in such a way as to emit energy into the resonator the latter is effectively pumped by the flow of charges. For sufficiently strong coupling the resonator can be driven into a state of self-sustained oscillation. Useful analogies can be made between the SSET-resonator system and quantum optical systems such as the laserBennett and Clerk (2006); Rodrigues et al. (2007a); Astafiev et al. (2007); Marthaler et al. (2008). In particular, there are a number of similarities between the predicted dynamics of a resonator driven by a SSET and a micromaser. In a micromaserScully and Zubairy (1997) a superconducting cavity interacts with a stream of two-level atoms which pass through it one at a time. The ordered flow of atoms which interact one at a time with the cavity leads to an interesting range of effects which are not seen in standard lasers such as the existence of a whole set of dynamical transitions and the regimes where the cavity is driven into non-classical statesScully and Zubairy (1997); Walther et al. (2006). In the SSET-resonator system the resonator interacts with only one pair of charges moving through the SSET island at any one time, but the SSET current is not independent of the resonator dynamics. Nevertheless, features similar to the micromaser such as the existence of a sequence of dynamical transitions and the possibility of driving the resonator into non-classical states are predicted to occur for the SSET-resonator systemRodrigues et al. (2007b, a); Marthaler et al. (2008).

Figure 1: (a) Schematic diagram of the system consisting of a superconducting island formed by two tunnel junctions and a gate capacitor incorporating a resonator. (b) The JQP cycle involving coherent Cooper pair oscillations between the and charge states and incoherent quasi-particles tunneling to go from the state to the state via the state.

In the micromaser system the state of the atoms emerging from the cavity provides the information about the cavity dynamicsWalther et al. (2006). Similarly, The SSET current provides a natural source of information about the dynamics of the resonatorAstafiev et al. (2007). However, as with many mesoscopic systemsBlanter and Buttiker (2000), the fluctuations in the current contain much more information about the dynamics of the system than the average current alone. A number of recent studiesNovotný et al. (2003); Pistolesi (2004); Armour (2004); Flindt et al. (2004, 2005); Haupt et al. (2006); Usmani et al. (2007); Merlo et al. (2008); Doiron et al. (2008); Bennett and Clerk (2008) have shown how the current noise of a nanoelectromechanical system can be used to infer quite a lot about the dynamics of the system. For example, it has been recognized that a bistability in the dynamics of the mechanical resonator can lead to an extremely large peak in the current noiseFlindt et al. (2005); Usmani et al. (2007). In the case of the SSET-resonator system it has already been shown that the onset of self-sustained oscillations in the resonator can be associated with strong features in the SSET current noiseBennett and Clerk (2006).

In this paper we present a systematic study of the current characteristics of a SSET when it is coupled to a resonator as a function of the SSET-resonator coupling strength, the choice of SSET bias point and the frequency of the resonator. In particular, we use a numerical approach based on the master equation to explore the relation between the resonator’s state (measured by the average resonator energy and associated variance) and the SSET current and zero-frequency current noise. For sufficiently strong coupling, strong correlations develop between fluctuations in the resonator energy and the current noise. This means that measurements of the current noise could provide a very useful probe of the resonator dynamics.

In addition to the full numerical calculations, we also use a series of simpler approximate methods to gain more insight into the coupled dynamics of the system. When the SSET-resonator coupling is sufficiently weak, the SSET acts on the resonator like an effective thermal bath. In this regime we use a mean field approach that allows us to include information about the resonator dynamics as well as the correlations between the SSET and the resonator progressively and hence discover how these affect the current noise without influencing the average current. In the strong coupling regime, where the resonator undergoes oscillations driven by the current, we use an eigenfunction expansion of the Liouville operator in the master equationBriegel and Englert (1993); Jakob and Stenholm (2004); Flindt et al. (2004); Brandes (2008) to understand the current noise. In the vicinity of a bistability the current noise is dominated by the slow switching of the system between the two effective states of the system which is manifested by one eigenvalue for the Liouvillian which is much smaller (in magnitude) than all the othersFlindt et al. (2005). Interestingly, we find elsewhere that the noise can also be approximated quite well using a single term in the eigenfunction expansion even when a wide separation between the smallest few eigenvalues does not exist.

The organization of this paper is as follows. In Sec. II we introduce the master equation we use to model the SSET resonator system. We also describe how the steady state properties of the resonator and the current noise can be calculated numerically. Then in Sec. III we present calculations of the SSET current and zero-frequency current noise together with the associated resonator energy and energy fluctuations for a wide range of system parameters. We then focus on the weak coupling regime in Sec. IV where we present details of how simple models based around the mean field equations of the system can be used to understand the current and noise in this regime. Then in Sec. V we consider the regime where the coupling is strong enough to drive the resonator into limit cycle states. We use eigenfunction expansions of the relevant Liouville operator to explore the extent to which the presence of a very slow timescale in the resonator motion affects the current noise. Finally in Sec. VI we draw our conclusions. Appendixes AC contain further details on certain aspects of the calculations.

Ii Master Equation Formalism

In the vicinity of the JQP resonance Averin and Aleshkin (1989); Choi et al. (2001a, 2003) the SSET island is confined by charging effects to one of three charge states as shown in Fig. 1. These states correspond to the presence on the island of no excess charges, , one Cooper pair , or one quasiparticle, . The master equation describing the SSET and resonator at the JQP resonance is given by Rodrigues et al. (2007a, b),

(1)

The first term describes the coherent evolution of the density matrix under the Hamiltonian while the second and third terms describe the dissipative effects of quasiparticle tunneling and the resonator’s environment respectively. The SSET operators are defined in terms of the three accessible charge states,

(2)

The Hamiltonian, , written in terms of these operators takes the form,

(3)

where is the detuning from the JQP resonance, is the Josephson energy and the resonator has frequency , mass , momentum operator and position operator . The final term represents the linear coupling of the resonator to the charge on the SSET island. The length scale is the shift in the resonator position due to the addition of a single electronic charge to the island. The coupling strength is conveniently expressed in terms of the dimensionless parameter , where is the drain source voltage and the electron charge. Note that here we use the language and notation appropriate for a nanomechanical resonator, but the Hamiltonian takes essentially the same form for a superconducting stripline resonatorAstafiev et al. (2007).

Quasiparticle decay at the right hand junction is described by the superoperator ,

(4)

where is the quasiparticle tunneling rate and is the anticommutator Rodrigues et al. (2007b); foo (a). For simplicity, we have neglected both the differences between the rates for the two quasiparticle decay processes and the (weak) dependence of the rates on the position of the resonator Rodrigues et al. (2007b). The final term in Eq. (1) represents the damping of the resonator by its external environment:

(5)

where is the damping rate and where is the temperature of the resonator’s surroundings.

The whole master equation can be represented by the single superoperator , known as the Liouvillian. The Liouvillian operates in Liouville space where a Hilbert space operator becomes a vector and both pre- (left) and post- (right) multiplication of the operator can be represented by an appropriate matrix multiplying  Blum (1996); Briegel and Englert (1993); Jakob and Stenholm (2004); Flindt et al. (2004); Brandes (2008). The inner product for two vectors in Liouville space is defined as . Using this notation Eq. (1) takes the form,

(6)

Since we are dealing with an open system, the Liouvillian is not Hermitian and hence has different right and left eigenvectors,

(7)
(8)

We choose to label the set of eigenvalues such that . Neglecting the possibility of degeneracy,Jakob and Stenholm (2004) we assume that the eigenvectors form a complete orthonormal set, . The solution to Eq. (6) can therefore be expanded in terms of the eigenvectors to give

(9)

where is the initial density matrix of the system. For a master equation with a well-defined steady state (such as the one we consider here) the lowest eigenvalue will be , a property which we used to obtain the second line above. The other eigenvalues must obeyJakob and Stenholm (2004) and the steady state density operator is . The normalization of is determined by , which gives , where is the identity operator (in Hilbert space). While corresponds to the steady state, the eigenvectors for each represent a change to the steady state density matrix that decays exponentially with rate .

The problem of finding the steady state density matrix is reduced to finding the right hand eigenvector of corresponding to the eigenvalue . By truncating the oscillator basis, Eq. (6) can be solved numerically to find the first few eigenvalues and eigenvectors of . Details of the numerical method and the approximations made are contained in Appendix A.

Our aim in this paper is to understand to what extent information about the dynamical state of the resonator becomes imprinted on the transport properties of the SSET. As well as calculating the current we also consider the zero frequency current noise, which is independent of the junction at which it is measured. We choose to calculate the noise at the junction at which the Cooper pair tunneling takes place (the left hand junction) as the current operator here is composed of system operators alone, which along with the Markovian nature of the master equation, allows the use of the quantum regression theorem Aguado and Brandes (2004a). (The results have also been calculated for the right hand junction and are in agreement.)

The symmetrized current noise at the left hand junction is defined as,

(10)

where the current operator at the left hand junction, , can be calculated by considering the flow of charge across the left hand junction Flindt et al. (2004). The operator is given by

(11)

For the current noise a symmetrized current operator can be defined in Liouville space,

(12)

Using the quantum regression theorem to evaluate the correlation function, the current noise can be written in terms of Liouville space operators as

(13)

In the zero frequency limit this has the solution Flindt et al. (2004)

(14)

where is the psuedoinverse of the Liouvillian and the projector . With this projection, the inversion takes place only in the space where is regular. The current noise is conveniently parametrized by the Fano factor, which is defined as

(15)

where is the Poissonian or shot noise limit corresponding to the current noise for a single tunnel junctionBlanter and Buttiker (2000).

Iii Features of the Current Noise

In this section we give a survey of the current and noise characteristics of the SSET and the corresponding resonator dynamics calculated numerically over a range of different parameters. In later sections we will provide a more detailed analysis of the most interesting regimes.

The SSET-resonator system is rather complex, even within the framework of the simple model we are using here. In particular, the behavior of the system depends on a rather large number of different parameters. Some of the system parameters such as the detuning and the SSET-resonator coupling can be changed during a particular experiment, while many of the other system parameters can be tuned by appropriate design of the device, including the frequency of the resonator which can be in the region of MHz for a mechanical deviceNaik et al. (2006) or of order GHz for a superconducting stripline resonatorAstafiev et al. (2007). We have not attempted a systematic survey of all feasible parameter regimes, but instead focus primarily on the effects of changing , and the resonator frequency.

We start by reviewing the characteristics of the SSET in the uncoupled limit . The current, , and Fano factor, , for a SSET tuned to the JQP resonance are given by Choi et al. (2001b, 2003):

(16)
(17)

The current has a peak at the center of the resonance , which has a width determined by and . Far from resonance the current Fano factor has a value of 2. This is because the rate at which the Cooper pairs tunnel onto the island is much slower than the quasiparticle decay rate and hence when a Cooper pair reaches the island it swiftly breaks up into quasiparticles. The two quasiparticle tunneling processes occur in quick successionChoi et al. (2003) (compared to the rate of Cooper pair tunneling) and hence the charge is effectively transferred in units of 2. However, close to the center of the resonance in the regime where (which we study here) there is a strong interplay between the coherent transfer of Cooper pairs and the quasiparticle tunneling which results in a suppression of the noise. This suppression is strongest at the center of the resonance where the coherent motion of Cooper pairs is most important.

Coupling a resonator to the SSET can significantly modify the behavior of the two individual systems. In particular, energy can be transferred between the SSET charges and the resonator. For low to moderate coupling the resonator reaches one of three types of steady-stateRodrigues et al. (2007a, b): a state in which the resonator fluctuates about a fixed point, a limit-cycle state where the resonator undergoes self-sustained oscillations and a ‘bistable’ state, where the two coexist. At larger couplings other states can be found for this system such as multiple limit cycles Rodrigues et al. (2007b), but these do not occur for the parameters used here. The different resonator states are readily identified from the steady state number state distribution of the resonator, , where is a Fock state of the resonator. The fixed-point state has a single peak in the distribution at , while the limit cycle has a peak at and we define the bistable state as having two peaks.

Figure 2: (color online). Average energy of the resonator as a function of the detuning from resonance and coupling strength for and . The dashed lines indicate transitions between dynamical states: for most of the range considered the resonator is in the fixed point state, but for large enough coupling a transition to the limit cycle state occurs close to the center of the resonance. The bistable region is the smallest and occurs for and .
Figure 3: (colour online). Average current () through the SSET as a function of the detuning from resonance and coupling strength. The dashed lines indicate transitions in the resonator’s state. The parameters are the same as in Fig. 2.
Figure 4: (color online). Resonator Fano factor, as a function of the detuning from resonance and coupling strength. The dashed lines indicate transitions in the resonator’s state. The parameters are the same as in Fig. 2 and the colors are on a log scale.
Figure 5: (color online). Current Fano factor, , as a function of the detuning from resonance and coupling strength. The dashed lines indicate transitions in the resonator’s state. The parameters are the same as in Fig. 2 and the colors are on a log scale.

The behavior of the resonator and the corresponding current characteristics of the SSET for a slow resonator, , are illustrated in Figs. 2-5. The average energy of the resonator is shown in Fig. 2 as a function of the detuning and the coupling for . We have chosen a junction resistance and the quasiparticle decay rate is taken to be . The Josephson energy is assumed to have a value so that throughout we will be in the regime where and hence the quasiparticles should be the dominant source of dephasing for the SSET island charge. The damping (which is somewhat higher than might be expected in experiment) is chosen to ensure numerical convergence and a small amount of thermal noise has been included .

The transitions between the three different dynamical states of the resonator are indicated by dashed lines in Fig. 2. For energy is transferred to the resonator and for strong enough coupling the resonator is driven into the limit cycle state which grows in size continuously as becomes more negative. However, for when is sufficiently negative () the resonator enters the bistable regime and then undergoes a transition back to the fixed point state in which the limit cycle disappears abruptlyRodrigues et al. (2007a). The corresponding behavior of the current is shown in Fig. 3. Although the current is clearly modified by the coupling to the resonator, it does not contain any clear signatures of the transitions in the resonator state.

An important measure of the resonator state is the Fano factor of the resonator occupation number, defined as where (where here is the number operator ). is plotted in Fig. 4. Unlike the average energy of the resonator, is strongly peaked around the transitions between the fixed point and limit cycle states, with the strongest feature occurring in the vicinity of the bistable regionRodrigues et al. (2007a).

The current Fano factor is plotted in Fig. 5. The behavior is rather complex, especially for relatively weak couplings, but overall it is clear that the current noise is a much better indicator of the presence of transitions in the resonator state than . The behavior is simplest for larger , well within the regime where the dynamical transitions occur and in this region we see well defined peaks in the noise at the transitions into and out of the fixed point state. The noise peak is particularly prominent in the case of the bistability. Although we have defined the bistable state on the basis of just the number of peaks in the distribution, rather than the coexistence of two well separated states it is certainly possible to find regimes where a true bistability in this sense exists (i.e. where the distribution not only has two peaks, but also has very small values for at least some range of the values between the peaks). For a well defined bistability the current noise can become extremely large, a phenomenon which has been shown to be due to the existence of a very slow timescale associated with the switching of the system between the metastable states of the system Isacsson and Nord (2004); Novotný et al. (2004); Usmani et al. (2007).

We now turn to consider the opposite regime of a fast resonator, . In this regime the coherent coupling between the resonator and the SSET is expected to give rise to well-defined features when the resonator frequency matches an eigenenergy of the SSET,

(18)

with an integer (for the relatively strong quasiparticle decay rates considered here means that ). Numerically we do indeed find that limit cycles occur at these resonance points, almost always via continuous transitions, with the higher order features () appearing at progressively larger couplings.

Figure 6 shows , and around the (one photon absorption) peak corresponding to with the dashed line marking the onset of the limit cycle state. This resonance was recently observed in experiments using a superconducting resonator Astafiev et al. (2007). It is clear that the peak in the current correlates well with the presence of the limit cycle, but this peak is present at the resonance even when the limit cycle is not formed, albeit with a very small size. However, the behavior of the noise shows a clear signature of the onset of the limit cycle with a distinctive peak structure forming along the transition lines for both and (Fig. 6b and c). Within the limit cycle regime both and show dips the size of which is an indication of how well defined the limit cycle is.

In practice it is also possible to build devices where the resistance of the junction where quasiparticle tunneling occurs is much larger than the quantum of resistanceAstafiev et al. (2007) (i.e. ). Increasing the resistance of the junction where quasiparticle tunneling occurs enhances the coherent interaction between the Cooper pairs and the resonator. Strictly speaking one would expect other effects which give rise to dephasing of the SSET charge (beyond just quasiparticle tunneling) to become relevant in this regime. Nevertheless, for our simple model we find that when and , the current noise at can become sub-Poissonian () and the resonator Fano factor can also drop below unity, indicating a non-classical resonator state. Interestingly, this regime is quite distinct from the case discussed in Ref. Rodrigues et al., 2007a (for which and ) where can also drop below unity. We note, however, that the thresholds for the sub-Poissonian regimes for and are not perfectly correlated: these two effects can occur at the same time or separately depending on the parameters chosenMerlo et al. (2008).

Figure 6: (color online). , and for , and ; the other parameters are the same as in Fig. 2. For the region within the dashed lines the resonator is in a limit cycle state and elsewhere it is in a fixed point state.

Iv Thermal Resonator

The simplest regime for the SSET-resonator system is that of very weak coupling where the resonator remains in the fixed point state for all values of the detuning . In this regime the effect of the SSET on the resonator dynamics is analogous to an additional thermal bath Clerk and Bennett (2005); Blencowe et al. (2005). In this section we consider in detail how the SSET current and noise are modified by the resonator in this regime and explore the extent to which the behavior can be understood in terms of simple models for the coupled SSET-resonator dynamics.

For sufficiently weak coupling the resonator’s steady state is a thermal (i.e. Gaussian) state to a good approximation Clerk and Bennett (2005); Blencowe et al. (2005). In this state the resonator position is determined by the average charge on the SSET island which in turn is proportional to the average (steady state) current flowing through the SSETBlencowe et al. (2005); Rodrigues et al. (2007b) (see also Appendix B),

(19)

The average occupation number of the resonator, , is given by a weighted sum of the contributions arising from the resonator’s thermal environment, , and the effective thermal bath it feels due to the SSET, ,

(20)

The weighting factors are the resonator damping rates due to the true thermal bath, , and the effective bath due to the SSET, . For a slow resonator (), both and are given by relatively simple analytic expressions Clerk and Bennett (2005); Blencowe et al. (2005); foo (b):

(21)
(22)

Both and are strongly dependent on the detuning . Furthermore, becomes negative for . However, for weak enough coupling the resonator is stabilized in a thermal state by the damping arising from the coupling to the external bath.

For weak SSET-resonator coupling the changes in the transport properties of the SSET due to the resonator are relatively small so it makes sense to examine just the difference between the values for the coupled and uncoupled cases. The change in the SSET current due to the resonator (calculated numerically) is shown in Fig. 7. We consider a slow resonator and very weak coupling so that although the SSET has quite a strong influence on the resonator state, the resonator nevertheless remains in a thermal state which is well described by Eqs. (20)–(22). From Fig. 7 we see that near the center of the resonance the current is suppressed by the resonator, but on either side of this there is an enhancement. The current noise is modified in a similar way to the current, but in the opposite sense, as shown in Fig. 8, thus there is an increase in the noise near to the resonance with a decrease on either side.

Figure 7: (color online). Change in current through the SSET as a function of for , , , and . The curves are labeled as num for the numerical results and fl for the change in the current calculated using Eq. (24). [Note that for these parameters varies from a value of 2 far from resonance to a peak value of 2.28 at . has maxima and minima of at .]
Figure 8: (color online). Change in the zero frequency Fano factor of the SSET due to the resonator. The curves are labeled as num for the numerical results, fl is obtained from Eq. (LABEL:eq:fl_FI), mean2 is calculated using the second order mean field equations and mean3 using the third order mean field equations. The parameters are the same as in Fig. 7.

The simplest way of including the influence of the resonator on the SSET is to include the effect of fluctuations in the position of the resonator on the current. Because the resonator acts as a gate for the SSET island, a shift of the position of the resonator leads to an effective change in the detuning energy (Eq. 3). Hence, when the resonator position fluctuates so will the detuning energy. We can incorporate the effect of the mechanical motion into the expressions for the current and noise, (Eqs. 16-17), by calculating them for a fixed position, making the replacement , and then averaging over the resonator state. For the current Eq. (16) becomes:

(23)

Assuming the shift term is small we can perform a Taylor expansion and then take the average over the resonator. Keeping terms up to order and using Eq. (19), we obtain

(24)

where and the averages are taken over the (Gaussian) steady state probability distribution for the resonator. The value of is calculated using Eq. (20).

For the current noise we naively replace to obtain

(25)

After expanding to second order in and taking the average over the resonator state we can then calculate the corresponding Fano factor,

(26)

where .

Looking first at the current (Fig. 7), it is clear that Eq. (24) accurately describes the modification due to the presence of the resonator. Thus in this weak coupling regime where the resonator remains in a thermal state, the modification of the current is simply due to the shift in the resonator’s position (which gives the asymmetric shape) and a smearing out of the JQP current peak due to fluctuations in the resonator position. In contrast, we can see from Fig. 8 that for the Fano factor Eq. (LABEL:eq:fl_FI) does not capture the behavior correctly. Although the qualitative shape is the same with a central peak with dips either side, the curves do not match and the asymmetry of the numerical curve is in the opposite direction to that predicted by the simple model.

The reason for the disagreement in the current noise is that the simple model of a fluctuating gate neglects both the correlations between the electrical and mechanical motion and the dynamics of the resonator. The current noise (in contrast to the average current) is sensitive to the correlations between the SSET charge and the resonator motion and hence to describe it accurately we need to include them in some way. A straightforward and systematic way to include correlations and information about the resonator dynamics can be found using the mean field equations of the system, namely the equations of motion for the expectation values of the SSET and resonator operators. The mean field equations are generated in turn by multiplying the master equation by an operator (or product of operators) and taking the trace over the full system Rodrigues et al. (2007b). The mean field equations for the SSET resonator system are given explicitly in Appendix B up to second order.

The set of mean field equations for the SSET resonator system never forms a closed set with equations of motion for the operators always including some higher order operators. However, progress can be made by making a semiclassical approximation, in which correlations between certain operators are neglected so that the system of equations then becomes closed. Here we retain at least some of the SSET-resonator correlations by only applying the semiclassical approximation to products of higher order than required. For the set of second order equations we approximate products of three system operators, thus we make the substitution . Crucially the correlations between products of any two operators are retained.

The resulting set of equations is closed, but is non-linear because of the terms generated by the semiclassical approximation. However, by again using what we know about the steady state of the resonator in this regime [i.e. Eqs. 19 and 20] to replace the expectation values involving only the resonator operators by their steady state values, we can recover a linear set of equations (full details are given in Appendix B). This set of equations is then solved to obtain the steady state moments of the system using the same approach as that employed to solve the master equation (namely solving for the null eigenvector of the matrix of coefficients). The current is then obtained directly from the moments of the SSET operators, . The current noise can be calculated using a slightly more elaborate calculation which introduces an electron counting variable Aguado and Brandes (2004b), full details of which are given in Appendix C.

The results from the second order mean field equations for the current agree very well with the numerical calculation (and Eq. 24) as one would expect. For the current noise the second order mean field calculation is qualitatively correct as can be seen in Fig. 8, capturing the asymmetry in the numerical results though quantitative agreement is still lacking. This is not surprising as the second order mean field calculation only partly includes the SSET-resonator correlations and does not describe the dynamics fully. However, the mean field approach is readily extended to third order (i.e. the semiclassical approximation is only applied to products of four operators), thereby including higher order correlations and more information about the resonator dynamics. We find that the third order calculation leads to quantitatively correct results as shown in Fig. 8. However, we also note that reducing the coupling reduces the importance of the higher order correlations which the second order mean field calculation neglects. Figure 9 provides a clear illustration of this as it shows that the second order calculation becomes accurate for low enough .

Figure 9: (color online). Change in the zero frequency Fano factor of the SSET due to the resonator for . All other parameters and labeling of curves are the same as Fig. 8.

V Time Scales of the System

In this section we discuss the signatures of the resonator dynamics in the current and current noise when the coupling is strong enough to drive the resonator into limit cycle states. In particular we investigate how simple approximations to the noise based on the eigenfunction expansion in Liouville space (Eq. 9) can be used to give insights into the connections between the fluctuations in the resonator state and the current noise.

Figure 10 shows the current calculated numerically as a function of for three very different values of the resonator frequency. For low resonator frequencies (), the current is slightly suppressed at the center of the JQP resonance and enhanced further away. This is qualitatively the same as was seen for weak coupling in Sec. IV, even though now the coupling is larger so the resonator is driven into a limit cycle state. For the case small peak features are seen for at points corresponding to the and resonances in Eq. (18). In the high frequency case (), the current is greatly enhanced at the resonance and relatively unchanged elsewhere.

Figure 10: (color online). Current as a function of for different resonator frequencies . In each case the values of and have been chosen to ensure that the system reaches the limit-cycle state for at least some values of whilst still remaining at low enough energies to allow a numerical calculation. For , and ; for , and ; and for , and . The other parameters are the same throughout: .

Figs. 1113 show the current noise calculated numerically for the same parameters. For and the two peaks in the current noise correspond in both cases to a continuous transition from a fixed point state to a limit cycle at and the presence of a region of bistability near the second (larger) peak in . In between these two peaks the system is in a limit cycle state. For the case the two peaks in both correspond to continuous transitions (from fixed point to limit cycle state) with the resonator in a limit cycle state between the peaks.

v.1 Bistability

As discussed in Sec. III, the current noise contains more information about the dynamical state of the resonator than the current alone. For example, the presence of a dip in the noise between two strong peaks gives a clear indication that the resonator is actually in a limit cycle stateNovotný et al. (2004). The current noise can also tell us about the types of fluctuations present in the system, and the time scales over which these fluctuations decay. This analysis is particularly clear in the case of a bistability.

Current noise peaks for bistable regions in nanoelectromechanical systems, such as the charge shuttle, have been studied extensively Isacsson and Nord (2004); Novotný et al. (2004); Flindt et al. (2005); Usmani et al. (2007). The current characteristics of a conductor coupled to a truly bistable system (i.e. one with only two accessible internal states) can be described by a model specified in terms of four parameters: the (different) currents associated with the two states , and the switching rates between them of . The current and current noise for this two-state model take the simple form Flindt et al. (2005); Usmani et al. (2007)

(27)
(28)

where is the variance in the current.

This two-state model can be applied to a more complex system in a bistable regime if the two metastable states are well separated so that the switching rate between the states is much slower than the other relevant time-scalesFlindt et al. (2005); Usmani et al. (2007). From Eq. (28) we can see how slow switching rates between the two states can lead to a large value for the current noise in this regime. However, we also note that when the two metastable states give rise to very different currents the large variance that results can also make an important contribution to the current noise.

For certain parameters the noise at the bistable transition in our system is very well described by this two state model (with the two metastable states being the fixed point and limit cycle). In practice this means that a single set of the four parameters can be found which allow us to fit the current and current noise to Eqs. (27) and (28) respectively. Furthermore, the same parameters can be used to calculate higher cumulants of the current which can also be compared with numerical results Flindt et al. (2005). The required parameters can be extracted as follows. The relative probabilities of the two states are obtained by inspection of the steady state probability distribution . Setting those elements of the steady state density matrix which correspond to just one of the two states to zero and recalculating the current then allows the currents and to be obtained. Finally, the sum of the rates can be determined by comparing the current noise (calculated numerically) with Eq. 28.

The two-state model can only be applied when a true bistability exists in the sense described in Sec. III (i.e. the distribution for the resonator steady state should have two peaks with a vanishingly small probability for some range of values in between) which we find generally occurs for . Using the methods described above, we found that the two state approximation can be used to describe the current and current noise for the bistable state seen around in Fig. 12 (), but not for the one around in Fig. 11 (), where there is significant overlap between the limit cycle and fixed point states. For the former case we obtained further confirmation that the two-state model could be applied by checking that the numerically calculated third cumulant agreed with that obtained using the two state model (an approach discussed in detail in Ref. [Flindt et al., 2005]) and also by checking that the smallest (non-zero) eigenvalue of the Liouvillian matched up well with the total rate (as we discuss below).

In an ideal experiment one would be able to monitor the current with sufficient time resolution to observe the slow switching between two distinct values of the current directly. However, measuring the current noise as well as the average current in a region where the theory predicts a bistability would provide convincing evidence if agreement was obtained. One could also make use of further generic predictions of the two-state modelUsmani et al. (2007), such as the presence of a Lorentzian peak in the finite frequency current noise (at zero frequency) with a width given by .

v.2 Eigenvector expansions

In general, we cannot describe our system in the vicinity of the dynamical transitions by a simple two-state model. As we have seen, even where the transition involves a region of coexistence between the limit cycle and fixed point states the states may not be well enough separated for a two-state model to apply. Near the continuous transitions between the limit cycle and fixed point states there are clearly not just two states involved. However, one element of the two state model which might be expected to apply more widely is the emergence of a single very slow timescale which dominates the current noise. In the case of the continuous transition such a slow timescale might result from the vanishing effective damping () of the system at the transition. In what follows we use the eigenvector expansion of the Liouvillian to investigate the extent to which the current noise at each of the dynamical transitions can be described by a single slow process.

We begin by rewriting Eq. (14) for the current noise in terms of the eigenvectors and eigenvalues of the Liouvillian,

(29)

This should be compared with a similar expansion for the variance in the current also for the left hand junction:

(30)

The variance is given by a sum over the same matrix elements as the current noise but this time unmodified by the eigenvalues, . Each of the eigenvectors of the Liouvillian describe a change to (or fluctuation away from) the steady state that decays with a purely exponential rate [see Eq. 9]. Thus, the matrix element can be thought of as the variance in the current due to a fluctuation of type . We then see that the current noise consists of a sum over the variances due to each type of fluctuation, each divided by the rate at which that fluctuation decays.

It is clear from Eq. (29) that if then we could expect the current noise to be dominated by the first term, which corresponds to the slowest timescale in the system. This is in indeed what happens when the system has a well-defined bistability. In this case an obvious connection can be made with an appropriate two state model (i.e. Eq. 28), with the relevant eigenvalue corresponding to the sum of the rates and the numerator gives the current variance, . More generally, it is not just a slow timescale that is important. For a single term in the eigenvector expansion to accurately describe the noise, the matrix element divided by the eigenvalue for must be much larger than for all .

In Figs. 1113 we compare the full current Fano factor with approximations using just the first term in Eq. (29). The peaks at the transitions are described quite well by just the first term in the eigenvector expansion. Away from the peaks, however, we find that the noise is not captured by the approximation based on the first eigenvector. It is particularly clear in Fig. 13 that something is missing from this approximation. The features that are simply due to the SSET alone are not captured, such as the dip at and the Fano factor of 2 far from resonance. We can understand this better by considering the meaning of the eigenvectors and eigenvalues of the LiouvillianBriegel and Englert (1993); Englert and Morigi (2002).

Figure 11: (color online). for , , , , and . The curve labeled num shows the numerical value of the noise, app is the approximate value of the noise using the first term in Eq. (29), and app5+ is the first five terms plus the noise for a SSET alone [Eq. 32].
Figure 12: (color online). for , , , , and . The curve labeled num shows the numerical value of the noise, app is the approximate value of the noise using the first term in Eq. (29), and app5+ is the first five terms plus the noise for a SSET alone [Eq. 32].
Figure 13: (color online). for , , , , and . The curve labeled num shows the numerical value of the noise, app is the approximate value of the noise using the first term in Eq. (29), and app5+ is the first five terms plus the noise for a SSET alone [Eq. 32].

In the limit the resonator-SSET system becomes uncoupled and the eigenvectors and eigenvalues of the system can be expressed in terms of those of the individual subsystems, namely the SSET and the resonator. When the resonator is decoupled from the SSET it still remains coupled to the external bath and its smallest (non-zero) eigenvalues are integer multiplesEnglert and Morigi (2002) of . Thus the smallest of these eigenvalues corresponds to the energy relaxation rate of the resonator, , and hence we can infer that the corresponding eigenvector describes fluctuations in the resonator’s energy (something which we will justify further below). There are also a set of eigenvectors (and corresponding eigenvalues) that describe fluctuations in the SSET charge state. In the uncoupled regime the current noise of the SSET can be obtained using Eq. (29), with the sum running over just the SSET eigenvalues, though we already know the result will be given by Eq. (17).

For the coupled SSET-resonator system we can still identify the eigenvalues and eigenvectors as corresponding to one or other of the subsystems by looking at their behavior for large detunings (i.e. large ) where the systems are effectively decoupled. The first few eigenvalues, which correspond to the resonator, are shown (for the slow resonator case ) in Fig. 14 as a function of . These first few eigenvalues indeed converge towards , , for large detunings. Thus at least for large detunings the first eigenstate, , should therefore represent fluctuations which change the resonator energy. This can be confirmed by comparing the resonator variance to an expansion in terms of the eigenvectors,

(31)

where . The full calculation of the energy variance is compared with approximations based on the first term in the eigenvector expansion in Fig. 15. It is clear that only the first eigenvector is needed to describe the energy fluctuations for large detunings as we expect. However, the approximation based on the first eigenvector also describes the energy fluctuations rather well at the peaks where the transitions occur, but not in between where the resonator is in a limit cycle state. However, Fig. 15 also shows that we can describe the energy fluctuations throughout by using more terms in the eigenvector expansion.

Figure 14: (color online). The four smallest (non-zero) eigenvalues as a function of the detuning, . The parameters are the same as in Fig. 11. The eigenvalues differ from each other by more less than one order of magnitude throughout and converge towards integer multiples of for large .
Figure 15: (color online). Energy fluctuations of the resonator, , as a function of . The three curves show the full numerical calculation, num, and approximations using just the first term, app and the first five terms, app5, of the eigenfunction expansion [Eq. 31], respectively. The parameters are the same as in Fig. 11.

We are now in a position to understand why the calculation of the current noise using just the first term of the eigenvector expansion works as well as it does and to see how and why this can easily be improved upon. Comparing Figs. 11 and 15 it is clear that the single-eigenvector approximation to the current noise matches the numerical results well around the two peaks marking the transitions (between the fixed point and limit cycle states) where the first term in the eigenvector expansion also describes the energy fluctuations in the resonator accurately. The fact that the first term in the eigenvector expansion does not capture the current noise far from resonance is not surprising as it only describes fluctuations in the resonator state and does not include the fluctuations of the SSET degrees of freedom. We can easily obtain better agreement for large detunings by extending our approximation to include the contribution of the uncoupled SSET, . Better agreement within the limit cycle regime can be attained by using sufficient eigenvectors in our approximation to ensure that the fluctuations in the resonator energy are described accurately. Thus we arrive at our final approximate expression for the current noise

(32)

where should be large enough so that the corresponding number of terms can be used to calculate accurately [via Eq. 31]. In this case we find is sufficient, and the current noise calculated this way agrees very well at almost all points as shown in Figs. 1113. The one area where good agreement is still lacking using Eq. 32 is within the limit cycle region for , shown in Fig. 13. This is because we have simply used the uncoupled contribution to the current noise arising from the SSET eigenvectors. In fact these SSET terms are strongly modified due to the resonant absorption of energy by the resonator from the Cooper pairs at this point.

From these approximations it is clear that in the vicinity of the resonator transitions the current noise is largely determined by the slow fluctuations in the energy of the resonator. This is because the current depends in the first instance on the resonator position and hence on the latter’s energy (as this is slowly changing compared to its period). Thus the current fluctuations depend strongly on the fluctuations of the resonator energy, rather than those of higher moments of the resonator. Thus when depends on more than one eigenvector, the current noise does too.

It is important to note that even in the regions where including just the first term in the eigenvector expansion describes the current noise fairly well this is not simply because the associated eigenvalue is very much smaller than all the others. We can see from Fig. 14 that (for these parameters) an overwhelming difference between the slowest two eigenvalues never develops and from Fig. 16, that the relative size of the corresponding matrix elements is important in causing the first term in the eigenvector expansion to dominate.

Figure 16: (color online). Matrix element corresponding to the two smallest magnitude eigenvalues, . The parameters are the same as in Fig. 11. Large differences between the magnitudes of the two matrix elements are visible in the same regions as the peaks in the corresponding plot of the current noise, Fig. 11.

Vi Conclusions

We have investigated the current and current noise of a SSET coupled to a resonator and how they relate to the latter’s dynamics. The steady state properties of the system and the zero-frequency current noise are readily calculated using a numerical approach based on a representation of the master equation as a matrix equation in Liouville space. Overall we found that the current noise varies widely depending on the precise choice of SSET bias point, the resonator frequency and the strength of the SSET-resonator coupling. For sufficiently strong couplings, the SSET current noise is strongly influenced by the fluctuations in the resonator energy. In particular, the resonator energy displays strong fluctuations in the regions where transitions between dynamical states occur and this behavior is reproduced in the current noise. This means that measuring the current noise could provide clear signatures of dynamical transitions in the resonator.

In addition to the full numerical calculations, we used a range of approximate methods to provide further insights into the coupled dynamics of the system. For very weak SSET-resonator couplings the SSET acts on the resonator like an effective thermal bath. We found that in this regime mean field equations for the system operators provided a convenient way of establishing the importance of the SSET-resonator correlations and the resonator dynamics in determining the SSET current and noise. For stronger couplings we used eigenfunction expansions of the Liouvillian matrix to demonstrate the strong connection between the energy fluctuations in the resonator and the current noise. In many cases the current noise is well approximated by just a few terms in the eigenfunction expansion, but we found that it cannot be approximated accurately without including all of the eigenvectors that are needed to describe the energy fluctuations in the resonator state.

Acknowledgements

We would like to thank C. Flindt, J. Imbers and F. Pistolesi for a series of very useful discussions. We acknowledge funding from the EPSRC-UK under grants EP/E03442X/1 and EP/D066417/1.

Appendix A Numerical Method

This appendix describes in more detail the numerical method used to solve the master equation and the approximations that are made. To find the steady state of the system numerically the basis of the resonator must be truncated. External damping sets a limit on the resonator energy. We therefore use a Fock state basis for the oscillator truncated to states, where is chosen to be large enough that the probability for the resonator to have an energy larger than is negligible. In Liouville space the density matrix is a vector and is a matrix with dimensions .

To obtain the steady state density matrix we need the eigenvector corresponding to the zero eigenvalue, or null eigenvector, of the Liouvillian. In Sec. V we also require the first few nonzero eigenvalues with the corresponding right and left eigenvectors. Due to the large dimension, complex nature and unsymmetric structure of the Liouvillian the inverse iteration method is used Wilkinson (1965); Golub and van Loan (1996). Starting from a random vector the iteration is

(33)

where we have expanded in terms of the eigenstates of the Liouvillian and is the identity in Liouville space. Repeating the iteration, the value of will converge to the eigenvector of closest to . To find the null eigenvector is set to so that the matrix to be inverted is not singular. To find the other eigenvalues and eigenvectors we must start from an initial guess for the eigenvalue to be found. This guess can be updated every few iterations until the solution is found. The efficiency of the algorithm depends on how close is to the exact eigenvalue in comparison to the next closest eigenvalue. If the eigenvalue is known to high accuracy then convergence is found in a single iterationWilkinson (1965).

In order to use the largest possible value of we make two approximations. These approximations are based on the knowledge that if certain elements of the density matrix are known to be zero in the steady state then they can be omitted from the calculation by removing the appropriate rows and columns of the Liouvillian. We note that these terms must further be very rapidly decaying for the current noise not to be affected by their omission.

The first approximation we make is to neglect the density matrix elements corresponding to the , , and operators. This is valid since there is no coherence between the state and the or states so these elements must decay to zero in the steady state. This approximation reduces the dimensions of the Liouvillian to .

The second approximation is made in terms of the oscillator basis. The coherence between resonator Fock states which are widely separated in energy is small. Elements of the oscillator density matrix far from the diagonal can therefore be neglected. To check that this is a valid approximation the largest value on the last included diagonal of the resonator density matrix is found. So long as this value is below the results are found to be indistinguishable from the exact solution. After making these approximations the problem can be solved for using an inverse iteration method on a desktop computer. The limiting factor is the memory required to calculate in the iteration procedure.

Appendix B Mean Field Equations

This appendix describes in detail the mean field equation approach used in Sec. IV to solve the SSET-resonator system in the weak coupling limit. The mean field equations up to first order in the system operators were calculated in Ref. Rodrigues et al., 2007b. Here we go further and work to second order in the first instance,

(34)
(35)
(36)
(37)
(38)
(39)
(40)
(41)
(42)
(43)
(44)
(45)
(46)
(47)
(48)
(49)
(50)

where here the averages imply a trace over the SSET and resonator weighted by the density operator.

The first thing to note is that Eqs. (34) and (35) can be used to obtain a simple approximation for the average resonator position. In the steady state we findRodrigues et al. (2007b), , but from Eq. (37) we see that . Using the fact that the average current is related to the average charge on the SSET island, we then readily find that .

In order to obtain a closed set of mean field equations at second order we need to make a semiclassical approximation to eliminate the third order terms (e.g. ). This is done by setting the third order cumulant Gardiner (2004) to zero. In this context we apply the semiclassical approximation to products of three operators , and by assuming that

provided , and all commute. Where the operators involved do not commute the expectation value should be symmetrized appropriately in order for the commutation relations to be preserved.

Applying the semiclassical approximation to the term , in Eq. (44), we make the replacement and similarly for in Eq. (45). The resulting approximate equations are given by:

(51)
(52)

These equations can be linearized by treating and as parameters. To a good approximation the average resonator position is given by Eq. (19) in the steady state (with the current evaluated in the limit) and as discussed in Sec. IV, within the weak coupling regime it is also possible to use Eq. (20) to obtain .

The other third order terms we need to consider are and , which arise in Eqs. (49) and (50) respectively. Since and do not commute we must first rewrite the expectation values in the following way before expansion so that the commutation relations are obeyed.

(53)

Performing the expansion as before we can make the replacement,

(54)

These equations are readily linearized by treating as a parameter and using the fact that when the resonator is in a thermal steady state. The same procedure can be followed for the term to give

(55)
(56)

which completes the set of linearized equations. In order to solve this set of equations, we write the moments as a vector so that the equations of motion can be rewritten in the form where is the matrix of coefficients. The steady state for the moments is then obtained from the null eigenvector of .

This method is readily extended to obtain the third order mean field equations by instead setting the fourth order cumulant to zero and following the same procedure.

Appendix C Current Noise Calculation in the Mean Field

This appendix describes the calculation of the current noise within the mean field model using the counting variable approach Aguado and Brandes (2004b); Armour (2004); Flindt et al. (2004). We carry this out in terms of the second order equations, but the method is easily extended to higher orders. First we write the master equation in a modified form, where the number of quasiparticles, , to have tunneled across the right-hand junction (since ) are included,

(57)

Mean field equations can be calculated for this master equation in the same manner as before, where we now include a subscript to indicate the number of quasi-particles to have tunneled. The majority of the equations are the same but with a subscript , with averages now defined by , with the trace taken over the system operators but not . The following equations have a modified form,

(58)