Sliding Mechanisms in Multilayered Hexagonal Boron Nitride and Graphene: The Effects of Directionality, Thickness, and Sliding Constraints

The interlayer sliding potential of multilayered hexagonal boron nitride (h-BN) and graphene is investigated using density-functional theory including many-body van der Waals (vdW) interactions. We find that interlayer sliding constraints can be employed to tune the contribution of electrostatic interactions and dispersive forces to the sliding energy profile, ultimately leading to different sliding pathways in these two materials. In this context, vdW interactions are found to contribute more to the interlayer sliding potential of polar h-BN than they do in nonpolar graphene. In particular, the binding energy, the interlayer distance, and the friction force are found to depend sensitively on the number of layers. By comparing with the experimental findings, we identify sliding pathways which rationalize the observed reduced friction for thicker multilayers and provide quantitative explanation for the anisotropy of the friction force.

Understanding nanoscale tribology is of fundamental importance for developing novel nanomaterials with desired (electro)mechanical properties [1,2].The main challenge for reliable modeling of tribological properties in nanostructured materials stems from the nonlocal anisotropic polarization that is accompanied by a rather strong contribution of many-body electronic correlations to their cohesion and dynamics [3][4][5].Promising and widely studied materials for nanoscale tribology include multilayered hexagonal boron nitride (h-BN) and graphene [6][7][8][9][10][11], due to their remarkable mechanical and electronic properties.In addition, bulk h-BN and graphite are widely used as solid lubricants.Therefore, the understanding of the mechanism of interlayer sliding of h-BN and graphene is important, both for elucidating their fundamental tribological properties and for applications of these materials in nanomechanics [12][13][14].
It is widely accepted that the interlayer potentials in h-BN and graphene are determined by electrostatic interactions and van der Waals (vdW) forces [9,10,15].However, the relative importance of different factors in determining the interlayer sliding potential is still under debate.By a combination of friction force microscopy and angle-resolved photoemission spectroscopy, Filleter et al. revealed that the friction of bilayer graphene is a factor of 2 smaller than a single layer of graphene on SiC, which was attributed, however, to electron-phonon coupling [16].In comparison, Lee et al. observed that friction is monotonically decreased as the number of layers increase from a single layer to four layers for h-BN and graphene, which was attributed to the local puckering of the substrate due to adhesion with the tip [7].
Many efforts have been made to understand the interlayer interactions of h-BN and graphene using density-functional approximations (DFA) with semilocal exchange-correlation functionals [17][18][19][20].However, these functionals miss the long-range vdW interactions for nonhomogeneous electron densities and cannot correctly describe multilayered h-BN and graphene.By adding a dispersion correction to DFA [21][22][23], this failure can be effectively cured.By applying the Tkatchenko-Scheffler (TS) vdW scheme to h-BN bilayer [23], Marom et al. proposed that vdW interactions play the role of fixing the interlayer distance, while the electrostatic forces determine the optimal mode and the interlayer sliding corrugation [9].In contrast, Constantinescu et al. have used the localized MP2 method to study bulk and bilayer h-BN, arguing that DFA augmented with a pairwise approximation for dispersion forces is unable to correctly describe the stacking of bilayer h-BN [10].Several other approaches based on DFA and manybody methods have been applied to describe dispersion forces in h-BN and graphite [4,[24][25][26][27][28][29][30].However, these studies did not address the interlayer potential of multilayered h-BN and graphene.
In this context, it is clear that a full and deep understanding of the sliding mechanism in multilayered materials has not yet been achieved.In this Letter, we systematically investigate this problem by employing DFA calculations including an accurate description of nonlocal many-body dispersion (MBD) interactions using the DFAþMBD method [31,32].We find that interlayer sliding constraints essentially determine the contribution of electrostatic interactions and dispersion forces to the sliding energy profile, leading to different sliding pathways.VdW interactions are found to contribute more to the interlayer sliding potential of polar bilayer h-BN than they do in nonpolar bilayer graphene.By comparing with the experimental findings, we identify the most likely sliding pathways which explain reduced friction for thicker vdW-bound layers and provide quantitative explanation for the observed anisotropy of the friction force.
All calculations in this work have been performed using the FHI-aims all-electron code with "tight" computational settings [33].We used a 6 × 6 × 1 supercell and a k-point grid of 4 × 4 × 1 throughout the work.To eliminate the interaction between periodic images, a vacuum width of 1000 Å is adopted in the direction perpendicular to the substrate surface.The geometries of multilayered h-BN and graphene were obtained using the Perdew-Burke-Ernzerhof (PBE) DFA functional [34] with the TS pairwise dispersion method [23].Subsequently, the energy was calculated at the PBEþMBD level of theory.Since the PBEþMBD binding energy for graphite obtained with the PBEþTS geometry deviates from the fully optimized PBEþMBD binding energy by 0.1 meV=atom, our approximation of computing PBE þ MBD energies using PBEþTS geometries is sufficiently reliable for all of the obtained conclusions.
As shown in Table I, the PBEþMBD method yields the interlayer distances, binding energies, and C 33 elastic constants of bulk h-BN and graphite in very good agreement with high-level calculations and experimental results.In particular, upon including many-body vdW interactions, there is a significant improvement in the obtained cohesive properties when compared with the pairwise TS approximation for these interactions.PBEþMBD calculations are especially good in reproducing the interlayer distances and elastic constants.The only apparent discrepancy with higher-level calculations is that the binding energy of bulk h-BN is slightly larger than that of graphite at the PBEþMBD level of theory.This is opposite to the trend obtained by random-phase approximation (RPA), which predicts graphite to be more stable than h-BN by 8 meV=atom [4,29].We note that RPA is prone to underestimating intermolecular vdW interactions [35], and this could explain part of the difference between PBEþMBD and RPA.Another reason for a possible overestimation of the binding energy in h-BN by PBEþMBD could be the ionic nature of bonding and the ensuing nonuniqueness of the electron density partitioning as discussed by Bučko et al. [36].We emphasize that here we study the relative interlayer energies involved in the sliding mechanism; thus, slight inaccuracies in the absolute binding energies would not influence the relevant findings of this paper.
Before proceeding to discuss sliding mechanisms, we first analyze the preferred stacking modes of bilayer and multilayer h-BN and graphene, and compare them with experimental observations.The interaction energies and interlayer distances are provided in Table SI of the Supplemental Material [45].In short, PBEþMBD predicts the AA 0 stacking to be the most stable configuration for bilayer h-BN, which is the stacking mode observed in experiments [46][47][48].As the number of layers increases from 2 to 5, multilayered h-BN structures consistently exhibit the same AA 0 stacking mode.Correspondingly, both the interaction energy and the interlayer distance asymptotically resemble the results of bulk h-BN: the binding energies increase from 25 to 43 meV=atom, while the interlayer distance decreases from 3.37 Å in bilayer to 3.33 Å in pentalayer.In the case of graphene, its bilayer exhibits AB stacking while trilayer prefers ABC stacking, which is the same order of stability found experimentally [49,50].As the number of graphene layers increases, multilayered graphenes again favor AB stacking modes (in tetralayer and pentalayer graphene).Meanwhile, the binding energies gradually increase from 23 meV=atom in bilayer to 39 meV=atom in pentalayer, while the interlayer distance decreases from 3.37 Å in bilayer to 3.35 Å in pentalayer, which also asymptotically resemble the results of graphite as found for multilayered h-BN.The larger change of the interlayer distance in h-BN compared with graphene is consistent with its smaller C 33 elastic constant, as shown in Table I.The magnitude of the change of the interlayer distance for different stacking modes will become crucial to understanding the different sliding mechanisms in h-BN and graphene, as we will explain below.
Starting from the most stable stacking modes, we now analyze the interlayer sliding pathways of h-BN and graphene.To repeat the commensurate stacking mode, there are two representative sliding directions: a short (S) pathway aligned along the metaposition of a sixmembered ring and a long (L) pathway along the paraposition of a six-membered ring [e.g.along N-N atoms and B-N atoms for h-BN in Fig. 1(a)].Our motivation to study these two pathways stems from recent experiments, which measured the friction of multilayered structures along commensurate directions [7].Moreover, intermediate pathways could be constructed as superpositions between S and L. The position of the top sliding layer is defined with respect to the bottom layer, which is always fixed.As shown in Fig. 1, by displacing the top layer along the S and L directions and relaxing all the other degrees of freedom (corresponding to S X and L Y pathways), we follow a minimum energy pathway for commensurate sliding.When the top layer slides along the S=L directions and only the Z direction of the top and intermediate layers is allowed to relax, we obtain the sliding energies of the partially constrained pathways S XY and L XY (see Fig. 1).Further fixing the Z direction of the top layer but fully relaxing the intermediate layers, the sliding profiles of the so-called S XYZ and L XYZ pathways are obtained.In short, the subscripts of the pathways denote the constrained sliding directions of the top layer.For graphene, the sliding pathways are defined in a similar manner.
For bilayer h-BN and graphene, the energetics for all of the discussed sliding pathways are shown in blue in Fig. 2. We denote the activation energy of interlayer sliding as E a , and separate the contribution from the dispersion energy (MBD) in red in Fig. 2. In our calculations, the L Y and L XY pathways turned out to be equivalent.Even upon manual displacement of the top layer in the X direction, the resulting L Y pathway for both graphene and h-BN followed a rigid displacement along the Y direction.This results in five different pathways shown in Fig. 2.
The sliding profiles in both materials can be discussed in terms of an interplay between electrostatics and vdW interactions (which include both vdW dispersive effects and Pauli repulsion).Since h-BN layer is polar [15], electrostatics are expected to play a larger role here than in nonpolar graphene.Even though carbon atoms in graphene do not carry a permanent charge or dipole moment, the charge density of graphene can be described by associating a permanent quadrupole moment to every carbon atom [51].Therefore, vdW interactions, rather than electrostatics, are expected to contribute more in graphene.
First, we note that there is a qualitative agreement between the different pathways in h-BN and graphene bilayers in Fig. 2. Clearly, vdW dispersive interactions make a crucial contribution to the S X pathways, although it has a different sign for h-BN and graphene.In the S XY and L XY pathways, the electrostatic interactions are found to play a more prominent role than in the S X pathway.Upon fully constraining the top layer (S XYZ and L XYZ ), the electrostatic interactions dominate the sliding energy profile.It is worth noticing that h-BN exhibits larger contribution from dispersion forces for all studied pathways than graphene does.Considering the polar nature of h-BN [15], one would expect electrostatics, rather than vdW dispersion, to have a more prominent role in the sliding mechanism.However, a deeper analysis of this issue reveals that the dispersion contribution depends sensitively on the variation of the interlayer distance upon sliding, which changes by a noticeable 0.07 Å in the S X pathway of h-BN compared to 0.04 Å in graphene.Finally, we note We now turn to analyze the thickness-dependent sliding mechanisms of multilayered h-BN and graphene, by computing the sliding energy profile of tri-, tetra-, and pentalayer h-BN and graphene (see Table SII and Fig. S2 in Ref. [45]).Overall, the complex interplay between different effects leads to a nonmonotonic dependence of E a on the thickness in the partially constrained S and L pathways.The relative contribution of dispersion interactions to sliding barriers either stays constant or increases for thicker multilayers of h-BN and graphene.As already found for bilayers, dispersion interactions have a more prominent role in the sliding profiles in multilayered h-BN than in graphene.In the case of fully constrained S XYZ and L XYZ pathways, the intermediate layers essentially serve the role of a lubricant.With more degrees of freedom involved in the sliding process, the sliding barriers are reduced with increasing thickness by up to 50% for h-BN and 45% for graphene.
We now discuss how our calculations can be used to rationalize existing experimental observations on friction of multilayered structures.Recent experiments on multilayers employed a sliding tip, and observed a monotonically decreasing friction for commensurate displacement in multilayered h-BN and graphene [7].We assume that the tip is attached to the substrate [1] and this fixes the vertical Z position of the top layer.In such an experimental situation, our computed S XYZ and L XYZ pathways provide a reasonable model.We approximate the friction of the S XYZ =L XYZ pathways by dividing the sliding barrier (Table SII) by the sliding distance (Table SIII), which is then normalized to the friction we obtain for the bilayer (see Fig. 3).The results demonstrate that the friction of each pathway (except the S XYZ pathway of h-BN) is monotonically decreased as the number of layers increases from bilayer to pentalayer for h-BN and graphene, providing a reasonable explanation for the experimental findings that are obtained by sliding a tip over 1-4 layers of h-BN and graphene [7].It is noteworthy that the experimental observations for h-BN [7] are reproduced quantitatively assuming sliding along the L XYZ pathway.In contrast, for multilayered graphene, the best agreement with experiment is achieved assuming sliding along the S XYZ pathway.In both cases, the monotonic decrease of friction is caused by a combination of decreasing sliding barrier (Table SII) and the change of sliding distance (Table SIII).The sliding distance associated with the sliding barrier can change by up to 33% depending on the thickness.
Since the friction of both the S XYZ and L XYZ pathways is essentially converged in terms of thickness for pentalayer graphene, this should thus be a good model for friction in bulk graphite.Therefore, the friction of graphite is calculated to be 5.05 pN=atom for the L 3 pathway and 1.92 pN=atom for the S XYZ pathway, to be compared to the experimental value of 2.11 pN=atom [6].These results suggest that the direction experiencing the peak friction force in experiments on graphite [6] is along the S XYZ pathway.We further calculated the friction anisotropy ratio of the L XYZ and S XYZ pathways for interlayer sliding of bilayer graphene.The value we obtain is 2.09, which is again in excellent agreement with the experimental value of 2.15 AE 0.08 that was measured for graphene with a sliding tip [52].
Finally, we calculated the friction of bilayer and bulk h-BN and compared it with graphene/graphite to motivate further experiments on this material.We obtained a friction of 5.47 pN=atom for the L XYZ pathway and 2.60 pN=atom for the S XYZ pathway in bulk h-BN; both values are similar to graphite.However, we observe much larger friction values in the the case of bilayer h-BN: 8.09 pN=atom for the L XYZ pathway and 5.28 pN=atom for the S XYZ pathway, leading to the friction anisotropy ratio equal to 1.53.Therefore, the magnitude of friction anisotropy in bilayer h-BN is smaller than in bilayer graphene, motivating further experimental studies to help understand this different behavior [53].
In conclusion, our PBEþMBD calculations provide fundamental insights into the mechanism of interlayer sliding of multilayered h-BN and graphene.We find that in both materials interlayer sliding constraints can be used to tune the contribution of electrostatic interactions and dispersion forces to the sliding energy profile, leading to rather different sliding pathways.As the number of layers increases, both the binding energy and the interlayer distance of multilayered h-BN and graphene asymptotically resemble the bulk results.Correspondingly, friction force is found to exhibit a monotonic decrease at a fixed top-bottom distance.These results enable us to provide explanations for several key experimental observations, including the reduced friction for thicker vdW-bound layers FIG. 3 (color online).Friction of the S XYZ and L XYZ pathways for multilayered h-BN and graphene along with experimental results [7].The friction is normalized to the value obtained for two layers.and the anisotropy of bulk and bilayer h-BN and graphene, and also make predictions about qualitatively different friction behavior of multilayered h-BN and graphene.This demonstrates the promise of recently developed first-principles methods for understanding nanomechanical behavior of low-dimensional materials.

FIG. 1 (
FIG. 1 (color online).Illustration of sliding pathways for bilayer h-BN.The top (bottom) panel refers to top view (side view) in each subfigure (a)-(f).The wavy lines indicate relaxation in the corresponding direction.