# Current-induced forces in nanosystems: A hierarchical equations of motion
approach

Samuel L. Rudge  Institute of Physics, University of Freiburg, Hermann-Herder-
Strasse 3, 79104 Freiburg, Germany  Yaling Ke  Institute of Physics,
University of Freiburg, Hermann-Herder-Strasse 3, 79104 Freiburg, Germany
Michael Thoss  Institute of Physics, University of Freiburg, Hermann-Herder-
Strasse 3, 79104 Freiburg, Germany

###### Abstract

A new approach to calculating current-induced forces in charge transport
through nanosystems is introduced. Starting from the fully quantum mechanical
hierarchical equations of motion formalism, a timescale separation between
electronic and vibrational degrees of freedom is used to derive a classical
Langevin equation of motion for the vibrational dynamics as influenced by
current-induced forces, such as the electronic friction. The resulting form of
the friction is shown to be equivalent to previously derived expressions. The
numerical exactness of the hierarchical equations of motion approach, however,
allows the investigation of transport scenarios with strong intrasystem and
system-environment interactions. As a demonstration, the electronic friction
of three example systems is calculated and analyzed: a single electronic level
coupled to one classical vibrational mode, two electronic levels coupled to
one classical vibrational mode, and a single electronic level coupled to both
a classical and quantum vibrational mode.

##  I Introduction

A comprehensive picture of nonequilibrium charge transport through
nanosystems, such as molecular junctions Tal  _et al._ (2008); Secker  _et
al._ (2011), nanoelectromechanical systems Craighead (2000); Ekinci and Roukes
(2005), or suspended carbon nanotubes Weig  _et al._ (2004); Sapmaz  _et al._
(2006), requires not only a description of charge effects but also the
coupling of electronic and vibrational degrees of freedom. Novel effects
arising from such electron-vibration interactions include the Franck-Condon
blockade Koch and von Oppen (2005); Koch  _et al._ (2006a, b); Leturcq  _et
al._ (2009), current-induced heating and cooling Härtle and Thoss (2011a);
Härtle  _et al._ (2018); Härtle and Thoss (2011b); Kaspar  _et al._ (2022),
and non-renewal statistics Koch and von Oppen (2005); Kosov (2017); Rudge and
Kosov (2019a); Schinabeck  _et al._ (2014). Including such interactions,
however, requires transport methods able to treat the quantum nature of the
nanosystem while simultaneously incorporating nonequilibrium effects from the
electrodes. Examples of such methods include nonequilibrium Green’s functions
(NEGFs) Mitra  _et al._ (2004); Frederiksen  _et al._ (2007); Galperin  _et
al._ (2007); Härtle  _et al._ (2008); Novotný  _et al._ (2011); Dash  _et al._
(2011); Agarwalla  _et al._ (2015); Ryndyk  _et al._ (2006); Martin-Rodero
_et al._ (2008); Entin-Wohlman  _et al._ (2009); Egger and Gogolin (2008);
Galperin  _et al._ (2004); Tikhodeev  _et al._ (2001); Paulsson  _et al._
(2005), which can treat strong molecule-lead couplings, and Born-Markov master
equations (BMMEs) Koch  _et al._ (2006a); Li  _et al._ (2005); Leijnse  _et
al._ (2010); Rudge and Kosov (2021), which can treat strong intra-system
interactions. In recent years, furthermore, many more sophisticated and
numerically exact methods have also been developed, including path integral
approaches Mühlbacher and Rabani (2008); Weiss  _et al._ (2008), the
multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) approach Wang
_et al._ (2011); Wang and Thoss (2016) or the hierarchical equations of motion
(HEOM) approach Schinabeck  _et al._ (2016, 2018); Schinabeck and Thoss
(2020); Jin  _et al._ (2008), although these can be computationally
challenging.

In the specific context of transport through molecular junctions, many
vibrational degrees of freedom, each corresponding to the motion of atomic
nuclei, are necessary for realistic simulations. In these scenarios, mixed
quantum-classical approaches, where the vibrational dynamics are treated
classically while subject to current-induced forces from the quantum
electronic degrees of freedom, are often more viable Tully (1990); Stock and
Thoss (1997); Martinez  _et al._ (1996). When the nuclear motion is slow or
limited to small displacements, furthermore, the vibrational degrees of
freedom follow a Langevin equation of motion. In such regimes, the current-
induced force splits into a zeroth-order term, corresponding to the Born-
Oppenheimer approximation, and a first-order correction linear in the
vibrational momenta or coordinates Dou and Subotnik (2018): the electronic
friction, which is the specific current-induced force of interest in this
study. Because the current-induced forces offer insight into the role of the
electron-vibration interaction on the dynamics Dzhioev and Kosov (2011);
Dzhioev  _et al._ (2013); Preston  _et al._ (2021); Mozyrsky  _et al._ (2006);
Kershaw and Kosov (2020); Bode  _et al._ (2012); Lü  _et al._ (2010, 2012);
Todorov  _et al._ (2011, 2014); Maurer  _et al._ (2016); Askerka  _et al._
(2016); Bustos-Marún  _et al._ (2013); Fernández-Alcázar  _et al._ (2019);
Calvo  _et al._ (2017), especially in the regime of vibrational instabilities
Sabater  _et al._ (2015); Ring  _et al._ (2020), various theories have been
developed for calculating them in the context of quantum transport through
nanostructures.

These approaches can be broadly categorized into two classes: those that
derive the electronic friction under a weak electron-vibration coupling, and
those that derive it using a timescale separation between electronic and
vibrational degrees of freedom. To the first class belongs an influence
functional approach that has been derived by several groups Hussein  _et al._
(2010); Metelmann and Brandes (2011); Lü  _et al._ (2012, 2010, 2011);
Brandbyge  _et al._ (1995); Caldeira and Leggett (1983); Schmid (1982). The
second class contains a wider variety of transport methods. These include the
approach introduced in the seminal paper by Head-Gordon and Tully Head‐Gordon
and Tully (1995); Stock and Thoss (1997); Meyer and Miller (1979), which has
since been extended to finite temperatures Dou  _et al._ (2017); Maurer  _et
al._ (2016) and has shown to be equivalent to a combined scattering theory and
NEGF approach Bode  _et al._ (2012); Dou and Subotnik (2017a). The connections
between these various approaches and assumptions have been explicitly explored
by Chen et al. Chen  _et al._ (2019a), who also used their unified formalism
to investigate interacting systems Chen  _et al._ (2019b). Although these
approaches explicitly treat the vibrational degrees of freedom as classical,
there have also been recent attempts to derive a quantum equivalent to the
classical electronic friction Martinazzo and Burghardt (2022a, b).

Under the timescale separation assumption, one can also formulate master
equation approaches, such as that proposed by Dou and Subotnik in Refs. Dou
_et al._ (2015); Dou and Subotnik (2016). Here, a quantum-classical Liouville
equation (QCLE) for the vibrational dynamics was embedded into a BMME for the
electronic dynamics. This approach has the benefit of allowing for strong
intra-system interactions, but is limited to weak molecule-lead couplings
under the Born-Markov assumption. This method is particularly relevant to the
purpose of this paper, in which a similar approach to calculating the
electronic friction is introduced, although with the numerically exact HEOM
approach instead of a BMME. Starting from the full HEOM method, in which the
vibrational degrees of freedom are originally treated quantum mechanically,
the QCLE is formed by applying a partial Wigner transform of the vibrational
degrees of freedom Kapral and Ciccotti (1999); İmre  _et al._ (1967). Similar
to Ref. Dou and Subotnik (2016), the QCLE is then transformed to a Fokker-
Planck equation in the adiabatic limit of slow vibrational motion, from which
the Markovian electronic friction of the corresponding Langevin equation can
be identified. The resulting electronic friction is then shown to be
equivalent to expressions derived by treating all electronic degrees of
freedom in the system and bath equally Dou and Subotnik (2018); Chen  _et al._
(2019a). Because it is derived from HEOM, which systematically incorporates
higher-order tunneling and non-Markovian effects for nonequilibrium charge
transport, the results are general in- and out-of-equilibrium for systems with
strong interactions and strong molecule-lead couplings.

In order to demonstrate the approach, three model systems are considered. In
the first system, a single electronic level is linearly coupled to a classical
vibrational degree of freedom. The electronic friction of this non-interacting
system has previously been calculated from NEGFs Dou and Subotnik (2018),
which the HEOM approach reproduces. In the second system, electron-electron
interactions are incorporated via an Anderson impurity, and the friction at
equilibrium is shown to be equivalent to the combined numerical
renormalization group (NRG) and NEGF method from Ref. Dou  _et al._ (2017).
The HEOM approach, however, can also investigate this model out-of-
equilibrium, and the corresponding friction at finite bias voltages is
analyzed. The final case study couples an electronic level to both a high-
frequency mode, treated quantum mechanically in the HEOM framework, and a low-
frequency mode, treated classically in the Langevin equation framework. This
system is similar to one treated recently in Ref. Chen  _et al._ (2019b),
where the high-frequency mode was used to model a strong light-matter
interaction. The additional vibrational degree of freedom adds extra structure
to the friction on the classical mode in nonequilibrium scenarios, even
turning it negative for particular parameters.

The paper is structured as follows. In Section II, the general model of a
molecular junction is outlined. Next, in Section III, the two methods used to
describe nonequilibrium charge transport are introduced, the HEOM and BMME
approaches. In Section IV, these are then Wigner-transformed to a QCLE and,
under the appropriate assumptions, to a Markovian Fokker-Planck equation, from
which the friction can be identified. In order to demonstrate the method,
Section V then presents results for the three model systems, with the
conclusions contained in Section VI. Relevant derivations can be found in
Appendices A, B, and C.

Throughout the paper we use units where $\hbar=e=k_{B}=1$.

##  II Model

The total Hamiltonian of a junction containing a nanosystem, $H$, can be
decomposed into a part describing the electrons in the leads, $H_{\text{B}}$,
a part containing the nanosystem, $H_{\text{S}}$, and a part containing the
interaction between them, $H_{\text{SB}}$:

$\displaystyle H$ $\displaystyle=H_{\text{S}}+H_{\text{B}}+H_{\text{SB}}.$ (1)

The system Hamiltonian contains several coupled electronic and vibrational
degrees of freedom, which can generally be written as

$\displaystyle H_{\text{S}}$
$\displaystyle=\sum_{mm^{\prime}}\varepsilon_{mm^{\prime}}(\hat{\mathbb{x}})d_{m}^{{\dagger}}d_{m^{\prime}}+\sum_{m,m^{\prime}<m}U_{mm^{\prime}}(\hat{\mathbb{x}})d_{m}^{{\dagger}}d_{m}d_{m^{\prime}}^{{\dagger}}d_{m^{\prime}}$
$\displaystyle\>\>\>\>\>+\sum_{i=1}^{N}\frac{\hat{p}^{2}_{i}}{2m_{i}}+V_{\text{vib}}(\hat{\mathbb{x}}),$
(2)

where
$\left(\hat{\mathbb{x}}=(\hat{x}_{1},\dots,\hat{x}_{N}),\hat{\mathbb{p}}=(\hat{p}_{1},\dots,\hat{x}_{N})\right)$
are vectors of displacement operators and their corresponding conjugate
momenta for each of the $N$ vibrational coordinates. The potential of the
unperturbed vibrational degrees of freedom is
$V_{\text{vib}}(\hat{\mathbb{x}})$. The vibrational motion is coupled to the
electronic degrees of freedom via
$\varepsilon_{mm^{\prime}}(\hat{\mathbb{x}})$, which is the energy of the
$m$th electronic energy level if $m=m^{\prime}$, or the direct hopping between
the $m$th and $m^{\prime}$th electronic energy levels if $m\neq m^{\prime}$,
and $U_{mm^{\prime}}(\hat{\mathbb{x}})$, which is the electron-electron
interaction between electrons in the $m$th and $m^{\prime}$th electronic
energy levels. The annihilation and creation operators corresponding to the
$m$th electronic energy level are $d_{m}$ and $d^{{\dagger}}_{m}$,
respectively. Note that the vibrational operators have been explicitly written
with operator notation to distinguish $\hat{\mathbb{x}}$ and
$\hat{\mathbb{p}}$ from their classical variable counterparts, $\mathbb{x}$
and $\mathbb{p}$, which will be introduced in Sec. IV.1.

Although the interest in current-induced forces and mixed quantum-classical
dynamics is often from the perspective of transport through molecular
junctions, the Hamiltonian in Eq.(2) and the theory that follows in Secs. III
and IV is general for transport through any generic nanosystem.

The leads are modeled as two reservoirs of non-interacting electrons, labeled
the left, $L$, and right, $R$, electrodes,

$\displaystyle H_{\text{B}}$
$\displaystyle=\sum_{\alpha\in\\{L,R\\}}\sum_{k_{\alpha}}\varepsilon_{k_{\alpha}}c^{{\dagger}}_{k_{\alpha}}c_{k_{\alpha}}.$
(3)

The $(k_{\alpha})$th level is formed through the corresponding annihilation
and creation operators, $c_{k_{\alpha}}$ and $c^{{\dagger}}_{k_{\alpha}}$, and
has energy $\varepsilon_{k_{\alpha}}$. Each electrode is held at local
equilibrium, with temperatures $T_{\alpha}$ and chemical potentials
$\mu_{\alpha}$. By introducing a temperature difference, $\Delta
T=T_{L}-T_{R}$, or voltage bias, $\Phi=\mu_{L}-\mu_{R}$, across the junction,
the system can be driven out of equilibrium, producing a non-zero total
current as electrons transport across the junction.

The interaction term in the total Hamiltonian, which is linear in both the
system and bath annihilation and creation operators, describes electron
transfer between the nanosystem and the electrodes,

$\displaystyle H_{\text{SB}}$
$\displaystyle=\sum_{\alpha,k_{\alpha}}\sum_{m}V_{k_{\alpha},m}\left(c^{{\dagger}}_{k_{\alpha}}d_{m}+d^{{\dagger}}_{m}c_{k_{\alpha}}\right),$
(4)

where the molecule-lead couplings, $V_{k_{\alpha},m}$, are assumed to be real
and independent of the vibrational positions. Although this assumption is
standard in the literature and greatly simplifies the quantum and semi-
classical transport theory, a more general approach would include
$V_{k_{\alpha},m}(\hat{\mathbf{x}})$, such as in Ref. Dou and Subotnik
(2017b), which will be incorporated in future work. Finally, the spectral
density of bath $\alpha$ is introduced,

$\displaystyle\Gamma_{\alpha,mm^{\prime}}(\varepsilon)$
$\displaystyle=2\pi\sum_{k_{\alpha}}V_{k_{\alpha},m}V_{k_{\alpha},m^{\prime}}\delta(\varepsilon-\varepsilon_{k_{\alpha}}),$
(5)

which describes the interaction strength between electronic levels in the
nanosystem and electrodes.

##  III Quantum transport theory

This section outlines the two transport methods used to simulate the
nonequilibrium quantum dynamics of the general Hamiltonian in Eq.(2). Although
both methods have been used extensively, and there exist several detailed
reviews, brief outlines are included here for self-consistency within the
paper.

###  III.1 Hierarchical equations of motion

The HEOM formalism is a powerful method able to simulate nonequilibrium
quantum transport with strong interactions and beyond the Born-Markov regime.
Although the general steps are outlined in this section, readers interested in
a comprehensive review are referred to the original work by Tanimura and Kubo
Tanimura and Kubo (1989), its extension to fermionic systems Tanimura (2006);
Härtle  _et al._ (2015); Jin  _et al._ (2007, 2008); Zheng  _et al._ (2009);
Yan (2014); Wenderoth  _et al._ (2016), and the review articles in Refs.
Tanimura (2020); Ye  _et al._ (2016).

As with all master equation methods, the central object of interest in the
HEOM approach is the reduced density matrix of the nanosystem, $\rho(t)$,
which is obtained by tracing out the bath degrees of freedom from the total
density matrix, $\rho(t)=\text{Tr}_{\text{B}}\left[\rho_{\text{T}}(t)\right]$.
One could obtain $\rho(t)$ by solving the corresponding Liouville-von Neumann
equation,

$\displaystyle\frac{\partial}{\partial t}\rho(t)$
$\displaystyle=-i\text{Tr}_{\text{B}}\left[H,\rho_{\text{T}}(t)\right].$ (6)

This, however, is a difficult task; the dynamics involves not only the
coherent time-evolution of the nanosystem, but non-trivial dissipation due to
interaction with the electrodes. By employing a coherent state representation
for the fermionic annihilation and creation operators, the HEOM approach
treats the system-bath interaction via the Feynman-Vernon influence
functional, which collects the bath effects into an effective action term: the
influence functional. In this work, the influence functional is derived under
the assumption of a factorized initial state,

$\displaystyle\rho_{\text{T}}(0)$
$\displaystyle=\rho(0)\otimes\rho_{\text{B}}(0),$ (7)

as initial correlations between the nanosystem and the bath require a more
complicated imaginary-time propagation approach Tanimura (2014); Ke  _et al._
(2022).

As long as the molecule-lead coupling is linear in the bath annihilation and
creation operators, one can then show that a cumulant expansion of the
influence functional terminates exactly after the second term. The bath
statistics are therefore Gaussian and their effect is exactly described with
the two-time bath correlation functions,

$\displaystyle C^{\sigma}_{\alpha,mm^{\prime}}(t-\tau)$
$\displaystyle=\sum_{k_{\alpha}}V_{k_{\alpha},m}V_{k_{\alpha},m^{\prime}}\text{Tr}_{\text{B}}\left[c^{\sigma}_{k_{\alpha}}(t)c^{\bar{\sigma}}_{k_{\alpha}}(\tau)\rho_{\text{B}}(0)\right],$
(8)

where the notation $\sigma\in\\{+,-\\}$, with
$c^{-(+)}_{k_{\alpha}}=c^{({\dagger})}_{k_{\alpha}}$, has been introduced.
Based on the spectral densities of the baths, the bath-correlation functions
can be rewritten as

