# Dissipative Tunneling Rates through the Incorporation of First-Principles
Electronic Friction in Instanton Rate Theory II: Benchmarks and Applications

Y. Litman<EMAIL_ADDRESS>MPI for the Structure and Dynamics of Matter,
Luruper Chaussee 149, 22761 Hamburg, Germany  E. S. Pós  MPI for the Structure
and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany  C. L.
Box  Department of Chemistry, University of Warwick, Coventry CV4 7AL, United
Kingdom  R. Martinazzo  Department of Chemistry, Università degli Studi di
Milano, Via Golgi 19, 20133 Milano, Italy  R. J. Maurer  Department of
Chemistry, University of Warwick, Coventry CV4 7AL, United Kingdom  M. Rossi
<EMAIL_ADDRESS>MPI for the Structure and Dynamics of Matter,
Luruper Chaussee 149, 22761 Hamburg, Germany

###### Abstract

In part I, we presented the ring-polymer instanton with explicit friction
(RPI-EF) method and showed how it can be connected to the ab initio electronic
friction formalism. This framework allows the calculation of tunneling
reaction rates that incorporate the quantum nature of the nuclei and certain
types of non-adiabatic effects (NAEs) present in metals. In this second part,
we analyze the performance of RPI-EF on model potentials and apply it to
realistic systems. For a 1D double-well model, we benchmark the method against
numerically exact results obtained from multi-layer multi-configuration time-
dependent Hartree calculations. We demonstrate that RPI-EF is accurate for
medium and high friction strengths and less accurate for extremely low
friction values. We also show quantitatively how the inclusion of NAEs lowers
the cross-over temperature into the deep tunneling regime, reduces the
tunneling rates, and in certain regimes, steers the quantum dynamics by
modifying the tunneling pathways. As a showcase of the efficiency of this
method, we present a study of hydrogen and deuterium hopping between
neighboring interstitial sites in selected bulk metals. The results show that
multidimensional vibrational coupling and nuclear quantum effects have a
larger impact than NAEs on the tunneling rates of diffusion in metals.
Together with part I, these results advance the calculations of dissipative
tunneling rates from first principles.

††preprint: AIP/123-QED

##  I Introduction

The transport kinetics of small molecules and atoms in metals play an
important role in technological applications in several areas, like fuel
cells, batteries, and nuclear reactors, among others Mohtadi and Orimo (2016);
Robertson  _et al._ (2015); Traidia, Chatzidouros, and Jouiad (2018);
Schlapbach and Züttel (2001). In particular, accurate measurements of hydrogen
diffusion are still experimentally challenging since diffusion coefficients
are very sensitive to the microstructure of the material and the composition
of the alloy Völkl and Alefeld (1978); Jothi, Croft, and Brown (2015); Fuk
(2005). As a consequence, reported values of diffusion constants by different
groups can be scattered over a few orders of magnitude.

First-principles calculations have the potential to provide an increased
understanding of the transport dynamics of light particles in well-defined
structures throughout a wide range of thermodynamic conditions (e.g.,
temperature and pressure), which can complement experiments and guide further
developments. Most of the atomistic calculations performed to study nuclear
dynamical processes rely on Newtonian dynamics on the ground-state adiabatic
potential energy surface (PES), as given by the Born-Oppenheimer approximation
(BOA). However, the coupling of the nuclear movement with electrons in the
metal can easily induce electronic excitations, which represent a breakdown of
the BOA and give rise to non-adiabatic effects (NAEs) Nienhaus (2002); Wodtke
(2016). Furthermore, light particles such as hydrogen and deuterium can
exhibit strong nuclear quantum effects (NQEs), which can either increase or
decrease its mobility through metals Kimizuka, Ogata, and Shiga (2018);
Kimizuka and Shiga (2021). As a consequence, an understanding of the interplay
of NQEs and NAEs in the transport process of light atoms in realistic systems
remains elusive.

Several exact methods to simulate non-adiabatic quantum dynamics have been
developed in the last decadesBeck  _et al._ (2000a); Tanimura and Kubo (1989);
Ilk and Makri (1994). However, due to the unfavorable scaling of these
theoretical approaches with the number of degrees of freedom, efficient but
accurate methods are required to capture NQEs and NAEs in high-dimensional
systems or to be used with costly ab initio potentials. In part I of this
paper, we presented the ring-polymer instanton with explicit friction (RPI-EF)
theory for the calculation of dissipative thermal tunneling rates, which has
the potential to fulfill these requirements. Briefly, RPI-EF allows an easy
and efficient incorporation of the electronic friction formalism, originally
proposed by Hellsing and Persson Hellsing and Persson (1984) and Head-Gordon
and TullyHead-Gordon and Tully (1995), into the the semi-classical ring-
polymer instanton (RPI) rate theoryRichardson (2018).

In this article, we show the merits and limitations of RPI-EF. We benchmark
the accuracy of RPI-EF rate predictions by comparing them with numerically
exact theories in model potentials. Subsequently, we examine the interplay of
the reaction barrier height and friction strengths on the tunneling rate in 1D
and 2D model potentials connected to a bath. Finally, the new approach is
applied to hydrogen and deuterium hopping reactions in bulk transition metals,
focusing on Pd, by employing Kohn-Sham density functional theory. These
calculations allow a quantitative analysis of the impact of different effects
on the rate constants, such as the dimensionality of the system, NQEs, and
NAEs.

Part II of this paper is structured as follows: In Section II the methods
employed to obtain the tunneling rates are briefly summarized and the
simulation details for each one are specified. In Section III, the model
potentials and the systems treated from first-principles are described.
Results of numerical simulations on low-dimensional models for position-
independent and position-dependent friction tensors are discussed in Section
IV.1 and IV.2, respectively. Finally, the first-principles results for the
hydrogen and deuterium hopping in metals are discussed in Section V. Section
VI concludes by summarizing the main results and giving an outlook to future
directions.

##  II Rate Calculation Methods

###  II.1 Ring-polymer instanton (with explicit friction) calculations

The RPI rate theory Richardson and Althorpe (2009a); Arnaldsson (2007) is a
semi-classical method that allows the calculation of tunneling rates. RPI
theory can be interpreted as the extension of Eyring transition state theory
(TST) Eyring (1935) into the deep tunneling regime, since the rates are
evaluated utilizing only a limited amount (often a single) special
configuration, circumventing the necessity of real-time sampling. While in TST
this special configuration is the first-order saddle-point connecting
reactants and products on the potential energy surface (PES), in RPI rate
theory the special configurations are found at the first-order saddle-points
of the extended space of the ring polymer (RP) potential. These trajectories
in imaginary time are known as instanton trajectories. The RPI rate expression
is analogous to the one proposed by TST and reads,

$\begin{split}k_{\text{inst}}(\beta)\propto
e^{-U_{P}^{\text{sys}}(\bm{\bar{q}})/\hbar},\end{split}$ (1)

where

$U_{P}^{\text{sys}}(\bm{q})=\\\
\sum_{k=1}^{P}\sum_{i=1}^{3N}m_{i}\frac{\omega_{P}^{2}}{2}(q_{i}^{(k)}-q_{i}^{(k+1)})^{2}+\\\
\sum_{k=1}^{P}V(q_{1}^{(k)},\dots q_{3N}^{(k)})$ (2)

