Spin-detection in a quantum electromechanical shuttle system

We study the electrical transport of a harmonically-bound, single-molecule endohedral fullerene shuttle operating in the Coulomb blockade regime, i.e. single electron shuttling. In particular we examine the dependance of the tunnel current on an ultra-small stationary force exerted on the shuttle. We derive a quantum master equation for the full shuttle system which includes the metallic contacts, the spatially dependent tunnel couplings to the shuttle, the electronic and motional degrees of freedom of the shuttle itself and a coupling of the shuttle's motion to a phonon bath. We analyse the resulting quantum master equation and find that, due to the exponential dependance of the tunnel probability on the shuttle-contact separation, the current is highly sensitive to very small forces. In particular we predict that the spin state of the endohedral electrons of the endohedral fullerene in a large magnetic gradient field can be distinguished from the resulting current signals within a few tens of nanoseconds. This effect could prove useful for the detection of the endohedral spin-state of individual paramagnetic molecules such as N@C60 and P@C60, or the detection of very small static forces acting on a molecular shuttle.


I. INTRODUCTION
The determination of the spin-state of a single electron spin has attracted much attention recently [1]. Following the ground breaking experiments of Park et al., in the observation of quantized motional transitions in the conductance of a single-molecule C 60 transistor [2], interest in quantum electromechanical systems, or QEMS, has blossomed [3]. In this work we wish to examine whether a QEMS, with the inclusion of spin degrees of freedom, or a Spin-QEMS, can be used to detect the spin-state of the included spin. We have in mind the particular case of Nitrogen or Phosphorous endohedrally doped C 60 , N@C 60 or P@C 60 . This material has previously been shown to possess unpaired endohedrally trapped electrons, in a quartet ground state, S = 3/2, and which display exceedingly long transverse relaxation (or T 2 ), times. This material has been considered for use in quantum information processing devices but may have uses in other areas related to classical spintronic memories.
As a first exploration one can consider the original experiment of Park et al. using the paramagnetic endohedral, like P@C 60 , instead of C 60 . It is known that the electronic spin degrees of freedom reside in the endohedral electrons which are concentrated isotropically at the centre of the C 60 cage in P@C 60 [4]. In the experiment of [2], the C 60 oscillates in a Harmonic potential at THz frequencies. For the Spin-QEMS system of a P@C 60 molecule as a movable island in a single molecule transistor, one must be capable of engineering a coupling between the state of the endohedral spin and the conductance properties of the transistor. This can be accomplished through the imposition of a large magnetic field gradient ∂B/∂x, along the plane of the transistor. The interaction between the endohedral spin and this gradient field will impart a small force on the molecule which will in turn, cause a slight shift in the equilibrium position within the damped Harmonic oscillator trapping potential, shifting the molecule towards one of the contacts. This small shift will lead to a change in conductance properties. It will be shown that if the configuration is asymmetric, i.e. the mobile island is closer to one lead than the other, the resulting current depends exponentially, on this small spatial displacement. In the case of a tight trapping, observed in Park et al., [2], the resulting spin-dependent change in conductance is extremely small and unobservable. If we instead configure for the largest achievable magnetic field gradient (obtained for instance in MRFM devices), by placing a nanoscopic permanent magnet near to the mobile island, we also find that due to the large difference between the energy scales of the endohedral electronic spin flip (GHz), and the island's harmonic motion in [2], (THz), the island's spin-dependent position shift is minute and is unobservable in the current.
Guided by these observations we instead consider the more "floppy" Spin-QEMS of a shuttle device where the restoring oscillator now operates in the GHz frequency regime (see Figure 1). MRFM cantilevers fabricated out of single Si crystals have recently been generated which have resonant frequencies in the GHz range [5]. The resulting device could be similar to that outlined by Isacsson [6]. It may also be possible to engineer non-conducting flexible linker molecules which will connect the shuttle molecule to the metallic contacts and have resonant frequencies in the GHz range.
We thus consider the system of a movable shuttle or island, which is harmonically bound to oscillate between two metallic contacts in the presence of significant motional damping. We further assume that we operate in the Coulomb blockade region where only one electron can reside on the shuttle at one time. For the particular case of an endohedral N{P}@C 60 island, as in [7], we make the assumption that this injected electron distributes itself isotropically around the surface of the C 60 cage thus averaging out any magnetic dipole coupling between the spin state of the injected electron and the spin state of the endohedral electrons. We also assume that the injected electrons are not spin polarised and the resulting force experienced by the shuttle due to the interaction between the spin of the injected electrons and the magnetic gradient will average out while the force on the shuttle due to the endohedral electrons interacting with the magnetic field remains static. We therefore neglect the random spin-force contribution of the injected electrons in what follows.
Before going into the model in detail, we can understand, in a simple but surprisingly robust manner the underlying reasons why we expect to achieve an exponential dependence of the current on the island's displacement from its equilibrium position. The current through a simple source-fixed island-drain, two-stage sequential tunneling process [8], in terms of the left and right tunnel rates γ L , γ R , is given by If we now take into account the exponential dependence of these tunnel rates on the position of the (now mobile) island, and letting δ (positive) be the small displacements of the island from its equilibrium position to the right (to the drain), we can set and thus the right tunnel rate, from the island to the drain, increases due to the smaller separation while the left tunnel rate decreases due to the larger separation. This is reversed if the motion is to the left (δ negative), or towards the source: Inserting these into (1), we easily get: where F = γ 0 R /γ 0 L . If F = 1 then Ω = 1, and there is no difference between the currents in the two cases. When F ≫ e +2λ|δ| , Ω ∼ e +2λ|δ| , while if F ≪ e −2λ|δ| , Ω ∼ e −2λ|δ| , and in either of these asymmetric cases the two currents can be very different from each other if λ|δ| is appreciable.
In the sections below, we analyse the classical and quantum models of this system. Following the above simple rationale, we can find a realistic parameter range where the right-hand tunnel current shows a large (exponential), dependance on small spin-dependent island position displacements. This is confirmed below using low-resolution numerical simulations of both the "semiclassical" and quantum systems, and through high-resolution quantum simulations. Besides demonstrating the two different current "signals" corresponding to the spin force, we must also determine the noise present within the system to estimate the signal to noise ratio or so-called Fano factor. Following standard methods we derive the current spectral noise density and from this the signal to noise properties of the coupled system. Using all this we finally deduce the measurement time required to distinguish the spin state from the current signals and can find device parameters where this acquisition time will be several tens of nanoseconds.

