Single molecular shuttle-junction: Shot noise and decoherence

Single molecular shuttle-junction is one kind of nanoscale electromechanical tunneling system. In this junction, a molecular island oscillates depending on its charge occupation, and this charge dependent oscillation leads to modulation of electron tunneling through the molecular island. This paper reviews recent development on the study of current, shot noise and decoherence of electrons in the single molecular shuttle-junction. We will give detailed discussion on this topic using the typical system model, the theory of fully quantum master equation and the Aharonov-Bohm interferometer.


Introduction
The experimental and theoretical study of molecular junctions has become an important topic of modern condensed matter research because of the potential applications for precision measurement and nano-electronic mechanical devices. A single molecular junction is one of the simplest structures, which can be fabricated by connecting two electrodes with a single molecular island. A bias voltage is loaded on the two electrodes, and a gate voltage is used to switch the energy level of the molecular island. The fabrication of molecular junctions utilizes a bottom-up approach. In this approach, elementary molecules are arranged into more complex self-assembled structures using chemical synthesis. In contrast, the topdown design approach would cut larger raw materials into desired devices using physical methods. Unfortunately, modern lithography tools cannot easily reach the scale associated with the synthesis of electronically functional molecules. One further advantage of the bottomup design approach is that nanoscale devices may be fabricated exactly. For example, it is possible to make devices with any degree of accuracy. The molecular islands in the junctions are generally nanometer scaled and have discrete electronic energy levels. Thereby, they work as quantum dots (QDs) with three-dimensional bounded barriers and discrete energy levels. The single molecular junction has another important degree of freedom, namely mechanical vibration of the center of mass of the molecular grain. This vibrational motion plays an important role in the transport of molecular junctions and leads to interesting phenomena.
The mechanical vibration in the molecular junction makes it possible to control electron transport through nanoscale devices in the single-electron regime [1]. Such system are also known as shuttle-junctions. In the semiclassical regime, time evolution of the shuttle-junction shows electron transfer one-by-one orderly [2,3]. Furthermore, the tunneling-position coupling in the electromechanical junction has significant applications in the detection of displacement [4][5][6][7], spin [8,9], charge [10,11], and mass [12]. The measurement correction can be better than the zero-point-motion uncertainty. In addition, vibrational junctions have abundant physical applications such as heat transfer and cooling based on nonequilibrium hot electron transport [13,14], subresonance inelastic electronic transport that is analogous to atomic laser-induced cooling [14,15] and tip-induced cooling [16]. The effective temperature of the molecular island is defined as the internal energy stored in the vibrational degrees of freedom. Moreover, electronphonon coupling enhances the charge transport and increases the temperature of the molecular junction. However, the charge transport leads to inelastic phonon absorption and decreasing temperature [17,18]. In particular, when the energy of incident electron is lower than vibrational level, the electron takes a phonon from the junction, inducing cooling at the molecular junction. The cooling effect can be very significant: resonant tunneling in a single-molecule device can generate sufficient heat to thermally decompose the molecular junction.
Theoretically, Green's functions are used to study the nanometer electromechanical junction in classical [19], semiclassical [20], and quantum regimes [21]. These approaches allow the calculation of transport properties on the molecular junction, such as weak to strong electronphonon coupling [21,22] and inelastic electron tunneling spectroscopy [23,24]. A real-time path-integral Monte Carlo approach can also be applied to study the transport process in the molecular vibrational junction over a large temperature and electron-phonon coupling range [25]. Quantum master equations, which stem from the Liouville-von Neumann equation [26,27] and the many-body Schrödinger equation [28], are widely applied to describe electromechanical vibrational junctions. In the junctions, charge-vibration (or electron-phonon) couplings play an important role, which can be treated in both incoherent [29][30][31][32][33] and coherent approaches [34,35]. In the incoherent approach, the off-diagonal couplings between electron tunneling and the vibrational states are not included. On the other hand, in the coherent approach, the off-diagonal couplings are considered. Using the master equation, fascinating electromechanical system features have been predicted, such as negative differential conductance [29,30], the shuttling effect [1,36], the super and sub-Poissonian Fano factor [2,32,34,37,38], spin-dependent transport [8,39,40], vibrational state instability [41], and negative damping instability [42]. In a fully quantum mechanical description of the electromechanical process [43], the coherent coupling is included at the cost of solving a large matrix in the mathematical treatment. The quantum mechanical model of the coherent dynamics is further developed to investigate the shuttling mechanism [3,36,37]. According to the charge-position (momentum) correlation, the motion of the quantum shuttle can be described by the shuttling and tunneling regimes as well as the coexistence of these regimes [2,36]. In the shuttling regime, the electron transport is highly deterministic and characterized by the extraordinary sub-Poissonian Fano factor [37]. By measuring the shot noise, the transition between tunneling and shuttling can be identified [3]. Moreover, single molecular vibrational junctions have recently been studied in strong correlation physics such as the Kondo effect [44][45][46][47][48], pair tunneling, and co-tunneling [49,50].
The molecular electromechanical junction in electromagnetic fields is also very interesting. Surface-enhanced Raman scattering is used to study the mechanisms in molecular junctions both experimentally [51][52][53] and theoretically [54][55][56]. The approaches based on surfaceenhanced Raman scattering have been widely applied in the study of equilibrium materials over the last 30 years [57][58][59][60][61][62][63]. Recently, these methods have been applied to the nonequilibrium charge transport process in molecular devices. Theoretically, the electron population at the molecular level and changes in the heat can be detected with this approach [54][55][56]. Under a magnetic field, the spin degree of freedom can be exploited. Based on the coupling between the transport of spin-polarized electrons and the mechanical degrees of freedom of the island, the shuttle instability is predicted to have two stationary domains such as vibronic and shuttling domains depending on the applied electric and magnetic fields, which suggests an extension to spintronics applications. Therefore, the coupling between the spin polarization of the current and the motion of the molecular island center of mass can be used to control the dynamics of the mechanical degrees of freedom using an external magnetic field [64]. The spin-dependent current in the external magnetic field has been recently studied [8,39,40]. In a magnetic shuttle device, the spin-exchange interaction causes spin-dependent tunneling rates, which leads to spin-polarized current [65]. This effect is known as the spintromechanical effect.
Experimentally, nanometer-scale electromechanical devices are fabricated by organizing molecular building blocks [66][67][68][69][70][71][72][73][74][75][76][77][78][79][80] and by optical or electron beam lithography methods [81][82][83]. Here, we briefly present two typical experiments on nanoscale electromechanical systems. In 2000, Park et al. fabricated a single molecular vibrational junction [69]. The vibrational molecule C 60 is connected to two gold electrodes (see Fig. 1). The gold electrodes were fabricated using electron-beam lithography, and the gap between the electrodes was created using the process of electromigration. The narrowest gap is about 1 nm, the lateral size of the electrodes near the gap was on the order of 100 nm, and the height of the electrodes was about 15 nm. The experimentally determined frequency of the C 60 mechanical oscillator is f = 1.2 THz at a temperature of 1.5k B T . Considering the constant = 4.135660 eVs and k B = 0.861739 eV·K −1 , the energy quanta of the oscillator is f = 4.9628 × 10 −3 eV, which is much larger than the heat energy k B T = 1.2926×10 −4 eV of its environment. At low temperatures, the molecular oscillator excites quantized energy levels, which leads to the step-structure in the current shown in Fig. 1. In 2009, Moskalenko et al. reported a fabrication method for an electromechanical shuttle-junction using gold grain [82,83]. Comparing with the C 60 molecule (1 nm order), the gold grain is about 20 nm in diameter. Therefore, it is easy to image the gold grain system. Using atomic force microscopy, three-dimensional images of the electrodes are obtained shown in Figs. 2(a), (b), and (c). The process for positioning the gold grain between the electrodes is illustrated in Fig. 2(d). The frequency of the nanomechanical oscillator is on the order of 10 8 Hz, so the energy quanta is 10 −7 eV. The energy is much smaller than the heat energy at room temperature Fig. 1 The I-V curves of the C 60 transistor at T = 1.5K, reproduced from Ref. [69], Copyright c 2000 Nature Publishing Group. Fig. 2 The process of embedding a 20 nm gold gain between two electrodes. The picture is obtained with atomic force microscopy. Reproduced from Ref. [82], Copyright c 2009 American Physical Society. k B T room = 2.5852 × 10 −2 eV, where T room represents room temperature. The experiment is performed at room temperature. Though the Coulomb blockade effect does not appear, it is difficult to achieve single electron processes. To realize single electron transport, the temperature should be decreased and the surface roughness of the gold grain and electrodes be removed. As shown in Fig. 3, for a low applied voltage, the gold grain is almost static, and the current is low. After the voltage exceeds a certain threshold value, the grain was driven to oscillate and promote the current.
There are two commonly used models for the study of nanoelectromechanical devices. In the first model, the electron tunneling amplitudes are considered to be independent of the nano-particle displacement [25,31,[84][85][86]. In the second model, the tunneling amplitudes are modulated by the displacement x of the center of mass of the nano-particle island [1,34,35,41,87]. For example, the modulation has an exponential form of T l e −x/λ and T r e x/λ , where T l and T r are bare tunneling rates on the left and right sides of the island, and λ is the characteristic length of electron tunneling. This latter model is known as the shuttle-junction model, which is the focus of the following discussion. In this paper, we review recent developments on the study of current, shot noise, and decoherence of electron transport through the electromechanical shuttle-junction using the quantum master equation. The Green's function, path integral, and master equation approaches are applied to study the model in which electron tunneling is independent of the island displacement. However, the shuttle-junction properties can only be described using the quantum master equation. The applications of this method for describing nanoelectromechanical junctions have been developed within the past 15 years. In this study, we develop a general master equation for the description of the electromechanical tunneling. This general equation considers the Fermi distribution functions of the electronic leads in which discrete energy levels of the vibrational modes will be included. When the bias voltage is low, the distribution functions are sensitive to the oscillator levels. In addition, both the diagonal and off-diagonal terms of the coupling between the oscillator dynamics and electron transfer are incorporated in the equation. Both the electron transfer and oscillator dynamics are shown to be very important for the device description. In a recent attempt of describing the electromechanical system, the position dependence is neglected, and only two oscillator modes are considered [85]. We will show that the position dependence of the tunneling rates and all vibrational modes of the harmonic oscillator can be treated within the master equation approach. Furthermore, the master equation can equally consider the tunneling and dissipation terms, unlike early derivations, which required treating these terms separately [36,43,85]. In the low temperature and high bias voltage regime, the presented master equation should agree with those given previously [3,36].
Previous studies on the nanoelectromechanical junction have primarily concentrated on the incoherent properties of electron transport. Both the coherence and dephasing character of the electron transfer will be discussed in this paper. By applying an Aharonov-Bohm (AB) interferometer with a QD embedded in one arm, the coherence character of electron tunneling through the QD can be studied in experiments [88][89][90]. These experiments show that a fixed QD supports coherent transport and causes the phase-shift of an electron (see Fig. 4). When a QD is allowed to mechanically oscillate around its equilibrium point, an electron transferring through the dot would be randomly accompanied by the absorption or emission of phonons. The phase property in mechanical vibration-assisted electron tunneling is still an open interesting question. The vibrational motion in the system can be approximated using a monochromatic oscillator. This vibrational motion differs from thermal fluctuating Bosonic baths, which cause decoherence to local electronic states of charge [91,92] and spin [93,94]. As recently reported, a single vibrational mode of a QD array enhances the electron transport and partially preserves its phase information [95]. The coherent transport of electrons in QDs is also sensitive to spin flip, electron-electron interaction, and external detectors [96][97][98][99][100][101][102][103][104]. A cantilever-based which-path charge detector has been previously studied [105,106]; this detector is based on dot-cantilever coupling, which causes remarkable dephasing to electrons. In the Armour et al. model [105,106], the dot-lead coupling does not depend on the oscillator position; therefore, this system differs from that considered in this study. Moreover, there are other excellent reviews [86,[107][108][109][110] presenting different features of the vibrational electronic junction.

Quantum master equation for single molecular shuttle-junction
Electron transport generally involves two electrodes separated by a conductor. The conductor has finite degrees of freedom, whereas the electrodes have infinite degrees of freedom. In this configuration, electron transport is similar to an open system process, where the conductor is the system and the electrodes are the environment. Typical examples of electron transport are transport through a QD and point contact. Therefore, the theory of an open system, such as the master equation, is suitable for describing electron transport. In fact, the master equation has been widely applied in the area of transport materials [3,28,[111][112][113][114]. In this section, we introduce a typical shuttle-junction model and develop a quantum master equation to discuss properties of single molecular shuttle-junctions.

Model
In theory, a single molecular shuttle-junction is modeled as a QD connected to two electronic leads with elastic materials. The QD center of mass is allowed to mechanically oscillate. In particular, the mechanical vibration of the dot can be regarded as a harmonic oscillator with an effective mass m [1,69,86,108]. At finite bias voltage V , the quantized energy of the mechanical oscillator plays an important role, especially when the bias voltage is low [34]. The energy structure of the conductor is shown in Fig. 5. With respect to the QD energy levels, the chemical potentials in the left and right leads are eV /2 and −eV /2, respectively, where e is absolute value of the electron charge. We assume that the capacitance of the QD is so small that its electron occupation number is 0 or 1. Initially, the QD is empty, and the harmonic oscillator is in its ground state. When an electron jumps onto the QD, the electric field caused by the bias voltage exerts a force on the charged dot and drives the mechanical oscillator into an excited state. As a consequence, the electron transfer is influenced by the vibration. The typical model of the single molecular shuttle-junction is [3,34,36,40,41,108]  which is sum of the electron tunneling Hamiltonian the mechanical oscillator Hamiltonian and the coupling between charge and the oscillator The first term in H tun is the energy of the QD with annihilation and creation operators c and c † , respectively. The second term in H tun describes the non-interacting electrons in the left (y = l) and right (y = r) leads. The operator d † yk (d yk ) creates (annihilates) an electron with momentum k in lead y. The spin degree of freedom is not involved here. The third term in H tun represents electron tunneling between the leads and the QD. The coupling strength is exponentially dependent on the dot position x, which will also be expressed using the second quantization form x = x 0 (a+a † ) in the remainder of this paper. Here, a and a † are the annihilation and creation operators of local phonons, respectively, and x 0 = 2mω0 is the zero-point position uncertainty with the Plank constant and the intrinsic oscillator frequency ω 0 . The characteristic tunneling length is given by λ. For simplicity, we define S l = −1, S r = 1, and α = x0 λ . The first term in H mech is the free evolution of the mechanical oscillator. The creation and annihilation operators, b † k and b k , respectively, characterize the Bosonic thermal bath. The last term in H mech represents oscillator dissipation due to the bath. We assume the electric field intensity is created from the bias voltage V over the distance between the two electrodes. Other electric environments of the QD may cause a bias-voltage-independent electric field. These environments also contribute to the chargeoscillator coupling, which primarily acts as the gate voltage by shifting the QD level and changing the number of vibrational states accessible for the electron transport [30]. We do not discuss this effect in detail here.

Equation of motion
Master equations are believed to be convenient and direct methods for describing complex motion. For a large bias voltage, the system can be studied with a semiclassical approach by treating the mechanical oscillator using a Newtonian equation and by treating electron tunneling with the master equation [1]. In this regime, the discrete energy levels of the oscillator are not obvious, and the system manifests a clear shuttling effect. The oscillator distribution in phase space is studied, and the results show that the current and shot noise are connected to the oscillator states [36,37]. For low bias voltage, only a few energy levels of the mechanical oscillator are involved in the transport, and quantum effects become more important. In the frame of the master equation, even incoherent models of the theoretical treatment can produce results that qualitatively fit those of experiments [29][30][31][32][33]. A general master equation, which can tune the bias voltage and temperature of the system and also includes coherent dynamics of the oscillator, has been developed to study electromechanical devices [3]. However, in Utami et al. [3], the Fermi-Dirac distribution functions do not consider the energy levels of the oscillator. Strictly speaking, oscillator levels are important for electron occupation, especially in areas of low applied voltage [29,30,85].
In what follows, we will show that a general master equation can be developed to describe the electromechanical tunneling [34,35]. In this approach, we employ the Fermi distribution functions of the electronic leads involving the discrete energy levels of the vibrational modes. In addition, both the diagonal and off diagonal terms of the coupling between oscillator dynamics and electron transfer are included in the equation as both terms are shown to be important for describing the device. Moreover, the position dependence of the tunneling rates and all vibrational modes of the harmonic oscillator are considered in the equation. This approach overcomes the shortcomings in a recent attempt to describe an electromechanical system where the position dependence is neglected and only two oscillator modes are considered [85]. Furthermore, the present master equation was obtained by treating the tunneling and the dissipation terms concurrently, unlike the separate prescription adopted in earlier derivations [36,43,85]. However, in the low temperature and high bias voltage limit, the present master equation should agree with those given previously [3,36]. Now, we proceed to derive the master equation using the Liouville-von Neumann equation: Here, ρ T (t) is the density matrix of the total configuration. Unlike traditional methods of deriving the quantum master equation [26,27], we divide the Hamiltonian into three parts instead of two parts. The three components are the free evolution Hamiltonian H 0 , the charge-oscillator interaction Hamiltonian H driv , and the interaction Hamiltonian H 1 : Note that H driv is given in the previous section. For convenience, the Hamiltonian H 1 is further divided into two parts as and In the interaction picture, according to O(t) = e itH0/ Oe −itH0/ for any operator O, these Hamiltonians become and where , which satisfies the following Liouville-von Neumann equation: Integrating Eq. (13) over time t 1 and substituting the result into the second commutator of Eq. (13), we obtain The first commutator in Eq. (13) is retained in Eq. (14), and the integral equation is substituted into the second commutator of Eq. (13). Note that this step is not performed in the traditional derivation of the master equation [26,27]. However, this step is valid because the Hamiltonian H driv is a measurable quantity, and together with the system free Hamiltonian, this Hamiltonian is strictly solvable. In other words, the following equation can be analytically solved: Therefore, we do need to iterate the first terms in the right side of Eq. (13). On the other hand, the integration of (13), which is substituted into the second commutator of the equation, is a standard approach in the master equation derivation. This process not only allows us to trace over the environment to obtain the non-zero terms and the equation of motion for the system but also creates measurable quantities such as the Fermi distribution function and tunneling rates. There are two environments for this system: the first environment is an electronic reservoir, such as that formed by the two leads of electrons, and the other environment is the Bosonic bath, which interacts with the mechanical oscillator. Under the Born approximation, the total density matrix can be factorized as ρ is the density matrix of the system composed of the QD and the mechanical oscillator. Both the leads and the Bosonic bath are always assumed to be in an equilibrium state described by the density matrix ρ L and ρ B , respectively. Therefore, if ρ T (0) = ρ(0)ρ L ρ B initially, then ρ T (t) = ρ(t)ρ L ρ B at time t. This approximation is valid as long as the coupling is weak and the environment is sufficiently large that the back action from the system to the environment is negligible. Performing a trace over the leads (tr L ) and bath (tr B ), we obtain the reduced density matrix for the system We further assume that two reservoirs are uncorrelated. Moreover, the reservoirs have no memory and do not preserve system information. As a result, the interaction between the system and the reservoirs is not affected by the system history. This effect causes an open system to lose coherence. In this limit, we can use the Markov approximation and replace ρ(t 1 ) by ρ(t) to arrive at the equation for the reduced density matrix:

Fig. 6 Time replacement in the integral
After evaluating the commutation relation (see Appendix 1) in Eq. (17), we can transform the equation to the Schrödinger picture using the unitary operator e −it(ε0c † c+ ω0a † a)/ . For convenience, we make a timedisplacement transformation τ = t − t 1 (see Fig. 6). Denoting the electron number in the left and right leads as w and v, respectively, we write the density matrix in the form ρ w,v (t). The density matrix is obtained from the following simple relation: for a closed, isolated equilibrium system, the traces tr L (d † rk d rk ρ L ) and tr L (d rk ρ L d † rk ) are the same. For the electronic leads in an open system, tr L (d † rk d rk ρ L ) and tr L (d rk ρ L d † rk ) have the same value, but the traces include different information about the electron number in the leads and QD. For a given electron number v in the right lead, we have In the first term, there is no change in the number of electron in the right lead; however, in the second term, one electron hops into the lead from the QD. Comparing tr . If one electron is annihilated (created) in the left lead, the density matrix becomes ρ w−1,v (t) (ρ w+1,v (t)), and if an electron is annihilated (created) in the right lead, the density matrix becomes ρ w,v+1 (t) (ρ w,v−1 (t)). Note that the density matrix satisfies the normalization condition ∞ w,v=0 ρ w,v = ρ. The counting method introduced here is similar to the counting approach presented in the theory of the many-body Schördinger equation [28]. Additionally, the Fermi distribution function in lead y is in- The Bose distribution function in the thermal bath of the mechanical oscillator is obtained from . Because the tunneling amplitudes exponentially depend on position of the molecular grain, there are exponential factors in the equation. We expand the exponential terms using the following equations: Considering the rapidly decaying electrons in the reservoirs, the time integral is extended to the infinite regime, i.e., t → ∞. The integrals over time τ are performed using the formula ∞ 0 dτ e ±ixτ = ±iP 1 x + πδ(x) for any variable x. The imaginary part ±iP 1 x corresponds to the Lamb shift in quantum optics [26,27]. Because we consider weak coupling between the system and the reservoirs (electronic leads and Bosonic thermal bath), the Lamb shift is very small and can be neglected. Subsequently, we define the density of states in lead y as N y (ξ yk ) and the density of states in the thermal bath as D(ω k ). By converting the summation to an integration k;y=l,r ∼ dξ yk N y (ξ yk ) and k ∼ dω k D(ω k ), we integrate over the electron and phonon variables, respectively. Under the wide band approximation, the tunneling rates are written as Γ y = 2π |T y | 2 N y (y = l, r) and the dissipation rate is γ = 2πDg 2 , which is independent of energy spectrum. Finally, we achieve the master equation: where we define and We now sum over all possible values: The commutator on the right side of the equation denotes the system evolution that is driven by an electric field. The second term represents phononassisted tunneling through the electromechanical junction in which the operators A − mznz and A + mznz are defined as Here, A − mznz indicates that an electron has been transported out of the dot along with the created phonon m z and annihilated phonon n z . Likewise, A + mznz indicates that an electron has been transported into the dot. Figure 7 illustrates these operators for any phonon number m and n. The Fermi distribution function is written in the form f y,m1n1 = In the present work, we assume the electronic leads and the thermal bath have the same temperature. However, it is straightforward to extend this to the situation in which temperature in the two reservoirs is not the same.

Numerical methods
Due to the large Hilbert space of the fully quantized system considered here, it is impossible to derive even an approximate analytical expression for the current vs. voltage characteristics. However, we can deduce an expression in terms of the density matrix. To understand the current vs. voltage characteristics, we project the master equation onto the Fock state bases spanning the Hilbert space of the system [115]. The base vectors of the Hilbert is the direct product, |0 and |1 are the eigenstates of the QD, indicating that the dot is occupied by zero and one electron, respectively, and |n is the eigenstate of the nth level of the mechanical oscillator. In the Hilbert space, the system density matrix element is written as ρ ij,mn = m| ⊗ i|ρ|j ⊗ |n for (i, j = 0, 1; m, n = 0, 1, 2, . . .). For any two given vibrational states |m , |n we have four density matrix elements ρ ij,mn (i, j, k, l = 0, 1) in terms of the electron states. However, we only consider the diagonal electron occupation states ρ 00,mn and ρ 11,mn because they are decoupled from the coherent terms ρ 01,mn and ρ 10,mn between the occupied and unoccupied electron states. Obviously, the density matrix contains both diagonal and off-diagonal terms for the vibrational states. In the following, we primarily discuss the stationary solutions of the master equation, and we apply the condition ∂ρ/∂t = 0 (24) in Eq. (20). Therefore, the projection onto the Fock states is performed as follows: where i = 0, 1. We then obtain a set of 2(N + 1) 2 linear equations, where N = 0, 1, 2, . . . is the number of excited vibrational levels considered. After the projection, the equations can be expressed in terms of the density matrix elements (see Appendix 2) associated with the normalization condition, We have directly solved the linear equations. The effective dimension of the Hilbert space is directly related to the excited states of the harmonic oscillator involved in practical transport. However, the master equation given in this paper is valid for a much wider parameter range than that considered here as long as an adequate number of vibrational states are utilized. As more vibrational levels are considered, the computational cost of the numerical implementation increases. The large memory and time requirements are weak points of the approach to directly solve the equation. Note that these problems could be circumvented in the iteration method in which preconditioning is necessary to ensure convergence [116]. The iteration reaches is completed when sum of the diagonal elements of the system density matrix closes to unity. In Fig. 8, the probabilities P n = ρ 00,nn + ρ 11,nn of the mechanical oscillator at different energy levels are plotted for N = 18. For the applied voltage eV = 11 ω 0 , the contribution from levels higher than N = 18 may be negligibly small. For the remainder of this paper, if is it not specifically stated otherwise, only 18 excited vibrational levels are included in the numerical calculation. The resulting Hilbert space is sufficiently large for understanding the system with a low bias voltage and weak dot-lead couplings.

Current through the single molecular shuttle-junction
In this section, we discuss the current as a function of the applied voltage and the oscillator damping rate. The current-voltage curves for different gate voltages are given for comparison with experimental results. Finally, we compare the contributions from the coherent and incoherent dynamics of mechanical oscillator to current.

Current formula
Considering the charge conservation, the stationary current can be calculated from the flow either in the left or the right lead. For convenience, we consider electrons moving in and out of the right lead. The probability of v electrons collected in the right lead is where tr mech is the trace over the occupation states of the nanomechanical oscillator, while tr char denotes the trace over the charge degrees of freedom in the QD. According to the counting theory in the electron transport system with the general formula [117] the current can be expressed by the system density matrix in the following form: In the above equation, all parameters are known, except for the system density matrices ρ 00 (t) and ρ 11 (t). After the density matrices are obtained from equation (20), the rate of electron transfer through the vibrational shuttlejunction can be calculated.

Current
We first focus on the mechanical oscillator. Figure 9(a) shows the averaged phonon number of the harmonic oscillator, which can be calculated from the formula n = ∞ n=0 n(ρ 00,nn + ρ 11,nn ). The mean phonon number is equivalent to the averaged oscillator energy. The energy increases suddenly when the bias voltage increases over an extra resonance level of the mechanical oscillator. However, the range between the two neighboring quantized energy levels is different from the original eigenenergy of the phonon. A similar average phonon number in an nanoelectromechanical system has also been discussed recently for the asymmetric coupling between the island and the two electronic leads [118]. However, their system is not a shuttle-junction because the tunneling length is independent of the oscillator position.
The current is plotted as a function of bias voltage in Fig. 9(b). These curves correspond to the quantized energy in Fig. 9(a) with the same parameters. We found nearly discontinuous transitions at low temperatures, and we found that no current is available at zero bias voltage. When the bias increases from the zero point to a small finite quantity, the current appears and increases sharply. However, when the voltage continues to increase, the current seems to be unchanged. In the voltage range of −2 ω 0 eV 2 ω 0 , the mechanical oscillator is found primarily in the ground state with zero mean phonon number [see Fig. 9(a)], but it still contributes to the current due to the zero-point fluctuation. In Fig.  9(b), we can easily verify that the first step is larger than the current of bare QD tunneling Γ l Γ r /(Γ l + Γ r ). The chemical potential in the left lead now lies between the ground and the first excited states of the mechanical oscillator (see Fig. 5). When the bias voltage increases to 2 ω 0 , the chemical potential reaches the first excited level of the oscillator, which results in a second jump of the current. More current steps emerge as the bias voltage increases. Each time the chemical potential reaches an additional level, an extra transport channel is opened, and a current step is observed. The energy spacing of the oscillator levels is reflected by the voltage range of the current steps. The height of the steps tends to decrease with the increase of the bias voltage. In fact, the height of the steps can be controlled by tuning the parameters Γ l , Γ r , Γ , and α. Comparing Figs. 9(a) and (b), every increase (decrease) of the current in the system is clearly accompanied by an emission (absorption) of energy by the oscillator. Therefore, the current is strongly correlated with the energy of the mechanical oscillator. The key of the theoretical model is that the mechanical oscillator introduces multiple modes into the dot conductor and quantizes the current. Discrete levels of the oscillator play the role of multiple sub-bands in a narrow conductor of a two-dimensional electron gas in which quantized conductance is observed [119,120]. As the left and right parts of the model are absolutely symmetric, the current appears to be antisymmetric for the positive and the negative bias. The antisymmetry could be broken if the left (Γ l ) and right (Γ r ) bare tunneling rates are not equal [31].
In the above analysis, the energy level of the QD is constant. We now study the influence of the QD resonance level on the voltage-current curves. This property has been shown in experiment [69], where levels of the QD can be shifted by tuning the gate voltage. We plot the current as a function of bias voltages in Fig. 10 for

108501-10
Wenxi Lai, Chao Zhang, and Zhongshui Ma, Front. Phys. 10, 108501 (2015) different dot levels. To compare the experimental results with our theoretical calculation, we use the experimental parameters were possible. The remaining parameters are selected with probable values. In Fig. 10 the curves are symmetric for both negative and positive bias voltage because the energy structure in our model is symmetric. However, in experiment, the current for the negative and positive bias voltage does not appear to be symmetric. This asymmetry may be caused by the hysteresis effect in the polarization of the real material. Figure 10 implies that, by adjusting the gate voltage, the current can be controlled with discrete quantities. As discussed in the previous section, the current of the system is correlated with the energy of the mechanical oscillator. Damping of the vibrational mode influences the mean energy of the oscillator as well as the electron  transfer. Figure. 11 illustrates this effect in the current by changing the dissipation rate γ. For small γ, the step structure of current is still clear; however, for large γ, the lifetime of the excited levels are very short, and the steps become unclear. Moreover, a sharp decrease in the current is observed in the area γ < 10Γ l , 10Γ r and eV > ω 0 . In the range eV ω 0 , the mechanical oscillator is primarily found in the ground state; therefore, the current is independent of the dissipation rate. The approximate calculation reveal that the currents for γ < 0.01ω 0 are further intensified but have finite quantities. These results are not shown here as they have been calculated in previous studies [2,36,37].

Coherent and incoherent dynamics
In the end of this section, we discuss the effects that result from the coherent coupling between the charge transport and dynamics of the mechanical resonator. We calculate the system current using an incoherent model in which only diagonal terms of the density matrix, such as ρ 00,nn andρ 11,nn , in Eq. (11) and the current formula are involved. As illustrated in Fig. 12, for a bias eV 4 ω 0 , the contribution from the off-diagonal terms is negligible small. However, in the case eV > 4 ω 0 , the current calculated from the coherent model, which includes both diagonal and off-diagonal terms of the density matrix, is lower than that achieved from the incoherent model.
The suppression of current in the coherent model may be due to the destructive interference between different transport channels. The incoherent model applied here is not absolutely the same as those considered previously Fig. 12 Comparison between the currents calculated from the coherent and incoherent models. The dashed (green) and dot-dashed (red) lines present the currents which only involve diagonal elements of the system density matrix. The solid (black) and dotted (blue) lines indicate the currents that calculated from both the diagonal and off diagonal terms of the density matrix. The corresponding parameters are x 0 /d = 0.005, α = 0.7, ε 0 = 0, because of the different derivation methods [29][30][31][32][33]. However, in this study, we only intend to clarify the importance of coherent coupling in the system.

Shot noise
In this section, the noise of electron transport is analyzed by calculating shot noise in the shuttle-junction. Especially, Fano factor spectrum as a function of bias voltage and temperature are analyzed and presented.

Noise in quantum devices
The shot noise in a classical device is given by Schottky's formula: In a quantum system, the shot noise can be higher than or lower than the value measured by the Fano factor, where ω is the frequency of spectrum. When the electron transfer is completely random, the shot noise is equal to S class . This is known as Poissonian noise, which indicates an uncorrelated transport process. Note that no correlation does not mean no noise is present. The quantum correlation occurs because of the stochastic features of the wave function. When a system has more than one channel and the channels are correlated due to the superposition of the wave function, the process would appear to bunch. In this case, the Fano factor satisfies F > 1, which is indicated as super-Poissonian. If these channels are not correlated, electrons are more likely to transfer independently in any individual channels, and the anti-bunching property of electrons would dominate the transport. In this case, the Fano factor appears to be sub-Poissonian with F < 1.
Concerning the probability of v electrons collected in the right lead, we employ the McDonald formula [121], to calculate shot noise of the system. It is convenient to use the Fourier transformation . Then, using our master equation, the shot noise can be described by the density matrix: Here, we define . From equation (20), it is easy to construct equations for Z 00 (ω), Z 11 (ω), Y 00 (ω), and Y 11 (ω). Using these values, we can calculate the shot noise.

Shot noise and mechanical oscillator
In the high bias voltage and low temperature limit, the shuttle-junction has a very simple distribution in the phase space of the mechanical oscillator [2,36]. According to the distribution, the system motion is described by shuttling, tunneling, and the coexistence of the shuttling and tunneling states. In the tunneling state, the Fano factor is close to 0.5, which is the same as that in tunneling through a single-level QD. In the coexistence regime, the Fano factor is very large [2,37]. This large noise is not hard to understand because oscillator is in two probable states of tunneling and shuttling at the same time. The two probable states increase the uncertainty of electron transport. In the shuttling regime, the oscillator is in the state of oscillation with a certain amplitude. This mechanical motion leads to shuttling effect and enables electron transport through the QD in a definitive manner. The shuttle mechanism "catches" an electron from one lead and transports it to the other lead. This process reduces the role of electron motion and noise. As a result, the Fano factor is much less than unity in the shuttling regime.
Now we consider the current noise in the low bias regime. Based on Eq. (15), Fig. 13 shows the Fano factor spectrum for the dissipation rate that satisfies γ < Γ l , Γ r . At the bias voltage eV = ω 0 , the channel of ground state is primarily open, and the probability of electrons passing through the system is concentrated in this chan-

108501-12
Wenxi Lai, Chao Zhang, and Zhongshui Ma, Front. Phys. 10, 108501 (2015) nel. Hence, the zero-frequency Fano factor appears to be sub-Poissonian, as illustrated in Fig. 13(a). This behavior is similar to a bare tunneling process; however, noise peaks appear at frequencies ±ω 0 . The presence of noise indicates the small probability of an additional channel contributing to the transfer. Figure 13(b) reveals that, at voltage eV = 2 ω 0 , the Fano factor appears to be super-Poissonian at frequencies that are integrals of ω 0 . At this voltage, the excited states of the mechanical oscillator begin to effectively contribute to the electron flow, which destroys the zero frequency sub-Poissonian statistics. This effect can be observed more clearly at higher voltages [see Figs. 13(c) and (d)] because more channels are open for transport. The super-Poissonian noise at zero frequency implies that the electron transfer through the ground state channel is interrupted by tunneling through the channels of the excited levels.
Due to the Coulomb blockade effect, there will be a competition between these channels with different tunneling probabilities, and correlation occurs among the transport processes of multiple channels. In the area of off-resonant frequencies, noise suppression appears in the sub-Poissonian statistics. This result was previously achieved using incoherent dynamics [38]. However, the previous system is a simple vibrational QD, not a real shuttle-junction, because the tunneling amplitudes are constant. In Fig. 14, Fano factor is shown for the situation γ Γ l , Γ r . Due to the fast damping effect, the contribution from the excited states of the mechanical resonator is very small. An electron transports with the dominant probability through the channel provided by the ground state of the system, causing suppression of the noise. When the bias voltage is increased to a finite quantity, the zero-frequency Fano factor is sub-Poissonian and approaches 0.5. In fact, the result is also true in the large bias voltage limit [37].
The above interpretation of the physical picture of the super-Poissonian statistics is consistent with the previous results. In a movable QD array [116], different current channels are formed due to different resonant quantum states connecting the neighboring dots in the co-tunneling regime. Switching between those channels gives rise to super-Poissonian noise in the small damping rate regime. In the semiclassical case, electron transport through the bistable coexistent shuttling and tunneling channels causes a super-Poissonian noise spectrum both at zero [37] and finite frequencies [122]. The bistability of the quantum shuttle is further illustrated with full counting statistics [123].
The relationship between the resonator state and charge transport is an interesting issue. By measuring the character of charge transport in the system, we expect to obtain the information about mechanical resonator. In Section 3, we showed that the current is sensitive to the mean phonon number of the resonator with varying bias voltage. Moreover, we demonstrate that the Fano factor spectrum of the charge transfer is dependent on the mechanical oscillator motion. Recently, the current noise together with the phonon statistics has been considered [124]. As illustrated in Ref. [124], a uniform relation of statistical characteristics between the localized phonon and electron current does not exist. The relation is determined by the system parameters. As an example, in the parameter regime where the charge-oscillator coupling is weak and two tunneling barriers are very asymmetric, we can expect a super-Poissonian current Fano factor associated with sub-Poissonian phononic population. The sub-Poissonian statistics of phonon distribution is also predicted in a system with a resonator coupled to a superconducting single-electron transistor [125]. Rodrigues et al. found that the system behaves as a micromaser and can generate number-squeezed resonator states. Because the phonon number distribution would be narrow in the squeezed state [26], the phonon noise is reduced to sub-Poissonian noise.

Shot noise and temperature
Next, we consider the effects related to the temperature of the environment around the system. As illustrated in Fig. 15, a critical low temperature is required to observe the quantized current. If we increase temperature, the step-like current structure becomes less clear. In particular, the plateaus disappear in the temperature range k B T > 0.3 ω 0 . This temperature dependent effect agrees with previous experimental [71] and theoretical [126,127] studies. For temperature k B T < 0.3 ω 0 at the molecular vibrational junction, the density matrix contains discrete peaks, and the quantized levels of the mechanical oscillation become important [128]. We now connect the temperature effects to the behavior of the zero-frequency current noise.
From Fig. 16, we know that due to high temperature, the Fermi surface in the electron leads span more than one level of the mechanical oscillator, which weakens the role of discrete levels in the junction, causing the current steps to disappear. The step disappearance can also be observed when the frequency-independent quality factor of the vibrational mode is very small [31]. Now, let us briefly discuss the influence of temperature on the zero-frequency current fluctuation. As illustrated in Fig. 17, the large Fano factor is predicted in the ω 0 /(k B T ) ∼ 1 case. The curves at different bias voltages have a common feature in that they do not obviously depend on temperature until ω 0 /(k B T ) decreases to

108501-14
Wenxi Lai, Chao Zhang, and Zhongshui Ma, Front. Phys. 10, 108501 (2015) about 3. When ω 0 /(k B T ) 3, the noise is dominated by thermal noise. We connect the Fano factor with the current-voltage curves at different temperatures in Fig. 15. The noise increase with temperature is accompanied by the disappearance of the current steps near ω 0 /(k B T ) = 3, which implies that the thermal noise that emerges due to the finite temperature removes the quantum mechanical characteristics of the current.

Vibrational excitation-induced dephasing
The electron coherence in the shuttle-junction has not been intensively studied. This section investigates decoherence of electrons induced by the electromechanical vibration in the single-molecular transistor. The decoherence is investigated by embedding a harmonically movable QD in one (target) arm of the AB interferometer and locating a fixed QD in the other (reference) arm.

Model of the AB interferometer
The schematic structure of our AB interferometer is illustrated in Fig. 18. It contains two single-level QDs coupled to two electronic leads in parallel. One QD with the energy ε 1 (QD1) is fixed in the upper arm and the other QD with energy ε 2 (QD2) is located in the lower arm. The two arms and the electrodes enclose a magnetic flux Φ, which passes through the loop-plane. Here, QD2 is assumed to be bounded in a harmonic potential, which consists of an electromechanical shuttle-junction, while QD1 provides reference path. We consider both the inter-dot and intra-dot Coulomb blockade limits in order to verify that electrons propagate through the two-path interferometer one-by-one. The spin degree of freedom is not involved in our approach. The Hamiltonian can be written in the form of [35] where Fig. 18 Two single-level quantum dots connect to two leads in parallel, which enclose a magnetic flux Φ for Aharonov-Bohm interference. The upper dot is fixed and the lower dot is bounded by an harmonic potential and movable in the horizontal direction. (33) describes the non-interacting electrons in the left (y = l) and right (y = r) leads. Here, d † yk and d yk are creation and annihilation operators of electrons with momentum k and energy ξ yk , respectively. In the following Hamiltonian, c † i (c i ) represents the creation (annihilation) operator of QDi (i=1,2): The last term is the contribution from the work of an electric field on charged QD2. The position x of the center of mass of the vibrational QD is defined in the same way as that in Section 2. Moreover, V denotes the bias voltage, d is the effective distance between the two electrodes, e is the absolute value of the electron charge, and x 0 is the zero-point position uncertainty /2mω 0 of the oscillator with frequency ω 0 and effective mass m. The nanomechanical vibration is treated in the quantum regime as where a (a † ) and b k (b † k ) are annihilation (creation) operators for the vibrational mode and its thermal bath, respectively, and ω k denotes the frequency of mode k in thermal bath that is coupled to the oscillator with coefficient g. Tunneling through the two QDs is represented by The tunneling coefficients between the two leads and QD1 are given by T 1y e iSyφ/4 (S l(r) = −1(+1)) and its complex conjugate, where the phase φ is related to the magnetic flux φ = 2πΦ/Φ 0 with the flux quantum Φ 0 = h/e. The coupling coefficient with respect to QD2 is written as T 2y e −iSyφ/4 e Syx/λ . For convenience, we define α = x 0 /λ, which is introduced in an earlier section. The state of total configuration is written by the density matrix ρ T (t), which satisfies the Liouville-von Neumann equation (5). Both the electronic leads and the thermal bath are assumed to be in the equilibrium state at all times and are described by the time-independent equilibrium density matrix ρ L and ρ B , respectively. Assuming that the initial state is ρ T (0) = ρ(0)ρ L ρ B , we can write the state at time t in the form ρ T (t) = ρ(t)ρ L ρ B under the Born approximation. Here, ρ(t) is the density matrix of the system composed of the two QDs and the mechanical oscillator. Using the Hamiltonian (32) and iterating Eq. (5) in the interaction picture to the second order and performing a trace over the leads (tr L ) and the bath (tr B ) variables, we obtain the master equation for the reduced density matrix of the system in the Markov approximation: On the right hand side of Eq. (37), v 0 denotes the evolution term of the system in which QD2 is coupled to the harmonic oscillator. Moreover, v 1 describes the contribution from direct tunneling by QD1 in the absence of QD2, v 2 is the right-hand side of the master equation in our previous work [34], representing the contribution from vibration assisted transfer through QD2 alone, v 12 is coherent term of the transport involving the two dots, and v d accounts for the dissipation of the vibrational mode. Explicit expressions of these terms are given by and The degrees of freedom in the electronic leads and the thermal bath are assumed to be continuous with densities of states N y (ξ yk ) and D(ω k ), respectively. The coefficients in the above equations, which correspond to particle hopping into or out of the system, are composed of integrals over these reservoir variables via Here, Γ ijy (ξ yk ) = 2πN y (ξ yk )T * iy T jy and γ(ω k ) = 2πD(ω k )g 2 for i, j = 1, 2. We have the Fermi-Dirac distribution function in lead y of f y (ξ yk ) = [e (ξ yk −μy)/(kBT ) + 1] −1 and the Bose-Einstein distribution function of the thermal bath of n B (ω k ) = e ω k /(kBT ) − 1 −1 , where T is temperature and k B is the Boltzmann constant. In the above equation, we have defined A − mn = c 2 (a † ) m (a) n and A + mn = c † 2 (a † ) m (a) n , where A + mn (A − mn ) describes an electron that hops into (out of) QD2 accompanied by the creation of m phonons and annihilation of n phonons. Moreover, v, v + y = v + (1 + S y )/2, and v − y = v − (1 + S y )/2 indicate the number of electrons accumulated in the right lead. These variables are determined in the following way. The traces tr L (d † rk d rk ρ L ) and tr L (d rk ρ L d † rk ) contain different information about the number of electrons in the right lead when the number is not infinite. Assuming v electrons are in the right lead, then the number of electrons in this lead can be expressed by ρ v f y (ξ rk ) = ρtr L (d † rk d rk ρ L ) and ρ v+1 f y (ξ rk ) = ρtr L (d rk ρ L d † rk ). In the same way, we have ρ v−1 f y (ξ rk ) = ρtr L (d † rk ρ L d rk ). The density matrix satisfies ∞ v=0 ρ v = ρ. The above method is equivalent to the counting approach in the many-body Schrödinger equation [28], representing how many par-

Numerical treatment and current formula
In the following numerical treatment, we consider the wide band approximation and apply the energyindependent transmission rates Γ ijy = 2πN y T * iy T jy (i, j = 1, 2 and y = l, r) and the damping rate γ = 2πDg 2 . We assume Γ 12y and Γ 21y are real and satisfy Γ 12y = Γ 21y = Γ 11y Γ 22y . The chemical potentials for the left and right electrodes are set to be μ l = eV /2 and μ r = −eV /2, respectively. Now, two QDs are considered in the AB ring. Therefore, the Hilbert space of the system is generated by the composite basis {|00 , |01 , |10 , |11 } ⊗ {|0 , |1 , . . . |n , . . .}, where ⊗ is the direct product. The state |ij represents i electrons in QD1 and j electrons in QD2, and |n is the eigenstate of the nth excited level of the mechanical oscillator. In the Hilbert space, the system density matrix element is written as ρ ijkl,mn = m| ⊗ ji|ρ|kl ⊗ |n (i, j, k, l = 0, 1; m, n = 0, 1, 2, . . .). For any two given vibrational states |m , |n , we have 16 density matrix elements ρ ijkl,mn (i, j, k, l = 0, 1) in terms of the electron states. However, six density matrix elements are sufficient to describe the transport process because they constitute a closed equation set for the system dynamics. These matrix elements are ρ ijij,mn and ρ jiij,mn (i, j = 0, 1).
The values of density matrix elements can be obtained by solving Eq. (37) under the condition (24). Similar to Eq. (25), we project the stationary term of Eq. (37) in the basis of the Hilbert space as and where i, j = 0, 1, excluding the case i = j = 1. Then, we can obtain a set of 5(N + 1) 2 linear equations, where N = 0, 1, 2, . . . is the number of excited vibrational levels. These equations can be solved with the use of the normalization condition (26). For the numerical treatment, we take N = 18. This approximation is valid for the low bias voltage, weak dot-lead couplings, and finite oscillator damping rate applied here because the contribution from the higher levels (n > N) of the vibrational mode is very small. For the case of strong inter-dot Coulomb interaction, we assume that the state of two-electron occupation is not inside the transport window. In other words, the bias voltage is so low that only one electron passes through the system at any time. As a consequence, the process involving the state ρ 1111,mn is not contained in our equations [28]. Substituting Eq. (37) into Eq. (27), we reach the following expression for the current: where and Here, I 1 is the current through QD1 alone, and I 2 is the current across the electromechanical junction in the absence of the reference arm. In fact, this is the same as the current directly derived from the master equation of the single-molecular junction [34]. The interference in terms of the off-diagonal density matrix for the electronic states is given by I 12 .

Phase relaxation
Due to interference of two waves, the visibility can be reduced not only by phase destruction of the waves but also by the difference in the amplitudes of the absolute values. The electromechanical systems substantially enhance the electron transport for certain bias voltages [29,30,69,82]. Therefore, to understand the net contribution of phase relaxation to the interference fringe, we balance the amplitudes of waves in the two paths by setting the bare transmission rates of QD2 to be less than those of QD1. To this end, the bare tunneling rates for the reference path are set to be Γ 11l = Γ 11r = 0.01ω and those for the target path taken as Γ 22l = Γ 22r = 0.001312ω 0 . Then, the current in equation nearly equals that in equation with the small difference I 1 − I 2 < 10 −5 /(eω 0 ). In this case, we suppose that the absolute values of the two amplitudes corresponding to the two paths is almost the same. The current, which plotted with respect to the magnetic flux, is indicated by the solid line in Fig. 19. The current shows AB oscillation with a period of 2π. The interference current does not vanish at its weakest points ((2n + 1)π, where n is an integer). Moreover, these results reveal that the coherence of the electron wave is influenced by electromechanical vibration. Using the same method for balancing the current amplitudes in the two paths, we provide three additional examples in Fig. 19 for different parameters. The low bias voltage eV = 3 ω 0 (red dotted line), high damping rate γ = 0.3ω 0 (green dashed line), and small tunneling length α = x 0 /λ = 0.3 (blue dot-dashed line) weaken the effects from vibrational mode. As a result, the interference fluctuation is enhanced. Figure 19 shows that there is no noticeable shift in the minimum and maximum values of the interference pattern under different parameters. Using this property, the current visibility can be easily calculated using I max I(φ = 0) and I min I(φ = π). This substitution works under the conditions of ε 1 = ε 2 = 0 and Ω ω 0 . The visibility of interference fringe is given by the formula V isibility = (I max − I min )/(I max + I min ). In Fig. 20(a), the substantial influence of the bias voltage to the interference visibility can be observed. For very low voltages eV < 2 ω 0 , there is no exited level contained in the transport window (eV ), and the visibility is close to unity. When the applied voltage is close to zero, the corresponding current approaches zero, causing a small drop of visibility near the zero voltage. The excited states of the mechanical oscillator play an important role in the phase relaxation of electrons. By increasing the bias voltage, the excited levels of the vibrational mode are involved in the transport, which suppresses the visibility.
In the low voltage area regime, a few discrete states of the vibration contribute to the transport, and the visibility displays a step profile. The mechanical oscillation is naturally coupled to the thermal bath, and it has an intrinsic lifetime that is the inverse of the damping rate γ. By increasing the damping rate, the visibility increase can be observed, as shown in Fig. 20(b), because the contribution from the mechanical motion would decrease in the case with a high damping rate. The visibility is no longer obviously enhanced for damping rate γ > 0.2ω 0 . These results agree with the transition of the electromechanical system from the so-called shuttling regime into the tunneling regime [36,37]. The visibility is not very high, even at the quality factor ω 0 /γ < 1, implying that the coherence of electron does not obviously depend on the intrinsic lift time of the mechanical oscillator for a large damping rate. In fact, the interference pattern is affected by the strength of electron-phonon interaction, which is determined by the parameter α = x 0 /λ. In Fig.  20(c), the visibility is plotted with respect to the coupling strength α. For a given oscillator with mass m and frequency ω 0 , the zero-point uncertainty x 0 is fixed, and the coupling strength is primarily related to the tunneling length λ. For an infinitely large tunneling length λ → ∞, we have α → 0. In this case, 2 in Eq. (37) is close to the form of 1 , and the effect of vibration in 2 and 12 vanishes. As a consequence, tunneling between the two electrodes and QD2 is almost independent of the dot displacement. Therefore, we obtain the visibility close to unity, as illustrated in Fig. 20(c). In general, the large current induced by the vibrational junction is the reason for the reduced visibility in the AB ring. For instance, when we takes the same bare tunneling rates for the two paths shown in Fig.  20, the probability of an electron passing the target arm is much larger than that of an electron propagating through the reference arm. There is another probable reason for the weak interference, namely the phase shift of electron waves. The propagation of an electron wave through the vibrational junction gives rise to many scattered excited states. In fact, at a sufficiently high applied voltage, the electron is in superposition of a large number of single-particle excited modes associated with the electron-phonon interaction. These excited states are characterized by phase acquirement related to the absorption and emission of phonons. We expect that the scattering in the space of positive phase shifts is symmetric with the space of negative phase shifts. As a consequence, the interference of all scattering waves does not exhibit a global phase shift between the two paths. This property is valid for the same dot levels, ε 1 = ε 2 and for the weak charge-field coupling Ω . In the next section, we discuss the case where ε 1 = ε 2 . The influence of the charge-field coupling strength to the electron coherence has been previously considered in a similar system [105,106].

Coherent phase shift
As we mentioned earlier, there is no global phase shift when an election propagates through the single-dot electromechanical system. However, it does not means there is no phase shift when one component of the electron wave transports through any individual level of the system. To observe the phase change of the propagating electron wave through the target system, we change the gate voltage in the reference arm to observe the variance of the AB interference oscillation. As illustrated in Fig.  21, the pattern of the interference oscillation shifts continuously in one direction when the resonant level of the reference arm is moving. The phase shift breaks the original symmetry of the interference fringe under φ ↔ −φ. This symmetry breaking is induced by detuning between the two QDs. In the AB interferometer of electron transport, the interference is strong only when the energy level in one path is close to that in another path [129]. In other words, the propagating waves in both paths must be oscillating in (at least nearly) the same frequency. Although QD1 is detuned from the electronic level of the molecular junction, the molecular system still has energy levels provided by the mechanical oscillator. Therefore, interference does not disappear in the case of the detuning, except there is some phase shift.
The propagating wave in the reference arm only interferes with the wave in the molecular junction whose resonant energy (ε 2 +Δε) is the same as the resonant energy (ε 1 ) of the reference arm. Here, Δε is defined as the  acquired or lost by an electron due to its inelastic scattering on the vibrating QD. In Fig. 21, the phase shift corresponding to the detuning ε 1 −ε 2 represents the phase change of the sub-transmission amplitude whose energy is ε 1 = ε 2 + Δε in the molecular junction. The total amplitude of electronic wave transferring through the molecular junction is, of course, the superposition of all sub-transmission amplitudes with different resonant energies.
By choosing the particular detunings ε 1 − ε 2 in Fig.  22, we analyze the quantitative phase shift, especially for the resonant levels. Without loss of generality, the zeropoint energy is taken at ε 2 = 0. The numerical results in Fig. 22 show the the phase shift Δθ of the transmission amplitude with energy ε 2 + Δε in the molecular junction roughly satisfies the relation We can write this relation in more compact form as Equation (50) connects the phase shift of an electron to its energy variance during the transport. When an electron gains or loses an integer number of phonons, its phase would be changed by the integer times π, suggesting that the phase difference of propagating waves corresponding to two adjacent vibrational levels is π. Equation (50) is an empirical formula based on the numerical results. This off-phase character is analogous to the phenomenon described by the Friedel sum rule [130][131][132]. The sum rule relates the phase shift of a scattering electron to the number of states in the energy interval due to scattering [133][134][135][136][137][138][139]. However, the electron number accumulated in the impurity, which is described by the general Friedel sum rule, is replaced by the phonon number involved in the electron transfer in our present system. The phase shift is very steep when the energy of incident electron sweeps over the resonant levels [140]. However, in our model, the tunneling is dependent of the position of the mechanical oscillator. Thus, the phase change is continuous and very smooth. The positiondependent tunneling causes an inelastic process, which improves the decay of the oscillating QD and broadens the energy levels of the system [141].
From Fig. 19, we know that the changes in the applied voltage, the electron-phonon coupling α, and the damping rate of the oscillator do not induce a global phase shift in the AB interferometer. Therefore, the definitive phase relation of π difference between two adjacent levels is independent of these parameters so long as they are properly considered such that the discrete levels of the mechanical vibration effectively contributes to the electron transport. For instance, on one hand, the applied voltage should be large enough that at least one excited level of the oscillator is included in the transport window. On the other hand, the voltage is not too large that the feature of discrete levels involved in the tunneling is not altered. In fact, the phase shift is related to the unit quanta of the mechanical oscillator as shown in Fig. 22.
According to the above analysis, the neighboring resonant levels in the molecular vibrational junction are offphase by π, which is the character of the one-dimensional quantum system. In the one-dimensional system, the upper energy level has one more wave function node than the lower level, and each node changes the phase of the transmission amplitude by π. This property may not be true if the system is not strictly one-dimensional [142]. In the AB interferometer experiment, where a fixed QD is embedded in one of the arms, the phase behaviors are the same for all resonant levels of the QD [88][89][90]. This is different from the present effect found in our electromechanical system, where all vibrational levels are coherently correlated with a definitive phase difference of π. The phase shift varies from zero to any large value, depending on the net number of phonons involved in an electron tunneling.
The reason for the decrease in visibility becomes more clear. In fact, an electron considers all channels of the discrete vibrational levels involved in the transport process. Therefore, interference not only occurs between the propagating waves in the two paths but also occurs among waves taking different channels of the vibrational junction. As we analyzed above, any two neighboring channels have a phase difference of π. The wave functions taking different vibrational levels interfere destructively because of the phase differences. This is the reason for the phase relaxation in the AB interference due to the vibrational junction (see Fig. 19). In a double-QD two-electron AB interference, two components of conductance oscillations with the same amplitude are canceled due to their phase difference of π [143]. Because two components of the conductance have the same amplitude, the final conductance disappears. In the present case, the electron occupation probabilities of the resonant levels of the electromechanical system are not the same. Therefore, there is a net current in the system, but it is not fully coherent. In fact, the interference between the different channels is also reflected in the direct transmission of charge through the electromechanical system [34]. The current calculated from the scheme considering both diagonal and off-diagonal couplings between the vibrational mode and the electron tunneling is remarkably lower than that obtained by the approach in which only diagonal terms are taken into account. This current suppression is related to the destructive interference between different transport channels.

Summary
The general master equation describing a single molecular shuttle-junction can be derived in the Born-Markovian approximation, including the position dependence of the tunneling rates and the Fermi distribution functions of the electronic leads. The equation of motion is numerically solved in the Hilbert space formed by the electron and phonon Fock state. In the high applied voltage limit, the mechanical oscillation of the island in the junction is semiclassical, which enables remarkable shuttling phenomena. However, in the low bias voltage regime, the contributions from discrete levels of the oscillator cannot be ignored. These contributions lead to current steps, super Poissonian noise, and phase shift of electrons, which is similar to the Friedel sum rule. The current steps, which are caused by quantized channels in the junction, tend to disappear when the temperature of the electron leads increases. The excited energy levels in the shuttle-junction lead to super Poissonian noise because these levels provide multiple paths to the electron and increase its uncertainty during the propagation. The general master equation introduced here suggests that the off-diagonal terms of the density matrix provide an important contribution to current. Electrons propagating through the single-molecular vibrational junction are dephased due to electron scattering by the excited levels of the vibrational mode. The transmission amplitudes corresponding to channels of the vibrational resonant levels are coherently correlated via neighboring channels and have a definitive phase difference of π. Because of the phase shifts between the resonant levels in the electromechanical junction, different branches of the transmission waves destructively interfere with each other. As a consequence, the electron tunneling through the system is not fully coherent. Over a wide range, the phase difference of π is independent of the bias voltage, tunneling length, and lifetime of the vibrational mode. The phase difference only depends on the frequency of the mechanical oscillator. When tunneling amplitudes are dependent on the displacement of the oscillator, the mathematical treatment becomes more complex. As we know, the master equation is primarily used to describe the shuttle-junction. The construction of a Green's function approach for such systems is a considerable future work. Furthermore, many phenomena such as the temperature of the junction, cooling, electromagnetic effects, and light spectrums are studied in molecular junctions in which the tunneling amplitudes are independent of the oscillator displacement. These problems have not yet been considered in a molecular shuttle-junction. Further research on the phase change and electron coherence in such systems is also required because the molecular shuttle-junction is related to nonlinear electron-phonon couplings and is important for the electron transport in real materials. 1 m 1 !n 1 ! ×[−e Syα(a † +a) (S y αa † ) m1 (S y αa) n1 cc † ρ w,v (t)f y (ξ yk )e i(ξ yk − 0−(m1−n1) ω)τ / −e Syα(a † +a) (S y αa † ) m1 (S y αa) n1 c † cρ w,v (t)(1 − f y (ξ yk ))e −i(ξ yk − 0+(m1−n1) ω)τ / +e Syα(a † +a) cρ w + y ,v (t)c † (S y αa † ) m1 (S y αa) n1 (1 − f y (ξ yk ))e i(ξ yk − 0−(m1 −n1) ω)τ / +e Syα(a † +a) c † ρ w − y ,v (t)c(S y αa † ) m1 (S y αa) n1 f y (ξ yk )e −i(ξ yk − 0+(m1−n1 ) ω)τ / +(S y αa † ) m1 (S y αa) n1 cρ w + y ,v (t)c † e Syα(a † +a) (1 − f y (ξ yk ))e −i(ξ yk − 0+(m1 −n1) ω)τ / where S y = −1 when (y = l), S y = +1 when (y = r), w + y = w +(S y −1)/2, and w − y = w −(S y −1)/2. Equation (20) can be obtained by performing the time integration and summing over wave vector k in the above equation.

Appendix 2: Projection of the master equation in Fock states
In Eq. (25), we encounter compound operators that will be projected. In this section, all possible projections are given, which allow us to arrive at the density matrix elements. The density matrix has been projected for the electron occupation states ρ 00 (t) = 0|ρ(t)|0 and ρ 11 (t) = 1|ρ(t)|1 . Consequently, the density matrix in the states of mechanical oscillator can be obtained. The compound operators in Eq. (25) are projected as follows: