#  Unraveling current-induced dissociation mechanisms in single-molecule
junctions

Yaling Ke  Institute of Physics, Albert-Ludwig University Freiburg, Hermann-
Herder-Strasse 3, 79104 Freiburg, Germany  André Erpenbeck  School of
Chemistry, The Raymond and Beverley Sackler Center for Computational Molecular
and Materials Science, Tel Aviv University, Tel Aviv 6997801, Israel  Uri
Peskin  Schulich Faculty of Chemistry, Technion-Israel Institute of
Technology, Haifa 32000, Israel  Michael Thoss  Institute of Physics, Albert-
Ludwig University Freiburg, Hermann-Herder-Strasse 3, 79104 Freiburg, Germany
EUCOR Centre for Quantum Science and Quantum Computing, Albert-Ludwig
University Freiburg, Hermann-Herder-Strasse 3, 79104 Freiburg, Germany

###### Abstract

Understanding current-induced bond rupture in single-molecule junctions is
both of fundamental interest and a prerequisite for the design of molecular
junctions, which are stable at higher bias voltages. In this work, we use a
fully quantum mechanical method based on the hierarchical quantum master
equation approach to analyze the dissociation mechanisms in molecular
junctions. Considering a wide range of transport regimes, from off-resonant to
resonant, non-adiabatic to adiabatic transport, and weak to strong vibronic
coupling, our systematic study identifies three dissociation mechanisms. In
the weak and intermediate vibronic coupling regime, the dominant dissociation
mechanism is stepwise vibrational ladder climbing. For strong vibronic
coupling, dissociation is induced via multi-quantum vibrational excitations
triggered either by a single electronic transition at high bias voltages or by
multiple electronic transitions at low biases. Furthermore, the influence of
vibrational relaxation on the dissociation dynamics is analyzed and strategies
for improving the stability of molecular junctions are discussed.

##  I Introduction

Current-induced rupture of chemical bonds is a major concern when single
molecules are being considered as electronic components in nano-scale devices.
The most widely studied architecture in this context is a molecular junction,
where a single molecule is bound to metal or semiconductor electrodes.
Molecular junctions represent a unique architecture to investigate molecules
in a distinct nonequilibrium situation and, in a broader context, to study
basic mechanisms of charge and energy transport in a many-body quantum system
at the nanoscale. Cuevas and Scheer (2010); Galperin, Ratner, and Nitzan
(2007); Bergfield and Ratner (2013); Aradhya and Venkataraman (2013); Bâldea
(2016); Su  _et al._ (2016); Thoss and Evers (2018); Evers  _et al._ (2020)

An important mechanism of bond rupture in molecular junctions is the coupling
of the transport electrons to the vibrations of the molecule, which gives rise
to current-induced vibrational excitation, i.e. local heating. While the level
of current-induced vibrational excitation is typically small for low voltages
in the off-resonant transport regime, it can be substantial for higher
voltages, in particular in the resonant transport regime. In that regime,
current-induced heating can cause mechanical instability of the junction and
may eventually result in bond rupture.Persson and Avouris (1997); Kim, Komeda,
and Kawai (2002); Koch  _et al._ (2006); Huang  _et al._ (2006, 2007); Schulze
_et al._ (2008); Ioffe  _et al._ (2008); Sabater, Untiedt, and van Ruitenbeek
(2015); Li  _et al._ (2015, 2016); Capozzi  _et al._ (2016); Schinabeck
(2018); Gelbwaser-Klimovsky  _et al._ (2018); Bi  _et al._ (2020); Peiris  _et
al._ (2020)

The process of current-induced bond rupture has recently been observed
experimentally in molecular junctions.Sabater, Untiedt, and van Ruitenbeek
(2015); Li  _et al._ (2015, 2016); Capozzi  _et al._ (2016) It is also known
from scanning tunneling microscopy studies of molecules at surfaces.Ho (2002);
Stipe  _et al._ (1997); Huang  _et al._ (2013) The understanding of the
underlying mechanisms of bond rupture and its implication for the stability in
molecular junctions is not only of fundamental interest, but is also crucial
for the design of molecular junctions, which are stable at higher
voltages.Gelbwaser-Klimovsky  _et al._ (2018); Härtle  _et al._ (2018);
Kuperman, Nagar, and Peskin (2020) Molecular junctions that are stable at
higher bias voltages are particularly relevant for possible nanoelectronic
applications. Furthermore, the understanding of current-induced bond rupture
is also crucial for current-induced chemistry and nano-scale chemical
catalysis.Li and Somorjai (2010); Kolasinski (2012); Seideman (2016)

The theoretical framework to study current-induced vibrational excitation in
molecular junctions is well established for models, which treat the
vibrational modes within the harmonic approximation. Galperin, Nitzan, and
Ratner (2006); Ryndyk, Hartung, and Cuniberti (2006); Benesch  _et al._
(2008); Härtle and Thoss (2011); Schinabeck, Härtle, and Thoss (2018);
Erpenbeck  _et al._ (2016) While such models have been used to investigate the
mechanical stability of molecular junctions,Härtle and Thoss (2011); Härtle
and Kulkarni (2015); Schinabeck, Härtle, and Thoss (2018) the study of bond
rupture requires to go beyond the harmonic approximation and use nuclear
potentials which can describe the dissociation process explicitly. This has
been achieved within a classical treatment of the nuclei based on the
Ehrenfest approachDzhioev and Kosov (2011); Dzhioev, Kosov, and Von Oppen
(2013); Pozner, Lifshitz, and Peskin (2014); Erpenbeck  _et al._ (2018a) or
using perturbative theories.Koch  _et al._ (2006); Foti and Vázquez (2018)
Similarly, the study of mechanical instabilities of molecular junctions under
the influence of non-conservative current-induced forces has so far been based
on classical treatments of the nuclei and/or used the harmonic approximation
for the description of the nuclear potentials.Lu, Brandbyge, and Hedegård
(2010); Lü, Hedegård, and Brandbyge (2011); Lü  _et al._ (2012); Preston,
Kershaw, and Kosov (2020); Preston, Gelin, and Kosov (2021) It is also noted
that the theoretical framework to study the related process of dissociative
electron attachment in the gas phase is well established, Domcke (1991);
Gertitschke and Domcke (1993); Čížek, Horáček, and Domcke (1999); Gallup and
Fabrikant (2011) but this problem is conceptually simpler because only a
single electron that is scattered from the molecule has to be considered.
Moreover, the processes of light-induced dissociation or desorption of
molecules at surfaces has also been studied in great detail
theoretically.Brandbyge  _et al._ (1995); Saalfrank (2006); Ho (2002); Kim
_et al._ (2015); Frederiksen, Paulsson, and Ueba (2014) In that scenario,
however, typically the system is only temporarily driven out of equilibrium by
a laser pulse, while in molecular junction transport, the electrical current
driven through the molecule results in a nonequilibrium steady state.

Recently, we have developed a fully quantum mechanical theoretical framework
to study current-induced bond rupture in molecular junctions, which takes
proper account of the many-electron, nonequilibrium nature of the transport
process and is not limited to weak coupling.Erpenbeck  _et al._ (2018a, 2020)
The method combines the hierarchical quantum master equation (HQME) approach
with a discrete variable representation (DVR) of the nuclear degrees of
freedom to facilitate the description of general potential energy surfaces
(PESs). The application to a model where the charged state is characterized by
a repulsive potential showed that the current-induced population of anti-
bonding states is another important mechanism, which can lead to fast bond
rupture in molecular junctions and can dominate over current-induced
heating.Erpenbeck  _et al._ (2020)

In the present work, we extend our previous study and consider a scenario,
where the potential energy surface of the charged molecule also supports bound
states. In this model, both bond-rupture via direct dissociation in the
continuum states of the charged state potential energy surface or via current-
induced heating are possible. The model also accounts for additional
dissociation channels via Feshbach resonances. Furthermore, it incorporates
vibrational relaxation processes induced by coupling of the dissociative
reaction mode to other inactive modes (intramolecular vibrational relaxation),
the phonons of the leads or a possible solution environment. The detailed
analysis of this extended model provides a rather comprehensive understanding
of the mechanisms of current-induced bond rupture in molecular junctions.

The remainder of the paper is organized as follows: In Sec. II, we introduce
the model and the HQME approach used to investigate current-induced reaction
dynamics. The results are presented and discussed in detail in Sec. III,
considering a broad range of different regimes and processes, comprising off-
resonant to resonant transport, weak to strong vibronic and molecule-lead
coupling, as well as vibrational relaxation due to coupling to a phonon bath.
Furthermore, time-dependent current-voltage characteristics are presented,
implications for experiments are addressed, and strategies for improving the
stability of molecular junctions are discussed. We close with a summary in
Sec. IV.

##  II Theory

###  II.1 Model

Figure 1:  Visualization of the molecular junction and the system energy
landscape. The molecular junction consists of a backbone coupled to electrodes
and a side-group. In case of dissociation, the side group detaches from the
backbone. The neutral [charged] state of the molecule is characterized by a
Morse potential $V_{g}(Q)$ [$V_{e}(Q)$] with equilibrium position $Q_{g}$
[$Q_{e}=Q_{g}+\Delta Q$] and fundamental oscillation frequency $\Omega_{0}$.
$E_{\rm C}$ denotes the charging energy, $E_{\rm D}$ the dissociation energy,
and $g_{\rm L/R}(Q)$ the coordinate-dependent molecule-lead coupling. The
energy levels of the bound states in the potentials $V_{g}(Q)$ $[V_{e}(Q)]$
are indicated by horizontal blue [red] lines. The gray dotted line and shaded
area indicate the complex absorbing potential $W(Q)$. The probability absorbed
by the complex absorbing potential is mapped to the representative auxiliary
grid point $Q_{\infty}$. The coupling to a phonon bath is depicted by black
wavy lines within the orange shaded round area.

We consider a molecular junction depicted in Fig. 1, where the molecule,
consisting of a backbone and a side group, is in contact with two macroscopic
electrodes and a phonon bath, which is described by the Hamiltonian

$H=H_{\rm mol}+H_{\rm leads}+H_{\rm mol-leads}+H_{\rm ph}+H_{\rm mol-ph}.$ (1)

For the molecule, a minimal model is adopted comprising a single vibrational
reaction mode, which describes the bonding of the side group, and a spinless
electronic level. The molecular Hamiltonian takes the form

$H_{\rm mol}=\frac{P^{2}}{2M}+V_{g}(Q)dd^{\dagger}+V_{e}(Q)d^{\dagger}d,$ (2)

where $Q$ and $P$ are the coordinate and momentum of the reaction mode,
respectively, and $M$ denotes the corresponding reduced nuclear mass. The
operator $d^{\dagger}$ creates an electron in the molecular electronic level,
and $d$ is its hermitian conjugate. $V_{g}(Q)$ and $V_{e}(Q)$ describe the
potential energy surfaces of the electronic ground state of the neutral and
charged molecule, respectively.

Specifically, the neutral state is described by a Morse potential,

$V_{g}(Q)=E_{\rm D}(1-e^{-a(Q-Q_{g})})^{2},$ (3)

where $Q_{g}$ denotes the equilibrium position, $E_{\rm D}$ the dissociation
energy, and $a=1.028a_{0}^{-1}$ the width parameter of the Morse potential
with $a_{0}$ being the Bohr radius. In the calculations reported below, the
parameters are chosen as $Q_{g}=1.78\textrm{ \AA}$, $a=1.028a_{0}^{-1}$,
$E_{\rm D}=2.38$ eV, and $M=1$ amu (atomic mass unit). The corresponding
fundamental frequency for small oscillations at the bottom of potential well
is $\hbar\Omega_{0}=a\sqrt{2E_{\rm D}/M}=274$ meV, which is within the typical
range of molecular vibrations. For such a high frequency, the nuclear quantum
effects are expected to be non-negligible even at the room temperature. There
are in total 16 bound vibrational states. The choice of these parameters is
motivated by the process of H2 desorption from metal surfaces.Halstead and
Holloway (1990) However, it is emphasized that the goal of this work is to
study the basic mechanisms of current-induced bond rupture and it does not
attempt to describe a specific molecule.

The potential energy surface of the charged state is assumed to be of the same
form as in the neutral state but with a shifted equilibrium position
$Q_{e}=Q_{g}+\Delta Q$,

$V_{e}(Q)=E_{\rm D}(1-e^{-a(Q-Q_{e})})^{2}+E_{\rm C}.$ (4)

Here, $E_{\rm C}$ denotes the charging energy, which is chosen as $E_{\rm
C}=1$ eV in the calculations reported below. The displacement $\Delta Q$
determines the electronic-vibrational (vibronic) coupling strength. It should
be emphasized that different from some other models of vibrationally-coupled
transport in molecular junctions using the harmonic approximation,Galperin,
Nitzan, and Ratner (2006); Härtle and Thoss (2011); Schinabeck  _et al._
(2016) the charging energy in the model employed here is independent of
$\Delta Q$. The potential energy surfaces are illustrated schematically in
Fig. 1.

The molecule is coupled to two leads ($\alpha=L/R$), which serve as electron
reservoirs and are modelled by non-interacting electrons,

$H_{\rm leads}=\sum_{\alpha}H_{\alpha}=\sum_{\alpha}\sum_{k}\epsilon_{\alpha
k}c_{\alpha k}^{\dagger}c_{\alpha k},$ (5)

where $c_{\alpha k}^{\dagger}(c_{\alpha k})$ creates (annihilates) an electron
in the $k$th state with energy $\epsilon_{\alpha k}$ in lead $\alpha$.

The molecule-lead interaction is described by

$H_{\rm mol-leads}=\sum_{\alpha}\sum_{k}g_{\alpha}(Q)(t_{\alpha k}c_{\alpha
k}^{\dagger}d+t^{*}_{\alpha k}d^{\dagger}c_{\alpha k}).$ (6)

Here, $t_{\alpha k}$ denotes the coupling strength between the electronic
state at the molecule and the $k$th state in lead $\alpha$. The dependence of
the molecule-lead interaction on the nuclear coordinate $Q$ allows the
modelling of situations where the conductance changes upon detachment of the
side group. This can, for example, result from a change of the
$\pi$-conjugation within the molecular backbone as a consequence of the side
group separating from the molecule. In this work, we employ a dependence on
the nuclear coordinate of the form

$g_{\rm
L/R}(Q)=\frac{1-q}{2}\left[1-\tanh\left(2(Q-\tilde{Q})/a_{0}\right)\right]+q,$
(7)

which is depicted by the green line in Fig. 1. The parameter $\tilde{Q}$
determines the region of the reaction mode, where the transition between
stronger and weaker molecule-lead coupling upon detachment of the side group
occurs, and is set below to $\tilde{Q}=4.0\textrm{~{}\AA}$. In this work, a
non-destructive dissociation is considered, that is, the molecule is still
conductive after the dissociation of the side group. The parameter $q$
determines the relative molecule-lead coupling strength for the dissociated
molecule and is chosen as $q=0.05$.

To model vibrational relaxation of the dissociative reaction mode induced by
coupling to other inactive modes (intramolecular vibrational relaxation), the
phonons of the leads or a possible solution environment, we include the
coupling of the reaction mode to a bosonic bath (in the following referred to
as phonon bath). The phonon bath is modelled by a collection of harmonic
oscillators,

$H_{\rm
ph}=\sum_{j}\frac{p_{j}^{2}}{2m_{j}}+\frac{m_{j}\omega_{j}^{2}q_{j}^{2}}{2},$
(8)

