Thermally induced hex-graphene transitions in 2D carbon crystals

Resourceful beyond-graphene two-dimensional (2D) carbon crystals have been proposed/synthesized; however, the fundamental knowledge of their melting thermodynamics remains lacking. Here, the structural and thermodynamic properties of nine contemporary 2D carbon crystals upon heating are investigated using firstprinciple-based ReaxFF molecular dynamics simulations. Those 2D carbon crystals show distinct evolution of energetic and Lindemann index that distinguish their thermal stabilities. There are two or three critical temperatures at which structural transformation occurs for non-hexagoncontained 2D carbon allotropes. Analysis of polygons reveals that non-hexagon-contained 2D carbon crystals show thermally induced hex-graphene transitions via mechanisms such as bond rotations, dissociation, and reformation of bonds. The study provides new insights into the thermodynamics and pyrolysis chemistry of 2D carbon materials, as well as structural transitions, which is of great importance in the synthesis and application of 2D materials in high-temperature processing and environment.


Introduction
Carbon is a unique element because it plays a crucial role in the chemistry of living things. It has an electron configuration of 1s 2 2s 2 2p 2 , enabling to form diverse sp n (n = 1, 2, 3)-hybrid bonds. As a result, a variety of carbon allotropes form, for example, three-dimensional (3D) diamond and graphite found since prehistoric times and low-dimensional carbon buckyballs [1], carbine [2], carbon nanotubes (CNTs) [3], and graphene [4] discovered in the last decades.
Remarkably, graphene ushers a new era of twodimensional (2D) materials. It is a quasi 2D monolayered structure of carbon atoms that has attracted particular attention due to its unique electronic properties [5], as well as its outstanding strength, elasticity, and flexibility [6]. As an example, due to the fact that rather high energy of around 5.0 eV is required to break a sp 2 -hybrid C-C bond of graphene [7], it exhibits excellent resistance to mechanical and thermal actions suggesting its extremely high fusing degree. This mainly stems from the arrangement of the sp 2 -hybrid carbon atoms in a hexagonal honeycomb structure [8]. Beyond graphene, there are other 2D materials such as 2D MoS 2 [9][10][11][12][13][14][15] and 2D SiC [16] that have also attracted great attentions due to their unique mechanical properties [12,17,18], electrical conductivity [9], and thermal stability [16,19].
Besides, the behaviors of carbon crystals at elevated temperatures are a crucial piece of knowledge for their formation mechanisms and applications; however, they are still very limited [29,35]. As is known, the thermal characteristics of materials are greatly dominated by the atomic structural arrangements [36]. Understanding the thermodynamics and structural transition of 2D carbon crystals is critical for developing novel carbon-based platforms for emerging nanotechnological applications in hash conditions such as high-temperature processing and environments. To this end, the objective of this work is to reveal the thermal properties of the abovementioned 2D carbon allotropes ( Figure 1) with distinct atomic arrangements in terms of melting and thermodynamics by first-principlebased ReaxFF molecular dynamics (MD) simulations. This study provides critical insights into the melting thermodynamics and pyrolysis chemistry of 2D carbon crystals, which sheds light on the synthesis and applications of 2D carbon allotropes and offers guidance for heat-resistant composite designs by 2D carbon crystals.

