Stiffness and Atomic-Scale Friction in Superlubricant MoS2 Bilayers

Molecular dynamics simulations, performed with chemically accurate ab initio machine-learning force fields, are used to demonstrate that layer stiffness has profound effects on the superlubricant state of two-dimensional van der Waals heterostructures. We engineer bilayers of different rigidity but identical interlayer sliding energy surface and show that a 2-fold increase in the intralayer stiffness reduces the friction by a factor of ∼6. Two sliding regimes as a function of the sliding velocity are found. At a low velocity, the heat generated by the motion is efficiently exchanged between the layers and the friction is independent of the layer order. In contrast, at a high velocity, the friction heat flux cannot be exchanged fast enough and a buildup of significant temperature gradients between the layers is observed. In this situation, the temperature profile depends on whether the slider is softer than the substrate.

The configurations are then obtained by random displacement of the atoms with a maximum amplitude of δ=0.05 of the strained lattices. DFT simulations are carried out with VASP 2 using the PBE exchange-correlation functional and projector augmented wave (PAW) pseudo-potentials. The energy cutoff is 500 eV and the k-point mesh is 18×18×1 for unit-cell calculations. The k-point mash converts to 6×6×1 for the 3×3 supercell.
The potential is then fitted by using the atomic coordinates and the total energy of the 720 configurations. In order to verify the accuracy of the SNAP potential, a molecular dynamics (MD) simulation is carried out for a 3×3 supercell at a temperature of 300 K.
Snapshots of the atoms are taken every 10 ps for 1 ns, after an equilibration of 0.5 ns. The resulting 100 configurations form the test set. The SNAP-predicted total energies and atomic forces of the test sets are compared against the DFT values, and the results are shown in Fig. S1. Although the forces are not used for the fitting, we compare the forces calculated by SNAP with the DFT ones on the training set. The root-mean-square errors of energies and forces are listed in Table 1. Table 1: The root-mean-square error (RMSE) of total energies and forces over the training and test sets. ∆E is the RMSE of energy and the unit is meV/atom. ∆f is the RMSE of the x, y and z components of atomic forces and the unit is meV/Å. The SNAP is only fitted with total energies.   Table 2.
The inter-layer interaction is computed through a Lennard-Jones (LJ) potential. The C 6 coefficients are extracted from the Tkatchenko-Scheffler van der Waals (TS-vdW) corrections with iterative Hirshfeld partitioning. FHI-AIMS 3 is used with the "tight" version of its   In order to find the possible supercell configurations of two identical hexagonal 2D lattices, as shown in yellow and purple in Fig. S4, the twist angle θ needs to be searched between 0 • and 30 • . If we constraint the supercell to be hexagonal and without in-plane strain, the smallest supercell areas will be found for 21.8 • ( √ 7 × √ 7). Figure S4: Sketch of supercell with twist angle 21.8 • .

Calculation of the thermal conductance across the layers
Since the out-of-plane thermal conductivity of a bilayer is ill-defined, we instead compute the interfacial thermal conductance, G. This simply relates the heat flux, q, with the temperature difference ∆T , q = G∆T , and it can be extracted from MD simulations. We first equilibrate the two layers at 400 K and 200 K, respectively. In the absence of external thermal reservoirs, the temperature difference between the two layers decays exponentially in time from its initial value, ∆T 0 =200 K (see Fig. S5), as where C v is the specific heat. G can then be extracted by monitoring the time-evolution of ∆T . In performing the fit C v =75.75 J/mol·K −1 is calculated from the total energy fluctuations. As expected G is found not to depend on the direction of the heat flux and the computed values (per unit area) are reported in Table I in the main text. Figure S5: (Color online) Time evolution of the temperature gap between the two layers of a bilayer, following equilibration at 400 K and 200 K, respectively. Each curve is the average over 10 different trajectories.

Temperature dependence of η
In order to study how the friction depends on temperature, we simulate the free sliding of the N-N bilayer, where the substrate is thermalized at different temperatures. The substrate temperatures are set to range from 150 K to 450 K in 50 K steps. The friction force is fitted to (1) in the main text]. Here k is set the value at 300 K, namely 0.83, so that η is the only parameter left to be determined. This approach eliminates numerical instabilities in the fit and leads to a better comparison between different temperatures. The η(T ) curve is then fitted again to a power-law η = C · T n , and we obtain n=1.6. Our results are shown in Fig. S4.