Introduction

The nucleation and growth of cracks in solids beyond their strength is controlled by the loading conditions that extend to the structural scale as well as the atomic-level fracture processes, which are localized at the crack tip. Forecasting the development of fracture remains to be challenging due to the multiscale nature and nonlinearity arising at crack tip, even with the recent advances in microstructural characterization and in situ instrumentation1,2,3,4. Neither continuum nor atomistic descriptions fail to fully rationalize the fracture processes, tackling the microstructural and physical complexities such as the interaction between the cracks and material imperfections (e.g., impurities, implanted defects, grain boundaries), and the intrinsic oscillatory dynamical instability5,6,7. By excluding the microstructural complexity, two-dimensional (2D) crystals such as graphene offer an ideal platform to explore the fracture processes. Their single-atom thickness allows each atom of the material to be exposed for experimental characterization. Events such as crack advancing, healing, structural destruction and reorganization can be identified using instruments such as the transmission electron microscope (TEM)8,9,10.

The theoretical foundation of fracture mechanics was laid by Griffith, who proposed that during crack extension, potential energy G is released to cleave new surfaces with the energy density Gc = 2γ, although in practice we usually have Gc > 2γ due to the presence of dissipative processes in fracture5. The original work was developed for brittle, isotropic solids, which could be extended to materials with anisotropy in the mechanical responses and fracture toughness. With anisotropy in the fracture toughness Gc(θ) accounted for, cracks do not necessarily follow the easy direction determined by G(θ)11. However, although the fracture strength and toughness of 2D crystals were reported for flaw-containing samples7,12,13,14,15, the anisotropy in Gc(θ) has only been discussed in theoretical studies16.

Graphene is a 2D crystal with the honeycomb lattice, which endows isotropic linear elastic responses at infinitesimal strain17. Material anisotropy in the nonlinear elastic responses arises under uniaxial tensile strain of ~15%, approaching the strain to failure16,18. The fracture toughness also demonstrates dependence on the orientation of cleaved edges19. However, most experimental studies on the elasticity are limited below ~5% strain because of the imperfections created during the growth process or sample preparation in prior to the tests15. Strain at this level cannot activate the nonlinearity in the regime of isotropic elasticity18. Graphene can thus be considered as isotropic and with elastic responses before fracture occurs. The failure behaviors of graphene were explored by loading the samples through nanoindentation, electron irradiation, or peeling20,21,22,23, showing that the cleaved edges are predominantly aligned to the armchair or zigzag edge, signaling anisotropy in the fracture toughness. Unfortunately, the stress state at the crack tip in these studies cannot be determined, not allowing further analysis with the effects of the loading conditions accounted for9,10,24,25,26. To address this issue, in situ tests of suspended graphene using a micromechanical device in a scanning electron microscope (SEM) were carried out, demonstrating the capability to measure the mechanical behaviors of 2D materials under well-defined loading conditions15,27.

The edge energy densities γ of the armchair and zigzag edges in graphene were studied theoretically by atomistic simulations using first-principles methods or empirical interatomic potentials. However, no data from experimental measurements is available to our knowledge. The values of γ or the fracture toughness Gc = 2γ in the Griffith theory reported in the literature depend on the choice of methods. The ratio between the energy density of armchair and zigzag edges (R = GA/GZ = γA/γZ) spans over the range of \(\left[0.74,1.07\right]\)12,19,20,22,28,29,30,31. The wide range of R across 1 makes it challenging to understand the fracture mechanics of graphene. The most notable fact is that first-principles calculations predict R < 1 due to electronic relaxation at the armchair edges19, in contrast to most of the empirical models where the edge energy density is largely determined by the atomic density along the edge. In this work, we conducted in situ uniaxial tensile tests of suspended graphene monolayers to constrain the value of R. The crack paths were characterized experimentally and analyzed through the extended finite element method (XFEM), resolving the microscopic processes of fracture in graphene such as crack deflection and formation of kinks.

