# Memory effects in non-adiabatic molecular dynamics at metal surfaces

###### Abstract

We study the effect of temporal correlation in a Langevin equation describing non-adiabatic dynamics at metal surfaces. For a harmonic oscillator the Langevin equation preserves the quantum dynamics exactly and it is demonstrated that memory effects are needed in order to conserve the ground state energy of the oscillator. We then compare the result of Langevin dynamics in a harmonic potential with a perturbative master equation approach and show that the Langevin equation gives a better description in the non-perturbative range of high temperatures and large friction. Unlike the master equation, this approach is readily extended to anharmonic potentials. Using density functional theory we calculate representative Langevin trajectories for associative desorption of N from Ru(0001) and find that memory effects lowers the dissipation of energy. Finally, we propose an ab-initio scheme to calculate the temporal correlation function and dynamical friction within density functional theory.

## I Introduction

Modern computational surface chemistry, as e.g. applied to heterogeneous catalysis, is largely based on the Born-Oppenheimer approximation and potential energy surfaces which are typically obtained using density functional theory (DFT).christensen () In the adiabatic approximation the electrons are assumed to follow the motion of the nuclei instantaneously and the dynamics thus becomes confined to the ground state potential energy surface. While the adiabatic approximation has certainly been successful in giving a detailed quantitative account of a range of chemical reactions on metal surfaces, it is still not clear under which general circumstances the approximation is reliable.luntz1 (); juaristi1 (); luntz3 (); juaristi2 () In particular, the role of non-adiabatic effects is often difficult to asses due to the inadequacy of low dimensional models of surface dynamics. For example, unusual sticking coefficients in the measured dissociative adsorption of N on Ru(0001) hints at strong non-adiabatic energy loss,diekhoner1 () but has been accounted for by multi-dimensional adiabatic dynamics.diaz1 (); diaz2 () For other reactions, such as associative desorption of N from Ru(0001), non-adiabatic effects still seem to be very importantmurphy (); diekhoner2 (); luntz1 () and multi-dimensional adiabatic simulations have not been able to account for large energy losses during desorption.diaz3 () Another example where adiabatic dynamics have failed is the dissociation of O on Al(111) where spin selection rules gives rise to highly non-adiabatic behavior.behler1 ()

Non-adiabatic dynamics for isolated molecules is usually handled by including the first few excited adiabatic potential energy surfaces and imposing some surface hopping scheme. When distinct adsorbate diabatic states are present there may be physical arguments why the adsorbate should remain in such a state during a reaction and the non-adiabatic dynamics can be evaluated by constraining the adsorbate to such a diabat.behler2 () This picture may then be improved by introducing surface hopping between diabats. However, for molecules adsorbed on metal surfaces there is an infinity of electronic excited states in the immediate vicinity of the ground state and surface hopping may not be the most practical scheme. Another popular and rather general method to handle non-adiabatic effects is through Langevin dynamics where electronic friction and stochastic forces account for dissipation and fluctuation as a result of coupling to excited electronic states.metiu (); schmid (); caldeira (); brandbyge () Usually, the so-called Markov approximation is employed where the fluctuating forces are not temporally correlated but can be related to the electronic friction through the fluctuation-dissipation theorem.tully () At high electronic temperatures, thermal excitations dominate the electronic system and the Markov approximation is good for describing chemical reactions mediated by hot electrons.luntz2 (); olsen5 () However, for non-adiabatic dynamics in general, the Markov approximation may fail and it then becomes important to take into account the ’memory’ of the system.

