The effect of temperature and velocity on superlubricity

Recent developments in surface scanning probes have stimulated a large interest in understanding and using the phenomenon of friction at the nanoscale [1, 2]. One of the most fascinating aspects of this field is the possibility of superlubricity, namely the possibility of sliding two surfaces onto each other with vanishingly small friction. Several experiments report such a low friction state [3–5]. The crucial ingredient of this phenomenon is the geometrical incommensurability of the sliding crystal faces. Other aspects like roughness of the sliding surfaces [6] and electronic excitations during sliding [7] are not discussed here. The key model for studying incommensurate systems is the Frenkel–Kontorova model (FK-model) [8] (shown in the top panel of figure 1), which describes a chain of atoms with period a on a potential profile with period b. The case where the ratio b/a is not a rational number identifies an incommensurate system with no finite global periodicity.

Zoom In Zoom Out Reset image size

Aubry has shown that a structural transition occurs at a critical coupling λ_c of the chain to the substrate potential [9, 10]. Below this critical coupling, the static friction vanishes because, while moving, the atoms of the chain always explore all possible positions on the substrate and thereby average out the effective corrugation and resulting friction force. One of the empirical laws that describes macroscopic friction states that kinetic friction is lower than static friction [1]. Therefore it is important to establish whether vanishing static friction does indeed imply vanishing kinetic friction also at the nanoscale. In the paper where they coined the term superlubricity, Shinjo and Hirano [11] studied the kinetic contribution to friction in the FK-model and proposed a phase diagram (bottom panel figure 1) with a frictional and a superlubric region as a function of the coupling λ and the velocity of the centre of mass of the chain. In the superlubric regime, once in motion, the chain would slide indefinitely without any kinetic friction but with a recurrent exchange of kinetic energy between the centre of mass (CM) and a single internal mode related to the periodic modulation due to the substrate. According to their findings, by increasing the velocity of the chain for weak coupling (λ < λ_c, no static friction), namely along a vertical line in the phase diagram of figure 1, the system would first acquire friction and then go back to a superlubric state at higher velocities. The proposed behaviour implies a surprising sudden change to a frictional state as soon as the system is set into motion.

This result was later criticized [12], by showing that the initial superlubric behaviour is destroyed with time, since the recurrent oscillation of the CM gives rise to the resonant parametric excitation of acoustical, long wavelength, vibrations in the chain. This phenomenon leads to a high friction regime where all the initial kinetic energy of the CM is eventually converted into internal motion that can be interpreted as heat. The resonant excitation of acoustical modes is however weaker at lower velocities and goes naturally towards the vanishing friction of the static limit. This result agrees well with the fact that most experimental observations of quasi-vanishing friction have been at the extremely low velocities of atomic force microscopes (AFM), of the order of μm s⁻¹. Indeed, the term superlubricity is currently used to indicate in general a low friction sliding rather than in the original meaning of [11].

Here we study how superlubricity is affected by temperature and velocity. First we show that the dynamical mechanism identified in [12] leads to a system in thermal equilibrium that satisfies energy equipartition over the internal vibrations, yielding a well defined temperature. We then show that at finite temperature the system displays a damped dynamics that allows one to calculate the damping coefficient from the microscopic description rather than by assuming its presence phenomenologically. We find that this damping is not a monotonic function of the sliding velocity, as it peaks close to specific phonon resonances. Away from these resonances and particularly in the limit of low velocity, the effective damping is finite, and thus at low velocities the friction vanishes, preserving superlubricity.

The paper is organized as follows. In section 2 we describe the model and its dynamics. In section 3 we review how kinetic friction arises and show how the system evolves into thermal equilibrium. In section 4 we show that, at thermal equilibrium, friction can be described as a viscous damping that can be calculated from the model. We also compare our results against damped driven models [13]. In section 5 we suggest experimental systems that could confirm the results presented in this paper and in section 6 we summarize our findings.

