Atomic Structure and Mechanical Properties of Twisted Bilayer Graphene

We studied the atomic structure and mechanical properties of twisted bilayer graphene with a different twist angle using molecular dynamic simulations. The two layers are corrugated after energy minimization. We found two different modes of corrugation. The mechanical properties are tested both in-plane and perpendicular to the plane. The in-plane properties are dominated by the orientation of graphene. The perpendicular properties depend on the twist angle, as the larger the twist angle, the higher the intrinsic strength.


Introduction
Graphene, as a special two-dimensional material exhibiting a number of intriguing properties, has attracted extensive attention worldwide for a long time [1][2][3][4][5][6][7][8][9]. The attention even shifted to graphene-related and graphene-like materials [10][11][12][13][14][15][16][17]. The stacking of two layers of graphene also brings some new properties to this material. When one of the stacking layers gets twisted with a certain angle relative to the other, so called twisted bilayer graphene (tBLG) is formed, which has been proved to have many new properties, such as unique electrical properties [18][19][20][21][22][23][24][25], optical properties [26][27][28][29], chemical properties [30] and phonon thermal properties [31][32][33][34][35][36]. One intriguing property is the emergence of a flat band in a small twist angle [18][19][20]. These flat bands exhibit insulating states at half-filling while superconducting the state upon electrostatic doping [21,22]. The unique twist-angle-dependent van Hove singularities (VHS) influence the absorption spectrum of tBLG, which may be useful in designing graphene photodetectors [26][27][28]. Besides, the enlarged unit cell of tBLG results in the reduction of the Brillouin zone size, which in turn leads to the emergence of new phonon branches [31]. This further decreases its thermal conductivity [32,33]. The specific heat of tBLG also depends on the twist angle, especially at low temperature, which means it is possible to control the thermodynamic properties of bilayer graphene by changing the twist angle [34,35]. Those unique properties make tBLG a promising material in future electrical and optical devices. Few people have paid attention to the mechanical properties of tBLG material which could play an important role in the thermo-mechanical stability, reliability and durability of applications. The knowledge about the mechanical properties of tBLG will be of great help in applications. Although Liu has simulated the mechanical properties of tBLG, he only focused on a twist angle of 25.04 • [37]. The mechanical properties of tBLG with other twist angles are still unknown. This work will focus on the mechanical properties of tBLG with different twist angles to find out their relationship with twist angles through molecular dynamical simulation.

Methods
The tBLG model was prepared according to Uchida's method [20]. Any tBLG can be labeled by a pair of integers (m, n) as explained below. A tBLG sample was obtained by rotate a vector V(m, n) to V (n, m) in one layer of the bilayer graphene with m, n the coordinates with respect to the basis vector a 1 ( √ 3 2 , − 1 2 ) and a 2 ( √ 3 2 , 1 2 ). The superperiodicity defined by (n, m) and (− m, n + m) vectors is the super unit cell of tBLG, with a cell length L equals: The twist angle can be calculated as: The x direction is set along (n − m, n + 2m) direction. Table 1 is the list of tBLGs calculated in this work. The simulations were carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software package. We combine two different potentials to describe the interaction between carbon atoms. For the interaction in the same layer, we choose the 2nd generation reactive empirical bond order (REBO) potential [38], which can precisely describe the in-plane elastic properties of graphite and diamond. Meanwhile for the interlayer interaction, we chose the Kolmogorov-Crespi interaction potential specially designed to describe interlayer interaction for graphitic system [39]. Such a combination of potential is capable of precisely describing the mechanical properties of twisted bilayer graphene.
We carried out a tensile test simulation in the x-y plane and compression test along the z direction using molecular dynamic methods to study the elastic properties of tBLG. Before the test, all of our samples were energy-minimized, and the stress in the boundary of the simulation box was set to zero. Then the sample was further relaxed under isothermal-isobaric ensemble for 100 ps. The temperature was set to be 0.01 kelvin in all the following simulations and the external press is set to zero. After that we carry out the tensile test under the strain rate of 4 × 10 −4 /ps along x direction and the compression test under the strain rate of 1 × 10 −3 /ps along y direction. The stress in the other direction remain unchanged during the entire time.