where $\omega_{j}$ denotes the frequency of the $j$th bath oscillator;
$q_{j}$, $p_{j}$, and $m_{j}$ are the corresponding coordinate, momentum, and
mass, respectively. The interaction between the reaction mode and the phonon
bath is given by

$H_{\rm mol-
ph}=-f(Q)\sum_{j}h_{j}q_{j}+\sum_{j}\frac{\left(h_{j}f(Q)\right)^{2}}{2m_{j}\omega_{j}^{2}}.$
(9)

Here, $h_{j}$ denotes the coupling strength between the $j$th bath oscillator
and the reaction mode. The second term in Eq. (9) is a counter term, which is
introduced to avoid an extensive effect of the bath on the potential of the
system.Weiss (2012) The coupling operator $f(Q)$ is taken in the following
form

$\begin{split}f(Q)=&\frac{(Q-Q_{g})}{a_{0}}e^{-\chi\left(\frac{Q-Q_{g}}{a_{0}}\right)^{2}}dd^{\dagger}\\\
&+\frac{(Q-Q_{e})}{a_{0}}e^{-\chi\left(\frac{Q-Q_{e}}{a_{0}}\right)^{2}}d^{\dagger}d,\end{split}$
(10)

where $\chi=\frac{\hbar\Omega_{0}}{4E_{\rm D}}$ denotes the anharmonicity of
the Morse potential.Ilk and Makri (1994) This specific form ensures that the
coupling of the reaction mode to the phonon bath vanishes for large values of
$Q$. In the harmonic limit, $\chi\rightarrow 0$, it reduces to the
paradigmatic linear form.Joutsuka and Ando (2011)

###  II.2 Method

To simulate the nonequilibrium dynamics of the molecular junction based on the
Hamiltonian in Eq. (1), we use the HQME method. The HQME approach, also
referred to as hierarchical equations of motion (HEOM), is a reduced density
matrix scheme, which describes the dynamics of a quantum system influenced by
an environment. In the case considered here, the molecule is the system and
the leads together with the phonon bath represent the environment. The HQME
approach generalizes perturbative quantum master equation methods by including
higher-order contributions as well as non-Markovian memory and allows for the
systematic convergence of the results. Härtle  _et al._ (2013, 2015); Xu  _et
al._ (2017); Trushechkin (2019) The method was originally developed by
Tanimura and Kubo to study relaxation dynamics. Tanimura and Kubo (1989);
Tanimura (2006) Later it was extended to fermionic charge transport,Jin  _et
al._ (2007); Jin, Zheng, and Yan (2008); Zheng  _et al._ (2009); Yan (2014);
Ye  _et al._ (2016); Härtle  _et al._ (2013, 2015); Wenderoth, Bätge, and
Härtle (2016) also including electronic-vibrational coupling.Schinabeck  _et
al._ (2016); Schinabeck, Härtle, and Thoss (2018); Dou  _et al._ (2018) For
more details about the developments and applications of the method, we refer
to the review in Ref. Tanimura, 2020.

In recent work, we have formulated and applied the HQME method to simulate
bond rupture in molecular junctions.Erpenbeck and Thoss (2019); Erpenbeck  _et
al._ (2020) Moreover, we have formulated the HQME method for an open quantum
system, which is coupled to multiple fermionic and bosonic environments, as is
the case in the current application of the method.Bätge  _et al._ (2021) In
the following, a brief recapitulation of the most important aspects of the
method is given.

We assume that the initial state is given by

$\rho(t=0)=\rho_{\rm s}(0)\rho_{\rm leads}^{\rm eq}\rho_{\rm ph}^{\rm eq},$
(11)

including the initial density matrices of the system, $\rho_{\rm s}(0)$, the
leads, $\rho_{\rm leads}^{\rm eq}$, and the phonon bath, $\rho_{\rm ph}^{\rm
eq}$, respectively. The latter are assumed to be in their respective
equilibrium state. Specifically, the leads are initially described by the
grand canonical distribution,

$\rho_{\rm leads}^{\rm eq}=\prod_{\alpha}\rho_{\alpha}^{\rm
eq}=\prod_{\alpha}\frac{e^{-\beta_{\alpha}(H_{\alpha}-\mu_{\alpha}N_{\alpha})}}{\mathrm{Tr}_{\alpha}\\{e^{-\beta_{\alpha}(H_{\alpha}-\mu_{\alpha}N_{\alpha})}\\}},$
(12)

where $\beta_{\alpha}=1/(k_{\rm B}T_{\alpha})$ denotes the inverse temperature
with Boltzmann constant $k_{\rm B}$, $\mu_{\alpha}$ is the chemical potential
and $N_{\alpha}=\sum_{k}c_{\alpha k}^{\dagger}c_{\alpha k}$ the occupation
number operator of lead $\alpha$, respectively. Moreover,
$\mathrm{Tr}_{\alpha}$ denotes the trace over all electronic degrees of
freedom in lead $\alpha$. The difference of the chemical potentials of the
left and right leads defines the bias voltage $\Phi$, which we assume to drop
symmetrically, i.e., $\mu_{\rm L}=-\mu_{\rm R}=e\Phi/2$. The initial
equilibrium state of the phonon bath at inverse temperature $\beta_{\rm
ph}=1/(k_{\rm B}T_{\rm ph})$ is given by

$\rho_{\rm ph}^{\rm eq}=\frac{e^{-\beta_{\rm ph}H_{\rm ph}}}{\mathrm{Tr}_{\rm
ph}\\{e^{-\beta_{\rm ph}H_{\rm ph}}\\}}.$ (13)

In the studies below, the temperatures of the leads and the phonon bath are
set to $T_{\rm L}=T_{\rm R}=T_{\rm ph}=300$ K.

Due to the Gaussian statistical properties of the electronic reservoirs and
the phonon bath, the influence of the environments on the system dynamics is
encoded in the two-time correlation functions $C_{\alpha}^{\pm}(t-\tau)$ and
$C_{\rm ph}(t-\tau)$ of the lead electrons and the bath phonons, respectively.

The two-time correlation function of the lead electrons is given by

$\begin{split}C_{\alpha}^{\sigma}(t-\tau)\equiv&\mathrm{Tr}_{\alpha}\left\\{\sum_{k}\tilde{c}_{\alpha
k}^{\sigma}(t)\tilde{c}_{\alpha k}^{\bar{\sigma}}(\tau)\rho^{\rm
eq}_{\rm\alpha}\right\\}\\\
=&\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{i\sigma\frac{\epsilon}{\hbar}(t-\tau)}f_{\alpha}^{\sigma}(\epsilon)\Gamma_{\alpha}(\epsilon)\mathrm{d}\epsilon,\end{split}$
(14)

where we have introduced the operators

$\tilde{c}_{\alpha
k}^{\sigma}(t)=e^{\frac{it}{\hbar}H_{\alpha}}t^{\sigma}_{\alpha k}c_{\alpha
k}^{\sigma}e^{-\frac{it}{\hbar}H_{\alpha}}$ (15)

and the notation $\sigma=\pm$, $\bar{\sigma}=-\sigma$, $c^{+}_{\alpha
k}=c^{\dagger}_{\alpha k}$, $c^{-}_{\alpha k}=c_{\alpha k}$, $t^{+}_{\alpha
k}=t^{*}_{\alpha k}$, $t^{-}_{\alpha k}=t_{\alpha k}$. It includes the Fermi
distribution function
$f_{\alpha}^{\sigma}(\epsilon)=1/(1+e^{\sigma\beta_{\alpha}(\epsilon-\mu_{\alpha})})$
and the coupling-weighted density of states of lead $\alpha$ (also called
level-width function),

$\Gamma_{\alpha}(\epsilon)=2\pi\sum_{k}\left|t_{\alpha
k}\right|^{2}\delta(\epsilon-\epsilon_{\alpha k}).$ (16)

For the scope of this work, the leads are described within the wide-band
limit, where the level-width function is energy-independent,
$\Gamma_{\alpha}=2\pi\left|t_{\alpha}\right|^{2}$.

The influence of the phonon bath on the reduced system dynamics is determined
by the correlation function

$\begin{split}C_{\rm ph}(t-\tau)=&\mathrm{Tr}_{\rm
ph}\left\\{\sum_{j}\tilde{q}_{j}(t)\tilde{q}_{j}(\tau)\rho^{\rm eq}_{\rm
ph}\right\\}\\\ =&\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{d}\omega
e^{-i\omega(t-\tau)}J(\omega)f_{\rm ph}(\omega),\end{split}$ (17)

where we have introduced the operators

$\tilde{q}_{j}(t)=e^{iH_{\rm ph}t/\hbar}h_{j}q_{j}e^{-iH_{\rm ph}t/\hbar},$
(18)

the Bose distribution function $f_{\rm ph}(\omega)=1/(e^{\beta_{\rm
ph}\hbar\omega}-1)$, and the spectral density of the phonon bath,

$J(\omega)=2\pi\hbar\sum_{j}\frac{h_{j}^{2}}{2m_{j}\omega_{j}}\delta(\omega-\omega_{j}).$
(19)

In this work, we assume that the spectral density of the phonon bath is of
Lorentz-Drude form,

$J(\omega)=2\hbar\lambda_{\rm
ph}\frac{\omega\omega_{c}}{\omega^{2}+\omega_{c}^{2}},$ (20)

where $\lambda_{\rm ph}$ characterizes the coupling strength and $\omega_{c}$
is the characteristic frequency, which is the inverse of the bath correlation
time $\tau_{c}$.

To obtain a formally closed set of hierarchical equations, a central step is
to represent the two-time correlation functions as sums of exponential
functions. To this end, various sum-over-poles decomposition schemes for the
Fermi and Bose distribution functions Hu, Xu, and Yan (2010); Hu  _et al._
(2011); Cui  _et al._ (2019); Abe, Yamashita, and Saalfrank (2003) or the
corresponding time correlation functions Tang  _et al._ (2015); Rahman and
Kleinekathöfer (2019); Erpenbeck  _et al._ (2018b) have been proposed. Here,
we employ the Padé decomposition scheme, which was found to be particularly
efficient at moderate temperatures. In this way, the correlation function in
Eq. (14) can be represented as

$C_{\alpha}^{\sigma}(t-\tau)\approx\hbar\pi\delta(t-\tau)+\sum_{l=1}^{\rm
L}\eta_{\alpha,l}e^{-\gamma_{\alpha,l}^{\sigma}(t-\tau)},$ (21)

with exponents
$\hbar\gamma_{\alpha,l}^{\sigma}=\xi_{l}/\beta_{\alpha}-i\sigma\mu_{\alpha}$
and coefficients $\eta_{\alpha,l}=-2i\pi\kappa_{l}/\beta_{\alpha}$. The Padé
parameters $\kappa_{l}$ and $\xi_{l}$ are obtained using the procedures
described in Ref. Hu  _et al._ , 2011. The delta-function term appears due to
the wide-band limit approximation.

Likewise, utilizing the Padé spectrum decomposition of the Bose distribution
function, the time-correlation function for the phonon bath is written as

$C_{\rm ph}(t-\tau)\approx\sum_{k=0}^{\rm
K}\eta_{k}e^{-\gamma_{k}(t-\tau)}+\sum_{k={\rm
K}+1}^{\infty}\frac{2\eta_{k}}{\gamma_{k}}\delta(t-\tau),$ (22)

with $\hbar\gamma_{0}=\omega_{c}$ and $\eta_{0}=\hbar\lambda_{\rm
ph}\omega_{c}\cot(\beta\hbar\omega_{c}/2)-i\hbar\lambda_{\rm ph}\omega_{c}$ as
well as $\hbar\gamma_{k}=\nu_{k}$ and $\eta_{k}=\frac{4\lambda_{\rm
ph}\xi_{k}}{\beta}\frac{\omega_{c}\nu_{k}}{\nu_{k}^{2}-\omega_{c}^{2}}$ for
$k>0$. The Padé parameters $\xi_{k}$ and $\nu_{k}$ for the Bose function are
all real numbers. In the above formula, it is assumed
$\gamma_{k}e^{-\gamma_{k}|t|}\approx 2\delta(t)$ when $\gamma_{k}$ is
significantly larger than the vibrational frequency of the system
$\Omega_{0}$.Ishizaki and Tanimura (2005)

Based on the exponential expansions of the correlation functions in Eqs. (21)
and (22), the HQME method uses a set of auxiliary density operators
$\rho_{\bm{m}}^{\bm{n}}$ to describe the dynamics of the open quantum
system.Jin, Zheng, and Yan (2008); Shi  _et al._ (2009); Härtle  _et al._
(2013); Schinabeck, Härtle, and Thoss (2018); Bätge  _et al._ (2021) The
subscript $\bm{m}$ is a bosonic index vector,
$\bm{m}=(m_{0},m_{1},\cdots,m_{K})$, where every element $m_{k}$ is a non-
negative integer. The superscript $\bm{n}$ denotes an ordered array of
fermionic multi-indices $\bm{a}_{j}=(\alpha_{j},\sigma_{j},l_{j})$, i.e.,
$\bm{n}=(\bm{a}_{1},\cdots,\bm{a}_{n})$, where $\alpha_{j}$ labels the left or
right lead, $\sigma_{j}=\pm$, and $l_{j}$ specifies the fermionic Padé pole.
The reduced density operator of the system, which is defined by tracing the
density operator of the complete system $\rho(t)$ over the environmental DOFs,

$\rho_{\rm s}(t)=\mathrm{Tr}_{{\rm leads+ph}}\\{\rho(t)\\},$ (23)

corresponds to the lowest tier auxiliary density operator, i.e., $\rho_{\rm
s}(t)=\rho_{0}^{0}(t)$.

The HQME is then given by the following equations of motion of the auxiliary
density operators,Jin, Zheng, and Yan (2008); Shi  _et al._ (2009); Bätge  _et
al._ (2021); Xu  _et al._ (2019)

$\begin{split}\frac{\partial\rho_{\bm{m}}^{\bm{n}}(t)}{\partial
t}=&-\frac{i}{\hbar}\left[H_{\rm
mol},\rho_{\bm{m}}^{\bm{n}}(t)\right]_{-}-\frac{1}{\hbar}\mathcal{L}_{\infty}\rho_{\bm{m}}^{\bm{n}}(t)-\left(\sum_{j=1}^{n}\gamma_{a_{j}}+\sum_{k=0}^{K}m_{k}\gamma_{k}\right)\rho_{\bm{m}}^{\bm{n}}(t)\\\
&-\sum_{\alpha\in{\rm
L/R},\sigma\in{\pm}}\frac{\pi}{2\hbar}\left[g_{\alpha}(Q)d^{\overline{\sigma}},\left[g_{\alpha}(Q)d^{\sigma},\rho_{\bm{m}}^{\bm{n}}(t)\right]_{(-)^{n+1}}\right]_{(-)^{n+1}}-\frac{i}{\hbar}\sum_{j=1}^{n}(-1)^{n-j}\mathcal{C}_{a_{j}}\rho_{\bm{m}}^{\bm{n}^{-}_{j}}(t)-\frac{i}{\hbar}\sum_{a_{n+1}}\mathcal{A}_{\alpha_{n+1}}^{\overline{\sigma}_{n+1}}\rho_{\bm{m}}^{\bm{n}^{+}}(t)\\\
&-\frac{\tilde{C}_{K}}{\hbar^{2}}\left[f(Q),\left[f(Q),\rho_{\bm{m}}^{\bm{n}}(t)\right]_{-}\right]_{-}-\frac{i}{\hbar}\sum_{k=0}^{K}\sqrt{(m_{k}+1)\left|\eta_{k}\right|}\mathcal{A}_{k}\rho_{\bm{m}_{k}^{+}}^{\bm{n}}(t)-\frac{i}{\hbar}\sum_{k=0}^{K}\sqrt{\frac{m_{k}}{\left|\eta_{k}\right|}}\mathcal{C}_{k}\rho_{\bm{m}_{k}^{-}}^{\bm{n}}(t).\end{split}$
(24)