Results

In situ uniaxial tensile tests

Displacement-controlled cyclic (loading, unloading) uniaxial tensile tests of single-crystalline graphene monolayers were performed until failure within SEM (Fig. 1, see details in Methods). Crack extension during the fracture process is manifested in the quasi-static regime, where the inertial effects are expected not to be significant due to the low loading rate and small sample size. In brittle materials, kinetic energies accumulated during crack advancing could modify the stress field at the crack tip. Moreover, crack propagation along a linear path becomes unstable above a critical velocity of about 0.73cR for linear, isotropic solids, where cR is the Rayleigh wave speed5,32. However, unstable features such as branching33,34 were not identified in our experiments, and the kinking of crack paths follows the crystallographic symmetry of the graphene lattice. These experimental evidences imply that the dynamical effects are not apparently relevant for the discussion on the crack edge orientation. To address this issue, resolving the dynamics of crack tip is necessary, which remains challenging from the experimental perspective and will be conducted in the future.

Fig. 1: In situ SEM tensile testing of single-crystalline graphene monolayers.
figure 1

a Experimental setup. Samples were placed on a push-to-pull (PTP) micromechanical device, which is actuated by an external quantitative pico-indenter15. b, c SEM images of Sample 1 at the initial (b) and fractured (c) states in the final loading-unloading cycle. d Stress-strain relations of Sample 1 (black) and Sample 2 (red), showing the brittle fracture features. The apparent nonlinearity arises mostly from the presence of cracks due to the limited size of samples. e, f SEM images of Sample 2 at the initial (e) and fractured (f) states in the final loading-unloading cycle. The dimensions and coordinates are measured on the initial state. The cleaved (upper and lower) edges along the crack paths (from left to right) were marked in (c), (f). The crystalline direction is characterized by TEM and shown in the insets of panels (c), (f). Panels (c), (f) show the non-planar or curling nature of the cleaved edges.

Full size image

Two samples were prepared and pre-stretched to be fully tightened before the tests. The size of Sample 1 is 6 × 3 μm. The stress-strain relation was recorded before brittle fracture at the strain of 4.6%, and the tensile strength is 44.3 GPa (Fig. 1d). The strain to failure and tensile strength are limited by the edge defects induced by focused ion beam (FIB) that was used to cut the sample (see the SEM image in Fig. 1b), as well as the effects of geometrical constraints at the clamping ends, or namely the clamping effects15,35. SEM characterization of the sample after fracture (Fig. 1c and Supplementary Movie 1) shows that crack advancing proceeds along the zigzag edge as characterized electron diffraction in TEM, which makes an inclination angle θZ = 11∘( ± 2∘) with the clamping ends. This result demonstrates the anisotropic fracture behaviors of graphene, which dominates over the effects of stress field that would lead to fracture along the direction in perpendicular to the maximum tension or following the principle of local symmetry11.

Sample 2 has a different geometry of 12 × 3 μm and θZ ∈ [− 5∘, 0], which was tested under three loading-unloading cycles before failure of the sample through brittle fracture. From in situ SEM characterization and the extracted stress-strain relation (Supplementary Movie 2), we identified a small, 0.84 μm-long crack during the second loading-unloading cycle at the strain of 1.33% (Fig. 1e). It was initiated at the edge of sample near the clamping end, where ripples or folds were formed after the sample-transfer and FIB-cutting processes. The strength of samples is thus reduced. The nucleated crack grows in the third loading cycle and gives rise to the fracture of sample at the strain of 1.8%. The tensile strength is 14.4 GPa (Fig. 1d). Compared to Sample 1, the crack first advances with an inclination angle of 12∘ with the clamping end instead of the low-index armchair (25∘) or zigzag (−5∘) edges. This primary crack is then deflected to be almost parallel to the clamping end, showing the dominance of the driving force G(θ) (Fig. 1f).

Modeling the fracture process