Atomic Corrugation
After the energy minimization, a periodical corrugation forms in the surface of two layers of graphene. The amplitude of the corrugation is as small as around 0.1 A, thus the surface is actually almost flat. However, we find two different modes of corrugation. When the twist angle is large, the corrugation is in the form of Figure 1b,c. In the AA stacking areas of the tBLG, the two layers of graphene buckle out-of-plane in the opposite direction. This result is inconsistent with the results obtained by Uchida et al. [20] and Dai et al. [40]. Dai et al. called this mode the breathing mode. The reason it forms may be that the interlayer spacing in AA stacking graphene is larger than the AB/AC stacking graphene. When the twist angle is smaller than 1.6 • , the corrugation mode is different, as shown in Figure 1e,f. In the AA stacking areas, the two layers still buckle in opposite directions, same as the breathing mode. While in the surrounding area of AA nodes, a new mode of corrugation emerges. The buckling forms three peaks alternating with three valleys around each AA stacking node in one layer, and two layers buckling in the same direction. This mode is similar to the bending mode mentioned by Dai et al. For convenience, we also refer to the second mode as the bending mode. all of our samples were energy-minimized, and the stress in the boundary of the simulation box was set to zero. Then the sample was further relaxed under isothermal-isobaric ensemble for 100 ps. The temperature was set to be 0.01 kelvin in all the following simulations and the external press is set to zero. After that we carry out the tensile test under the strain rate of 4 10 /ps along x direction and the compression test under the strain rate of 1 10 /ps along y direction. The stress in the other direction remain unchanged during the entire time.

Atomic Corrugation
After the energy minimization, a periodical corrugation forms in the surface of two layers of graphene. The amplitude of the corrugation is as small as around 0.1 A, thus the surface is actually almost flat. However, we find two different modes of corrugation. When the twist angle is large, the corrugation is in the form of Figure 1b,c. In the AA stacking areas of the tBLG, the two layers of graphene buckle out-of-plane in the opposite direction. This result is inconsistent with the results obtained by Uchida et al. [20] and Dai et al. [40]. Dai et al. called this mode the breathing mode. The reason it forms may be that the interlayer spacing in AA stacking graphene is larger than the AB/AC stacking graphene. When the twist angle is smaller than 1.6°, the corrugation mode is different, as shown in Figure 1e,f. In the AA stacking areas, the two layers still buckle in opposite directions, same as the breathing mode. While in the surrounding area of AA nodes, a new mode of corrugation emerges. The buckling forms three peaks alternating with three valleys around each AA stacking node in one layer, and two layers buckling in the same direction. This mode is similar to the bending mode mentioned by Dai et al. For convenience, we also refer to the second mode as the bending mode. For tBLGs with twist angle of 1.30°, 3.89° and 4.41°, the same simulations are done with different box size respectively. The results show that the box size has a slight effect on the atomic corrugation. Another difference between breathing mode and bending mode is their buckling height, which is defined as the largest difference in z coordinate of atoms in one layer. As shown in Figure 2, the buckling height of breathing mode is sustained at around 0.11 A in a small twist angle. It starts to decrease at the angle of 7° and finally becomes zero when the twist angle reaches 20°. The bending mode only arises when the twist angle is smaller than 1.6°. The buckling height of bending mode increases rapidly with the decrease of the twist angle. It reaches 1.4 A when the angle is as small as 0.71°. Such a large buckling height may bring something new to the mechanical properties of tBLG Another difference between breathing mode and bending mode is their buckling height, which is defined as the largest difference in z coordinate of atoms in one layer. As shown in Figure 2, the buckling height of breathing mode is sustained at around 0.11 A in a small twist angle. It starts to decrease at the angle of 7 • and finally becomes zero when the twist angle reaches 20 • . The bending mode only arises when the twist angle is smaller than 1.6 • . The buckling height of bending mode increases rapidly with the decrease of the twist angle. It reaches 1.4 A when the angle is as small as 0.71 • . Such a large buckling height may bring something new to the mechanical properties of tBLG in small twist angle. The corrugations also serve as defects in the graphene layers, which may enhance the defect scattering and boundary scattering of phono and lead to the decrease in thermal conductivity of tBLG [36,41]. This is inconsistent with the results obtained by Li et. al [32] and Cocemasov et. al [35].
in small twist angle. The corrugations also serve as defects in the graphene layers, which may enhance the defect scattering and boundary scattering of phono and lead to the decrease in thermal conductivity of tBLG [36,41]. This is inconsistent with the results obtained by Li et. al [32] and Cocemasov et. al [35].