is the RP potential. In the previous equation, $q^{(k)}_{i}$ is the position
of the $i$-th degree of freedom of the $k$-th bead of the ring polymer,
$m_{i}$ is the mass of the $i$-th degree of freedom, $N$ is the number of
atoms, $P$ is the number of beads (replicas), $\bm{q}$ is an abbreviated
notation to represent all the degrees of freedom, $\bm{\bar{q}}$ denotes the
instanton geometry, and $\omega_{P}=(\beta_{P}\hbar)^{-1}$ with
$\beta_{P}=\frac{1}{k_{B}PT}$.

In part I, we showed how RPI rate theory can be extended to compute tunneling
rates for systems connected to a harmonic bath which simulates a dissipative
environment. Irrespective of the dissipative mechanism, and assuming that the
environment degrees of freedom adjust adiabatically to the system position,
one can fully characterize the system-environment coupling by a position
($\bm{q}$) and frequency ($\lambda$) dependent friction tensor,
$\tilde{\eta}(\bm{q},\lambda)$. Moreover, when the frequency and position
dependence are decoupled, we proved that the RP potential that enters Eq. 1 is
renormalized, adopting the form

$\begin{split}U^{\text{MF}}_{P}=&U_{P}^{\text{sys}}+\\\
&\sum_{l=-P/2+1}^{P/2}\sum_{i=1}^{3N}\frac{\omega_{l}}{2}\bigg{[}\sum_{k=1}^{P}C_{lk}\bigg{(}\int_{\bm{q}^{\text{ref}}}^{\bm{q}^{(k)}}\tilde{\eta}(\bm{q}^{\prime},\omega_{k})^{1/2}\cdot
d\bm{q}^{\prime}\bigg{)}\bigg{]}^{2},\end{split}$ (3)

where $\omega_{l}=2\omega_{P}\sin(|l|\pi/P)$ are the free RP normal mode
frequencies and $\bm{C}$ is the transformation matrix between the RP normal
modes and Cartesian coordinates. In the limiting case of a position-
independent friction, the previous expression can be simplified to

$\begin{split}U^{\text{MF}}_{P}=U_{P}^{\text{sys}}+\sum_{l=-P/2+1}^{P/2}\sum_{i=1}^{3N}\frac{\tilde{\eta}(\omega_{l}){\omega_{l}}}{2}(Q_{i}^{(l)})^{2},\\\
\end{split}$ (4)

where $\bm{Q}^{(l)}$ represent the free RP normal mode coordinates.

The RPI calculations were performed using the i-PI Kapil  _et al._ (2019)
code. The forces and energies required by the algorithm were passed to i-PI
from an external code through an interface based on internet sockets. The RPI-
EF calculations required an extension of the i-PI communication protocol. We
added to the existing communication of the energy, forces, and stresses from
external codes to i-PI, the possibility to pass additional information as
JSON-formatted strings. In this way, it is possible to communicate any type of
data and, more importantly, when the data is numeric, it becomes available to
be used by any implemented algorithm within i-PI. This enables the use of
quantities beyond energies and forces, that change along the simulation,
within different types of dynamics.

The FHI-aims codeBlum  _et al._ (2009) and an in-house python code were used
in connection to i-PI for the DFT and model calculations, respectively.
Transition-state geometries were obtained either by the string methodE, Ren,
and Vanden-Eijnden (2007) combined with the climbing image techniqueHenkelman,
Uberuaga, and Jónsson (2000) or using a minimum-mode-following algorithm
Nichols  _et al._ (1990), as implemented in i-PI. The RPI calculations were
initialized after finding the transition state, by stretching the transition-
state geometry along the mode with imaginary frequency using a number of
replicas between 10 and 16. Optimizations were started at a temperature of 10
K below the corresponding cross-over temperature, $T_{c}^{\circ}$. After
converging the instanton pathway for the first calculation, successive steps
of temperature decrease and RP interpolation to increase the number of beads
were performed until the target temperature was reached. If required, further
calculations with more beads were performed to guarantee that, in all cases
and for all temperatures, the final rates were converged within a 10% error
Beyer  _et al._ (2016); Litman (2020). See more details in section I of the
Supplemental Information (SI).

###  II.2 Multi-configuration time-dependent Hartree Calculations

Numerically exact results for selected models in this paper were obtained with
the multi-layer variantWang and Thoss (2003); Manthe (2009); Vendrell and
Meyer (2011) of the multi-configuration time-dependent Hartree methodMeyer,
Manthe, and Cederbaum (1990); Beck  _et al._ (2000b); Meyer, Gatti, and Worth
(2009) (ML-MCTDH), as implemented in the Heidelberg packageWorth  _et al._ .
MCTDH is a variational method that relies on optimal, time-dependent basis
functions to alleviate the exponential scaling problem of standard methods
based on direct expansions on time-independent basis. In the MCTDH _ansatz_
the wave function $\Ket{\Psi(t)}$ is expanded on orthogonal configurations
$\Ket{\Phi_{J}(t)}$, which in turn are products of ‘single particle’ functions
(SPFs),

$\begin{split}\Ket{\Psi(t)}=\sum_{J}A_{J}(t)\Ket{\Phi_{J}(t)}=\sum_{J}A_{J}(t)\prod_{k}\ket{\phi_{j_{k}}^{(k)}}(t),\end{split}$
(5)

and both the expansion coefficients ($A_{J}$) and the SPFs
($\ket{\phi_{j}(t)}$) are variationally optimized. Here,
$J=(j_{1},..j_{k},..j_{F})$ is a multi-index, the index $k=1,2,..F$ runs over
the single-particles, and $j_{k}=1,..n_{k}$ labels the SPFs used for the
$k^{\textrm{th}}$ mode. Conventional MCTDH uses single-particles for each
degree of freedom or small groups thereof, and represents their SPFs by a
direct expansion on a grid/basis-set (the so-called primitive grid) designed
for the single particle at hand. This limits its capability of handling large
systems. In contrast, in ML-MCTDH the single particles are high-dimensional
modes whose SPFs are described by further MCTDH expansions employing lower
dimensional SPFs. The procedure can be indefinitely iterated till reasonably
small single-particles are defined that can be described on primitive grids.
This recipe, similarly to tensor networks and matrix-product statesCirac  _et
al._ (2021), endows the wave function with a hierarchical, flexible structure
that allows the treatment of considerably larger systems. In particular, ML-
MCTDH has been successfully applied to the calculation of thermal rate
constants in condensed-phase problemsWang, Skinner, and Thoss (2006a); Craig,
Thoss, and Wang (2007a) in the framework of the reactive flux-side approach,
where $k(T)$ is given by the long-time limit of the equilibrium flux-side
time-correlation functionMiller, Schwartz, and Tromp (1983); Miller (1998).

In this work, we followed closely the original work by Wang and ThossWang,
Skinner, and Thoss (2006b); Craig, Thoss, and Wang (2007b), who introduced an
importance sampling technique to recast the trace expression of the flux-side
correlation function as an accessible ensemble-average over time-evolving
wavepackects. In a nutshell, the evaluation of the rate constant is reduced
to: i) a preparation step where the wavepackets are initialized by combining
(system) Boltzmannized flux eigenvectors with bath states drawn from the
canonical ensemble of the uncoupled bath, ii) an equilibration step where
imaginary-time dynamics introduces the correlations present in the coupled
system, and iii) a propagation step where the real-time dynamics is followed
up to the onset of the kinetic regime (the long-time limit alluded to above).
Details about the calculations, including an overview of the flux-side
approach, the Monte Carlo sampling, the tree structure of the ML-MCTDH
wavefunction, the number of SPFs and the primitive grids used, are provided in
section II of SI. Converged calculations were obtained with 50 bath modes and
using 128-256 realizations for each value of $T$ and coupling strength.

