We investigated the effect of the interlayer twisting angle on the frictional performance of four new heterostructures. In Figure 3, the interlayer friction of the heterostructure is clearly shown to vary with the rotation angle at a constant normal force of 0.02 nN/atom and a temperature of 1 K. The number of atoms of the upper sheet is denoted as N. For the graphene/MoS₂ and graphene/NbSe₂ systems, the friction coefficient shows a periodic variation with the rotation angle. This is consistent with the previous results of Song et al. (Song et al., 2018). The curve takes 60° as the period and is determined by the periodic hexagonal lattice structure. When the rotation angle is 0° and 60°, the friction coefficient reaches the maximum, and the heterostructure is in the aligned contact state. The friction coefficient is at its lowest at a rotation angle of 30°. Furthermore, we considered the effect of temperature on the interlayer friction performance for various rotation angles, and the results are shown in Supplementary Figure S2. The increase in temperature leads to an increase in the thermal vibration of the atoms, which results in a decrease in the interlayer frictional resistance, and then the dependence of the friction coefficient on the rotation angle is weakened. As shown in Figures 3B–F, we compared the friction coefficients obtained by sliding sheets of different sizes on the substrate, and we can perceive that the smaller the size of the sheet is, the more difficult it is to achieve superlubricity (see Effect of Torsion Angle on Friction for a detailed discussion). The larger the size of the sheet is, the larger the range of angles to achieve superlubricity. Comparing these two different systems, the friction coefficient of graphene/NbSe₂ is lower than that of graphene/MoS₂ under the same conditions. This is due to the greater lattice mismatch of the graphene/NbSe₂ heterostructure. In comparison with the graphene/graphene homostructure system (Ru et al., 2020), it can be seen that the friction coefficient of the homogeneous structure system is higher than that of the heterostructure in the aligned contact state. In addition, the homostructure system is more sensitive to the rotation angle, while the friction coefficient of the heterostructure decreases slowly after the change in the rotation angle. This is consistent with previous studies (Dienwiebel et al., 2004; Filippov et al., 2008; Wang et al., 2019a; Ru et al., 2020).

By combining the analysis with the potential energy surface (PES), we further explored the relationship between the rotation angle and the friction force. At the microscopic scale, potential energy fluctuations at the interface are the origin of kinetic friction (Guo et al., 2011). Therefore, we mapped the potential energy surface = profile of the heterostructure at different rotation angles, as shown in Figure 4. The interlayer interaction energy between the top graphene sheet and the substrate was calculated using the LJ potential as a function of the position coordinates of the center mass of the sheet (X- and Y-axes chosen along the armchair and zig-zag directions, respectively) and the relative twist angle θ. From Figure 4, we can clearly see that for graphene/MoS₂ and graphene/NbSe₂ heterostructures, the PES fluctuations at the rotation angle θ = 0° are much larger than those at other angles because the lattice of the graphene sheet and the substrate form an aligned contact. The corresponding PES fluctuations at other rotation angles are smoother. The smoother the PES fluctuations are, the less energy consumed in the sliding process, and the system structure is more likely to enable structural superlubricity.

In the same way, we investigated the interlayer sliding behavior of α-Te/MoS₂ and α-Te/NbSe₂ heterostructures under the same simulation conditions, as shown in Figure 5. From the calculation results in Figures 5A–C, it can be found that for the α-Te/MoS₂ heterostructure, the interlayer friction shows an obvious periodicity with the change in rotation angles. However, for the α-Te/NbSe₂ heterostructure, the periodicity of the interlayer friction is not significant when the slider size is small. After increasing the size, the periodicity gradually appears. Since the lattice periodicity is broken at the boundaries, this leads to a reduction in the binding of the boundary atoms. As a result, the boundary atoms are more active than the central atoms. A smaller slider size means more pronounced edge effects, which we discuss in detail in Effect of Torsion Angle on Friction. To further elucidate the reason for the variation in friction with angle for the α-Te/MoS₂ and α-Te/NbSe₂ systems, we plotted the sliding potential energy surface configurations of these two systems, as shown in Figure 6. We can clearly see that the degree of undulation of the potential energy surface corresponds to the magnitude of the friction coefficient. A small undulation of potential energy means that the system has a small energy loss during sliding, i.e., a small sliding resistance. For both heterostructures, α-Te/MoS₂ and α-Te/NbSe₂, the potential energy surface is flatter when the interlayer twist angle is 15° than when the twist angles are 0° and 30°. In Figure 5, these two heterostructures do have lower friction coefficients at a rotation angle of 15°, and both are consistent with the results from the potential energy surface calculation. Meanwhile, we found that for α-Te/MoS₂ and α-Te/NbSe₂ heterostructures, the rotation angle regulates the interlayer friction of the system, but superlubricity cannot occur (the coefficient of friction is greater than one thousandth of an order of magnitude). The α-Te/MoS₂ and α-Te/NbSe₂ heterostructures require more potential barriers to be overcome in the slip path than the graphene/MoS₂ and graphene/NbSe₂ heterostructures due to the dominance of the frictional mechanism of interfacial sliding resistance in tellurene heterostructures. However, for graphene heterojunctions, friction is dominated by boundary effects, as detailed in the analysis in Effect of Torsion Angle on Friction.

