On the Origin of Frictional Energy Dissipation

The origin of the friction between sliding bodies establishes an outstanding scientific problem. In this article, we demonstrate that the energy loss in each microscopic slip event between the bodies readily follows from the dephasing of phonons that are generated in the slip process. The dephasing mechanism directly links the typical timescales of the lattice vibrations with those of the experienced energy ‘dissipation’ and manifests itself as if the slip-induced motion were close to critically damped.


Introduction
Friction continues to fascinate engineers and scientists in spite of its seeming simplicity. We are taught that sliding friction involves the conversion of mechanical energy into heat [1,2] and that this conversion necessarily is irreversible [3]. The microscopic picture that is usually associated with this conversion involves phonons [4][5][6] that are excited within the sliding bodies, typically through mechanical instabilities, such as stick-slip events. Due to phonon coupling and the associated, finite phonon lifetimes [7][8][9][10], the energy that specific phonon modes acquire from these instabilities is thought to be quickly redistributed over the full spectrum of possible phonons.
The coupling between a specific degree of freedom, e.g., associated with the motion of a sliding body, and the other degrees of freedom of a system, such as the amplitudes and phases of all vibrational eigenmodes of the body itself and those of the solid over which the body is forced to move, can be treated theoretically as the coupling to a bath of harmonic oscillators. This formulation introduces the combination of all other degrees of freedom as a thermal or harmonic bath. From this description, one readily derives the generalized Langevin equation and the fluctuation-dissipation theorem that relates the strength of thermal fluctuations to the dissipation rate [11,12].
In this article, we concentrate on the bath of harmonic oscillators, i.e., on the full spectrum of vibrational modes in either of the two bodies that are sliding over each other, and ask the question which modes are really coupled directly to a specific slip event and how they manage to carry away the energy invested in them in such an event so efficiently that it becomes effectively irretrievable for re-use in later slip events. The tempting answer to the latter question seems that again coupling between phonons, i.e., finite phonon lifetimes, would be at play, which would re-establish thermal equilibrium after each slip event. Instead, we will demonstrate that there is a more natural explanation.
Finite phonon lifetimes derive from the anharmonicity of the interaction between atoms in solids and from the presence of impurities and other structural defects. Whereas anharmonicity is an intrinsic material property, defect densities and impurity concentrations often vary over orders of magnitude, usually without a strong influence on friction. This should be taken as a serious indication that the internal redistribution of mechanical energy is not primarily due to the thermalization of phonons. Energy line widths measured by inelastic neutron scattering on bulk phonons [7][8][9][10] and by inelastic helium atom or electron scattering on surface phonons [13] show that in spite of all anharmonicities and structural non-idealities, these mechanical eigenmodes 'live' for tens of vibrational periods or more. In order to redistribute the energy rapidly enough to cause significant friction, however, phonons would have to be close to critically damped, i.e., convert their energy into other phonons on a timescale of one or just a few vibrational periods. Such strong damping is predicted theoretically only for the vibrations of isolated adsorbate atoms on a semi-infinite solid [2], but is not characteristic for the surface or interior of the solid itself. Near-critical damping [14][15][16] also forms an essential ingredient in the interpretation of friction force microscopy (FFM) images, in which atomic stick-slip patterns are routinely observed [15][16][17], with occasional slips of the FFM tip over two or more lattice spacings [17][18][19]. Again, details of the materials and their structural perfection, both for the tip and the substrate, seem not to be of critical importance.
In the following, we present an alternative non-thermodynamic view on the redistribution of energy within a sliding body, accompanying a slip event. We treat the dynamics in terms of the dephasing of phonons that are excited in the slip event. While being consistent with the formulation of damping due to the coupling to a thermal bath, our description leaves out phonon coupling altogether. In the calculations presented in this article, near-critically damped motion nevertheless emerges and we show that this is a natural consequence of our description of slip events as the simultaneous excitation of a large number of vibrational eigenmodes. This result invites us to speculate about new methods to modify friction, simply by manipulating the spectrum of available, vibrational eigenmodes.