Forcefield
In this study, the first-principle-based ReaxFF potential [37] is used to describe the interatomic interactions of 2D carbon allotropes. This ReaxFF forcefield comprises of three potential terms, namely bond-order-dependent covalent interaction, nonbonded standard Coulomb and Morse interaction terms. As a result, such forcefield is capable of mimicking chemical reactions of simulation systems during MD calculations, such as dissociation and formation of covalent bonds. Moreover, previous studies showed that the ReaxFF potential could accurately predict the chemical and mechanical behavior of various carbon-based structures [8,[38][39][40][41][42]. Here, the version of ReaxFF C-2013 potential developed by Srinivasan et al. [43] is adopted to predict the thermal properties of graphene and its allotropes. The ReaxFF C-2013 was developed on the basis of DFT data for equations of state of sp 2 -graphite and sp 3 -diamond, as well as formation energies of a variety of defects in graphene, fullerene, and amorphous carbon phases [43]. As a result, the ReaxFF C-2013 is capable of mimicking the chemistry and dynamics of carbon condensed phases, for example, thermal decompositions and C-C bond formation of hydrocarbons like graphene and fullerene. The similarity between graphene and its allotropes (pentagraphene, penta-hexoctite, HOP-graphene, H-haeckelite, O-haeckelite, R-haeckelite, T-graphene, and S-graphene) considered in this study indicates that the ReaxFF C-2013 potential should be capable of predicting their thermal and chemical properties.

MD simulations
Initially, each 2D carbon crystal is quasi-statically relaxed to a local minimum-energy configuration through the conjugate gradient method, where the convergence tolerances of energy and force are 1.0 × 10 −4 kcal/mol and 1.0 × 10 −4 kcal/(mol Å), respectively. Subsequently, the as-optimized sample is relaxed for 10 ps at zero pressure and temperature of 10 K under NPT (constant number of particles, constant volume, and constant temperature) ensemble. Lastly, the as-relaxed sample is gradually heated to an extremely high temperature of 7,000 K from 10 K by 7,000,000 MD steps, namely a constant heating rate of around 0.00986 K/fs. Periodic boundary conditions (PBCs) are imposed in the two planar directions of 2D carbon crystals to mimic infinite sheet, while non-PBC is applied in the off-plane direction. Such settings avoid any spurious boundary effects. The velocity-Verlet integration algorithm with a small timestep of 0.1 fs is used to integrate Newton's motions in the MD simulations. All MD calculations are carried out using the Large-scale Atomic-Molecular Massively Parallel Simulator software package.

Lindemann index (LI)
It is well-known that the root-mean-square relative bondlength variance, termed as LI [44,45], is a reasonable and effective method to estimate the melting temperature of diverse materials, such as 3D bulk materials, nanoclusters [46], and 2D systems [47]. In this study, the distancefluctuation of the LI is adopted to identify the melting temperature of our investigated 2D carbon crystals. For a system composed of N atoms, the local LI for the ith atom in the system is calculated as follows [48,49]: and the average LI of the system is determined by where r ij is the distance between the ith and jth atoms, N is the number of atoms, and r ij T 〈 〉 represents the global average distance of the ith and jth atoms at given temperature.

Polygon statistics
For structural characterization of 2D carbon crystals upon heating, their microstructure features of carbon polygons can be examined by exact geometric methods for capturing structural transformations. In this study, by means of the "shortest path ring" algorithm developed by Franzblau [50], a variety of carbon polygons ranging from trigon to decagon in our investigated 2D carbon structures upon heating are recognized on the assumption that the minimum and maximum carbon-carbon bond distances are 1.2 and 1.9 Å [8,38], respectively.
3 Results and discussion 3.1 Energetics Figure 2 shows the evolution of potential energy per atom (E) of the nine 2D carbon crystals as a function of temperature (T). Apparently, all 2D carbon structures show unique E-T curves. For graphene, the curve is described by that E increases linearly at low T region, but starts to positively deviate from the linear behavior approaching T of around 5,000 K, and finally rises sharply when over around 5,600 K. For structures of penta-hexoctite, the haeckelite family, and T-graphene, they show similar nonlinear characteristics of E-T curves. Within low T regions, E increases linearly with increasing T. Upon heating at intermediate T of around 2,000-3,500 K (depending on the type of 2D carbon structures), E starts to negatively deviate the linear behavior, conflicting with the case of graphene. When T is increased to critical values, there is a turning point in E. Above which, E of H-and R-haeckelites becomes much less increased with increasing T, while penta-hexoctite, O-haeckelite, and T-graphene show a sharp decrease in the E within finite T regimes. Those indicate the occurrence of large-scale structural transformations. Thereafter, E increases nonlinearly for those structures within finite T regimes, followed by rapid increases in the E at high T. As for the 2D carbon crystal of HOP-graphene, its curve is primarily characterized by a number of turning points within the T region of around 2,400-4,600 K, implying step-by-step structural transformations. With regard to 2D carbon crystals of penta-graphene and S-graphene who possess the highest E at ground state; however, their E-T curves are highlighted by double drops of E at two critical T. This signifies the thermally induced occurrence of two obvious structural transitions in penta-graphene and S-graphene.

LI
LI is a critical indicator in revealing the microstructural changes of materials. Figure 3 shows the average LI of the nine 2D carbon allotropes as a function of T. Based on the characteristics of LI-T curves, three groups of 2D carbon allotropes can be roughly classified. The first group can be represented by graphene, HOP-graphene and H-haeckelite. For this group, the LI-T curves are smooth throughout and are described by that the LI increases linearly with low T region, but then positively deviates from the linear behavior and eventually increases steeply. In terms of the intermediate deviation region, they are sorted as H-haeckelite > HOP-graphene > graphene. The second group consists of 2D carbon crystals of penta-hexoctite, O-haeckelite, R-haeckelite, T-and S-graphene. The LI-T curves of this group are primarily characterized by that there are two critical T at which the increase in LI becomes much more significant. This suggests that those 2D carbon crystals undergo two distinct thermally induced structural changes during the whole heating process. Remarkably, S-graphene shows fluctuation in the LI within the intermediate T region, indicating complex structural transformations. The last group is represented by a 2D carbon crystal of penta-graphene. The LI-T curve of penta-graphene is characterized by four distinct stages in the change of LI. In stage I, LI linearly decreases within the extremely low T region. In stage II, reduction in the LI becomes less pronounced. In stage III, LI suddenly increases, indicating the occurrence of large-scale structural transformations. In the last stage, LI increases nonlinearly and the increase in the LI becomes more significant with increasing T, which is similar to the whole curves of the first group.

Critical temperatures for structural destabilizations
As a result of the unique composition and arrangement of carbon polygons in the 2D carbon crystals, 1-3 inflection points that correspond to critical temperatures (T c ) can be identified from the nonlinear E-T and LI-T curves depending on the structural type. The critical temperature (T c ) of the 2D carbon crystals is quantitatively determined by the intersection T of the E-T and LI-T linear-part curves, as indicated by the dash lines in Figures 2 and 3. In this study, the last T c is referred to as the melting temperature (T m ) of 2D carbon crystals. Previous studies showed that the LI of nanoparticles and homopolymers at melting temperature (T m ) varies around 0.03-0.05 [48]. In our study of hexagondominated graphene, the value of LI is found to be 0.03 at T m , similar to hexagon-dominated one-dimensional CNTs [51], indicating the accuracy of T m by LI-T curve. Table 1 lists the T c of the nine 2D carbon crystals obtained from the E-T and LI-T curves. As is seen, graphene shows only one T c (T m ) of around 5,628 and 5,568 K using the LI-T and E-T curves, respectively. Our predicted T m of graphene falls in the regime determined by previous studies [52][53][54] that reported the minimum and maximum T m of around 4,510 and 7,750 K, respectively, suggesting the reliability in predicting T m of 2D carbon crystals. As for penta-hexoctite, HOP-graphene, and H-, O,-and R-haeckelite, they exhibit two T c . The first T c varies from around 2,300 to 4,500 K, and in terms of the first T c , they are ranked as HOP-graphene > H-haeckelite > R-haeckelite > penta-hexoctite > O-haeckelite > T-graphene. This implies that tetragon in 2D carbon crystals is critical to reducing thermal stability. Whereas in terms of the second T c , they are sorted as O-haeckelite > penta-hexoctite > R-haeckelite More intriguingly, penta-and S-graphene present three T c , suggesting their complex structural transformations upon heating. In a nutshell, non-hexagon-contained 2D carbon crystals show low T m than hex-graphene and complex pyrolysis chemistry.

Thermally induced hex-graphene transitions
To quantitively understand their thermal instability, the number of various polygons in those 2D carbon crystals subjected to the heating process is countered, as shown in Figure 4. For graphene as reference (Figure 4a), the ring number -T curve is mainly featured with initial long stable plateau and then rapid deep drop of the hexagon number at high T, accompanied by slight increases of non-hexagons. This indicates that graphene is thermally destabilized primarily by the dissociation of hexagons. For other 2D carbon crystals, there are unique variations in the number of polygons with temperature, strikingly differing from that in graphene. For penta-graphene (Figure 4b), its ring number -T curves are mainly characterized by that, as the heating T is over around 1,570 K, there is a sudden drop of the pentagon, accompanied by a sharp rising in the number of hexagon and heptagon, signifying instantaneous large-scale structural transformations from pentagon to hexagon + heptagon. With a further increase of T, there is a crossover in the number of hexagons, while the numbers of pentagon and heptagon monotonically decrease.
For penta-hexoctite (Figure 4c), it is observed from the curves that above critical T, the numbers of pentagon and octagon monotonically reduce, while the numbers of hexagon and heptagon initially increase but then decrease at high temperatures. Note that the transition in the change of the number of hexagon and heptagon occurs at different T. With regard to HOP-graphene (Figure 4d), it is uniquely found from the curves that above around 2,400 K, a two-step reduction in the number of pentagon is identified, whereas for the number of octagon, it initially increases, but then monotonically decreases. Concerning the haeckelite-based structures (Figure 4e-g) subjected to heating of high T, they are primarily involved in the nonlinear changes of the number of hexagon, pentagon, and heptagon, similar to the case of penta-graphene. Moreover, they show almost synchronous changes in the number of pentagon and heptagon, indicating the formation of pentagon-heptagon pair defects.
Remarkably, as displayed in Figure 4h and i, as a result of their composition of energetically unfavorable tetragon and octagon, T-and S-graphene show more complex changes in the polygons than other 2D carbon crystals. For example, above the critical temperature, besides the significant increases in the number of hexagons, pentagons, and heptagons, there is a considerable formation of trigon and nonagon as the rapid reduction in the number of native tetragon and octagon. Additionally, the number of the newly formed pentagon is larger than that of the newly formed heptagon, indicating that they are not purely pentagon-heptagon pair defects.
In a nutshell, excluding graphene, as indicated by the maximum peak of hexagons in the curves of Figure 4, nonhexagons contained 2D carbon crystals subjected to critical elevated T are dominated by the formation of hexagons (over 80%). This clearly demonstrates that thermally induced hex-graphene transitions are a key roadmap in the course of the pyrolysis of non-hexagon-contained 2D carbon crystals.

Atomistic origins in the thermallyinduced structural transformations
To reveal their atomic origins in the thermally induced structural transformations, the evolution of molecular configurations of those 2D carbon crystals subjected to heating is captured. Based on the formation mechanisms of local polygons and the characteristics of microstructural changes, those 2D carbon crystals can be classified into several groups. Figure 5 shows the top-viewed snapshots of graphene at different T. Upon heating at the low T region, graphene is thermally characterized by local out-of-plane displacements, resulting in wrinkling morphology. At 5,000 K, graphene is an imperfect crystal with nucleation of a variety of non-hexagonal defects varying from tetragon to decagon. As the heating T reaches 5,400 K, more non-hexagons are generated and well-distributed in graphene. Once the heating T increases, a number of voids nucleate at non-hexagonal defects and further grow, resulting in the melting of graphene. Figure 6a presents the top views of penta-graphene at eight different T. Note that at an extremely low T of around 34 K, penta-graphene is globally structurally transformed via sudden changes in bond configurations, resulting in global in-plane expansion and off-plane contraction in penta-graphene. At 1,580 K, sp 2 + sp 3 mixed penta-graphene is thermally destabilized by several local nucleations of sp 2 -polygons (mainly sp 2 -pentagon, sp 2hexagon, and sp 2 -heptagon). As illustrated in Figure 6b, nucleation of an sp 2 -hexagon occurs by dissociation of C(1)-C(2), C(3)-C(4), and C(1)-C(5) bonds and formation of C(1)-C(3) and C(1)-C(4) bonds, while an sp 2 -heptagon nucleates via breaking C(6)-C(7) and C(8)-C(9) bonds and formation of C(7)-C(8) bond. However, nucleation of a local polygon rapidly proceeds to cover all the surface of penta-graphene as the heating T approaches 1,730 K, achieving a second large-scale structural transformation to form an sp 2 -carbon sheet, consistent with the previous study [55]. With further increase of the heating T to around 4,900 K, the structural change is dominated by large-scale structural transitions of topological motifs from pentagon and heptagon to hexagon through the roadmap of dissociation of three bonds and formation of two bonds involved in a pentagon as indicated in Figure 6b. Finally, over 4,900 K, similar to graphene, the destabilization of the as-formed sheet through nucleation of voids and their growth due to thermal disintegration, resulting in the melting of penta-graphene. Figure 7a and Figure S1a show snapshots of pentagon + hexagon + octagon contained penta-hexoctite and HOPgraphene subjected to heating, respectively. As is observed, both structures are thermally destabilized by the following stages. Initially, well-distributed local nucleation of nonnative polygons occurs. However, the critical heating T of HOP-graphene responsible for such an event of microstructural change is higher than that of penta-hexoctite. In this stage, as a result of the structural changes, thermally induced buckling occurs in penta-hexoctite, as indicated by the inset of the snapshot at 3,500 K. Second, the thermally induced structural changes are dominated by the formation of hexagon accompanied by the local nucleation of non-native polygons. Third, as-nucleated pentagon and heptagon co-rearrange to form hexagons. In those above stages, "bond rotation" is the key mechanism in the largescale structural transitions, as illustrated in Figure 7b and Figure S1b. As an example, the structural transitions of penta-hexoctite result from stepwise rotations of bonds shared by pentagon-pentagon, heptagon-octagon, and octagon-octagon, whereas for HOP-graphene, they are mainly from stepwise rotation of bonds shared by hexagon-hexagon and hexagon-heptagon. Lastly, as-formed sheets primarily composed of hexagons are thermally disintegrated at several sites.  a higher rapid rate in the nucleation of hexagons than H-haeckelite during the heating process. As the nucleation of hexagons terminates, O-and R-haeckelites contain over 95% hexagon in the structures, indicating excellent thermally induced hex-graphene transition for their 2D crystals. Figure 8b and Figures (S2b and S3b) show their representative local structural changes during the heating process. Similar to penta-hexoctite and HOP-graphene, the microstructural transitions in the family of haeckelites are mainly dominated by the "bond rotation" mechanisms. For example, in O-haeckelite, the structural transformations result from two-step rotation of bonds shared by hexagon-heptagon and heptagon-octagon, whereas in H-/R-haeckelite, they are mainly from one-step rotation of bonds shared by hexagon-heptagon/heptagon-heptagon. Such difference stems from the distinct arrangement of pentagons, hexagons, and heptagons in the sheets. Figure 9a and Figure S4a display snapshots of tetragon + octagon dominated S-and T-graphene at different T, respectively. Because of the inhomogeneous arrangement of tetragons and octagons in the sheet, S-graphene shows strong thermally induced buckling morphology at low T regions, while T-graphene remains flat morphology at low T region due to its homogeneous arrangement of tetragons and octagons. As a result of their energetically unfavorable configurations of tetragon and octagon, they show scattered nucleation in the sheets at relatively low T. With increasing T, in comparison with other non-hexagon contained 2D carbon crystals, more complex polygons forms. At high T, both structures tend to form hexagon-dominated stable sheets, followed by thermal disintegration, similar to other 2D carbon crystals.   Figure S4b show the representative snapshots of local structures of S-and T-graphene during the structural changes, respectively. As is illustrated, for S-graphene, the initial structural transformations are primarily described by rotation and dissociation of bonds shared by tetragon-tetragon and bond formation in the octagon, resulting in the formation of trigon, pentagon, hexagon, and nonagon in the sheets. As a result, the microstructures in S-and T-graphene becomes close to other 2D carbon crystals. Therefore, the following structural transformations of S-and T-graphene are analogous to multi-polygon-dominated 2D carbon crystals.

Conclusion
In summary, the melting thermodynamic behaviors of graphene and its contemporary allotropes are comprehensively investigated by conducting first-principle-based ReaxFF MD simulations. All graphene allotropes show unique nonlinearity in the LI-T and E-T curves that reflect their distinct thermal properties. Differing from that of hexagon-dominated graphene, the melting behaviors of investigated graphene allotropes are primarily characterized by multiple phase transformations at critical temperatures that are greatly dictated by the polygonal compositions and arrangement. At low heating T regions, graphene, H-, O-, and R-haeckelites that are mainly composed of pentagon, hexagon, and heptagon are thermally stable structures, while penta-graphene, penta-hexoctite, HOP-, T-, and S-graphene show apparent local structural transformations. Upon heating at moderate T regions, all graphene allotropes containing non-hexagonal defects transform into hexagondominated structures via distinct mechanisms that are also controlled by their unique microstructural features such as polygonal compositions and arrangement. This indicates that thermally induced hex-graphene transitions are a key roadmap in the course of the pyrolysis of graphene allotropes. As they are heated at high T, similar to graphene, the thermodynamic behaviors of graphene allotropes are characterized by the nucleation and growth of nanovoids as a result of thermal disintegration. Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.

Conflict of interest:
Authors state no conflict of interest.