In the present paper we explore the consequences of Langevin dynamics with and without the Markov approximation at various temperatures. We follow the approach of Brandbyge et al. brandbyge () and base the analysis on a model Hamiltonian from which the electronic friction and correlation function can be derived explicitly. We start by modelling the internal stretch mode of CO adsorbed on Cu(100) by a considering a harmonic oscillator coupled to a thermal reservoir of electrons and compare the results to those obtained with a master equation approach. The general trend we see is that at low temperatures the Markov approximation overestimates the effect of dissipation. This is because the fluctuating forces are of thermal origin in the Markov approximation while the dissipative terms originate from non-thermal excitations and the relative effect of dissipation compared to fluctuations is thus increased. We also show that non-Markovian dynamics is needed in order to ensure energy conservation at low temperatures and thus maintains detailed balance between fluctuations and dissipation. Associative desorption of N from Ru(0001) is then studied using the Langevin equation on representative trajectories and again, memory effects are shown to reduce the dissipation of energy. Finally, we comment on a possible method to obtain the full correlation function and thus include memory effects in an ab-initio setting based on DFT.trail () In appendix A, we review the connection between the reduced density matrix and the Langevin equation and emphasize the probabilistic interpretation of the correlation function. In appendix B, we review how correlated stochastic forces can be randomly sampled given a discretized version of the correlation function.

## Ii Model

A commonly used electronic Hamiltonian describing of atoms or molecules adsorbed on metals surfaces is the Newns-Anderson model,newnsanderson1 (); newnsanderson2 () where the adsorbate is described by a single adsorbate state which hybridizes with metallic states and thus acquires a broadening in energy. A very simple non-adiabatic extension of this model is obtained by coupling the resonant states to an adsorbate degree of freedom and extending the Hamiltonian with a nuclear kinetic energy and adiabatic potential. Assuming a quadratic nuclear potential and linear coupling to the resonance, the Hamiltonian becomes

(1) | ||||

where is the nuclear momentum, the adsorbate effective mass, and and are creation operators for adsorbate and metallic electronic states respectively. The Hamiltonian can be thought of as a harmonic oscillator which is shifted when the resonance becomes occupied and the coupling constant is the force felt by adsorbate in this state. We will furthermore restrict ourselves to the wide band approximation in which the metallic band of electrons is assumed to be much wider than the width of the resonant state. The resonance projected density of states is then a Lorentzian:

(2) |

with full width at half maximum given by

(3) |

Due to the non-adiabatic coupling in Eq. 1, the adsorbate may exchange energy with the electronic system via the resonant state . However, usually we are only interested in the nuclear degrees of freedom and it is then convenient to trace out the electronic degrees of freedom from the full dynamics. This is accomplished by the reduced time dependent density matrix:

(4) |

where means the trace over electronic states and is the full density matrix at . Choosing an adsorbate basis , the diagonal elements of the reduced density matrix give the probabilities that the adsorbate is in a particular state at time .

### ii.1 Non-Markovian master equation

We will consider the time-dependent probability of being in a particular energy eigenstate . The equation governing these probabilities is known as a master equation and can be derived by taking the trace of the Liouville equation for the full density matrix. The result is

(5) |

where the influence functional depends on the complete history of the full density matrix. Gao gao97 () has shown how to evaluate using the Hamiltonian (1) within the self-consistent Born approximation. Furthermore, imposing the Markov approximation, where it is assumed that the influence functional only depends on the instantaneous value of the density matrix, and taking the diagonal elements of Eq. (5) led to an explicit expression for the master equation. However, using the formalism of Gao gao97 (), it is straightforward to generalize the results to a non-Markovian master equation. For completeness we state the result here which is

(6) | ||||

with the differential rates given by

(7) | ||||

where , is an eigenstate of with eigenvalue , , is the Fermi-Dirac distribution, is the projected density of states (2) and .

The Markov approximation is obtained by extending the temporal integration to infinity which is justified when is much larger than some electronic correlation time . The probabilities are then assumed to depend on t rather than and integrating over yields the usual golden rule type expression for the transition rates.gao97 () The master equation was derived assuming that off-diagonal elements of the density matrix are not important. Allthough it is straightforward to generalize the result (6) to a coherent master equation which takes off-diagonal elements into account, it has been shown that, if the initial state is not coherent, the off-diagonal elements will have very little influence on the diagonal elements.gao97 () In terms of multiple inelastic scattering events, neglecting coherency corresponds to associating a probability distribution to each scattering event.olsen1 (); olsen3 ()

In Fig. 1 we show the differential transition rate at three different temperatures. The time dependence only depends on the properties of the electronic system at the given temperature and the non-adiabatic coupling simply gives an overall scaling. Since the probabilities typically change on timescales of and the approach zero within a few femtoseconds, the Markov approximation is expected to be very good for the master equation in a large range of temperatures.

### ii.2 Non-Markovian Langevin dynamics

If we calculate the diagonal of the reduced density matrix in a basis of position eigenstates, a Langevin equation emerges. This is achieved by writing Eq. (4) as a path integral and using the Feynman-Vernon formalism of influence functionals to obtain an effective action to second order in the frictional coupling .feynman-vernon (); caldeira (); schmid () The result is given in Eq. (25). This approach is not perturbative in the same sense as the master equation where the derivation is based on a direct second order expansion of the reduced density matrix. Rather, the second order expansion of the action leads to a density matrix which contains all orders of the frictional coupling. As explained in appendix A, the result can be interpreted as a sum over classical Langevin trajectories with initial conditions sampled from the Wigner distribution of the initial state and the equation of motion is

(8) |

where the dynamical electronic friction and is a Gaussian distributed stochastic force specified by its correlation function:

(9) |

With the model Hamiltonian (1) it is possible to evaluate the friction and correlation function explicitly. The result is:brandbyge ()

(10) |

with

(11) | ||||

(12) |

and is the projected density of states Eq. (2). The correlation function is

(13) |

In Fig. 2 we show the correlation function for three different temperature. The structure and typical correlation time is very similar to the differential rate shown in Fig. 1. However, in contrast to the master equation the timescale of motion in the Langevin equation is on the order of and correlation effects may be very important at low temperatures. In general the stochastic forces at a given time depends on the state of the electronic system which again depends on the path taken by the adsorbate. This gives rise to correlation between forces at different times and this ”memory” is the price one pays for tracing out the electronic degrees of freedom. The range of memory in the system depends strongly on temperature since high temperatures tend to rapidly destroy coherence in the state of the electronic system. In the high temperature limit where , one obtains the well known expression

(14) |

where . This is the Markov approximation in which there is no correlation between forces at different times. Taking CO on Cu(100) as an example, the smallest timescale is the period of vibrational motion which is . With a standard Verlet integration one needs a timestep of to describe ground state vibrations and a first estimate of the validity of the Markov approximation is obtained as: . To get a better quantitative estimate of the validity of the Markov approximation we can consider the correlation time given by

(15) |

It should be noted that the correlation time is only a function of the electronic system and does not depend on the non-adiabatic coupling . We have calculated as a function of temperature and the result is shown in Fig. 3. For molecular dynamics requiring a time step no larger than , we see that the correlation time becomes larger than this when the temperature comes below . Thus below this temperature non-Markovian processes play an important role in the dynamics.

## Iii Results

Before we test the role of non-Markovian effects on a generic non-adiabatic surface reaction, we will compare Markovian and non-Markovian Langevin dynamics for a harmonic oscillator potential with results obtained from a master equations approach.

### iii.1 Quadratic potential

It is easy to see that the fluctuating force in the Langevin equation has to vanish within the Markov approximation Eq. (14) when . With a quadratic potential and it is then possible to solve the Langevin equation analytically which gives the time-dependent energy

(16) |

