Reactions involving adsorbates on metallic surfaces and impurities in bulk metals are ubiquitous in a wide range of technological applications. The theoretical modeling of such reactions presents a formidable challenge for theory because nuclear quantum effects (NQEs) can play a prominent role and the coupling of the atomic motion with the electrons in the metal gives rise to important non-adiabatic effects (NAEs) that alter atomic dynamics. In this work, we derive a theoretical framework that captures both NQEs and NAEs and, due to its high efficiency, can be applied to first-principles calculations of reaction rates in high-dimensional realistic systems. More specifically, we develop a method that we coin ring polymer instanton with explicit friction (RPI-EF), starting from the ring polymer instanton formalism applied to a system–bath model. We derive general equations that incorporate the spatial and frequency dependence of the friction tensor and then combine this method with the ab initio electronic friction formalism for the calculation of thermal reaction rates. We show that the connection between RPI-EF and the form of the electronic friction tensor presented in this work does not require any further approximations, and it is expected to be valid as long as the approximations of both underlying theories remain valid.

Metallic systems lack an energy gap between unoccupied and occupied electronic states. Thus, low energy excitation and deexcitation of electron–hole pairs can easily exchange energy with nuclear vibrations, representing a violation of the Born–Oppenheimer principle, where electrons are assumed to adjust adiabatically to the position of the nuclei. This type of non-adiabatic effect (NAE) has been verified experimentally numerous times in the past1–3 and found to be particularly important for hot-electron-induced reactions,4 surface scattering,5,6 and vibrational relaxation lifetimes.7–9 Many theoretical approaches have been developed to account for NAEs in these contexts.10–13 Among them, the method coined “molecular dynamics with electronic friction”13 (MDEF) has the advantage of being a method that can currently be coupled to ab initio electronic structure theory without resorting to model parameterization.14,15

MDEF describes the motion of classical nuclei through a Langevin equation, where the friction forces and the corresponding random noise embody the effects of the electronic excitations, therefore, approximately including NAEs to an otherwise classical nuclear dynamics. The electronic friction tensor lies at the core of the definition of the frictional force in MDEF and can be understood as a first-order correction to the Born–Oppenheimer approximation in the presence of a manifold of fast relaxing electronic states.16 As a consequence, MDEF is expected to break down when strong non-adiabatic effects are present1,17 and in cases where charge transfer mechanisms are dominant.11 Despite these shortcomings, MDEF has proven to be a useful approach for several realistic systems and conditions, wherein high-level simulations can explain well-controlled experiments.6,9,18–21

The classical-nuclei approximation in the context of MDEF is appropriate for many situations. However, when studying the motion of light atoms, such as hydrogen, deuterium, and lithium, the quantum nature of the nuclei can lead to significant quantum effects, including tunneling, isotope effects, and non-Arrhenius temperature dependence of reaction rates.22–25 Indeed, Fang et al.26 showed that depending on the shape of the barrier for hydrogen diffusion on metallic surfaces, a coexistence of deep tunneling through the barrier and classical hopping above the barrier can take place even at elevated temperatures, while Kimizuka et al. exposed a strong interplay between the relevance of the quantum fluctuations and the lattice strain of interstitial H diffusion in metals.27 

In this work, we propose a new approach based on the ring polymer instanton (semi-classical) rate theory,28 which includes NAEs through the electronic friction formalism initially proposed by Hellsing and Persson29 and Head-Gordon and Tully.13 Because the instanton theory relies on an imaginary-time propagation, we arrive at a formulation in which the friction modifies the ring polymer potential but no fluctuation-force term appears, in contrast to real-time dynamics methods. As we will discuss later, the description of the electrons as a harmonic bath of non-interacting particles is the key consideration that makes the connection between ring polymer instanton with explicit friction (RPI-EF) and electronic friction possible. We show how both theories can be combined naturally and include the full tensorial nature, frequency and position dependence of the electronic friction.

This paper is organized as follows: In Sec. II, we present a brief review of the ring polymer instanton (RPI) theory. In Sec. III, we introduce and discuss our proposed theory coined ring polymer instanton with explicit friction (RPI-EF). In Sec. IV, we present and discuss the connection of RPI-EF and an ab initio electronic friction, and finally conclude in Sec. V. In Paper II, we present benchmarks of our method for model systems and an application to hydrogen hopping in Pd, employing Kohn–Sham density-functional theory.

The RPI approximation30,31 is a semi-classical method based on the path integral formulation of quantum mechanics and allows the evaluation of reaction rates in the deep tunneling regime.28,32 The theory assumes that only a couple of well-defined reactant and product states are sufficient to describe the reactive process under consideration. This condition is met for gas-phase reactions and atomic diffusion on (or in) solids,26 but it is rarely satisfied for liquid environments.33 The RPI approximation replaces the quantum mechanical time propagator in the flux-side time correlation function with its semi-classical counterpart and allows one to express the reaction rate in terms of the dominant stationary trajectory in imaginary time that connects reactants and products, i.e., the instanton trajectory. This special trajectory can be interpreted as the main tunneling pathway at a given temperature and offers an intuitive picture of the process under study. The RPI rate theory is based on the “Im F” premise, where the rate is related to the imaginary part of the free energy.34 It has been shown by Althorpe35 to be equivalent to the original derivation by Miller32 expressed in terms of the stability parameters of the instanton orbit, and, more recently, it has been derived by Richardson from quantum scattering theory.36 

To find the instanton pathway, it is numerically convenient to discretize the closed trajectories and represent them by ring polymers.28 The discretized Euclidean action for a given trajectory in imaginary time, SP, is related to the potential energy of the ring polymer, UP, by

SP/=βPUP
(1)

