Registry-Dependent Peeling of Layered Material Interfaces: The Case of Graphene Nanoribbons on Hexagonal Boron Nitride

Peeling of layered materials from supporting substrates, which is central for exfoliation and transfer processes, is found to be dominated by lattice commensurability effects in both low and high velocity limits. For a graphene nanoribbon atop a hexagonal boron nitride surface, the microscopic peeling behavior ranges from stick-slip, through smooth-sliding, to pure peeling regimes, depending on the relative orientation of the contacting surfaces and the peeling angle. The underlying mechanisms stem from the intimate relation between interfacial registry, interlayer interactions, and friction. This, in turn, allows for devising simple models for extracting the interfacial adhesion energy from the peeling force traces.


Convergence of the quasi-static simulations with respect to the force criterion
All quasi-static simulations presented in the main text were performed using the FIRE algorithm 1 with a force criterion of 10 -3 eV/Å. To verify the adequacy of this criterion we repeated some of the quasi-static peeling simulations of the aligned ( = 0°, Figure S1a) and misaligned ( = 90°, Figure S1b) contacts formed between a graphene nanoribbon (GNR) and a hexagonal boron nitride (h-BN) surface with a peeling angle of = 90° using a tighter force criterion of 10 -6 eV/Å. As can be seen from Figure S1, both force criteria give similar peeling traces, indicating the convergence of the simulation results with respect to this parameter. Figure S1. Convergence of the quasi-static peeling force trace with respect to the force criterion used in the FIRE algorithm 1 for two different relative orientations of the GNR and h-BN lattices: (a) = 0°, and (b) = 90°. A peeling angle of = 90° is used for both systems.

Effect of vertical distance criterion on the calculation of the contact length
As is common in atomic-scale simulations, the contact area between a surface and a substrate can be defined as the area covered by surface atoms that feel an overall repulsive force from the substrate.
This corresponds to distances smaller than the equilibrium distance ( eq ) . 2 In the GNR/h-BN system studied in this work, eq = 3. ) while recalculating the contact length using the same approach. The results of both calculations are compared in Figure S2 with those obtain using the cutoff distance applied in the main text (4 Å). From this figure we can conclude that choosing a value of 3.3 Å results in pronounced fluctuations and that increasing cut to 8 Å induced only a slight shift in the estimated value of the contact length with no effect on its qualitative behavior during the entire simulated peeling process. This justifies the used of the value cut = 4 Å. According to the model suggested by Gigli et al, 3 the shape of a GNR during peeling from a rigid substrate can be described by the following equation: , (S1) where = + sin( ) and = eq + [1 − cos( )]. This equation describes the GNR as three connected sections: (i) a flat horizontal section in contact with the substrate, (ii) a circular section at the region of detachment, and (iii) a leading inclined straight section. To check the applicability of Eq. S1 to the studied GNR/h-BN systems, we extracted snapshots of the peeling process and fitted them to Eq. S1, the results are illustrated in Figure S3. Figure S3. The shape of a 20 nm long GNR at various times ( ) during quasi-static peeling simulations for the aligned ( = 0°, upper panels) and misaligned ( = 90°, lower panels) GNR/h-BN contacts, performed using the simulation protocol presented in the main text. Black squares and red lines correspond to the simulation data and fitting curves based on Eq. S1, respectively. The values of the fitting parameters are listed in each panel. The height of the GNR ( ) is calculated by averaging the atomic height along the narrow dimension of the GNR (excluding the passivating hydrogen atoms) and over a single axial unit-cell. The position of each axial unit-cell along the GNR ( ) is calculated as the distance of its center-of-mass from the ribbon's trailing edge.
Interestingly, although Eq. S1 is proposed only for superlubric contacts with smooth-sliding peeling behavior, this formula is found to be suitable also for the aligned GNR/h-BN contacts ( = 0°), where stick-slip peeling behavior is observed. The radius of curvature of the circular GNR section is found to decrease during the peeling process. For the misaligned ( = 90°) contact, the radius of curvature saturates at around 0.9 nm, which is slightly larger than the value reported by Gigli et al. (~0.8 nm) for a GNR peeled-off of a gold substrate. 3 We attribute this small difference to the different simulation protocols and substrates. To support this point, we performed additional simulations with the same protocol used in Ref. 4 (see Figure 5 in the main text). The corresponding GNR snapshots during the peeling process appear in Figure S4, which indeed shows even better agreement between the obtained saturated radius of curvature and that reported by Gigli et al. 3 for misaligned contacts ( = 90°).
Overall, it should be noted that the exact shape of the GNR (and the corresponding elastic energy) during the peeling process has a minor effect on the obtained peeling forces. Figure S4. The shape of a 20 nm long GNR at various times ( ) during quasi-static peeling simulation for the aligned ( = 0°, upper panels) and misaligned ( = 90°, lower panels) GNR/h-BN contacts performed using the simulation protocol of Ref. 4 (see Figure 5 in the main text). Black squares and red lines correspond to the simulation data and fitting curves based on Eq. S1, respectively. The values of the fitting parameters are listed in each panel. The height of the GNR ( ) is calculated by averaging the atomic height along the narrow dimension of the GNR (excluding the passivating hydrogen atoms) and over a single axial unit-cell. The position of each axial unitcell along the GNR ( ) is calculated as the distance of its center-of-mass from the ribbon's trailing edge.