##  III Simulation details and parameters

###  III.1 1D and 2D double well models

We analyze the performance of RPI-EF on a double-well model similar to the
ones usually employed to study quantum dynamics in system-bath models Topaler
and Makri (1994); Craig and Manolopoulos (2005). The potential energy surface
of the system is given by

$\begin{split}V_{\text{DW}}(q)=-\frac{1}{2}m\omega^{{\ddagger}2}(q-q_{0})^{2}+\frac{m^{2}\omega^{{\ddagger}4}}{16V_{0}}q^{4},\end{split}$
(6)

where, unless otherwise specified, we set $m$ as the mass of atomic hydrogen
and $\omega^{\ddagger}=500$ cm-1. The coupling between the system and the bath
can be made position-dependent according toStraus, Gomez Llorente, and Voth
(1993)

$\begin{split}f(q)=q[1+\epsilon_{1}\exp(-\Delta
q^{2}/2)+\epsilon_{2}\tanh(\Delta q)],\end{split}$ (7)

where $\Delta q=(q-q^{\ddagger})/\delta$, $\delta$ determines the length-scale
of the nonlinear couplings, and $\epsilon_{1}$ and $\epsilon_{2}$ the
magnitude of its symmetric and anti-symmetric components, respectively. The
calculations with position-independent friction were obtained by setting
$\epsilon_{1}=\epsilon_{2}=0$. Naturally, the position-dependent couplings of
real systems do not generally adopt such simple forms. However, as it will be
shown in Sec. IV.2, this simplified form will prove sufficient to expose the
importance of including this position-dependence in the tunneling rates.

We consider an Ohmic (linear) spectral density multiplied by an exponential
cutoff, leading to

$\begin{split}\tilde{\eta}_{il}(\bm{q},\lambda)=\frac{2}{\pi}\int_{0}^{+\infty}d\omega\tilde{\eta}_{0}\bigg{(}\frac{\partial
f(\bm{q})_{i}}{\partial
q_{l}}\bigg{)}^{2}\frac{\lambda}{\omega^{2}+\lambda^{2}}e^{-\omega/\omega_{c}}\end{split}$
(8)

where $\tilde{\eta}_{0}$ is the static friction coefficient and, unless
otherwise specified, we set $\omega_{c}=500$ cm-1.

To illustrate the tensorial nature of the friction we also considered a two-
dimensional double double-well (DDW) potential given by

$\begin{split}V_{\text{DDW}}(q_{1},q_{2})=V_{\text{DW}}(q_{1})+V_{\text{DW}}(q_{2})+Cq_{1}q_{2},\end{split}$
(9)

where $C$ is a constant to be specified. Additionally, for the sake of
simplicity we considered that the bath couples to each degree of freedom
independently, and the coupling function is given by

$\begin{split}f(q_{1},q_{2})=(f(q_{1}),f(q_{2}))\end{split}$ (10)

where $f(q_{i})$ is given by Eq. 7.

###  III.2 Fcc metals

We also perform first-principles atomistic simulations to showcase the
methodology developed in this paper. More specifically, we focus on hydrogen
and deuterium hopping reactions within interstitial sites in different bulk
fcc metals: Pd, Pt, Cu, and Ag. The bulk systems were modelled by 2 $\times$ 2
$\times$ 2 cubic supercells containing one hydrogen or deuterium atom and 32
metal atoms. Energies and forces were computed employing density-functional
theory (DFT) using the FHI-aims Blum  _et al._ (2009) code and the Perdew,
Burke, and Ernzerhof (PBE) exchange-correlation functional Perdew, Burke, and
Ernzerhof (1996). Geometries were relaxed using the standard light settings
(in the case of Pd, increasing the radial multiplier to 2) from FHI-aims until
all forces were below 10-3 eV/Å. Minimum-energy pathways (MEP) were obtained
with the string methodE, Ren, and Vanden-Eijnden (2007) combined with the
climbing image techniqueHenkelman, Uberuaga, and Jónsson (2000) as implemented
in the aimsChain package provided with the FHI-aims code. The BFGS algorithm
was used as the optimization procedure and the residual forces converged below
10-3 eV/Å. Unless specified otherwise, a 6 $\times$ 6 $\times$ 6 k-point
sampling was used. This setup ensures that errors in relative energies are
below 1 meV/atom. We obtained lattice constants of 3.95 Å, 3.97 Å, 3.63 Å, and
4.16 Å, for Pd, Pt, Cu and Ag, respectively, in good agreement with Ref. Haas,
Tran, and Blaha (2009). See section III A of the SI for more details regarding
convergence tests.

As explained in part I of this paper, the electronic friction tensor was
computed as

$\begin{split}\tilde{\eta}_{ij}(\textbf{q},\lambda)&=\hbar\sum_{\nu,\nu^{\prime}}\braket{\psi_{\nu}}{\partial_{i}\psi_{\nu^{\prime}}}\braket{\psi_{\nu^{\prime}}}{\partial_{j}\psi_{\nu}}(f(\epsilon_{\nu})-f(\epsilon_{\nu^{\prime}}))\\\
&\times\frac{\lambda\Omega_{\nu\nu^{\prime}}}{\lambda^{2}+\Omega_{\nu\nu^{\prime}}^{2}},\end{split}$
(11)

where $f(\epsilon)$ is the state occupation given by the Fermi-Dirac
distribution, $\psi_{\nu}$ and $\epsilon_{\nu}$ are the KS electronic orbitals
and orbital energies of the $\nu$-th level, $i$ and $j$ label the nuclear
degrees of freedom, $\partial_{i}=\partial/\partial q_{i}$, and
$\Omega_{\nu\nu^{\prime}}=(\epsilon_{\nu}-\epsilon_{\nu^{\prime}})/\hbar$. The
calculation of the non-adiabatic coupling elements was obtained through a
finite-difference approach, as currently implemented in the FHI-aims code
Maurer  _et al._ (2016). We stress that Eq. 11 is different from the
expression commonly used in the literature Dou, Miao, and Subotnik (2017);
Maurer  _et al._ (2016); Head-Gordon and Tully (1995); Luntz and Persson
(2005); Monturet and Saalfrank (2010). We used a step length of 0.001 Å for
the finite-difference evaluation, and a 16 $\times$ 16 $\times$ 16 k-point
sampling for the friction tensor. We only calculated the tensor components
related to the hydrogen or deuterium atoms. More details regarding convergence
tests can be found in section III. B of the SI.

##  IV Numerical results on model potentials

###  IV.1 Position-independent friction

Figure 1:  Reaction rate constants for the DW potential ($V_{0}=258$ meV,
$\omega_{c}=500$ cm-1 and $q_{0}=0$) for temperatures between 50 K and 100 K
for different friction values computed with RPI-EF (solid lines with filled
circles) and ML-MCTDH (dashed lines with empty circles).