with

UP(q)=k=1Pi=13NmiωP22(qi(k)qi(k+1))2+k=1PV(q1(k),,q3N(k)).
(2)

Here, qi(k) is the position of the ith degree of freedom of the kth replica, mi is the mass of the ith degree of freedom, N is the number of atoms, P is the number of replicas, q is an abbreviated notation to represent all the degrees of freedom, and ωP=(βP)1 with βP=1kBPT, where kB is the Boltzmann constant and T is the temperature. As a result of Eq. (1) and the fact that the instanton pathway is a stationary trajectory in the imaginary time, the instanton geometry represents a stationary point on UP. Moreover, it constitutes a first-order saddle point and can easily be found by standard saddle-point search algorithms.37 

The tunneling rate can be expressed as28 

kinst(β)=2βImF=limP1ZPr(β)2βImeSP(q)/dq,
(3)

where F is the system’s complex free energy and ZPr is the reactant canonical partition function. The evaluation of the integral is performed by a steepest-descent integration around the instanton geometry q̄ for all modes with positive eigenvalues, while the mode with negative eigenvalue and the mode with zero eigenvalue, which corresponds cyclic permutation of beads, require special care. As a result, the instanton rate reads

kinst(β)=limP1ZPr(β)1βPBP(q̄)2πβP2ZPv(β)eSP(q̄)/,
(4)

with

BP(q̄)=i=13Nk=1Pmi(q̄i(k+1)q̄i(k))2
(5)

and

ZPv=k1βP|λk|.
(6)

In the expression above, λk represent the P × 3N eigenvalues of the mass scaled ring polymer Hessian defined by the second derivatives of UP with respect to the replica positions, and the prime indicates that the product is taken over all modes except the ones with zero eigenvalue. The contribution of translational and rotational degrees of freedom has been discarded in Eq. (4) since they are not relevant to the present work. The accuracy of the tunneling rates defined by this theory is limited mainly due to two reasons: (i) the fluctuations orthogonal to the reactive direction are considered to be harmonic and, (ii) due to lack of real-time information, the recrossing effects are completely neglected. Despite these shortcomings, the RPI approximation constitutes a valuable and practical method due to its favorable trade-off between accuracy and computational cost, and the method has been successfully applied to systems containing up to several hundreds of degrees of freedom, by using an ab initio description of their electronic structure.38–40 

Instanton trajectories only exist below a critical temperature known as cross over temperature, Tc°, which, in most cases, can be estimated by a parabolic barrier approximation as

kBTc°=ω2π,
(7)

where ω represents the imaginary frequency at the barrier top between reactants and products. At temperatures below Tc°, the reactive process is dominated by tunneling, while above Tc°, the classical “over-the-barrier hopping” mechanism represents the major contribution, with nuclear tunneling playing a minor role.41 Due to the lack of real-time information, the RPI approach is sometimes presented as a “thermodynamic” method,42 and therefore, can be seen as an extension of the Eyring transition state theory into the deep tunneling regime (i.e., extension for temperatures below Tc°).

We consider a system coupled to a harmonic bath, which leads to the following modified RP potential:

UPsb=UPsys+j=1Nbk=1PμjωP22(xj(k)xj(k+1))2+i=13Nk=1Pμjωj22xj(k)fij(q(k))μjωj22,
(8)

where UPsys is the system RP potential given by Eq. (2), q(k)={q1(k),q2(k),,qN(k)}; Nb is the total number of bath modes; xj(k)represents the coordinate of the kth replica of the jth bath mode; and μj and ωj are its mass and frequency, respectively. The function fij determines the coupling between the jth bath mode and the ith system degree of freedom, even though it is, in principle, a function of all the system degrees of freedom.43 

The harmonic bath can be completely characterized by a second-rank tensor known as spectral density whose components are given by

Jil(q,ω)=π2j=1Nbfij(q)ql21μjω(δ(ωωj)+δ(ω+ωj)),
(9)

and the time- and position-dependent friction tensor, which will be an important quantity in this paper, can be expressed in terms of the spectral density as

ηil(q,t)=1πdωJil(q,ω)ωcos(ωt),
(10)

where it is understood that the previous equation is valid only for t ≥ 0.42 

For reasons that will become clear later, we write the Laplace transform of ηil(t) as

η̃il(q,λ)=0eλtηil(q,t)dt=1πdωJil(q,ω)ωλω2+λ2.
(11)

The derivation of the expression for the effective ring polymer potential involves a coordinate transformation from the Cartesian representation to the normal modes of the free RP and a later Gaussian integral as detailed below. The quantum canonical partition function, Z, for such a system can be related to a classical partition function, ZP, as44 

Z=limPZP,
(12)

with

ZP=m2πβP2P/2jNbμj2πβP2P/2×dq(1)dq(P)jNbdxj(1)dxj(P)eβPUPsb,
(13)

where UPsb is the RP potential of Eq. (8).

We now perform a unitary transformation into the RP normal modes space45 to get

UPsb=l=P/2+1P/2i=13N12miωl2(Qi(l))2+k=1PV(q1(k),,q3N(k))+j=1Nbl=P/2+1P/212μjωl2(Xj(l))2+i=13N12μjωj2Xj(l)fij(l)μjωj22,
(14)

where Xj(l)and Qi(l)represent the coordinates in the normal mode space and fij(l)=k=1PClkfij(q(k))represents the RP transformed system–bath coupling with C being the RP normal mode transformation matrix (see  Appendix A) and ωl = 2ωP sin(|l|π/P) . In the previous expression, an even number of replicas, P, has been assumed. It is straightforward to treat an odd number of beads, but more involved and not necessary for the present derivation.