GNR
In the main text, we presented peeling force traces of a GNR from a h-BN surface obtained at a peeling angle of = 90°. For completeness, here we present complementary peeling force traces for several other peeling angles obtained using the quasi-static peeling protocol described in the main text. Figure S5- Figure S7 show the peeling force (left axis) and the contact length (right axis) as a function of peeling distance along the peeling path with four peeling angles ( = 5°, 30°, 75°, and 135°) and for three orientations ( = 0°, 90°, and 45°) of the GNR with respect to the h-BN substrate.
Interestingly, for very low peeling angles (see Figure S5a) of the aligned contact ( = 0°), the peeling force (red trace) is found to be initially independent of the peeling distance. This results from the fact that at large contact length, the static friction force is approaching a constant value. 5 At larger peeling angles, this effect is eliminated since already at the initial stages of the peeling process, the contact length considerably reduces resulting in a decrease of the static friction. For longer GNRs, the contact length during the early stages of the peeling process exceeds the threshold value below which the static friction becomes length dependent. Hence, a similar constant force regime is expected to be observed in this case also for larger peeling angle scenarios. Above a peeling angle of > 90° stick-slip peeling is completely eliminated. Instead, pure peeling with no interfacial sliding motion is obtained due to the high interfacial static friction (see Supplementary  For the misaligned ( = 90°, Figure S6) contact the overall peeling forces are found to be an order of magnitude smaller than those of the aligned contact. After initial transient dynamics a steadystate, characterized by force oscillations, is obtained for all peeling angles considered.
At a misfit angle of = 45° (see Figure S7) a more complex behavior is obtained involving intermediate peeling forces with pronounced oscillations and complex stick-slip motion, for peeling angles below = 90°. This is associated with the reported reorientation of the GNR during the mixed peeling and sliding process.

Quasi-static peeling
In the main text, all peeling simulations were performed with a 20 nm long GNR. To check the effect of the GNR's length on its peeling behavior, we performed additional quasi-static simulations for a 30 nm long GNR using the FIRE algorithm 1 with a force criterion of 10 -6 eV/Å. Both aligned (θ = 0°) and misaligned (θ = 90°) GNR/h-BN contacts were considered with two peeling angles of ϕ = 60° and 90°. The peeling force traces and the corresponding evolution of the contact length are presented in Figure S8, which shows a qualitatively similar peeling behavior as that for a 20 nm long GNR.

Constant velocity peeling
In this section, we present representative peeling force traces of a 30 nm long GNR on an h-BN substrate, obtained using the constant velocity peeling protocol with a driving velocity of 1 m/s, for various peeling angles. Figure S9 and Figure S10 show the peeling force (left axis) and the contact length (right axis) as a function of peeling distance along the peeling path for four peeling angles (ϕ = 5°, 30°, 90°, and 135°) and two orientations (θ = 0° and 90°) of the GNR with respect to the h-BN substrate. Comparing to the quasi-static results presented above, we find that the qualitative behavior of the peeling force and contact length, obtained using the two protocols, is similar. The magnitudes of the peeling forces are also similar with some differences in the amplitude of the force and the contact area oscillations.  In Figure 2 of the main text, we found that the peeling force traces for the aligned GNR/h-BN contact ( = 0°) show stick-slip behavior with a continuous overall reduction of the peeling force during the detachment of the GNR from the h-BN substrate with no apparent steady state. We attributed this behavior to the contact length dependence of the interfacial friction (and hence the peeling force) in the aligned GNR/h-BN contact, as presented in our recent study. 5 To further demonstrate this, we plot in Figure S11,

Theoretical model for extracting the adhesion energy from smooth-sliding peeling
In the main text, we showed that existing steady-state peeling models fail to describe the peeling behavior of incommensurate GNR/h-BN contact, where the steady-state peeling force is independent of the peeling angle. To describe this unexpected behavior, we develop a simple model using the method proposed by Lin et al. 6 In this model, the vertical force acting on the driving end of the GNR at steady-state is balanced by the vdW interaction between the substrate and the detached section of the GNR, which is defined by a vertical separation exceeding the equilibrium distance, . The remaining GNR section in contact with the substrate is assumed to reside at the equilibrium distance and hence does not contribute to the vertical force. Under these assumptions, the vertical force can be written as follows: where CC = To proceed, one needs to provide an explicit expression for the potential ( ). In our numerical simulations we used the dedicated interlayer potential (ILP) 7-10 as detailed in the main text. Since this potential does not allow for an analytical solution of Eq. S2 and since this equation considers large separations between the peeled GNR section and the h-BN substrate, it is sufficient to use a Lennard-Jones (LJ) type potential to describe the long-range interaction appearing in Eq. S2.
Adopting the standard 6-12 LJ potential functional form, the interaction between a carbon atom and an infinite flat substrate takes the following form 11 : In the peeled region, the repulsive energy is negligible, thus we can approximate ( ) as: where the parameter mimics the effect of taper function and the short-range Fermi−Dirac type damping function in the attractive long-range interaction term of the ILP, 9 which are replaced here by the LJ potential.
To further simplify Eq. S2, we assume that ( ) is approximately constant (see Section 3 above) and given by the angle at the middle of the peeled region, which is denoted by m . With this, Eq.
S2 can be integrated analytically, yielding: To obtain an expression for sin( m ), we use continuum elastic beam theory, where the radius of (S10) Substituting Eq. S10 into Eq. S7, we have: From the simulation results, we found that the peeled GNR section is nearly straight, thus the projected length, , can be approximated as: Substituting Eq. S12 into Eq. S11 and using the trigonometric relation sin 2 ( m ) = The solution of Eq. S13 is: Substituting Eq. S14 in Eq. S7 and using the same trigonometric relation as above, the vertical force as a function of height z can be expressed as: .
(S15) Figure S12 compares results obtained using the above model with the simulation results for the misaligned ( = 90°) contact, which exhibits smooth-sliding peeling. Eq. S15 is fitted to the simulations results (red lines in Figure S12) using a single parameter, namely the binding energy, bind , appearing in Eq. S5. The equilibrium distance eq and the constant are fixed at 3.3 Å and 0.8, respectively. A very good agreement between the simulation results and the model predictions is obtained, where the fitted binding energies agree well with the value calculated from the ILP (54.3 meV/C atom). We note that using = 1.0 results in binding energies that are ~10% lower than that predicted by the ILP. Importantly, the peeling force predicted by Eq. S15 is found to be weakly dependent on the peeling angle for the misaligned ( = 90°) contact. Figure S12. Normal peeling force ( ) as a function of vertical peeling distance (z) for the misaligned ( = 90°) contact between a 20 nm long GNR and an h-BN substrate at various peeling angles ( ) . The red and blue lines represent simulation and model results (Eq. S15), respectively.
As seen in our simulations, in the case of stick-slip peeling, steady-state is not achieved. Therefore, the model presented in section 7 above, is not suitable to treat this case. Therefore, in order to extract the adhesion energy from the stick-slip peeling force traces observed for the aligned ( = 0°) GNR/h-BN contact, we rely on energy conservation considerations rather than force balancing. 12,13 Here, the variation of the adhesion energy between two consecutive stick events is given by: where 1 , 2 and 12 are the surface energy densities of the GNR, the h-BN substrate, and the GNR/h-BN contact, respectively, ≡ 1 + 2 − 12 , Δ = Δ , and Δ is the change in contact length between the two stick events (see Figure S13a). Assuming a dissipationless process, this energy change should be balanced by the sum of changes in the elastic energy stored in the GNR, Δ el , and the external work performed by the peeling force, δW: According to our simulations Δ el is smaller than δW (see Figure S13b). Therefore, for simplicity we discard it in Eq. S17, yielding: In the linear regime, where the peeling force depends almost linearly on the stage displacement, the work during each stick phase can be approximated as: where and Δ are the compliance of the system and the change in the peeling force, respectively. Therefore, the adhesion energy density in each stick step is given by: Eq. S20 allows for calculating the adhesion energy density directly from the peeling force traces.
Due to variance between different stick events, a more reliable estimation may be obtained by averaging the result of Eq. S20 over many stick events along the peeling force trace. Figure S13c shows the calculated averages and the corresponding error bars extracted from the variance. Figure S13. Extraction of the adhesion energy of the aligned GNR/h-BN contact ( = 0°) under stick-slip conditions. Panels (a) and (b) correspond to the peeling force (red) and energies as a function of peeling distance at a peeling angle of = 90°, respectively. Panel (a) also shows the contact length (right vertical axis, green) as a function of peeling distance. The parameters labeled in panel (a), i.e., , Δ and Δ correspond to the system compliance, the change in the peeling force, and the change in the contact length, respectively, during a given stick event. Panel (b) shows the profiles for the total potential energy ( tot , red line), the interlayer energy ( ILP , blue line), the intralayer (elastic) energy ( REBO , black line), and the work of applied force (stored in the driving springs) along the lateral pulling, x, ( sx , green line) and vertical, z, ( sz , orange line) directions. Zero energy corresponds to the state where the GNR is fully detached from the h-BN substrate. Panel (c) shows the estimated adhesion energy as a function of peeling angle using Eq. S20.
Eq. S20 can be rewritten as: 13 showing that the adhesion energy can be extracted by recording Δ values as a function of √2Δ / from many peeling traces of various peeling angles, as illustrated in Figure 3 in the main text. A linear fit of this data gives a reasonable estimation of the adhesion energy for the GNR/h-BN heterojunction.