The FK-model [8] consists of a harmonic linear chain of particles of mass m with period a₀, subjected to an external periodic potential with period a_s and amplitude U₀/2π. The spring constant, between the particles, is K. The ratio of the depth of the potential to the spring constant of the chain defines a coupling parameter: $\lambda \equiv {U}_{0}/(K{a}_{\mathrm{s}}^{2})$. In the context of friction the particles represent the atoms of a sliding body and the potential represents a rigid solid substrate. The model is depicted schematically in the top panel of figure 1.

If the ratio of the parameters a₀/a_s is an irrational number, the system is called incommensurate; otherwise the system is called commensurate. The incommensurate case is the most interesting, because for two arbitrary surfaces in contact a common periodicity is not likely.

The ground state is determined by the competition of the harmonic interaction between the particles and the interaction of the particles with the periodic substrate, which is controlled by λ. For incommensurate systems there exists a structural phase transition [9, 10]: for λ < λ_c the system is in a so called floating phase and there is no static friction; if λ > λ_c the system is in the pinned phase and exhibits static friction. The value of λ_c is largest when the ratio between the length scales is equal to the golden mean, τ_g: ${a}_{0}/{a}_{\mathrm{s}}={\tau }_{\mathrm{g}}=(1+\sqrt{5})/2=1.6 1 8 \hspace{0.167em} 0 3 3 \hspace{0.167em} 9 8 9\ldots$ [14], for which λ_c = 0.1546... [15]. In the numerical implementation, one may approximate the golden mean by a ratio of successive Fibonacci numbers, τ_g = lim_n→∞F_n+1/F_n. Here, all calculations are done for n = 12, i.e. N = 144 particles are distributed onto 233 periods of the substrate potential.

The model can be enriched by considering the dynamics of the particles [11, 12]. We include the kinetic energy of the particles by defining τ = t/τ₀, where t is the time and ${\tau }_{0}\equiv \sqrt{m/K}$. If we express all energies in units of $K{a}_{\mathrm{s}}^{2}$, times in units of τ₀ and lengths in units of a_s, the Hamiltonian of the system becomes

$$\begin{eqnarray} &&\fl \mathcal{H}=\sum _{i=1}^{N}\biggl [\frac{1}{2}\left (\frac{\mathrm{d}{u}_{i}}{\mathrm{d}\tau }\right )^{2}+\frac{1}{2}({u}_{i+1}-{u}_{i}-{a}_{0})^{2}\nonumber\\ &&+~\frac{\lambda }{2 \pi }(1-\cos (2 \pi {u}_{i}))\biggr ], \end{eqnarray} \tag{ 1 }$$

where u_i is the position of particle i of the chain with the periodic boundary condition u_i+N = u_i. The resulting equations of motion can be written as a dynamical system of 2N first-order ordinary differential equations:

$$\begin{eqnarray} &&\frac{\mathrm{d}{u}_{i}}{\mathrm{d}\tau }={y}_{i}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&\frac{\mathrm{d}{y}_{i}}{\mathrm{d}\tau }={u}_{i+1}+{u}_{i-1}-2{u}_{i}-\lambda \sin (2 \pi {u}_{i}). \end{eqnarray} \tag{ 3 }$$

We solve these equations numerically by using the fourth-order Runge–Kutta (RK) method with a time step of τ₀/150 and calculate also the position $Q=(1/N){\mathop{\sum }\nolimits }_{i=1}^{N}{u}_{i}$ and velocity $P=(1/N){\mathop{\sum }\nolimits }_{i=1}^{N}{y}_{i}$ of the CM. As the system is Hamiltonian (no heuristic damping) the total energy is conserved. Nevertheless, energy can be transferred from the centre of mass to the internal degrees of freedom, leading to the arrest of the chain in time. This effect can be interpreted as an effective friction.

If there is no coupling between the chain and the substrate, λ = 0, the dispersion relation of the phonon modes is that of a harmonic chain, ω_k = 2|sin(k/2)|, with ω_k the frequency of the phonons and k = 2πi/N the wavevectors, in which i∈(−N/2,N/2]. When all the particles move with the same velocity p₀, the atoms slide over the periodic potential with the washboard frequency Ω = 2πp₀. Strong resonances occur if Ω is close to the frequency of the phonon mode with the wavevector q = 2πa₀/a_s or to its harmonics nq, multiples of the wavevector q. If λ is small, these frequencies can be approximated by the frequencies ω_nq of the phonons of the harmonic chain. The resonance phenomenon therefore occurs for CM velocities [13, 12]

