Electronic transport in a shuttle with a discrete energy spectrum

We present a multiscale approach for a molecular nanoelectromechanical shuttle, which takes into account the discrete energy levels and realistic orbitals of a specific molecular system. We deploy a Cu13 cluster as a central island and describe the electrodes within the jellium model. The shuttling features are then examined for a wide spectrum of bias voltages. The microscopic data of the cluster are obtained with density functional theory (DFT) while the macroscopic part combines stochastic charge dynamics, based on the microscopically evaluated tunnelling rates, with Newtonian dynamics. We find four transport regimes within the analysed bias voltage interval: a blocked regime (Vbias⩽ 2 V), a pure shuttling regime (3 V ⩽Vbias⩽ 4 V), a mixed regime (5 V ⩽Vbias⩽ 6 V) and a direct tunnelling regime (Vbias⩾7 V). In the two latter cases, the shuttling motion leads to a reduced current as compared to an island with a fixed position. The transport characteristics are strongly coupled to the availability of energetically allowed channels and the system therefore depends strongly on the microscopic details and bias voltages. For the setup at hand, however, the shuttling transport of one electron is ubiquitous and stable for all but the blocked regime.

. Setup. The electrodes consist of two semi-infinite jellium slabs (r s = 3a 0 ). The electronic characteristics of the mobile cluster (Cu 13 ) are obtained with DFT. The cluster can move in the direction perpendicular to the electrode surfaces, and the bias voltage is chosen so that the left and right electrodes act as the source and drain, respectively. Note that the geometric surfaces z b and the effective surfaces z 0 are non-equivalent. The former represents the edges of the positive charge background whereas the latter is seen as the effective locations of the metal surfaces. All parameters used in the computations are calculated for the effective surface. The variables L and Z c are the effective gap and the cluster centre-of-mass relative to the (left) electrode surface, respectively. The distance between the geometrical and effective surfaces is determined to be z = 0.34 Å.
several transport regimes and behaviours. Our results offer insight into the complexity of the molecular shuttle.
We focus on the setup that is most favourable for a molecular shuttle, and consider only electrode-shuttle pairs that do not exhibit significant chemisorption. In this weak coupling approximation, the electronic structures of the shuttling cluster and the metal electrodes can be regarded as independent of each other and can be evaluated separately. The electronic structure of the shuttling cluster and its effective potential is computed with density functional theory (DFT), while transition probabilities are calculated with the transfer Hamiltonian method. Finally, coupling between the mechanical and charge degrees of freedom are handled by a dynamic Monte Carlo scheme. For the surface-cluster interactions, we use a Born-Mayer shortrange repulsion and a classical image charge potential modified for close-range interactions.
The text is organized in three sections. In section 2, we present the model we use for the molecular shuttle. In section 3, we present our main findings for the physical observables, in particular discussing the I-V characteristics and the properties of the shuttling motion. Finally, we summarize our observations in section 4.

