Nonadiabatic dynamics of molecules interacting with metal surfaces: Extending the hierarchical equations of motion and Langevin dynamics approach to position-dependent metal–molecule couplings

Mäck, Martin; Thoss, Michael; Rudge, Samuel L.

doi:10.1063/5.0222076

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,⁵⁵ path integrals,^63,64 scattering theory,⁵¹ nonequilibrium Green’s functions (NEGFs),⁵⁶ 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⁶⁹ 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⁵³ 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 $\hat{x}$ and $\hat{p}$ ⁠, while their classical counterparts will be denoted by x and p. Moreover, in this work, we use units where e = ℏ = 1.

II. MODEL

In this section, we introduce the general model considered within this work. The total setup comprises a molecule that is attached to one or more metal leads or surfaces. The entire setup is therefore described by the Hamiltonian

H = H_{mol} + H_{met} + H_{met - mol},

(1)

where 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

\begin{align} H_{mol} = & \sum_{m m'} ϵ_{m m'} (\hat{x}) d_{m}^{†} d_{m'} + \sum_{m, m' < m} U_{m m^{'}} (\hat{x}) d_{m}^{†} d_{m} d_{m'}^{†} d_{m^{'}} \\ + \sum_{i} \frac{{\hat{p}}_{i}^{2}}{2 M_{i}} + U_{vib} (\hat{x}) . \end{align}

(2)

Here,

\hat{x} = ({\hat{x}}_{1}, {\hat{x}}_{2}, \dots)

and

\hat{p} = ({\hat{p}}_{1}, {\hat{p}}_{2}, \dots)

are vectors of the position and momentum operators of the vibrational coordinates with corresponding mass M_i, respectively. Moreover,

U_{vib} (\hat{x})

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′,

ϵ_{m m^{'}} (\hat{x})

represents the energy of the mth 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

d_{m}^{†}

and d_m, respectively. The electron–electron interaction between two electrons in states m and m′ is considered via the matrix

U_{m m^{'}} (\hat{x})

⁠. 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.⁷²

The metal surfaces or leads are modeled as reservoirs of noninteracting electrons,

H_{met} = \sum_{α} H_{met, α},

(3)

with

H_{met, α} = \sum_{k} ϵ_{k α} c_{k α}^{†} c_{k α},

(4)

where ϵ_kα is the energy of state k in lead or surface α, and

c_{k α}^{†}

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.

The interaction between the molecule and the metal is described by the metal–molecule interaction Hamiltonian

H_{met - mol} = \sum_{k, α, m} g_{α} (\hat{x}) t_{k α, m} (c_{k α}^{†} d_{m} + d_{m}^{†} c_{k α}),

(5)

where

g_{α} (\hat{x})

is a real-valued function describing the spatial dependency of the metal–molecule coupling, and 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.

We note that, unlike previous calculations of electronic friction using HEOM, the metal–molecule coupling depends explicitly on the vibrational coordinates

\hat{x}

⁠. This means that the level-width function of surface α

Γ_{α, m m^{'}} (ϵ, \hat{x}) = 2 π \sum_{k} g_{α}^{2} (\hat{x}) t_{k α, m} t_{k α, m^{'}} δ (ϵ - ϵ_{k α}),

(6)

also depends on

\hat{x}

⁠.

Initially, the molecule and the surfaces are assumed to be uncoupled, such that the total density matrix of the setup, ρ_T, factorizes

ρ_{T} (0) = ρ_{mol} (0) ρ_{met} (0) .

(7)

Here, ρ_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,

ρ_{met} (0) = \prod_{α} \frac{e^{- (H_{met, α} - μ_{α}) / k_{B} T}}{{tr}_{met, α} [e^{- (H_{met, α} - μ_{α}) / k_{B} T}]} .

(8)

The interaction between the molecule and the metal surfaces is completely described by the correlation function

\begin{align} C_{α, m m^{'}}^{σ} (\hat{x}, t - τ) = & \sum_{k} g_{α}^{2} (\hat{x}) t_{k α, m} t_{k α, m^{'}} \\ \times {tr}_{met} [c_{k α}^{σ} (t) c_{k α}^{\bar{σ}} (τ) ρ_{met} (0)], \end{align}

(9)

where σ = ± and

\bar{σ} = \mp

⁠,

c_{k α}^{-} = c_{k α}

and

c_{k α}^{+} = c_{k α}^{†}

⁠, and

c_{k α}^{σ} (t) = e^{i H_{met} t} c_{k α}^{σ} e^{- i H_{met} t}

⁠.

The correlation and level-width functions are connected by Fourier transform via

C_{α, m m^{'}}^{σ} (\hat{x}, t) = \frac{1}{2 π} \int d ϵ e^{i σ ϵ t} Γ_{α, {m m}^{'}} (ϵ, \hat{x}) f_{α}^{σ} (ϵ),

(10)

where

f_{α}^{σ} (ϵ) = \frac{1}{1 + e^{σ (ϵ - μ_{α}) / k_{B} T}},

(11)

is the Fermi–Dirac distribution of surface α. The correlation functions of the metal can be expanded as a sum of exponential functions,

C_{α, m m^{'}}^{σ} (t) = g_{α}^{2} (\hat{x}) t_{α, m} t_{α, m^{'}} \sum_{ℓ = 0}^{ℓ_{\max}} η_{α, σ, ℓ, m} e^{- κ_{α, σ, ℓ, m} t},

(12)

which is an essential step in the HEOM approach. To calculate the coefficients η_α,σ,ℓ,m and κ_α,σ,ℓ,m, we decompose the Fermi–Dirac distribution via a sum-over-poles approach based on the AAA method,⁷³ or a Padé decomposition,⁷⁴ 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. .

Although such decomposition methods can be applied to general

Γ_{α, m m^{'}} (ϵ, \hat{x})

⁠, in this work we exclusively consider the energy dependence of the level-width function to be a Lorentzian,

Γ_{α, m m^{'}} (ϵ, \hat{x}) = 2 π t_{α, m} t_{α, m^{'}} g_{α}^{2} (\hat{x}) \frac{W_{α}^{2}}{W_{α}^{2} + {(ϵ - μ_{α})}^{2}},

(13)

where W_α is the bandwidth. The constant Γ_α,mm′ = 2πt_α,mt_α,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

The dynamics of the entire molecule–metal system are described by the Liouville–von Neumann equation

\frac{\partial}{\partial t} ρ_{T} = - i [H, ρ_{T}],

(14)

where ρ_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,

ρ_{mol} = {tr}_{met} [ρ_{T}] .

(15)

The dynamics of the molecular density matrix is then given by

\frac{\partial}{\partial t} ρ_{mol} = - i {tr}_{met} ([H, ρ_{T}]) .

(16)

Although we are only interested in the dynamics of the molecule, the interaction of the molecule with the metal still influences the dynamics of the molecule itself, which makes the solution of Eq. challenging to obtain.

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.⁷⁷

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, $ρ_{j}^{(n)}$ ⁠. In this formalism, the zeroth tier ADO corresponds to the molecular density matrix, ρ⁽⁰⁾ = ρ_mol, while higher tier ADOs contain information about metal–molecule interactions and non-Markovian effects in the metal surfaces.²⁸ The ADOs are connected by a hierarchy of first-order coupled differential equations.

Given the form of metal–molecule coupling and the factorized initial state outlined in Sec. II, the equation of motion for an nth-tier ADO is

\begin{align} \frac{\partial}{\partial t} ρ_{j}^{(n)} = & - i [H_{mol}, ρ_{j}^{(n)}] - (\sum_{r = 1}^{n} κ_{j_{r}}) ρ_{j}^{(n)} \\ - i \sum_{r = 1}^{n} {(- 1)}^{n - r} C_{j_{r}} ρ_{j^{-}}^{(n - 1)} - i \sum_{j} A_{α, m}^{\bar{σ}} ρ_{j^{+}}^{(n + 1)}, \end{align}

(17)

with the coupling-up superoperator

A_{α, m}^{\bar{σ}}

and the coupling-down superoperator

C_{j_{r}}

defined by their action on

ρ_{j}^{(n)}

⁠,

A_{α, m}^{\bar{σ}} ρ_{j}^{(n)} = g_{α} (\hat{x}) t_{α, m} d_{m}^{\bar{σ}} ρ_{j}^{(n)} + {(- 1)}^{n} ρ_{j}^{(n)} g_{α} (\hat{x}) t_{α, m} d_{m}^{\bar{σ}},

(18)

and

C_{j_{r}} ρ_{j}^{(n)} = η_{j_{r}} g_{α} (\hat{x}) t_{α, m} d_{m}^{σ} ρ_{j}^{(n)} - {(- 1)}^{n} ρ_{j}^{(n)} η_{{\bar{j}}_{r}}^{*} g_{α} (\hat{x}) t_{α, m} d_{m}^{σ} .

(19)

Equations – use the superindex j = [j_n, …, j₁],

j_{r} = (α_{j_{r}}, σ_{j_{r}}, ℓ_{j_{r}}, m_{j_{r}})

⁠, with

{\bar{j}}_{r} = (α_{j_{r}}, {\bar{σ}}_{j_{r}}, ℓ_{j_{r}}, m_{j_{r}})

⁠. Moreover, the superindex j₋ = [j_n, …, j_r+1, j_r−1, …, j₁] is formed by removing the index j_r, and the superindex j₊ = [j, j_n, …, j₁] 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

We start by considering the partial Wigner transform of an operator O with respect to N vibrational DOFs,

O^{W} (x, p) = \frac{1}{{(2 π)}^{3 N}} \int d y e^{i p y} ⟨x - \frac{y}{2} |O| x + \frac{y}{2}⟩,

(20)

where x = (x₁, …, x_N), p = (p₁, …, 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.

The Wigner transform of the product of two operators, O and Q, is⁸¹

{(O Q)}^{W} = O^{W} e^{- i Λ / 2} Q^{W},

(21)

where

Λ = {\overset{⃖}{\nabla}}_{p} {\vec{\nabla}}_{x} - {\overset{⃖}{\nabla}}_{x} {\vec{\nabla}}_{p},

(22)

and the arrows indicate in which direction the derivatives act. A classical approximation is obtained by expanding Eq. up to first order in Λ, which results in

{(O Q)}^{W} \approx O^{W} Q^{W} - i \{O^{W}, Q^{W}\},

(23)

with the Poisson-bracket

\{O^{W}, Q^{W}\} = \sum_{i = 1}^{N} (\frac{\partial O^{W}}{\partial x_{i}} \frac{\partial Q^{W}}{\partial p_{i}} - \frac{\partial O^{W}}{\partial p_{i}} \frac{\partial Q^{W}}{\partial x_{i}}) .

(24)

Moreover, we introduce the notation

{\{O^{W}, Q^{W}\}}_{a} = \frac{1}{2} (\{O^{W}, Q^{W}\} - \{Q^{W}, O^{W}\}) .

(25)

Within this classical approximation for the vibrational DOFs, the time evolution of the Wigner-transformed total density matrix,

ρ^{W} (x, p, t)

⁠, follows the quantum–classical Liouville equation (QCLE),^78,82

\frac{\partial}{\partial t} ρ^{W} = - i [H^{W}, ρ_{T}^{W}] + {\{H^{W}, ρ_{T}^{W}\}}_{a} .

(26)

From this, the phase-space probability density (PSPD) of the now classical variables is introduced,

A (x, p, t) = {tr}_{{mol}^{qu} + met} [ρ_{T}^{W} (x, p, t)],

(27)

where we have used the notation

{tr}_{{mol}^{qu} + met} [\dots]

⁠, which denotes the trace over the electronic DOFs in the metal and the remaining quantum DOFs in the molecule after the Wigner transformation. We note that, in the case of a constant metal–molecule coupling, Eq. simplifies to the expression in Ref. 66,

A (x, p, t) = {tr}_{{mol}^{qu}} [ρ_{mol}^{W} (x, p, t)],

(28)

because all operators depending on x and p are contained only in the molecular DOFs.

The time evolution of the PSPD from Eq. can be obtained by tracing out the quantum DOFs from the QCLE in Eq. , which yields

\frac{\partial}{\partial t} A = {tr}_{{mol}^{qu} + met} [{\{H_{mol}^{W} + H_{met - mol}^{W}, ρ_{T}^{W}\}}_{a}] .

(29)

Although we have simplified the problem somewhat via the classical approximation for the vibrational DOFs, the metal–molecule interaction Hamiltonian

H_{met - mol}^{W}

still contains products of operators in the molecule and metal subspaces. Hence, performing the trace over the joint Hilbert space of molecule and metal is still a difficult task. We now show, therefore, how to express the time evolution of the PSPD in terms of the HEOM formalism.

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_{α} (\hat{x})$ affects the partial Wigner transform with respect to the vibrational DOFs.

In particular, the coupling-up,

A_{α, m}^{\bar{σ}} (\hat{x})

⁠, and coupling-down,

C_{j_{r}} (\hat{x})

⁠, superoperators now depend on

\hat{x}

⁠. Therefore, applying the partial Wigner transformation to the HEOM from Eq. gives additional terms compared to Ref. 66 and reads

\begin{aligned} \frac{\partial}{\partial t} ρ_{j}^{(n), W} & = - i [H_{mol}^{W}, ρ_{j}^{(n), W}] + {\{H_{mol}^{W}, ρ_{j}^{(n), W}\}}_{a} - (\sum_{r = 1}^{n} κ_{j_{r}}) ρ_{j}^{(n), W} \\ - i \sum_{r = 1}^{n} {(- 1)}^{n - r} (C_{j_{r}}^{W} ρ_{j^{-}}^{(n - 1), W} - \frac{i}{2} \sum_{i = 1}^{N} {(C_{j_{r}}^{W})}_{i}^{'} \frac{\partial}{\partial p_{i}} ρ_{j^{-}}^{(n - 1), W}) \\ - i \sum_{j} (A_{α, m}^{\bar{σ}, W} ρ_{j^{+}}^{(n + 1), W} - \frac{i}{2} \sum_{i = 1}^{N} {(A_{α, m}^{\bar{σ}, W})}_{i}^{'} \frac{\partial}{\partial p_{i}} ρ_{j^{+}}^{(n + 1), W}), \end{aligned}

(30)

where we have used the notation

\begin{align} {(A_{α, m}^{\bar{σ}, W})}_{i}^{'} \frac{\partial}{\partial p_{i}} ρ_{j}^{(n), W} = & \frac{\partial}{\partial x_{i}} g_{α} (x) t_{α, m} d_{m}^{\bar{σ}} \frac{\partial}{\partial p_{i}} ρ_{j}^{(n), W} \\ - {(- 1)}^{n} \frac{\partial}{\partial p_{i}} ρ_{j}^{(n), W} \frac{\partial}{\partial x_{i}} g_{α} (x) t_{α, m} d_{m}^{\bar{σ}}, \end{align}

(31)

and

\begin{aligned} {(C_{j_{r}}^{W})}_{i}^{'} ρ_{j}^{(n), W} & = η_{j_{r}} \frac{\partial}{\partial x_{i}} g_{α} (x) t_{α, m} d_{m}^{σ} \frac{\partial}{\partial p_{i}} ρ_{j}^{(n), W} \\ + {(- 1)}^{n} \frac{\partial}{\partial p_{i}} ρ_{j}^{(n), W} η_{{\bar{j}}_{r}}^{*} \frac{\partial}{\partial x_{i}} g_{α} (x) t_{α, m} d_{m}^{\bar{σ}} \end{aligned} .

(32)

If we unfold all ADOs into a vector,

ρ^{W} = {[\begin{matrix} ρ^{(0), W}, ρ_{j_{1}}^{(1), W}, \dots, ρ_{j}^{(n), W}, \dots \end{matrix}]}^{T},

(33)

then the Wigner-transformed HEOM from Eq. can be written in the joint Liouville space of all ADOs,

\begin{aligned} \frac{\partial}{\partial t} ρ^{W} & = - L^{W} ρ^{W} + {\{\{H_{mol}^{W}, ρ^{W}\}\}}_{a} \\ - \frac{i}{2} \sum_{i = 1}^{N} ({\tilde{C}}_{\partial x_{i}}^{W} + {\tilde{A}}_{\partial x_{i}}^{W}) \frac{\partial}{\partial p_{i}} ρ^{W}, \end{aligned}

(34)

where we have introduced the superoperator

L^{W} ρ^{W} = i [H_{mol}^{W}, ρ^{W}] + (\tilde{κ} - {\tilde{C}}^{W} - {\tilde{A}}^{W}) ρ^{W},

(35)

which comprises all terms of the zeroth-order expansion of the Wigner-transformed HEOM acting on

ρ^{W}

⁠. Moreover, we introduced the new superoperators

{\tilde{C}}^{W}

⁠,

{\tilde{A}}^{W}

⁠, and

\tilde{κ}

⁠, which are defined by their action on the vector of ADOs. The coupling-down superoperator in the joint Liouville space, for example, is

{\tilde{C}}^{W} [\begin{matrix} ρ^{(0), W} \\ ρ_{j_{1}}^{(1), W} \\ ⋮ \\ ρ_{j}^{(n), W} \\ ⋮ \end{matrix}] = [\begin{matrix} 0 \\ - i C^{W} ρ^{(0), W} \\ ⋮ \\ - i \sum_{r = 1}^{n} {(- 1)}^{n - r} C_{j_{r}}^{W} ρ_{j^{-}}^{(n - 1), W} \\ ⋮ \end{matrix}] .

(36)

The next term in Eq. applies a symmetrized Poisson bracket to each Wigner-transformed ADO

{\{\{H_{mol}^{W}, ρ^{W}\}\}}_{a} = {[{\{H_{mol}^{W}, ρ_{mol}^{W}\}}_{a}, {\{H_{mol}^{W}, ρ_{j_{1}}^{(1), W}\}}_{a}, \dots]}^{T},

(37)

and arises from the first-order term in the Wigner transform of the commutator with the molecular Hamiltonian in the HEOM. In a similar manner, the Wigner transform of the coupling-up and coupling-down superoperators also yields first-order terms,

{(A_{α, m}^{\bar{σ}, W})}_{i}^{'}

and

{(C_{j_{r}}^{W})}_{i}^{'}

⁠. In the joint Liouville space, these are collected into

{\tilde{C}}_{\partial x_{i}}^{W}

and

{\tilde{A}}_{\partial x_{i}}^{W}

⁠, which are also defined by their action on the vector of ADOs. The

{\tilde{C}}_{\partial x_{i}}^{W}

superoperator, for example, is

{\tilde{C}}_{\partial x_{i}}^{W} [\begin{matrix} ρ^{(0), W} \\ ρ_{j_{1}}^{(1), W} \\ ⋮ \\ ρ_{j}^{(n), W} \\ ⋮ \end{matrix}] = [\begin{matrix} 0 \\ - i {(C^{W})}_{i}^{'} ρ^{(0), W} \\ ⋮ \\ - i \sum_{r = 1}^{n} {(- 1)}^{n - r} {(C_{j_{r}}^{W})}_{i}^{'} ρ_{j^{-}}^{(n - 1), W} \\ ⋮ \end{matrix}] .

(38)

Explicit equations for all other superoperators are given in Appendix A.

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

To express the time evolution of the PSPD introduced in Eq. in terms of the Wigner-transformed HEOM, we follow a similar procedure as in Ref. 25. By comparing the time evolution of the molecular density matrix, governed by the Wigner-transformed Liouville equation, and the time evolution of the zeroth tier ADO in the Wigner-transformed HEOM from Eq. , we show in Appendix B that the time evolution of A(x, p, t) in the HEOM formalism can be expressed as

\begin{align} \frac{\partial A}{\partial t} = & {tr}_{{mol}^{qu}} [{\{\{H_{mol}^{W}, ρ^{W}\}\}}_{a}] \\ + {tr}_{{mol}^{qu}} [\sum_{i = 1}^{N} H_{met - mol, i}^{W} \frac{\partial ρ^{W}}{\partial p_{i}}] . \end{align}

(39)

Here, we have introduced the new superoperator

H_{met - mol, i}^{W}

⁠, which is defined via its action on the vector of all ADOs

H_{met - mol, i}^{W} ρ^{W} = [\begin{matrix} \sum_{j} \frac{\partial g_{α}}{\partial x_{i}} t_{α, m} d_{m}^{\bar{σ}} ρ_{j}^{(1), W} \\ 0 \\ 0 \\ ⋮ \end{matrix}],

(40)

and where the trace now selects and traces over only the zeroth-tier

{tr}_{{mol}^{qu}} [ρ^{W}] = {tr}_{{mol}^{qu}} [ρ^{(0), W}] .

(41)

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.

First, we write the vector of all ADOs with the ansatz⁵⁶

ρ^{W} (x, p, t) = A (x, p, t) σ_{ss}^{W} (x) + B (x, p, t),

(42)

where

σ_{ss}^{W} (x)

is the electronic steady state frozen at one vibrational frame,

L^{W} σ_{ss}^{W} = 0 .

(43)

One can view this term as the electronic steady state in the purely adiabatic limit of Eq. . The contribution B simply contains the difference between the exact vibronic dynamics and this adiabatic limit.

The key step is to recognize that, under the timescale separation assumption, 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. with Eq. , the time evolution of A in Eq. is given by the Fokker–Planck equation

\begin{align} \frac{\partial}{\partial t} A & = - \sum_{i} \frac{p_{i}}{m_{i}} \frac{\partial A}{\partial x_{i}} + \sum_{i} F_{i}^{ad} (x) \frac{\partial A}{\partial p_{i}} \\ + \sum_{i j} γ_{i j} (x) \frac{\partial}{\partial p_{i}} (\frac{p_{j}}{m_{j}} A) + \sum_{i j} D_{i j} (x) \frac{\partial^{2} A}{\partial p_{i} \partial p_{j}} . \end{align}

(44)

The dynamics of Eq. are governed by the adiabatic contribution to the mean force

F_{i}^{ad} (x)

⁠, the electronic friction tensor, γ_ij(x), and the correlation function of the stochastic force, D_ij(x). Although Eq. 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.

The adiabatic contribution to the mean force, for example, is now

F_{i}^{ad} (x) = F_{i, mol}^{ad} (x) + F_{i, met - mol}^{ad} (x),

(45)

where there are contributions originating not only from the molecular Hamiltonian, but also the metal–molecule coupling,

F_{i, mol}^{ad} (x) = - {tr}_{{mol}^{qu}} [\frac{\partial H_{mol}^{W}}{\partial x_{i}} σ_{ss}^{W}],

(46)

F_{i, met - mol}^{ad} (x) = - {tr}_{{mol}^{qu}} [H_{met - mol, i}^{W} σ_{ss}^{W}] .

(47)

Similarly, the electronic friction reads

γ_{i j} (x) = γ_{i j, mol} (x) + γ_{i j, met - mol} (x),

(48)

where

γ_{i j, mol} (x) = - {tr}_{{mol}^{qu}} [\frac{\partial H_{mol}^{W}}{\partial x_{i}} {(L^{W})}^{- 1} \frac{\partial σ_{ss}^{W}}{\partial x_{j}}],

(49)

γ_{i j, met - mol} (x) = - {tr}_{{mol}^{qu}} [H_{met - mol, i}^{W} {(L^{W})}^{- 1} \frac{\partial σ_{ss}^{W}}{\partial x_{j}}] .

(50)

For the correlation function of the random force, we obtain

\begin{align} D_{i j} (x) = & {tr}_{{mol}^{qu}} [(\frac{\partial H_{mol}^{W}}{\partial x_{i}} + H_{met - mol, i}^{W}) {(L^{W})}^{- 1} \\ \times ((- \frac{i}{2} ({\tilde{C}}_{\partial x_{j}}^{W} + {\tilde{A}}_{\partial x_{j}}^{W}) σ_{ss}^{W}) \\ + \frac{1}{2} (\frac{\partial H_{mol}^{W}}{{\partial x}_{j}} σ_{ss}^{W} - σ_{ss}^{W} \frac{\partial H_{mol}^{W}}{{\partial x}_{j}}) + F_{j}^{ad} (x) σ_{ss}^{W})] . \end{align}

(51)

In Eqs. –,

{(L^{W})}^{- 1}

refers to a formal inversion of Eq. .

The Langevin equation corresponding to the Fokker–Planck equation in Eq. is given by

m_{i} {\ddot{x}}_{i} = F_{i}^{ad} (x) + \sum_{j} γ_{i j} (x) {\dot{x}}_{j} + f_{i} (t),

(52)

where f_i(t) is a Gaussian random force with white noise

⟨ f_{i} (t) f_{j} (t^{'}) ⟩ = 2 D_{i j} (x) δ (t - t^{'}) .

(53)

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

First, we consider a single electronic level linearly coupled to a harmonic vibrational mode, with molecular Hamiltonian

H_{mol} = (ϵ_{0} + λ \hat{x}) d^{†} d + \frac{ω}{2} ({\hat{x}}^{2} + {\hat{p}}^{2}),

(54)

where we have implicitly used dimensionless coordinates,

\hat{x} \to \hat{x} \sqrt{m ω}

and

\hat{p} \to \hat{p} / \sqrt{m ω}

⁠. Moreover, ϵ₀ 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

H_{mol} = U_{g} (\hat{x}) d d^{†} + U_{c} (\hat{x}) d^{†} d + \frac{ω}{2} {\hat{p}}^{2},

(55)

which allows us to identify the potential energy surfaces of the neutral,

U_{g} (\hat{x}) = \frac{1}{2} ω {\hat{x}}^{2}

⁠, and charged,

U_{c} (\hat{x}) = U_{g} (\hat{x}) + λ \hat{x} + ϵ_{0}

⁠, states of the molecule, respectively.

Once the vibrational coordinate is treated classically, the electronic part of this Hamiltonian is noninteracting. Therefore, the electronic friction and other forces can be calculated exactly not only within the HEOM approach but also using NEGFs.⁵³ Following Ref. 53 and rewriting the electronic friction in our notation, one obtains the following formula at equilibrium in the wideband limit:

γ (x) = - \frac{1}{4 π} \int d ϵ {(\frac{\partial h}{\partial x} + \frac{ϵ - h (x)}{{(g (x))}^{2}} \frac{\partial {(g)}^{2}}{\partial x})}^{2} A^{2} (ϵ, x) \frac{\partial f (ϵ)}{\partial ϵ},

(56)

where

h (x) = U_{c} (x) - U_{g} (x) = ϵ_{0} + λ x

and

A (ϵ, x) = \frac{Γ \cdot {(g (x))}^{2}}{{(ϵ - h (x))}^{2} + {(Γ \cdot {(g (x))}^{2} / 2)}^{2}},

(57)

is the spectral function. Out of equilibrium, Eq. can be generalized as in Ref. 83. In the following, we will not only investigate the electronic forces of this model but also verify the extended HEOM-LD approach by comparing the electronic friction as calculated from Eq. to the expression in Eq. . To mimic the wideband approximation, we choose a bandwidth of 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_{α} (\hat{x}) = 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.

Figure 1(a) shows the potential energy surfaces corresponding to the two diabatic potential energy surfaces in the molecular Hamiltonian, as well as the adiabatic potential energy surface calculated from the adiabatic contribution to the mean force,

U_{ad} (x) = \int_{x_{0}}^{x} d x^{'} F^{ad} (x^{'}) .

(58)

The diabatic potential energy surfaces of the neutral and charged states cross at

x_{ε_{F}} = - (ϵ_{0} - ε_{F}) / λ,

(59)

which corresponds to the coordinate at which the electronic energy equals the Fermi energy of the leads. In the vicinity of

x_{ε_{F}}

⁠, the adiabatic potential shows a typical avoided crossing, which is due to the nonadiabatic coupling between the charged and uncharged diabats via the interaction with the metal. For large negative 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.

FIG. 1.

View large Download slide

Potential energy surfaces for the vibronic model (a) and friction and correlation function calculated with HEOM and NEGFs (b). The parameters are ϵ₀ = 0.05 eV, λ = 0.0075 eV, ω = 0.001 eV, and Γ_α = 0.0049 eV. The temperature of the metal is T = 300 K.

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_{ε_{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_{ε_{F}}$ ⁠.

Moreover, Fig. 1(b) shows that the friction and the correlation function of the stochastic force not only match between the HEOM and NEGF approaches, but that they also satisfy the fluctuation–dissipation theorem,⁸⁴

D (x) = k_{B} T γ (x) .

(60)

This reflects that the relaxation of vibrational modes via coupling to EHPs is balanced by the reverse process at equilibrium.

Next, we add a position-dependent metal–molecule coupling to the model. Although one can obtain such couplings from a first-principles calculation, here we choose a simple form of the coupling that is designed to reproduce the basic features one would expect in a molecular nanojunction. In particular, if the vibrational coordinate represents center-of-mass motion between the two leads, then the coupling

g_{α} (\hat{x})

increases and decreases as the molecule approaches and leaves lead α, respectively. The spatial part of the metal–molecule coupling, therefore, takes the form

g_{L} (\hat{x}) = [\frac{1 - q}{2} (1 - \tanh (\frac{\hat{x}}{a}) + q)],

(61)

and

g_{R} (\hat{x}) = - g_{L} (\hat{x}) + g_{c},

(62)

where we choose 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

Γ_{α} = 2 π t_{α}^{2} = 0.01 eV

⁠, since now for

\hat{x} \to \pm \infty

the molecule is essentially only coupled to one surface.

FIG. 2.

View large Download slide

Spatial part of the metal–molecule couplings g_L/R(x) for the vibronic model in Eq. (61). The parameters entering g_L/R(x) are q = 0.1 and a = 1.

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_{ε_{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, $F_{i, met - mol}^{ad} (x)$ ⁠, is proportional to the spatial derivative of the metal–molecule couplings, $\frac{\partial g_{α}}{\partial x}$ ⁠.

FIG. 3.

View large Download slide

Potential energy surfaces for the system with the position-dependent metal–molecule coupling from Fig. 2 (a), and friction and correlation function calculated with HEOM and NEGFs (b). Apart from setting Γ_α = 0.01 eV, all parameters are the same as in the previous model from Fig. 1.

At x = 0, both $\frac{\partial g_{L}}{\partial x}$ and $\frac{\partial g_{R}}{\partial x}$ 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).

FIG. 4.

View large Download slide

Contributions to the friction from the molecular Hamiltonian, γ_mol(x), and the metal–molecule coupling, γ_met–mol(x), as well as the friction for a constant coupling, γ_c(x), from the previous model in Fig. 1(b).

As in the constant coupling case, the full electronic friction also displays a peak at $x_{ε_{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 $\frac{\partial σ_{ss}^{(0), W}}{\partial x}$ ⁠, 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 $σ_{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 $\frac{\partial σ_{ss}^{(0), W}}{\partial x}$ and $\frac{\partial g_{α}}{\partial x}$ 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 $\frac{\partial g_{α}}{\partial x}$ or $\frac{\partial H_{mol}^{W}}{\partial x}$ ⁠, 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 $\frac{\partial σ_{ss}^{(0), W}}{\partial x}$ and $\frac{\partial σ_{ss, j}^{(1), W}}{\partial x}$ ⁠.

Next, we drive the nanojunction out of equilibrium by applying a bias voltage, Φ = 1 V. The corresponding electronic friction and correlation function of the stochastic force are shown in Fig. 5(a). Out of equilibrium, EHP creation and annihilation processes are resonant when the electronic energy crosses μ_α. Consequently, the main equilibrium friction peak originating from

\frac{\partial H_{mol}^{W}}{\partial x}

splits into two peaks at positions

x_{μ_{α}} = - (ϵ_{0} - ε_{F} - μ_{α}) / λ,

(63)

which has already been discussed in Ref. 66. Additionally, in Fig. 5(a), we observe that there is still a peak close to x_mm, which is now significantly more pronounced than the equilibrium case.

FIG. 5.

View large Download slide

Friction and correlation function out of equilibrium (a), and different contributions to the friction (b). Apart from setting Φ = 1 V, the parameters are the same as in Fig. 3.

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 $\frac{\partial σ_{ss}^{(0), W}}{\partial x}$ ⁠.

Finally in this subsection, we note that the fluctuation–dissipation theorem is broken out of equilibrium, such that γ(x) ≠ D(x)/k_BT. 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_{μ_{R}}$ and $x_{μ_{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_BT > γ(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

Here, we consider a strongly interacting system, for which one cannot calculate the electronic forces exactly using NEGFs. In contrast, the HEOM-LD approach to electronic forces remains numerically exact even for interacting systems, demonstrating the usefulness of the extended HEOM-LD theory. Specifically, we consider a molecular model in which an electronic level is linearly coupled to one low-frequency vibrational mode, which we will treat classically, and one high-frequency vibrational mode, which we will treat quantum mechanically alongside the electronic DOFs. The molecular Hamiltonian is given by

\begin{align} H_{mol} = & (ϵ_{0} + λ_{cl} {\hat{x}}_{cl}) d^{†} d + \frac{1}{2} ω_{cl} ({\hat{x}}_{cl}^{2} + {\hat{p}}_{cl}^{2}) \\ + λ_{qu} {\hat{y}}_{qu} d^{†} d + \frac{1}{2} ω_{qu} ({\hat{y}}_{qu}^{2} + {\hat{p}}_{qu}^{2}), \end{align}

(64)

where ϵ₀ 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,⁶⁰ or a molecule with multiple modes operating on different timescales coupled to the transport.⁶⁷ 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, ϵ₀, such that the energy after small polaron shift is $ϵ = ϵ_{0} - \frac{λ_{qu}^{2}}{ω_{qu}} = 0.05 e V$ ⁠, 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.

FIG. 6.

View large Download slide

Potential energy surfaces of the classical vibrational DOFs for the interacting system (a) and friction and correlation function (b). The parameters are ϵ = 0.05 eV, λ_cl = 0.0075 eV, ω_cl = 0.001 eV, λ_qu = 0.55 eV, and ω_qu = 0.3 eV.

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_{ε_{F}} = - (ϵ - ε_{F}) / λ_{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_{μ_{L}}$ and $x_{μ_{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.

Figure 7 shows additional features, which can only be observed for the interacting system. These additional extrema in the friction appear at positions⁶⁶

x_{μ_{α}, n} = \frac{μ_{α} - (ϵ_{0} - ε_{F} - λ_{qu}^{2} / ω_{qu} + n ω_{qu})}{λ_{cl}} .

(65)

At these positions, EHP creation and annihilation processes mediated by the exchange of 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

FIG. 7.

View large Download slide

Friction and the correlation function for the interacting system out of equilibrium. The bias voltage is Φ = 1 V. All other parameters are the same as in Fig. 6.

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

The superoperators

{\tilde{A}}^{W}

⁠,

{\tilde{A}}_{\partial x_{i}}^{W}

⁠, and

\tilde{κ}

introduced in Eq. are defined by their action on the vector of all ADOs,

{\tilde{A}}^{W} [\begin{matrix} ρ^{(0), W} \\ ρ_{j_{0}}^{(1), W} \\ ⋮ \\ ρ_{j}^{(n), W} \\ ⋮ \end{matrix}] = [\begin{matrix} - i \sum_{j} A_{α, m}^{\bar{σ}, W} ρ_{j}^{(1), W} \\ - i \sum_{j} A_{α, m}^{\bar{σ}, W} ρ_{j j_{0}}^{(2), W} \\ ⋮ \\ - i \sum_{j} A_{α, m}^{\bar{σ}, W} ρ_{j^{+}}^{(n + 1), W} \\ ⋮ \end{matrix}],

(A1)

{\tilde{A}}_{\partial x_{i}}^{W} [\begin{matrix} ρ^{(0), W} \\ ρ_{j_{0}}^{(1), W} \\ ⋮ \\ ρ_{j}^{(n), W} \\ ⋮ \end{matrix}] = [\begin{matrix} - i \sum_{j} {(A_{α, m}^{\bar{σ}, W})}_{i}^{'} ρ_{j}^{(1), W} \\ - i \sum_{j} {(A_{α, m}^{\bar{σ}, W})}_{i}^{'} ρ_{j j_{0}}^{(2), W} \\ ⋮ \\ - i \sum_{j} {(A_{α, m}^{\bar{σ}, W})}_{i}^{'} ρ_{j^{+}}^{(n + 1), W} \\ ⋮ \end{matrix}],

(A2)

and,

\tilde{κ} [\begin{matrix} ρ^{(0), W} \\ ρ_{j_{1}}^{(1), W} \\ ⋮ \\ ρ_{j}^{(n), W} \\ ⋮ \end{matrix}] = [\begin{matrix} 0 \\ κ_{j_{1}} ρ_{j_{1}}^{(1), W} \\ ⋮ \\ \sum_{r = 1}^{n} κ_{j_{r}} ρ_{j}^{(n), W} \\ ⋮ \end{matrix}],

(A3)

where we have used the notation introduced in Eqs. and .

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.

We start by considering the time evolution of the molecular density matrix, given by the Wigner-transformed Liouville–von Neumann equation

\begin{align} \frac{\partial}{\partial t} ρ_{mol}^{W} = & {tr}_{met} (- i [H_{mol}^{W}, ρ_{T}^{W}] - i [H_{met - mol}^{W}, ρ_{T}^{W}] \\ + {\{H_{mol}^{W}, ρ_{T}^{W}\}}_{a} + {\{H_{met - mol}^{W}, ρ_{T}^{W}\}}_{a}) . \end{align}

(B1)

This is equivalent to the time evolution of the Wigner-transformed zeroth-tier ADO,

\begin{align} \frac{\partial}{\partial t} ρ^{(0), W} = & - i [H_{mol}^{W}, ρ^{(0), W}] + {\{H_{mol}^{W}, ρ^{(0), W}\}}_{a} \\ - i {(\sum_{j} A_{m, α}^{\bar{σ}} ρ_{j}^{(1)})}^{W}, \end{align}

(B2)

where the Wigner transform of the coupling-up superoperator is evaluated explicitly as

\begin{align} {(\sum_{j} A_{m, α}^{\bar{σ}} ρ_{j}^{(1)})}^{W} & = \sum_{j} t_{α, m} (g_{α} (\hat{x}) d_{m}^{\bar{σ}} ρ_{j}^{(1)} - {ρ_{j}^{(1)} g_{α} (\hat{x}) d_{m}^{\bar{σ}})}^{W} \\ = \sum_{j} g_{α} (x) t_{α, m} (d_{m}^{\bar{σ}} ρ_{j}^{(1), W} - ρ_{j}^{(1), W} d_{m}^{\bar{σ}}) \\ - \frac{i}{2} \sum_{j} \sum_{i = 1}^{N} \frac{\partial g_{α}}{\partial x_{i}} t_{α, m} (d_{m}^{\bar{σ}} \frac{\partial ρ_{j}^{(1), W}}{\partial p_{i}} - \frac{\partial ρ_{j}^{(1), W}}{\partial p_{i}} d_{m}^{\bar{σ}}) . \end{align}

(B3)

Now, we relate the different terms in Eq. to their respective counterparts in Eq. ,

- i {(\sum_{j} A_{m, α}^{\bar{σ}} ρ_{j}^{(1)})}^{W} = {tr}_{met} (i [H_{met - mol}^{W}, ρ_{T}^{W}] + {\{H_{met - mol}^{W}, ρ_{T}^{W}\}}_{a}) .

(B4)

We write out the {⋯}_a brackets on the right hand side of Eq. and match the terms on both sides of Eq. depending on

\frac{\partial}{\partial p_{i}} ρ_{j}^{(1), W}

and

\frac{\partial}{\partial p_{i}} ρ_{T}^{W}

⁠,

\begin{align} \sum_{j} \sum_{i = 1}^{N} \frac{\partial g_{α}}{\partial x_{i}} t_{α, m} (d_{m}^{\bar{σ}} \frac{\partial ρ_{j}^{(1), W}}{\partial p_{i}} - \frac{\partial ρ_{j}^{(1), W}}{\partial p_{i}} d_{m}^{\bar{σ}}) \\ = {tr}_{met} [\sum_{i = 1}^{N} (\frac{\partial H_{met - mol}^{W}}{\partial x_{i}} \frac{\partial ρ_{T}^{W}}{\partial p_{i}} - \frac{\partial ρ_{T}^{W}}{\partial p_{i}} \frac{\partial H_{met - mol}^{W}}{\partial x_{i}})] . \end{align}

(B5)

Writing out the explicit form of

H_{met - mol}^{W}

and comparing both sides in Eq. then implies

\sum_{j} \frac{\partial g_{α}}{\partial x_{i}} t_{α, m} d_{m}^{\bar{σ}} \frac{\partial ρ_{j}^{(1), W}}{\partial p_{i}} = \sum_{α, m} \frac{\partial g_{α}}{\partial x_{i}} t_{α, m} {tr}_{met} (c_{α}^{σ} d_{m}^{\bar{σ}} \frac{\partial ρ_{T}^{W}}{\partial p_{i}}) .

(B6)

Equation allows us to express the trace over the metal DOFs in terms of the first-tier ADOs.

We now consider again the expression for the time evolution of the PSPD from Eq. . The first term on the right hand side of Eq. can be written as

{tr}_{{mol}^{qu} + met} [{\{H_{mol}^{W}, ρ_{T}\}}_{a}] = {tr}_{{mol}^{qu}} [{\{H_{mol}^{W}, ρ_{mol}^{W}\}}_{a}] .

(B7)

By exploting the connection between the trace over the metal DOFs and the first tier ADOs from Eq. , we can write the second term on the right hand side of Eq. as

\begin{aligned} {tr}_{{mol}^{qu} + met} [{H_{met - mol}^{W}, ρ_{T}^{W}}_{a}] = {tr}_{{mol}^{qu}} [\sum_{α, m} \sum_{i = 1}^{N} \frac{\partial g_{α} (x)}{\partial x_{i}} t_{α, m} \\ \times (d_{m}^{†} {tr}_{met} [c_{α} \frac{\partial ρ_{T}^{W}}{\partial p_{i}}] + d_{m} {tr}_{met} [c_{α}^{†} \frac{\partial ρ_{T}^{W}}{\partial p_{i}}])] \\ = {tr}_{{mol}^{W}} [\sum_{α, σ, l, m} \sum_{i = 1}^{N} \frac{\partial g_{α} (x)}{\partial x_{i}} t_{α, m} d_{m}^{\bar{σ}} \frac{\partial}{\partial p_{i}} ρ_{j}^{(1), W}] . \end{aligned}

(B8)

Finally, adding Eqs. and , the time evolution of the PSPD can be written as

\frac{\partial A}{\partial t} = {tr}_{{mol}^{qu}} [{\{H_{mol}^{W}, ρ_{S}^{W}\}}_{a}] + {tr}_{{mol}^{qu}} [\sum_{j} \sum_{i = 1}^{N} \frac{\partial g_{α}}{\partial x_{i}} t_{α, m} d_{m}^{\bar{σ}} \frac{\partial ρ_{j}^{(1), W}}{\partial p_{i}}] .

(B9)

Considering the vector of all ADOs introduced in Eq. and the superoperator

H_{met - mol, i}^{W}

⁠, which has been introduced in Eq. , we obtain the expression in Eq. .

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.

We assume that the dynamics of the quantum mechanical DOFs and the classical DOFs take place on timescales τ_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. . Recalling the time evolution of the PSPD in terms of the HEOM from Eq. and using Eq. , we obtain

\frac{\partial A}{\partial t} = {tr}_{{mol}^{qu}} [{\{\{H_{mol}^{W}, ρ^{W}\}\}}_{a}] + {tr}_{{mol}^{qu}} [\sum_{i = 1}^{N} H_{met - mol, i}^{W} \frac{\partial ρ^{W}}{\partial p_{i}}],

(C1)

which, when we insert the ansatz from Eq. , becomes

\begin{align} \frac{\partial A}{\partial t} = & {tr}_{{mol}^{qu}} [{\{\{H_{mol}^{W}, A σ_{ss}^{W}\}\}}_{a}] + {tr}_{{mol}^{qu}} [{\{\{H_{mol}^{W}, B\}\}}_{a}] \\ + {tr}_{{mol}^{qu}} [\sum_{i = 1}^{N} H_{met - mol, i}^{W} \frac{\partial}{\partial p_{i}} A σ_{ss}^{W}] \\ + {tr}_{{mol}^{qu}} [\sum_{i = 1}^{N} H_{met - mol, i}^{W} \frac{\partial}{\partial p_{i}} B] . \end{align}

(C2)

We now write out the antisymmetrized Poisson brackets explicitly,

\begin{align} \frac{\partial A}{\partial t} = & - \sum^{N} i = 1 (\frac{p i}{m_{i}} \frac{\partial A}{\partial x_{i}} + {tr}_{{mol}^{qu}} [\frac{\partial H_{mol}^{W}}{\partial x_{i}} σ_{ss}^{W}] \frac{\partial A}{\partial p i} \\ + {tr}_{{mol}^{qu}} [\frac{\partial H_{mol}^{W}}{\partial x_{i}} \frac{\partial B}{\partial p_{i}}] + {tr}_{{mol}^{qu}} [H_{met - mol, i}^{W} \frac{\partial}{\partial p_{i}} A σ_{ss}^{W}] \\ + {tr}_{{mol}^{qu}} [H_{met - mol, i}^{W} \frac{\partial}{\partial p_{i}} B]) . \end{align}

(C3)

Next, we rearrange Eq. to obtain the time evolution of B,

\frac{\partial}{\partial t} B = \frac{\partial}{\partial t} ρ^{W} - \frac{\partial}{\partial t} A σ_{ss}^{W},

(C4)

which we can express in terms of the Wigner-transformed HEOM from Eq. ,

\begin{align} \frac{\partial}{\partial t} B = & - L^{W} B + {\{H_{mol}^{W}, (A σ_{ss}^{W} + B)\}}_{a} \\ - \frac{i}{2} \sum_{i = 1}^{N} ({\tilde{C}}_{\partial x_{i}}^{W} + {\tilde{A}}_{\partial x_{i}}^{W}) \frac{\partial}{\partial p_{i}} (A σ_{ss}^{W} + B) - \frac{\partial}{\partial t} A σ_{ss}^{W}, \end{align}

(C5)

where we have also used Eq. . We now formally solve Eq. for B, obtaining

\begin{align} B = & {(L^{W})}^{- 1} (- \frac{\partial}{\partial t} B + {\{H_{mol}^{W}, (A σ_{ss}^{W} + B)\}}_{a} \\ - \frac{i}{2} \sum_{i = 1}^{N} ({\tilde{C}}_{\partial x_{i}}^{W} + {\tilde{A}}_{\partial x_{i}}^{W}) \frac{\partial}{\partial p_{i}} (A σ_{ss}^{W} + B) - \frac{\partial}{\partial t} A σ_{ss}^{W}) . \end{align}

(C6)

We note that the pseudoinverse of the joint Liouville operator,

{(L^{W} (x))}^{- 1}

⁠, can be calculated via the Laplace transform

{(L^{W} (x))}^{- 1} = \lim_{δ \to 0_{+}} \int_{0}^{\infty} d t e^{- (L^{W} (x) + δ) t} .

(C7)

Now we express

\frac{\partial}{\partial t} A

in Eq. through Eq. and insert the whole expression for B into Eq. . Thereby, we assume

\frac{\partial}{\partial t} B ≪ \frac{\partial}{\partial t} A σ_{ss}^{W}

⁠,

\frac{\partial}{\partial x_{i}} B ≪ \frac{\partial}{\partial x_{i}} (A σ_{ss}^{W})

⁠, and

\frac{\partial}{\partial p_{i}} B ≪ \frac{\partial}{\partial p_{i}} A σ_{ss}^{W}

⁠, and that, therefore, derivatives of B can be neglected. Under this assumption, it becomes apparent that only four terms contribute to B, so that

\begin{align} B = & - \sum_{i = 1}^{N} \frac{p_{i}}{m_{i}} {(L^{W})}^{- 1} \frac{\partial σ_{ss}^{W}}{\partial x_{i}} A \\ + \frac{1}{2} {(L^{W})}^{- 1} \sum_{i = 1}^{N} (\frac{\partial H_{mol}^{W}}{\partial x_{i}} σ_{ss}^{W} - σ_{ss}^{W} \frac{\partial H_{mol}^{W}}{\partial x_{i}}) \frac{\partial A}{\partial p_{i}} \\ - \frac{i}{2} {(L^{W})}^{- 1} \sum_{i = 1}^{N} ({\tilde{C}}_{\partial x_{i}}^{W} + {\tilde{A}}_{\partial x_{i}}^{W}) σ_{ss}^{W} \frac{\partial A}{\partial p_{i}} \\ + {(L^{W})}^{- 1} \sum_{i = 1}^{N} (- \frac{p_{i}}{m_{i}} \frac{\partial A}{\partial x_{i}} + {tr}_{{mol}^{qu}} [\frac{\partial H_{mol}^{W}}{\partial x_{i}} σ_{ss}^{W}] \frac{\partial A}{\partial p_{i}} \\ + {tr}_{{mol}^{qu}} [H_{met - mol, i}^{W} \frac{\partial}{\partial p_{i}} A σ_{ss}^{W}]) σ_{ss}^{W} . \end{align}

(C8)

The first two terms on the right hand side of Eq. come from the antisymmetrized Poisson bracket in Eq. and also appear in the constant metal–molecule coupling case. The third term comes from the inclusion of a x-dependency in the metal–molecule coupling, and the last term comes from inserting Eq. into Eq. and neglecting derivatives of B. Inserting these contributions in the last line of Eq. , we obtain the Fokker–Planck equation from Eq. .

REFERENCES

1.

R. J.

Maurer

,

B.

Jiang

,

H.

Guo

, and

J. C.

Tully

,

Phys. Rev. Lett.

118

,

256001

(

2017

).

https://doi.org/10.1103/physrevlett.118.256001

Google Scholar

Crossref

PubMed

2.

M.

Askerka

,

R. J.

Maurer

,

V. S.

Batista

, and

J. C.

Tully

,

Phys. Rev. Lett.

116

,

217601

(

2016

).

https://doi.org/10.1103/physrevlett.116.217601

Google Scholar

Crossref

PubMed

3.

Y.

Zhang

,

R. J.

Maurer

,

H.

Guo

, and

B.

Jiang

,

Chem. Sci.

10

,

1089

(

2019

).

https://doi.org/10.1039/c8sc03955k

Google Scholar

Crossref

PubMed

4.

Y.

Zhang

,

C. L.

Box

,

T.

Schäfer

,

A.

Kandratsenka

,

A. M.

Wodtke

,

R. J.

Maurer

, and

B.

Jiang

,

Phys. Chem. Chem. Phys.

24

,

19753

(

2022

).

https://doi.org/10.1039/d2cp03312g

Google Scholar

Crossref

PubMed

5.

R. B.

Gerber

,

Chem. Rev.

87

,

29

(

1987

).

https://doi.org/10.1021/cr00077a003

Google Scholar

Crossref

6.

M.

Brandbyge

,

P.

Hedegård

,

T. F.

Heinz

,

J. A.

Misewich

, and

D. M.

Newns

,

Phys. Rev. B

52

,

6042

(

1995

).

https://doi.org/10.1103/physrevb.52.6042

Google Scholar

Crossref

7.

Y.

Kim

,

T.

Komeda

, and

M.

Kawai

,

Phys. Rev. Lett.

89

,

126104

(

2002

).

https://doi.org/10.1103/physrevlett.89.126104

Google Scholar

Crossref

PubMed

8.

A.

von Hippel

,

Science

123

,

315

(

1956

).

https://doi.org/10.1126/science.123.3191.315

Google Scholar

Crossref

PubMed

9.

A.

Aviram

and

M. A.

Ratner

,

Chem. Phys. Lett.

29

,

277

(

1974

).

https://doi.org/10.1016/0009-2614(74)85031-1

Google Scholar

Crossref

10.

J. C.

Cuevas

and

E.

Scheer

,

Molecular Electronics

(

World Scientific

,

2010

).

Google Scholar

Crossref

11.

N. A.

Zimbovskaya

and

M. R.

Pederson

,

Phys. Rep.

509

,

1

(

2011

).

https://doi.org/10.1016/j.physrep.2011.08.002

Google Scholar

Crossref

12.

J. P.

Bergfield

and

M. A.

Ratner

,

Phys. Status Solidi B

250

,

2249

(

2013

).

https://doi.org/10.1002/pssb.201350048

Google Scholar

Crossref

13.

M.

Thoss

and

F.

Evers

,

J. Chem. Phys.

148

,

030901

(

2018

).

https://doi.org/10.1063/1.5003306

Google Scholar

Crossref

14.

M.

Galperin

,

M. A.

Ratner

, and

A.

Nitzan

,

J. Phys.: Condens. Matter

19

,

103201

(

2007

).

https://doi.org/10.1088/0953-8984/19/10/103201

Google Scholar

Crossref

15.

R.

Härtle

and

M.

Thoss

,

Phys. Rev. B

83

,

125419

(

2011

).

https://doi.org/10.1103/physrevb.83.125419

Google Scholar

Crossref

16.

A.

Erpenbeck

,

C.

Schinabeck

,

U.

Peskin

, and

M.

Thoss

,

Phys. Rev. B

97

,

235452

(

2018

).

https://doi.org/10.1103/physrevb.97.235452

Google Scholar

Crossref

17.

A.

Erpenbeck

,

Y.

Ke

,

U.

Peskin

, and

M.

Thoss

,

Phys. Rev. B

102

,

195421

(

2020

).

https://doi.org/10.1103/physrevb.102.195421

Google Scholar

Crossref

18.

Y.

Ke

,

A.

Erpenbeck

,

U.

Peskin

, and

M.

Thoss

,

J. Chem. Phys.

154

,

234702

(

2021

).

https://doi.org/10.1063/5.0053828

Google Scholar

Crossref

19.

Y.

Ke

,

J.

Dvořák

,

M.

Čížek

,

R.

Borrelli

, and

M.

Thoss

,

J. Chem. Phys.

159

,

024703

(

2023

).

https://doi.org/10.1063/5.0155290

Google Scholar

Crossref

20.

B.

Persson

and

P.

Avouris

,

Surf. Sci.

390

,

45

(

1997

).

https://doi.org/10.1016/s0039-6028(97)00507-4

Google Scholar

Crossref

21.

M.

Ring

,

D.

Weber

,

P.

Haiber

,

F.

Pauly

,

P.

Nielaba

, and

E.

Scheer

,

Nano Lett.

20

,

5773

(

2020

).

https://doi.org/10.1021/acs.nanolett.0c01597

Google Scholar

Crossref

PubMed

22.

C.

Sabater

,

C.

Untiedt

, and

J. M.

van Ruitenbeek

,

Beilstein J. Nanotechnol.

6

,

2338

(

2015

).

https://doi.org/10.3762/bjnano.6.241

Google Scholar

Crossref

PubMed

23.

Y.

Tanimura

and

R.

Kubo

,

J. Phys. Soc. Jpn.

58

,

101

(

1989

).

https://doi.org/10.1143/jpsj.58.101

Google Scholar

Crossref

24.

C.

Schinabeck

and

M.

Thoss

,

Phys. Rev. B

101

,

075422

(

2020

).

https://doi.org/10.1103/physrevb.101.075422

Google Scholar

Crossref

25.

J.

Jin

,

X.

Zheng

, and

Y.

Yan

,

J. Chem. Phys.

128

,

234703

(

2008

).

https://doi.org/10.1063/1.2938087

Google Scholar

Crossref

26.

R.

Härtle

,

G.

Cohen

,

D. R.

Reichman

, and

A. J.

Millis

,

Phys. Rev. B

88

,

235426

(

2013

).

https://doi.org/10.1103/physrevb.88.235426

Google Scholar

Crossref

27.

R.

Härtle

,

C.

Schinabeck

,

M.

Kulkarni

,

D.

Gelbwaser-Klimovsky

,

M.

Thoss

, and

U.

Peskin

,

Phys. Rev. B

98

,

081404

(

2018

).

https://doi.org/10.1103/physrevb.98.081404

Google Scholar

Crossref

28.

C.

Schinabeck

,

A.

Erpenbeck

,

R.

Härtle

, and

M.

Thoss

,

Phys. Rev. B

94

,

201407

(

2016

).

https://doi.org/10.1103/physrevb.94.201407

Google Scholar

Crossref

29.

H.

Wang

and

M.

Thoss

,

J. Chem. Phys.

131

,

024114

(

2009

).

https://doi.org/10.1063/1.3173823

Google Scholar

Crossref

30.

H.

Wang

,

I.

Pshenichnyuk

,

R.

Härtle

, and

M.

Thoss

,

J. Chem. Phys.

135

,

244506

(

2011

).

https://doi.org/10.1063/1.3660206

Google Scholar

Crossref

31.

H.

Wang

and

M.

Thoss

,

J. Chem. Phys.

145

,

164105

(

2016

).

https://doi.org/10.1063/1.4965712

Google Scholar

Crossref

32.

H.-H.

Teh

,

W.

Dou

, and

J. E.

Subotnik

,

Phys. Rev. B

104

,

L201409

(

2021

).

https://doi.org/10.1103/physrevb.104.l201409

Google Scholar

Crossref

33.

H.-H.

Teh

,

W.

Dou

, and

J. E.

Subotnik

,

Phys. Rev. B

106

,

184302

(

2022

).

https://doi.org/10.1103/physrevb.106.184302

Google Scholar

Crossref

34.

J.-T.

Lü

,

M.

Brandbyge

, and

P.

Hedegård

,

Nano Lett.

10

,

1657

(

2010

).

https://doi.org/10.1021/nl904233u

Google Scholar

Crossref

PubMed

35.

T. N.

Todorov

,

D.

Dundas

, and

E. J.

McEniry

,

Phys. Rev. B

81

,

075416

(

2010

).

https://doi.org/10.1103/physrevb.81.075416

Google Scholar

Crossref

36.

J.-T.

Lü

,

T.

Gunst

,

P.

Hedegård

, and

M.

Brandbyge

,

Beilstein J. Nanotechnol.

2

,

814

(

2011

).

https://doi.org/10.3762/bjnano.2.90

Google Scholar

Crossref

PubMed

37.

R. J.

Preston

,

T. D.

Honeychurch

, and

D. S.

Kosov

,

Phys. Rev. B

106

,

195406

(

2022

).

https://doi.org/10.1103/physrevb.106.195406

Google Scholar

Crossref

38.

A. P.

Horsfield

,

D. R.

Bowler

,

A. J.

Fisher

,

T. N.

Todorov

, and

C. G.

Sánchez

,

J. Phys.: Condens. Matter

16

,

8251

(

2004

).

https://doi.org/10.1088/0953-8984/16/46/012

Google Scholar

Crossref

39.

A. P.

Horsfield

,

D. R.

Bowler

, and

A. J.

Fisher

,

J. Phys.: Condens. Matter

16

,

L65

(

2004

).

https://doi.org/10.1088/0953-8984/16/7/l03

Google Scholar

Crossref

40.

C.

Verdozzi

,

G.

Stefanucci

, and

C.-O.

Almbladh

,

Phys. Rev. Lett.

97

,

046603

(

2006

).

https://doi.org/10.1103/physrevlett.97.046603

Google Scholar

Crossref

PubMed

41.

A.

Kartsev

,

C.

Verdozzi

, and

G.

Stefanucci

,

Eur. Phys. J. B

87

,

14

(

2014

).

https://doi.org/10.1140/epjb/e2013-40905-5

Google Scholar

Crossref

42.

J. E.

Subotnik

,

J. Chem. Phys.

132

,

134112

(

2010

).

https://doi.org/10.1063/1.3314248

Google Scholar

Crossref

43.

T. N.

Todorov

,

D.

Dundas

,

J.-T.

Lü

,

M.

Brandbyge

, and

P.

Hedegård

,

Eur. J. Phys.

35

,

065004

(

2014

).

https://doi.org/10.1088/0143-0807/35/6/065004

Google Scholar

Crossref

44.

B.

Cunningham

,

T. N.

Todorov

, and

D.

Dundas

,

Beilstein J. Nanotechnol.

6

,

2140

(

2015

).

https://doi.org/10.3762/bjnano.6.219

Google Scholar

Crossref

PubMed

45.

N.

Bellonzi

,

A.

Jain

, and

J. E.

Subotnik

,

J. Chem. Phys.

144

,

154110

(

2016

).

https://doi.org/10.1063/1.4946810

Google Scholar

Crossref

46.

J. C.

Tully

and

R. K.

Preston

,

J. Chem. Phys.

55

,

562

(

1971

).

https://doi.org/10.1063/1.1675788

Google Scholar

Crossref

47.

J. C.

Tully

,

J. Chem. Phys.

93

,

1061

(

1990

).

https://doi.org/10.1063/1.459170

Google Scholar

Crossref

48.

W.

Ouyang

,

W.

Dou

, and

J. E.

Subotnik

,

J. Chem. Phys.

142

,

084109

(

2015

).

https://doi.org/10.1063/1.4908032

Google Scholar

Crossref

49.

N.

Shenvi

,

S.

Roy

, and

J. C.

Tully

,

Science

326

,

829

(

2009

).

https://doi.org/10.1126/science.1179240

Google Scholar

Crossref

PubMed

50.

M.

Head-Gordon

and

J. C.

Tully

,

J. Chem. Phys.

103

,

10137

(

1995

).

https://doi.org/10.1063/1.469915

Google Scholar

Crossref

51.

N.

Bode

,

S. V.

Kusminskiy

,

R.

Egger

, and

F.

von Oppen

,

Beilstein J. Nanotechnol.

3

,

144

(

2012

).

https://doi.org/10.3762/bjnano.3.15

Google Scholar

Crossref

PubMed

52.

J.-T.

Lü

,

M.

Brandbyge

,

P.

Hedegård

,

T. N.

Todorov

, and

D.

Dundas

,

Phys. Rev. B

85

,

245444

(

2012

).

https://doi.org/10.1103/physrevb.85.245444

Google Scholar

Crossref

53.

W.

Dou

and

J. E.

Subotnik

,

J. Chem. Phys.

146

,

092304

(

2016

).

https://doi.org/10.1063/1.4965823

Google Scholar

Crossref

54.

W.

Dou

and

J. E.

Subotnik

,

J. Chem. Phys.

145

,

054102

(

2016

).

https://doi.org/10.1063/1.4959604

Google Scholar

Crossref

55.

R. J.

Maurer

,

M.

Askerka

,

V. S.

Batista

, and

J. C.

Tully

,

Phys. Rev. B

94

,

115432

(

2016

).

https://doi.org/10.1103/physrevb.94.115432

Google Scholar

Crossref

56.

W.

Dou

,

G.

Miao

, and

J. E.

Subotnik

,

Phys. Rev. Lett.

119

,

046001

(

2017

).

https://doi.org/10.1103/physrevlett.119.046001

Google Scholar

Crossref

PubMed

57.

W.

Dou

and

J. E.

Subotnik

,

Phys. Rev. B

96

,

104305

(

2017

).

https://doi.org/10.1103/physrevb.96.104305

Google Scholar

Crossref

58.

W.

Dou

and

J. E.

Subotnik

,

J. Chem. Phys.

148

,

230901

(

2018

).

https://doi.org/10.1063/1.5035412

Google Scholar

Crossref

59.

F.

Chen

,

K.

Miwa

, and

M.

Galperin

,

J. Phys. Chem. A

123

,

693

(

2018

).

https://doi.org/10.1021/acs.jpca.8b09251

Google Scholar

Crossref

60.

F.

Chen

,

K.

Miwa

, and

M.

Galperin

,

J. Chem. Phys.

150

,

174101

(

2019

).

https://doi.org/10.1063/1.5095425

Google Scholar

Crossref

61.

R. J.

Preston

,

M. F.

Gelin

, and

D. S.

Kosov

,

J. Chem. Phys.

154

,

114108

(

2021

).

https://doi.org/10.1063/5.0045652

Google Scholar

Crossref

62.

A.

Nocera

,

C. A.

Perroni

,

V.

Marigliano Ramaglia

, and

V.

Cataudella

,

Phys. Rev. B

83

,

115420

(

2011

).

https://doi.org/10.1103/physrevb.83.115420

Google Scholar

Crossref

63.

S.

Weiss

,

J.

Eckel

,

M.

Thorwart

, and

R.

Egger

,

Phys. Rev. B

77

,

195316

(

2008

).

https://doi.org/10.1103/physrevb.77.195316

Google Scholar

Crossref

64.

L.

Mühlbacher

and

E.

Rabani

,

Phys. Rev. Lett.

100

,

176403

(

2008

).

https://doi.org/10.1103/physrevlett.100.176403

Google Scholar

Crossref

PubMed

65.

W.

Dou

,

A.

Nitzan

, and

J. E.

Subotnik

,

J. Chem. Phys.

143

,

054103

(

2015

).

https://doi.org/10.1063/1.4935713

Google Scholar

Crossref

66.

S. L.

Rudge

,

Y.

Ke

, and

M.

Thoss

,

Phys. Rev. B

107

,

115416

(

2023

).

https://doi.org/10.1103/physrevb.107.115416

Google Scholar

Crossref

67.

S. L.

Rudge

,

C.

Kaspar

,

R. L.

Grether

,

S.

Wolf

,

G.

Stock

, and

M.

Thoss

,

J. Chem. Phys.

160

,

184106

(

2024

).

https://doi.org/10.1063/5.0204307

Google Scholar

Crossref

68.

A.

Erpenbeck

and

M.

Thoss

,

J. Chem. Phys.

151

,

191101

(

2019

).

https://doi.org/10.1063/1.5128206

Google Scholar

Crossref

69.

A.

Erpenbeck

,

Y.

Ke

,

U.

Peskin

, and

M.

Thoss

,

Nanoscale

15

,

16333

(

2023

).

https://doi.org/10.1039/d3nr02176a

Google Scholar

Crossref

PubMed

70.

L. Y.

Gorelik

,

A.

Isacsson

,

M. V.

Voinova

,

B.

Kasemo

,

R. I.

Shekhter

, and

M.

Jonson

,

Phys. Rev. Lett.

80

,

4526

(

1998

).

https://doi.org/10.1103/physrevlett.80.4526

Google Scholar

Crossref

71.

W.

Lai

,

C.

Zhang

, and

Z.

Ma

,

Front. Phys.

10

,

59

(

2015

).

https://doi.org/10.1007/s11467-014-0443-z

Google Scholar

Crossref

72.

A.

Erpenbeck

,

R.

Härtle

,

M.

Bockstedte

, and

M.

Thoss

,

Phys. Rev. B

93

,

115421

(

2016

).

https://doi.org/10.1103/physrevb.93.115421

Google Scholar

Crossref

73.

X.

Dan

,

M.

Xu

,

J. T.

Stockburger

,

J.

Ankerhold

, and

Q.

Shi

,

Phys. Rev. B

107

,

195429

(

2023

).

https://doi.org/10.1103/physrevb.107.195429

Google Scholar

Crossref

74.

T.

Ozaki

,

Phys. Rev. B

75

,

035123

(

2007

).

https://doi.org/10.1103/physrevb.75.035123

Google Scholar

Crossref

75.

Y.

Tanimura

,

J. Phys. Soc. Jpn.

75

,

082001

(

2006

).

https://doi.org/10.1143/jpsj.75.082001

Google Scholar

Crossref

76.

L.

Ye

,

X.

Wang

,

D.

Hou

,

R.-X.

Xu

,

X.

Zheng

, and

Y.

Yan

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

6

,

608

(

2016

).

https://doi.org/10.1002/wcms.1269

Google Scholar

Crossref

77.

Y.

Tanimura

,

J. Chem. Phys.

153

,

020901

(

2020

).

https://doi.org/10.1063/5.0011599

Google Scholar

Crossref

78.

R.

Kapral

and

G.

Ciccotti

,

J. Chem. Phys.

110

,

8919

(

1999

).

https://doi.org/10.1063/1.478811

Google Scholar

Crossref

79.

E.

Wigner

,

Phys. Rev.

40

,

749

(

1932

).

https://doi.org/10.1103/physrev.40.749

Google Scholar

Crossref

80.

J. E.

Moyal

,

Math. Proc. Cambridge Philos. Soc.

45

,

99

(

1949

).

https://doi.org/10.1017/s0305004100000487

Google Scholar

Crossref

81.

K.

İmre

,

E.

Özizmir

,

M.

Rosenbaum

, and

P. F.

Zweifel

,

J. Math. Phys.

8

,

1097

(

1967

).

https://doi.org/10.1063/1.1705323

Google Scholar

Crossref

82.

G.

Stock

and

M.

Thoss

,

Classical Description of Nonadiabatic Quantum Dynamics

,

Advances in Chemical Physics

(

John Wiley and Sons, Ltd.

,

2005

), Chap. 5, pp.

243

–

375

.

Google Scholar

83.

R. J.

Preston

and

D. S.

Kosov

,

J. Chem. Phys.

158

,

224106

(

2023

).

https://doi.org/10.1063/5.0153000

Google Scholar

Crossref

84.

A.

Einstein

,

Ann. Phys.

322

,

549

(

1905

).

https://doi.org/10.1002/andp.19053220806

Google Scholar

Crossref

Nonadiabatic dynamics of molecules interacting with metal surfaces: Extending the hierarchical equations of motion and Langevin dynamics approach to position-dependent metal–molecule couplings

I. INTRODUCTION

II. MODEL

III. QUANTUM TRANSPORT THEORY

IV. ELECTRONIC FRICTION AND LANGEVIN DYNAMICS

A. Partial Wigner transform

B. Wigner-transformed HEOM

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

V. RESULTS

A. Simple vibronic model

B. Interacting vibronic model

VI. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

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

APPENDIX C: DERIVATION OF THE FOKKER–PLANCK EQUATION

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Nonadiabatic dynamics of molecules interacting with metal surfaces: Extending the hierarchical equations of motion and Langevin dynamics approach to position-dependent metal–molecule couplings

I. INTRODUCTION

II. MODEL

III. QUANTUM TRANSPORT THEORY

IV. ELECTRONIC FRICTION AND LANGEVIN DYNAMICS

A. Partial Wigner transform

B. Wigner-transformed HEOM

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

V. RESULTS

A. Simple vibronic model

B. Interacting vibronic model

VI. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

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

APPENDIX C: DERIVATION OF THE FOKKER–PLANCK EQUATION

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only