Current-induced rotations of molecular gears

Downsizing of devices opens the question of how to tune not only their electronic properties, but also of how to influence ‘mechanical’ degrees of freedom such as translational and rotational motions. Experimentally, this has been meanwhile demonstrated by manipulating individual molecules with e.g. current pulses from a Scanning Tunneling Microscope tip. Here, we propose a rotational version of the well-known Anderson-Holstein model to address the coupling between collective rotational variables and the molecular electronic system with the goal of exploring conditions for unidirectional rotation. Our approach is based on a quantum–classical description leading to effective Langevin equations for the mechanical degrees of freedom of the molecular rotor. By introducing a time-dependent gate to mimic the influence of current pulses on the molecule, we show that unidirectional rotations can be achieved by fine tuning the time-dependence of the gate as well as by changing the relative position of the potential energy surfaces involved in the rotational process.


Introduction
Gaining control over molecular-scale collective mechanical degrees of freedom, such as translations and rotations, poses a big challenge to current state of the art nanoscale manipulation techniques. It therefore represents a major advance in the field that individual molecular rotors as well as collective rotations in molecular assemblies have been experimentally demonstrated [1][2][3][4][5][6][7][8]. The propagation of angular momentum in such assemblies may open the door, e.g. to the implementation of molecular scale analog computing devices, such as the Pascaline or more recent mechanical computers [9].
A crucial condition for realizing single-molecule gears is the ability to induce unidirectional rotations, which can be achieved, e.g. by using current pulses [6,[10][11][12], voltage pulses [13,14] or mechanical way [15,16] in a scanning tunneling microscope (STM) [17][18][19]. Meanwhile, diverse milestones have been achieved such as stepby-step molecular rotation [20], controlling the rotational direction of a molecular rotor [21], and collective rotation effects [22,23], among others [24][25][26][27]. However, the underlying physical mechanisms leading to unidirectional molecular rotation are not well understood, since they involve in general terms a delicate interplay between collective mechanical and electronic degrees of freedom. On the theoretical side, a combination of model-based approaches [28][29][30][31][32][33][34], catching the basic physics of the problem with more advanced first-principles methodologies able to provide atomistic, system-specific information is required. Some of the problems here include, e.g. the computation of the potential energy surface(s) (PES) necessary to describe the excitation of molecular motion and the definition of one or more collective degrees of freedomthrough an appropriate coarse-graining procedure-to describe the rotational dynamics [35].
In this study, we approach the problem of unidirectional rotation getting inspiration from the well-known Anderson-Holstein model (AHM) [36][37][38], which describes an electronic system linearly coupled to a harmonic optical phonon mode. In general terms, changes in the occupation of the relevant electronic states lead in the AHM to a linear shift of the equilibrium position of the vibrational mode PES. We aim at extending the AHM The electronic part of the Hamiltonian is written, within a single particle picture, using molecular orbitals and is thus diagonal in this basis.  { } n are the corresponding molecular orbital energies and the operator † d n (d n ) creates (annihilates) an electron in the nth molecular orbital (including spin degrees of freedom). These orbitals in general depend on the individual coordinates of the constitutive N a -atoms a , meaning that intramolecular distortions can modify the orbital energies. The nuclear part is treated classically in the spirit of the Born-Oppenheimer approximation by invoking the large mass difference between electrons and nuclei [46]. The variables P i and M i denote the momentum and mass of the ith nucleus, respectively. In order to address the rotational dynamics of a molecule deposited on a substrate, we need to compute its potential energy surface in terms of the full set of nuclear coordinates a and include the influence of the surface; this requires, however, an atomistic approach to the problem, which goes beyond the scope of the current study. To further proceed, we assume that the substrate is acting only as a conformational constraint, so that the dynamics of the molecule on the substrate can be described in terms of only three collective variables: the center of mass coordinates (implying that the internal relative motion is neglected) and two collective angular variables. Thus, we write the Hamiltonian as: where R CM and M mol are the center of mass coordinate of the molecule and the total molecular mass, respectively, θ is the azimuthal angle in a plane parallel to the substrate, and ψ is a polar angle or tilt angle with respect to the normal to the substrate (see figure 2). Note that these angles can always be clearly defined with respect to the specific orientation of the molecule in the state with lowest energy. The vectors L is angular momentum and I n(θ, ψ) is the moment of inertia with respect to the principle axis q y ( ) n , . Clearly, a more systematic approach would imply an explicit coarse-graining of the full set of molecular coordinates {R i } (N a degrees of freedom) down to a set of few collective variables, going beyond q y R , , CM [35]. Here, however, we assume that this has been already carried out and based our choice of the collective coordinates on physical intuition. To further simplify our model, we consider situations where the molecule can not be tilted with respect to the normal to the substrate, thus removing the angular variable ψ from our description. If the center of mass is also fixed by the interaction with the substrate, we are then left with a single angular degree of freedom θ. We can separate the electronic Hamiltonian into a contribution arising from occupied states -up to the highest-occupied molecular orbital (HOMO)-and a contribution from the unoccupied states -beginning with the lowest-unoccupied molecular orbital (LUMO): In general, only orbitals below the LUMO are filled up. Therefore, the electronic operators are acting on a subspace of the Fock space in which all the occupation numbers below HOMO are equal to 1. As a result, the first term in the previous equation becomes a scalar function of the collective variables q y , which defines the ground state potential energy surface of the molecule. If now additional electrons are added to the molecule, the Hamiltonian in this subspace becomes: where U n =ò n +U 0 represents now the potential energy surface with an occupied nth-orbital.
Once this minimal molecular Hamiltonian has been introduced, we can easily extend it to include the coupling to electronic degrees of freedom describing the electronic systems of the substrate and the STM tip: The molecule is contacted by the STM tip (α=T) and deposited on the substrate (α=S) with the energy spectra ò kα , where k stand for the corresponding wave vector. The electronic matrix element a T k n describes the coupling between the levels in the molecule and the reservoirs. The operators a † c k and c kα are creation and annihilation operators of an electron in level ò kα , respectively. The third term in equation (6) mimics the local electrostatic gating effects coming from the action of the STM tip on the molecule such as consecutive probing; the overall effect is included in a time-dependent gating Δ n (t).

Rotational Anderson-Holstein model (RAHM)
Combining equations (3) and (6), we arrive at the following general Hamiltonian: nuc . Now,H e represents the effective electronic Hamiltonian relevant for the dynamics andH m is the unperturbed molecular Hamiltonian with effective ground state potential energy surfaceŨ 0 . The interaction term -V e m , which is crucial in our approach, describes the possibility that the molecular conformational state can change from its ground state PESŨ 0 to any other PESŨ n in dependence on the occupation of electronic states with energies    n LUMO . Note that the interaction -Ũ U n 0 is in general a non-linear function in the angular variables θ and Ψ. If the angular distortions and the separation of the two PES minima can be considered as small, then a Taylor expansion can be used to obtain a linear coupling between the angles and the electronic degrees of freedom. In this case, we recover the AHM. However, in our case, where the focus is the possibility of inducing global rotations of the molecule, both the full non-linear potential and the non-linear electron-rotation coupling must be taken into account. Therefore, our model can be viewed as a generalized version of the AHM. For the sake of simplicity, we will drop all the tildes from now on.
To extend the model given by equation (7), one can include the first-order corrections due to the fast nuclear dynamics in the Hamiltonian, which yields off-diagonal coupling terms. In general, one may also introduce an angle dependence of the molecule-lead coupling strength. This dependence will be determined by the details of the setup, e.g. the geometry of molecule and the tip.

Langevin equation
Using Hamilton's equations of motion, we can derive a Langevin equation for the rotational dynamics of the classical variable θ as: where s denotes the reduced electronic density matrix. We have further introduced a matrix U with elements given by U mn =U n δ mn . In this equation of motion, the first term on the right hand side will be denoted as current-induced torque. It is given by the expectation value of an operator valued torque. The second term, x , which is a stochastic operator, quantifies the deviation of the torque from the mean value. If one assumes that the time scales for the rotational dynamics of the molecule and for the electron transfer from the tip to the molecule are well separated, then the electronic system can always reach a stationary state according to the corresponding molecular configuration. This represents the so called non-equilibrium Born-Oppenheimer (NEBO) approximation [40], and one can show that the noise term in the adiabatic limit is always delta-correlated in time. In order to account for the noise, the operator x is often replaced by a classical stochastic torque ξ(t) with an appropriate correlation function [40,42,43,45]. Note that the damping is implicitly hidden in the first term. One can perform the adiabatic expansion of the density matrix up to first-order to get an explicit expression for the damping, which can be shown to fulfill the fluctuation-dissipation theorem in the limiting case of thermal equilibrium. On the other hand, according to the equation (8) the dynamics of the relevant nuclear degrees of freedom are mainly determined by an ensemble-averaged PES given on the right-hand side if the torque noise is sufficiently weak. The opposite limit is captured by considering a single-trajectory dynamics, where the switching between potential energy surfaces is purely stochastic [35].

Electronic dynamics
To actually solve equation (8), we first need to know the reduced density matrix, which can be obtained by solving the following equation of motion in the time domain [47,48]: . 9 Here, ( ) H t is the matrix with and P a are so-called current matrices, which can be expressed in terms of Green's functions and self-energies as: Here  G and  S are lesser/greater Green's functions and self-energies, respectively [49]. The current matrices are closely related to the current by: t  e  t  2 Re Tr  ,  11 with e denoting the elementary charge. However, in the general case the calculation of the current matrices is very challenging. To overcome this, we consider the wide-band limit [50], where the real part of the retarded selfenergy S R is vanishing and the imaginary part is a constant independent of the energy. Then, we can solve the following system of differential equations to get the reduced density matrix without evaluating the complicated convolution integral in equation (10): where G a denotes the broadening matrix with G G S = å =a a 2 Im R , and P a ( ) t p are auxiliary current matrices, which are related to the current matrix through the relation: The coefficients R p and poles c a + p are determined by the Matsubara expansion [51] with coefficients R p =1 and poles c m where μ α is the chemical potential in reservoir α and 1/β=k B T represents the thermal energy. For a better convergence, we choose Karrasch's approach [52] to obtain the coefficients R p and poles z p . To get additional insight into the problem, we perform now an adiabatic expansion of the reduced density matrix and the current matrices The superscript (n) denotes the n-th order correction with respect to the time derivative. The zeroth-order term is the instantaneous solution of the time-dependent Hamiltonian and the first-order term will be related below to the damping and the influence of external driving. We therefore limit or discussion to these two contributions in the expansion and provide analytic expressions for them in the appendix A. A similar expansion can be carried out for the current J α (t), giving:

A single electronic level coupled to the rotational motion
We apply the formalism presented above to the simple case of a planar molecule with N-fold symmetry placed on top of a metallic substrate and we consider only one relevant electronic level (LUMO) as relevant for the rotational dynamics. In addition, we focus on the current-induced torque only [54], so that the stochastic term is not considered (Ehrenfest approximation). Since the ground state PES should display the molecular symmetry, we make the Ansatz: where τ 0 gives the amplitude of the angle-dependent potential. We also need to define the PES U 1 (θ) corresponding to the excited state with non-zero electron occupancy. The simplest way to describe it is by introducing a phase shift between the two PES [35]: Once we specify the molecular potential, we can write down the adiabatic corrections to the reduced density matrix by using the adiabatic expansion; details of the calculations can be found in appendix B. Inserting the obtained reduced density matrix into equation (8), we arrive at an equation of motion for the rotational degree of freedom of the molecule: The three terms on the right-hand side are the current-induced mean torque, current-induced damping, and the external-driving torque, respectively, which are given by the expressions: c o s c o s , .  (20); the solutions will be divided into two cases, depending whether the orbital shift is time-dependent or not.

Time-independent case
In this case, we consider a planar molecule with a three-fold symmetry (N=3), which experiences a constant orbital shift Δ(t)=const. One can immediately find that the external-driving torque vanishes according to equation (B.2). To solve the equation of motion, we use the following parameters which are representative of typical experimental conditions [14-16, 35, 55]: T=5 K, f=π/2, I=10 −42 kg·m 2 , Δμ=μ T −μ S = 10 meV, τ 0 =10 meV and Γ=0.1 meV. For the initial condition, we set θ(0)=−π/4 and q = ( ) / 0 0 rad. ns. By solving the adiabatic Langevin equation, we can obtain the phase-space trajectory of the solution, which is shown in figure 3(a). The red (black) line represents the solution with the orbital shift Δ=1 eV (−1 eV) and the solution exhibits a typical damped oscillation behavior with corresponding fixed-points at * * q q q = and L R (close to −π/3 and −π/6 respectively).
To explain this oscillatory behavior, we rewrite equation (20) as follows: . From equation (24), one can understand that the molecule is rotating according to the effective potential with an angle-dependent damping γ(θ), which are shown in figure 3(b). In terms of the effective potential, it is easy to explain why the fix-points are at * * q q q = and L R : when the orbital shift is high (Δ=1 eV), the corresponding electron density on the molecule is nearly zero. According to our previous discussion when introducing the RAHM, we know that in this case U eff (θ)≈U 0 (θ). As a result, the molecule starts to change the orientation until it reaches the nearest local minimum on the effective potential surface, which is simply at * q R (fixed-point on the right). On the contrary, if the orbital shift is low (Δ=−1 eV), the average electron occupation approaches unity, which implies U eff (θ)≈U 1 (θ), so that the nearest local minimum is at * q L in this case.

Time-dependent case
In the previous section, we have already seen that the location of the fixed-points depends on the orbital shift. This suggests that a time-dependent orbital shift may be considered a way to manipulate the location of the fixed-points, i.e. to control the rotational behavior of the molecule. Here, we assume the following orbital shift (see figure 4): where Δis the initial orbital shift; k is the switching rate of the orbital shift and t 1 is the switching time.
In the following calculations, we use the parameters f=−11π/12, k=100 GHz, Δ=1 eV, I=10 −41 kg·m 2 , Δμ=μ L −μ R =10 meV, τ 0 =10 meV, Γ=0.1 eV and t 1 =100 ps [14-16, 35, 55]. For the initial conditions, we choose θ(0)=θ * , σ(0)=σ (0) (θ * ) and , which means that the molecule starts from one of the fixed-points. One can clearly see in figure 4(b) that as Δ(t) decreases the electron occupation increases. The time-dependency of the total current flow into the molecule J L +J R is also consistent with the change of occupation. Since the occupation is changed, the fixed-points are also shifted. On the other hand, the trajectory in figure 4(d) shows a similar oscillatory behavior as in the time-independent case. From the pattern of the phase-space trajectory, it is straightforward to consider a consecutive switching of orbital shift to see whether an open trajectory is possible, which would mean that the molecule is rotating uni-directionally instead of oscillating around a local minimum. To implement consecutive switching, we apply the following time-dependent orbital shift:   Figures (a)-(d) show the time-dependency of orbital shift, electron density, total current flow into the molecule and phase-space trajectory for single-switching with switching timing t 1 =100 ps, respectively. In the bottom, figures (e)-(h) illustrate the same quantities with two switching timing t 1 =100 ps and t 2 =400 ps.
with t 1 =100 ps and t 2 =400 ps. Due to the second switch, the instantaneous change of the fixed-point can possibly make the molecule rotate one-way. To achieve this kind of rotation, the switching timing has to be in coherent with the instantaneous angular velocity; otherwise, the molecule is not able to cross the nearby potential maximum and it will hence display an oscillatory behavior. In figure 4(h), since the second switch is coherent, the trajectory becomes open, which means the rotation is uni-directional. Note that, in our approach, the rotation is driven by manipulating the electrostatic gating effect explicitly instead of using a stochastic driving torque [35].

Response to different phase-shifts
An interesting issue is whether a phase shift f between the two PES may be found, such that the molecule will rotate with the largest angular velocity. To answer this question, we have used the same parameters as for the time-dependent case while adjusting the phase-shift from −π to π to see the behavior of the average angular velocity 6 w . In figure 5, one can clearly see that there exist certain windows in the values of the phase shift where the angular velocity of the molecule is in resonance with the external switching. We therefore denote this behavior as resonant rotation. The response is fully anti-symmetric with respect to f=0, since the ground state PES t q = U sin 3 0 0 is also anti-symmetric. In figure 5 there are peaks close to ±π and ±π/2. The direction of the rotation for f close to −π can be explained by considering the effective potential U 1 . Suppose t<t 1 , then the molecule is always staying on a minimum of the PES. When t 1 <t<t 2 it is moving to a nearby minimum of , which is on the right-hand side of the original location. For t>t 2 , the potential is then switched back to U 0 , but the angular velocity is large enough such that the molecule can further rotate in the same direction, which is similar to the scenario in figure 4(h). For the peaks near f≈−π/2, the mechanism is similar. The only difference is that the second switching takes place before the molecule crosses the first local maximum in the U 0 surface for f≈−π, whereas the second switching happens after the molecule passes the first local maximum in the U 0 surface for f≈−π/2.

Conclusion
We have proposed a rotational version of the Anderson-Holstein model to describe the molecular rotational dynamics in a generic setup consisting of a single molecule adsorbed on a metallic substrate and electrically addressed by an STM tip. By applying an adiabatic expansion to the time-evolution of the reduced density matrix, the damping and external driving can be related to the first-order corrections of the reduced density matrix. To demonstrate our approach, a planar molecule with N-fold spatial symmetry was considered, where a single rotational variable couples non-linearly to an electronic level, the occupation of the latter being controlled Figure 5. Phase response of molecular average angular velocity w with parameters I=10 −41 kg·m 2 , Δμ=μ L −μ R =10 meV, τ 0 =1 meV and Γ=0.1 eV with two switching timing t 1 =100 ps and t 2 =400 ps. by a local, in general time-dependent gate mimicking the influence of current pulses coming from the STM tip. We have shown that unidirectional rotation can be achieved by specific tuning of the time-dependent gate as well as of the relative phase difference of the potential energy surfaces. Our framework can be systematically extended to include multi-level electronic systems as well as more than one collective variable. It thus opens the possibility to make contact with atomistic simulations able to provide quantitative information e.g. on the shape of the potential energy surfaces involved in the rotational process [11] and on the strength of the coupling between mechanical and electronic degrees of freedom. On this basis, a further going step would be the study of the mechanisms leading to the transmission of angular momentum in coupled molecular rotor arrays, which represents a major issue in designing molecule-based mechanical devices.
In order to further simply the problem, we consider the adiabatic approximation, the reduced density matrix can be expanded by the following: Similarly, this also holds for the auxiliary current matrices: As in zeroth-order case, the equations above enable one to get the first-order corrections.