II. MAGNETIC FIELD GRADIENT
We now estimate the magnitude of the maximum magnetic field gradient we can achieve at the location of the shuttle. We follow a similar derivation in Magnetic Resonance Force Microscopy (MRFM) [9]. The Hamiltonian for the shuttle's endohedral spin, S, interacting with a magnetic field which varies along the x−axis (see Figure 2), B(x), is H = − µ · B = −µ B g e B(x)S z , while the resulting force on the shuttle is given by F spin = µ B g e S z dB(x)/dx. If we consider placing a nanoscopic ferromagnetic spherical particle of radius R a distance d from the shuttle, one has where µ 0 = 4π × 10 −7 H/m, m is the magnetic moment of the ferromagnetic particle, r is the vector connecting the particle to the shuttle and r = | r| = R + d. For a sphere the magnitude of the magnetic moment is | m| = 4/3πR 3 M , where M is the magnetization of the ferromagnetic material. Following this we have B F erro = −µ 0 R 3 M/r 3ẑ . As the shuttle moves along the x−axis, the separations d and r, between the shuttle and the magnet change. Setting the distance between the shuttle and the magnet's centre to be r ≡ R + d 0 + x =x + x, where x is the small oscillation of the shuttle, d 0 is the equilibrium separation, andx = R + d 0 , one can derive For an iron ferromagnet, µ 0 M = 2.2Tesla, and choosing d 0 = 5nm, R = 5nm,x = 10nm, S z = ±3/2 (for P@C 60 ), one obtains F spin ∼ ±1 × 10 −15 N. Using this we can re-examine the spin Hamiltonian and expand to first order in x about the equilibrium shuttle position, where H 0 is a constant which we drop. From (6), above we set whereχ The Hamiltonian for the Coloumb-energy of the charged island/shuttle can be written as where c represents the operator which annihilates an electron on the shuttle, and c † c is the number of injected electrons on the shuttle, which in the blockade situation is either 0 or 1.