where is the initial energy. However, as shown in Ref. olsen5, , the Langevin equation is quantum mechanically exact for a harmonic potential if the initial conditions are accounted for correctly and the total energy should thus not be allowed to decay below . The problem is that the Markov approximation neglects all non-thermal excitations of the electronic system and leads to pure dissipation at . In reality, an oscillating adsorbate will induce (non-thermal) excitations of the electron gas which may then influence the propagation of the adsorbate. In general, it is therefore expected that the Markov approximation tends to underestimate the influence of the electronic system on the adsorbate. This non-Markovian effect should vanish at high temperature where the thermal excitations of the electronic system dominate. In Fig. 4 we show the time evolution of the average energy of a harmonic oscillator interaction with a thermal reservoir of electrons at six different temperatures. The average energy is calculated using the full non-Markovian correlation function as described in appendix B and within the Markov approximation. The initial state was chosen as the vibrational ground state and included exactly by phase space sampling the Wigner distribution.olsen5 () The parameters used were chosen to match the internal vibrational mode of CO adsorbed on Cu(100)olsen2 (); olsen5 () and we have thus taken , , , and Å. The failure of the Markov approximation and resulting decay of the average energy is clearly seen at low temperatures. In particular, at the Markov approximation gives rise to exponentially decaying energy whereas the energy remains nearly fixed at when memory effects are included. For high temperatures () thermal excitations dominate and the Markov approximation becomes reliable. In all calculations we have converged the results by decreasing the time steps.

In Fig. 5 we show the average energy of the harmonic oscillator after of interaction with a thermal reservoir of electrons. The energy is calculated with non-Markovian Langevin dynamics, Markovian Langevin dynamics, the master equation with rates obtained from perturbation theory and the master equation with exact (non-perturbative) rates.gao97 () Using the full non-Markovian master equation Eqs. (6)-(7) does not change the results. The non-Markovian Langevin approach matches the exact non-perturbative Master equation approach, whereas the perturbative Master equation fails at high temperatures and the Markovian Langevin approach fails at low temperatures. It may be surprising that the non-Markovian Langevin equation reproduces the exact and not the perturbative master equation. However, as mentioned above, the perturbative derivation of the master equation is based on a direct evaluation of the reduced density matrix to second order in the non-adiabatic coupling,gao97 () while the Langevin equation is derived by constructing an effective action to second order in the non-adiabatic coupling.brandbyge () Thus, while the reduced density matrix calculated from the effective action only becomes exact in the small friction limit, it does contain terms to all orders in the non-adiabatic coupling and is a much better approximation for large frictional coupling and high temperatures than the direct perturbative derivation leading to the master equation Eqs. (6)-(7).

It may seem like a complete overkill to apply a non-adiabatic Langevin dynamics to a harmonic potential when the results are readily obtainable from the master equation approach. However, for anharmonic potentials it is not possible to derive transition rates for the master equation exactly and the best approximation is then the non-Markovian Langevin equation. This was also concluded in Ref. olsen5, where a perturbative master equation approach underestimated transition rates in a Morse potential compared to a Markovian Langevin approach.

### iii.2 Associative desorption of N from Ru(0001)

As an example of a potential where the master equation approach is not readily applicable, we consider the well-known example of associative desorption of N from Ru(0001). This system has been subject to extensive experimentaldiekhoner1 (); diekhoner2 (); murphy () and theoreticalmurphy (); diaz1 (); diaz2 (); diaz3 (); luntz1 () studies and much evidence points to a non-adiabatic dissipation of energy during associative desorption.

The Langevin equation can be generalized to an arbitrary potential by a semiclassical expansion of the potential and the excited state forces acting on the adsorbate.brandbyge () The potential is included by making the substitution in the Hamiltonian Eq. (1) and the friction arises from an excited resonant state with potential energy and is included by the generalizations and which in the wide band limit leads to a position dependent resonance width . When multiple coordinates are considered , , and in Eqs. (10)-(13) become tensors and the Langevin equations for each coordinate become coupled through terms like . In addition to temporal correlation, the stochastic forces acting on different coordinates become correlated through the off-diagonal terms in the correlated function:

(17) |