Method
The shuttle system in this work consists of two electrodes and a mechanically active part as depicted in figure 1. Note that in contrast to e.g. [6], a mechanical restoring force is not 4 included. The electrodes are modelled as two semi-infinite (in the z-direction) jellium slabs with a 20 Å gap between their geometrical surfaces. The mobile part, a Cu 13 cluster, is situated in the electrode gap and is free to move in the z-direction perpendicular to the electrode surfaces. In all instances, the intra-cluster atomic positions correspond to those in a relaxed neutral cluster, which has a highly symmetric icosahedral structure. The orientation of the cluster is chosen so that its highest symmetry axis coincides with the z-axis (see figure 1). This rigidity of the molecule is introduced in order to reduce computational complexity and increase overall method efficiency. However, it should be noted that non-centre-of-mass vibration-assisted electron transport has been shown to be of some importance for similar setups. A more quantitative approach should therefore also include intra-molecular degrees of freedom [27].
Our approach is based on well-established approximations such as the adiabatic approximation, weak coupling between the electrodes and the central island and the treatment of the cluster as a classical particle [7]. The weak coupling allows us to regard the electrons as being localized either on the leads or on the cluster. The electrons in the different subsystems (the cluster and the electrodes) are considered to be in equilibria and the transitions between leads and cluster are regarded as instantaneous. Key molecular features, such as the effective potential and electronic structure, are obtained numerically from the ground-state DFT calculations. (We use the plane-wave DFT code DaCapo with a PW91 exchange correlation functional.) For simplicity, we consider a system operating at zero temperature.
Previous studies have shown that, except for cluster positions very close to the leads, first-order perturbation theory adequately describes the influence of the electrodes on a qualitative level [26], supporting the assumption of weak coupling. Accordingly, the full system Kohn-Sham eigenvalues and wave functions can be written as and respectively, where v pert is a change in the one-electron potential for the cluster due to the vicinity of the metal electrode, and DFT j and ψ DFT j are the Kohn-Sham eigenvalues and orbitals for an isolated cluster. The subscript j denotes the orbital and Z c is the centre-of-mass position of the cluster. The electronic structure calculations for the cluster and electrodes can thus be separated and the number of DFT simulations is reduced to one for each integer charge state of the isolated cluster.
The perturbation potential comprises the bias voltage and Coulomb forces from the induced surface charges on the electrodes. Short-range Pauli repulsion, which is a many-body effect, is included in the forces but not in the one-electron perturbation potential. The implementation of van der Waals forces and chemisorption are outside the scope of this work. The former gives a small contribution compared to other system interactions and the latter violates the assumption of weak coupling.
The one-electron perturbation potential for z L 0 < z < z R 0 is given by The first term in (3) is the potential due to the first order of images of the atoms on both sides of the gap. The sum is over the number of atoms, N, andr L/R is the location of the image charge. Near the electrode surface the image potential deviates from its classical (divergent) form, and a damping factor (1 − exp[ − λ(z − z 0 )]) is introduced following [28]. The two parameters λ = 1.6 Å −1 and z 0 = 0.34 Å have been determined by fitting the effective potentials in [29] to the one-electron image potential suggested in [28]. The parameter z 0 can be associated with the effective surface of the jellium slab for the image charge and electrical field calculations [30]. The electron distribution around each nucleus, (4), is approximated by a Gaussian, with the width parameter β = 1/2σ 2 determined by fitting the electron density of DFT calculations of single Cu atoms to a Gaussian distribution. The effective nuclear charge as seen by the valence electrons is Q v = 11e, while Q e = Q v − Q/N is the valence electron charge on the copper atoms and Q is the net charge on the cluster. The first part of the second line in (3) is the image potential contribution from the higher order terms in the image charge expansion, i.e. images of images. To simplify the calculations, we approximate the higher order terms that arise from image charge quite far from the cluster by their values at the centre-of-mass of the cluster. This allows their closed form expression in terms of the digamma function ψ(z). The symbol γ denotes Euler's constant and L is the effective gap. The final term in (3) is due to the external DC bias voltage. Inside the electrodes, z < z L 0 and z > z R 0 , the perturbation potential is modified to balance the long-range Coulomb forces from the cluster, giving a total of zero potential (plus bias voltage). The eigenvalue corrections are calculated for a set of centre-of-mass cluster positions, while values for intermediate locations are determined by linear interpolation.
Data concerning the electrodes are extrapolated from a jellium model with Wigner-Seitz radius equal to 3a 0 . For parameters beyond the jellium model, such as lattice constants and electrode atom sizes, we use the properties of gold.
In the weak coupling limit, transition rates can be calculated with the transfer Hamiltonian method [31,32]. Since we consider only lowest order perturbative effects, the two electrode-cluster barriers are treated separately. For the left electrode, the Hamiltonian is where H g is the Hamiltonian for the free cluster and H e is the Hamiltonian for a free-standing electrode. The corresponding transition rate for one cluster orbital is In (6), the expectation value is taken over the effective gap. The occupation of orbital j is n j , while V g and V e are the potentials associated with H g and H e , respectively. The electronic electrode-vacuum barrier is V 0 . The electrode wave functions ψ k (k denotes the wave vector) are calculated using a semi-infinite square potential [33]. Due to the complexity of (6), it is computationally inefficient to repeat the calculation of j for every cluster position. Instead we use the phenomenological relation j = 0 exp[ − z/ ] [34]. The parameters and 0 are determined for each orbital j by computing (6) for two 6 different cluster positions. The total transition rate is found by summing over all energetically allowed orbitals for a position Z c . A number of low-energy unoccupied Kohn-Sham states are included in the summation. These are regarded as excited states of the system [35], and give rise to several channels for the tunnelling of charge onto the cluster. The mechanical path of the cluster is described with a dynamic Monte Carlo scheme where each simulation gives one possible trajectory {Q(t), r(t)}. The time evolution of the cluster charge Q(t) is determined by stochastic tunnelling events with tunnelling rates that depend on the position and charge of the cluster. The time evolution of the cluster position r(t) is determined by the forces acting on the cluster which include interactions with induced surface charges, short-range repulsion and a dissipation source. The first of these forces are where W is the electrostatic energy deduced from the image potential used in (3) and integrated over the charge distribution of (4). The electrostatic energy from the left lead is thus where the last term V h.o. (Z c ) is the electrostatic potential due to the higher order terms in the image charge expansion.
For the short-range repulsion we use a Born-Mayer potential. Integrating over the electrode surface and summing over the cluster atoms gives The above equation is derived by treating the surface as a continuum except for one atom coinciding with the axis of motion of atom i. We do this to maximize wall repulsion, consequently minimizing the coupling between the systems. The pair-potential-specific Born-Mayer parameters A r and b r are calculated consistently with [36] for the gold-copper interaction. Lattice spacing, a, is taken from a (100) fcc gold crystal surface. Dissipation is necessary for the stability of shuttling. For a many-atom specimen impinging on the surface, in the expected range of velocities, phononic emission dissipation dominates loss processes [37]. We employ a simple model of phonon emission, which captures many of the qualitative features of a real system. We model the electrode surface as a mass on a spring and couple it to the cluster motion following standard mechanical analysis (see figure 2). We use the equations where the first row of (9) describes the cluster motion and the second equation relates to the surface motion. The variable F diss corresponds to the first-order dissipative force on the cluster, The parameter Z is the distance between the surface and the cluster. of the mobile part of the surface, M in the second equation, is a function of the cluster-surface distance due to the finite range of the cluster-surface interactions, while k is the spring constant for the single gold atoms that make up the electrodes. The methods we have used to calculate these surface parameters are presented in [26] and will not be discussed here. In contrast to viscous damping, often used in more phenomenological studies, the damping in our model is explicitly attributed to surface-cluster interactions, and is most effective when the cluster is near the electrode surface. Note that the intra-cluster vibrational degrees of freedom not included in this work can play a significant part in energy transfer from the centre-of-mass translational motion.