To perform an integration over the bath degrees of freedom, it is convenient to rewrite the previous equation as follows:

UPsb=l=P/2+1P/2i=13N12miωl2(Qi(l))2+k=1PV(q1(k),,q3N(k))+l=P/2+1P/2j=1Nbi=13N12(fij(l))2ωl2μj(ωj2+ωl2)+l=P/2+1P/2j=1Nbi=13N12μj(ωj2+ωl2)Xj(l)fij(l)μj(ωj2+ωl2)2.
(15)

Introducing Eq. (15) into Eq. (13) and integrating over the bath degrees of freedom yield

ZP=ZPbathm2πβP2P/2dq(1)dq(P)eβPUPeff,
(16)

where

ZPbath=jNbl=P/2+1P/21βP(ωj2+ωl2),
(17)

which converges to the harmonic oscillator partition function jN(2sinh(βωj/2))1 in the limit of P → ∞,46 and the effective RP potential is given by

UPeff=UPsys+l=P/2+1P/2i=13Nj=1Nb12ωl2μj(ωj2+ωl2)ωj2(f̃ij(l))2.
(18)

From this point, an expression of the rate can be obtained in analogy to Eq. (4) by replacing the SP/ by SPeff/=βPUPeff and computing ZPv and ZPr consistently, by using the same effective RP potential.

We shall refer to this formulation as RPI with explicit friction (RPI-EF). The discretized ring polymer formulation that we present in this paper is theoretically equivalent to previous formulations proposed in Refs. 42, 43, and 47; however, it presents several advantages: it allows for a more intuitive analysis, it is mathematically simpler, and it is computationally more efficient. We note that related methodologies have been proposed earlier in the literature.48–51 

Equation (18) exactly shows how, according to quantum mechanics, the effect of the bath modifies time-independent equilibrium properties, while in the classical limit (i.e., P = 1 and ωl=0 = 0), the bath contribution to these properties becomes zero. This formulation does not account for the dynamical effects of the bath on the system, which makes it particularly suitable to be combined with the ring polymer instanton method. We note also that even though the random force is a crucial element in the MDEF approach,52 rooted in the second fluctuation dissipation theorem,53 it does not appear in the RPI-EF theory due to the lack of real-time trajectories. Next, we consider different possible forms of the system–bath coupling.

The first type of coupling considered is a linear coupling given by

fij(q)=cjqi.
(19)

Using this expression, we can write Eq. (18) as

UPeff=UPsys+l=P/2+1P/2i=13Nj=1Nb12ωk2cj2μj(ωj2+ωk2)ωj2(Qi(l))2=UPsys+l=P/2+1P/2i=13N1πdωJii(ω)ωωl2ω2+ωl2(Qi(l))2=UPsys+l=P/2+1P/2i=13Nη̃ii(ωl)ωl2(Qi(l))2,
(20)

where we used Eq. (9) in the second line, and Eqs. (11) and (19) in the last line. Importantly, as a consequence of Eq. (9), a linear coupling function results in a friction tensor that is position independent.

In Fig. 1, we show a cartoon representation of the ring polymer corresponding to one free particle and to the same particle coupled to a harmonic bath. The second term on the right-hand side of Eq. (20), the “friction springs,” couple the beads beyond the nearest neighbors, representing a term that is non-local in the imaginary time (see further discussion in  Appendix B). For a positive-definite friction tensor, these friction spring-terms increase the effective coupling among the beads, causing the system to behave more classically when compared to the free-particle system. Indeed, in the low friction limit, the zero-point energy (ZPE) of a damped harmonic oscillator with frequency ω, EZPEω,η, decreases as

EZPEω,ηEZPEω,0η̃0ωD4ω+η̃(0)2πlnωDω<0,
(21)
FIG. 1.

Schematic representing the ring polymer of (a) a free particle and (b) a particle in contact with a harmonic bath. The standard harmonic springs are depicted in black, while the springs emerging from the system–bath coupling are depicted in orange. These two couplings correspond to the second and the last term on the right-hand side of Eq. (20), respectively. Note that, for the sake of clarity, we only draw the latter term for a single bead, but it is present between all pairs of beads.

FIG. 1.

Schematic representing the ring polymer of (a) a free particle and (b) a particle in contact with a harmonic bath. The standard harmonic springs are depicted in black, while the springs emerging from the system–bath coupling are depicted in orange. These two couplings correspond to the second and the last term on the right-hand side of Eq. (20), respectively. Note that, for the sake of clarity, we only draw the latter term for a single bead, but it is present between all pairs of beads.

Close modal

where a spectral density with the Drude cutoff JDrude(ω)=ηωωD2/(ωD2+ω2) was assumed in the derivation.54 

Next, we show that when the friction is position independent, we can derive an extension of the Grote–Hynes approximation for the reaction rate in the deep-tunneling regime.

1. Extension of the Grote–Hynes approximation into the deep-tunneling regime

Grote, Hynes,55,56 and Pollak57 showed that for intermediate to strong friction values, the classical reaction rate of the system–bath model can be written in terms of the system rate and η̃(ω)as

kTST(η̃)=kTSTsysη̃(ω0)2mω2+11/2η̃(ω0)2mω,
(22)

where kTSTsys and kTST are the transition state rates for the system and system bath, respectively, mis the system mass, ω is the imaginary frequency of the system at the barrier top, and ω0 is given by the following relation:

ω0ω=η̃(ω0)2mω2+11/2η̃(ω0)2mω.
(23)

