Atomic Insights into Fracture Characteristics of Twisted Tri-Layer Graphene

Graphene twistronics have recently gained significant attention due their superconductive behavior as a consequence of their tunable electronic properties. Although the electronic properties of twisted graphene have been extensively studied, the mechanical properties and integrity of twisted trilayer graphene (tTLG) under loading is still elusive. We investigated the fracture mechanics of tTLG with a twist angle of ±1.53° utilizing molecular dynamics simulation. This twist angle was chosen because it is known to exhibit highly superconductive behavior. The results indicate that tTLG does not preserve the excellent mechanical properties typically associated with graphene, with toughness and fracture strain values much lower in comparison. The Young’s modulus was an exception with values relatively close to pristine graphene, whereas the tensile strength was found to be roughly half of the intrinsic strength of graphene. The fracture toughness, fracture strain and strength converge as the crack length increases, reaching 0.26 J/m3, 0.0217 and 39.9 GPa at a crack length of 8 nm, respectively. The Griffth critical strain energy is 19.98 J/m2 and the critical stress intensity factor Kc is 4.47 MPa M1/2, in good agreement with that of monolayer graphene in the experiment. Our atomic insights might be helpful in the material design of twisted trilayer graphene-based electronics.


Introduction
Graphene is a carbon allotrope made up of single layer co-planar atoms arranged in a hexagonal lattice. Graphene is mainly known for its excellent mechanical properties coupled with its extraordinary thermal and electrical conductivity. Monolayer graphene was designated as the strongest material ever tested, with a Young's modulus of 1 TPa [1]. Moreover, graphene has great potential in the electronics industry. Twisted layered graphene, or the field of "twistronics", consists of layers of graphene that have been rotated relative to each other, as illustrated in Figure 1 [2]. The emergence of Moiré patterns at specific angles known as "magic angles" changes the material's electronic properties [3,4] due to the lower energy requirements for free electron tunneling between the graphene layers [5]. Jarillo-Herrero et al. discovered that at 1.1°, the band of the electronic structure became flat near zero Fermi energy, causing unconventional superconductive behavior [6]. This discovery could open the possibility of superconducting transistors to be constructed where the level of superconductivity could be varied by manipulating the intensity of the applied electric field. The discovery of superconductivity in twisted bilayer graphene was groundbreaking in the field of material science and was followed by a sudden rise in research concerning twistronic materials [7][8][9][10][11]. The scope of this discovery also extended to other Moirépatterned 2D materials that exhibit superconductivity, such as twisted bilayer tungsten diselenide (WSe 2 ) and boron nitride (BN) between layers of graphene [12,13].
In 2019, Ashvin Vishwanath of Harvard university predicted that multilayer graphene with alternating twist angles would also exhibit superconductivity. This was very recently confirmed, and it was found that at an average angle of ±1.56°, twisted trilayer graphene exhibited extraordinary superconductive properties and considering its low charge density, it had the strongest electron pairing ever observed [14]. tTLG was noted As superior to bilayer graphene because it was superconductive at a higher temperature and had better tunability in terms of its electronic structure and superconducting characteristics than bilayer graphene [15]. The electrical properties of tTLG are already well established, however, there is a research gap concerning the mechanical behavior of this material [16][17][18].
Hence, the purpose of this research was to study the fracture mechanics of tTLG. As demonstrated in our earlier work concerning graphene [19][20][21][22][23], this would allow us to better understand the material's viability in terms of its potential applications. To gain atomistic insights, we used molecular dynamics simulations, which are powerful and insightful for various atomic-level large-scale finite-temperatures applications [24]. The remainder of this paper is structured as follows: in Section 2, the model and technique are provided; Section 3 contains the findings and discussions; and Section 4 includes the conclusions.

Initial Structure
To develop a twisted trilayer structure, we opted to rotate the top and bottom layers by ±1.53°, respectively, to test the material when it exhibits superconductive behavior [25]. The 1 × 1 primitive cell structure contains 8326 atoms equally divided into the three layers, which is sufficient to map out a complete repeating cell of such a low twist angle. The X and Y axes had periodic boundary conditions whereas the Z axis was fixed to accurately model the 2D graphene structure. The size of the primitive cell was 9.05 nm and 7.84 nm along the X and Y directions, respectively.
The cracked structure was obtained by removing atoms from every layer of pristine tTLG within a rectangular prism with dimensions of 2a 0 × 0.5 nm, where 2a 0 is the crack length. As mentioned in the literature, the crack length is perpendicular to the X axis which underwent tensile stress and is denoted by the armchair direction. The pristine and cracked structures are displayed in Figure 2 with the top and isometric views, respectively.

