Electronic friction and Langevin dynamics is a popular mixed quantum–classical method for simulating the nonadiabatic dynamics of molecules interacting with metal surfaces, as it can be computationally more efficient than fully quantum approaches. In this work, we extend the theory of electronic friction within the hierarchical equations of motion formalism to models with a position-dependent metal–molecule coupling. We show that the addition of a position-dependent metal–molecule coupling adds new contributions to the electronic friction and other forces, which are highly relevant for many physical processes. Our expressions for the electronic forces within the Langevin equation are valid both in and out of equilibrium and for molecular models containing strong interactions. We demonstrate the approach by applying it to different models of interest.

## I. INTRODUCTION

Investigating the dynamics of molecules interacting with metal surfaces is a crucial task for understanding many phenomena in physics and chemistry. This broad research area describes a plethora of scenarios, such as the scattering of molecules off surfaces^{1–4} or chemical reactions and catalytic processes.^{5–7} When multiple surfaces are considered, moreover, this also includes nonequilibrium transport in molecular nanojunctions.^{8–13}

However, simulating the dynamics of such systems is a challenging task. One must not only describe the interaction between the molecule and the metal surfaces accurately, but also the interaction between the electronic and vibrational degrees of freedom (DOFs) within the molecule itself. In the context of molecular nanojunctions, for example, the electronic–vibrational interaction can lead to interesting transport phenomena such as current rectification, negative differential resistance, and current induced bond rupture.^{14–22} Due to the continuum of electronic states within the metal surfaces, such electronic–vibrational interactions often result in highly nonadiabatic dynamics that requires treatment beyond the Born–Oppenheimer approximation.

Although the most accurate approaches, such as the hierarchical equations of motion (HEOM)^{23–28} and the multilayer multiconfigurational time-dependent Hartree method,^{29–31} treat all DOFs quantum mechanically, they can be numerically costly. This can occur especially in regimes where one needs a large vibrational basis, such as for heavy molecules or large-amplitude vibrational motion. To circumvent this problem, one can use a mixed quantum–classical method, in which vibrational DOFs within the molecule are treated classically while influenced by quantum mechanical electronic DOFs, both in the molecule and the surfaces. Beyond numerical efficiency, such treatments have the additional benefit that they often give valuable insight into the electronic forces acting on the molecular vibrations.^{32,33} This information can be especially insightful in the context of vibrational instabilities^{34–36} and can be difficult to extract from a fully quantum mechanical description of the system. Within the family of such mixed quantum–classical approaches, popular methods include Ehrenfest dynamics,^{16,35,37–45} surface hopping,^{46–49} and Langevin dynamics.^{50–62}

Within the Langevin dynamics approach, one first calculates the electronic friction and other forces via some quantum mechanical method, which, given the complexity of molecular models, is a nontrivial task. Many approaches to calculating electronic friction have been developed, such as using first-principles calculations,^{55} path integrals,^{63,64} scattering theory,^{51} nonequilibrium Green’s functions (NEGFs),^{56} and quantum master equations.^{54,65} In this context, we have developed the HEOM-LD approach, which combines the numerically exact HEOM method with Langevin dynamics (LD).^{66,67} Since the HEOM method allows for numerically exact simulation of systems with strong metal–molecule interactions, non-Markovian effects in the surfaces, and strong electron–electron or electronic–vibrational interactions within the molecule, it is a particularly promising approach for calculating electronic friction.

So far, however, the HEOM-LD approach has incorporated the electronic–vibrational coupling only within the molecule itself, neglecting any metal–molecule coupling that depends on the position of vibrational coordinates, which can strongly influence the vibronic dynamics. Natural examples include scattering^{1,3,4} and desorption^{2,16,17,68} processes at metal surfaces, where the metal–molecule coupling depends sensitively on the distance to the surface. Additionally, in the context of transport through molecular nanojunctions, position-dependent metal–molecule couplings are crucial in understanding vibrational instabilities^{69} and quantum shuttling.^{70,71} We note that there exist other expressions for the electronic friction of systems with position-dependent metal–molecule coupling, both for noninteracting^{53} and interacting^{59,60} systems. Practical calculations for interacting systems, however, often rely on further approximations, such as a mean-field approach. The expression obtained via HEOM, while formally identical to that found, for example, in Ref. 59, has the benefit of being easily computationally accessible, even for interacting systems. For this reason, in this paper, we extend the HEOM-LD approach to systems with a position-dependent metal–molecule coupling.