Friction is independent of the macroscopic contact area because the actual area of contact on the atomic scale is a fraction of the total surface area. However, at the atomic scale, it is crucial to understand the relationship between contact size and friction to control friction. Mandelli et al. found that for all contact sizes considered, the commensurate homogeneous interface exhibits highly dissipative viscous slip motion, leading to size-independent static and dynamic friction (Mandelli et al., 2017). In contrast, incommensurate interfaces have significant size effects. In this section, we investigate the relationship between contact size and friction for heterogeneous structures that are under aligned contact.

Herein, we discuss the variation of the frictional force with size for three different shapes of sliders using the parameter, i.e., the perimeter-to-area ratio (P/A), which has been used in reference (Belghachi and Khelifi, 2006). The higher the value of P/A is, the larger the proportion of atoms at the edge of the slider. As shown in Figure 7, for the four heterostructures, the friction of the heterostructures tends to decrease with increasing P/A, which is consistent with the conclusion of previous calculations (Wang et al., 2019c). In addition, Figure 7 shows that the friction decreases at a faster rate when P/A < 0.2, with an approximately linear decrease. However, when P/A>0.2, the rate of decline gradually slows down. As mentioned before, in our MD simulation, normal loading is achieved by loading 0.02 nN to each atom of the outermost layer of the slider. A larger slider implies a larger normal load. Due to the different structures of graphene and α-Te, the same contact area corresponds to different amounts of normal load. Therefore, we redrew the relationship curve between the friction coefficient and P/A, as shown in Supplementary Figure S3 (Supplementary Information S1). The coefficient of friction vs. specific perimeter (P/A) for the heterostructures by graphene and α-Te exhibited different variation trends. For the graphene/MoS₂ and graphene/NbSe₂ heterostructures, the friction coefficient increases linearly with increasing P/A, which means that the smaller the slider size is, the higher the frictional resistance. A smaller slider size also means a higher percentage of atoms at the edge. However, for the α-Te/MoS₂ and α-Te/NbSe₂ heterostructures, the friction coefficient decreases linearly with increasing P/A, implying that the larger the slider size is, the larger the frictional resistance. A larger slider size also means a larger contact area at the interface. This indicates that they have two completely different friction mechanisms.

To further investigate the two different interlayer friction mechanisms, we plotted the variation in the friction force (average to each atom of the slider) with a specific perimeter (P/A). The ratio of edge atoms to internal atoms decreases as the slider size increases (as shown in Supplementary Figure S4, Supplementary Material S1), and the contribution of the edge atoms to the system friction decreases. In the graphene/MoS₂ and graphene/NbSe₂ heterostructures, as shown in Figures 8A,C, the single-atom friction shows a linear increase with increasing P/A of the graphene slider. This phenomenon is related to the fact that the friction mainly comes from the contact interface edge pinning. The same conclusion was shown by Qu et al. (Qu et al., 2020), where they separated the inner area of the slider contact interface from the edge area, calculated the friction force of both separately, and confirmed experimentally that the friction force mainly comes from the edge. In addition, edge atoms can exhibit large distortions due to their low confinement, which can lead to large friction (Guo et al., 2011; Mandelli et al., 2017).

