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.

## I. INTRODUCTION

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 past^{1–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 present^{1,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 Persson^{29} 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.

## II. RING POLYMER INSTANTON RATE THEORY

The RPI approximation^{30,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 Althorpe^{35} to be equivalent to the original derivation by Miller^{32} 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, *S*_{P}, is related to the potential energy of the ring polymer, *U*_{P}, by

with

Here, $qi(k)$ is the position of the *i*th degree of freedom of the *k*th replica, *m*_{i} is the mass of the *i*th 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 $\omega P=(\beta P\u210f)\u22121$ with $\beta P=1kBPT$, where

*k*

_{B}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

*U*

_{P}. 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 as^{28}

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\u0304$ 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

with

and

In the expression above, *λ*_{k} represent the *P* × 3*N* eigenvalues of the mass scaled ring polymer Hessian defined by the second derivatives of *U*_{P} 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, *T*_{c}°, which, in most cases, can be estimated by a parabolic barrier approximation as

where *ω*^{‡} represents the imaginary frequency at the barrier top between reactants and products. At temperatures below *T*_{c}°, the reactive process is dominated by tunneling, while above *T*_{c}°, 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 *T*_{c}°).

## III. RING POLYMER INSTANTON RATE THEORY WITH EXPLICIT FRICTION

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

where $UPsys$ is the system RP potential given by Eq. (2), $q(k)={q1(k),q2(k),\u2026,qN(k)}$; *N*_{b} is the total number of bath modes; $xj(k)$represents the coordinate of the *k*th replica of the *j*th bath mode; and *μ*_{j} and *ω*_{j} are its mass and frequency, respectively. The function *f*_{ij} determines the coupling between the *j*th bath mode and the *i*th 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

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

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

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, *Z*_{P}, as^{44}

with

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

We now perform a unitary transformation into the RP normal modes space^{45} to get

where $Xj(l)$and $Qi(l)$represent the coordinates in the normal mode space and $fij(l)=\u2211k=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:

where

which converges to the harmonic oscillator partition function $\u220fjN(2sinh(\beta \omega j/2))\u22121$ in the limit of *P* → ∞,^{46} and the effective RP potential is given by

From this point, an expression of the rate can be obtained in analogy to Eq. (4) by replacing the *S*_{P}/*ℏ* by $SPeff/\u210f=\beta 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.

### A. Position-independent friction

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

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

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\omega ,\eta $, decreases as

where a spectral density with the Drude cutoff $JDrude(\omega )=\eta \omega \omega D2/(\omega D2+\omega 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 Pollak^{57} 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 $\eta \u0303(\omega )$as

where $kTSTsys$ and *k*_{TST} are the transition state rates for the system and system bath, respectively, *m*is the system mass, *ω*^{‡} is the imaginary frequency of the system at the barrier top, and *ω*_{0} is given by the following relation:

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 *k*_{inst}(*β*, *η*) and *k*_{inst}(*β*, *η* = 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 *T*_{a} on $UPsys$, at which temperature *T*_{b} will an instanton obtained on $UPeff$ present (approximately) the same geometry?” Mathematically, given *β*_{a} = 1/*k*_{B}*T*_{a} and *q*_{a} the solution to

we aim to find *β*_{b} such that *q*_{a} is an approximate solution to

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

and

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

where $Qi(l)\u22600$ 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 $\omega l=2\pi |l|\beta \u210f$, so we can solve the quadratic equation for *T*_{b} to obtain

where *T*_{c}° is given by Eq. (7) and *l* ≠ 0. The previous equation cannot be fulfilled for all values of *l* unless $\eta \u0303$ 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 $\omega lb$ depends on *T*_{b}. 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 *T*_{b}. Thus, it can be interpreted as a generalization of the GH equation into the deep tunneling regime.

### B. Position-dependent friction

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

leading to the following spectral density:

This coupling function is equivalent to assuming that the zero-frequency value of the friction tensor, $\eta \u0303(q,\omega =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

where we consider that $qk(s):R\u2192R3N$ is a parameterization, such that **q**^{k}(0) = **q**^{ref} and **q**^{k}(1) = **q**^{k}. In the last line, we used that $\u2211k=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,

and finally, by rearranging the terms, we obtain

with $\eta \u0303i$ being the *i*th 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

### C. Renormalization of cross over temperature

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

where *ω*_{l} = 2*ω*_{P} sin(|*l*|*π*/*P*) are the free RP normal mode frequencies, *iω*^{‡} is the imaginary frequency at the barrier top, and $\eta \u0303$ has been evaluated at the barrier top. In the limit of large *P*, *ω*_{l} = 2*π*|*l*|/*βℏ* and the lowest three frequencies are

where *ω*_{1} refers to the first Matsubara frequency,^{59} which depends on the temperature. The cross over temperature is the temperature, $\beta csb=1/kBTcsb$, below which *λ*_{±1} becomes imaginary (i.e., $\lambda \xb112$ 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 $\beta csb$, one obtains

where *ω*_{1} is evaluated at $Tcsb$. Finally, identifying *T*_{c}° [Eq. (7)] in the equation above leads to

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 *T*_{c}°. 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 $\eta \u0303(\omega 1)$ by $\eta \u0303(q\u2021,\omega 1)$, where *q*^{‡} refers to the transition state geometry. Interestingly, Eq. (39) results from Eq. (29) for *T*_{a} = *T*_{c}°, so the former equation can be interpreted as a special case of the latter.

## IV. *AB INITIO* ELECTRONIC FRICTION

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:

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} = *∂*/*∂q*_{i}. A Fourier transform of the expression above leads to the usual expression employed in Refs. 13, 14, and 60 and reads

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.

### A. Electronic spectral density of non-interacting electrons

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:

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.

## V. CONCLUSIONS

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.

## ACKNOWLEDGMENTS

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

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

### APPENDIX A: FREE RING POLYMER NORMAL MODES

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}