III. ELECTROSTATIC FORCE
In most studies of electron shuttles [3], one considers the shuttle to be driven primarily by the electrostatic force exerted on the shuttle through the interaction of the shuttle charge and an electric field present between the source and drain metallic contacts created by a potential difference V = V L − V R , across the inter-contact gap D. This gives an electric field of strength E = V /D.
In a standard treatment of electron shuttling [10], the shuttle acquires the electron at the average displacement of zero slightly towards the source and continues to move towards the source. It then undergoes electrostatic force to be propelled back towards the drain and loses the electron at a slightly displaced position towards the drain while continuing to move towards it. Again under electrostatic force the shuttle moves back towards the source and repeats the process. In this case no harmonic restoring force is required to sustain the shuttle instability. However in our case, in the situation where our shuttle is operating in the Coulomb blockade regime, the shuttle can acquire a single excess negative charge near the source and thus be forced electrostatically to the drain, but once the shuttle dumps its electron at a displaced position near the drain there is no electrostatic force which can bring this shuttle back to the source contact to repeat the shuttle sequence. Thus in our case a Harmonic restoring force is required to sustain the shuttling action.
In more detail,

IV. HARMONIC OSCILLATOR
The Hamiltonian of the Harmonic oscillator is dealt with in the standard manner. Setting the annihilation operator for motional quanta to be, √ 2a = x/x 0 + ip/p 0 , where x 0 = /(mω 0 ) and p 0 = √ mω 0 are the extents in x and p, of the ground state wave function of the harmonic oscillator, we have

V. COMPLETE HAMILTONIAN
We are now in a position to write out the complete Hamiltonian for the drain, source and motional, spin and electronic degrees of freedom of the shuttle: where χ ≡χx 0 , and η ≡ηx 0 . In the above: (12) is the Coloumb self-energy of the shuttle, (13) is the spin-dependent potential energy, (14) is the electrostatic potential energy, (15) is the Harmonic oscillator self-energy, (16) and (17) are the self-energies of the Fermi baths in the source and drain contacts, (18) is the self-energy and coupling of a phonon bath to the motion of the shuttle, (19) is the tunneling interaction between the source contact and the shuttle, while (20) is the tunneling interaction between the shuttle and the drain contact. We note the operatorsÊ L/R (x), in the tunneling terms. The amplitude for tunneling depends exponentially on the spatial separation between the source/drain and the shuttle. One can assume the formŝ where x is the displacement of the shuttle from its equilibrium position. In the following we assume the source and drain are identical in shape and material and thusλ L =λ R =λ, where for Gold contacts,λ −1 ∼ 3 Angstroms. As above we set λ ≡λx 0 to getÊ VI. DISPLACEMENT PICTURE As we mentioned above, the small static spin force exerted on the shuttle, in the absence of the electrostatic force, will cause a slight shift in the shuttle's equilibrium position within the Harmonic trapping potential. We now move to a frame of reference where we shift this displacement back to the origin. This can be done by considering the transformation where D(α) ≡ exp(αa † − α * a), is the standard displacement operator for operators satisfying [a † , a] = 1. By choosing α = −(χS z )/( √ 2 ω 0 ), which corresponds to the physical displacement of the shuttle's equilibrium position due to the spin force F spin , ofδ we can apply the displacement transformation (23), to the above full Hamiltonian to get where nowÊ where δ ≡δ/x 0 .

VII. MASTER EQUATION
Assuming that the time scales for the couplings to the source, drain and phonon baths are faster than those associated with the Harmonic oscillator and electrostatic potential, one can derive the following master equation [10]: where the super-operator D is defined to be with a dependence on the bias voltage (through the chemical potential). Notice the dependency of the Fermi factor on the spin force. If the displacement factor δ is large enough, the Fermi energies corresponding to the two spin directions will become significantly different and large enough to be observable at low temperatures. If one supposes that the bias voltage can be tuned to be between these two energies such that only one Fermi factor (corresponding to spin up say), is near unity while the other is nearly vanishing, it should be possible to differentiate the two spins based on the detection or absence of current. However with such low bias, we believe the noise in the system will be overwhelming and we thus proceed by assuming that the system is in large bias such that this effect is negligible. Now at milli-Kelvin temperatures, k B T ∼ 0.2µeV, while ω I ∼ 5meV. So if V > 10meV then f L → 1, while f R → 0, giving the simpler form for the Master equation:

VIII. EXPECTATION VALUES
From the master equation (41), we can derive the following differential equations: These equations can be further tidied by going to natural time units associated with the Harmonic oscillator, i.e. τ = ω 0 t, and by definingγ we get The above equations cannot be solved as they are not closed. One can break the correlations in the semiclassical regime of high damping. When this is done we finally obtain: while the complete master equation in these rescaled parameters takes the form, where In what follows we take the phonon bath to be at zero temperatureN = 0, and that the typical phonon bath excitation will have energies much smaller than that of the oscillator, i.e. below 1GHz . We will, in particular, be interested in the "semiclassical" and fully quantum steady-state behaviour of the right tunnel current, or the expected value of the operator: wherex ≡ x/x 0 ,p ≡ p/p 0 , andī R , is the re-scaled average current flowing from the moving island to the drain. We wish to determine whether ī R (t → ∞) = ī ∞ R , is significantly different in the two cases when the endohedral spin is up or down. However it is not enough that the two currents (endohedral spin up/down), are different but we must also then determine the steady-state (or DC or zero-frequency), quantum noise associated with these two steady-state currents. This DC noise will allow us to determine how distinguishable the two current signals are. This then will determine the minimum length of signal acquisition time one can average over to achieve a signal to noise ratio which will allow the endohedral spin states to be distinguished to a high level of confidence. In the following subsections (IX B) we discuss the initial "semi-classical" and quantum mechanical steady-state for the coupled electronic-vibronic shuttle system, paying particular regard to the change in right-tunnel-current with endohedral spin-state. This is done in a parameter range of the system where we are certain that our numerics are very precise and faithfully simulate the system. We find a parameter regime where ī ↑ R / ī ↓ R can be substantial, indicating that the endohedral spin-state measurement could be possible. In subsection (IX C), we expand the size of the numerical simulation to cover a much larger range of parameters (see parameter set (B) in Table I), which models a much larger electrostatic drive forceχ ES . This is done using the method of quantum trajectories [11,12,13]. This type of simulation is very costly and only a small number of such simulations were performed. We found that when the electrostatic driving force is increased, the ratio ī ↑ R / ī ↓ R remains significant. Further the highly accurate quantum trajectory simulations corresponds very well with the lower-resolution steady-state results indicating that our steady-state numerics could be trusted to quite large values of the electrostatic force. In subsection (IX D), we derive the analytical expressions for the noise spectra and numerically derive the DC-noise component over a limited (but precisely modeled), range of parameters. We reproduced the expected values for the DC noise of a standard, spatially fixed quantum dot, in the case when λ = 0, i.e. no coupling between the electronic and vibration degrees of freedom. We further found a slight decrease in this noise with the amount of electrostatic driving and degree of the exponential coupling. Finally, in subsection (IX E), we gather together all this information to estimate the measurement time required to distinguish between the current signals due to endohedral spin up(down).

A. Physical Parameters
We now determine the typical values of the dimensionless parameters appearing in the system's dynamics. We will take, m C60 = 1.2 × 10 −24 kg, ω 0 = 2π × 1GHz, the potential difference between the contacts to be V = 10meV, and the spatial separation between the contacts D = 10nm. From this we have x 0 ∼ 1.2Å, δ ∼ 0.18, giving the physical displacementδ ∼ .2Å. The dimensionless displacement of the harmonically bound shuttle due to the electrostatic force alone in units of x 0 is δ ES ∼ |F ES |/(mω 2 0 x 0 ) ∼ 28, while λ ∼ 0.4. From symmetry conditions one might suspect that we cannot gain any advantage from the exponential dependence of the tunneling amplitudes on the shuttle-contact separation if the shuttle's equilibrium position is exactly midway between the contacts (effectively whenγ L =γ R ). Thus we chooseγ L ∼ 0.1, whileγ R ∼ 0.9. We chooseκ ∼ 2, thus setting the system close to critical damping. We also setN = 0, or assume a zero temperature vibronic bath, throughout. The example re-scaled parameter values that we study below are summarized in Table I Thus the fixed point must satisfy: Following the derivation outlined in [10], the stability of the fixed point can be found by looking at the eigenvalues derived from a linearised dynamics. From [10], the system possesses one real eigenvalue and two complex conjugate eigenvalues which indicate the existence of a limit cycle. This limit cycle appears at a bifurcation of the fixed point which happens when the eigenvalues are purely imaginary. This condition appear at a critical value of χ ES : where A * = γ L e −2λx∓2λδ + γ R e 2λx±2λδ . Interestingly, the fixed point for the two spins are different and will result in a slightly different critical χ h values. However we choose to operate away from this regime in order to again avoid any possibility of noise upsetting the distinguishability of the resulting signals.
A more straight forward method is to keep system such that it falls within the fixed point regime and proceed in calculating the difference in the steady state current. This is possible by a careful choice of values such that the coupling χ ES are always below this critical values for both of the spins. This puts the system well in the supercritical regime. We can thus proceed by looking at the steady state at the fixed point and thus eliminating the possible error related to the existence of a stable limit cycle at the outer position of the doped fullerene's oscillations.
From (54)-(56), we see that when λ = 0, there is obviously no dependence of the current on the endohedral spin-state whenγ R =γ L . This is not the case when the coupling is non-vanishing i.e. λ > 0. This reflects the inherent asymmetry of the current transport from the source to the drain which is now mediated by an island which has, essentially, steadystate tunneling probabilities to either contact which are now effected by the endohedral spin dependent equilibrium position of the island. Taking, the steady-state semiclassical right tunnel current to be ī R