$\displaystyle C^{\sigma}_{\alpha,mm^{\prime}}(t)$
$\displaystyle=\frac{1}{2\pi}\int^{\infty}_{-\infty}d\omega\>e^{\sigma i\omega
t}\Gamma_{\alpha,mm^{\prime}}(\omega)f^{\sigma}_{\alpha}(\omega),$ (9)

with the Fermi-Dirac functions,

$\displaystyle f^{\sigma}_{\alpha}(\omega)$
$\displaystyle=\frac{1}{1+e^{\sigma\beta(w-\mu_{\alpha})}},$ (10)

and $\beta=1/k_{B}T$.

The next key step in deriving the HEOM is to assume that the bath-correlation
functions can be expanded in a series of exponential functions,

$\displaystyle C^{\sigma}_{\alpha,mm^{\prime}}(t)$
$\displaystyle=V_{\alpha,m}V_{\alpha,m^{\prime}}\sum_{\ell=0}^{\ell_{\max}}\eta_{\alpha,\sigma,\ell,m}e^{-\kappa_{\alpha,\sigma,\ell,m}t}.$
(11)

By incorporating a sum-over-poles decomposition scheme for the Fermi-Dirac
functions and choosing a specific form of the bath spectral densities, the
coefficients, $\eta_{\alpha,\sigma,\ell,m}$, and exponents,
$\kappa_{\alpha,\sigma,\ell,m}$, can then be calculated directly with residue
theory from Eq.(9). In this paper, the Padé decomposition is exclusively used,
which performs much bettter than a standard Matsubara decomposition Ozaki
(2007). All results below are also calculated under the wide-band limit, such
that
$\Gamma_{\alpha,mm^{\prime}}(\varepsilon)=\Gamma_{\alpha,mm^{\prime}}=2\pi
V_{\alpha,m}V_{\alpha,m^{\prime}}$, where the constant density of states has
been absorbed into the definition of $V_{\alpha,m}$. Although the wideband
limit is a physically reasonable assumption for many systems of interest, such
as molecules attached to metal electrodes Verzijl  _et al._ (2013), it is not
actually necessary for the HEOM approach. The decomposition of the bath-
correlation functions in Eq.(9) can be performed for a wide variety of
$\Gamma_{\alpha,mm^{\prime}}(\varepsilon)$, particularly with the recent
development of numerical schemes able to treat density of state functions
without an explicit analytic decomposition Xu  _et al._ (2022); Dan  _et al._
(2022).

By repeatedly differentiating the bath-correlation functions, a hierarchy of
first-order differential equations coupling the reduced system density matrix
to a series of auxiliary density operators (ADOs), $\rho^{(n)}_{\mathbf{j}}$,
is formed,

$\displaystyle\frac{\partial}{\partial t}\rho^{(n)}_{\mathbf{j}}(t)$
$\displaystyle=-i\left[H_{\text{S}},\rho^{(n)}_{\mathbf{j}}\right]-\left(\sum_{r=1}^{N}\kappa_{j_{r}}\right)\rho^{(n)}_{\mathbf{j}}$
$\displaystyle-i\sum_{r=1}^{n}(-1)^{n-r}\mathcal{C}_{j_{r}}\rho^{(n-1)}_{\mathbf{j}^{-}}-i\sum_{j}\mathcal{A}^{\bar{\sigma}}_{\alpha,m}\rho^{(n+1)}_{\mathbf{j}^{+}}.$
(12)

The notation here uses a vector of condensed super-indices,
$\mathbf{j}=\left(j_{n}\dots j_{1}\right)$ with
$j_{r}=\\{\alpha_{j_{r}},\sigma_{j_{r}},\ell_{j_{r}},m_{j_{r}}\\}$. The ADO
tier, which is given the label $n$, is then the number of indices required to
form the $\rho^{(n)}_{\mathbf{j}}$ ADOs. In this formulation, the $0$th-tier
ADO is then just the reduced density matrix of the nanosystem,
$\rho^{(0)}(t)=\rho(t)$. The HEOM also contain two more super-indices,
$\mathbf{j}^{-}=\left(j_{n},\dots,j_{r+1},j_{r-1},\dots,j_{1}\right)$, which
is formed by removing an ADO, and
$\mathbf{j}^{+}=\left(j,j_{n},\dots,j_{1}\right)$, which is formed by adding
one. The superoperators $\mathcal{C}_{j_{r}}$ and
$\mathcal{A}^{\bar{\sigma}}_{\alpha,m}$, meanwhile, are defined by their
action on the ADOs,

$\displaystyle\mathcal{C}_{j_{r}}\rho^{(n)}_{\mathbf{j}}(t)$
$\displaystyle=V_{\alpha,m}\left(\eta_{j_{r}}d^{\sigma}_{m}\rho^{(n)}_{\mathbf{j}}(t)-(-1)^{n}\eta^{*}_{j_{r}}\rho^{(n)}_{\mathbf{j}}d^{\sigma}_{m}(t)\right)$
(13)
$\displaystyle\mathcal{A}^{\bar{\sigma}}_{\alpha,m}\rho^{(n)}_{\mathbf{j}}(t)$
$\displaystyle=V_{\alpha,m}\left(d^{\bar{\sigma}}_{m}\rho^{(n)}_{\mathbf{j}}(t)+(-1)^{n}\rho^{(n)}_{\mathbf{j}}d^{\bar{\sigma}}_{m}(t)\right).$
(14)

Eq.(12) represents a large set of coupled $1$st-order differential equations.
In order to generate the hierarchy, one first constructs the equation of
motion for the reduced density matrix of the nanosystem, which couples
according to Eq.(14) to the $1$st-tier ADOs, $\rho^{(1)}_{j}(t)$. The
respective $1$st-tier equations of motion are then generated from Eq.(12),
which couple according to Eqs.(13) and (14) to the reduced density matrix and
$2$nd-tier ADOs, $\rho^{(2)}_{j_{2}j_{1}}$, and reduced density matrix,
respectively. This process continues until the hierarchy naturally terminates
at some, typically large, tier. In practice, the hierarchy must be truncated
before this point due to numerical constraints, but at a large enough tier
such that the dynamics are converged. In this work, a simple scheme is used
that truncates the hierarchy directly at a maximum tier, $N$, such that the
$\\{\rho^{(N)}_{\mathbf{j}}(t)\\}$ couple only to ADOs of the same tier and
ADOs of the tier below, $\\{\rho^{(N-1)}_{\mathbf{j}^{-}}(t)\\}$. The
resulting coupled set of $1$st-order differential equations, with the
appropriate initial condition, is solved with a $4$th-order Runge-Kutta
method.

###  III.2 Born-Markov master equation

In order to compare to the method in Ref. Dou and Subotnik (2016), a Born-
Markov master equation is also used to calculate the electronic friction.
Unlike the numerically exact HEOM method, Born-Markov theory expands the
Liouville-von Neumann equation to second-order in the system-bath interaction
and explicitly assumes that bath-correlations decay on a timescale quicker
than the dynamics of the nanosystem Li  _et al._ (2005); Rudge and Kosov
(2021). The resulting equation of motion for the reduced system density
matrix,

$\displaystyle\frac{\partial}{\partial t}\rho(t)$
$\displaystyle=-i\left[H_{\text{S}},\rho(t)\right]$
$\displaystyle\>\>\>\>\>\>-\int^{\infty}_{0}d\tau\>\text{Tr}_{\text{B}}\left[H_{\text{SB}},\left[H_{\text{SB}}(\tau),\rho(t)\otimes\rho^{\text{eq}}_{\text{B}}\right]\right],$
(15)

is closed and, therefore, easier to treat numerically, but it is only valid in
the weak molecule-lead coupling regime. In Eq.(15), the system-bath
interaction follows the time-evolution
$H_{\text{SB}}(\tau)=e^{-i(H_{\text{S}}+H_{\text{B}})\tau}H_{\text{SB}}e^{i(H_{\text{S}}+H_{\text{B}})\tau}$.
In all calculations, the BMME will be further simplified by excluding the off-
diagonal elements of the density matrix from the transport, yielding
essentially a rate equation between pure states of the nanosystem Rudge and
Kosov (2019b).

##  IV Quantum-classical transport theory

###  IV.1 Quantum-classical Liouville formulation of the HEOM

In this subsection, the fully quantum HEOM approach is transformed into a
quantum-classical Liouville equation. First, the reduced system density matrix
is partially Wigner transformed with respect to the vibrational degrees of
freedom Kapral and Ciccotti (1999); Dou and Subotnik (2016),

$\displaystyle\rho^{\text{el}}(\mathbf{x},\mathbf{p};t)$
$\displaystyle=(2\pi)^{-3N}\int
d\mathbf{y}\>e^{i\mathbf{p}\cdot\mathbf{y}}\langle\mathbf{x}-\frac{\mathbf{y}}{2}|\rho(t)|\mathbf{x}+\frac{\mathbf{y}}{2}\rangle.$
(16)

Note that the Wigner transformed quantities, such as
$H^{\text{el}}_{S}(\mathbf{x},\mathbf{p})$, are evaluated at positions and
momenta, $(\mathbf{x},\mathbf{p})$, but are still operators in the electronic
subspace of the nanosystem.

Partially Wigner transforming the first of the coupled differential equations
in the HEOM yields

$\displaystyle\frac{\partial}{\partial
t}\rho^{\text{el}}(\mathbf{x},\mathbf{p},t)$
$\displaystyle=-i\left(\left[H^{\text{el}}_{\text{S}},\rho(t)\right]\right)^{\text{el}}(\hat{\mathbf{x}},\hat{\mathbf{p}})$
$\displaystyle\>\>\>\>\>-i\sum_{j}\mathcal{A}^{\bar{\sigma}}_{K}\rho^{(1),\text{el}}_{j}(\mathbf{x},\mathbf{p};t),$
(17)

where the partial Wigner transform of an $n$th-tier ADO is

$\displaystyle\rho^{(n),\text{el}}_{\mathbf{j}}(\mathbf{x},\mathbf{p};t)$
$\displaystyle=(2\pi)^{-3N}\int
d\mathbf{y}\>e^{i\mathbf{p}\cdot\mathbf{y}}\langle\mathbf{x}-\frac{\mathbf{y}}{2}|\rho^{(n)}_{\mathbf{j}}(t)|\mathbf{x}+\frac{\mathbf{y}}{2}\rangle.$
(18)

The crucial step that transforms this into a quantum-classical equation comes
when evaluating the partial Wigner transform of a product of operators İmre
_et al._ (1967), which is

$\displaystyle(AB)^{\text{el}}(\mathbf{x},\mathbf{p})$
$\displaystyle=A^{\text{el}}(\mathbf{x},\mathbf{p})e^{\hbar\Lambda/2i}B^{\text{el}}(\mathbf{x},\mathbf{p}),$
(19)

with the Poisson bracket operator

$\displaystyle\Lambda$
$\displaystyle=\overleftarrow{\nabla}_{\mathbf{p}}\cdot\overrightarrow{\nabla}_{\mathbf{x}}-\overleftarrow{\nabla}_{\mathbf{x}}\cdot\overrightarrow{\nabla}_{\mathbf{p}}.$
(20)

The direction in which each operator acts is indicated by the arrows. In
expanding the exponential, a classical approximation for the vibrational
degrees of freedom is obtained by keeping only terms up to linear in $\hbar$,
such that the commutator in Eq.(17) becomes

$\displaystyle\left(\left[H^{\text{el}}_{\text{S}},\rho(t)\right]\right)^{\text{el}}(\hat{\mathbf{x}},\hat{\mathbf{p}})\approx$
$\displaystyle\hskip
14.22636pt\left[H^{\text{el}}_{\text{S}}(\mathbf{x},\mathbf{p}),\rho^{\text{el}}(\mathbf{x},\mathbf{p};t)\right]+i\left\\{H^{\text{el}}_{\text{S}}(\mathbf{x},\mathbf{p}),\rho^{\text{el}}(\mathbf{x},\mathbf{p};t)\right\\}_{a},$
(21)

where
$\\{A_{1},A_{2}\\}_{a}=\frac{1}{2}\left(\\{A_{1},A_{2}\\}-\\{A_{2},A_{1}\\}\right)$
and the Poisson bracket is

$\displaystyle\\{A_{1},A_{2}\\}$
$\displaystyle=\sum_{i=1}^{N}\left(\frac{\partial A_{1}}{\partial
x_{i}}\frac{\partial A_{2}}{\partial p_{i}}-\frac{\partial A_{1}}{\partial
p_{i}}\frac{\partial A_{2}}{\partial x_{i}}\right).$ (22)

Overall, then, Eq.(17) becomes

$\displaystyle\frac{\partial}{\partial
t}\rho^{\text{el}}(\mathbf{x},\mathbf{p};t)$
$\displaystyle=\left\\{H^{\text{el}}_{\text{S}}(\mathbf{x},\mathbf{p}),\rho^{\text{el}}(\mathbf{x},\mathbf{p};t)\right\\}_{a}$
$\displaystyle\>\>\>\>\>-i\left[H^{\text{el}}_{\text{S}}(\mathbf{x},\mathbf{p}),\rho^{\text{el}}(\mathbf{x},\mathbf{p};t)\right]$
$\displaystyle\>\>\>\>\>-i\sum_{j}\mathcal{A}^{\bar{\sigma}}_{\alpha,m}\rho^{(1),\text{el}}_{j}(\mathbf{x},\mathbf{p};t).$
(23)

Each of the coupled equations in the hierarchy is then partially Wigner
transformed, yielding in general

$\displaystyle\frac{\partial}{\partial
t}\rho^{(n),\text{el}}_{\mathbf{j}}(\mathbf{x},\mathbf{p};t)$
$\displaystyle=\left\\{H^{\text{el}}_{\text{S}}(\mathbf{x},\mathbf{p}),\rho^{(n),\text{el}}_{\mathbf{j}}\right\\}_{a}$
$\displaystyle\>\>\>\>\>\>\>-i\left[H^{\text{el}}_{\text{S}}(\mathbb{x}),\rho^{(n),\text{el}}_{\mathbf{j}}\right]-\left(\sum_{r=1}^{N}\kappa_{j_{r}}\right)\rho^{(n)}_{\mathbf{j}}$
$\displaystyle\>\>\>\>\>\>\>-i\sum_{r=1}^{n}(-1)^{n-r}\mathcal{C}_{j_{r}}\rho^{(n-1),\text{el}}_{\mathbf{j}^{-}}$
$\displaystyle\>\>\>\>\>\>\>-i\sum_{j}\mathcal{A}^{\bar{\sigma}}_{\alpha,m}\rho^{(n+1),\text{el}}_{\mathbf{j}^{+}}.$
(24)

In Eq.(24), the partial Wigner transform of the part connecting the $n$th-tier
ADO to the $(n+1)$th- and $(n-1)$th-tier ADOs only affects the ADOs
themselves, as the superoperators $\mathcal{A}$ and $\mathcal{C}$ do not
depend on the vibrational coordinates. It is also noted that, although the
momenta enter $H^{\text{el}}_{S}(\mathbf{x},\mathbf{p})$ through the kinetic
energies of the nuclei, they disappear in the commutator with the electronic
density matrix. This means that the HEOM part of Eq.(24), that is, everything
except for the Poisson bracket, depends only on the vibrational coordinates,
$\mathbf{x}$.

The reduced density matrix and all ADOs are next unfolded into vector format
and collected into one joint vector,
$\bm{\tilde{\rho}}^{\text{el}}(\mathbf{x},\mathbf{p};t)=\left(\rho^{\text{el}},\rho^{(1),\text{el}}_{\mathbf{j}},\dots,\rho^{(N),\text{el}}_{\mathbf{j}}\right)$.
Everything can then be rewritten in the joint Liouville space of all the ADOs,

$\displaystyle\frac{\partial}{\partial
t}\bm{\tilde{\rho}}^{\text{el}}(\mathbf{x},\mathbf{p};t)$
$\displaystyle=\\{\mskip-3.0mu\mskip-3.0mu\\{H^{\text{el}}_{\text{S}}(\mathbf{x},\mathbf{p}),\bm{\tilde{\rho}}^{\text{el}}(\mathbf{x},\mathbf{p};t)\\}\mskip-3.0mu\mskip-3.0mu\\}_{a}$
$\displaystyle\>\>\>\>\>\>\>-\mathcal{L}^{\text{el}}(\mathbf{x})\bm{\tilde{\rho}}^{\text{el}}(\mathbf{x},\mathbf{p};t).$
(25)

The first term in Eq.(25) contains the symmetrized Poisson bracket of the
Wigner-transformed Hamiltonian with each Wigner-transformed ADO in
$\bm{\tilde{\rho}}^{\text{el}}(\mathbf{x},\mathbf{p};t)$,

$\displaystyle\\{\mskip-3.0mu\mskip-3.0mu\\{H^{\text{el}}_{\text{S}}(\mathbf{x},\mathbf{p}),\bm{\tilde{\rho}}^{\text{el}}(\mathbf{x},\mathbf{p};t)\\}\mskip-3.0mu\mskip-3.0mu\\}_{a}=$
$\displaystyle\>\>\>\>\>\left(\\{H^{\text{el}}_{\text{S}},\rho^{(0),\text{el}}\\}_{a},\\{H^{\text{el}}_{\text{S}},\rho^{(1),\text{el}}_{j_{1}}\\}_{a},\\{H^{\text{el}}_{\text{S}},\rho^{(1),\text{el}}_{j_{2}}.\\}_{a},\dots\right),$
(26)

which is then also unfolded into vector format. The second term,
$\mathcal{L}^{\text{el}}(\mathbf{x})$, contains all dynamics of the electronic
HEOM for a fixed vibrational frame, $\mathbf{x}$. Although not shown here, the
same procedure can be applied to the BMME in Eq.(15) Dou and Subotnik (2016).

The procedure above involves in general a classical approximation for the
vibrational dynamics. If the electron-vibration coupling is linear in
$\mathbf{x}$ and the vibrational potential is harmonic, as it is in all
transport scenarios considered in this work, then the expansion in Eq.(19)
naturally terminates after the first term and Eq.(21) is exact. For such a
system, therefore, the mixed quantum-classical approach is just a formally
exact rewriting of the original HEOM. In the next section, a further
assumption will be introduced that will make the resulting equations of motion
approximate, but for now, they remain formally exact.