Tensile test along x direction
The energy minimized sample is further relaxed under isothermal-isobaric ensemble for 100 ps with T = 0.01K and P = 1 KPa before tensile test. Considering that tBLGs with different twist angles have different interlayer distances, and in order to ignore the debate about the thickness of graphene layers, all the stress in this test is defined as the force applied to the unit length along the y direction of the sample, which is N/m. To verify the validity of our method, we firstly carry out tensile test on single layer graphene. All of the elastic responses obey the nonlinear relationship. The function describing such a nonlinear relationship can be expressed as [8]:

= +
where is the axial stress, is the strain, E is young's modulus and D is the third order elastic modulus. The young's modulus and the third order elastic modulus can be obtained by fitting the stress-strain curve with quadratic polynomial. The young's modulus of single layer graphene along the zigzag and armchair direction are simulated to be 294 N/m and 353 N/m respectively, which match well with the experimental results [8].
Then we carry out tensile test for samples with different twist angle respectively. Graphene is an anisotropic material, the mechanical properties along different orientation is different. The results is as Figure 3a. In order to find out how much the orientation can affect the mechanical properties of tBLG, we also carry out tensile simulation to AA stacking bilayer graphene (BLG) under the same conditions. Therefore, we can compare the elastic properties of tBLG with BLG, which has the same orientation. Since each layer of tBLG with a twist angle has an orientation of /2 relative to x direction, the corresponding BLG should also has an orientation of /2. For convenience, the angle in the horizontal axis is the twist angle of tBLG, which is actually twice the orientation of BLG. The Young's modulus of tBLG and BLG both decrease with the increase of twist angle, the third order elastic modulus of tBLG and BLG both increase with the increase of twist angle. Figure 3d is the intrinsic strength obtained by = − /4 , which is the highest stress the system can support before failure. The curve is almost flat in a small twist angle, and even in the large angle region, the difference is small. The blue dot in Figure 3b and Figure 3d is the simulation results obtained by Liu [37]. Our results are approximately the same, and the differences may be caused by a different fitting process.

Tensile Test Along x Direction
The energy minimized sample is further relaxed under isothermal-isobaric ensemble for 100 ps with T = 0.01K and P = 1 KPa before tensile test. Considering that tBLGs with different twist angles have different interlayer distances, and in order to ignore the debate about the thickness of graphene layers, all the stress in this test is defined as the force applied to the unit length along the y direction of the sample, which is N/m. To verify the validity of our method, we firstly carry out tensile test on single layer graphene. All of the elastic responses obey the nonlinear relationship. The function describing such a nonlinear relationship can be expressed as [8]: where σ is the axial stress, is the strain, E is young's modulus and D is the third order elastic modulus. The young's modulus and the third order elastic modulus can be obtained by fitting the stress-strain curve with quadratic polynomial. The young's modulus of single layer graphene along the zigzag and armchair direction are simulated to be 294 N/m and 353 N/m respectively, which match well with the experimental results [8].
Then we carry out tensile test for samples with different twist angle respectively. Graphene is an anisotropic material, the mechanical properties along different orientation is different. The results is as Figure 3a. In order to find out how much the orientation can affect the mechanical properties of tBLG, we also carry out tensile simulation to AA stacking bilayer graphene (BLG) under the same conditions. Therefore, we can compare the elastic properties of tBLG with BLG, which has the same orientation. Since each layer of tBLG with a twist angle θ has an orientation of θ/2 relative to x direction, the corresponding BLG should also has an orientation of θ/2. For convenience, the angle in the horizontal axis is the twist angle of tBLG, which is actually twice the orientation of BLG. The Young's modulus of tBLG and BLG both decrease with the increase of twist angle, the third order elastic modulus of tBLG and BLG both increase with the increase of twist angle. Figure 3d is the intrinsic strength obtained by σ int = −E 2 /4D, which is the highest stress the system can support before failure. The curve is almost flat in a small twist angle, and even in the large angle region, the difference is small. The blue dot in Figures 3b,d is the simulation results obtained by Liu [37]. Our results are approximately the same, and the differences may be caused by a different fitting process. The same tendency and similarity of the three curves means that the elastic properties of tBLG are mainly caused by the orientation of the graphene layer. The twist only slightly retards the change of the properties. In the small angle region where the atomic corrugation becomes apparent, the modulus and intrinsic strength of tBLG is almost the same as AA stacking BLG. Above all, we can conclude that the atomic corrugation has little effect in the mechanical properties of tBLG. Therefore, it is safe to bring in tBLG in most of applications without worrying about the change of mechanical properties compared with graphene.