Here,

$\tilde{C}_{\rm K}=\sum_{k={\rm
K}+1}^{\infty}\frac{\eta_{k}}{\gamma_{k}}=\frac{2\lambda}{\beta\omega_{c}}-\sum_{k=0}^{K}\frac{\eta_{k}}{\gamma_{k}}$
(25)

and a few shorthand notation are used, including the bosonic index arrays
$\bm{m}_{k}^{+}=(m_{0},\cdots,m_{k}+1,\cdots m_{K})$ and
$\bm{m}_{k}^{-}=(m_{0},\cdots,m_{k}-1,\cdots,m_{K})$, the fermionic index
arrays
$\bm{n}_{j}^{-}=(\bm{a}_{1},\cdots,\bm{a}_{j-1},\bm{a}_{j+1},\cdots,\bm{a}_{n})$
and $\bm{n}^{+}=(\bm{a}_{1},\cdots,\bm{a}_{n},\bm{a}_{n+1})$ as well as the
commutator $[\mathcal{O},\rho_{\bm{m}}^{\bm{n}}]_{-}$ and anticommutator
$[\mathcal{O},\rho_{\bm{m}}^{\bm{n}}]_{+}$. Note that the first term in the
second (third) line of Eq. (24) stems from the delta-function approximation in
Eq. (21) (Eq. (22)). We found that the numerical performance benefits
significantly from these two approximations, because the stiffness of the
equation can be significantly reduced as compared to considering a very large
bandwidth or a very high Padé frequency. The superoperators
$\mathcal{A}_{\alpha}^{\sigma}$, $\mathcal{C}_{a}$, $\mathcal{A}_{k}$, and
$\mathcal{C}_{k}$ connect the different tiers of the hierarchy and are given
by

$\mathcal{A}_{\alpha}^{\sigma}\rho_{\bm{m}}^{\bm{n}}(t)=g_{\alpha}(Q)d^{\sigma}\rho_{\bm{m}}^{\bm{n}}(t)+(-1)^{n}\rho_{\bm{m}}^{\bm{n}}(t)d^{\sigma}g_{\alpha}(Q),$
(26)
$\mathcal{C}_{a}\rho_{\bm{m}}^{\bm{n}}(t)=\eta_{\alpha,l}^{\sigma}g_{\alpha}(Q)d^{\sigma}\rho_{\bm{m}}^{\bm{n}}(t)-(-1)^{n}\eta_{\alpha,l}^{\bar{\sigma}*}\rho_{\bm{m}}^{\bm{n}}(t)d^{\sigma}g_{\alpha}(Q),$
(27)
$\mathcal{A}_{k}\rho_{\bm{m}}^{\bm{n}}(t)=f(Q)\rho_{\bm{m}}^{\bm{n}}(t)-\rho_{\bm{m}}^{\bm{n}}(t)f(Q),$
(28)
$\mathcal{C}_{k}\rho_{\bm{m}}^{\bm{n}}(t)=\eta_{k}f(Q)\rho_{\bm{m}}^{\bm{n}}(t)-\rho_{\bm{m}}^{\bm{n}}(t)f(Q)\eta^{*}_{k}.$
(29)

The numerically exact description provided by the HQME approach employs an
infinite hierarchy of auxiliary density operators, which needs to be truncated
in a suitable manner. For a detailed discussion, we refer to Refs. Tanimura,
2006; Ye  _et al._ , 2016.

Within the HQME formalism outlined above, all (auxiliary) density operators
are also acting on the nuclear DOF of the system, the reaction coordinate. In
order to allow for a description of nuclear DOFs with generic PESs, we employ
a discrete variable representation (DVR),Colbert and Miller (1992); Echave and
Clary (1992); Seideman and Miller (1992) which represents the nuclear reaction
coordinate effectively by a finite set of grid points $Q_{i}$.

In order to avoid finite size effects, a complex absorbing potential (CAP),
$W(Q)$, is introduced, which absorbs the parts of the density reaching the
boundary of the DVR grid.Riss and Meyer (1996) In the calculations reported
below, a power-law form of the CAP is used,Erpenbeck  _et al._ (2020)

$W(Q)=\zeta(Q-Q_{\rm CAP})^{4}\cdot\Theta(Q-Q_{\rm CAP}),$ (30)

with the Heaviside step function $\Theta$. The parameters of the CAP, $Q_{\rm
CAP}=4.0\textrm{ \AA}$ and $\zeta=5$ eV/$\textrm{\AA}^{4}$, were determined by
test calculations to ensure that the observables obtained do not depend on the
CAP.

As discussed previously,Selstø and Kvaal (2010); Kvaal (2011); Prucker  _et
al._ (2018); Erpenbeck and Thoss (2019); Erpenbeck  _et al._ (2020) the
introduction of the CAP may result in problems associated with the
conservation of the particle number. In particular, the action of the CAP
causes a decrease of the trace of the density matrix and consequently an
artificial loss of the number of electrons. To avoid these problems, we
compensate for the loss of populations due to the action of the CAP by
introducing an additional Lindblad-like source term into the HQME,Erpenbeck
and Thoss (2019); Erpenbeck  _et al._ (2020)

$\begin{split}\mathcal{L}_{\infty}\rho_{\bm{m}}^{\bm{n}}(t)=&2C_{\infty}(Q)\rho_{\bm{m}}^{\bm{n}}(t)C^{\dagger}_{\infty}(Q)\\\
&-\left[C^{\dagger}_{\infty}(Q)C_{\infty}(Q),\rho_{\bm{m}}^{\bm{n}}(t)\right]_{+}\end{split}$
(31)

with

$C_{\infty}(Q)=\sqrt{W(Q)}\left|Q_{\infty}\right\rangle\left\langle Q\right|.$
(32)

This source term maps the probability absorbed by the CAP to an auxiliary grid
point $Q_{\infty}$ (see Fig. 1 and Ref. Erpenbeck  _et al._ , 2020).

###  II.3 Observables of interest

We briefly comment on the calculation of observables of interest. Any system
observable can be obtained directly from the reduced density matrix
$\rho_{s}(t)=\rho_{0}^{0}(t)$. Several system observables are considered
below. In particular, we are interested in the dissociation probability, which
can be calculated as

$P_{\infty}(t)=\mathrm{Tr}_{s}\left\\{\left|Q_{\infty}\right\rangle\left\langle
Q_{\infty}\right|\rho_{s}(t)\right\\},$ (33)

where $\mathrm{Tr}_{s}$ denotes the trace over electronic and nuclear DOF of
the system. The remaining population $1-P_{\infty}(t)$, which corresponds to
the portion of non-dissociated molecules, is called the survival probability
in the following. Assuming an exponential kinetics of the dissociation process
in the long-time limit, the dissociation rate is given by

$k_{\mathrm{diss}}=-\lim_{t\rightarrow\infty}\frac{d\ln(1-P_{\infty}(t))}{dt}.$
(34)

Another important observable to characterize the dynamics is the population of
the vibrational states in the neutral and charged state of the molecule, given
by

$P^{g}_{v}=\mathrm{Tr}_{s}\left\\{dd^{\dagger}|\psi_{v}^{g}\rangle\langle\psi^{g}_{v}|\rho_{s}(t)\right\\}$
(35)

and

$P^{e}_{v}=\mathrm{Tr}_{s}\left\\{d^{\dagger}d|\psi_{v}^{e}\rangle\langle\psi^{e}_{v}|\rho_{s}(t)\right\\}.$
(36)

Here, $|\psi^{g/e}_{v}\rangle$ denotes the $v$th vibrational eigenfunction in
the neutral/charged state of the molecule. Due to dissociation, the population
of the undissociated molecule decays over time. In the analysis below,
vibrational populations are renormalized by the survival probability
$(1-P_{\infty}(t))$ to obtain a stationary distribution in the long-time limit
despite the dissociation. The renormalized average vibrational excitation can
be obtained as

$\langle
n_{\mathrm{vib}}\rangle=\sum_{v=0}v\frac{\mathrm{Tr}_{s}\left\\{|\psi_{v}^{g}\rangle\langle\psi^{g}_{v}|\rho_{s}\right\\}}{1-P_{\infty}}.$
(37)

For stable molecular junctions, the notion of an effective temperature can be
introduced to quantify the non-thermal vibrational excitation.Preston,
Kershaw, and Kosov (2020); Wang, Nian, and Lü (2020); Zhang, Zheng, and Di
Ventra (2019)

Bath-related observables of interest can be obtained from the auxiliary
density operators. For example, the electronic current from lead $\alpha$ to
the molecule is expressed in terms of the zeroth and first fermionic tier
auxiliary density operators,

$\begin{split}I_{\alpha}(t)=&-2e\frac{\mathrm{d}}{\mathrm{d}t}\left\langle\sum_{k}c_{\alpha
k}^{\dagger}c_{\alpha k}\right\rangle\\\
=&\frac{2ie}{\hbar}\sum_{l=1}\mathrm{Tr}_{s}\left\\{g_{\alpha}(Q)\left(d\rho^{(\alpha,+,l)}_{\bm{0}}(t)-d^{\dagger}\rho_{\bm{0}}^{(\alpha,-,l)}(t)\right)\right\\}\\\
&+\frac{2\pi
e}{\hbar}\mathrm{Tr}_{s}\left\\{g_{\alpha}(Q)g_{\alpha}(Q)\left[dd^{\dagger}\rho_{\bm{0}}^{0}(t)-d^{\dagger}d\rho_{\bm{0}}^{0}(t)\right]\right\\}\end{split}$
(38)

where $e$ denotes the electron charge and where the spin degeneracy is taken
into account. The current passing through the molecule is given by
$I(t)=(I_{L}(t)-I_{R}(t))/2$. The contribution of the zeroth tier auxiliary
density operator in the above current expression is due to the wide-band
approximation. The details of the derivation are provided in the supplementary
material.

###  II.4 Numerical details

Here, we provide some details of the numerical calculations.

The initial state of the system is chosen as the vibrational ground state of
the neutral molecule,

$\rho_{s}(0)=dd^{\dagger}\left|\psi^{g}_{v=0}\right\rangle\left\langle\psi_{v=0}^{g}\right|.$
(39)

More details about the influence of the initial state on the dynamics can be
found in the supplementary material. In short, the dissociation rate and the
quantities rescaled by the survival probability at long times are insensitive
to the initial preparation.

The HQME, Eq. (24), is solved using the propagation scheme proposed in Ref.
Wilkins and Dattani, 2015, which is based on the power series expansion of the
propagator and gives a more efficient use of memory. The numerical
calculations are performed on GPUs and the shared memory parallel programming
technique is employed for further speed-up. Furthermore, it should be
emphasized that a large percent of auxiliary density matrices are exactly zero
due to the exclusion principle. Therefore, a sparse matrix multiplication
algorithm is used, where the sparsity of all auxiliary density matrices and
the Hamiltonian are checked before propagation.

For all data presented below, we have tested the convergence of the
observables with respect to the number of DVR grid points, the number of Padé
poles used to represent Fermi and Bose function, the time step of the
integrator, and the truncation tier of the HQME. At least 64 DVR grid-points
for the reactive coordinate are required. The calculations employ 20 Padé
poles for the Fermi function and two for the Bose function. We adopt a
hierarchy truncation scheme, where all auxiliary density matrices beyond the
specified truncation tier are set to zero. The converged hierarchical
truncation tier depends on the molecule-lead coupling strength, bias voltage,
vibration-phonon bath interaction. Stronger molecule-lead coupling or a lower
bias voltage requires a higher fermionic hierarchical truncation tier.
Likewise, a stronger phonon bath coupling requires a higher bosonic hierarchy
tier. For most parameter sets chosen in this work, converged results are
obtained at the second or third fermionic and bosonic tier of the hierarchy.

##  III Results and Discussion

In this section, we study the current-induced dissociation dynamics in single-
molecule junctions based on the methods and model introduced above. To provide
a comprehensive analysis of the underlying mechanisms, we consider in Sec.
III.1 \- Sec. III.3 a broad range of different regimes and processes,
comprising off-resonant to resonant transport, weak to strong vibronic and
molecule-lead coupling, as well as vibrational relaxation due to coupling to a
phonon bath. In Sec. III.4, time-dependent current-voltage characteristics are
presented and the implications for experiments are addressed. Furthermore,
strategies for improving the stability of molecular junctions are discussed.

###  III.1 Overview of dissociation mechanisms

As a basis for the subsequent detailed analysis, we first give an overview of
the most important dissociation mechanisms.

Fig. 2 summarizes the basic vibrational heating and cooling processes in a
molecular junction. When current-induced vibrational heating exceeds heat
dissipation, the energy accumulated in a certain bond can reach the
dissociation threshold and the bond ruptures. Depending on the specific bond
affected, the functionality of the junction may be destroyed. The underlying
current-induced dissociation mechanism is determined by the applied bias
voltage and intrinsic molecular properties. The bias voltage dictates whether
a process is energetically possible in principle, whereas the kinetics of a
process is controlled by specific properties of the molecular junctions, such
as molecule-lead and vibronic coupling.

A schematic diagram of relevant dissociation mechanisms is shown in Fig. 3.
When vibronic coupling is weak, stepwise vibrational ladder climbing (cf. Fig.
3, process M1) is the dominant dissociation mechanism over the whole bias
range. In this regime, the probability of inelastic transport processes, where
the vibrational mode is excited by multiple vibrational quanta, is small. In
the case of stronger vibronic coupling, multi-quantum vibrational excitations
are favored. When the bias voltage exceeds the dissociation threshold,
$e\Phi>2E_{\rm D}$, the energy of an incoming electron is sufficient to
directly excite the molecule into unbound electronic states, leading to
ultrafast direct dissociation (M2). If the bias voltage is lower, dissociation
can be induced by multiple electronic transitions via multi-quantum vibration
excitations (M3). Similar mechanisms have been invoked to explain molecule
desorption from metal surfaces using scanning tunneling microscopy.Stipe  _et
al._ (1997); Lee, Sorescu, and Deng (2011); Tan  _et al._ (2011); Zhao  _et
al._ (2013); Chen  _et al._ (2019) In this context, it was found that the
dissociation rate shows a linear (power-law) dependence on the tunneling
current when the desorption is induced by a single (multiple) electronic
transition(s).Salam, Persson, and Palmer (1994); Lee, Sorescu, and Deng
(2011); Ueba (2003); Tikhodeev and Ueba (2004); Persson and Avouris (1997)

Figure 2: Vibrational heating and cooling processes in molecular junctions.
(a), (b): Transport-related heating (a) and cooling (b) processes, where the
transport of an electron from one lead to the other is accompanied by the
emission or absorption of vibrational energy. (c): Electron-hole pair creation
process, where an electron is transferred from the valence band of one lead
onto the molecule, absorbs energy from the vibrations and then returns to the
conduction band of the same lead. (d): Vibrational relaxation process due to
the coupling to a phonon bath. Figure 3:  Sketch of current-induced
dissociation mechanisms in molecular junctions under various conditions. The
horizontal green dotted line marks the boundary between off-resonant and
resonant transport regimes. Red dashed lines indicate the transition region
between dominant dissociation mechanism.