We start by analyzing the linear coupling case and benchmarking the RPI-EF
results against the ML-MCTDH results. In Fig. 1 the rate constants calculated
for the DW potential (Eq. 6) at different temperatures and friction values are
shown. The RPI-EF results are in good agreement with the exact calculations
for $\tilde{\eta}_{0}/m\omega^{\ddagger}>0.1$ at all the temperatures
considered. We note that the lowest temperature, 50 K, represents less than
half of the cross-over temperature evaluated without friction. Even at
considerably higher friction values, where the RPI approach has been predicted
to be inadequate Richardson and Althorpe (2009b), the agreement is quite
remarkable, showing that real-time dynamical effects, such as recrossing, play
a minor role in these cases. At lower friction values, specially below
$\tilde{\eta}_{0}/m\omega^{\ddagger}=0.05$, the agreement deteriorates. The
poor performance of the RPI-EF method in the weak coupling regime is not
surprising since RPI predicts a finite value for the rate even at
$\tilde{\eta}_{0}=0$, where the dynamics in such a 1D model would be described
by a Rabi oscillation and, strictly speaking, a rate process cannot be
defined. Indeed, RPI-EF reaches a plateau at
$\tilde{\eta}_{0}/m\omega^{\ddagger}\approx 0.05$ which might be interpreted
as “an intrinsic” dissipation inherent to that theory, as a consequence of
considering trajectories that bounce only once in the evaluation of the
imaginary time kernel Richardson (2018). When $\eta$ goes to zero, the exact
results approach the limit of coherent tunneling dynamics, yielding ML-MCTDH
rates that are larger than the ones calculated with RPI. In passing, we note
that both methods present a minimum of the rate at around 70 K in the weak-
friction regime (see inset in Fig. 1), which differs from the low-temperature
power law observed in metastable systems Grabert, Weiss, and Hanggi (1984);
Grabert, Olschowski, and Weiss (1987). A deeper study of this subtle but
interesting quantum effect is beyond the scope of the current work and will be
the subject of future research.

Figure 2:  Tunneling enhancement factors
$\log_{10}[\kappa(\beta,\tilde{\eta}_{0})]$ as a function of energy barrier
height and friction strength, on the DW potential with $q_{0}=0.0$ Å
(symmetric reaction profile) at $T=0.7T^{\circ}_{c}=80$ K. The colour-scale is
logarithmic and contour lines are drawn for isosurfaces spaced by 2
logarithmic units.

One way to evaluate the relative importance of tunneling to the total rate is
to analyze the tunneling enhancement factor, $\kappa^{\text{tun}}$, defined as

$\begin{split}\kappa^{\text{tun}}(\beta,\tilde{\eta})=\frac{k_{\text{inst}}(\beta,\tilde{\eta})}{k_{\text{TST}}(\beta,\tilde{\eta})}\end{split}$
(12)

with $k_{\text{TST}}(\beta,\eta)$ being the TST rateBeyer  _et al._ (2016).
Figs. 2 and 3 show the calculated tunneling enhancement factors for different
barrier heights and friction values for symmetric and asymmetric reactions,
respectively. For the range of parameters considered here, it can be observed
that the the tunneling enhancement factor calculated with and without friction
can differ up to almost ten orders of magnitude, that it increases with the
increase of the barrier height, and that it decreases with increasing friction
strength. For barriers larger than 500 meV, the isosurfaces are approximately
straight lines with slope one, meaning that barrier heights and friction
strengths have a comparable but opposite effect on tunneling.

Figure 3:  Same as Fig. 2 with $q_{0}=0.08$Å (asymmetric reaction profile).
The exoergic reaction is considered.

We now proceed to discuss why the friction strength and the barrier height
impact the tunneling contribution to the rate. The friction value determines
the system-bath coupling strength. As evidenced by $\kappa^{\text{tun}}\sim 1$
at higher friction values, the stronger the coupling, the more classical the
system behaves. However, the reason why the impact of the friction becomes
more relevant at higher barrier heights is less straightforward to understand
and requires the analysis of the instanton pathways. In Fig. 4, we show the
decomposition of the instanton geometry into the free RP normal mode basis. We
first consider the case without coupling to the bath. For symmetric barriers,
Fig. 4a, the instanton pathway expands only along the odd RP normal modes due
to the symmetry of the underlying potential, and the first two degenerate RP
normal modes ($l=\pm 1$) contribute with more than 99$\%$ to the path. For
asymmetric barriers, Fig. 4c, even though all normal modes are in principle
allowed by symmetry, the first two degenerate RP normal modes exhibit the
highest contribution, with the centroid mode ($l=0$) presenting a non-
negligible contribution as well. For both barrier shapes, the population of
the $l=\pm 1$ modes increases with the barrier height, simply because the
pathway from reactants to product becomes longer. The same trends are observed
for calculations with intermediate system-bath coupling, Fig. 4c and 4d, where
the only difference is an overall smaller population of the RP normal modes
due to the shorter instanton pathway. Thus, higher barriers correlate with a
relatively larger impact of the friction on the rate. This is due to an
increase in the RP normal mode population, which leads to an increase of the
last term in Eq. 4.

Figure 4: Decomposition of instanton geometry into the free RP normal mode
basis at a representative temperature of $T=0.70T^{\circ}_{c}=80$ K for (a)
symmetric DW with $\tilde{\eta}_{0}/m\omega^{\ddagger}=0.00$, (b) symmetric DW
with $\tilde{\eta}_{0}/m\omega^{\ddagger}=0.50$, (c) asymmetric DW
($q_{0}=0.08$Å) with $\tilde{\eta}_{0}/m\omega^{\ddagger}=0.00$, and (d)
asymmetric DW ($q_{0}=0.08$Å) with $\tilde{\eta}_{0}/m\omega^{\ddagger}=0.50$.
Five different barrier heights were considered: 125 meV (purple), 258 meV
(green), 500 meV (red), 750 meV (light blue), and 1000 meV (gray).
Coefficients are ordered and grouped by their corresponding RP normal mode
(NM) index (l).

Grote-Hynes (GH) theory Grote and Hynes (1980, 1981); Pollak (1986) defines a
relationship between reaction rates obtained with a finite friction strength
and those obtained with vanishing friction strength in the classical limit. In
part I of this article, we showed an extension of GH theory to the deep
tunneling regime for the case of position-independent friction. Briefly, we
proposed that tunneling rates with finite friction strength – RPI-EF rates –
at a given temperature, can be related to tunneling rates without friction –
RPI rates – performed at a scaled temperature. The scaling relation for the
temperatures is given by

$\begin{split}T_{b}/T_{c}^{\circ}=\sqrt{\bigg{(}\frac{\tilde{\eta}(\omega_{l}^{b})}{2m{\omega^{\ddagger}}}\bigg{)}^{2}\frac{1}{l^{2}}+\left(\frac{T_{a}}{T_{c}^{\circ}}\right)^{2}}-\frac{\tilde{\eta}(\omega_{l}^{b})}{2m{\omega^{\ddagger}}}\frac{1}{l}\end{split}$
(13)

where $T_{c}^{\circ}$ is the cross-over temperature without friction, $T_{b}$
the target temperature at which the RPI-EF result is desired, $T_{a}$ the
temperature at which the RPI calculation must be performed, and
$\omega_{l}^{b}$ the $l$ free RP normal mode frequency at $T_{b}$. Since
$\omega^{b}_{l}$ depends on $T_{b}$, this equation has to be solved self-
consistently. Even though Eq. 13 must be fulfilled for all $l$, from Fig. 4
one can expect that considering only $\omega_{l=\pm 1}$ should be a good
assumption.

