On the connection between sliding friction and phonon lifetimes: Thermostat induced thermolubricity effects in molecular dynamics simulations

TL;DR

This study investigates how thermostat parameters influence phonon behavior and friction in molecular dynamics simulations of nanotribology, revealing the roles of phonon lifetime, substrate temperature, and sliding speed in thermolubricity effects.

Contribution

The paper clarifies the relationship between thermostat tuning, phonon lifetimes, and friction, providing insights into thermolubricity in molecular dynamics simulations with harmonic and anharmonic interactions.

Findings

Friction maximum occurs at optimal damping parameter GAMMA.

Substrate temperature TSUB shows a pronounced minimum at high sliding speeds.

Thermolubricity effects differ between harmonic and anharmonic substrate interactions.

Abstract

A typical nanotribology simulation setup is the semi-infinite substrate, featuring a sliding bead on top, and with the lower substrate layers thermostatted to control temperature. A challenge is dealing with phonons that backreflect from the substrate lower boundary, as these will artificially reduce the friction FR acting on the sliding bead. One proposed solution is to use a Langevin thermostat, operating at temperature TLAN, and with the corresponding damping parameter, GAMMA, optimally tuned such that FR is maximized [Benassi etal, PHYSREVB, V82, P081401, 2010]. In this paper, the method is revisited, and related to the substrate phonon lifetime, the substrate temperature TSUB, and the sliding speed. At low sliding speed, where the time between stick-slip events is large compared to the phonon lifetime, we do not observe much dependence of FR on GAMMA, and here thermostat tuning is not required. At high sliding speed, upon varying GAMMA, we confirm the aforementioned friction maximum, but also observe a pronounced minimum in TSUB, which here deviates from TLAN. For substrate particle interactions that are strongly anharmonic, the variation of FR with GAMMA can be understood as a manifestation of thermolubricity, backreflections being essentially unimportant. In contrast, for harmonic interactions, where phonon lifetimes become very long, FR is strongly affected by backreflecting phonons, though not enough to overturn thermolubricity.

Equations

$${\vec{F}}_{\mathrm{lan}}=-\gamma m\vec{v}+\sqrt{(24k_B{T}_{\mathrm{lan}}\gamma m/\delta t)}\vec{h},$$

$$u_{\mathrm{afm}}(r)=c_1ϵ\left[{\left(\frac{c_2}{r}\right)}^{12}-{\left(\frac{c_2}{r}\right)}^6\right]+c_{3}r+c_4,$$

$$\tilde{D}(\vec{k})=\frac{4{\alpha }_2}{m}\sum _{\nu =0,1,2}[1-\cos (a\vec{k}\cdot {\hat{n}}_{\nu })]{\hat{n}}_{\nu }\otimes {\hat{n}}_{\nu },$$

$$\begin{array}{rl}{A}_p(\vec{k},t)& =\\ e^{ı{\omega }_p(\vec{k})t}& \sqrt{\frac{N}{2{\omega }_p(\vec{k})}}\left({\omega }_p(\vec{k})Q_p(\vec{k},t)+ı{\dot{Q}}_p(\vec{k},t)\right).\end{array}$$

Full text

I Introduction

Molecular Dynamics (MD) simulations are an established tool in the field of (nano)tribology. Relatively straightforward to setup, these simulations can be used to study the dynamics of a contact of some kind (e.g. the tip of an atomic force microscope) sliding across a surface. A practical problem facing all these simulations is how to remove the excess energy pumped into the system at the sliding contact, for, if this energy is not removed, the system temperature would steadily rise. The problem has received considerable attention Benassi et al. (2010); Toton et al. (2010); Benassi et al. (2012); Kantorovich (2008); Li and E (2007); Meyer and Reuter (2014) and one practical solution that has emerged is to remove the excess energy via a thermostat acting on the lower part of the substrate, as depicted schematically in Fig. 1(a). The thermostat will, of course, alter the true dynamics of the system in the lower region, but if the simulation cell is large enough, one might hope that the dynamics in the upper region remains unaffected, facilitating meaningful friction measurements.

Alas, simulations using the geometry of Fig. 1(a) unambiguously show that the measured friction force $F_{\rm fr}$ depends quite sensitively on thermostat details Benassi et al. (2010, 2012); Apostoli et al. (2017). In these works, this dependence is mainly attributed to acoustic phonons that get “backreflected” from the lower simulation box boundary. In very simple terms, as the bead is dragged across the surface, acoustic phonons are generated in the underlying substrate. The fixed atoms at the bottom of the simulation cell act as a “mirror” reflecting these phonons back toward the bead, potentially allowing the bead to regain some of its energy, implying that the measured friction force will be “too small”, compared to what it would be in the infinite system without backreflections. The key point to note next is that, by tuning the thermostat, which acts on the particles directly above the frozen bottom layer, the degree of phonon backreflection can be regulated to some extent Toton et al. (2010). Hence, it seems logical to tune the thermostat such that the measured friction force is maximized, since, under this condition, the artificial reduction of friction due to phonon backreflections must necessarily be minimized. Indeed, taking for concreteness a Langevin thermostat, where the adjustable parameter is the damping coefficient $\gamma$, the existence of an “optimal” damping coefficient, where the friction force reaches a maximum, was strikingly confirmed Benassi et al. (2010, 2012); Apostoli et al. (2017).

The purpose of this paper is to relate these findings to the phonon properties of the underlying substrate, specializing to the regime where the motion of the bead is of the stick-slip type. Quite surprisingly, we confirm the previous findings Benassi et al. (2010, 2012); Apostoli et al. (2017) in certain dynamic regimes, but not in others. To be precise: the friction maximum upon variation of the thermostat damping parameter $\gamma$ is only observed at high sliding speeds; at low sliding speeds, the friction force does not reveal any systematic $\gamma$ dependence. The defining criterion turns out to be the typical phonon lifetime in the substrate: When the latter is small compared to the time between single stick-slip events, the influence of the thermostat on friction vanishes. In the opposite limit, friction depends quite sensitively on thermostat details, but the extent to which this impedes meaningful measurements depends sensitively on the degree of anharmonicity in the substrate particle interactions. For anharmonic substrates, the dependence of friction on thermostat details is, for the most part, explained by thermolubricity Riedo et al. (2003); Jansen et al. (2010).

We will, in what follows, present the results of MD simulations leading to these conclusions. Our simulations are based on a two-dimensional (2D) model system, described in Section II, which is similar in spirit to that of Ref. Benassi et al. (2010). In Section III, we present the corresponding friction measurements, as well as the substrate phonon properties required to explain these measurements. The relation between friction and phonon properties, in particular the phonon lifetime, is discussed in Section IV. We end with a summary and some recommendations in Section V.

II Model, Methods, and Definitions

We simulate the 2D setup of Fig. 1(a). The substrate is a hexagonal lattice, whose primitive cell is spanned by the vectors $\vec{a}_{1}=a\hat{x}$ and $\vec{a}_{2}=\frac{a}{2}(\hat{x}+\sqrt{3}\hat{z})$, each cell containing exactly one particle of mass $m$, with $a$ the lattice constant. The cell is replicated $\sqrt{N}$ times in both directions $\vec{a}_{1,2}$, $N=1600$ being the particle number (the simulation box is thus triclinic). We apply periodic boundary conditions in the (horizontal) $x$-direction, but not in the (vertical) $z$-direction. Between nearest neighboring particles, springs are attached, the energy of a single spring being $u_{\rm spr}(r)=\sum_{n=2}^{4}\alpha_{n}(r-a)^{n}$, where $r$ is the spring length. To anchor the substrate, the lowest row of particles have their positions fixed.

A Langevin thermostat Schneider and Stoll (1978); DÜNWEG and PAUL (1991) is applied to the row of particles directly above the fixed layer. To the forces acting on these particles, damping and stochastic terms are added:

[TABLE]

in addition to the forces arising from the springs. Here, $\gamma$ is the damping parameter, $\delta t$ the MD integration time step, $k_{B}$ the Boltzmann constant, $\vec{v}$ the particle velocity, and $\vec{h}$ a two-dimensional vector with components drawn uniformly from $[-0.5:0.5]$. For each thermostatted particle, at each MD step, a new vector $\vec{h}$ is to be generated! Note that the Langevin limits $\gamma\to 0$ and $\gamma\to\infty$ are equivalent, both corresponding to a system without thermostat. For $\gamma=0$ this is intuitively clear, since here $\vec{F}_{\rm lan}=0$. For $\gamma\to\infty$, the stochastic term, $\propto\sqrt{\gamma}$, becomes negligible compared to the viscous term, $\propto\gamma$, leading to an infinite viscous force, which will impede any motion of the thermostat particles. Hence, $\gamma\to\infty$ is similar to $\gamma=0$, but with an extra layer of frozen particles (direct simulations in the regime $\gamma\to\infty$ are, however, numerically challenging, since a progressively smaller MD timestep $\delta t$ is then required).

The temperature $T_{\rm lan}$ appearing in Eq. (1) is the thermostat temperature. In thermal equilibrium, this is also the temperature the substrate will adopt, $T_{\rm sub}=T_{\rm lan}$, with $T_{\rm sub}$ computed directly from the particle kinetic energy, $k_{B}T_{\rm sub}=(2/d)\langle K\rangle$, with $d$ the spatial dimension, and where $\langle K\rangle$ is the average kinetic energy of the substrate particles (when computing $\langle K\rangle$, we exclude the bead particle, as well as the frozen and thermostatted substrate layers). We already announce here that, in non-equilibrium situations, $T_{\rm sub}$ and $T_{\rm lan}$ can be very different!

On top of the substrate, a bead particle is placed, having the same mass $m$ as the substrate particles. The bead is subjected to a vertical load, $\vec{F}_{L}=-L\hat{z}$, pressing it down onto the underlying substrate. In addition, the bead is attached to one end of a harmonic spring (spring constant $k$, zero rest length), while the other end of the spring is dragged to the right with constant velocity $v$. The spring is assumed to act only in the lateral $x$-direction, that is, it induces a force $\vec{F}_{\rm spr}=k[X(t)-vt]\hat{x}$ onto the bead, where $t$ is the time, and $X(t)$ the $x$-coordinate of the AFM bead at time $t$. The average friction force magnitude is then obtained by time averaging the lateral spring force: $F_{\rm fr}=k\langle X(t)-vt\rangle_{t}$.

The interaction energy between the AFM bead and the underlying substrate particles is taken to be a sum of (short-ranged) pair potentials, $U_{\rm afm}=\sum_{s}^{\prime}u_{\rm afm}(R_{s})$, with the sum over all substrate particles $s$ whose distance $R_{s}$ from the AFM bead is smaller than a specified cutoff distance: $R_{s}<R_{c}=2.5a$. The pair potential is of the $(12,6)$ Lennard-Jones (LJ) form, with added linear term:

[TABLE]

where $\epsilon$ sets the energy scale, and with the constants $c_{i}$ chosen such that the minimum of the pair potential is located at $r=a$, with corresponding value $u_{\rm afm}(a)=-0.6\epsilon$, while at the cutoff $u_{\rm afm}(R_{c})=u_{\rm afm}^{\prime}(R_{c})=0$, i.e. there is no force discontinuity.

We adopt LJ units throughout: $a=\epsilon=m=k_{B}\equiv 1$. The vertical load is set to $L=10$, the AFM spring constant to $k=5$. The time evolution of the system is obtained via standard molecular dynamics, using the velocity-verlet algorithm as implemented in LAMMPS Plimpton (1995), with integration time step $\delta t=0.001$. The sliding velocity $v$ will be varied, but we take care to remain in the “stick-slip” regime, the AFM bead thus spending most of its time in low-energy positions on the substrate, the transitions between such positions being rapid. Since our primary interest is friction on solid supports, we consider low temperatures only, typically choosing $T_{\rm lan}=0.035$, which is well below melting Wierschem and Manousakis (2011).

In what follows, we shall refer to harmonic and anharmonic substrates. For the harmonic substrate, the spring energy parameters $\alpha_{2}=36,\alpha_{3}=\alpha_{4}=0$; for the anharmonic substrate $\alpha_{2}=36,\alpha_{3}=-252,\alpha_{4}=1113$. These parameters correspond to a Taylor expansion of a $(12,6)$ LJ potential around its minimum, with well-depth $\epsilon=1$, and minimum located at $r=1$. We emphasize that both substrates are strictly speaking anharmonic, since, even for the harmonic version, the potential energy is not quadratic in the particle displacements. However, owing to the explicit absence of 3rd and 4th order terms, anharmonicity effects should be much weaker, which suffices for our purposes.

II.1 Phonon properties

In our analysis, we will relate our friction measurements to the vibrational properties of the substrate using the language of phonons. Each phonon is characterized by a wavevector $\vec{k}=(k_{x},k_{z})$ and polarization $p$. In a finite system, the number of wavevectors inside the first Brillouin zone (FBZ) equals the number of primitive unit cells, which here equals the number of particles $N$, since the hexagonal lattice has a 1-atom basis. In Fig. 1(b), we show the FBZ of the hexagonal lattice, properly oriented, i.e. corresponding to the primitive vectors $\vec{a}_{1,2}$ of the lattice in real space. The center of the FBZ marks the $\Gamma$ point: $\vec{k}=(0,0)$. Each wavevector yields two polarizations, longitudinal-acoustic ($p=\rm LA$) and transverse-acoustic ($p=\rm TA$), the respective speeds of sound being $c_{\rm TA}=(\sqrt{3}/2)a\sqrt{\alpha_{2}/m}$ and $c_{\rm LA}=(3/2)a\sqrt{\alpha_{2}/m}$.

For each phonon propagating in the direction $\vec{k}$, we introduce its lifetime $\tau(\vec{k},p)$, and energy $E(\vec{k},p)$, with $p\in\rm TA,LA$. These quantities can be obtained from the substrate particle positions and velocities in the MD trajectory (details in Appendix). Phonons with short lifetimes $\tau(\vec{k},p)$ are dissipating, since these phonons quickly distribute their energy over other phonon modes, leading to rapid thermalization of the entire phonon population. Analogously, phonons with long lifetimes $\tau(\vec{k},p)$ are non-dissipating, since these can store their energy over longer times, delaying thermalization. In applications, thermalization of the substrate, i.e. the conversion of the bead’s kinetic energy into heat, is typically undesirable, and so to understand how the sliding motion of the bead couples to the phonons of the underlying substrate, is very important.

III Results

III.1 Friction measurements

We first measure how the friction force, $F_{\rm fr}$, depends on the Langevin thermostat parameter, $\gamma$, as well as on the sliding speed, $v$. The temperature of the Langevin thermostat $T_{\rm lan}=0.035$. However, this being a non-equilibrium situation, the actual temperature of the substrate, $T_{\rm sub}$, may well deviate from $T_{\rm lan}$. For each friction measurement, $\sim 20\cdot 10^{6}$ MD steps were applied to bring the system into a steady state, followed by production runs of at least $500\cdot 10^{6}$ MD steps, during which the average friction force was measured. In Fig. 2(a), we show, for the harmonic substrate, the dependence of $F_{\rm fr}$ on $\gamma$ for three values of the sliding speed, $v$, as indicated. The lower panel, Fig. 2(b), shows the variation of $T_{\rm sub}$ with $\gamma$. For the anharmonic substrate, the analogous analysis is presented in Fig. 3.

At first sight, the data of Fig. 2 and Fig. 3 look very similar. For both the harmonic and anharmonic substrate, at the highest sliding speed $v=0.01$ considered (which remains well below the speed of sound), the curve of $F_{\rm fr}$ versus $\gamma$ reveals a maximum. The existence of the friction maximum agrees with previous studies Benassi et al. (2010, 2012); Apostoli et al. (2017), and is attributed to phonon backreflections, which are maximally suppressed at the friction maximum. However, still considering the highest sliding speed, the data also reveal that the substrate temperature, $T_{\rm sub}$, varies with $\gamma$, reaching, in fact, a minimum at the friction maximum (with the value of $T_{\rm sub}$ then being close to $T_{\rm lan}$ of the thermostat). Since friction is known to depend on temperature as well, an effect called thermolubricity Riedo et al. (2003); Jansen et al. (2010), our data indicate that the observed friction maximum could be the manifestation of two effects, namely, suppression of phonon backreflections and reducing the substrate temperature, both of which would increase friction.

Interestingly, lowering the sliding speed $v$, these effects gradually vanish, $F_{\rm fr}$ and $T_{\rm sub}$ then being essentially independent of $\gamma$, with $T_{\rm sub}$ close to $T_{\rm lan}$ of the thermostat. For both the harmonic and anharmonic substrate, reducing the sliding speed leads to lower friction. Note also that, for the anharmonic substrate, friction significantly exceeds that of the harmonic substrate.

III.2 Effect of sliding on substrate phonons

We now identify the phonon modes most excited by the sliding bead, at the highest considered sliding speed $v=0.01$, thermostat temperature $T_{\rm lan}=0.035$, for various values of $\gamma$, and for both substrate types (harmonic, anharmonic). To this end, we measure the time-averaged energy $E(\vec{k},p)$ of the phonon characterized by the wavevector $\vec{k}$ and polarization $p\in\rm TA,LA$. As measure for the degree of excitation, we use the quantity $\kappa(\vec{k},p)\equiv E(\vec{k},p)/{\rm Mdn}[E(\vec{k},p)]$, where ${\rm Mdn}[E(\vec{k},p)]$ is the median (not mean) of the obtained $E(\vec{k},p)$ values. In thermal equilibrium, at low temperature, one may assume equipartition approximately holds, in which case all phonon modes should have the same energy: $E(\vec{k},p)={\rm Mdn}[E(\vec{k},p)]=k_{B}T_{\rm sub}$, implying $\kappa(\vec{k},p)=1$. Under driving, we still expect a substantial fraction of the phonon population to be thermalized, with the exception of modes that couple strongly to the sliding motion of the bead, whose energy should then exceed the thermal value, implying $\kappa(\vec{k},p)>1$. To confirm these ideas, we identify, for each polarization separately, the set $S_{n}$ of $n=10$ phonons with the largest value of $\kappa(\vec{k},p)$, and mark the corresponding wavevectors in the FBZ. For these phonons, we additionally compute the average lifetime, $\tau_{n}=(1/n)\sum^{\prime}\tau(\vec{k},p)$, and the average excitation, $\kappa_{n}=(1/n)\sum^{\prime}\kappa(\vec{k},p)$, where the sum is over all phonons in the set $S_{n}$. For the harmonic (anharmonic) substrate, the result of this procedure is shown in Fig. 4 [Fig. 5].

For the harmonic substrate, $\tau_{n}$ varies strongly with $\gamma$ [Fig. 4]. For both polarizations, by increasing $\gamma$, $\tau_{n}$ first reaches a maximum, then decreases to a minimum at $\gamma\sim 10-100$, which roughly corresponds to the friction maximum [Fig. 2(a)]. In addition, for small $\gamma\leq 0.1$, we observe a vertical band around the $\Gamma$ point of LA phonons with very long lifetimes; by increasing $\gamma$, this band disappears. Apparently, the effect of increasing $\gamma$ is to drastically reduce the lifetime of long wavelength, vertically propagating, LA phonons. Also of interest is the variation of ${\rm Mdn}[E(\vec{k},p)]$ with $\gamma$. As can be seen by comparing to Fig. 2(b), the median rather closely follows the substrate temperature, $k_{B}T_{\rm sub}$, confirming that many phonons are still thermalized. However, phonons inside the sets $S_{n}$, i.e. those which couple most strongly to the bead, have energies that exceed the thermal value by factors typically $\kappa_{n}\sim 3$, the exception being LA phonons at small $\gamma$, for which the enhancement is significantly larger, $\kappa_{n}>6$.

For the anharmonic substrate, we observe a much weaker dependence of $\tau_{n}$ on $\gamma$, whose value, in comparison to the harmonic substrate, is now much smaller [Fig. 5]. In contrast to the harmonic substrate, vertical bands in the FBZ persist for all values of $\gamma$, i.e. a preferred suppression of the lifetime of vertically propagating LA phonons upon increasing $\gamma$ does not take place. The typical phonon enhancement $\kappa_{n}\sim 3$ for all values of $\gamma$. As before, the median closely follows the substrate temperature, $k_{B}T_{\rm sub}$, see Fig. 3(b). There is, however, one subtle difference: For the harmonic substrate, ${\rm Mdn}[E(\vec{k},p)]$ for TA and LA modes differs by about $\pm 0.001$, while for the anharmonic substrate, the difference is only $\pm 0.0001$, i.e. ten times smaller. This indicates that the deviation from thermal equilibrium is largest for the harmonic substrate.

III.3 Typical phonon lifetimes

Finally, we still measure the typical substrate phonon lifetime, $\langle\tau_{p}\rangle$, where typical means an average over all wavevectors (details in Appendix). When computing $\langle\tau_{p}\rangle$, we disregard whether the phonon modes being averaged over are actually excited by the bead, which is the main difference from $\tau_{n}$ defined previously. In Fig. 6, we show $\langle\tau_{p}\rangle$, for both substrate types, and for both polarizations $p$. Results are shown for the equilibrium case ($v=0$: circles), and under sliding ($v=0.01$: squares) using again $T_{\rm lan}=0.035$. For the equilibrium case, after initial equilibration, active thermostatting is not required to maintain the temperature. In this case, we can “turn off” the Langevin thermostat, yielding the star symbols in Fig. 6.

For the harmonic substrate, Fig. 6 reveals a strong dependence of $\langle\tau_{p}\rangle$ on $\gamma$, reaching a minimum at $\gamma\sim 10$, which again is close to the friction maximum [cf. Fig. 2(a)]. In addition, upon sliding, $\langle\tau_{p}\rangle$ is significantly reduced from its equilibrium value. Apparently, the sliding bead introduces extra “noise” into the substrate, in addition to that of the Langevin stochastic term, promoting phonon mixing. We also observe that, in the limits $\gamma\to 0$ and $\gamma\to\infty$, the equilibrium values (circles) approach the star symbols obtained without Langevin thermostat.

For the anharmonic substrate, compared to the harmonic one, Fig. 6 reveals an overall much smaller phonon lifetime: Only for $\gamma\sim 10$ and under sliding, do the values become somewhat close. Furthermore, there are pronounced qualitative differences. For LA modes, neither $\gamma$ nor sliding systematically affect $\langle\tau_{p}\rangle$, the corresponding values remaining close (in absolute terms) to those obtained without Langevin thermostat (stars). For TA modes, a systematic dependence on $\gamma$ remains absent, but sliding does appear to slightly reduce $\langle\tau_{p}\rangle$. Still, in absolute terms, also for TA modes, the observed lifetimes are all rather similar, remaining close to those obtained without Langevin thermostat.

IV Discussion

IV.1 Friction at low sliding speeds

At the lowest considered sliding speed, $v=10^{-4}$, for both substrate types, neither $F_{\rm fr}$ nor $T_{\rm sub}$ show any appreciable dependence on the Langevin damping parameter $\gamma$ [Fig. 2 and Fig. 3]. This can be understood from the typical phonon lifetime $\langle\tau_{p}\rangle$. At $v=10^{-4}$, the time between stick-slip events $\tau_{\rm slip}=a/v=10000$ LJ time units, which far exceeds $\langle\tau_{p}\rangle$, for all values of $\gamma$, and for both substrate types [Fig. 6]. The reduction of friction due to phonon backreflections thus cannot occur since any coherence between the bead and the substrate phonons generated during slip, will long have decayed by the time of the next slip event. Consecutive slip events are thus uncorrelated, and, at the start of each such event, the underlying substrate in a state of thermal equilibrium. Friction simulations in this regime are thus relatively straightforward, as essentially any value of $\gamma$ suffices, always yielding a friction value corresponding to the temperature of the thermostat, which then equals that of the substrate: $T_{\rm lan}=T_{\rm sub}$.

At the intermediate sliding speed, $v=10^{-3}$, $\tau_{\rm slip}=1000$ LJ time units. For the anharmonic substrate, this still far exceeds $\langle\tau_{p}\rangle$. In line with our argumentation, the corresponding friction data do not show any appreciable dependence on $\gamma$ [Fig. 3(a)]. For the harmonic substrate, at small $\gamma$, $\langle\tau_{p}\rangle$ is still below $\tau_{\rm slip}$, but not that much lower, which means there could be some influence of backreflections. Consistent with this interpretation, a slight decrease of $F_{\rm fr}$ at small $\gamma$ is visible then [Fig. 2(a)].

IV.2 Friction at high sliding speed

At high sliding speed, $v=0.01$, the behavior is far more subtle, since now $F_{\rm fr}$ and $T_{\rm sub}$ both depend on $\gamma$. In addition, the time between stick-slip events is now much shorter: $\tau_{\rm slip}=100$ LJ time units. For the anharmonic substrate, $\tau_{\rm slip}$ still amply exceeds the lifetime $\tau_{n}$ of even the longest living phonon modes [Fig. 5]. We thus do not expect phonon backreflections to be important here. Instead, for the anharmonic substrate, the variation of friction with $\gamma$ should largely be due to the changing substrate temperature, $T_{\rm sub}$, via thermolubricity. In contrast, for the harmonic substrate, $\tau_{\rm slip}$ is comparable to $\tau_{n}$, or even exceeds it [Fig. 4]. In this case, friction may well be affected by phonon backreflections, in addition to thermolubricity effects. To verify, we choose a fixed value of $\gamma$, then vary the thermostat temperature $T_{\rm lan}$, and monitor how $F_{\rm fr}$ changes with $T_{\rm sub}$ (this thus requires additional simulations to be performed, extending the range of $T_{\rm lan}$). For the harmonic substrate, for each value of $\gamma$, a different curve “$F_{\rm fr}$ versus $T_{\rm sub}$” is obtained [Fig. 7(a)]. In contrast, for the anharmonic substrate, the data for different $\gamma$ collapse onto a single curve [Fig. 7(b)].

Fig. 7 strikingly confirms the expectations. For the anharmonic substrate, friction is set by the substrate temperature, $T_{\rm sub}$, the value of $\gamma$ being irrelevant here [Fig. 7(b)]. In addition, the decrease of $F_{\rm fr}$ with $T_{\rm sub}$ agrees with thermolubricity Riedo et al. (2003); Jansen et al. (2010), the data being well described by $F_{c}-F_{\rm fr}\propto T^{2/3}$ Dong et al. (2012); Sang et al. (2001), indicated by the dashed curve in Fig. 7(b). In contrast, for the harmonic substrate, Fig. 7(a) shows that $F_{\rm fr}$ depends on both $T_{\rm sub}$ and $\gamma$. If one increases $T_{\rm sub}$ keeping $\gamma$ fixed, $F_{\rm fr}$ decreases, similar to thermolubricity. Changing $\gamma$ affects both the strength of phonon backreflections, and the substrate temperature $T_{\rm sub}$, leading to a shift of the “$F_{\rm fr}$ versus $T_{\rm sub}$” curves.

Even though, for the harmonic substrate, backreflections are thus important, they are not strong enough to overturn thermolubricity. As Fig. 2 shows, a strong anti-correlation between $F_{\rm fr}$ and $T_{\rm sub}$ remains clearly visible. In contrast, following Fig. 4, the correlation between $F_{\rm fr}$ and $\tau_{n}$ is much weaker (if backreflections were to dominate, $F_{\rm fr}$ should reach a minimum at $\gamma\sim 0.1-1$, since here $\tau_{n}$ is maximal). Hence, at least for the present data, the occurrence of the friction maximum with $\gamma$ appears to be mainly due to the reduction of $T_{\rm sub}$, rather than to the elimination of phonon backreflections.

IV.3 Role of thermostat

We still mention one subtle point concerning the use of a (Langevin) thermostat. For the harmonic substrate, phonon lifetimes are extremely sensitive to $\gamma$ [cf. Fig. 4 and Fig. 6(a,c)]. These lifetimes thus mainly reflect the influence of the thermostat, i.e. they are “biased”. An “unbiased” estimate could be obtained by removing the thermostat, which obviously can only be done in equilibrium, in which case the lifetime becomes maximally large [Fig. 6(a,c): star symbols]. Alas, under driving, removing the thermostat is not possible, and here one unavoidably simulates in a regime where the phonon dynamics of the entire substrate is largely dictated by the thermostat, which is against the notion of the thermostat being just a local perturbation, affecting only the lowest few substrate layers.

In contrast, for the anharmonic substrate, the phonon lifetime does not strongly depend on $\gamma$ [cf. Fig. 5 and Fig. 6(b,d)]. In this case, the lifetime is determined by the substrate particle interactions, thus reflecting a “true” (unbiased) material property. One could say that the anharmonic substrate is “self-thermalizing”: Due to the 3rd and 4th order terms in the particle interaction, phonon scattering is strong, and so the energy injected by the sliding bead at the top of the substrate, is quickly converted into heat. In order to control the temperature, all the thermostat has to do is remove this excess heat, but there is no need for the thermostat to actively generate heat, i.e. a stochastic term is not needed. In fact, for the anharmonic substrate, even a simple velocity-rescaling thermostat Bussi et al. (2007) suffices, which simply rescales the velocities of the particles being thermostatted every so often to match some chosen thermostat temperature $T_{\rm lan}$. In this case, varying now the frequency at which rescaling is performed, one still finds that $T_{\rm sub}$ deviates from $T_{\rm lan}$, but the curve of $F_{\rm fr}$ versus $T_{\rm sub}$ coincides with the one obtained using the Langevin thermostat [Fig. 7(b): hexagon symbols].

V Conclusions

In this work, we have studied the problem of how to thermostat a nanotribology MD simulation, following-up on previous works Benassi et al. (2010, 2012); Apostoli et al. (2017). The main insight of this work has been to relate thermostat effects to the lifetime of phonons in the underlying substrate. One finding is that, provided the sliding speed is sufficiently low, thermostat effects become negligible, the defining criterion being the typical phonon lifetime, which must be small compared to the time between stick-slip events. This result is convenient for modeling (clean) single-asperity AFM contacts, since here the use of a lower sliding speed would also better resemble experimental conditions Li et al. (2011).

The opposite regime, where the phonon lifetime is comparable to or larger than the time between stick-slip events, could occur with single molecule adsorbates, due to the now much larger possible sliding velocity Buldum et al. (1999). In this regime, we confirm previous findings Benassi et al. (2010, 2012); Apostoli et al. (2017), namely, that the variation of friction with the Langevin damping parameter, $\gamma$, reveals a maximum. In addition, we also find that the substrate temperature, $T_{\rm sub}$, depends strongly on $\gamma$ then, reaching a minimum at the friction maximum. In this situation, we recommend plotting $F_{\rm fr}$ versus $T_{\rm sub}$. In case the latter yields a single curve, the $\gamma$ dependence could be a manifestation of thermolubricity Riedo et al. (2003); Jansen et al. (2010), which one can verify by comparing to theoretical predictions Dong et al. (2012); Sang et al. (2001). The data may then be deemed to reliably describe friction, but at the substrate temperature, $T_{\rm sub}$, which here deviates from $T_{\rm lan}$ of the thermostat, typically exceeding it.

In contrast, when measurements of $F_{\rm fr}$ versus $T_{\rm sub}$ do not yield a single curve, one likely has long-lived backreflecting phonons in the system, which will strongly reduce friction. In this case, one could try to minimize them, by tuning $\gamma$ to the friction maximum, which follows the recommendation of earlier works Benassi et al. (2010, 2012); Apostoli et al. (2017). The tuning of $\gamma$ has two effects, namely, (1) a reduction of the lifetime of backreflecting phonons, as shown by the vanishing of the vertical LA bands in Fig. 4, making these phonons more dissipating, and (2) a reduction of the substrate temperature, $T_{\rm sub}$, which increases friction due to thermolubricity. It is the combination of both effects that gives rise to the friction maximum.

Our phonon analysis also yields interesting insights into the underlying frictional mechanisms for a single sliding contact, at least for the present 2D model. Upon sliding, the energy of the bead appears to be mainly transferred into vertically propagating phonon modes (the “patterns” in Fig. 4 and Fig. 5 have a clear tendency to align in the $z$-direction). Depending now on how long these phonons live, this energy gets dissipated quickly to other phonon modes (short lifetime $\to$ large friction), or not (long lifetime $\to$ lower friction). Hence, frictional control tactics could be aimed at tuning the lifetimes of vertically propagating phonons. In situations where the phonon lifetime is long (i.e. our harmonic substrate) using backreflections in this way, could potentially reduce friction by over 20% [Fig. 7(a)].

We conclude with a word about system size effects. The present study used substrates containing $N=1600$ particles. Increasing $N$ will introduce new phonon modes into the system, which, depending on how these couple to the sliding bead, will likely affect friction (in fact, finite size effects in friction have been reported Benassi et al. (2010); Kajita et al. (2009)). Further investigation of the interplay between system size and friction, indeed, how this might be exploited as control tactic, could be a topic for future work.

Acknowledgements.

We acknowledge support by the German research foundation (SFB-1073, TP A01).

Appendix A Computation of phonon properties

Central in the analysis of phonon properties is the dynamical matrix, which, for the 2D hexagonal lattice, is a $2\times 2$ matrix Cserti and Tichy (2004); Jones and March (1985)

[TABLE]

with unit vectors $\hat{n}_{\nu}=(\cos\theta_{\nu},\sin\theta_{\nu})$, $\theta_{\nu}=2\pi\nu/3$, $m$ the substrate particle mass, and $\alpha_{2}$ the coefficient of the quadratic term in the bond energy. We choose the wavevectors on a grid inside the FBZ, $\vec{k}(n_{1},n_{2})=(n_{1}\vec{b}_{1}+n_{2}\vec{b}_{2})/n_{\rm max}+{\rm mod}(l_{s}\vec{d}_{i})$, integers $0\leq n_{i}<n_{\rm max}$, $n_{\rm max}=\sqrt{N}$, with $\vec{b}_{i}$ the reciprocal lattice vectors (which follow trivially from the real space primitive vectors $\vec{a}_{i}$ given in the main text). The modulo operation corresponds to repeated translations of length $l_{s}=4\pi/\sqrt{3}a$ (which is the height of the FBZ hexagon) in the directions $\vec{d}_{i}=\pm(\cos\phi_{i},\sin\phi_{i})$, $\phi_{i}\in\{-\pi/6,\pi/6,\pi/2\}$, until $\vec{k}$ lies inside the FBZ [Fig. 1(b)]. For each wavevector $\vec{k}$ inside the FBZ, diagonalization of $\tilde{D}(\vec{k})$ yields two eigenvalues, $\lambda_{p}(\vec{k})$, with corresponding eigenvectors $\hat{e}_{p}(\vec{k})$, normalized to unit length, where $p$ denotes the polarization. The eigenvalue yields the mode frequency, $\omega_{p}^{2}(\vec{k})=\lambda_{p}(\vec{k})$; upon inspecting the inner product, $\vec{k}\cdot\hat{e}_{p}(\vec{k})$, one finds that the mode with the lowest frequency is predominantly transversal, the other longitudinal, and so we set $p\in\rm TA,LA$ in our notation. During the MD simulations, we record, for each particle $i=1,\ldots,N$, the displacement $\vec{u}_{i}(t)$ from its (time-averaged) position, and its velocity $\vec{v}_{i}(t)$, as function of time $t$. For this analysis, $\sim 20000$ time measurements were taken, each one separated by 10000 (1000) MD timesteps for the harmonic (anharmonic) substrate. From these data, we compute, for each wavevector $\vec{k}$ and polarization $p$, the normal mode coordinates at time $t$, $Q_{p}(\vec{k},t)=(\sqrt{m}/N)\sum_{i=1}^{N}\vec{u}_{i}(t)\cdot\hat{e}_{p}(\vec{k})\,e^{-\imath\vec{k}\cdot\vec{R}_{i}}$, with an analogous expression for $\dot{Q}_{p}(\vec{k},t)$, where one replaces $\vec{u}_{i}(t)\to\vec{v}_{i}(t)$. In these equations, $\vec{R}_{i}$ is the perfect hexagonal lattice position of particle $i$ at the start of the simulation (i.e. does not depend on $t$). The normal mode coordinates are then converted to mode amplitudes:

[TABLE]

Note that the normal mode coordinates and amplitudes are generally complex numbers. The instantaneous phonon energy is given by $E_{p}(\vec{k},t)=\omega_{p}(\vec{k})|A_{p}(\vec{k},t)|^{2}$, which can be time-averaged over the MD trajectory to obtain $E(\vec{k},p)$. The phonon lifetime $\tau(\vec{k},p)$ is obtained from the autocorrelation function $\chi_{p}(\vec{k},t)$ of the corresponding amplitude time series $|A_{p}(\vec{k},t)|$. We use normalization $\chi_{p}(\vec{k},0)=1$; the time-averaged value of $|A_{p}(\vec{k},t)|$ is subtracted from the time series; fast Fourier transforms are used to compute $\chi_{p}(\vec{k},t)$. Precise estimates of $\tau(\vec{k},p)$ are difficult to obtain, since the functions $\chi_{p}(\vec{k},t)$ can be quite noisy. To this end, we first averaged $\chi_{p}(\vec{k},t)$ over all wavevectors, $\bar{\chi}_{p}(t)=(\sum_{\vec{k}}w_{\vec{k}}\chi_{p}(\vec{k},t))/(\sum_{\vec{k}}w_{\vec{k}})$, with “weights” $w_{\vec{k}}=1/|\vec{k}|$. The function $\bar{\chi}_{p}(t)$ decays to zero with increasing $t$ rather smoothly. We take the corresponding rate of the decay as a measure of the typical phonon lifetime, $\langle\tau_{p}\rangle$, shown in Fig. 6, which we compute using the integral measure: $\langle\tau_{p}\rangle=\int\bar{\chi}_{p}(t)\,dt$. To obtain the lifetime $\tau(\vec{k},p)$ of individual modes, we use the same integral measure, but apply it to the envelope function of $\chi_{p}(\vec{k},t)$, which one obtains via a Hilbert transform. The resulting estimates are scaled by a factor $f$ afterward, such that their (weighted) average matches the typical lifetime: $\langle\tau_{p}\rangle=(f\sum_{\vec{k}}w_{\vec{k}}\tau(\vec{k},p))/(\sum_{\vec{k}}w_{\vec{k}})$, with $w_{\vec{k}}$ as given above.