$$\begin{eqnarray} &&n{p}_{0}\sim \frac{{\omega }_{n\mathrm{q}}}{2 \pi }\tqs n\in \mathbb{N}\quad n\geq 1. \end{eqnarray} \tag{ 4 }$$

When these resonances occur, after an initial recurrent energy exchange [11, 12], the CM translational kinetic energy is eventually totally transformed into internal kinetic and potential energy of vibration in the chain. This energy transfer occurs because the CM motion with velocity P induces a modulation with the washboard frequency Ω that leads to parametric resonances with exponential growth of phonon modes with wavevector k whenever [13, 12, 17]

$$\begin{eqnarray} &&\Omega \simeq \frac{\omega (k)+\omega (m q-k)}{m} \end{eqnarray} \tag{ 5 }$$

for some integer m. From the above relation, one can derive a maximum value P_max = Ω_max/(2π) = 4|cos(πτ_g/2)|/(2π) ≈ 0.5254. Above this value, the transfer of energy from the CM to the phonons is blocked and the motion would be really superlubric.

Lastly we note that, although the total energy is conserved, entropy increases when the energy is dispersed over the large number of internal degrees of freedom of the chain. Finite systems can show temporary decreases in the entropy, and, in principle, due to Poincaré recurrence, the system could return close to the initial state with all energy in the CM momentum. In practice, however, for systems with a sufficiently large number of internal degrees of freedom, such a large decrease in entropy is highly improbable and occurs only on extremely long timescales. Indeed, we have never observed this reversal in the length of our simulations.

In this section, we show the evolution of the CM velocity P for initial velocities close to the first four resonances ω_nq and show how this process leads to thermal equilibrium. We take λ = 0.05 ≈ λ_c/3 for which there is no static friction, namely we consider values of P on the vertical dashed line in figure 1. For the initial positions of the particles we take a ground state of the system, calculated by minimizing the potential energy. We give all the particles the same initial velocity in order to have the CM velocity P equal to the velocity corresponding to the first four resonances (table 1).

Table 1.  For each resonance n we give the initial velocity per particle p₀. The τ-range gives the time range in which the probability density function of P and the averages of the phonon amplitudes have been calculated. The absolute value of the averages (μ) of the Gaussian fits of the histograms of P is given in the fourth column in units of 10⁻⁴a_s/τ₀. These are a measure of the deviation of the system from a perfect thermal equilibrium. The last three columns are the temperatures obtained from the initial kinetic energy (T_eq), the distribution of $P~({T}_{\mathrm{eq}}^{\mathrm{CM}})$ and the phonon modes (${T}_{\mathrm{eq}}^{\mathrm{phon}}$).

n	p0 (as/τ0)	τ-range (τ0)	|μ|	Teq  (λ)	${T}_{\mathrm{eq }}^{\mathrm{CM }}~(\lambda )$	${T}_{\mathrm{eq }}^{\mathrm{phon }}~(\lambda )$
1	0.2966	(5–405) × 104	6.5 ± 2	0.88	0.89 ± 0.01	0.91 ± 0.04
2	0.1075	(25–425) × 104	1.5 ± 0.7	0.12	0.13 ± 0.01	0.13 ± 0.01
3	0.046 92	(4–16) × 106	3.1 ± 1	0.022	0.026 ± 0.002	0.030 ± 0.004
4	0.079 27	(4–8) × 106	2.0 ± 0.7	0.063	0.055 ± 0.002	0.074 ± 0.004