In fact, since the friction tensor is well approximated by with , has only one non-zero eigenvalue. This implies that there is a single (coordinate dependent) mode on which the stochastic force acts and the random forces can thus be regarded as completely spatially correlated at any given time.

We have studied associative desorption of N from Ru(0001) using the code gpaw,gpaw (); mortensen () which is a real-space Density Functional Theory (DFT) code that uses the projector augmented wave method.blochl1 (); blochl2 () The Ru(0001) substrate was modelled by a three layer slab where the top layer was relaxed. We used a grid spacing of Å and the calculations were performed in a (2x4) supercell sampled by a (4x6) grid of -points using the RPBEhammer () exchange-correlation functional. The friction is assumed to be dominated by the orbital which is only partly occupied in the ground state. To calculate the excited state potential energy we applied a generalization of the -self-consistent field method where the resonant state is expanded in a basis of Kohn-Sham orbitals and the resulting resonant density is added to the density in each iteration step. Thus for each adsorbate position we calculate the energy resulting from forcing an electron into a orbital. For details on the method and comparison with experiments we refer to Ref. gavnholt, We have restricted the analysis to the two-dimensional desorption process considered in Ref. murphy, where the two N atoms are adsorbed at adjacent hcp hollow sites and desorbs by moving perpendicular to the bridge towards the fcc hollow while changing the center of mass coordinate. While a two-dimensional analysis is almost certainly not sufficient to obtain quantitative results,luntz1 (); diaz2 (); diaz3 () we do expect to draw some qualitative conclusions about the validity of the Markov approximation for this system. The calculated ground and excited state potential energy surfaces are shown in Fig. 6. To obtain we have fitted the width of the projected density of states of the orbital along the minimum reaction to an exponential and obtained , Å and is the center of mass position at the transition state. The frictional force in the internal mode become large in the exit channel and gives rise to large dissipation of the internal energy while the molecule desorbs. However due to the rapid decay of the friction tensor essentially vanishes at Å. The amount of dissipated energy thus largely depends on the time spend in the exit channel in the immediate vicinity of the transition state.

To examine the impact of non-adiabatic dissipation of energy and, in particular, the validity of the Markov approximation, we have considered four representative initial conditions leading to desorption. All four are initially at the transition state with a kinetic energy of . The kinetic energy is then concentrated in positive or negative center of mass momentum or positive or negative internal momentum. We have taken the surface and thus the electronic temperature to be .murphy (); luntz1 (). Table 1 displays the average energy loss in a desorption event of the four initial conditions with and without the Markov approximation. The reason for the large difference is due to the average time spend in the exit channel which for initial negative internal momentum is and for initial positive center of mass momentum is . The shift to lower dissipation in non-Markovian dynamics is what we would expect from the conclusions in Sec. III.1 and Fig 4. In general, memory effects tends to increase the importance of fluctuating forces and thus decrease the overall dissipation of energy.

The present analysis should in no way be taken as a quantitative study of non-adiabatic effect in associative desorption of N from Ru(0001). In such a study one would need to sample a thermal distribution of initial configurations at the transition state and include all 6 molecular degrees of freedom.diaz3 () Furthermore, the present study assumes that the electronic friction originates from a single resonance () and it is well described by the wide band approximation. That this is not the case has already been establishedluntz1 () and a complete ab-initio scheme for the electronic friction is needed. Such a scheme based on DFT has already been suggested and put to use within the Markov approximationtrail (); luntz1 (); luntz2 (), but need to be generalized slightly to take memory effects into account. In section IV we will propose such a generalization.

Mode | z- | z+ | d- | d+ |
---|---|---|---|---|

Markovian | 0.31 | 0.49 | 0.079 | 0.10 |

non-Markovian | 0.13 | 0.42 | 0.053 | 0.10 |