However, in the α-Te/MoS₂ and α-Te/NbSe₂ heterostructures, as shown in Figures 8B,D, the abnormal phenomenon was observed, i.e., a trend of decreasing friction (averaged to each atom) with increasing P/A. The increase in P/A means that the percentage of boundary atoms increases. The friction shows a decreasing trend, indicating that for this heterostructure, the interlayer friction is not dominated by edge pinning but due to the in-plane interface sliding resistance mainly contributing to the friction. Liao et al. defined the shear strength S and edge-pinning strength E of a finite size slider as follows (Liao et al., 2021): S = F r A (1) E = F r P (2) where F r is the frictional force, A is the area of the domain, and P is the perimeter of the slider. When the friction is mainly influenced by the sliding resistance at the interface, the shear strength S is constant with respect to the slider area. When the friction force is mainly from edge pinning, E is constant with respect to the slider perimeter. As shown in Figure 9, we plotted the shear strength S versus domain area for the four systems. In addition, the inset shows the trend of the corresponding edge pinning strength E with the domain perimeter.

In the graphene/MoS₂ and graphene/NbSe₂ heterostructures, as shown in Figures 9A,C, the edge pinning strength E is constant relative to the perimeter of the slider. We believe that the friction force mainly comes from the edge. In the α-Te/MoS₂ and α-Te/NbSe₂ heterostructures, as shown in Figures 9B,D, the shear strength S is constant relative to the area of the slider. We believe that in these two systems, the friction force is dominated by the sliding resistance in the interface.

According to the Prandtl-Tomlinson model (Tomlinson, 1929), when the upper moving atom falls into the potential well formed by the static atoms, if the upper atom moves, it is necessary to overcome the potential barrier of the potential well, resulting in increased resistance. However, if the upper atoms do not fall into the potential well, the resistance is very small when moving. This atomic-scale resistance is friction in the macro view, which is also the origin of friction. For the systems mentioned above, that is, a two-dimensional material slider sliding on the two-dimensional material plane, we can regard the two-dimensional material matrix as a potential field. From the atomic level, the contribution of each atom of the slider to the friction force depends on the extent to which the atom falls into the matrix potential well. The deeper it falls into the potential well, the greater its contribution to the friction force, and the shallower it falls into the potential well, the smaller the contribution to friction. From the slider, each atom on the contact surface of the same slider is the same, but the position of the atom is different. The range of atomic motion in the sliding process is different, which leads to the fact that if the range of atomic motion is large, the atom easily falls into the potential well of the matrix. The range of atomic motion can be expressed by the mean-square displacement (MSD) of the atom (Kubečka et al., 2016). By using the built-in commands in the Lammps software, we can easily calculate the mean square displacement of the system atoms with the corresponding defining equation shown below: M S D = 〈 1 N ∑ i = 1 N [ r i ( t 0 + t ) − r i ( t 0 ) ] 2 〉 (3) where N is the number of particles, t represents time, r (t ₀ + t)-r (t ₀) is the vector distance passed by a given particle in a period of time, and <...> is the mean of the system synthesis after equilibrium.

In other words, the contribution of the slider atom to the friction is related to the MSD of the atom. If the MSD of the atom is large, the contribution to the friction is large; otherwise, the contribution to the friction is small.

As shown in Figure 10, we distinguish the slider into two parts: the edge and the inner area. We further calculated the average MSD of the atoms in the two areas separately and averaged over the whole MD trajectory. The final calculation results are presented in Table 3.

For the graphene/MoS₂ heterostructure, the average atomic displacement in the edge area is significantly larger than the average displacement in the central area, which indicates that the edge atoms have a relatively larger range of motion, and thus, the edge atoms are more likely to be trapped in the energy minimum of the substrate, leading to the edge pinning effect. Thus, the edge atoms cause more friction than the central atoms, and for the total friction of the slider, the contribution of the edge atoms dominates.

For the α-Te/MoS₂ heterostructure, the displacements of the atoms in the boundary region remain larger than the displacements of the atoms in the inner region, which can be explained by the fact that edge effects are always available for a finite size slider. Differently from the graphene/MoS₂ system, the atoms in the internal region of the α-Te slider also have a large atomic shift, which indicates that the internal atoms are equally prone to fall into the potential wells formed by the substrate, and the frictional force increases with the increase of the interfacial contact area; therefore, for the α-Te/MoS₂ heterostructure, the interfacial sliding resistance is the main source of frictional force.