Model System: Harmonium
In order to reduce frictional contact dynamics to its very essence, we performed a combined theoretical and numerical study of the simplest possible model system in which friction might arise, namely that of two initially static slabs of an idealized, completely regular solid material, made of atoms that interact with each other through short-ranged, exclusively harmonic forces. By construction, our system contains neither anharmonicities, nor defects, impurities, or adsorbate layers. This renders the lifetime of all phonons infinite. Our calculations nevertheless indicate that slip events are followed by behavior that is best described as very rapid, near-critical damping. We show that this damping originates from the progressive destructive interference of the phonons that are excited in the slabs by the slip events.
We refer to our harmonic model substance as harmonium, (Hr). The Hr is organized in a body-centered cubic (bcc) lattice, with a dimensionless cubic lattice constant of unity. In order to ensure stability of the bcc lattice, harmonic forces are introduced between nearest-and next-nearest-neighbor atom pairs with dimensionless spring coefficients of 2 and 1 respectively (see Appendix 3). We imagine mechanical contact to be established between the (001) surfaces of two identical, infinite slabs of this material, touching each other only via two individual Hr atoms, one on each of the two surfaces. In the calculation, each slab is represented by a periodically repeated rectangular block, containing N atoms. We keep the situation completely symmetric between the two contacting slabs, so that we can concentrate fully on one of the two. One stick-slip cycle then corresponds to a sequence with a stick-phase, in which the slab is first deformed via an external force exerted on the contacting Hr atom at the surface, followed by a slip event in which the external force is suddenly reduced to zero. We follow the resulting motion of all Hr atoms in the slab, paying specific attention to the characteristic time scales of their response. Note, that our geometry also represents the situation in which an FFM tip is responsible for the initial surface deformation in the slab. In that case, our calculations describe the motion in the slab, following the slip event of the tip. Figure 1a shows the configuration in the harmonium slab, prior to the slip event. In this example, a (dimensionless) lateral force of 0.1 (see Appendix 3) is exerted along the [100] surface direction on the central surface atom, which is displaced as a consequence. An accompanying deformation pattern is also visible within the slab that decays with distance to that atom. The configuration shown in Fig. 1a is the equilibrium displacement pattern for the specific lateral force exerted on the central surface atom. This is a stationary pattern; the atoms are all standing still.

Numerical Results
At t = 0, the external force is removed-the 'slip' eventand we numerically evaluate the equations of motion of all atoms in the system. Our first observation is that the deformation pattern rapidly accelerates back, as is illustrated by the movie (available online). The solid black curve in Fig. 1b shows the displacement of the central surface atom along the [100] direction as a function of time. We see that the atom overshoots the zero position and goes through a rapidly damped oscillation. The occurrence of a small number of zero-crossings shows that the motion is slightly underdamped.
As strong as the damping may seem in Fig. 1b, we have performed our calculation without any explicit damping on the motion of any of the atoms. In fact, the motion for each atom was obtained by simply integrating Newton's equations along the x-, y-, and z-directions, i.e., along [100], [010] and [001], in response to the forces exerted on the atom by its direct and next-nearest neighbors. No velocity-related terms entered this description [14]. Also, we have left out all other 'hidden' forms of damping, for example via a thermostat in the calculation or via absorbing boundary conditions [20,21] (see Appendix 3 for more details). We observe nearcritical damping, even in the complete absence of an explicit damping mechanism.
At this point we have to acknowledge that there is one mild, implicit form of anharmonicity, which we cannot avoid in our calculations. It reflects the simple fact that a displacement of an atom along one direction changes the directions of most of the nearest-and next-nearest-neighbor bonds that the atom is involved in. This leads to higher-order contributions to the forces on the atoms that make the response of the lattice deviate increasingly from perfectly harmonic with increasing amplitudes of displacement. By repeating our calculations for various values of the initial displacement amplitude (or, equivalently, for various values of the initial external force), we could easily verify that this higher-order effect is not causing the damping observed in Fig. 1b.

Lattice Dynamics Results
We now turn to the dashed blue curve in Fig. 1b, which indicates the result of an alternative computation of the response of the distorted lattice to the removal of the external force at t = 0. Underlying this curve is a calculation of the complete set of 3N − 6 ≈ 3N phonons, the mechanical eigenmodes of our harmonium slab. The result of the phonon calculation is represented by the dispersion curves in Fig. 1c and in Fig. 3 in Appendix 3. They show the angular frequencies of the lattice vibrations k x , k y as a function of the parallel wave vector k x , k y along the three symmetry directions of the reciprocal surface unit cell. The curves display the typical combination of bulk phonon bands and a small number of surface modes, associated with the two free surfaces of the slab. This calculation is completely harmonic; also the subtle anharmonicity is absent that was present in the solid black curve in Fig. 1b. For each vibrational mode k x , k y , our calculation also provides the polarization vector that contains the relative amplitudes, directions and phases with which all atoms in the slab participate in that vibration. Since the eigenmodes form a complete set, each configuration of the slab can be expressed as a unique combination of amplitudes and phases of the 3N phonons. If we perform this projection operation on the initial configuration of Fig. 1a, we obtain a complete picture of the phonons that are excited by the external, lateral force on the central atom. This is indicated in Fig. 1c, which displays the same dispersion curves as Fig. 3, but also shows their relative amplitudes. From the perspective of the phonons, the only change at t = 0, when the external force is removed, is that they all commence their periodic time evolution. This fully harmonic time evolution is nearly indistinguishable from the result of the integrated Newton equations of motion, as is illustrated by the close match between the dashed blue curve in Fig. 1b and the solid black curve. Again, we observe near-critical damping of the motion of the central surface atom, this time for a rigorously harmonic system.
We should stress that the outcome of our calculations does not depend on the strength of the springs or the value of the atomic mass. When these parameters are changed by arbitrary factors, the horizontal and vertical axes of Fig. 1b and the vertical axis of Fig. 1c are rescaled, but apart from this, both figures remain completely unchanged and the near-critically damped character of the motion is not affected.

Phonon Dephasing
There is a natural reason for the observed damping, which should be regarded as an inherent property of solids. The spatial localization of the initial deformation pattern, i.e., the mere fact that the frictional contact is local, necessarily makes that the pattern is composed of contributions from a variety of phonons with different wavelengths. Each individual phonon is a collective vibrational mode of the entire slab and is fully delocalized over the slab. It is only by combining a large number of phonons that a localized deformation can emerge at all. This local concentration of displacement occurs exclusively at the place and time where the phonons are largely in phase, so that their displacements add up constructively. Everywhere else the phase relation between the phonons is 'random' or sufficiently close to random that no significant displacement results. This special combination of place and time is that of the central surface atom that is exposed to the external force and the time origin t = 0, at which the displacement is at its maximum and slip starts. At this point all phonons start to evolve in time, each one with its own angular frequency k x , k y . Because these frequencies are all different, the phonons rapidly run out of phase with each other, thus quickly reducing the displacement amplitude of the central atom. The observed damping is the direct result of the dephasing of the excited phonons with respect to each other. If all frequencies between 0 and the maximum phonon frequency max were represented equally strongly in this process, we should expect the central atom to oscillate effectively with a frequency in the order of 1∕2 max and a dephasing rate in the same order of magnitude, which would render the oscillation critically damped. Figure 1c shows that the excited phonons are not The red and black data are for two different threshold levels of the kinetic energy, namely the value of 10 −14 (red), illustrated in panel (a), and a significantly higher value of 10 −10 (black). The differences are minor. The blue line is a linear fit to the red data points for times above t = 2. It shows that the wave front expands radially at a constant speed of 1.2 distributed completely evenly over all available frequencies, which explains why the observed oscillation is somewhat underdamped.
Since the phonons are all completely delocalized over the entire system, one might expect that the potential energy that is invested by the external force would also be delocalized. This is certainly not the case. It is the region directly around the central surface atom that carries the largest distortions and hence the highest potential energy per deformed bond. After t = 0, it is the same region where the atoms develop the highest velocities, i.e., the highest kinetic energies. When the phase matching is lost progressively for the central surface atom, a wave front travels outwards. Figure 2a shows snapshots of the surface of this front, obtained from our calculations as the outer contour of atoms with kinetic energies above a certain, low threshold value, in this case 10 −14 (see Appendix 4). At each point in time, the maximum of the outgoing wave resides somewhat inside this shell and corresponds to the surface at which the phases match best. The wave travels out at a group velocity that is, like the dephasing rate of the central surface atom, determined by the frequency differences between the excited phonons and, in addition, also by their wavelength differences. One should expect polarization effects to be visible in the form of anisotropy of the wave velocity, for example with the wave running out faster along the [100] axis of the initial surface displacement, due to the stronger longitudinal character of the wave in that direction, and slower along the two perpendicular directions, [010] and [001], due to the more transverse character along those. However, the size of our periodically repeated block was too modest to measure this anisotropy sufficiently accurately. Therefore, we have only determined the orientational average of the velocity of the outgoing wave, as is illustrated in Fig. 2b. We find that the front expands radially with a constant velocity of 1.2. This velocity is in the order of the average value of ∕|k| for the excited phonons, which should be regarded as the appropriate, effective speed of sound for this wave.

Larger Contact Sizes
The single-atom contact that we have considered so far forms an extreme case that can be addressed experimentally only through special instruments such as an atomic force microscope or a friction force microscope. Importantly, our results carry a much more generic character and are also relevant for larger contacts. This is illustrated in Fig. 4 (see Appendix 3), in which the central 3 × 3 surface atoms were subjected to an initial collective displacement along the [100] direction. Qualitatively, the time evolution of the displacement of the central surface atom is rather similar to the single-atom case (Fig. 1a), but the motion proceeds a factor ~ 2 more slowly. Also the dephasing is slower, by approximately the same factor, so that the motion is again somewhat underdamped. Figure 4 shows that the larger contact is associated with a clear selection of the wavelengths of the phonons that are excited; these are concentrated near the Γ-point and correspond to wavelengths of approximately 3 lattice spacings and larger. Again, the phonon-based time evolution is nearly identical to the result of the direct integration of the equations of motion.
The change from the single-atom contact to the 3 × 3-atom contact illustrates an element of inherent scaling that we expect to hold even up to typical tribological contacts with micrometer-size asperities and larger. The size of the contact is a measure for the spatial scale, both along and perpendicular to the surface, of the elastic deformation patterns that result from forces on that contact. The motion induced by the slip event of a macroscopic contact can therefore be viewed as a coarse-grained version of the response that we have followed for the single-atom contact and the 3 × 3-atom contact. The coarse-graining involves the effective volume and mass of the regions in the solid that are set into relative motion, the effective spring coefficients that describe their interactions with each other and the resulting wavelengths and frequencies of the vibrational eigenmodes with which they predominantly move. Even though each of these quantities scales in its own way with the contact size, the qualitative feature remains unchanged, that the slip motion is composed of a superposition of eigenmodes; they result in a damped oscillation with a frequency that is some average over these eigenfrequencies and with a dephasing rate that is determined by the typical difference between these eigenfrequencies. As the 3 × 3-atom contact illustrates, the effective frequency is lower for larger contacts and the dephasing rate is lower in the same manner. As the effective frequency and the dephasing rate are intimately related, their ratio cannot change much with contact size, which renders the slip motion close to critically damped for all contact sizes.

Summary and Discussion
The main conclusion from this work is that the dephasing of excited phonons forms a natural 'recipe' for damping. The essentially new and non-trivial element is that this mechanism occurs in purely linear systems, even though its consequences may seem similar to those of the dynamic stochastization, well-known for nonlinear systems [22,23]. Our observation has been made here in the context of friction, but it applies to all cases where the wave packet of phonons that is excited in a process contains a sufficiently large number of phonons with different frequencies. Other examples of surface phenomena that should be expected to obey similar 'rules' are surface diffusion, the adsorption of atoms and molecules on surfaces and surface chemical reactions [24]. Damping of motion involved in phenomena inside three-dimensional materials, for example bulk diffusion or internal structural changes, or the deposition of energy and momentum by an impinging ion, should behave similarly. In the examples presented in this article the resulting dephasing rate is close to the maximum of critical damping that is possible in this way. The dephasing takes place on the timescale of a small number of vibrational periods, i.e., well before the finite lifetime of the excited phonons would become noticeable. Due to this dephasing, the energy and momentum that are invested in the initial displacements are irretrievably 'lost' on a rapid timescale that is fully decoupled from the slow thermalization of the excess energy by conversion of the excited phonon wave packet in the appropriate thermal (Bose-Einstein) distribution of phonons. Interestingly, the thermalization depends on a multitude of subtle properties of the solids involved, such as the anharmonicity of the interatomic potentials and the character and density of defects and impurities that can act as scattering centers for the phonons. The inherent nature of the phonon dephasing makes the resulting damping mechanism quite robust with respect to these subtleties and therefore very similar even for widely different materials.
The dephasing mechanism invites us to speculate about approaches to modify friction, for example via the contact geometry. One possibility lies in the dimensionality of the materials, such as for graphene and other layered materials with the strongly two-dimensional nature of their phononic eigenstates. Another possibility is offered by the prospect of nanostructuring contacting surfaces into geometries that strongly confine those phonons that are excited during the stick-slip process. Both elements are explicitly present in the spectacular reduction in friction, recently observed by Wada et al. [25].
Acknowledgements This work has been carried out at ARCNL, a public-private partnership of UvA, VU, NWO and ASML, with financial support from the ERC Advanced Grant Project "Science F(r)iction" and from FOM-program Fundamental Aspects of Friction.

The Friction Force Microscope
The instrument that seems most suitable to search for the microscopic origin of friction is the so-called friction force microscope (FFM). In an FFM experiment, a sharp tip is forced to slide over a counter-surface. The normal force, with which the tip is pressed against the surface, is held constant at a pre-selected, adjustable value. The tip is made to perform a two-dimensional scanning motion, while the lateral forces are recorded, thus building up a two-dimensional lateral force map. In such experiments, one routinely obtains lateral force patterns with atomic periodicities [15-17, 26, 27]. Usually, these patterns display a saw-tooth character that is referred to as 'stick-slip' motion; the interaction with the countersurface holds the tip in a local (atomic-scale) potential energy minimum until the lateral force is high enough to make it slip into the next minimum. The mere fact that this process repeats itself from one potential energy minimum to the next seems to indicate that the energy that is released during the slip event is quickly dissipated, so that the tip has to start again at or near the bottom of the next potential energy minimum. What is seen on the atomic-scale in the FFM also happens on longer length scales in macroscopic friction geometries, as we are familiar with from the squeaking of hinges and the stick-slip motion in earthquake dynamics.

Appendix 2 Dissipation
A popular model that captures the essence of this stick-slip behavior has been introduced independently by Tomlinson and Prandtl in the nineteen twenties [28,29]. The Tomlinson-Prandtl model assumes that the excess mechanical energy that is released during the slip phase in the stick-slip motion is removed from the system 'instantaneously'. Interestingly, this assumption works well on all length scales, which makes this a popular model also for the description of FFM experiments [16,26].
A more general theoretical approach to motion with dissipation is given by the familiar Langevin equation, Here, the motion is described by the time-dependent coordinate U(t), is the position-dependent potential energy and is the damping coefficient. Thermal fluctuations contribute a random force R(t) . In the simple form, given here, the equation contains a damping term that corresponds to a force that points against the direction of motion and is proportional to the instantaneous velocity. The Langevin equation is successful in describing a wide variety of phenomena in which damping plays a role, such as atomic-scale friction [16,17], surface chemical reactions [30] and surface diffusion [31], provided that the damping rate is made sufficiently high. In fact, the damping rate has to be in the same order of magnitude as the relevant vibrational frequencies in the system, which makes the system close to critically damped [16,17]. This rapid damping is consistent with the crude assumption of instantaneous dissipation in the Tomlinson-Prandtl model, which should be read as 'fast enough to look like instantaneous'. Various attempts have been made to provide a thermodynamic and atomistic foundation of the semi-phenomenological Langevin equation, leading to a Generalized Langevin equation with memory [32][33][34][35]. To our knowledge, this first-principle-based approach has not reached applications yet.
We have indicated that critical/nearly critical damping is needed in the interpretation of FFM experiments, but remark that this seems to also set the stage for friction on other length scales. On each length scale, the dissipation exhibits a rate that is consistent with the dephasing rate of the vibrational modes that are relevant for that length scale.

General Configuration and Setup
We have performed the calculations in this article for supercells of a body-centered-cubic (bcc) crystal of a hypothetical material, 'harmonium (Hr)', that we constructed from periodically repeated copies of a simple cubic unit cell with a twoatom basis. The results shown here were for supercells with dimensions of 20 × 20 simple cubic unit cells along the x-and y-directions, i.e., the [100] and [010]-directions, parallel to the (001) surface, and a thickness of 15 of these unit cells along the z-direction ([001] direction), corresponding to a total of N = 12.000 Hr atoms. We have conducted our calculations also for other sizes, in order to verify that the reported results were not affected by the finite size of the 20 × 20 × 15 unit cell system. In order to give the supercell the character of an infinite slab, we applied periodic boundary conditions along the two surface directions x and y, while no periodicity was imposed along the z-direction; the top and bottom were free surfaces. The dimensionless lattice constant was set to a = 1.0. The dimensionless mass of the Hr atom was chosen to be m = 1.0 . The atomic interactions were modeled as harmonic bonds between nearest neighbor and next-nearest neighbor pairs of Hr atoms. As dimensionless spring coefficients for these bonds, we chose C 1 = 2.0 and C 2 = 1.0, respectively; this combination ensured the stability of the bcc structure under externally applied, linear and shear forces.
Identical initial configurations were chosen for two types of calculation of the dynamics (see below). We have used a variety of initial conditions. Those that were illustrated by Figs. 1 and 2 were obtained for a situation in which an initial (dimensionless) force of 0.1 parallel to the [100] surface axis was applied to a single surface atom. Below, we also present results for calculations in which the same force was equally distributed over a group of 3 × 3 atoms at the surface. In both cases, an equally large but oppositely oriented force was distributed over all atoms of the bottom surface, to keep the total external force zero and thus avoid acceleration of the supercell as a whole. We started by computing the distorted equilibrium configuration under the influence of these forces, which we obtained by minimizing the potential energy stored in all bonds in the system. This was achieved by numerically integrating Newton's equations of motion for all N Hr atoms in the harmonic system and applying an appropriate (high) damping force on each atom, proportional and opposed to the atom's velocity: In this equation, we define U j (t) as the displacement vector of atom j with respect to the atom's equilibrium position at time t. The atom's mass m carries no subscript, as all masses are chosen equal in our calculation. F j is the initial, external force on atom j , which is zero for most atoms and 0.1 along [100] for the central surface atom (or distributed over the central group of 3 × 3 atoms, depending on the specific starting conditions) and − 0.1 along [100] distributed over all atoms in the bottom surface. C jj ′ is the force constant of the interaction between atoms j and j ′ , which is either C 1 or C 2 , depending on whether the atoms are nearest or nextnearest neighbors (see above) and ΔR jj � is the change in distance between atoms j and j ′ from the equilibrium distance, due to their displacements U j and U j ′ . The sum runs over all nearest and next-nearest neighbors j ′ of atom j . The elastic forces are all oriented along the bonds, hence the bond unit vectors r jj ′ in Eq. 2. The last term in Eq. 2 is the damping term, characterized by the damping constant ′ , which we chose equal to 0.1 in order to obtain rapid convergence. The numerical integration was performed using the Verlet algorithm [36] with discrete time interval dt 1 = 0.002 , and the calculation was continued until the total (dimensionless) kinetic energy of the entire supercell was less than 10 −12 . The resulting, distorted configuration was regarded as sufficiently converged to serve as the static starting configuration for the two subsequent dynamic calculations.

Slip Event
We used the distorted configurations, obtained by the procedure described above, as the starting configuration for the slip event. The slip event was introduced by the sudden removal of the external forces, both on the atom (or atoms) in the top surface and on the entire bottom layer of the supercell. We define the time of this sudden removal of the forces as t = 0 and the initial, distorted configuration, discussed in the previous section, is then given by the set of U j (t) | | |t=0 . The dynamics evolving as a consequence of the removal of the forces were followed by computing U j (t) for times t > 0 in two types of calculations in the complete absence of any damping. The first type is a direct numerical integration, while the second is completely analytical.

Calculation 1: Numerical Integration of Equations of Motion
The first type of calculation consisted of the integration of Newton's equations of motion for all atoms in the supercell, starting with the distorted, initial configuration at t = 0. To this end, we numerically integrated the set of equations Note that this set of equations is identical to that of Eq. 2, but without the damping terms and without the external force(s).
Using the Verlet algorithm, the accelerations, velocities and coordinates were updated with a time interval dt 2 = 2 × 10 −4 . Snapshots were generated every 100 time steps.

Calculation 2: Analytical Solution Via Lattice Dynamics
The coupled set of Eq. 3 can also be solved analytically. Due to the harmonic nature of the interactions and the periodicity of the lattice along the x-and y-directions, the free motion of the Hr lattice is composed of plane waves, with wave vectors in the xy-plane. They are the phonons of the Hr system and, here, we have ignored their quantum-mechanical nature. Together, these 3N − 6 ≈ 3N (in this case 36.000) so-called normal modes form a complete set, which means that linear combinations of these modes can be used to describe the complete vibrational response of the N-atom Hr supercell to any combination of initial positions and velocities of the atoms in the cell. Here, the initial conditions were taken according to the distorted configuration described earlier, characteristic for the situation at the start of the slip event, and with all initial velocities equal to zero. The time-dependent displacement of atom j then takes the following form.
Here, the label k runs over all 3N normal modes and each mode is characterized by a combination of a two-dimensional wave vector q k and an angular frequency k . The polarization vector kj has the x-, y-and z-components of the (complex) amplitude, i.e., the amplitude and phase, with which atom j participates in mode k when the mode has a normalized amplitude of unity, while A k denotes the complex amplitude with which mode k is excited by the initial conditions (see below). The equilibrium position of atom j is denoted by r 0 j . The normalization factor N ′ is the number of wavevectors that enter the calculation as a result of the discrete Fourier transform over the supercell. In our calculation, a three-dimensional periodic supercell with N atoms, would result in N � = N∕2 , due to the two-atom basis of the cubic unit cell. For two-dimensional periodic boundary conditions, along the x-and y-axes, the wavevectors are restricted to two dimensions, q = q x , q y , and N � = N qx N qy , where N qx and N qy are the numbers of allowed wavevectors along the x-and y-directions of the reciprocal supercell, corresponding to the numbers of cubic unit cells along these directions of a large, real-space supercell.
We started by finding the 3N normal modes, for which we used the traditional slab method [37][38][39] for surface phonons. Due to the geometry of our supercell, there were 20 × 20 different choices for the wave vector q k and for each of them 3 × 2 × 15 solutions, i.e., combinations of k and combined polarization vectors k = k1 , k2 , … , kN , the 3 deriving from the dimensionality, the 2 from the two-atom basis of the bcc unit cell of the Hr lattice and the 15 from the thickness of our slab. All k and k were obtained as solutions of the well-known eigenvalue equation In this equation, D q k is the dynamical matrix for wave vector q k [37,38], defined as the (N × N) combination of (3 × 3) submatrices of the type 1 Here, j and j ′ denote two of the N atoms in the supercell. C jj ′ is the (3 × 3) matrix of force constants between the x-, y-, and z-displacements of these atoms, which is obtained readily from the (scalar) force constant C jj ′ , introduced above, and the x-, y-and z-components of the vector r 0 j − r 0 j � , con- necting the equilibrium positions of atoms j and j ′ . For each of the 20 × 20 wave vectors q k , Eq. 5 is solved by finding the eigenvalues 2 k of the matrix D q k and, for each eigenvalue, the corresponding eigenvector k . The result of this calculation is shown as the phonon dispersion relations 2 k q k in Fig. 3.
Once the 3N normal modes were known, we calculated the complex amplitude A k with which each of them was excited/populated by the initial conditions of the static, distorted starting configuration. To this end, we projected the set of initial displacements U j (t) | | |t=0 on the displacement patterns associated with each of the modes, according to where * kj is the complex conjugate of kj . The total energy contributed by mode k is thus given by Note that the specific choice of initial conditions, discussed so far, made the system start exclusively with potential energy, As all initial velocities were zero (to within 10 −12 ), so was the initial kinetic energy in each of the modes: Finally, we have used Eq. 4 with the set of complex amplitudes A k to compute the displacements of all N atoms in the supercell at each point in time, t > 0, after the slip event.
The result was compared in Fig. 1 and in Fig. 4 with the direct numerical integration of the equations of motion. As explained in Sect. 3.1, the modest differences between the calculation results derive from the small degree of anharmonicity inherent in the integration results. The phonon-based analytical calculation should be considered as truly harmonic.
In cases where the initial conditions were not static, i.e., for which the initial velocities were non-zero (second case in next section), the projection operation of Eq. 7 was complemented with an equivalent velocity-based contribution. Of course, in such cases, the kinetic energies of the modes (Eq. 10) were non-zero.

Dependence on Initial Conditions
As explained before, the examples shown in Figs. 1 and 2 were computed for two specific sets of starting conditions. In the first, an external, lateral force was exerted on the central atom in the top free surface. In the second case, an external, lateral force was shared by a group of 3 × 3 atoms in the top surface. In both cases, an equally large, but oppositely oriented force was evenly distributed over the entire bottom surface of the supercell. In both cases the system was allowed to relax completely under the influence of these forces, which resulted in static starting configurations, characteristic of the situation just prior to or at the very beginning of a single slip event in a stick-slip sliding sequence. The results for the 3 × 3 atom case are shown in Figs. 4 and 5. Even though the damping that emerged for the singleatom and 3 × 3 atom cases was qualitatively the same, one might wonder to what extent these results depend critically on the initial conditions. This is why we explored two alternative sets of starting conditions. In the first alternative configuration, we simply displaced the central atom in the top surface laterally, along the [100] azimuth by 4% of a lattice spacing with respect to its equilibrium position. None of the other atoms were given the freedom (or time) to relax their positions in response to this local distortion.
In the second alternative situation, we kept all atoms precisely in their equilibrium positions, but provided the central atom in the top surface with an initial velocity.
While the original starting configurations can be viewed as characteristic for the start of a slip event in case the sliding velocity is so low that the system can relax completely between subsequent slip events, the first alternative configuration captures some of the non-equilibrium character that should be expected at higher sliding velocities. The second alternative configuration can be viewed as more appropriate for situations in which a sudden transfer of momentum occurs, as could be the case when an atom, molecule or ion hits a surface or when kinetic energy is released in a chemical reaction at a surface. Figure 6 compares the time dependence of the displacement along the [100] direction of the central atom in the top surface for the two alternative starting configurations with the original one (the one with the initial force exerted only on a single surface atom). These results were obtained with the numerical integration method. What is evident from the comparison is that the damping behavior exhibited in these three strongly different cases is very similar. Even though the difference in the degree of localization of the initial distortion seems very large between the relaxed (original) and non-relaxed (first alternative) starting configurations, significant amplitudes are already required for the normal modes with the shortest wavelengths (at the X-and M-points in Fig. 1) to properly describe the relaxed case. As a result, the variation in the contributing normal mode frequencies is very similar for the relaxed and non-relaxed cases. This makes the average frequency of the excited modes, which dictates the effective frequency with which the central atom moves, and the dephasing rate, which determines the damping rate of the central atom's motion, for the two cases also very similar. The qualitative difference between these two results and the second alternative, with an initial velocity for the central atom in the top surface, is a time shift by a quarter of the period of the central atom's damped oscillation. The frequency and damping rate are very similar to those for the first two cases.
This comparison between three extreme cases for the initial conditions, (i) locally distorted, fully relaxed and static, (ii) locally distorted, completely un-relaxed and static, and (iii) completely undistorted and locally dynamic, demonstrates that the emerging oscillation frequency and damping rate are rather robust with respect to these characteristics. The only feature of the starting conditions that has a direct influence is the shortest length scale of the initial excitation pattern, irrespective of whether this pattern is defined by displacements or velocities. For the three cases that we compared in Fig. 6, this length scale is equal to that of a single atomic spacing. After all, we dictated either the displacement or the velocity of a single atom. It is only when we change this, as we illustrate in Fig. 4, that we change the emerging oscillation period and damping rate. By increasing the region of influence from a single surface atom to a 3 × 3-atom area, we increased the minimum length scale in the excitation pattern by a factor 3, as is illustrated convincingly by panel (c) of Fig. 4, which shows significant amplitudes only for normal modes with wavelengths of 3 unit cell sizes and longer (wave numbers of 1/3 and less of the Brillouin zone size). When we compare panels (b) of Figs. 1 and 4, we recognize that this large change in excitation pattern has reduced the emerging oscillation and damping rates by approximately a factor 2, while Fig. 5 shows that the 'wave front' moves out at the same speed as in the single-atom case (cf. Fig. 2).
Finally, we illustrate that the damping that we report in Fig. 1 and in Figs. 4 and 5 is not a hidden artifact of our calculations. We use our analytical phonon-based calculation for this purpose. The left panel of Fig. 7 repeats the Γ-X part of the phonon dispersion curves that were shown already in Fig. 3. As before, we use Eq. 4 to compute the time-dependent displacement along the [100] direction of the central atom in the top surface. In this case, we restrict the amplitudes of the excited normal modes to just a single one. In other words, only for a single mode k, the A k -value is made non-zero. In the phonon dispersion plot of Fig. 7 (left panel), we picked a 'random' mode, which is indicated by the blue dot. As Eq. 4 prescribes, the resulting motion is a single standing wave for the entire lattice, in which each atom performs a straightforward sinusoidal oscillation in time without any damping. The right panel of Fig. 7 shows that this is indeed the case. It is the destructive interference between simultaneously active modes with different frequencies that is responsible for the emergent damping, reported in this article and illustrated in Figs. 1, 4 and 5.