The deflection of cracks observed in the experiments signals the competition between the stress state at the crack tip and the material anisotropy in the fracture toughness, which can be measured by the directional dependence of G(θ) and Gc(θ), respectively. To analyze the crack patterns identified in the experiments, we extend the XFEM to include the anisotropy in Gc(θ). Compared to atomistic simulations, the generalized continuum mechanics framework can model the actual sample size in our experimental setup across a few micrometers. The stress intensity factors (SIFs) were calculated through the interaction integrals36. This approach was validated in our previous work, which nicely captures the crack deflection profiles in anisotropic solids37,38.

The initiation of fracture is controlled by the Griffith criterion, G > Gc, and the direction of crack advancing is determined by the generalized maximum energy release rate (MERR) criterion, that is, crack occurs at37,38,39,40

$$\theta =\arg \mathop{\max }\limits_{\theta }\frac{G(\theta )}{{G}_{{{{{{{{\rm{c}}}}}}}}}(\theta )}.$$
(1)

The direction of kinking can be selected in analogy to the Wulff construction for crystal growth, where the curve of G(θ)−1 approaches that of the \({G}_{{{{{{{{\rm{c}}}}}}}}}^{-1}(\theta )\) in polar coordinates (Fig. 2a, b)11. G(θ) is calculated from the SIFs of deflected or kinked cracks, \({\tilde{K}}_{{{{{{{{\rm{I}}}}}}}}}(\theta )\) and \({\tilde{K}}_{{{{{{{{\rm{II}}}}}}}}}(\theta )\), the values of which are derived from KI and KII defined in the local coordinates of the original crack (Fig. 2a, see details in Methods)41.

Fig. 2: XFEM models and simulation results.
figure 2

a Sample geometry, loading conditions, local energy release rate (G(θ)), and the anisotropy in fracture toughness (\({G}_{{{{{{{{\rm{c}}}}}}}}}^{-1}(\theta )\), plotted for two values across R = 0.87) at the crack tip. The stress intensity factors KI, KII and \({\tilde{K}}_{{{{{{{{\rm{I}}}}}}}}},{\tilde{K}}_{{{{{{{{\rm{II}}}}}}}}}\) (see Methods for details). b-e Crack-path predictions and standard deviations measured from the experimental data for Sample 1 (b, c) and Sample 2 (d, e), respectively. The results are obtained by varying R = GA/GZ in the range of [0.74, 1.07].

Full size image

In the honeycomb lattice of graphene, the anisotropic edge energy density can be captured in an analytical formalism19

$$\gamma (\alpha )=2{\gamma }_{{{{{{{{\rm{A}}}}}}}}}\sin (\alpha )+2{\gamma }_{{{{{{{{\rm{Z}}}}}}}}}\sin \left(3{0}^{\circ }-\alpha \right),$$
(2)

where α = θ − θZ ∈ [0, 30∘] is the angle between the cleaved edge and the zigzag motif (Fig. 2a). According to Eq. (2), the anisotropy can be classified into 3 categories as indicated by the value of R. Specifically, we have γZ ≥ γ ≥ γA for R ≤ 0.87, γ ≥ γZ ≥ γA for 0.87 ≤ R ≤ 1.0, and γ ≥ γA ≥ γZ for 1.0 ≤ R ≤ 1.15. For R > 1.15 (γA ≥ γ ≥ γZ), the zigzag direction has the lowest edge energy. However, this range of R is excluded in our study following the reported values of \(\gamma \in \left[0.74,1.07\right]\) from theoretical calculations. It is interesting to note that the edge with the lowest energy density, or the easy direction for fracture as determined by Gc(θ), is either armchair or zigzag.

A 2D finite element model was developed following the samples studied in experiments (Fig. 2a) with a mesh size of 60 nm. A mesh-independent piecewise linear crack model was implemented, as explained in our previous work37. An isotropic, linear elastic constitutive relation was used due to the limited strain range before fracture of the samples in the experimental tests. The Young’s modulus and Poisson’s ratio are 1 TPa and 0.169, respectively, and the nominal thickness is 0.335 nm42. By choosing R values in the range of \(\left[0.74,1.07\right]\), and γZ = 5.6Jm−2, γA = RγZ, we obtained a series of predicted crack paths from XFEM modeling. The results were compared with the experimental data to extract the best-fitting value of R.

The crystalline orientation (the zigzag direction) is set as θZ = 11∘ and − 5∘ for Samples 1 and 2, following the experimental setup. Experimental evidences show that the embryo cracks were nucleated from the edges due to the presence of material imperfections, and may be deflected before subsequent growth. In our model, the pre-crack was created along with the direction after the first deflection characterized in the experiment, to avoid the local effects from the edge defects created during sample preparation. The length and orientation of the pre-cracks are 2.2 μm, 11∘ (Sample 1) and 1.4 μm, 12∘ (Sample 2), respectively. The pre-cracks are located near the clamping ends and in directions not perpendicular to the loading direction. The load is applied by displacing the clamping ends, and the cracks advance by 0.03 μm at each XFEM step.

Material anisotropy in the fracture toughness

For Sample 1 (Fig. 2b, c), the pre-crack aligns with the zigzag direction. All cracks simulated based on the anisotropic MERR criterion deviate slightly from the path predicted for the isotropic model, indicating a weak anisotropy in the fracture toughness Gc(θ) (Fig. 3a). For R < 0.87, the simulated armchair direction (θZ ± 30∘) has the lowest edge energy density, which deviates from that of the pre-crack. The crack is deflected to the direction in perpendicular to the maximum tension with kinks following the armchair motif (α = ± 30∘), which highlights the effect of G(θ). This instability between neighboring armchair directions in crack advancing is attributed to the non-convex feature of \({G}_{{{{{{{{\rm{c}}}}}}}}}^{-1}(\theta )\) in the marginal condition Gc(θ)″ + Gc(θ) = 011 (Figs. 2a and 3b). For 0.87 < R < 1.0, after a few kinks between the armchair and zigzag edges that correspond to the local angular minima of force or maxima of G(θ)/Gc(θ)11 (Figs. 2a and 3b), the crack ends up with a straight line along the zigzag edge, even the armchair edge has a lower edge energy density. The alignment between the XFEM predictions and the experimental observation is improved as R increases in this range, approaching the isotropic limit. For R > 1, the zigzag edge has the lowest Gc, the crack proceeds with straight and clean edges as a result. The standard deviation decreases with R in this range (Fig. 2c). However, as this zigzag edge aligns closely to the direction in perpendicular to the load, we can only conclude that the anisotropy in the fracture toughness is weak (R ≈ 1) from the data obtained for Sample 1.

Fig. 3: Mechanisms of kinking during crack advancing in Sample 2.
figure 3

a Mechanical energy release rate and fracture toughness plotted as functions of the orientation angle θ with KII/KI = 0 and 0.25. The star marks the direction of kinks to be activated. b Kink angles (the stars) determined by the local values of KII/KI = 0 and R. The direction changes between two armchair edges (±30∘) for R = 0.85, and between the armchair (±30∘) and zigzag edges (0∘) for R = 0.94. c The next-preferred edge cleavage plotted as a function of R and the value of KII/KI at the crack tip, which changes as the crack advances due to the finite size of clamped samples. The solid line is the result obtained with R = 0.94. The open or filled symbols indicate that the crack kinks or not in the next step of advancing. d The sequences of kinking during crack advancing (from left to right), the kink angles, and the changes induced in the ratio KII/KI.

Full size image