SC
in the limit of largeγ R /γ L (Γ R /Γ L ≫ e +4λδ ), and large λδ. Thus we might expect an exponential dependance of the resulting current ratio on the endohedral spin-state.
We can obtain the steady-state solution to (51), by converting the resulting Master equation into a 4N 2 × 4N 2 , (where N is the size of the truncation of the Fock-state representation of the vibronic Hilbert space), complex matrix and searching for the eigenstate (ρ ∞ ), with vanishing eigenvalue. We use the inverse power method to find the eigenvector with smallest eigenvalue [14]. In the following we mostly choose N = 36, giving matrices of size 4096 × 4096 ∼ 256MBytes in memory. As the memory resources increase with N 4 , we cannot realistically go beyond N ∼ 60. However, if the resulting large matrix is sparse then routines which only use matrix vector products can be used to determine the steady-state solution [15]. In Figure 3, we see that although the matrix has many small elements, it is unclear that it is highly sparse. As the current (53), is exponential inx, the expectation value ī R , could depend on very small matrix elements. To ensure convergence we will work with the full matrix in what follows and choose to truncate the bosonic vibrational Hilbert space at N = 16, or N = 32.
In Figure4, we plot the "semiclassical" evolution of (50), and quantum steady-state expectation values whenχ ES = 3,κ = 1,γ L = .1,γ R = .9, and λ = δ = .4. We see that the different endohedral spin states lead to very different phase space dynamics, differing island populations n , and currents ī R ∞ . One further notices a remarked difference between the magnitudes of the quantum and semiclassical expectation values for the currents in Fig. (4.c). This difference grows with the size ofχ ES , and is primarily due to the extended nature of the quantum state in the vibronic phase space. With the availability of the full quantum steady-state solution we are also able to visualise the quantum state via quasi-distribution functions and in Figure 5, we see that the quantum state is localised near the phase space origin. Moreover, from the compactness of the resulting Wigner or Husimi representations we can be relatively certain that the quantum steady-state solution is well represented with a Fock state truncation level set at N = 32. Asχ ES is increased towards the more physically realistic regime ofχ ES ∼ 20, one is justifiably concerned whether a truncation level of N ∼ 32, would remain sufficient to accurately represent the resulting quantum steadystate solution. We will see below however, that N ∼ 32, remains fairly accurate whenχ ES = 10. In the section below we are primarily interested in knowing with relative certainty that our quantum simulations are highly accurate. In Figure XI, we superimpose the semiclassical phase space evolution on to the Wiger function of the quantum steadystate solution for the two cases of endohedral spin up and down. We note that due to a slight assymetry in the Wigner function the quantum expectation value for the steady-state equilibrium solution is not located at the peak of the Wigner function.
We can see from this analysis that there is a significant difference between ī ↑ R , and ī ↓ R . To characterize this difference more systematically we numerically evaluate the ratio ī ↑ R / ī ↓ R , as a function ofγ L andγ R . We choosē χ ES = 3,κ = 2, λ = δ = .4, and range through values ofγ L andγ R , computing the resulting quantum steady-state current ratios for the two cases of endohedral spin up and spin down. The results are shown in Figure XI, and in the limit of largeγ R /γ L (Γ R /Γ L ≫ e +4λδ ), the current ratio appears to asymptote, as it should to exp(4λδ) ∼ 1.9. By repeating these simulations at the reduced truncation N = 16, and finding that the resulting current ratios are within 10 −9 of those computed at N = 32, indicates that our quantum steady-state simulations are extremely precise.

