Size-dependent commensurability and its possible role in determining the frictional behavior of adsorbed systems

Recent nanofriction experiments of xenon on graphene revealed that the slip onset can be induced by increasing the adsorbate coverage above a critical value, which depends on temperature. Moreover, the xenon slippage on gold is much higher than on graphene in spite of the same physical nature of the interactions. To shed light on these intriguing results we have performed molecular dynamics simulations relying on ab initio derived potentials. By monitoring the interfacial structure factor as a function of coverage and temperature, we show that the key mechanism to interpret the observed frictional phenomena is the size-dependence of the island commensurability. The latter quantity is deeply affected also by the lattice misfit, which explains the different frictional behavior of Xe on graphene and gold.


I. INTRODUCTION
Submonolayer islands of rare gas atoms adsorbed on crystal surfaces offer an excellent platform to address friction at crystalline interfaces. In the submonolayer range (0 < θ < 1, where θ is the coverage) and at low temperatures, adsorbate phase diagrams are well known to display phase-separated two-dimensional (2D) solid islands, usually incommensurate with the surface lattice, coexisting with the 2D adatom vapour [1]. The inertial sliding friction of these islands can be accurately measured using a quartz crystal microbalance (QCM) [2]. The condensation of a film on the QCM electrodes is signaled by a decrease in the resonance frequency. Any dissipation taking place at the solid-film interface is instead detected by a decrease in the corresponding resonance amplitude. From the shifts in the resonance frequency and amplitude of the QCM one can calculate the slip time τ s [3]: it represents the time constant of the exponential film velocity decrease due to a hypothetical sudden stop of the oscillating substrate. Very low τ s means high interfacial viscosity and, in the case of a film rigidly locked to the substrate, τ s goes to zero.
Many different friction phenomena of sub-monolayer islands have been studied using the QCM, including: the dynamic depinning of Kr on gold [4], the impact of substrate corrugation on the sliding friction levels of adsorbed films [5], the existence of an onset coverage for sliding of Ne islands on lead [6], the dynamical sticking of solid helium on graphite [7], the nanofriction of various adsorbates on superconducting lead [8][9][10], and the superlubricity of Xe islands on copper [11].
The nanotribological studies based on the QCM, which involve clean, well defined interfaces, represent an ideal situation where simulations can combine with experiments to unravel fundamental aspects of friction. Recent * mcrighi@unimore.it findings include the scaling behavior of the island edge contribution to static friction, [12] the effect of non linear dynamics, [13] the influence of the adsorbate-substrate interaction on friction, [14] the size dependence of static friction, [15] the importance of the domain nucleation, [16] the nature of thermally activated island creep [17] and, more in general, the energy dissipation in monolayer films sliding along substrates. [18][19][20] The nanotribological studies based on the QCM, which involve clean, well defined interfaces, represent an ideal situation where simulations can combine with experiments to unravel fundamental aspects of friction. Recent findings include the scaling behavior of the island edge contribution to static friction [12], the effect of non-linear dynamics [13], the influence of the adsorbate-substrate interaction on friction [14], the size dependence of static friction [15], the importance of the domain nucleation [16], the nature of thermally activated island creep [17] and, more in general, the energy dissipation in monolayer films sliding along substrates [18][19][20].
Nanoparticle manipulation [21] is also emerging as an invaluable tool to investigate the frictional properties of nanoaggregates on surfaces, both experimentally and by theoretical simulations. Recent findings, complementary and closely related to the QCM ones, include the sizescaling of structural lubricity [22][23][24], the influence of contact aging [25], and thickness promoted lubricity for metallic clusters [26].
Recently, Pierno et al. have measured the sliding time of Xe on graphene using a QCM technique between 25 and 50 K [27] and found that the Xe monolayer is pinned to the graphene surface at temperatures below 30 K. At T = 35 K, the Xe film starts to slide for a coverage of 0.45 monolayer (ML). The critical coverage for sliding decreases with increasing temperature afterwards. By comparison, under similar conditions, the slip time of Xe on the bare gold surface is about twice higher than that on graphene, showing a lower depinning temperature of 25 K. This is counterintuitive since it is believed that the arXiv:1610.06323v2 [cond-mat.mes-hall] 25 Jan 2019 hexagonal graphene layer is smoother for adsorbate than gold surfaces. The understanding of the physical mechanisms behind this puzzling observation is the subject of this theoretical investigation, in which we perform a comparative study of xenon on graphene and gold substrates. By means of classical molecular dynamics, relying on ab initio derived potentials, we simulate the particle diffusion on the two substrates and monitor the commensurability of the growing islands as a function of coverage and temperature.

II. SYSTEMS AND METHODS
At full monolayer (ML) coverage and temperature below 70 K, the xenon adatoms form a √ 3 × √ 3 R30 • sublattice on the graphene layer [28] (GL) (left panel of Fig. Fig. 1). The mismatch between the lattice parameters of an isolated Xe monolayer and the graphene √ 3 × √ 3 R30 • sublattice is very small, equal to −2.8%. The same adsorption configuration for xenon-on-gold (Xe/Au) (right panel of Fig. 1) is characterized by a much larger mismatch of about +13.9%.
We sample the potential energy surface (PES), which describes the interaction between the Xe adlayer and the substrate as a function of their relative position, by means of ab initio calculations based on density functional theory (DFT). The results are used to optimize the parameters of the interaction potentials subsequently employed in molecular dynamics simulation systems with a large number of atoms. We perform DFT calculations using the Quantum Espresso package [29], initially within the local density approximation (LDA) for the exchange correlation functional. Afterwards, to better describe the van der Waals (vdW) interactions, the numerical approach has been improved by using the rVV10 density functional [30]. rVV10 denotes a more efficient version of the original VV10 scheme [31] and is based on a nonlocal correlation functional which provides an accurate description of van der Waals effects. rVV10 is expected to perform better than, for instance, PBE-D [32] which is instead a scheme where the semilocal PBE functional is modified by simply adding an empirical potential. This is particularly true for metal and semimetal systems, such as the substrates considered in our study, where the electronic charge density is relatively delocalized, so that empirical, atom-based vdW corrections turn out to be inadequate.
The ionic species are described by pseudopotentials and the electronic wavefunctions are expanded in plane waves. For the Xe/GL (Xe/Au) system a kinetic energy cutoff of 30 Ry (25 Ry) is used to truncate the expansion on the basis of preliminary calculations optimizing the lattice parameter of graphene (Au). We use periodic supercells with ( √ 3 × √ 3) in-plane size and vertical edge 42Å (46Å) long. The k-point sampling of the Brillouin zone is realized with a 8×8×1 (6×6×1) Monkhorst-Pack grid. [33] The interaction energy per particle, E, for a Xe layer on the substrate is calculated as , where E tot is the total energy of the adsorbate system, E sub is the energy of the isolated substrate and E Xe is the energy per particle in the isolated xenon layer. All these systems are described with the same ( The supercell used to model the Xe/GL system contains N = 1 Xe atom, while the supercell used to model the Xe/Au system contains N = 2 Xe atoms adsorbed on the two opposite surfaces of the slab at sites with the same symmetry on the opposite sides. This reduces spurious dipole-dipole interactions between periodically repeated images. The Au slab is 7 layer thick. Classical molecular dynamics simulations (MD) are performed by using the LAMMPS [34] computer code. The adopted force fields are constructed by putting together the following set of interactions. The xenon-xenon interaction is described by the pair potential proposed by Tang and Toennies [35] where r is the distance between two xenon atoms, f 2n is a damping function for the attractive part of the potential that can be expressed in terms of the incomplete gamma . The adopted numerical values for the A T , b T and C 2n coefficients are the same as in the original paper [35]. This potential provides an improved description, with respect to the Lennard Jones potential, of both the repulsive part of the interaction (with the use of the exponential function instead of the too steep r −12 term) and the attractive part (including the higher order dispersion r −8 and r −10 terms in addition to the standard r −6 one).
The xenon-carbon interaction is described by the func-tional form of the Buckingham potential [36] where the A B , b B and C B coefficients are optimized in order to reproduce the ab initio PES, as detailed in Section III A. A distance cutoff of 12Å is used for both the Tang-Toennies and the Buckingham potentials.
The interactions among the carbon atoms of graphene are modelled by the second generation REBO potential [37,38]. To avoid large out-of-plane deformations of the graphene layer and mimic the experimental conditions, we model the presence of a (rigid) substrate underlying the graphene layer by means of the following analytic potential: where (x, y, z) are the coordinates of a carbon atom on the substrate, The numerical values for the a and C min/max i coefficients are reported in Ref. 39. Finally, the interaction of xenon with the Au(111) surface is modeled by the three-dimensional periodic function [40] V (x, y, z) y) for i = 1, 2, with u(x, y) = 3 − g cos g · r and the summation running over the first three g vector of the reciprocal lattice and A 0 (x, y) is given by (5) This function is able to accurately describe the peculiar features of the PES for rare gases on metals [40,41]. In particular, its "anticorrugation", which originates from the on-top site preference for rare gas adsorption rather than the hollow site [42,43]. This feature cannot be reproduced by pair-wise potentials, such as the Lennard-Jones potential, which favor adsorption at the highest coordinated sites. The numerical parameters entering the A i (x, y) periodic functions are determined by fitting ab initio adsorption energies calculated for several configurations of the xenon atom on the metal surface, as described in the next section.
MD calculations are performed using periodic boundary conditions. The orthorhombic cell used to model the Xe/GL (Xe/Au) system has 15.25 nm × 14.67 nm × 2.00 nm (18.17 nm × 17.48 nm × 2.004 nm) dimensions. We consider three different Xe coverages on both the GL and Au(111) surfaces. The adatom coverage is calculated assuming that y = 100% corresponds to a complete Xe ML with an areal density of 5.94 atoms per nm 2 , i.e., considering a nearest neighbor distance among Xe atoms of 0.44 nm. This choice is consistent with that adopted in the analysis of the QCM experiments [27]. The simulated coverages y = 11%, 22%, and 44% correspond, thus, to 147, 294, and 588 adparticles in the Xe/GL system, and to 208, 416, and 832 particles in the Xe/Au system, 1328 and 1886 being the total numbers of atoms necessary to realize a 100% coverage of the two-dimensional cells used in the two systems.
The MD simulations are started from initial conditions constructed by assigning to the xenon atoms random, non-overlapping, positions sampled from a uniform distribution on the surface plane in the MD cell, and random velocities sampled from the Maxwell distribution at the temperature T o = 25 K. The system is then thermalized at T = 30 K with a constant-volume, constanttemperature (NVT) run, 10 ns long. A second, higher temperature value, T h = 50 K, is also considered to match the experimental conditions, where both T and T h are considered. The initial configurations at T h are generated from those at the end of the T runs, by increasing sharply the temperature and letting the systems equilibrate for further 10 ns. The temperature is controlled by means of Langevin thermostats, [44] with two separate thermostats used for the Xe/GL system, one for each atomic species, and only one for the Xe/Au system. The same integration time step δt = 1 fs is used in all calculations.
The MD simulations are started from the initial conditions constructed by assigning to the xenon atoms random, non-overlapping, positions on the surface, and random velocities sampled from the Maxwell distribution at the temperature T o = 25 . The system is then thermalized at T = 30 K with a constant-volume, constanttemperature (NVT) run, 10 ns long. A second, higher temperature value, T h = 50 K, is also considered to match the experimental conditions, where both T and T h are considered. The initial configurations at T h are generated from those at the end of the T runs, by increasing sharply the temperature and letting the systems equilibrate for further 10 ns. The temperature is controlled by means of Langevin thermostats [44], with two separate thermostats used for the Xe/GL system, one for each atomic species, and only one for the Xe/Au system. The same integration time step δt = 1 fs is used in all calculations. We calculate the Xe adsorption energy on graphene and gold substrates for different relative positions, obtaining in this way a sampling of the PESes for the two adsorbate systems shown in Fig. 1. We firstly perform standard DFT calculations with the exchange correlation functional described by LDA. Then, we take into account the vdW interactions by means of the rVV10 method [30], as motivated in the previous section. The results obtained within rVV10 are used as a benchmark to derive the parameter values for the xenon-carbon (Eq. 2) and xenon-Au(111) (Eq. 4) potentials used in the MD simulations.
The calculated interaction energy, E, between the Xe adatom and the substrate (see Section II for the definition of E) is represented as a function of the adatom-substrate separation, z, in Fig. 2 and b. The Xe atom position relative to the substrate corresponds to the most favorable one in both the cases. By comparing the red and green curves, it appears evident that the DFT-LDA description underestimates the depth of the potential wells, a problem which is solved by the rVV10 method that increases the strength and widens the range of the attractive part of the interactions. The equilibrium distance, z eq , and the corresponding energy in the minimum, i.e., the adsorption energy E ads , are reported in Table I both for the Xe/GL and Xe/Au systems. It can be seen that the absolute value of the adsorption energy of Xe on the gold substrate is significantly higher, about 100 meV per atom, compared with that on the graphene substrate; this is most likely due to the partial hybridization of Xe 5p orbitals with Au d states [45,46]. The equilibrium adsorption energies and distances are in good agreement with previous experimental and theoretical data. In particular, the Xe adsorption energy on graphene, E GL ads = −197 meV per atom, is within the range of theoretical data available in the literature (-128.6 meV per atom by DFT/vdW-WF [47], -142.9 meV per atom by MP2 [48], -204 meV per atom by all-electron full-potential linearized augmented plane wave plus local orbitals method [49] and -209.7 meV per atom by the DFT/vdW-DF method [46]), and the adsorption distance, z GL eq = 3.61Å, is in excellent agreement with the experimental value of 3.59 ± 0.04Å. In the case of the Xe/Au system, our adsorption energy, E Au ads = −289 meV per atom, is in good agreement with previous theoretical calculations (-262.3 meV per atom by DFT/vdW-DF [46]) and within the experimental range [50][51][52][53], while the equilibrium distance, z Au eq = 3.43Å, is slightly shorter than the experimental (3.62Å) and theoretical (3.58Å) values reported in the literature [46,50].
We then calculate the Xe-substrate interaction energy for different relative lateral positions. In particular, we consider the configurations along the direction (referred to as the y direction) passing through the high symmetry sites labelled in Fig. 1. The Xe-substrate distance is optimized at each location. The energy profile obtained for the Xe/GL system, represented in Fig. 2c, presents the minima (maxima) at hollow (on-top) sites. In contrast, the energy profile obtained for the Xe/Au system presents the minima (maxima) at on-top (hollow) sites, in agreement with experimental observations. The energy difference between the maxima and the minima (∆E), referred to as the PES corrugation, is reported in Table I. The PES corrugation calculated within both the DFT schemes is slightly higher for the Xe/Au system than for Xe/GL. The order of magnitude of the PES corrugations is the same as those estimated from experiments. This is consistent with the fact that Xe presents a stronger binding on gold than on graphene. It is, in fact, often observed that the energy barriers to displace an adsorbate (or a countersurface) along a substrate increase with the strength of the adsorbate-substrate interaction.
The results obtained within the DFT-rVV10 scheme for the adsorption of Xe on both graphene and gold substrates appear to be accurate enough to be used as the reference dataset for tuning the unknown parameters appearing in the force field for the Xe-C (Eq. (2)) and Xe-Au (Eq. (4)) interactions. By fitting the data shown in Fig. 2a and b using a nonlinear least-squares Marquardt-Levenberg algorithm [56], we obtain the parameter values reported in Table II. The fitting is very satisfactory as can be seen in Fig. 2a and b by comparing the Xe adsorption energy as a function of the adatom-substrate separation calculated within the MD scheme (in black color) with the ab initio rVV10 data (in green color); in particular the equilibrium energies and distances are in excellent agreement. The energy variation as a function of the adatom lateral position, shown in Fig. 2c and d, is well reproduced by the adopted force fields. An ac-curate description of the lateral PES is essential for a reliable simulation of tribological systems. The shape and the corrugation of the PES determine, in fact, the frictional forces. It is typically difficult to correctly reproduce the shape of the PES for rare gases on metals by analytical force fields. Pair wise potentials, like the Lennard-Jones potential, favor, in fact, the sites with the highest coordination. Experiments reveal, instead, that for many rare gas atoms on metals, such as Xe on gold, the adsorption occurs on-top. Such hard substrates can be accurately modeled by an external potential representing the interaction with a fixed and rigid triangular lattice frame [12,40]. The potential we adopted to model the Xe/Au system correctly reproduces the PES shape (Fig. 2), with minima at on-top sites and maxima at hollow sites. Not only the potential corrugation is in excellent agreement with that obtained by the rVV10 method. A small underestimation, less than 1 meV, is instead observed in the PES corrugation for the Xe/GL system. The uncertainty in the estimation of the PES corrugation is not expected to affect the results of the MD simulations presented in the next section. The different lattice match of the Xe/Au and Xe/GL plays, in fact, a major role in determining the degree of commensurability of the rare gas islands in the two systems.
In conclusion, the parametric potentials tuned on the DFT-rVV10 data provide an accurate description of the adsorbate systems. Their use in MD simulations will allow us to study the evolution of a system relevant for tribology, containing a large number of xenon atoms over long time scales, typically inaccessible by ab initio MD [57,58].

