# Emergence of negative viscosities and colored noise under current-driven
Ehrenfest molecular dynamics

Riley J. Preston  College of Science and Engineering, James Cook University,
Townsville, QLD, 4811, Australia  Thomas D. Honeychurch  College of Science
and Engineering, James Cook University, Townsville, QLD, 4811, Australia
Daniel S. Kosov  College of Science and Engineering, James Cook University,
Townsville, QLD, 4811, Australia

###### Abstract

Molecules in molecular junctions are subject to current-induced forces that
can break chemical bonds, induce reactions, destabilize molecular geometry,
and halt the operation of the junction. Theories behind current-driven
molecular dynamics simulations rely on a perturbative time-scale separation
within the system with subsequent use of nonequilibrium Green’s functions
(NEGF) to compute conservative, non-conservative, and stochastic forces
exerted by electrons on nuclear degrees of freedom. We analyze the
effectiveness of this approximation, paying particular attention to the
phenomenon of negative viscosities. The perturbative approximation is directly
compared to the nonequilibrium Ehrenfest approach. We introduce a novel time-
stepping approach to calculate the forces present in the Ehrenfest method via
exact integration of the equations of motion for the nonequilibrium Green’s
functions, which does not necessitate a time-scale separation within the
system and provides an exact description for the corresponding classical
dynamics. We observe that negative viscosities are not artifacts of a
perturbative treatment but also emerge in Ehrenfest dynamics. However, the
effects of negative viscosity have the possibility of being overwhelmed by the
predominantly positive dissipation due to the higher-order forces unaccounted
for by the perturbative approach. Additionally, we assess the validity of the
white-noise approximation for the stochastic forces, finding that it is
justifiable in the presence of a clear time-scale separation and is more
applicable when the current-carrying molecular orbital is moved outside of the
voltage window. Finally, we demonstrate the method for molecular junction
models consisting of one and two classical degrees of freedom.

##  I INTRODUCTION

The prototypical building block of nanoscale electronics technology is the
molecular junction, in which a chosen molecule is bonded between two
conducting leads. Current densities in nanoscale electronic devices far
outweigh that of their macroscopic counterparts[1]. This results in joule
heating[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13] whereby energy from the
tunneling electrons is transferred to the vibrational degrees of freedom of
the molecule. This can result in a plethora of technologically unwanted but
physically interesting effects such as large-scale conformational changes of
the molecule which affect the conduction properties of the system[14, 15], as
well as total bond rupture and device breakdown[16, 17, 18, 19, 20, 21, 3].
Clearly, it is desirable for us to be able to predict and even manipulate such
behaviours to our advantage, emphasizing the need for reliable theoretical
models which accurately describe the interplay between electronic and
vibrational degrees of freedom.

The prevailing theory can be partitioned into two main archetypes; methods
which opt to treat the entire system on a quantum footing, and quasi-classical
methods which utilise a classical description of nuclei interacting with a
quantum environment. Purely quantum models of the molecular dynamics are often
limited to harmonic vibrations[22, 23, 24, 25, 26] or other generic
vibrational potentials[27, 28] which are unaffected by the electronic
environment. Other quantum methods which are able to accurately account for
the large-scale motion required to study device breakdown are generally
limited to highly simplistic systems with minimal relevance to experiment or
are computationally infeasible to simulate[29, 30, 31]. It is instead
convenient to adopt a quasi-classical approach in which the nuclei are assumed
to behave according to classical equations of motion under the effects of the
quantum tunneling electrons. This enables the description of arbitrary
molecular potentials for systems along with being able to capture
experimentally observed phenomena such as current-induced heating[32, 33, 34,
35, 36, 37] along with bond rupture and electronically induced chemical
reactions[38, 39, 40, 41], at the cost of a fully quantum description.

A common accompanying assumption is the clear separation of time-scales
between the classically described nuclei and the electronic environment, which
through the use of a perturbative approximation allows for a Langevin
description of the classical particle[32, 41, 33, 34, 42, 43, 44, 35, 45, 46,
47, 36, 48]. In this regime, the particle dynamics are described by stochastic
differential equations in which nuclei evolve under the influence of a
frictional force which acts to subdue nuclear vibrations, and a stochastic
force which delivers energy to the nuclei; the balance of these two forces
yielding the effective temperature of the molecule. In addition, the influence
of the electronic environment on the molecular potential is captured by an
additional, typically non-conservative[49, 50, 51], renormalization term which
can naturally describe switching phenomena and bond-stability[32, 35, 44, 41].
The Langevin approach has had great success in the description and simulation
of large-scale systems corresponding to real molecular electronic junction
configurations, such as for graphene nanoribbons[52, 53].

One such theoretically predicted phenomenon which leads to device breakdown is
the notion of a negative dissipation; under a Langevin approach, this implies
that the viscosity coefficient becomes negative[32, 44, 36, 47, 54, 55, 56].
In this regime, energy is applied to rather than dissipated from the classical
coordinate until it reaches unsustainable temperatures for the device.
Negative dissipations have also been predicted using purely quantum mechanical
methods in which a population inversion in the quantized phonons leads to
ever-increasing temperatures[46, 57, 58, 59, 60, 61, 62]. The physicality of
such theoretical results have often been called into question, thought to
arise as an artifact of assumptions applied in the theory rather than a
demonstration of real behaviour[63].

In this paper, we extend the quasi-classical formalism by utilising a
numerically exact method to calculate the forces acting on our classical
coordinates without resorting to the assumption of a time-scale separation
within the system. This is commonly known as the Ehrenfest approach to
classical dynamics and it has garnered significant interest in the modelling
of nuclear motion in nanoscale systems[64, 56, 65, 66, 67, 38, 68, 49, 51, 50,
55]. This approach treats the electron-nuclear interaction on a mean-field
level; it therefore does not provide a full description of the inelastic
scattering between electrons and nuclei and cannot fully capture the effects
of joule heating within the system[69, 70, 66, 71, 65]. Efforts have been made
to remedy these caveats via perturbative corrections to the electron-nuclear
correlations with some success[70, 71]. However, in this study we seek only to
utilise Ehrenfest dynamics as a basis for comparison for predictions made by
the perturbative Langevin approach. In doing so, we once again observe the
emergence of negative dissipations and validate the results of the
perturbative Langevin approximation, depending on the parameter range applied.

Our approach additionally allows for the exploration of the autocorrelations
in the stochastic force. In general, the stochastic force is correlated at
different times dependent on the electronic structure of the system
considered; this corresponds to coloured noise[72, 45, 73]. The correlations
are additionally dependent on the non-adiabatic motion of the classical
coordinate. However, it is often computationally infeasible to account for
coloured noise and it becomes necessary to employ a white-noise approximation.
Under this approximation, the coloured-noise diffusion is replaced by a
Markovian, white-noise equivalent which attempts to produce the same dynamical
behaviour[32, 33, 41, 44, 42, 39]. We apply our method to investigate the
effects of a time-scale separation on the diffusion coefficient along with
assessing the validity of the white-noise approximation to the diffusion under
a variety of regimes.

In section II we introduce the theory which describes the classical dynamics.
This involves the introduction of the Ehrenfest force along with its
corresponding perturbative Langevin approximation, both expressed in terms of
non-equilibrium Green’s functions. We also discuss our iterative method for
evaluating the Ehrenfest force through time. In section III.2, we demonstrate
the accuracy of the perturbative approximation for a simple system consisting
of a single classical degree of freedom, while in III.3 we then apply
Ehrenfest dynamics to a system with two classical degrees of freedom. Finally,
in section III.4 we utilise our method to evaluate the diffusion coefficient
and the suitability of the white-noise approximation for the stochastic force
over a range of parameters. We use atomic units for all parameters and results
presented.

##  II THEORY

###  II.1 Hamiltonian

Our system is described by a generic tunnelling Hamiltonian as per

$H(t)={H}_{M}(t)+{H}_{L}+{H}_{R}+{H}_{LM}+{H}_{MR}+H_{\text{cl}}(t).$ (1)

The total system Hamiltonian is partitioned into the following components; the
molecular Hamiltonian ${H}_{M}(t)$, the left and right leads Hamiltonians
${H}_{L}$ and ${H}_{R}$, the leads-molecule coupling Hamiltonians ${H}_{LM}$
and ${H}_{MR}$ which describe the coupling between the electronic states on
the central molecule and the left and right leads, respectively, and the
classical Hamiltonian $H_{\text{cl}}$ which describes the time-evolving
molecular geometry.

The molecular Hamiltonian takes the form:

${H}_{M}(t)=\sum_{ij}h_{ij}(\mathbf{x}(t))d^{{\dagger}}_{i}d_{j},$ (2)

where the operators $d^{{\dagger}}_{i}$ and $d_{j}$ denote the creation and
annihilation operators for the atomic electronic states whose energies are
given by the diagonal Hamiltonian matrix elements $h_{ii}(\mathbf{x}(t))$
while the hopping amplitudes are given by $h_{ij}(\mathbf{x}(t))$. Note the
explicit time dependence here, which arises as a result of the time evolution
of the multi-dimensional classical coordinate $\mathbf{x}$. We will usually
not show the time dependence of $\mathbf{x}$ explicitly.

The leads Hamiltonian is taken in the standard form:

${H}_{L}+{H}_{R}=\sum_{k\alpha}\epsilon_{k\alpha}d^{\dagger}_{k\alpha}d_{k\alpha},$
(3)

where the creation and annihilation operators are given by
$d^{\dagger}_{k\alpha}$ and $d_{k\alpha}$, and the subscript $k\alpha$ denotes
the operator acting on state $k$ in the $\alpha$ lead which has energy
$\epsilon_{k\alpha}$.

The system-lead coupling Hamiltonians $H_{LM}$ and $H_{MR}$ are given by:

${H}_{LM}+{H}_{MR}=\sum_{k\alpha i}\Big{(}t_{k\alpha
i}d^{\dagger}_{k\alpha}d_{i}+\text{h.c.}\Big{)}.$ (4)

The matrix elements $t_{k\alpha i}$ (and their conjugates) describe the
tunnelling amplitudes between lead states $k\alpha$ and the molecular orbitals
$i$. We take $t_{k\alpha i}$ to be independent of $\mathbf{x}$.

Finally, the classical Hamiltonian is given by

$H_{\text{cl}}(t)=\sum_{\nu}\frac{p_{\nu}^{2}}{2m_{\nu}}+U(\mathbf{x}),$ (5)