As elegantly proved by Pollak,57 Eq. (23) can be interpreted as a renormalized effective barrier frequency due to dissipation. In a similar spirit, we would like to derive a relation between the RPI rates of the system without dissipation and the RPI-EF rates, which include dissipation. Moreover, it would be desirable to obtain such a relation without resorting to any assumption on the potential energy surface. The last condition forbids any direct relation between kinst(β, η) and kinst(β, η = 0) at the same temperature, since both instanton pathways will have different extensions and, therefore, will be affected by different regions of the potential energy surface. Another possibility to tackle this problem is to ask the following question: “Given an instanton obtained at Ta on UPsys, at which temperature Tb will an instanton obtained on UPeff present (approximately) the same geometry?” Mathematically, given βa = 1/kBTa and qa the solution to

0=SPsys(βa,qa)qi(k),k=1,,Pi=1,,3N,
(24)

we aim to find βb such that qa is an approximate solution to

0=SPeff(βb,qa)qi(k),k=1,,P,i=1,,3N.
(25)

In the previous equations, the sub-indices a and b refer to the inverse temperatures βa and βb, respectively. Using Eq. (1) and going to the RP normal mode representation, we look for βb that simultaneously satisfies

0=Ssys(βa,qa)/Qi(l)=βPaV(qa)Qi(l)+miωla2Qi(l),0=V(qa)Qi(l)+miωla2Qi(l),
(26)

and

0=Seff(βb,qa)/Qi(l)=βPbV(qb)Qi(l)+miωlb2Qi(l)+η̃(ωl)Qi(l),0=V(qa)Qi(l)+miωlb2Qi(l)+η̃(ωlb)ωlbQi(l),
(27)

for l = −P/2 + 1, …, P/2 and i = 1, …, 3N.

We combine Eqs. (26) and (27) to get

ωlb2+η̃(ωlb)miωlbωla2=0,
(28)

where Qi(l)0 was assumed since in that case Eqs. (26) and (27) became identical and, therefore, the latter is trivially satisfied when the former is. In the limit of P → ∞, we have ωl=2π|l|β, so we can solve the quadratic equation for Tb to obtain

Tb/Tc°=η̃(ωlb)2miω21l2+TaTc°2η̃(ωlb)2miω1l,
(29)

where Tc° is given by Eq. (7) and l ≠ 0. The previous equation cannot be fulfilled for all values of l unless η̃ has a very specific frequency dependence. In Paper II, we will see that we can exploit the fact that the normal modes with |l| = 1 dominate the instanton pathways, and thus, we can solve Eq. (29) for |l| = 1. Note that this equation has to be solved self-consistently, since ωlb depends on Tb. Equation (29) allows one to compute the tunneling rates of a system coupled to a bath by finding the instanton pathways of the uncoupled system at the scaled temperature Tb. Thus, it can be interpreted as a generalization of the GH equation into the deep tunneling regime.

The simplest system–bath coupling function that results in a position dependence of the friction is given by

fij(q)=cjgi(q),
(30)

leading to the following spectral density:

Jil(q,ω)=gi(q)ql2π2j=1Nbcj2μjω(δ(ωωj)+δ(ω+ωj)).
(31)

This coupling function is equivalent to assuming that the zero-frequency value of the friction tensor, η̃(q,ω=0), is position-dependent and its frequency dependence is identical for all positions.58 Thus, it is sometimes referred to as “separable coupling,” and can be shown to yield a lower limit for the tunneling rate.43 The RP potential in this scenario becomes

UPeff=UPsys+l=P/2+1P/2i=13Nj=1Nb12cj2ωl2μj(ωj2+ωl2)ωj2k=1PClkg(q(k))i2=UPsys+l=P/2+1P/2i=13Nj=1Nb12cj2ωl2μj(ωj2+ωl2)ωj2×k=1PClk01dsdg(qk(s))idsk=1PClkg(qref)i2=UPsys+l=P/2+1P/2i=13Nj=1Nb12cj2ωl2μj(ωj2+ωl2)ωj2×k=1PClk01dsdg(qk(s))ids2,
(32)

where we consider that qk(s):RR3N is a parameterization, such that qk(0) = qref and qk(1) = qk. In the last line, we used that k=1PClkg(qref)i only contributes to the l = 0 term and, since ω0 = 0, its contribution vanishes. As a consequence, the reference position is a free parameter that does not affect the results.

We continue by applying the chain rule,

UPeff=UPsys+l=P/2+1P/2i=13Nk=1PClk01dsrqrk(s)s×j=1Nb12cj2μjωj2g(qk(s))iqrk2ωl2(ωj2+ωl2)1/22,
(33)

and finally, by rearranging the terms, we obtain

UPeff=UPsys+l=P/2+1P/2i=13Nωl2×k=1PClkr=13Nqrefq(k)η̃ir(q,ωk)1/2dqr2=UPsys+l=P/2+1P/2i=13Nωl2×k=1PClkqrefq(k)ηĩ(q,ωk)1/2dq2,
(34)

with η̃i being the ith row of the friction tensor. We note that a straightforward extension of the Grote–Hynes approximation is not possible in this case. For completeness, for a one-dimensional system, the previous equation simplifies to

UPeff=UPsys+l=P/2+1P/2ωl2k=1PClkqrefq(k)dqη̃(q,ωl)1/22.
(35)

Naturally, the coupling of the bath to the system impacts the nuclear tunneling. One can study, for example, how the tunneling cross over temperature is modified by the coupling to the bath. A trivial stationary point on the extended phase space of the ring polymer in the pathway that connects reactants and products can be found by locating all the beads at the top of the barrier. For a 1D system with position-dependent or independent friction under the parabolic barrier approximation, one can write

λl=ωl2+η̃(ωl)mωl(ω)2,
(36)

