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 one-dimensional (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 AFM-experiments. (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 PT-model 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.

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 ) d t 2 + m γ d r ( t ) d t + ∇ V ( r , t ) = ζ ( t ) , (1) with r (t) being the time-dependent position of the AFM-tip of mass m and γ 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 tip-drive subsystem, respectively, then V ( r , t ) = U ( r ) + 1 2 k c v s t − r 2 , (2) 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 v s = v s i cos Φ + j sin Φ , (3) 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 ζ ( t ) ζ ( t ′ ) = 2 m γ k B T δ t − t ′ , (4) 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 ζ ( t ) = ζ ( t ) i cos ⁡ θ + j sin ⁡ θ (5) with θ being a uniformly distributed random angle.

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, d ξ ( t ) d t = v ( t ) d v ( t ) d t = 1 m ζ ( t ) − γ v ( t ) − 1 m ∂ V ( ξ , t ) ∂ ξ , (6) which, once numerically solved, yields simultaneously the position ξ(t) and velocity v(t) of the AFM-tip along the sliding path for ∀t ≥ t ₀, if the initial position ξ ₀ and velocity v ₀ of the AFM-tip are both known at the starting time t ₀, i.e., ξ ( t 0 ) = ξ 0 v ( t 0 ) = v 0 . (7)

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 t i = t n + ∑ j = 1 i − 1 a i j Δ t ξ i = ξ n + ∑ j = 1 i − 1 a i j k i , ξ v i = v n + ∑ j = 1 i − 1 a i j k i , v w i = G i q i Δ t W ⇒ k i , ξ = Δ t v i k i , v = − Δ t γ v i + 1 m ∂ V ( ξ i , t i ) ∂ ξ + Δ t w i m , to finally get ξ n + 1 = ξ n + b 1 k 1 , ξ + b 2 k 2 , ξ + b 3 k 3 , ξ + b 4 k 4 , ξ v n + 1 = v n + b 1 k 1 , v + b 2 k 2 , v + b 3 k 3 , v + b 4 k 4 , v + O ( Δ t 4 ) , ∀ n = 0,1,2 , … , where all the involved real-valued coefficients as derived in Ref. (Kasdin, 1995), G _i (i = 1, …, 4) are Gaussian random numbers and W = 2mγ 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 21 a 31 a 32 a 41 a 42 a 43 b 1 b 2 b 3 b 4 = 1 2 0 1 2 0 0 1 1 6 1 3 1 3 1 6 . (8)

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 a 1 = a 1,0 a 2 = a 1 2 , 3 2 , (9) 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 ₁ a ₁ + n ₂ a ₂, with n 1 , n 2 ∈ Z and correspond to a 2D hexagonal Bravais lattice with an angle γ = π/3 rad between a ₁ and a ₂, since   cos γ = a ₁ ⋅ a ₂ = 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, b 1 = 2 π a 1 , − 1 3 b 2 = 2 π a 0 , 2 3 (10) and hence the corresponding translation vector within the 2D reciprocal space reads K _m = m ₁ b ₁ + m ₂ b ₂ ( m 1 , m 2 ∈ Z ).

Now, from the Fourier series expansion of the 2D periodic potential U( r ) = U( r + T _n ), U ( r ) = ∑ K m U K m ⁡ exp ( i K m r ) , 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 R e   U ( r ) ± 10 U ± 10 + R e   U ( r ) 0 ± 1 U 0 ± 1 = cos 2 π a ξ cos ⁡ φ − sin ⁡ φ 3 + cos 2 π a ξ 2 ⁡ sin ⁡ φ 3 , 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 ₀₀. In this manner, one could have any desired sinusoidal corrugation as defined by φ, for example, U ( ξ ) = − U 0 2 cos 2 π a ξ cos ⁡ φ − sin ⁡ φ 3 , (11) providing the temperature independent force acting on the AFM-tip in the form of ∂ V ( ξ , t ) ∂ ξ = π a U 0 cos ⁡ φ − sin ⁡ φ 3 sin 2 π a ξ cos ⁡ φ − sin ⁡ φ 3 − k c v s t − ξ . (12)

Indeed, when φ = 0, the most applied and simplest sinusoidal corrugation of period equal to the lattice constant a results, U ( ξ ) φ = 0 = − U 0 2 cos 2 π a ξ a n d ∂ V ( ξ , t ) ∂ ξ φ = 0 = π a U 0 ⁡ sin 2 π a ξ − k c v s t − ξ . (13)

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 high-symmetry 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.6 eV, and hence correspond to fcc Au(111). However, for the large-parametric study presented in Seciton 6, both U ₀ 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, η = π 2 U 0 k c a 2 2 (14) 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.

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, F i = 1 N t i ∑ j = 1 N t i F t j (15) 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 t i 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⁴. 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 〈 F 〉 n = 1 n ∑ i = 1 n F i (16) either frictional or reaction one (the latter only within the 2D PT-model), and its standard deviation σ n = 1 n ∑ i = 1 n F i − 〈 F 〉 n 2 (17) are both calculated using numerical results from n = 10⁴ 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 〈 F 〉 = lim n → ∞ 〈 F 〉 n σ F = lim n → ∞ σ n . (18)

As can be directly seen from Figure 3A, the distribution of time-averaged friction force within a PT-model is a normal (Gaussian) one, N F | 〈 F 〉 , σ F 2 = 1 2 π σ F 2 exp − F − 〈 F 〉 2 2 σ F 2 , (19) of an almost perfect shape in view of its third (skewness) and fourth (kurtosis) central moments, s k n = 1 σ n 3 1 n ∑ i = 1 n F i − 〈 F 〉 n 3 k u n = 1 σ n 4 1 n ∑ i = 1 n F i − 〈 F 〉 n 4 . (20)

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 time-averaged 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, r = 1 σ F σ q 1 n ∑ i = 1 n F i − 〈 F 〉 q i − 〈 q 〉 , (21) 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 ₀-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 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⁴ 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⁴. 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 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³, 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., s k n ≤ ϵ G & k u n − 3 ≤ ϵ G . (22)

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⁴, 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 〉 1 0 4 taken as reference when n = 10⁴. 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.

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

When U ₀ 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 non-convergence 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 stick-slip 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 ₀ 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 ₀ 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.

Apart from identifying the frictional regime according to the material and system properties, the 1D and 2D PT-models also enable us to study the temperature and velocity dependence of friction and reaction forces in both stick-slip and superlubric frictional regimes. In Figures 9, 10 only the calculated friction force for a large range of temperature and sliding velocity values is given, because in magnitude much smaller reaction force behaves pretty similarly. It can be readily observed that the increase of temperature leads to a monotonic decrease of friction forces within both frictional regimes, as the increasing thermal motion aids the AFM-tip to overcome the interaction energy barriers and slide with a lower mechanical energy. On the other hand, the increase of sliding velocity results in increased friction forces strongly depending on the frictional regime, probably due to different impact of the damping. A linear dependence of the friction and reaction force (although the latter is not shown explicitly here) on the logarithm of the sliding velocity can be observed within the stick-slip regime, a 1 k B T + a 0 ln v s v s ( r e f ) , (23) of which slope a ₁ k _B T + a ₀ is linearly depending on the absolute temperature T. In this Eq. 23 k _B stands for the Boltzmann constant and v s ( r e f ) denotes the selected reference sliding velocity, here and in the following, the minimum of all considered sliding velocities. Within the superlubric frictional regime, a similar functional dependence as in Eq. 23 holds only when U ₀ is varied. When k _c is varied, then the quadratic form a 2 ( k B T ) 2 + a 1 k B T + a 0 ln v s v s ( r e f ) 2 , (24) provides the temperature and sliding velocity dependence of the friction and reaction forces. For the values of the polynomial coefficients a _p (p = 0, 1, 2) and further computational details, see the Supplementary Material. All these are suggesting at least a power-law dependence of friction reaction force on the sliding velocity with a temperature-dependent exponent.

In summary, the PT-model was considered in one and two dimensions for the numerical description of an AFM-tip sliding over a periodic and atomically flat surface. A system of SDEs was considered and implemented for the inclusion of the effects of temperature on the motion of the tip. The conditions for the derivation of statistically accurate friction and reaction forces in a numerically efficient way were defined. The proposed statistically proper scheme can be applied for sliding on both commensurate and incommensurate paths, while it allows the study of the impact of both system and environmental conditions on the resulting friction forces.

By the numerical study of the AFM-tip sliding using the PT-models, it was concluded that 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, can uniquely determine the type of frictional regime, either stick-slip or superlubric between the AFM-tip and any semi-infinite crystalline substrate. In addition, the friction force in both regimes is predicted to be a function of temperature and sliding velocity.

All these findings can support the design of appropriate experimental studies for the realization of structural superlubricity, while they can also shed light into the design of novel systems showing ultra-low friction which can be beneficial in several engineering applications.