Sample 2 (Fig. 2d, e) has a larger length-to-width ratio of 4, and the effects of stress state are expected to be enhanced. XFEM modeling using the isotropic MERR criterion shows that the crack deviates continuously with a curved path and ends up in the direction with the maximum driving force, that is, in perpendicular to the loading direction (KII = 0, following the principle of local symmetry11). The direction of pre-crack (12∘) deviates significantly from the low-index directions, featuring a mixed armchair and zigzag nature (Fig. 3b). Compared to Sample 1, this deviation in Sample 2 allows us to pinpoint the value of R by matching the crack path observed in the experiments. It should be noted that in the determination of crack paths from SEM images, the non-planar and curling effects should be excluded (Fig. 1c, f). The modeling results can well capture the crack deflection at x = ~6 μm as observed in the experiments for specific R values. The best-fitting value of R = 0.94 is finally determined from the location of this deflection along the loading direction, which downshifts toward the clamping end as R increases (Fig. 2d, e).

Discussion

The driving force G(θ) increases with the tensile load, and the crack starts to advance after the Griffith criterion is reached. This criterion characterizes the activation of deflection and kinks as crack advances, which is determined by the minimum value of Gc(θ) in the condition that the distribution of G(θ) over θ is wide enough (Fig. 3a). If the envelope of G(θ) is enclosed by that of Gc(θ), the crack has to proceed along with the armchair or zigzag edge (the one with the lowest edge energy density, from Eq. (2)). Compared to Sample 1, Sample 2 provides data that is more suitable to extract the value of R since in Sample 1 the orientation of the pre-crack aligns closely to the zigzag direction.

The ratio of local SIFs KII/KI changes as the crack advances, which may lead to the formation of kinks along the crack path and history dependence of the fracture process. The kink angle can be determined from this ratio and Gc(θ). Both armchair and zigzag kinks could be activated in the range of 0.87 ≤ R ≤ 1.07. Out of this window, only armchair or zigzag edge could occur (Fig. 3b). For Sample 2 (Fig. 2d), the path of the deflected crack is composed mainly of armchair (25∘) and zigzag (−5∘) edges, which have positive and negative values of KII, respectively. To explore the effects of local SIFs for this sample, the next-preferred kink angles are summarized as a function of KII/KI and R through the anisotropic MERR criterion (Fig. 3c). The crack that remains straight could be deflected at the critical values of KII/KI that are marked as the boundaries between the regions with cleaved armchair or zigzag edges. The critical values of KII/KI are also controlled by the anisotropy of fracture toughness, R (Fig. 3c). For R = 0.94, a crack tip in the armchair direction advances with a kink to the closest zigzag direction if the local value of KII/KI (>0) increases to 0.2. KII/KI then turns to be negative. For KII/KI < −0.07, the crack tip in the zigzag direction tilts to the closest armchair direction, and the factor KII/KI becomes positive again. These theoretical arguments are confirmed by the XFEM results (Fig. 3c, d). Our model thus allows the prediction of crack deflection and advances with the knowledge of the pre-crack direction, the crystalline orientation of the samples, and the loading conditions.

To highlight the competition between the orientational dependence in the elastic driving force and fracture toughness, we consider a model problem with a pre-crack and under the pure mode-I loading condition. The pre-crack is perpendicular to the loading direction. With the condition of R = 0.94 and that the pre-crack aligns with the zigzag motif (θZ = 0∘), the crack will advance along the zigzag direction (Fig. 4a). However, this does not necessarily mean that the edge energy density in the zigzag direction is the lowest one12. Consequently, one cannot determine the anisotropy of fracture toughness without knowledge of the stress state. Changing the crystalline orientation of the sample or the value of θZ, we find that for ∣θZ∣ < 10∘, the crack follows the zigzag motif (Fig. 4a, b). Beyond this value of θZ, the crack will be deflected into the armchair direction (θZ ± 30∘). The critical value of KII/KI to change the crack direction can be estimated theoretically as well (Fig. 4c), allowing us to extract the lattice orientation from R and the crack path, or determine the value of R by characterizing the direction of cracks with the crystalline orientation known in prior. Similar conclusions are obtained for the pre-crack aligns with the armchair direction (Supplementary Fig. 1).