###  IV.2 Fokker-Planck equation

The goal of this section is to derive a Fokker-Planck equation from Eq.(25).
The approach is very similar to Ref. Dou and Subotnik (2016), which is from
the perspective of Born-Markov theory, but requires minor changes to suit
HEOM.

First, the steady state solution,
$\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}(\mathbf{x})$, is defined as the
joint vector of all ADOs at vibrational frame $\mathbf{x}$ that satisfies
$\mathcal{L}^{\text{el}}(\mathbf{x})\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}(\mathbf{x})=0$,
which is subject to the normalization condition that
$\text{Tr}_{\text{el}}\left[\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}(\mathbf{x})\right]=1$.
Here, the electronic trace is an operation that returns the trace of only the
reduced density matrix. Its action on the joint Wigner transform of all ADOs,
for example, is

$\displaystyle\text{Tr}_{\text{el}}\left[\bm{\tilde{\rho}}^{\text{el}}(\mathbf{x},\mathbf{p};t)\right]$
$\displaystyle=\text{Tr}_{\text{el}}\left[\rho^{(0),\text{el}}(\mathbf{x},\mathbf{p};t)\right]=\text{Tr}_{\text{el}}\left[\rho^{\text{el}}(\mathbf{x},\mathbf{p};t)\right].$
(27)

Next, the phase-space quasi-probability distribution is defined as

$\displaystyle A(\mathbf{x},\mathbf{p};t)$
$\displaystyle=\text{Tr}_{\text{el}}\left[\bm{\tilde{\rho}}^{\text{el}}(\mathbf{x},\mathbf{p};t)\right]=\text{Tr}_{\text{el}}\left[\rho^{\text{el}}(\mathbf{x},\mathbf{p};t)\right],$
(28)

for which a Fokker-Planck equation will be derived.

The key step is to write
$\bm{\tilde{\rho}}^{\text{el}}(\mathbf{x},\mathbf{p};t)$ as the probability to
be at the steady state for a fixed vibrational frame, plus some difference
$\mathbf{\tilde{B}}^{\text{el}}(\mathbf{x},\mathbf{p};t)$,

$\displaystyle\bm{\tilde{\rho}}^{\text{el}}(\mathbf{x},\mathbf{p};t)$
$\displaystyle=A(\mathbf{x},\mathbf{p};t)\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}(\mathbf{x})+\mathbf{\tilde{B}}^{\text{el}}(\mathbf{x},\mathbf{p};t),$
(29)

and consider the implications if
$\mathbf{\tilde{B}}^{\text{el}}(\mathbf{x},\mathbf{p};t)$ is small.

Eq.(29) and Eq.(25) are then used to write equations of motion for
$A(\mathbf{x},\mathbf{p};t)$ and
$\mathbf{\tilde{B}}^{\text{el}}(\mathbf{x},\mathbf{p};t)$. Under the Markovian
assumption that the electronic degrees of freedom relax to equilibrium
immediately for each vibrational coordinate, $\mathbf{x}$, the equation of
motion for $\mathbf{\tilde{B}}^{\text{el}}(\mathbf{x},\mathbf{p};t)$ can be
solved and substituted into that of $A(\mathbf{x},\mathbf{p};t)$ Dou and
Subotnik (2016). This results in a Markovian Fokker-Planck equation,

$\displaystyle\frac{\partial A(\mathbf{x},\mathbf{p};t)}{\partial t}$
$\displaystyle=-\sum_{i}\frac{p_{i}}{m_{i}}\frac{\partial A}{\partial
x_{i}}+\sum_{ij}\gamma_{ij}(\mathbf{x})\frac{\partial}{\partial
p_{i}}\left(\frac{p_{j}}{m_{j}}A\right)$
$\displaystyle\>\>\>\>\>\>+\sum_{i}F_{i}(\mathbf{x})\frac{\partial A}{\partial
p_{i}}+\sum_{ij}D_{ij}(\mathbf{x})\frac{\partial^{2}A}{\partial p_{i}\partial
p_{j}},$ (30)

with details found in Appendix A. Here, the average electronic force is

$\displaystyle F_{i}(\mathbf{x})$
$\displaystyle=-\text{Tr}_{\text{el}}\left[\frac{\partial
H^{\text{el}}_{\text{S}}(\mathbf{x},\mathbf{p})}{\partial
x_{i}}\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}(\mathbf{x})\right],$ (31)

the Markovian electronic friction tensor is

$\displaystyle\gamma_{ij}(\mathbf{x})$ $\displaystyle=-\lim_{\eta\rightarrow
0^{+}}\int^{\infty}_{0}dt\>\text{Tr}_{\text{el}}\left[\frac{\partial
H^{\text{el}}_{\text{S}}}{\partial
x_{i}}e^{-(\mathcal{L}^{\text{el}}(\mathbf{x})+\eta)t}\frac{\partial\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}}{\partial
x_{j}}\right],$ (32)

and the Markovian diffusion tensor is

$\displaystyle D_{ij}(\mathbf{x})$
$\displaystyle=\frac{1}{2}\lim_{\eta\rightarrow
0^{+}}\int^{\infty}_{0}dt\>\text{Tr}_{\text{el}}\Big{[}\delta
F_{i}(\mathbf{x})e^{-(\mathcal{L}^{\text{el}}(\mathbf{x})+\eta)t}$
$\displaystyle\>\>\>\>\>\hskip 56.9055pt\left(\delta
F_{j}(\mathbf{x})\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}(\mathbf{x})+\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}(\mathbf{x})\delta
F_{j}(\mathbf{x})\right)\Big{]}.$ (33)

In Eq.(33), the operator containing the difference between the electronic
force and the average electronic force has been introduced,

$\displaystyle\delta F_{j}(\mathbf{x})$ $\displaystyle=\frac{\partial
H^{\text{el}}_{\text{S}}}{\partial
x_{i}}-\text{Tr}_{\text{el}}\left[\frac{\partial
H^{\text{el}}_{\text{S}}}{\partial
x_{i}}\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}(\mathbf{x})\right].$ (34)

Note that the Markovian assumption used to derive Eq.(30) and the
corresponding current-induced forces is completely separate to the assumption
made in deriving the BMME in Sec. III.2, in which the bath electronic degrees
of freedom relaxed quickly with respect to the system degrees of freedom.

The Markovian Langevin equation for the $i^{\text{th}}$ vibrational
coordinate, which corresponds to the Fokker-Planck equation in Eq.(30), is

$\displaystyle m_{i}\ddot{x}_{i}$
$\displaystyle=F_{i}(\mathbf{x})-\sum_{j}\gamma_{ij}(\mathbf{x})\dot{x}_{j}+f_{i}(t),$
(35)

where $f_{i}(t)$ is a Gaussian random force with white noise,

$\displaystyle\langle f_{i}(t)f_{j}(t^{\prime})\rangle$
$\displaystyle=2D_{ij}(\mathbf{x})\delta(t-t^{\prime}).$ (36)

###  IV.3 Comments on the current-induced forces

In this section, the important properties of the current-induced forces will
be discussed, as well as relevant comments on their derivation and
computation.

The form of the friction in Eq.(32) is very similar to the expressions derived
in Refs. Dou and Subotnik (2017a); Dou  _et al._ (2017); Chen  _et al._
(2019a). These approaches, however, use the total electronic density matrix
and the full time-evolution due to some general electronic system and do not
explicitly include the splitting of electronic degrees of freedom into system
and electrodes. A natural question, then, is whether the current-induced
forces in Eqs.(31)-(33) are equivalent to those in Refs. Dou and Subotnik
(2017a); Dou  _et al._ (2017); Chen  _et al._ (2019a)?

In Appendix B, this question is answered and it is shown that the current-
induced forces derived from the HEOM approach are, in fact, exactly the same
quantities as have previously been reported in Refs. Dou and Subotnik (2017a);
Dou  _et al._ (2017); Chen  _et al._ (2019a). The only condition for this is
that the vibrational degrees of freedom are restricted to the reduced system,
an assumption that already underpins the mixed quantum-classical results of
the previous section. The HEOM form of the current-induced forces, however, is
an extremely useful representation, as it now provides access to a wide
variety of transport scenarios, including systems with strong molecule-lead
couplings and strong intra-system interactions, both in- and out-of-
equilibrium.

In Refs. Dou and Subotnik (2017a); Dou  _et al._ (2017); Chen  _et al._
(2019a), furthermore, several important properties of the current-induced
forces were also derived. First, the friction and diffusion tensors are
positive definite and symmetric. Second, the fluctuation-dissipation theorem
is satisfied at equilibrium,

$\displaystyle D_{ij}(\mathbf{x})=k_{B}T\gamma_{ij}(\mathbf{x}).$ (37)

Finally, it was also shown that, under these conditions, the renormalized
vibrational force is conservative, enabling one to define an effective
potential of the mean force, $U_{\text{pmf}}(\mathbf{x})$, as

$\displaystyle U_{\text{pmf}}(\mathbf{x})$
$\displaystyle=\sum_{i}\int^{\mathbf{x}}_{\mathbf{x}_{0}}\>F_{i}(\mathbf{x}^{\prime})d\mathbf{x}^{\prime}.$
(38)

Given that the current-induced forces in Eqs.(31)-(33) are equivalent to those
in Refs. Dou and Subotnik (2017a); Dou  _et al._ (2017); Chen  _et al._
(2019a), they must also satisfy these properties.

Because a Markovian approximation has been made in deriving Eq.(35), such that
the electronic degrees of freedom are assumed to relax immediately to the
steady state for each vibrational frame, the friction and diffusion tensors
depend only on the vibrational coordinates at each time, $t$, and not on the
full history of the trajectory. If the explicit Markovian assumption is
removed, such that there was still a separation of timescales between
electronic and vibrational degrees of freedom, but
$\frac{\partial\mathbf{\tilde{B}}^{\text{el}}}{\partial t}\neq 0$, then the
derivation would yield a friction tensor with a simple form of memory; details
are given in Appendix C. To include a the full non-Markovianity of the
friction tensor, however, would require an approach such as that outlined in
Ref. Chen  _et al._ (2019a).

The friction tensor in Eq.(32) is also conceptually similar to the one derived
from Born-Markov theory in Ref. Dou and Subotnik (2016),

$\displaystyle\gamma_{ij}^{\text{BM}}(\mathbf{x})$
$\displaystyle=-\lim_{\eta\rightarrow
0^{+}}\int^{\infty}_{0}dt\>\text{Tr}_{\text{el}}\left[\frac{\partial
H^{\text{el}}_{\text{S}}}{\partial
x_{i}}e^{-(\mathcal{L}^{\text{el}}_{\text{BM}}(\mathbf{x})+\eta)t}\frac{\partial\sigma^{\text{el}}_{\text{ss}}}{\partial
x_{j}}\right].$ (39)

In Eq.(39), however, the spatial derivative acts only on the steady state
reduced system density matrix, $\sigma^{\text{el}}_{\text{ss}}(\mathbf{x})$.
Similarly, the time-evolution defined in the
$\mathcal{L}^{\text{el}}_{\text{BM}}(\mathbf{x})$ superoperator is also due to
Born-Markov theory, and not the numerically exact HEOM method.

So far the theory has been developed under the assumption that the vibrational
degrees of freedom are restricted to the reduced nanosystem. There are many
interesting systems and transport scenarios, such as current-induced bond
rupture Erpenbeck  _et al._ (2020a), where position-dependent molecular-lead
couplings are necessary for a complete description. Although such position-
dependent molecule-lead couplings can be treated within the HEOM formalism
Erpenbeck  _et al._ (2020a); Ke  _et al._ (2022), they would substantially
complicate the theory of electronic friction considered here. Eq.(24), for
example, would contain not only a Poisson bracket from the partial Wigner
transform of $H_{S}\rho^{(n)}_{\mathbf{j}}$, but also from the
$\mathcal{C}_{j_{r}}(\mathbf{x})\rho^{(n-1)}_{\mathbf{j}^{-}}$ and
$\mathcal{A}^{\bar{\sigma}}_{\alpha,m}(\mathbf{x})\rho^{(n+1)}_{\mathbf{j}^{+}}$
terms, as these superoperators now also depend on $\mathbf{x}$. Similar to
Ref.Dou and Subotnik (2017c), then, the friction would contain additional
terms arising from these extra contributions.

Finally, although the Langevin equation in Eq.(35) relates the vibrational
forces linearly to the vibrational velocities, it includes the full non-linear
dependence of the vibrational coordinates. This is different to several other
treatments of the friction Hussein  _et al._ (2010); Metelmann and Brandes
(2011); Chen  _et al._ (2019a, b); Lü  _et al._ (2012, 2010, 2011), in which
an expansion is made in the electron-vibration coupling, or, equivalently,
small vibrational displacements are assumed.

##  V Results

In order to assess the performance of the HEOM approach when calculating
current-induced forces nanosystems out-of-equilibrium, in this section, the
electronic friction of three model systems is analyzed.

In all calculations, the Fermi energy of the electrodes is set to
$\epsilon_{F}=0\text{ eV}$.

###  V.1 One electronic level coupled to one vibrational mode

Figure 1: Electronic friction as a function of the vibrational coordinate
$\tilde{x}$ and for three different values of the molecule-lead coupling
strength, $\Gamma$. The results are calculated at equilibrium, $\Phi=0\text{
V}$, and the other parameters are chosen to match Ref. Dou and Subotnik
(2018): $\varepsilon_{0}=0.15\text{ eV}$, $g=0.02\text{ eV}$,
$k_{B}T=0.026\text{ eV}$, and $\hbar\omega=0.003\text{ eV}$. Note that there
is no regime in which exists a distinguishable difference between the exact
friction and that calculated from the HEOM approach.

The first, well-studied system is a single electronic level linearly coupled
to a single vibrational coordinate. Using the dimensionless coordinate,
$\tilde{x}=x\sqrt{m\omega/\hbar}$, and momentum,
$\tilde{p}=p/\sqrt{m\omega\hbar}$, the corresponding Hamiltonian is

$\displaystyle H^{\text{el}}_{\text{S}}(\tilde{x},\tilde{p})$
$\displaystyle=\varepsilon(\tilde{x})d^{{\dagger}}d+\frac{1}{2}\omega\left(\tilde{p}^{2}+\tilde{x}^{2}\right),$
(40)

where $\varepsilon(\tilde{x})=\varepsilon_{0}+g\sqrt{2}\tilde{x}$. The
vibrational degree of freedom moves according to a harmonic potential,
$U_{0}(\tilde{x})=\frac{1}{2}\omega\tilde{x}^{2}$, when the electronic level
is unoccupied, and a shifted harmonic potential,
$U_{1}(\tilde{x})=U_{0}(\tilde{x})+\varepsilon(\tilde{x})$, when occupied.

Under the wide-band limit and because the system without the vibrational mode
is non-interacting, the Markovian electronic friction tensor under a timescale
separation assumption can be analytically derived from NEGFs, as in Refs. Dou
and Subotnik (2018) and Preston  _et al._ (2022),

$\displaystyle\gamma(\tilde{x})$
$\displaystyle=\frac{\Gamma}{2k_{B}T}\frac{m\omega}{\hbar}\left(\frac{\partial\varepsilon(\tilde{x})}{\partial\tilde{x}}\right)^{2}$
$\displaystyle\>\times\int\frac{d\omega}{4\pi}\frac{\Gamma_{L}f^{+}_{L}(\omega)f^{-}_{L}(\omega)+\Gamma_{R}f^{+}_{R}(\omega)f^{-}_{R}(\omega)}{\left(\left[\omega-\varepsilon(\tilde{x})\right]^{2}+[\Gamma/2]^{2}\right)^{2}},$
(41)

with $\Gamma=(\Gamma_{L}+\Gamma_{R})$. In Fig. 1, the analytic equilibrium
friction in Eq.(41) is compared to the friction calculated from the two
transport methods presented in Sec. III and for three different molecule-lead
couplings. For this system, the only feature of the electronic friction is a
peak at $\tilde{x}_{0}=-\varepsilon_{0}/g\sqrt{2}$. This is the vibrational
coordinate for which $\varepsilon(\tilde{x})$ crosses the chemical potential
of both electrodes and electron-hole pair creation is the dominant relaxation
process, corresponding to the well-known mechanism of electronic friction Dou
and Subotnik (2018). One can see either from Fig. 1 or Eq.(41) that the
strength of the electronic friction is inversely related to the strength of
the molecule-lead coupling. In this regime of slow vibrational motion, an
increased $\Gamma$ means that the electrons do not spend long enough on the
molecule to interact with the vibrational mode, so that both the damping rate
and rate of energy injection are reduced.

Fig. 1 also demonstrates the accuracy of the HEOM approach. For such a non-
interacting electronic system, it has previously been shown that the
corresponding hierarchical equations of motion terminates exactly at $2$nd-
tier Zheng  _et al._ (2008), such that the friction here is equivalent to the
NEGF result in Eq.(41). The Born-Markov approach, in contrast, proves to have
a limited regime of validity. When $\Gamma=10\text{ meV}$, corresponding to
the black lines in Fig. 1, all three approaches agree; as expected, Born-
Markov theory performs well when $\Gamma\ll k_{B}T$. As the molecule-lead
coupling increases to $\Gamma\approx 2k_{B}T$ and $\Gamma\approx 4k_{B}T$,
corresponding to the red and blue lines, respectively, Born-Markov theory
deviates from the exact analytic and HEOM results, producing a peak that is
narrower and taller due to its incorrect treatment of level-broadening effects
Dou and Subotnik (2018).

The HEOM approach works equally well for nonequilibrium scenarios, which is
shown in Fig. 2. At a finite bias voltage, the friction displays two peaks at
$\tilde{x}_{\alpha}=-(\varepsilon_{0}-\mu_{\alpha})/g\sqrt{2}$, where the
renormalized energy level crosses the chemical potential of each electrode.
The magnitude is also halved, as now there are only electron-hole pair
creation processes associated with one electrode. Between the two peaks,
$\gamma(\tilde{x})$ has a minimum at $\tilde{x}_{0}$, as the dominant
transport process at this point is resonant tunneling, which leads to the
well-understood mechanism of current-induced heating Erpenbeck  _et al._
(2020b); Härtle and Thoss (2011a).

Figure 2: Electronic friction as a function of the vibrational coordinate
$\tilde{x}$ and at non-zero bias voltages. Apart from the voltage, which is
applied symmetrically such that $\mu_{S}=-\mu_{D}=\Phi/2$, all parameters are
the same as in Fig. 1. The molecule-lead coupling is $\Gamma=0.1\text{ eV}$.
Note that, even for $\Phi\neq 0$, there is no regime in which exists a
distinguishable difference between the exact friction and that calculated from
the HEOM approach.

###  V.2 Anderson impurity coupled to one vibrational mode

In this section, a model describing an Anderson impurity coupled to one
vibrational mode is considered, which, due to the presence of electron-
electron interactions, cannot be solved exactly. The system Hamiltonian now
contains a spin-degenerate level linearly coupled to the vibrational
coordinate of a harmonic vibrational mode, with a non-zero, coordinate-
independent Coulomb repulsion between electrons of opposite spin,

$\displaystyle H^{\text{el}}_{\text{S}}(\tilde{x},\tilde{p})$
$\displaystyle=\varepsilon(\tilde{x})\sum_{s\in\\{\uparrow,\downarrow\\}}d^{{\dagger}}_{s}d_{s}+Ud^{{\dagger}}_{\uparrow}d_{\uparrow}d^{{\dagger}}_{\downarrow}d_{\downarrow}$
$\displaystyle\>\>\>\>\>\>+\frac{1}{2}\omega\left(\tilde{p}^{2}+\tilde{x}^{2}\right),$
(42)

with the bath annihilation and creation operators also now carrying a spin
index,

$\displaystyle H_{\text{B}}$
$\displaystyle=\sum_{\alpha\in\\{L,R\\}}\sum_{k_{\alpha},s}\varepsilon_{k_{\alpha},s}c^{{\dagger}}_{k_{\alpha},s}c_{k_{\alpha},s}$
(43) $\displaystyle H_{\text{SB}}$
$\displaystyle=\sum_{\alpha,k_{\alpha}}\sum_{s}V_{k_{\alpha},s}\left(c^{{\dagger}}_{k_{\alpha},s}d_{s}+d^{{\dagger}}_{s}c_{k_{\alpha},s}\right).$
(44)

Figure 3: Electronic friction as a function of vibrational coordinate
$\tilde{x}$ for an Anderson impurity and at a range of temperatures.
Calculations are performed at equilibrium: $\Phi=0\text{ V}$. To compare
directly with Ref. Dou  _et al._ (2017), parameters were chosen as
$\varepsilon_{0}=-0.05\text{ eV}$, $U=0.1\text{ eV}$, $g=0.0075\text{ eV}$,
and $\Gamma=0.01\text{ eV}$.

As with the single electronic level model, the Anderson impurity model has
already been investigated in the context of electronic friction Dou  _et al._
(2017); Askerka  _et al._ (2016); Dzhioev  _et al._ (2013); Plihal and
Langreth (1999, 1998). In Ref. Dou  _et al._ (2017), for example, the authors
diagonalize the electronic part of the total Hamiltonian with a numerical
renormalization group (NRG) technique and then use NEGFs to calculate
$\gamma(\tilde{x})$ at equilibrium, demonstrating the superiority of NRG to a
mean-field theory approach Askerka  _et al._ (2016); Dzhioev  _et al._ (2013),
which does not correctly recover the doubly-peaked behavior of
$\gamma(\tilde{x})$. The HEOM method developed in this paper, in comparison,
is able to perform a numerically exact calculation of the friction for this
model out-of-equilibrium.

To compare with the NEGF approach, Fig. 3 shows, for a range of temperatures,
the electronic friction as a function of the vibrational coordinate for the
same parameters as in Fig. 1a from Ref. Dou  _et al._ (2017). At equilibrium,
the friction displays two peaks around the points where there is an electron-
hole pair creation resonance, $\varepsilon(\tilde{x})=0$ and
$\varepsilon(\tilde{x})+U=0$. At high to moderate temperatures, the HEOM
approach yields the same results as NEGFs. At low temperatures, however, so
many poles are required for the Padé decomposition of the Fermi-Dirac function
that the HEOM become numerically intractable, a well-known problem that limits
the temperature range of the method. Recent proposals Chen  _et al._ (2022);
Xu  _et al._ (2022); Nakatsukasa  _et al._ (2018) for a more efficient pole
expansion appear promising in extending the HEOM approach to low temperatures.

Figure 4: Electronic friction as a function of vibrational coordinate
$\tilde{x}$ for an Anderson impurity model and at a range of bias voltages.
Calculations are performed for $k_{B}T=\Gamma=0.01\text{ eV}$, and all other
parameters are the same as in Fig. 3.

At finite bias voltage, which is shown in Fig. 4, the two friction peaks split
into four, as there is now a resonance whenever an electron-hole pair creation
process is resonant with the chemical potential of an electrode:
$\varepsilon(\tilde{x})=\mu_{\alpha}$ and
$\varepsilon(\tilde{x})+U=\mu_{\alpha}$. The peaks at
$\varepsilon(\tilde{x})=\mu_{\alpha}$ are larger than those at
$\varepsilon(\tilde{x})+U=\mu_{\alpha}$, because in the first case the
degenerate level reaches $\mu_{\alpha}$ and two electron-hole pair creation
processes become available, whereas in the second it is only the doubly-
occupied process that activates. At $\Phi=0.1\text{ V}$, there is only one
extra peak at $\tilde{x}=0$, because $\varepsilon(0)+U=\mu_{S}$ and
$\varepsilon(0)=\mu_{D}$.

###  V.3 One electronic level coupled to one classical and one quantum
vibrational mode

Figure 5: Electronic friction in equilibrium as a function of vibrational
coordinate $\tilde{x}$ for a single electronic level coupled to a high-
frequency quantum vibrational mode and low-frequency classical vibrational
mode. Results are shown for two different molecule-lead couplings,
$\Gamma=0.01\text{ eV}$ (red) and $\Gamma=0.026\text{ eV}$ (blue), and from
the HEOM and Born-Markov approaches, corresponding to the solid and dashed
lines, respectively. Other parameters are $g_{\text{cl}}=0.02\text{ eV}$,
$k_{B}T=0.026\text{ eV}$, $\hbar\omega_{\text{qm}}=0.2\text{ eV}$, and in (a)
$\lambda_{\text{qm}}=0.1\text{ eV}$, while in (b)
$\lambda_{\text{qm}}=0.4\text{ eV}$. The energy, $\varepsilon_{0}$, is chosen
so that the renormalized energy level after the small polaron transformation
is always
$\varepsilon_{0}-\frac{\lambda_{\text{qm}}^{2}}{\omega_{\text{qm}}}=0.3\text{
eV}$.

Figure 6: Electronic friction at finite bias voltages and as a function of
vibrational coordinate $\tilde{x}$ for an electronic level coupled to a high-
frequency quantum vibrational mode as well as a low-frequency classical
vibrational mode. The molecule-lead coupling is $\Gamma=0.026\text{ eV}$. The
upper row corresponds to $\lambda_{\text{qm}}=0.1\text{ eV}$, while the lower
row corresponds to $\lambda_{\text{qm}}=0.4\text{ eV}$. The left column is
held at a voltage of $\Phi=0.5\text{ V}$, while the right column is held at a
voltage of $\Phi=1.0\text{ V}$. All other parameters are the same as in Fig.
5.

In the final example, a single electronic level is coupled to two vibrational
modes. One, which is treated classically, is subject to the limit of slow
vibrational driving, with a low frequency $\hbar\omega_{\text{cl}}\ll\Gamma$,
and one, which remains in the fully quantum mechanical HEOM treatment, has a
high frequency $\hbar\omega_{\text{qm}}>\Gamma$. Although this degree of
freedom is vibrational and not electronic, it can still be included in
$H^{\text{el}}_{\text{S}}(x_{\text{cl}},p_{\text{cl}})$ as long as it follows
the timescale separation necessary for the adiabatic approximation. As
discussed in Ref. Chen  _et al._ (2019b), in which the electronic friction for
a similar model was calculated under different assumptions, such a model could
describe a single-molecule in an optical cavity, where the high-frequency mode
describes a cavity mode with a strong light-matter interaction Kongsuwan  _et
al._ (2018); Chikkaraddy  _et al._ (2016). An alternative experimental
realization would be a molecule with two bonds of distinctly different
strengths.

Similar to the previous two case studies, the electronic level is linearly
coupled to both the classical and quantum vibrational coordinates, such that
the resulting Hamiltonian is

$\displaystyle
H^{\text{el}}_{\text{S}}(\tilde{x}_{\text{cl}},\tilde{p}_{\text{cl}})$
$\displaystyle=\varepsilon(\tilde{x}_{\text{cl}})d^{{\dagger}}d+\omega_{\text{qm}}b^{{\dagger}}b+\lambda_{\text{qm}}(b^{{\dagger}}+b)d^{{\dagger}}d$
$\displaystyle\>\>\>\>\>\>+\frac{1}{2}\omega_{\text{cl}}\left(\tilde{p}_{\text{cl}}^{2}+\tilde{x}_{\text{cl}}^{2}\right).$
(45)

The bosonic annihilation and creation operators are related to the
dimensionless coordinate,
$\hat{\tilde{x}}_{\text{qm}}=\hat{x}_{\text{qm}}\sqrt{m_{\text{qm}}\omega_{\text{qm}}/\hbar}$,
and momentum,
$\hat{\tilde{p}}_{\text{qm}}=\hat{p}_{\text{qm}}/\sqrt{m_{\text{qm}}\omega_{\text{qm}}\hbar}$,
of the quantum vibrational mode via

$\displaystyle b$
$\displaystyle=\frac{1}{\sqrt{2}}(\hat{\tilde{x}}_{\text{qm}}+i\hat{\tilde{p}}_{\text{qm}})\hskip
7.11317pt\text{ and }\hskip
7.11317ptb^{{\dagger}}=\frac{1}{\sqrt{2}}(\hat{\tilde{x}}_{\text{qm}}-i\hat{\tilde{p}}_{\text{qm}}).$
(46)

Unlike the previous two case studies, however, it is difficult to investigate
this model with other transport methods, as it contains a strong electron-
phonon interaction. This model serves, therefore, as an excellent example of
the new variety of nanosystems now accessible by semi-classical Langevin
equations. The quantum part of the Hamiltonian can be diagonalized via a small
polaron transformation Lang and Firsov (1963), which yields

$\displaystyle
H^{\text{el}}_{\text{S}}(\tilde{x}_{\text{cl}},\tilde{p}_{\text{cl}})$
$\displaystyle=\tilde{\varepsilon}(\tilde{x}_{\text{cl}})d^{{\dagger}}d+\omega_{\text{qm}}b^{{\dagger}}b+\frac{1}{2}\omega_{\text{cl}}\left(\tilde{p}_{\text{cl}}^{2}+\tilde{x}_{\text{cl}}^{2}\right),$
(47)

a Hamiltonian similar to Eq.(40), except that the electronic energy level has
been renormalized to

$\displaystyle\tilde{\varepsilon}(\tilde{x}_{\text{cl}})$
$\displaystyle=\varepsilon_{0}-\frac{\lambda_{\text{qm}}^{2}}{\omega_{\text{qm}}}+\sqrt{2}g_{\text{cl}}\tilde{x}_{\text{cl}}.$
(48)

The quantum transport under the small polaron transformation occurs between
eigenstates of the system Hamiltonian
$|q,\nu\rangle\rightarrow|q^{\prime},\nu^{\prime}\rangle$, where the quantum
numbers $q=0,1$ and $\nu=0,1,2,\dots$ refer to charge and phonon occupancy,
respectively. Transporting electrons can exchange quanta of vibrational energy
with the quantum harmonic oscillator, with the transition probability between
vibrational states of the occupied and unoccupied system given by the overlap
of the wave functions of the original and displaced quantum harmonic
oscillator, which is the Franck-Condon matrix,

$\displaystyle X_{\nu,\nu^{\prime}}$
$\displaystyle=\langle\nu|X|\nu^{\prime}\rangle=\langle\nu|e^{\frac{\lambda_{\text{qm}}}{\omega_{\text{qm}}}(b-b^{{\dagger}})}|\nu^{\prime}\rangle.$
(49)

In this transport picture, the ratio $\lambda_{\text{qm}}/\omega_{\text{qm}}$
determines the structure of the Franck-Condon matrix. For small
$\lambda_{\text{qm}}$, the $X_{\nu,\nu^{\prime}}$ are close to diagonal and
the transitions between low-lying vibrational states are favored, while the
opposite is true for strong $\lambda_{\text{qm}}$ Koch and von Oppen (2005);
Koch  _et al._ (2006a); Härtle and Thoss (2011a, b).

Due to the strong interaction between the electronic level and the quantum
vibrational mode, Eq.(47) is a model well-suited to the HEOM approach. Under
the appropriate high-temperature, weak-coupling limit, one could also use a
Born-Markov master equation, as this formulation is also able to treat strong
intra-system interactions exactly. To test this, Fig. 5(a) and Fig. 5(b) show
the Born-Markov friction against the HEOM friction at equilibrium and for two
different molecule-lead couplings, $\Gamma$, and electron-phonon couplings,
$\lambda_{\text{qm}}$. As expected, in the strong-coupling regime of
$\Gamma=k_{B}T=0.026\text{ eV}$, Born-Markov theory does not reproduce the
numerically exact HEOM result. Even in the weak-coupling regime of
$\Gamma=0.01\text{ eV}<k_{B}T$, furthermore, there are significant differences
between $\gamma$ calculated from the HEOM and BMME approaches, indicating that
the electronic friction for this model requires the superior HEOM approach
even for moderate temperatures and molecule-lead couplings.

Apart from the differences between the Born-Markov and HEOM approaches, the
qualitative friction behavior strongly resembles that of the single electronic
level in equilibrium, in that there is a single peak at

$\displaystyle\tilde{x}_{0}$
$\displaystyle=-\frac{\varepsilon_{0}-\lambda_{\text{qm}}^{2}/\omega_{\text{qm}}}{g_{\text{cl}}\sqrt{2}}.$
(50)

In Fig. 5(b), however, the height of the friction peak is an order of
magnitude larger than it is in the non-interacting case, which originates from
the additional electron-hole pair creation processes introduced by the quantum
vibrational mode, as was discussed in Ref. Härtle and Thoss (2011a). At
equilibrium, an electron-hole pair can be created with an effective cooling
process when an electron tunnels in from electrode $\alpha$ and absorbs one or
more vibrational quanta, allowing it to then tunnel out to an unoccupied state
in the same electrode; see Fig.(2) in Ref. Härtle and Thoss (2011a) for a
schematic example. Because the transition probability from the vibrational
ground to excited states,
$|X_{0,\nu}|^{2}=\frac{1}{\nu!}(\frac{\lambda_{\text{qm}}}{\omega_{\text{qm}}})^{2\nu}e^{-(\lambda_{\text{qm}}/\omega_{\text{qm}})^{2}}$,
is not suppressed for $\lambda_{\text{qm}}/\omega_{\text{qm}}=2$, these
electron-hole pair creation processes contribute significantly to the
friction. For the weaker electron-phonon coupling in Fig. 5(a), where
$\lambda_{\text{qm}}/\omega_{\text{qm}}=0.5$, the transition probabilities are
suppressed and the friction is also quantitatively similar to the non-
interacting system in Fig. 1.

At finite bias voltages, the electronic friction displays a much richer
structure, which is shown in Fig. 6. In Fig. 6(c) and Fig. 6(d), where
$\lambda_{\text{qm}}/\omega_{\text{qm}}=2$, the friction displays
characteristic peaks at points where

$\displaystyle\tilde{x}_{\alpha,\nu}$
$\displaystyle=\frac{\mu_{\alpha}-(\varepsilon_{0}-\lambda_{\text{qm}}^{2}/\omega_{\text{qm}}+q\omega_{\text{qm}})}{g_{\text{cl}}\sqrt{2}}:$
(51)

that is, points where the renormalized energy level plus some integer of
vibrational quanta crosses the chemical potential of electrode $\alpha$. At
each $\tilde{x}_{\alpha,\nu}$, the electron-hole pair creation process in
electrode $\alpha$ associated with the absorption of $\nu$ vibrational quanta
is enhanced. These side peaks are not visible at equilibrum, however, as only
the ground to ground state transition is energetically allowed Koch and von
Oppen (2005).

In contrast, the friction minima occur when the chemical potential of
electrode $\alpha$ is located between two such quasi-levels and the system
experiences current-induced heating. The peaks also decrease as the
renormalized level is pushed further into the center of the bias window; the
further
$\varepsilon_{0}-\lambda_{\text{qm}}^{2}/\omega_{\text{qm}}+g_{\text{cl}}\sqrt{2}\tilde{x}$
is away from $\mu_{\alpha}$, the more vibrational quanta are required for
electron-hole pair creation, and these processes are generally less probable.
Finally, the main difference between Figs. 6(c) and Fig. 6(d) is that the
higher voltage in Fig. 6(d) energetically allows more vibrational transitions
to contribute to electron-hole pair creation, causing the extra two friction
peaks.