In figure 2 we show the time evolution P(τ) for the second resonance. We see an initial oscillatory phase, enlarged in the inset of figure 2(a), where P periodically goes back to the initial value. After some time this recurrence mechanism breaks down and P decays with time. This effect, the onset of friction, was not reported in [11], possibly because of the use of a too small number of particles or too short integration times [12]. The mechanism for this behaviour, parametric excitation of acoustic phonons, has been analyzed in [12]. A similar behaviour has also been found for sliding concentric nanotubes [16].

Zoom In Zoom Out Reset image size

When the decay has finished, P jiggles around 0. From figure 2(c) we see that the CM kinetic energy has been transformed into equal amounts of kinetic and potential energy as expected classically for vibrations, a first indication that the system has reached thermal equilibrium. For the other resonances, the evolution of P in time is qualitatively similar but with different timescales: the higher the order of resonance the longer the time before the start of the decay of P. For the fourth resonance, we had to use a slightly perturbed ground state to speed up this process. The increasingly long time that it takes for the parametric growth of the phonon excitations to lead to friction makes it very difficult to study numerically the low velocity limit. We cannot rule out superlubric motion at very low velocities where the parametric growth of internal vibrations could either occur on timescales too long for the simulations or not occur at all. Notice that these results do not support the phase diagram [11] reported in figure 1.

If, after the decay of P, the system has reached thermal equilibrium, equipartition should apply, namely the initial CM kinetic energy should be distributed over the N oscillators, each having k_BT_eq. We use two independent methods to ascertain whether this is the case.

At thermal equilibrium, the probability density function of P(τ) should be described by a Maxwell–Boltzmann (MB) distribution. In one dimension, the MB distribution coincides with a Gaussian distribution around 0, with variance ${\sigma }^{2}={k}_{\mathrm{B}}{T}_{\mathrm{eq}}^{\mathrm{CM}}/M$, where M = Nm is the mass of the chain, k_B is the Boltzmann constant and ${T}_{\mathrm{eq}}^{\mathrm{CM}}$ is the equilibrium temperature. We calculate the probability density function from a normalized histogram of P, divided into 100 uniform bins, in the time interval in which P jiggles around 0. The result, shown for the second resonance in figure 2(b), is indeed described well by a Gaussian with mean μ = 1.5 × 10⁻⁴ a_s/τ₀ and variance ${\sigma }^{2}=4.5\times 1{0}^{-5}~{a}_{\mathrm{s}}^{2}/{\tau }_{0}^{2}$. The estimated values for the equilibrium temperature are given in table 1.

Another estimate for the temperature can be given by the amplitudes, x_k, of the phonon modes [17]: ${\mathop{\sum }\nolimits }_{k=1}^{N}\frac{1}{2}m{\omega }_{k}^{2}\langle {x}_{k}^{2}\rangle =\frac{1}{2}{k}_{\mathrm{B}}{T}_{\mathrm{eq}}^{\mathrm{phon}}$, where the brackets denote a time average and ${T}_{\mathrm{eq}}^{\mathrm{phon}}$ is the equilibrium temperature corresponding to the distribution of the phonon amplitudes. In figure 2(d) we show, for the second resonance, that the amplitudes of the phonon modes are proportional to $1/{\omega }_{k}^{2}$ as expected at equilibrium.

The behaviour of the distribution of P and the phonon modes is qualitatively the same for the other resonances. In table 1 we compare the temperatures T_eq with the values obtained from the distribution of $P,{T}_{\mathrm{eq}}^{\mathrm{CM}}$, and from the phonon amplitudes, ${T}_{\mathrm{eq}}^{\mathrm{phon}}$. In order to define high and low temperature regimes, all temperatures are given in units of the value of λ used in this work, 0.05. In this way we can distinguish high (T ≈ 1 λ) and low temperature (T ≪ 1 λ) regimes. In the next section, we also study the generic value T = 0.16 λ, between the first and second resonance, for which the distribution of P was checked and found to be a Gaussian. The agreement between the different estimates demonstrates that the system has reached thermal equilibrium in the studied range of velocities around the first four resonances.

