# Statistical mechanics of rate-independent stick-slip on a corrugated surface composed of parabolic wells

Stefano Giordano

► To cite this version:

Stefano Giordano. Statistical mechanics of rate-independent stick-slip on a corrugated surface composed of parabolic wells. Continuum Mechanics and Thermodynamics, 2022, 34 (5), pp.1343-1372. 10.1007/s00161-022-01129-0. hal-03747657

HAL Id: hal-03747657

https://hal.science/hal-03747657v1

Submitted on 13 Sep 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Statistical mechanics of rate-independent
stick-slip on a corrugated surface composed
of parabolic wells

Stefano Giordano

Univ. Lille, CNRS, Centrale Lille, Univ. Polytechnique
Hauts-de-France, UMR 8520 - IEMN - Institut d'Electronique de
Microélectronique et de Nanotechnologie, F-59000 Lille, France.

Contributing authors: stefano.giordano@univ-lille.fr;

Abstract

The stick-slip phenomenon, at the basis of friction, is crucial for several
applications ranging from nanotechnology and biophysics to mechan-
ics and geology. Deep understanding of friction mechanisms and, in
particular, the methodologies for its reduction must be sought in its
nanoscopic nature, where atomic interactions and stick-slip processes
play a crucial role. At this scale, thermal fluctuations clearly have a
major effect on the physics of the problem. Hence, we develop here
a theory for rate-independent stick-slip, based on equilibrium statisti-
cal mechanics. In particular, we introduce suitably modified Prandtl-
Tomlinson and Frenkel-Kontorova models in order to study the system
with one particle and the chain with *N* particles, respectively. The
adopted corrugated substrate is composed of a sequence of quadratic
wells. Interestingly, the calculation of corresponding partition func-
tions shows a conceptual link with the theory of Jacobi and Riemann
theta functions, allowing an efficient determination of the average static
frictional force and other relevant quantities. We show some applica-
tions including the study of structural lubricity and thermolubricity.

**Keywords:** static nanofriction, Prandtl-Tomlinson and Frenkel-Kontorova models, equilibrium statistical mechanics, thermolubricity
# 1 Introduction

The stick-slip phenomenon is the basic feature of friction processes, which play a crucial role in many physical systems covering a wide range of scales, from molecular or nanoscopic to mesoscopic and macroscopic [1–4]. It rules the physics of small contacts in nanoscience and nanotechnology [5–9], the sliding of macromolecules in biological and soft structures [10–13], the mechanics of machines with interacting surfaces [14–17], and the evolution of geophysical systems [18–20]. Moreover, stick-slip and friction are also at the basis of understanding plastic phenomena in solid materials, being able to control the nucleation of dislocations and fractures and to regulate the shear transformations and the ductile-to-brittle failure transition [21–28].

From the historical point of view, the empirical laws governing the frictional force occurring between a slider and a substrate have been investigated and established by Leonardo da Vinci (1452-1519), Guillaume Amontons (1663-1705), and Charles Augustin de Coulomb (1736-1806) [29]: (i) frictional force is independent of the apparent contact area between slider and substrate, (ii) the frictional force is proportional to the loading normal force, and (iii) kinetic friction is independent of the sliding velocity and is typically smaller than static friction. These macroscopic observations are purely empirical and ultimately derive from the complex physical processes at the interface between slider and substrate, involving the attachment and detachment of asperities and the resulting variation of the real contact area [30, 31]. The nature of friction at the nanoscale has begun to be studied more systematically with the advent of the atomic force microscope [32], the surface force apparatus [33], and the development of efficient molecular dynamics simulations and multiscale models [34–37]. These approaches to nanotribology have enabled the experimental and numerical observation of the characteristic stick-slip motion of the slider atoms, which interact with the periodically corrugated potential energy representing the atoms of the substrate. This result allowed the validation of an older conceptual model, which is usually referred to as the Prandtl-Tomlinson model [38, 39]. This is the most efficient model currently known for describing stick-slip motion or nanoscale friction and, for this reason, has been largely investigated and compared with experimental and numerical results [40, 41]. It consists in a single point mass moving on a one-dimensional periodic potential (typically sinusoidal) pulled by a linear spring characterized by its elastic constant. The analysis of this model reveled different regimes of temperature and velocity and these theoretical achievements explained several atomic force microscope experiments [42–52]. The most natural generalization of the Prandtl-Tomlinson model considers a one-dimensional elastic chain of interacting atoms (or particles) moving over the periodic potential mimicking the substrate. This scheme is known as Frenkel-Kontorova model and perfectly describes the mutual sliding of two crystalline interfaces. It has been first introduced by Dehlinger [53], and then adopted by Kontorova and Frenkel to study the dislocations motion and explain the plastic deformation in crystals [54–56]. The Frenkel-Kontorova model is able
to describe a large range of nonlinear phenomena, such as solitons, chaos, friction, dislocations, incommensurate phases and glass-like systems. Because of its simplicity and universality, it has been widely studied, eventually becoming one of the most popular models used in low-dimensional nonlinear physics [57–59]. We remark that the Frenkel-Kontorova model is sometimes named discrete sine-Gordon model, since the continuous limit of the Frenkel-Kontorova model exactly yields the sine-Gordon differential equation. Interestingly, this equation was originally obtained by Edmond Bour in 1862 while conducting researches on surfaces of constant negative curvature [60], and only rediscovered by Frenkel and Kontorova in 1939 in their study on dislocations [54–56]. The tribological importance of the Frenkel-Kontorova model and of the sine-Gordon equation is that their dynamics is characterized by topological solitons named kinks (in case of the chain compression) and anti-kinks (in the case of the chain extension), which move in the chain easily than particles, controlling the overall chain motion. Indeed, the shift of the whole chain on the left or the right by one period of the corrugated substrate corresponds to the propagation of a kink (or anti-kink) along the chain. The energy barrier to be crossed to generate the kink movement is referred to as Peierls-Nabarro barrier (introduced in the study of dislocations motion) and is smaller than the amplitude of energetic oscillations of the corrugated substrate potential [61–63]. In this context, Aubry proved an important transition for an infinite incommensurate chain between a pinned state and an unpinned (frictionless) state as its elastic constant increases [64–66] (see Section 3 for details). This transition has been much studied in the past years, generalized to chains of finite length and observed experimentally [67–75].

Although the absence of friction can be dangerous in some cases, facilitating for example slipping on ice [76] or a banana peel [77], a great scientific challenge consists in reducing its value in order to limit energy consumption in technological process [16]. In this regard, the superlubricity phenomenon, corresponding to a state of very low friction between two surfaces, has been observed and thoroughly investigated [78, 79]. This very special condition between surfaces can be achieved by various techniques, the most important of which are the structural superlubricity, and the thermolubricity. In a structural superlubric state, the two interacting surfaces have different or incommensurate lattices whose mismatch leads to a nearly vanishing frictional force [80–82]. This effect has been experimentally observed by scanning tunneling microscopy for W(011) tip on Si(001) surfaces [83], by frictional force microscopy for metal tips over graphite [84, 85] and by noncontact atomic force microscopy for graphene nanoribbons on gold surfaces [86], just to name a few. These results have been complemented by computational studies [87, 88]. The effect of temperature on friction and the thermolubricity phenomenon are well described by the Prandtl-Tomlinson minimalistic model introduced above [38, 39]. Actually, in the early pioneering article by Prandtl [38], one already finds that friction decreases with temperature since thermal energy fosters the crossing of energy barriers and eventually enables slip [42–52]. Once again, experimental
4 *Statistical mechanics of rate-independent stick-slip*

evidence for this principle has been obtained for various materials using the
atomic force microscope [89, 90].

In order to have a comprehensive view of the effect of temperature on friction, we develop here a theory for the rate-independent stick-slip phenomenon, based on equilibrium statistical mechanics. To do this, we implement the so-called spin variable approach useful to deal with nonconvex potential energies [91]. This technique has already been successfully applied to different problems such as the physics of muscles [92, 93], the folding of macromolecules [94–97], the adhesion processes [98, 99] and the phase transformations in solids [100, 101]. We remark that this approach is complementary to the methodologies commonly used to study the behavior of physical systems with multiple stable and metastable states [102–105]. Here, we modify the Prandtl-Tomlinson and Frenkel-Kontorova models by substituting the sinusoidal energy profile of the substrate with a sequence of quadratic potentials, each of which represents a substrate well. It is important to remark that the sinusoidal profile and the sequence of quadratic wells have a qualitatively similar behavior due to their periodicity but different features from the quantitative point of view because of the quite different shape. In this sense, the quadratic substrate should not be seen as an approximation of the sinusoidal one, but rather as an independent model with application, e.g., to frictional motion of defects in solids or sliding of macromolecules on biological structures [12, 28]. Our approach allows a simpler calculation of the partition functions for both (modified) Prandtl-Tomlinson and Frenkel-Kontorova systems. The spin variables approach requires the introduction of a discrete variable (spin) for each particle of the system, useful to identify the substrate well locally explored by the particle itself. The application of this methodology to the (modified) Prandtl-Tomlinson and Frenkel-Kontorova models leads to explicit partition functions written in terms of Jacobi or Riemann theta functions. This allows an efficient calculation of all thermodynamic variables, including the average frictional force (when the energy barrier between the wells is sufficiently larger than $K_B T$). Eventually, this force depends on temperature and other geometrical and physical variables, thus explaining possible structural and thermal lubricity phenomena. In particular, we proved the existence for finite chains of specific commensurate values of the lattice constant ratio between slider and substrate, which are able to sensibly reduce the friction. We finally studied the effect of the temperature on this friction reduction. The models introduced, although rather simple, offer the future possibility of taking into account more complex phenomena not accessible to current theories. We refer for example to periodic but complex energy profiles and corrugated but deformable substrates where the slider affects the shape of the substrate itself (this is relevant in soft matter applications).

The paper is structured as follows. In Section 2, we develop the theory for
the modified Prandtl-Tomlinson model, and we discuss the results in Section
3. Then, in Section 4, we elaborate the statistical mechanics for the modified
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Scheme of the modified Prandtl-Tomlinson model, characterized by a sequence of quadratic potentials with stiffness $k_0$ and placed at positions $jd, \forall j \in \mathbb{Z}$. The energy barrier $W$ characterizes the corrugated substrate and must be sufficiently larger than $K_B T$. The particle at point P is linked by a first linear spring (constant $k_b$) to the origin of the x-axis, and by a second linear spring (constant $k$) to the position $\ell$, representing the sliding device.</figcaption></figure>

Frenkel-Kontorova model, and we show some applications in Section 5. The conclusions (Section 6) and a mathematical Appendix close the paper.

## 2 One-dimensional stick-slip model

We introduce a simple model concerning the one-dimensional stick-slip behavior of a single particle in contact with a substrate described by a periodic potential energy. In particular, we propose here a modification of the classical Prandtl-Tomlinson model [38–41], able to consider the effect of thermal fluctuations on the static system behavior. Indeed, these fluctuations can significantly modify the friction experienced by the particle deposited on the substrate. While the interplay between temperature and velocity in the Prandtl-Tomlinson model has been largely studied by means of several different approaches [42–51], the static stick-slip behavior is directly investigated here through statistical mechanics. In our approach, we assume to be close to the thermodynamic equilibrium and therefore we consider rate independent processes observed with the slider at rest or at sufficiently low speed. The characteristic sinusoidal profile of the Prandtl-Tomlinson substrate potential is substituted here by a periodic sequence of quadratic wells with elastic constant $k_0$ and spacing $d$ (see Fig.1). The energy barrier between adjacent wells is given by $W = k_0d^2/8$ (it must be much larger than $K_B T$, as discussed below). Of course, we do not consider the natural transition rates between adjacent wells (Kramers rate theory [106]) and the interaction of this dynamical processes with the (low) speed of the traction device. These effects can be taken into account only by means of the out of equilibrium statistical mechanics, and this analysis has not been performed here. In order to identify the well occupied by the particle, we have to introduce a new discrete variable $n \in \mathbb{Z}$ within the phase space of the system, which represents the number corresponding to the concerned substrate site. It is important to underline that the energy barrier of the corrugated substrate (quantified by $W$ or $k_0$) has been directly related to the normal load applied between slider and substrate [52]. This gives still more importance to the Prandtl-Tomlinson model, whose features can be easily put in relation to the experimental parameters.
6 Statistical mechanics of rate-independent stick-slip

In our model, the manipulation of the particle is performed by a first linear spring (constant $k_b$) tethering the particle at the origin of the x-axis, and a second linear spring (constant $k$) pulling the particle at position $\ell$ (see Fig.1). These springs represent the device used to investigate the stick-slip behavior. Of course, if we want to simulate a non-constrained particle interacting with the substrate, we can let $k_b = 0$, thus eliminating the bonding at the origin. Moreover, $\ell$ must be considered slowly variable to allow the particle to explore quasi-statically the substrate. These premises enable us to write the potential energy of the system as

$$U(x, n) = \frac{1}{2}k_0(x-dn)^2 + \frac{1}{2}k_b x^2 + \frac{1}{2}k(\ell-x)^2, \quad (1)$$

with $x \in \mathbb{R}$ and $n \in \mathbb{Z}$. We suppose that the system is in contact with a thermal reservoir at temperature $T$ and we denote the Boltzmann constant by $K_B$. Hence, to introduce the partition function of the canonical distribution, we have to sum over the discrete variable $n$ and to integrate over the continuous one $x$. Importantly, since we integrate the variable $x$ for any potential well from $-\infty$ to $+\infty$, we have to consider an energy barrier $W$ sufficiently larger than the thermal energy $K_B T$. So doing, the effects generated by the overlapping of quadratic potentials are negligible. For example, if we consider a silicon cantilever sliding on a NaCl single crystal, we can observe that a normal load in the range between 0 and 6 nN corresponds to an energy barrier $W$ between 0 and 0.6 eV [52]. Then, it is not difficult to prove that the relation $W \gg K_B T$ (we can say $W = 10K_B T$) is verified when $W > 0.3$ eV, which means a normal load larger than 3 nN. This proves the applicability of our theory to real cases. In addition, we work within the Helmholtz ensemble of statistical mechanics since the distance $\ell$ is prescribed [107, 108]. These remarks lead to the expression

$$Z_H(\ell) = \sum_{n=-\infty}^{+\infty} \int_{-\infty}^{+\infty} e^{-\frac{U(x,n)}{K_B T}} dx, \quad (2)$$

where we use Eq.(1). The integral can be now performed by means of the classical gaussian integration

$$\int_{-\infty}^{+\infty} e^{-\alpha x^2} e^{\beta x} dx = \sqrt{\frac{\pi}{\alpha}} e^{\frac{\beta^2}{4\alpha}}, \quad (\alpha > 0), \quad (3)$$

which is able to deliver the following result. To write it in compact form, we define the parameters

$$\varphi = \frac{k}{2K_B T} \frac{k_0 + k_b}{k_t}, \quad \eta = \frac{k_0}{2K_B T} \frac{k_b + k}{k_t}, \quad \xi = \frac{1}{2K_B T} \frac{k_0 k}{k_t}, \quad (4)$$
where $k_t = k_0 + k_b + k$. We therefore observe that the partition function can be written as follows

$$Z_H(\ell) = e^{-\varphi\ell^2} \sum_{n=-\infty}^{+\infty} e^{-\eta d^2 n^2} e^{2\xi d\ell n} = e^{-\varphi\ell^2} \vartheta_3\left(-i\xi d\ell, \frac{i}{\pi}\eta d^2\right), \quad (5)$$

where we have neglected the non influential multiplicative constant and being $\vartheta_3(z, \tau)$ the third Jacobi theta function defined below

$$\vartheta_3(z, \tau) = \sum_{n=-\infty}^{+\infty} e^{\pi i n^2 \tau} e^{2inz} = 1 + 2 \sum_{n=1}^{+\infty} e^{\pi i n^2 \tau} \cos(2nz). \quad (6)$$

Here $z = -i\xi d\ell \in \mathbb{C}$ is the argument, and $\tau = i\eta d^2/\pi \in \mathbb{C}$ is the so-called lattice parameter satisfying the condition $\Im(\tau) > 0$ [109–113]. Sometimes, also the nome $q = e^{i\pi\tau} = e^{-\eta d^2}$ is introduced with the assumption $|q| < 1$, assuring the series convergence. The introduced function satisfies the two relations

$$\vartheta_3(z + \pi, \tau) = \vartheta_3(z, \tau), \quad (7)$$

$$\vartheta_3(z + \pi\tau, \tau) = \frac{e^{-2iz}}{q} \vartheta_3(z, \tau), \quad (8)$$

stating that the function is completely determined in the entire complex plane by the values it assumes in the parallelogram identified by the four points $z_0, z_0 + \pi\tau, z_0 + \pi + \pi\tau$ and $z_0 + \pi \forall z_0 \in \mathbb{C}$ (the fundamental domain) [109–113]. There is a wide variety of notations for the theta functions. Here, rather than the original Jacobi notation, we prefer the more modern notation introduced in Refs. [109–113]. Interestingly, theta functions, historically used for the study of elliptic functions and various problems in number theory, are also able to completely describe the statistical mechanics of rate-independent friction. From the thermodynamics of the system we know that the average value of the force necessary to maintain the position $\ell$ of the device is given by

$$\langle F_s \rangle = -K_B T \frac{\partial}{\partial \ell} \log Z_H(\ell) = -K_B T \frac{1}{Z_H(\ell)} \frac{\partial Z_H(\ell)}{\partial \ell}. \quad (9)$$

By using the following relation giving the logarithmic derivative of the third theta function [109–113]

$$\frac{\vartheta'_{3}(z, \tau)}{\vartheta_{3}(z, \tau)} = 4 \sum_{n=1}^{+\infty} (-1)^{n} \frac{q^{n}}{1 - q^{2n}} \sin(2nz), \quad (10)$$
8 Statistical mechanics of rate-independent stick-slip

where $\vartheta'_3 = \frac{\partial\vartheta_3}{\partial z}$, $q = e^{i\pi\tau}$, and which is valid for $|\Im(m(z))| < \frac{\pi}{2}\Im(m(\tau))$, we easily obtain

$$ \langle F_s \rangle = \frac{dk_0k}{k_0 + k_b + k} \left[ \frac{\ell}{d} \frac{k_0 + k_b}{k_0} + 2 \sum_{n=1}^{+\infty} (-1)^n \frac{e^{-n\eta d^2}}{1 - e^{-2n\eta d^2}} \sinh(2n\xi d\ell) \right], (11) $$

with validity in the range $|\xi| < \eta/2$. For the simplest case in which the particle is not bound to the origin of the reference system ($k_b = 0$), we obtain

$$ \langle F_s \rangle = \frac{dk_0k}{k_0 + k} \left[ \frac{\ell}{d} + 2 \sum_{n=1}^{+\infty} (-1)^n \frac{e^{-n\eta d^2}}{1 - e^{-2n\eta d^2}} \sinh(2n\ell d\eta) \right], \quad (12) $$

with $\eta = \frac{1}{2K_B T} \frac{k_0 k}{k_0 + k}$ and $|\ell| < d/2$.

Concerning the average value of the discrete variable $n$, directly from Eq.(5), we can write

$$ \langle n \rangle = \frac{\sum_{n=-\infty}^{+\infty} n e^{-n\eta d^2} e^{2\xi d\ell n}}{\sum_{n=-\infty}^{+\infty} e^{-n\eta d^2} e^{2\xi d\ell n}}, \quad (13) $$

or, equivalently

$$ \langle n \rangle = \frac{1}{2i} \frac{\vartheta'_3 \left(-i\xi d\ell, \frac{i}{\pi}\eta d^2\right)}{\vartheta_3 \left(-i\xi d\ell, \frac{i}{\pi}\eta d^2\right)}. \quad (14) $$

Using again Eq.(10), we get

$$ \langle n \rangle = 2 \sum_{n=1}^{+\infty} (-1)^{n+1} \frac{e^{-n\eta d^2}}{1 - e^{-2n\eta d^2}} \sinh(2n\xi d\ell), \quad (15) $$

which is valid if $|\xi| < \eta/2$. If $k_b = 0$, this result can be easily simplified and gives

$$ \langle n \rangle = 2 \sum_{n=1}^{+\infty} (-1)^{n+1} \frac{e^{-n\eta d^2}}{1 - e^{-2n\eta d^2}} \sinh(2n\eta d\ell), \quad (16) $$

with $\eta = \frac{1}{2K_B T} \frac{k_0 k}{k_0 + k}$ and $|\ell| < d/2$.

Although correct and rather interesting, Eqs.(11)-(12) and (15)-(16) are not very useful to give a physical interpretation of the phenomenon under study. To get a new form for these solutions, we can introduce the following Jacobi functional identity for the theta function [112, 113]

$$ \vartheta_3(z, \tau) = \frac{1}{\sqrt{-i\tau}} e^{\frac{z^2}{\pi i \tau}} \vartheta_3\left(\frac{z}{\tau}, -\frac{1}{\tau}\right), \quad (17) $$
where the square root is to be interpreted as the principal value (see the Appendix for details). If we apply this identity to our partition function given in Eq.(5), we eventually obtain

$$Z_H(\ell) = e^{-\varphi\ell^2} \frac{1}{d} \sqrt{\frac{\pi}{\eta}} e^{\frac{\xi^2}{\eta}\ell^2} \vartheta_3\left(-\pi\frac{\xi}{\eta}\frac{\ell}{d}, i\frac{\pi}{\eta d^2}\right), \quad (18)$$

which is a new form completely equivalent to the previous one. From Eq.(9), by using Eq.(10) we can obtain a new expression for the static friction force as follows

$$\langle F_s \rangle = \frac{k_b k \ell}{k_b + k} + \frac{4\pi K_B T k}{d(k_b + k)} \sum_{n=1}^{+\infty} \frac{(-1)^{n+1} e^{-\frac{\pi^2 n}{\eta d^2}}}{1 - e^{-\frac{2\pi^2 n}{\eta d^2}}} \sin \left( 2\pi n \frac{\xi}{\eta} \frac{\ell}{d} \right). \quad (19)$$

After the use of the Jacobi functional identity, this expression, contrary to Eq.(11), is always convergent and is composed of a linear term in $\ell$ added to a Fourier series in the same variable $\ell$. The convergence everywhere of the new form in Eq.(19) depends on the fact that we are now dealing with a Fourier series, whereas before in Eq.(11) we had an arbitrary function series with limited convergence. The first term represents the effect of the bonding spring $k_b$ whereas the second one comes from the periodic behavior of the corrugated substrate. The physical interpretation of Eq.(19) is therefore evident and transparent. When $k_b = 0$, the linear term vanishes and we obtain

$$\langle F_s \rangle = \frac{4\pi K_B T}{d} \sum_{n=1}^{+\infty} (-1)^{n+1} \frac{e^{-\frac{\pi^2 n}{\eta d^2}}}{1 - e^{-\frac{2\pi^2 n}{\eta d^2}}} \sin \left( 2\pi n \frac{\ell}{d} \right), \quad (20)$$

which is the periodic rate-independent friction force. This expression allows a simple analysis of the low temperature regime. Indeed, by using the limit $\lim_{x \to 0} x/(1 - e^{-ax}) = 1/a$, Eq.(20) delivers

$$\lim_{T \to 0} \langle F_s \rangle = \frac{k_0 k d}{2\pi(k_0 + k)} \sum_{n=1}^{+\infty} (-1)^{n+1} \frac{2}{n} \sin \left( 2\pi n \frac{\ell}{d} \right), \quad (21)$$

where we can recognize the Fourier series of a sawtooth wave

$$x = \sum_{n=1}^{+\infty} (-1)^{n+1} \frac{2}{n} \sin nx, \quad x \in (-\pi, +\pi), \quad (22)$$

thus obtaining $\lim_{T \to 0} \langle F_s \rangle = k_0 k \ell / (k_0 + k)$, for any $\ell$ such that $-1/2 < \ell/d < 1/2$. This simple result exactly corresponds to the purely mechanical force (without temperature effects) necessary to maintain the particle at position $\ell$ when it is attracted by the well centered at $x = 0$ (elastic constant $k_0$)
and pulled by the spring with constant $k$. The theory is therefore perfectly consistent with the pure mechanics at $T = 0$.

We can now use Eq. (14) combined with Eq. (18) to find the new expression for the average value $\langle n \rangle$ of the discrete or spin variable $n$. The result is

$$ \langle n \rangle = \frac{\ell}{d} \frac{k}{k_b + k} + \frac{\pi}{2\eta d^2} \sum_{n=1}^{+\infty} \frac{(-1)^n e^{-\frac{n^2 n}{\eta d^2}}}{1 - e^{-\frac{2\pi^2 n}{\eta d^2}}} \sin \left( 2\pi n \frac{\xi \ell}{\eta d} \right), \quad (23) $$

which is again always convergent and composed of a linear term and a periodic one. When $k_b = 0$, we finally obtain

$$ \langle n \rangle = \frac{\ell}{d} + \frac{\pi}{2\eta d^2} \sum_{n=1}^{+\infty} (-1)^n \frac{e^{-\frac{\pi^2 n}{\eta d^2}}}{1 - e^{-\frac{2\pi^2 n}{\eta d^2}}} \sin \left( 2\pi n \frac{\ell}{d} \right), \quad (24) $$

where the linear term now remains present because it represents the increase of $\langle n \rangle$ during the sliding, i.e. the sequential transitions between the different wells of the substrate potential. In conclusion, the latest results obtained through the Jacobi functional relation are represented by Fourier series that always converge, converge quickly, and have an immediate physical interpretation.

## 3 Results for the modified Prandtl-Tomlinson model

All obtained expressions can be numerically implemented and some examples can be found in Fig. 2. In panels (a) and (b) of Fig. 2 one can find the behavior of the static force $\langle F_s \rangle / (k_0 d)$ and the spin variable $\langle n \rangle$ versus $\ell/d$ for a system with $k_b/k_0 = 0$ (without bonding spring). It means that the device slowly spans some periods of the corrugated substrate and the main physical quantities are calculated for different temperatures. By observing the energy profile of the substrate, we see that the transitions of the particle between different wells occur at $\ell/d = 1/2, 3/2, \dots$, and so on. Then, in correspondence to these positions, we have peaks in the force [see panel (a) of Fig. 2] and switches between two adjacent integers in the average spin variable [see panel (b) of Fig. 2]. This is coherent with the Fourier series in Eqs. (20) and (24). Indeed, since $k_b/k_0 = 0$, the curve representing $\langle F_s \rangle / (k_0 d)$ in panel (a) is periodic in $\ell/d$ with period 1. In panels (c) and (d) of Fig. 2 we can find the behavior of the same quantities $\langle F_s \rangle / (k_0 d)$ and $\langle n \rangle$ versus $\ell/d$ for a system with $k_b/k_0 = 0.2$ (with a spring tethering the particle at the axis origin). In this case, the periodicity of the friction force is evidently lost, since a larger device force will be necessary moving away from the origin of the axis, as described by Eq. (19), with $k_b > 0$. It is important to remark that, in both panels (a) and (c) of Fig. 2, we observe a reduction of the friction force with an increasing temperature, clearly explaining the thermolubricity phenomenon [78, 79]. In practice, the increase in temperature promotes thermal fluctuations and thus
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Thermal stick-slip behavior of the modified Prandtl-Tomlinson model. Panels a) and b): average adimensional static force $\langle F_s \rangle / (k_0 d)$ and average spin variable $\langle n \rangle$ versus $\ell/d$ for a system with $k_b/k_0 = 0$. We used Eq.(12) or (20) for the force and Eq.(16) or (24) for the average value of $n$. Panels c) and d): the same quantities $\langle F_s \rangle / (k_0 d)$ and $\langle n \rangle$ versus $\ell/d$ for a system with $k_b/k_0 = 0.2$. We used Eq.(11) or (19) for the force and Eq.(15) or (23) for the average value of $n$. In all panels a)-d) we adopted $k/k_0 = 1$ and $K_B T / (k_0 d^2)$ assumes 8 values in the range (0.01, 0.12). Panel e): behavior of the maximum static friction force as function of $k/k_0$ and $K_B T / (k_0 d^2)$ with $k_b/k_0 = 0$. Here, $k/k_0 = 5^\zeta$ with ten equispaced values of $\zeta$ in $(-1, 1)$, and $W > K_B T$.</figcaption></figure>

the wider exploration of substrate energy basins. This then allows energy barriers between wells to be more easily overcome and ultimately reduces friction. We also note that in both in panels (b) and (d) the average value $\langle n \rangle$ exhibits a sharp behavior with net transitions between the adjacent wells for low values of the temperature, and a more blurred response for high temperatures. This is consistent with Eqs.(23) and (24).

To conclude, in panel (e), we show the behavior of the friction in terms of the temperature and of the elastic constants ratio $k/k_0$. To do this, we quantify the static friction force by means of the maximum value $\max \{\langle F_s \rangle / (k_0 d), 0 < \ell/d \le 1\}$ and we represent it as function of $K_B T / (k_0 d^2)$
12 Statistical mechanics of rate-independent stick-slip

and $k/k_0$. It is important to recognize that this maximum force value corresponds exactly to the force that the device must be able to apply to allow the particle to advance from one period to another of the substrate corrugated potential. This concept is in close analogy with the Peierls-Nabarro stress introduced for the study of the onset of movement of dislocations [61–63]. We observe that the maximum value of $\langle F_s \rangle / (k_0 d)$ in the period is decreasing with an increasing adimensional temperature $K_B T / (k_0 d^2)$ and is increasing with an increasing ratio $k/k_0$. Importantly, we remark that the validity of the results is in the range with $W$ sufficiently larger than $K_B T$. In Fig.2, we explored the behavior of the system up to the limit $W \simeq K_B T$ in order to better show the trend of the obtained results. The plot in panel (e) is useful to study the behavior of the system when $k$ is variable and $k_0$ is constant. But we can show that the dual behavior can be deduced directly from the previous one. Indeed, if we define the functions

$$G(x, y) = \frac{x}{1+x} \left[ \frac{\ell}{d} + 2 \sum_{n=1}^{+\infty} \frac{(-1)^n e^{-\frac{n}{2y} \frac{x}{1+x}}}{1 - e^{-\frac{n}{y} \frac{x}{1+x}}} \sinh \left( \frac{n \frac{\ell}{d} \frac{x}{y}}{1+x} \right) \right], \quad (25)$$

$$F(x, y) = 4\pi y \sum_{n=1}^{+\infty} (-1)^{n+1} \frac{e^{-2\pi^2 n y \frac{1+x}{x}}}{1 - e^{-4\pi^2 n y \frac{1+x}{x}}} \sin \left( 2\pi n \frac{\ell}{d} \right), \quad (26)$$

we can rewrite Eqs.(12) and (20) as

$$\frac{\langle F_s \rangle}{k_0 d} = G \left( \frac{k}{k_0}, \frac{K_B T}{k_0 d^2} \right) = F \left( \frac{k}{k_0}, \frac{K_B T}{k_0 d^2} \right), \quad (27)$$

which represents the plot shown in panel (e) of Fig.2. Similarly, from Eqs.(12) and (20) we can also easily obtain the dual result

$$\frac{\langle F_s \rangle}{kd} = G \left( \frac{k_0}{k}, \frac{K_B T}{kd^2} \right) = F \left( \frac{k_0}{k}, \frac{K_B T}{kd^2} \right), \quad (28)$$

which can be used when we consider $k_0$ variable and $k$ constant. We remark that the function $G$ of Eq.(25) has been obtained with the original partition function in Eq.(5), while the function $F$ of Eq.(26) with the partition function in Eq.(18), modified by means of the Jacobi functional relation. Since the behavior of Eqs.(27) and (28) is exactly the same, we deduce that in the panel (e) of Fig.2 we can substitute $\langle F_s \rangle / (k_0 d)$, $k/k_0$ and $K_B T / (k_0 d^2)$ with $\langle F_s \rangle / (kd)$, $k_0/k$ and $K_B T / (kd^2)$ and the curves remain unchanged. Therefore, the maximum value of $\langle F_s \rangle / (kd)$ is decreasing with an increasing adimensional temperature $K_B T / (kd^2)$ and is increasing with an increasing ratio $k_0/k$. This result is typically observed in several experimental measurements where the friction is an increasing function of the relative corrugation $k_0/k$ [44, 52].
<figure><img src="image_4.png" /><figcaption>Fig. 3 Scheme of the modified Frenkel-Kontorova model, characterized by a sequence of quadratic potentials with stiffness $k_0$ and placed at positions $jd, \forall j \in \mathbb{Z}$. The energy barrier $W$ characterizes the corrugated substrate and must be sufficiently larger than $K_B T$. The mass-spring chain (elastic constant $k_m$, equilibrium length $s$), placed at $x_1, ..., x_N$, is linked by a first linear spring (constant $k_1$) to the position $\ell_1$ (left device), and by a second linear spring (constant $k_N$) to the position $\ell_N$ (right device).</figcaption></figure>

# 4 Thermally modified Frenkel-Kontorova model

We consider a one-dimensional chain composed of $N$ particles linked by linear springs with elastic constant $k_m$ and equilibrium length $s$ (see Fig.3). This chain interacts with a substrate lattice composed by a sequence of pinning sites with stiffness $k_0$ and spaced by a uniform distance $d$ (see Fig.3). This arrangement corresponds to the Frenkel-Kontorova model, where a competition between the intrinsic chain spacing $s$ and the substrate lattice spacing $d$ is introduced [53–59]. At the same time, the elastic potential energy of the moving chain competes with the pinning potential energy of the substrate. At the extremity of the chain, we also introduce two more linear springs, with elastic constant $k_1$ and $k_N$, which are mimicking the devices adopted to manipulate the chain at the nanoscale. They try to keep the two ends in the positions $\ell_1$ and $\ell_N$, respectively to the left and right of the chain (see Fig.3). Actually, in order to study the stick-slip phenomenon at the nanoscale, it is sufficient to use only one device, let’s say at the right of the chain, but we keep in the formalism the two lateral springs for symmetry reasons. In the original Frenkel-Kontorova model, the lattice substrate is represented by a periodic potential energy [53–59].

$$U_{FK}(\vec{x}) = -\frac{1}{2}W \sum_{j=1}^{N} \cos(\frac{2\pi x_j}{d}), \qquad (29)$$

where $W$ is the energy barrier between adjacent sites and $x_j$ are the particle positions. As before, we define here the substrate potential energy by introducing a sequence of quadratic wells (elastic constant $k_0$). Then, we introduce thermal fluctuations in the system and we work at thermodynamic equilibrium, corresponding to rate-independent or static friction. In order to identify the pinning site concerned for each particle, we add $N$ discrete or spin variables $n_1, ..., n_N$ to the phase space of the system. It means that the value of $n_j$ corresponds to the number of the site occupied by the $j$-th particle. To draw a comparison with the original Frenkel-Kontorova model, we observe that the
14 Statistical mechanics of rate-independent stick-slip

energy barrier in our modified approach is given by $W = k_0d^2/8$, being defined by the intersection of two adjacent quadratic wells. As previously discussed, $W$ must be sufficiently larger than the thermal energy $K_B T$. These assumptions lead to total energy

$$ U(\vec{x}, \vec{n}) = \sum_{j=1}^{N} \frac{1}{2} k_0 (x_j - dn_j)^2 + \sum_{j=1}^{N-1} \frac{1}{2} k_m (x_{j+1} - x_j - s)^2 \\ + \frac{1}{2} k_1 (x_1 - \ell_1)^2 + \frac{1}{2} k_N (\ell_N - x_N)^2, \qquad (30) $$

where $\vec{x} = (x_1, ..., x_N) \in \mathbb{R}^N$ is the particle positions vector and $\vec{n} = (n_1, ..., n_N) \in \mathbb{Z}^N$ is the spin variables vector. We stress again that both those vectors belong to the phase space of the system. Then, the Helmholtz partition function can be calculated as

$$ Z_H = \sum_{\vec{n} \in \mathbb{Z}^N} \int_{\mathbb{R}^N} e^{-\frac{U(\vec{x}, \vec{n})}{k_B T}} d\vec{x}, \qquad (31) $$

where we integrated over $\vec{x}$ and summed over $\vec{n}$. To simplify the calculations, it is not difficult to rewrite the total energy $U(\vec{x}, \vec{n})$ in the following compact form

$$ U(\vec{x}, \vec{n}) = \frac{1}{2} k_m \vec{x} \cdot A \vec{x} + sk_m \vec{v} \cdot \vec{x} - k_0 d \vec{n} \cdot \vec{x} + C. \qquad (32) $$

Here, we introduced the tridiagonal matrix $A$ defined below

$$ A = \begin{bmatrix} a_1 & -1 & 0 & \dots & 0 \\ -1 & a_2 & -1 & \dots & \vdots \\ 0 & -1 & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & a_{N-1} & -1 \\ 0 & \dots & 0 & -1 & a_N \end{bmatrix} \in \mathcal{M}_{N,N}(\mathbb{R}), \qquad (33) $$

with all the subdiagonal and superdiagonal elements equal to -1 and the diagonal elements given by

$$ \left\{ \begin{array}{l} a_1 = 2 + \frac{k_1+k_0}{k_m}, \\ a_2 = \dots = a_{N-1} = 2 + \frac{k_0}{k_m}, \\ a_N = 1 + \frac{k_N+k_0}{k_m}. \end{array} \right. \qquad (34) $$

Moreover, we defined the vector $\vec{v}$ as

$$ \vec{v} = \left( 1 - \frac{k_1 l_1}{k_m s}, 0, \dots, 0, -1 - \frac{k_N l_N}{k_m s} \right) \in \mathbb{R}^N, \qquad (35) $$
and the scalar quantity (independent of $\vec{x}$)

$$C = \frac{1}{2}k_m(N-1)s^2 + \frac{1}{2}k_0d^2\vec{n} \cdot \vec{n} + \frac{1}{2}k_1\ell_1^2 + \frac{1}{2}k_N\ell_N^2. \quad (36)$$

These premises allow us to rewrite the partition function in the form

$$Z_H = \sum_{\vec{n} \in \mathbb{Z}^N} e^{-\frac{C}{k_B T}} \int_{\mathbb{R}^N} e^{-\frac{k_m s^2}{2k_B T} \vec{x} \cdot A \vec{x}} e^{-\frac{k_0 d \vec{n} \cdot \vec{x} - sk_m \vec{x} \cdot \vec{n}}{k_B T}} d\vec{x}, \quad (37)$$

that can be simplified using the gaussian integration

$$\int_{\mathbb{R}^N} e^{-\frac{1}{2} \vec{y} \cdot \mathcal{N} \vec{y}} e^{\vec{b} \cdot \vec{y}} d\vec{y} = \sqrt{\frac{(2\pi)^N}{\det \mathcal{N}}} e^{\frac{1}{2} \vec{b} \cdot \mathcal{N}^{-1} \vec{b}}, \quad (38)$$

holding for any symmetric and positive definite matrix $\mathcal{N}$. This expression represents the multidimensional generalization of Eq. (3). The integration leads to the result

$$Z_H = \sqrt{\frac{(2\pi K_B T)^N}{k_m^N \det \mathcal{A}}} e^{\frac{k_m s^2}{2K_B T} \vec{v} \cdot A^{-1} \vec{v}} e^{-\frac{k_m s^2}{2K_B T} (N-1)} e^{-\frac{k_1 \ell_1^2}{2K_B T}} e^{-\frac{k_N \ell_N^2}{2K_B T}} \\ \times \sum_{\vec{n} \in \mathbb{Z}^N} e^{-\frac{k_0 d^2}{2K_B T} \vec{n} \cdot (\mathcal{I} - \frac{k_0}{k_m} A^{-1}) \vec{n}} e^{-\frac{k_0 s d}{K_B T} \vec{n} \cdot A^{-1} \vec{v}}, \quad (39)$$

where $\mathcal{I} \in \mathcal{M}_{N,N}(\mathbb{R})$ is the identity matrix. The partition function can be finally written in terms of the so-called Riemann theta function [112, 113]

$$\Theta(\vec{z}|\Omega) = \sum_{\vec{n} \in \mathbb{Z}^N} e^{2\pi i (\frac{1}{2}\vec{n} \cdot \Omega \vec{n} + \vec{n} \cdot \vec{z})}, \quad (40)$$

which is defined when the matrix $\Omega$ is in the following set: $\mathcal{H}_N = \{\Omega \in \mathcal{M}_{N,N}(\mathbb{C}) : \Omega = \Omega^T, \Im(\Omega) > 0\}$. This is the set of symmetric square matrices whose imaginary part is positive definite. The set $\mathcal{H}_N$ is called the Siegel upper half-space and is the multi-dimensional analog of the upper half-plane in $\mathbb{C}$. Here, $\vec{z} \in \mathbb{C}^N$ is an arbitrary complex vector. It is known that the Riemann theta function converges absolutely and uniformly on compact subsets of $\mathbb{C}^N \times \mathcal{H}_N$ [112, 113]. It can be seen as the multidimensional generalization of the function $\vartheta_3(z, \tau)$ introduced in Eq. (6) to study the single particle friction. Through the identifications

$$\Omega = \frac{k_0 d^2 i}{2\pi K_B T} \left( \mathcal{I} - \frac{k_0}{k_m} A^{-1} \right), \quad (41)$$

$$\vec{z} = \frac{k_0 s d i}{2\pi K_B T} A^{-1} \vec{v}, \quad (42)$$
we obtain the partition function as

$$Z_H = \sqrt{\frac{(2\pi K_B T)^N}{k_m^N \det \mathcal{A}}} e^{\frac{k_m s^2}{2K_B T} \vec{v} \cdot \mathcal{A}^{-1} \vec{v}} e^{-\frac{k_m s^2}{2K_B T} (N-1)} e^{-\frac{k_1 \ell_1^2}{2K_B T}} e^{-\frac{k_N \ell_N^2}{2K_B T}} \\ \times \Theta \left( \frac{k_0 s d i}{2\pi K_B T} \mathcal{A}^{-1} \vec{v} \left| \frac{k_0 d^2 i}{2\pi K_B T} \left( \mathcal{I} - \frac{k_0}{k_m} \mathcal{A}^{-1} \right) \right| \right). \quad (43)$$

Since it is possible to prove that the real matrix $\mathcal{I} - \frac{k_0}{k_m} \mathcal{A}^{-1}$ is positive definite, $\Omega$ defined in Eq.(41) belongs to $\mathcal{H}_N$, and then $\Theta(\vec{z}|\Omega)$ in Eq.(43) is always convergent. The Riemann theta function is able therefore to fully describe the complexity of the rate-independent stick-slip phenomena modeled by our thermally modified Frenkel-Kontorova approach. In particular, it perfectly handles the geometric and energetic compromise between chain and substrate.

We can give a physical interpretation of the exponential term summed in Eq.(39), i.e. summed within the Riemann theta function. If we look at the total energy of the system given in Eq.(32), we can minimize it with respect of the particle positions $\vec{x}$ by stating that

$$\frac{\partial U(\vec{x}, \vec{n})}{\partial \vec{x}} = 0. \quad (44)$$

This condition delivers the following relation between the minimizing vector $\vec{x}_0$ and $\vec{n}$

$$\vec{x}_0 = s \mathcal{A}^{-1} \left( \frac{k_0}{k_m} \frac{d}{s} \vec{n} - \vec{v} \right). \quad (45)$$

Now, we can substitute this value $\vec{x}_0$ in $U(\vec{x}, \vec{n})$, by obtaining

$$\begin{aligned} U(\vec{x}_0, \vec{n}) &= \frac{1}{2} k_m \vec{x}_0 \cdot \mathcal{A} \vec{x}_0 + s k_m \vec{v} \cdot \vec{x}_0 - k_0 d \vec{n} \cdot \vec{x}_0 + C \\ &= \frac{1}{2} k_0 d^2 \vec{n} \cdot \left( \mathcal{I} - \frac{k_0}{k_m} \mathcal{A}^{-1} \right) \vec{n} + k_0 s d \vec{n} \cdot \mathcal{A}^{-1} \vec{v} + \dots, \end{aligned} \quad (46)$$

where we have omitted all terms not depending on $\vec{n}$. We can note that the $\vec{n}$-dependent part of $U(\vec{x}_0, \vec{n})$ in Eq.(46) represents the exponent of the terms summed in Eq.(39) (or within the Riemann theta function), except for the Boltzmann factor $-1/(K_B T)$. This means that everything goes as if the Boltzmann statistic were applied only to the discrete spin variables after minimizing the energy with respect to the continuous coordinate variables. This is possible only because all interactions are based on quadratic energies and therefore on a linear response with respect to the applied forces. The nonlinearity of the system emerges only thanks to the jumps of the spin variables, which make the overall potential energy non-convex.
Similarly to the one-particle case, also for the N-particle chain we can obtain a second mathematical form for $Z_H$. To this aim, it is possible to introduce the following reciprocal formula or modular transformation (sometimes also called Jacobi identity) for the Riemann theta function [112–114]

$$ \Theta(\vec{z}|\Omega) = \frac{1}{\sqrt{\det(-i\Omega)}} e^{-\pi i \vec{z} \cdot \Omega^{-1} \vec{z}} \Theta(\Omega^{-1} \vec{z}| - \Omega^{-1}), \quad (47) $$

which is valid for any $\Omega \in \mathcal{H}_N$ (see Appendix for details). Here, the square root assumes its principal value. This property can be applied to Eq. (43) and allows to obtain a second form for the partition function of our system, as follows

$$ Z_H = \frac{(2\pi K_B T)^N}{\sqrt{k_m^N k_0^N d^{2N} \det \mathcal{B}}} e^{\frac{k_m s^2}{2K_B T} \vec{v} \cdot \mathcal{B}^{-1} \vec{v}} e^{-\frac{k_m s^2}{2K_B T} (N-1)} \\ \times e^{-\frac{k_1 \ell_1^2}{2K_B T}} e^{-\frac{k_N \ell_N^2}{2K_B T}} \Theta \left( \frac{s}{d} \mathcal{B}^{-1} \vec{v} \left| \frac{2\pi K_B T}{k_0 d^2} i \mathcal{A} \mathcal{B}^{-1} \right. \right), \quad (48) $$

where $\mathcal{B} = \mathcal{A} - \frac{k_0}{k_m} \mathcal{I}$. The possibility of having two different forms for the partition function is useful from a numerical point of view. In fact, if we look at the structure of the theta function, we observe that its convergence is faster if all the eigenvalues of $2\pi\Im(\Omega)$ are much larger than one [113]. As a matter of fact, if this condition is verified, the Gaussian shape summed within the theta function is very narrow and therefore only a few values of $\vec{n}$ are needed to find the whole sum accurately. So, when $2\pi\Im(\Omega)$, defined in Eq. (41), has large eigenvalues (larger than one) it is better to use Eq. (43), and when $2\pi\Im(\Omega)$ has small eigenvalues (smaller than one) it is better to use Eq. (48). We have observed numerically that frequently (but not always) the matrix $2\pi\Im(\Omega)$ is well-conditioned in the sense that the eigenvalues are typically either all larger than one or all between zero and one. We can also affirm that Eq. (43) is a partition function expression rapidly convergent in the low temperature regime whereas Eq. (48) is more adapted in the high temperature regime (always with $K_B T \ll W$). Anyway, numerical algorithms for the calculation of the Riemann theta function can be found in the literature [115].

The knowledge of the partition function in both forms given in Eqs. (43) and (48) allows the determination of the force applied by the device to hold the chain in a given position. In the most general case, we fix the left device at position $\ell_1$ and the right device at position $\ell_2$. Consequently, the two applied forces can be determined as

$$ \langle F_1 \rangle = -K_B T \frac{\partial \log Z_H(\ell_1, \ell_N)}{\partial \ell_1}, \quad (49) $$

$$ \langle F_N \rangle = -K_B T \frac{\partial \log Z_H(\ell_1, \ell_N)}{\partial \ell_N}, \quad (50) $$
18 Statistical mechanics of rate-independent stick-slip

$-K_B T \log Z_H(\ell_1, \ell_N)$ being the Helmholtz free energy of the system. Two particular cases are of interest: (i) if the first extremity is left free, which means that $k_1 = 0$, then the partition function depends only on $\ell_N$ and the right force reads

$$ \langle F_s^{(i)} \rangle = -K_B T \frac{\partial \log Z_H(\ell)}{\partial \ell}, \qquad (51) $$

where we have considered $\ell_N = \ell$ to simplify the notation; (ii) if the two extremities are rigidly connected, we can write $\ell_1 = \ell - (N-1)s$ and $\ell_N = \ell$ [which means that the distance $\ell_N - \ell_1$ is maintained equal to the equilibrium value $(N-1)s$] and the total force is

$$ \langle F_s^{(ii)} \rangle = -K_B T \frac{\partial \log Z_H(\ell - (N-1)s, \ell)}{\partial \ell}, \qquad (52) $$

representing the sum of the two left and right forces. Examples of application of Eqs.(51) and (52) are given below.

To determine the average value $\langle \vec{n} \rangle$ of the spin variables vector, we can simply use the expression of the partition function given in Eq.(39), and we eventually get

$$ \langle \vec{n} \rangle = \frac{\sum_{\vec{n} \in \mathbb{Z}^N} \vec{n} e^{-\frac{k_0 d^2}{2K_B T} \vec{n} \cdot (\mathcal{I} - \frac{k_0}{k_m} \mathcal{A}^{-1})} \vec{n} e^{-\frac{k_0 s d}{K_B T} \vec{n} \cdot \mathcal{A}^{-1} \vec{v}}}{\sum_{\vec{n} \in \mathbb{Z}^N} e^{-\frac{k_0 d^2}{2K_B T} \vec{n} \cdot (\mathcal{I} - \frac{k_0}{k_m} \mathcal{A}^{-1})} \vec{n} e^{-\frac{k_0 s d}{K_B T} \vec{n} \cdot \mathcal{A}^{-1} \vec{v}}}. \qquad (53) $$

This expression can be written in terms of the theta function and its derivatives as follows. First of all, we can write the gradient of the Riemann theta function with respect of the vector variable $\vec{z}$ in the form

$$ \vec{\nabla} \Theta(\vec{z} | \Omega) = \frac{\partial \Theta(\vec{z} | \Omega)}{\partial \vec{z}} = 2\pi i \sum_{\vec{n} \in \mathbb{Z}^N} \vec{n} e^{2\pi i (\frac{1}{2} \vec{n} \cdot \Omega \vec{n} + \vec{n} \cdot \vec{z})}. \qquad (54) $$

Then, we can rewrite Eq.(53) as

$$ \langle \vec{n} \rangle = \frac{\vec{\nabla} \Theta \left( \frac{k_0 s d i}{2 \pi K_B T} \mathcal{A}^{-1} \vec{v} \middle| \frac{k_0 d^2 i}{2 \pi K_B T} \left( \mathcal{I} - \frac{k_0}{k_m} \mathcal{A}^{-1} \right) \right)}{2 \pi i \Theta \left( \frac{k_0 s d i}{2 \pi K_B T} \mathcal{A}^{-1} \vec{v} \middle| \frac{k_0 d^2 i}{2 \pi K_B T} \left( \mathcal{I} - \frac{k_0}{k_m} \mathcal{A}^{-1} \right) \right)}, \qquad (55) $$

which is the multi-dimensional generalization of Eq.(14). Unfortunately, there are no explicit formulas for the logarithmic derivatives of the Riemann theta function and therefore we cannot obtain more explicit expressions for $\langle \vec{n} \rangle$ or other quantities as we had done for a single particle system. By deriving
Eq.(38) with respect to $\vec{b}$, we obtain the integral expression

$$ \int_{\mathbb{R}^N} e^{-\frac{1}{2}\vec{y}\cdot\mathcal{N}\vec{y}} e^{\vec{b}\cdot\vec{y}} \vec{y} d\vec{y} = \sqrt{\frac{(2\pi)^N}{\det \mathcal{N}}} e^{\frac{1}{2}\vec{b}\cdot\mathcal{N}^{-1}\vec{b}} \mathcal{N}^{-1}\vec{b}, \quad (56) $$

which can be used to determine the average value $\langle \vec{x} \rangle$ of the particles position vector. Indeed, by definition we have

$$ \langle \vec{x} \rangle = \frac{\sum_{\vec{n} \in \mathbb{Z}^N} e^{-\frac{C}{K_B T}} \int_{\mathbb{R}^N} \vec{x} e^{-\frac{k_m}{2K_B T} \vec{x} \cdot \vec{A} \vec{x}} e^{-\frac{k_0 d \vec{n} \cdot \vec{x} - s k_m \vec{v} \cdot \vec{x}}{K_B T}} d\vec{x}}{\sum_{\vec{n} \in \mathbb{Z}^N} e^{-\frac{C}{K_B T}} \int_{\mathbb{R}^N} e^{-\frac{k_m}{2K_B T} \vec{x} \cdot \vec{A} \vec{x}} e^{-\frac{k_0 d \vec{n} \cdot \vec{x} - s k_m \vec{v} \cdot \vec{x}}{K_B T}} d\vec{x}}, \quad (57) $$

and by identifying $\mathcal{N}$ with $k_m A/(K_B T)$ and $\vec{b}$ with $(k_0 d \vec{n} - s k_m \vec{v})/(K_B T)$, we easily obtain the simple relation

$$ \langle \vec{x} \rangle = s \mathcal{A}^{-1} \left( \frac{k_0}{k_m} \frac{d}{s} \langle \vec{n} \rangle - \vec{v} \right). \quad (58) $$

Remarkably, it means that the knowledge of $\langle \vec{n} \rangle$ is sufficient to determine the vector $\langle \vec{x} \rangle$. Moreover, it is also interesting to remark that Eq.(58) has exactly the same mathematical form of Eq.(45). It means that the minimizing of the total energy with respect to $\vec{x}$ delivers correct results even if the thermal fluctuations act on the system.

In order to evaluate the effect of the substrate on the sliding chain, we can define this sort of adimensional order parameter

$$ \Psi = \frac{1}{d^2} \sum_{j=1}^{N-1} (\langle x_{j+1} \rangle - \langle x_j \rangle - s)^2, \quad (59) $$

which assumes the value zero if the chain is at equilibrium and larger values for an heterogeneous non-equilibrium spacing. Of course, $\Psi$ depends on $\ell$ during the sliding of the chain and therefore we can also introduce its mean value over a period of the corrugated substrate

$$ \langle \Psi \rangle = \frac{1}{d} \int_{0}^{d} \Psi(\ell) d\ell = \frac{1}{d^3} \sum_{j=1}^{N-1} \int_{0}^{d} (\langle x_{j+1} \rangle - \langle x_j \rangle - s)^2 d\ell, \quad (60) $$

which is useful to have an overall quantification of the trade-off between chain and substrate energies.

# 5 Results for the modified Frenkel-Kontorova model

The theoretical results previously shown have been numerically implemented in order to explore the behavior of the thermally modified Frenkel-Kontorova
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Static force $\langle F_s \rangle / (k_0 d)$, average particles positions $\langle \vec{x} \rangle / d$ and order parameter $\Psi$ for a chain with $N = 4$ (left column) and $N = 5$ (right column). Both the cases (i) with the device acting only on the second chain extremity (red curves), and (ii) with the device that keeps the distance between the ends fixed at $s(N-1)$ (blue curves) are considered. The adimensional physical quantities are plotted versus the sliding ratio $\ell/d$. In the panels showing $\langle \vec{x} \rangle / d$, the shaded areas are useful to distinguish the substrate wells. We adopted the parameters $K_B T / (k_0 d^2) = 0.002$, $k_m / k_0 = 0.5$, $s/d = 0.39$. Moreover, we imposed $k_1 / k_0 = 0$ and $k_N / k_0 = 1$ for the case (i) and $k_1 / k_0 = k_N / k_0 = 1$ for the case (ii).</figcaption></figure>

model. All examples have been studied under the assumption $W \gg K_B T$ to ensure that we stay within the validity of the theory. We explore here the following issues: the Frenkel-Kontorova stick-slip with sliding boundary devices, the effect of thermal fluctuations, the effect of the elastic contrast between chain and substrate, the effect of the boundary devices intrinsic elasticity and, finally, the effect of the geometric contrast between chain and substrate.

## 5.1 Frenkel-Kontorova stick-slip with sliding boundary devices

A first example can be found in Fig.4, where we show the evolution of the static friction force, the particles position and the order parameter as function of the position $\ell$ of the device for the two cases (i) and (ii) defined above. It means that we considered the case (i) with the device acting only on the right of the chain and the case (ii) where the device acts on both extremities, keeping
their distance fixed at $s(N - 1)$. In both cases, the device spans one period of the corrugated substrate and the previously mentioned physical quantities are recorded for two chains of length $N = 4$ and $N = 5$. In this first example, all parameters are fixed as declared in the caption of Fig.4. We observe that in both cases (i) and (ii), to keep the device at $\ell/d = 0$, we have to apply a given force since this configuration is not at equilibrium. Of course, this issue depends in a complex way on the ratios $s/d$ and $k_m/k_0$. Indeed, we can observe that there are some values of $\ell/d$ for which the force is zero, and then they correspond to the equilibrium configurations. Anyway, in the curves of $\langle F_s \rangle/(k_0d)$, we can see some jumps corresponding to the transitions of the particles between adjacent substrate wells. We adopted a low temperature and a low ratio $k_m/k_0$ to amplify and then easily identify this effect. This behavior is consistent with the evolution of $\langle \vec{x} \rangle/d$, where we can observe the transitions of each single particle (stick-slip processes). It is interesting to see that a transition of one particle between two adjacent substrate wells may induce a step behavior (smoothed by the temperature) in the position of other particles not transiting between the wells. This is due to the strong energetic coupling of all variables within the Frenkel-Kontorova model. Furthermore, we can observe that at transitions there is a significant decrease in the order parameter $\Psi$. In fact, it represents a measure of the equilibrium state of the chain and, therefore, must decrease as a result of the transitions that occur due to the geometric and energetic compromise between chain and substrate. In other words, transitions occur because they are energetically favored and thus bring the system into a configuration closer to equilibrium, with a lower $\Psi$. From Fig.4 we can also deduce that the friction force patterns are more complicated in the case of the boundary condition (ii) since this configuration is more constraining for the system, thus generating a more complex sequence of transitions.

We are studying the geometric and energetic competition between a mass-spring chain and a corrugated substrate. In this context, Aubry proved an interesting transition originated by this competition in the Frenkel-Kontorova model with an infinite ($N \to \infty$) and incommensurate ($s/d$ irrational) chain without temperature effects [64–66]. For values of the ratio $k_0/k_m$ below a certain threshold, chain elasticity prevails over the energy barrier of the substrate and we observe an unpinned sliding phase. In this case the static friction vanishes and we say that the sliding is superlubric (an infinitesimal force can move the chain). We remark that even if the static friction is zero, the kinetic friction is nonzero since the positive speed induces the excitation of phonons in the chain, requiring a certain amount of energy. On the other hand, above the critical value of $k_0/k_m$ the particles are pinned by the corrugated substrate since the substrate energy is predominant over the chain elastic energy. In this case, we have a non-zero static friction to move the chain since it must be able to overcome the Peierls-Nabarro barrier (and in this case the static friction force is analogous to the Peierls-Nabarro stress) [61–63]. The threshold or critical value of $k_0/k_m$ that discriminates unpinned and pinned states
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Hull functions for the chains considered in Fig.4. We considered the chain lengths $N = 4$ (first column) and $N = 5$ (second column) and the two boundary conditions (i) (first row) and (ii) (second row). All parameters are indicated in the caption of Fig.4.</figcaption></figure>

depends in a very intricate way on the incommensurability ratio $s/d$ [67]. In particular, it reaches its highest value when $s/d$ is equal to the golden ratio $(\sqrt{5}-1)/2$ [59]. The shift between the unpinned sliding phase to the pinned condition with stick-slip motion generates the so-called analyticity breaking of the particle positions. It means that the curves generated by plotting the actual particle position $x_j$ versus the unperturbed position $x_{j,0}$ (without substrate) exhibit discontinuities or gaps, induced by the jump of the particle between substrate wells. In general, the relation between $\langle \vec{x} \rangle/d$ and $\langle \vec{x}_0 \rangle/d$ (with $\ell$ variable) is called hull function and shows discontinuities in the pinned Aubry state. Then, in this pinned state the hull function analyticity breaking, the crossing of the Peierls-Nabarro barrier, and the stick-slip motion are three different aspects of the same phenomenology. Although the Aubry transition exists only for infinite and incommensurable chains at $T=0$, it is interesting to observe the shape of the hull functions for our system as well. In our case, the reference positions are given by $x_{j,0} = \ell - (N-j)s \ \forall j = 1, ..., N$ where $\ell$ is the variable position of the sliding device. Of course, these references positions can be used for both boundary conditions (i) and (ii), previously introduced. In Fig.5, one can find the hull functions for the chains considered in Fig.4. We considered the chain lengths $N = 4$ (first column) and $N = 5$ (second column) and the two boundary conditions (i) (first row) and (ii) (second row). We can easily recognize the jumps in the hull functions but they are not real discontinuities because of the effect of the thermal fluctuations that are able to smooth out their shape. We observe that our finite-size finite-temperature Frenkel-Kontorova model is always statically pinned, even for an irrational value of $s/d$ since we always have a positive static friction corresponding to a positive
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Temperature behavior of the static force $\langle F_s \rangle / (k_0 d)$, the average particles positions $\langle \bar{x} \rangle / d$ and the order parameter $\Psi$ for a chain with $N = 4$ with both boundary conditions (i) (left column) and (ii) (right column). The adimensional physical quantities are plotted versus the sliding ratio $\ell/d$ and parametrized by the temperature $K_B T / (k_0 d^2) = 0.001, 0.003, 0.005, 0.007, 0.009, 0.011$. We adopted the parameters $k_m / k_0 = 0.5$ and $s/d = 0.39$. Moreover, we imposed $k_1 / k_0 = 0$ and $k_N / k_0 = 1$ for the case (i) and $k_1 / k_0 = k_N / k_0 = 1$ for the case (ii).</figcaption></figure>

Peierls-Nabarro barrier. It is also interesting to remember that an Aubry-like transition, based a symmetry-breaking phenomenon, has been identified for finite chain at zero temperature [68–70]. However, in our finite chain at finite temperature all transitions are smooth as one can see in Fig.5, and the response is controlled by the compromise between pinning intensity and thermal fluctuations that can facilitate the chain sliding. Of course, the slope of the jumps in our hull functions is larger for lower temperatures. Besides, from Fig.5 we deduce that the jumps of the transitions are larger for the boundary condition (i) since the second condition (ii) is more constraining, limiting the positions discontinuity. However, under condition (ii), the height of the jumps is smaller but the pattern of transitions is more intricate, as seen in Fig.4. Interestingly, Aubry-type transitions have been experimentally observed in different finite systems, as discussed in the recent literature for one-dimensional [71–74] and two-dimensional [75] structures.
<figure><img src="image_8.png" /><figcaption>**Fig. 8** Influence of $k_m$ and $T$ on the behavior of the maximum static friction force $\max \{|\langle F_s \rangle/(k_0d)|, 0 < \ell/d \le 1\}$ and on the average order parameter $\langle \Psi \rangle$ for a chain with $N=4$ with both boundary conditions (i) (first column) and (ii) (second column). We adopted the parameters $s/d = 0.39$ and $K_B T/(k_0 d^2) = 0.001, 0.006, 0.01, 0.015, 0.02$. Moreover, we imposed $k_1/k_0 = 0$ and $k_N/k_0 = 1$ for the case (i) and $k_1/k_0 = k_N/k_0 = 1$ for the case (ii).</figcaption></figure>

the temperature fosters these transitions allowing a more complete exploration of the phase space, even where the total energy is higher. This behavior is directly reflected in the graph of $\Psi$, which decreases significantly as the temperature increases. In particular, the minimum points of the $\Psi$ function are in correspondence of the transitions, occurring as soon as they are energetically favored and, therefore, bringing the system closer to equilibrium. As the chain is already closer to equilibrium as a result of these transitions, temperature has a lesser effect in reducing the value of $\Psi$ at these points. This can be seen in the graphs of $\Psi$ versus $\ell/d$ in Fig. 6, where, at the transitions the curves at different temperatures are very closed together as opposed to the regions without transitions where the curves are more spaced.

## 5.3 Effect of the elastic contrast between chain and substrate

The effects induced by varying the chain elastic constant $k_m$ can be seen in Fig. 7. First of all, from the plots of $\langle F_s \rangle/(k_0d)$ versus $\ell/d$, we deduce that the static friction force rises as the elastic constant $k_m$ increases. This behavior is similar to the one obtained for a single particle in Fig. 2. However, in the present case with $N=4$ particles, we also observe that the transitions occur at positions that depend significantly on the elastic constant $k_m$ [for both boundary conditions (i) and (ii)]. This is easily explained by the fact that $k_m$ controls the forces in the system and thus modifies the energy competition between chain and substrate. The transitions shift with $k_m$ can be clearly seen in both plots of $\langle \vec{x} \rangle/d$ versus $\ell/d$ and $\Psi$ versus $\ell/d$. The latter was represented
<figure><img src="image_9.png" /><figcaption>**Fig. 9** Effect of the boundary devices on the transitions pattern: we varied the elastic constant $k_N$ in a system with $N = 4$ under of the condition (i). The static friction force $\langle F_s \rangle / (k_0 d)$, the average particles positions $\langle \vec{x} \rangle / d$ and the order parameter $\Psi$ are plotted versus $\ell/d$. The hull function $\langle \vec{x} \rangle / d$ versus $\langle \vec{x}_0 \rangle / d$ is shown as well. The dashed ellipses show the emergence of a transition for the largest values of $k_N$. We adopted the parameters $s/d = 0.39$, $k_m/k_0 = 1$ and $K_B T/(k_0 d^2) = 0.001$. Moreover, we imposed $k_1/k_0 = 0$ and $k_N/k_0 = 0.5, 2, 8, 32$ corresponding to the condition (i).</figcaption></figure>

on a semi-logarithmic scale because the effects of varying $k_m$ on $\Psi$ are very pronounced. It is in fact evident that a strong increase in $k_m$ brings the system very close to the equilibrium configuration with the particles equispaced at distance s.

The effect of varying the temperature $T$ and the elastic constant $k_m$ can be summarized by means of the plots in Fig.8. Here, one can find the influence of $k_m$ and $T$ on the behavior of the maximum static friction force $\max \{|\langle F_s \rangle / (k_0 d)|, 0 < \ell/d \le 1\}$ and on the average order parameter $\langle \Psi \rangle$ for a chain with $N = 4$ with both boundary conditions (i) (first column) and (ii) (second column). In this analysis we considered a fixed value of $s$ but we spanned an entire period of the corrugated substrate by varying $\ell$ in the interval $(0, d)$. This allowed us to determine the maximum friction force and the average value of the order parameter, as defined in Eq.(60). The final result is that the maximum of friction force $\langle F_s \rangle / (k_0 d)$ increases with $k_m/k_0$ and decreases with $K_B T/(k_0 d^2)$ (for each value of $N$ and for both boundary conditions (i) and (ii)). Similarly, the average order parameter $\langle \Psi \rangle$ decreases with both $k_m/k_0$ and $K_B T/(k_0 d^2)$. This analysis is useful when $k_0$ is constant and $k_m$ variable. The dual situation with $k_0$ variable and $k_m$ constant is not represented in Fig.8 but, as discussed for the system with one particle, its overall behavior can be summarized as follows. The maximum value of $\langle F_s \rangle / (k_m d)$ is decreasing with an increasing temperature $K_B T/(kd^2)$ and is increasing with an increasing ratio $k_0/k_m$. At the same time, $\langle \Psi \rangle$ decreases with $K_B T/(k_m d^2)$ and increases with $k_0/k_m$.
## 5.4 Effect of the boundary devices intrinsic elasticity

To further show the complexity of the system under investigation, we show in Fig.9 the effect of the boundary devices on the transitions pattern. We considered a chain with length $N = 4$ under boundary condition (i), i.e. with a single device pulling the chain from the right extremity. This device is here modified by changing its elastic constant $k_N$. It is interesting to observe that for the first two values adopted $k_N/k_0 = 0.5, 2$ we have no transitions for $\ell/d < 0.5$ while for the other values $k_N/k_0 = 8, 32$ we observe the emergence of one transition around $\ell/d \approx 0.3$. This transition origination can be easily identified in all panels of Fig.9 (see dashed ellipses). We give a direct interpretation of this result as follows. The initial positions of the ideal chain particles for $\ell/d = 0$ would be $x_1/d \approx -1.2$, $x_2/d \approx -0.8$, $x_3/d \approx -0.4$ and $x_4/d \approx 0$, without interaction with the corrugated substrate. Hence, at the beginning, the first two particles are within the well with $n = -1 (-3/2 < x/d < -1/2)$, and the others in the well with $n = 0 (-1/2 < x/d < 1/2)$. So, in an ideal situation, we should expect a first transition for $\ell/d \approx 0.3$, when $x_2/d \approx -0.5$, i.e. when the second particle reaches the transition threshold between the two neighboring wells. However, this transition is not observed for low values of $k_N$ since in this case the device, instead of generating chain sliding, generates its own elongation due to its weak elasticity. In this case the chain pinning is stronger than the device pulling. Instead, for higher values of $k_N$, the device deforms much less and the chain has to move, generating the transition of the second particle when $\ell/d \approx 0.3$. This example explains that there is not only a geometric and energetic trade-off between chain and substrate but also a competition between device elasticity and chain elasticity. This brings additional complexity to the friction phenomenon we are studying.

## 5.5 Effect of the geometric contrast between chain and substrate

Up to this point, we have always considered the $s/d$ ratio to be constant and therefore assumed that we were working with the same geometric competition between chain and substrate. Now it is interesting to explore the effects of the ratio $s/d$ on the stick-slip process and on the static friction. The results of this analysis can be found in Figs.10, 11, 12, and 13, where we considered chains of length $N = 2, N = 3, N = 4$ and $N = 5$, respectively. In each of these figures, we plotted the static friction surface $(F_s)/(k_0d)$ as function of $\ell/d$ and $s/d$ and the corresponding contour lines, for three different values of the chain to substrate ratio $k_m/k_0$ of elastic constants (with fixed temperature and fixed device elastic constant). For the sake of brevity, we only considered the boundary condition (i) in this analysis. These figures show that the friction pattern over the period $0 < \ell < d$ is strongly sensible to the variations of $s/d$ and this influence is stronger for large values of $k_m/k_0$. Obviously the pattern of the friction response is increasingly complicated as $N$ and the ratio $k_m/k_0$ increase. For a given value of $s/d$ the transitions observable in the
<figure><img src="image_10.png" /><figcaption>**Fig. 10** Static friction $\langle F_s \rangle / (k_0 d)$ as function of $\ell/d$ and $s/d$ for $N = 2$. We represented the surface and the contour lines of $\langle F_s \rangle / (k_0 d)$ for $k_m/k_0 = 0.5$ (first column), $k_m/k_0 = 1$ (second column), and $k_m/k_0 = 10$ (third column). We adopted the parameters $K_B T / (k_0 d^2) = 0.01$ and we analyzed only the case (i) with $k_1/k_0 = 0$ and $k_N/k_0 = 1$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 11** Static friction $\langle F_s \rangle / (k_0 d)$ as function of $\ell/d$ and $s/d$ for $N = 3$. We represented the surface and the contour lines of $\langle F_s \rangle / (k_0 d)$ for $k_m/k_0 = 0.5$ (first column), $k_m/k_0 = 1$ (second column), and $k_m/k_0 = 10$ (third column). We adopted the parameters $K_B T / (k_0 d^2) = 0.01$ and we analyzed only the case (i) with $k_1/k_0 = 0$ and $k_N/k_0 = 1$.</figcaption></figure>

curve $\langle F_s \rangle / (k_0 d)$ versus $\ell/d$ (the friction pattern) can be identified by tracing straight lines parallel to the $\ell/d$ axis on the plane with axes $s/d$ and $\ell/d$. This is especially evident when considering the plane where the contour lines are drawn. In fact, the highly visible height jumps noticeable in this plane are precisely the transitions in the friction responses.
<figure><img src="image_12.png" /><figcaption>**Fig. 12** Static friction $\langle F_s \rangle/(k_0 d)$ as function of $\ell/d$ and $s/d$ for $N = 4$. We represented the surface and the contour lines of $\langle F_s \rangle/(k_0 d)$ for $k_m/k_0 = 0.5$ (first column), $k_m/k_0 = 1$ (second column), and $k_m/k_0 = 10$ (third column). We adopted the parameters $K_B T/(k_0 d^2) = 0.01$ and we analyzed only the case (i) with $k_1/k_0 = 0$ and $k_N/k_0 = 1$.</figcaption></figure>

The results shown in Figs.10, 11, 12, and 13 concerns only one value of temperature. We performed the same analysis by changing the temperature over the interval $K_B T/(k_0 d^2) = 0.01j, \forall j = 1, ..., 10$, and we collected the outcomes in order to show the maximum static friction force $\max\{|\langle F_s \rangle/(k_0 d)|, 0 < \ell/d \le 1\}$ as function of $s/d$ and parametrized by the temperature. The result of this operation can be found in Fig.14, where we considered the chains of length $N = 2, N = 3, N = 4$ and $N = 5$ in the first, second, third and fourth row respectively. Moreover, we used the chain to substrate elastic constants ratio $k_m/k_0 = 1$ in first column, $k_m/k_0 = 10$ in the second column, and $k_m/k_0 = 50$ in the third column. The collection of these graphs shows well the effect of geometric compromise between chain and substrate as the temperature increases. In fact, the first thing we can observe is that the maximum static friction is always decreasing with temperature regardless of the structural parameter $s/d$. As far as the geometric ratio $s/d$ is concerned, we note that there are particular values of this parameter that significantly reduce the maximum friction force when $k_m$ is sufficiently large. For $N = 2$ this happens when $s/d = 1/2$, for $N = 3$ when $s/d = 1/3$ and $2/3$, for $N = 4$ when $s/d = 1/4, 1/2$, and $3/4$, and finally for $N = 5$ when $s/d = 1/5, 2/5, 3/5$, and $4/5$. On the contrary, for $s/d = 0$ and $s/d = 1$ the friction is always higher than for all other values of the structural ratio $s/d$. In general, we can summarize this result by means of the two following properties: (I) if $s/d = 0, 1 \forall N$ then the friction assume the value $\max_{0<s/d \le 1} \max_{0<\ell/d \le 1} |\langle F_s \rangle/(k_0 d)|$; (II) if $s/d = j/N \forall j = 1, ..., N - 1$ then the friction assume the value $\min_{0<s/d \le 1} \max_{0<\ell/d \le 1} |\langle F_s \rangle/(k_0 d)|$.

We can explain this behavior as follows. We suppose to work with $k_m$ and $k_N$ sufficiently large and therefore we consider here a sliding finite chain at
<figure><img src="image_13.png" /><figcaption>Fig. 13 Static friction $\langle F_s \rangle/(k_0 d)$ as function of $\ell/d$ and $s/d$ for $N = 5$. We represented the surface and the contour lines of $\langle F_s \rangle/(k_0 d)$ for $k_m/k_0 = 0.5$ (first column), $k_m/k_0 = 1$ (second column), and $k_m/k_0 = 10$ (third column). We adopted the parameters $K_B T/(k_0 d^2) = 0.01$ and we analyzed only the case (i) with $k_1/k_0 = 0$ and $k_N/k_0 = 1$.</figcaption></figure>

equilibrium, corresponding to the positions $x_j = \ell - s(N - j) \forall j = 1, \dots, N$. It means that we introduce the geometrical competition between chain and substrate but we neglect the corresponding energetic trade-off. Similarly, we also neglect the thermal fluctuations. We also consider in this calculation the original Frenkel-Kontorova potential stated in Eq. (29). So, we can write

$$U_{FK} = -\frac{W}{2} \sum_{j=1}^{N} \cos \left( \frac{2\pi\ell}{d} - \frac{2\pi sN}{d} + \frac{2\pi sj}{d} \right) = -\frac{W}{2} \sum_{j=1}^{N} \Re \{ e^{i\eta} e^{-i\xi N} e^{i\xi j} \}, \quad (61)$$

where we defined $\eta = 2\pi\ell/d$ and $\xi = 2\pi s/d$. Now, we can easily develop this expression as follows

$$U_{FK} = -\frac{W}{2} \Re \left\{ e^{i\eta} e^{-i\xi N} \sum_{j=1}^{N} e^{i\xi j} \right\} = -\frac{W}{2} \Re \left\{ e^{i\eta} \frac{e^{-i\xi(N-1)} - e^{i\xi}}{1 - e^{i\xi}} \right\}, (62)$$

where we used the sum $\sum_{j=1}^{N} x^j = (x - x^{N+1})/(1 - x)$.

We can now multiply numerator and denominator of the internal fraction by $1 - e^{-i\xi}$ and we get

$$\begin{aligned} U_{FK} &= -\frac{W}{4} \frac{\Re \{ e^{i\eta} [e^{-i\xi(N-1)} - e^{i\xi}] (1 - e^{-i\xi}) \}}{1 - \cos \xi} \\ &= -\frac{W}{4} \left\{ \frac{\cos[\eta - (N-1)\xi] - \cos(\eta + \xi)}{1 - \cos \xi} - \frac{\cos(\eta - N\xi) - \cos \eta}{1 - \cos \xi} \right\}. (63) \end{aligned}$$
<figure><img src="image_14.png" /><figcaption>**Fig. 14** Behavior of the maximum static friction force max {$|\langle F_s \rangle/(k_0d)|, 0 < l/d \le 1$} of the Frenkel-Kontorova modified model as function of $N$, $k_m/k_0$, $s/d$ and the temperature $T$. We consider the chains of length $N = 2, N = 3, N = 4$ and $N = 5$ in the first, second, third and fourth row respectively. Moreover, we show the results for $k_m/k_0 = 1$ (first column), $k_m/k_0 = 10$ (second column), and $k_m/k_0 = 50$ (third column). All curved are plotted versus $s/d$ and parametrized by the temperature $K_B T/(k_0 d^2) = 0.01j, \forall j = 1, \dots, 10$. We also imposed $k_1/k_0 = 0$ and $k_N/k_0 = 10$ for simulating the boundary condition (i).</figcaption></figure>

Here, we can use a prosthaphaeresis formula obtaining the simplified result

$$ U_{FK} = -\frac{W}{2} \frac{\sin \frac{N\xi}{2} \left[ \sin \left( \eta + \xi - \frac{N\xi}{2} \right) - \sin \left( \eta - \frac{N\xi}{2} \right) \right]}{1 - \cos \xi}. \quad (64) $$

To conclude, we can use another prosthaphaeresis formula to get the final result

$$ U_{FK} = -W \frac{\sin \frac{N\xi}{2} \sin \frac{\xi}{2}}{1 - \cos \xi} \cos \left( \eta + \frac{\xi}{2} - \frac{N\xi}{2} \right). \quad (65) $$

This expression represents the Frenkel-Kontorova potential energy of a non-deformable chain interacting with the corrugated substrate.

The force necessary to keep the chain at $l/d = \eta/(2/\pi)$ is given by $F_s = \partial U_{FK}/\partial l = (2\pi/d)\partial U_{FK}/\partial \eta$, or more explicitly by

$$ F_s = \frac{2\pi W}{d} \frac{\sin \frac{N\xi}{2} \sin \frac{\xi}{2}}{1 - \cos \xi} \sin \left( \eta + \frac{\xi}{2} - \frac{N\xi}{2} \right). \quad (66) $$
<figure><img src="image_15.png" /><figcaption>Fig. 15 Behavior of the quantity $\max_{\eta \in (0, 2\pi)} d|F_s|/(2\pi W)$ given in Eq. (67) versus $\xi/(2\pi) = s/d$ for different values of $N$. The zeros of this function exactly correspond to the minima observed in Fig. 14.</figcaption></figure>

The maximum value of this force over the period of the corrugated substrate is then simply obtained as

$$ \max_{\eta \in (0, 2\pi)} |F_s| = \frac{2\pi W}{d} \left| \frac{\sin \frac{N\xi}{2} \sin \frac{\xi}{2}}{1 - \cos \xi} \right| \quad (67) $$

This result perfectly explains the behavior of the friction force observed in Fig. 14. Indeed, the numerator of Eq. (67) vanishes when $N\xi/2 = j\pi$, $\forall j = 0, ..., N$, which means $s/d = \xi/(2\pi) = j/N$, $\forall j = 0, ..., N$. However, the two extremes values $s/d = 0$ and $s/d = 1$ must not be considered because they are zeros of the denominator and they in fact correspond to the maximum value of Eq. (67). This scenario exactly corresponds to the properties (I) and (II) stated above, as one can also see in Fig. 15. Finally, the particular values of $s/d$ that significantly reduce the maximum friction force are explained by the pure geometrical competition between chain and substrate. It is interesting to remark that for an infinite chain the Aubry sliding phase is observed with incommensurate structures while for a finite chain the maximum reduction of friction can be obtained for commensurate structures. This is not a contradiction since the incommensurate Aubry chain is elastically deformable whereas we considered here a rigid chain to explain the friction reduction. Moreover in our analysis of Fig. 14 we have considered the complete behavior of the system including elasticity and thermal fluctuations.

# 6 Conclusions

In this work we have considered the equilibrium statistical mechanics of stick-slip processes, with application to the nanoscale friction. In other words, we can say that we studied the effect of thermal fluctuations on the rate-independent nanoscale interactions between a slider and a corrugated substrate. We have conveniently modified the one-particle Prandtl-Tomlinson model and the N-particle Frenkel-Kontorova model so that we can easily calculate the partition
function of both systems. In particular, instead of considering the classical sinu-
soidal energy profile of the corrugated substrate, we introduced a sequence of
quadratic potentials, each of which represents a substrate well. It is important
to remark that the sinusoidal and the quadratic profiles have a qualitatively
similar but quantitatively different behavior. Indeed, both profiles are periodic
but their shapes are quite different. Therefore, the scheme with the sequence
of parabolic wells should be considered as an independent model with specific
applications. First of all, it is interesting because of the possibility of a com-
prehensive analytical development, and moreover, it is well suited to the study
of the propagation of dislocations and fractures in solids or the sliding of bio-
logical structures [12, 28]. Of course, the quadratic potentials are integrable
and then allows a rather simple analytic treatment of the problem. Follow-
ing the classical spin variables approach, we had to add to the phase space
a discrete variable for each particle useful to identify the well explored by
the particle itself. In this context, it is important to remark that the calcula-
tion of the partition function based on the spin approach assumes that for all
substrate well, all possible particle positions are considered (also outside the
well extension). This effectively corresponds to a multivalued energy function
(see the superpositions of the quadratic curves in Figs.1 and 3). As shown in
Refs.[91, 96], the effect of this approximation is statistically negligible since
the potential energy in the regions outside the actual well extension is sensibly
higher than the energy in the internal region. Of course, this is true if the bar-
rier W between two adjacent wells is sufficiently larger than K<sub>B</sub>T. On the one
hand, the application of this technique to the (modified) one-particle Prandtl-
Tomlinson model leads to a partition function which can be written in terms
of the third Jacobi theta function. Then, the application of the classical Jacobi
functional identity allows us to obtain the closed form Fourier series describ-
ing the static friction force experienced by the single particle sliding on the
corrugated periodic substrate through a stick-slip motion. Importantly, this
force depends on the temperature and therefore the obtained results explain
clearly the thermolubricity phenomenon. On the other hand the spin variable
method applied to the (modified) N-particle Frenkel-Kontorova model leads
to a partition function written in terms of the N-dimensional Riemann theta
function. It is interesting to underline that the Riemann theta function is able
to fully describe the complexity of the rate-independent stick-slip phenomenon
controlled by the geometric and energetic compromise between chain and sub-
strate. In fact, in this case, we have the geometric compromise induced by the
difference of the characteristic lattice lengths of chain and substrate, and the
energetic compromise induced by the difference between the elastic constants
of chain and substrate. We discussed the effect of the elastic constants and of
the temperature of the stick-slip and friction responses. Moreover, we plotted
the hull functions for the particles of the chain in order to compare the tran-
sitions between the wells in our system and in the infinite incommensurate
Aubry configuration. We have shown that in our case with finite chains and
nonzero temperature the transitions are smooth, in contrast to those observed
by Aubry for infinite incommensurable chains at zero temperature, where one finds the so-called analyticity breaking of the particle positions. To conclude we explored the effects of the geometric trade-off induced by the ratio $s/d$. In particular, we proved the existence for finite chains of specific commensurate values of the ratio $s/d$, which are able to strongly reduce the friction experienced by the slider. We investigated the effect of the temperature on this friction reduction. To conclude, we can say that our proposed method to take into account the temperature effect in friction phenomena can be adopted for several problems, including the application to the motion of defects in solid materials or the sliding of biological structures [12, 28]. The originality of this approach will allow the study of a number of problems that are for the moment only partially or not yet theoretically modeled. These perspectives concern the generalization of this approach to the two-dimensional case [120], to more complex corrugated substrates having for example two or more energy levels for the intercalated wells [121], to corrugated soft substrates with elastically moving wells [122], to the interplay between adhesion and friction with application to bio-systems [123], to the kinetic or dynamic case by means of the Langevin and/or Fokker-Planck methodologies [124, 125], and so forth. In particular, specific techniques for studying the dynamics of systems with multiwell energy landscapes will be expolited [106, 126, 127].

## Declarations

* The author certifies that he have no affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this manuscript.

* No funding was received to assist with the preparation of this manuscript.

* All data generated or analyzed during this study are included in this published article.

## Appendix A Reciprocal relation for Jacobi and Riemann theta functions

In this Appendix, we briefly discuss the conceptual relationship between direct and reciprocal Bravais lattices, Dirac comb, Poisson summation formula and reciprocal relation for the Riemann theta function and third Jacobi theta function. We consider a Bravais lattice composed of points $\vec{r} = n_1\vec{a}_1 + \dots + n_N\vec{a}_N \in \mathbb{R}^N$ where $n_j \in \mathbb{Z}$ and $\vec{a}_j \in \mathbb{R}^N$ are the primitive vectors. We introduce a periodic sampling function $f(\vec{x})$ on this Bravais lattice constructed by delta functions (Dirac comb)

$$f(\vec{x}) = \sum_{\vec{n} \in \mathbb{Z}^N} \delta(\vec{x} - n_1\vec{a}_1 - \dots - n_N\vec{a}_N). \qquad (A1)$$
To simplify the notation, we define a matrix $\mathcal{C} \in \mathcal{M}_{N,N}(\mathbb{R})$ where the columns represent the vectors $\vec{a}_j$. It means that $\mathcal{C} = [\vec{a}_1|...|\vec{a}_N]$ and we can write

$$f(\vec{x}) = \sum_{\vec{n} \in \mathbb{Z}^N} \delta(\vec{x} - \mathcal{C}\vec{n}), \qquad (A2)$$

with $\vec{n} = (n_1, ..., n_N)^T$. We assume that $\mathcal{C}$ is not singular. Then, if we define $\vec{\eta}$ such that $\vec{x} = \mathcal{C}\vec{\eta}$, we have

$$f(\vec{\eta}) = \sum_{\vec{n} \in \mathbb{Z}^N} \delta(\mathcal{C}(\vec{\eta} - \vec{n})). \qquad (A3)$$

Now, the function $f(\vec{\eta})$ is multi-periodic with period one along all directions $\eta_1, ..., \eta_N$. Hence, it can be developed in Fourier series

$$f(\vec{\eta}) = \sum_{\vec{m} \in \mathbb{Z}^N} c_{\vec{m}} e^{2\pi i \vec{m} \cdot \vec{\eta}}, \qquad (A4)$$

where the coefficients are given by

$$c_{\vec{m}} = \int_{[0,1]^N} f(\vec{\eta}) e^{-2\pi i \vec{m} \cdot \vec{\eta}} d\vec{\eta}. \qquad (A5)$$

In the set $[0, 1]^N$ we have $f(\vec{\eta}) = \delta(\mathcal{C}\vec{\eta})$ and therefore we get

$$c_{\vec{m}} = \int_{[0,1]^N} \delta(\mathcal{C}\vec{\eta}) e^{-2\pi i \vec{m} \cdot \vec{\eta}} d\vec{\eta} = \frac{1}{\det \mathcal{C}}, \qquad (A6)$$

as we can easily prove by applying the substitution $\vec{y} = \mathcal{C}\vec{\eta}$. Coming back to the variable $\vec{x}$, we finally obtain

$$f(\vec{x}) = \frac{1}{\det \mathcal{C}} \sum_{\vec{m} \in \mathbb{Z}^N} e^{2\pi i \vec{m} \cdot \vec{C}^{-1} \vec{x}}, \qquad (A7)$$

which is the Fourier series of the Dirac comb defined in Eq. (A1) or (A2). We can also determine the Fourier transform of the same function

$$F(\vec{k}) = \int_{\mathbb{R}^N} f(\vec{x}) e^{-i\vec{k} \cdot \vec{x}} d\vec{x} = \frac{1}{\det \mathcal{C}} \sum_{\vec{m} \in \mathbb{Z}^N} \int_{\mathbb{R}^N} e^{2\pi i \vec{m} \cdot \vec{C}^{-1} \vec{x}} e^{-i\vec{k} \cdot \vec{x}} d\vec{x}. \qquad (A8)$$

Since $\delta(\vec{q}) = \int_{\mathbb{R}^N} e^{i\vec{q} \cdot \vec{x}} d\vec{x}/(2\pi)^N$, we easily get

$$F(\vec{k}) = \frac{(2\pi)^N}{\det \mathcal{C}} \sum_{\vec{m} \in \mathbb{Z}^N} \delta(\vec{k} - 2\pi \mathcal{C}^{-T} \vec{m}). \qquad (A9)$$
It means that if $\vec{r} = C\vec{n}$ $\forall \vec{n} \in \mathbb{Z}^N$ is the direct Bravais lattice, then $\vec{k} = 2\pi C^{-T}\vec{m}$ $\forall \vec{m} \in \mathbb{Z}^N$ is the reciprocal Bravais lattice. In other words, the Dirac comb on the direct lattice is Fourier transformed into a Dirac comb on the reciprocal lattice. This property can be applied to find the so-called N-dimensional Poisson summation formula. We consider an arbitrary function $\phi(\vec{x})$ for which the Fourier transform $\Phi(\vec{k})$ exists. We define the replication or periodic summation

$$g(\vec{x}) = \sum_{\vec{n} \in \mathbb{Z}^N} \phi(\vec{x} - C\vec{n}) = \phi * f, \qquad (A10)$$

where $f$ is defined in Eq. (A2) and $*$ means convolution. Since the Fourier transform of the convolution is the product of the two Fourier transforms (convolution theorem), we have

$$G(\vec{k}) = \frac{(2\pi)^N}{\det C} \sum_{\vec{m} \in \mathbb{Z}^N} \Phi(\vec{k})\delta(\vec{k} - 2\pi C^{-T}\vec{m}) \qquad (A11) \\ = \frac{(2\pi)^N}{\det C} \sum_{\vec{m} \in \mathbb{Z}^N} \Phi(2\pi C^{-T}\vec{m})\delta(\vec{k} - 2\pi C^{-T}\vec{m}).$$

Now, we can remember that the Fourier transform of $e^{i\vec{x} \cdot \vec{v}}$ is given by $(2\pi)^N \delta(\vec{k} - \vec{v})$, and then from Eq. (A11) we come back to the original function $g(\vec{x})$, eventually obtaining

$$\sum_{\vec{n} \in \mathbb{Z}^N} \phi(\vec{x} - C\vec{n}) = \frac{1}{\det C} \sum_{\vec{m} \in \mathbb{Z}^N} \Phi(2\pi C^{-T}\vec{m})e^{2\pi i \vec{x} \cdot C^{-T}\vec{m}}, \qquad (A12)$$

which is the Fourier series of the periodic summation. For $\vec{x} = 0$, this result delivers the Poisson summation formula

$$\sum_{\vec{n} \in \mathbb{Z}^N} \phi(C\vec{n}) = \frac{1}{\det C} \sum_{\vec{m} \in \mathbb{Z}^N} \Phi(2\pi C^{-T}\vec{m}). \qquad (A13)$$

We take now into account the following particular function $\phi(\vec{x})$ with its Fourier transform $\Phi(\vec{k})$

$$\phi(\vec{x}) = e^{2\pi i (\frac{1}{2}\vec{x} \cdot \mathcal{T}\vec{x} + \vec{x} \cdot \vec{\sigma})}, \qquad (A14)$$

$$\Phi(\vec{k}) = \frac{1}{\sqrt{\det(-i\mathcal{T})}} e^{-\pi i (\vec{\sigma} - \frac{\vec{k}}{2\pi}) \cdot \mathcal{T}^{-1} (\vec{\sigma} - \frac{\vec{k}}{2\pi})}, \qquad (A15)$$

where $\mathcal{T} \in \mathcal{M}_{N,N}(\mathbb{C})$ with $\mathcal{T} = \mathcal{T}^T$, $\mathfrak{Im}(\mathcal{T}) > 0$ and $\vec{\sigma} \in \mathbb{C}^N$. To conclude, we substitute Eqs. (A14) and (A15) in the Poisson summation formula given in Eq. (A13) and we get, after straightforward calculations, the N-dimensional
Jacobi reciprocal relation for the Riemann theta function

$$ \Theta(\bar{z}|\Omega) = \sqrt{\frac{1}{\det(-i\Omega)}} e^{-\pi i \bar{z} \cdot \Omega^{-1} \bar{z}} \Theta(\Omega^{-1} \bar{z} | -\Omega^{-1}), \quad (A16) $$

where we have identified $\Omega = C^T T C$ and $\bar{z} = C^T \vec{\sigma}$. This proves Eq. (47) of the main text. Other more refined properties of $\Theta(\bar{z}|\Omega)$ can be found in the literature [112, 113]. The result given in Eq. (A16) can be specialized to the case with $N = 1$ by obtaining the original Jacobi identity

$$ \vartheta_3(z, \tau) = \frac{1}{\sqrt{-i\tau}} e^{\frac{z^2}{\pi i \tau}} \vartheta_3\left(\frac{z}{\tau}, -\frac{1}{\tau}\right), \quad (A17) $$

where we have conveniently compared the definitions of $\vartheta_3(z, \tau)$ and $\Theta(\bar{z}|\Omega)$ with $N = 1$ [112, 113]. This finally proves Eq. (17) of the main text as well.

## References

[1] J. Gao, W. D. Luedtke, D. Gourdon, M. Ruths, J. N. Israelachvili, and Uzi Landman, Frictional Forces and Amontons' Law: From the Molecular to the Macroscopic Scale, J. Phys. Chem. B 108, 3410-3425 (2004).

[2] M. Urbakh, J. Klafter, D. Gourdon, and J. Israelachvili, The nonlinear nature of friction, Nature 430, 525-528 (2004).

[3] A. Vanossi, N. Manini, M. Urbakh, S. Zapperi, and E. Tosatti, Colloquium: Modeling friction: From nanoscale to mesoscale, Rev. Mod. Phys. 85, 529-552 (2013).

[4] A.I. Vakis, V.A. Yastrebov, J. Scheibert, L. Nicola, D. Dini, C. Minfray, A. Almqvist, M. Paggi, S. Lee, G. Limbert, J. F. Molinari, G. Anciaux, R. Aghababaei, S. Echeverri Restrepo, A. Papangelo, A. Cammarata, P. Nicolini, C. Putignano, G. Carbone, S. Stupkiewicz, J. Lengiewicz, G. Costagliola, F. Bosia, R. Guarino, N.M. Pugno, M. H. Müser, M. Ciavarella, Modeling and simulation in tribology across scales: An overview, Tribology International 125, 169-199 (2018).

[5] M. Ternes, C. P. Lutz, C. F. Hirjibehedin, F. J. Giessibl, and A. J. Heinrich, The Force Needed to Move an Atom on a Surface, Science 319, 1066-1069 (2008).

[6] I. Szlufarska, M. Chandross, and R. W. Carpick, Recent advances in single-asperity nanotribology, J. Phys. D: Appl. Phys. 41 123001 (2008).

[7] Y. Mo, K. T. Turner, and I. Szlufarska, Friction laws at the nanoscale, Nature 457, 1116 (2009).
[8] S. Y. Krylov, and J. W. M. Frenken, The physics of atomic-scale friction: Basic considerations and open questions, Phys. Status Solidi **251**, 711-736 (2014).

[9] N. Manini, G. Mistura, G. Paolicelli, E. Tosatti, and A. Vanossi, Current trends in the physics of nanoscale friction, Advances in Physics: X **2**, 569-590 (2017).

[10] V. Bormuth, V. Varga, J. Howard, and E. Schäffer, Protein Fric- tion Limits Diffusive and Directed Movements of Kinesin Motors on Microtubules, Science **325**, 870 (2009).

[11] R. Sahli, G. Pallares, C. Ducottet, I. E. Ben Ali, S. Al Akhrass, M. Guibert, and J. Scheibert, Evolution of real contact area under shear and the value of static friction of soft materials, PNAS **115**, 471-476 (2018).

[12] P. Sens, Stick-slip model for actin-driven cell protrusions, cell polariza- tion, and crawling, PNAS **117**, 24670-24678 (2020).

[13] E. Liamas, S. D. Connell, S. N. Ramakrishna, and A. Sarkar, Probing the frictional properties of soft materials at the nanoscale, Nanoscale **12**, 2292 (2020).

[14] K. Holmberg, A. Matthews, and H. Ronkainen, Coatings tribology - contact mechanisms and surface design, Tribology International **31**, 107-120 (1998).

[15] V. L. Popov, Contact Mechanics and Friction: Physical Principles and Applications (Springer-Verlag, Berlin, 2010).

[16] K. Holmberg, A. Erdemir, Influence of tribology on global energy consumption, costs and emissions, Friction **5**, 263-284 (2017).

[17] M. Białas, J. Maciejewski, S. Kucharski, Friction coefficient of solid lubri- cating coating as a function of contact pressure: experimental results and microscale modeling, Continuum Mech. Thermodyn. **33**, 1733-1745 (2021).

[18] C. H. Scholz, Earthquakes and friction laws, Nature (London) **391**, 37 (1998).

[19] C. Marone, The effect of loading rate on static friction and the rate of fault healing during the earthquake cycle, Nature **391**, 69 (1998).

[20] E. G. Daub, J. M. Carlson, Friction, Fracture, and Earthquakes, Annual Review of Condensed Matter Physics **1**, 397-418 (2010).
[21] E. Fried, Energy release, friction, and supplemental relations at phase interfaces, Continuum Mech. Thermodyn. **7**, 111-121 (1995).

[22] D. Amitrano, J. R. Grasso, D. Hantz, From diffuse to localised damage through elastic interaction, Geophysical Research Letters **26**, 2109-2112 (1999).

[23] E. Gerde and M. Marder, Friction and fracture, Nature **413**, 285 (2001).

[24] O. Kresse, L. Truskinovsky, Lattice friction for crystalline defects: from dislocations to cracks, J. Mech. Phys. Sol. **52**, 2521-2543 (2004).

[25] F. Gimbert, D. Amitrano, and J. Weiss, Crossover from quasi-static to dense flow regime in compressed frictional granular media, EPL **104**, 46001 (2013).

[26] P. Biscari, M.F. Urbano, A. Zanzottera, G. Zanzotto, Intermittency in Crystal Plasticity Informed by Lattice Symmetry, Journal of Elasticity **123**, 85-96 (2016).

[27] K. Karimi, D. Amitrano, and J. Weiss, From plastic flow to brittle fracture: Role of microscopic friction in amorphous solids, Phys. Rev. E **100**, 012908 (2019).

[28] N. Gorbushin, G. Mishuris, and L. Truskinovsky, Frictionless Motion of Lattice Defects, Phys. Rev. Lett. **125**, 195502 (2020).

[29] D. Dowson, *History of Tribology* (Wiley, New York, 1998).

[30] T. W. J. de Geus, M. Popovic, W. Ji, A. Rosso, and M. Wyart, How collective asperity detachments nucleate slip at frictional interfaces, PNAS **116**, 23977-23983 (2019).

[31] T. D. B. Jacobs, A. Martini, Measuring and Understanding Contact Area at the Nanoscale: A Review, Applied Mechanics Reviews **69**, 060802 (2017).

[32] G. Binnig, C. F. Quate, and C. Gerber, Atomic Force Microscope, Phys. Rev. Lett. **56**, 930 (1986).

[33] J. N. Israelachvili, Adhesion forces between surfaces in liquids and condensable vapours, Surf. Sci. Rep. **14**, 109 (1992).

[34] Y. Dong, Q. Li, A. Martini, Molecular dynamics simulation of atomic friction: A review and guide, J. Vac. Sci. Technol. A **31**, 030801 (2013).

[35] F. Bonelli, N. Manini, E. Cadelano, and L.Colombo, Atomistic simulations of the sliding friction of graphene flakes, Eur. Phys. J. B **70**,
449–459 (2009).

[36] B. Q. Luan, S. Hyun, J. F. Molinari, N. Bernstein, and M. O. Robbins, Multiscale modeling of two-dimensional contacts, Phys. Rev. E **74**, 046710 (2006).

[37] M. Wolloch, G. Levita, P. Restuccia, and M. C. Righi, Interfacial Charge Density and Its Connection to Adhesion and Frictional Forces, Phys. Rev. Lett. **121**, 026804 (2018).

[38] L. Prandtl, Ein Gedankenmodell zur kinetischen Theorie der festen Körper, Z. Angew. Math. Mech. **8**, 85 (1928).

[39] G. A. Tomlinson, A molecular theory of friction, Philos. Mag. **7**, 905 (1929).

[40] V. L. Popov and J. A. T. Gray, Prandtl-Tomlinson model: History and applications in friction, plasticity, and nanotechnologies, Z. Angew. Math. Mech. **92**, 683-708 (2012).

[41] U. D. Schwarz and H. Hölscher, Exploring and Explaining Friction with the Prandtl-Tomlinson Model, ACS Nano **10**, 38-41 (2016).

[42] Y. Sang, M. Dubé, and M. Grant, Thermal Effects on Atomic Friction, Phys. Rev. Lett. **87**, 174301 (2001).

[43] E. Riedo, E. Gneco, R. Bennewitz, E. Meyer, and H. Brune, Interaction Potential and Hopping Dynamics Governing Sliding Friction, Phys. Rev. Lett. **91**, 084502 (2003).

[44] S. Yu. Krylov, K. B. Jinesh, H. Valk, M. Dienwiebel, and J. W. M. Frenken, Thermally induced suppression of friction at the atomic scale, Phys. Rev. E **71**, 065101(R) (2005).

[45] S. Yu Krylov, and J. W. M. Frenken, The crucial role of temperature in atomic scale friction, J. Phys.: Condens. Matter **20**, 354003 (2008).

[46] K. B. Jinesh, S. Yu. Krylov, H. Valk, M. Dienwiebel, and J. W. M. Frenken, Thermolubricity in atomic-scale friction, Phys. Rev. B **78**, 155440 (2008).

[47] L. Jansen, H. Hölscher, H. Fuchs, and A. Schirmeisen, Temperature Dependence of Atomic-Scale Stick-Slip Friction, Phys. Rev. Lett. **104**, 256101 (2010).

[48] D. Perez, Y. Dong, A. Martini, and A. F. Voter, Rate theory description of atomic stick-slip friction, Phys. Rev. B **81**, 245415 (2010).
[49] Martin H. Müser, Velocity dependence of kinetic friction in the Prandtl-Tomlinson model, Phys. Rev. B **84**, 125419 (2011).

[50] Y. Dong, A. Vadakkepatt, A. Martini, Analytical Models for Atomic Friction, Tribol. Lett. **44**, 367-386 (2011).

[51] P. C. Torche, T. Polcar, and O. Hovorka, Thermodynamic aspects of nanoscale friction, Phys. Rev. B **100**, 125431 (2019).

[52] A. Socoliuc, R. Bennewitz, E. Gnecco, and E. Meyer, Transition from Stick-Slip to Continuous Sliding in Atomic Friction: Entering a New Regime of Ultralow Friction, Phys. Rev. Lett. **92**, 134301 (2004).

[53] U. Dehlinger, Zur Theorie der Rekristallisation reiner Metalle, Ann. Phys. (Leipzig) **394**, 749 (1929).

[54] T. A. Kontorova and Ya. I. Frenkel, On the theory of the plastic deformation and twinning, Zh. Eksp. Teor. Fiz. **8**, 89-95 (1938).

[55] T. A. Kontorova and Ya. I. Frenkel, On the theory of plastic deformation and twinning. II, Zh. Eksp. Teor. Fiz. **8**, 1340-1348 (1938).

[56] T. A. Kontorova and Ya. I. Frenkel, On the theory of plastic deformation and twinning. III, Zh. Eksp. Teor. Fiz. **8**, 1349-1358 (1938).

[57] V. L. Pokrovskij and A. L. Talapov, Theory of incommensurate crystals (Harwood Academic Publishers, New York, 1984)

[58] O. M. Braun and Yu. S. Kivshar, Nonlinear dynamics of the Frenkel-Kontorova model, Physics Reports. **306**, 1-108 (1998).

[59] O. M. Braun and Yu. S. Kivshar, The Frenkel-Kontorova Model: Concepts, Methods, and Applications (Springer, Berlin, 2004).

[60] E. Bour, Théorie de la déformation des surfaces, Journal de l'École Impériale Polytechnique **19**, 1-48 (1862).

[61] R. E. Peierls, The size of a dislocation, Proc. Phys. Soc. **52**, 34 (1940).

[62] F. R. N. Nabarro, Dislocations in a simple cubic lattice, Proc. Phys. Soc. **59**, 256 (1947).

[63] J. P. Hirthand, and J. Lothe, Theory of Dislocations (Wiley, New York, 1982).

[64] S. Aubry, The twist map, the extended Frenkel-Kontorova model and the devil's staircase, Phys. D **7**, 240-258 (1983).
[65] S. Aubry, Devil's staircase and order without periodicity in classical condensed matter, J. Physique **44**, 147-162 (1983).

[66] M. Peyrard and S. Aubry, Critical behaviour at the transition by breaking of analyticity in the discrete Frenkel-Kontorova model, J. Phys. C: Solid State Phys. **16**, 1593-1608 (1983).

[67] O. Biham and D. Mukamel, Global universality in the Frenkel-Kontorova model, Phys. Rev. A **39**, 5326-5335 (1989).

[68] Y. Braiman, J. Baumgarten, J. Jortner, and J. Klafter, Symmetry-breaking transition in finite Frenkel-Kontorova chains, Phys. Rev. Lett. **65**, 2398 (1990).

[69] S. R. Sharma, B. Bergersen, and B. Joos, Aubry transition in a finite modulated chain, Phys. Rev. B **29**, 6335 (1984).

[70] T. Pruttivarasin, M. Ramm, I. Talukdar, A. Kreuter and H. Häffner, Trapped ions in optical lattices for probing oscillator chain models, New J. Phys. **13**, 075012 (2011).

[71] A. Benassi, A. Vanossi, and E. Tosatti, Nanofriction in cold ion traps, Nature Comm. **2**, 236 (2011).

[72] A. Bylinskii, D. Gangloff, I. Counts and V. Vuletić, Observation of Aubry-type transition in finite atom chains via friction, Nature Mat. **15** 717 (2016).

[73] J. Kiethe, R. Nigmatullin, D. Kalincev, T. Schmirander, and T.E. Mehlstäubler, Probing nanofriction and Aubry-type signatures in a finite self-organized system, Nature Comm. **8**, 15364 (2017).

[74] D. A. Gangloff, A. Bylinskii, and V. Vuletić, Kinks and nanofriction: Structural phases in few-atom chains, Phys. Rev. Research, **2**, 013380 (2020).

[75] T. Brazda, A. Silva, N. Manini, A. Vanossi, R. Guerra, E. Tosatti, and C. Bechinger, Experimental Observation of the Aubry Transition in Two-Dimensional Colloidal Monolayers, Phys. Rev. X **8**, 011050 (2018).

[76] R. Rosenberg, Why Is Ice Slippery?, Phys. Today **58**, 50-55 (2005).

[77] K. Mabuchi, K. Tanaka, D. Uchijima, R. Sakai, Frictional coefficient under banana skin, Tribol. Online **7**, 147-151 (2012).

[78] M. Z. Baykara, M. R. Vazirereshk, and A. Martini, Emerging superlubricity: A review of the state of the art and perspectives on future research, Appl. Phys. Rev. **5**, 041102 (2018).
[79] A. Kumar, Advancements in emerging superlubricity: A review of the atomistic models, simulation techniques and their applications to explore the state of ultra-low friction, Materials Today: Proceedings **42**, 884-892 (2021).

[80] M. Hirano and K. Shinjo, Atomistic locking and friction, Phys. Rev. B **41**, 11837 (1990).

[81] K. Shinjo and M. Hirano, Dynamics of friction: superlubric state, Surf. Sci. **283**, 473 (1993).

[82] M. H. Muser, Structural lubricity: Role of dimension and symmetry, Europhys. Lett. **66**, 97 (2004).

[83] M. Hirano, K. Shinjo, R. Kaneko, and Y. Murata, Observation of Super-lubricity by Scanning Tunneling Microscopy, Phys. Rev. Lett. **78**, 1448 (1997).

[84] M. Dienwiebel, G. S. Verhoeven, N. Pradeep, and J. W. M. Frenken, J. A. Heimberg, H. W. Zandbergen, Superlubricity of Graphite, Phys. Rev. Lett. **92**, 126101 (2004).

[85] E. Koren, E. Lörtscher, C. Rawlings, A. W. Knoll, U. Duerig, Adhesion and friction in mesoscopic graphite contacts, Science **348** 679 (2015).

[86] S. Kawai, A. Benassi, E. Gncco, H. Soede, R. Pawlak, X. Feng, K. Muellen, D. Passerone, C. A. Pignedoli, P. Ruffieux, R. Fasel, and E. Meyer, Sup erlubricity of graphene nanoribbons on gold surfaces, Science **351**, 957 (2016).

[87] P. Restuccia, M. Ferrario, and M. C. Righi, Monitoring water and oxygen splitting at graphene edges and folds: Insights into the lubricity of graphitic materials, Carbon **156**, 93-103 (2020).

[88] G. Losi, P. Restuccia, and M. C. Righi, Superlubricity in phosphorene identified by means of *ab initio* calculations, 2D Mater. **7**, 025033 (2020).

[89] X. Zhao, S. R. Phillpot, W. G. Sawyer, S. B. Sinnott, and S. S. Perry, Transition from Thermal to Athermal Friction under Cryogenic Conditions, Phys. Rev. Lett. **102**, 186102 (2009).

[90] I. Barel, M. Urbakh, L. Jansen and A. Schirmeisen, Multibond Dynamics of Nanoscale Friction: The Role of Temperature, Phys. Rev. Lett. **104**, 066104 (2010).

[91] S. Giordano, Spin variable approach for the statistical mechanics of folding and unfolding chains, Soft Matter **13**, 6877-6893 (2017).
44 Statistical mechanics of rate-independent stick-slip

[92] M. Caruel and L. Truskinovsky, Statistical mechanics of the Huxley-Simmons model, Phys. Rev. E **93**, 062407 (2016).

[93] M. Caruel, L. Truskinovsky, Physics of muscle contraction, Rep. Prog. Phys. **81**, 036602 (2018).

[94] M. Benedito, S. Giordano, Thermodynamics of small systems with conformational transitions: the case of two-state freely jointed chains with extensible units, J. Chem. Phys. **149**, 054901 (2018).

[95] M. Benedito, S. Giordano, Isotensional and isometric force-extension response of chains with bistable units and Ising interactions, Phys. Rev. E **98**, 052146 (2018).

[96] G. Florio, G. Puglisi, Unveiling the influence of device stiffness in single macromolecule unfolding, Sci. Rep. **9** 4997 (2019).

[97] L. Bellino, G. Florio, and G. Puglisi, The influence of device handles in single-molecule experiments, Soft Matter **15**, 8680-8690 (2019).

[98] G. Florio, G. Puglisi, and S. Giordano, Role of temperature in the decohesion of an elastic chain tethered to a substrate by onsite breakable links, Phys. Rev. Research **2**, 033227 (2020).

[99] A. Cannizzo, G. Florio, G. Puglisi, and S. Giordano, Temperature controlled decohesion regimes of an elastic chain adhering to a fixed substrate by softening and breakable bonds, J. Phys. A: Math. and Theor. **54**, 445001 (2021).

[100] L. Bellino, G. Florio, S. Giordano, and G. Puglisi, On the competition between interface energy and temperature in phase transition phenomena, Appl. Eng. Sci. **2**, 100009 (2020).

[101] A. Cannizzo, L. Bellino, G. Florio, G. Puglisi, S. Giordano, Thermal control of nucleation and propagation transition stresses in discrete lattices with non-local interactions and non-convex energy, Eur. Phys. J. Plus **137**, 569 (2022).

[102] A. Prados, A. Carpio, and L. L. Bonilla, Sawtooth patterns in force-extension curves of biomolecules: An equilibrium-statistical-mechanics theory, Phys. Rev. E **88**, 012704 (2013).

[103] L. L. Bonilla, A. Carpio, and A. Prados, Theory of force-extension curves for modular proteins and DNA hairpins, Phys. Rev. E **91**, 052712 (2015).

[104] D. De Tommasi, N. Millardi, G. Puglisi, and G. Saccomandi, An energetic model for macromolecules unfolding in stretching experiments, J.
R. Soc. Interface **10**, 20130651 (2013).

[105] I. Benichou and S. Givli, Structures undergoing discrete phase transformation, J. Mech. Phys. Sol. **61**, 94 (2013).

[106] H. A. Kramers, Brownian motion in a field of force and the diffusion model of chemical reactions, Physica (The Hague) **7**, 284 (1940).

[107] F. Manca, S. Giordano, P. L. Palla, R. Zucca, F. Cleri, and L. Colombo, Elasticity of flexible and semiflexible polymers with extensible bonds in the Gibbs and Helmholtz ensembles, J. Chem. Phys. **136**, 154906 (2012).

[108] J. H. Weiner, *Statistical Mechanics of Elasticity* (Dover Publications, New York, 1983).

[109] E. T. Whittaker, G. N. Watson, *A Course of Modern Analysis* (Cambridge University Press, Cambridge UK, 2021).

[110] I. S. Gradshteyn and I. M. Ryzhik, *Table of Integrals, Series and Products* (Academic Press, San Diego, 1965).

[111] M. Abramowitz, I.A. Stegun, *Handbook of Mathematical Functions* (Dover Publication, New York, 1970).

[112] F.W.J. Olver, D.W. Lozier, R.F. Boisvert, C.W. Clark, *NIST Handbook of Mathematical Functions* (National Institute of Standards and Technology and Cambridge University Press, New York, 2010).

[113] R. Bellman, *A brief introduction to theta functions* (New York, Holt, Rinehart and Co., 1961).

[114] R. Bellman and R. S. Lehman, The Reciprocity Formula for Multidimensional Theta Functions, Proceedings of the American Mathematical Society **12**, 954-961 (1961).

[115] B. Deconinck, M. Heil, A. Bobenko, M. Van Hoeij, and M. Schmies, Computing Riemann Theta Functions, Mathematics of Computation **73**, 1417-1442 (2003).

[116] B. V. Petukhov, M. Bartsch, and U. Messerschmidt, Temperature dependence of the flow stress and the strain rate sensitivity at the transition from the Peierls mechanism to pinning by localized obstacles, Eur. Phys. J. Appl. Phys. **9**, 89-95 (2000).

[117] P. Cordier, J. Amodeo, and P. Carrez, Modelling the rheology of MgO under Earth's mantle pressure, temperature and strain rates, Nature **481**, 177 (2012).
[118] Y. Kamimura, K. Edagawa, S. Takeuchi, Experimental evaluation of the Peierls stresses in a variety of crystals and their relation to the crystal structure, Acta Materialia **61**, 294-309 (2013).

[119] J. Amodeo, S. Merkel, C. Tromas, P. Carrez, S. Korte-Kerzel, P. Cordier, and J. Chevalier, Dislocations and Plastic Deformation in MgO Crystals: A Review, Crystals **8**, 240 (2018).

[120] J. F. Curry, A. R. Hinkle, T. F. Babuska, M. A. Wilson, M. T. Dugger, B. A. Krick, N. Argibay, and M. Chandross, Atomistic Origins of Temperature-Dependent Shear Strength in 2D Materials, ACS Appl. Nano Mater. **1**, 5401-5407 (2018).

[121] O. Y. Fajardo and J. J. Mazo, Effects of surface disorder and temperature on atomic friction, Phys. Rev. B **82**, 035435 (2010).

[122] Y. Peng, C. M. Serfass, A. Kawazoe, Y. Shao, K. Gutierrez, C. N. Hill, V. J. Santos, Y. Visell, and L. C. Hsiao, Elastohydrodynamic friction of robotic and human fingers on soft micropatterned substrates, Nat. Mater. **20**, 1707-1711 (2021).

[123] J. C. Mergel, J. Scheibert, R. A. Sauer, Contact with coupled adhesion and friction: Computational framework, applications, and new insights, J. Mech. Phys. Sol. **146**, 104194 (2021).

[124] H. Risken, *The Fokker-Planck equation* (Springer Verlag, Berlin, 1989).

[125] W. T. Coffey, Yu. P. Kalmykov, and J. P. Waldron, *The Langevin equation* (World Scientific, Singapore, 2004).

[126] I. Benichou, S. Givli, Rate Dependent Response of Nanoscale Structures Having a Multiwell Energy Landscape, Phys. Rev. Lett. **114**, 095504 (2015).

[127] I. Benichou, Y. Zhang, O. K. Dudko, S. Givli, The rate dependent response of a bistable chain at finite temperature, J. Mech. Phys. Sol. **95**, 44 (2016).