While our approach is able to describe molecules interacting with metal surfaces in general, we specifically focus on molecular junctions for the examples that we show in this paper later on. By comparing the electronic friction calculated via HEOM to that calculated via NEGFs, which can be used for noninteracting systems, we verify the updated HEOM approach. Following this, we calculate the electronic friction of an interacting vibronic model via HEOM. Here, we demonstrate the usefulness of our extended approach since the HEOM remains numerically exact even for interacting systems.

The paper is structured as follows. In Sec. II, we introduce a general model of a molecule interacting with metal surfaces. In Sec. III, we give a short overview of the HEOM approach. Following this, in Sec. IV, we generalize the HEOM approach to electronic friction to systems with position-dependent metal–molecule couplings. In Sec. V, we apply the approach to different models of molecular nanojunctions.

Since we have to differentiate between classically and quantum mechanically treated coordinates and momenta, we will explicitly denote the vectors of position and momentum operators as $x\u0302$ and $p\u0302$, while their classical counterparts will be denoted by **x** and **p**. Moreover, in this work, we use units where *e* = ℏ = 1.

## II. MODEL

*H*

_{mol}is the Hamiltonian of the molecule,

*H*

_{met}is the Hamiltonian of the metal surface(s), and

*H*

_{met–mol}is the metal–molecule interaction Hamiltonian. In general, the molecular Hamiltonian is given by

*M*

_{i}, respectively. Moreover, $Uvib(x\u0302)$ describes the potential energy surface of a reference electronic state of the molecule, e.g., the electronic ground state of the neutral molecule. When

*m*=

*m*′, $\u03f5mm\u2032(x\u0302)$ represents the energy of the

*m*th electronic level of the molecule, with the coupling between states

*m*and

*m*′ given by

*m*≠

*m*′. The corresponding fermionic creation and annihilation operators for level

*m*are $dm\u2020$ and

*d*

_{m}, respectively. The electron–electron interaction between two electrons in states

*m*and

*m*′ is considered via the matrix $Umm\u2032(x\u0302)$. Although spin has not been written explicitly within

*H*

_{mol}, it can be included with no loss of generality in the later theory. Note that the electronic energies, hopping between electronic states, and electron–electron interactions can all depend on the vibrational coordinates, such that these terms already incorporate the electronic–vibrational interaction.

^{72}

*ϵ*

_{kα}is the energy of state

*k*in lead or surface

*α*, and $ck\alpha \u2020$ and

*c*

_{kα}are the respective creation and annihilation operators. The surfaces are held at local equilibrium, with a well-defined chemical potential,

*μ*

_{α}, and temperature,

*T*. Note that, since we exclusively investigate molecular nanojunctions in this work, we will consider two surfaces, such that

*α*∈ {

*L*,

*R*}, where

*L*and

*R*refer to the left and right leads of the nanojunction, respectively.

*t*

_{kα,m}describes the coupling strength between state

*m*in the molecule and state

*k*in lead

*α*. Here, we explicitly assumed that the position dependency of the metal–molecule coupling affects all states

*k*equally.

*α*

*ρ*

_{T}, factorizes

*ρ*

_{mol}is the initial state of the molecule, and

*ρ*

_{met}(0) is the initial state of the surfaces, which is a tensor product of Gibbs’ states,

*σ*= ± and $\sigma \u0304=\u2213$, $ck\alpha \u2212=ck\alpha $ and $ck\alpha +=ck\alpha \u2020$, and $ck\alpha \sigma (t)=eiHmettck\alpha \sigma e\u2212iHmett$.