where the summation is taken over all classical degrees of freedom. The
classical momentum for the $\nu$ degree of freedom is given by $p_{\nu}$,
while $m_{\nu}$ is the corresponding mass and $U(\mathbf{x})$ is the multi-
dimensional potential.

###  II.2 Forces on the Classical Coordinate

The force operator for the force acting on the classical coordinates due to
the quantum environment is given by

$\displaystyle\mathbf{f}(t)$ $\displaystyle=-\nabla H(t),$ (6)
$\displaystyle=-\nabla H_{M}(\mathbf{x})-\nabla U(\mathbf{x}),$ (7)

where $\nabla=[\partial_{\nu_{1}},\partial_{\nu_{2}},\ldots]^{T}$ and
$\partial_{\nu_{1}}$ is the partial derivative with respect to the classical
coordinate $\nu_{1}$. We have retained only the Hamiltonian components which
depend on the classical coordinates. By taking the quantum average of the
force operator, we introduce the so-called Ehrenfest force according to

$\mathbf{F}^{\text{ehr}}(t)=-\langle\nabla H_{M}(\mathbf{x})\rangle.$ (8)

This is the main object of the present study. As discussed, the Ehrenfest
force alone is ill-equipped to capture the entirety of the dynamics of the
system. In order to completely describe the classical motion, an additional
stochastic force term is required according to

$f_{\nu}(t)=-\partial_{\nu}U+F^{\text{ehr}}_{\nu}(t)+\delta f_{\nu}(t),$ (9)

where the subscript $\nu$ denotes the force acting on the $\nu$ classical
degree of freedom. This equation can be obtained by taking the classical limit
of the Feynman-Vernon influence functional for a system interacting with a
quantum environment and applying a Hubbard-Stratonovich transformation, in
which the dynamical evolution is split into a deterministic (Ehrenfest)
component and a stochastic component[74, 45]. For the purposes of our
simulations, we disregard the stochastic component in (9) and instead focus
solely on the Ehrenfest force along with the classical force arising from $U$
since an accurate description of joule heating is not the aim of this study.
We separately consider the influence of the stochastic force in Section III.4.
We choose to model the classical potential for any degree of freedom according
to a harmonic potential as given by

$U(\mathbf{x})=\sum_{\nu}\frac{1}{2}k_{\nu}(\nu-\nu_{0})^{2},$ (10)

where $k_{\nu}$ is the spring constant for the $\nu$ degree of freedom, and
$\nu_{0}$ is a parameter which shifts the potential along the $\nu$
coordinate. The model is is not restricted to this choice of potential; it is
capable of dealing with arbitrary classical potentials. Our primary tool for
evaluating the forces is the non-equilibrium Green’s function. For our
Hamiltonian, the Ehrenfest force can be represented in terms of the lesser
component of exact non-equilibrium Green’s functions as

$F_{\nu}^{\text{ehr}}(t)=i\text{Tr}\Big{\\{}\partial_{\nu}h(t)G^{<}(t,t)\Big{\\}},$
(11)

where $h$ contains the elements of the molecular Hamiltonian and the elements
of $G^{<}$ are defined according to

$G_{ij}^{<}(t,t^{\prime})=i\langle
d_{j}^{\dagger}(t^{\prime})d_{i}(t)\rangle.$ (12)

Each are square matrices with dimension corresponding to the number of
molecular electronic energy levels.

###  II.3 Perturbative Expansion of the Ehrenfest Force

By enforcing a time-scale separation between the classical degrees of freedom
and the electronic environment, a perturbative expansion of the Ehrenfest
force can be taken, from which we can derive the Langevin description of
dynamics [32, 41, 33, 44, 42, 34, 36, 47]. Here, our small parameter is
$\Omega/\Gamma$, where $\Omega$ is the characteristic classical vibrational
frequency and $1/\Gamma$ is associated with the electron tunnelling time in
molecular electronic junctions[1]. This amounts to assuming that the classical
coordinate vibrates over much longer time-scales than the time an electron
spends in the central region. This interpretation is appropriate for the
systems considered in our paper. $\Gamma$ is generally a highly tunable
experimental parameter while $\Omega$ depends on the vibrational mode
considered. In order to identify the separate time-scales, we introduce Wigner
coordinates for the time. These are the central time ($T$) and relative time
($\tau$) which are defined according to

$T=\frac{t+t^{\prime}}{2},\;\;\;\;\;\;\;\tau=t-t^{\prime}.$ (13)

We then define an auxilliary force $\mathcal{F}(t,t^{\prime})$ which is equal
to $F^{\text{ehr}}(t)$ when $t=t^{\prime}$. Thus, we have

$\mathcal{F}(t,t^{\prime})=i\text{Tr}\Big{\\{}\partial_{\nu}h(t)G^{<}(t,t^{\prime})\Big{\\}}.$
(14)

Application of the Wigner transform along with its inverse results in

$\mathcal{F}(t,t^{\prime})=i\text{Tr}\Big{\\{}\int\frac{d\omega}{2\pi}e^{-i\omega\tau}\partial_{\nu}h(T)\widetilde{G}^{<}(T,\omega)\Big{\\}},$
(15)

where $G$ denotes a Green’s function in the time domain, and $\widetilde{G}$
denotes the Wigner space defined as

$\widetilde{G}(T,\omega)=\int
d(t-t^{\prime})e^{i\omega(t-t^{\prime})}G(t,t^{\prime}).$ (16)

We adopt the convention that if the integral limits are not shown explicitly,
they range from $-\infty$ to $\infty$. We will often subdue the $(T,\omega)$
indices for brevity. We now propose the ansatz,

$\widetilde{G}^{<}=\widetilde{G}_{(0)}^{<}+\widetilde{G}_{(1)}^{<}+\widetilde{G}_{(2)}^{<}+...,$
(17)

in which $\widetilde{G}_{(n)}^{<}$ is of $n^{th}$ order in our small
parameter. Here, $\widetilde{G}_{(0)}^{<}$ is adiabatic and corresponds to the
Born-Oppenheimer approximation, while the higher order terms go beyond this
and account for the non-zero motion of the classical coordinates. Applying
this ansatz and letting $t=t^{\prime}$ finally yields

$\displaystyle F^{\text{ehr}}_{\nu}(t)$
$\displaystyle=i\text{Tr}\Big{\\{}\int\frac{d\omega}{2\pi}\partial_{\nu}h(t)\Big{(}\widetilde{G}_{(0)}^{<}(t,\omega)+\widetilde{G}_{(1)}^{<}(t,\omega)+...\Big{)}\Big{\\}},$
(18) $\displaystyle=\underbrace{F_{\nu,(0)}+F_{\nu,(1)}}_{\text{retained in
Langevin approach}}+F_{\nu,(2)}+....$ (19)

Retaining only the zeroth (adiabatic) and first order contributions to the
Ehrenfest force yields the deterministic forces present in the Langevin
description, where the adiabatic force is purely position dependent while the
first order force is the correction which accounts for non-zero velocities.
The viscosity coefficient matrix $\hat{\xi}$ is obtained via

$\mathbf{F_{(1)}}=-\hat{\xi}\mathbf{v},$ (20)

where $\hat{\xi}$ is an $n\times n$ matrix; $n$ corresponding to the number of
classical degrees of freedom in the system.

###  II.4 Solving for the Adiabatic and First Order Green’s Functions

The adiabatic and first order Green’s functions in the frequency domain can be
solved for via the Keldysh-Kadanoff-Baym equations, expressed in the Wigner
space as [44, 41, 32]

$\Big{(}\omega+\frac{i}{2}\partial_{T}-e^{\frac{1}{2i}\partial_{\omega}^{G}\partial_{T}^{h}}h(T)\Big{)}\widetilde{G}^{R/A}=I\\\
+e^{\frac{1}{2i}\left(\partial_{T}^{\Sigma}\partial_{\omega}^{G}-\partial_{\omega}^{\Sigma}\partial_{T}^{G}\right)}\widetilde{\Sigma}^{R/A}\widetilde{G}^{R/A},$
(21)

and

$\Big{(}\omega+\frac{i}{2}\partial_{T}-e^{\frac{1}{2i}\partial_{\omega}^{G}\partial_{T}^{h}}h(T)\Big{)}\widetilde{G}^{</>}=e^{\frac{1}{2i}\left(\partial_{T}^{\Sigma}\partial_{\omega}^{G}-\partial_{\omega}^{\Sigma}\partial_{T}^{G}\right)}\\\
\times\Big{(}\widetilde{\Sigma}^{R}\widetilde{G}^{</>}+\widetilde{\Sigma}^{</>}\widetilde{G}^{A}\Big{)},$
(22)

where we have shown the retarded/advanced and the lesser/greater collectively.
Here, we adopt the convenient notation for derivatives, $\partial_{T}^{G}$,
which denotes a partial derivative acting on the $G$ term with respect to $T$
and so on. Again, each quantity here is a square matrix with dimensions
corresponding to the number of electronic levels in the central region. Since
$T$ is associated with the motion of the classical coordinate, terms
containing $(\partial_{T})^{n}$ will be of $n^{\text{th}}$ in our small
parameter. Thus, we apply our ansatz in (17), along with applying an expansion
of the exponentials in orders of $\partial_{T}$. Truncating after the zeroth
order and solving for $\widetilde{G}_{(0)}$ yields the standard adiabatic
Green’s functions as follows:

$\widetilde{G}_{(0)}^{R/A}=\Big{(}\omega\hat{I}-h-\widetilde{\Sigma}^{R/A}\Big{)}^{-1},$
(23)
$\widetilde{G}_{(0)}^{</>}=\widetilde{G}_{(0)}^{R}\widetilde{\Sigma}^{</>}\widetilde{G}_{(0)}^{A}.$
(24)

An equivalent approach for the first order yields[44, 41, 32]

$\widetilde{G}_{(1)}^{R/A}=\frac{1}{2i}\widetilde{G}_{(0)}^{R/A}\Big{[}\widetilde{G}_{(0)}^{R/A},\partial_{T}h\Big{]}\widetilde{G}_{(0)}^{R/A},$
(25)
$\widetilde{G}_{(1)}^{</>}=\widetilde{G}_{(0)}^{R}\widetilde{\Sigma}^{</>}\widetilde{G}_{(1)}^{A}+\widetilde{G}_{(1)}^{R}\widetilde{\Sigma}^{</>}\widetilde{G}_{(0)}^{A}\\\
+\frac{1}{2i}\widetilde{G}_{(0)}^{R}\Big{(}\partial_{T}h\widetilde{G}_{(0)}^{R}\partial_{\omega}\widetilde{\Sigma}^{</>}+\widetilde{G}_{(0)}^{</>}\partial_{T}h+h.c\Big{)}\widetilde{G}^{A}_{(0)}.$
(26)