In Fig. 6(a) and Fig. 6(b), the electron-phonon coupling is again weak,
$\lambda_{\text{qm}}/\omega_{\text{qm}}=0.5$, and electron-hole pair creation
via absorption of vibrational quanta is suppressed with increasing voltage. At
$\tilde{x}_{\alpha,0}$, the resonant electron-hole pair creation process still
yields a friction peak because the vibrational ground to excited state
transition has a strong Franck-Condon factor. At each of the other
$\tilde{x}_{\alpha,\nu\neq 0}$, however, the electronic friction experiences a
local or global minimum, as the $|X_{0,\nu>0}|^{2}$ are suppressed for
increasing $\nu$ and $\lambda_{\text{qm}}/\omega_{\text{qm}}<1$. In this
regime, the electronic friction actually becomes negative. This implies that
the loss of cooling by electron-hole pair creation not only heats the quantum
vibrational mode, but that energy is also being pumped into the classical mode
by the usually dissipative processes, driving it to instability. The
identification of negative friction is of high interest in transport through
nanostructures, as it provides an easy insight into regimes of instability
where the junction can break Preston  _et al._ (2022); Lü  _et al._ (2012).

##  VI Conclusion

In this work, the HEOM approach to nonequilibrium charge transport through
nanosystems was used to develop a new method for calculating current-induced
forces. Similar to other master equation approaches to electronic friction Dou
and Subotnik (2016), the fully quantum HEOM is first transformed to a QCLE and
then, under an adiabatic approximation for the vibrational degrees of freedom,
to a Markovian Fokker-Planck equation, from which the average electronic
force, Markovian friction tensor, Markovian diffusion tensor, and
corresponding Langevin equation are identified. These quantities were then
shown to be equivalent to previous derivations of the current-induced forces
under a timescale separation, but in a more computationally useful form. The
electronic friction, for example, can now be calculated for junctions with
both strong molecule-lead couplings and strong intra-system interactions. The
resulting friction and diffusion tensors, furthermore, are completely general
in- and out-of-equilibrium.

The HEOM approach to calculating electronic friction was applied to three
models of different complexity. First, a single non-interacting electronic
level coupled linearly to one classical vibrational mode was considered, a
system for which exact results can be calculated; these were used to
demonstrate the numerical accuracy of the HEOM approach. In the second model,
an Anderson impurity coupled to one classical vibrational mode was considered.
Because adding extra levels and interactions does not complicate the HEOM
theory, this system can be treated easily both in- and out-of-equilibrium.
Finally, a system with an additional interaction between the electronic level
and a quantum vibrational mode was considered, which yielded a richer
structure in the friction of the classical mode. In particular, it was shown
that for a simple one level, two mode system, heating of the quantum
vibrational mode caused negative friction in the classical mode, indicating
that it is also being driven to instability.

##  VII Acknowledgements

This work was supported by the German Research Foundation (DFG) through
FOR5099. S.L.R thanks the Alexander von Humboldt Foundation for the award of a
Research Fellowship. Furthermore, support by the state of Baden-Württemberg
through bwHPC and the DFG through Grant No. INST 40/575-1 FUGG (JUSTUS 2
cluster) is gratefully acknowledged.

##  Appendix A On the derivation of the Markovian Fokker-Planck equation

In this appendix, the Markovian Fokker-Planck equation in Eq.(30) will be
explicitly derived. First, the equations of motion for
$A(\mathbf{x},\mathbf{p};t)$ and
$\mathbf{\tilde{B}}^{\text{el}}(\mathbf{x},\mathbf{p};t)$ are

$\displaystyle\frac{\partial A}{\partial t}$
$\displaystyle=\text{Tr}_{\text{el}}[\\{\mskip-3.0mu\mskip-3.0mu\\{H^{\text{el}}_{\text{S}},A\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\\}\mskip-3.0mu\mskip-3.0mu\\}_{a}]+\text{Tr}_{\text{el}}[\\{\mskip-3.0mu\mskip-3.0mu\\{H^{\text{el}}_{\text{S}},\mathbf{\tilde{B}}^{\text{el}}\\}\mskip-3.0mu\mskip-3.0mu\\}_{a}]$
$\displaystyle=-\sum_{i}\frac{p_{i}}{m_{i}}\frac{\partial A}{\partial
x_{i}}+\sum_{i}\text{Tr}_{\text{el}}\left[\frac{\partial
H^{\text{el}}_{\text{S}}}{\partial
x_{i}}\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\right]\frac{\partial
A}{\partial p_{i}}$
$\displaystyle\>\>\>\>\>\>+\sum_{i}\text{Tr}_{\text{el}}\left[\frac{\partial
H^{\text{el}}_{\text{S}}}{\partial
x_{i}}\frac{\partial\mathbf{\tilde{B}}^{\text{el}}}{\partial
p_{i}}\right]\hskip 28.45274pt\text{ and }$ (52)
$\displaystyle\frac{\partial\mathbf{\tilde{B}}^{\text{el}}}{\partial t}$
$\displaystyle=\\{\mskip-3.0mu\mskip-3.0mu\\{H^{\text{el}}_{\text{S}},\mathbf{\tilde{B}}^{\text{el}}\\}\mskip-3.0mu\mskip-3.0mu\\}_{a}-\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\text{Tr}_{\text{el}}[\\{\mskip-3.0mu\mskip-3.0mu\\{H^{\text{el}}_{\text{S}},\mathbf{\tilde{B}}^{\text{el}}\\}\mskip-3.0mu\mskip-3.0mu\\}_{a}]-\mathcal{L}^{\text{el}}\mathbf{\tilde{B}}^{\text{el}}$
$\displaystyle\>\>\>\>\>\>+\\{\mskip-3.0mu\mskip-3.0mu\\{H^{\text{el}}_{\text{S}},A\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\\}\mskip-3.0mu\mskip-3.0mu\\}_{a}-\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\text{Tr}_{\text{el}}[\\{\mskip-3.0mu\mskip-3.0mu\\{H^{\text{el}}_{\text{S}},A\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\\}\mskip-3.0mu\mskip-3.0mu\\}_{a}].$
(53)

where the $(\mathbf{x},\mathbf{p};t)$ notation has been suppressed for brevity
and
$\text{Tr}_{\text{el}}\left[\mathcal{L}^{\text{el}}\bm{\tilde{\rho}}^{\text{el}}\right]=0$
has been used. This relation can be seen if one considers that
$\text{Tr}_{\text{el}}\left[\dots\right]$ only traces over the $0$th-tier ADO,
and the corresponding part of the hierarchy in $\mathcal{L}^{\text{el}}$
contains only commutators,

$\displaystyle\text{Tr}_{\text{el}}\left[\mathcal{L}^{\text{el}}\bm{\tilde{\rho}}^{\text{el}}\right]$
$\displaystyle=-i\text{Tr}_{\text{el}}\left(\left[H^{\text{el}}_{\text{S}},\rho^{\text{el}}\right]\right)$
$\displaystyle\>\>\>\>-i\sum_{j}V_{\alpha,m}\text{Tr}_{\text{el}}\left(\left[d^{\bar{\sigma}},\bm{\tilde{\rho}}^{(1),\text{el}}_{j}\right]\right)$
(54) $\displaystyle=0.$ (55)

The first few terms in Eq.(53) involve derivatives of
$\mathbf{\tilde{B}}^{\text{el}}$ with respect to $t$, $x_{i}$, and $p_{i}$,
due to the Poisson bracket. As has been discussed in Refs. Dou  _et al._
(2017, 2015), the limit of slow vibrational degrees of freedom in comparison
to fast electronic degrees of freedom allows these terms to be neglected, as
$\mathbf{\tilde{B}}^{\text{el}}$ should be much smaller than $A$. Eq.(53)
becomes, therefore,

$\displaystyle\mathcal{L}^{\text{el}}\mathbf{\tilde{B}}^{\text{el}}$
$\displaystyle=\\{\mskip-3.0mu\mskip-3.0mu\\{H^{\text{el}}_{\text{S}},A\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\\}\mskip-3.0mu\mskip-3.0mu\\}_{a}-\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\text{Tr}_{\text{el}}[\\{\mskip-3.0mu\mskip-3.0mu\\{H^{\text{el}}_{\text{S}},A\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\\}\mskip-3.0mu\mskip-3.0mu\\}_{a}]$
(56)
$\displaystyle=-\sum_{j}\frac{p_{j}}{m_{j}}\frac{\partial\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}}{\partial
x_{j}}A-\sum_{j}\text{Tr}_{\text{el}}\left[\frac{\partial
H^{\text{el}}_{\text{S}}}{\partial
x_{j}}\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\right]\frac{\partial
A}{\partial p_{j}}\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}$
$\displaystyle\>\>\>\>\>\>+\frac{1}{2}\sum_{j}\left(\frac{\partial
H^{\text{el}}_{\text{S}}}{\partial
x_{j}}\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}+\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\frac{\partial
H^{\text{el}}_{\text{S}}}{\partial x_{j}}\right)\frac{\partial A}{\partial
p_{j}}.$ (57)

With the $(\mathbf{x},\mathbf{p};t)$ notation now explicitly written, and
using the identity

$\displaystyle[\mathcal{L}^{\text{el}}(\mathbf{x})]^{-1}$
$\displaystyle=\lim_{\eta\rightarrow
0^{+}}\int^{\infty}_{0}dt\>e^{-(\mathcal{L}^{\text{el}}(\mathbf{x})+\eta)t},$
(58)

$\mathcal{L}^{\text{el}}(\mathbf{x})$ is formally inverted and
$\mathbf{\tilde{B}}^{\text{el}}(\mathbf{x},\mathbf{p};t)$ can be directly
calculated. This is then substituted into the equation of motion for
$A(\mathbf{x},\mathbf{p};t)$, yielding the Markovian Fokker-Planck equation in
Eq.(30).

##  Appendix B Comparison to previous approaches

The goal of this section is to show that the current-induced forces derived
via the HEOM approach in the main text are equivalent to other approaches that
treat all electronic degrees of freedom without an explicit system-bath
partitioning. See, for example, Eq.(11), Eq.(18), and Eq.(19) in Ref. Dou and
Subotnik (2018) or the Markovian version of Eq.(27) in Ref. Chen  _et al._
(2019a). In these approaches, the average force is

$\displaystyle F_{i}(\mathbf{x})$
$\displaystyle=-\text{Tr}_{\text{tot},\text{el}}\left[\frac{\partial
H^{\text{el}}}{\partial x_{i}}\rho^{\text{el}}_{\text{tot},\text{ss}}\right],$
(59)

the Markovian electronic friction is

$\displaystyle\gamma_{ij}(\mathbf{x})$
$\displaystyle=-\int^{\infty}_{0}dt\>\text{Tr}_{\text{tot},\text{el}}\left[\frac{\partial
H^{\text{el}}}{\partial
x_{i}}e^{-iH^{\text{el}}t}\frac{\partial\rho^{\text{el}}_{\text{tot},\text{ss}}}{\partial
x_{j}}e^{iH^{\text{el}}t}\right],$ (60)

and the Markovian diffusion tensor is

$\displaystyle D_{ij}(\mathbf{x})$
$\displaystyle=\int^{\infty}_{0}dt\>\text{Tr}_{\text{tot},\text{el}}\left[\delta
F_{i}e^{-iH^{\text{el}}t}\left(\delta
F_{j}\rho^{\text{el}}_{\text{tot},\text{ss}}\right.\right.$
$\displaystyle\hskip
76.82234pt\left.\left.+\rho^{\text{el}}_{\text{tot},\text{ss}}\delta
F_{j}\right)e^{iH^{\text{el}}t}\right],$ (61)

with

$\displaystyle\delta F_{i}$ $\displaystyle=\frac{\partial
H^{\text{el}}}{\partial x_{i}}-F_{i}(\mathbf{x}).$ (62)

Here, $\text{Tr}_{\text{tot},\text{el}}(\dots)$ and
$\rho^{\text{el}}_{\text{tot},\text{ss}}$ are the trace over and the steady
state density matrix of all electronic degrees of freedom in the system and
bath, respectively. This is in contrast to Eqs.(31)-(33) in the main text,
where the trace is only over the system electronic degrees of freedom and
$\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}$ contains the steady state system
electronic density matrix as well as all steady state ADOs.

The first question, then, is whether the information one extracts from
incorporating the bath effects via the ADOs is the same as that from treating
the full electronic density matrix of system and bath. This is, however, not
the only difference between the approaches. Eq.(60) and Eq.(61) contain the
time-evolution of all electronic degrees of freedom according to the full
Hamiltonian, $H^{\text{el}}$, whereas Eq.(60) and Eq.(61) contain the time-
evolution according to a HEOM derived under the assumption that at the initial
state of the object being propagated, $A_{T}(0)$, can be factorized into
system and bath components, $A_{T}(0)=A\otimes A_{B}$. Since the quantities
$\frac{\partial\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}}{\partial x_{j}}$
and $\delta F_{j}(\mathbf{x})\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}$
evidently cannot be factorized, this raises another question; is the HEOM
propagation defined by $\mathcal{L}^{\text{el}}(\mathbf{x})$ correct for these
quantities? In the remainder of this section, both of these questions will be
answered and the equivalence of the two forms of the current-induced forces
will be shown.

###  B.1 Average electronic force

As in the main text, the key assumption is that all vibrational degrees of
freedom are restricted to the system only, such that

$\displaystyle\frac{\partial H^{\text{el}}(\mathbf{x})}{\partial x_{i}}$
$\displaystyle=\frac{\partial H^{\text{el}}_{\text{S}}(\mathbf{x})}{\partial
x_{i}}$ (63)

To start, consider the average electronic force, $F_{i}(\mathbf{x})$. Because
the spatial derivative now only involves system operators, it can be
simplified to

$\displaystyle F_{i}(\mathbf{x})$
$\displaystyle=-\text{Tr}_{\text{el}}\left[\frac{\partial
H^{\text{el}}_{S}}{\partial
x_{i}}\text{Tr}_{B}\left\\{\rho^{\text{el}}_{\text{tot},\text{ss}}\right\\}\right]=-\text{Tr}_{\text{el}}\left[\frac{\partial
H^{\text{el}}_{S}}{\partial x_{i}}\sigma^{\text{el}}_{\text{ss}}\right],$ (64)

which is exactly the result one obtains from the HEOM approach in Eq.(31). The
system electronic steady state density matrix obtained by HEOM is, by
definition, the same as that obtained by the time-evolution of all electronic
degrees of freedom to the steady state and then tracing out the bath, because
HEOM is numerically exact and all systems considered have a unique steady
state.

###  B.2 Friction tensor

Figure 7: Schematic of the transformation of the forward paths into one joint
path. In the left plot, the paths $(\bm{\xi_{a}},\bm{\xi^{\prime}_{a}})$ and
$(\bm{\xi_{b}},\bm{\xi^{\prime}_{b}})$ are joined at $t=0$ by the midpoint
$(\bm{\xi_{m}},\bm{\xi^{\prime}_{m}})$. In Eq.(73), however, this midpoint is
integrated over and the two paths combine into one, albeit with different
actions before and after $t=0$, which is shown in the right plot.

Next, consider the more difficult quantities defined in Eq.(60) and Eq.(61).
In what follows, the equivalence will be shown explicitly for the friction, as
the diffusion follows a very similar process. Similar to the previous
subsection, Eq.(60) can be rewritten as

$\displaystyle\gamma_{ij}(\mathbf{x})$
$\displaystyle=-\int^{\infty}_{0}dt\>\text{Tr}_{\text{el}}\left[\frac{\partial
H^{\text{el}}_{S}}{\partial x_{i}}A(\mathbf{x},t)\right],$ (65)

where

$\displaystyle A(\mathbf{x},t)$
$\displaystyle=\text{Tr}_{B}\left[e^{-iH^{\text{el}}(\mathbf{x})t}\frac{\partial\rho^{\text{el}}_{\text{tot},\text{ss}}}{\partial
x_{j}}e^{iH^{\text{el}}(\mathbf{x})t}\right],$ (66)

has been defined; it is evidently this quantity that now must be shown to be
equivalent to the HEOM approach. Note that the time-evolution has been
explicitly written only as a function of the vibrational coordinates,
$\mathbf{x}$, because the classical vibrational momenta do not couple directly
to the electronic degrees of freedom.

The next step, then, is to introduce the formal definition of the spatial
derivative,

$\displaystyle\frac{\partial\rho^{\text{el}}_{\text{tot},\text{ss}}}{\partial
x_{j}}$ $\displaystyle=\lim_{dx_{j}\rightarrow
0^{+}}\frac{\rho^{\text{el}}_{\text{tot},\text{ss}}(\mathbf{x}+dx_{j})-\rho^{\text{el}}_{\text{tot},\text{ss}}(\mathbf{x)}}{dx_{j}}.$
(67)

Because the electronic systems under consideration have unique steady states,
$\rho^{\text{el}}_{\text{tot},\text{ss}}(\mathbf{x}+dx_{j})$ and
$\rho^{\text{el}}_{\text{tot},\text{ss}}(\mathbf{x})$ can be found by starting
the total system at some point in the distant past, where they are assumed to
be uncorrelated,
$\lim\limits_{\tau\rightarrow-\infty}\rho^{\text{el}}_{\text{tot}}(\tau)=\rho^{\text{el}}(\tau)\rho^{\text{el}}_{B}(\tau)$,
and evolving them up until time $t=0$, such that

$\displaystyle\rho^{\text{el}}_{\text{tot},\text{ss}}(\mathbf{x}+dx_{j})$
$\displaystyle=\lim_{\tau\rightarrow-\infty}e^{iH^{\text{el}}(\mathbf{x}+dx_{j})\tau}\rho^{\text{el}}(\tau)\rho^{\text{el}}_{B}(\tau)e^{-iH^{\text{el}}(\mathbf{x}+dx_{j})\tau}$
(68) $\displaystyle\rho^{\text{el}}_{\text{tot},\text{ss}}(\mathbf{x})$
$\displaystyle=\lim_{\tau\rightarrow-\infty}e^{iH(x)\tau}\rho^{\text{el}}(\tau)\rho^{\text{el}}_{B}(\tau)e^{-iH(x)\tau}.$
(69)

Eq.(66), therefore, will contain two terms, each with two sets of time-
evolution. Consider the first term,

$\displaystyle A_{\mathbf{x}+dx_{j}}(\mathbf{x},t)$
$\displaystyle=\text{Tr}_{B}\left[e^{-iH^{\text{el}}(\mathbf{x})t}\rho^{\text{el}}_{\text{tot},\text{ss}}(\mathbf{x}+dx_{j})e^{iH^{\text{el}}(\mathbf{x})t}\right]$
(70)
$\displaystyle=\lim_{\tau\rightarrow-\infty}\text{Tr}_{B}\left[e^{-iH^{\text{el}}(\mathbf{x})t}e^{iH^{\text{el}}(\mathbf{x}+dx_{j})\tau}\right.$
$\displaystyle\hskip
42.67912pt\left.\rho^{\text{el}}(\tau)\rho^{\text{el}}_{B}(\tau)e^{-iH^{\text{el}}(\mathbf{x}+dx_{j})\tau}e^{iH^{\text{el}}(\mathbf{x})t}\right],$
(71)

which can now be evaluated using HEOM. Note that the subscript notation
containing $\mathbf{x}+dx_{j}$ specifies the vibrational coordinates as the
electronic degrees of freedom evolve to the steady state from $\tau$ to $t=0$.
It is only this part that needs to be specified, because, from Eq.(65), one
can see that the evolution after $t=0$ for any quantity is always at
$\mathbf{x}$.

Not all steps of the standard HEOM derivation will be reproduced here, as they
can be found in Refs. Tanimura and Kubo (1989); Tanimura (2006); Härtle  _et
al._ (2015); Jin  _et al._ (2007, 2008); Zheng  _et al._ (2009); Yan (2014);
Wenderoth  _et al._ (2016); Tanimura (2020); Ye  _et al._ (2016). First, the
basis of fermionic coherent states, which are eigenvectors of the fermionic
annihilation and creation operators, is introduced

$\displaystyle d_{i}|\xi\rangle$ $\displaystyle=\xi_{i}|\xi\rangle\hskip
14.22636pt\text{ and }\hskip
14.22636pt\langle\xi|d^{{\dagger}}_{i}=\langle\xi|\xi_{i}^{*}.$ (72)

The eigenvalues $\xi_{m}$, which are often collected into vectors,
$\bm{\xi}=\left(\xi_{1},\xi_{2},\dots\right)$ and
$\bm{\xi^{*}}=\left(\xi^{*}_{1},\xi^{*}_{2},\dots\right)$, are not ordinary
complex numbers, but rather Grassmann variables satisfying the anti-
commutation relation between fermionic annihilation and creation operators. In
this basis, $A_{\mathbf{x}+dx_{j}}(\mathbf{x},t)$ is

$\displaystyle
A_{\mathbf{x}+dx_{j}}(\bm{\xi_{f}},\bm{\xi^{\prime}_{f}};\mathbf{x},t)$
$\displaystyle=\lim_{\tau\rightarrow-\infty}\int(d\bm{\xi_{i}^{*}}d\bm{\xi_{i}})\>(d\bm{\xi_{i}^{*\prime}}d\bm{\xi_{i}^{\prime}})\>(d\bm{\xi_{i_{B}}^{*}}d\bm{\xi_{i_{B}}})\>(d\bm{\xi_{i_{B}}^{*\prime}}d\bm{\xi^{\prime}_{i_{B}}})e^{-\bm{\xi_{i}^{*}}\bm{\xi_{i}}}e^{-\bm{\xi_{i}^{*\prime}}\bm{\xi^{\prime}_{i}}}e^{-\bm{\xi_{i_{B}}^{*}}\bm{\xi_{i_{B}}}}e^{-\bm{\xi_{i_{B}}^{*\prime}}\bm{\xi^{\prime}_{i_{B}}}}$
$\displaystyle\hskip
36.98866pt\int(d\bm{\xi_{m}^{*}}d\bm{\xi_{m}})\>(d\bm{\xi_{m}^{*\prime}}d\bm{\xi^{\prime}_{m}})(d\bm{\xi_{m_{B}}^{*}}d\bm{\xi_{m_{B}}})\>(d\bm{{\xi_{m_{B}}^{*\prime}}}d\bm{\xi^{\prime}_{m_{B}}})e^{-\bm{\xi_{m}^{*}}\bm{\xi_{m}}}e^{-\bm{\xi_{m}^{*\prime}}\bm{\xi^{\prime}_{m}}}e^{-\bm{\xi_{m_{B}}^{*}}\bm{\xi_{m_{B}}}}e^{-\bm{\xi_{m_{B}}^{*\prime}}\bm{\xi^{\prime}_{m_{B}}}}$
$\displaystyle\hskip
36.98866pt\text{Tr}_{B}\Big{[}G_{\mathbf{x}}(\bm{\xi_{f_{B}}},\bm{\xi_{f}},t;\bm{\xi_{m_{B}}},\bm{\xi_{m}},0)G_{\mathbf{x}+dx_{j}}(\bm{\xi_{m_{B}}},\bm{\xi_{m}},0;\bm{\xi_{i_{B}}},\bm{\xi_{i}},\tau)$
$\displaystyle\hskip
56.9055pt\langle\xi_{i_{B}},\xi_{i}|\rho^{\text{el}}(\tau)\rho^{\text{el}}_{B}(\tau)|\xi^{\prime}_{i},\xi^{\prime}_{i_{B}}\rangle
G^{*}_{\mathbf{x}+dx_{j}}(\bm{\xi^{\prime}_{m_{B}}},\bm{\xi_{m}^{\prime}},0;\bm{\xi^{\prime}_{i_{B}}},\bm{\xi^{\prime}_{i}},\tau)G^{*}_{\mathbf{x}}(\bm{\xi^{\prime}_{f_{B}}},\bm{\xi^{\prime}_{f}},t;\bm{\xi_{m_{B}}^{\prime}},\bm{\xi^{\prime}_{m}},0)\Big{]}.$
(73)

In Eq.(73), coherent states with a $B$ subscript denote bath variables, while
those without a subscript denote system variables. The object,

$\displaystyle
G_{\mathbf{x}}(\bm{\xi_{f_{B}}},\bm{\xi_{f}},t;\bm{\xi_{m_{B}}},\bm{\xi_{m}},0)$
$\displaystyle=\langle\xi_{f_{B}},\xi_{f}|e^{-iH^{\text{el}}(\mathbf{x})t}|\xi_{m},\xi_{m_{B}}\rangle,$
(74)

for example, denotes forward time propagation from coherent states
$|\xi_{m_{B}},\xi_{m}\rangle$ at time $t=0$ to coherent states
$|\xi_{f_{B}},\xi_{f}\rangle$ at time $t$, all while the system is held at
vibrational position $\mathbf{x}$.

Since the total density matrix factorizes in the infinite past,
$\rho^{\text{el}}(\tau)$ can be brought outside the bath trace, and the
expression can be written with a propagator,

$\displaystyle
A_{\mathbf{x}+dx_{j}}(\bm{\xi_{f}},\bm{\xi^{\prime}_{f}};\mathbf{x},t)$
$\displaystyle=\lim_{\tau\rightarrow-\infty}\int(d\bm{\xi_{i}^{*}}d\bm{\xi_{i}})\>(d\bm{\xi_{i}^{*\prime}}d\bm{\xi_{i}^{\prime}})(d\bm{\xi_{m}^{*}}d\bm{\xi_{m}^{*}})\>(d\bm{\xi_{m}^{*\prime}}d\bm{\xi_{m}^{\prime}})e^{-\bm{\xi_{i}^{*}}\bm{\xi_{i}}}e^{-\bm{\xi_{i}^{*\prime}}\bm{\xi^{\prime}_{i}}}e^{-\bm{\xi_{m}^{*}}\bm{\xi_{m}}}e^{-\bm{\xi_{m}^{*\prime}}\bm{\xi_{m}^{\prime}}}$
$\displaystyle\hskip
42.67912pt\times\mathcal{J}_{\mathbf{x}+dx_{j}}\left(\bm{\xi_{f}},\bm{\xi^{\prime}_{f}},t;\bm{\xi_{m}},\bm{\xi^{\prime}_{m}},0;\bm{\xi_{i}},\bm{\xi^{\prime}_{i}},\tau\right)\rho^{\text{el}}(\bm{\xi}_{i},\bm{\xi}^{\prime}_{i};\tau).$
(75)

The propagator is then decomposed into four path integrals over the system
coherent states, with path $(\bm{\xi_{a}},\bm{\xi^{*}_{a}})$ from $\tau$ to
the steady state at time $t=0$, and path $(\bm{\xi_{b}},\bm{\xi^{*}_{b}})$
from $t=0$ to time $t$,

$\displaystyle\mathcal{J}_{\mathbf{x}+dx_{j}}\left(\bm{\xi_{f}},\bm{\xi^{\prime}_{f}},t;\bm{\xi_{m}},\bm{\xi^{\prime}_{m}},0;\bm{\xi_{i}},\bm{\xi^{\prime}_{i}},\tau\right)$
$\displaystyle=\int_{\bm{\xi_{b}}(0)=\bm{\xi_{m}}}^{\bm{\xi_{b}^{*}}(t)=\bm{\xi^{*}_{f}}}\mathcal{D}\left(\bm{\xi_{b}},\bm{\xi_{b}^{*}}\right)\int_{\bm{\xi^{\prime}_{b}}(0)=\bm{\xi^{\prime}_{m}}}^{\bm{\xi^{*\prime}_{b}}(t)=\bm{\xi^{*\prime}_{f}}}\mathcal{D}\left(\bm{\xi_{b}^{\prime}},{\bm{\xi^{*\prime}_{b}}}\right)\int_{\bm{\xi_{a}}(\tau)=\bm{\xi_{i}}}^{\bm{\xi_{a}^{*}}(0)=\bm{\xi^{*}_{m}}}\mathcal{D}\left(\bm{\xi_{a}},\bm{\xi_{a}^{*}}\right)\int_{\bm{\xi^{\prime}_{a}}(\tau)=\bm{\xi^{\prime}_{i}}}^{\bm{\xi^{*\prime}_{a}}(0)=\bm{\xi^{*\prime}_{m}}}\mathcal{D}\left(\bm{\xi_{a}}^{\prime},{\bm{\xi^{*\prime}_{a}}}\right)$
$\displaystyle\>\>\>\>\>\times
e^{i\mathcal{S}_{S}\left[\bm{\xi^{*}_{b}},\bm{\xi_{b}}\>;\>\mathbf{x},t,0\right]}e^{i\mathcal{S}_{S}\left[\bm{\xi^{*}_{a}},\bm{\xi_{a}}\>;\>\mathbf{x}+dx_{j},0,\tau\right]}\mathcal{F}[\bm{\xi_{a}},\bm{\xi^{\prime}_{a}},\bm{\xi^{*}_{a}},\bm{\xi^{*\prime}_{a}},\bm{\xi_{b}},\bm{\xi^{\prime}_{b}},\bm{\xi^{*}_{b}},\bm{\xi^{*\prime}_{b}};t,\tau]$
$\displaystyle\>\>\>\>\>\times
e^{-i\mathcal{S}^{*}_{S}\left[\bm{\xi^{*\prime}_{a}},\bm{\xi^{\prime}_{a}}\>;\>\mathbf{x}+dx_{j},0,\tau\right]}e^{-i\mathcal{S}^{*}_{S}\left[\bm{\xi^{*\prime}_{b}},\bm{\xi^{\prime}_{b}}\>;\>\mathbf{x},t,0\right]}.$
(76)

A schematic of these four path integrals for the forward propagation is shown
in the left plot of Fig. 7, where paths $(\bm{\xi_{a}},\bm{\xi^{*}_{a}})$ and
$(\bm{\xi_{b}},\bm{\xi^{*}_{b}})$ are joined by the midle point,
$(\bm{\xi_{m}},\bm{\xi^{*}_{m}})$. The effect of the bath on each path is
contained in the Feynman-Vernon influence functional,
$\mathcal{F}[\bm{\xi_{a}},\bm{\xi^{\prime}_{a}},\bm{\xi^{*}_{a}},\bm{\xi^{*\prime}_{a}},\bm{\xi_{b}},\bm{\xi^{\prime}_{b}},\bm{\xi^{*}_{b}},\bm{\xi^{*\prime}_{b}};t,\tau]$.
Notice that, because the vibrational coordinates enter only through
$H^{\text{el}}_{S}$, the influence-functional is the same for path
$(\bm{\xi_{a}},\bm{\xi^{*}_{a}})$ and path $(\bm{\xi_{b}},\bm{\xi^{*}_{b}})$.
The difference between the two sets of paths, then, resides exclusively in the
system actions,

$\displaystyle\mathcal{S}_{S}\left[\bm{\xi^{*}_{b}},\bm{\xi_{b}}\>;\>\mathbf{x},t,0\right]$
$\displaystyle=-i\bm{\xi^{*}_{b}}(t),\bm{\xi_{b}}(t)+\int^{t}_{0}d\tau_{1}\left[i\bm{\xi^{*}_{b}}(\tau_{1})\frac{\partial\bm{\xi_{b}}(\tau_{1})}{\partial\tau_{1}}-H^{\text{el}}_{S}\left(\bm{\xi_{b}}(\tau_{1}),\bm{\xi^{*}_{b}}(\tau_{1})\>;\>\mathbf{x}\right)\right]$
(77)
$\displaystyle\mathcal{S}_{S}\left[\bm{\xi^{*}_{a}},\bm{\xi_{a}}\>;\>\mathbf{x}+dx_{j},0,\tau\right]$
$\displaystyle=-i\bm{\xi^{*}_{a}}(0),\bm{\xi_{a}}(0)+\int^{0}_{\tau}d\tau_{1}\left[i\bm{\xi^{*}_{a}}(\tau_{1})\frac{\partial\bm{\xi_{a}}(\tau_{1})}{\partial\tau_{1}}-H^{\text{el}}_{S}\left(\bm{\xi_{a}}(\tau_{1}),\bm{\xi^{*}_{a}}(\tau_{1})\>;\>\mathbf{x}+dx_{j}\right)\right].$
(78)

This greatly simplifies the theory, as now the two distinct paths can be
joined by the integral over the middle point in Eq.(73). The resulting joint
path is shown in the right plot in Fig. 7, with the total propagator and
influence functional now

$\displaystyle\mathcal{J}_{\mathbf{x}+dx_{j}}\left(\bm{\xi_{f}},\bm{\xi}^{\prime}_{f},t;\bm{\xi}_{i},\bm{\xi}^{\prime}_{i},\tau\right)$
$\displaystyle=\int_{\bm{\xi}(\tau)=\bm{\xi_{i}}}^{\bm{\xi^{*}}(t)=\bm{\xi^{*}_{f}}}\mathcal{D}\left(\bm{\xi},\bm{\xi^{*}}\right)\int_{\bm{\xi^{\prime}}(\tau)=\bm{\xi^{\prime}_{i}}}^{\bm{\xi^{*\prime}}(t)=\bm{\xi^{*\prime}_{f}}}\mathcal{D}\left(\bm{\xi^{\prime}},{\bm{\xi^{*\prime}}}\right)$
$\displaystyle\>\>\>\>\>e^{i\mathcal{S}_{S}\left[\bm{\xi^{*}},\bm{\xi}\>;\>\mathbf{x},t,0\right]}e^{i\mathcal{S}_{S}\left[\bm{\xi^{*}},\bm{\xi}\>;\>\mathbf{x}+dx_{j},0,\tau\right]}\mathcal{F}[\bm{\xi},\bm{\xi^{\prime}},\bm{\xi^{*}},\bm{\xi^{*\prime}};t,\tau]e^{-i\mathcal{S}^{*}_{S}\left[\bm{\xi^{*\prime}},\bm{\xi^{\prime}}\>;\>\mathbf{x}+dx_{j},0,\tau\right]}e^{-i\mathcal{S}^{*}_{S}\left[\bm{\xi^{*\prime}},\bm{\xi^{\prime}}\>;\>\mathbf{x},t,0\right]}$
(79)
$\displaystyle\mathcal{F}[\bm{\xi},\bm{\xi^{\prime}},\bm{\xi^{*}},\bm{\xi^{*\prime}};t,\tau]$
$\displaystyle=\int(d\bm{\xi_{i_{B}}^{*}}d\bm{\xi_{i_{B}}})\>(d\bm{\xi_{i_{B}}^{*\prime}}d\bm{\xi_{i_{B}}^{\prime}})\>e^{-\bm{\xi_{i_{B}}^{*}}\bm{\xi_{i_{B}}}}e^{-\bm{\xi_{i_{B}}^{*\prime}}\bm{\xi^{\prime}_{i_{B}}}}$
$\displaystyle\>\>\>\>\>\text{Tr}_{B}\Big{[}G_{B}(\bm{\xi_{f_{B}}},\bm{\xi_{f}},t;\bm{\xi_{i_{B}}},\bm{\xi_{i}},\tau)\langle\xi_{i_{B}}|\rho_{B}(0)|\xi^{\prime}_{i_{B}}\rangle
G_{B}^{*}(\bm{\xi^{\prime}_{f_{B}}},\bm{\xi^{\prime}_{f}},t;\bm{\xi^{\prime}_{i_{B}}},\bm{\xi^{\prime}_{i}},\tau)\Big{]}.$
(80)

Here, $G_{B}(\bm{\xi_{f_{B}}},\bm{\xi_{f}},t;\bm{\xi_{i_{B}}},\bm{\xi_{i}},0)$
is the propagator for the bath density matrix while the system evolves along
path $(\bm{\xi^{*}},\bm{\xi})$.

In the bath interaction picture, where
$H_{SB}^{I}(t)=e^{-iH_{B}t}H_{SB}e^{iH_{B}t}$, this time-evolution can be
alternatively expressed with time-ordered exponentials,

$\displaystyle\mathcal{F}$
$\displaystyle=\text{Tr}_{B}\left\\{\mathcal{T}e^{-i\int^{t}_{\tau}d\tau_{1}H_{SB}^{I}(\tau_{1}\>;\>\bm{\xi}(\tau_{1}),\bm{\xi^{*}}(\tau_{1}))}\rho^{\text{el}}_{B}(\tau)\right.$
$\displaystyle\hskip
28.45274pt\left.\times\mathcal{T}^{-1}e^{i\int^{t}_{\tau}d\tau_{1}H_{SB}^{I}(\tau_{1}\>;\>\bm{\xi^{\prime}}(\tau_{1}),{\bm{\xi^{*\prime}}}(\tau_{1}))}\right\\},$
(81)

where the coherent state notation has been suppressed from $\mathcal{F}$ for
brevity.

The influence functional in Eq.(81) is similar to that found in standard HEOM
theory, except now the initial time is at some point in the distant past,
$\tau$. Following the usual HEOM approach Tanimura and Kubo (1989); Tanimura
(2006); Härtle  _et al._ (2015); Jin  _et al._ (2007, 2008); Zheng  _et al._
(2009); Yan (2014); Wenderoth  _et al._ (2016); Tanimura (2020); Ye  _et al._
(2016) allows the influence functional to be expressed in terms of the two-
time bath-correlation functions, $C^{\sigma}_{\alpha,mm^{\prime}}(t)$, and
their expansion from Eq.(11) in the main text,

$\displaystyle\mathcal{F}$
$\displaystyle=\exp\Big{(}-i\int^{t}_{\tau}d\tau_{1}\sum_{j}\mathcal{A}^{\bar{\sigma}}_{\alpha,m}\left[\tau_{1};\bm{\xi},\bm{\xi}^{\prime}\right]\mathcal{B}_{j}\left[\tau_{1},\tau;\bm{\xi},\bm{\xi}^{\prime}\right]\Big{)}.$
(82)

Here, two new Grassmann variables have been introduced,

$\displaystyle\mathcal{B}_{j}\left[\bm{\xi},\bm{\xi}^{\prime};t,\tau\right]$
$\displaystyle=-iV_{\alpha,m}\left(\eta_{j}\int^{t}_{\tau}d\tau_{1}\>e^{-\kappa_{j}\tau_{1}}{\xi}^{\sigma}_{m}(t-\tau_{1})\right.$
$\displaystyle\hskip
56.9055pt\left.-\eta^{*}_{j}\int^{t}_{\tau}d\tau_{1}\>e^{-\kappa_{j}\tau_{1}}{\xi}^{\sigma\prime}_{m}(t-\tau_{1})\right)$
(83)
$\displaystyle\mathcal{A}^{\sigma}_{\alpha,m}\left[\bm{\xi},\bm{\xi}^{\prime};t\right]$
$\displaystyle=V_{\alpha,m}\left[\xi^{\sigma}_{m}(t)+\xi^{\sigma\prime}_{m}(t)\right],$
(84)

along with the super-index $j=(\alpha,\sigma,\ell,m)$ from the main text. One
can now construct the standard hierarchical equations of motion for
$\mathcal{F}$,

$\displaystyle\frac{\partial}{\partial t}\mathcal{F}^{(n)}_{\mathbf{j}}$
$\displaystyle=-\left(\sum_{r=1}^{N}\kappa_{j_{r}}\right)\mathcal{F}^{(n)}_{\mathbf{j}}-i\sum_{r=1}^{n}(-1)^{n-r}\mathcal{C}_{j_{r}}\mathcal{F}^{(n-1)}_{\mathbf{j}^{-}}$
$\displaystyle\hskip
14.22636pt-i\sum_{j}\mathcal{A}^{\bar{\sigma}}_{\alpha,m}\mathcal{F}^{(n+1)}_{\mathbf{j}^{+}},$
(85)

where the $n$th-tier influence functional is formed by application of $n$
Grassmann variables
$\mathcal{B}_{j}\left[\bm{\xi},\bm{\xi^{\prime}};t,\tau\right]$,

$\displaystyle\mathcal{F}^{(n)}_{\mathbf{j}}$
$\displaystyle=\mathcal{B}_{j_{n}}\dots\mathcal{B}_{j_{1}}\mathcal{F},$ (86)

such that the $0$th tier corresponds to the original influence functional,
$\mathcal{F}^{(0)}=\mathcal{F}$.

Now explicitly including the coherent state notation, a corresponding $n$th
tier propagator and auxiliary operator (AO) can be respecetively defined,

$\displaystyle\mathcal{J}^{(n)}_{\mathbf{j},\mathbf{x}+dx_{j}}\left(\bm{\xi_{f}},\bm{\xi}^{\prime}_{f},t;\bm{\xi}_{i},\bm{\xi}^{\prime}_{i},\tau\right)$
$\displaystyle=\int_{\bm{\xi}(\tau)=\bm{\xi_{i}}}^{\bm{\xi^{*}}(t)=\bm{\xi^{*}_{f}}}\mathcal{D}\left(\bm{\xi},\bm{\xi^{*}}\right)\int_{\bm{\xi^{\prime}}(\tau)=\bm{\xi^{\prime}_{i}}}^{\bm{\xi^{*\prime}}(t)=\bm{\xi^{*\prime}_{f}}}\mathcal{D}\left(\bm{\xi^{\prime}},{\bm{\xi^{*\prime}}}\right)e^{i\mathcal{S}_{S}\left[\bm{\xi^{*}},\bm{\xi}\>;\>\mathbf{x},t,0\right]}e^{i\mathcal{S}_{S}\left[\bm{\xi^{*}},\bm{\xi}\>;\>\mathbf{x}+dx_{j},0,\tau\right]}$
$\displaystyle\>\>\>\>\>\times\mathcal{F}^{(n)}_{\mathbf{j}}[\bm{\xi},\bm{\xi^{\prime}},\bm{\xi^{*}},\bm{\xi^{*\prime}};t,\tau]e^{-i\mathcal{S}^{*}_{S}\left[\bm{\xi^{*\prime}},\bm{\xi^{\prime}}\>;\>\mathbf{x}+dx_{j},0,\tau\right]}e^{-i\mathcal{S}^{*}_{S}\left[\bm{\xi^{*\prime}},\bm{\xi^{\prime}}\>;\>\mathbf{x},t,0\right]}$
(87) $\displaystyle
A^{(n)}_{\mathbf{j},\mathbf{x}+dx_{j}}(\bm{\xi_{f}},\bm{\xi^{\prime}_{f}};\mathbf{x},t)$
$\displaystyle=\lim_{\tau\rightarrow-\infty}\int(d\bm{\xi_{i}^{*}}d\bm{\xi_{i}})\>(d\bm{\xi_{i}^{*\prime}}d\bm{\xi_{i}^{\prime}})e^{-\bm{\xi_{i}^{*}}\bm{\xi_{i}}}e^{-\bm{\xi_{i}^{*\prime}}\bm{\xi^{\prime}_{i}}}\mathcal{J}^{(n)}_{\mathbf{j},\mathbf{x}+dx_{j}}\left(\bm{\xi_{f}},\bm{\xi^{\prime}_{f}},t;\bm{\xi_{i}},\bm{\xi^{\prime}_{i}},\tau\right)\rho^{\text{el}}(\bm{\xi_{i}},\bm{\xi^{\prime}_{i}};\tau).$
(88)

Consider now the case at $t=0$. In this instance, two of the four system
actions vanish,
$\mathcal{S}_{S}\left[\bm{\xi^{*}},\bm{\xi}\>;\>\mathbf{x},0,0\right]=0$, and,
after moving to a more computationally suitable basis, the AOs propagate
according to

$\displaystyle\frac{\partial}{\partial
t^{\prime}}{A}^{(n)}_{\mathbf{j},\mathbf{x}+dx_{j}}(t^{\prime})$
$\displaystyle=i\left[H_{\text{S}}(\mathbf{x}+dx_{j}),A^{(n)}_{\mathbf{j},\mathbf{x}+dx_{j}}\right]-\left(\sum_{r=1}^{N}\kappa_{j_{r}}\right)A^{(n)}_{\mathbf{j},\mathbf{x}+dx_{j}}+i\sum_{r=1}^{n}(-1)^{n-r}\mathcal{C}_{j_{r}}A^{(n-1)}_{\mathbf{j}^{-},\mathbf{x}+dx_{j}}+i\sum_{j}\mathcal{A}^{\bar{\sigma}}_{\alpha,m}A^{(n+1)}_{\mathbf{j}^{+},\mathbf{x}+dx_{j}},$
(89)

where $\tau<t^{\prime}\leq 0$. Eq.(89) is solved subject to the initial
condition $A^{(0)}_{\mathbf{x}+dx_{j}}(\tau)=\rho^{\text{el}}(\tau)$, which is
the initial state of the system electronic degrees of freedom in the distant
past, and $A^{(n>0)}_{\mathbf{j},\mathbf{x}+dx_{j}}(\tau)=0$. If one solves
these equations of motion for

$\displaystyle\bm{A}(\mathbf{x}+dx_{j},0)$
$\displaystyle=\left[A^{(0)}_{\mathbf{x}+dx_{j}}(0),A^{(1)}_{j_{1},\mathbf{x}+dx_{j}}(0),\dots\right],$
(90)

then one is simply calculating

$\displaystyle\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}(\mathbf{x}+dx_{j})$
$\displaystyle=\left[\sigma^{(0),\text{el}}_{\text{ss}}(\mathbf{x}+dx_{j}),\sigma^{(1),\text{el}}_{j_{1},\text{ss}}(\mathbf{x}+dx_{j}),\dots\right],$
(91)

the steady of all ADOs for vibrational frame $\mathbf{x}+dx_{j}$. This entire
process can be repeated for the second term in Eq.(67),
${A}_{\mathbf{x}}(\mathbf{x},t)$, which just yields
$\bm{\tilde{\sigma}}^{\text{el}}_{ss}(\mathbf{x})$.

With these results, the propagation from time $t=0$ to time $t$ can be easily
treated. First, note that both ${A}_{\mathbf{x}}(\mathbf{x},t)$ and
${A}_{\mathbf{x}+dx_{j}}(\mathbf{x},t)$ will follow the same equation of
motion for vibrational frame $\mathbf{x}$ during this interval. Since the AOs
at $t=0$ are known from the previous section, furthermore,

$\displaystyle\lim_{\tau\rightarrow-\infty}A^{(n)}_{\mathbf{j},\mathbf{x}+dx_{j}}(0,\tau)$
$\displaystyle=\sigma^{(n),\text{el}}_{\mathbf{j},\text{ss}}(\mathbf{x}+dx_{j})$
$\displaystyle\lim_{\tau\rightarrow-\infty}A^{(n)}_{\mathbf{j},\mathbf{x}}(0,\tau)$
$\displaystyle=\sigma^{(n),\text{el}}_{\mathbf{j},\text{ss}}(\mathbf{x}),$
(92)

both can just be propagated from this point to time $t$. After putting both
the $\mathbf{A}_{\mathbf{x}+dx_{j}}(\mathbf{x},t)$ and
$\mathbf{A}_{\mathbf{x}}(\mathbf{x},t)$ back together, then, Eq.(66) becomes

$\displaystyle
A^{(n)}_{\mathbf{j},\mathbf{x}+dx_{j}}(\bm{\xi_{f}},\bm{\xi^{\prime}_{f}};\mathbf{x},t)$
$\displaystyle=\lim_{\tau\rightarrow-\infty}\int(d\bm{\xi_{i}^{*}}d\bm{\xi_{i}})\>(d\bm{\xi_{i}^{*\prime}}d\bm{\xi_{i}^{\prime}})e^{-\bm{\xi_{i}^{*}}\bm{\xi_{i}}}e^{-\bm{\xi_{i}^{*\prime}}\bm{\xi^{\prime}_{i}}}\mathcal{J}^{(n)}_{\mathbf{j},\mathbf{x}}\left(\bm{\xi_{f}},\bm{\xi^{\prime}_{f}},t;\bm{\xi_{i}},\bm{\xi^{\prime}_{i}},0\right)\frac{\partial\sigma^{(n),\text{el}}_{\mathbf{j},\text{ss}}(\bm{\xi_{i}},\bm{\xi^{\prime}_{i}};\mathbf{x})}{\partial
x_{j}}.$ (93)

Again, returning to a more computationally suitable basis yields the
corresponding HEOM,

$\displaystyle\frac{\partial}{\partial t}A^{(n)}_{\mathbf{j}}(\mathbf{x},t)$
$\displaystyle=-i\left[H_{\text{S}}(\mathbf{x}),A^{(n)}_{\mathbf{j}}\right]-\left(\sum_{r=1}^{N}\kappa_{j_{r}}\right)A^{(n)}_{\mathbf{j}}$
$\displaystyle\>\>\>\>\>-i\sum_{r=1}^{n}(-1)^{n-r}\mathcal{C}_{j_{r}}A^{(n-1)}_{\mathbf{j}^{-}}$
$\displaystyle\>\>\>\>\>-i\sum_{j}\mathcal{A}^{\bar{\sigma}}_{\alpha,m}A^{(n+1)}_{\mathbf{j}^{+}},$
(94)

an equation of motion that is solved subject to the initial condition

$\displaystyle A^{(n)}_{\mathbf{j}}(\mathbf{x},0)$
$\displaystyle=\frac{\partial\bm{\tilde{\sigma}}^{\text{el}}_{ss}(\mathbf{x})}{\partial
x_{j}}.$ (95)

Note that this is exactly the same HEOM one would obtain for the time-
evolution of a reduced system density matrix based on a factorized initial
condition, except that now it time-evolves the spatial derivative of the
steady state of all ADOs.

Evidently, therefore, both questions from the start of this section have been
answered. First, when the vibrational degrees of freedom are only in the
system Hamiltonian, $H^{\text{el}}_{S}$, one extracts the same information
from $\bm{\tilde{\sigma}}^{\text{el}}_{ss}(\mathbf{x})$ as from
$\rho^{\text{el}}_{\text{tot},\text{ss}}(\mathbf{x})$. Second, the resulting
time-propagation is under the same HEOM as that derived for a factorized
initial state. Finally, to fully connect this with the form of the friction of
the main text, consider the HEOM in Eq.(94) in a joint Liouville space, where
it has the solution

$\displaystyle\bm{A}(\mathbf{x},t)$
$\displaystyle=e^{-\mathcal{L}^{\text{el}}(\mathbf{x})t}\frac{\partial\bm{\tilde{\sigma}}^{\text{el}}_{ss}(\mathbf{x})}{\partial
x_{j}}.$ (96)

Here, $\mathcal{L}^{\text{el}}(\mathbf{x})$ is a superoperator containing the
numerically exact dynamics in Eq.(94) and is, by definition, the same quantity
as in Eq.(25). The friction in Eq.(60), consequently, is now exactly the same
expression as that derived in the main text.

###  B.3 Diffusion tensor

The equivalence of the diffusion tensors can be shown through a similar
process. The only difference is that an operator in the system space, $\delta
F_{j}(\mathbf{x})$, is applied at time $t=0$, instead of a classical spatial
derivative, which is the same problem that the authors already considered in
Ref. Ke  _et al._ (2022).

##  Appendix C Non-Markovian friction tensor

The friction and diffusion tensors derived in the previous section are
Markovian, in that they depend only on the vibrational position at time $t$.
In general, introducing non-Markovian dynamics for this problem is nontrivial;
however, a simple form of memory can be introduced by setting
$\frac{\partial\mathbf{\tilde{B}}^{\text{el}}}{\partial t}\neq 0$, as in Ref.
Dou  _et al._ (2017). This still keeps the timescale separation between
electronic and vibrational degrees of freedom, in that the time-evolution of
$\mathbf{\tilde{B}}^{\text{el}}(\mathbf{x},\mathbf{p};t)$ and
$A(\mathbf{x},\mathbf{p};t)$ is still for a fixed vibrational frame, but now
some relaxation effects from the electronic degrees of freedom are retained.
Again suppressing the $(\mathbf{x},\mathbf{p};t)$ notation, this can be seen
in the equation of motion for
$\mathbf{\tilde{B}}^{\text{el}}(\mathbf{x},\mathbf{p};t)$,

$\displaystyle\frac{\partial\mathbf{\tilde{B}}^{\text{el}}}{\partial t}$
$\displaystyle=\\{\mskip-3.0mu\mskip-3.0mu\\{H^{\text{el}}_{\text{S}},A\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\\}\mskip-3.0mu\mskip-3.0mu\\}_{a}-\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\text{Tr}_{\text{el}}[\\{\mskip-3.0mu\mskip-3.0mu\\{H^{\text{el}}_{\text{S}},A\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\\}\mskip-3.0mu\mskip-3.0mu\\}_{a}]-\mathcal{L}^{\text{el}}\mathbf{\tilde{B}}^{\text{el}}.$
(97)

After solving this first-order differential equation and assuming that the
system starts in the electronic steady state, such that
$\mathbf{\tilde{B}}^{\text{el}}(\mathbf{x},\mathbf{p};0)=0$, time-dependent
friction and diffusion tensors,

$\displaystyle\bar{\gamma}_{ij}(t-\tau)$
$\displaystyle=-\text{Tr}_{\text{el}}\left[\frac{\partial}{\partial
x_{i}}e^{-\mathcal{L}^{\text{el}}(\mathbf{x})(t-\tau)}\frac{\partial\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}(\mathbf{x})}{\partial
x_{j}}\right]$ (98) $\displaystyle\bar{D}_{ij}(t-\tau)$
$\displaystyle=\frac{1}{2}\text{Tr}_{\text{el}}\Big{[}\delta
F_{i}e^{-\mathcal{L}^{\text{el}}(\mathbf{x})(t-\tau)}\left(\delta
F_{j}\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}+\bm{\tilde{\sigma}}^{\text{el}}_{\text{ss}}\delta
F_{j}\right)\Big{]},$ (99)

are derived. They form part of the resulting non-Markovian Fokker-Planck
equation,

$\displaystyle\frac{\partial A(t)}{\partial t}$
$\displaystyle=-\sum_{i}\frac{p_{i}}{m_{i}}\frac{\partial A(t)}{\partial
x_{i}}+\sum_{i}F_{i}(\mathbf{x})\frac{\partial A(t)}{\partial
p_{i}}+\sum_{ij}\int^{t}_{0}d\tau\bar{\gamma}_{ij}(t-\tau)\frac{\partial}{\partial
p_{i}}\left(\frac{p_{j}}{m_{j}}A(\tau)\right)+\sum_{ij}\int^{t}_{0}d\tau\bar{D}_{ij}(t-\tau)\frac{\partial^{2}A(\tau)}{\partial
p_{i}\partial p_{j}},$ (100)

