Structural lubricity of physisorbed gold clusters on graphite and its breakdown: Role of boundary conditions and contact lines

The sliding motion of gold slabs adsorbed on a graphite substrate is simulated using molecular dynamics. The central quantity of interest is the mean lateral force, that is, the kinetic friction rather than the maximum lateral forces, which correlates with the static friction. For most setups, we find Stokesian damping to resist sliding. However, velocity-insensitive (Coulomb) friction is observed for finite-width slabs sliding parallel to the armchair direction if the bottom-most layer of the three graphite layers is kept at zero stress rather than at zero displacement. Although the resulting kinetic friction remains much below the noise produced by the erratic fluctuations of (conservative) forces typical for structurally lubric contacts, the nature of the instabilities leading to Coulomb friction could be characterized as quasi-discontinuous dynamics of the Moiré patterns formed by the normal displacements near a propagating contact line. It appears that the interaction of graphite with the second gold layer is responsible for the symmetry break occurring at the interface when a contact line moves parallel to the armchair rather than to the zigzag direction.


Introduction
The origin of friction between surfaces in relative sliding motion has led to speculations since the times of Coulomb (1785). He recognized a possible friction mechanism that surface asperities must deform so that they can slide past each other or that surface atoms, due to their proximity, can form a coherence, which needs to be overcome to initiate sliding. In an extreme point of view, individual atoms could be seen as the smallest possible asperities and the molecular generation of coherence as a (incommensurate) superstructure. When irregularities of a counterbody make surface atoms deflect from their ideal lattice positions, vibrations are excited. As long as they can be treated as a small perturbation of the crystal lattice, the ensuing dissipation can be approximately described within a linear-response theory (Adelman and Doll (1976)); that is, friction would be Stokes-like. However, once the system is driven to an instability point, where, at a given moment of time, the closest stable atomic position is a finite distance away from the previous one, an atom-or a collection of atoms-quasidiscontinuously pops toward the next available equilibrium site (Prandtl (1928)). As already demonstrated by Prandtl himself, such a process leads to Coulomb-like friction, unless the system is thermal and driven so slowly that the relevant degree of freedom can sample both the old and new potential energy minima even before the ultimate instability point, in which case Stokes friction can arise again. Interestingly, the transition between the Coulomb and Stokes regimes in the Prandtl model, which consists of a simple damped spring pulled past a sinusoidal potential, can be well described by the Carreau-Yasuda equation, which Yasuda (1979) derived to describe the shear thinning of polystyrene [Müser (2020)]. Hirano and Shinjo (1990) were arguably the first to come to the conclusion that ideally flat but incommensurate (iron) solids placed on top of each other do not interlock to a degree leading to static friction, implying the absence of kinetic friction at a (infinitesimally) small sliding velocity. They later called the phenomenon superlubricity [ Shinjo and Hirano (1993)]. Their prediction of super-low friction was confirmed experimentally on MoS 2 -coated solids in ultrahigh vacuum: Martin et al. (1993) measured friction coefficients as small as 10 -3 . In theory, the small friction originates from the systematic annihilation of lateral forces, where some interfacial atoms are pushed to the right and others to the left. Thus, the more incommensurate the interface, the smaller the friction. This expectation was confirmed in experiments by Dienwiebel et al. (2004), who rotated two graphite flakes against each other. Friction reduces dramatically with increasing angular misfit [Filippov et al. (2008)]. To distinguish ultrasmall friction between two (incommensurate) solids in direct mechanical contact from that due to an intervening fluid layer, Müser (2004) introduced the term structural lubricity, which can be seen as a special form of superlubricity [Erdemir and Donnet (2006); Baykara et al. (2018)].
If either the substrate or the physisorbed cluster is amorphous, simple scaling arguments based on the law of large numbers suggest root-mean square lateral or friction forces F to increase with A √ rather than linearly with the area of the cluster's basal plane A as long as elasticity keeps the upper hand over interfacial forces [Müser et al. (2001);Müser (2004)]. For incommensurate surfaces, the increase of friction with contact area is generally weaker and may even depend sensitively on the cluster's shape as well as on the sliding direction [de Wijn (2012)]. However, in certain situations, even incommensurate structures can yield (maximum) lateral forces that scale with A √ , for example, when the cluster has a sharp edge that is perfectly aligned with the substrate and the sliding direction is normal to that edge. For a random contact line, Koren and Duerig (2016) from scaling arguments since the annihilation of lateral forces in incommensurate contacts is quite systematic and randomness comes from the outer rim only, while spherical adsorbed solids (noble-gas monolayers) reveal yet another scaling with approximately A 0.37 [Varini et al. (2015)]. Experiments on metal clusters on graphite [Dietzel et al. (2013[Dietzel et al. ( , 2018; Hartmuth et al. (2019)] as well as computer simulations [Chen and Gao (2020)] seem to confirm the scaling arguments. Müser et al. (2001) showed that the scaling arguments can also be motivated from a continuum description of the repulsive forces between interacting surfaces. Marom et al. (2010) introduced a similar parameter for discrete systems, which was called the registry index. Hod (2013) demonstrated that it provides quantitative descriptions for various layered compounds. The abovementioned scaling arguments assume the solid bodies to be rigid. Thus, an important question to be addressed is what makes two crystals with plane surfaces pin and/or what controls kinetic friction. Analysis of low-dimensional models, such as the Frenkel-Kontorova model, can only provide qualitative insight as they ignore the long-range nature of elastic restoring forces [Shinjo and Hirano (1993); Müser (2004)]. Candidates to foster pinning between two incommensurate solids with flat surfaces are adsorbed layers or boundary lubricants [He et al. (1999); Müser et al. (2001); Dietzel et al. (2008);Feldmann et al. (2014)] and (extended) lattice defects. For adsorbed layers to act as pinning agents, normal pressures have to be large enough since the interfacial layer would simply form an intervening viscous medium between the surfaces otherwise; that is, depending on the ratio of energy activation barriers and thermal energy, the response of a molecular layer can range from Stokes-to Coulomb-like friction, as discussed, for example, by Müser (2020). In fact, Özoğul et al. (2017) found superlubric states for metal clusters adsorbed on graphite despite conducting their experiments under ambient conditions. Roughness can also lead to pinning, for example, through many small contact patches that carry relatively little normal load but, due to their small size, exert relatively large frictional shear stresses [Müser (2019)]. Finally, when interfacial interactions dominate the ones inside the bulk, which happens when chemical bonds form across the interface [Dietzel et al. (2017)], two solids or clusters have no choice but to pin.
In this article, we will focus on ultrahigh vacuum (UHV) conditions and idealized systems with perfectly smooth surfaces. Pinning, or rather the onset of pinning under such clean conditions, is a competition between elastic restoring forces and interfacial interactions. The question to be addressed then is which force "wins" at large length scales. For threedimensional amorphous systems, Müser (2004) predicted elastic restoring forces to scale linearly with the linear dimension of a contact so that the details determine which of the two effects keeps the upper hand since lateral and interfacial Frontiers in Chemistry frontiersin.org forces frequently obey the same scaling. This prediction is consistent with the results by Sharp et al. (2016) and Monti and Robbins (2020), who found that large, disordered clusters pin, while small ones do not. Also, Wang et al. (2019)'s results that the friction of thin graphite sheets crosses over from a F ∝ A √ to F ∝ A scaling are consistent with the scaling arguments since the flakes are two-dimensional so that a cross-over to linear scaling is expected.
Despite much progress in the presence or rather the absence of static friction and the absence of instabilities leading to Coulomb-like friction, it is surprisingly unexplored what parameters affect the prefactors to viscous-like damping in a contact satisfying the requirements for structural lubricity. In recent work, Lodge et al. (2016) compared simulations to experimental results on the slip time of gold nanocrystals sliding past the graphene substrate using the quartz crystal microbalance. However, their work does not reveal how damping depends on the velocity or on the geometric features describing the contact.
One aspect, which could be particularly relevant to dissipation at the small scale, is the sliding velocity relative to the contact line. From a continuum perspective, stress gradients are the largest near a contact line so that moving parallel to a contact line would be expected to yield small dissipation and moving normal to it would yield large dissipation. Quantifying this effect is one purpose of this article. In addition, we would like to explore how the elasticity of the solids affects dissipation. Most solids of practical interest can be described as semi-infinite. However, simulations assume a few layers only and keep the bottom layer fixed. This can become problematic when the contact radius is larger than the height of the simulated slab. To assess the role of boundary conditions, we explore the two extreme limits of keeping the bottom-most layer completely flexible or completely rigid, thereby providing lower and upper bounds for the true elastic response of the substrate. Through this manipulation, we also effectively change the dimensionality of the objects. A stiff bottom layer resembles a Winkler foundation; that is, the elastic properties of the solid are those of a high-dimensional object, while a soft foundation mimics the response of a two-dimensional sheet, each time, of course, for undulations at wavelengths clearly exceeding the height of the object.

Model
Using molecular dynamics (MD), we model a six-layer gold slab with the (111) surface sliding against a three-layer highly oriented pyrolytic graphite (HOPG) substrate in the absence of contaminants. Both the contact surfaces are atomically flat and defect-free. As shown in Figure 1, three modeling scenarios are investigated, termed as (1) full coverage, (2) x-stripe, and (3) ystripe. Periodic boundary conditions are applied in the xy-plane so that contact lines appear only in (2) and (3). The in-plane dimensions of the simulation cell are 11.5 × 11.6 nm 2 , to which the strains of the gold slab in (1) as well as those in (2) and (3) along the longitudinal directions are less than 0.3% with respect to its minimum-energy (bulk) configuration. The widths of the stripes in (2) and (3) are approximately one half of the corresponding cell length.
In each case, the center of mass (COM) of the HOPG's bottom layer is kept fixed, while the COM of the gold top layer is constrained to move at a constant in-plane velocity ranging typically from 20~120 m/s at different angles relative to the stripe direction. However, the COM of the gold cluster's top layer is unconstrained in the direction normal to the interface so that it moves under zero external normal load. Atoms in the outermost layers are either constrained ("rigid") relative to the COM of that layer or free to move ("flexible") relative to it according to Newton's equations of motion. From an elasticity point of view, a flexible HOPG bottom layer can be loosely associated with a freely suspended three-layer thick graphite solid. The Snapshots of three simulation models at rest. Yellow arrows indicate the sliding direction, that is, if motion is parallel (||) or perpendicular (⊥) to the contact line. L x and L y are 11.5 and 11.6 nm, respectively.
Frontiers in Chemistry frontiersin.org elasticity of a semi-infinite graphite would be softer than that of the rigid bottom layer but stiffer than the flexible layer. The rigidity assumption for gold is made to mimic the constraining effect that an atomic-force microscope tip moving a cluster over the surface might have on that cluster. Again, reality will be somewhere in between the ideal, limiting cases of perfect rigidity and flexibility. Neither forces nor torques or any other external constraint other than those on the COM velocities of the two outermost layers are imposed. The interatomic interactions of graphite are described by the AIREBO potential developed by Stuart et al. (2000), while the interactions of gold are described by an embedded-atom potential proposed by Zhou et al. (2004). The cross interactions between the two substances are described by a Morse potential [de la Rosa-Abad et al. (2016)] with a cutoff of 10 Å. To improve the signal-to-noise ratio, thermal noise is kept small by setting the target temperature to 1 K using a Langevin thermostat. It is only applied to the mid-graphite layer in the z-direction with a coupling time constant of 1/γ = 100 fs. The simulation timestep is 1 fs. All the simulations are carried out using the open-source MD code LAMMPS [Plimpton (1995)].
Shear stresses and (lateral) forces, whose averages are frictional forces, were measured in various ways and different locations in the system, that is, the total force acting on either the bottom-most graphite layer or the top-most gold layer as well as the interfacial forces acting between the gold and the carbon atoms. All these forces must be identical on an average during steady-state sliding, except for their sign, due to Newton's third law. We confirmed this to hold within the small scatter, which is due to finite run time effects. Averages were typically accumulated over a sliding distance of 2 μm. In addition, the sliding-induced dissipated power P γ i∈mid (m〈v 2 iz 〉 − k B T) adsorbed by the thermostat was averaged and also successfully correlated with the directly measured friction forces by using the equation P = 〈F〉·v. Here, i ∈ mid refers to carbon atoms in the middle graphite layer and m is their mass, while 〈. . .〉 indicates a time average during steady-state sliding.
The shear stress was obtained by dividing the mean force in the sliding direction through the area of the adsorbate, for which each atom in the bottom gold layer was assigned the same area. Statistical error bars, ΔO 2 , of an observable O are deduced from the integral over its time autocorrelation function C OO (t)≔〈O (t′) O (t′ + t)〉 via ΔO 2 2 τ τ 0 C OO (t)dt, where τ is the simulated time. A representative measurement of the lateral force, which was conducted on a y-stripe sliding in the x-direction, is shown in Figure 2. It can be seen that the mean friction force is only a small fraction of typical instantaneous forces, which thus are predominantly conservative in nature and expected for structurally lubric contacts.

Reference cluster and preliminary considerations
Before discussing the results on sliding, it is useful to analyze how the interaction between "indenter", that is, the cluster, and the substrate deforms the solids because similar deformations are dragged along as the cluster is moved with respect to the substrate, whereby surface vibrations are excited, which ultimately propagate toward the bulk and get adsorbed into a heat sink, in our case into the thermostat. Expectations from generic continuum considerations would be as follows: the substrate can increase the interaction with the indenter by moving (near-surface) atoms below the indenter. Such a process would increase the number of atoms below the indenter, which would make the substrate raise up below the indenter and go down right outside the contact line-assuming the Poisson's ratio of the substrate to be positive. For our system, the situation is qualitatively different. Radial displacements are small because of the stiff in-plane bonds in graphite. This reduces the propensity of carbon atoms to be pulled below the indenter. Moreover, due to the interactions being body rather than surface forces, the second graphite layer is attracted toward the gold cluster. As a consequence, the lattice contracts below the indenter, as is evidenced in Figure 3, panels (b) and (d). The displacements in the normal (z) direction turn out to be a factor of 3 larger than in the in-plane (y) direction, as revealed in Figure 3B.
Although a continuum description of the x-stripe and ystripe would be identical, a clear difference between their normal displacements shows up. Specifically, the z-displacements for the y-stripe reveal Moiré patterns, which clearly violate the inversion symmetry, while all other displacements obey the symmetry Typical dependence of the instantaneous force/stress (black lines) on slid distance for the y-stripe at a sliding velocity of v ⊥ = 50 m/s. Averages were only taken over the last 2 μm.
Frontiers in Chemistry frontiersin.org  Frontiers in Chemistry frontiersin.org 05 expectations, except for minor fluctuations, which are unavoidable in discrete/atomic systems. However, the symmetry breaking for the z-displacements in the y-stripe do not necessarily reflect broken ergodicity since the contact is not mirror-symmetric about the yz-plane. In fact, it is the second gold layer, which breaks the mirror symmetry for a y-stripe as can be seen in Figure 4A. The broken symmetry is also revealed in the normal displacements, as depicted in panel (b).

Sliding simulations
Scaling arguments for the magnitude of instantaneous or maximum lateral forces usually focus on the contact area for amorphous clusters and on the circumference for incommensurate interfaces between crystals. However, this does not mean that the cluster height or the boundary conditions or constraints on the outermost layers are irrelevant. Assuming two clusters to have identical basal planes, the thinner cluster is more compliant (or "flexible") than the thicker one, which is effectively more "rigid". This regards both in-plane and out-of-plane deformations. A larger in-plane compliance in the direction of sliding (i.e., in the longitudinal direction) implies increased friction since the cluster can interlock more easily with the counterbody. However, an increased out-of-plane-as well as in-plane, transverse-compliance allows the surface atoms to deflect away from irregularities on the counterface so that a larger cluster height can also increase friction. To investigate which of the effects is more relevant for metal clusters adsorbed on graphite, Figure 5 depicts the shear stress as a function of velocity for (a) perfectly rigid and (b) flexible outermost layers for a ystripe sliding in the x-direction. Figure 5A reveals that the motion of both stripes and the full layer show Stokes-like damping, roughly linear in velocity, when both the outermost layers are rigid. The damping coefficient for the x-stripe roughly equals that of the full layer, while the y-stripe is damped a little less than twice as strongly. The increased damping for the y-stripe can be easily rationalized since its velocity is perpendicular to the contact line. From a continuum perspective of brittle fracture, a (small) hysteresis between a closing crack at the leading edge and an opening crack at the trailing edge must be expected, which adds to the bulk or areal dissipation far away from the crack Popov et al. (2021).
After releasing the rigidity constraint at both outermost layers, the lateral force opposing the motion of both the xstripe and full layer parallel to x still is Stokesian with a  Interestingly, the y-stripe reveals a much increased resistance to sliding and a rather weak velocity dependence of the kinetic friction on sliding velocity at small v. This increased friction is related predominantly to the increased compliance of the flexible graphite substrate, as can be seen from Figure 6. Figure 7 highlights the interplay of sliding direction and mean lateral force in the case of flexible outermost layers, which we find to be the most interesting case, since it can deviate from Stokes damping. The full-coverage layer does not show any detectable direction dependence. The x-stripe reveals a damping coefficient, which increases with increasing angle between the contact line and sliding direction, while the ystripe transitions between Stokes-and Coulomb-like friction at an angle between 30°and 60°.
Instabilities are required in order for Coulomb-like friction to occur [ Prandtl (1928); Tomlinson (1929); Müser (2002)]. Movies of the displacements were produced to characterize the instabilities, that is, animated versions of all panels shown in Figure 3. They clearly reveal that the normal displacement of the y-stripe moving in the x-direction displays quasidiscontinuous dynamics, while all other displacement fields evolve continuously with time. Snapshots of these movies are shown in Figure 8, from where it becomes apparent that the observed Moiré patterns only changed marginally in the relatively large time periods separating panels (a) from (b), (c) from (d), and (e) from (f), but quite substantially during the brief time periods separating configuration (b) from (c) and (d) from (e).
The dynamics loosely resembles that of a crack front, which discontinuously advances by a distance of the order of one atomic  spacing (Wang et al., 2020). However, for the studied case, different types of instabilities occur, where the corrugated pattern does not only advance by a small multiple integer of the graphite lattice constant but simultaneously, for example, in the transition from Figure 8 panel (b) to (c), but also undergoes a phase shift of 180°from (d) to (e). As is always the case for Coulomb friction, the energy dissipated in the process is approximately the difference of the closest energy minimum of the configuration just before and just after the pop. It is relatively insensitive to the rate with which locally released energy or heat is transported away from the interface, unless the heat conductance is so small that the interface heats up substantially. Finally, the size dependence of the Coulomb shear stress is investigated in Figure 9. It decreases with the length of the y-slab much more quickly than simply with 1/L x . Specifically, the decay is consistent with an exponential dependence on the relatively small investigated domain. This means that the contribution of the near-contact line zone cannot be simply assigned a unique value, but it appears to be a non-local contribution.

Conclusion
In this article, we studied the kinetic friction of the adsorbed, periodically continued gold stripes, sliding past a graphite substrate under their mutual adhesive attraction. We found the systems to exhibit structural lubricity for the most part in that mean lateral forces disappeared quickly, that is, linearly with velocity. When atoms in the outermost layers of the graphite substrate and the gold adsorbate were constrained to zero displacement relative to their equilibrium position in the respective outermost layer, Stokes friction was observed in all cases. When sliding a full-coverage cluster kn, any direction or a finite-width slab parallel to the direction of the stripe, Stokes stress, defined as Stokes friction per unit area, turned out to be roughly σ 20 m/s = 1.3 kPa at a sliding velocity of 20 m/s. This is the same as the shear stress that would be obtained with a~1.5 μm thick film of the lubricant as water under ambient conditions, η = 0.1 mPa s, and stick boundary conditions. In these "Stokesian cases" of our study, dissipation can be related to quasi-elastic deformations of the two solids in contact and the unavoidable hysteresis, which is consistent with linear harmonic theory, as proposed, for example, in the pioneering study by Adelman and Doll (1976) or the more modern formulation by Kajita et al. (2010) using the timedependent elastic Green's functions.
Allowing the outermost layers to be flexible increased the friction in all cases, but most so for the y-stripe when sliding parallel to the armchair direction of the graphite lattice, in which case Coulomb friction occurred. Although the shear stress also increased for the x-stripe, in fact by a factor of 4 when sliding perpendicular rather than parallel to the stripe, it remained Stokesian. This latter behavior can be rationalized to be a consequence of a slowly moving adhesive crack pair, in which the crack closure on the leading edge does not fully recuperate the energy needed to open the crack at the closing edge. When elastic or phonon relaxation times are small, this process leads to a hysteresis linear rather than algebraic in velocity [Müser and Persson (2022)]. When moving parallel to the slab, no contact-line contribution was observed, which, however, can be due to the fact that both the solids had a very smooth and unjagged edge.
For flexible outermost layers, which in the case of the modeled graphite can be associated with a suspended substrate, the Coulomb friction was induced by the added compliance. This result is in agreement with the pioneering study of Lee et al. (2010) on the friction of suspended layered solids, including graphite, where friction was observed to continuously decrease with width, which is in line with the argument that long-range elastic restoring forces are a key element for superlubricity [Shinjo and Hirano (1993); Müser (2004)]. Lee et al. (2010) related the increased friction to the easier out-of-plane puckering for thin sheets. We have little reason to object to that conclusion but refine the mechanism for our system in that puckering occurs not only in front of the indenter but also below it. Moreover, the induced undulations below the indenter near the contact line have Moiré pattern characteristics and move continuously for the most part with occasional quasi-discontinuities at isolated moments of time. These Moiré patterns differ from those considered earlier in the context of friction, for example, Figure 1 in Ref. [He et al. (1999)] or the Moiré pattern analysis by Müser et al. (2001) in Area dependence of shear stress for the xand y-stripe geometries and flexible outermost layers at v ⊥ = 50 m/s. Here, the size of the substrate was kept constant and the length of the slab was varied parallel to the sliding velocity direction. The blue dashed line shows an exponential dependence with L ⊥ with σ(L ⊥ → ∞) constrained to the shear stress obtained at full coverage (≈ 5 kPa).
Frontiers in Chemistry frontiersin.org 08 that normal rather than in-plane displacements form the Moiré pattern. Placing our results into a historical context, our Moiré pattern instabilities are rather a Fourier version of the discontinuous dynamics suggested by Tomlinson (1929), who considered instabilities normal to the interface, rather than those by Prandtl (1928), who focused implicitly more on motion within atoms within their planes.
Despite some commonality with experiments, the shear stresses in our system are extremely small, that is, about 5 kPa when the linear size of the stripe extends over 25 nm, which is roughly 1 order of magnitude less than typical peaks in lateral forces or static friction. In most other situations, instabilities are required to lead to detectable friction; that is, the friction in our systems would still have to be classified as superlubric by all means, yet we would not label them as structurally lubric due to the presence of multistability at scales smaller than the width of the stripes. While it remains to be seen to what extent our puckering-Moiré-pattern mechanism matters in real (UHV) systems, our study certainly revealed that utmost care has to be taken when modeling the friction of graphite layers. It is also a paradigm support for Coulomb's ingenious insight of the origin of friction, reiterated in the very beginning of our introduction, except that in our case, the roughness that needs to be overcome to continue sliding is not pre-existing but the consequence of the coherence that graphite and gold want to form due to their proximity.

Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions
HG set up, run, and analyzed the simulations and edited the manuscript. MM designed the study, assisted in the analysis of the simulations, and wrote the manuscript.

Funding
This research was founded by the German Research Foundation (DFG) through grant number GA 3059/2-1.