where ωl = 2ωP sin(|l|π/P) are the free RP normal mode frequencies, is the imaginary frequency at the barrier top, and η̃ has been evaluated at the barrier top. In the limit of large P, ωl = 2π|l|/βℏ and the lowest three frequencies are

λ0=iω,λ±1=4π2β22+2πη̃(ω1)βm(ω)2,
(37)

where ω1 refers to the first Matsubara frequency,59 which depends on the temperature. The cross over temperature is the temperature, βcsb=1/kBTcsb, below which λ±1 becomes imaginary (i.e., λ±12 becomes negative) and the location of the first-order saddle point is not at the top of the barrier, i.e., a non-trivial instanton pathway becomes possible. By taking λ±1 = 0 in the previous equation and solving the quadratic equation for βcsb, one obtains

kBTcsb=1βcsb=4πη̃(ω1)2m2+4ω2η̃(ω1)2m2,
(38)

where ω1 is evaluated at Tcsb. Finally, identifying Tc° [Eq. (7)] in the equation above leads to

Tcsb=T°c×η̃(ω1)2mω2+1η̃(ω1)2mω.
(39)

Since ω1 depends on Tcsb, Eq. (39) has to be solved self-consistently. The number between square brackets is always positive and less than 1, so Tcsb is always lower than Tc°. Moreover, it is straightforward to see that the stronger the friction, the lower Tcsb becomes, and tunneling becomes less important at a given temperature. If the friction tensor is position-dependent, Tcsb can be calculated by using Eq. (39) replacing η̃(ω1) by η̃(q,ω1), where q refers to the transition state geometry. Interestingly, Eq. (39) results from Eq. (29) for Ta = Tc°, so the former equation can be interpreted as a special case of the latter.

For systems in which the ground electronic state can be approximated by effectively independent quasi-particles such as the ones obtained with Kohn–Sham (KS) density-functional theory (DFT), the adiabatic electronic friction tensor can be obtained from first-principles simulations assuming non-interacting electrons and, as shown in  Appendix C, adopts the following form for t > 0:

ηijel(q,t)=ν,νψν|iψνψν|jψνΩνν×(f(εν)f(εν))cos(Ωννt),
(40)

where f(ɛ) is the state occupation given by the Fermi–Dirac occupation function, Ωνν = (ɛνɛν)/; ψν and ɛν are the KS electronic orbitals and orbital energies of the νth level, respectively; i and j label the nuclear degrees of freedom; and i = /∂qi. A Fourier transform of the expression above leads to the usual expression employed in Refs. 13, 14, and 60 and reads

η̂ijel(q,ω)=πν,νψν|iψνψν|jψν×(f(εν)f(εν))×Ωννδ(Ωννω),
(41)

where the k-point dependence has been omitted.

Most of the applications of the MDEF approach only consider an electronic friction tensor that is local in time to avoid the complexities of handling a non-instantaneous memory kernel and normally invoke the Markov approximation. This limit is also often referred to in the literature as the quasi-static limit since the Markov approximation is normally realized by taking the ω → 0 limit in Eq. (41).29 In the cases where the system presents a constant density of states (DOS) around the Fermi level, an equivalent derivation is possible by applying the constant coupling approximation.13 The quasi-static limit involves the evaluation of the friction tensor in Eq. (41) for excitations infinitesimally close to the Fermi level. In practical calculations with finite k-point grids, it is numerically challenging to accurately describe the DOS at the Fermi energy. This is typically circumvented by introducing a finite width for the delta function in Eq. (41). The choice of width depends on the system and, in the literature, values between 0.01 and 0.60 eV can be found.14,61–63

The connection of the ab initio electronic friction and the RPI-EF rate theory might seem a trivial substitution of Eq. (41) into Eq. (34). However, as we shall show in Sec. IV A, this is not the case and in order to obtain a better connection between the electronic friction and the system–bath model used in the formulation of RPI-EF, a different expression should be employed.

Starting from Eq. (40), and in a similar spirit to Eq. (10), we perform a Laplace transform to get

η̃ijel(q,λ)=0dteλtηijel(q,t)=ν,νψν|iψνψν|jψν(f(εν)f(εν))Ωννλλ2+Ωνν2.
(42)

The equation above adopts the same limit for λ → 0 as Eq. (41). However, for λ > 0, instead of the δ function, we obtain a sum of Lorentzian functions of width 2λ. By comparingEqs. (42) and (11), we can identify the equivalent of the spectral density in RPI-EF as follows:

Jij(q,ω)=πν,νψν|iψνψν|jψνω2×(f(εν)f(εν))δ(ωΩνν),
(43)

which provides a seamless connection between RPI-EF and electronic friction.

In RPI-EF, the spectral density of electronic friction shown in Eq. (43) is evaluated simultaneously at the ring polymer normal mode frequencies. Thus, we have derived viable expressions to combine RPI-EF with an electronic friction formulation that can be calculated from first-principles, without any further approximations, except for the assumption of separable coupling, which we shall examine for real systems in Paper II. We note that, as discussed in Ref. 13, the electronic-friction formalism itself is only well-defined under the assumption that the system–bath coupling is linear in the coupling parameter, making the separable coupling considered here a rather general form that makes both RPI-EF and the electronic friction applicable. We expect that the connection of these two theories will be suitable as long as the approximations of both underlying theories remain valid.

We have presented an extension of the ring polymer instanton rate theory to describe a system coupled to a bath of harmonic oscillators through the definition of an effective friction tensor that enters the instanton ring polymer potential energy expression, and therefore, we refer to this method as the RPI-EF approach. The theory is rather general and allows the inclusion of frequency and position dependence in the system–bath coupling for the calculation of thermal tunneling rates within the instanton approximation. For the case of linear coupling, we derived an approximation that allows one to predict RPI-EF reaction rates by using only RPI calculations. The approximation can be understood as an extension of the Grote–Hynes approximation to the deep tunneling regime. This may be useful to estimate whether it is necessary to carry out full RPI-EF calculations for a particular reaction.