*α*. The correlation functions of the metal can be expanded as a sum of exponential functions,

*η*

_{α,σ,ℓ,m}and

*κ*

_{α,σ,ℓ,m}, we decompose the Fermi–Dirac distribution via a sum-over-poles approach based on the AAA method,

^{73}or a Padé decomposition,

^{74}which performs better than a Matsubara decomposition. For a given spectral density, one can then calculate

*η*

_{α,σ,ℓ,m}and

*κ*

_{α,σ,ℓ,m}by applying residue theory to Eq. (10).

*W*

_{α}is the bandwidth. The constant Γ

_{α,mm′}= 2

*πt*

_{α,m}

*t*

_{α,m′}reflects the maximum of the energy-dependent part of the level-width function, which is independent of the individual bath states, and will be referred to from now on as the metal–molecule coupling strength.

## III. QUANTUM TRANSPORT THEORY

*ρ*

_{T}(

*t*) is the total density matrix at time

*t*, including the metal surfaces. However, one is usually only interested in the explicit dynamics of the molecule itself, which is obtained by tracing out the metal DOFs,

There are a variety of methods to circumvent this issue and to describe the influence of the metal on the molecule, some of which have already been mentioned in the introduction. In this work, we use the numerically exact HEOM approach based on the Feynman–Vernon influence functional. Here, we only give a short overview over the HEOM for the sake of comprehensibility of the derivations and calculations in this paper. For a more detailed description and derivation of the HEOM, we refer to Refs. 23, 25, 28, 75, and 76, and the recent review of Tanimura.^{77}

The main idea of the HEOM approach is to describe the dynamics of the molecular density matrix as coupled to so-called auxiliary density operators (ADOs) of different tiers, $\rho j(n)$. In this formalism, the zeroth tier ADO corresponds to the molecular density matrix, *ρ*^{(0)} = *ρ*_{mol}, while higher tier ADOs contain information about metal–molecule interactions and non-Markovian effects in the metal surfaces.^{28} The ADOs are connected by a hierarchy of first-order coupled differential equations.

*n*th-tier ADO is

**j**= [

*j*

_{n}, …,

*j*

_{1}], $jr=(\alpha jr,\sigma jr,\u2113jr,mjr)$, with $j\u0304r=(\alpha jr,\sigma \u0304jr,\u2113jr,mjr)$. Moreover, the superindex

**j**

_{−}= [

*j*

_{n}, …,

*j*

_{r+1},

*j*

_{r−1}, …,

*j*

_{1}] is formed by removing the index

*j*

_{r}, and the superindex

**j**

_{+}= [

*j*,

*j*

_{n}, …,

*j*

_{1}] is formed by adding an additional index

*j*.

## IV. ELECTRONIC FRICTION AND LANGEVIN DYNAMICS

It is often more efficient to use mixed quantum–classical approaches, in which the vibrational DOFs of the molecule are treated classically, as influenced by quantum mechanical electronic DOFs. In this section, we will show how to obtain such a mixed quantum–classical description.

### A. Partial Wigner transform

*O*with respect to

*N*vibrational DOFs,

**x**= (

*x*

_{1}, …,

*x*

_{N}),

**p**= (

*p*

_{1}, …,

*p*

_{N}), and the superscript $W$ indicates that an operator has been Wigner-transformed.

^{78–80}For the sake of better readability later on, we will no longer explicitly write down the

**x**and

**p**dependency of the Wigner-transformed operators except when introducing new quantities for which the dependence is important.

*O*and

*Q*, is

^{81}

^{78,82}

**x**and

**p**are contained only in the molecular DOFs.

### B. Wigner-transformed HEOM

In this section, we apply the partial Wigner transform to Eq. (17), deriving an analog to the quantum–classical Liouville equation in the HEOM formalism. The Wigner-transformed HEOM for the case of a position-independent metal–molecule coupling has already been introduced in Ref. 66, with similar discussions for bosonic reservoirs also in Ref. 77. The focus in this section, therefore, will be on how the inclusion of position-dependent $g\alpha (x\u0302)$ affects the partial Wigner transform with respect to the vibrational DOFs.