In general, the boundary pinning effect exists in all studied sliders. For graphene/MoS₂ and graphene/NbSe₂ heterostructure systems, this boundary pinning is the main source of friction. For α-Te/MoS₂ and α-Te/NbSe₂ heterostructures, the in-plane atomic pinning effect is more obvious; that is, boundary pinning is not the main source of friction, and in-plane interface friction plays a leading role.

Monolayer graphene has been reported to have higher friction than multilayer graphene and graphite (Lee et al., 2009); however, the underlying mechanism remains to be discussed. Earlier studies attributed the increased friction of thinner two-dimensional samples to puckering and local pinning (Lee et al., 2010; Li et al., 2010; Ye et al., 2012; Li et al., 2016). In our simulations, graphene sheets with different layer thicknesses were chosen to slide on the substrate with the number of graphene layers (M) varying from 1 to 10. Meanwhile, we also considered three shapes: triangular, rectangular, and hexagonal sheets. For different shapes, the number of atoms per layer (N) is 1980, 1860 and 1,260. The initial spacing between adjacent layers of graphene is set to 0.34 nm. Figure 11 shows the relationship between the friction force and the number of layers of graphene with different shapes. In general, a smaller number of layers corresponds to a larger frictional force. The friction coefficient shows a decreasing trend as the number of graphene layers increases, and the decreasing trend slows down when the number of layers M > 4.

According to the theory given in reference (Liao et al., 2021), the interlayer friction comes from two parts: one is caused by adhesion, and the other is caused by increasing load. When the load is constant, the magnitude of friction changes with the change in adhesion. When we simulate the relationship between the number of layers and friction, the load is also constant. According to the literature (Pourzand et al., 2013), the interlayer adhesion decreases with increasing layer thickness (number of layers). That is, as the number of layers of the slider increases, the adhesion between the slider and the substrate decreases, resulting in a decrease in friction, as shown in Figure 11.

For a further comparative study of the effect of temperature on the interlayer friction of different systems, simulations were carried out for four heterostructures at temperatures of 1 K, 50 K, 100 K, 150 K, 200 K, 250 and 300 K. Throughout our simulations, the normal load and slip velocity were set to 0.02 nN/atom and 1 Å/ps, respectively. The sliding direction of the slider is set to slide along the X-direction. The simulation results are shown in Figure 12. We find that the variation trend of the friction force of the four systems with temperature is almost the same; in other words, the friction first decreases with increasing temperature and then remains essentially constant, which is consistent with previous reports (Zhao et al., 2007; Smolyanitsky, 2015; Gongyang et al., 2020).

The influence of temperature on the interlayer friction of the heterostructure can be explained as follows: as mentioned above, the interlayer friction of the two-dimensional van der Waals heterostructure can be regarded as being composed of two parts (Liao et al., 2021): one is from the adhesion between layers (interlayer interaction), and the other is the friction caused by loading (related to the structure of the two-dimensional material). To study the effect of temperature on friction, we have taken the same load in the simulation. In this way, we only need to study the relationship between interlayer adhesion and temperature. Deng et al. show theoretical studies that (Deng and Berry, 2016), with the increase of temperature, atomic thermal vibration intensifies, and nanoscale ripples will inevitably appear in two-dimensional materials, resulting in the reduction of effective contact area and the reduction of adhesion force. The conclusion that the interlayer adhesion force of two-dimensional materials decreases with increasing temperature has been confirmed by experiments (Polfus et al., 2021).

We calculated the coordinates of the upper S atoms of the base material MoS₂ in the Z-axis direction, which is intended to indicate the degree of fluctuation of the sliding contact interface. (see Supplementary Figure S5). It can be clearly seen that the contact interface between the slider and substrate fluctuates greatly, while it remains flat in other noncontact regions. In addition, as shown in Supplementary Figure S5B, the kinetic energy distribution of the contact interface is higher than that of other regions. This result shows that it is the nanoscale ripple in the contact interface that directly leads to the reduction of the real contact area. In other words, the interlayer friction of two-dimensional materials decreases with increasing temperature. It is worth mentioning that the interlayer friction changes rapidly with temperature and soon reaches a constant value, which is similar for the four systems, indicating that the nanoscale ripple no longer increases after reaching a certain temperature, the adhesion force reaches a constant value, and the interlayer friction also reaches a constant value.

