Multiscale modelling of a nanoelectromechanical shuttle

In this paper, we report a theoretical analysis of a nanoelectromechanical shuttle based on a multiscale model that combines microscopic electronic structure data with macroscopic dynamics. The microscopic part utilizes a (static) density functional description to obtain the energy levels and orbitals of the shuttling particle together with the forces acting on the particle. The macroscopic part combines stochastic charge dynamics that incorporates the microscopically evaluated tunnelling rates with a Newtonian dynamics. We have applied the multiscale model to describe the shuttling of a single copper atom between two gold-like jellium electrodes. We find that energy spectrum and particle surface interaction greatly influence shuttling dynamics; in the specific example that we studied the shuttling is found to involve only charge states Q = 0 and Q = +e. The system is found to exhibit two quasi-stable shuttling modes, a fundamental one and an excited one with a larger amplitude of mechanical motion, with random transitions between them.


Introduction
Nanoelectromechanical systems that combine electrical and mechanical functionalities on the nanometre scale have in recent years attracted a great deal of theoretical and experimental interest [1,2]. The nanoelectromechanical shuttle is a structure that resembles a single electron transistor but incorporates mechanical motion of the central island. Previous theoretical studies on the shuttle have shown that in the presence of a dc applied bias the charge and velocity of the central island are correlated, Q(t)Ż(t) = 0, which implies that the shuttle absorbs energy from the dc field and converts it into mechanical motion. The shuttle motion facilitates charge transfer through the system, and signatures of mechanical motion can be seen both in the current-voltage characteristics and in the noise properties of the device [3]- [5].
Several theoretical studies have been carried out for different set-ups of the shuttle since the first description of this phenomenon. The theoretical studies cover different size regimes of the shuttle, featuring coherent [6] or sequential [7,8] tunnelling and quantum mechanically [9,10] or classically [11,12] described mechanical motion. The studies have shown that the shuttle instability strongly depends on the bias voltage and the system set-up. This sensitivity also renders the shuttle behaviour dependent on the precise description of the problem.
Experimental evidence of coupling between vibrational degrees of freedom and electron transfer has been found for both microscopic [13,14] and macroscopic [15] systems. In particular, the experiment by Park et al [16], using a C 60 molecule between gold electrodes, has demonstrated the type of coupling that has been considered by many theoretical studies and has increased the interest for a molecular shuttle [6,9,10,12,17].
In the shuttle geometry the mechanical motion is on a nearly macroscopic time scale, typically from picoseconds for small molecules to nanoseconds for large molecules such as carbon nanotubes. The motion is excited due to tunnelling events between the mobile object and the stationary electrodes, which have a typical time scale of femtoseconds that is determined by the electronic structures of the mobile molecule and the electrodes. Hence, a theoretical description of the shuttle system naturally calls for multiscale methods that combine the fast electronic time scales with the slower mechanical ones. Thus far research has concentrated on the slow degrees of freedom (long length scales) while dealing with the fast ones (short length scales) in a phenomenological approximation.
In this study, we examine some of the issues that are important for a molecular shuttle system at zero temperature. We focus on the thus far unexplored effects that emerge in the shuttling of quantum mechanical objects with realistic electronic structure comprising many orbitals with different symmetries. A key feature of the system is the interaction between the shuttling object and the stationary electrodes which includes a purely energetic part associated with the potential energy surfaces of the surface-shuttle system, and a charge transfer part describing electron transfer between the different subsystems (shuttling object and electrodes). The latter issue is problematic from an electron structure point of view since the charge equilibration rate strongly depends on the separation between the shuttling particle and the electrode, and at large separations the subsystems are independent of each other as far as their electronic structures are concerned while at near proximity the electronic structures must be treated jointly. In this study, we concentrate on the situation when the hybridization between the shuttling object and the surfaces is small and charge transfer can be described by stochastic tunnelling events between unhybridized electronic states-systems in which chemisorption plays an important role are hence excluded from this study. Since the impact of the small shuttling object on the electronic structure of the electrodes is quite small in the absence of chemisorption, we have opted for modelling the electrodes as two semi-infinite jellium slabs. This allows us to restricted the numerical electronic structure calculations to the shuttling object alone and incorporate the effects of the electrodes as external surface potentials. This decoupling also eliminates spurious effects such as hybridization between two degenerate orbitals that are separated by a large distance.
A priori, it is unclear if the predominantly attractive surface forces are so strong that shuttling of a nanoscale object in general is prohibited, or if the adhesive forces can be overridden by reasonable electric fields.Addressing this issue is one of the main motivations of the present study. The surface forces are included as external potentials in the static density functional description of the shuttling object [18,19], which provides information on the energy spectrum of the island as well as structure of the relevant orbitals. The atomistic description of the central island is hence coupled to an effective description of the electrode surfaces. The mechanical motion of the shuttle is described using classical dynamics that allows us to treat macroscopic time scales. The electronic and mechanical descriptions are connected by the forces that are determined by the atomistic model and by stochastic tunnelling events that describe charge transfer between the electrodes and the shuttle. The mechanical model also incorporates a dissipative term which prevents catastrophic runaway by transferring mechanical energy from the shuttle to lattice vibrations in the electrodes.
A more accurate treatment of this problem in a future study would be to incorporate a time-dependent density functional theory (TD-DFT) [20] module that describes both the island and the electrode during the crucial parts of the shuttle cycle; at present, however, that type of description is prohibitively expensive from a computational point of view.

Set-up
For computational efficiency, we have chosen to focus on the simplest possible system where the central island comprises just one atom. However, the methods and qualitative results should be applicable also to more complex systems. The system we consider includes two electrodes 15 Å apart, described as semi-infinite jellium slabs. The central island is a copper atom that can move in a direction normal to the electrode surfaces. For the electrodes, the Wigner-Seitz radius is set to 3 au and the electrode work function W is set to 3.5 eV [21]. The Fermi energy, F , is calculated from the Wigner-Seitz radius while other material specific electrode parameters are taken from gold. A bias voltage of 3 V is applied over the gap, the potential dropping from left to right. The region described by the DFT module consists of the region between the two electrodes plus a buffer region of 2.5 Å inside each of the electrodes as shown in figure 1. The buffer regions are needed so that the plane-wave-based code can better describe orbitals that are localized in the inter-electrode gap.

Electronic structure calculations
The electrodes and the space between them are treated separately: the electronic structure of the shuttling object is described using a density functional code that is limited to the interelectrode space while the electrodes themselves are described analytically within the jellium approximation. For the island electronic structure calculation we use the DaCapo code [22] with the Perdew-Wang 91 exchange correlation functional [23] and employ the adiabatic approximation to separate the electronic structure calculation of the island from the dynamics description. The adiabatic approximation cannot, however, be invoked to separate charge transfer between the island and the electrodes from the dynamics due to the fact that the relevant charge transfer rates depend exponentially on the island-electrode separation.
The effect of the electrodes on the island is modelled as a one-electron potential and inserted directly into the effective potential in the Kohn-Sham equations. The potential is divided into two parts: the interaction between an electron and a metal surface, and the interaction between an electron and the induced charge arising from the remaining charges in the space between the electrodes. For the former part we use the saturated image barrier suggested by Jones and Jennings [24]. Here, the parameters are related to the work function W and Fermi energy F of the electrodes where the parameters z 0 and λ are found by fitting (1) with DFT calculated tabulated data on the effective potential, v eff [25]. The values were determined to be λ = 1.6 Å −1 and which multiplies the classical image potential, tends towards zero as the mirrored charge closes in on the surface. The image plane, z 0 , is used as the effective surface of the material, in compliance with Lang and Kohn [26]. The value of z 0 derived for (1) differs from the value derived by Lang and Kohn, for e.g., image charge potential [27]. However, for the interactions with the induced charge, the individual contributions from the electrons and the nucleus should cancel out far from the surface and it is therefore suitable to keep only one parameter for the effective surface. For the second term, the potential at r due to a charge at position r , we use a form v(r; r ) that satisfies Poisson's equation in the region outside the electrodes and reduces to the saturated form (1) if the mirrored charge is at the same point where the potential is measured, The resulting potential for interaction with the core is where Z is the position of the island and Q v is the core charge of the island, i.e., the number of protons minus the number of core electrons. The model for the interaction with the induced charge from the other electrons is obtained similarly. Here, we assume that the electron distribution of the island is well approximated with a Gaussian, as is typically the case, and integrate over the mirrored charge to obtain where is a form function and Q e is the number of valence electrons for the central island. The width σ(Q) is calculated by fitting a Gaussian to the unperturbed electron density obtained by DFT for different charge states, Q denoting the number of extra charge units on the central island.
Assuming localization of electrons to the island or the leads (no chemisorption), Q is strictly an integer. The added bias voltage is where the electric field is E = V bias / gap and gap is the distance between the electrodes. Finally, the Pauli repulsion that confines electrons within the gap is for the electronic structure calculations described as a repulsive square potential wall placed at z b . The repulsion term is important for the separation of electronic structure calculation and charge transfer mediated by tunnelling, however, the form of the repulsion is less important as most of the tunnelling events take place when the island is relatively far from the surface on an atomic scale. Owing to the time-consuming and time-independent character of DFT, it is not practical to determine the system properties continuously. Instead, the simulations are performed on a number of positions and charge states. A continuous description is produced with interpolation.