Figure 5: Error on the RPI rate values obtained by the scaling relation shown
in Eq. 13 ($k_{\text{RPI}}^{\text{GH}}$), compared to RPI-EF rates
($k_{\text{RPIEF}}$). The error is reported as the logarithm of the ratio
between these rates in the symmetric DW model. Temperatures of
0.70$T^{\circ}_{c}$ (squares) and 0.55$T^{\circ}_{c}$ (circles), and reaction
barriers of 258 meV (red), 500 meV (blue), and 1000 meV (black) are shown. An
analogous plot for an asymmetric barrier is presented in the SI.

In Fig. 5, the error obtained by computing the rate, using only $\omega_{l=\pm
1}$ in Eq. 13, for different temperatures and coupling strengths in the
symmetric DW potential is presented. The estimated RPI rates at the scaled
temperatures are within one order of magnitude from the full RPI-EF rates for
all friction strengths, but they are in better agreement for
$\tilde{\eta}_{0}/m\omega^{{\ddagger}}<0.5$. Similar accuracy is observed for
an asymmetric DW model with this approximation, even though the $l=0$ mode is
appreciably activated (see section IV in the SI).

###  IV.2 Position-Dependent Friction

Figure 6: Instanton pathways obtained at 40 K using the DDW model with
$V_{0}=258$ meV, $q_{0}=0$ Å, $C=97.1$ meV/Å2 ( 0.001 a.u.), $\epsilon_{1}$ =
0, $\epsilon_{2}=-0.8$, and $\Delta=1.0$ for
$\tilde{\eta}_{0}/m\omega^{\ddagger}$= 0, 0.10, 0.25, and 0.50 represented by
full black, red, orange, and blue lines, respectively. The pathways are shown
on top of a) a heat map representing the underlying potential energy surface
b) a heat map that helps visualize the position-dependence of the friction
tensor. The map is computed as the sum of the diagonal elements of
$\tilde{\eta}$ for $\tilde{\eta}_{0}/m\omega^{\ddagger}=0.50$ (see Eq. 8 and
10 in the main text).

We now consider the case where the coupling between the system and the bath
depends on the position of the system coordinate, i.e. a position-dependent
friction. Fig. 6 shows instanton pathways obtained for different system-bath
coupling strengths on a DDW potential, as described by Eq. 9, at a temperature
considerably lower than $T_{c}^{\circ}$ in this model. The parameters of the
model were chosen such that the PES presents two global minima, two local
minima, four first order saddle points that connect each global minimum with
the closest local minimum, and one second order saddle point at
$(q_{1},q_{2})=(0,0)$, as shown in Fig. 6a. The optimal tunneling pathway in
the absence of dissipation is represented by the black curve in Fig. 6 and it
is a linear trajectory that connects the two global minima by crossing the
second order saddle point. For $\tilde{\eta}>0$, the anisotropy of the
friction (see Fig. 6b) results in a modification of the instanton pathway,
which bends towards regions of lower friction values. The magnitude of the
bending of this path increases as the strength of the friction becomes larger.
This shows that the optimal dissipative tunneling pathway is a compromise
between the path with the shortest length, the path with the lowest potential
energy, and the path with lowest friction. Indeed, for
$\tilde{\eta}_{0}/m\omega^{\ddagger}>0.5$, the dissipation is so strong close
to the second-order saddle point that no instanton pathway that connects
directly the two global minima can be found.

##  V Hopping of Hydrogen and Deuterium in Bulk Metals

###  V.1 Minimum energy paths, barrier heights and friction strengths

Having characterized the performance of RPI-EF in model potentials, we now
address the interplay of NAEs and NQEs on the hopping reaction of H within
bulk metals, which we calculate from first-principles electronic-structure
simulations.

We first analyze the MEP of the reactions in Pd, Pt, Cu, Ag and Al. We note
that we focus on the hopping reaction between neighboring
octahedral$\rightarrow$tetrahedral interstitial sites. Even though in perfect
solids these reactions determine the diffusion rate, in real materials other
mechanisms might become the rate determining step of the diffusion process Di
Stefano, Mrovec, and Elsässer (2015); Pedersen and Jónsson (2009). In Table 1,
we report the reaction energy, reaction barrier, and electronic-friction
values along the MEP for the reactions considered in this work. As shown in
column 3, the energy barriers are in the 100-300 meV range, in accordance with
previous studies Ferrin  _et al._ (2012).

System |  $E_{\text{T-O}}$ (meV) |  $E_{\text{TS-O}}$ (meV) |  $\tilde{\eta}$ (ps-1) |  $T_{c}^{\circ}$ (K) |  $\omega^{\ddagger}$ (cm-1)
---|---|---|---|---|---
H$@$Pd | 43 | 148 | 0.7 - 2.7 | 115 | 501
H$@$Pt | -35 | 44 | 0.8 - 2.8 | 96 | 420
H$@$Cu | 188 | 300 | 0.7 - 1.1 | 140 | 612
H$@$Ag | 52 | 160 | 0.7 - 1.0 | 116 | 504
H$@$Al | -71 | 88 | 1.8 - 3.1 | 84 | 365
Table 1:  Reaction energy, $E_{\text{T-O}}=E_{\text{T}}-E_{\text{O}}$, and
energy barrier heights, $E_{\text{TS-O}}=E_{\text{TS}}-E_{\text{O}}$, for the
different fcc metals considered in this work. $E_{\text{T}}$, $E_{\text{O}}$,
and $E_{\text{TS}}$ refer to the potential energy corresponding to structures
where the H atom is located at the tetrahedal (T), octahedral (O) and
transition state (TS) sites, respectively. Values are reported without ZPE
corrections. Minimum and maximum values adopted by the electronic friction,
$\tilde{\eta}$, along the MEP are presented in column four. Values are
evaluated at 54 meV (which corresponds to the first non-zero ring-polymer
normal mode frequency at 100 K) and considering a projection of $\tilde{\eta}$
on the reaction coordinate. Columns 5 and 6 show the crossover temperature,
$T_{c}^{\circ}$, and imaginary frequency at the TS, $\omega^{\ddagger}$,
respectively.

We continue by analyzing the electronic friction tensor
$\tilde{\eta}_{il}(\bm{q},\lambda)$ on the hydrogen atom along the MEPs. In
Fig. 7a, we present the electronic friction values evaluated at the first non-
zero ring-polymer normal mode frequency at 100 K, projected onto the direction
parallel to the reaction coordinate for the case of Pd. The friction values
vary up to almost an order of magnitude, indicating the necessity of having a
rate theory that takes into account such position dependence. The electronic
friction along the MEP for the other metals shows a strong position dependence
as well (see section III B of the SI). As shown in column 4 of Table 1, the
values vary from 0.7 ps-1 to 3.1 ps-1. The magnitude of the friction
coefficients is large enough to impact vibrational lifetimes Maurer  _et al._
(2016), adsorption mechanisms Bünermann  _et al._ (2015), and scattering
experiments Kandratsenka  _et al._ (2018); Bünermann  _et al._ (2015).
However, the dimensionless coefficient $\tilde{\eta}/m\omega^{\ddagger}$
yields at most a value of 0.05 (for Al), which, given the relatively low
barrier heights, would result in a reduction of the tunneling rates by less
than a factor of 5, according to our study on model potentials presented in
the previous sections (see Figs. 2 and 3).