Interatomic Potential
The potential that portrays the microscopic interactions between the carbon atoms in the graphene sheet is the adaptive intermolecular reactive bond order (AIREBO) [26]. This force-field allows breaking and forming carbon-carbon covalent bonds, enabling the precise visualization of the atomic interactions within severe environments, including fracture. Therefore, this potential was suitable for investigating the mechanical properties of monolayer graphene [27,28]. In previous literary work, multilayer graphene was generally modeled by the adaptive intermolecular reactive bond order potential (AIREBO) [29][30][31]. However, it is known that while it accurately models covalent bond breaking and forming, it is not suitable for studying fracture behavior and interlayer van der Waals interactions at high pressures. The work of Mattsson et al. on Polyethylene demonstrated these limitations [32]. In order to obtain accurate results, it is critical to consider van der Waals interactions when modeling layered graphene systems [33,34]. Therefore, we utilized the modified AIREBO potential (AIREBO-M) where the LJ term was replaced with the nondivergent Morse potential. The inclusion of the Morse potential was shown to be valid for bilayer graphene over a greater pressure range while preserving performance at ambient conditions [35]. The potential consists of three terms and is described by the following: where the REBO term has the same functional form as the potential generated by Brenner et al. [36]. The coefficients in AIREBO are the same as Brenner's potential. The MORSE potential covers the long-range interactions and depends on the given cutoff. The TORSION term refers to a four-body potential that describes the different dihedral angle inclinations in a hydrocarbon arrangement. We opted to use a cutoff of 3.0 to precisely model the carbon-carbon interactions. The inner and outer cutoff distances between carbon atoms were set to 0.192 nm and 0.20 nm, respectively. This parameter was well examined to prevent artificial strengthening near the fracture [25,37]. The van der Waals interactions between layers were modeled in the form of Morse potential, which can adequately capture the characteristics of interlayer interactions of multilayer graphene systems.

Molecular Dynamics Simulations
The large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS, stable version of 3-Mar-2020) was used for all the molecular dynamics (MD) simulations. We first carried out the potential energy minimization via the conjugate gradient algorithm, followed by the geometry optimization via the NPT ensemble with a temperature of T = 300 K and a pressure of 101,325 Pa (1 atm) for 50 ps. The MD timestep of the projector integration was 0.001 ps. We carried out the tensile test simulations with a strain rate of 1 × 10 9 s −1 .

Equilibrium
The system was first equilibrated by performing an isothermal-isobaric ensemble (NPT) simulation for 50 ps at a timestep of 0.001 ps with T = 300 K and P = 0.0001 GPa to mimic ambient conditions on pristine tTLG. Figure 3 reveals that both pressure and temperature fluctuate over constant values, indicating that equilibrium has been reached. The model was then subjected to a tensile load during the simulation with a strain rate of 1 × 10 9 s −1 to observe the stress-strain characteristics of our model. L(t) = L 0 exp (n × t) is the equation for elongation, where L 0 is the starting length, n is the strain rate, in the unit of s −1 and L(t) is the length of the simulation box at time t [38]. The stress-strain curve plays a major role in characterizing the mechanical properties of a structure and will be used to validate the model by comparing the Young's modulus with the experiment.