C. Quantum Trajectories
As mentioned above, going much beyond a Fock state truncation of N ∼ 60, with the full matrix representation of the Liouvillian is not possible without exotic computational resources. Another method of simulating the dynamical evolution of the full quantum master equation (51), is to use the method of quantum trajectories [11,12,13]. Briefly, in this method one represents the quantum master equation (51), as an average of a stochastic Schröedinger equation over some stochastic noise. The temporal evolution of an initial state determined by this stochastic Schröedinger equation is known as a quantum trajectory. One can determine the time evolution of quantum expectation values of an observable by taking the expectation values of that observable within a single quantum trajectory and then averaging the expectation values found over numerous stochastic repetitions of the quantum trajectory. Since the method primarily involves the Schröedinger evolution of pure quantum states, the quantum evolution of an individual quantum trajectory can be simulated with much higher Fock state truncations N , than in the above steady-state method. The down side is that to obtain precise expectation values one must average over a large number of repetitions of the quantum trajectory.
The evolution of the stochastic Schrodinger equation corresponding to (51), whenN = 0, is conditioned on the appearance of random quantum jumps. The information gained by the environment, in a sense, yields information which conditions the quantum state following the quantum jump. In (51), there are three types of quantum jump: (I) |ψ → c † E L |ψ (unnormalised), which corresponds to an electron tunneling onto the island from the Fermi bath of the left electrode with some associated monitoring (via E L ), of the island's position state by the vibrational bath, (II) |ψ → cE R |ψ , corresponding to an electron moving off the island into the Fermi bath of the right electrode with some associated monitoring (via E R ), of the island's position state by the vibrational bath, and (III) |ψ → a|ψ , where a motional quanta escape's from the island's position state into the vibrational bath. In Figure 8, we show an example of a single quantum trajectory of the coupled electronic-vibronic shuttle system corresponding to (51). We choose no motional damping and set most of the parameters to vanish so we can illustrate the effects of the quantum jumps clearly. At time points (A) in Figure 8, a single electron jumps onto the island. Since this is more likely for negativē x, due to the exponential dependance of E L on position, we see a jump in the island's position towards the source contact located in the −x region. At time points (B), this electron now jumps off the island into the drain. Since again this is more likely closer to the drain due to the form of E R , we see a quantum jump in the island's position to larger values of +x. While the island is either empty/occupied, the non-Hermitian evolution suffered by the quantum trajectory causes the system to move towards the situation where the occupied/empty state becomes more and more probable. As a consequence of the exponential dependances of E L and E R , this means that the oscillations grow towards larger +x, if the island is occupied, while they grow towards −x, if the island is empty. The net effect is to pump the vibrational motion of the shuttle during the single trajectory and even though this growth is interspersed with dips due to the movement on or off of an electron. Once averaged over many quantum trajectories, the resulting effect is to pump the average energy of the oscillation in an unbounded manner (if there is no motional damping).
We now consider the full-blown quantum trajectory simulation of the system using the parameter values χ ES = 10,κ = 2,γ L = .1,γ R = .8, and λ = δ = .4. We now have chosen a far larger value for the electrostatic force and are able to expand the Fock basis truncation to N = 100 using quantum trajectories. To determine the steady state quantum solution using a full matrix representation would now require over 26GBytes of memory! The resulting simulation averages over 2000 quantum trajectories and required many days to execute on a 2GHz Operton Dual processor 64-bit computer. In Figure 9, we plot the right tunnel current ī R , as computed via quantum trajectories for the two endohedral spin states. We also plot the "semiclassical" dynamical evolution and the steady-state quantum expectation values for the currents where the latter is computed at the necessarily lower Fock truncation of N = 40. We can clearly see that the (incredibly computationally expensive), quantum trajectory simulation remarkably agrees very well with the values obtain from the steady-state numerics. It is not possible to plot out the Wigner function from the quantum trajectory simulation as this method can only yield the values of a small number of expectation values without running into the difficulty of requiring huge amounts of storage. However, the results shown in Figure 9, strongly indicates that the technique of directly determining the quantum steady-state via the eigenvector of the Liouvillian matrix with vanishing eigenvalue, can extend to very large values of the electrostatic driving forceχ ES , without any significant lack of precision. As the physically relevant regime for χ ES ∼ 20, the behaviours we obtain here withχ ES ∼ 3 − 10, should therefore also hold in the physically relevant regime.

D. Estimating the Steady-State Quantum Current Noise
There are a number of methods to deduce the noise power spectral density S(ω). We will use the methods outlined in [16].
We first represent the instantaneous right tunnel current as i R (t) = edN/dt, where dN (t), is a classical point process which represents the number (either zero or one), of tunneling events seen in an infinitesmal time, dt, and e is the electric charge. We can consider dN (t) to be the increment in the number of electrons N (t), in the drain electrode in a time dt. The steady-state current is computed by , represents an average over the stochastic process describing the tunneling events. In our case this is given by, where the jump superoperator J [A]B ≡ ABA † . The fluctuations in the current are quantified by the two-time correlation function From the theory of open quantum systems [12], one can show that where J R = √ γ R c exp(−λx) (see (21)), and the Liouvilian evolution, exp(Lt), is according to the master equation (41). In the case of a decoupled system (λ = 0), this leads to a normally ordered correlation function while one obtains an antinormally ordered correlation function if one were studying the noise of the left tunnel current [17]. For the coupled system the interpretation of the operator ordering is more complicated.
The current noise spectral density is given by In terms of our rescaled variables and parameters, τ = ω 0 t, and (46), we setτ = ω 0t , to obtain, where now the Liouvillian evolution is given by (51). The spectral density is now given by From this we can derive the Fano factor, F = S(ω)/(2ei ∞ R ), to be The Fano factor F gives information on the statistics of the tunnel current noise. If F = 1 then the noise is completely uncorrelated/Poissonian, i.e. white noise. If F > 1, then the noise is known as super-Poissonian and the tunnel events are bunched. If F < 1, then the noise is known as sub-Poissonian and the tunneling events are anti-bunched, i.e. there is little chance that a second tunnel event will closely follow a previous event. Quantum correlations within the coupled system are primarily responsible for F = 0. In the case of the decoupled system, In Figure XI, we compare the results for the decoupled system (λ = 0), with the standard results for a two-state sequential tunneling process [8] In Figure 11, we display the current noise spectra for various frequencies, i.e.
and see that typically the noise is sub-Poissonian. In Figure 12, we contour plot the DC (or zero frequency), Fano factors for the decoupled and coupled system and observe that the coupled system displays stronger anti-bunching whenγ L ∼γ R , than in the decoupled case. However in the parameter region where we are interested,γ R ≫γ L , the DC Fano factors of the decoupled and coupled system are very similar.

E. Measurement Time
The final quantity we must determine, given the DC current signals and current noise associated with the two endohedral spin-states, is the measurement time τ m , required to differentiate the two currentsī ↑ ∞ R andī ↓ ∞ R . This is given by the expression [18], where S 1 , and S 2 , are the DC current noise spectral densities for the two cases of spin up/down while ∆I ≡ | ī ↑ ∞ R − ī ↓ ∞ R |. In Figure 13, we plot log 10 τ m , for variousγ L ,γ R . We see clearly that whenγ R ∼γ L , τ m must be very large to discriminate between the very small difference in steady-state currents. In the parameter regions we are concerned with, i.e.γ R /γ L ≫ 1, this measurement time is τ m ∼ 10 1 −10 2 , in the natural units of the oscillator (ω 0 = 2π ×1GHz), or a few tens of nanoseconds. Thus in this parameter region, the measurement time is very realistic.

X. SUMMARY
In summary, we have shown that the right-hand tunnel current of a strongly damped single electron shuttle operating in the Coulomb blockade regime has very high sensitivity to small equilibrium position displacements of the shuttle. We considered a N@C 60 molecular shuttle whose endohedral spin can exert a force on the molecule when placed in a large magnetic field gradient. By placing a large nanoscopic magnet nearby, a tiny ∼ ±10 −15 N force is suffered by the molecular island arising from the endohedral spin state S z = ±3/2. This small force alters the equilibrium position of the island in the presence of the electrostatic driving force, motional damping, and harmonic restoring force. Through the extremely sensitive (exponential), dependence of the conductance on the tunnel separations, the shift in the equilibrium position of the island profoundly influences the current flow. We determined the current noise spectral density and the measurement time required to distinguish between the two DC steady-state current signals (spin up and down), in the presence of steady-state quantum noise.
The results we have found would strongly indicate that there are realistic parameter regimes where the spin state dependent currents are distinguishable within several tens of nanoseconds. Thus this device should thus be capable of spin-detection. Alternatively, similar devices (without an endohedral spin), should be capable of discriminating small, ∼ ±10 −15 N , static forces on a C 60 island.