Forces
The forces on the atomic core are calculated using the DFT code and later inserted in the dynamic module. In the Kohn-Sham single-particle equations all charges but one electron are treated as a mean-field static charge distribution. We implement the potential due to surface interactions as an external field. This results in proper (mean field) orbitals and energy eigenvalues but incorrect forces and total energies. Therefore, a correction term in the form of ∂V DFT ∂Z − ∂V el ∂Z is added to the forces calculated by the electronic structure code. Here, V DFT is the induced charge potential used in the DFT calculation, the sum of (1), (3) and (4). The actual potential, V el is obtained as a variational derivative of the full electronic energy as where g(r) comprises all terms independent of Z which do not contribute to the forces. A smooth many-body short-range repulsion is added to the total force. We have chosen a Born-Mayer type pair potential [28] and integrated over the electrode surface which yields the force For a Cu-atom outside a gold surface, the effective radii are σ Cu = 0.77 Å and σ Au = 1.37 Å, n s = 2a −2 is the surface density of atoms and a = 4.078 Å is the lattice constant for a face-centred cubic [100] gold surface [29,30]. The closest electrode lattice site is considered separately (in order to maximize the localization) while the other lattice sites are handled as a continuum. The values f 0 and ρ are 6.95 10 −11 N and 0.3 Å respectively [31]- [34].

Transition rates
Charge transfer rates between the central island and electrodes are calculated with the transfer Hamiltonian method [35,36] using the overlap between the atomic and electrode wavefunctions. The transition rates to and from a specific atomic orbital are where Here, ψ p,a is a Kohn-Sham orbital (localized on the island) with eigenenergy E p , ψ q,m is a metal wavefunction and n f is the relevant number of states available for tunnelling in the (often degenerate) orbital ψ p,a . The potentials V a and V m are the effective potentials for ψ p,a and ψ q,m respectively. The atomic orbitals ψ p,a and the effective gap potential V a are extracted from the DFT simulations. The transition rates are calculated numerically for all energetically allowed transitions determined by comparing the electrode chemical potentials with the Kohn-Sham eigenvalues.
The electrode wavefunctions are calculated analytically from a square potential of length L with a finite barrier of height ( F + W) at the electrode surface and a hard wall at −L. The electrode wavefunctions in the z-direction are with κ z = (2m( F + W)/h 2 − k 2 z ) and δ z = arcsin (hk z / √ 2m( F + W)). A small offset between the physical and the geometrical surface is implemented to attain charge neutrality of the jellium slabs [21]. Periodic boundary conditions are assumed parallel to the surface.

Dynamics
A dynamic Monte Carlo approach is used to calculate the shuttle dynamics, in the spirit of the method proposed by Tully [37]. Input parameters to this module consist of forces and transition rates as functions of the island position. The motion of the central island is described classically by where F ext are the core forces given by the previous calculations and F dissip is a dissipation term. The island position x(t) is calculated by numerically integrating the equation of motion, while the island charge Q(t) is allowed to change stochastically using the tunnel rates determined above. This results in a coupled stochastic dynamics for the mechanical and electrical degrees of freedom.

Dissipation.
Since the shuttle absorbs energy from the bias voltage, a dissipation term is essential for the stability of the island motion [3,38]. Earlier theoretical study has mainly used viscous damping, −ηẋ [3,12], [39]- [41]. While viscous damping is often used in macroscopic systems, its justification is less clear on the nanoscale where the damping arises from the interactions between the small subsystem (the shuttling island) and the large system that is assumed to stay in a near thermal equilibrium (the electrodes). The damping rate is therefore largest when the interaction between the subsystems is strongest, i.e., when the shuttling island is close to the electrodes. The model we employ is based on mechanical damping motivated by phonon emission into the electrodes. An additional electronic dissipation mechanism is due to the creation of electronhole pairs in the electrodes. The relative importance of these mechanisms scales roughly as (m s /m)(πh/mvd) where m s is the mass of a substrate atom, v the impact velocity and d the distance that the shuttle impinges into the surface [42]. For our case, m ≈ 0.5 m s and the electronhole pair creation contribution is only a few percent. Upon impact, the shuttling island exerts a force on the electrode surface which sets a part of the surface in motion (phonon emission). The size of the surface area that is appreciably affected by interactions with the shuttle changes with the island-surface separation. We incorporate these physical effects in a simple phenomenological dissipation model where the surface is described by a single degree of freedom X(t). The mass M(x − X) associated with this degree of freedom is determined so that the acceleration of the single surface degree of freedom subjected to the force felt by the entire electrode surface equals the acceleration of the surface atom that experiences the maximum force (i.e. the surface atom closest to the shuttle island). The resulting equation of motion for the island (x, m) and the surface (X, M) degrees of freedom reads where K = kM(X)/m is the effective elastic spring constant for the surface degree of freedom and the -function restricts the flow of energy so that island motion is dissipated. The reverse energy flow, from an oscillating electrode surface to the shuttling island, would correspond to the island absorbing phonons from the electrode. Since we assume that the electrode surfaces equilibrate rapidly to their ground states (T = 0), the absorption processes are excluded. The atomic spring constant k is given by k = 8mv 2 s /a 2 where a is the lattice constant and v s is the sound velocity. All parameters for the electrode surface atoms have been chosen to correspond to those of gold.
The choice of the dissipation term has rather a small influence on the behaviour of the system, and it mostly changes details such as the threshold voltage for onset of shuttling. It also renders the threshold dependent on the initial conditions. Note that for separations beyond the range of surface forces our dissipation model implies free propagation of the shuttle, while at very small separations the total force is dominated by the maximum force so that M ≈ m and K ≈ k. At intermediate separations the damping coefficient goes smoothly to zero with increasing separation.