For anharmonic potentials where the friction tensor acquires a position dependence, non-Markovian Langevin dynamics is rather time consuming, since one has to calculate and diagonalize the correlation function in each time step. In the present simulation, the time required for a single time step in the dynamics was increased by a factor of compared to Markovian dynamics, but the computational time is, however, vanishing compared to that required for a full DFT calculation at a given position. For N on Ru(0001) the memory effects are clearly seen but probably not important compared to neglecting four degrees of freedom. The calculated energy dissipation is not large enough to account for the vibrational damping observed in Ref. diekhoner2, , but it is very likely that inclusion of more degrees of freedom would result in a larger amount of time spend in the exit channel and thus a much larger dissipation of internal energy.

## Iv Non-Markovian friction and fluctuating forces from density functional theory

Within linear response theory, it is possible to derive an expression for the electronic friction for a general non-adiabatic Hamiltonian if one assumes classical adsorbate motion.trail () The result depends on the response function of the electronic system as well as the derivative of the electron-vibron coupling with respect to adsorbate coordinates. Replacing the true response function with a Kohn-Sham response function and the coupling Hamiltonian by a Kohn-Sham potential, give the result for the electronic friction

(18) | ||||

where are Kohn-Sham orbitals with eigenenergies . The result is valid within the Markov approximation, since the memory in the Kohn-Sham potential has been neglected. Generalizing this result to include non-Markovian dynamics would require a non-adiabatic exchange-correlation potential, which is presently out of reach. However, since the result is equivalent to that obtained within the reduced density matrix formalism and the Hamiltonian (1) we can impose a very simple generalization which reduces to the adiabatic result (18) in the Markov approximation and to (10)-(13) in the case of non-interacting electrons. Indeed, it is easy to verify that (18) reduces to the Markovian limit of (10)-(13) if is replaced with the Hamiltonian (1) and one is led to a non-Markovian Langevin equation based on DFT which is given by (8)-(11) and (13), but with

(19) | ||||

This result for the dynamical friction was also derived in Ref. trail, , however, the dominating non-Markovian effect is the correlation function (13) which follows from the reduced density matrix formalism in conjunction with (19).

## V Summary

From a fundamental point of view it is important to realize that the Langevin equation gives an exact description of a harmonic oscillator interacting with a reservoir of electrons if the initial quantum state is taken correctly into account.olsen5 () However, it is easy to see that the fluctuating forces vanish at low temperatures in the Markov approximation, which then results in a exponentially decaying energy of the oscillator. This of course contradicts the quantum description of the oscillator and the problem can be traced to the Markov approximation which does not take non-thermal electronic excitations into account. In Fig. 4, we have shown explicitly how memory effects ’saves’ the quantum behavior of the oscillator and conserves the energy of the vibrational ground state at low electronic temperatures.

Another way of handling dissipative systems is using the master equation. This approach is based on a perturbative calculation of the reduced density matrix in basis of energy eigenstates. While this approach is fast and intuitively appealing, the method breaks down at high temperature or large friction due to the perturbative nature of the method. In contrast, the Langevin equation is based on an effective action Eq. (21)-(22) giving a non-perturbative flavor. Furthermore, the master equation requires quantization of the potential energy surface and becomes impractical for complicated potentials with many bound states.

As an example of such a potential we have considered the associative desorption of N from Ru(0001). We have not tried to perform a quantitative two-dimensional study of this system as done by Luntz et al.luntz1 (), but rather examined the effect of temporal correlation in two representative trajectories. As expected, the effect is an increased significance of the fluctuating forces leading to lower dissipation when memory is included. While this is may be of qualitative interest, the effect of including all degrees of freedom and performing an ab-initio calculation of the friction tensor, would almost certainly lead to corrections which are quantitatively much more important.diaz3 ()

Finally, we have provided an expression for the correlation function within an ab-initio DFT scheme. The result follows naturally by combining the usual DFT based friction tensortrail () with the relationship between the friction and fluctuating forces in a non-adiabatic Newns-Anderson Model.brandbyge () In principle, this scheme allows one to model non-adiabatic dynamics at metal surfaces by Langevin dynamics with ab-initio non-Markovian friction and fluctuating forces.