System Size Effects
The periodic boundary conditions in our simulations cause replicate image cells to interact with the simulation cell. This can drastically skew the results if our simulation size is not large enough relatively to the crack size. To analyze the effects of the simulation size, we constructed six different structures of varying sizes ranging from 2 × 2 (32,843 atoms) to 10 × 10 (832,133 atoms) while keeping a constant crack size of a 0 = 4 nm (so the full crack length is 2a 0 = 8 nm). Figure 4 illustrates the stress-strain relation in the different cell sizes. The graph indicates that sizes larger than 2 × 2 exhibit similar fracture behavior with a linear increase in stress initially followed by plastic behavior until the material fractures. The 2 × 2 system behaves erratically and we suspect that the structure is not large enough to accurately depict the tensile fracture of the tTLG structure. Figure 5 indicates the effects of the system size on Young's modulus (GPa), strength (GPa), toughness (J/m 3 ) and fracture strain. The Young's modulus seems to be converging as the system size becomes larger, increasing from 584 GPa to 964 GPa, resulting in a 65% increase over our system size range. The toughness and fracture strain follow a general linear trend with the system size, increasing from 0.214 (J/m 3 ) to 0.409 (J/m 3 ), and 0.021 to 0.027, respectively. This signifies a 91.1% increase in toughness and a 28.6% increase in fracture strain for this system size range. The strength also seems to generally increase from 36.0 GPa to 50.1 GPa with a 39.1% increase. However, there is a sharp drop in strength for the 6 × 6 structure which can be treated as an anomaly. Overall, the mechanical properties tend to improve as the system size increases, which is expected as the crack size occupies a smaller space relative to the total simulation size and is projected to be closer to the values observed under pristine conditions.

Crack Size Effect
The size of the crack plays an important role in the characteristics of a material [39][40][41][42]. In real applications, fracture properties carry great significance as cracks and defects are common in materials. Cracks also exist in various different sizes and affect mechanical properties to varying degrees. Therefore, by simulating such features, we will have a better understanding of the engineering qualities of this material. We varied the half length a 0 of the pre-crack. We generated a pre-crack in the shape of a rectangular prism with dimensions 2a 0 × 0.5 nm, where the half crack size a 0 varied from 0.5 nm to 4 nm. To evaluate the material's Mode-I stress intensity factor, the length of the crack was perpendicular to the axis of the strain [43,44]. The tensile test was then performed uni-axially in the zig-zag direction to assess the stress-strain response of the pre-cracked structures. Our results of the stress-strain relationships of the eight structures are displayed in Figure 6 to examine the crack size effect.
The eight different structures share the same stress-strain pattern and trend: initially, the stress values linearly rise with strain, plastically deform and then abruptly fracture after reaching the maximum tensile strength. The plastic component of the stress-strain curve reduces in range as the crack length increases and the structure starts behaving in a more brittle manner. Interestingly, the plastic behavior starts at around the same strain values across the different systems, irrespectively of the crack length. Figure 6. Crack size effect on the stress-strain relationship: A rectangular crack is generated in tTLG before tensile simulations under the conditions of temperature T = 300 K, pressure P = 0.0001 GPa, system size 4 × 4, and strain rate 1 × 10 9 s −1 .
The relations between the crack length and the toughness J/m 3 , fracture strain, strength GPa and Young's modulus GPa are demonstrated in Figure 7. The results are compared and benchmarked against pristine monolayer graphene. The toughness, fracture strain and strength appear to be converging as the crack length increases, reaching 0.26 J/m 3 , 0.0217 and 39.9 GPa at a crack length of 8 nm, respectively. The drop in these mechanical properties quantifies a reduction of 83.46% in toughness, 35.03% in strain and 84.06% in strength values over the crack length range. The toughness and fracture strain are quite low compared to pristine monolayer graphene with fracture strain values quite comparable to the values found in highly brittle materials such as silicate glasses [45]. These findings can potentially limit this material's applicability in cases where graphene usually excels. The Young's modulus decreases in a roughly linear manner as the crack length increases, from 981 GPa to 857 GPa, representing a 12.64% decrease. These values are quite comparable to pristine graphene's Young's modulus of 1005 GPa.