After the system has reached thermal equilibrium, as described in the previous section, we give the CM an initial velocity P(τ = 0) and study how it evolves with time. If the friction is purely viscous, F_fric =− ηMP, as is often assumed in the study of frictional dynamics [13, 18, 19], the sliding velocity decays exponentially with time with a timescale determined by the friction coefficient η. This type of viscous behaviour can be expected on general theoretical grounds for many-particle systems [20, 21].

We start from initial conditions obtained by choosing many (typically a hundred) configurations of the system in thermal equilibrium at temperature T_eq in the τ-range given in table 1. We increase the CM velocity P by adding a uniform initial sliding velocity p⁺ to all the particles. The systems are then evolved with time and the time dependence of the averaged CM velocity 〈P(τ)〉 is obtained by averaging over all initial configurations. This averaging is necessary as, in such small systems, fluctuations are large. In figure 3 we show 〈P(τ)〉 for p⁺ = 0.09 a_s/τ₀, averaged over 500 configurations at T = 0.88 λ. The average CM velocity decays exponentially, 〈P(τ)〉(a_s/τ₀) = 0.0923exp(−ητ) with a friction coefficient η = 1.94 × 10⁻³(1/τ₀).

Zoom In Zoom Out Reset image size

In this system, it is the microscopic dynamics that give rise to the friction and determine the value of the friction coefficient, without the need for imposing a phenomenological friction term. A related result for damping of phonon modes in relation to thermostats in crystalline materials was obtained and tested in the context of friction [22].

The time range used for the calculation of η needs to be not too short because of the reliability of the fit, nor too long because of the change of temperature of the system which occurs during the time range. We choose for the high temperature case, the time range τ = 0 − 500 τ₀ and for the low temperature case, the time range τ = 0 − 1000 τ₀, for determining η and its error bar. Therefore, we do not consider the effect of different time ranges for the fit, although we have checked that the behaviour of η remains qualitatively the same.

The behaviour of η as a function of the sliding velocity is shown in figure 4. For all T, in the limit p⁺ → 0,η becomes constant (as indicated by the arrows in figure 4), implying that the viscous friction, F_fric, goes to zero and superlubricity is conserved. The same conclusion has been drawn in [23] for the Tomlinson model. The friction coefficient η grows with temperature. This is not surprising, because thermal fluctuations distort the incommensurate contact between the surfaces. It is known that friction decreases with temperature, particularly in the stick–slip regime [24, 23] where the 'slips' from one minimum of the potential energy to the other can be activated by thermal energy. The attempt frequency of these processes is determined by the damping η, which is usually taken to be temperature independent.

Zoom In Zoom Out Reset image size

The vertical dashed lines in figure 4 are the velocities corresponding to the first four resonances of the system (equation (4)). For the high temperature case in figure 4(a), the friction coefficient has a peak at frequencies just above the first resonance, ω_q, and no pronounced features at the lower three resonances. For the low temperature case, figure 4(b), the influence of the lower resonances on the friction characteristics becomes apparent. Again the peaks in η are shifted to the right. In addition, the influence of temperature on the broadening of the peaks is visible: the higher the temperature, the broader the peak. For the lowest resonance studied, the fourth, there is no clear peak of η at any of the temperatures studied, whereas for the third resonance it can be seen at the lowest temperature studied. Nevertheless, at sliding velocities near the resonances, superlubricity is suppressed for all temperatures.

In section 2 we have shown that at T = 0 there is a maximum resonance velocity, P_max ≈ 0.5254. For initial velocities sufficiently far above this value, there can be no effective damping at T = 0. At finite temperature, thermal fluctuations distort the incommensurate geometry, and nonzero effective damping can still exist even for high velocities. We have indeed observed a slow decay of velocity in our simulations at p⁺ = 0.75 and T = 0.88 λ.

Our results can be compared to the velocity–force characteristics calculated by Strunz and Elmer [13] for a driven and damped FK-model at T = 0. The inherent friction in their system, calculated by subtracting the imposed damping term from the total force, shows an effect of the resonances similar to our results of figure 4. A difference between our results and those of Strunz and Elmer is the location of the peaks in the friction characteristics. In [13] these peaks are situated at the values of the resonances calculated from equation (4), while in our situation the peaks are shifted towards higher values of the sliding velocity. A possible explanation for this shift towards the right in our undamped case is the influence of the use of the approximation λ = 0 for the resonances equation (4). Indeed, for smaller coupling, λ = 0.005 and T = 8.8 λ, we have observed a smaller shift for the first resonance.