Results
In accordance with Coulomb blockade arguments, the energies of specific orbitals grow linearly with the number of electrons on the system, in a similar fashion for all occupied orbitals. The relative distribution of the occupied eigenvalues within a charge state therefore changes only marginally. For the less bonded excited states, however, the width of the eigenvalue distribution decreases with an increasing number of electrons, as shown in figure 3. This generates an asymmetry between positively and negatively charged clusters. Another asymmetry in the charge spectra is due to the perturbative image potential in (3). Since the potential is repulsive for Q > 0 and attractive for Q < 0, the eigenvalues as functions of the cluster position behave differently for positive and negative charges. Finally, the energies of the Kohn-Sham orbitals are distributed in more or less dense bands, so that a small change in the external potential can cause the opening or closing of a large number of transport channels. This sensitivity together with the aforementioned charge asymmetries are the main contributors to the versatility of the system characteristics. These features are typically seen neither in the extreme quantum limit of a few electron levels in the mechanically active element, nor in the case of macroscopic metallic grains, but are a typical mesoscopic phenomenon.
The form of the wave functions also plays an important role in the system characteristics. Atypical orbitals that extend far from the atomic nuclei along the z-axis can give a significant contribution to the transition rates even when the cluster's centre-of-mass is far from the surface. The activation of such a level can change the system properties substantially in comparison to a somewhat smaller or larger bias voltage in which the specific level is energetically unavailable. Consequently, an interesting feature of the system is a significant probability for long-distance tunnelling. In the left half, specifically for Q = 1e (for larger biases also for Q = 0e), a large number of excited states fall below the chemical potential µ R of the right electrode, and are therefore available for electrons tunnelling into the cluster from the right. These orbitals are intrinsically quite extended and therefore give a considerable possibility of a charge transfer event when energetically allowed. A similar but weaker effect takes place in the right half, where for Q = −1e the occupied states are very high above µ L . The highest occupied orbital (HOMO) is considerably extended in the z-direction, resulting in a substantial transition probability also far from the left electrode. This is an example of how the detailed characteristics of the electronic structure can determine the traits of the entire system, showing the importance of atomistic modelling of each specific experimental system. The main transition rates and their dependence on charge and position are given in figure 4.
An important feature of the system is the very low probability for transitions to charge states Q = ±2e. These transitions are often energetically allowed, even for small voltages, but in fact, even for V bias = 7 V, there is at least an order of magnitude higher probability of transferring  to Q = 0e before reaching Q = −2e for all island locations so that the highly charged states are dynamically rather than energetically blocked. For lower bias voltages this difference is several orders of magnitude. Transitions to Q = 2e are seen first at V bias = 8 V. Due to the much higher probability of reaching charge neutrality, the number of higher charge events are so few as to be qualitatively negligible for the considered range of bias voltages. For an interval close to the right lead, the cluster acts as a bridge for electron tunnelling between the leads (Q = −1e ↔ Q = 0e). Nevertheless, the shuttling motion is very stable and the shuttling transport from the previous regime is continuously upheld. The contribution from the bridging events that involve tunnelling over long distances increases current fluctuations. (d) The direct tunnelling regime, V bias = 7 V. A second interval where the cluster acts like an effective bridge opens up in the left part of the gap. Higher charge states (Q = −2e) are possible but uncommon and of little relevance for the transport properties. The motion is still stable, except for small, localized variations not seen in the previous regimes. The ubiquitous shuttle transport is present but contributes negligibly to the total current. Instead, the number of charges transferred per revolution varies significantly. Interestingly, in the mixed and direct tunnelling regimes, the shuttling motion reduces the total current through the system as compared with a fixed cluster located in one of the bridging positions. The force on the central island shows an image charge potential well followed by a sharp repulsion as the cluster closes in on the surface (see figure 5). The attractive well is somewhat sharper for positive island charges. This is partly due to their more compact form deferring the effects of the repulsion, but there is also a contribution from the dispersion of electrons as compared to the point-like localization of the nuclei, decreasing the interaction between the surface and the negative charge. As a consequence, a more shallow potential well is formed for negative Q while the positive charge well deepens. However, these effects are more pronounced only for island charges |Q| 2e, and are therefore of little concern for the Cu 13 setup.
For the chosen nanogap geometry and the investigated range of bias voltages, we identify four separate transport regimes. The differences in paths and transfer characteristics are presented in figure 6. In the first regime, 1 V V bias 2 V, the cluster moves back and forth in the potential well, 'bouncing off' the short-range repulsion close to the leads. The current is predominantly zero, but sporadic shuttling events, in the case of V bias = 2 V, result in an occasional transport of charge through the system. For all higher bias voltages, one positive charge (Q = 0e ↔ Q = 1e) is shuttled for every revolution. The second regime, 3 V V bias 4 V, represents an (almost) pure shuttling regime with minimal current fluctuation. In this regime, the mesoscopic shuttle behaves very much like the phenomenological models predict [7]. In the third regime, 5 V V bias 6 V, a large part of the electronic transport involves a negatively charged state of the cluster, Q = 0e ↔ Q = −1e. For these events, both charging and de-charging takes place in the right half of the gap. Rather than carrying charge mechanically across the gap, the cluster divides the gap between the electrodes into two shorter segments, effectively operating as a bridge. The number of tunnelling events varies significantly from revolution to revolution, introducing a substantial amount of noise into the system. It is worth noticing that the noise mainly affects the current, while the fluctuations in the mechanical motion are extremely small. For the fourth regime, V bias 7 V, a similar bridging process appears also in the left half of the gap. For this regime, non-shuttling transitions dominate transport. Fluctuations in current are large while the mechanical motion start to show some, if still small, variations. This kind of direct tunnelling is a feature of the nanoscopic gap size we use in our setup. Even though larger gaps can be used in theoretical studies to focus on the shuttling process, nanoscale gaps are closer to the experimental reality of molecular electronics [38]. The I-V characteristics and number of electrons transferred per cycle are given in figures 7 and 8.
Except for small variations in the highest and lowest regime, the mechanical motion is impressively stable. For 3 V V bias 7 V the median amplitudes are nearly identical, showing only a minuscule increase with bias voltage (A ≈ 11.4 Å with an offset of 0.1 Å towards the right lead). This implies that the location of the turning point is more determined by the sharpness of the short-range repulsion than an effect of the bias voltage acting on the system charge. The oscillation period, instead, varies noticeably also within the Monte Carlo paths (T ∼ 4 ps): the current increases due to increasing frequency of the mechanical motion rather than due to larger number of charges transported per cycle. Even though the overall trend shows a clear decrease in period as the bias voltage increases, the period is very sensitive to the amount of time the island stays charged during a revolution. Since the latter time scale is strongly dependent on the exact features of the transition rates, it is difficult to make a quantitative prediction about the period from one bias voltage to another.
While the ubiquitous stability of motion also has a stabilizing effect on the current, a minimum of charge transport is needed in order to balance dissipation and stabilize the motion, as can be seen in figure 6(a).