Fig. 4: Fracture behaviors of a pre-crack-containing sample under mode-I loading (KII = 0), where the pre-crack aligns with the zigzag direction.
figure 4

a Crack paths predicted for different values of θZ with R = 0.94. b The first crack deflection θ in panel a at x = 20 μm. The results show that the angle θ − θZ is a constant in each of the three regions I–III, which are spaced by 30∘. c The next-preferred edge cleavage after the first deflection, plotted as a function of R and KII/KI at the crack tip. The conditions corresponding to crack advancing along the zigzag direction and without kinks in panel (a) are annotated as the white segment of the solid line plotted for R = 0.94.

Full size image

In our XFEM modeling, the loading condition is well-defined, which, however, is difficult to achieve in the experiments. Imperfect clamping boundaries would lead to deviation from the uniaxial loading conditions. XFEM simulations show that, changing the uniform loading condition to a gradient one has a minor effect on the results (Supplementary Fig. 2a). Adding shear to the displacement loading conditions could modify the stress field and crack paths, while the crack tilts (anti-)clockwise under negative (positive) shear (Supplementary Fig. 2b). The non-planar or curling nature of graphene edges (Fig. 1c, f) also brings uncertainties, as well as the determination of crystalline orientations from TEM diffraction (Fig. 1f). Modifying the zigzag direction of Sample 2 to θZ = −5∘ ± 3∘ for R = 0.94 shows similar effects as those from the shear (Supplementary Fig. 3a). The crack rotates anti-clockwise as θZ decreases. These results suggests that the experimental design on the lattice orientation and pre-crack is of critical importance to extract the degree of anisotropy in a feasible way.

In the XFEM model, samples were analyzed at the micrometer scale, but the edge energy density is defined by the arrangement of atoms at the nanoscale (Eq. (2)). Each step of the crack advance thus may accommodate a series of kinks, which cannot be assessed from the continuum-based simulations. The simulation results show that as the step size increases, the coarse-level direction of the crack path remains the same. Only armchair and zigzag edges were identified in the simulation results, in consistency with the analysis on G(θ) and Gc(θ) (Fig. 3b). However, the local features of the kinks are different, and the critical events of crack deflection may not be fully resolved (Supplementary Fig. 3b). In other words, the macroscopic orientation of crack edges observed in experiments can be different from the actual local ones. However, our simulation results suggest that the effects of the local arrangement of atoms or microstructures, although crucial for the microscale motifs of cleaved edges, are suppressed at length scales much larger than that of the microstructural features. It should also be noted that even the chiral edges in the graphene lattice can always be decomposed into short segments of armchair and zigzag edges at the atomic scale, and the bridging between the continuum and atomistic descriptions remain as an interesting topic to explore, possibly through the discussion on the energy- and strength-based criteria for fracture11,31. Experimental characterization with nanoscale resolution in combination with theoretical development could yield more insights into this multiscale problem. Unfortunately, resolving the edge structures of a crack remains experimentally challenging. SEM allows the characterization of the crack with a resolution of the edge morphology at microscale, but the local lattice orientation and the edge structures can only be identified by high-resolution techniques such as atomic force microscope (AFM) and TEM, which cannot cover a large area of samples.

Conclusion

In brief, we characterize crack development in free-standing single-crystalline graphene monolayers through in situ tensile tests and modeling. Brittle fracture initiated from damaged lattice sites at the edges and advances across the whole sample. XFEM modeling by considering the anisotropy in the fracture toughness reproduces the fracture behaviors observed in the experiments, showing the competition between the orientational dependence of elastic driving force and fracture toughness. The micromechanical approach to the uniaxial tensile tests here offers well-defined loading conditions that allow quantitative analysis of the fracture processes, in contrast to the previous studies20,21,22,23. The weak anisotropy (R ≈ 1) in the edge energy density is identified by matching the theoretical predictions to experimental evidence. The concluded value of R = 0.94 is less than unit (GA < GZ), agreeing with the first-principles predictions of γA < γZ in the Griffith scenario19,28,30,31, but contradicts with that obtained from the empirical potentials such as the adaptive intermolecular reactive empirical bond order (AIREBO) potential12,19,31. This result also aligns with the experimental fact that zigzag edges dominate in the kinetic regime of graphene growth43,44, reporting one of the key material parameters in understanding the fracture of crystals. The combined experimental and theoretical approach presented here can be generalized to other 2D crystal, and provides a platform to explore fracture mechanics behaviors (e.g., deflection, kinking) in situ and with the atomic-scale resolution.

Methods

Sample preparation

Monolayer graphene was synthesized on the copper foil by chemical vapor deposition (CVD). The film were coated by PMMA with thickness <100 nm, and transferred to the PTP device through chemical etching. First, to release graphene from the copper foil, the sample was immersed into a Cu etchant. Secondly, graphene with PMMA was fished up with a piece of polyethylene terephthalate (PET) and released into distilled water for three times to remove Cu etchant thoroughly. Thirdly, graphene with PMMA coating was fished onto the PTP device, and then left overnight in a dry cabinet to assure that graphene is attached firmly with the PTP device. Fourthly, PMMA was dissolved by acetone in a critical point dryer, which reduces the surface tension of acetone and protects the suspended graphene on the gap in the PTP device. Finally, the suspended graphene was cut into a ribbon shape by the FIB using a low electric current to reduce the beam-induced damage in the sample. The remaining part covers the surface of the PTP device, providing robust clamping for the tensile tests.

In situ SEM tensile tests

Graphene on the PTP device was tested using a Hysitron pico-indenter (PI85) inside a FEI Quanta 450 SEM. The PTP device loaded with graphene was firstly mounted onto the PI85 testing platform. The platform was then mounted into the SEM chamber for in situ tensile testing. Displacement loading was applied at a constant rate of 1 nms−1, which corresponds to the strain rate of 3.3 × 10−4s−1 to our samples. The load and displacement were recorded with the indenter transducer. Videos were recorded in situ by SEM that operates at 5–20 kV to reduce the electron beam effects. Tensile strain can be well controlled and measured directly from the SEM image sequences in the videos. The stress was calculated using the area computed from the width in the middle of the samples and the nominal thickness of 0.335 nm of graphene15.

Material characterization

The graphene samples were characterized by Raman spectroscopy with a 514-nm laser. The laser power was set as 1.7 mW to prevent thermal heating and damage in graphene. TEM characterization of the graphene samples before and after testing was carried out using a JEOL 2100F TEM. The crystalline structure of graphene ribbons suspended on the gap was characterized by TEM diffraction analysis. The selected area electron diffraction (SAED) shows only one set of the hexagonal patterns, which confirms the single-crystalline nature of the tested samples.

Energy release rate calculations

To calculate the values of G(θ), the SIFs at kinks in mode-I and mode-II were calculated as \({\tilde{K}}_{{{{{{{{\rm{I}}}}}}}}}(\theta )\) and \({\tilde{K}}_{{{{{{{{\rm{II}}}}}}}}}(\theta )\), respectively, which are measured at an angle of θ measured in the local coordinates defined by the original crack (Fig. 1a). Amestoy and Leblond41 applied complex analysis to this problem and reported the results in series expansion. The elastic driving force of fracture is

$$G\left(\theta ;{K}_{{{{{{{{\rm{I}}}}}}}}},{K}_{{{{{{{{\rm{II}}}}}}}}}\right)=\frac{1-{\nu }^{2}}{E}\left[{\tilde{K}}_{{{{{{{{\rm{I}}}}}}}}}^{2}(\theta )+{\tilde{K}}_{{{{{{{{\rm{II}}}}}}}}}^{2}(\theta )\right].$$
(3)