where the *P* × *P* matrix *C*^{(P)} is defined as

### APPENDIX B: EFFECTIVE RING POLYMER POTENTIAL IN CARTESIAN REPRESENTATION FOR SPATIALLY INDEPENDENT COUPLING

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

with

and

where a linear (Ohmic) spectral density was considered. In Fig. 2, we show a graphical representation of the spring coupling matrices, *O*_{k,k′} and *D*_{k,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.

### APPENDIX C: ARRIVING AT EQ. (40)

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:

where *h*_{qp}(**q**) is the general notation to represent that matrix elements might depend on the nuclear degrees of freedom, **q**, and $d\u0302p+$ and $d\u0302q$ 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 as^{60}

where $L\u0302\u0302$ is the Liouvillian superoperator, tr_{e} implies tracing over the electronic degrees of freedom, $\rho \u0302ss$ is the steady state electronic density matrix, and *i* and *j* represent two nuclear degrees of freedom.

By recalling that $e\u2212iLt\u0302\u0302(\u22c5)=e\u2212ih\u0302t/\u210f(\u22c5)eih\u0302t/\u210f$, the invariance of the trace under cyclic permutations, and considering the quadratic Hamiltonian presented above, we find

By noting that $eih\u0302t/\u210fd\u0302m+e\u2212ih\u0302t/\u210f=\u2211a(eih\u0302t/\u210f)mad\u0302a+$ and $eih\u0302t/\u210fd\u0302ne\u2212ih\u0302t/\u210f=\u2211b(e\u2212ih\u0302t/\u210f)bnd\u0302b$ (see the supplementary material in Ref. 60), and defining $\sigma abss=tre(d\u0302a+d\u0302b\rho \u0302ss)$, Eq. (C3) can be expressed as

where *∑*_{n} represents a sum over electronic orbitals. If we take the basis in which $h\u0302$ is diagonal,

Equation (C4) simplifies to

At equilibrium, we can write $\sigma \u0302ss$ as

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

leads to

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}