###### Acknowledgements.

This work was supported by the Danish Center for Scientific Computing. The Center for Individual Nanoparticle Functionality (CINF) is sponsored by the Danish National Research Foundation.## Appendix A Path integral representation of the reduced density matrix

In this appendix we give a path integral representation of the reduced density matrix Eq. (4) in a coordinate basis. We will focus on the probabilistic interpretation of the path integral which leads to a Gaussian distributed classical Langevin equation.

In a coordinate representation the reduced density matrix is

(20) |

and the diagonal elements give the probabilities of finding the adsorbate at a particular position regardless of the state of the electronic system. As shown in Ref. [brandbyge, ] the reduced density matrix of the Hamiltonian (1) can be represented as a double path integral:

(21) |

with the effective action given by

(22) | ||||

where , and and are given by Eqs. (10) and (13) respectively. With a quadratic potential the non-interacting action is given by . Changing to coordinates and and performing a partial integration on the kinetic term then gives for the diagonal part of the density matrix:

(23) |

where is the Wigner distribution of ,

(24) |

and have the additional constraint that . For a non-quadratic potential , one is forced to make a semiclassical expansion of the potential to second order and the exponential in Eq. (23) would contain additional terms of order . Similarly, with a nonlinear interaction in (1) one can perform a semiclassical expansion of the frictional terms which leads to a position dependence in Eq. (12).

Without the quadratic term in , the density matrix (23) would give a delta functional in and the dynamics would be governed by a classical equation of motion with dynamical friction function . However the last term in the exponential of (23) gives rise to a Gaussian broadening of the classical path. To see this explicitly we ”complete” the square in the exponential and perform the path integral in which gives

(25) |

where solves

(26) |

The exponential in (25) can be interpreted as the probability density of taking the path given the endpoints and and the initial velocity . It has a maximum at corresponding to the classical dynamics and the classical path is broadened by . However, it will be more convenient to consider the probability density of which obviously has dimensions of a force. It is then necessary to change the path integral measure from to and it can be shown that the Jacobian of this transformation is independent of .schmid () The two-point correlation function of can then be calculated by

(27) |

This is the most compact way of specifying the statistical properties of and Eq. (24) can be regarded as a classical equation motion with a stochastic Gaussian distributed force .

## Appendix B Discretization of the correlation function

To sample a correlated ”force path” we need to discretize the correlation function and diagonalize the resulting correlation matrix. For a set of Gaussian distributed stochastic variables with probability distribution

(28) |

The correlation matrix can be assumed symmetric without loss of generality. Hence, there exist a diagonal basis of uncorrelated variables which can be sampled from independent normalized Gaussians. The transformation can be obtained by a Cholesky decomposition of such that , where .

The stochastic force appearing in the Langevin equation can be regarded as an infinite number of stochastic variables; one for each point in time from to . Thus, to obtain an expression for the fluctuation force in a time interval , we need the statistical properties of the integrals

(29) |

From the theory of multivariate Gaussian distributions it is readily shown that the set of these integrals are Gaussian distributed with the correlation matrix:

(30) |

and this is the expression used when calculating molecular trajectories using non-Markovian Langevin dynamics.

## References