The SIFs can be calculated at θ = mπ from their values at θ = 0 (KI, KII), that are

$${\tilde{K}}_{{{{{{{{\rm{I}}}}}}}}}(m\pi )={F}_{11}(m){K}_{{{{{{{{\rm{I}}}}}}}}}+{F}_{12}(m){K}_{{{{{{{{\rm{II}}}}}}}}},\\ {\tilde{K}}_{{{{{{{{\rm{II}}}}}}}}}(m\pi )={F}_{21}(m){K}_{{{{{{{{\rm{I}}}}}}}}}+{F}_{22}(m){K}_{{{{{{{{\rm{II}}}}}}}}}.$$
(4)

The coefficients expanded up to the 20th order are

$${F}_{11}(m)= \; 1-\frac{3{m}^{2}{\pi }^{2}}{8}+\left({\pi }^{2}-\frac{5{\pi }^{4}}{128}\right){m}^{4}+\left(\frac{{\pi }^{2}}{9}-\frac{11{\pi }^{4}}{72}+\frac{119{\pi }^{6}}{15360}\right){m}^{6}\\ +5.0779{m}^{8}-2.88312{m}^{10}-0.0925{m}^{12}+2.996{m}^{14}\\ -4.059{m}^{16}+1.63{m}^{18}+4.1{m}^{20},$$
(4a)
$${F}_{12}(m)= -\frac{3m\pi }{2}+\left(\frac{10\pi }{3}+\frac{{\pi }^{3}}{16}\right){m}^{3}+\left(-2\pi -\frac{133{\pi }^{3}}{180}+\frac{59{\pi }^{5}}{1280}\right){m}^{5}\\ +12.3139{m}^{7}-7.32433{m}^{9}+1.5793{m}^{11}+4.0216{m}^{13}\\ -6.915{m}^{15}+4.21{m}^{17}+4.56{m}^{19},$$
(4b)
$${F}_{21}(m)= \; \frac{m\pi }{2}-\left(\frac{4\pi }{3}+\frac{{\pi }^{3}}{48}\right){m}^{3}+\left(-\frac{2\pi }{3}+\frac{13{\pi }^{3}}{30}-\frac{59{\pi }^{5}}{3840}\right){m}^{5}\\ -6.17602{m}^{7}+4.44112{m}^{9}-1.534{m}^{11}-2.07{m}^{13}\\ +4.684{m}^{15}-3.95{m}^{17}-1.32{m}^{19},$$
(4c)
$${F}_{22}(m)= \; 1-\left(4+\frac{3{\pi }^{2}}{8}\right){m}^{2}+\left(\frac{8}{3}+\frac{29{\pi }^{2}}{18}-\frac{5{\pi }^{4}}{128}\right){m}^{4}+\left(-\frac{32}{15}-\frac{4{\pi }^{2}}{9}-\frac{1159{\pi }^{4}}{7200}+\frac{119{\pi }^{6}}{15360}\right){m}^{6}\\ +10.5825{m}^{8}-4.78511{m}^{10}-1.8804{m}^{12}+7.28{m}^{14}\\ -7.591{m}^{16}+0.25{m}^{18}+26.6681{m}^{20}.$$
(4d)

Statistical analysis

The graphene sample was pre-stretched to be fully tightened, marked as the zero-strain state. For analysis, we assume that the clamping ends align with the x-axis, while the loading direction aligns with y. The corner is set as the origin. After failure with the nucleation of cracks, the indenter load declines. Two edges were cleaved, and the strain was released. The crack path predicted by XFEM simulations is \(\left\{{x}_{i},{y}_{i}\right\}\), which is sampled at a set of N points (i = 1, N). The standard deviation with the experimental data is calculated as \(\sqrt{{\sum }_{i = 1,N}{({y}_{i}-{y}_{i,\exp })}^{2}/(N-1)}\), where \({y}_{i,\exp }\) is the path interpolated from experimental measurements.