Note that this non-Markovian friction is not equivalent to that derived in
Ref. Chen  _et al._ (2019a), where all non-Markovian effects are rigorously
included and the full history of the classical trajectory is used.

## References

  * Tal  _et al._ (2008) O. Tal, M. Krieger, B. Leerink, and J. M. van Ruitenbeek, Phys. Rev. Lett. 100, 196804 (2008).
  * Secker  _et al._ (2011) D. Secker, S. Wagner, S. Ballmann, R. Härtle, M. Thoss, and H. B. Weber, Phys. Rev. Lett. 106, 136807 (2011).
  * Craighead (2000) H. G. Craighead, Science 290, 1532 (2000).
  * Ekinci and Roukes (2005) K. L. Ekinci and M. L. Roukes, Rev. Sci. Instrum. 76, 061101 (2005).
  * Weig  _et al._ (2004) E. M. Weig, R. H. Blick, T. Brandes, J. Kirschbaum, W. Wegscheider, M. Bichler, and J. P. Kotthaus, Phys. Rev. Lett. 92, 046804 (2004).
  * Sapmaz  _et al._ (2006) S. Sapmaz, P. Jarillo-Herrero, Y. M. Blanter, C. Dekker, and H. S. J. van der Zant, Phys. Rev. Lett. 96, 026801 (2006).
  * Koch and von Oppen (2005) J. Koch and F. von Oppen, Phys. Rev. Lett. 94, 206804 (2005).
  * Koch  _et al._ (2006a) J. Koch, F. von Oppen, and A. V. Andreev, Phys. Rev. B 74, 205438 (2006a).
  * Koch  _et al._ (2006b) J. Koch, M. Semmelhack, F. von Oppen, and A. Nitzan, Phys. Rev. B 73, 155306 (2006b).
  * Leturcq  _et al._ (2009) R. Leturcq, C. Stampfer, K. Inderbitzin, L. Durrer, C. Hierold, E. Mariani, M. G. Schultz, F. von Oppen, and K. Ensslin, Nat. Phys. 5, 327 (2009).
  * Härtle and Thoss (2011a) R. Härtle and M. Thoss, Phys. Rev. B 83, 115414 (2011a).
  * Härtle  _et al._ (2018) R. Härtle, C. Schinabeck, M. Kulkarni, D. Gelbwaser-Klimovsky, M. Thoss, and U. Peskin, Phys. Rev. B 98, 081404(R) (2018).
  * Härtle and Thoss (2011b) R. Härtle and M. Thoss, Phys. Rev. B 83, 125419 (2011b).
  * Kaspar  _et al._ (2022) C. Kaspar, A. Erpenbeck, J. Bätge, C. Schinabeck, and M. Thoss, Phys. Rev. B 105, 195435 (2022).
  * Kosov (2017) D. S. Kosov, J. Chem. Phys. 147, 104109 (2017).
  * Rudge and Kosov (2019a) S. L. Rudge and D. S. Kosov, Phys. Rev. B 99, 115426 (2019a).
  * Schinabeck  _et al._ (2014) C. Schinabeck, R. Härtle, H. B. Weber, and M. Thoss, Phys. Rev. B 90, 075409 (2014).
  * Mitra  _et al._ (2004) A. Mitra, I. Aleiner, and A. J. Millis, Phys. Rev. B 69, 245302 (2004).
  * Frederiksen  _et al._ (2007) T. Frederiksen, M. Paulsson, M. Brandbyge, and A.-P. Jauho, Phys. Rev. B 75, 205413 (2007).
  * Galperin  _et al._ (2007) M. Galperin, M. A. Ratner, and A. Nitzan, J. Phys. Condens. Matter 19, 103201 (2007).
  * Härtle  _et al._ (2008) R. Härtle, C. Benesch, and M. Thoss, Phys. Rev. B 77, 205314 (2008).
  * Novotný  _et al._ (2011) T. Novotný, F. Haupt, and W. Belzig, Phys. Rev. B 84, 113107 (2011).
  * Dash  _et al._ (2011) L. K. Dash, H. Ness, and R. W. Godby, Phys. Rev. B 84, 085433 (2011).
  * Agarwalla  _et al._ (2015) B. K. Agarwalla, J.-H. Jiang, and D. Segal, Phys. Rev. B 92, 245418 (2015).
  * Ryndyk  _et al._ (2006) D. A. Ryndyk, M. Hartung, and G. Cuniberti, Phys. Rev. B 73, 045420 (2006).
  * Martin-Rodero  _et al._ (2008) A. Martin-Rodero, A. Levy Yeyati, F. Flores, and R. C. Monreal, Phys. Rev. B 78, 235112 (2008).
  * Entin-Wohlman  _et al._ (2009) O. Entin-Wohlman, Y. Imry, and A. Aharony, Phys. Rev. B 80, 035417 (2009).
  * Egger and Gogolin (2008) R. Egger and A. O. Gogolin, Phys. Rev. B 77, 113405 (2008).
  * Galperin  _et al._ (2004) M. Galperin, M. A. Ratner, and A. Nitzan, J. Chem. Phys. 121, 11965 (2004).
  * Tikhodeev  _et al._ (2001) S. Tikhodeev, M. Natario, K. Makoshi, T. Mii, and H. Ueba, Surf. Sci. 493, 63 (2001).
  * Paulsson  _et al._ (2005) M. Paulsson, T. Frederiksen, and M. Brandbyge, Phys. Rev. B 72, 201101(R) (2005).
  * Li  _et al._ (2005) X. Q. Li, J. Luo, Y. G. Yang, P. Cui, and Y. J. Yan, Phys. Rev. B 71, 205304 (2005).
  * Leijnse  _et al._ (2010) M. Leijnse, M. R. Wegewijs, and K. Flensberg, Phys. Rev. B 82, 045412 (2010).
  * Rudge and Kosov (2021) S. L. Rudge and D. S. Kosov, Phys. Rev. B 104, 165307 (2021).
  * Mühlbacher and Rabani (2008) L. Mühlbacher and E. Rabani, Phys. Rev. Lett. 100, 176403 (2008).
  * Weiss  _et al._ (2008) S. Weiss, J. Eckel, M. Thorwart, and R. Egger, Phys. Rev. B 77, 195316 (2008).
  * Wang  _et al._ (2011) H. Wang, I. Pshenichnyuk, R. Härtle, and M. Thoss, J. Chem. Phys. 135, 244506 (2011).
  * Wang and Thoss (2016) H. Wang and M. Thoss, J. Chem. Phys. 145, 164105 (2016).
  * Schinabeck  _et al._ (2016) C. Schinabeck, A. Erpenbeck, R. Härtle, and M. Thoss, Phys. Rev. B 94, 201407(R) (2016).
  * Schinabeck  _et al._ (2018) C. Schinabeck, R. Härtle, and M. Thoss, Phys. Rev. B 97, 235429 (2018).
  * Schinabeck and Thoss (2020) C. Schinabeck and M. Thoss, Phys. Rev. B 101, 075422 (2020).
  * Jin  _et al._ (2008) J. Jin, X. Zheng, and Y. Yan, J. Chem. Phys. 128, 234703 (2008).
  * Tully (1990) J. C. Tully, J. Chem. Phys. 93, 1061 (1990).
  * Stock and Thoss (1997) G. Stock and M. Thoss, Phys. Rev. Lett. 78, 578 (1997).
  * Martinez  _et al._ (1996) T. J. Martinez, M. Ben-Nun, and R. D. Levine, J. Phys. Chem. 100, 7884 (1996).
  * Dou and Subotnik (2018) W. Dou and J. E. Subotnik, J. Chem. Phys. 148, 230901 (2018).
  * Dzhioev and Kosov (2011) A. A. Dzhioev and D. S. Kosov, J. Chem. Phys. 135, 074701 (2011).
  * Dzhioev  _et al._ (2013) A. A. Dzhioev, D. S. Kosov, and F. von Oppen, J. Chem. Phys. 138, 134103 (2013).
  * Preston  _et al._ (2021) R. J. Preston, M. F. Gelin, and D. S. Kosov, J. Chem. Phys. 154, 114108 (2021).
  * Mozyrsky  _et al._ (2006) D. Mozyrsky, M. B. Hastings, and I. Martin, Phys. Rev. B 73, 035104 (2006).
  * Kershaw and Kosov (2020) V. F. Kershaw and D. S. Kosov, J. Chem. Phys. 153, 154101 (2020).
  * Bode  _et al._ (2012) N. Bode, S. V. Kusminskiy, R. Egger, and F. von Oppen, Beilstein J. Nanotechnol. 3, 144 (2012).
  * Lü  _et al._ (2010) J.-T. Lü, M. Brandbyge, and P. Hedegård, Nano Lett. 10, 1657 (2010).
  * Lü  _et al._ (2012) J.-T. Lü, M. Brandbyge, P. Hedegård, T. N. Todorov, and D. Dundas, Phys. Rev. B 85, 245444 (2012).
  * Todorov  _et al._ (2011) T. N. Todorov, D. Dundas, A. T. Paxton, and A. P. Horsfield, Beilstein J. Nanotechnol. 2, 727 (2011).
  * Todorov  _et al._ (2014) T. N. Todorov, D. Dundas, J.-T. Lü, M. Brandbyge, and P. Hedegård, Eur. J. Phys. 35, 065004 (2014).
  * Maurer  _et al._ (2016) R. J. Maurer, M. Askerka, V. S. Batista, and J. C. Tully, Phys. Rev. B 94, 115432 (2016).
  * Askerka  _et al._ (2016) M. Askerka, R. J. Maurer, V. S. Batista, and J. C. Tully, Phys. Rev. Lett. 116, 217601 (2016).
  * Bustos-Marún  _et al._ (2013) R. Bustos-Marún, G. Refael, and F. von Oppen, Phys. Rev. Lett. 111, 060802 (2013).
  * Fernández-Alcázar  _et al._ (2019) L. J. Fernández-Alcázar, H. M. Pastawski, and R. A. Bustos-Marún, Phys. Rev. B 99, 045403 (2019).
  * Calvo  _et al._ (2017) H. L. Calvo, F. D. Ribetto, and R. A. Bustos-Marún, Phys. Rev. B 96, 165309 (2017).
  * Sabater  _et al._ (2015) C. Sabater, C. Untiedt, and J. M. van Ruitenbeek, Beilstein J. Nanotechnol. 6, 2338 (2015).
  * Ring  _et al._ (2020) M. Ring, D. Weber, P. Haiber, F. Pauly, P. Nielaba, and E. Scheer, Nano Lett. 20, 5773 (2020), pMID: 32589039.
  * Hussein  _et al._ (2010) R. Hussein, A. Metelmann, P. Zedler, and T. Brandes, Phys. Rev. B 82, 165406 (2010).
  * Metelmann and Brandes (2011) A. Metelmann and T. Brandes, Phys. Rev. B 84, 155455 (2011).
  * Lü  _et al._ (2011) J.-T. Lü, P. Hedegård, and M. Brandbyge, Phys. Rev. Lett. 107, 046801 (2011).
  * Brandbyge  _et al._ (1995) M. Brandbyge, P. Hedegård, T. F. Heinz, J. A. Misewich, and D. M. Newns, Phys. Rev. B 52, 6042 (1995).
  * Caldeira and Leggett (1983) A. Caldeira and A. Leggett, Phys. A: Stat. Mech. Appl 121, 587 (1983).
  * Schmid (1982) A. Schmid, J. Low Temp. Phys. 49, 609 (1982).
  * Head‐Gordon and Tully (1995) M. Head‐Gordon and J. C. Tully, J. Chem. Phys. 103, 10137 (1995).
  * Meyer and Miller (1979) H. Meyer and W. H. Miller, J. Chem. Phys. 70, 3214 (1979).
  * Dou  _et al._ (2017) W. Dou, G. Miao, and J. E. Subotnik, Phys. Rev. Lett. 119, 046001 (2017).
  * Dou and Subotnik (2017a) W. Dou and J. E. Subotnik, Phys. Rev. B 96, 104305 (2017a).
  * Chen  _et al._ (2019a) F. Chen, K. Miwa, and M. Galperin, J. Phys. Chem. A 123, 693 (2019a).
  * Chen  _et al._ (2019b) F. Chen, K. Miwa, and M. Galperin, J. Chem. Phys. 150, 174101 (2019b).
  * Martinazzo and Burghardt (2022a) R. Martinazzo and I. Burghardt, Phys. Rev. A 105, 052215 (2022a).
  * Martinazzo and Burghardt (2022b) R. Martinazzo and I. Burghardt, Phys. Rev. Lett. 128, 206002 (2022b).
  * Dou  _et al._ (2015) W. Dou, A. Nitzan, and J. E. Subotnik, J. Chem. Phys. 143, 054103 (2015).
  * Dou and Subotnik (2016) W. Dou and J. E. Subotnik, J. Chem. Phys. 145, 054102 (2016).
  * Kapral and Ciccotti (1999) R. Kapral and G. Ciccotti, J. Chem. Phys. 110, 8919 (1999).
  * İmre  _et al._ (1967) K. İmre, E. Özizmir, M. Rosenbaum, and P. F. Zweifel, J. Math. Physics 8, 1097 (1967).
  * Dou and Subotnik (2017b) W. Dou and J. E. Subotnik, J. Chem. Phys. 146, 092304 (2017b).
  * Tanimura and Kubo (1989) Y. Tanimura and R. Kubo, J. Phys. Soc. Jpn. 58, 101 (1989).
  * Tanimura (2006) Y. Tanimura, J. Phys. Soc. Jpn 75, 082001 (2006).
  * Härtle  _et al._ (2015) R. Härtle, G. Cohen, D. R. Reichman, and A. J. Millis, Phys. Rev. B 92, 085430 (2015).
  * Jin  _et al._ (2007) J. Jin, S. Welack, J. Luo, X.-Q. Li, P. Cui, R.-X. Xu, and Y. Yan, J. Chem. Phys. 126, 134113 (2007).
  * Zheng  _et al._ (2009) X. Zheng, J. Jin, S. Welack, M. Luo, and Y. Yan, J. Chem. Phys. 130, 164708 (2009).
  * Yan (2014) Y. Yan, J. Chem. Phys. 140, 054105 (2014).
  * Wenderoth  _et al._ (2016) S. Wenderoth, J. Bätge, and R. Härtle, Phys. Rev. B 94, 121303(R) (2016).
  * Tanimura (2020) Y. Tanimura, J. Chem. Phys. 153, 020901 (2020).
  * Ye  _et al._ (2016) L. Ye, X. Wang, D. Hou, R.-X. Xu, X. Zheng, and Y. Yan, WIREs Comput. Mol. Sci. 6, 608 (2016).
  * Tanimura (2014) Y. Tanimura, J. Chem. Phys. 141, 044114 (2014).
  * Ke  _et al._ (2022) Y. Ke, C. Kaspar, A. Erpenbeck, U. Peskin, and M. Thoss, J. Chem. Phys. 157, 034103 (2022).
  * Ozaki (2007) T. Ozaki, Phys. Rev. B 75, 035123 (2007).
  * Verzijl  _et al._ (2013) C. J. O. Verzijl, J. S. Seldenthuis, and J. M. Thijssen, J. Chem. Phys. 138, 094102 (2013).
  * Xu  _et al._ (2022) M. Xu, Y. Yan, Q. Shi, J. Ankerhold, and J. T. Stockburger, Phys. Rev. Lett. 129, 230601 (2022).
  * Dan  _et al._ (2022) X. Dan, M. Xu, J. T. Stockburger, J. Ankerhold, and Q. Shi, arXiv:2211.04089 (preprint) (2022).
  * Rudge and Kosov (2019b) S. L. Rudge and D. S. Kosov, J. Chem. Phys. 151, 034107 (2019b).
  * Erpenbeck  _et al._ (2020a) A. Erpenbeck, Y. Ke, U. Peskin, and M. Thoss, Phys. Rev. B 102, 195421 (2020a).
  * Dou and Subotnik (2017c) W. Dou and J. E. Subotnik, J. Chem. Phys. 146, 092304 (2017c).
  * Preston  _et al._ (2022) R. J. Preston, T. D. Honeychurch, and D. S. Kosov, Phys. Rev. B 106, 195406 (2022).
  * Zheng  _et al._ (2008) X. Zheng, J. Jin, and Y. Yan, J. Chem. Phys. 129, 184112 (2008).
  * Erpenbeck  _et al._ (2020b) A. Erpenbeck, Y. Ke, U. Peskin, and M. Thoss, Phys. Rev. B 102, 195421 (2020b).
  * Plihal and Langreth (1999) M. Plihal and D. C. Langreth, Phys. Rev. B 60, 5969 (1999).
  * Plihal and Langreth (1998) M. Plihal and D. C. Langreth, Phys. Rev. B 58, 2191 (1998).
  * Chen  _et al._ (2022) Z.-H. Chen, Y. Wang, X. Zheng, R.-X. Xu, and Y. Yan, J. Chem. Phys. 156, 221102 (2022).
  * Nakatsukasa  _et al._ (2018) Y. Nakatsukasa, O. Sète, and L. N. Trefethen, SIAM J. Sci. Comput. 40, A1494 (2018).
  * Kongsuwan  _et al._ (2018) N. Kongsuwan, A. Demetriadou, R. Chikkaraddy, F. Benz, V. A. Turek, U. F. Keyser, J. J. Baumberg, and O. Hess, ACS Photonics 5, 186 (2018).
  * Chikkaraddy  _et al._ (2016) R. Chikkaraddy, B. de Nijs, F. Benz, S. J. Barrow, O. A. Scherman, E. Rosta, A. Demetriadou, P. Fox, O. Hess, and J. J. Baumberg, Nature 535, 127 (2016).
  * Lang and Firsov (1963) I. G. Lang and Y. A. Firsov, J. Exp. Theor. Phys. 16, 1301 (1963).