Results
The different parts of the one-electron potential are depicted in figure 2. The small widths of the form function (σ ∼ 0.24-0.27 Å for an unperturbed pseudopotential, Q = 1 · · · − 1  and ∼0.30-0.35 Å for the double junction potential) imply that sufficiently far from the surface the point charge approximation would be quite accurate: for distances 2 Å from the surface, the spatial distribution of charge has little effect on the potential. For island positions close to the electrode surfaces, the spread in the valence electron distribution and the rapid saturation effectively bare the core image making the effective potential strongly repellent.
The resulting Kohn-Sham eigenvalues, depicted in figure 3, are used as energy spectra for the central island. Comparison between the eigenvalues and the electrode chemical potentials gives the possible transitions. Full relaxation into the N/2 lowest bands is assumed to be instantaneous, where N is the number of valence electrons (11 for the used Cu generalized gradient approximation pseudopotential, Q = 0). Higher bands are treated as excitations [43].  The temperature is taken to be zero in the treatment of tunnelling events; however, in the DFT calculation a finite temperature is used to improve convergence.
A small correction is needed for some calculations in order to use the equilibrium DFT calculations within a dynamic picture. For some island positions and Q 0, it is energetically favourable to place some of the extra charge in the surface potential well outside the positive electrode surface instead of on the central island. However, the time scale for this direct equilibration between leads is very long. In order to find the energy spectra and orbitals that are relevant for the dynamic evolution, the surface well near the left electrode surface is manually suppressed for core positions near the right electrode; due to the polarity of the applied bias, similar problems do not arise for core positions near the left electrode. The possibility of transitions directly between the electrodes is kept, but the transition rates are small enough to be of no importance for the results.
For the tunnelling rate calculation the Kohn-Sham eigenvalues are regarded as electron energies, which is known to be rigorously correct for the ionization potential involving the highest occupied molecular orbital level [43], and believed to be reasonably accurate for the other levels as well [44]. We assume that tunnelling rates are sufficiently low so that the island fully relaxes between each tunnelling event to the configuration determined by time-independent DFT. The time scale for this relaxation is typically in the femtosecond range which is fast compared to the tunnelling rates except for core positions very close to the electrode surfaces; however, It is interesting to notice the asymmetry of attainable charge states that arises from the asymmetry of the energy spectrum of the Cu atom: the dynamical evolution only involves charge states Q = 0 and Q = 1 but not Q = −1. The symmetric expression for charging energy E = Q 2 /2C used for larger metallic grains is only justified if the level spacing on the island is small enough so that the electrostatic energy scales dominate. This asymmetry implies that the shuttling is asymmetric also in the sense that energy is absorbed from the DC field only during half a period which makes the system more sensitive to dissipative mechanisms.
The occupied Kohn-Sham eigenfunctions are identified as d-and s-orbitals in accordance with the expected electron configuration of Cu, 3d 10 4 s 1 . Close to the electrode surfaces, the orbitals deform against the repulsion wall. For all but the closest position to the electrodes the 4s-orbital gives the widest electron distribution and the largest contribution to the transition rates.
The core forces for the central positions are strongly dependent on the delocalization of the electron distribution of the central island ( figure 4). For the positive ion with Q = 1 the dominant force is the electrical bias while the potential of Q = −1 is a nearly symmetric image charge potential. For Q = 0, the sign depends on the description of the surface interactions, and with the interaction model we have chosen the neutral atom feels a slight net force towards the negatively charged right electrode. The repulsion from the surface is dominant for the two outermost positions on either side giving a physisorption minimum between 3 and 5 Å from the electrode surface. The transition rates are much less sensitive to the surface description than the forces (figure 5). Their distance dependence is approximately exponential as assumed by effective theories [45] with slight saturation for core positions nearest to the electrodes with a tunnelling length that is approximately 0.4 Å with some variation for the different allowed transitions. Near the electrodes the energetics considerations inhibit tunnelling, as seen in figure 3, which can be viewed as a molecular equivalent of a Coulomb blockade.
In the dynamics simulations, the calculated forces and transition rates are joined. The result is indeed a stable shuttling regime where Q(t)Ż(t) = 0. We have performed dynamical simulations starting from a variety of initial states, and seen that for most starting conditions the results are quite similar: as a rule, the model does shuttle electrons (figure 6). However, for some initial configurations such as Z(t = 0) ≈ 10 Å and Q(t = 0) = 0 the applied bias of 3 V is not sufficient to initiate shuttling.
One of the differences between our results and previous studies is the complexity of the forces acting on the shuttle, particularly near the electrode surfaces. The results are quite sensitive to the details of the surface interaction. A slightly different potential may render the forces on the neutral atom positive over a larger range of positions, and the range of initial conditions that would result in shuttling would be smaller, implying that the threshold voltage for shuttling depends sensitively on the model for surface-island interactions and on the initial conditions.
The detailed structure of the forces of the middle positions is less important after shuttling is well established. The main forces become the close-range exponential forces of the electrodes and the applied electric field. For the asymmetric shuttle, energy is absorbed by the charged shuttle from the field during half a cycle while during the other half-cycle, after an elastic collision with the electrode surface, a neutral shuttle moves nearly freely in the opposite direction. The energy loss during the shuttle-electrode collision cannot exceed the energy absorbed from the field if a stable periodic motion is to be established.  The distribution of positions at which tunnelling events take place depends on the transition rates, and the width of the distribution is connected to the spatial derivatives of the rates, i.e., on the tunnelling lengths. The steeper 1 → 0 gives a more compact distribution as seen in figure 7.
The shuttle reaches a stable shuttling motion quickly with a current of ∼0.19 µA and an amplitude of ∼11.1 Å, (see figure 8). The random character of the transition processes influences the turning points very little. There is, however, a possibility for the system to undergo semistable excitations due to randomness of transfer events and the position dependence of the energy spectrum near the right lead ( figure 9). Very close to the negative (right) lead, there is a possibility of a process in which an electron first tunnels from the electrode to an initially positively charged (Q = 1) shuttle that continues to move towards the right lead, followed by tunnelling against the bias back into the lead, and finally a new tunnelling event after the shuttle has changed its direction of motion. During the time that the charged shuttle spends near the electrode surface after the second tunnelling event, it experiences a larger force than a neutral shuttle would, which allows it to absorb more energy from the potential and results in an enhanced shuttling amplitude. The increase in amplitude enhances the possibility for this sequence of three tunnelling events to take place also in the next period. The excitation lasts until a transfer without reverse tunnelling takes place near the negative lead. For the system we have studied, the amplitude of this excited cycle is about 0.3 Å larger than that of the simple cycle, and the current level is increased by approximately 20% to 0.23 µA. The possibility of two stable shuttling amplitudes has recently been discussed by Usmani et al [46].

