Oscillatory motion in layered materials: graphene, boron nitride, and molybdenum disulfide

Figures 2(a), (c), and (e) show the two-dimensional potential-energy surfaces for sliding calculated using DFT. The minimum-energy configuration of graphene has the carbon atoms staggered, while the minima for BN and MoS₂ correspond to a situation with maximum overlap of unlike atoms: B over N and S over Mo. The binding energies per atom are 24.78 (graphene), 27.68 (BN), and 39.46 (MoS₂) meV. The dispersion contributions to these are 31.98, 35.22, and 48.18 meV, respectively, indicating that dispersion interactions are dominant in inter-layer binding. The potential energy surface and the binding energy of BN are similar to graphene, which is no surprise given their similar structures and polarizabilities.

Zoom In Zoom Out Reset image size

The energy per unit area along the long diagonal of the contour plots is shown in figures 2(b), (d), and (f). The maximum energy per area (energy barrier) of MoS₂ is more than three times greater than that of graphene or BN. This may be explained by the corrugated nature of the MoS₂ surface, which increases both the dispersion binding and the energy barriers, and also by the sulfur anions being more polarizable than the first-row atoms in graphene and BN. Also shown on these figures is the energy per area predicted by the empirical potentials used in the MD simulations. Comparison to the DFT reference data reveals that there is excellent agreement for all three materials, indicating that the MD simulations are able to capture the inter-layer energetics accurately.

MD simulations were first used to characterize the dynamic response of the laterally offset materials for ideally flat layers (i.e. no externally imposed surface roughness). figures 3(a), (b), and (c) show the oscillatory and self-retraction motion of laterally offset graphene, BN and MoS₂. The graphene and BN layers exhibit oscillatory motion while MoS₂ retracts to the maximally overlapping position without oscillation. Since all three materials are ideally flat, and the commensurability and temperature conditions are the same, the observed difference in the MoS₂ dynamics is attributed to the large energy barrier (see figure 2) arising from the corrugated nature of the surface. A higher energy barrier is associated with more resistance to sliding (greater friction), which can in turn be expected to damp the oscillations more quickly. This suggests that self-retraction and oscillation behavior is affected by the inherent resistance of the material to sliding (as measured by the static energy barriers), as well as external factors (e.g. temperature, roughness) as reported previously [4, 11]. While maximizing dispersion binding provides a driving force for self-retraction for all layered materials, small sliding barriers are needed to observe oscillation.

Zoom In Zoom Out Reset image size

To further analyze the effect of the potential energy barriers on the dynamics, we varied the inter-layer interaction strength by changing the depth of the potential well () in the LJ potential for each material in the MD simulations. The results are shown in figure 3. MoS₂ exhibits self-retraction, i.e. the motion is overdamped, for all potential well depths studied. However, for BN and graphene, increasing the interaction strength increases the rate at which the oscillation amplitude decreases over time. This effect is quantified in terms of the decay rate obtained by fitting the maximum position versus time data to an exponential expression. Figure 3(d) shows that the decay rate is faster for stronger interactions, depending exponentially on the barrier.

We also considered the dependence of the motion on the σ parameter. While modifying the LJ parameter changes only the energy barrier, modifying σ affects both the equilibrium inter-layer distance (i.e. the layer corrugation) and the energy barrier. Illustrative results for graphene and MoS₂ are shown in figure 4. For MoS₂, doubling σ, which increases the interlayer distance, results in the onset of oscillatory motion. The minimum-energy separation between layers becomes larger relative to the intra-layer atomic distances, making the interlayer contact less corrugated. This facilitates sliding because both the energy barrier is smaller and the corrugated nature of the surface impedes the relative motion of the layers, similar to the effect of roughness-induced friction. For graphene, which shows oscillatory motion with the original LJ parameters, we decreased σ by half. As shown in figure 4, instead of oscillating, the upper flake became trapped in a local minimum on the potential-energy surface, which is similar to the behavior seen for MoS₂ with 0.5 ${{\epsilon }_{0}}$ in figure 3(c). This occurs because, with decreased interlayer distance, the sliding barrier is so high that it cannot be surmounted and the upper layer cannot retract back to the minimum-energy position. Therefore, changing the potential-energy surface can qualitatively alter the self-retraction/oscillation behavior. In particular, a relatively flat (i.e. uncorrugated) surface with low energy barriers seems to be a necessary condition in order to observe oscillation.

Zoom In Zoom Out Reset image size

Our previous study showed that the decay rate of the oscillatory motion of graphene increases linearly with imposed surface roughness and a reduced-order model was developed to describe that motion [11]. The model is briefly reviewed here. The potential energy, E, as a function of the dimensionless offset ratio $x=L/{{L}_{0}}$, is written as

$$\begin{eqnarray}&&E={{k}_{1}}{{A}_{0}}|x|-{{k}_{2}}{{A}_{0}}(1-|x|){\rm cos} \left( \frac{2\pi {{L}_{0}}x}{p} \right).\end{eqnarray} \tag{ 2 }$$

The force can be obtained as

$$\begin{eqnarray}&&F=-\frac{{\rm d}E}{{\rm d}x}\frac{1}{{{L}_{0}}}+c{{L}_{0}}\frac{{\rm d}x}{{\rm d}t}\end{eqnarray} \tag{ 3 }$$

and the dynamics can be described by solving equation of motion:

$$\begin{eqnarray}&&F=m\ddot{x}.\end{eqnarray} \tag{ 4 }$$

In these expressions, A₀ is the area (in this case $L_{0}^{2}$), k₁ is the binding energy per contact area, k₂ is the amplitude of the oscillatory energy profile, p is the periodicity of the material structure, and c is a coefficient related to the roughness, which controls the friction. Increasing roughness increases friction which, in turn, linearly increases the decay rate of the oscillations [11].

We now examine the effect of increasing roughness on the self-retraction motion of both graphene and BN (since MoS₂ does not oscillate even for the ideally flat layers, we do not expect any change in this behavior by increasing the roughness). The results in figure 5 show that graphene and BN each exhibit a linear relationship between decay rate and roughness with similar slopes, but BN has a slightly higher decay rate at zero roughness, likely reflecting the larger energy barrier (figures 2(b), (d)). We also considered the effect of roughness on the decay rate for graphene and BN with the interaction strengths and sliding barriers artificially increased by taking $\varepsilon =2{{\varepsilon }_{0}}$. The results show that the linear relationship is retained and increasing the interaction strength causes an effectively constant shift in the decay rate to higher values. This suggests that oscillatory motion is affected by both the potential energy barriers (intercept) and the roughness (slope), at least for materials with a structure similar to graphene.

Zoom In Zoom Out Reset image size

To describe these observations, the reduced-order model presented in our previous paper [11] must be generalized. In our previous work, we observed that c depends linearly on the decay rate and then fit the dissipative friction coefficient c using the linear increase in decay rate with roughness predicted by the MD simulations. However, we observe in figure 5 that there is a finite decay rate, even for ideally flat surfaces, which is associated with the surface structure and energy barriers, as discussed in the preceding section. To capture the effects of both roughness and energy barrier on the motion, the definition of the coefficient, c, must be made more flexible. We now propose that this coefficient have the general form

$$\begin{eqnarray}&&c=A{{R}_{q}}+B,\end{eqnarray} \tag{ 5 }$$

where A and B are determined by the slope and intercept of the decay-rate versus roughness plot (figure 5), and account for the effects of imposed roughness and energy, respectively.

The A and B parameters, as well as k₁ and k₂, can be obtained from static DFT and MD data. In this work, the surface energy, k₁, was obtained from periodic DFT calculations and found to be 0.01888 eV Å⁻² (24.78 meV/atom) for graphene and 0.02030 eV Å⁻² (27.68 meV/atom) for BN. The value of k₂ was calculated as the amplitude of the oscillatory energy profile that results from subtracting a V-shaped curve (given by ${{k}_{1}}{{A}_{0}}|x|$) from the MD potential energy versus the mismatch ratio. From our MD simulations, k₂ was calculated to be 2.5 × 10⁻⁴ eV Å⁻¹and 3.2 × 10⁻⁴ eV Å⁻¹ for graphene and BN, respectively.

Table 2 shows the values of A and B obtained from linear fits to the MD data for graphene and BN with original and artificially increased inter-layer interaction strengths. The differences between all values of A (related to the slopes in figure 5) are on the same order as the standard error of the linear fit, indicating that the relationship between imposed roughness and decay rate is nearly independent of the material and sliding energetics. This contrasts with the B parameter that depends strongly on both material and interaction strength, and corresponds to the decay rate for ideally flat surfaces. We showed in figure 3(d) that the decay rate increases exponentially with the energy barrier. Hence, the coefficient c in equation (3) depends linearly on the roughness (with slope A) and exponentially on the energy barrier (through B), and these two effects are independent of each other.

The generalized analytical model can, in principle, be used to predict offset-driven motion for any material with any roughness. For example, figure 6 shows the predicted oscillation behavior of graphene and BN with the original and artificially increased interaction strengths, which are consistent with the results obtained from the MD simulations. However, the more useful feature of the model is its ability to predict offset-driven motion for cases where MD simulations are not possible. For example, the roughness that can be imposed in MD is limited by the tractable simulation size and is typically smaller than roughnesses measured experimentally. However, the reduced-order model can be used to extrapolate the behavior of graphene and BN to experimental roughness values in order to predict the self-retraction dynamics of real materials. Figure 6 also shows the oscillation behavior predicted by our model extrapolated to a roughness of 5 Å, which is consistent with the range of experimentally measured roughnesses of graphene [24] and BN [25]. For both BN and graphene, the predicted behavior is oscillatory but heavily damped, with at most one observable oscillation. This result is consistent with the previously-reported experiments on graphene [3–5].

Zoom In Zoom Out Reset image size

Our simple reduced-order model appears to capture the energetics- and roughness-dependent oscillation behavior of BN and graphene accurately. In the case of MoS₂, the overdamped behavior of the ideally flat system at any interaction strength (figure 3(c)) prevents calculation of the decay rate and, consequently, the A and B parameters for this material. The value of B is likely one or more orders of magnitude higher than those for graphene and BN. It is important to note that the energy barrier for MoS₂ is only three times larger than those of graphene and BN, which means that this is not the only factor contributing to B, and that the layer geometry (and in particular its corrugation, as illustrated in figure 1(f)) also affects this coefficient. Nevertheless, since the decay rate always increases with roughness, we do not expect this material (or others with similarly-corrugated layers) to present oscillation experimentally in any case.