Conclusions
Many of the predictions of more phenomenological studies have proved valid also for the mesoscopic system at hand, in particular in the pure shuttling regime. Central features such as the limitation of charge on the system and the charge dependence of the eigenvalues seem to be comprehensively explained by classical Coulomb blockade arguments, although we find that sometimes states that would be energetically allowed are dynamically blocked by more probable transitions. We find that the self-sustaining shuttling motion is ubiquitous for all but the lowest voltages. From V bias = 3 V to 6 V, a stable shuttling regime develops, where electron shuttling plays a more (low biases) or less (high biases) important role in the transport of charges through the system. However, for the highest voltages, the cluster works more effectively as a bridge than a shuttle; hence, in this regime, the shuttle motion reduces the current that is dominated by bridging events, cf [39]. In experimental implementations, further complications may arise through effects that are beyond the weak coupling approximation such as chemisorption, or more complicated cluster-surface interactions. Such effects are likely to be quite geometry and materials specific, and were not the focus of our study.
The key mechanism of the electronic transport is the energetic availability of transfer channels as a function of position and charge. This dependence renders the system extremely sensitive to small changes in the energetics-potential landscape, applied voltages and electronic spectrum. Together with the complexity of molecular spectra and surface interactions, this sensitivity implies that microscopic energetics must be carefully considered when analysing charge transport in a molecular shuttle system. The details of the electronic structure manifest themselves in, for instance, an asymmetric behaviour of the shuttle system, both in charge and position. Experimental observation of shuttling in this simplest implementation is therefore quite challenging as the shuttling phenomena are sensitive to the details of the system and the shuttle current may be dominated by other current carrying mechanisms. The same technique applied to a single Cu-atom also yields a stable shuttling regime for small voltages above the critical voltage and a noisy regime for higher bias voltages. The atomic high-noise regime, however, is due to other phenomena than the direct tunnelling regimes of the cluster shuttle [40]. This again implies a strong dependence of transport on setup details.
The mechanical motion, on the other hand, shows surprising stability, with almost identical amplitude for all voltages. Specifically the latter is due to the sharpness of the short-range repulsion potential, determining the turning points almost independently of the kinetic energy of the shuttle. Again, a precise treatment of cluster-electrode interactions is fundamental for a precise description of the shuttle. Since the regular mechanical motion persists even when the total current is dominated by other mechanisms, it may be advisable to search for shuttling phenomena using probes that are more directly sensitive to the mechanical motion rather than the overall current.