### C. Time evolution of the phase-space probability density and Fokker–Planck equation

*A*(

**x**,

**p**,

*t*) in the HEOM formalism can be expressed as

Although we now have an expression for time evolution of the PSPD in terms of the HEOM formalism, it is still in the form of a quantum–classical Liouville equation. In order for electronic friction to appear, we need to further apply a limit of weak nonadiabaticity, which we do by assuming a timescale separation between fast electronic relaxation and slow vibrational dynamics. In the following, we explore the effects of such an assumption.

^{56}

**simply contains the difference between the exact vibronic dynamics and this adiabatic limit.**

*B***provides only a small nonadiabatic correction. By following a similar derivation as in Ref. 66, we show in Appendix C that, under this assumption and by combining Eq. (42) with Eq. (39), the time evolution of**

*B**A*in Eq. (39) is given by the Fokker–Planck equation

*γ*

_{ij}(

**x**), and the correlation function of the stochastic force,

*D*

_{ij}(

**x**). Although Eq. (44) has the same form as was derived in Ref. 66, the inclusion of a position dependence in the metal–molecule coupling adds additional terms to these electronic forces.

*f*

_{i}(

*t*) is a Gaussian random force with white noise

## V. RESULTS

In this section, we apply the extended HEOM-LD formalism introduced in Sec. IV C to several example scenarios, investigating how the inclusion of a position-dependent metal–molecule coupling affects the electronic forces entering the Langevin dynamics. Although the extended HEOM-LD approach is applicable to the general problem of molecules interacting with metal surfaces, we will consider the specific scenario of charge transport through a molecular nanojunction, in which we have two metal surfaces, the left, *α* = L, and right, *α* = R, leads. These can either be held at equilibrium, such that the chemical potentials and temperatures of both leads are equal, or driven out of equilibrium via a voltage bias, Φ = *μ*_{L} − *μ*_{R}, applied symmetrically around the Fermi level, which is always set to zero, *ɛ*_{F} = 0, such that *μ*_{L} = −*μ*_{R} = *e*Φ/2.

The first scenario is nonequilibrium transport in a simple vibronic model, in which the molecule contains a single electronic level linearly coupled to a single harmonic vibrational mode. The second system introduces a strong electronic–vibrational interaction, such that the molecular model contains a single electronic level coupled linearly to both a low-frequency mode, which is treated classically via Langevin dynamics, and a high-frequency mode, which must be treated quantum mechanically alongside the electronic DOFs within the HEOM. For both systems, we will first investigate adiabatic and nonadiabatic effects at equilibrium before discussing the nonequilibrium situation.

### A. Simple vibronic model

*ϵ*

_{0}is the energy of the electronic level without coupling to the vibrational mode,

*λ*is the electronic–vibrational coupling strength, and

*ω*is the frequency of the vibrational mode. We can also rewrite the molecular Hamiltonian as

^{53}Following Ref. 53 and rewriting the electronic friction in our notation, one obtains the following formula at equilibrium in the wideband limit:

*W*

_{α}= 30 eV for all HEOM calculations.

To help identify the specific effect of the position-dependent metal–molecule coupling, we first investigate the electronic forces of a constant coupling, $g\alpha (x\u0302)=1$. Although similar systems have already been analyzed previously, for the sake of self-completeness, we will introduce the main features of the electronic forces, starting in the equilibrium situation, Φ = 0 V.

*x*values, the adiabatic potential follows the diabat of the charged state, since the electronic level at these positions is so far below

*ɛ*

_{F}that it is mostly occupied. In contrast, for large positive

*x*values, the energy level is so far above

*ɛ*

_{F}that the electronic level is mostly unoccupied and the adiabatic potential follows the neutral state diabat.