B. Monitoring the commensurability of the adsorbed monolayers by molecular dynamics
Classical molecular dynamics simulation runs, each 20 ns long, for three coverage values and two different temperatures have been analyzed in parallel for the Xe/GL and Xe/Au systems. The final atomic configurations of all the simulated systems are pictured in Fig. 3 using a range of gray shades from black to white to illustrate the different level of commensurability between the xenon layer and the substrates below: black color means that the adatom is commensurate with respect to the substrate (i.e. resides in a minimum of the potential); white color means that the adatom is fully incommensurate (i.e. is located in a maximum of the potential), and gray levels represent intermediate conditions.
In both the studied systems, it is possible to notice that the xenon atoms at low temperature and low coverage gather into many separated islands of small size with the characteristic hexagonal pattern. In the Xe/GL case there is a predominance of black as the Xe atoms tend to occupy the PES minima; in this case the mismatch of the lattices is small and the number of neighbors of each Xe atom is not enough to provide the energy gain TABLE I: Optimized adsorption distance, z eq (Å), adsorption energy, E ads (meV per atom), and potential corrugation, ∆E (meV), for xenon on graphene and gold. The DFT results obtained within the LDA and rVV10 approaches are compared to other theoretical and experimental data present in the literature.   to promote incommensurability between the Xe and GL lattices. For the Xe/Au system, instead, while the size of the island is comparable, the presence of black Xe atoms is much less relevant.
With increasing temperature, one can see immediately that the size of these islands increases in all cases, as the increase of thermal energy makes the Xe atoms to easily overcome the energy barriers of the PES and Xe diffusion on the surface is correspondingly enhanced. The grayscale patterns do not change significantly for the Xe/Au systems, but for the Xe/GL systems one can immediately notice the presence of a lesser number of neighboring black Xe atoms.
As expected, increasing the coverage leads to an increase of the size of the clusters. In particular, while at 11% and 22% coverages one can count several isolated ones, at 44% for both systems Xe atoms are arranged in a unique connected cluster, i.e., the size of the island reaches the dimension of the MD cell. In the Xe/Au case there are mainly isolated black Xe atoms occupying the PES minima while for Xe/GL several, somewhat smaller, black patches remain.
Overall one sees that the size of these xenon clusters increases with both coverage and temperature and, correspondingly, their structure becomes more and more incommensurate with the substrate. The change in the range of temperature and coverage explored is more significant for the Xe/GL system, where the mismatch between the Xe and the substrate is smaller, favoring the existence of small commensurate islands at lower temperatures and coverages. This qualitative behavior of the adatoms can be confirmed quantitatively by computing the structure factor S(G) of the xenon atoms where G is the summation of the two reciprocal lattice basis vectors of the substrate cell and r j is the real space position of the j-th xenon atom. S(G) crosses over from ∼ 0, for fully incommensurability, to 1 for fully commensurability. In panel (a) of Fig. 4 we report the behaviors of each individual S(G) as a function of time, together with the averaged values of S(G) over the last 10 ns of the simulations for both systems. It is possible to notice that, by varying the coverage and the temperature of the systems, S(G) changes significantly; in particular, at low coverage and low temperature, S(G) attains its maximum, instead at high coverage and temperature, S(G) decreases significantly in both systems. Moreover, S(G) is always higher in the Xe/GL system than in the Xe/Au one. These data allow us to confirm quantitatively the qualitative analysis based on the final atomic configurations. Our results complemented with a previous study, where we described the close relation between the commensurability and static friction of adsorbed islands [15], can provide an explanation for the QCM experimental observations reported in Ref. [27]. In particular, (i) the existence of pinning forces in nominally incommensurate systems at low coverage can be accounted for by the result that small adsorbed islands are commensurate with the substrate even in the presence of lattice misfit. (ii) The existence of a critical coverage necessary for the film slipping is related to the critical size that the growing Xe islands need to reach for the depinning process to be activated at the considered temperature. (iii) Lower critical coverages are observed for Xe/Au than for Xe/GL in spite of the larger PES corrugation of the former system because the lattice misfit has a larger impact on the commensurability of islands. We showed, in fact, that the critical size for adsorbed islands, R C ∼ 1 e 2 ∆E , depends on the ratio between the potential corrugation, ∆E, and the interparticle interaction strength , and is dominated by the lattice misfit, e [15]. We expect that the above result can be generalized to solid clusters on crystalline substrates: clusters will interlock with substrates if their size shrinks to a critical value. How small is this value depends on the bulk modulus, interfacial corrugation and lattice mismatch. An example of this situation may be provided by AFM tips that quite often display a stickand-slip behavior, due to interlocking, independently of the nominal incommensurability of the tip and the substrate materials.

IV. CONCLUSIONS
In a recent experiment [27], the nanofriction of Xe monolayers deposited on graphene has been measured by means of a quartz crystal microbalance. At low temperatures, the adsorbate was fully pinned to the substrate, and it started sliding as soon as a critical coverage was reached. The critical coverage was found to depend on the temperature, in particular it decreases with the temperature. Similar measurements repeated on bare gold showed an enhanced slippage of the Xe films and a decrease of the depinning temperature.
To shed light on the atomistic mechanism governing the above described experimental observations, we have performed a comparative study of xenon dynamics on graphene and gold substrates at different coverages and temperatures. The first part of our work is devoted to the accurate theoretical description of the Xe-graphene and Xe-Au(111) interactions. Our DFT calculations show that the rVV10 method [30], which allows to treat nonlocal van der Waals interactions from first principles, improves the LDA description and provides adsorption energies and distances in very good agreement with the available experimental data. The DFT rVV10 results are used to optimize the force fields employed in the MD simulations. An excellent fit of the PES for the Xe/GL system is obtained by using the Buckingham pairwise potential for the Xe-C interactions, while the Xe-Au(111) interaction is accurately described by means of the analytic function proposed in Ref. 40 for rare gas on metals.
The MD simulations performed in the second part of our work show that at low temperatures and coverages the Xe atoms deposited on the graphene cluster in small islands. The island size increases with the temperature due to an increased particle diffusion on the substrate. The island size increases also with coverage until the Xe layer percolates and a unique patch of coalesced islands is formed. We monitor the commensurability of islands during their growth by both visual and quantitative methods based on the calculated structure factor. The results uncover the existence of a close correlation between the island size and commensurability: small islands are in register with the substrate, while larger islands are less commensurate. The simulations repeated under equivalent conditions on gold reveal a similar trend, but a much lower commensurability is found for every considered temperature and coverage.
Considering the close relation between the static friction and the interfacial commensurability [15], our results can explain the existence of a critical coverage for the depinning transition and its dependence on temperature as observed in QCM experiments [27]. Furthermore, they confirm the theory on the size dependence of static friction, according to which nominal incommensurate interfaces becomes commensurate below a critical size [15]. The critical size, that depends on the interparticle interaction strength and the potential corrugation, is dominated by the lattice misfit. The latter dependence accounts for the different frictional behavior observed for the Xe/GL and Xe/Au systems.

V. ACKNOWLEDGMENTS
We acknowledge the CINECA consortium for the availability of high performance computing resources and support through the ISCRA-A "Lubric" project.