Thermal effects and spontaneous frictional relaxation in atomically thin layered materials

We study the thermal effects on the frictional properties of atomically thin sheets. We simulate a simple model based on the Prandtl-Tomlinson model that reproduces the layer dependence of friction and strengthening effects seen in AFM experiments. We investigate sliding at constant speed as well as reversing direction. We also investigate contact aging: the changes that occur to the contact when the sliding stops completely. We compare the numerical results to analytical calculations based on Kramers rates. We find that there is a slower than exponential contact aging that weakens the contact and that we expect will be observable in experiments. We discuss the implications for sliding as well as aging experiments.


Introduction
During the past decade, atomically thin sheets consisting of single or multiple layers of a 2D material have been the subject of much investigation [1,2] and are an important component in nanotechnology in general [3,4]; as well as a promising candidate as friction reducing agents in microscopic as well as macroscopic applications [5]. The frictional properties of atomically thin sheets have been found to exhibit several interesting phenomena. Specifically, friction decreasing with increasing number of layers, and an initial friction strengthening. These were first observed in AFM experiments over a decade ago [6,7]. To further the understanding of these materials, atomistic MD simulations have also been performed by a number of groups, (such as [8,9]). However, significant developments were recently made in this area using a simple model [10].
While much understanding has been gained of the strengthening and layer-number-dependence of friction of layered-materials, one aspect that is yet to be understood or investigated is how friction, in such systems, is impacted by thermal effects. Thermal fluctuations play an important role in experimental nano-scale friction, especially for systems like these, where the energy barriers against sliding tend to be relatively low [11,12]. Conversely, thermal effects in general in layered materials have been the subject of much scrutiny (see, for example [13,14]).
In this work, we analyse the thermal behaviour of the phenomena of strengthening and layer dependence of the friction of atomically thin sheets using the simple modified Prandtl-Tomlinson (PT) model [15][16][17] that has previously been shown to capture and explain these phenomena [10]. We add thermal noise to this model, and investigate its effect on the distortion of the sheets and the friction. Furthermore we investigate how contact age- * Electronic address: astrid.dewijn@ntnu.no ing results from thermally activated relaxations in the system. Contact ageing -when the strength of the contact changes during sticking periods -is an established phenomenon that occurs on wide range of length scales from the atomic (see, for example [18][19][20]) to the geological [21,22]. Ageing usually leads to increasingly strong contacts, and thus plays a role in the phenomenological (macroscopic) observation that static friction is typically larger than dynamic friction. Simple extensions to the PT-model has been used before to investigate the effect of aging on friction, see e.g. [23], the novelty in the present work however is that the model geared towards atomically thin layered materials.
We will begin by introducing and quickly reviewing the model proposed in [10]. Following this, we will discuss how the thermal fluctuations affect the dynamics during sliding. The main body of this report will be devoted to studying thermal relaxations after sliding has stopped: an AFM tip that has been scanning and deforming a sheet for some distance is stopped; then, due to thermally activated decay, the system slowly re-equilibrates as the sheet relaxes. We combine this with analytical calculations of thermal activation using Kramers theory, which we show can calculate the most probable decay paths in the complex potential landscape, and thus predict estimated relaxation times.
An interesting consequence of the structure of the potential landscape is that relaxation times appear over many orders of magnitude, which implies that it should be possible to observe these effects on experimental time scales. It also shows that the contact ageing is slower than exponential, which is typically what is observed in experiments (see e.g. [18]). What is however more unusual is that the ageing actually weakens the contact. tip support substrate sheet Figure 1: Illustration of the system and model setup [10]. The model is similar to the PT model, with a support pulling the tip. A single additional variable is introduced to describe the state of distortion of the entire sheet. Figure from [ The model we consider here was first presented in Ref. [10], and is an extension of the traditional onedimensional PT model. It is based around the setup shown in Fig. 1.
As in the PT model, the tip (position x) is pulled via a spring by a support moving at a constant velocity, v. The tip is subject to viscous damping and the system has a potential energy that describes the interaction between the tip and the material it slides on. To describe the atomically thin sheets between the tip and the substrate, this extended model introduces an extra degree of freedom, q, which describes the distortion in the sheet. Here, distortion means any change in the sheet or sheets; this can be bending, shearing, shifting, or very local displacements of individual atoms. The total potential energy of the system is then described as a function of both the tip and sheet via The first term in the potential describes the stiffness of the tip and AFM cantilever with a single spring constant k. The other terms in the total potential energy are: the internal potential energy of the sheet V sheet ; the interaction between the tip and the sheets V tip−sheet ; and the interaction between the tip and the substrate through the sheet V tip−substrate . The terms in the potential energy related to the sheet and substrate can be written in a general way by including leading order, and in some cases next-to-leading order, terms in an expansion in q [10]. Higher orders were found to give no extra insight or relevance to the problem. This yields where ν 2 and ν 4 are parameters describing the potential energy from distortion of the sheet; V 1 and V 2 are the potential-energy corrugation of the tip on an undistorted sheet and the contribution of the substrate respectively; κ 1 and κ 2 account for changes in the corrugation; and a and b are the lattice parameters of the sheets and substrate respectively. The extra degree of freedom for distortion, q, has both a global and local effect on the potential landscape. The local effect is created by the local interaction with the tip through the cosine term in equation 3; whereas all other instances of q give rise to global effects in the potential landscape of the system. This model captures a number of typical frictional behaviours exhibited by layered materials in AFM experiments. Specifically, it displays the strengthening that appears as an initial increase in the force until the regular steady-state stick-slip pattern begins, as well as the dependence of the friction on the thickness of the sheet. These characteristics of the friction of layered materials have been reproduced many times in both AFM experiments and molecular dynamic simulations [6,8].
A crucial role in the dynamics of this model is played by the position of the energy minima in the potentialenergy landscape. In the low-velocity limit, the system is quasi-static, and the tip is pulled from minimum to minimum, initially distorting the sheet and then later sliding along it. For more details, see Ref. [10].
In Fig. 2 we show a baseline summary of how the model behaves in the absence of thermal noise. It shows the strengthening, and in the xq plot one can see how this is related to the structure of the energy landscape. After the strengthening, during steady-state sliding, the tip remains in a small range of q and stays near one row of minima in the potential-energy landscape. This is explored in some detail in [10], and is reflected in Fig. 2. This results in a periodic steady-state, that looks very similar to the one-dimensional PT model.
In this work, we introduce thermal fluctuations into this model using Langevin dynamics, i.e. in addition to the viscous damping with damping constants η x and η q in x and q respectively, we include Gaussian random forces ξ x (t) and ξ q (t). The equations of motion then read where m x and m q are the effective masses of the tip and sheet distortion respectively. The thermal noise is given superimposed on the periodic terms in the potential-energy landscape. The rest of the parameter values are given in the text. The behaviour of the system is controlled by the minima in the potential-energy landscape. Figure from [10], shared under the Creative Commons Attribution 4.0 International License [24].
by the amplitudes A x,q = 2k B T η x,q and Gaussian uncorrelated noises ξ x (t) and ξ q (t) with standard deviation 1 and 0 mean. Unless otherwise mentioned, the parameter values we use are the same as in our previous work [10]. They were estimated from detailed molecular-dynamics simulations that showed strengthening and layer-dependence [8] and that modelled an AFM tip sliding over layers of graphene. In some cases of missing information, we have used AFM experiments as a starting point and picked our values based on this. The parameter value are as follows. The masses are m x = 10 −23 kg and m q = 3.67 × 10 −24 kg, which are in the order of magnitude thousands of carbon atoms, similar to that used by previous authors in this context [8]. The damping coefficients are η x = 18.75ps −1 and η q = 42.86 ps −1 , which are typical values for an AFM. For the tip-support interaction: v = 1.0 ms −1 , k = 2.0 Nm −1 = 12.5 eVnm −2 and a = 2.5 A; the support velocity is similar to that used in MD simulations, the spring constant is on the order of magnitude of an ordinary cantilever, and the lattice spacing is that of graphene when scanning in the geometrically regular, zig-zag, direction. To keep matters simple, we use a commensurate combination of lattice parameters, a/b = 1. The corrugation parameters are V 1 = 0.31 eV and V 2 = 0.15 eV, which are common energy levels use d in the ordinary PT model. The corrugation coefficients: κ 1 = 0.375 eVnm −2 and κ 2 = 0.1875 eVnm −2 , were picked on the same grounds as the corrugation parameters. The distortion energies are chosen to be ν 2 = 2.39eVnm −2 and ν 4 = 3.64eVnm −4 , these represent binding energies between the sheet and the substrate as well as the internal bonding energies in the sheet respectively, and their values were picked accordingly.
The exact parameter values are not critical for the model to qualitatively reproduce friction dynamics of layered materials. Rather the model works for a wide range of parameter values. We will show however how the V and κ energy parameters change the barrier heights of the potential energy landscape, and how ν parameters restricts the accessibility of the corresponding potential minima. Crucially, ν 4 needs to be large enough to stabilize the q variable, while still allowing some freedom for q to change.

Numerical simulations
The equations of motion [Eqs. (5) and (6)] were integrated using a fourth order Runge-Kutta algorithm implemented in C++ with timestep 15 fs. This timestep size was verified by checking energy conservation, as well as giving a reliable sampling frequency even during slips.
For sliding simulations we choose as initial conditions an undistorted sheet (q = 0) and the tip and support both at x = vt = 0, with tip velocityẋ = 0. Along with studying the trajectories, we investigate the thermal relaxation after sliding has stopped. In this case, as v = 0, we refer to the support position as r in this stationary case. For this purpose, we collect statistics from a large number of simulations. To simplify the analysis, we start each of these simulations from the same initial conditions. We use as initial conditions the zero-temperature slipping point in the steady state sliding regime. This point was used as initial state for the relaxation simulations. This corresponds to the following values: x = 1.03 nm, x = 1.00 ms −1 ,ẍ = 1.86 × 10 13 ms −2 , q = 0.74 nm, q = 0.51 ms −1 ,q = 2.17 × 10 13 ms −2 and r = 2.448 nm. While this is not a procedure that can be repeated in an experiment, as we will discuss later on, the aspects of the decay that would be experimentally detectable are not affected.
We identify decay points using a simple procedure. To eliminate the risk of large, short-range thermal fluctuations being picked up as decays, we first reduce the noise by smoothing it with a running mean over 1000 timesteps (15ns). From the potential energy we can calculate all the minima in the potential landscape. Furthermore, we know by construction that the tip starts in the minimum labeled m1. We can define a proximity measure: of the tip to each minimum with some coordinates (x 0 , q 0 ). The algorithm at every timestep checks if the tip is sufficiently close to any minimum (d < 0.025 nm).
If that is the case, then it checks if this minimum is the same minimum that it was assigned to previously. If that is not the case it further checks if the tip has been closer to this minimum than any other for ten consecutive timesteps. If that is the case, the algorithm determines that a decay has occurred into the new minimum at the time all these conditions were first met.

Simulation results
The thermal fluctuations in this distortion can influence the friction and other behaviour of the system in several ways. In this section, we first investigate numerically the effects during normal sliding, and then move onto the phenomena that appear when the sliding is reversed or stopped.

Sliding
In Fig. 3, we show a typical example of the evolution of a force and q trace, and a path taken in the underlying potential-energy landscape by the tip. We observe that besides the normal irregularity in the force trace, which one would expect due to thermal noise, there is an additional pattern of larger short-time decreases in the force corresponding to temporarily lower values in q: when the tip falls into a minimum on a row closer to q = 0. This corresponds to lower corrugation of the sheet and thus lower lateral force.
The model displays a thermal behaviour that is similar to the one-dimensional PT model [25]. The random thermal noise kicks the tip over the barrier into the next minimum of the substrate earlier than it would have moved without the noise. This leads to slips at a lower lateral force, and thus lowers friction. This phenomenon is called thermolubricity, see [11] for a formal treatment, and [12] for experimental support within the context of a one-dimensional PT system. Since the present model has one extra degree of freedom compared to the PT model, the thermal behaviour is more complex. There superimposed on the periodic terms in the potential-energy landscape. The rest of the parameter values are given in the text. We see an irregular stick-slip pattern due to thermal noise, and the system sometimes randomly changes to a lower value of q.
are, as in the one-dimensional PT model, thermally activated slips in the x-direction. However, the tip can also overcome the barrier separating two rows of minima, and briefly switch to a row of minima with a lower distortion, q. During the next slip, the strengthening process once again pulls it up to its native row. This is responsible for a very typical dip in the lateral force and is a signature of the presence of the sheet that might be observable in experiments.
In the one-dimensional PT model with nonzero temperature, thermal activation of slips can be handled using Kramers rate theory, resulting in a friction that depends on the corrugation of the substrate, the velocity, and temperature through [11] Figure 4: The friction as a function of (a) velocity and (b) temperature, along with theoretical estimates based on a constant value of q and Sang's Kramers rate treatment [11] of the one-dimensional PT model, given in Eqn. (8).
where F c is the lateral force in the inflection point of the substrate potential in the absence of thermal noise, v * is a dimensionless velocity that depends on the support velocity as well as temperature, and ∆F is a constant. Both F c and ∆F involve the corrugation of the substrate. For a formal treatment, as well as full expressions for these constants we refer to [11]. We can estimate the friction based on these onedimensional expressions by reducing our system to onedimension. This is done by approximating q as a constant and obtaining an approximate corrugation from that. Previously, we showed that q in the steady state could be estimated reasonably well from its maximum attainable value, q max , given by the sole real root to the polynomial: 2ν 2 q max + 4ν 4 q 3 max + 2κ 1 q max − 2π a (V 1 + κ 1 q 2 max ) = 0 [10]. From Eqn. (1), we can see that the corrugation potential then has two terms, one from the sheet and one from the substrate, with corrugations U 1 = V 1 + κ 1 q 2 max and U 2 = V 2 + κ 2 q 2 max respectively. The tip slips between minima of both terms, and must cross barriers inbetween, which lie on the ridge of maximum tip-sheet energy. We can estimate the effective corrugation by considering the energy difference between the minima (where both terms are at a minimum) and the points on the ridges where the tip-sheet term has a maximum, but the tip-substate term can have any value between −U 2 and +U 2 . This gives lower and upper bounds for the effective corrugation of Figure 4 shows the friction as a function of velocity and temperature computed from the simulations, as well as the estimated lower and upper bounds for the friction.
In the simulations, the friction was obtained as the average lateral force from 120 ns simulations, excluding the strengthening period. The observed friction is within the estimated limits, but the range between our upper and lower bound is quite wide, and the quantitative behaviour is also somewhat different. This is not surprising given that there is a plethora of new dynamics in the q degree of freedom that is unaccounted for, including, for example, the hopping between different rows of minima that can be seen in Fig. 3.

Reversing direction
We now briefly consider the transient effects when the sliding direction changes. This is common in experiments to measure friction loops, where the direction of sliding is reversed halfway through. A simulation of such a force trace is shown in Fig. 5. The parameter values in this figure are as with sliding. The sliding direction is reversed after 2.5 ns, after which time the support retraces the same path it came. The distortion of the sheet q does not immediately disappear when the tip changes sliding direction, and so signatures of the distortion generated by the first half of the loop are visible during the beginning of the second half in the form of a double strengthening regime. This friction evolution is in good agreement with that found in detailed MD simulations and experiments on layered materials [6][7][8].

Relaxation after stopping
We now turn to the ageing effects that appear after the sliding stops. An example of this is shown in Fig. 6. Once the tip has stopped, the force decays by finite increments, as the distortion also decreases and the sheet relaxes. This effect can again be understood by considering the potential-energy landscape. At this point, the distortion variable, q, is non-vanishing, and thus the system is not in a global energy minimum, but only in a local minimum. The thermal noise provides random kicks to  the tip, allowing it to overcome potential energy barriers by thermal activation, and move into a minimum with a lower potential energy. This decay is a type of contact ageing.
Since we are interested in understanding what signatures of this might be visible in experiments through changes in the lateral force we investigate the characteristics of this process of relaxation in more detail in terms of the path that the tip is likely to take across the landscape and the lifetime of each decay out of each minimum on its journey.
In the standard one-dimensional PT model, the system would relax down to a potential minimum dependent on the position of the support. However, as can be seen from Fig. 6, the added degree of freedom, q, also relaxes. In general, the global minimum does not correspond to a completely relaxed sheet (q = 0), but the sheet remains somewhat distorted. Only for very specific support po-sitions does the global minimum correspond to a fully relaxed sheet. We have chosen a position for stopping the support that is typical, and does not produce any special geometry in the potential-energy landscape.
The precise path that the tip takes through the minima depends on the locations of the minima and barriers between them. While the spring force makes obvious slips towards the lowest minima, the initial values of q remain largely unchanged, as can be seen in Fig. 6. This will be further examined below in subsection 4 .1. Another immediately notable feature of Fig. 6 is that a slip in the force is often accompanied by a quickly reversed slip in q. This will be explained in further detail later in section 4 .2.
100 30 ns simulations of relaxation with the same initial conditions were also generated and examined. Two of the tip trajectories from this data set can be seen in Fig. 7 and will be discussed in section 5 .1.

Theoretical considerations and calculations of the decay rates
We now look more closely at the decay and discuss analytical approaches for better understanding the ageing. The computational simulations of the model can show us what happens at the beginning of the decay, but due to limitations in computation time, we cannot determine numerically what will happen to the system after more than around a µs of waiting, and we cannot explore the entire space of system-specific parameters. We can however obtain insight from the analytical calculations discussed in this section.

Geometry of the potential-energy landscape
The ageing of the system is controlled by thermally activated decay and by the available metastable minima in the potential-energy landscape. We must therefore consider the general features of the geometry of these minima. An example of the total potential-energy landscape for the parameter values we have used in the simulations is shown in Fig. 8, with black contour lines signifying level curves. To understand its structure, we must look at the various contributions to the potential energy separately.
As was discussed in the previous work [10] and can be seen, for example, in Fig. 2, the minima of the tipsheet-substrate system sit more or less in a grid defined by the cosine terms in Eqns. 3 and 4 for V tip−sheet and V tip−substrate . However the V sheet term of the potential energy (Eqn. 1) parabolicly warps the potential landscape around q = 0, while the support term in a similar fashion warps the landscape around the support position (indicated by r). This introduces a global 2d potential well centred on x = r, q = 0. This removes any minima that are too far away from this point and shifts the rest Relaxation measured for 30 ns at 300 K. (b) The corresponding path (grey line) of the same data in relation to minima in the potential energy landscape. Contours are spaced every 0.4 eV to give a better guide of the shape of the landscape.
towards it. Due to the higher corrugation for large absolute values of q, more minima survive further away from q = 0. For each minimum there are at most eight nearby minima. Double decays may be possible. However, overcoming of the initial potential barrier to escape the minimum is the same regardless of how many minima or saddle points the tip goes through during the rest of the decay.
The position of the stopped support in the lattice unit cell impacts the position of the global minimum and structure of the other local minima around it, and thus the behaviour of the system in the final stages of relaxation. The support may be positioned such that there is a simultaneous minimum for the tip-support interaction and the tip-sheet interaction for an undistorted sheet.  This is the case when the support position is almost equal to an integer, n, multiple of the lattice constant, r ≈ na. In this situation, the system has a possibility of nearly relaxing fully to where the tip position aligns with the support's position and the sheets have relaxed to their lowest energy so that the final state is close to x = r, q = 0. Conversely, the support may have stopped close to a maximum in the tip-sheet interaction for an undistorted sheet. In this case, the global minima to which the system may relax will in general be further away from x = r or q = 0, and the tip and sheet will not be able to relax completely simultaneously. The most extreme case of this is the most uncongruous support position, r ≈ (n + 1 2 )a.

Calculations of the relaxation rates using Kramers theory
The relaxation path once the tip has stopped is long and consists of many steps. Since the simulations can only show us what happens in the beginning of the decay, in order to understand better the relaxation and the path that the system takes, we calculate the relaxation rates for each possible step along the path. We investigate the escape time of the tip from each metastable minimum through nearby saddle points during relaxation.
The decay rate for each step is determined by the local structure of the potential-energy landscape and the temperature. Here, we use these to calculate the escape rate from one minimum through a saddle point using Kramers' theory [26][27][28].
The Kramers escape rate is given by: where E m and E s are the Hessian matrices of the potential energy in the (metastable) minimum and saddle point respectively. The energy difference between the minimum and saddle point is given by ∆U , and λ + is the exponential growth rate of a small deviation from the saddle point, and T is the temperature of the system. We note that the mass m x only appears in λ + , and the functional form of this dependence is analytically known. We numerically determine the minima and saddle points by minimising the expression ||(∂V /∂x, ∂V /∂q)|| 2 , and extracting any points that simultaneously meet the condition of both the partial derivatives in x and q are zero, ∂V /∂x = 0, ∂V /∂q = 0. They are indicated in Fig. 8. We then compute the 2 × 2 Hessian matrices. The growth of deviations from the saddle point λ + is obtained as the largest positive eigenvalue of the 4 × 4 dynamic matrix M in the saddle point, where M is defined by d(x, q,ẋ,q)/dt = M · (x, q,ẋ,q). Applying these values to Eqn. (11), we then obtain the Kramers escape rates.
What remains is the question of where the system goes after passing through a saddle point. As can be seen from Fig. 8, in many cases the saddle point is located almost directly between two nearby minima. In those cases there is little question which minimum the system will decay into after passing through the saddle point. However, there are some saddle points where this is less obvious, especially near the higher potential-energy minima. This can be seen in Fig. 8 for example for m3 decaying via s4. For these decays, as well as for some others, for example from m2 through s2, the reaction path is not a straight line. In Fig. 6, this can be seen as a reversed movement in q.
Various minima, saddle points and their energy differences, along with the corresponding Kramers rates calculated for 300 K, are given in table I. The most probable decay path and a sample of the decay times are indicated in Fig. 9.

Comparison of decay times and paths between the analytical calculations and simulations
We now compare the geometric arguments and Kramers rate calculations detailed in the previous section to results from simulations as far as the available computing power allows. It requires less statistics to determine the tendencies in the path of the tip than the relaxation time, so we first focus on this.

Complexities in the potential-energy landscape and consequences for the decay of the tip
Examining 100 30 ns simulations starting from identical initial conditions described above and a fixed support at r = 2.448 nm, we discovered several unsurprising, but nonetheless odd, paths the tip had taken. These differed from the theoretically expected paths outlined in Fig. 9 and demonstrated in Fig. 6. These included 13 of which decayed from one minimum and missed several neighbouring minima before occupying another. One of these interesting paths can be seen in Fig. 7(a) where the tip moves from minimum m1 to minimum m4. Other long trajectories that bypassed minima occurred between m1 and m3, and m2 and m5.
Another interesting decay path that occurred more frequently is a trajectory between m2 and m4 where the tip appears to occupy a minimum that does not exist, as seen in Fig. 7(b). 34 of the 100 paths examined had a trajectory similar to the one in the figure between these two minima. Even though there is no true minimum present here, the potential landscape is sufficiently shaped to allow the tip to maintain its locality for some time before finding its next minimum. This is, of course, a symptom of the geometry of the potential landscape and the support's position, as discussed in section 4 .1. In comparison, 33 of the paths had similar trajectories between m2 and m3 as seen in Fig. 6(b).
Worth noting in Fig. 9, that may not be explicit in the image, is the predicted path is indicating repeated movements between some minima, indicated by two arrows. This prediction of a jiggle between two minima arrives from the calculated life times. For example, minimum m8's shortest life time, as seen in Fig. 9 (and table I), is towards m10, while minimum m10's shortest life time is towards m8. This motion would repeat until the tip travels in the next most probable direction, i.e. either m11 from m8, or m11 from m10. This behaviour is due to the structure of the landscape and should appear regardless of the precise parameters. It is therefore a signature of the sheet dynamics that could be observed in experiments.

Lifetimes
To obtain more detailed information about the lifetimes, 10000 120 ns simulations were performed in the same manner. The decay detection algorithm described in section 2 was used to extract the time stamps and modes of decay. In figure 10 we show the frequency of decay as a function of the length of time that the tip remained in a specific minimum. The decay is exponential, as expected from any process with short-time correlations.
We fit exponential functions to a 9 point moving average (9PMA) and extract the decay rates and life times. The results of this fitting are indicated in the figure. The observed and calculated lifetimes are very close. We can use the results for the Kramers decay rates for rest of the minima, that we cannot reach in the simulations, to understand what happens as the tip moves closer to the support position: the decay times increase. As can be seen from the more complete table of calculated life  times (table I), the decay times can become quite long, in the order of s. These long times are not achievable in the simulations, and therefore we cannot observe the full decay in simulations.

Implications for ageing in experiments
Signatures of this ageing and these decay rates might be visible in experiments. While the specific decay rates and lifetimes are linked to the precise potential-energy landscape, there are general qualitative behaviours that appear regardless of the quantitative energies.
The number of decays is much larger and the path is more complex than without the presence of the sheet. For our parameters, if the sheet is removed, there would be a maximum of 3 minima before the tip reaches the support position. With the sheet included there are an order of magnitude more energy minima at high q and positions further away from the support. This complexity gives rise to such things as the hopping back and forth that we predict will happen (for example m11 and m14 in Fig. 9). This behaviour is linked to the geometry, and does not qualitatively depend on the precise parameters. It is, therefore, likely to also appear in experiments.
Another signature of these decays is the broad range of time scales involved, from picoseconds to milliseconds (see table I). This is a general property, and is a result of the structure of the energy landscape and the wide range of energy barriers. This means that in experiments, where longer time scales are accessible (as opposed to in our simulations) these slower decays might be observ- able. It would take the shape of a slower than exponential contact ageing. The ageing is different from the type of contact ageing found in more complex contacts [19], in that it weakens the contact and decreases the restarting friction when sliding begins again, rather than increasing it. Although a higher number of dimensions may be truer to the complex system expected with contact ageing, it appears as though two dimensions offers enough complexity to gain insight into some of the behaviour that may appear in contact ageing.

Conclusions
We have investigated thermal effects in friction on atomically thin sheets in a simple model based on the Prandtl-Tomlinson model that reproduces the strength-ening effect seen in AFM experiments. We show that at finite temperatures there is thermolubricity, as in the normal Prandtl-Tomlinson model. We also identify a number of thermal effects that could be observed in further experiments.
During sliding, the tip can slip to a lower distortion as well as in the sliding direction. This can be observed as a larger than normal change in the force during slipping followed by a repeat strengthening, and a very typical shape of the force trace.
After sliding stops, the contact ages slowly due to the complex potential-energy landscape combined with thermal decay. Due to the presence of the sheet, there are many decays with different barriers and decay rates that would not have appeared in a plain PT model. We have performed analytical calculations to estimate the decay times that are not accessible in simulations. The initial decays are fast, but as the tip approaches the position of the support, the decay slows down. The calculated and observed lifetimes show that there is a wide range of decay times, from several picoseconds to seconds, giving hope to the possibility that the phenomenon of relaxation can be observed in AFM experiments. In addition, there are other signatures that might be observable in experiments, such as the hopping back and forth between sets of minima on the same slanted line, which is a general feature. We also note that while there is no way to directly measure q in an experiment, it is possible to obtain it indirectly by simply starting to slide again and observing the amount of strengthening that occurs.
The wide range of lifetimes and the sequence of the decays raise an interesting possibility to understand contact ageing in a general way. While this model was designed for friction on thin sheets, it is so general that it could also be used to investigate other systems with extra degrees of freedom. In realistic systems, many different degrees of freedom inside an asperity could play a similar role; including, but not limited to elastic deformation of the material and adsorbed layers from the atmosphere. The model shows that even one such degree of freedom can give rise to a sequence of decay times of many different orders of magnitude, leading to a nontrivial, non-exponential ageing of the contact.  Table I: Energy barriers and lifetimes for decay from various minima via specific saddle points as estimated using Kramers' escape rate theory. The temperature is T = 300 K. The lifetimes range from picoseconds to a second. The minima, saddle points, and their labels are indicated in Fig. 8. We do not include any decay times longer than one day. The values used in figure 9 are in bold in this table.