a) $\Delta Q=0.025$Å

b) $\Delta Q=0.3$Å

Figure 4:  Average vibrational excitation $\langle n_{\mathrm{vib}}\rangle$ as
well as dissociation rate $k_{\mathrm{diss}}$ as a function of bias voltage
for $\Delta Q=0.025\textrm{\AA}$ (a) and $\Delta Q=0.3\textrm{ \AA}$ (b),
respectively. The close-up of the dissociation rate in the low-bias regime is
shown in the respective inset. Solid and dotted lines represent results
without and with the coupling to a thermal phonon bath, i.e., $\lambda_{\rm
ph}=$ 0 and $\lambda_{\rm ph}=$ 0.05 eV, respectively. Other parameters are
$\Gamma_{\rm L}=\Gamma_{\rm R}=$ 0.05 eV and $\omega_{c}=3\Omega_{0}$.

To demonstrate the characteristics of the different dissociation mechanisms,
Fig. 4 displays the renormalized average vibrational excitation $\langle
n_{\mathrm{vib}}\rangle$ as well as dissociation rate $k_{\mathrm{diss}}$ as a
function of bias voltage for values of the displacement between the neutral
and charged PES of $\Delta Q=0.025\textrm{\AA}$ and $\Delta
Q=0.3\textrm{\AA}$, respectively. These two values of $\Delta Q$ represent
cases of weak and strong vibronic coupling, respectively. The molecule is
coupled symmetrically to both leads and a weak molecule-lead coupling strength
of $\Gamma_{\rm L}=\Gamma_{\rm R}=0.05$ eV is chosen to minimize the influence
of electronic level broadening.

We first discuss the case of weak vibronic coupling ($\Delta Q=0.025\textrm{
\AA}$), in which stepwise vibrational ladder climbing (process M1 in Fig. 3)
is the dominant dissociation mechanism. The results in the left panel of Fig.
4 (a) show that the vibrational excitation $\langle n_{\rm vib}\rangle$
increases from negligible values in the off-resonant transport regime to
moderate values in the resonant regime upon increase of the bias voltage. A
similar behavior is observed for the dissociation rate $k_{\rm diss}$ depicted
in the right panel of Fig. 4 (a). It is noted, though, that the dissociation
rate only increases at voltages of about $\Phi=2.2$ V, which is already above
the onset of resonant transport at 2 V. Additional vibrational relaxation
induced by coupling to a phonon bath ($\lambda_{\rm ph}=0.05$ eV) results in a
substantial reduction of both the vibrational excitation and the dissociation
rate.

In the off-resonant regime, electron transport is dominated by cotunneling and
thus current-induced heating is slow and ineffective. In principle, once the
bias voltage exceeds the threshold $e\Phi>E_{10}$, where $E_{10}=0.258$ eV is
the energy gap between the vibrational ground and first excited state, the
molecule can be excited by a succession of inelastic electrons tunneling
processes to high-lying vibrational bound states in a stepwise manner. In
practice, however, efficient dissociation also requires that vibrational
heating [cf. Fig. 2 (a)] must be faster than vibrational cooling [cf. Fig. 2
(b), (c) and (d)]. As a result, in the off-resonant transpot regime, the
molecule remains preferentially in the vibrational ground state and the
junction is stable.

In the resonant transport regime, the current is significantly larger and
heating becomes significant. For bias voltages in the vicinity of the onset of
resonant transport (in the present model at 2 V), however, cooling effects due
to electron-hole pair creation [cf. Fig. 2 (c)] counteract current-induced
heating in the course of stepwise vibrational ladder climbing, such that the
onset of dissociation appears at a somewhat higher bias voltage. In the high
bias voltage regime, where electron-hole pair creation processes are fully
blocked, harmonic models predict an extremely large vibrational excitation
(vibrational instability). Mitra, Aleiner, and Millis (2004); Härtle and Thoss
(2011); Schinabeck, Härtle, and Thoss (2018) For the more realistic
dissociative model considered here, however, $\langle n_{\mathrm{vib}}\rangle$
saturates to a moderate value. This is the combined result of stepwise heating
and the presence of a dissociation threshold. Dissociation happens before the
vibrational instability can set in. Finally, when energy dissipation to a
phonon bath is efficient and fast enough to compete with current-induced
heating, dissociation is almost completely suppressed.

Next, we consider the case of strong vibronic coupling depicted in Fig. 4 (b)
for the example of $\Delta Q=0.3\textrm{ \AA}$. The results differ from those
obtained in the weak coupling regime in three aspects. First, in contrast to
the saturation observed in the high bias regime for $\Delta Q=0.025\textrm{
\AA}$, both $\langle n_{\mathrm{vib}}\rangle$ and $k_{\mathrm{diss}}$ increase
continuously with increasing bias voltage within the sampled bias window.
Second, the dissociation rate is already significant at the onset of resonant
transport, as shown in the inset of Fig. 4 (b). For example, at a bias voltage
of 1.8 V, the dissociation rate is about 0.5 $\rm{ns}^{-1}$. Furthermore, the
coupling to a phonon bath only slightly reduces the vibrational excitation and
the dissociation rate.

These distinct features point to a different dissociation mechanism. Strong
vibronic coupling allows direct vibrational excitation from low-lying bound
states to continuum states. For higher-bias voltages in the resonant transport
regime this facilitates direct dissociation (process M2 in Fig. 3). Because
dissociation in this case is very fast, vibrational relaxation due to
electron-hole pair creation or coupling to a phonon bath is inefficient.
Moreover, dissociation in the off-resonant lower-voltage regime is possible
because energy exchange with only a few electrons is sufficient for the
molecule to overcome the dissociation threshold, as inelastic transport
processes are accompanied by multi-quantum vibrational excitations (process M3
in Fig. 3).

###  III.2 Detailed analysis of current-induced dissociation

In the following, we analyze in detail the mechanisms of current-induced
dissociation in molecular junctions. To simplify the discussion, the coupling
to a phonon bath is neglected throughout this section. We will study the
influence of the phonon bath in Sec. III.3.2.

####  III.2.1 High bias voltage regime

a)

b)

c)

Figure 5: (a) Franck-Condon factors for the system with $\Delta Q=0.01\textrm{
\AA}$. (b) Population dynamics at $\Phi=8$ V. The upper and lower panel
correspond to the charged and neutral state, respectively. In each panel, the
renormalized populations of vibrational bound states
$P^{e/g}_{v}/(1-P_{\infty})$ ($v=0-15$) are shown in color solid lines and the
summation over the populations of the continuum states is shown in black
dotted lines. The renormalization factor is the survival probability
$1-P_{\infty}(t)$, with $P_{\infty}(t)=P_{\infty}^{g}+P_{\infty}^{e}(t)$.
$P^{g}_{\infty}(t)$ and $P^{e}_{\infty}(t)$ are plotted in green dashed lines.
(c) Schematic illustration of the stepwise vibrational heating. Other
parameters are $\Gamma_{\rm L/R}=0.05eV$.

We first consider the high bias voltage regime, exemplified by $\Phi=8$ V. In
this regime, electron-hole pair creation processes are fully blocked and,
thus, only transport-related vibrational heating and cooling processes are
active. Furthermore, electron transport takes place resonantly and direct
dissociation (process M2 in Fig. 3) is energetically possible. We note that
for such high bias voltages, in realistic molecular junctions more than a
single electronic state included in the model may enter the bias window. This
can result in additional phenomena,Härtle, Benesch, and Thoss (2009) which
will be studied in future work.

We start the analysis with the weak vibronic coupling case, $\Delta
Q=0.01\textrm{ \AA}$. The Franck-Condon transition matrix depicted in Fig. 5
(a) shows that in this regime inelastic transport processes are dominated by
single-quantum vibrational excitation and deexcitation. Fig. 5 (b) depicts the
population dynamics in the vibrational state manifold of the neutral and the
charged state for a symmetric molecule-lead coupling scenario. Starting in the
initially prepared vibrational ground state, inelastic transport processes
excite the molecule to the first vibrationally excited state after a cycle of
charging and discharging, as illustrated by the black arrows in Fig. 5 (c). In
a sequence of such inelastic transport processes higher vibrationally excited
states are sequentially populated, as shown by the short-time dynamics in Fig.
5 (b). On a timescale of tens of picoseconds, a quasi-steady distribution is
established. The population of the vibrational continuum states above the
dissociation threshold (black dotted lines) rises on the same timescale as the
population of the highest bound state, which confirms the stepwise vibrational
heating mechanism.

a)

b)

c)

Figure 6: Same as Fig. 5, but for $\Delta Q=0.3\textrm{ \AA}$.  Figure 7:
Dissociation rate as a function of current $I_{c}$ for $\Delta Q=0.3\textrm{
\AA}$ and $\Phi=8$ V. The data are obtained by varying $\Gamma_{\rm L/R}$ from
0.001 to 0.01 eV. $I_{c}$ is the plateau value of $I(t)/(1-P_{\infty}(t))$, as
shown as an example for $\Gamma_{\rm L/R}=0.005$ eV in the inset.

In the case of strong vibronic coupling ($\Delta Q=0.3$Å), depicted in Fig. 6,
excitation and deexcitation processes involving multiple vibrational quanta
are favored, as confirmed by the Franck-Condon matrix (see Fig. 6 (a)). The
population dynamics of the vibrational states displayed in Fig. 6 (b) show
that a dissociation probability of 100% (sum of green dashed lines) is reached
within one picosecond and the rescaled population of continuum states (black
dotted lines) is rather high. These observations point to the existence of two
direct dissociation pathways that are induced by a single tunneling electron.
At a voltage of $\Phi=$ 8 V, the incoming electron can excite the molecule
directly from the vibrational ground state or a low-lying vibrationally
excited state into a continuum state, as schematically illustrated by path 1
in Fig. 6 (c). Alternatively, in path 2, the molecule is first charged and
excited by an incoming electron into a high-lying vibrational bound state
($6\leq n_{v}\leq 15$), and then proceeds resonantly into a continuum state
upon discharging, corresponding to dissociation mediated by a Feshbach
resonance.Brisker and Peskin (2008)

It is known from studies of molecular desorption from metal surfaces that the
reaction rate depends linearly on the tunneling current when the desorption is
induced by a single electronic transition. Fig. 7 shows the dissociation rate
$k_{\mathrm{diss}}$ as a function of the quasi steady-state current of the
undissociated molecule $I_{c}$. Here, $I_{c}$ is obtained in our model by
taking the plateau value of $I(t)/(1-P_{\infty}(t))$, as illustrated in the
inset of Fig. 7. The data are obtained by varying the molecule-lead coupling
$\Gamma_{L/R}$ from 0.001 to 0.01 eV. The least-squares fitting of
$k_{\mathrm{diss}}\propto I_{c}^{n}$ yields $n=0.98$, which corroborates our
hypothesis that in the case of strong vibronic coupling and high bias voltage
($e\Phi>2E_{\rm D}$), dissociation is induced by a single tunneling electron.
We note that deviations from the linear dependence are observed when
$\Gamma_{L/R}$ is larger than 0.02 eV. The influence of stronger molecule-lead
coupling will be discussed in Sec. III.3.1.

a)

b)

Figure 8: (a) Dissociation rate $k_{\mathrm{diss}}$ as a function of $\Delta
Q$ at $\Phi=8$ V. The molecule-lead coupling is $\Gamma_{\rm L/R}=$ 0.05 eV.
(b) Franck-Condon matrix elements
$S_{01}=\left|\langle\psi^{e}_{1}|\psi^{g}_{0}\rangle\right|^{2}$ and the
summation of the transition probabilities from the vibrational ground state to
the sixth and higher vibrationally excited states, $S_{0c}=\sum_{v\geq
6}\left|\langle\psi^{e}_{v}|\psi^{g}_{0}\rangle\right|^{2}$.

Having distinguished different dissociation mechanisms in the limiting cases
of weak and strong vibronic coupling, we next analyze the intermediate regime.
To this end, Fig. 8 (a) shows the dissociation rate as a function of the
displacement $\Delta Q$. The dissociation rate exhibits a non-monotonous
dependence on the displacement with a first increase up to about $\Delta
Q=0.16\textrm{ \AA}$, followed by a slight decrease in the range of 0.16
$\textrm{ \AA}<\Delta Q<0.2\textrm{ \AA}$ and then a further increase. The
slope of the second increase is much smaller than that of the first one.

In order to explain these findings, we examine the dependence of Franck-Condon
matrix elements
$S_{vv^{\prime}}=\left|\langle\psi^{e}_{v^{\prime}}|\psi^{g}_{v}\rangle\right|^{2}$
on $\Delta Q$ in Fig. 8 (b). The excitation rate of the $v=0\rightarrow
v^{\prime}=1$ transition is proportional to $S_{01}$. Moreover,
$S_{0c}=\sum_{v\geq
6}\left|\langle\psi^{e}_{v}|\psi^{g}_{0}\rangle\right|^{2}$ is a measure for
the transition probability from the vibrational ground state to the sixth and
higher vibrationally excited states, which is related to the direct
dissociation mechanism. $S_{01}$ exhibits a turnover at $\Delta Q=0.16\textrm{
\AA}$ and $S_{0c}$ increases monotonically and exceeds $S_{01}$ at $\Delta
Q=0.23\textrm{ \AA}$. This is in line with the transitions observed in the
dissociation rate and suggests the following mechanisms in the different
regimes: In the weak to intermediate vibronic coupling regime ($0<\Delta
Q<0.16\textrm{ \AA}$), dissociation is dominated by stepwise vibrational
ladder climbing. In this mechanism, which comprises multiple consecutive
heating steps, the dissociation rate scales non-linearly with the single-step
excitation rate, $k_{s}$,Salam, Persson, and Palmer (1994); Ueba (2003) i.e.,
$k_{\mathrm{diss}}\propto k_{s}^{n}$ ($n\gg 1$). Thus, as $k_{s}$ ($\propto
S_{01}$) increases, the dissociation rate rises also. For $\Delta Q>0.16\text{
\AA}$, $S_{01}$ starts to decrease and multi-quantum vibrational excitations
are preferred. The transition of the dominant dissociation mechanism takes
place in this region. As such, the dissociation rate drops only slightly and
increases again when the dissociation is mainly induced by a single electronic
transition, which is the case for $\Delta Q>0.22\text{ \AA}$, where $S_{0c}$
is close to or larger than $S_{01}$.

####  III.2.2 Intermediate bias voltage regime

a) SYMM ($\Gamma_{\rm R}=\Gamma_{\rm L}$)

b) SYMM ($\Gamma_{\rm R}=\Gamma_{\rm L}$)

c) ASYMM ($\Gamma_{\rm R}=0.1\Gamma_{\rm L}$)

Figure 9: Dissociation rate $k_{\mathrm{diss}}$ as a function of $\Delta Q$
for the model SYMM in (a) and (b) and for the model ASYMM in (c). Different
lines correspond to different bias voltages. The coupling of the molecule to
the left lead is fixed at $\Gamma_{\rm L}=0.05$ eV.  Figure 10: Energy-level
scheme illustrating the suppression of electron-hole pair creation processes
with increasing bias voltage for different $\Delta Q$. Shown is the Fermi
distribution of the electrons in the left lead (yellow) for lower and higher
bias voltage, the molecular energy level as well electron-hole pair creation
processes including the absorption of up to three vibrational quanta.

