Nonadiabatic Dynamics of Molecules Interacting with Metal Surfaces: A Quantum-Classical Approach Based on Langevin Dynamics and the Hierarchical Equations of Motion

A novel mixed quantum-classical approach to simulating nonadiabatic dynamics of molecules at metal surfaces is presented. The method combines the numerically exact hierarchical equations of motion approach for the quantum electronic degrees of freedom with Langevin dynamics for the classical degrees of freedom, namely, low-frequency vibrational modes within the molecule. The approach extends previous mixed quantum-classical methods based on Langevin equations to models containing strong electron-electron or quantum electronic-vibrational interactions, while maintaining a nonperturbative and non-Markovian treatment of the molecule-metal coupling. To demonstrate the approach, nonequilibrium transport observables are calculated for a molecular nanojunction containing strong interactions.


I. INTRODUCTION
The dynamics of molecules interacting with metal surfaces is a highly relevant topic in both physics and chemistry, with many technological applications.It covers a broad range of physical scenarios, including the scattering of molecules off surfaces [1,2], reactive and catalytic processes [3][4][5], and charge transport through STM setups [6,7] or molecular nanojunctions [8][9][10][11].Particularly interesting in these systems is understanding how energy is transferred between the continuum of electronic states in the metal and the molecule, and then dispersed among the molecule's vibrational degrees of freedom.Such electronic and vibrational relaxation processes are often nonadiabatic and require theoretical treatments that go beyond the Born-Oppenheimer approximation.
In Langevin dynamics, one must first calculate the electronic forces, such as the friction, from a suitable quantum method.Although this can be done from first-principles [39], in recent years, a number of new approaches based on quantum transport methods have arisen [40][41][42][43][44][45].In Ref. [46], for example, two of the authors demonstrated how to calculate the electronic friction within the HEOM formalism, which generalizes methods such as nonequilibrium Green's functions [47][48][49][50], scattering theory [51][52][53][54], and path integrals [4,[55][56][57][58] by including strong interactions within the molecule, all while treating the molecule-metal coupling in a nonperturbative and non-Markovian manner.Consequently, the HEOM approach to electronic forces represents a significant step forward in the range of molecular models available for Langevin dynamics.However, in Ref. [46], the authors restricted their investigation exclusively to the forces themselves.
In this work, to our knowledge for the first time, the electronic forces calculated from HEOM are used in the corresponding Langevin equation, introducing the novel HEOM-LD approach to nonadiabatic dynamics of molecules interacting with metal surfaces.First, a brief introduction to HEOM is presented, alongside a summary of how one obtains electronic forces from it.Next, it is shown how to incorporate these forces into the corresponding Langevin equation and how to then solve for the dynamics.Finally, the new approach is demonstrated via several illustrative examples, including molecular models with strong electron-electron and quantum electronic-vibrational interactions.While the HEOM-LD method is completely applicable to transient molecular dynamics, in this work, all results will focus on one of the most challenging problems in context of open quantum systems, namely calculating steady-state nonequilibrium transport observables in molecular nanojunctions.Alongside these mixed quantum-classical results, numerically exact transport calculations from the fully quantum HEOM are also provided, verifying the accuracy of the HEOM-LD method in appropriate parameter regimes.These quantum calculations are performed via a novel implementation of the HEOM designed to treat lowfrequency vibrational modes, which is an extension of the method presented in Ref. [59].arXiv:2402.13161v3[cond-mat.mes-hall]3 May 2024 Overall, the example applications of HEOM-LD presented in this work opens up calculations of molecules interacting with metal surfaces to models that were previously inaccessible.For example, many molecules have multiple vibrational modes participating in the transport with very different frequencies, for which fully quantum treatments can be prohibitively expensive and mixed quantum-classical treatments are inaccurate.An approach such as HEOM-LD, in contrast, allows one to systematically separate those parts of the nanosystem requiring a quantum treatment and those that will submit to a classical one.
The paper is structured as follows.The model describing a molecular nanojunction is introduced in Section II.Section III contains the quantum transport theory used for the numerically exact comparisons, while in Section IV an overview of HEOM-LD is presented.This includes a summary of how one calculates the electronic forces and a description of the numerical algorithm used to propagate the Langevin equation.The results of the direct comparison between quantum HEOM and HEOM-LD are contained in Section V and the conclusions are presented in Section VI.
Throughout the paper, units are used such that ℏ = e = m e = 1.

II. MODEL
In this section, the general model describing nonadiabatic dynamics of molecules interacting with metal surfaces is introduced.While the theory is broadly applicable to any of the scenarios mentioned in the introduction, such as scattering or desorption problems, the specific focus in Sec.V will be on nonequilibrium transport in molecular nanojunctions.With this in mind, the electronic degrees of freedom within the metal will be referred to as the leads, with Hamiltonian H leads .These are coupled to the molecule, H mol , via the molecule-lead interaction, H mol-leads .The total Hamiltonian is The molecular model contains a series of N el spinindependent electronic states linearly coupled to N vib harmonic vibrational modes, with Hamiltonian The mth electronic level is described by its annihila- In transport approaches that treat the vibrational degrees of freedom quantum mechanically, it will prove useful to introduce bosonic annihilation and creation operators, b i and b † i , which are related to the ith vibrational mode via Although the molecular Hamiltonian in Eq.( 56) assumes harmonic oscillators and linear electronicvibrational couplings, the HEOM-LD approach is not restricted in complexity.The theory presented in Sec.IV is completely applicable to molecules containing multiple spin-dependent electronic states, anharmonic vibrational modes with nonlinear couplings, and even nonadiabatic electronic-vibrational interactions.
The left, α = L, and right, α = R, metallic leads are modeled with noninteracting electrons and are assumed to be initially in local equilibrium, defined by temperatures T α and chemical potentials µ α .In the junctiontype setups considered within Sec.V, a voltage bias is applied symmetrically between the left and right leads: where the operators c kα and c † kα annihilate and create an electron in lead α with energy ε kα , respectively.These states are linearly coupled to the mth electronic state in the molecule via the molecule-lead interaction Hamiltonian, which introduces V m,kα as the coupling strength between electronic level m in the molecule and electronic state k in lead α.The molecule-lead coupling is also characterized by the level-width function of lead α, In this work, all results are calculated for a Lorentzian density of states, such that with a bandwidth of W = 10 eV.Throughout the paper, the parameter Γ m,α = 2π|V m,α | 2 will be referred to as the molecule-lead coupling.Here, the V m,α now represent a constant coupling strength between state m in the molecule and lead α, with the dependence on states in the leads, k, being described solely by the Lorentzian density of states.
In the HEOM approach, it is necessary to describe the influence of the leads on the time-evolution of the molecular degrees of freedom.For the noninteracting leads and linear molecule-lead coupling introduced in Eqs.( 5) and ( 6), this influence is entirely described by two-time lead-correlation functions [22,60], In Eq.( 9), the compact notation σ = ± and σ = ∓ with c − kα = c kα and c + kα = c † kα has been introduced, alongside the interaction picture of the leads, c σ kα (t) = e iH leads t c σ kα e −iH leads t .Additionally, it has been implicitly assumed that the initial state of the junction factorizes, ρ total (0) = ρ(0)ρ leads (0), (10) where ρ(0) is the initial state of the molecule and ρ leads (0) is a Gibbs state: Computationally, the lead-correlation functions can be difficult objects to include.In order to obtain a closed set of equations in the HEOM approach, for example, one expands the C σ m,α (t − τ ) as a series of exponential functions, which is generated via a pole decomposition of the Fermi-Dirac function using the AAA algorithm [61][62][63].All results in this work were calculated at T = 300K, such that the decomposition of the bath-correlation functions converged at ℓ max = 8.

III. QUANTUM TRANSPORT THEORY
This section introduces the numerically exact HEOM approach as the quantum transport method underpinning all results in the paper.This includes both the standard formulation, in which molecular vibrations are treated within the molecular density matrix, and the reservoir formulation, in which some vibrational degrees of freedom in the molecule are treated as an extra bath.To finish the section, a discussion on how to calculate relevant transport observables, such as the electric current and vibrational excitation, is presented.

A. Hierarchical Equations of Motion
The HEOM transport method, which derives from the Feynman-Vernon influence functional, incorporates the effect of the leads on the time-evolution of the molecular density matrix by coupling the molecular density matrix, ρ(t), to a series of auxiliary density operators (ADOs), ρ (n) j (t), in a hierarchical fashion.Since the complete derivation is lengthy, and the approach used in this work closely follows previous formulations, this section will just present a brief overview and leave the details to Refs.[21-23, 59, 60, 64-69].
Considering the setup and underlying assumptions introduced in Sec.II, the equation of motion for an nth tier ADO, ρ Each nth tier ADO is described by an n-dimensional vector of super-indices, j = (j n , . . ., j 1 ), where each index j r = {α jr , σ jr , ℓ jr , m jr } can be thought of as a mapping of the essentially continuous interactions between lead and molecule to a series of finite ones.The reduced density matrix of the molecule is simply the 0th-tier ADO, ρ (0) (t) = ρ(t).Within the HEOM, there are also two more relevant super-indices, j − = (j n , . . ., j r+1 , j r−1 , . . ., j 1 ) and j + = (j, j n , . . ., j 1 ), which are formed by removing and adding a molecule-lead interaction from and to an ADO, respectively.They relate directly to the superoperators C jr and A σ α,m , which couple different tiers of the hierarchy, Eq.( 14) has also introduced the super-index jr = {α jr , σjr , ℓ jr , m jr }.Within the HEOM approach, one needs to converge the dynamics with respect to both the fermionic tier, which roughly corresponds to how many higher-order interactions with the lead are relevant to the transport, and with respect to the size of the vibrational basis.All results in this work were converged at a fermionic tier of either 2 or 3 and a vibrational basis size of 15 to 25.

B. Reservoir Treatment of the Vibrational Degrees of Freedom
The HEOM formulation in Eq.( 13) treats all degrees of freedom within the molecular Hamiltonian on an equal footing within ρ(t).In regimes of vibrational instability, where the vibrational excitation is large, it can be more numerically efficient to treat the corresponding vibrational degrees of freedom as an extra reservoir.This idea has already been explored in Ref. [59,[70][71][72][73], albeit under the assumption that all vibrational degrees of freedom would be treated within the reservoir formulation.In scenarios where some vibrational degrees of freedom exhibit instability while some do not, perhaps due to different electronic-vibrational couplings or frequencies, it would be useful to have a framework that allows stronglycoupled vibrational modes to be kept within the reduced density matrix and weakly-coupled modes to receive the reservoir treatment.
This essentially amounts to rewriting the junction Hamiltonian as where N sys and N res are the number of modes treated in the system and reservoir, respectively.Note that Eqs. ( 17) and ( 18) exclude the zero point energies of the harmonic oscillators, which, since they just add a constant to the Hamiltonian, do not affect the dynamics.
Because the vibrational modes receiving the reservoir treatment are harmonic and the corresponding electronic-vibrational coupling is linear in xk , the effect of the kth vibrational mode in the reservoir on the molecule is also described by a two-time correlation function, where Application of the bath interaction picture now means that vibrational operators treated in the reservoir are time- In order to treat these modes within the HEOM framework, the reservoir vibrational initial state must be chosen so that they are in thermal equilibrium with some temperature, T vib , and that the total initial state still factorizes as in Eq. (10).Under these assumptions, the bath-correlation functions describing the vibrational reservoir decompose analytically into where For steady-state observables, the specific choice of T vib is not particularly important, although it is highly relevant for transient dynamics.
The corresponding HEOM now has two hierarchies arising from the original fermionic electrodes and the new bosonic baths, with the new vibrational superoperators The ADOs of qth bosonic tier and nth fermionic tier now contain not just information about interactions with the metallic leads, but also qth-order interactions with vibrational modes in the reservoir.It is important to note that the truncation tier, q, and super-indices, q, do not refer to specific operators.Rather, each qthtier ADO is described by a vector of super-indices, q = (q − 0 , q + 0 , q − 1 , . . ., q + Nres ), where each index q v k is an integer whose value refers to the number of b v k contained within this ADO.Alternatively, one can think of q v k as the number of interactions of type v with mode k that this ADO describes.The bosonic tier, q, is then calculated by summing over all of these interactions: q = k (q + k + q − k ).For a more detailed formulation, we refer to Ref. [59], specifically Eq. ( 24).
This updated reservoir formulation of the HEOM was implemented via an extension of the solver introduced in Ref. [74], which uses an efficient iterative scheme to solve the HEOM directly for the steady state.Similar to the fermionic hierarchy, one needs to converge the new bosonic hierarchy with respect to the tier, q.All results in this work were converged with a bosonic tier between 10 and 20.

C. Transport Observables from the Quantum HEOM Approach
The ADOs are not just mathematical objects required for the exact time-evolution of ρ(t).Rather, they also provide direct information on bath properties, such as heat and charge flow.The average electric current through electrode α, for example, can be directly obtained from ADOs that are 1st-tier in the fermionic hierarchy [23].When all vibrational degrees of freedom are treated within the molecular density matrix, as in Sec.III A, this yields In the reservoir treatment introduced in Eq.( 23), the current is 1st-tier in the fermionic hierarchy and 0th-tier in the bosonic hierarchy [59], The steady-state average vibrational excitation of mode i, , on the other hand, is calculated in a very different manner between the two treatments.When mode i is treated as part of the molecular Hamiltonian, the excitation is When mode i is treated as a bosonic reservoir, b † i b i is now an observable 2nd order in bath operators, so it can be reconstructed from ADOs 2nd-tier in the bosonic hierarchy [59], where Ni (0) is the initial excitation due to choosing a thermalized initial state.

IV. ELECTRONIC FRICTION AND LANGEVIN DYNAMICS: THE HEOM-LD APPROACH
In Ref. [46], two of the authors formulated a mixed quantum-classical version of the HEOM, which, for selfcompleteness, is outlined here.From it, one can calculate quantum electronic forces acting on classical vibrational degrees of freedom.This section also describes how one uses these electronic forces in the corresponding Langevin dynamics to form the complete HEOM-LD method.

A. Derivation of the HEOM-LD approach
Often, vibrational modes in the molecule can be separated into N qu high-frequency modes that need to be treated quantum mechanically, ( xqu , pqu ), and N cl lowfrequency modes that can be well-approximated with a classical treatment, (x cl , p cl ).As discussed in Ref. [46], this amounts to Wigner transforming the HEOM with respect to the low-frequency vibrational degrees of freedom and truncating the resulting expansion for terms linear in ℏ and higher.In the joint Liouville space of all ADOs, this yields where Õ is the Wigner transform of operator Ô and . The leading-order contribution in the expansion of the Wigner-transformed HEOM corresponds to the first term in Eq. (30), which is proportional to 1/ℏ and contains the HEOM dynamics frozen at one classical vibrational frame, L(x cl ).The equation of motion for the nth-tier ADO described by this superoperator is Here, the ADOs are written with σ and not ρ so as to distinguish between the objects following the time-evolution in Eq.( 31) and the objects following the time-evolution in Eq. (30).
The second term in Eq.( 30), meanwhile, arises from the next-to-leading-order term in the expansion of the Wigner-transformed HEOM, which is 0th-order in ℏ.Similar to a quantum-classical Liouville equation, this contains a symmetrized Poisson bracket with each ADO [45,75], where A Langevin equation of motion for the classical vibrational modes can then be derived by directly applying a timescale separation between quantum and classical degrees of freedom to the equation of motion for A(x cl , p cl ; t) = Tr mol,qu [ρ(x cl , p cl ; t)], which is the phasespace probability density of the classical vibrational modes.Note that the Tr mol,qu [. . .] is now a trace over the quantum degrees of freedom remaining in the molecular Hamiltonian after the Wigner transform and that these degrees of freedom may include any high-frequency vibrational modes requiring a quantum mechanical treatment.One expects this to be justified if the timescale of electronic relaxation is much faster than the timescale of the vibrational dynamics.In scenarios where the moleculelead coupling is larger than the effective vibrational frequency, Γ, k B T ≫ ω, such as the dynamics of heavy molecules, the HEOM-LD should perform well.At high temperatures, furthermore, not only are quantization effects negligible, which justifies a classical treatment of the vibrations, but the Markovian approximation of electronic dynamics is also valid.
After these two approximations, the classical vibrational modes in the molecule follow a Markovian Langevin equation of motion, which for the ith coordinate is where f i (t) is a Gaussian random force with white noise, Here, D ij (x cl ) is the strength of the correlation function of the stochastic force between vibrational coordinates i and j.
The adiabatic contribution to the mean force, F ad i (x cl ), as well as the force difference operator, δ F ad i (x cl ), are Here, ∂ i is shorthand for a derivative with respect to the ith classical vibrational coordinate, ∂ i = ∂/∂x cl,i , and σss (x cl ) is the steady-state joint density operator of the HEOM frozen at one classical vibrational frame, Eq.( 34) also introduces the electronic friction tensor as well as the correlation function of the stochastic force, At equilibrium, the friction tensor is always symmetric and positive-definite, such that it has an overall damping effect on the classical vibrational motion [46,76].It is also related to the correlation function of the stochastic force via the fluctuation-dissipation relation, Out of equilibrium, however, the fluctuation-dissipation theorem and the positive-definiteness of the friction tensor are not guaranteed.At finite bias voltages, for example, molecular vibrations still have excitation and relaxation mechanism from coupling to electron-hole pairs near the chemical potentials of the leads, but are also excited by the nonequilibrium electric current, which manifests as Although the theory here focuses on electronic friction, the calculation of friction and diffusion tensors is also well established for classical molecular simulations [77][78][79][80].
To end this subsection, note that one could also solve the corresponding Fokker-Planck equation for decomposition of the propagator defined by the Fokker-Planck equation [81].Restricting Eq.( 43) to one dimension for notational simplicity, the time-evolution operator of the classical phase-space probability density, L cl , can be split as The classical time-evolution is defined by the propagator e −∆tL cl , but one can construct a family of approximations via the Trotter decomposition [82], where j ∈ {A, B, O}.Since the propagators e −(∆t/2)Lj do not commute, different orderings within Eq.( 45) yield different approximations.In generating the results in Sec.V, both the ABOBA [83] and OBABO [82] algorithms were tested, where the letters denote the order within Eq.( 45).Drawing from our experience with numerical integration of Markovian Langevin equations, we know that, for large-dimensional systems, the OBABO algorithm has been shown to be generally more accurate because it splits the stochastic part of the propagation [81,84].In all the one-dimensional examples discussed in Sec.V, however, there was no discernible difference in accuracy or convergence speed between the two algorithms.

C. Transport Observables from the HEOM-LD Approach
Consider a general operator treated within the quantum part of the system in the HEOM-LD approach, O, which does not explicitly depend on classical coordinates.At a particular phase-space coordinate, (x cl , p cl ), it will have an instantaneous average, Here, ρtotal is the Wigner transform of the total density matrix with respect to the classical vibrational degrees of freedom and the subscript qu indicates an average over only the quantum part of the system.One calculates the actual expectation value, ⟨O⟩(t), by averaging ⟨O⟩ qu over the classical part of the system, In this work, the instantaneous quantum averages will be approximated adiabatically, such that ρtotal (x cl , p cl ; t) is replaced by σ ss (x cl ), the adiabatic steady-state density matrix frozen at vibrational coordinate x cl : see Eq. (38).For example, the instantaneous average current and electronic occupations are approximated as which is now only a function of x cl .If the quantum part of the system contains vibrational mode i, then this same approach would yield the following approximation for the instantaneous vibrational excitation, In the following section, only steady-state quantities will be discussed.These will be denoted with an ss superscript and can be obtained from the above approach by taking the long-time limit of the phase-space probability density.Since in this work Langevin and not Fokker-Planck dynamics are used, the classical averages are calculated by averaging over many trajectories, rather than directly evaluating Eq.(47).To generate the trajectories and sample points, the following algorithm was used.
First, the electronic forces and instantaneous steadystate expectation values of quantum operators are calculated via the procedure outlined in Sec.IV A. A large number of Langevin trajectories, N traj , are initialized.The initial coordinates and momenta for each trajectory are sampled from the Wigner function of the initial quantum state, which provides the initial distribution for the classical calculation [38].Since all systems considered in this work contain harmonic oscillators and are onedimensional in the classical coordinates, this is always chosen to be the ground state of the harmonic oscillator of the uncharged electronic state, which in dimensionless coordinates is Because the Fokker-Planck equation for the systems treated here has a unique steady-state, the exact choice of ρ W (x cl , p cl ; 0) is, in this work, relatively unimportant.However, if one were interested in transient dynamics, then the choice of ρ(0) would be crucial.Each trajectory is propagated to a long-enough time, t ss , such that A(x cl , p cl ; t ss ) → A ss (x cl , p cl ).In practice, this is converged by first running a low number of trajectories and observing the time-evolution of the average kinetic and potential energies.Once a trajectory reaches t ss , it is propagated further and sampled N sample times, with some interval time, t int , between each sample point to ensure independence.
Finally, these N total = N traj × N sample points are used to calculate the steady-state expectation values.For quantum operators, this is Expectation values of classical observables are calculated in the same manner.Of particular interest is the steadystate excitation of the ith classical vibrational mode.For a harmonic oscillator, this can be approximated from E i , the total energy contained in the mode: In the quantum case, if one were to define the vibrational excitation via the energy, then there would be a contribution from the zero point energy, which is something that the classical approach cannot reproduce.However, as will become evident in the next section, the contribution from the zero point energy is generally small when the vibrational frequency is lower than the temperature and molecule-lead coupling: that is, in the parameter regimes for which HEOM-LD performs well.Finally, note that all results presented in Sec.V have been calculated with N traj = N sample = 500.

V. RESULTS
In this section, the steady-state transport properties of three reduced models for molecules in molecular junctions are explored.The first system, discussed in Sec.V A, contains a single electronic level linearly coupled to a harmonic oscillator, which will be named the onelevel, one-mode (1L1M) model and will be used to explore appropriate parameter regimes in the HEOM-LD approach.The second system extends this simple model by considering two electronic levels that are linearly coupled to a harmonic oscillator and have a Coulomb repulsion between them.This will be named the two-level, one-mode (2L1M) model and is analyzed in Sec.V B. In Sec.V C, a single electronic level is linearly coupled to low-frequency oscillator, which is treated classically, and a high-frequency harmonic oscillator, which is treated quantum mechanically within the HEOM.

A. Basic Vibronic Model
In the mass-and frequency-transformed coordinates, the molecular Hamiltonian is Given that there is only one electronic and vibrational degree of freedom and there are thus only a few parameters, it is an excellent testbed in which to assess the quality of the HEOM-LD approach.This is shown in Fig. (1), where the steady-state electric current and vibrational excitation are plotted as a function of bias voltage for a range of vibrational parameters.Before starting the analysis, note that in Fig. (1), the electronic energy level has been chosen such that the energy level after the small polaron shift is always ε 0 = ε0 − λ 2 ω = 150 meV.This ensures that the onset of resonant transport occurs at the same voltage for all parameter regimes [85].The electronic current displays a single step at the opening of the elastic channel, when ε enters the bias window.Since k B T ≫ ω, the vibrational mode behaves essentially classically and its excitation increases monotonically.Due to the small frequency and electronic-vibrational couplings, the mode experiences a vibrational instability at relatively low voltages, which arises from the well-known effect of current-induced heating [59].Since Γ > ω, the timescale of vibrational motion is too slow for electrons to see anything but fixed nuclei during relaxation and the HEOM-LD approach exactly reproduces the steady-state observables calculated from the full quantum approach for all electronic-vibrational couplings.
In the intermediate regime, shown in Fig. (1c) and Fig. (1d), the vibrational frequency is on the same order of magnitude as the temperature and molecule-lead coupling: ω = 30 meV.Although the timescale separation is formally violated, the HEOM-LD approach still performs quite well.This can be attributed to two factors.First, since k B T ∼ ω, the vibrational dynamics are still essentially classical, such that there is again only one step in the current at eΦ/2 = ε 0 .Second, whether a vibrational mode experiences adiabatic or nonadiabatic dynamics depends not just on the relative timescales of electronic and vibrational motion, but also strongly on the electronicvibrational coupling.In Holstein-type models, this manifests via the separation between the ground and excited state diabats: the small polaron shift.Even in this intermediate regime, the largest shift is λ 2 /ω = 30 meV ∼ Γ, which ensures essentially adiabatic dynamics [86].Finally, note that the mode also experiences a vibrational instability in this regime, although this occurs at higher voltages than in the classical regime.
In the nonadiabatic regime, where Γ, k B T < ω, quantization of the vibrational energy becomes relevant to transport processes.Since the electronic transport occurs on a timescale slower than vibrational motion, the vibrational motion has sufficient time to respond during electronic relaxation, such that electrons occupying the molecule form polarons with vibrational phonons.Under a small polaron transformation, it can be shown that this rescales the molecule-lead coupling to a much smaller effective value, Γ eff ∼ Γe −(λ/ω) 2 [85].As shown in Figs.
(1e) and (1f), where ω = 300 meV, this translates to additional steps in the average current and vibrational excitation, corresponding to the opening of phonon-assisted transport channels entering the bias window, ε + νω with ν = 0, 1, 2, . . .[87,88].For λ/ω = 1, the current and vibrational excitation are suppressed at low bias voltages due to the Franck-Condon blockade [9,89,90].While the mixed quantum-classical approach roughly predicts the quantitative behavior of the electronic current, especially at high bias voltages, it fails to reproduce the extra steps, which are a purely quantum effect.
In comparison to the parameter regimes explored above, many semiclassical treatments restrict ω to the meV range to ensure that the molecule-lead coupling or the temperature is much larger than the effective frequency [33,34,44,51,55,91].Considering that vibrational modes of actual molecules often have effective frequencies in the range of 0.01 − 0.4eV [92] and that realistic molecule-lead couplings also fall within this range, the results of this section demonstrate that the HEOM-LD approach can potentially be used in a wider parameter regime than what has been previously reported, at least for models in which the quantum part of the system is noninteracting.

B. Influence of Electron-Electron Interaction
In this subsection, the simple vibronic model is extended by adding an extra electronic level and an electron-electron interaction, forming the two-level, onemode (2L1M) model, with Hamiltonian In all calculations, both energy levels are coupled equally to both leads: Γ α,m = Γ/2.In addition, parameters are chosen such that the energies after the small polaron shift are The inclusion of the extra electronic level and electronelectron interaction has several effects on the steady-state transport properties, which are shown in Fig. (2).First, the electric current now displays four steps at the voltages where the ε 1 , ε 2 , ε 1 + U , and ε 2 + U transport channels open.In purely electronic systems without a vibrational mode, the current plateaus between steps, but in Figs.
(2a) and (2c), the current increases gradually with the voltage.This is due to the interaction with the vibrational mode, which, since it has such a low ω and λ, can exchange an essentially continuous amount of energy with the electronic degrees of freedom.
These steps are also visible in Figs.(2b) and (2d), which contain the corresponding steady-state vibrational excitation.At low voltages, only the first electronic transport channel is open and ⟨N ⟩ ss increases with voltage in a similar manner to the 1L1M model shown in Fig. (1d).As the second transport channel opens, at eΦ = 2ε 2 = 1.2 eV, the vibrational excitation starts decreasing with increasing voltage.Such local cooling arises from a competition between the already open transport channel and the channel that is currently opening, and has been reported for a similar system in Ref. [87].When only the first transport channel is open, increasing the bias voltage increases the excess electronic energy available to be transferred to the vibrational mode, hence the current-induced excitation.As the second transport channel opens, high-energy electrons can also transport almost adiabatically through ε 2 , exchanging only a small amount of energy with the vibrational mode, hence the decrease in vibrational excitation.This process continues as the four transport channels successively start to open and then saturate.There is a particularly large peak in ⟨N ⟩ ss between the opening of ε 2 and ε 1 + U because the strong Coulomb repulsion induces a large energetic gap between these charging energies.If one were to increase the voltage beyond Φ = 5 V, then ⟨N ⟩ ss would increase monotonically, as all transport channels are fully open.The electron-electron interaction, therefore, has an overall stabilizing effect on the nanojunction in comparison to the one-level, one-mode model.
Finally, note that in this interacting system, the size of the molecule-lead coupling in comparison to the vibrational frequency is crucial for the accuracy of the semiclassical HEOM-LD approach.In the top row, where Γ < ω, the timescale separation assumption is violated, and the HEOM-LD approach clearly does not approximate the steady-state observables as well as it does in the bottom row, where Γ > ω.This contrasts the steady-state observables of the one-level, one-mode system, where the HEOM-LD approach performed well even outside the timescale separation assumption.The key difference is that the electron-electron interaction slows electronic relaxation, such that the molecule-lead coupling must be larger to maintain a Markovian approximation of the electronic forces.

C. Influence of a Quantum Electronic-Vibrational Interaction
In a previous publication, two of the authors investigated the electronic friction of a system containing a single electronic level coupled to both a low-frequency nuclear mode, which could be treated in the semiclassical HEOM-LD framework, and a high-frequency nuclear mode, which must be treated quantum mechanically [46].For brevity, this model will be referred to as the one-level, two-mode (1L2M) in the following.The motivation for such a model could be a molecule with some high-frequency internal vibration accompanied by low-frequency center-of-mass motion, or coupling of the molecule to a cavity mode [50].The purpose of this subsection is to explore how the addition of the quantum high-frequency mode affects the dynamics of the classical low-frequency mode.
The molecular Hamiltonian is where the "lf" and "hf" indices represent low-and highfrequency, respectively.In what follows, only the lowfrequency mode will be treated classically, such that (p lf , xlf ) → (p lf , x lf ), while the high-frequency mode will be treated quantum mechanically alongside the electrons.In a similar manner, two modes operating on such different timescales require a unique treatment in the fully quantum HEOM, in that the low-frequency mode will be treated as an extra reservoir while the high-frequency mode remains within ρ(t): this is the splitting outlined in Eqs.( 17)- (19).Note that, to satisfy the requirement of fast relaxation within the quantum part of the system, the electronic-vibrational coupling of the high-frequency mode is set quite large: λ hf = 1.5ω hf = 450 meV.
In order to investigate how the addition of the highfrequency mode affects the stability of the low-frequency  3), the electronic forces acting on the lowfrequency mode and adiabatic quantum transport observables are plotted for the 1L2M model at one example voltage.Shown alongside are also the electronic forces and adiabatic transport observables for the 1L1M model with the same parameters, just without the highfrequency mode.As reported in Ref. [46], the electronic friction acting on the low-frequency mode is generally much stronger than in the one-mode system, as there are now extra electron-hole pair creation and annihilation processes via the exchange of energy with the highfrequency mode.These extra processes are also the origin for the extra peaks in the friction for the 1L2M model at ε(x lf ) + νω hf = eΦ/2, where ε(x lf ) = ε − λ lf √ 2 x lf and ν ∈ N. Since the electronic forces are plotted for a finite bias voltage, the fluctuation-dissipation theorem is not satisfied and the correlation function of the stochastic force displays a broad peak when the level sits within the bias window while the friction does not, which is a symptom of current-induced heating.
In addition to the extra structure in the electronic friction and correlation function, the instantaneous steadystate electric current also exhibits typical features of Franck-Condon blockade.As shown in Fig. (3b), instead of a single plateau when ε(x lf ) lies within the bias window, the current exhibits steps at the same voltages as peaks in the friction and plateaus inbetween.Likewise, the instantaneous steady-state excitation of the high-frequency mode, ⟨N ⟩ qu (x lf ), displays similar steplike behavior.As the classical coordinate decreases and the level is pushed below the bias window, the level becomes totally occupied and ⟨N ⟩ qu (x lf ) reaches some saturation value.Evidently, introducing a high-frequency mode to the problem significantly changes the electronic forces and adiabatic transport observables.A natural question, therefore, is how the additional structure shown in First, note that ⟨I L ⟩ ss displays steps at the emergence of each phonon-assisted transport channel, eΦ/2 = νω hf .This is quantitatively and qualitatively similar behavior to what a single-mode model with ω hf and λ hf predicts, which is shown by the solid green line.This is due to the relative strengths of the electronic-vibrational interactions, which differ by an order of magnitude between the high-and low-frequency modes.Unlike the single-mode model, however, ⟨I L ⟩ ss does not plateau between steps, instead increasing monotonically due to the essentially continuous exchange of energy with the low-frequency mode.On the other hand, the single-mode model containing ω lf and λ lf predicts only one step in the current and saturates to a much higher value, as the suppression of the current due to the Franck-Condon blockade is not present.
A similar difference between the two models can be seen in Fig. (4b), in which the excitation of the highfrequency mode is plotted against the voltage.In comparison to the single-mode case, the 1L2M model does not plateau between steps of ⟨N hf ⟩ ss .At high bias voltages, the excitation is lower for the 1L2M model, as there is now competition for the energy of transporting electrons with the low-frequency mode.Since λ hf ≫ λ lf , the excitation is reduced by only a small amount.This is in stark contrast to the excitation of the low-frequency mode, which is shown in Fig. (4c).
At low voltages, Φ < 0.6 V, electrons do not have enough incident energy to interact efficiently with the high-frequency mode, and the resulting excitation of the low-frequency mode, ⟨N ⟩ ss lf , follows the same monotonic increase that the corresponding 1L1M model with ω lf and λ lf displays.After some threshold voltage, ⟨N ⟩ ss lf exhibits negative differential vibrational excitation, with maxima just before the phonon-assisted transport channel at eΦ/2 = νω hf becomes available, and minima inbetween.When these channels are fully open, energy dissipated in the junction is distributed unevenly between the two vibrations due to the large difference in λ lf and λ hf , which strongly favors excitation of the high-frequency mode.Adding a strong quantum electronic-vibrational interaction, therefore, has a similar local cooling effect to that already shown in Sec.V B.
As a final comment, one should be careful when using the semiclassical approach with additional quantum modes, as the timescale of electronic motion can be strongly affected by the quantum electronic-vibrational coupling.A strong λ hf , for example, leads to the formation of a polaron and an effective molecule-lead coupling of Γ eff = Γe −(λ/ω) 2 that is orders of magnitude smaller than ω lf .At room temperature, the HEOM-LD approach is still reasonably accurate in comparison to the full quantum HEOM, and it can always be improved by increasing the molecule-lead coupling.Increasing the molecule-lead coupling in the quantum HEOM reservoir approach, on the other hand, is much more difficult, as one needs more electronic tiers within the hierarchy and calculations for the 1L2M model become prohibitively expensive.
This speaks to the overall usefulness of the HEOM-LD approach, in that it works best in regimes where it is most difficult to apply the quantum HEOM.Another key restriction to the quantum HEOM approach is that the vibrational degrees of freedom have been modeled as harmonic oscillators with linear electronicvibrational couplings, which enables the highly efficient reservoir treatment to be used for small λ.Since realistic molecular models include anharmonic modes and nonlinear electronic-vibrational couplings, for which such an approach cannot be used, regimes of high vibrational excitation become very difficult to model.In a discrete variable representation (DVR), for example, large amplitude motion or the motion of heavy nuclei requires many grid points and, correspondingly, a large basis of vibrational states in ρ(t) [93].In these scenarios, mixed quantum-classical approaches are significantly more efficient than their quantum counterparts.The HEOM-LD approach, furthermore, performs well in these regimes even when the quantum part of the nanosystem contains strong interactions.

VI. CONCLUSION
This work introduced the HEOM-LD approach as a novel mixed quantum-classical method for simulating the nonadiabatic dynamics of molecules under the influence of quantum electronic forces from one or more metal surfaces.The method combines the numerically exact HEOM approach for electronic degrees of freedom with a Langevin equation for the vibrational degrees of freedom within the molecule.Within the new approach, strong interactions within the quantum part of the system can be treated while simultaneously incorporating nonadiabatic effects in the vibrational dynamics, such as relaxation due to electron-hole pair creation.
To demonstrate the efficacy of HEOM-LD, steadystate transport observables for a molecular nanojunction model were presented, alongside exact calculations from the fully quantum HEOM approach.For a simple vibronic model, where the electronic part of the system is noninteracting, HEOM-LD performs excellently even when the timescale separation is formally violated.Only in regimes where vibrational quantization becomes important, k B T ≪ ω, are there differences between the quantum and HEOM-LD methods.Under the influence of an electron-electron or quantum electronic-vibrational interaction, however, the electronic relaxation is slowed and fulfilling the timescale separation assumption is crucial for accurate dynamics.Such strong interactions also introduce a stabilizing effect to the low-frequency mode, which manifests via negative differential vibrational excitation.
These calculations highlight the strong potential of HEOM-LD as a method for modeling the nonadiabatic dynamics of molecules interacting with metal surfaces, with a particular focus on systems containing strong interactions.Possible applications include investigating the influence of the electron-electron interaction on the desorption or scattering of a molecule from a metal surface, light-induced processes, and vibrational instabilities in molecular nanojunctions arising from nonconservative electronic forces.

B.
FIG. 1: Steady-state transport observables as a function of bias voltage for the 1L1M model.Left column: Steadystate electronic current, ⟨I L ⟩ ss .Right column: Corresponding steady-state vibrational excitation, ⟨N ⟩ ss .The top row shows results for a vibrational frequency of ω = 3 meV, the middle a vibrational frequency of ω = 30 meV, and the bottom a vibrational frequency of ω = 300 meV.For each vibrational frequency, results for a range of electronic-vibrational couplings are also plotted, which are indicated by color in the plots.Other parameters are Γ = 20 meV and k B T = 25.8 meV.The electronic energy is chosen such that the energy level after the small polaron shift is always ε 0 = ε0 − λ 2 ω = 150 meV.

FIG. 3 : 2 lf ω lf − λ 2 hfω
FIG. 3: Electronic friction and correlation function of the stochastic force (a), instantaneous steady-state electronic current in the adiabatic approximation (b), and vibrational excitation of the high-frequency mode (c), as a function of the low-frequency vibrational coordinate and at finite bias voltage: Φ = 1.5 V. Solid lines correspond to the 1L2M model and dashed lines correspond to the 1L1M model evaluated at the same parameters but without the high-frequency mode.Parameters are Γ = 20 meV, k B T = 25.8 meV, ω lf = 30 meV, λ lf = 10 meV, ω hf = 300 meV, λ hf = 450 meV.Again, the electronic energy level is chosen such that the energy after the small polaron shift(s) is ε 0 = 150 meV.
Fig. (3) manifests itself in the steady state transport observables, which are shown in Fig. (4).

FIG. 4 :
FIG. 4: Blue: steady state transport observables for the 1L2M model as a function of bias voltage.Electric current is shown in (a), excitation of the high-frequency mode in (b), and excitation of the low-frequency mode in (c).Results calculated from the fully quantum HEOM are given solid lines, while those calculated from the HEOM-LD approach are given dashed lines.Also plotted are the corresponding transport observables for the corresponding single-mode model with just the low-frequency (red) and just the high-frequency (green) mode.All parameters are the same as in Fig. (3) U mm ′ , between electrons.The vibrational degrees of freedom, meanwhile, are described by coordinates x = {x 1 , . . ., xN vib } and momenta p = {p 1 , . . ., pN vib }.These degrees of freedom can represent center-of-mass motion of atomic nuclei or bond vibrations.Vibrational operators have been written with explicit operator notation to distinguish between the quantum and classical treatments.Under a transformation to dimensionless coordinates, xi → xi √ m i ω i and pi → pi / √ m i ω i , the molecular Hamiltonian simplifies to tion and creation operators, d m and d † m , with some Coulomb repulsion,