RPI-EF is a method tailored for the description of tunneling rates and based on imaginary-time trajectories. Therefore, it cannot be applied for the simulation of vibrational relaxation or scattering experiments,2,21,64 where some kinds of NQEs and NAEs could interplay strongly. It would be interesting to write similar extensions to approaches based on path integral molecular dynamics,65–68 since they would yield efficient approximations to model these situations. We hope that the derivations presented in this work stimulate further theoretical developments in this area and allow new phenomena to be explained in situations that we have not yet explored.

Y.L., E.S.P., and M.R. acknowledge financial support 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 Ph.D. studentship. R.J.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 (No. RPG-2019-078) and the UKRI Future Leaders Fellowship program (Grant No. MR/S016023/1).

The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article.

The free ring polymer potential is given by setting V = 0 in Eq. (2). The resulting potential is harmonic; however, due to the presence of degenerate eigenvalues, there is no unique transformation to diagonalize it. Assuming P is even, one possibility is the following orthogonal coordinate transformation:45,69

Qi(l)=k=1PClk(P)qiki=1,,3N,l=P/2+1,,P/2,
(A1)

where the P × P matrix C(P) is defined as

Clk(P)=1P,l=0,2Pcos2πkln,1lP/21,1P(1)j,l=P/2,2Psin2πklP,P/2+1l1.

The mean-field RP potential in the Cartesian representation for the linear coupling case is obtained by introducing Eq. (A1) into (20) leading to

UPeff=l=1PV(q1(l),,q3N(l))+k=1k=1Pi=13N12qi(k)Oi,k,kqi(k)+12qi(k)Di,k,kqi(k),
(B1)

with

Ok,k=miωP2(2δk,kδk,k1δk,k+1)
(B2)

and

Dk,k=η0ωPl=1(P1)/24Psin(πl/P)cos(2πl(kk)/P)2P(1)k+k,
(B3)

where a linear (Ohmic) spectral density was considered. In Fig. 2, we show a graphical representation of the spring coupling matrices, Ok,k and Dk,k. The latter has small but non-zero matrix elements outside the tridiagonal entries present in the former, representing coupling beyond nearest-neighbor beads. Moreover, the matrix elements decay rapidly with the increase of the bead index distances when periodic boundary conditions are considered.

FIG. 2.

Spring coupling matrices Ok,k (left) and Dk,k (right) at 25 K for 15 beads and η0/0 = 1.0. The values have been scaled by P to ease visual comparison.

FIG. 2.

Spring coupling matrices Ok,k (left) and Dk,k (right) at 25 K for 15 beads and η0/0 = 1.0. The values have been scaled by P to ease visual comparison.

Close modal

The steps presented in this appendix closely follow Ref. 60, and we repeat them here merely for completeness. We consider a quadratic electronic Hamiltonian of the following form:

ĥ=pqhqp(q)d̂p+d̂q,
(C1)

where hqp(q) is the general notation to represent that matrix elements might depend on the nuclear degrees of freedom, q, and d̂p+ and d̂q are the electronic creation and annihilation operators, respectively.

Starting from the quantum–classical Liouville equation, the electronic friction tensor in the adiabatic limit with nuclei fixed at position q can be written as60 

ηijel(q,t)=treiĥeiL̂̂tjρ̂ss,
(C2)

where L̂̂ is the Liouvillian superoperator, tre implies tracing over the electronic degrees of freedom, ρ̂ss is the steady state electronic density matrix, and i and j represent two nuclear degrees of freedom.

By recalling that eiLt̂̂()=eiĥt/()eiĥt/, the invariance of the trace under cyclic permutations, and considering the quadratic Hamiltonian presented above, we find

ηijel(q,t)=treeiĥtiĥeiĥtjρ̂ss=mnihnmtreeiĥtd̂m+d̂neiĥtjρ̂ss=mnihnmtreeiĥtd̂m+eiĥteiĥtd̂neiĥtjρ̂ss.
(C3)

By noting that eiĥt/d̂m+eiĥt/=a(eiĥt/)mad̂a+ and eiĥt/d̂neiĥt/=b(eiĥt/)bnd̂b (see the supplementary material in Ref. 60), and defining σabss=tre(d̂a+d̂bρ̂ss), Eq. (C3) can be expressed as

ηijel(q,t)=mnabihnm(eiĥt/)ma(eiĥt/)bntred̂a+d̂bjρ̂ss=mnabihnm(eiĥt/)majσabss(eiĥt/)bn=niĥeiĥt/jσsseiĥt/nn,
(C4)

where n represents a sum over electronic orbitals. If we take the basis in which ĥ is diagonal,

ĥ=νενψνψν.
(C5)

Equation (C4) simplifies to

ηijel(q,t)=ν,νψνiĥψνeiενt/ψνjσssψνeiενt/.
(C6)

At equilibrium, we can write σ̂ss as

σ̂ss=νf(εν)ψνψν,
(C7)

and, therefore, the second matrix element in Eq. (C6) can be evaluated as

ψνiσ̂ssψν=iενf(εk)ενδνν+(f(εν)f(εν))ψν|iψν.
(C8)

Introducing Eq. (C8) into Eq. (C6) and using

iψνĥψν=iενδνν=ψνiĥψν+ψν|iψν(ενεν)
(C9)

leads to