Compression test along z direction
Now we consider about the mechanical properties along the z direction. We created two walls, one beneath the lower layer of graphene and one upon the upper layer. The energy of wall-particle interaction is given by the following equation: where r is the distance between particle and wall, = 1, = 1, = 0.7. Then the upper wall moves down to a certain speed to simulate the compress process. The energy minimization and relaxation under isothermal-isobaric ensemble is also executed before the compression. The results are shown in Figure 4. The strain-stress curve doesn't exhibit a typical polynomial relationship but rather like an exponential relationship. The stress increase faster with the increase of strain, then the system breaks down in certain strain when two layers of graphene finally contact. When twist angle is larger than 3°, the intrinsic strength of tBLG is larger than AA stacking (the black curve in Figure 4a) and AB stacking (the red curve in Figure 4a) bilayer graphene, which is around 346 GPa (the red dashed line in Figure 4b), and increases with twist angle. While for tBLGs with twist angle smaller than 3°, the strength is smaller than 346 GPa, which means that the tBLGs are much "softer" in small twist angle. The same tendency and similarity of the three curves means that the elastic properties of tBLG are mainly caused by the orientation of the graphene layer. The twist only slightly retards the change of the properties. In the small angle region where the atomic corrugation becomes apparent, the modulus and intrinsic strength of tBLG is almost the same as AA stacking BLG. Above all, we can conclude that the atomic corrugation has little effect in the mechanical properties of tBLG. Therefore, it is safe to bring in tBLG in most of applications without worrying about the change of mechanical properties compared with graphene.

Compression Test along z Direction
Now we consider about the mechanical properties along the z direction. We created two walls, one beneath the lower layer of graphene and one upon the upper layer. The energy of wall-particle interaction is given by the following equation: where r is the distance between particle and wall, = 1, σ = 1, r c = 0.7. Then the upper wall moves down to a certain speed to simulate the compress process. The energy minimization and relaxation under isothermal-isobaric ensemble is also executed before the compression. The results are shown in Figure 4. The strain-stress curve doesn't exhibit a typical polynomial relationship but rather like an exponential relationship. The stress increase faster with the increase of strain, then the system breaks down in certain strain when two layers of graphene finally contact. When twist angle is larger than 3 • , the intrinsic strength of tBLG is larger than AA stacking (the black curve in Figure 4a) and AB stacking (the red curve in Figure 4a) bilayer graphene, which is around 346 GPa (the red dashed line in Figure 4b), and increases with twist angle. While for tBLGs with twist angle smaller than 3 • , the strength is smaller than 346 GPa, which means that the tBLGs are much "softer" in small twist angle. Therefore, the strength of tBLG in z direction can be tailored by change the twist angle ranging from 325 GPa to 460 GPa while leaving the elastic properties in x-y plane unchanged. This may be very useful in ballistic protection applications [42].

Conclusions
We studied the atomic corrugation and mechanical properties of twisted bilayer graphene with a different twist angle. After energy minimization, the surface of tBLG becomes corrugated. We found two different modes of corrugation. One is a breathing mode that has two layers of graphene buckling in the opposite direction. The other is a bending mode that has two layers of graphene buckling in the same direction and only exists in a small twist angle. The buckling height of bending becomes significantly large in a very small twist angle. However, the tensile test shows that such an atomic corrugation has little effect on the elastic properties of tBLG in x-y plane. The elastic properties of tBLG is dominated by the orientation of graphene. The twist only slightly retards the change of these properties. While in the z direction, the twist has a great effect on the strength, as the strength increases with the twist angle. Understanding these mechanical properties is helpful for the application of tBLG. Such properties make tBLG a potential material in ballistic protection applications.