A closer look at the expression used to compute the electronic friction
tensor, Eq. 11, allows us to rationalize the reasons behind such small
coefficients Forsblom and Persson (2007). Large friction values will arise in
the case of a high DOS close to the Fermi level and due to the presence of
hydrogen states close to the Fermi level. The former contributes to Eq. 11 via
the Fermi-Dirac factors, the latter contributes via the strength of the
nonadiabatic coupling. While the former depends mainly on the metal at hand,
the latter is affected by both the impurity and the metalLecroart  _et al._
(2021). We analyzed the atomic projected DOS (see section III C in the SI) and
confirmed that Pd and Pt are, at the same time, the systems that present the
largest electronic friction values and the highest DOS at the Fermi level
among the transition metals. Surprisingly, Al presents slightly larger
friction values, without a high DOS at the Fermi energy. This might suggest
that the nonadiabatic couplings are comparatively large in this case. However,
in all cases, the hydrogen atom neither creates new states nor affects the DOS
appreciably in the vicinity of the Fermi level, which ultimately leads to
rather small friction coefficients along the MEP for these systems.

Figure 7:  a) Minimum energy pathway (MEP) and friction along the reaction
coordinate for the H hopping reaction in Pd. The energy is set to zero at the
reactant geometry. b) Frequency dependence of the friction tensor projected on
the reaction coordinate at the reactant (black), transition state (red), and
product (blue) states. To ease visual comparison, all curves in panel b have
been scaled to adopt the value of 1 at the highest friction value. The first
three non-zero ring-polymer normal modes frequencies at 100K are depicted as
vertical dashed gray lines.

Up to this point, we have considered values of the friction tensor at a single
frequency. However, the calculation of the tunneling rates require the
evaluation of the friction tensor at all the RP normal modes frequencies (Eq.
3), and more importantly in the derivation of RPI-EF with a position-dependent
friction tensor, we have assumed that a separable coupling ansatz is valid.
This ansatz is equivalent to assuming that the frequency dependence of the
friction tensor remains the same at all (relevant) positions. As an
illustrative example, we present in Fig. 7b the frequency dependence of the
friction tensor at the stationary points of the MEP (reactant, transition and
product state) for the hydrogen hopping reaction in Pd. The frequency
dependence shows a non-monotonic profile with a maximum around 0.6 eV and
remains fairly similar along the MEP suggesting that the non-adiabatic
couplings in bulk metals are, to a great extent, well described by a
“separable coupling". This observation is equally valid for the other metals
(see SI section III B). A different type of coupling might be observed in
scattering reactions, where atoms or molecules transition from vacuum to
electron-rich environmentsBünermann  _et al._ (2015); Maurer  _et al._ (2017).
As discussed in section IV.1, only the first few ring-polymer normal modes are
appreciably activated. For this reason, only a relatively small fraction of
spectral density contributes to the rates (see Fig. 7b).

###  V.2 Tunneling rates: The case of Pd

We performed full-dimensional instanton calculations on Pd in order to gauge
the predictive power of our studies on low-dimensional models and static
estimators presented in the previous sections. We selected Pd because it
presents high values of friction along the MEP and the diffusion of H in Pd
has been well studied theoretically and experimentally before. In order to
reduce the computational cost, we performed instanton calculations on a fcc
cubic cell containing 4 Pd atoms and 1 H or D atom. These calculations were
performed using a 12 $\times$ 12 $\times$ 12 k-point sampling. The relatively
small size of the unit cell induces an effective increase of the barrier when
compared to larger unit cells and the new $T_{c}^{\circ}$ increases to 136 K
(see section I in the SI). Since larger barriers magnify the effect of the
friction on the rates, these calculations can be considered an upper-limit
estimation of the impact of the friction on these rates.

In Fig. 8, the reaction rates for the hopping reaction of H and D in Pd from
the octahedral to the tetrahedral site using TST, RPI theory and RPI-EF theory
are presented. At temperatures below $T_{c}^{\circ}$, the tunneling effects,
evidenced by the difference between the RPI and TST rates, become increasingly
important, enhancing the rate by several orders of magnitude. The comparison
of the TST predictions for H and D indicates an inverse kinetic isotope effect
(KIE). The inverse KIE can be traced back to the softening of the normal modes
orthogonal to the reaction pathway at the transition-state geometry Meisner
and Kästner (2016). Since this effect is mainly due to ZPE, it is also present
at temperatures above $T_{c}$ and has been reported experimentally Völkl and
Alefeld (1978). Below $T_{c}$, the emergence of tunneling creates a
competition between ZPE and tunneling effects as already reported by Shiga et
al. Kimizuka, Ogata, and Shiga (2019). Moreover, the similarity of the
reaction rates for both isotopologues around 80 K is in agreement with the
results of Ref. Kimizuka, Ogata, and Shiga (2019) and with experiments.
However, the absolute values reported here are considerably smaller as a
consequence of the rather small unit-cell employed in the calculations.

Figure 8: Reaction rates for octahedral$\to$tetrahedral-site hopping reaction
of H (black) and D (orange) in Pd, calculated by TST (dashed lines), and RPI
rate theory (solid lines). RPI-EF rates with position-independent friction,
fixed at a value of 4 ps-1, are shown by black and orange circles for H and D,
respectively, while the RPI-EF rates with position-dependent friction is
depicted as a red cross. RPI rate calculations with the Pd atoms fixed at
their reactant, transition and product states are presented by diamond, square
and circle gray symbols, respectively. $T_{c}^{\circ}$ for H and D are marked
by vertical black and orange dotted lines, respectively.

We now inspect the impact of the lattice on the reaction rates. For this
purpose, we performed RPI calculations where the Pd atoms were fixed at the
reactant, product and transition state positions at 115 K. The rate obtained
when the Pd atoms are fixed at their reactant positions is 7 and 4 orders of
magnitude slower than the RPI and TST estimates, respectively. This confirms
the significant contribution of the lattice fluctuations to the reactive
process Kimizuka, Ogata, and Shiga (2018) and highlights the multidimensional
nature of the reaction. In contrast, the rate estimates obtained from the
calculations with the Pd atoms fixed at their transition-state or product
positions are within an order of magnitude of the RPI estimates. The former
result is expected since we are analyzing the rates at only 20 K below
$T_{c}^{\circ}$, such that the instanton pathway lies very close to the TS
geometry. The latter result, instead, shows that the lattice relaxation
between the transition state and the product is of comparatively lower
relevance.

Finally, we consider the effect of the electronic friction on the rate. At 100
K, we performed a calculation with on-the-fly estimations of the friction
tensor along the instanton pathway (also included in the instanton
optimization). In order to gauge the importance of the spatial dependence of
the friction, we also performed, at the same temperature, an instanton
calculation using a constant value of 4 ps-1 for the friction (slightly higher
than the maximum value reported in Table 1). The results were numerically
indistinguishable within the accuracy of our calculations, which may seem an
unexpected result, at first. We proceeded to perform calculations with a
constant and spatially-independent friction value (which is less
computationally demanding) at several other temperatures. In Figure 8, we show
these results and, as predicted by our earlier assessment, the friction
coefficients are not large enough to produce a significant effect on the
rates. This explains why the spatial dependence also does not appreciably
change the rate constants.

##  VI  Conclusions