ηijel(q,t)=ν,νψν|iψν(ενεν)+iενδνν×jενf(εν)ενδνν+(f(εν)f(εν))ψν|jψν×ei(ενεν)t/.
(C10)

Finally, upon noticing that the factors in front of the exponential are invariant under the exchange of orbital labels such that the imaginary part of the complex exponential vanishes, and the fact that for ν = ν′, the expression is zero, one arrives to Eq. (40) presented in the main text. This expression agrees with an alternative recent derivation based on the exact factorization of the electronic–nuclear wavefunction.70 

1.
C.
Bartels
,
R.
Cooper
,
D. J.
Auerbach
, and
A. M.
Wodtke
,
Chem. Sci.
2
,
1647
(
2011
).
2.
A. M.
Wodtke
,
Chem. Soc. Rev.
45
,
3641
(
2016
).
3.
D. J.
Auerbach
,
J. C.
Tully
, and
A. M.
Wodtke
,
Nat. Sci.
1
(
2021
).
4.
B.
Schindler
,
D.
Diesing
, and
E.
Hasselbrink
,
J. Chem. Phys.
134
,
034705
(
2011
).
5.
O.
Cohen
,
G.
Bartal
,
H.
Buljan
,
T.
Carmon
,
J. W.
Fleischer
,
M.
Segev
, and
D. N.
Christodoulides
,
Nature
433
,
503
(
2005
).
6.
C. L.
Box
,
Y.
Zhang
,
R.
Yin
,
B.
Jiang
, and
R. J.
Maurer
,
J. Am. Chem. Soc. Au
1
,
164
(
2021
).
7.
M.
Persson
and
B.
Hellsing
,
Phys. Rev. Lett.
49
,
662
(
1982
).
8.
J. C.
Tully
,
M.
Gomez
, and
M.
Head‐Gordon
,
J. Vac. Sci. Technol. A
11
,
1914
(
1993
).
9.
S. P.
Rittmeyer
,
J.
Meyer
,
J. I.
Juaristi
, and
K.
Reuter
,
Phys. Rev. Lett.
115
,
046102
(
2015
).
10.
W.
Dou
and
J. E.
Subotnik
,
J. Chem. Phys.
144
,
024116
(
2016
).
11.
N.
Shenvi
,
S.
Roy
, and
J. C.
Tully
,
J. Chem. Phys.
130
,
174107
(
2009
).
12.
I. G.
Ryabinkin
and
A. F.
Izmaylov
,
J. Phys. Chem. Lett.
8
,
440
(
2017
); arXiv:1611.06525.
13.
M.
Head-Gordon
and
J. C.
Tully
,
J. Chem. Phys.
103
,
10137
(
1995
).
14.
R. J.
Maurer
,
M.
Askerka
,
V. S.
Batista
, and
J. C.
Tully
,
Phys. Rev. B
94
,
115432
(
2016
).
15.
M.
Blanco-Rey
,
J. I.
Juaristi
,
R.
Díez Muiño
,
H. F.
Busnengo
,
G. J.
Kroes
, and
M.
Alducin
,
Phys. Rev. Lett.
112
,
103203
(
2014
).
16.
W.
Dou
and
J. E.
Subotnik
,
J. Chem. Phys.
148
,
230901
(
2018
).
17.
A. J.
Coffman
and
J. E.
Subotnik
,
Phys. Chem. Chem. Phys.
20
,
9847
(
2018
).
18.
Y.
Zhang
,
R. J.
Maurer
,
H.
Guo
, and
B.
Jiang
,
Chem. Sci.
10
,
1089
(
2019
).
19.
P.
Spiering
,
K.
Shakouri
,
J.
Behler
,
G.-J.
Kroes
, and
J.
Meyer
,
J. Phys. Chem. Lett.
10
,
2957
(
2019
).
20.
A.
Kandratsenka
,
H.
Jiang
,
Y.
Dorenkamp
,
S. M.
Janke
,
M.
Kammler
,
A. M.
Wodtke
, and
O.
Bünermann
,
Proc. Natl. Acad. Sci. U. S. A.
115
,
680
(
2018
).
21.
O.
Bünermann
,
H.
Jiang
,
Y.
Dorenkamp
,
A.
Kandratsenka
,
S. M.
Janke
,
D. J.
Auerbach
, and
A. M.
Wodtke
,
Science
350
,
1346
(
2015
).
22.
K. M.
Forsythe
and
N.
Makri
,
J. Chem. Phys.
108
,
6819
(
1998
).
23.
H.
Kimizuka
,
S.
Ogata
, and
M.
Shiga
,
Phys. Rev. B
100
,
024104
(
2019
).
24.
M.
Rossi
,
M.
Ceriotti
, and
D. E.
Manolopoulos
,
J. Phys. Chem. Lett.
7
,
3001
(
2016
).
25.
H.
Kimizuka
and
M.
Shiga
,
Phys. Rev. Mater.
5
,
065406
(
2021
).
26.
W.
Fang
,
J. O.
Richardson
,
J.
Chen
,
X.-Z.
Li
, and
A.
Michaelides
,
Phys. Rev. Lett.
119
,
126001
(
2017
).
27.
H.
Kimizuka
,
S.
Ogata
, and
M.
Shiga
,
Phys. Rev. B
97
,
014102
(
2018
).
28.
J. O.
Richardson
,
Int. Rev. Phys. Chem.
37
,
171
(
2018
).
29.
B.
Hellsing
and
M.
Persson
,
Phys. Scr.
29
,
360
(
1984
).
30.
J. O.
Richardson
and
S. C.
Althorpe
,
J. Chem. Phys.
131
,
214106
(
2009
).
31.
A.
Arnaldsson
, “
Calculation of quantum mechanical rate constants directly from ab initio atomic forces
,” Ph.D. thesis,
University of Washington
,
2007
.
32.
W. H.
Miller
,
J. Chem. Phys.
62
,
1899
(
1975
).
33.
E.
Grifoni
,
G.
Piccini
, and
M.
Parrinello
,
Proc. Natl. Acad. Sci. U. S. A.
116
,
4054
(
2019
).
34.
I.
Affleck
,
Phys. Rev. Lett.
46
,
388
(
1981
).
35.
S. C.
Althorpe
,
J. Chem. Phys.
134
,
114104
(
2011
).
36.
J. O.
Richardson
,
J. Chem. Phys.
144
,
114106
(
2016
).
37.
J. B.
Rommel
,
T. P. M.
Goumans
, and
J.
Kästner
,
J. Chem. Theory Comput.
7
,
690
(
2011
).
38.
J. B.
Rommel
,
Y.
Liu
,
H.-J.
Werner
, and
J.
Kästner
,
J. Phys. Chem. B
116
,
13682
(
2012
).
39.
Y.
Litman
,
J. O.
Richardson
,
T.
Kumagai
, and
M.
Rossi
,
J. Am. Chem. Soc.
141
,
2526
(
2019
).
40.
Y.
Litman
and
M.
Rossi
,
Phys. Rev. Lett.
125
,
216001
(
2020
).
41.
M. J.
Gillan
,
J. Phys. C: Solid State Phys.
20
,
3621
(
1987
).
42.
U.
Weiss
,
Quantum Dissipative Systems
, 3rd ed. (
World Scientific
,
2008
).
43.
A. O.
Caldeira
and
A. J.
Leggett
,
Ann. Phys.
149
,
374
(
1983
).
44.
R. P.
Feynman
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
, International Series in Pure and Applied Physics (
McGraw Hill
,
New York
,
1965
).
45.
T. E.
Markland
and
D. E.
Manolopoulos
,
J. Chem. Phys.
129
,
024105
(
2008
).
46.
H.
Kleinert
,
Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets
, 5th ed. (
World Scientific
,
2009
).
47.
P.
Hänggi
,
P.
Talkner
, and
M.
Borkovec
,
Rev. Mod. Phys.
62
,
251
(
1990
).
48.
S.
Ranya
and
N.
Ananth
,
J. Chem. Phys.
152
,
114112
(
2020
).
49.
J. E.
Lawrence
,
T.
Fletcher
,
L. P.
Lindoy
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
151
,
114119
(
2019
).
50.
J. E.
Lawrence
and
D. E.
Manolopoulos
,
J. Chem. Phys.
152
,
204117
(
2020
).
51.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
105
,
6856
(
1996
).
52.
N.
Hertl
,
R.
Martin-Barrios
,
O.
Galparsoro
,
P.
Larrégaray
,
D. J.
Auerbach
,
D.
Schwarzer
,
A. M.
Wodtke
, and
A.
Kandratsenka
,
J. Phys. Chem. C
125
,
14468
(
2021
).
53.
R.
Kubo
,
Rep. Prog. Phys
29
,
255
(
1966
).
54.
G.-L.
Ingold
, “
Path integrals and their application to dissipative quantum systems
,” in
Coherent Evolution in Noisy Environments
, edited by
A.
Buchleitner
and
K.
Hornberger
(
Springer
,
Berlin, Heidelberg
,
2002
), pp.
1
53
.
55.
R. F.
Grote
and
J. T.
Hynes
,
J. Chem. Phys.
73
,
2715
(
1980
).
56.
R. F.
Grote
and
J. T.
Hynes
,
J. Chem. Phys.
74
,
4465
(
1981
).
57.
E.
Pollak
,
J. Chem. Phys.
85
,
865
(
1986
).
58.
J. B.
Straus
,
J. M.
Gomez Llorente
, and
G. A.
Voth
,
J. Chem. Phys.
98
,
4082
(
1993
).
59.
T.
Matsubara
,
Prog. Theor. Phys.
14
,
351
(
1955
).
60.
W.
Dou
,
G.
Miao
, and
J. E.
Subotnik
,
Phys. Rev. Lett.
119
,
046001
(
2017
).
61.
D.
Novko
,
J. C.
Tremblay
,
M.
Alducin
, and
J. I.
Juaristi
,
Phys. Rev. Lett.
122
,
016806
(
2019
).
62.
R. J. M.
Connor
,
L.
Box
, and
W. G.
Stark
, arXiv:2112.00121 (
2021
).
63.
A. M.
Shipley
,
M. J.
Hutcheon
,
M. S.
Johnson
,
R. J.
Needs
, and
C. J.
Pickard
,
Phys. Rev. B
101
,
224511
(
2020
).
64.
H.
Jiang
,
X.
Tao
,
M.
Kammler
,
F.
Ding
,
A. M.
Wodtke
,
A.
Kandratsenka
,
T. F.
Miller
, and
O.
Bünermann
,
J. Phys. Chem. Lett.
12
,
1991
(
2021
).
65.
I. R.
Craig
and
D. E.
Manolopoulos
,
J. Chem. Phys.
121
,
3368
(
2004
).
66.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
100
,
5106
(
1994
).
67.
M.
Rossi
,
M.
Ceriotti
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
140
,
234116
(
2014
).
68.
G.
Trenins
,
M. J.
Willatt
, and
S. C.
Althorpe
,
J. Chem. Phys.
151
,
054109
(
2019
).
69.
I. R.
Craig
, “
Ring polymer molecular dynamics
,” Ph.D. thesis,
University of Oxford
,
2006
.
70.
R.
Martinazzo
and
I.
Burghardt
, arXiv:2108.02622 (
2021
).