- (1) C. H. Christensen and J. K. Nørskov, J. Chem. Phys. 128, 182503 (2008).
- (2) A. C. Luntz and M. Persson, J. Chem. Phys. 123, 074704 (2005).
- (3) J. I. Juaristi, M. Alducin, R. D. Muiño, H. F. Busnengo, and A. Salin, Phys. Rev. Lett. 100, 116102 (Mar 2008).
- (4) A. C. Luntz, I. Makkonen, M. Persson, S. Holloway, D. M. Bird, and M. S. Mizielinski, Phys. Rev. Lett. 102, 109601 (Mar 2009).
- (5) J. I. Juaristi, M. Alducin, R. D. Muiño, H. F. Busnengo, and A. Salin, Phys. Rev. Lett. 102, 109602 (Mar 2009).
- (6) L. Diekhöner, H. Mortensen, A. Baurichter, and A. C. Luntz, J. Chem. Phys. 115, 3356 (2001).
- (7) C. Díaz, J. K. Vincent, G. P. Krishnamohan, R. A. Olsen, G. J. Kroes, K. Honkala, and J. K. Nørskov, Phys. Rev. Lett. 96, 096102 (2006).
- (8) C. Díaz, J. K. Vincent, G. P. Krishnamohan, R. A. Olsen, G. J. Kroes, K. Honkala, and J. K. Nørskov, J. Chem. Phys. 125, 114706 (2006).
- (9) M. J. Murphy, J. F. Skelly, A. Hodgson, and B. Hammer, J. Chem. Phys. 110, 6954 (1999).
- (10) L. Diekhöner, L. Horneær, H. Mortensen, E. Jensen, A. Baurichter, V. V. Petrunin, and A. C. Luntz, J. Chem. Phys. 117, 5018 (2002).
- (11) C. Díaz, A. Perrier, and G. J. Kroes, Chem. Phys. Lett. 434, 231 (2007).
- (12) J. Behler, B. Delley, S. Lorenz, K. Reuter, and M. Scheffler, Phys. Rev. Lett. 94, 036104 (Jan 2005).
- (13) J. Behler, K. Reuter, and M. Scheffler, Phys. Rev. B 77, 115421 (Mar 2008).
- (14) H. Metiu and G. Schön, Phys. Rev. Lett. 53, 13 (1984).
- (15) A. Schmid, J. Low Temp. Phys. 49, 609 (1982).
- (16) A. O. Caldeira and A. J. Leggett, Physica A. 121, 587 (1983).
- (17) M. Brandbyge, P. Hedegård, T. F. Heinz, J. A. Misewich, and D. M. Newns, Phys. Rev. B 52, 6042 (1995).
- (18) J. C. Tully, M. Gomez, and M. Head-Gordon, J. Vac. Sci. Technol. A 11, 1914 (1993).
- (19) A. C. Luntz, M. Persson, S. Wagner, C. Frischkorn, and M. Wolf, J. Chem. Phys. 124, 244702 (2006).
- (20) T. Olsen and J. Schiøtz, Submitted , arxiv:1003.2318.
- (21) J. R. Trail, D. M. Bird, M. Persson, and S. Holloway, J. Chem. Phys. 119, 4539 (2003).
- (22) P. W. Anderson, Phys. Rev. 124, 41 (1961).
- (23) D. M. Newns, Phys. Rev. 178, 1123 (1969).
- (24) S. Gao, Phys. Rev. B 55, 1876 (1997).
- (25) T. Olsen, J. Gavnholt, and J. Schiøtz, Phys. Rev. B 79, 035403 (2009).
- (26) T. Olsen and J. Schiøtz, Phys. Rev. Lett. 103, 238301 (2009).
- (27) R. P. Feynman and F. L. Vernon, Ann. Phys. (N.Y.) 24, 118 (1963).
- (28) T. Olsen, Phys. Rev. B 79, 235414 (2009).
- (29) The gpaw code is available as a part of the CAMPOS software: www.camd.dtu.dk/Software.
- (30) J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen, Phys. Rev. B 71, 035109 (2005).
- (31) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
- (32) P. E. Blöchl, C. J. Först, and J. Schimpl, Bull. Mat. Sci. 26, 33 (2003).
- (33) B. Hammer, L. B. Hansen, and J. K. Nørskov, Phys. Rev. B 59, 7413 (1999).
- (34) J. Gavnholt, T. Olsen, M. Engelund, and J. Schiøtz, Phys. Rev. B 78, 075441 (2008).