We have benchmarked the RPI-EF method and showed its performance in model
potentials and first-principles calculations. By performing numerically exact
simulations in 1D model systems including a spatially-independent friction, we
showed that RPI-EF yields accurate rates for all but very small friction
coefficients at a much reduced computational cost. A systematic analysis of 1D
and 2D double-well potentials allowed us to determine the magnitude of the
decrease of the tunneling rates caused by friction, as a function of the
coupling strength and barrier height. We found that the suppression of
tunneling is promoted by high coupling strengths and high energy barriers. We
were also able to demonstrate that for a spatially-dependent friction tensor,
the instanton pathway can be considerably deformed towards low-friction
regions, when compared to the “non-dissipative” path (without friction). In
comparison to previous similar approachesCaldeira and Leggett (1983); Weiss
(2008), the RPI-EF method is advantageous because it is highly-efficient, more
intuitive and mathematically simpler.

In the context of reactions involving atoms and molecules in metallic
environments, RPI-EF allows the inclusion of NQEs and NAEs as described by an
effective electronic friction. While here we used an electronic friction
formulation that disregards electronic correlation Head-Gordon and Tully
(1995); Dou, Miao, and Subotnik (2017); Maurer  _et al._ (2016), the RPI-EF
approach is rather general and can be combined with other flavours of
electronic friction that go beyond the independent quasi-particle picture. As
a consequence of the relatively low computational cost of both RPI and the
employed ab initio electronic friction formalism, RPI-EF allows the study of
high-dimensional systems with on-the-fly ab initio evaluation of the forces
and the electronic friction tensor.

In this work, we presented calculations of hydrogen and deuterium hopping
between nearest interstitial sites of selected fcc metals, employing density-
functional theory calculations. By evaluating the impact that different
factors have on the reactions rates of this reaction in bulk Pd, we
established that nuclear tunneling and lattice relaxation play a larger role
in determining the magnitude of the rate than electronic friction. The latter
turned out to have a negligible impact on the reaction rates of these systems.
This negative result, however, answers an important theoretical question
regarding the interplay between NAEs (modeled by electronic friction) and
NQEsDou and Subotnik (2018), and also validates the rate constants currently
used in this context for multiscale modelingHüter  _et al._ (2018).
Nonetheless, we anticipate that for impurities or adsorbates that present
electronic levels in the vicinity of the Fermi energy of the metal, NAEs might
play a more prominent role in connection with tunneling. We also expect to
observe a larger effect for lighter particles such as Muons Teichler (1977),
or for surface reactions with higher energy barriers del V. Gómez  _et al._
(2017).

###### Acknowledgements.

Y.L., E.S.P. and M.R. acknowledge financing from the Max Planck Society and
computer time from the Max Planck Computing and Data Facility (MPCDF). Y.L and
M.R. thank Jeremy Richardson, Aaron Kelly, and Stuart Althorpe for a critical
reading of the manuscript. C.L.B. acknowledges financial support through an
EPSRC-funded PhD studentship. R.M. acknowledges Unimi for granting computer
time at the CINECA HPC center. R.J.M. acknowledges financial support through a
Leverhulme Trust Research Project Grant (RPG-2019-078) and the UKRI Future
Leaders Fellowship programme (MR/S016023/1).