To investigate the nonadiabatic effects, we turn to Fig. 1(b), which shows the electronic friction for this system calculated via NEGFs as well as the HEOM-LD approach. The friction shows a peak at $x\epsilon F$, as the friction represents a nonadiabatic correction to the adiabatic force, and it is exactly at this position that there is an avoided crossing in *U*_{ad}(*x*). In particular, the friction describes the relaxation of the vibrational mode via coupling to electron–hole pairs (EHPs), in that electrons may transfer from a lead to the molecule, absorb energy from the vibrational mode, and transfer back to a higher energetic state in the lead. At equilibrium, such processes are resonant exactly at the position where the electronic level crosses the chemical potential, $x\epsilon F$.

^{84}

*α*, respectively. The spatial part of the metal–molecule coupling, therefore, takes the form

*g*

_{c}so that

*g*

_{L}(±∞) =

*g*

_{R}(∓∞). Similar metal–molecule couplings have been used in Refs. 16–18 and 68, and recently in Ref. 69. The position dependency of the metal–molecule coupling is shown in Fig. 2. To compare with the case of a constant metal–molecule coupling, we set the metal–molecule-coupling strength to $\Gamma \alpha =2\pi t\alpha 2=0.01eV$, since now for $x\u0302\u2192\xb1\u221e$ the molecule is essentially only coupled to one surface.

Figure 3(a) shows the potential energy surfaces for the vibronic model with the position-dependent metal–molecule coupling from Eq. (61). As in the constant coupling case, *U*_{ad}(*x*) displays an avoided crossing at $x\epsilon F$. However, in the vicinity of *x* = 0, two minima appear, which originate from *g*_{α}(*x*). First, we note that the contribution to the adiabatic force arising from the metal–molecule coupling, $Fi,met\u2212molad(x)$, is proportional to the spatial derivative of the metal–molecule couplings, $\u2202g\alpha \u2202x$.

At *x* = 0, both $\u2202gL\u2202x$ and $\u2202gR\u2202x$ reach their extrema. Since this is close to the point at which the energy level crosses the Fermi level of the leads, any change in the vibrational coordinate at this point is accompanied by a strong change in the electronic occupancy. Therefore, in the vicinity of *x* = 0, the contribution to the adiabatic force from the metal–molecule coupling is maximized. Hence, the classical vibrational DOF will move away from these positions in either a positive or negative direction to positions where the forces acting on it are minimized, which is why there are two minima close to *x* = 0.

This contribution also manifests itself in the electronic friction and correlation function, shown in Figs. 3(b) and 4. First, we note in Fig. 3(b) that the electronic forces calculated from the HEOM and NEGFs methods still match exactly, even with the addition of a position-dependent metal–molecule coupling. We also note that the fluctuation–dissipation theorem is still satisfied because Φ = 0 V. To understand the behavior of these forces, we turn to Fig. 4, which plots the full friction, *γ*(*x*), as well as the contributions from the molecular Hamiltonian, *γ*_{mol}(*x*), the metal–molecule interaction, *γ*_{met–mol}(*x*), and the friction for a constant metal–molecule coupling from the previous model, *γ*_{c}(*x*).

As in the constant coupling case, the full electronic friction also displays a peak at $x\epsilon F$, which originates from the same EHP creation processes that were discussed in Fig. 1(b). However, there is now an additional side peak at a point that we will name *x*_{mm}, which arises directly from the position dependency of the metal–molecule coupling. To understand the origin of this side peak, we consider the individual contributions to the total friction. First, we note that *γ*_{mol}(*x*) also contains a side peak at *x*_{mm} and does not just reproduce the friction in the constant coupling case. Although *γ*_{mol}(*x*) does not depend explicitly on the metal–molecule coupling, it depends on it implicitly via $\u2202\sigma ss(0),W\u2202x$, which is the origin of this side peak. At position *x*_{mm}, the metal–molecule couplings are changing rapidly and the molecular electronic energy is close to the Fermi level, such that $\sigma ss(0),W$ is also changing rapidly. When one includes the *γ*_{met–mol}(*x*) contribution, which has a peak at *x*_{mm}, we recover the total friction. From a physical point of view, electronic friction represents a nonadiabatic vibrational relaxation process via coupling to EHP creation. Such processes are dominant at points where $\u2202\sigma ss(0),W\u2202x$ and $\u2202g\alpha \u2202x$ are large, as these are points in space at which electrons can easily jump back and forth between metal and molecule.

This discussion also helps to connect the HEOM-LD expression for the electronic friction to previous expressions, such as that shown for noninteracting electronic systems from NEGFs. In Eq. (56), for example, the friction comprises four terms: two that contain only $\u2202g\alpha \u2202x$ or $\u2202HmolW\u2202x$, and two cross terms. In the HEOM-LD expression, although it appears there are only two components to the electronic friction, these cross terms are contained within $\u2202\sigma ss(0),W\u2202x$ and $\u2202\sigma ss,j(1),W\u2202x$.

*μ*

_{α}. Consequently, the main equilibrium friction peak originating from $\u2202HmolW\u2202x$ splits into two peaks at positions

*x*

_{mm}, which is now significantly more pronounced than the equilibrium case.

This can be explained considering the specific form of the metal–molecule coupling for this system. At *x*_{mm}, the metal–molecule couplings are equal and the electronic level sits well within the bias window, such that the electronic level is half filled. From Fig. 2, we see that a small positive displacement in *x* reduces the coupling to the left lead while increasing the coupling to the right lead. Given the direction of the voltage bias, the molecule immediately decharges. Similarly, a small negative displacement results in coupling only to the left lead, such that the molecule immediately charges. These two features combine for a new set of EHP creation processes, which manifest as a large $\u2202\sigma ss(0),W\u2202x$.

Finally in this subsection, we note that the fluctuation–dissipation theorem is broken out of equilibrium, such that *γ*(*x*) ≠ *D*(*x*)/*k*_{B}*T*. While friction describes vibrational relaxation via coupling to EHP creation, the stochastic force and, therefore, its correlation function describe vibrational excitation not only via coupling to EHP relaxation but also to the electric current. For this reason, *D*(*x*) displays a broad peak between $x\mu R$ and $x\mu L$, as for these coordinates the electronic level is within the bias window and there is a nonzero electric current. For most of this coordinate range, *D*(*x*)/*k*_{B}*T* > *γ*(*x*), which is a symptom of Joule heating. Around *x*_{mm}, however, the strong peak in the friction is larger than the corresponding peak in the correlation function of the stochastic force, indicating that the vibrational relaxation processes at this point outweigh the heating originating from the electric current.

### B. Interacting vibronic model

*ϵ*

_{0}is the energy of the electronic level,

*ω*

_{cl}and

*ω*

_{qu}are the frequencies of the classical and quantum vibrational modes, respectively, and

*λ*

_{cl}and

*λ*

_{qu}are the corresponding electronic–vibrational couplings.

Such a model can describe a situation where a molecule is coupled not only to the electrodes but also to a high-frequency cavity mode,^{60} or a molecule with multiple modes operating on different timescales coupled to the transport.^{67} For example, one could consider a molecular nanojunction consisting of a large molecule with a small side group attached to it. The center-of-mass motion is slow due to the large effective mass, while the frequency of the side group bond is much higher.

To connect with the noninteracting models considered in Sec. V A, we again choose the metal–molecule coupling given in Eq. (61). Furthermore, in this subsection, we use a bandwidth of *W*_{α} = 10 eV and set the electronic energy, *ϵ*_{0}, such that the energy after small polaron shift is $\u03f5=\u03f50\u2212\lambda qu2\omega qu=0.05eV$, which is the same as in Sec. V A. To ensure that the timescale separation assumption is satisfied, we choose *ω*_{qu} = 300 meV = 10 *ω*_{cl} and *λ*_{qu} = 550 meV. All other parameters are the same as in Subsection V A.

We consider first the equilibrium case, Φ = 0 V. Figure 6(a) shows the potential energy surfaces of the neutral state, *U*_{g}(*x*), charged state, *U*_{c}(*x*), and the adiabatic potential, *U*_{ad}(*x*) for this system. We observe that the addition of a quantum vibrational mode for this system barely influences the adiabatic potential. In contrast, the nonadiabatic forces are strongly influenced by the addition of a quantum vibrational mode, which can be seen in Fig. 6(b), which shows the friction and the correlation function for this system.

We note that the fluctuation–dissipation theorem from Eq. (60) is also fulfilled for the interacting system. Again, we observe the double peak structure of the friction, which can be explained in the same way as for the noninteracting model in Subsection V A. The peak at $x\epsilon F=\u2212(\u03f5\u2212\epsilon F)/\lambda cl$ is at the position where the electronic level crosses the chemical potential. In contrast, the second peak originates from the position dependency of the metal–molecule coupling. Here, at *x*_{c}, the spatial derivatives of the metal–molecule couplings reach their extrema, resulting in the peak in the friction at this position.

We note that the electronic friction of the interacting system in Fig. 6(b) is significantly larger than the electronic friction of the noninteracting system for the entire range of vibrational coordinates [cf. Fig. 3(b)]. Even at equilibrium, the exchange of quanta with the high-frequency mode allows additional lead electrons to transfer to the molecular junction, absorb vibrational energy from the classical mode, and transfer back to a higher energetic state in the leads. Therefore, the addition of a quantum vibrational mode greatly enhances the number of EHP creation processes, resulting in a larger electronic friction for the interacting system compared to the noninteracting system.^{15,66}

However, while for the interacting system the peak at *x*_{c} in Fig. 6(b) has a smaller value than the peak at *x*_{μ}, for the corresponding peaks in the noninteracting system in Fig. 3(b) the opposite is true and can be explained as in the following. The enhancement of the electron–hole pair creation process is especially dominant in the vicinity of *x*_{μ}, where the electronic level crosses the chemical potential. Therefore, it influences the friction at *x*_{μ} more than at *x*_{mm}. Hence, for the interacting system in Fig. 6(b), the peak in the friction at around *x*_{mm} is less pronounced compared to the peak at around *x*_{μ}, in contrast to the noninteracting case in Fig. 3(b).

We now consider the nonequilibrium case and apply a bias voltage Φ = 1 V to the leads. Figure 7 shows the friction and the correlation function for the interacting system. Out of equilibrium, the peak in the friction in Fig. 6(b) at *x*_{μ} splits into two peaks at $x\mu L$ and $x\mu R$. Again, we observe that the peak in the friction close to *x*_{mm} is significantly more pronounced in comparison to the other peaks than in the equilibrium situation. This can be explained in the same way as in Subsection V A, considering the specific form of the metal–molecule coupling for this system.

^{66}

*n*vibrational quanta with the high-frequency mode are resonant, resulting in additional extrema of the friction compared to the noninteracting case. Most interestingly, for this interacting system, such processes do not always result in a cooling of the low-frequency mode, as the electronic friction is negative for a significant portion of the coordinate range, which indicates areas of space in which local heating effects dominate. The appearance of such negative friction could connect to vibrational instabilities within molecular nanojunctions.

^{34–36}

## VI. CONCLUSION

In this work, we extended the theory of electronic friction within the HEOM formalism to systems with a position-dependent metal–molecule coupling. The incorporation of such couplings is critical for a realistic description of many scenarios, such as reaction dynamics at metal surfaces, quantum shuttling, and vibrational instabilities in molecular nanojunctions. Consequently, the updated theory represents one of the most general mixed quantum–classical approaches to the nonadiabatic dynamics of molecules interacting with metal surfaces currently available.

Generally, our theory showed that incorporating a position-dependent metal–molecule coupling adds extra terms to the adiabatic contribution to the mean force, the electronic friction, and the correlation function of the stochastic force. We next demonstrated that such extra terms can critically affect both the adiabatic and nonadiabatic forces. We first applied the updated HEOM-LD approach to a nonequilibrium charge transport through a simple vibronic system containing a harmonic vibrational mode linearly coupled to an electronic level. Since the electronic part of this system is noninteracting, we were also able to compare our updated theory with electronic friction obtained from NEGFs, verifying our extended approach. We next demonstrated the power of the HEOM-LD approach by calculating the electronic forces of an interacting two-mode system, where one of the modes was treated classically and the other quantum mechanically.

In both scenarios, we demonstrated that incorporating a position-dependent metal–molecule coupling significantly affects the electronic forces at positions where the spatial derivatives of the metal–molecule couplings reach their extrema. This was particularly relevant for the electronic friction and stochastic force, which are critical for understanding nonadiabatic processes at metal surfaces. Beyond the reduced models considered here, we foresee great applicability of the updated HEOM-LD approach to more realistic models, as the HEOM treats the quantum mechanical part of the problem in a numerically exact manner, while the Langevin dynamics is capable of handling anharmonic vibrational degrees of freedom with nonlinear electronic–vibrational interactions.

## ACKNOWLEDGMENTS

This work was supported by the German Research Foundation (DFG) through Research Unit FOR5099. S.L.R. acknowledges the Alexander von Humboldt Foundation for support via a research fellowship. Moreover, the authors acknowledge Riley Preston for helpful discussions. 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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Martin Mäck**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (supporting); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **Michael Thoss**: Conceptualization (equal); Formal analysis (supporting); Funding acquisition (lead); Methodology (supporting); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (equal). **Samuel L. Rudge**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Supervision (equal); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: DEFINITION OF WIGNER-TRANSFORMED SUPEROPERATORS IN THE HEOM FORMALISM

### APPENDIX B: TIME EVOLUTION OF THE PHASE-SPACE PROBABILITY DENSITY IN THE HEOM FORMALISM

In the following, we show a detailed derivation of the time evolution of the PSPD from Eq. (39) in terms of the HEOM. We follow a similar procedure as in Ref. 25 and compare the time evolution of the molecular density matrix with the time evolution of the zeroth tier ADO. In contrast to Ref. 25, however, we consider the time evolution of the Wigner-transformed operators.

_{a}brackets on the right hand side of Eq. (B4) and match the terms on both sides of Eq. (B4) depending on $\u2202\u2202pi\rho j(1),W$ and $\u2202\u2202pi\rho TW$,

### APPENDIX C: DERIVATION OF THE FOKKER–PLANCK EQUATION

In the following, we show that by assuming a timescale separation between the dynamics of classical and quantum DOFs, the time evolution of the PSPD in Eq. (39) is given by a Fokker–Planck equation.

*τ*

_{qu}and

*τ*

_{cl}, respectively. We make the same assumption as in Ref. 66, that the dynamics of the classical DOFs are considerably slower than the quantum dynamics, so that

*τ*

_{cl}≫

*τ*

_{qu}. Hence, on a timescale determined by

*τ*

_{cl}, the quantum DOFs will immediately reach their steady state. Moreover, within the timescale determined by

*τ*

_{qu}, the classical vibrational DOFs will barely influence the quantum DOFs. Therefore, we can write the ADOs as in Eq. (42). Recalling the time evolution of the PSPD in terms of the HEOM from Eq. (39) and using Eq. (42), we obtain

**B**,

**B**, obtaining

**B**into Eq. (C3). Thereby, we assume $\u2202\u2202tB\u226a\u2202\u2202tA\sigma ssW$, $\u2202\u2202xiB\u226a\u2202\u2202xiA\sigma ssW$, and $\u2202\u2202piB\u226a\u2202\u2202piA\sigma ssW$, and that, therefore, derivatives of

**B**can be neglected. Under this assumption, it becomes apparent that only four terms contribute to

**B**, so that

**x**-dependency in the metal–molecule coupling, and the last term comes from inserting Eq. (C3) into Eq. (C6) and neglecting derivatives of

**B**. Inserting these contributions in the last line of Eq. (C3), we obtain the Fokker–Planck equation from Eq. (44).

## REFERENCES

*Molecular Electronics*

*Classical Description of Nonadiabatic Quantum Dynamics*

*Advances in Chemical Physics*