Discussion
The energy spectra and the (Kohn-Sham) orbitals of small molecules near metal surfaces exhibit a rich structure and vary substantially as a function of the molecule-metal surface separation. The transition rates are largely exponential functions of the tunnelling distance as assumed in phenomenological theories, but the allowed transitions are determined by the energy spectrum, and in particular near the surfaces certain transitions are forbidden by energy considerations. This results in an asymmetry in possible charge states and in asymmetric shuttling where energy is absorbed only during one half cycle of the periodic motion. The average current for a path with system excitation due to backtunnelling close to the negative lead (see inset). The increase of amplitude caused by the larger Q = 1 forces enhances the possibility for the tunnelling triplet to take place also in the next period. The transition is sharp with an amplitude increase of ∼0.3 Å and a current of ∼0.23 µA. The excitation lasts until a period without reverse tunnelling takes place.
The forces in the system are highly sensitive to the description of the electrode-molecule interaction and to the electronic structure of the shuttling object. For the particular system that we studied we found that the attractive surface forces can be overcome by reasonable electric fields and shuttling of a nanoscale object can be established. In the stable shuttling regime the island velocity is large enough to overcome physisorption allowing the island to bounce between the electrodes. The details of the forces near the middle of the gap are less important than the balance between dissipation and short range surface forces.
The shuttle excitations depicted in figure 9 are an example of effects that arise due to the details of the energy spectrum of a small system. For a more complicated spectrum and larger bias voltage more phenomena of the same type can be expected; for a slightly different model of interactions between the shuttle and the electrodes we have even observed that the regular shuttling motion may pass into a more chaotic behaviour. Therefore, it is likely that a microscopic picture of both forces and energy levels is paramount for both quantitative and qualitative predictions of molecule-sized shuttles.