## References

  * Mohtadi and Orimo (2016) R. Mohtadi and S.-i. Orimo, Nat. Rev. Mater. 2, 16091 (2016).
  * Robertson  _et al._ (2015) I. M. Robertson, P. Sofronis, A. Nagao, M. L. Martin, S. Wang, D. W. Gross, and K. E. Nygren, Metall. Mater. Trans. B 46, 1085 (2015).
  * Traidia, Chatzidouros, and Jouiad (2018) A. Traidia, E. Chatzidouros, and M. Jouiad, Corros. Rev. 36, 323 (2018).
  * Schlapbach and Züttel (2001) L. Schlapbach and A. Züttel, Nature 414, 353 (2001).
  * Völkl and Alefeld (1978) J. Völkl and G. Alefeld, “Diffusion of hydrogen in metals,” in  _Hydrogen in Metals I: Basic Properties_, edited by G. Alefeld and J. Völkl (Springer Berlin Heidelberg, Berlin, Heidelberg, 1978) pp. 321–348.
  * Jothi, Croft, and Brown (2015) S. Jothi, T. Croft, and S. Brown, Int. J. Hydrog. Energy 40, 2882 (2015).
  * Fuk (2005) “Diffusion,” in  _The Metal-Hydrogen System: Basic Bulk Properties_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2005) pp. 303–400.
  * Nienhaus (2002) H. Nienhaus, Surf. Sci. Rep. 45, 1 (2002).
  * Wodtke (2016) A. M. Wodtke, Chem. Soc. Rev. 45, 3641 (2016).
  * Kimizuka, Ogata, and Shiga (2018) H. Kimizuka, S. Ogata, and M. Shiga, Phys. Rev. B 97, 014102 (2018).
  * Kimizuka and Shiga (2021) H. Kimizuka and M. Shiga, Phys. Rev. Mater. 5, 065406 (2021).
  * Beck  _et al._ (2000a) M. Beck, A. Jäckle, G. Worth, and H.-D. Meyer, Phys. Rep. 324, 1 (2000a).
  * Tanimura and Kubo (1989) Y. Tanimura and R. Kubo, J. Phys. Soc. Japan 58, 101 (1989).
  * Ilk and Makri (1994) G. Ilk and N. Makri, J. Chem. Phys. 101, 6708 (1994).
  * Hellsing and Persson (1984) B. Hellsing and M. Persson, Phys. Scr. 29, 360 (1984).
  * Head-Gordon and Tully (1995) M. Head-Gordon and J. C. Tully, J. Chem. Phys. 103, 10137 (1995).
  * Richardson (2018) J. O. Richardson, Int. Rev. Phys. Chem. 37, 171 (2018).
  * Richardson and Althorpe (2009a) J. O. Richardson and S. C. Althorpe, J. Chem. Phys. 131, 214106 (2009a).
  * Arnaldsson (2007) A. Arnaldsson,  _Calculation of quantum mechanical rate constants directly from ab initio atomic forces_ , Ph.D. thesis, University of Washington (2007).
  * Eyring (1935) H. Eyring, J. Chem. Phys. 3, 107 (1935).
  * Kapil  _et al._ (2019) V. Kapil, M. Rossi, O. Marsalek, R. Petraglia, Y. Litman, T. Spura, B. Cheng, A. Cuzzocrea, R. H. Meißner, D. M. Wilkins, B. A. Helfrecht, P. Juda, S. P. Bienvenue, W. Fang, J. Kessler, I. Poltavsky, S. Vandenbrande, J. Wieme, C. Corminboeuf, T. D. Kühne, D. E. Manolopoulos, T. E. Markland, J. O. Richardson, A. Tkatchenko, G. A. Tribello, V. V. Speybroeck, and M. Ceriotti, Comp. Phys. Comm. 236, 214 (2019).
  * Blum  _et al._ (2009) V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, Comp. Phys. Comm. 180, 2175 (2009).
  * E, Ren, and Vanden-Eijnden (2007) W. E, W. Ren, and E. Vanden-Eijnden, J. Chem. Phys. 126, 164103 (2007).
  * Henkelman, Uberuaga, and Jónsson (2000) G. Henkelman, B. P. Uberuaga, and H. Jónsson, J. Chem. Phys. 113, 9901 (2000).
  * Nichols  _et al._ (1990) J. Nichols, H. Taylor, P. Schmidt, and J. Simons, J. Chem. Phys. 92, 340 (1990).
  * Beyer  _et al._ (2016) A. N. Beyer, J. O. Richardson, P. J. Knowles, J. Rommel, and S. C. Althorpe, J. Phys. Chem. Lett. 7, 4374 (2016).
  * Litman (2020) Y. Litman,  _Tunneling and Zero-Point Energy Effects in Multidimensional Hydrogen Transfer Reactions: From Gas Phase to Adsorption on Metal Surfaces_ , Ph.D. thesis, Freie Universität Berlin (2020).
  * Wang and Thoss (2003) H. Wang and M. Thoss, J. Chem. Phys. 119, 1289 (2003).
  * Manthe (2009) U. Manthe, J. Chem. Phys. 130, 054109 (2009).
  * Vendrell and Meyer (2011) O. Vendrell and H.-D. Meyer, J. Chem. Phys. 134, 044135 (2011).
  * Meyer, Manthe, and Cederbaum (1990) H.-D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. Phys. Lett. 165, 73 (1990).
  * Beck  _et al._ (2000b) M. H. Beck, A. Jäckle, G. A. Worth, and H.-D. Meyer, Phys. Rep. 324, 1 (2000b).
  * Meyer, Gatti, and Worth (2009) H.-D. Meyer, F. Gatti, and G. A. Worth, eds.,  _Multidimensional Quantum Dynamics: MCTDH Theory and Applications_ (Wiley-VCH, Weinheim, 2009).
  * (34) G. A. Worth, M. H. Beck, A. Jäckle, O. Vendrell, and H.-D. Meyer, The MCTDH Package, Version 8.5 (2013), http://mctdh.uni-hd.de/.
  * Cirac  _et al._ (2021) J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Rev. Mod. Phys. 93, 045003 (2021).
  * Wang, Skinner, and Thoss (2006a) H. Wang, D. E. Skinner, and M. Thoss, J. Chem. Phys. 125, 174502 (2006a).
  * Craig, Thoss, and Wang (2007a) I. R. Craig, M. Thoss, and H. Wang, J. Chem. Phys. 127, 144503 (2007a).
  * Miller, Schwartz, and Tromp (1983) W. H. Miller, S. D. Schwartz, and J. W. Tromp, The Journal of Chemical Physics 79, 4889 (1983).
  * Miller (1998) W. H. Miller, Journal of Physical Chemistry A 102, 793 (1998).
  * Wang, Skinner, and Thoss (2006b) H. Wang, D. E. Skinner, and M. Thoss, The Journal of Chemical Physics 125, 174502 (2006b).
  * Craig, Thoss, and Wang (2007b) I. R. Craig, M. Thoss, and H. Wang, Journal of Chemical Physics 127 (2007b), 10.1063/1.2772265.
  * Topaler and Makri (1994) M. Topaler and N. Makri, J. Chem. Phys. 101, 7500 (1994).
  * Craig and Manolopoulos (2005) I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 122, 0 (2005).
  * Straus, Gomez Llorente, and Voth (1993) J. B. Straus, J. M. Gomez Llorente, and G. A. Voth, J. Chem. Phys. 98, 4082 (1993).
  * Perdew, Burke, and Ernzerhof (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Haas, Tran, and Blaha (2009) P. Haas, F. Tran, and P. Blaha, Phys. Rev. B 79, 085104 (2009).
  * Maurer  _et al._ (2016) R. J. Maurer, M. Askerka, V. S. Batista, and J. C. Tully, Phys. Rev. B 94, 115432 (2016).
  * Dou, Miao, and Subotnik (2017) W. Dou, G. Miao, and J. E. Subotnik, Phys. Rev. Lett. 119, 046001 (2017).
  * Luntz and Persson (2005) A. C. Luntz and M. Persson, J. Chem. Phys. 123, 074704 (2005).
  * Monturet and Saalfrank (2010) S. Monturet and P. Saalfrank, Phys. Rev. B 82, 075404 (2010).
  * Richardson and Althorpe (2009b) J. O. Richardson and S. C. Althorpe, J. Chem. Phys. 131, 214106 (2009b).
  * Grabert, Weiss, and Hanggi (1984) H. Grabert, U. Weiss, and P. Hanggi, Phys. Rev. Lett. 52, 2193 (1984).
  * Grabert, Olschowski, and Weiss (1987) H. Grabert, P. Olschowski, and U. Weiss, Phys. Rev. B 36, 1931 (1987).
  * Grote and Hynes (1980) R. F. Grote and J. T. Hynes, J. Chem. Phys. 73, 2715 (1980).
  * Grote and Hynes (1981) R. F. Grote and J. T. Hynes, J. Chem. Phys. 74, 4465 (1981).
  * Pollak (1986) E. Pollak, J. Chem. Phys. 85, 865 (1986).
  * Di Stefano, Mrovec, and Elsässer (2015) D. Di Stefano, M. Mrovec, and C. Elsässer, Acta Mater. 98, 306 (2015).
  * Pedersen and Jónsson (2009) A. Pedersen and H. Jónsson, Acta Mater. 57, 4036 (2009).
  * Ferrin  _et al._ (2012) P. Ferrin, S. Kandoi, A. U. Nilekar, and M. Mavrikakis, Surf. Sci. 606, 679 (2012).
  * Bünermann  _et al._ (2015) O. Bünermann, H. Jiang, Y. Dorenkamp, A. Kandratsenka, S. M. Janke, D. J. Auerbach, and A. M. Wodtke, Science 350, 1346 (2015).
  * Kandratsenka  _et al._ (2018) A. Kandratsenka, H. Jiang, Y. Dorenkamp, S. M. Janke, M. Kammler, A. M. Wodtke, and O. Bünermann, Proc. Natl. Acad. Sci. USA 115, 680 (2018).
  * Forsblom and Persson (2007) M. Forsblom and M. Persson, J. Chem. Phys. 127, 154303 (2007).
  * Lecroart  _et al._ (2021) L. Lecroart, N. Hertl, Y. Dorenkamp, H. Jiang, T. N. Kitsopoulos, A. Kandratsenka, O. Bünermann, and A. M. Wodtke, J. Chem. Phys. 155, 034702 (2021).
  * Maurer  _et al._ (2017) R. J. Maurer, B. Jiang, H. Guo, and J. C. Tully, Phys. Rev. Lett. 118, 256001 (2017).
  * Meisner and Kästner (2016) J. Meisner and J. Kästner, Angew. Chem. Int. Ed. 55, 5400 (2016).
  * Kimizuka, Ogata, and Shiga (2019) H. Kimizuka, S. Ogata, and M. Shiga, Phys. Rev. B 100, 024104 (2019).
  * Caldeira and Leggett (1983) A. Caldeira and A. Leggett, Ann. Phys. 149, 374 (1983).
  * Weiss (2008) U. Weiss,  _Quantum Dissipative Systems_, 3rd ed. (WORLD SCIENTIFIC, 2008).
  * Dou and Subotnik (2018) W. Dou and J. E. Subotnik, J. Chem. Phys. 148, 230901 (2018).
  * Hüter  _et al._ (2018) C. Hüter, P. Shanthraj, E. McEniry, R. Spatschek, T. Hickel, A. Tehranchi, X. Guo, and F. Roters, Metals 8 (2018).
  * Teichler (1977) H. Teichler, Phys. Lett. A 64, 78 (1977).
  * del V. Gómez  _et al._ (2017) E. del V. Gómez, S. Amaya-Roncancio, L. B. Avalle, D. H. Linares, and M. C. Gimenez, Appl. Surf. Sci. 420, 1 (2017).