Interfacial Friction Anisotropy in Few-Layer Van der Waals Crystals

Friction anisotropy is one of the important friction behaviors for two-dimensional (2D) van der Waals (vdW) crystals. The effects of normal pressure and thickness on the interfacial friction anisotropy in few-layer graphene, h-BN, and MoSe2 under constant normal force mode have been extensively investigated by first-principle calculations. The increase of normal pressure and layer number enhances the interfacial friction anisotropy for graphene and h-BN but weakens that for MoSe2. Such significant deviations in the interfacial friction anisotropy of few-layer graphene, h-BN and MoSe2 can be mainly attributed to the opposite contributions of electron kinetic energies and electrostatic energies to the sliding energy barriers and different interlayer charge exchanges. Our results deepen the understanding of the influence of external loading and thickness on the friction properties of 2D vdW crystals.


Introduction
Layered van der Waals (vdW) materials such as graphite, boron nitride, and transition metal dichalcogenides (TMDs) which layers bind with weak interlayer vdW interactions, have been widely used as solid lubricants in engineering technology to reduce friction and wear. When size goes down to nanoscale, two-dimensional (2D) vdW crystals exhibit exceptional and excellent friction properties and have attracted numerous scientific interests . Lee et al. [1] used atomic force microscopy (AFM) technique to characterize the microscopic friction characteristics of monolayer and multilayer graphene, MoS 2 , h-BN, and NbSe 2 that were mechanically peeled off from weakly adherent SiO 2 substrates, and reveal that the friction force decreases with the increase of the number of layers. Zhang et al. studied the friction behavior of a graphene flake sliding on a supported graphene substrate and showed that the friction force increases exponentially with the decreasing stiffness [2][3][4][5]. Guo et al. [6] find that the interlayer friction force of graphene sheets at the commensurate state increases sharply with the decrease of the interlayer distance while the incommensurate contact maintains an ultra-low friction state. By using a pressurized bubble loading device, Wang et al. [7] found the existence of ultra-low friction characteristics in the incommensurate graphene/graphene interface at the microscopic scale. Negative friction coefficients are found at the interfaces of layered graphene-hexagonal boron nitride (h-BN) heterojunctions due to load-induced suppression of out-of-plane distortions [8], and the interfaces between graphene sheets and the tip of atomic force microscope tip [9,10]. The novel friction properties emerging in 2D vdW crystals make them possessing a promising future in the application of functional devices and nanoelectromechanical systems.
Due to surface morphology, crystal lattice, or structural deformation, the friction along various directions is usually different. Friction anisotropy is one of the important friction behaviors and has been extensively studied for solid interfaces and surfaces by theoretical and experimental methods. For 2D vdW crystals, the friction anisotropy also becomes obvious and remarkable [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46]. Nevertheless, most previous first-principles studies on the interlayer friction for two 2D vdW crystals were conducted with constant interlayer distance mode rather than constant normal force mode. Empirical potential based molecular dynamics simulations are hard to accurately describe interlayer coulomb interaction and charge exchange under compression. It is necessary to use more accurate methods to investigate the interfacial friction behaviors of 2D vdW crystals in the presence of ideal constant normal force. On the other hand, mechanical loading is an effective way to modify the physical and chemical properties of low-dimensional materials. Monolayer and few-layer vdW crystals usually possess different friction behaviors. However, the effects of normal pressure and thickness on the interlayer friction anisotropy of few-layer 2D vdW crystals have been seldom studied and remains elusive.
In this work, the interfacial friction properties of few-layer graphene, h-BN, and MoSe 2 under constant normal force mode have been extensively studied by using first-principles calculations. It is found that the deviations in the energy barriers between the interlayer maximum energy sliding paths and minimum energy sliding paths increase for graphene and h-BN but decrease for MoSe 2 with the applied normal pressure increases. The friction anisotropy in graphene and h-BN is therefore enhanced but reduced in MoSe 2 by the normal pressure. Moreover, the increase of layer number increases the friction anisotropy in graphene and h-BN and decreases that in MoSe 2 . The significant difference in the influence of normal pressure and layer number on the friction anisotropy among these three kinds of 2D vdW crystals can be attributed to different interlayer charge exchanges and the opposite changes in electron kinetic energies and electrostatic energies. Our results deepen the understanding of the effects of normal pressure and thickness on the friction behaviors of 2D vdW crystals.

Methods
In our model, we established 2 to 8-layer graphene, h-BN and 2H MoSe 2 with interlayer AB stacking in the rhombus unit cells, where graphene, h-BN and MoSe 2 monolayers consist of two C atoms, 1 B, and 1 N atom, and 1 Mo and 2 Se atoms, respectively. Figure 1 shows the atom structures of six-layer graphene, h-BN and MoSe 2 . In the unit cells, a vacuum region larger than 15 Å is in the direction perpendicular to the atomic plane. All computations were performed within the framework of density-functional theory (DFT) as implemented in the FHI-aims code with "tight" computational settings [47] in which the Perdew-Burke-Ernzerhof (PBE) [48] was employed. The influence of vdW interactions was considered by using a many body dispersion (MBD) vdW model [49,50]. A k point grid of 15 × 15 × 1 was used throughout the work. First the whole systems were relaxed by PBE + MBD until the force on each atom was less than 0.01 eV/Å. The optimized lattice constants of graphene, h-BN and MoSe 2 unit cells are 2.467, 2.507, and 3.28 Å, respectively.
To simulate interlayer sliding, the whole a 1 -a 2 planes of the unit cells were equally divided into 81 positions where the nearest translational positions were separated by 0.274, 0.279, and 0.364 Å for graphene, h-BN, and MoSe 2 , respectively. The top parts of few-layer graphene, h-BN and MoSe 2 including 1, 2, 3, or 4 layers were transversely moved as a whole part with respect to the bottom parts, as shown in Figure 1d, and shifted relatively to different divided positions in the a 1 -a 2 planes. Different interlayer stacking were achieved when the top parts shifted to different divided positions. The interlayer distances d of the few-layer 2D crystals were modified by changing the z-direction coordinates of the atoms in the top layers. For graphene and h-BN, 1 C and 1 B atom in the top layers and 1 C and 1 B atom in the bottom layers were fully fixed at each shifted position. For MoSe 2 , only the Se atoms at the top and bottom surfaces were fully fixed. Then those systems were relaxed again by PBE + MBD. After relaxation, the normal forces F n at different interlayer distances and shifted positions were calculated by summing the z-direction forces of all atoms in the uppermost layers, and the total energies E total of the few-layer 2D crystals at different shifted positions and interlayer distances were calculated. Next, the total energies at a given shifted position were fitted with the corresponding normal pressure P n (P n = F n /A, A is the unit cell area) by a 3-order polynomial, as shown in Figure 2.
According to the fitting curves, the total energies at a given normal pressure was obtained. Through this way, the interlayer sliding simulation under constant normal force mode was realized. Then the potential energy surfaces (PES) for the structures with all the shifted positions under constant normal force mode were constructed by ∆E = E total − E min , where E min is the lowest total energy. To simulate interlayer sliding, the whole a1-a2 planes of the unit cells were equally divided into 81 positions where the nearest translational positions were separated by 0.274, 0.279, and 0.364 Å for graphene, h-BN, and MoSe2, respectively. The top parts of few-layer graphene, h-BN and MoSe2 including 1, 2, 3, or 4 layers were transversely moved as a whole part with respect to the bottom parts, as shown in Figure 1d, and shifted relatively to different divided positions in the a1-a2 planes. Different interlayer stacking were achieved when the top parts shifted to different divided positions. The interlayer distances d of the few-layer 2D crystals were modified by changing the z-direction coordinates of the atoms in the top layers. For graphene and h-BN, 1 C and 1 B atom in the top layers and 1 C and 1 B atom in the bottom layers were fully fixed at each shifted position. For MoSe2, only the Se atoms at the top and bottom surfaces were fully fixed. Then those systems were relaxed again by PBE + MBD. After relaxation, the normal forces Fn at different interlayer distances and shifted positions were calculated by summing the z-direction forces of all atoms in the uppermost layers, and the total energies of the fewlayer 2D crystals at different shifted positions and interlayer distances were calculated. Next, the total energies at a given shifted position were fitted with the corresponding normal pressure Pn (Pn = Fn/A, A is the unit cell area) by a 3-order polynomial, as shown in Figure 2. According to the fitting curves, the total energies at a given normal pressure was obtained. Through this way, the interlayer sliding simulation under constant normal force mode was realized. Then the potential energy surfaces (PES) for the structures with all the shifted positions under constant normal force mode were constructed by ∆ = − , where is the lowest total energy.

Results and Discussion
As shown in Figures 3 and 4, for all of these 2D crystals, the AA stacking has the highest energy while the AB stacking has the lowest energy at a given normal pressure. There are two typical sliding paths P1 and P2, see Figures 3 and 4: P1 is from AB to AB via an apex point AH in which path the energy barrier ∆ (∆ = − ) is the lowest compared with the other paths, another P2 is AB to AA to AB in which path the energy barrier ∆ (∆ = − ) is the highest. The apex point AH with the highest energy along the sliding path P1 is defined as AH stacking. Figure 5 shows the

Results and Discussion
As shown in Figures 3 and 4, for all of these 2D crystals, the AA stacking has the highest energy while the AB stacking has the lowest energy at a given normal pressure. There are two typical sliding paths P 1 and P 2 , see Figures 3 and 4: P 1 is from AB to AB is the lowest compared with the other paths, another P 2 is AB to AA to AB in which path the energy barrier ∆E P 2 max (∆E P 2 max = E AA total − E AB total ) is the highest. The apex point AH with the highest energy along the sliding path P 1 is defined as AH stacking. Figure 5 shows the variations of ∆E P 1 max and ∆E P 2 max with normal pressure for 2L-and 6L-layer graphene, h-BN and MoSe 2 . For bilayers, both ∆E P 1 max and ∆E P 2 max increase with the pressure increases. When the layer number becomes 6, ∆E P 1 max of graphene first increases and then decreases with normal pressure. In contrast, ∆E P 2 max of MoSe 2 decreases with the normal pressure increases. The interlayer sliding energy barriers are affected by the change in layer number. On the other hand, the deviation ∆E P 2 max − ∆E P 1 max between the energy barriers ∆E P 1 max and ∆E P 2 max is directly related to the interfacial friction anisotropy of the considered few-layer 2D vdW crystals. The higher ∆E P 2 max − ∆E P 1 max means the stronger friction anisotropy while the lower means the weaker friction anisotropy. It can be seen from Figure 6 that the deviations ∆E P 2 max − ∆E P 1 max of graphene and h-BN increase with the pressure increases, indicating that normal loading enhances the interfacial friction anisotropy. Moreover, ∆E P 2 max − ∆E P 1 max of graphene and h-BN is further increased with the increase of layer number. When the thicknesses of few-layer graphene and h-BN becomes thicker, the interfacial friction anisotropy increases as well. On the contrary, the deviation ∆E P 2 max − ∆E P 1 max of MoSe 2 decreases with the pressure increases so that normal loading weakens the interfacial friction anisotropy. Meanwhile, the interfacial friction anisotropy of few-layer MoSe 2 decreases with the increase of layer number. Obviously, the normal pressure and thickness increasing impose completely different influence on the interfacial friction anisotropy of few-layer graphene, h-BN, and MoSe 2 .         In our DFT calculations, the total energy consists of kinetic, electrostatic, exchange-correlation, and vdW energies, namely = + + + . To better In our DFT calculations, the total energy E total consists of kinetic, electrostatic, exchangecorrelation, and vdW energies, namely E total = E K + E H + E XC + E vdW . To better understand the different friction anisotropy in these few-layer 2D crystals, the energy differences in exchange-correlation energy ∆E XC , electron kinetic energy ∆E K , electrostatic energy ∆E H and vdW energy ∆E vdW between AH or AA and AB stacking have been calculated by ∆E where X = K, H, XC, and vdW representing electron kinetic, electrostatic, exchange-correlation, and vdW energies, respectively. Then Figure 7 shows the variations of ∆E P 2 X − ∆E P 1 X with normal pressure for 2L-and 6L-layer graphene, h-BN and MoSe 2 . For graphene, ∆E P 2 K − ∆E P 1 K increases but ∆E P 2 H − ∆E P 1 H decreases as the pressure increases, while ∆E P 2 XC − ∆E P 1 XC and ∆E P 2 vdW − ∆E P 1 vdW slightly changes with normal pressure. Similarly, H decreases with the increase of normal pressure. In contrast, ∆E P 2 K − ∆E P 1 K of MoSe 2 decreases but ∆E P 2 H − ∆E P 1 H increases as the pressure increases. Therefore, the opposite contributions and cancelling effects of the electron kinetic energies and electrostatic energies lead to such quite different friction anisotropy between graphene, h-BN and MoSe 2 . It is shown from Figure 7 that the changes in pressure and layer number slightly influence ∆E P 2 vdW − ∆E P 1 vdW . The increase in the layer number of the considered 2D crystals modifies the values of ∆E P 2 X − ∆E P 1 X but will not change the total variation trends.
Therefore, the opposite contributions and cancelling effects of the electron kinetic energies and electrostatic energies lead to such quite different friction anisotropy between graphene, h-BN and MoSe2. It is shown from Figure 7 that the changes in pressure and layer number slightly influence ∆ − ∆ . The increase in the layer number of the considered 2D crystals modifies the values of ∆ − ∆ but will not change the total variation trends. In order to further elucidate the mechanism of interfacial friction anisotropy, the interlayer charge density differences ∆ for 6L-layer graphene, h-BN, and MoSe 2 at AB, AA, and AH stacking have been calculated under the normal pressures of 0 and 2.5 GPa by ∆ = − − , where is the total charge density, and are the charge densities of top three layers and bottom three layers, respectively. As shown in Figure 8, for graphene and h-BN the interlayer charge depletions increase with the pressure increases, and more charges move to the interfacial C atoms and N atoms. For MoSe2, the interlayer charge accumulations and depletions are relatively weak and charge exchange mainly occurs around the interfacial Se atoms. As a result, the effect of pressure on the interfacial friction anisotropy in few-layer MoSe2 is different from that of few-layer graphene and h-BN. Moreover, we also considered a 6L-layer graphene In order to further elucidate the mechanism of interfacial friction anisotropy, the interlayer charge density differences ∆ρ for 6L-layer graphene, h-BN, and MoSe 2 at AB, AA, and AH stacking have been calculated under the normal pressures of 0 and 2.5 GPa by ∆ρ = ρ 6L−layer − ρ top−part − ρ bottom−part , where ρ 6L−layer is the total charge density, ρ top−part and ρ bottom−part are the charge densities of top three layers and bottom three layers, respectively. As shown in Figure 8, for graphene and h-BN the interlayer charge depletions increase with the pressure increases, and more charges move to the interfacial C atoms and N atoms. For MoSe 2 , the interlayer charge accumulations and depletions are relatively weak and charge exchange mainly occurs around the interfacial Se atoms. As a result, the effect of pressure on the interfacial friction anisotropy in few-layer MoSe 2 is different from that of few-layer graphene and h-BN. Moreover, we also considered a 6L-layer graphene with initial ABC stacking, h-BN with initial AA' stacking, and MoSe 2 with initial AA' stacking. The energy barriers between the interlayer maximum energy sliding paths and minimum energy sliding paths were calculated by the same method and procedure. ∆E P 2 max − ∆E P 1 max of 6L-layer graphene with ABC stacking and h-BN with AA' stacking increases with the normal pressure increases, which is consistent with that of 6L-layer graphene and h-BN with initial AB stacking. On the contrary, ∆E P 2 max − ∆E P 1 max of MoSe 2 with AA' stacking decreases with the normal pressure increases, which is also consistent with that of 6L-layer MoSe 2 with initial AB stacking. Therefore, the initial stacking mode will not change the qualitative variations of friction anisotropy with normal pressure.

∆
− ∆ of 6L-layer graphene with ABC stacking and h-BN with AA' stacking increases with the normal pressure increases, which is consistent with that of 6L-layer graphene and h-BN with initial AB stacking. On the contrary, ∆ − ∆ of MoSe2 with AA' stacking decreases with the normal pressure increases, which is also consistent with that of 6L-layer MoSe2 with initial AB stacking. Therefore, the initial stacking mode will not change the qualitative variations of friction anisotropy with normal pressure.

Conclusions
In summary, we show by comprehensive DFT calculations that under constant normal force mode the interfacial friction anisotropy in few-layer 2D vdW crystals is significantly modified by normal pressure and layer number. For graphene and h-BN, the friction anisotropy increases with increasing the normal pressure and layer number. On the contrary, the friction anisotropy for MoSe2 decreases with increasing the normal pressure and layer number. The opposite contributions of electron kinetic energies and electrostatic energies to the sliding energy barriers and different interlayer charge exchanges in fewlayer graphene, h-BN, and MoSe2 lead to such remarkable deviation in interfacial friction anisotropy. These results deepen our understanding of friction behaviors of 2D vdW crystals and provide some new insights into the application of 2D vdW crystals as solid lubricants.

Conclusions
In summary, we show by comprehensive DFT calculations that under constant normal force mode the interfacial friction anisotropy in few-layer 2D vdW crystals is significantly modified by normal pressure and layer number. For graphene and h-BN, the friction anisotropy increases with increasing the normal pressure and layer number. On the contrary, the friction anisotropy for MoSe 2 decreases with increasing the normal pressure and layer number. The opposite contributions of electron kinetic energies and electrostatic energies to the sliding energy barriers and different interlayer charge exchanges in few-layer graphene, h-BN, and MoSe 2 lead to such remarkable deviation in interfacial friction anisotropy. These results deepen our understanding of friction behaviors of 2D vdW crystals and provide some new insights into the application of 2D vdW crystals as solid lubricants.