a)

b)

Figure 11: (a) Dissociation rate $k_{\mathrm{diss}}$ as a function of the
current $I_{c}$ for model SYMM with $\Gamma_{\rm L/R}=0.05$ eV. The colored
lines varying from dark purple to yellow correspond to $\Delta Q$ from $0.025$
to $0.3\textrm{ \AA}$ with the increment of $0.025\textrm{ \AA}$. The gray
dotted lines represent a least-squares fitting. (b) Fitting parameter $n$ in
$k_{\mathrm{diss}}\propto I_{c}^{n}$ for different $\Delta Q$.

Next, we turn to the intermediate bias voltage regime, 2 V $\lesssim\Phi<4.76$
V $=2E_{\rm D}/e$. In this regime, electron transport is still resonant,
however, direct dissociation (process M2 in Fig. 3) is no longer possible.
Furthermore, electron-hole pair creation processes are active. In order to
distinguish between transport induced vibrational mechanisms and the influence
of electron-hole pair creation processes, it is insightful to also consider
systems that display different symmetries with respect to the coupling to the
left and the right lead. Fig. 9 depicts the dissociation rates in this regime
for a symmetric (SYMM, panels (a) and (b)) and an asymmetric (ASYMM, panel
(c)) molecule-lead coupling scenario, respectively. The following analysis
focuses exclusively on positive $\Delta Q$. Results for negative $\Delta Q$
will be discussed in Sec. III.3.3.

Overall, the results reveal a more complex dependence of the dissociation rate
on the vibronic coupling $\Delta Q$, as compared to the high bias voltage
regime in Fig. 8. For moderate voltages (e.g. 4 V, depicted by the blue line
in Fig. 9 (a)), the dissociation rate rises for small $\Delta Q$, then
exhibits a turnover, followed by a local minimum and a further increase with
additional structures for larger $\Delta Q$. Lowering the voltage (e.g. to 3 V
or 2.6 V, depicted by the orange and green line in Fig. 9 (a), respectively),
the dissociation rate overall decreases and the additional structures at
larger $\Delta Q$ disappear. For even lower voltage, (e.g. 2.2 V or 1.8 V
shown in Fig. 9 (b)), the dissociation rate increases monotonically with
$\Delta Q$. A similar behavior is observed for the case of asymmetric
molecule-lead coupling depicted in Fig. 9 (c).

We first analyze these observations in the weak to intermediate vibronic
coupling regime. As mentioned before, dissociation in this regime is dominated
by the stepwise vibrational ladder climbing mechanism, which is particularly
sensitive to relaxation effects. In the model considered here, vibrational
relaxation is caused by transport related deexcitation [cf. Fig. 2 (b)] and
electron-hole pair creation [cf. Fig. 2 (c)] processes. The efficiency of
electron-hole pair creation processes depends sensitively on the applied bias
voltage.Härtle and Thoss (2011); Härtle, Peskin, and Thoss (2013); Schinabeck,
Härtle, and Thoss (2018) As shown previously,Härtle and Kulkarni (2015);
Nitzan and Galperin (2018) cooling due to electron-hole pair creation
processes is particularly effective at the onset of resonant transport, but
becomes successively blocked for larger voltages. As a consequence, the
dissociation rate is significantly reduced with decreasing bias voltage, where
pair creation processes become less suppressed. Furthermore, it should be
noted that with increasing $\Delta Q$, vibrational deexcitation processes with
more vibrational quanta are enabled and become dominant in the electron-hole
pair creation processes, as illustrated in Fig. 10. The first turnover of the
dissociation rate with increasing $\Delta Q$ indicates where the electron-hole
pair creation processes are no longer largely blocked, and it is shifted to a
larger $\Delta Q$ with increasing bias voltage.

The importance of electron-hole pair creation processes can be confirmed by
considering the case of asymmetric molecule-lead coupling (model ASYMM,
depicted in Fig. 9 (c)). As is known from the study of asymmetric models with
harmonic vibrations, the cooling efficiency of electron-hole pair creation
processes depends sensitively on the bias polarity.Härtle and Thoss (2011);
Härtle, Peskin, and Thoss (2013); Schinabeck, Härtle, and Thoss (2018) Fig. 9
(c) shows that, while the dissociation rate is independent on bias polarity
for a large bias voltage of 8 V when electron-hole pair creation processes are
blocked, it does depend on the bias polarity for moderate and small voltages.
For example, a bias voltage of +4 V results in a significant smaller
dissociation rate than a voltage of -4 V. This is because for positive bias
voltage, electron-hole pairs are generated predominantly with respect to the
left lead and for negative bias voltage mostly with respect to the right lead.
As in model ASYMM the molecule is stronger coupled to the left lead, the
cooling due to electron-hole pair creation processes is more effective for
positive bias voltage and thus the dissociation rate is smaller. This cooling
effect becomes stronger with increasing $\Delta Q$, because more electron-hole
pair creation processes accompanied by multi-quantum vibrational deexcitations
are kinetically allowed. As a consequence, a pronounced negative slope is
observed in the intermediate vibronic coupling regime, which was also found in
Ref. Koch  _et al._ , 2006.

Next, we consider the intermediate to strong vibronic coupling regime. In this
regime, the dominant dissociation mechanism undergoes a transition from
vibrational ladder climbing to multi-quantum vibrational excitations to the
continuum states. The crossover also depends on the current.Salam, Persson,
and Palmer (1994) As mentioned before, at a high bias voltage of 8 V, the
dominance of stepwise vibrational ladder climbing extends to $\Delta
Q=0.16\textrm{ \AA}$. At lower bias voltages, the current is decreased and,
simultaneously, vibrational relaxation is more effective. As a result, the
average time between electron tunneling events becomes comparable or longer
than the vibrational relaxation time. In this situation, as shown by Salam and
coworkers,Salam, Persson, and Palmer (1994) vibrational ladder climbing is
less effective and multi-quantum vibrational excitation to the continuum
states becomes more important. As a consequence, the transition of the
dominant dissociation mechanism can take place at a smaller $\Delta Q$ for a
smaller bias voltage. The dissociation rate increases with increasing vibronic
coupling because of the increased transition probability to continuum states
and the decreased sensitivity to vibrational relaxation effects. The latter is
due to the fact that the dissociation after excitation above the dissociation
threshold is faster than the cooling rates.

In order to verify the above assertion, Fig. 11 (a) shows the dissociation
rate $k_{\mathrm{diss}}$ as a function of the quasi steady-state current of
the undissociated molecule, $I_{c}$, for various $\Delta Q$. The data were
obtained by varying the bias voltage in the range of 2.6 - 3.2 V. It is known
from studies of molecular desorption from metal surfaces that when the
reaction is induced by multiple electronic transitions, the reaction rates
exhibit a power-law dependence on the tunneling current,
$k_{\mathrm{diss}}\propto I_{c}^{n}$, and the exponent $n$ was interpreted as
the number of electrons involved to induce the reaction.Salam, Persson, and
Palmer (1994); Ueba (2003); Tikhodeev and Ueba (2004); Persson and Avouris
(1997) The results in Fig. 11 indicate indeed a power-law relation with values
of $n$ varying from $n=60$ at $\Delta Q=0.025\text{ \AA}$ to $n=7$ at $\Delta
Q=0.3\text{ \AA}$, thus confirming the multiple-electron nature of the
mechanism. The dependence of the power $n$ on $\Delta Q$ depicted in Fig. 11
(b) reveals a pronounced decrease at about $\Delta Q=0.1\text{ \AA}$, which is
in line with the onset of multi-phonon inelastic processes, where fewer
electronic transitions are needed to reach the dissociation threshold. This
confirms a transition of the dominant dissociation mechanism from stepwise
vibrational ladder climbing to multi-quantum vibrational excitations. When
multi-quantum vibrational excitations become dominant, increasing the vibronic
coupling leads to the increased probability of the direct transition from the
lower-lying vibrational states to the continuum states. Therefore, fewer
electrons are involved to induce the dissociation and the dissociation rate is
increased.

####  III.2.3 Lower bias voltage regime

Figure 12: Population dynamics for $\Delta Q=0.3\text{\AA}$ at 1.8 V. The
molecule-lead coupling is $\Gamma_{\rm L/R}=$ 0.05 eV. The upper and lower
panel corresponds to the charged and neutral state, respectively. For details,
see Fig. 5.

Finally, we consider the lower bias voltage regime before the onset of
resonant transport. In many cases studied experimentally, single-molecule
junctions are found to be rather stable at low bias voltage.Su  _et al._
(2016) But a recent work has reported that bond rupture can also take place in
the off-resonant regime.Li  _et al._ (2016)

We have performed a series of calculations for bias voltages $\Phi<2$ V and a
range of $\Delta Q$. At bias voltages $\Phi\lesssim 1.6$ V, the dissociation
rate is negligible over the time scale accessible with our simulations, which
employ direct time propagation. In this regime, the use of different reaction
rate approaches, based, e.g., on the flux correlation function formalism is
necessary, which will be the subject of future work.

Fig. 9 (b) shows the dissociation rate at a bias voltage of $\Phi=1.8$ V,
which is just below the threshold for resonant transport. In this regime the
dissociation rate is overall very small, but shows significant values for
$\Delta Q>0.25\text{ \AA}$ In order to analyze the dissociation mechanism in
this case, the corresponding population dynamics are shown in Fig. 12 for
$\Delta Q=0.3\text{ \AA}$. It can be seen that the molecule remains
predominantly in the neutral state in this regime. At short times (before 1
ps), the dissociation probability $P^{g}_{\infty}$ (green dashed line in the
lower panel) exhibits a step-like increase, which reflects the sudden switch-
on of bias voltage and molecule-lead coupling at the initial moment. After
this transient dynamics, through co-tunneling transport-induced vibrational
heating, a quasi-steady distribution of vibrationally excited states is formed
at about one picosecond. This steady distribution is independent on the
initial preparations. After the establishment of the quasi-steady
distribution, the dissociation probability $P^{g}_{\infty}$ increases steadily
over time. In this case of low current and strong vibrational relaxation due
to the electron-hole pair creation processes, dissociation is dominated by the
mechanism induced by multi-quantum vibrational excitations to the continuum
states. The energetic analysis suggests that at a bias voltage of $\Phi=1.8$
V, the lowest vibrationally excited state allowed to be excited into the
continuum states (via Feshbach resonance) is $\nu=7$. The required energy is
$\Delta E=E_{\rm D}-E^{g}_{7}=0.77$ eV, which is smaller than the chemical
potential of the left lead, $\mu_{\rm L}=0.9$ eV.

###  III.3 Further aspects

In this section, we study further aspects that are important to unravel the
mechanisms of current-induced dissociation in molecular junctions. This
includes the effects of the molecule-lead coupling strength, the influence of
vibrational relaxation induced by coupling to a phonon bath, as well as
situations, where the considered bond shortens upon charging of the molecule,
i.e. $\Delta Q<0$.

####  III.3.1 Influence of molecule-lead coupling strength

a) $\Phi=$ 4 V

b) $\Delta Q=0.2$Å

c) $\Phi=$ 4 V, $\Delta Q=0.2$Å

Figure 13: Dissociation rate $k_{\mathrm{diss}}$ as a function of the
molecule-lead coupling strength $\Gamma$ for different displacement $\Delta Q$
(a), bias voltages (b), and molecular vibrational frequencies (c). We assume a
symmetric molecule-lead coupling scenario, $\Gamma_{\rm L}=\Gamma_{\rm
R}=\Gamma$. The molecular vibrational frequency $\hbar\Omega_{0}=137$ meV in
panel (c) is obtained by setting the nuclear mass to $M=4$ amu but keeping the
Morse potential parameters $E_{\rm D}=2.38$ eV and $a=$1.028 $a_{0}^{-1}$
unchanged. For all other cases shown, the molecular frequency is
$\hbar\Omega_{0}=$ 274 meV.

So far, we have considered molecular junctions with weak coupling to the
electrodes. The molecule-lead coupling strength depends on the specific
molecule, the anchoring group and geometry as well as the electrode
material.Xin  _et al._ (2019) For example, it has been reported that graphene
electrodes can provide strong covalent bonding to molecules.Sun  _et al._
(2018); Leitherer, Papior, and Brandbyge (2019) In the following, we analyze
the influence of the molecule-lead coupling strength on dissociation dynamics
in the symmetric coupling scenario, $\Gamma_{\rm L}=\Gamma_{\rm R}$.

Fig. 13 shows the dissociation rate $k_{\mathrm{diss}}$ as a function of the
molecule-lead coupling strength $\Gamma_{\rm L}$ for different potential
surface displacements $\Delta Q$ and various bias voltages. In addition,
results for different fundamental vibrational frequencies $\hbar\Omega_{0}$,
obtained by varying the reduced mass $M$, are shown. The range of the
molecule-lead coupling $\Gamma_{\rm L/R}$ covers the transition from non-
adiabatic ($\Gamma_{\rm L/R}\ll\hbar\Omega_{0}$) to adiabatic ($\Gamma_{\rm
L/R}\gg\hbar\Omega_{0}$) transport. A turnover of the dissociation rate is
observed for all parameter sets.

Increasing the molecule-lead coupling can affect the dissociation dynamics in
different ways. On the one hand, increasing $\Gamma_{\rm L/R}$ leads to a
shorter resonance state lifetime and a larger current. On the other hand, it
also increases the adiabaticity of transport dynamics and facilitates
electron-hole pair creation processes.

The dissociation rate is determined by the number of electrons passing through
the molecule per unit time (i.e. the current) and the amount of energy
transferred per electron. When $\Gamma_{\rm L/R}\ll\hbar\Omega_{0}$, i.e. in
the non-adiabatic regime, the rate-determining factor is the current, which
increases with the molecule-lead coupling. In the adiabatic regime,
$\Gamma_{\rm L/R}\gg\hbar\Omega_{0}$, however, electrons tunnel through the
molecule much faster than the vibrational period such that the energy exchange
efficiency is very low. In this case, notwithstanding the increased current,
the transmitted energy per tunneling electron to the molecule reduces with
increasing $\Gamma_{\rm L/R}$. The trade-off between these two counteracting
effects leads to the turnover of the dissociation rate. The slight shift of
the turnover to a smaller $\Gamma$ with increasing displacement and decreasing
bias voltage is due to the electron-hole pair creation processes, which are
facilitated by stronger molecule-lead coupling.

As a side note, the results obtained for larger molecule-lead coupling
$\Gamma$ can be used to evaluate the range of validity of low-order kinetic
schemes.Härtle and Kulkarni (2015); Foti and Vázquez (2018) We compared the
results obtained using different hierarchical truncation tiers for the
parameters shown in Fig. 13 and found that when $\Gamma>\hbar\Omega_{0}/2$,
the first-tier results (equivalent to a second order perturbative treatment of
molecule-lead coupling) always overestimate the dissociation rate.

####  III.3.2 Vibrational relaxation due to coupling to a phonon bath

a) $\omega_{c}=3\Omega_{0}$

b) $\Delta Q=0.025$Å

c) $\Delta Q=0.3$Å

Figure 14:  Dissociation rate $k_{\mathrm{diss}}$ for a scenario which
includes vibrational relaxation due to the coupling to a phonon bath. (a)
Dissociation rate $k_{\mathrm{diss}}$ as a function of $\Delta Q$ for
$\lambda_{\rm ph}=0$ and 0.1 eV. For comparison, results without coupling to a
phonon bath are shown by the dotted line. (b,c) Dissociation rate
$k_{\mathrm{diss}}$ as a function of $\lambda_{\rm ph}$ for $\Delta
Q=0.025\textrm{ \AA}$ (b) and $\Delta Q=0.3\textrm{ \AA}$ (c), respectively.
The molecule-lead coupling is fixed at $\Gamma_{\rm L}=\Gamma_{\rm R}=$ 0.05
eV and the bias voltage is $\Phi=4$ V.

Vibrational relaxation in molecular junctions can also be induced by coupling
of the considered reaction mode to other inactive modes (intramolecular
vibrational relaxation), the phonons of the leads or a possible solution
environment. Here, we consider the effect of such additional relaxation
processes on current-induced dissociation by coupling the reaction mode to a
bath of phonons characterized by Lorentz-Drude spectral density function as
described in Sec. II.2. The influence of this coupling is determined by two
parameters, the coupling strength $\lambda_{\rm ph}$ and the characteristic
frequency of the phonon bath $\omega_{c}$. In the overview, Fig. 4 (a), it was
already shown that the coupling to a phonon bath strongly quenches heating and
dissociation. Here, we analyze this effect in more detail based on the data
depicted in Fig. 14.

Fig. 14 (a) shows that the vibrational relaxation process induced by coupling
to the phonon bath causes a pronounced reduction of the dissociation rate, in
particular for small vibronic coupling $|\Delta Q|<0.05\textrm{ \AA}$. To
analyze this effect in more detail, Fig. 14 (b) and (c) depict the
dissociation rate as a function of the bath coupling parameter $\lambda_{\rm
ph}$ for $\Delta Q=0.025\textrm{ \AA}$ (b) and $\Delta Q=0.3\textrm{ \AA}$
(c). Two cut-off frequencies of the phonon bath are chosen, which are
smaller/larger than the vibrational frequency $\hbar\Omega_{0}$.

In the case of weak vibronic coupling, as shown in Fig. 14 (b), the
dissociation rate drops quickly to very small values upon increasing
$\lambda_{\rm ph}$, especially for a fast phonon bath,
$\omega_{c}=3\Omega_{0}$. We emphasize that, different from the relaxation
effect due to electron-hole pair creation processes, this behavior is found to
be independent on the bias voltage. The remarkably effective suppression of
the dissociation is mainly due to the fact that in this regime, dissociation
is dominated by stepwise vibrational ladder climbing, which involves many
consecutive steps. Even if the rate of every step $k_{s}$ is only slightly
reduced by vibrational relaxation, the dissociation rate
$k_{\mathrm{diss}}\propto k_{s}^{n}$ can be many orders of magnitude smaller.

For strong vibronic coupling, as shown in Fig. 14 (c) for $\Delta
Q=0.3\textrm{ \AA}$, the dissociation rate is reduced with increasing
$\lambda_{\rm ph}$ and $\omega_{c}$, but not as pronounced as in the weak
vibronic coupling case. We recall that dissociation in this regime is induced
by multi-quantum vibrational excitations. At the intermediate bias voltage of
4 V considered in Fig. 14 (c), more than a single electronic transition is
required for the molecule to finally reach the dissociation energy. As long as
the vibrational state of the molecule has not yet reached the dissociative
continuum, it is sensitive to vibrational relaxation, which reduces the
dissociation rate. For high bias voltages (e.g. $\Phi=8$ V, data not shown),
direct dissociation becomes the dominating process and the influence of
vibrational relaxation is even smaller.

####  III.3.3 Molecular junctions with negative displacement $\Delta Q$

a)

b)

c) $\Delta Q=-0.1\text{\AA}$

d) $\Delta Q=0.1\text{\AA}$

Figure 15: Franck-Condom factors for the system with $\Delta Q=-0.1\textrm{
\AA}$ (a) and $\Delta Q=0.1\textrm{ \AA}$ (b), respectively. Population
dynamics for $\Delta Q=-0.1\textrm{ \AA}$ (c) and $\Delta Q=0.1\textrm{ \AA}$
(d), respectively. The bias voltage is $\Phi=2.6$ V and $\Gamma_{\rm L/R}=$
0.05 eV.

The charging of the molecule leads a reorganization of the nuclear geometry.
Depending on the specific situation, this may result in bond stretching
($\Delta Q>0$) or compression ($\Delta Q<0$). In the above sections, we
focused on the case with positive $\Delta Q$. Here, we provide a brief
comparative discussion of the case of negative $\Delta Q$. To facilitate the
comparison, the dissociation rates for negative $\Delta Q$ are depicted in the
same plot as for positive $\Delta Q$ in Fig. 9.

It should be noted that within a harmonic model for the vibrational degrees of
freedom, the dynamics is independent on the sign of $\Delta Q$. Thus, the
dependence of the dynamics on the sign of $\Delta Q$ discussed below is also a
manifestation of the anharmonicity of the model. This was also found in a
previous study.Brisker and Peskin (2006)

The results in Fig. 9 show that the dissociation rates are rather insensitive
to the sign of $\Delta Q$ for larger bias voltages. For small to moderated
bias voltages, however, the dissociation rates for negative $\Delta Q$ are
found to be in general significantly larger than for positive $\Delta Q$. For
example, at bias voltage $\Phi=$2.6 V, the dissociation rate for $\Delta
Q=-0.1\textrm{ \AA}$ is more than three orders of magnitude larger than for
$\Delta Q=0.1\textrm{ \AA}$.

To explain the underlying mechanism in more detail, Fig. 15 (a) and (b) depict
the Franck-Condon matrices for $\Delta Q=-0.1\textrm{ \AA}$ and $0.1\textrm{
\AA}$, respectively. Both exhibit the characteristic pattern of anharmonic
models, i.e., an asymmetry with respect to the diagonal line. As an example
for a specific elementary process consider the charging of the molecule
accompanied by excitation to the sixth vibrationally excited state and a bias
voltage of $\Phi=$ 2.6 V. At this voltage, electron-hole pair creation
processes with respect to the left lead can be accompanied by the vibrational
deexcitation processes $n^{e}_{6}\rightarrow n^{g}_{v=0-4}$
($E^{e}_{6}-E^{g}_{v=0-4}>1.3$ eV). The transition probabilities of these
vibrational deexcitation processes for $\Delta Q=-0.1\textrm{ \AA}$ are
smaller than their counterparts for $\Delta Q=0.1\textrm{ \AA}$, as shown in
Fig. 15 (a) and (b). Moreover, the dissociation pathways mediated by Feshbach
resonances from the low-lying vibrationally excited states in the potential
surface of the charged state to the continuum states of the neutral molecule
are available for $\Delta Q=-0.1\textrm{ \AA}$, but are blocked for $\Delta
Q=0.1\textrm{ \AA}$. Thus, the significantly larger increase of the
dissociation rate for negative $\Delta Q$, compared to positive $\Delta Q$, is
caused by less efficient electron-hole pair creation processes and an
increased transition probability to the continuum states.

The above reasoning is confirmed by the population dynamics depicted in Fig.
15 (c) and (d). For $\Delta Q=-0.1\textrm{ \AA}$, due to the weaker electron-
hole pair creation cooling effect, a broad distribution among the vibrational
states is quickly formed and the dissociation is almost completed at 1 ps. In
contrast, for $\Delta Q=0.1\textrm{ \AA}$, the dissociation rate at ten
picoseconds is less than 0.5%. In addition, we observe that for $\Delta
Q=-0.1\textrm{ \AA}$, the population of states above the dissociation barrier
in the neutral state, $P^{g}_{\rm continuum}$, shown as the black dotted line
in the lower panel of Fig. 15 (c), is higher than the population of high-lying
vibrational bound states $(n^{g}_{v}>8)$. This is a result of the population
transfer from the low-lying vibrationally excited states in the potential
surface of the charged molecule to continuum states of the neutral molecule,
which then leads to ultrafast dissociation. But for $\Delta Q=0.1\textrm{
\AA}$, this shortcut is blocked.

The above analysis clearly shows that in more realistic anharmonic models of
molecular junctions, in addition to the strength of the vibronic coupling, its
sign also plays an important role in the current-induced dissociation
dynamics.

###  III.4 Time-dependent current-voltage characteristics and implications for
experiments

a) $\lambda_{\rm ph}=0$ eV

b) $\lambda_{\rm ph}=0.05$ eV

c) $\Delta Q=0.025$Å

d) $\Delta Q=0.1$Å

e) $\Delta Q=0.3$Å

Figure 16: Time-dependent current-voltage characteristics. The parameters
$\Delta Q$ and $\lambda_{\rm ph}$ are given in the plot, other parameters of
the calculations are $\Gamma_{\rm L}=\Gamma_{\rm R}=$ 0.05 eV and
$\omega_{c}=3\Omega_{0}$. The vertical dashed lines in (c) - (e) indicate the
onset of resonant transport.

Current-voltage characteristics serve as an important tool for acquiring
information about charge transport in molecular junctions. For a molecular
junction which undergoes bond rupture, or a different structural change, the
conductance may vary over the time of the reactive process, as was observed in
a recent experiment on the picosecond timescale.Arielly  _et al._ (2017) The
details depend on the specific process. Here, we analyze the change of the
current-voltage characteristic for a molecular junction, where the bond to a
side group ruptures and the conductance changes upon dissociation of the side
group, as described by the model introduced in Sec. II.1.

Fig. 16 shows time-dependent current-voltage characteristics for the cases of
weak and strong vibronic coupling, without (a) and with (b) coupling to a
phonon bath.

First, we consider cases without coupling to a phonon bath, depicted in Fig.
16 (a). For weak vibronic coupling, $\Delta Q=0.025\textrm{ \AA}$, transport
is dominated by elastic processes and the current-voltage characteristics for
short times ($t=20$ fs) resembles a typical IV curve of a resonant level
model. For longer times, the current at higher bias voltages decreases,
resulting in a negative differential resistance feature in the current-voltage
characteristic. This decrease of the current is a result of the dissociation
process.

For strong vibronic coupling, $\Delta Q=0.3\textrm{ \AA}$, Franck-Condon
blockade reduces the current at the onset of the resonant transport
regime.Koch, Von Oppen, and Andreev (2006); Schinabeck  _et al._ (2014) Due to
broadening and the nonequidistant level structure of the vibrational energies,
individual vibronic steps are not seen in the current-voltage characteristic.
Again, the current at higher bias voltages decreases for longer times. The
resulting negative differential resistance is much more pronounced than for
weak vibronic coupling. This is a result of the fast dissociation process in
this regime, which is almost complete within a few picoseconds.

Taking into account additional vibrational relaxation due to the coupling to a
phonon bath (Fig. 16 (b)), dissociation is suppressed for weak vibronic
coupling ($\Delta Q=0.025\textrm{ \AA}$) and thus the negative differential
resistance feature disappears. For strong vibronic coupling ($\Delta
Q=0.3\textrm{ \AA}$), on the other hand, the time-dependent current-voltage
characteristics are very similar to the case without coupling to a phonon
bath.

The behavior of the time-dependent current voltage characteristics in the
voltage region at the onset of resonant transport is depicted in more detail
in Fig. 16 (c)-(e). For weak and strong vibronic coupling, the results show
again the decrease of the current at longer times for higher bias voltages due
to the dissociation process. Interestingly, this is not observed for
intermediate vibronic coupling. The reason for this different behavior is the
effective suppression of the dissociation due to electron-hole pair creation
processes in this parameter regime, as explained in Sec. III.2.2.

These results indicate that strategies, which facilitate electron-hole pair
creation processes will be helpful to further increase the stability of
molecular junctions at moderate bias voltages.Gelbwaser-Klimovsky  _et al._
(2018); Härtle  _et al._ (2018) Other strategies to increase the stability
include the devise of junctions with efficient coupling of the molecular
vibrations to electrode phonons or a solution environment and the use of
anchoring groups that provide strong molecule-lead coupling such that the
junction operates in the adiabatic transport regime.

These findings may also be interesting in the context of recent experimental
studies which showed the change of the transport characteristics related to
bond rupture and structural changes.Li  _et al._ (2015); Capozzi  _et al._
(2016); Fung  _et al._ (2019); Zang  _et al._ (2020) Although the time-scales
observed in the experiments are significantly longer than in our study,Fung
_et al._ (2019) the basic relation between a structural change of the molecule
and the time-dependent change of the current-voltage characteristic should
appear similar.

##  IV Conclusion

We have investigated current-induced bond rupture in single-molecule junctions
employing a fully quantum mechanical method based on the HQME approach.
Extending our previous work,Erpenbeck  _et al._ (2018a, 2020) we have
considered a model, which includes more general potential energy surfaces,
accounting for both bound and continuum states of the charged molecule, as
well as vibrational relaxation processes induced by coupling of the
dissociative reaction mode to other inactive modes, the phonons of the leads
or a possible solution environment. The model also accounts for additional
dissociation channels via Feshbach resonances. Based on this model, we have
analyzed current-induced dissociation dynamics in a broad range of different
regimes, comprising off-resonant to resonant transport, weak to strong
vibronic coupling as well as non-adiabatic to adiabatic transport.

The study provides a comprehensive analysis of the reaction mechanisms
prevailing in the different regimes. Specifically, we found that for weak to
intermediate vibronic coupling, dissociation is induced by current-induced
stepwise vibrational ladder climbing. In this case, dissociation is sensitive
to vibrational relaxation. For strong vibronic coupling, multi-quantum
vibrational excitations are favored. When the applied bias voltage is high
enough, the molecule can be directly excited into a continuum state and
dissociates. Otherwise, dissociation is induced by a few electronic
transitions. Because of fast dissociation in the continuum states,
dissociation is less sensitive to vibrational relaxation in this regime.

The analysis also revealed a turnover of the dissociation rate upon increase
of molecule-lead coupling, which arises mainly from the transition from non-
adiabatic to adiabatic transport. This shows that strong molecule-lead
coupling can stabilize a molecular junction. Moreover, the results showed that
the dissociation dynamics is affected by the sign of vibronic coupling, i.e.
it exhibits different characteristics depending on whether the charging of the
molecule leads to bond stretching or compression.

Finally, it is noted that the presented method can also be used to study other
processes of current-induced reaction dynamics, such as proton transfer or
isomerization,Hofmeister, Coto, and Thoss (2017); Weckbecker, Coto, and Thoss
(2017) which are important for the realization of molecular switches, diodes
or transistors. With further extensions, it may also be useful to investigate
more complex processes of current-induced chemistry in realistic systems. For
instance, the extension of the current model to higher dimensional systems is
possible by utilizing the low-storage matrix product state representation of
the hierarchical approach.Shi  _et al._ (2018); Borrelli (2019); Yan  _et al._
(2021)

## Data Availability Statement

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

## Acknowledgements

We thank C. Kaspar and J. Bätge for helpful discussions. This work was
supported by the German Research Foundation (DFG). Y.K. gratefully
acknowledges a Research Fellowship of the Alexander von Humboldt Foundation.
A.E. was supported by the Raymond and Beverly Sackler Center for Computational
Molecular and Materials Science, Tel Aviv University. U. P. wishes to
acknowledge the Freiburg Institute for Advanced Studies for support and
stimulating atmosphere, Prof. Thoss and his group for the warm and lovely
hospitality and collaboration during his sabbatical stay at Freiburg, and the
Israel Science Foundation and the Israeli ministry of science and education
for supporting this research. Furthermore, the authors acknowledge support by
the High Performance and Cloud Computing Group at the Zentrum für
Datenverarbeitung of the University of Tübingen, the state of Baden-
Württemberg through bwHPC and the German Research Foundation (DFG) through
grant no INST 37/935-1 FUGG.