The self-energy components in the frequency domain are given by

$\widetilde{\Sigma}^{R}_{\alpha}=-\frac{i}{2}\Gamma_{\alpha},\;\;\;\;\;\;\widetilde{\Sigma}^{A}_{\alpha}=\frac{i}{2}\Gamma_{\alpha},$
(27)
$\widetilde{\Sigma}^{<}_{\alpha}=if_{\alpha}(\omega)\Gamma_{\alpha},\;\;\;\;\;\;\widetilde{\Sigma}^{>}_{\alpha}=-i[1-f_{\alpha}(\omega)]\Gamma_{\alpha}.$
(28)

Here, $f_{\alpha}(\omega)$ is the Fermi-Dirac distribution which describes the
$\alpha$ lead. In solving for the above, we have applied the wide-band
approximation to the leads such that
$\Gamma_{\alpha}(\omega)=\Gamma_{\alpha}$. The explicit expression for an
element of $\Gamma_{\alpha}$ is

$\Gamma_{c,c^{\prime},\alpha}=2\pi\rho_{\alpha}t^{*}_{c,\alpha}t_{c^{\prime},\alpha},$
(29)

where $\rho_{\alpha}$ is the now constant density of states for the $\alpha$
lead, and $t_{c,\alpha}$ is the coupling strength between any state in the
$\alpha$ lead and the central region electronic state $c$. We take
$\Gamma_{\alpha}$ as a parameter for our model. In equations (27) to (29), the
full quantity can be obtained by taking a sum over the leads $\alpha$. For
example,

$\Gamma_{c,c^{\prime}}=\sum_{\alpha}\Gamma_{c,c^{\prime},\alpha}.$ (30)

###  II.5 Time-Stepping Approach to the Evolution of the Green’s Function

Clearly, if we are able to calculate $G^{<}(t,t)$ at each time, then (11)
allows us to readily calculate the Ehrenfest force without need to resort to a
perturbative approach involving time-scale separation. Our method is based off
of the seminal work of Jauho et al. in reference [75], from which we borrow
the main defined quantities. A similar time-stepping solution to Ehrenfest
dynamics was presented in reference [68]. However, our application of the
approach as detailed in this section is more sophisticated and allows for
better convergence and stability. In approaching these calculations, we borrow
the main defined quantities from reference [75]. The lesser and greater
Green’s function in the time domain are given by the Keldysh equation,

$G^{</>}(t,t^{\prime})=\int dt_{1}\int
dt_{2}G^{R}(t,t_{1})\Sigma^{</>}(t_{1},t_{2})G^{A}(t_{2},t^{\prime}).$ (31)

For the sake of clarity, in this section we will use $(...)$ to denote a
functional dependence, whereas $[...]$ will denote a term in the equation. Now
we will take (31) and express the $\Sigma^{</>}$ term according to the inverse
Wigner transform of (28). In the lesser case, we find

$\displaystyle G^{<}(t,t^{\prime})$ $\displaystyle=i\sum_{\alpha}\int
dt_{1}\int
dt_{2}G^{R}(t,t_{1})\int\frac{d\omega}{2\pi}e^{-i\omega[t_{1}-t_{2}]}f_{\alpha}(\omega)\Gamma_{\alpha}G^{A}(t_{2},t^{\prime}),$
(32) $\displaystyle=i\sum_{\alpha}\int dt_{1}\int
dt_{2}G^{R}(t,t_{1})\int\frac{d\omega}{2\pi}e^{-i\omega[t_{1}-t_{2}]}f_{\alpha}(\omega)\Gamma_{\alpha}G^{A}(t_{2},t)e^{i\omega[t-t]}e^{i\omega[t^{\prime}-t^{\prime}]},$
(33)
$\displaystyle=i\sum_{\alpha}\int\frac{d\omega}{2\pi}f_{\alpha}(\omega)e^{i\omega[t^{\prime}-t]}\int
dt_{1}e^{i\omega[t-t_{1}]}G^{R}(t,t_{1})\Gamma_{\alpha}\int
dt_{2}e^{-i\omega[t^{\prime}-t_{2}]}G^{A}(t_{2},t^{\prime}),$ (34)
$\displaystyle=i\sum_{\alpha}\int\frac{d\omega}{2\pi}e^{-i\omega\tau}f_{\alpha}(\omega)A(\omega,t)\Gamma_{\alpha}A^{\dagger}(\omega,t^{\prime}).$
(35)

An equivalent derivation can be applied in the greater case to find

$G^{>}(t,t^{\prime})=-i\sum_{\alpha}\int\frac{d\omega}{2\pi}e^{-i\omega\tau}[1-f_{\alpha}(\omega)]A(\omega,t)\Gamma_{\alpha}A^{\dagger}(\omega,t^{\prime}).$
(36)

Here we have defined the quantity,

$A(\omega,t)=\int dt_{1}e^{i\omega[t-t_{1}]}G^{R}(t,t_{1}).$ (37)

The evolution of $G^{R}$ is found via the Kadanoff-Baym equation in the time
domain:

$\Big{(}i\frac{\partial}{\partial
t}-h(t)\Big{)}G^{R}(t,t^{\prime})=\delta(t-t^{\prime})+\int
dt_{1}\Sigma^{R}(t,t_{1})G^{R}(t_{1},t^{\prime}).$ (38)

By utilising the fact that

$\Sigma^{R}(t,t^{\prime})=-\frac{i}{2}\Gamma\delta(t-t^{\prime}),$ (39)

we can solve for the general solution of (38) as

$G^{R}(t,t^{\prime})=-i\Theta(t-t^{\prime})\hat{T}\text{exp}\Big{\\{}-i\int_{t^{\prime}}^{t}dt_{1}h(t_{1})-\frac{1}{2}\Gamma[t-t^{\prime}]\Big{\\}}.$
(40)

Here, $\hat{T}$ is the time-ordering operator and $\Theta$ is the Heaviside
step function. By substituting (40) into (37), we arrive at

$A(\omega,t)=\\\
-i\int^{t}_{-\infty}dt_{1}\hat{T}\text{exp}\Big{\\{}i\omega[t-t_{1}]-i\int_{t_{1}}^{t}dt_{2}h(t_{2})-\frac{1}{2}\Gamma[t-t_{1}]\Big{\\}}.$
(41)

This is the quantity which we will iteratively evolve forwards in time, from
which we can then extract $G^{<}$ via (35). We begin by applying a step
$\Delta t$ forwards in time to (41)),

$\displaystyle A(\omega,t+\Delta t)$ $\displaystyle=-i\int^{t+\Delta
t}_{-\infty}dt_{1}\hat{T}\text{exp}\Big{\\{}i\omega[t+\Delta
t-t_{1}]-i\int_{t_{1}}^{t+\Delta t}dt_{2}h(t_{2})-\frac{1}{2}\Gamma[t+\Delta
t-t_{1}]\Big{\\}},$ (42)
$\displaystyle=-i\int^{t}_{-\infty}dt_{1}\hat{T}\text{exp}\Big{\\{}i\omega[t+\Delta
t-t_{1}]-i\int_{t_{1}}^{t+\Delta t}dt_{2}h(t_{2})-\frac{1}{2}\Gamma[t+\Delta
t-t_{1}]\Big{\\}}$ (43) $\displaystyle-i\int^{t+\Delta
t}_{t}dt_{1}\hat{T}\text{exp}\Big{\\{}i\omega[t+\Delta
t-t_{1}]-i\int_{t_{1}}^{t+\Delta t}dt_{2}h(t_{2})-\frac{1}{2}\Gamma[t+\Delta
t-t_{1}]\Big{\\}},$ (44) $\displaystyle=A^{A}(\omega,t+\Delta
t)+A^{B}(\omega,t+\Delta t).$ (45)

Here we have split $A$ into $A^{A}$ which reflects how the history before $t$
informs the system at $t+\Delta t$, and $A^{B}$ which contains the effects of
the system from $t$ to $t+\Delta t$. $A^{A}$ can be further simplified by
partitioning the inner integral as follows:

$\displaystyle A^{A}(\omega,t+\Delta t)$
$\displaystyle=-i\int^{t}_{-\infty}dt_{1}\hat{T}\text{exp}\Big{\\{}i\omega\Delta
t-i\int_{t}^{t+\Delta t}dt_{2}h(t_{2})-\frac{1}{2}\Gamma\Delta t\Big{\\}}$
(46)
$\displaystyle\times\hat{T}\text{exp}\Big{\\{}i\omega[t-t_{1}]-i\int_{t_{1}}^{t}dt_{3}h(t_{3})-\frac{1}{2}\Gamma[t-t_{1}]\Big{\\}},$
(47) $\displaystyle=\hat{T}\text{exp}\Big{\\{}i\omega\Delta
t-i\int_{t}^{t+\Delta t}dt_{2}h(t_{2})-\frac{1}{2}\Gamma\Delta
t\Big{\\}}A(\omega,t).$ (48)

We will make the assumption that for a sufficiently small $\Delta t$, $h(t)$
can be approximated as being piecewise constant. Thus, at some time-step
$t_{n}$, we assume that for $t_{n}<t<t_{n+1}$, we have $h(t)=h(t_{n})$. This
entails that our time step size should be sufficiently small in the
computational implementation when calculating classical trajectories. This
enables us to directly calculate the integrals over $h(t)$ appearing in (48)
and (44) along with removing the time-ordering operators. We thus re-express
$A^{A}(\omega,t+\Delta t)$ and $A^{B}(\omega,t+\Delta t)$ according to

$A^{A}(\omega,t+\Delta t)=\text{exp}\Big{\\{}i\Delta
t\Big{(}\omega-h(t)+\frac{i}{2}\Gamma\Big{)}\Big{\\}}A(\omega,t),$ (49)

and

$\displaystyle A^{B}(\omega,t+\Delta t)$ $\displaystyle=-i\int^{\Delta
t}_{0}dt_{2}\text{exp}\Big{\\{}it_{2}\Big{(}\omega-h(t)+\frac{i}{2}\Gamma\Big{)}\Big{\\}},$
(50) $\displaystyle=-i\Big{(}e^{\Lambda\Delta t}-I\Big{)}\Lambda^{-1},$ (51)