Interlayer friction anisotropy in two-dimensional materials has been studied for a long time, and there is evidence that interlayer friction anisotropy in two-dimensional materials is caused by different lattice orientations. Here, we compare the variation in interlayer friction in four heterostructures when the slider moves along different sliding directions. We selected the upper slider to slide along the 0°, 15°, 30°, 45°, 60°, 75° and 90° directions, where 0° exactly corresponds to the sliding direction along the X-direction. The normal load was set to a constant 0.02 nN/atom, and the sliding speed was set at 1 Å/ps.

To further explicitly describe the effect of slip direction on interlayer friction, the system temperature was set to 1 K for the simulation. As shown in Figure 13, we can find that the variation in the friction force with the slip direction is relatively small for the heterostructures; that is, the heterostructures have lower friction anisotropy than the homostructures, which is consistent with the conclusions in previous research works (Ru et al., 2020). Due to the natural lattice mismatch of the heterogeneous structure, even when the heterogeneous structure is in aligned contact, the upper slider can slide smoothly regardless of the sliding direction, thus exhibiting a small frictional anisotropy. The friction anisotropy still shows variability for different heterogeneous structural systems. As shown in Figure 13B, both α-Te/MoS₂ and α-Te/NbSe₂ heterostructures have greater anisotropy than the graphene/MoS₂ and graphene/NbSe₂ systems and exhibit a significant 60° periodicity. Compared with the latter, the potential energy surface of the former fluctuates greatly during the sliding process, showing more obvious anisotropy.

We further studied the velocity-dependent interlayer friction. The normal load was kept at 0.02 nN/atom, the sliding direction was fixed along the X-direction, and the system temperature was controlled at 1 K. The detailed results are shown in Figure 14.

As seen from Figure 14A, for the graphene/MoS₂ and graphene/NbSe₂ heterostructures, the friction force remains essentially constant at low sliding speeds, while it increases significantly with increasing speed. The simulation results are the same as those reported in the previous literature (Zwörner et al., 1998; Sheehan and Lieber, 2017). Additionally, Zwörner demonstrated experimentally that the linear increase in friction with sliding speed at large sliding speeds is due to the strong damping effect (Zwörner et al., 1998). However, we obtained completely different results in the α-Te/MoS₂ and α-Te/NbSe₂ heterostructures. As the sliding velocity increases from 0.1 to 0.5 Å/ps, the friction force increases approximately linearly. However, when the sliding speed exceeds 0.5 Å/ps, there is a decreasing tendency for the friction, as shown in Figure 14B. We find that for α-Te/MoS₂ and α-Te/NbSe₂ heterostructures, as the sliding speed increases, the temperature in the contact part will also increase, i.e., the local temperature increases significantly with increasing sliding velocity. The increase in system temperature comes from the increased atomic thermal vibration during high-speed sliding, as mentioned in Effect of temperature on friction, which leads to an increase in friction in the α-Te/MoS₂ and α-Te/NbSe₂ systems. It has been reported that surface-to-surface contact contributes to thermal activation. Thermal activation, in turn, is closely related to the sliding speed (Smolyanitsky, 2015). When the sliding velocity is relatively fast, thermal activation may occur, resulting in a smaller slip resistance, thus showing a decrease in friction, as shown in Figure 14B.

We studied the normal load versus interlayer friction for the four heterostructures, while the chosen normal loads were 0.02, 0.1, 0.2 and 0.5 nN/atom. The simulation results are shown in Figure 15.

It is found that the normal load is approximately linearly related to the friction in the heterostructures studied, agreeing with the results of previous studies (Miura and Kamiya, 2002; Lee et al., 2009; Filleter and Bennewitz, 2010; van Wijk et al., 2013; Levita et al., 2014; Baykara et al., 2018; Wang et al., 2019a). Large normal loads tend to cause deformation of the plane. In addition, for a single layer slider, an increase in normal load means a decrease in layer spacing. The effective contact area of the system interface increases, leading to the stick-slip phenomenon in the system. As a result, the interaction between atoms becomes stronger, leading to a linear increase in friction.