## References

  * Cuevas and Scheer (2010) J. C. Cuevas and E. Scheer,  _Molecular electronics: an introduction to theory and experiment_ (World Scientific, Singapore, 2010).
  * Galperin, Ratner, and Nitzan (2007) M. Galperin, M. A. Ratner, and A. Nitzan, “Molecular transport junctions: vibrational effects,” J. Phys.: Condens. Matter 19, 103201 (2007).
  * Bergfield and Ratner (2013) J. P. Bergfield and M. A. Ratner, “Forty years of molecular electronics: Non-equilibrium heat and charge transport at the nanoscale,” physica status solidi (b) 250, 2249–2266 (2013).
  * Aradhya and Venkataraman (2013) S. V. Aradhya and L. Venkataraman, “Single-molecule junctions beyond electronic transport,” Nat. Nanotechnol. 8, 399 (2013).
  * Bâldea (2016) I. Bâldea,  _Molecular Electronics: An Experimental and Theoretical Approach_ (CRC Press, 2016).
  * Su  _et al._ (2016) T. A. Su, M. Neupane, M. L. Steigerwald, L. Venkataraman, and C. Nuckolls, “Chemical principles of single-molecule electronics,” Nat. Rev. Mater. 1, 16002 (2016).
  * Thoss and Evers (2018) M. Thoss and F. Evers, “Perspective: Theory of quantum transport in molecular junctions,” J. Chem. Phys. 148, 030901 (2018).
  * Evers  _et al._ (2020) F. Evers, R. Korytár, S. Tewari, and J. M. van Ruitenbeek, “Advances and challenges in single-molecule electron transport,” Rev. Mod. Phys. 92, 035001 (2020).
  * Persson and Avouris (1997) B. Persson and P. Avouris, “Local bond breaking via stm-induced excitations: the role of temperature,” Surf. Sci. 390, 45–54 (1997).
  * Kim, Komeda, and Kawai (2002) Y. Kim, T. Komeda, and M. Kawai, “Single-molecule reaction and characterization by vibrational excitation,” Phys. Rev. Lett. 89, 126104 (2002).
  * Koch  _et al._ (2006) J. Koch, M. Semmelhack, F. Von Oppen, and A. Nitzan, “Current-induced nonequilibrium vibrations in single-molecule devices,” Phys. Rev. B 73, 155306 (2006).
  * Huang  _et al._ (2006) Z. Huang, B. Xu, Y. Chen, M. D. Ventra, and N. Tao, “Measurement of current-induced local heating in a single molecule junction,” Nano Lett. 6, 1240–1244 (2006).
  * Huang  _et al._ (2007) Z. Huang, F. Chen, R. D’agosta, P. A. Bennett, M. Di Ventra, and N. Tao, “Local ionic and electron heating in single-molecule junctions,” Nat. Nanotechnol. 2, 698 (2007).
  * Schulze  _et al._ (2008) G. Schulze, K. J. Franke, A. Gagliardi, G. Romano, C. Lin, A. Rosa, T. A. Niehaus, T. Frauenheim, A. Di Carlo, A. Pecchia, _et al._ , “Resonant electron heating and molecular phonon cooling in single c 60 junctions,” Phys. Rev. Lett. 100, 136801 (2008).
  * Ioffe  _et al._ (2008) Z. Ioffe, T. Shamai, A. Ophir, G. Noy, I. Yutsis, K. Kfir, O. Cheshnovsky, and Y. Selzer, “Detection of heating in current-carrying molecular junctions by raman scattering,” Nat. Nanotechnol. 3, 727–732 (2008).
  * Sabater, Untiedt, and van Ruitenbeek (2015) C. Sabater, C. Untiedt, and J. M. van Ruitenbeek, “Evidence for non-conservative current-induced forces in the breaking of au and pt atomic chains,” Beilstein J. Nanotechnol. 6, 2338–2344 (2015).
  * Li  _et al._ (2015) H. Li, T. A. Su, V. Zhang, M. L. Steigerwald, C. Nuckolls, and L. Venkataraman, “Electric field breakdown in single molecule junctions,” J. Am. Chem. Soc. 137, 5028–5033 (2015).
  * Li  _et al._ (2016) H. Li, N. T. Kim, T. A. Su, M. L. Steigerwald, C. Nuckolls, P. Darancet, J. L. Leighton, and L. Venkataraman, “Mechanism for si–si bond rupture in single molecule junctions,” J. Am. Chem. Soc. 138, 16159–16164 (2016).
  * Capozzi  _et al._ (2016) B. Capozzi, J. Z. Low, J. Xia, Z.-F. Liu, J. B. Neaton, L. M. Campos, and L. Venkataraman, “Mapping the transmission functions of single-molecule junctions,” Nano Lett. 16, 3949–3954 (2016).
  * Schinabeck (2018) C. Schinabeck,  _Hierarchical quantum master equation approaches to nonequilibrium charge transport through single-molecule junctions_ , Ph.D. thesis, Universität Erlangen-Nürnberg (2018).
  * Gelbwaser-Klimovsky  _et al._ (2018) D. Gelbwaser-Klimovsky, A. Aspuru-Guzik, M. Thoss, and U. Peskin, “High voltage assisted mechanical stabilization of single-molecule junctions,” Nano Lett. (2018).
  * Bi  _et al._ (2020) H. Bi, C.-A. Palma, Y. Gong, K. Stallhofer, M. Nuber, C. Jing, F. Meggendorfer, S. Wen, C. Yam, R. Kienberger, _et al._ , “Electron–phonon coupling in current-driven single-molecule junctions,” J. Am. Chem. Soc. 142, 3384–3391 (2020).
  * Peiris  _et al._ (2020) C. R. Peiris, S. Ciampi, E. M. Dief, J. Zhang, P. J. Canfield, A. P. Le Brun, D. S. Kosov, J. R. Reimers, and N. Darwish, “Spontaneous s–si bonding of alkanethiols to si (111)–h: towards si–molecule–si circuits,” Chem. Sci. 20, 5246–5256 (2020).
  * Ho (2002) W. Ho, “Single-molecule chemistry,” J. Chem. Phys. 117, 11033–11061 (2002).
  * Stipe  _et al._ (1997) B. Stipe, M. Rezaei, W. Ho, S. Gao, M. Persson, and B. Lundqvist, “Single-molecule dissociation by tunneling electrons,” Phys. Rev. Lett. 78, 4410 (1997).
  * Huang  _et al._ (2013) K. Huang, L. Leung, T. Lim, Z. Ning, and J. C. Polanyi, “Single-electron induces double-reaction by charge delocalization,” J. Am. Chem. Soc. 135, 6220–6225 (2013).
  * Härtle  _et al._ (2018) R. Härtle, C. Schinabeck, M. Kulkarni, D. Gelbwaser-Klimovsky, M. Thoss, and U. Peskin, “Cooling by heating in nonequilibrium nanosystems,” Phys. Rev. B 98, 081404 (2018).
  * Kuperman, Nagar, and Peskin (2020) M. Kuperman, L. Nagar, and U. Peskin, “Mechanical stabilization of nanoscale conductors by plasmon oscillations,” Nano Lett. 20, 5531–5537 (2020).
  * Li and Somorjai (2010) Y. Li and G. A. Somorjai, “Nanoscale advances in catalysis and energy applications,” Nano Lett. 10, 2289–2295 (2010).
  * Kolasinski (2012) K. W. Kolasinski,  _Surface science: foundations of catalysis and nanoscience_ (John Wiley & Sons, 2012).
  * Seideman (2016) T. Seideman,  _Current-driven phenomena in nanoelectronics_ (CRC Press, 2016).
  * Galperin, Nitzan, and Ratner (2006) M. Galperin, A. Nitzan, and M. A. Ratner, “Resonant inelastic tunneling in molecular junctions,” Phys. Rev. B 73, 045314 (2006).
  * Ryndyk, Hartung, and Cuniberti (2006) D. Ryndyk, M. Hartung, and G. Cuniberti, “Nonequilibrium molecular vibrons: An approach based on the nonequilibrium green function technique and the self-consistent born approximation,” Phys. Rev. B 73, 045420 (2006).
  * Benesch  _et al._ (2008) C. Benesch, M. Cizek, J. Klimeš, I. Kondov, M. Thoss, and W. Domcke, “Vibronic effects in single molecule conductance: First-principles description and application to benzenealkanethiolates between gold electrodes,” J. Phys. Chem. C 112, 9880–9890 (2008).
  * Härtle and Thoss (2011) R. Härtle and M. Thoss, “Resonant electron transport in single-molecule junctions: Vibrational excitation, rectification, negative differential resistance, and local cooling,” Phys. Rev. B 83, 115414 (2011).
  * Schinabeck, Härtle, and Thoss (2018) C. Schinabeck, R. Härtle, and M. Thoss, “Hierarchical quantum master equation approach to electronic-vibrational coupling in nonequilibrium transport through nanosystems: Reservoir formulation and application to vibrational instabilities,” Phys. Rev. B 97, 235429 (2018).
  * Erpenbeck  _et al._ (2016) A. Erpenbeck, R. Härtle, M. Bockstedte, and M. Thoss, “Vibrationally dependent electron-electron interactions in resonant electron transport through single-molecule junctions,” Phys. Rev. B 93, 115421 (2016).
  * Härtle and Kulkarni (2015) R. Härtle and M. Kulkarni, “Effect of broadening in the weak-coupling limit of vibrationally coupled electron transport through molecular junctions and the analogy to quantum dot circuit qed systems,” Phys. Rev. B 91, 245429 (2015).
  * Dzhioev and Kosov (2011) A. A. Dzhioev and D. Kosov, “Kramers problem for nonequilibrium current-induced chemical reactions,” J. Chem. Phys. 135, 074701 (2011).
  * Dzhioev, Kosov, and Von Oppen (2013) A. A. Dzhioev, D. S. Kosov, and F. Von Oppen, “Out-of-equilibrium catalysis of chemical reactions by electronic tunnel currents,” J. Chem. Phys. 138, 134103 (2013).
  * Pozner, Lifshitz, and Peskin (2014) R. Pozner, E. Lifshitz, and U. Peskin, “Charge transport-induced recoil and dissociation in double quantum dots,” Nano Lett. 14, 6244–6249 (2014).
  * Erpenbeck  _et al._ (2018a) A. Erpenbeck, C. Schinabeck, U. Peskin, and M. Thoss, “Current-induced bond rupture in single-molecule junctions,” Phys. Rev. B 97, 235452 (2018a).
  * Foti and Vázquez (2018) G. Foti and H. Vázquez, “Origin of vibrational instabilities in molecular wires with separated electronic states,” J. Phys. Chem. Lett. 9, 2791–2796 (2018).
  * Lu, Brandbyge, and Hedegård (2010) J.-T. Lu, M. Brandbyge, and P. Hedegård, “Blowing the fuse: Berry’s phase and runaway vibrations in molecular conductors,” Nano letters 10, 1657–1663 (2010).
  * Lü, Hedegård, and Brandbyge (2011) J.-T. Lü, P. Hedegård, and M. Brandbyge, “Laserlike vibrational instability in rectifying molecular conductors,” Phys. Rev. Lett. 107, 046801 (2011).
  * Lü  _et al._ (2012) J.-T. Lü, M. Brandbyge, P. Hedegård, T. N. Todorov, and D. Dundas, “Current-induced atomic dynamics, instabilities, and raman signals: Quasiclassical langevin equation approach,” Phys. Rev. B 85, 245444 (2012).
  * Preston, Kershaw, and Kosov (2020) R. J. Preston, V. F. Kershaw, and D. S. Kosov, “Current-induced atomic motion, structural instabilities, and negative temperatures on molecule-electrode interfaces in electronic junctions,” Phys. Rev. B 101, 155415 (2020).
  * Preston, Gelin, and Kosov (2021) R. J. Preston, M. F. Gelin, and D. S. Kosov, “First-passage time theory of activated rate chemical processes in electronic molecular junctions,” J. Chem. Phys. 154, 114108 (2021).
  * Domcke (1991) W. Domcke, “Theory of resonance and threshold effects in electron-molecule collisions: The projection-operator approach,” Phys. Rep. 208, 97–188 (1991).
  * Gertitschke and Domcke (1993) P. Gertitschke and W. Domcke, “Time-dependent wave-packet description of dissociative electron attachment,” Phys. Rev. A 47, 1031 (1993).
  * Čížek, Horáček, and Domcke (1999) M. Čížek, J. Horáček, and W. Domcke, “Associative detachment, dissociative attachment, and vibrational excitation of hcl by low-energy electrons,” Phys. Rev. A 60, 2873–2881 (1999).
  * Gallup and Fabrikant (2011) G. A. Gallup and I. I. Fabrikant, “Vibrational feshbach resonances in dissociative electron attachment to uracil,” Phys. Rev. A 83, 012706 (2011).
  * Brandbyge  _et al._ (1995) M. Brandbyge, P. Hedegård, T. Heinz, J. Misewich, and D. Newns, “Electronically driven adsorbate excitation mechanism in femtosecond-pulse laser desorption,” Phys. Rev. B 52, 6042 (1995).
  * Saalfrank (2006) P. Saalfrank, “Quantum dynamical approach to ultrafast molecular desorption from surfaces,” Chem. Rev. 106, 4116–4159 (2006).
  * Kim  _et al._ (2015) Y. Kim, K. Motobayashi, T. Frederiksen, H. Ueba, and M. Kawai, “Action spectroscopy for single-molecule reactions–experiments and theory,” Prog. Surf. Sci. 90, 85–143 (2015).
  * Frederiksen, Paulsson, and Ueba (2014) T. Frederiksen, M. Paulsson, and H. Ueba, “Theory of action spectroscopy for single-molecule reactions induced by vibrational excitations with stm,” Phys. Rev. B 89, 035427 (2014).
  * Erpenbeck  _et al._ (2020) A. Erpenbeck, Y. Ke, U. Peskin, and M. Thoss, “Current-induced dissociation in molecular junctions beyond the paradigm of vibrational heating: The role of antibonding electronic states,” Phys. Rev. B 102, 195421 (2020).
  * Halstead and Holloway (1990) D. Halstead and S. Holloway, “The influence of potential energy surface topologies on the dissociation of h2,” J. Chem. Phys. 93, 2859–2870 (1990).
  * Schinabeck  _et al._ (2016) C. Schinabeck, A. Erpenbeck, R. Härtle, and M. Thoss, “Hierarchical quantum master equation approach to electronic-vibrational coupling in nonequilibrium transport through nanosystems,” Phys. Rev. B 94, 201407 (2016).
  * Weiss (2012) U. Weiss,  _Quantum dissipative systems_ , Vol. 13 (World scientific, 2012).
  * Ilk and Makri (1994) G. Ilk and N. Makri, “Real time path integral methods for a system coupled to an anharmonic bath,” J. Chem. Phys. 101, 6708–6716 (1994).
  * Joutsuka and Ando (2011) T. Joutsuka and K. Ando, “Vibrational spectroscopy and relaxation of an anharmonic oscillator coupled to harmonic bath,” J. Chem. Phys. 134, 204511 (2011).
  * Härtle  _et al._ (2013) R. Härtle, G. Cohen, D. Reichman, and A. Millis, “Decoherence and lead-induced interdot coupling in nonequilibrium electron transport through interacting quantum dots: A hierarchical quantum master equation approach,” Phys. Rev. B 88, 235426 (2013).
  * Härtle  _et al._ (2015) R. Härtle, G. Cohen, D. Reichman, and A. Millis, “Transport through an anderson impurity: Current ringing, nonlinear magnetization, and a direct comparison of continuous-time quantum monte carlo and hierarchical quantum master equations,” Phys. Rev. B 92, 085430 (2015).
  * Xu  _et al._ (2017) M. Xu, L. Song, K. Song, and Q. Shi, “Convergence of high order perturbative expansions in open system quantum dynamics,” J. Chem. Phys. 146, 064102 (2017).
  * Trushechkin (2019) A. Trushechkin, “Higher-order corrections to the redfield equation with respect to the system-bath coupling based on the hierarchical equations of motion,” Lobachevskii J. Math. 40, 1606–1618 (2019).
  * Tanimura and Kubo (1989) Y. Tanimura and R. Kubo, “Time evolution of a quantum system in contact with a nearly gaussian-markoffian noise bath,” J. Phys. Soc. Jpn. 58, 101–114 (1989).
  * Tanimura (2006) Y. Tanimura, “Stochastic liouville, langevin, fokker–planck, and master equation approaches to quantum dissipative systems,” J. Phys. Soc. Jpn. 75, 082001 (2006).
  * Jin  _et al._ (2007) J. Jin, S. Welack, J. Luo, X.-Q. Li, P. Cui, R.-X. Xu, and Y. Yan, “Dynamics of quantum dissipation systems interacting with fermion and boson grand canonical bath ensembles: Hierarchical equations of motion approach,” J. Chem. Phys. 126, 134113 (2007).
  * Jin, Zheng, and Yan (2008) J. Jin, X. Zheng, and Y. Yan, “Exact dynamics of dissipative electronic systems and quantum transport: Hierarchical equations of motion approach,” J. Chem. Phys. 128, 234703 (2008).
  * Zheng  _et al._ (2009) X. Zheng, J. Jin, S. Welack, M. Luo, and Y. Yan, “Numerical approach to time-dependent quantum transport and dynamical kondo transition,” J. Chem. Phys. 130, 164708 (2009).
  * Yan (2014) Y. Yan, “Theory of open quantum systems with bath of electrons and phonons and spins: Many-dissipaton density matrixes approach,” J. Chem. Phys. 140, 054105 (2014).
  * Ye  _et al._ (2016) L. Ye, X. Wang, D. Hou, R.-X. Xu, X. Zheng, and Y. Yan, “Heom-quick: a program for accurate, efficient, and universal characterization of strongly correlated quantum impurity systems,” WIREs Comput Mol Sci 6, 608–638 (2016).
  * Wenderoth, Bätge, and Härtle (2016) S. Wenderoth, J. Bätge, and R. Härtle, “Sharp peaks in the conductance of a double quantum dot and a quantum-dot spin valve at high temperatures: A hierarchical quantum master equation approach,” Phys. Rev. B 94, 121303 (2016).
  * Dou  _et al._ (2018) W. Dou, C. Schinabeck, M. Thoss, and J. E. Subotnik, “A broadened classical master equation approach for treating electron-nuclear coupling in non-equilibrium transport,” J. Chem. Phys. 148, 102317 (2018).
  * Tanimura (2020) Y. Tanimura, “Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (heom),” J. Chem. Phys. 153, 020901 (2020).
  * Erpenbeck and Thoss (2019) A. Erpenbeck and M. Thoss, “Hierarchical quantum master equation approach to vibronic reaction dynamics at metal surfaces,” J. Chem. Phys. 151, 191101 (2019).
  * Bätge  _et al._ (2021) J. Bätge, Y. Ke, C. Kaspar, and M. Thoss, “Nonequilibrium open quantum systems with multiple bosonic and fermionic environments: A hierarchical quantum master equation approach,” arXiv preprint arXiv:2102.09484 (2021).
  * Hu, Xu, and Yan (2010) J. Hu, R.-X. Xu, and Y. Yan, “Communication: Padé spectrum decomposition of fermi function and bose function,” J. Chem. Phys. 133, 101106 (2010).
  * Hu  _et al._ (2011) J. Hu, M. Luo, F. Jiang, R.-X. Xu, and Y. Yan, “Padé spectrum decompositions of quantum distribution functions and optimal hierarchical equations of motion construction for quantum open systems,” J. Chem. Phys. 134, 244106 (2011).
  * Cui  _et al._ (2019) L. Cui, H.-D. Zhang, X. Zheng, R.-X. Xu, and Y. Yan, “Highly efficient and accurate sum-over-poles expansion of fermi and bose functions at near zero temperatures: Fano spectrum decomposition scheme,” J. Chem. Phys. 151, 024110 (2019).
  * Abe, Yamashita, and Saalfrank (2003) A. Abe, K. Yamashita, and P. Saalfrank, “Stm and laser-driven atom switch: An open-system density-matrix study of h/si (100),” Phys. Rev. B 67, 235411 (2003).
  * Tang  _et al._ (2015) Z. Tang, X. Ouyang, Z. Gong, H. Wang, and J. Wu, “Extended hierarchy equation of motion for the spin-boson model,” J. Chem. Phys. 143, 224112 (2015).
  * Rahman and Kleinekathöfer (2019) H. Rahman and U. Kleinekathöfer, “Chebyshev hierarchical equations of motion for systems with arbitrary spectral densities and temperatures,” J. Chem. Phys. 150, 244104 (2019).
  * Erpenbeck  _et al._ (2018b) A. Erpenbeck, C. Hertlein, C. Schinabeck, and M. Thoss, “Extending the hierarchical quantum master equation approach to low temperatures and realistic band structures,” J. Chem. Phys. 149, 064106 (2018b).
  * Ishizaki and Tanimura (2005) A. Ishizaki and Y. Tanimura, “Quantum dynamics of system strongly coupled to low-temperature colored noise bath: Reduced hierarchy equations approach,” J. Phys. Soc. Jpn. 74, 3131–3134 (2005).
  * Shi  _et al._ (2009) Q. Shi, L. Chen, G. Nan, R.-X. Xu, and Y. Yan, “Efficient hierarchical liouville space propagator to quantum dissipative dynamics,” J. Chem. Phys. 130, 084105 (2009).
  * Xu  _et al._ (2019) M. Xu, Y. Liu, K. Song, and Q. Shi, “A non-perturbative approach to simulate heterogeneous electron transfer dynamics: Effective mode treatment of the continuum electronic states,” J. Chem. Phys. 150, 044109 (2019).
  * Colbert and Miller (1992) D. T. Colbert and W. H. Miller, “A novel discrete variable representation for quantum mechanical reactive scattering via the s-matrix kohn method,” J. Chem. Phys. 96, 1982–1991 (1992).
  * Echave and Clary (1992) J. Echave and D. C. Clary, “Potential optimized discrete variable representation,” Chem. Phys. Lett. 190, 225–230 (1992).
  * Seideman and Miller (1992) T. Seideman and W. H. Miller, “Calculation of the cumulative reaction probability via a discrete variable representation with absorbing boundary conditions,” J. Chem. Phys. 96, 4412–4422 (1992).
  * Riss and Meyer (1996) U. V. Riss and H. Meyer, “Investigation on the reflection and transmission properties of complex absorbing potentials,” J. Chem. Phys. 105, 1409–1419 (1996).
  * Selstø and Kvaal (2010) S. Selstø and S. Kvaal, “Absorbing boundary conditions for dynamical many-body quantum systems,” J. Phys. B 43, 065004 (2010).
  * Kvaal (2011) S. Kvaal, “Multiconfigurational time-dependent hartree method to describe particle loss due to absorbing boundary conditions,” Phys. Rev. A 84, 022512 (2011).
  * Prucker  _et al._ (2018) V. Prucker, M. Bockstedte, M. Thoss, and P. Coto, “Dynamical simulation of electron transfer processes in self-assembled monolayers at metal surfaces using a density matrix approach,” J. Chem. Phys. 148, 124705 (2018).
  * Wang, Nian, and Lü (2020) T. Wang, L.-L. Nian, and J.-T. Lü, “Nonthermal vibrations in biased molecular junctions,” Phys. Rev. E 102, 022127 (2020).
  * Zhang, Zheng, and Di Ventra (2019) D. Zhang, X. Zheng, and M. Di Ventra, “Local temperatures out of equilibrium,” Phys. Rep. 830, 1–66 (2019).
  * Wilkins and Dattani (2015) D. M. Wilkins and N. S. Dattani, “Why quantum coherence is not important in the fenna–matthews–olsen complex,” J. Chem. Theory Comput. 11, 3411–3419 (2015).
  * Lee, Sorescu, and Deng (2011) J. Lee, D. C. Sorescu, and X. Deng, “Electron-induced dissociation of co2 on tio2 (110),” J. Am. Chem. Soc. 133, 10066–10069 (2011).
  * Tan  _et al._ (2011) S. Tan, Y. Zhao, J. Zhao, Z. Wang, C. Ma, A. Zhao, B. Wang, Y. Luo, J. Yang, and J. Hou, “Co 2 dissociation activated through electron attachment on the reduced rutile tio 2 (110)-1$\times$ 1 surface,” Phys. Rev. B 84, 155418 (2011).
  * Zhao  _et al._ (2013) A. Zhao, S. Tan, B. Li, B. Wang, J. Yang, and J. Hou, “Stm tip-assisted single molecule chemistry,” Phys. Chem. Chem. Phys. 15, 12428–12441 (2013).
  * Chen  _et al._ (2019) C. Chen, L. Kong, Y. Wang, P. Cheng, B. Feng, Q. Zheng, J. Zhao, L. Chen, and K. Wu, “Dynamics of single-molecule dissociation by selective excitation of molecular phonons,” Phys. Rev. Lett. 123, 246804 (2019).
  * Salam, Persson, and Palmer (1994) G. P. Salam, M. Persson, and R. E. Palmer, “Possibility of coherent multiple excitation in atom transfer with a scanning tunneling microscope,” Phys. Rev. B 49, 10655–10662 (1994).
  * Ueba (2003) H. Ueba, “Motions and reactions of single adsorbed molecules induced by vibrational excitation with stm,” Surf. Rev. Lett. 10, 771–796 (2003).
  * Tikhodeev and Ueba (2004) S. Tikhodeev and H. Ueba, “Relation between inelastic electron tunneling and vibrational excitation of single adsorbates on metal surfaces,” Phys. Rev. B 70, 125414 (2004).
  * Mitra, Aleiner, and Millis (2004) A. Mitra, I. Aleiner, and A. Millis, “Phonon effects in molecular transistors: Quantal and classical treatment,” Phys. Rev. B 69, 245302 (2004).
  * Härtle, Benesch, and Thoss (2009) R. Härtle, C. Benesch, and M. Thoss, “Vibrational nonequilibrium effects in the conductance of single molecules with multiple electronic states,” Phys. Rev. Lett. 102, 146801 (2009).
  * Brisker and Peskin (2008) D. Brisker and U. Peskin, “Charge-transport-induced dissociation in donor-bridge-acceptor complexes,” J. Chem. Phys. 129, 244709 (2008).
  * Härtle, Peskin, and Thoss (2013) R. Härtle, U. Peskin, and M. Thoss, “Vibrationally coupled electron transport in single-molecule junctions: The importance of electron–hole pair creation processes,” Physica Status Solidi (b) 250, 2365–2377 (2013).
  * Nitzan and Galperin (2018) A. Nitzan and M. Galperin, “Kinetic schemes in open interacting systems,” J. Phys. Chem. Lett. 9, 4886–4892 (2018).
  * Xin  _et al._ (2019) N. Xin, J. Guan, C. Zhou, X. Chen, C. Gu, Y. Li, M. A. Ratner, A. Nitzan, J. F. Stoddart, and X. Guo, “Concepts in the design and engineering of single-molecule electronic devices,” Nature Reviews Physics 1, 211–230 (2019).
  * Sun  _et al._ (2018) H. Sun, Z. Jiang, N. Xin, X. Guo, S. Hou, and J. Liao, “Efficient fabrication of stable graphene-molecule-graphene single-molecule junctions at room temperature,” ChemPhysChem 19, 2258–2265 (2018).
  * Leitherer, Papior, and Brandbyge (2019) S. Leitherer, N. Papior, and M. Brandbyge, “Current-induced atomic forces in gated graphene nanoconstrictions,” Phys. Rev. B 100, 035415 (2019).
  * Brisker and Peskin (2006) D. Brisker and U. Peskin, “Vibrational anharmonicity effects in electronic tunneling through molecular bridges.” J. Chem. Phys. 125, 111103 (2006).
  * Arielly  _et al._ (2017) R. Arielly, N. Nachman, Y. Zelinskyy, V. May, and Y. Selzer, “Picosecond time resolved conductance measurements of redox molecular junctions,” J. Chem. Phys. 146, 092306 (2017).
  * Koch, Von Oppen, and Andreev (2006) J. Koch, F. Von Oppen, and A. Andreev, “Theory of the franck-condon blockade regime,” Phys. Rev. B 74, 205438 (2006).
  * Schinabeck  _et al._ (2014) C. Schinabeck, R. Härtle, H. Weber, and M. Thoss, “Current noise in single-molecule junctions induced by electronic-vibrational coupling,” Phys. Rev. B 90, 075409 (2014).
  * Fung  _et al._ (2019) E.-D. Fung, D. Gelbwaser, J. Taylor, J. Low, J. Xia, I. Davydenko, L. M. Campos, S. Marder, U. Peskin, and L. Venkataraman, “Breaking down resonance: Nonlinear transport and the breakdown of coherent tunneling models in single molecule junctions,” Nano letters 19, 2555–2561 (2019).
  * Zang  _et al._ (2020) Y. Zang, E.-D. Fung, T. Fu, S. Ray, M. H. Garner, A. Borges, M. L. Steigerwald, S. Patil, G. Solomon, and L. Venkataraman, “Voltage-induced single-molecule junction planarization,” Nano Lett. 21, 673–679 (2020).
  * Hofmeister, Coto, and Thoss (2017) C. Hofmeister, P. B. Coto, and M. Thoss, “Controlling the conductance of molecular junctions using proton transfer reactions: A theoretical model study,” J. Chem. Phys. 146, 092317 (2017).
  * Weckbecker, Coto, and Thoss (2017) D. Weckbecker, P. Coto, and M. Thoss, “Controlling the conductance of a graphene–molecule nanojunction by proton transfer,” Nano Lett. 17, 3341–3346 (2017).
  * Shi  _et al._ (2018) Q. Shi, Y. Xu, Y. Yan, and M. Xu, “Efficient propagation of the hierarchical equations of motion using the matrix product state method,” J. Chem. Phys. 148, 174102 (2018).
  * Borrelli (2019) R. Borrelli, “Density matrix dynamics in twin-formulation: An efficient methodology based on tensor-train representation of reduced equations of motion,” J. Chem. Phys. 150, 234102 (2019).
  * Yan  _et al._ (2021) Y. Yan, M. Xu, T. Li, and Q. Shi, “Efficient propagation of the hierarchical equations of motion using the tucker and hierarchical tucker tensors,” J. Chem. Phys. 154, 194104 (2021).