where $I$ is the identity matrix and for cleanliness, we have used
$\Lambda=i\Big{(}\omega-h(t)+\frac{i}{2}\Gamma\Big{)}$. Here we have used the
matrix exponential identity:

$\int^{T}_{0}dte^{Xt}=\Big{(}e^{XT}-I\Big{)}X^{-1}.$ (52)

Now we have everything we need to be able to perform an Ehrenfest simulation
of the classical dynamics. The algorithm functions as follows. First, we
assume an adiabatic initial condition, such that (37) can be calculated
through the use of (23). The initial force can then be computed according to
(35) and (11) given the chosen initial conditions for the classical
coordinate. The calculated Ehrenfest force feeds into an external algorithm
which updates the nuclear dynamics to the next time step, at which point we
use (49) and (51) to non-adiabatically compute $A$ at the next time step. From
here the process repeats, yielding a wholly non-adiabatic, classical
trajectory in time.

##  III RESULTS

###  III.1 Purely Positive Viscosity in the Single-Level Case

In the case where there is only a single electronic level in the central
region and only the molecular Hamiltonian contains classical dependence, the
diagonal components of the viscosity matrix $\xi$ can be given by

$\xi_{\nu\nu}=\frac{\Gamma(\partial_{\nu}h)^{2}}{4\pi}\int
d\omega\frac{\frac{f_{L}}{T_{L}}(1-f_{L})\Gamma_{L}+\frac{f_{R}}{T_{R}}(1-f_{R})\Gamma_{R}}{\Big{(}(\omega-h)^{2}+\Gamma^{2}/4\Big{)}^{2}},$
(53)

where $T_{L}$ and $T_{R}$ are the temperatures for the left and right lead,
respectively. In our model, we set $T_{L}=T_{R}=0.001$. A rigorous derivation
of (53) has been relegated to Appendix A. Thus, we observe that the diagonal
viscosity elements are purely positive as each term in (53) is positive. From
a Langevin standpoint, theory is known to produce negative viscosities in non-
equilibrium in the case of multiple electronic levels[44, 54, 47], as well as
resulting from non-Condon effects in which the geometry at the leads-interface
becomes relevant[32]. Due to the cumbersome nature of the equations in the
latter case, we focus only on the former in this paper.

###  III.2 Single Classical Degree of Freedom

Here we present results for a two-level system with a single classical degree
of freedom. The matrix elements of the molecular Hamiltonian take the form

$H_{M}=\begin{pmatrix}-\lambda x&\eta\\\ \eta&\lambda x\end{pmatrix},$ (54)

where the classical coordinate $x$ has the effect of linearly shifting the two
atomic orbitals and $\eta$ dictates the overlap. We choose the left lead to be
coupled to the first state and the right lead to be coupled to the second
state which gives our $\Gamma$ matrices the following form:

$\Gamma_{L}=\begin{pmatrix}\gamma&0\\\
0&0\end{pmatrix},\;\;\;\;\;\;\Gamma_{R}=\begin{pmatrix}0&0\\\
0&\gamma\end{pmatrix}.$ (55)

In this section, the common parameters are $\lambda=\eta=0.1$, $\gamma=0.05$,
$k=1$, $x_{0}=-0.22$ and $\mu_{L}=-\mu_{R}=0.12$, where $\mu_{\alpha}$ is the
chemical potential of the $\alpha$ lead. In contrast to the single-level case,
this model allows us to observe negative viscosities. An example of this is
shown in Figure 1 where the coordinate dependent viscosity coefficient (now a
scalar function in the case of a single classical degree of freedom) becomes
negative in a small region around $x\approx-0.2$.

Figure 1: The viscosity coefficient as a function of the classical coordinate.
Atomic units are used for quantities here and in all further plots in the
paper.

By utilising our time-stepping algorithm to evolve $G^{<}$ in time as
presented in II.5 we can calculate the Ehrenfest force as a function of time,
from which we can then computationally simulate the dynamics of the classical
coordinate. In Figure 2, we calculate a classical Ehrenfest trajectory and
record $F^{\text{ehr}}$ at each time step. We also record the corresponding
forces which the classical coordinate would experience under the perturbative
truncation of $F^{\text{ehr}}$ after the first order. We note that the methods
clearly differ in the transient regime depending on the chosen initial
conditions, before each settles into periodicity as the coordinate oscillates
in time. The cause of the differences between the calculated Ehrenfest and
perturbative forces are made clearer in Figure 3, where a direct comparison
has been made of the magnitude of the first order force relative to the sum of
all higher order forces ($F_{(2)}+F_{(3)}+...$) for different values of
$\Omega/\Gamma$. This is achieved by varying the effective mass of the
classical coordinate while keeping all other parameters constant. We can then
estimate $\Omega$ by assuming that the oscillations are approximately
harmonic. We observe that when $\Omega/\Gamma$ is small in (a), the higher
order forces have only a small relative contribution. In (b) however,
$\Omega/\Gamma$ is no longer small and thus our perturbative assumption is no
longer satisfied, meaning the higher order forces have become increasingly
relevant.

Figure 2: Trajectory in time of $F^{\text{ehr}}(t)$ (Ehrenfest) and
$F_{(0)}(t)+F_{(1)}(t)$ (perturbative). $\Omega/\Gamma\approx 0.32$.

Figure 3: Comparison of the first order force against the sum of all higher
order forces. (a) $\Omega/\Gamma\approx 0.14$, (b) $\Omega/\Gamma\approx 3.2$.

Calculation of the time-dependent forces allows us to then simulate the phase-
space trajectory of the classical coordinate as a function of time. Figure 4
shows time dependent trajectories of the classical coordinate $x$ for three
different values of the small parameter $\Omega/\Gamma$. We also calculate the
instantaneous power of the classical coordinate at each time step according to
$(F^{\text{ehr}}-F_{(0)})\cdot\frac{p}{m}$, which includes only the effects of
the excitational and dissipative forces. The classical coordinate is
intentionally confined to the region of negative viscosity shown in Figure 1
via our choice of $U(x)$ such that under the perturbative assumption, the
instantaneous power will be positive at all time-steps. We observe that when
$\Omega/\Gamma\ll 1$ as in (a); rendering the perturbative approximation as
valid, the instantaneous power is overwhelmingly positive in agreeance with
the perturbative approximation. This entails that the oscillations will
increase until reaching regions of positive viscosity, whereby they will form
a limit cycle. This is further illustrated in Figure 6 via long trajectories
for the same parameters, where the amplitude of oscillations increases in both
the Ehrenfest and perturbative approaches. This result demonstrates that even
in the absence of time-scale separation, negative viscosities emerge and
dictate the behaviour of the system. However, we observe the resultant
trajectories from the two methods to diverge from each other over longer time-
scales, demonstrating the importance of the higher order forces even for these
parameters where we expect they can be reasonably neglected. We note that
Ehrenfest dynamics has been used to observe the effects of negative
viscosities in different regimes in the literature[68, 55, 56, 62]. In Figure
4 (b) and (c) where $\Omega/\Gamma>1$ and our perturbative assumption is no
longer adequate, we observe that instantaneous power is more often negative
wherein the electronic environment is taking energy away from the classical
coordinate. This means that the higher order terms (2nd order and above) in
our perturbative expansion for $F^{\text{ehr}}$ have become more relevant and
the dissipative nature of these forces are overwhelming the negative viscosity
produced by $F_{(1)}$. These results are summarised in Figure 5 in which the
average power input to the classical coordinate over a period of oscillation
is calculated and classified as positive (blue) or negative (red). We observe
that the average power input to the classical coordinate is positive far
beyond when the perturbative assumption is valid. We anticipate that this cut-
off between positive and negative average power input is highly dependent on
the model and parameters at hand. However, this demonstrates that while
negative viscosities will still emerge under a numerically exact approach, it
can be dominated by dissipative higher order forces for large values of
$\Omega/\Gamma$. These results suggest that the inclusion of the higher order
forces which emerge through Ehrenfest dynamics acts to further subdue
classical vibrations within the system.

Figure 4: Samples of Ehrenfest trajectories of the classical coordinate
against time for (a) $\Omega/\Gamma\approx 0.14$, (b) $\Omega/\Gamma\approx
2.2$, (c) $\Omega/\Gamma\approx 3.2$. Colour shows the instantaneous power of
the classical coordinate. Figure 5: Regimes of negative (blue) or positive
(red) average power input to the classical coordinate over a period of
oscillation. The green points refer to the corresponding plots in Figure 4.
Figure 6:  Comparison of long trajectories of the classical coordinate in time
using the Ehrenfest approach and the truncated perturbative approach. The
viscosity coefficient is negative between the red dashed lines and positive
outside. $\Omega/\Gamma\approx 0.14$.

###  III.3 Two Classical Degrees of Freedom

The algorithm for calculating $F^{\text{ehr}}$ can be readily extended to
account for many classical degrees of freedom. In this section, we consider
another 2-level model consisting now of two classical degrees of freedom; a
stretching component and an angular component. As a result, in contrast to the
case of a single vibrational degree of freedom, $F_{(0)}$ now becomes
typically non-conservative. The stretching coordinate $q$ modulates the
hopping amplitude between electronic sites in the central region, while the
angular coordinate $\theta$ accounts for the shift in electric levels due to
the electric field produced by the applied voltage bias. Thus, the molecular
Hamiltonian now takes the form

$H_{M}=\begin{pmatrix}h_{L}(\theta)&\eta(q)\\\
\eta(q)&h_{R}(\theta)\end{pmatrix}.$ (56)

The hopping amplitude is given by

$\eta(q)=\eta e^{-q}\Big{(}1+q+q^{2}/3\Big{)}.$ (57)

This is a generic $1s$-orbital overlap scaled by some constant $\eta$[76].
Meanwhile, the atomic levels are given by

$h_{L/R}(\theta)=\pm Eq_{0}\cos(\theta),$ (58)

where $q_{0}$ is the constant bond-length parameter which appears in our
classical potential for $q$, while $E$ is the linear approximation to the
electric field across the junction as determined by

$E=\frac{\mu_{L}-\mu_{R}}{L},$ (59)

and $L$ is the junction length. For all results in this section, we set
$\eta=0.1$, $q_{0}=2$, $\theta_{0}=0$, $L=4$, $k_{q}=0.1$, $k_{\theta}=0.05$,
and $\mu_{L}=-\mu_{R}=0.1$. All other components of the full Hamiltonian take
the same form as in section III.2.

In Figure 7, we calculate an Ehrenfest trajectory and record the Ehrenfest and
perturbative forces on each degree of freedom separately. Given that we have
more than one classical degree of freedom, this also includes the so-called
non-equilibrium anti-symmetric forces in which the motion of $\theta$ induces
a force on $q$ and vice versa; the net anti-symmetric force being
perpendicular to the motion of the classical coordinate. Figure 8 performs a
direct comparison between the time dependent trajectories calculated via the
Ehrenfest force and the perturbative approximation. To do so, we first
calculate a classical Ehrenfest trajectory of our two coordinates. From this
trajectory at some point in time after the unusual transient behaviour in
$F^{\text{ehr}}$ has subsided, we extract the initial conditions for a
comparative trajectory where the forces are calculated according to the
perturbative approximation. In Figure 8 (a), the trajectories are compared for
large effective masses for each coordinate such that the perturbative
assumption should be satisfied. We observe only small differences between the
trajectories resulting from the different methods in this case. Contrarily,
Figure 8 (b) demonstrates for small effective masses that the trajectories
deviate away from each other almost instantly and undergo largely different
oscillatory behaviour. This is because the perturbative assumption is no
longer valid since $\Omega/\Gamma$ is not sufficiently small.

Figure 7: The force trajectories in time for the Ehrenfest force (blue) and
perturbative force (black) acting on (a) $q$ and (b) $\theta$.
$\Omega_{q}/\Gamma\approx 0.32$, $\Omega_{\theta}/\Gamma\approx 0.22$.

Figure 8: Short classical trajectories in the coordinate-space comparing the
Ehrenfest method and perturbative method. Green(red) points denote the
start(end) of the trajectory. The vector field shows $F^{cl}+F_{(0)}$. (a)
$\Omega_{q}/\Gamma\approx 0.32$, $\Omega_{\theta}/\Gamma\approx 0.22$, and a
trajectory length of $600$, and (b) $\Omega_{q}/\Gamma=1$,
$\Omega_{\theta}/\Gamma\approx 0.7$, and a trajectory length of $200$.

###  III.4 Evaluating the Diffusion Coefficient

The stochastic force term which we have thus far avoided in (9) is a
stochastic Gaussian process entirely defined by the following[45]:

$\langle\delta f(t)\rangle=0,\;\;\;\;\;\;\ \langle\delta f(t)\delta
f(t^{\prime})\rangle=D(t,t^{\prime}),$ (60)

where we will refer to $D(t,t^{\prime})$ as the exact diffusion coefficient.
We will consider a system consisting of a single electronic level coupled to a
single classical degree of freedom. Our molecular Hamiltonian is given by

$H_{M}=h_{0}+\lambda x,$ (61)

where $x$ represents any generic classical degree of freedom. We will let
$h_{0}=0$ and our $\Gamma_{\alpha}$ become scalar inputs in the single level
case. The exact diffusion coefficient in (60) for our model can then be
expressed in terms of nonequilibrium Green’s functions according to

$D(t,t^{\prime})=\lambda^{2}G^{>}(t,t^{\prime})G^{<}(t^{\prime},t).$ (62)

A derivation of this result is presented in Appendix B. Our $D(t,t^{\prime})$
here is not exclusively real. Since our stochastic force is classical, we
symmetrize (62) by taking the real component such that our expression for
$D(t,t^{\prime})$ is now in correspondence with reference [45] as per

$D(t,t^{\prime})=\text{Re}\left\\{\lambda^{2}G^{>}(t,t^{\prime})G^{<}(t^{\prime},t)\right\\}.$
(63)

$D(t,t^{\prime})$ accounts for the effects of the random fluctuations about
the mean-field of the electronic environment on the classical degrees of
freedom along with accounting for the feedback of the classical coordinate on
the electronic environment due to its motion. However, generally a time-scale
separation within the system is utilised in order to produce a perturbative
solution to $D(t,t^{\prime})$ in which the non-adiabatic feedback due to the
motion of the classical coordinate is not included; in other words, the
classical coordinate evolves adiabatically. We can solve for the adiabatic
form of $D(t,t^{\prime})$ by considering the Wigner transform of (62) and
retaining only the adiabatic terms. With some work which has been relegated to
Appendix C, the Wigner transform of (62) can be expressed as

$\tilde{D}(\omega,T)=\lambda^{2}\int\frac{d\omega^{\prime}}{2\pi}\tilde{G}^{>}(\omega+\omega^{\prime},T)\tilde{G}^{<}(\omega^{\prime},T).$
(64)

We note that (64) is a result of the phonon self-energy derived in reference
[47]. Now, the adiabatic diffusion is found by taking $\tilde{G}^{</>}$ to be
our adiabatic Green’s functions as in (24). The subsequent application of the
inverse Wigner transform then yields the adiabatic diffusion in the time
domain as

$D_{(0)}(\tau,T)=\lambda^{2}\int\frac{d\omega}{2\pi}e^{-i\omega\tau}\int\frac{d\omega^{\prime}}{2\pi}\tilde{G}_{(0)}^{>}(\omega+\omega^{\prime},T)\tilde{G}_{(0)}^{<}(\omega^{\prime},T).$
(65)

The adiabatic diffusion will serve as a base of comparison with the exact
diffusion in order to analyse the effects of a time-scale separation within
the system. We additionally analyse the validity of the white-noise
approximation for different parameters. The white-noise approximation is a
method of coarse-graining the diffusion coefficient which allows for a simpler
mathematical and computational treatment. This is done by assuming that the
coloured-noise exact diffusion coefficient can instead be replaced by a white-
noise equivalent according to

$D(t,t^{\prime})\equiv D^{w}(T)\delta(t-t^{\prime}),$ (66)

where we have introduced the white-noise diffusion coefficient $D^{w}$. Thus,
the aim is to accurately reproduce the effects of the stochastic force under
the assumption that it is entirely uncorrelated in time. In order to replicate
the correct dynamics using this approximation, it can be shown that $D^{w}$
must then take the form

$D^{w}(t_{f})=\int^{t_{f}}_{-t_{f}}d\tau D(t_{f},\tau),$ (67)

where to avoid ambiguity, $t_{f}$ denotes a specific point along the
trajectory. This is demonstrated in Appendix D. (67) is just the Wigner
transform of $D(t,t^{\prime})$ where $\omega=0$ such that $D^{w}$ is
independent of $\omega$; hence the name “white-noise” diffusion. Thus, the
white-noise diffusion coefficient contains information about the correlations
in the stochastic force but applies that information in a Markovian manner. A
time-scale separation can be applied to the white-noise diffusion in a similar
way as previously, whereby the adiabatic white-noise diffusion is given by

$D^{w}_{(0)}(t_{f})=\int^{t_{f}}_{-t_{f}}d\tau D_{(0)}(t_{f},\tau).$ (68)

The adiabatic white-noise diffusion is the most commonly used approximation to
the diffusion coefficient[32, 41, 44] and we note in passing that it satisfies
a fluctuation-dissipation relation with (53) when $\partial_{\nu}h=\lambda$.
The white-noise approximation is valid when correlations in the exact
diffusion coefficient decay over time-scales in which the effects of other
forces present in the system (in our case, the Ehrenfest and the classical
forces) are negligible. In order to quantify its validity, we introduce the
correlation time defined according to

$t_{\text{corr}}(t_{f})=\frac{\int^{t_{f}}_{-t_{f}}d\tau
D(t_{f},\tau)}{D(t_{f},\tau=0)}.$ (69)

This is a measure of the persistence of correlations in the stochastic force,
independent of their strength. The numerator of (69) is our expression for the
white-noise diffusion while the denominator is the variance in the stochastic
force at time $t_{f}$. With all quantities now defined, we can begin to
discuss the results.

(a)

(b)

(c)

(d)

Figure 9: Exact diffusion coefficient as a function of $\tau$ for where time
$t$ occurs at: (a)/(b) $x=0$, (c)/(d) a turning point. Left: $V=0$, Right:
$V=0.2$. Dashed line is the corresponding $D_{(0)}(t,t^{\prime})$ for the same
parameters.

We once again apply the time-stepping algorithm presented in Section II.5 to
find classical Ehrenfest trajectories in time for our classical coordinate. In
doing so, we store $A(\omega,t)$ (41) at each time-step. The two-time lesser
and greater Green’s functions can then be computed according to (35) and (36)
by inputting $A(\omega,t)$ at different points in the stored trajectory and
numerically integrating over $\omega$ for each possible value of $\tau$. We
then calculate the exact diffusion coefficient as a function of $\tau$
according to (62), where $D(\tau=0)$ corresponds to the variance in the
stochastic force at the end point of the trajectory and $D(\tau>0)$ is the
correlation in the stochastic force between the times $t$ and $t-\tau$. To
clarify the method, the stochastic force is not included in the simulation. We
instead calculate the trajectory using Ehrenfest dynamics as a means to assess
the behaviour of the exact diffusion coefficient. The common parameters in
this section are $\lambda=0.1$, $k=1$ and $x_{0}=0$.

In Figure 9 (a)-(d) we observe the exact diffusion as in (62) for different
values of our small parameter (solid line) plotted against the corresponding
adiabatic diffusion as per (65) (dashed line). Here, time $t$ corresponds to
the end point of the trajectory where in (a)/(b) the trajectory is ended at a
point when $x=0$, while in (c)/(d) the trajectory is ended at a turning point.
The left plots are calculated in equilibrium while the right plots are
calculated at a non-zero voltage. The differences between the solid and dashed
lines give an indication of the effects of the feedback on the electronic
environment due to the motion of the classical coordinate. To our knowledge,
the effects of non-adiabatic motion on the diffusion coefficient have not been
directly observed previously. As expected, this feedback becomes especially
important when the small parameter $\Omega/\Gamma$ becomes larger, such that
the perturbative truncation of (62) is no longer satisfied. However, we
observe that an increase to the voltage nullifies the effects of the feedback.
This is justified by the knowledge that the electronic tunnelling time-scale
is faster at higher voltages[1], meaning that the electronic environment can
more readily equilibrate to any changes in the classical geometry. Thus, the
time-scale separation becomes increasingly justifiable further from
equilibrium, even despite $\Omega/\Gamma$ being large. We also find that the
correlations in the stochastic force are dissipated over shorter time-scales
when $\Omega/\Gamma$ is smaller such that $D(\tau)$ approaches a shape more
reminiscent of a delta function. Finally, we note that the time-scale
separation appears less satisfactory at the turning points of the trajectory.
The acceleration of the classical coordinate is largest around the turning
points and we posit that acceleration dependent terms become important here,
which are otherwise unaccounted for under the assumption of adiabatic motion.

The effectiveness of the time-scale separation for different ending positions
along the trajectory is quantified in Figure 10, in which we have calculated
$D^{w}/D^{w}_{(0)}$ for different end positions along a period of the
trajectory when $\Omega/\Gamma=1$. For these parameters, the time-scale
separation is ineffective as $D^{w}$ is over twice as large as $D^{w}_{(0)}$
at a minimum. We observe the adiabatic assumption to be at its weakest in the
vicinity of the turning points.

Figure 10: An example trajectory of the classical coordinate. Colour scale
shows the calculated $D^{w}/D^{w}_{(0)}$ at that point in the trajectory.
Parameters: $\Omega/\Gamma=1$, $V=0$.

In Figure 11, we use $t_{\text{corr}}/T_{p}$ as a measure of the suitability
of the white-noise approximation, where $T_{p}$ is the period of oscillations
in the classical coordinate. Here, $t_{\text{corr}}$ contains any information
about the electronic forces acting on the coordinate, while $T_{p}$ contains
information on the classical force. We observe that an increase to
$\Omega/\Gamma$ results in a corresponding increase to
$t_{\text{corr}}/T_{p}$, implying that the white-noise approximation is more
valid at smaller $\Omega/\Gamma$ where the time-scale separation is more well-
defined. However, we note that the validity decreases upon increasing the
voltage which would ordinarily serve to further increase the time-scale
separation. For this system, a larger voltage results in more persistent
correlations in the stochastic force which emerges in the numerator of (69),
while the change in the denominator is comparatively irrelevant.

Figure 11: The ratio of the correlation time to the period of classical
oscillation, as a function of $\Omega/\Gamma$ for equilibrium and non-
equilibrium cases. Trajectory length = 300.

This counter-intuitive result also indirectly emerges in Figure 12 where we
observe the white-noise approximation to become more applicable outside of the
voltage window and away from the chemical potentials of the leads. In this
region, electrons tunnel more slowly through the central region. We do not yet
have a convincing explanation for this.

Figure 12: The ratio of correlation time to the period of classical
oscillation for shifted values of $h_{0}$. The dashed blue lines show the
edges of the voltage window. Parameters: $\Omega/\Gamma=1$.

##  IV CONCLUSIONS

In this paper, we have introduced a novel time-stepping algorithm for
evaluating the exact lesser Green’s function at equal times, which allows us
to simulate the trajectory of multiple classical coordinates simultaneously
via an Ehrenfest approach. We use this to benchmark the commonly used Langevin
approach which necessitates the use of a time-scale separation between
classical and electronic coordinates within the system. We observe that
despite our avoidance of a time-scale separation within our system, we observe
negative dissipations (positive power input to the classical coordinate) as
predicted by the Langevin approach. We also note that these negative
dissipations can be overwhelmed by positive dissipations due to higher order
forces unaccounted for in the Langevin approach. We also apply our method to
the calculation of the diffusion coefficient in which we observe the effect of
the non-adiabatic feedback of the classical coordinate onto the electronic
environment and its behaviour under a time-scale separation. Additionally, we
assess the validity of the white-noise approximation for the diffusion
coefficient for a range of parameters and find that it is most applicable
under a clear time-scale separation within the system and is better applied
when the transport channel energy is shifted outside of the voltage window.

##  Appendix A Viscosity Coefficient in the Single-Level Case

For the case of a single electronic level and any number of classical degrees
of freedom, the first order force for an arbitrary degree of freedom $\nu$ in
(19) is given by

$F_{\nu,(1)}=\int\frac{d\omega}{2\pi}i\partial_{\nu}h\widetilde{G}_{(1)}^{<},$
(70)

where $\widetilde{G}_{(1)}^{<}$ in the single-level case can be simplified to

$\widetilde{G}_{(1)}^{<}=\frac{i}{2}\partial_{T}h\widetilde{G}_{(0)}^{R}\widetilde{G}_{(0)}^{A}\Big{(}\widetilde{G}_{(0)}^{A}-\widetilde{G}_{(0)}^{R}\Big{)}\partial_{\omega}\widetilde{\Sigma}^{<}.$
(71)

Through the use of (23), (28), and the fact that
$\partial_{\omega}f_{\alpha}=-\frac{1}{T_{\alpha}}f_{\alpha}(1-f_{\alpha})$,
some simplification yields

$\widetilde{G}_{(1)}^{<}=\frac{i\partial_{T}h\Gamma[\frac{f_{L}}{T_{L}}(1-f_{L})\Gamma_{L}+\frac{f_{R}}{T_{R}}(1-f_{R})\Gamma_{R}]}{2\Big{[}(\omega-h)^{2}+\frac{\Gamma^{2}}{4}\Big{]}^{2}}.$
(72)

We obtain the simplified form of $F_{\nu,(1)}$ by substituting (72) into (70):

$F_{\nu,(1)}=-\partial_{\nu}h\partial_{T}h\Gamma\int\frac{d\omega}{4\pi}\frac{\frac{f_{L}}{T_{L}}(1-f_{L})\Gamma_{L}+\frac{f_{R}}{T_{R}}(1-f_{R})\Gamma_{R}}{\Big{[}(\omega-h)^{2}+\frac{\Gamma^{2}}{4}\Big{]}^{2}}.$
(73)

We apply the chain rule such that
$\partial_{T}h=\sum_{\nu^{\prime}}v_{\nu^{\prime}}\partial_{\nu^{\prime}}h$
where $v_{\nu^{\prime}}$ is the velocity of the $\nu^{\prime}$ coordinate and
the summation is over all classical degrees of freedom. The diagonal
components of the force then correspond to when $\nu^{\prime}=\nu$, for which
we will utilise an additional subscript:

$F_{\nu\nu,(1)}=-v_{\nu}(\partial_{\nu}h)^{2}\Gamma\int\frac{d\omega}{4\pi}\frac{\frac{f_{L}}{T_{L}}(1-f_{L})\Gamma_{L}+\frac{f_{R}}{T_{R}}(1-f_{R})\Gamma_{R}}{\Big{[}(\omega-h)^{2}+\frac{\Gamma^{2}}{4}\Big{]}^{2}}.$
(74)

Finally, the diagonal components of the viscosity coefficient can then be
found according to $\xi_{\nu\nu}=-\frac{F_{\nu\nu,(1)}}{v_{\nu}}$, which
yields

$\xi_{\nu\nu}=\frac{\Gamma(\partial_{\nu}h)^{2}}{4\pi}\int
d\omega\frac{\frac{f_{L}}{T_{L}}(1-f_{L})\Gamma_{L}+\frac{f_{R}}{T_{R}}(1-f_{R})\Gamma_{R}}{\Big{[}(\omega-h)^{2}+\Gamma^{2}/4\Big{]}^{2}}.$
(75)

##  Appendix B Exact Diffusion Coefficient

Here we derive (62) for the case of a single electronic level and single
classical degree of freedom $x$. We express the quantum force acting on the
classical coordinate due to the electronic environment according to a mean
component which contains the deterministic forces, and a stochastic component
which captures the probabilistic fluctuations about the mean:

$\hat{f}(t)=f(t)+\delta\hat{f}(t).$ (76)

By multiplying this equation by itself at a different time $t^{\prime}$ and
utilising the fact that $\langle\delta\hat{f}(t)\rangle=0$, we find

$\langle\delta\hat{f}(t)\delta\hat{f}(t^{\prime})\rangle=\langle\hat{f}(t)\hat{f}(t^{\prime})\rangle-f(t)f(t^{\prime}),$
(77)

where the quantum force operator according to (7) is given by

$\hat{f}(t)=-\partial_{x}U-\partial_{x}hd^{\dagger}d.$ (78)

After inputting this expression into (77) and applying Wick’s theorem to
decompose the strings of creation and annihilation operators, we observe that
the terms in $f(t)f(t^{\prime})$ cancel exactly with part of the decomposition
of $\langle\hat{f}(t)\hat{f}(t^{\prime})\rangle$. We are then left with

$\displaystyle\langle\delta\hat{f}(t)\delta\hat{f}(t^{\prime})\rangle$
$\displaystyle=(\partial_{x}h)^{2}\langle
d^{\dagger}(t)d(t^{\prime})\rangle\langle d(t)d^{\dagger}(t^{\prime})\rangle,$
(79) $\displaystyle=\lambda^{2}G^{>}(t,t^{\prime})G^{<}(t^{\prime},t).$ (80)

Here we have used the fact that $\partial_{x}h=\lambda$ for our model and we
have utilised the definitions for the lesser and greater Green’s functions.

##  Appendix C Wigner Transform of (62)

Direct application of the Wigner transform to (62) yields

$\tilde{D}(\omega,T)=\lambda^{2}\int d\tau
e^{i\omega\tau}G^{>}(t,t^{\prime})G^{<}(t^{\prime},t).$ (81)

Now, we replace $G^{<}$ in the time domain by the inverse Wigner transform of
$\tilde{G}^{<}$ and simplify to find

$\displaystyle\tilde{D}(\omega,T)$ $\displaystyle=\frac{\lambda^{2}}{2\pi}\int
d\tau e^{i\omega\tau}G^{>}(t,t^{\prime})\int
d\omega^{\prime}e^{-i\omega^{\prime}(-\tau)}\tilde{G}^{<}(\omega^{\prime},T),$
(82) $\displaystyle=\frac{\lambda^{2}}{2\pi}\int d\omega^{\prime}\int d\tau
e^{i(\omega+\omega^{\prime})\tau}G^{>}(t,t^{\prime})\tilde{G}^{<}(\omega^{\prime},T),$
(83) $\displaystyle=\frac{\lambda^{2}}{2\pi}\int
d\omega^{\prime}\tilde{G}^{>}(\omega+\omega^{\prime},T)\tilde{G}^{<}(\omega^{\prime},T),$
(84)

where we have used the definition of the Wigner transform.

##  Appendix D Equation for the white-noise diffusion

In (66), we introduced the white-noise diffusion $D^{w}$ as a means of
representing the colored-noise diffusion in a Markovian manner. We are then
tasked with finding an equation for $D^{w}$ whose application will accurately
reproduce the dynamics produced by the exact colored-noise diffusion. To do
so, we consider a pedagogical example in which we neglect all forces except
the stochastic force. The governing equation of motion is given by

$m\frac{dv}{dt}=\delta f(t).$ (85)

This is solved for the velocity at time $t$ according to

$v(t)=v(0)+\frac{1}{m}\int^{t}_{0}dt^{\prime}\delta f(t^{\prime}).$ (86)

The change in kinetic energy of the classical coordinate is then found by
squaring the above and taking an average over the fluctuations such that we
obtain

$\Delta
KE=\frac{1}{2m}\int^{t}_{0}dt^{\prime}\int^{t}_{0}dt^{\prime\prime}D(t^{\prime},t^{\prime\prime}).$
(87)

Now, we implement the white noise approximation according to (66) and solve
for an equation for $D^{w}$ which yields the same change in kinetic energy.
Making the transformation from $(t^{\prime},t^{\prime\prime})$ to $(T,\tau)$
yields

$\displaystyle\Delta KE$
$\displaystyle=\frac{1}{2m}\int^{t}_{0}dT\int^{t}_{-t}d\tau
D^{w}(T)\delta(\tau),$ (88)
$\displaystyle=\frac{1}{2m}\int^{t}_{0}dTD^{w}(T).$ (89)

By enforcing that (89) and (87) are equal, we find

$\int^{t}_{0}dTD^{w}(T)=\int^{t}_{0}dT\int^{t}_{-t}d\tau D(T,\tau).$ (90)

Application of $\frac{d}{dt}$ to both sides and using the symmetry property,
$D(T,\tau)=D(T,-\tau)$, then yields

$D^{w}(t)=\int^{t}_{-t}d\tau D(t,\tau)+2\int^{t}_{0}dTD(T,t).$ (91)

Finally, we assume that time $t$ is sufficiently large such that $D(T,t)$ is
approximately zero. This amounts to assuming that the correlations in the
stochastic force decay to zero over time intervals of $t$ or longer. The
second term disappears and we are left with our final expression for the
white-noise diffusion coefficient:

$D^{w}(t)=\int^{t}_{-t}d\tau D(t,\tau).$ (92)

Clearly, the white noise approximation in this case will produce the same
observable change in kinetic energy over any time-scale. However, the
inclusion of a frictional force and an external potential will limit the
validity of the approximation to finite time-scales dependent on their
respective strengths. This is because these forces may produce an appreciable
effect on the dynamics over the time-scales for which the stochastic force is
correlated.

## References

  * Cuevas and Scheer [2017] J. C. Cuevas and E. Scheer,  _Molecular Electronics_, World Scientific Series in Nanoscience and Nanotechnology, Vol. 15 (WORLD SCIENTIFIC, 2017).
  * Ward  _et al._ [2011] D. R. Ward, D. A. Corley, J. M. Tour, and D. Natelson, Vibrational and electronic heating in nanoscale junctions, Nature Nanotechnology 6 (2011).
  * Sabater  _et al._ [2015] C. Sabater, C. Untiedt, and J. M. van Ruitenbeek, Evidence for non-conservative current-induced forces in the breaking of Au and Pt atomic chains, Beilstein J. Nanotechnol. 6, 2338 (2015).
  * Tsutsui  _et al._ [2008] M. Tsutsui, M. Taniguchi, and T. Kawai, Local Heating in Metal-Molecule-Metal Junctions, Nano Letters 8, 3293 (2008).
  * Schulze  _et al._ [2008a] G. Schulze, K. J. Franke, A. Gagliardi, G. Romano, C. S. Lin, A. L. Rosa, T. A. Niehaus, T. Frauenheim, A. Di Carlo, A. Pecchia, and J. I. Pascual, Resonant Electron Heating and Molecular Phonon Cooling in Single ${\mathrm{C}}_{60}$ Junctions, Phys. Rev. Lett. 100, 136801 (2008a).
  * Ioffe  _et al._ [2008] Z. Ioffe, T. Shamai, A. Ophir, G. Noy, I. Yutsis, K. Kfir, O. Cheshnovsky, and Y. Selzer, Detection of heating in current-carrying molecular junctions by Raman scattering, Nat Nano 3, 727 (2008).
  * Franke and Pascual [2012] K. J. Franke and J. I. Pascual, Effects of electron–vibration coupling in transport through single molecules, Journal of Physics: Condensed Matter 24, 394002 (2012).
  * Mukherjee  _et al._ [2013] S. Mukherjee, F. Libisch, N. Large, O. Neumann, L. V. Brown, J. Cheng, J. B. Lassiter, E. A. Carter, P. Nordlander, and N. J. Halas, Hot Electrons Do the Impossible: Plasmon-Induced Dissociation of H2 on Au, Nano Letters 13, 240 (2013).
  * Huang  _et al._ [2006] Z. Huang, B. Xu, Y. Chen, M. D. Ventra, and N. J. Tao, Measurement of Current-Induced Local Heating in a Single Molecule Junction, Nano Letters 6, 1240 (2006).
  * de Leon  _et al._ [2008] N. P. de Leon, W. Liang, Q. Gu, and H. Park, Vibrational Excitation in Single-Molecule Transistors: Deviation from the Simple Franck-Condon Prediction, Nano Letters 8, 2963 (2008).
  * Hüttel  _et al._ [2009] A. K. Hüttel, B. Witkamp, M. Leijnse, M. R. Wegewijs, and H. S. J. van der Zant, Pumping of Vibrational Excitations in the Coulomb-Blockade Regime in a Suspended Carbon Nanotube, Physical Review Letters 102, 225501 (2009).
  * Stettenheim  _et al._ [2010] J. Stettenheim, M. Thalakulam, F. Pan, and Others, A macroscopic mechanical resonator driven by mesoscopic electrical back-action., Nature 466, 86 (2010).
  * Schulze  _et al._ [2008b] G. Schulze, K. J. Franke, and J. I. Pascual, Resonant heating and substrate-mediated cooling of a single C 60 molecule in a tunnel junction, New Journal of Physics 10, 065005 (2008b).
  * Hahn and Ho [2005] J. R. Hahn and W. Ho, Chemisorption and dissociation of single oxygen molecules on Ag(110), The Journal of Chemical Physics 123, 214702 (2005).
  * Roy  _et al._ [2013] S. Roy, V. Mujica, and M. A. Ratner, Chemistry at molecular junctions: Rotation and dissociation of O2 on the Ag(110) surface induced by a scanning tunneling microscope, The Journal of Chemical Physics 139, 74702 (2013).
  * Hla  _et al._ [2000] S.-W. Hla, L. Bartels, G. Meyer, and K.-H. Rieder, Inducing All Steps of a Chemical Reaction with the Scanning Tunneling Microscope Tip: Towards Single Molecule Engineering, Phys. Rev. Lett. 85, 2777 (2000).
  * Lauhon and Ho [2000] L. J. Lauhon and W. Ho, Control and Characterization of a Multistep Unimolecular Reaction, Phys. Rev. Lett. 84, 1527 (2000).
  * Stipe  _et al._ [1997] B. C. Stipe, M. A. Rezaei, W. Ho, S. Gao, M. Persson, and B. I. Lundqvist, Single-Molecule Dissociation by Tunneling Electrons, Phys. Rev. Lett. 78, 4410 (1997).
  * Li  _et al._ [2015] H. Li, T. A. Su, V. Zhang, M. L. Steigerwald, C. Nuckolls, and L. Venkataraman, Electric Field Breakdown in Single Molecule Junctions, Journal of the American Chemical Society 137, 5028 (2015).
  * Li  _et al._ [2016] H. Li, N. T. Kim, T. A. Su, M. L. Steigerwald, C. Nuckolls, P. Darancet, J. L. Leighton, and L. Venkataraman, Mechanism for Si–Si Bond Rupture in Single Molecule Junctions, Journal of the American Chemical Society 138, 16159 (2016).
  * Pobelov  _et al._ [2017] I. V. Pobelov, K. P. Lauritzen, K. Yoshida, A. Jensen, G. Mészáros, K. W. Jacobsen, M. Strange, T. Wandlowski, and G. C. Solomon, Dynamic breaking of a single gold bond, Nature Communications 8, 15931 (2017).
  * Ding  _et al._ [2016] G.-H. Ding, B. Xiong, and B. Dong, Transient currents of a single molecular junction with a vibrational mode, Journal of Physics: Condensed Matter 28, 65301 (2016).
  * Erpenbeck  _et al._ [2016] A. Erpenbeck, R. Härtle, M. Bockstedte, and M. Thoss, Vibrationally dependent electron-electron interactions in resonant electron transport through single-molecule junctions, Physical Review B 93, 115421 (2016).
  * Schinabeck  _et al._ [2014] C. Schinabeck, R. Härtle, H. B. Weber, and M. Thoss, Current noise in single-molecule junctions induced by electronic-vibrational coupling, Phys. Rev. B 90, 75409 (2014).
  * Haupt  _et al._ [2010] F. Haupt, T. c. š. Novotný, and W. Belzig, Current noise in molecular junctions: Effects of the electron-phonon interaction, Phys. Rev. B 82, 165441 (2010).
  * Galperin  _et al._ [2009] M. Galperin, K. Saito, A. V. Balatsky, and A. Nitzan, Cooling mechanisms in molecular conduction junctions, Phys. Rev. B 80, 115427 (2009).
  * Erpenbeck and Thoss [2019] A. Erpenbeck and M. Thoss, Hierarchical quantum master equation approach to vibronic reaction dynamics at metal surfaces, Journal of Chemical Physics 151, 191101 (2019).
  * Koch  _et al._ [2006] J. Koch, M. Semmelhack, F. von Oppen, and A. Nitzan, Current-induced nonequilibrium vibrations in single-molecule devices, Phys. Rev. B 73, 155306 (2006).
  * Wilner  _et al._ [2014] E. Y. Wilner, H. Wang, M. Thoss, and E. Rabani, Nonequilibrium quantum systems with electron-phonon interactions: Transient dynamics and approach to steady state, Phys. Rev. B 89, 205129 (2014).
  * Wang  _et al._ [2011] H. Wang, I. Pshenichnyuk, R. Härtle, and M. Thoss, Numerically exact, time-dependent treatment of vibrationally coupled electron transport in single-molecule junctions, J. Chem. Phys. 135, 244506 (2011).
  * Schinabeck  _et al._ [2016] C. Schinabeck, A. Erpenbeck, R. Härtle, and M. Thoss, Hierarchical quantum master equation approach to electronic-vibrational coupling in nonequilibrium transport through nanosystems, Phys. Rev. B 94, 201407 (2016).
  * Preston  _et al._ [2020a] R. J. Preston, V. F. Kershaw, and D. S. Kosov, Current-induced atomic motion, structural instabilities, and negative temperatures on molecule-electrode interfaces in electronic junctions, Phys. Rev. B 101, 155415 (2020a).
  * Preston  _et al._ [2020b] R. J. Preston, T. D. Honeychurch, and D. S. Kosov, Cooling molecular electronic junctions by AC current, Journal of Chemical Physics 153, 121102 (2020b).
  * Kershaw and Kosov [2020] V. F. Kershaw and D. S. Kosov, Non-adiabatic effects of nuclear motion in quantum transport of electrons: A self-consistent Keldysh–Langevin study, The Journal of Chemical Physics 153, 154101 (2020).
  * Micchi  _et al._ [2016] G. Micchi, R. Avriller, and F. Pistolesi, Electromechanical transition in quantum dots, Phys. Rev. B 94, 125417 (2016).
  * Lu  _et al._ [2010] J.-T. Lu, M. Brandbyge, and P. Hedegard, Blowing the Fuse: Berry’s Phase and Runaway Vibrations in Molecular Conductors, Nano Letters 10, 1657 (2010).
  * Lü  _et al._ [2015] J.-T. Lü, R. B. Christensen, J.-S. Wang, P. Hedegård, and M. Brandbyge, Current-Induced Forces and Hot Spots in Biased Nanojunctions, Physical Review Letters 114, 096801 (2015).
  * Erpenbeck  _et al._ [2018] A. Erpenbeck, C. Schinabeck, U. Peskin, and M. Thoss, Current-induced bond rupture in single-molecule junctions, Phys. Rev. B 97, 235452 (2018).
  * Dzhioev  _et al._ [2013] A. A. Dzhioev, D. S. Kosov, and F. von Oppen, Out-of-equilibrium catalysis of chemical reactions by electronic tunnel currents, J. Chem. Phys. 138, 134103 (2013).
  * Dzhioev and Kosov [2011] A. A. Dzhioev and D. S. Kosov, Kramers problem for nonequilibrium current-induced chemical reactions, J. Chem. Phys. 135, 74701 (2011).
  * Preston  _et al._ [2021] R. J. Preston, M. F. Gelin, and D. S. Kosov, First-passage time theory of activated rate chemical processes in electronic molecular junctions, The Journal of Chemical Physics 154, 114108 (2021).
  * Dou  _et al._ [2017] W. Dou, G. Miao, and J. E. Subotnik, Born-Oppenheimer Dynamics, Electronic Friction, and the Inclusion of Electron-Electron Interactions, Phys. Rev. Lett. 119, 46001 (2017).
  * Dou and Subotnik [2018] W. Dou and J. E. Subotnik, Universality of electronic friction. II. Equivalence of the quantum-classical Liouville equation approach with von Oppen’s nonequilibrium Green’s function methods out of equilibrium, Physical Review B 97, 10.1103/PhysRevB.97.064303 (2018), arXiv:1801.06108 .
  * Bode  _et al._ [2012] N. Bode, S. V. Kusminskiy, R. Egger, and F. von Oppen, Current-induced forces in mesoscopic systems: A scattering-matrix approach, Beilstein Journal of Nanotechnology 3, 144 (2012).
  * Chen  _et al._ [2019a] F. Chen, K. Miwa, and M. Galperin, Current-Induced Forces for Nonadiabatic Molecular Dynamics, The Journal of Physical Chemistry A 123, 693 (2019a).
  * Chen  _et al._ [2019b] F. Chen, K. Miwa, and M. Galperin, Electronic friction in interacting systems, The Journal of Chemical Physics 150, 174101 (2019b).
  * Lü  _et al._ [2012] J.-T. Lü, M. Brandbyge, P. Hedegård, T. N. Todorov, and D. Dundas, Current-induced atomic dynamics, instabilities, and Raman signals: Quasiclassical Langevin equation approach, Physical Review B 85, 245444 (2012).
  * Nocera  _et al._ [2011] A. Nocera, C. A. Perroni, V. Marigliano Ramaglia, and V. Cataudella, Stochastic dynamics for a single vibrational mode in molecular junctions, Physical Review B 83, 115420 (2011).
  * Dundas  _et al._ [2009] D. Dundas, E. J. McEniry, and T. N. Todorov, Current-driven atomic waterwheels, Nature Nanotechnology 4, 99 (2009).
  * Dundas  _et al._ [2012] D. Dundas, B. Cunningham, C. Buchanan, A. Terasawa, A. T. Paxton, and T. N. Todorov, An ignition key for atomic-scale engines, Journal of Physics: Condensed Matter 24, 402203 (2012).
  * Cunningham  _et al._ [2014] B. Cunningham, T. N. Todorov, and D. Dundas, Nonconservative dynamics in long atomic wires, Physical Review B 90, 115430 (2014).
  * Lü  _et al._ [2011a] J.-T. Lü, T. Gunst, P. Hedegård, and M. Brandbyge, Current-induced dynamics in carbon atomic contacts, Beilstein Journal of Nanotechnology 2, 814 (2011a).
  * Christensen  _et al._ [2016] R. B. Christensen, J.-T. Lü, P. Hedegård, and M. Brandbyge, Current-induced runaway vibrations in dehydrogenated graphene nanoribbons, Beilstein Journal of Nanotechnology 7, 68 (2016).
  * Hussein  _et al._ [2010] R. Hussein, A. Metelmann, P. Zedler, and T. Brandes, Semiclassical dynamics of nanoelectromechanical systems, Physical Review B 82, 165406 (2010).
  * Hopjan  _et al._ [2018] M. Hopjan, G. Stefanucci, E. Perfetto, and C. Verdozzi, Molecular junctions and molecular motors: Including Coulomb repulsion in electronic friction using nonequilibrium Green’s functions, Physical Review B 98, 041405 (2018).
  * Kartsev  _et al._ [2014] A. Kartsev, C. Verdozzi, and G. Stefanucci, Nonadiabatic Van der Pol oscillations in molecular transport, The European Physical Journal B 87, 14 (2014).
  * Foti and Vázquez [2018] G. Foti and H. Vázquez, Origin of Vibrational Instabilities in Molecular Wires with Separated Electronic States, The Journal of Physical Chemistry Letters 9, 2791 (2018).
  * Gunst  _et al._ [2013] T. Gunst, J.-T. Lü, P. Hedegård, and M. Brandbyge, Phonon excitation and instabilities in biased graphene nanoconstrictions, Physical Review B 88, 161401 (2013).
  * Simine and Segal [2012] L. Simine and D. Segal, Vibrational cooling, heating, and instability in molecular conducting junctions: full counting statistics analysis, Physical Chemistry Chemical Physics 14, 13820 (2012).
  * Härtle and Thoss [2011] R. Härtle and M. Thoss, Vibrational instabilities in resonant electron transport through single-molecule junctions, Phys. Rev. B 83, 125419 (2011).
  * Lü  _et al._ [2011b] J.-T. Lü, P. Hedegård, and M. Brandbyge, Laserlike Vibrational Instability in Rectifying Molecular Conductors, Physical Review Letters 107, 046801 (2011b).
  * Rizzi  _et al._ [2016] V. Rizzi, T. N. Todorov, J. J. Kohanoff, and A. A. Correa, Electron-phonon thermalization in a scalable method for real-time quantum dynamics, Physical Review B 93, 024306 (2016).
  * Nitzan and Galperin [2018] A. Nitzan and M. Galperin, Kinetic Schemes in Open Interacting Systems, The Journal of Physical Chemistry Letters 9, 4886 (2018).
  * Hyldgaard [2012] P. Hyldgaard, Nonequilibrium thermodynamics of interacting tunneling transport: variational grand potential, density functional formulation and nature of steady-state forces, Journal of Physics: Condensed Matter 24, 424219 (2012).
  * Todorov  _et al._ [2010] T. N. Todorov, D. Dundas, and E. J. McEniry, Nonconservative generalized current-induced forces, Phys. Rev. B 81, 75416 (2010).
  * Verdozzi  _et al._ [2006] C. Verdozzi, G. Stefanucci, and C.-O. Almbladh, Classical Nuclear Motion in Quantum Transport, Phys. Rev. Lett. 97, 46603 (2006).
  * Brandbyge  _et al._ [2003] M. Brandbyge, K. Stokbro, J. Taylor, J. L. Mozos, and P. Ordejon, Origin of current-induced forces in an atomic gold wire: A first-principles study, Phys. Rev. B 67, 193104 (2003).
  * Metelmann and Brandes [2011] A. Metelmann and T. Brandes, Adiabaticity in semiclassical nanoelectromechanical systems, Physical Review B 84, 155455 (2011).
  * Horsfield  _et al._ [2004a] A. P. Horsfield, D. R. Bowler, A. J. Fisher, T. N. Todorov, and M. J. Montgomery, Power dissipation in nanoscale conductors: classical, semi-classical and quantum dynamics, Journal of Physics: Condensed Matter 16, 3609 (2004a).
  * McEniry  _et al._ [2007] E. J. McEniry, D. R. Bowler, D. Dundas, A. P. Horsfield, C. G. Sánchez, and T. N. Todorov, Dynamical simulation of inelastic quantum transport, Journal of Physics: Condensed Matter 19, 196201 (2007).
  * Horsfield  _et al._ [2004b] A. P. Horsfield, D. R. Bowler, A. J. Fisher, T. N. Todorov, and C. G. Sánchez, Beyond Ehrenfest: correlated non-adiabatic molecular dynamics, Journal of Physics: Condensed Matter 16, 8251 (2004b).
  * Nitzan [2006] A. Nitzan,  _Chemical Dynamics in Condensed Phases_ (Oxford University Press, 2006).
  * Todorov  _et al._ [2014] T. N. Todorov, D. Dundas, J.-T. Lü, M. Brandbyge, and P. Hedegård, Current-induced forces: a simple derivation, European Journal of Physics 35, 065004 (2014).
  * Lü  _et al._ [2019] J.-T. Lü, B.-Z. Hu, P. Hedegård, and M. Brandbyge, Semi-classical generalized Langevin equation for equilibrium and nonequilibrium molecular dynamics simulation, Progress in Surface Science 94, 21 (2019).
  * Jauho  _et al._ [1994] A. P. Jauho, N. S. Wingreen, and Y. Meir, Time-dependent transport in interacting and noninteracting resonant-tunneling systems, Phys. Rev. B 50, 5528 (1994).
  * McQuarrie [2008] D. A. McQuarrie,  _Quantum Chemistry_ , 2nd ed. (University Science Books, 2008).