Fracture Toughness
We studied the fracture toughness of tTLG in terms of critical stress intensity factor K c under varying crack sizes. The stress intensity factor enabled us to gauge the properties and trends of tTLG under different conditions. The fracture toughness is described as an intrinsic property of a material to withstand crack growth under stress. The Griffith theory of brittle fracture was demonstrated to be applied to graphene, and thus is suitable for analyzing the fracture in our system [46]. The Griffith criterion is expressed as Equation (2) and describes the critical stress of the inception of rapid fracture. Previously, we evaluated the stress-strain response of the different structures containing a pre-crack. The results from those findings were used in the equations below to determine the critical stress intensity factor: where γ s is the surface energy density taken as 8.0 J/m 2 [46], E is the Young's modulus found to be approximately 1 TPa from pristine tTLG simulations, and a 0 = L/2 is the crack slit's half-length. The stress intensity factor, which is commonly used to describe fracture toughness, is represented by the following equation: The theoretical value for K c was calculated as 4 MPa M 1/2 while the value obtained from our MD simulations for the longest crack length of 8 nm was found to be 4.47 MPa M 1/2 , which signifies a difference of 11.75%. The theoretical value was calculated using Equation (2) and the simulation value using Equation (3). The stress intensity factor was evaluated under varying crack sizes and the results were compared with the theoretical value in Figure 8. The results suggest that K c is approaching the theoretical value of √ 2γ s E as the cracks become longer. This was expected as a longer crack more closely resembles the ideal Griffith crack which is characterized as long and sharp. These findings demonstrate that the Griffith theory can be applied to multilayered graphene structures. The corresponding critical strain energy release rate of a fracture, G c = σ 2 c πa o from MD is 19.98 J/m 2 , which is in relative agreement with an experimental study conducted on monolayer graphene, stating G c as 15.9 J/m 2 [46]. This is caused by the layers comprising of much stronger C-C covalent bonds than the weak van der Waals interactions observed between the layers. Hence, the addition of layers causes negligible effects to the fracture response. Coupling numerical data with visualization is required for the complete analysis of a material's fracture mechanics. Hence, the 8 nm cracked tTLG structure is displayed at different stages of tensile loading in Figure 9. After equilibrium, the initial structure (a) comprises of interlayer bonds around the crack boundaries as the carbon atoms reach a more stable configuration due to the removal of atoms forming the crack. Wrinkle formation is observed in (c) which coincides with the start of the material's plastic behavior. This topological wrinkling phenomenon is well documented to also exist in other graphene systems [47][48][49][50][51]. The amplitude of the wrinkles increase with strain as the crack propagates in the direction of the crack length until the material fractures at a strain of 0.022 (d). This trend of strain-correlated corrugation was also observed in monolayer and bilayer graphene, and has been demonstrated to influence their electronic and chemical properties [52][53][54][55]. The same fracture mechanism is observed in all the other samples tested consisting of different crack lengths, system sizes and strain rates.
The interlayer interactions play a key role in the characteristics of the material. The presence of sp 3 bonding between the layers has been shown to have a minor negative impact on the Young's modulus, ultimate tensile strength and the facture strain, while significantly increasing the interlayer shear modulus and load transfer rate, resulting in higher axial compression stability [56]. Mirparizi et al. demonstrated that interlayer shear interactions have also a considerable effect on the vibrational behavior of bilayer graphene and increase the material's natural frequencies by restricting lateral displacements [57]. Figure 10 shows carbon atoms near the crack site after equilibrium (a and b) and just before failure (c and d) from a top-down and side view. Interlayer sp 3 bonding was observed near the crack site and increased in spatial density as the structure approached failure. The ductile fracture behavior identified in our system is similarly connected to a higher interlayer bonding density [19].

Conclusions
We investigated the mechanical properties of pre-crack twisted trilayer graphene via molecular dynamics simulations. The pre-crack tTLG graphene systems exhibited a rapid brittle fracture behavior with toughness and fracture strain values which were orders of magnitude lower than that of monolayer graphene and on par with those observed in highly brittle silicate glass materials. Similarly, the fracture stress was found to be roughly half of that the intrinsic strength of graphene, whereas the Young's modulus was favorable with values relatively close to pristine graphene. The mechanical properties seem to be converging at longer crack lengths apart from the Young's modulus. The fracture toughness described by the critical stress intensity factor K c was consistent with the theoretical results with a difference of only 11.75% using the Griffith equations. Interestingly, the associated critical strain energy release rate of fracture G c was found to be in good agreement with G c values of monolayer graphene [46]. Overall, this study demonstrated the applicability of Griffith criteria for multilayer graphene systems and characterized the fracture response of tTLG, which could be helpful in determining the viability of the material for various applications. Data Availability Statement: The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.