Also Strunz and Elmer used this approximation, but they studied the so called uniform sliding state in a damped system, whereas we solve the full dynamics. Our results confirm the influence of the resonances on the friction found in [13], and show that this aspect of the model is not limited to T = 0 or to uniform sliding states of a damped driven system.

Our 1D model is highly simplified and focuses on the geometrical incommensurability which is essential for superlubricity. In real systems, other aspects like roughness [6] and electronic friction [7, 25, 26] may also play a role. Nevertheless, we attempt to find possible experimental situations in which our results could be confirmed. An example of a well studied system with incommensurate contacting surfaces is a layer of Xe atoms sliding on an Ag(111) substrate [26–28]. The ratio between the lattice parameters of the Xe monolayer and the substrate is a₀/a_s = (4.55 × 10⁻¹⁰ m/2.892 × 10⁻¹⁰ m) = 1.57 ≈ τ_g. The mass of the Xe atoms is m = 2.16 × 10⁻²⁵ kg. The spring constant between the Xe atoms can be obtained from a Taylor expansion of the Lennard-Jones potential with parameters for Xe–Xe around its minimum [26], K = 1.103 J m⁻². Consequently, the typical scale of the velocity is a_s/τ₀ = 6.5 × 10² m s⁻¹. In a similar way U₀ is estimated to be in the range 1–5 × 10⁻²¹ J [26, 29], which gives 0.01 < λ < 0.05 < λ_c close to the value of λ that we used in the simulations.

The velocity where the effective friction becomes a strongly nonlinear function of velocity is rather high, orders of magnitude too high for AFM and at the limit of reachable velocities for QCM. Therefore it is important to realize that in the low velocity limit, friction grows smoothly with velocity from the vanishing static friction, contrary to the Shinjo–Hirano model of figure 1. It would be interesting however to explore also the transition from low viscous friction to high friction in correspondence with the resonances described before. The friction characteristics for the third resonance could possibly be seen at sliding velocities of approximately 30 m s⁻¹, for temperature in the region of about 8 K for Xe sliding on Ag(111). The second resonance, at around 70 m s⁻¹, would be visible up to higher temperatures, at least around 50 K.

Another interesting experimental setting for the FK-model is in colloidal suspensions, where atomic-scale systems can be reproduced on larger scales, with colloidal particles playing the role of atoms. Recently, the static FK-model was studied experimentally in such systems [30]. However, as colloidal crystals have overdamped decay of phonon modes [31], the present results cannot easily be scaled up in the same way.

We have numerically calculated effective friction coefficients as a function of temperature and velocity in the FK-model without any driving forces or phenomenological damping terms. We have critically reviewed previous findings and shown how the initial kinetic energy of the centre of mass is converted into heat by excitation of the internal motion of the chain, bringing the system to thermal equilibrium. We then studied the temperature and velocity dependence of the friction and found that it can be described as an effective viscous damping. The friction coefficient however, shows a peculiar velocity dependence with peaks at sliding velocities close to resonances between the phonons of the chain and wavevectors related to the modulation induced by the incommensurate periodic potential. We argue that these peaks in the friction coefficient could be observed experimentally. Moreover, we find that for incommensurate contacts, the effective damping increases with temperature.

Despite the fact that both sliding velocity and temperature can suppress superlubricity in incommensurate systems, the fact that the friction coefficient remains finite at low velocities means that superlubricity is preserved in the important limit of very low velocity probed by AFM experiments.

ASdW's work was financially supported by a Veni grant of the Netherlands Organization for Scientific Research (NWO) and by an Unga Forskare grant from the Swedish Research Council. This work is part of the research programme of the Stichting voor Fundamenteel Onderzoek der Materie (FOM), which is financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO).