Statistics of Sliding on Periodic and Atomically Flat Surfaces

Among the so-called analytical models of friction, the most popular and widely used one, the Prandtl-Tomlinson model in one and two dimensions is considered here to numerically describe the sliding of the tip within an atomic force microscope over a periodic and atomically flat surface. Because in these PT-models, the Newtonian equations of motion for the AFM-tip are Langevin-type coupled stochastic differential equations the resulting friction and reaction forces must be statistically correctly determined and interpreted. For this, it is firstly shown that the friction and reaction forces as averages of the time-resolved ones over the sliding part, are normally (Gaussian) distributed. Then based on this, an efficient numerical scheme is developed and implemented to accurately estimate the means and standard deviations of friction and reaction forces without performing too many repetitions for the same sliding experiments. The used corrugation potential is the simplest one obtained from the Fourier series expansion of the two-dimensional (2D) periodic potential, e.g., for an fcc(111) surface, which permits sliding on both commensurate and incommensurate paths. In this manner, it is proven that the PT-models predict both frictional regimes, namely the structural superlubricity and stick-slip along (in)commensurate sliding paths, if the ratio of mean corrugation and elastic energies is properly set.


INTRODUCTION
Frictional experiments in an atomic force microscope can be straightforwardly simulated by using either molecular dynamics (MD) or ab-initio (first-principles) methods, see for example, Refs. (Kobayashi et al., 2016;Wolloch et al., 2015). From a computational point of view, however, these numerical calculations are extremely demanding and time consuming. Therefore, the so-called onedimensional (1D) and two-dimensional (2D) analytical models of friction, such as the Prandtl-Tomlinson and Frenkel-Kontorova models, (Dong et al., 2011) are still remaining to be highly efficient alternatives to the atomistic models and hence very popular because of their simplicity with which they are capturing the main undergoing mechanisms within nanotribological AFMexperiments. (Pawlak et al., 2016) The Prandtl-Tomlinson model (PT-model), in its original formulation, (Prandtl, 1928) was introduced by Prandtl while developing a kinetic theory of solids assuming a periodic interaction between the bodies. Curiously, in the other paper ascribed to the PT-model, (Tomlinson, 1929) Tomlinson was not dealing at all with these concepts, since he was outlining a molecular theory of friction based on the findings by Lennard-Jones. Despite this historical inaccuracy, (Popov and Gray, 2012) within a 1D or 2D PT-model, one can clearly distinguish between various frictional regimes by considering a single parameter, namely the ratio between the average interaction and elastic energies. (Socoliuc et al., 2004) Traditionally, if the coefficient-of-friction (CoF) observed for a tribological system is less than or equal to 0.01 that system is said to evolve within a superlubric frictional regime. (Martin and Erdemir, 2018) Superlubricity was predicted theoretically in 1990 by Hirano and Shinjo, (Hirano and Shinjo, 1990) and its occurrence was shown experimentally 1 year later. (Hirano et al., 1991) Nowadays, beyond this structural superlubricity discovered in the nineties as caused by the incommensurability of interacting surfaces, one speaks also about thermolubricity (Krylov et al., 2005) and chemolubricity (Erdemir and Martin, 2007) which are two other forms of superlubricity caused by the temperature acting as lubricant and due to chemistry, respectively. For further information on superlubricity, the reader is recommended to consult the excellent reviews in the newest book edited by J.-M. Martin and A. Erdemir. (Erdemir et al., 2021) In this contribution, the structural superlubricity will be addressed numerically by using the most popular and hence widely applied analytical model of friction, namely the PTmodel in one-dimension and two-dimensions, in its well-known form from the kinetic theory of solids and solid mechanics, see Section 2. For this, the periodic corrugation of interest is explicitly derived in Section 4 and directly introduced into the corresponding 1D/2D PT-model. Since the resulting Newtonian equations of motion are forming a set of stochastic differential equations (SDEs), they are solved numerically by applying an adequate fourth-order Runge-Kutta method as presented in Section 3. Due to the randomness of the thermally induced force mimicking the impact of temperature on the sliding, both frictional regimes, either stick-slip or superlubric, must be statistically evaluated and interpreted, for example, in terms of distributions as shown in Secion 5. Beyond the development of a statistically proper numerical scheme for comparing the calculated frictional performance along various paths on periodic atomic surfaces, the ultimate goal of this contribution is to elucidate in which circumstances, i.e., range of various parameters such as materials, elastic coupling of AFM-tip to the drive, sliding velocity or absolute temperature, one could detect the superlubric frictional regime with high accuracy, see Section 6.

PRANDTL-TOMLINSON MODELS
In an atomic force microscope, the movement of the tip on a corrugated surface of a substrate in dry contact conditions is described within the kinetic theory of solids by the Langevin equation of motion, (Filippov et al., 2010) m d 2 r(t) dt 2 + mc with r(t) being the time-dependent position of the AFM-tip of mass m and c is characterizing the damping proportional to the velocity of the AFM-tip. Here the first force on the left-hand side is that which accelerates the AFM-tip, the second one acts against the former and it is proportional to the velocity of the AFM-tip, whereas the last, third force on the left-hand side is due to the interaction of the AFM-tip with the surface and with the drive. If the total potential energy consists of the time independent corrugation potential U(r) and the elastic potential of the tipdrive subsystem, respectively, then where k c stands for the effective stiffness of the linear coupling (e.g., via spring) between the AFM-tip and drive. Note that the drive is supposed to move with a constant sliding velocity v s , such that with Φ defining the orientation of the sliding path with respect to the Cartesian x-axis, whereas i and j being the in-plane unity vectors along the Cartesian axes. The relative direction of v s also sets the commensurability/incommensurability of the sliding path versus the periodicity of corrugation U(r), for more details see Secion 4 below. For example, for a 2D corrugation one can consider the same functional form as for the 1D, such that the latter 1D corrugation provides a periodic atomic chain according to Φ, whereas in the 2D case it represents the simplest 2D periodic potential corresponding to the crystalline orientation of substrate's surface. In Eq. 1, a normal (Gaussian) distribution is assumed for the amplitude of the thermal-noise-induced random force ζ(t), thus for its ensemble averaged auto-correlation function holds that where δ denotes the Dirac delta-function, k B is the Boltzmann constant and T stands for the absolute temperature, and its random orientation is given by with θ being a uniformly distributed random angle.

RUNGE-KUTTA METHODS
In the following, some details on the implementation are given, which can be directly used for a 1D PT-model as well as for a 2D PT-model where ξ(t) is measured along the sliding path. In this case, the second-order stochastic differential equation in Eq. 1 is equivalent to a system of two coupled first-order differential equations, which, once numerically solved, yields simultaneously the position ξ(t) and velocity v(t) of the AFM-tip along the sliding path for ∀t ≥ t 0 , if the initial position ξ 0 and velocity v 0 of the AFM-tip are both known at the starting time t 0 , i.e., The Runge-Kutta (RK) methods for solving differential equations numerically were independently developed by Runge (Runge, 1895) and Kutta (Kutta, 1901) more than a century ago. During their long history, (Butcher and Wanner, 1996) they were continuously improved, (Kalogiratou et al., 2014) such that also RK-methods exist for solving SDEs, one of which will be here applied as follows. The fourth-order Runge-Kutta method applied for the coupled SDEs in Eq. 6, i.e., when T > 0, means that one knowing ξ n ξ(t n ) and v n v(t n ) for t n n Δt (∀n 0, 1, 2, . . .), computes for i 1, 2, 3, 4 where all the involved real-valued coefficients as derived in Ref. (Kasdin, 1995), G i (i 1, . . ., 4) are Gaussian random numbers and W 2mc k B T in accordance with Eq. 5. Due to this latter randomness at finite absolute temperature, Eq. 6 admits not only a single exact solution, but a set of differently valued solutions which, however, are all statistically equivalent down to T 1 mK-as we found. On the other hand, in accordance with Eq. 5, the thermal-noise-induced random force ζ(t) is vanishing at absolute zero temperature, since W 0 when T 0 K, and hence the first-order SDE in Eq. 6 turns at T 0 K into an ordinary differential equation (ODE) and therefore the 'classical' fourth-order Runge-Kutta method for ODEs can be immediately applied. Formally, this means that one considers the former expressions with w i 0 (i 1, . . ., 4) and the following tableaux of Runge-Kutta coefficients, (Press et al., 2007) a (8)

CORRUGATION POTENTIAL AND FRICTIONAL REGIMES
In the following, the corrugation potential of fcc(111)-the conceptual model considered in this contribution-is derived from the Fourier series expansion of the 2D translational invariant potential. For this, consider the 2D primitive lattice vectors for an fcc(111) crystalline surface as being with a a fcc / 2 √ denoting the 2D lattice constant, where a fcc is the lattice constant of the 3D fcc crystal. Accordingly, the 2D translation vectors are given by T n n 1 a 1 + n 2 a 2 , with n 1 , n 2 ∈ Z and correspond to a 2D hexagonal Bravais lattice with an angle c π/3 rad between a 1 and a 2 , since cos c a 1 ·a 2 1/2. The 2D reciprocal primitive vectors b j (j 1, 2), on the other hand, are directly obtained from Laue's equations, a i ·b j 2πδ ij (i, j 1, 2), where δ ij is the Kronecker symbol, such that by using Eq. 9, and hence the corresponding translation vector within the 2D reciprocal space reads K m m 1 b 1 + m 2 b 2 (m 1 , m 2 ∈ Z). Now, from the Fourier series expansion of the 2D periodic potential U (r) U (r + T n ), written for r x i + y j ξ( i cos φ + j sin φ) and by also using Eq. 10, the real part of the "first" Fourier expansion terms can be immediately considered as a proper model of the interaction potential between the AFM-tip and fcc (111) surface, e.g., by assuming for sake of simplicity that U ±10 U 0±1 U 00 . In this manner, one could have any desired sinusoidal corrugation as defined by φ, for example, providing the temperature independent force acting on the AFMtip in the form of Indeed, when φ 0, the most applied and simplest sinusoidal corrugation of period equal to the lattice constant a results, Considering in these latter expressions ξ Γ cos Φ, via Φ one can easily select any commensurate/incommensurate path for sliding, by letting Φ either to show or not to point along a highsymmetry direction on fcc(111), for example, when Φ 0 or Φ π/12 15 • .
The values for all the parameters entering Eqs 1-4, 13 are from Ref. (Dong et al., 2011), e.g., a 2.88 Å, and U 0 0.6 eV, and hence correspond to fcc Au(111). However, for the largeparametric study presented in Seciton 6, both U 0 and k c were continuously varied to determine the impact of the substrate's material and that of the coupling strength on the friction. When k c is kept unaltered, then its constant value is set to 1 N/m. As already mentioned, it turned out that independently of the sliding path, whether it is commensurate or incommensurate, the realization of a normal stick-slip as illustrated in Figure 1, or a superlubric frictional regime provided by Figure 2 indeed is only depending on a single quantity, namely on the ratio between the mean interaction energy of the AFM-tip with the surface and the elastic energy of spring coupling the AFM-tip to the drive, (Socoliuc et al., 2004) see also Section 6.

STATISTICS OF FORCES
From Eq. 13 it becomes immediately clear that the average of the lateral force F(t) k c [v s t − ξ(t)] taken over the sliding time or distance, for example, is a proper measure for the friction force and completely in accordance with the experimental practice. (Schwarz and Hölscher, 2016;Peng et al., 2020) Note, however, that accordingly the calculated friction force could be either positive or negative valued, depending on whether the drive is more times in advance with respect to the AFM-tip or the other way around during the entire sliding period. In this Eq. 15, N (i) t denotes the number of time moments considered during the ith repetition of the same sliding experiments, normally, a couple of thousands time moments up to 10 4 . Here and in the following for producing reference data, at a non-vanishing absolute temperature between 1 mK and the room-temperature (RT) of 300 K, for example, the mean of a force either frictional or reaction one (the latter only within the 2D PTmodel), and its standard deviation are both calculated using numerical results from n 10 4 times repeated "identical" sliding experiments, since it is assumed that after such an amount of repetitions 〈F〉 n and σ n are fairly approximating their exact values As can be directly seen from Figure 3A, the distribution of time-averaged friction force within a PT-model is a normal (Gaussian) one, of an almost perfect shape in view of its third (skewness) and fourth (kurtosis) central moments, Furthermore, by inspecting the normal probability plot given in Figure 3B which per definition identifies graphically the departure of random data from normality, i.e., from the Gaussian behavior, it can be observed that the ordered timeaveraged friction force shows against the theoretical quantiles q i (i 1, . . ., n − 1) indeed a perfect linear correlation in terms of the Pearson's coefficient, where 〈q〉 and σ q denote the mean and standard deviation of q i s. Note that theoretical quantiles are grid-points dividing the range of a probability density function, in our case the Gaussian one, into adjacent intervals of equal probabilities, thus the q i s are values of the inverse cumulative Gaussian density function at i/ n. This normal distribution of friction forces persists in all frictional regimes, i.e., in both stick-slip and superlubric regimes, for any material pairing (arbitrary U 0 -value), and any strength of the elastic coupling (arbitrary k c -value) between the AFM-tip and drive as detailed in the Supplementary Material. Even more, also the reaction force aiming the AFM-tip perpendicularly to the friction force to follow the straight sliding path is perfectly normally (Gaussian) distributed, recall Figure 4. Probably, all these statistical features are less surprising for somebody knowing that averages of random numbers of an arbitrary distribution are always FIGURE 2 | As in Figure 1, for an η-value leading to a superlubric sliding of the AFM-tip on an fcc(111) substrate. normally distributed. Accordingly, as further detailed in the Supplementary Material when applying the 1D PT-model for sake of simplicity, independently of the frictional regime, also the block-averaged time-dependent friction force is normally (Gaussian) distributed. Formally, this is still a direct consequence of the averaging and not necessarily of the fact that the time-dependent friction force is normally distributed too due to the Gaussian-shape of the thermal-noise-induced random force in Eq. 4.
Since our main goal is to numerically study the occurrence of structural superlubricity as function of various parameters, such as temperature and velocity, for example, a 10 4 times repetition of the sliding for a single parametric configuration is too demanding computationally and therefore one should develop a numerical scheme to achieve accurate approximations of F and σ F , recall Eq. 18, by performing much fewer repetitions n ≪ 10 4 . In a first attempt, it was taken advantage on the central limit theorem stating that the composite variable z n n √ (〈F〉 n − 〈F〉) / σ F , with 〈F〉 n as introduced in Eq. 16, tends to be normally distributed with zero mean and unity standard deviation, i.e., N (z n | 0, 1) in accordance with Eq. 19, when the number of repetitions n goes to infinity. This in terms of distribution means that the probability P(z n ≤ z sup ) tends to the cumulative Gaussian density function at z sup , i.e., the supremum of z n , when n → ∞. Indeed, by setting the accuracy for fulfilling the latter identity, one could have reasonable 〈F〉 n and σ n , recall Eq. 17, in comparison with 〈F〉 and σ F for n in range of thousands as presented in Supplementary Material. Unfortunately, for a large-parametric study of structural superlubricity, thousands of repetitions per configuration are still too high.
In a next attempt, it was iteratively attempted to achieve the convergence of the mean force. 1) In a first step, one increases step-by-step the number of repetitions n until the mean 〈F〉 n becomes zero or positive-valued. 2) Then the repetitions are further increased step-by-step until the mean 〈F〉 n+m will be again zero or positive. 3) If the difference of these means |〈F〉 n − 〈F〉 n+m | is below a given positive threshold ϵ F , then it is assumed that the convergence is achieved, and one considers (〈F〉 n + 〈F〉 n+m ) / 2 as being the accurate mean. 4) Otherwise, the FIGURE 4 | As in Figure 3, for the reaction force. FIGURE 3 | In the top panel, the histogram of the time-averaged friction force (blue) for the same sliding experiment repeated 10 4 times together with the corresponding normal (Gaussian) distribution (solid red line) and the mean (dashed red line). Here the skewness (sk) and kurtosis (ku), i.e., the third and fourth central moments, are measuring the high quality of the normal distribution for the time-averaged friction force. This is additionally also quantified by Pearson's linear correlation coefficient r 2 in the bottom panel, where the probability plot of the time-averaged friction force is given.
Frontiers in Mechanical Engineering | www.frontiersin.org October 2021 | Volume 7 | Article 742684 6 algorithm is continued from 1) by setting n to n + m. Unfortunately, in this iterative scheme too the number of repetitions n remains to be of order of 10 3 , see depicted in Supplementary Material.
Finally, it was decided to realize the convergence of the mean friction and/or reaction force by setting the accuracy ϵ G > 0 with which the normality of forces should be achieved. Namely, increasing either step-by-step or blockwise the number of repetitions, it is continuously verified whether 〈F〉 n ≥ 0 and simultaneously the skewness together with the kurtosis, recall Eq. 20, are accurate enough, i.e., If these criteria are satisfied, the mean of friction and reaction forces are considered to be converged and the normality of their distributions is additionally quantified by Pearson's linear correlation coefficient, recall Eq. 21. For the case that one of the criteria in Eq. 22 and/or 〈F〉 n ≥ 0 could not be fulfilled, a predefined maximum number of repetitions n max , e.g., 10 4 , will be performed by assuming that this n max is large enough for obtaining good mean values for the forces. What one observes in Figure 5 is that by decreasing ϵ G for the skewness and kurtosis, the normality of friction force is increasing as quantified by these moments and also by Pearson's (linear) correlation coefficient. However, after a while a further decease of ϵ G does not really improve the positive mean 〈F〉 n in comparison with 〈F〉 10 4 taken as reference when n 10 4 . From Figure 6 it becomes evident that a relatively large ϵ G 0.2 is recommended for using within the stick-slip regime leading to a couple of tenths of repetitions and hence to a fast convergence of the mean values for forces, whereas a ten times smaller ϵ G 0.02 is suitable to achieve highly accurate means of forces with the superlubric regime by performing a couple of hundreds of repetitions.

LARGE-PARAMETRIC STUDIES
Having fixed the convergence criteria with ϵ G 0.02, large parametric studies of sliding along (in)commensurate paths were performed first using the 1D PT-model. By varying the system and environment properties, we aim to gain understanding on the dependence of the friction forces and thus the frictional regime on any of these parameters. As illustrated in Figure 7, the temperature values of T 0 K and T 300 K were considered when solving the corresponding ODEs and SDEs respectively. In these calculations, the η-dependence of structural superlubricity was studied by keeping all parameters entering Eq. 1 constant, except the corrugation potential amplitude U 0 and the spring constant of the coupling of the AFM-tip to the drive k c , respectively, which were varied independently. A different behaviour is observed in the case of each variation.
FIGURE 5 | Convergence of the friction force skewness (blue) and kurtosis (green) within the 1D PT-model, together with the number n of repetitions needed (red dashed line) to achieve a given accuracy ε G , and to this accuracy corresponding Pearson's linear correlation coefficient r 2 n . In addition, the mean of the friction force 〈 F 〉 n is compared with the reference mean 〈 F 〉 10 4 obtained by repeating 10 4 times one and the same numerical AFM-experiment.
Frontiers in Mechanical Engineering | www.frontiersin.org October 2021 | Volume 7 | Article 742684 When U 0 is kept constant, one can observe a short region of very low friction forces, followed by a quick increase with a diminishing slope as the ratio η increases. The introduction of thermal noise at non-zero temperatures leads to the reduction of the friction forces in the stick-slip regime, however, the region of the superlubric regime is not significantly affected. At very low η values, which correspond to very high spring constants k c , negative friction force values of low magnitude are observed. This can be due to the nonconvergence of the average force, as it was already observed that a larger number of statistically independent realisations is typically required in the superlubric regime. On the other hand, it could pinpoint to an incompatibility of the very high elastic coupling energy for the given interaction energy of the AFM-tip with the surface, where the assumptions of the PT-model may not be valid.
In the opposite case where the spring constant k c is kept fixed and the interaction potential is varied, only positive friction forces are observed. The superlubric regime holds for larger values of η, while the increase of forces is less abrupt and shows a linear trend. Here, the introduction of thermal noise also leads to a reduction of the friction forces, however, it also leads to a significant expansion of the range of η values which lead to a superlubric regime. It becomes clear, therefore, from the analysis of Figure 7 that the value of η uniquely determines the type of frictional regime, either stickslip or superlubric between the AFM-tip and any semi-infinite crystalline substrate.
Note that all these numerical studies evidence significantly larger intervals for η where the superlubric frictional regime could occur than that of η < 1 assumed in the literature for both U 0 and k c kept constant, recall for example Ref. (Socoliuc et al., 2004). This means that by varying the material and system properties, the occurrence of superlubricity could be further tuned and even extended on (in)commensurate sliding paths. Furthermore, it seems also that the superlubricity for η < 1, recall Figure 7, it is more pronounced than over its extension for 1 ≤ η < 20, probably because the former is energetically and the latter symmetrically more favorable. Even small negative values for friction force are obtained when η < 1 and U 0 is varied, which shows that probably a sliding distance of only 3.5 nm which was considered here, is not enough to get in average the drive in front of the AFM-tip.
All the same holds for the friction forces calculated using the 2D PT-model, cf. Figure 7 with the left panels of Figure 8. What concerns the η-dependence of the reaction force, this differs not too much from that of the friction force, as one observes on the right column of Figure 8, and shows that keeping the sliding on an incommensurate path needs in magnitude higher reaction forces, when η increases. Interestingly, the incommensurate paths are not equivalent, since the η-dependence of the friction and reaction forces is Φ-dependent too. Furthermore, both forces, i.e., the frictional in both 1D and 2D PT-models, and the reaction force within the 2D PT-model, feature a strong temperature dependence for a given sliding velocity. Figure 5, for the friction (top row) and reaction force (bottom row) within a stick-slip (left panels) and a superlubric frictional regime (right panels) obtained by applying the 2D PT-model.

FIGURE 6 | As in
Frontiers in Mechanical Engineering | www.frontiersin.org October 2021 | Volume 7 | Article 742684 8 FIGURE 7 | η-dependence of structural superlubricity as predicted by the 1D PT-model at 0 K (solid line) and at RT of 300 K (full circles) along various incommensurate sliding paths Φ 5, 10, 15, 20, 25 and 30 deg (depicted in violet as specified in the legend), recall Figures 1, 2, when either U 0 (top row) or k c (bottom row) is kept constant in Eq. 14. Figure 7, for three two-dimensional incommensurate sliding paths Φ 5, 10 and 15 deg, when using the 2D PT-model. Frontiers in Mechanical Engineering | www.frontiersin.org

FIGURE 8 | As in
October 2021 | Volume 7 | Article 742684 9 FIGURE 9 | Temperature and sliding velocity dependence of the friction force as predicted by the 2D PT-model within the stick-slip regime along the incommensurate path for Φ 15 deg, when either U 0 (top row) or k c (bottom row) is kept constant in Eq. 14.
FIGURE 10 | As in Figure 9, but within the superlubric frictional regime.