Thermal Stability and Phase Transition Behavior in In-Plane Graphene/hBN/Graphene Heterostructures: A Molecular Dynamics Study

Thermal Stability and Phase Transition Behavior in In-Plane Graphene/hBN/Graphene Heterostructures: A Molecular Dynamics Study | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Thermal Stability and Phase Transition Behavior in In-Plane Graphene/hBN/Graphene Heterostructures: A Molecular Dynamics Study Hang T. T. Nguyen, Nghia Minh Dong, Dan Lin Lieu, Tin H. Nguyen, and 1 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-8822403/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract The lateral hybridization of graphene is studied using molecular dynamics simulation. The goal is to understand the structural and thermal properties of the in-plane hybrid graphene/hBN/graphene heterostructure, focusing on the dependence of the phase transition on the relative amount of hBN to graphene in the in-plane graphene/hBN/graphene hybrid heterostructure. The Tersoff potential is used to describe covalent bonds between carbon, boron, and nitrogen atoms. Four in-plane graphene/hBN/graphene hybrid configurations were constructed, in which the numbers of atoms in the left graphene, central hBN, and right graphene regions were 32600/25600/32600, 35860/25600/35860, 39120/25600/39120, and 45640/25600/45640, respectively. These configurations were used to examine the effect of the relative graphene-to-hBN ratio on the phase transition behavior. All systems exhibited first-order phase transitions, with broader transition ranges and higher melting points observed as the graphene content increased. Heat capacity peaks shifted from 6120 K to 6650 K, correlating with increased carbon concentration. Coordination number analysis revealed earlier structural breakdown in hBN compared to graphene, reflecting weaker B–N bonding. Angular distribution studies and bond orientational order parameters confirmed that higher graphene content delays structural distortion and enhances thermal resilience. These results demonstrate the critical role of graphene in preserving structural integrity under extreme thermal conditions, offering valuable insights for the design of high-temperature 2D heterostructures. molecular dynamics simulations in-plane Graphene/hBN heterostructure Thermal stability Melting behavior BOO parameter Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 Figure 11 Introduction In recent years, heterostructures have garnered increasing attention from researchers. By definition, the interface of two layers or regions of materials of different compositions is called a heterojunction [ 1 ]. For two-dimensional (2D) materials, heterostructures can be classified into two types: stacked 2D heterostructures and in-plane heterostructures [ 2 ]. The structure consisting of one or many heterojunctions is known as a heterostructure. As this novel structure allows bandgap engineering, heterostructures have become a critical area of research. Graphene, a 2D material with a honeycomb structure, exhibits numerous intriguing characteristics [ 3 ], making it a promising candidate for various applications. Many attempts have been made to incorporate graphene with other compatible materials to create a heterostructure with desired properties. For this reason, other materials with the same honeycomb structure have also been intensely investigated. Silicene and the group-III-nitride family are the most well-known groups of hexagonal-structure materials that are frequently researched. In the group-III-nitride family, the hexagonal Boron Nitride (hBN) is analogous to graphene with a small mismatch (about 1.6%) [ 4 ]. It has the potential to form van der Waals (vdW) interactions with graphene and produce heterostructures with novel characteristics. hBN interactions with living cells were also investigated. It was reported that at a nontoxic concentration, hBN can relieve molecular stress induced by drugs [ 5 ]. The influence of hBN nanoparticles of different doses was evaluated in Wistar albino rats, in which it was shown that only high doses of hBN nanoparticles (1600 and 3200 µg/kg) caused damage in the liver, kidney, heart, spleen, and pancreas, while lower doses (from 50 to 80 µg/kg) were considered as non-toxic [ 6 ]. hBN was also reported to increase dendritic cell maturation without inducing downstream T cell polarization [ 7 ]. In the biomedical field, hBN was suggested to be used in several areas: wound healing [ 8 ], prostate cancer treatment [ 9 ], and anticancer drug carriers [ 10 ]. Due to their relatively small lattice mismatch (about 1.6%) [ 4 ] and the same crystal structure, hBN and graphene can be combined by vertically stacking together. The major three stacking patterns in most studies are: (1) AA stacking, which means C atoms are directly aligned above the B and N atoms; (2) AB stacking, which means one C atom is positioned above the N atom, while the other C atom is positioned above the hBN ring; and (3) AB stacking, meaning that a carbon atom is placed above the B, and the other one is above the hexagonal hBN ring [ 11 ]. The other type of heterostructure of hBN and graphene is the in-plane graphene/hBN heterostructure. The electronic band structures and properties of this lateral graphene/hBN heterostructure were expected to be different from those of pristine graphene and hBN, giving potential for the development of band gap-engineered applications in electronics and optoelectronics [ 12 ]. The in-plane graphene/hBN heterostructure can be used to study the DNA. Bingquan Luan and colleagues introduced an ssDNA stretching platform using a two-dimensional in-plane heterostructure composed of graphene and hBN, demonstrating that ssDNA could be stretched on an hBN nanostripe sandwiched between two adjacent graphene domains ("nanochannel") [ 13 ]. Fabio and his group explored the sensitivity and selectivity of a graphene/hBN-based nanopore architecture for detecting and distinguishing synthetic Hachimoji nucleobases [ 14 ]. Ali Kiakojouri and his team proposed an in-plane graphene/hBN/graphene heterostructure, which exhibited higher sensitivity to different nucleotides compared to a similar structure made of pristine graphene, making it a promising candidate for DNA sequencing applications [ 15 ]. To further investigate how the phase transition behavior of such in-plane graphene/hBN/graphene heterostructures depends on the relative composition of graphene and hBN, we constructed four atomic graphene/hBN/graphene configurations in which the number of hBN atoms was kept constant while the number of graphene atoms was varied such as 32600/25600/32600, 35860/25600/35860, 39120/25600/39120, and 45640/25600/45640 atom arrangements. This approach allows for a detailed assessment of how the graphene content influences the thermal properties and structural stability of the hybrid system. The computational methodology is described in Section 2, followed by the presentation of results in Section 3. Section 4 provides a detailed discussion, and the final section summarizes the main conclusions of the study. Methods and Implementation The role of the interaction potential is critical in molecular dynamics (MD) simulations as it governs how particles behave and move within the simulated system. In MD simulations, the key aim is to accurately replicate real-world interactions between atoms or molecules in a computationally feasible way. The interaction potential encompasses the forces and energies of various atomic interactions, including bonded and non-bonded interactions. In systems with covalent bonds, where chemical bonds connect atoms, the interaction potential dictates the stretching, bending, and torsional forces that occur as atoms undergo movement and rearrangement. This capability allows the simulation to accurately depict processes like bond breaking and formation, which are essential for studying chemical reactions and structural transitions. Various empirical potentials, including those representing two- and three-body interactions in different forms, are available [ 16 – 20 ]. The Tersoff potential specifically addresses the interatomic forces and energies involved in covalent bonding in materials, making them highly suitable for studying systems characterized by strong chemical interactions [ 21 , 22 ]. What distinguishes the Tersoff potential is its ability to model complex phenomena such as bond breaking and reformation, anisotropic and directional bonding, and interactions at both short and long ranges. This broad representation extends their usefulness from basic elements like carbon to a diverse range of covalent materials such as semiconductors and heterostructures. In this study, the Tersoff potential [ 23 ] is used for the interactions between Boron (B), Nitrogen (N), and Carbon (C) in the initial configurations, which are written as below: $$\:{E}_{b}=\frac{1}{2}{\sum\:}_{i\ne\:j}^{}{f}_{c}\left({r}_{ij}\right)\left[{f}_{R}\right({r}_{ij})+{f}_{a}({r}_{ij}\left)\right]$$ 1 . where, \(\:{r}_{ij}\) is the distance from the atom i to atom j. The repulsive \(\:{f}_{R}\left({r}_{ij}\right)\) and the attractive \(\:{f}_{a}\left({r}_{ij}\right)\) terms are based on the Morse potential as proposed by Brenner [ 23 ]. The \(\:{f}_{c}\left({r}_{ij}\right)\) term represents the cutoff function used for calculating the number of neighbors as well as making the potential zero outside the interaction shell. To perform the calculations, the software package Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) is used [ 24 ]. First, to prepare the size of the hBN layer, we based on a previous study on the influence of the armchair/zigzag edge ratio on the melting process of free-standing hBN [ 25 ]. Accordingly, the configuration with an armchair/zigzag edge ratio of 10.745098 and the number of atoms in the model of 25600 is selected for investigation in this study. Then, the hBN configuration will be hybridized with graphene to form an in-plane graphene/hBN/graphene heterostructure (C–B bridge). The interaction at the edge of the graphene and the hBN layer can be seen in Fig. 1 . To study the dependence on the relative amount of hBN to graphene of the melting process of in-plane hybrid graphene/hBN/graphene, four different in-plane graphene/hBN/graphene configurations are considered by changing the number of atoms of the graphene layer (32600/25600/32600, 35860/25600/35860, 39120/25600/39120, and 45640/25600/45640 atom arrangements). After that, the hybrid in-plane graphene/hBN/graphene configurations are relaxed again for \(\:6\times\:1{0}^{5}\) MD steps to ensure that the system reaches mechanical and thermal equilibrium before the heating process begins. Finally, the hybrid graphene/hBN/graphene configuration is heated from 50 K to 8000 K to ensure that the graphene/hBN/graphene configuration is in a liquid state. The change in temperature upon heating obeys \(\:{T}_{0}+vt\) , in which \(\:{T}_{0}=50\) K, \(\:v\) is the heating rate \(\:5\times\:1{0}^{12}\) K/s, and t is the time required for heating. All configurations at the chosen temperatures are relaxed for \(\:6\times\:1{0}^{5}\) MD steps to ensure their stability before studying structural evolution or 2D visualization. The time interval for each step is 0.0001 picoseconds. We use Visual Molecular Dynamics (VMD) [ 26 ] to visualize atomic configurations. Results and Discussions Generally, the atoms of the solid structure are packed tightly together in an orderly arrangement, meaning that atoms vibrate about a fixed point. The fixed position of atoms minimizes the total energy per atom. Upon heating, the total energy increases as kinetic energy is introduced to the system, causing the atoms to vibrate rapidly. When the total energy per atom reaches a critical threshold, the amplitude of these vibrations becomes sufficient to overcome the intermolecular forces operating within the system, initiating a phase transition from the crystalline solid to the liquid state. Therefore, analyzing the total energy per atom provides informative insights into the stable configuration, energy barriers, and phase transitions, such as the solid-to-liquid transition. To investigate how varying the hBN:graphene ratio affects melting in an in-plane graphene/hBN/graphene hybrid, the energy per atom of the four configurations is calculated and shown in Fig. 2 . Overall, the four configurations exhibit a first-order phase transition from the crystalline to the liquid phase (Fig. 2 ). The pattern of first order type includes the three separate regions: i) at low temperatures (below T 1 ), the total energy per atom increases linearly as the atom vibrates around the equilibrium position; ii) when increasing to a critical temperature point (from T 1 to T 3 ), there are sudden jumps in total energy per atom, reflecting the breakdown of the crystalline lattice and the onset of disorder. The disorder arrangement proves a phase transition from the crystalline state to the liquid state of the system; and iii) after the phase transition temperature (from T 3 ), the total energy per atom of the system again increases linearly with temperature, indicating the dominance of translational motion in the liquid state. Although the four configurations exhibit a similar tendency to the first-order type, the phase transition temperature increases with the number of atoms per graphene layer. The 32600/25600/32600 configuration comprises the smallest amount of graphene and has the lowest melting temperature range, from 5850 K to 6320 K, as shown in Fig. 2 (line with crosses). The 35860/25600/35860 configuration (line with hollow circles in Fig. 2 ) and the 39120/25600/39120 configuration (line with crossed diamonds in Fig. 2 ) exhibit higher phase transition temperatures, ranging from 6050 K to 6450 K and from 6240 K to 6550 K, respectively. The highest phase transition temperature observed in the 45640/25600/45640 configuration (line with black squares in Fig. 2 ) is 6790 K. One can see that an upward shift in the phase transition temperature with increasing graphene content in the in-plane hybrid graphene/hBN/graphene systems can be attributed to several points: i) the extended graphene domains, composed of strong sp²-hybridized carbon networks, enhance the rigidity and bonding strength of the structure, making it more resistant to thermal disruptions; ii) in addition, the outer graphene layers function as mechanical anchors, restricting vibrational amplitudes and stabilizing the intermediate hBN region. This structural reinforcement also promotes uniform heat distribution, reducing local fluctuations that might otherwise initiate disorder. As a result, the hybrid system exhibits greater thermal tolerance, requiring higher temperatures to transition from a crystalline to a disordered phase. While the energy per atom provides a general overview of the melting process, analyzing heat capacity identifies the melting points of four configurations of hybrid graphene/hBN/graphene materials. The heat capacity is calculated according to the expression: $$\:{C}_{p}=\:\frac{dE}{dT}$$ 2 During a solid-to-liquid phase transition, atoms gain mobility, and the system undergoes significant reorganization, leading to a sharp change in energy-absorption behavior. This is manifested as a distinct peak in the heat capacity curve. Consequently, the temperature at which this peak occurs can be reliably used to identify the onset of melting. For ease of observation, we will present the heat capacity in Fig. 3 and Fig. 4 . As shown in Fig. 3 and Fig. 4 , the peak in heat capacity corresponds to the disruption of C–C, B–N, and N–B bonds, marking the transition from the crystalline to the liquid phase. The peak heat capacity of four configurations 32600/25600/32600, 35860/25600/35860, 39120/25600/39120, and 45640/25600/45640 is defined at 6120 K, 6270 K, 6390 K, and 6650 K, respectively (T 2 in Fig. 3 and Fig. 4 ). The disparity in temperature points between the four configurations indicates the influence of the number of graphene atoms on the melting process of the hybrid graphene/hBN/graphene configuration. The rising melting temperature associated with increasing graphene size across these configurations suggests a relationship between graphene size and melting behavior. The visualization of four graphene/hBN/graphene configurations was presented in Fig. 5 using VMD. Figure 5 illustrates the atomic structures of four graphene/hBN/graphene configurations at two temperature stages, before and after the melting point. As the temperature increases, significant structural disorder is observed, indicating the breakdown of C–C, B–N, and N–B bonds. The transition becomes more evident with increasing graphene content, which shifts the onset of melting to higher temperatures. This trend highlights the role of graphene size in enhancing the thermal stability of the heterostructures. The structural evolution was visualized using VMD. The coordination number (CN) provides a sensitive probe into the local atomic environment during heating and is therefore highly effective in characterizing the melting process of the graphene/hBN/graphene hybrid heterostructures. In the initial crystalline state at lower temperatures, the graphene and hBN regions form a honeycomb lattice in which each atom is threefold coordinated (CN = 3). Nearly all C, B, and N atoms satisfy this structural requirement, except for a negligible fraction at the interfaces and boundaries. This is consistent with the stability of the sp² bonding framework in both graphene and hBN. Upon heating, thermal vibrations increase progressively, leading to bond stretching and breaking. The destruction of C–C, B–N, and N–B bonds reduces the number of atoms maintaining CN = 3, accompanied by the appearance of atoms with lower coordination number (CN = 2, 1, or 0). At a certain temperature, the CN plot exhibits sudden collapse, corresponding to the melted point of the material. By monitoring the evolution of coordination numbers for B, N, and C atoms as functions of temperature, the phase transition can be identified and correlated with other thermodynamic indicators such as the total energy per atom and the heat capacity. The transition temperatures and melting points vary systematically with the graphene:hBN ratio. As the number of graphene atoms increases, the critical temperatures shift to higher values, showing that graphene layers act as stabilizing components to the heterostructure. In the 32600/25600/32600 configuration, at early temperature points (3750 K to 4750 K), the CN plots of B, N, and C atoms decrease at steady rate (Fig. 6 a). At 3750 K, the CN plots of B and N are both approximately 50% (line with black squares for B, line with red circles for N in Fig. 6 a), much higher than that of C (less than 25% at this temperature point, line with blue triangles for C in Fig. 6 a). At about 4750 K, the CN plots of B and N collapse sharply (line with black squares for B, line with red circles for N in Fig. 6 a), marking the initial disordering of the hBN sublattice. The collapse of C occurs between 6100–6250 K (line with blue triangles for C in Fig. 6 a), overlapping with the identified melting point at about 6120 K in the heat capacity curve (Fig. 3 a). The fact that the CN plot of C degrades later than that of B and N reflects the delayed disordering of graphene relative to hBN. At higher temperature points (from 6500 K), both B and N atoms reach near-zero coordination (line with black squares for B, line with red circles for N in Fig. 6 a), while C atoms retain a small fraction, but observably greater than 0, of CN = 3 (line with blue triangles for C in Fig. 6 a). After the steep decline of C, the CN plot of C is consistently higher than that of B and N, which is opposite to the observed trend at lower temperatures. In the 35860/25600/35860 configuration, the CN plots of B, N, and C have a similar appearance to the CN plots of the 32600/25600/32600 configuration. However, at 3750 K, while the percentage of CN = 3 of C is roughly 25% (line with blue triangles for C in Fig. 6 b), which is about the same that of C in the 32600/25600/32600 configuration, the percentage of CN = 3 of B and N is only less than 45% (line with black squares for B, line with red circles for N in Fig. 6 b). In the 4750 K − 5000 K range, the CN plots of B and N drop dramatically from 30% to nearly 0% (line with black squares for B, red circle for N in Fig. 6 b). The collapse of C’s CN plot happens at a higher temperature, about 6100 K (line with blue triangles for C in Fig. 6 b). In the 39120/25600/39120 configuration, at 3750 K, the CN plots of B and N are around 25–30% (line with black squares for B, red circle for N in Fig. 6 c). Meanwhile, the plot of C is in the vicinity of 20% (line with blue triangles for C in Fig. 6 c). Like the plots of other discussed configurations, at lower temperature points, the plots of B and N are observably higher than that of C. The C atoms in the 39120/25600/39120 configuration also show stronger resilience against heat compared to B and N atoms: the plot of C degrades at higher temperature, about 6250 K (line with blue triangles for C in Fig. 6 c), while the plots of B and N deteriorate at roughly 4750 K (line with black squares for B, line with red circles for N in Fig. 6 c). At temperatures higher than 6500 K, the plot of C continuously decreases as the temperature increases, but still remains well above 0% (line with blue triangles for C in Fig. 6 c); whereas the plots of B and N are close to 0% at these temperature points (line with black squares for B, line with red circles for N in Fig. 6 c). In the 45640/25600/45640 configuration, the CN plots of B, N, and C demonstrate the most delayed transition among all examined configurations. At 3750 K, the CN plot of C does not indicate great differences from the plots of B and N like other configurations (all plots are approximately 8–12% at this temperature point; line with blue triangles for C, black squares for B, and red circles for N in Fig. 6 d). The CN plots of B and N begin to collapse above 6300 K, while the CN plot of C plummets at nearly 6700 K (black square for B, red circle for N in Fig. 6 d). The persistence of C coordination underscores the dominant role of graphene in stabilizing the hybrid heterostructure. Thus, larger graphene layers act as a thermal backbone, resisting complete collapse and ensuring superior stability. At temperatures higher than 6500 K, while the plots of B and N are approximately 0% (line with black squares for B, line with red circles for N in Fig. 6 d), the plot of C remains at about 3–4% (line with blue triangles for C in Fig. 6 d), constantly higher than that of B and N. A closer comparison among atomic species reveals that B and N atoms lose coordination earlier and more rapidly at lower temperatures than those of carbon. This behavior can be attributed to the relative weakness of B–N bonds compared to C–C bonds. As temperature rises, B–N bonds break first, destabilizing the hBN regions and initiating local disorder. Graphene, with its stronger C–C bonds, resists this disorder longer, maintaining partial threefold coordination and thereby delaying the global phase transition. The result is a progressive shift of the melting point with increasing graphene content. The coordination number analysis is also consistent with energy and heat capacity calculations. In summary, the coordination number analysis demonstrates that melting in graphene/hBN/graphene hybrid heterostructures proceeds via a first-order phase transition. The critical transition temperatures and melting points systematically increase with the graphene:hBN ratio, highlighting the stabilizing role of graphene. The preferential survival of C–C coordination to higher temperatures further confirms graphene as the primary contributor to thermal resilience. These results underline the importance of atomic composition in dictating the structural integrity of 2D in-plane hybrid graphene/hBN/graphene heterostructures under extreme thermal conditions. To analyze structural deformation and clarify the stabilizing role of graphene layers in planar graphene/hBN/graphene heterostructures during heating, the angular distribution of C–C–C and B–N–B was evaluated for four configurations (32600/25600/32600, 35860/25600/35860, 39120/25600/39120, and 45640/25600/45640). The analysis was performed at two characteristic temperatures: the melting start temperature, determined from the initial abrupt change in total energy per atom (T 1 in Fig. 2 ), and the temperature point just above the melting point (T 3 in Fig. 2 ). For ease of reference, the results are presented in Fig. 7 to Fig. 10 . First, for the configuration with the lowest carbon content (the 32600/25600/32600 configuration), at the pre-melting temperature (5850 K), the angular distribution of both C–C–C and B–N–B remains narrow, exhibiting distinct peaks around 120° (approximately 0.8% for graphene and 0.9% for hBN). This demonstrates that the hexagonal lattice is largely preserved in both regions (Fig. 7 a) because at this temperature, thermal vibrations are insufficient to significantly disrupt the sp² linkage network. As the temperature increases to 6320 K, the angular distribution of C–C–C expands significantly (spanning approximately 85°–155°) while retaining a distinct peak at 120° with a decrease in intensity (≈ 0.6%) (line with red circles for C–C–C in Fig. 7 b). The persistence of the peak at 120° suggests that the short-range hexagonal order in graphene is not completely destroyed at the melting point. Meanwhile, the B–N–B angular distribution becomes significantly less concentrated around 120°. Specifically, the peak intensity decreases to approximately 0.7% and the shape is noticeably distorted (line with black squares for B–N–B in Fig. 7 b). This marked expansion reflects a more severe angular disorder in the hBN sublattice. When increasing the number of atoms in the graphene layer to 35860 atoms per layer (35860/25600/35860), the melting start temperature increased, implying that the stability of this configuration is significantly improved. Similar to the previous configuration, we will investigate the pre-melting temperature. At 6050 K, the peak at 120° of the C–C–C angle distribution is narrow and sharp (peak at approximately 0.9%), whereas the B–N–B angle distribution is flatter with more ripples at the top (peaks of these ripples are approximately 0.6%; line with red circles for C–C–C and line with black squares for B–N–B in Fig. 8 a). The B–N–B angle distribution overlaps significantly. Upon reaching 6450 K, the C–C–C peak at 120° is broadened, but the peak remains visible (peak reduction to 0.7%, line with red circles for C–C–C in Fig. 8 b), while the B–N–B begins to significantly deform (peak reduction to 0.5%, line with black squares in Fig. 8 b), indicating a higher degree of structural disorder. Compared to the 32600/25600/32600 configuration, due to the higher amount of C in this configuration, the loss of integrity of the hBN layer occurs more slowly. Continuing to increase the carbon content in the graphene layer to 39120 atoms per layer (39120/25600/39120), the melting temperature begins to rise higher than the two previous configurations. At the pre-melting temperature of 6240 K, the angular distribution of C–C–C shows a clear and concentrated peak (peak at approximately 0.8%, line with red circles for C–C–C in Fig. 9 a). The angular distribution of B–N–B is more dispersed, with the distribution at 120° accounts for about 0.5% and the highest peak is about 0.6% (line with black squares for B–N–B in Fig. 9 a). Continuing to increase the temperature to 6550 K, the angular distribution of C–C–C becomes slightly more dispersed (peak decreases to 0.7%, line with red circles for C–C–C in Fig. 9 b). The B–N–B angle distribution is flattened (line with black squares for B–N–B in Fig. 9 b) but is significantly less affected by temperature than the previous 32600/25600/32600 and 35860/25600/35860 configurations. This implies that the abundance of C molecules in the graphene layer better inhibits hBN structural deformation. Finally, when the number of atoms in the graphene layer is 45640 atoms, the best heat resistance is observed. Specifically, the pre-melting temperature is 6430 K. At the initial melting temperature of 6430 K, the angular distribution of C–C–C is sharp and narrow (the angular distribution center is approximately 120°, with the peak at over 0.8% (line with red circles for C–C–C in Fig. 10 a). On the other hand, B–N–B shows low structural integrity (line with black squares for B–N–B in Fig. 10 a). When the temperature increases to 6790 K (melting point), the C–C–C angular distribution is affected but still maintains a peak at 120° (peak at approximately 0.7%, line with red circles for C–C–C in Fig. 10 b). The B–N–B angle distribution is wide, but its variation is the least compared to other configurations (line with black squares for B–N–B in Fig. 10 b). This configuration has the highest melting temperature, suggesting that the amount of C in the graphene layer has the strongest impact on hBN stabilization. One can see that graphene contributes significantly to the stabilization of the hBN layer. In all configurations, the change in the B–N–B angle distribution is the most sensitive sign reflecting the structural degradation of the hBN layer. As the carbon content increases, this change is delayed, and the intensity is reduced, indicating the protective effect of graphene. At the same time, the C–C–C angle distribution still maintains a clear peak at high temperatures, demonstrating that the graphene network has better thermal resistance than hBN. It can be concluded that the angle distribution analysis provides microscopic evidence that graphene not only stabilizes itself but also plays an important role in maintaining the structure of the neighboring hBN layer. Therefore, adjusting the carbon content is an effective strategy to control the thermal stability of 2D heterojunction material systems, aiming at applications in electronic devices or materials operating in high-temperature environments. To characterize the structural evolution of the hybrid graphene/hBN/graphene configuration under increasing temperature, this study employs the Bond-Orientational Order (BOO) parameter, which quantifies the degree of local orientational symmetry in solid-state systems. BOO provides a measure of the bond alignment around each particle, thereby offering insights into phase transitions and structural disorder that may not be fully captured by simple radial distribution functions. In particular, for hexagonal lattices commonly found in two-dimensional materials such as graphene and hBN, the sixth-order BOO parameter (Q₆) is especially relevant, as it reflects the inherent hexagonal symmetry of these systems. The value of Q₆ ranges from 0 to 1, with crystalline phases exhibiting values close to unity, while reductions in magnitude indicate enhanced thermal agitation, amorphization, or melting transitions. In our simulations, monitoring Q₆ as a function of temperature reveals the onset of structural instabilities, providing quantitative evidence of how elevated thermal energy alters the orientational order of bonds within the hybrid graphene/hBN/graphene configuration. Overall, across all four configurations (Fig. 11 ), the BOO values associated with B, N, and C atoms at the initial temperature of 50 K are close to unity, confirming that the configurations possess nearly perfect orientational symmetry consistent with a hexagonal lattice. However, the BOO parameter subsequently decreases with rising temperature. This trend is consistent with the physical expectation that thermal fluctuations intensify atomic vibrations, thereby displacing atoms from their original positions. Nevertheless, the degree of reduction is not identical across configurations, nor is it uniform among different atomic species within the same configuration. For each configuration, the variations can be described as below. For the 32600/25600/32600 configuration (Fig. 11 a), at 50 K, the BOO values of B and N are nearly identical, both around 0.99, while the corresponding value for C is slightly lower at 0.96. Overall, both the hBN and graphene layers exhibit a high degree of hexagonal orientational symmetry at the initial temperature (50 K). With respect to B and N, their evolution during heating follows almost the same trend: from 50 K to 3000 K, the BOO decreases minimally (reaching 0.96 for both elements at 3258 K); beyond 3500 K, the rate of decrease gradually accelerates until around 6472 K, where the drop becomes most significant with the BOO falling to 0.63, indicating a rapid structural transformation of the hBN layer at this temperature. Thereafter, the decline rate slows down, making the plot reaching a final value of 0.36 at 8000 K. In contrast, C shows a distinct three-stage behavior: in the first stage, from 50 K to 5600 K, the BOO decreases linearly from 0.96 to 0.66; in the second stage, between 5600 K and 6300 K, the parameter drops rapidly to 0.34, with the sharpest change observed at 6110 K, reflecting a clear distortion in the original hexagonal structure of the graphene layer, characteristic of the material’s phase transition. Finally, in the third stage, the curve continues to decrease linearly, reaching as low as 0.20. For the 35860/25600/35860 configuration (Fig. 11 b), at 50 K, the BOO values for B and N are 0.99, while that for C is 0.96. Similar to the first configuration, the BOO of B and N decreases very slowly from 50 K to 3000 K, after which the rate of decline gradually accelerates until 6343 K, where it reaches 0.55, indicating a rapid structural transformation characteristic of the phase transition. The BOO then decreases more slowly to 0.29 at 8000 K. For C, the reduction occurs in three distinct stages: in the first stage, from 50 K to 5700 K, the BOO decreases linearly from 0.96 to 0.66; in the second stage, from 5700 K to 6300 K, the curve drops sharply to 0.38, confirming the occurrence of a structural phase transition most evident at 6250 K; and in the final stage, it continues to decline linearly, reaching 0.22 at the end of the simulation. In the 39120/25600/39120 configuration (Fig. 11 c), the BOO values for both elements in the hBN layer are 0.99 at the initial temperature, then slightly decrease to 0.96 at 2500 K. From 2500 K to 6443 K, the rate of decrease accelerates, with the BOO reaching 0.45 at the end of this interval, representing the phase transition stage in the hBN layer and subsequently slowing down to a final value of 0.26 at 8000 K. For C, the decline follows a pattern similar to that observed in the first two configurations and can be divided into three distinct stages: stage 1, from 50 K to 6100 K, the BOO decreases linearly from 0.96 to 0.62; stage 2, from 6100 K to 6500 K, the parameter drops rapidly, reaching 0.35 by the end of this stage, indicating the temperature range where the phase transition occurs as the hexagonal structure becomes distorted and loses its initial symmetry, with the most pronounced jump observed at 6342 K; stage 3, from 6500 K to 8000 K, the decline returns to a linear trend, resulting in a final value of 0.22. However, unlike the two previous configurations, where the BOO of B and N remained consistently higher than that of C throughout heating, the third configuration exhibits a reversed trend between 4900 K and 6400 K, during which the BOO of C surpasses that of the other two elements. This behavior arises because the structural degradation in the hBN layer initiates earlier at lower temperatures, whereas the rapid decline in C occurs at higher temperatures compared to the preceding configurations. In the 45640/25600/45640 configuration (Fig. 11 d), the BOO curves of B and N display earlier structural changes compared to the previous configurations. Initially, the BOO is 0.99, then the rate of decrease accelerates rapidly up to around 6000 K, at which point the BOO drops to 0.62, representing a rapid structural deformation characteristic of this phase-transition stage. Thereafter, the decrease slows progressively, reaching 0.24 at 8000 K. The BOO of C continues to follow a three-stage pattern: stage 1, from 50 K to 6200 K, the BOO decreases linearly from 0.96 to 0.65; stage 2, from 6200 K to 6800 K, where a sharp drop is observed particularly at 6648 K to 0.36 at the end of this stage, confirming the phase transition in the graphene layer; stage 3, beyond this point, the curve resumes a linear decrease, reaching a final value of 0.24 at 8000 K. In this configuration, a clear distinction can be observed: for most of the simulated temperature range (from 3800 K to 8000 K), the BOO of C remains higher than that of the other two elements. In general comparison, the transformation stages across all four configurations exhibit similar trends. However, as the number of atoms in the graphene layer increases, the BOO variation in the hBN layer occurs at lower temperatures, whereas the BOO variation in the graphene layer shifts to higher temperatures. This leads to noticeable changes in the BOO-temperature distributions among the configurations: in the first configuration, the BOO of the hBN layer consistently remains higher than that of graphene, resulting in two clearly separated curves; in the second configuration, the two curves gradually converge; and in the last two configurations, a distinct temperature interval appears where the BOO of graphene exceeds that of hBN. This trend is most evident in the progressive increase of the temperature at which the sharp BOO drop in graphene occurs, from 6110 K in configuration 32600/25600/32600, to 6250 K in configuration 35860/25600/35860, and 6342 K and 6648 K in configurations 39120/25600/39120 and 45640/25600/45640, respectively. These values are consistent with the preceding peak positions of the Cp curves (Fig. 3 and Fig. 4 ), which are predicted to correspond to the phase transition points of the configurations at approximately 6120 K, 6270 K, 6390 K, and 6650 K, respectively. Although the hBN and graphene layers exhibit opposite trends as the number of Carbon atoms increases, the average BOO values of the two layers across all four configurations remain within the range of 0.20–0.35. This indicates that a significant portion of the initial hexagonal orientational structure is still preserved within the configuration and has not been destroyed even at the final simulation temperature of 8000 K. Conclusions The melting process of 32600/25600/32600, 35860/25600/35860, 39120/25600/39120, and 45640/25600/45640 configurations of the in-plane hybrid graphene/hBN/graphene heterostructure is investigated using MD simulations. The main results are: An analysis of the energy per atom reveals that the four hybrid configurations of graphene/hBN/graphene exhibit a first-order phase transition. Configurations with more graphene atoms exhibit a phase transition at higher temperature ranges, indicating greater thermal tolerance. Additionally, examining the heat capacity of the hybrid system provides insights into how the melting point varies with the number of graphene atoms. The results indicate that as the number of graphene atoms increases, the peak heat capacity also rises, with peaks observed at 6120 K, 6270 K, 6390 K, and 6650 K. This finding highlights the relationship between graphene size and the melting behavior of the hybrid configurations. The coordination number analysis reveals a clear first-order melting transition in the graphene/hBN/graphene hybrid heterostructures under heating. All four configurations exhibit a gradual reduction in threefold atomic coordination at lower temperatures, followed by an abrupt collapse at higher temperatures corresponding to the melting point. The B and N atoms in hBN lose their coordination earlier than the C atoms, confirming weaker B–N bonding compared to C–C bonding. Increasing the graphene fraction systematically shifts the transition to higher temperatures, demonstrating that graphene enhances structural stability under extreme thermal conditions. The persistence of CN = 3 in carbon atoms to higher temperatures highlights the dominant role of graphene as a thermally resilient backbone. These results are consistent with heat capacity behavior and energy trends, emphasizing that atomic composition strongly governs the thermal resistance of graphene/hBN hybrid systems. Through observations of the angular distributions at the melting onset temperature and the temperature point just above the melting point of four configurations, one can conclude that increasing the carbon content in the graphene layer helps to improve the thermal stability of the graphene/hBN/graphene heterostructure. As the carbon content increases, the change in the B–N–B angle distribution is delayed and the peak value of the distribution is less changed, indicating the protective effect of graphene. To comprehensively evaluate the orientational order and the preservation of the hexagonal framework during the heating process of both hBN and graphene layers within the hybrid material, this study employs the BOO parameter. Across all four configurations, the BOO values decrease with increasing temperature due to enhanced atomic vibrations that gradually distort the ideal hexagonal arrangement. As the number of C atoms in the graphene layer increases, an opposite trend arises between the two layers: the hBN layer exhibits a faster BOO reduction, indicating earlier structural transformations, whereas the graphene layer shows a slower variation with higher carbon concentration. This contrasting behavior strongly affects the predicted phase transition temperature, as the abrupt BOO drop in graphene coincides closely with the Cp peaks identified previously. Declarations Author Contribution Hang T. T. Nguyen conceived of the presented idea, performed the computations, analyzed the data, supervised the findings of this work, and wrote the manuscript. Dan Lin Lieu analyzed the data, wrote the manuscript, and submitted the manuscript. Nghia Minh Dong analyzed the data and wrote the manuscript. Yen Khanh Nguyen analyzed the data and wrote the manuscript. Tin H. Nguyen analyzed the data and wrote the manuscript. Acknowledgement We acknowledge Ho Chi Minh City University of Technology (HCMUT), VNU-HCM for supporting this study. Data Availability Data sets generated during the current study are available from the corresponding author on reasonable request. References N. F. Andrade Neto, A. B. Lima, R. R. Y. O. V. Wilson, T. C. N. Nicacio, M. R. D. Bomio, and F. V. Motta, “Heterostructures obtained by ultrasonic methods for photocatalytic application: A review,” Materials Science in Semiconductor Processing , vol. 139, p. 106311, Mar. 2022, doi: 10.1016/j.mssp.2021.106311 . H. Wang, F. Liu, W. Fu, Z. Fang, W. Zhou, and Z. Liu, “Two-dimensional heterostructures: fabrication, characterization, and application,” Nanoscale , vol. 6, no. 21, pp. 12250–12272, Aug. 2014, doi: 10.1039/C4NR03435J . N. Sharma, R. Dev Gupta, R. Chandmal Sharma, S. Dayal, and A. Singh Yadav, “Graphene: An overview of its characteristics and applications,” Materials Today: Proceedings , vol. 47, pp. 2752–2755, 2021, doi: 10.1016/j.matpr.2021.03.086 . L. H. Li and Y. Chen, “Atomically Thin Boron Nitride: Unique Properties and Applications,” Adv Funct Materials , vol. 26, no. 16, pp. 2594–2608, Apr. 2016, doi: 10.1002/adfm.201504606 . I. C. Taskin, O. Sen, M. Emanet, M. Culha, and B. Yilmaz, “Hexagonal boron nitrides reduce the oxidative stress on cells,” Nanotechnology , vol. 31, no. 21, p. 215101, May 2020, doi: 10.1088/1361-6528/ab6fdc . F. Kar et al., “In Vivo Assessment of the Effect of Hexagonal Boron Nitride Nanoparticles on Biochemical, Histopathological, Oxidant and Antioxidant Status,” J Clust Sci , vol. 32, no. 2, pp. 517–529, Mar. 2021, doi: 10.1007/s10876-020-01811-w . H. Lin et al., “2D Materials and Primary Human Dendritic Cells: A Comparative Cytotoxicity Study,” Small , vol. 18, no. 20, p. 2107652, May 2022, doi: 10.1002/smll.202107652 . Ö. Şen, M. Emanet, and M. Çulha, “Stimulatory Effect of Hexagonal Boron Nitrides in Wound Healing,” ACS Appl. Bio Mater. , vol. 2, no. 12, pp. 5582–5596, Dec. 2019, doi: 10.1021/acsabm.9b00669 . M. Emanet Ciofani, Ö. Şen, and M. Çulha, “Hexagonal Boron Nitride Nanoparticles for Prostate Cancer Treatment,” ACS Appl. Nano Mater. , vol. 3, no. 3, pp. 2364–2372, Mar. 2020, doi: 10.1021/acsanm.9b02486 . S. M. Sharker, “Hexagonal Boron Nitrides (White Graphene): A Promising Method for Cancer Drug Delivery,” Int J Nanomedicine , vol. 14, pp. 9983–9993, 2019, doi: 10.2147/IJN.S205095 . J. Wang, F. Ma, and M. Sun, “Graphene, hexagonal boron nitride, and their heterostructures: properties and applications,” RSC Adv. , vol. 7, no. 27, pp. 16801–16822, 2017, doi: 10.1039/C7RA00260B . H. Wang, F. Liu, W. Fu, Z. Fang, W. Zhou, and Z. Liu, “Two-dimensional heterostructures: fabrication, characterization, and application,” Nanoscale , vol. 6, no. 21, pp. 12250–12272, Aug. 2014, doi: 10.1039/C4NR03435J . B. Luan and R. Zhou, “Spontaneous ssDNA stretching on graphene and hexagonal boron nitride in plane heterostructures,” Nat Commun , vol. 10, no. 1, p. 4610, Oct. 2019, doi: 10.1038/s41467-019-12584-w . A. Kiakojouri, I. Frank, and E. Nadimi, “In-plane graphene/h-BN/graphene heterostructures with nanopores for electrical detection of DNA nucleotides,” Phys. Chem. Chem. Phys. , vol. 23, no. 44, pp. 25126–25135, 2021, doi: 10.1039/D1CP03597E . F. A. L. De Souza, G. Sivaraman, M. Fyta, R. H. Scheicher, W. L. Scopel, and R. G. Amorim, “Electrically sensing Hachimoji DNA nucleotides through a hybrid graphene/ h -BN nanopore,” Nanoscale , vol. 12, no. 35, pp. 18289–18295, 2020, doi: 10.1039/D0NR04363J . G. C. Abell, “Empirical chemical pseudopotential theory of molecular and metallic bonding,” Phys. Rev. B , vol. 31, no. 10, pp. 6184–6196, May 1985, doi: 10.1103/PhysRevB.31.6184 . D. W. Brenner, “Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films,” Phys. Rev. B , vol. 42, no. 15, pp. 9458–9471, Nov. 1990, doi: 10.1103/PhysRevB.42.9458 . M. Büttiker, Y. Imry, R. Landauer, and S. Pinhas, “Generalized many-channel conductance formula with application to small rings,” Phys. Rev. B , vol. 31, no. 10, pp. 6207–6215, May 1985, doi: 10.1103/PhysRevB.31.6207 . J. Tersoff, “New empirical approach for the structure and energy of covalent systems,” Phys. Rev. B , vol. 37, no. 12, pp. 6991–7000, Apr. 1988, doi: 10.1103/PhysRevB.37.6991 . F. H. Stillinger and T. A. Weber, “Computer simulation of local order in condensed phases of silicon,” Phys. Rev. B , vol. 31, no. 8, pp. 5262–5271, Apr. 1985, doi: 10.1103/PhysRevB.31.5262 . J. Tersoff, “Modeling solid-state chemistry: Interatomic potentials for multicomponent systems,” Phys. Rev. B , vol. 39, no. 8, pp. 5566–5568, Mar. 1989, doi: 10.1103/PhysRevB.39.5566 . A. Kınacı, J. B. Haskins, C. Sevik, and T. Çağın, “Thermal conductivity of BN-C nanostructures,” Phys. Rev. B , vol. 86, no. 11, p. 115410, Sept. 2012, doi: 10.1103/PhysRevB.86.115410 . D. W. Brenner, “Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films,” Phys. Rev. B , vol. 42, no. 15, pp. 9458–9471, Nov. 1990, doi: 10.1103/PhysRevB.42.9458 . S. Plimpton, “Fast Parallel Algorithms for Short-Range Molecular Dynamics,” Journal of Computational Physics , vol. 117, no. 1, pp. 1–19, Mar. 1995, doi: 10.1006/jcph.1995.1039 . H. T. T. Nguyen, “Edge effects on the melting process of two-dimensional hexagonal boron nitride,” J Nanopart Res , vol. 26, no. 8, p. 199, Aug. 2024, doi: 10.1007/s11051-024-06108-x . W. Humphrey, A. Dalke, and K. Schulten, “VMD: Visual molecular dynamics,” Journal of Molecular Graphics , vol. 14, no. 1, pp. 33–38, Feb. 1996, doi: 10.1016/0263-7855(96)00018-5 . Additional Declarations No competing interests reported. Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-8822403","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":600269662,"identity":"66e81848-bec5-43e1-86bc-8ff7a99b9ada","order_by":0,"name":"Hang T. T. Nguyen","email":"","orcid":"","institution":"Ho Chi Minh City University of Technology (HCMUT)","correspondingAuthor":false,"prefix":"","firstName":"Hang","middleName":"T. T.","lastName":"Nguyen","suffix":""},{"id":600269664,"identity":"002eba15-4b16-464c-b68f-acbeb4b34b1c","order_by":1,"name":"Nghia Minh Dong","email":"","orcid":"","institution":"Ho Chi Minh City University of Technology (HCMUT)","correspondingAuthor":false,"prefix":"","firstName":"Nghia","middleName":"Minh","lastName":"Dong","suffix":""},{"id":600269667,"identity":"2d8a69fd-49db-4603-9268-daec28eab6c9","order_by":2,"name":"Dan Lin Lieu","email":"","orcid":"","institution":"Ho Chi Minh City University of Technology (HCMUT)","correspondingAuthor":false,"prefix":"","firstName":"Dan","middleName":"Lin","lastName":"Lieu","suffix":""},{"id":600269669,"identity":"ec0079e9-9d8f-4f93-9b8e-bbec65d00c3c","order_by":3,"name":"Tin H. Nguyen","email":"","orcid":"","institution":"Ho Chi Minh City University of Technology (HCMUT)","correspondingAuthor":false,"prefix":"","firstName":"Tin","middleName":"H.","lastName":"Nguyen","suffix":""},{"id":600269670,"identity":"07c714ba-9b5a-4d3a-a39c-7ca7d6f88ea3","order_by":4,"name":"Yen Khanh Nguyen","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABIElEQVRIiWNgGAWjYHADNjYgcUDOAMplbCBWizFYywFStCRuIKTF4HjvAeaCX4cT+6Xb0h783HEnfTv76TTpDww2shsO8Jg9wKblzLkE5pl9hxNnzjl23LD3zLPcnT252yQOMKQZA7WYG2DTciPHgJm353DihhvpbRK8bYdzN9zgBWkBihxgS5PApuX+G4iW/UAtkn/bDqcbQLT8x63lBo8BM88PoJkSacekgbYkQLUAw+EA8zFsWiTP5Bgc5m1IN55xIy1NWrbtsOGGM7mbLc4YJBvPPIxdC9/xM4aPef5Yy/bPSDOTfNt2WN7g+NmNNyoq7GT7jje2YdOicAAYCYxtDI4NaA4GYmYs6oFAHqz0D4M9dulRMApGwSgYBUAAAFd6c8gamo4vAAAAAElFTkSuQmCC","orcid":"","institution":"Ho Chi Minh City University of Technology (HCMUT)","correspondingAuthor":true,"prefix":"","firstName":"Yen","middleName":"Khanh","lastName":"Nguyen","suffix":""}],"badges":[],"createdAt":"2026-02-08 14:53:14","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-8822403/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-8822403/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":104020964,"identity":"8747dd3e-5431-4b3c-a5e1-ae7c7fb3711b","added_by":"auto","created_at":"2026-03-05 18:34:31","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":66969,"visible":true,"origin":"","legend":"\u003cp\u003eThe initial in-plane graphene/hBN/graphene heterostructure (C–B bridge). B, N, and C atoms are colored in ochre, blue, and cyan, in turn\u003c/p\u003e","description":"","filename":"Fig.1.jpg","url":"https://assets-eu.researchsquare.com/files/rs-8822403/v1/ea57d7a24c5dd75b408cc9a3.jpg"},{"id":104403013,"identity":"b7a216a6-c584-4b0d-9a05-ab3ebcc51035","added_by":"auto","created_at":"2026-03-11 12:17:10","extension":"jpeg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":421967,"visible":true,"origin":"","legend":"\u003cp\u003eThe energy per atom of the four configurations: a) The 32600/25600/32600 configuration (line with crosses), b) The 35860/25600/35860 configuration (line with hollow circles), c) The 39120/25600/39120 configuration (line with crossed diamonds), and d) The 45640/25600/45640 configuration (line with black squares)\u003c/p\u003e","description":"","filename":"Fig.2.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-8822403/v1/41c891a49e8f13ca0e842b93.jpeg"},{"id":104020966,"identity":"d4cf94d4-862e-4e40-98f6-7dc24ae37e32","added_by":"auto","created_at":"2026-03-05 18:34:31","extension":"jpeg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":653971,"visible":true,"origin":"","legend":"\u003cp\u003eEnergy per atom and heat capacity for two configurations: a) The 32600/25600/32600 configuration (line with squares), b) The 35860/25600/35860 configuration (line with circles)\u003c/p\u003e","description":"","filename":"Fig.3.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-8822403/v1/e536aa786690f3f6cfd151af.jpeg"},{"id":104020967,"identity":"f651cf6c-e8cb-4215-a1e1-7c18b6ef1c5d","added_by":"auto","created_at":"2026-03-05 18:34:31","extension":"jpeg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":659299,"visible":true,"origin":"","legend":"\u003cp\u003eEnergy per atom and heat capacity for two configurations: c) The 39120/25600/39120 configuration (line with up-pointing triangles), and d) The 45640/25600/45640 configuration (line with down-pointing triangles)\u003c/p\u003e","description":"","filename":"Fig.4.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-8822403/v1/64e0404965c15f0e4e7ed1ba.jpeg"},{"id":104403009,"identity":"508e735a-575d-4ce2-8d4b-d1037ed3ed84","added_by":"auto","created_at":"2026-03-11 12:17:10","extension":"jpeg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":123789,"visible":true,"origin":"","legend":"\u003cp\u003eThe visualization of the four graphene/hBN/graphene configurations at the beginning and at the end of the melting process: a) The 32600/25600/32600 configuration, b) The 35860/25600/35860 configuration, c) The 39120/25600/39120 configuration, and d) The 45640/25600/45640 configuration\u003c/p\u003e","description":"","filename":"Fig.5.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-8822403/v1/520d731c6d2bef2ab5c6b966.jpeg"},{"id":104020975,"identity":"001be5ef-62ee-48d8-9820-397569ffb953","added_by":"auto","created_at":"2026-03-05 18:34:32","extension":"jpeg","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":43756,"visible":true,"origin":"","legend":"\u003cp\u003eThe coordination number of 3 for C (line with blue triangles), B (line with red circles), and N (line with black squares) in the in-plane hybrid graphene/hBN/graphene heterostructure upon heating: a) The 32600/25600/32600 configuration, b) The 35860/25600/35860 configuration, c) The 39120/25600/39120 configuration, and d) The 45640/25600/45640 configuration\u003c/p\u003e","description":"","filename":"Fig.6.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-8822403/v1/21ceed00544757d3f0352519.jpeg"},{"id":104402661,"identity":"6ccc6e2e-9e49-4f09-bca3-373a4599fadb","added_by":"auto","created_at":"2026-03-11 12:16:01","extension":"jpeg","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":577041,"visible":true,"origin":"","legend":"\u003cp\u003eAngular distribution of the 32600/25600/32600 configuration of B–N–B (line with black squares) and C–C–C (line with red circles) at a) 5850 K and b) 6320 K\u003c/p\u003e","description":"","filename":"Fig.7.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-8822403/v1/231bf6fe6dc3b3a8b525d44e.jpeg"},{"id":104020968,"identity":"8998ddc1-c036-477d-a9ee-6fe2a380f12c","added_by":"auto","created_at":"2026-03-05 18:34:31","extension":"jpeg","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":563186,"visible":true,"origin":"","legend":"\u003cp\u003eAngular distribution of the 35860/25600/35860 configuration of B–N–B (line with black squares) and C–C–C (line with red circles) at a) 6050 K and b) 6450 K\u003c/p\u003e","description":"","filename":"Fig.8.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-8822403/v1/b2c4ce9cb997d4574abec7d9.jpeg"},{"id":104403299,"identity":"171c61f5-8842-402d-a81b-96f5ffe27f3c","added_by":"auto","created_at":"2026-03-11 12:17:58","extension":"jpeg","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":560945,"visible":true,"origin":"","legend":"\u003cp\u003eAngular distribution of the 39120/25600/39120 configuration of B–N–B (line with black squares) and C–C–C (line with red circles) at a) 6240 K and b) 6550 K\u003c/p\u003e","description":"","filename":"Fig.9.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-8822403/v1/f46004be9039df712b65a184.jpeg"},{"id":104020971,"identity":"634ed81b-879b-4efc-8e36-c9cc2ad075e4","added_by":"auto","created_at":"2026-03-05 18:34:32","extension":"jpeg","order_by":10,"title":"Figure 10","display":"","copyAsset":false,"role":"figure","size":571937,"visible":true,"origin":"","legend":"\u003cp\u003eAngular distribution of the 45640/25600/45640 configuration of B–N–B (line with black squares) and C–C–C (line with red circles) at a) 6430 K and b) 6790 K\u003c/p\u003e","description":"","filename":"Fig.10.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-8822403/v1/1b76ed8a44c77468cd619ffa.jpeg"},{"id":104020973,"identity":"dadff546-1bb1-499b-ae4d-be36428d63d3","added_by":"auto","created_at":"2026-03-05 18:34:32","extension":"jpeg","order_by":11,"title":"Figure 11","display":"","copyAsset":false,"role":"figure","size":1421128,"visible":true,"origin":"","legend":"\u003cp\u003eThe graph of the average BOO values for B (line with black squares), N (line with red circles), and C (line with blue triangles) in four configurations: a) The 32600/25600/32600 configuration, b) The 35860/25600/35860 configuration, c) The 39120/25600/39120 configuration, d) The 45640/25600/45640 configuration\u003c/p\u003e","description":"","filename":"Fig.11.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-8822403/v1/639d4f1c748de041ed53d7e5.jpeg"},{"id":104409021,"identity":"23364452-9318-4ed0-8ef4-3290c35acb97","added_by":"auto","created_at":"2026-03-11 12:43:57","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":6159950,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-8822403/v1/d083418b-e23a-4084-a146-1689441930c0.pdf"}],"financialInterests":"No competing interests reported.","formattedTitle":"Thermal Stability and Phase Transition Behavior in In-Plane Graphene/hBN/Graphene Heterostructures: A Molecular Dynamics Study","fulltext":[{"header":"Introduction","content":"\u003cp\u003eIn recent years, heterostructures have garnered increasing attention from researchers. By definition, the interface of two layers or regions of materials of different compositions is called a heterojunction [\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e]. For two-dimensional (2D) materials, heterostructures can be classified into two types: stacked 2D heterostructures and in-plane heterostructures [\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e]. The structure consisting of one or many heterojunctions is known as a heterostructure. As this novel structure allows bandgap engineering, heterostructures have become a critical area of research.\u003c/p\u003e \u003cp\u003eGraphene, a 2D material with a honeycomb structure, exhibits numerous intriguing characteristics [\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e], making it a promising candidate for various applications. Many attempts have been made to incorporate graphene with other compatible materials to create a heterostructure with desired properties. For this reason, other materials with the same honeycomb structure have also been intensely investigated. Silicene and the group-III-nitride family are the most well-known groups of hexagonal-structure materials that are frequently researched.\u003c/p\u003e \u003cp\u003eIn the group-III-nitride family, the hexagonal Boron Nitride (hBN) is analogous to graphene with a small mismatch (about 1.6%) [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e]. It has the potential to form van der Waals (vdW) interactions with graphene and produce heterostructures with novel characteristics. hBN interactions with living cells were also investigated. It was reported that at a nontoxic concentration, hBN can relieve molecular stress induced by drugs [\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e]. The influence of hBN nanoparticles of different doses was evaluated in Wistar albino rats, in which it was shown that only high doses of hBN nanoparticles (1600 and 3200 \u0026micro;g/kg) caused damage in the liver, kidney, heart, spleen, and pancreas, while lower doses (from 50 to 80 \u0026micro;g/kg) were considered as non-toxic [\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e]. hBN was also reported to increase dendritic cell maturation without inducing downstream T cell polarization [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e]. In the biomedical field, hBN was suggested to be used in several areas: wound healing [\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e], prostate cancer treatment [\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e], and anticancer drug carriers [\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eDue to their relatively small lattice mismatch (about 1.6%) [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e] and the same crystal structure, hBN and graphene can be combined by vertically stacking together. The major three stacking patterns in most studies are: (1) AA stacking, which means C atoms are directly aligned above the B and N atoms; (2) AB stacking, which means one C atom is positioned above the N atom, while the other C atom is positioned above the hBN ring; and (3) AB stacking, meaning that a carbon atom is placed above the B, and the other one is above the hexagonal hBN ring [\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e]. The other type of heterostructure of hBN and graphene is the in-plane graphene/hBN heterostructure. The electronic band structures and properties of this lateral graphene/hBN heterostructure were expected to be different from those of pristine graphene and hBN, giving potential for the development of band gap-engineered applications in electronics and optoelectronics [\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e]. The in-plane graphene/hBN heterostructure can be used to study the DNA. Bingquan Luan and colleagues introduced an ssDNA stretching platform using a two-dimensional in-plane heterostructure composed of graphene and hBN, demonstrating that ssDNA could be stretched on an hBN nanostripe sandwiched between two adjacent graphene domains (\"nanochannel\") [\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e]. Fabio and his group explored the sensitivity and selectivity of a graphene/hBN-based nanopore architecture for detecting and distinguishing synthetic Hachimoji nucleobases [\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e]. Ali Kiakojouri and his team proposed an in-plane graphene/hBN/graphene heterostructure, which exhibited higher sensitivity to different nucleotides compared to a similar structure made of pristine graphene, making it a promising candidate for DNA sequencing applications [\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eTo further investigate how the phase transition behavior of such in-plane graphene/hBN/graphene heterostructures depends on the relative composition of graphene and hBN, we constructed four atomic graphene/hBN/graphene configurations in which the number of hBN atoms was kept constant while the number of graphene atoms was varied such as 32600/25600/32600, 35860/25600/35860, 39120/25600/39120, and 45640/25600/45640 atom arrangements. This approach allows for a detailed assessment of how the graphene content influences the thermal properties and structural stability of the hybrid system. The computational methodology is described in Section 2, followed by the presentation of results in Section 3. Section 4 provides a detailed discussion, and the final section summarizes the main conclusions of the study.\u003c/p\u003e"},{"header":"Methods and Implementation","content":"\u003cp\u003eThe role of the interaction potential is critical in molecular dynamics (MD) simulations as it governs how particles behave and move within the simulated system. In MD simulations, the key aim is to accurately replicate real-world interactions between atoms or molecules in a computationally feasible way. The interaction potential encompasses the forces and energies of various atomic interactions, including bonded and non-bonded interactions. In systems with covalent bonds, where chemical bonds connect atoms, the interaction potential dictates the stretching, bending, and torsional forces that occur as atoms undergo movement and rearrangement. This capability allows the simulation to accurately depict processes like bond breaking and formation, which are essential for studying chemical reactions and structural transitions. Various empirical potentials, including those representing two- and three-body interactions in different forms, are available [\u003cspan additionalcitationids=\"CR17 CR18 CR19\" citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e]. The Tersoff potential specifically addresses the interatomic forces and energies involved in covalent bonding in materials, making them highly suitable for studying systems characterized by strong chemical interactions [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e, \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e]. What distinguishes the Tersoff potential is its ability to model complex phenomena such as bond breaking and reformation, anisotropic and directional bonding, and interactions at both short and long ranges. This broad representation extends their usefulness from basic elements like carbon to a diverse range of covalent materials such as semiconductors and heterostructures.\u003c/p\u003e \u003cp\u003eIn this study, the Tersoff potential [\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e] is used for the interactions between Boron (B), Nitrogen (N), and Carbon (C) in the initial configurations, which are written as below:\u003cdiv id=\"Equ1\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ1\" name=\"EquationSource\"\u003e\n$$\\:{E}_{b}=\\frac{1}{2}{\\sum\\:}_{i\\ne\\:j}^{}{f}_{c}\\left({r}_{ij}\\right)\\left[{f}_{R}\\right({r}_{ij})+{f}_{a}({r}_{ij}\\left)\\right]$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e1\u003c/div\u003e\u003c/div\u003e.\u003c/p\u003e \u003cp\u003ewhere, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{r}_{ij}\\)\u003c/span\u003e\u003c/span\u003e is the distance from the atom i to atom j. The repulsive \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{f}_{R}\\left({r}_{ij}\\right)\\)\u003c/span\u003e\u003c/span\u003e and the attractive \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{f}_{a}\\left({r}_{ij}\\right)\\)\u003c/span\u003e\u003c/span\u003e terms are based on the Morse potential as proposed by Brenner [\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e]. The \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{f}_{c}\\left({r}_{ij}\\right)\\)\u003c/span\u003e\u003c/span\u003e term represents the cutoff function used for calculating the number of neighbors as well as making the potential zero outside the interaction shell. To perform the calculations, the software package Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) is used [\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eFirst, to prepare the size of the hBN layer, we based on a previous study on the influence of the armchair/zigzag edge ratio on the melting process of free-standing hBN [\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e]. Accordingly, the configuration with an armchair/zigzag edge ratio of 10.745098 and the number of atoms in the model of 25600 is selected for investigation in this study.\u003c/p\u003e \u003cp\u003eThen, the hBN configuration will be hybridized with graphene to form an in-plane graphene/hBN/graphene heterostructure (C\u0026ndash;B bridge). The interaction at the edge of the graphene and the hBN layer can be seen in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eTo study the dependence on the relative amount of hBN to graphene of the melting process of in-plane hybrid graphene/hBN/graphene, four different in-plane graphene/hBN/graphene configurations are considered by changing the number of atoms of the graphene layer (32600/25600/32600, 35860/25600/35860, 39120/25600/39120, and 45640/25600/45640 atom arrangements).\u003c/p\u003e \u003cp\u003eAfter that, the hybrid in-plane graphene/hBN/graphene configurations are relaxed again for \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:6\\times\\:1{0}^{5}\\)\u003c/span\u003e\u003c/span\u003e MD steps to ensure that the system reaches mechanical and thermal equilibrium before the heating process begins.\u003c/p\u003e \u003cp\u003eFinally, the hybrid graphene/hBN/graphene configuration is heated from 50 K to 8000 K to ensure that the graphene/hBN/graphene configuration is in a liquid state.\u003c/p\u003e \u003cp\u003eThe change in temperature upon heating obeys \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{T}_{0}+vt\\)\u003c/span\u003e\u003c/span\u003e, in which \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{T}_{0}=50\\)\u003c/span\u003e\u003c/span\u003e K, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:v\\)\u003c/span\u003e\u003c/span\u003e is the heating rate \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:5\\times\\:1{0}^{12}\\)\u003c/span\u003e\u003c/span\u003e K/s, and t is the time required for heating. All configurations at the chosen temperatures are relaxed for \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:6\\times\\:1{0}^{5}\\)\u003c/span\u003e\u003c/span\u003e MD steps to ensure their stability before studying structural evolution or 2D visualization. The time interval for each step is 0.0001 picoseconds. We use Visual Molecular Dynamics (VMD) [\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e] to visualize atomic configurations.\u003c/p\u003e"},{"header":"Results and Discussions","content":"\u003cp\u003e \u003c/p\u003e \u003cp\u003eGenerally, the atoms of the solid structure are packed tightly together in an orderly arrangement, meaning that atoms vibrate about a fixed point. The fixed position of atoms minimizes the total energy per atom. Upon heating, the total energy increases as kinetic energy is introduced to the system, causing the atoms to vibrate rapidly. When the total energy per atom reaches a critical threshold, the amplitude of these vibrations becomes sufficient to overcome the intermolecular forces operating within the system, initiating a phase transition from the crystalline solid to the liquid state. Therefore, analyzing the total energy per atom provides informative insights into the stable configuration, energy barriers, and phase transitions, such as the solid-to-liquid transition.\u003c/p\u003e \u003cp\u003eTo investigate how varying the hBN:graphene ratio affects melting in an in-plane graphene/hBN/graphene hybrid, the energy per atom of the four configurations is calculated and shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e. Overall, the four configurations exhibit a first-order phase transition from the crystalline to the liquid phase (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). The pattern of first order type includes the three separate regions: i) at low temperatures (below T\u003csub\u003e1\u003c/sub\u003e), the total energy per atom increases linearly as the atom vibrates around the equilibrium position; ii) when increasing to a critical temperature point (from T\u003csub\u003e1\u003c/sub\u003e to T\u003csub\u003e3\u003c/sub\u003e), there are sudden jumps in total energy per atom, reflecting the breakdown of the crystalline lattice and the onset of disorder. The disorder arrangement proves a phase transition from the crystalline state to the liquid state of the system; and iii) after the phase transition temperature (from T\u003csub\u003e3\u003c/sub\u003e), the total energy per atom of the system again increases linearly with temperature, indicating the dominance of translational motion in the liquid state.\u003c/p\u003e \u003cp\u003eAlthough the four configurations exhibit a similar tendency to the first-order type, the phase transition temperature increases with the number of atoms per graphene layer. The 32600/25600/32600 configuration comprises the smallest amount of graphene and has the lowest melting temperature range, from 5850 K to 6320 K, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e (line with crosses). The 35860/25600/35860 configuration (line with hollow circles in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e) and the 39120/25600/39120 configuration (line with crossed diamonds in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e) exhibit higher phase transition temperatures, ranging from 6050 K to 6450 K and from 6240 K to 6550 K, respectively. The highest phase transition temperature observed in the 45640/25600/45640 configuration (line with black squares in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e) is 6790 K.\u003c/p\u003e \u003cp\u003eOne can see that an upward shift in the phase transition temperature with increasing graphene content in the in-plane hybrid graphene/hBN/graphene systems can be attributed to several points: i) the extended graphene domains, composed of strong sp\u0026sup2;-hybridized carbon networks, enhance the rigidity and bonding strength of the structure, making it more resistant to thermal disruptions; ii) in addition, the outer graphene layers function as mechanical anchors, restricting vibrational amplitudes and stabilizing the intermediate hBN region. This structural reinforcement also promotes uniform heat distribution, reducing local fluctuations that might otherwise initiate disorder. As a result, the hybrid system exhibits greater thermal tolerance, requiring higher temperatures to transition from a crystalline to a disordered phase.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWhile the energy per atom provides a general overview of the melting process, analyzing heat capacity identifies the melting points of four configurations of hybrid graphene/hBN/graphene materials. The heat capacity is calculated according to the expression:\u003cdiv id=\"Equ2\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ2\" name=\"EquationSource\"\u003e\n$$\\:{C}_{p}=\\:\\frac{dE}{dT}$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e2\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eDuring a solid-to-liquid phase transition, atoms gain mobility, and the system undergoes significant reorganization, leading to a sharp change in energy-absorption behavior. This is manifested as a distinct peak in the heat capacity curve. Consequently, the temperature at which this peak occurs can be reliably used to identify the onset of melting.\u003c/p\u003e \u003cp\u003eFor ease of observation, we will present the heat capacity in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e and Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e. As shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e and Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e, the peak in heat capacity corresponds to the disruption of C\u0026ndash;C, B\u0026ndash;N, and N\u0026ndash;B bonds, marking the transition from the crystalline to the liquid phase. The peak heat capacity of four configurations 32600/25600/32600, 35860/25600/35860, 39120/25600/39120, and 45640/25600/45640 is defined at 6120 K, 6270 K, 6390 K, and 6650 K, respectively (T\u003csub\u003e2\u003c/sub\u003e in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e and Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). The disparity in temperature points between the four configurations indicates the influence of the number of graphene atoms on the melting process of the hybrid graphene/hBN/graphene configuration. The rising melting temperature associated with increasing graphene size across these configurations suggests a relationship between graphene size and melting behavior. The visualization of four graphene/hBN/graphene configurations was presented in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e using VMD. Figure\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e illustrates the atomic structures of four graphene/hBN/graphene configurations at two temperature stages, before and after the melting point. As the temperature increases, significant structural disorder is observed, indicating the breakdown of C\u0026ndash;C, B\u0026ndash;N, and N\u0026ndash;B bonds. The transition becomes more evident with increasing graphene content, which shifts the onset of melting to higher temperatures. This trend highlights the role of graphene size in enhancing the thermal stability of the heterostructures. The structural evolution was visualized using VMD.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eThe coordination number (CN) provides a sensitive probe into the local atomic environment during heating and is therefore highly effective in characterizing the melting process of the graphene/hBN/graphene hybrid heterostructures. In the initial crystalline state at lower temperatures, the graphene and hBN regions form a honeycomb lattice in which each atom is threefold coordinated (CN\u0026thinsp;=\u0026thinsp;3). Nearly all C, B, and N atoms satisfy this structural requirement, except for a negligible fraction at the interfaces and boundaries. This is consistent with the stability of the sp\u0026sup2; bonding framework in both graphene and hBN.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eUpon heating, thermal vibrations increase progressively, leading to bond stretching and breaking. The destruction of C\u0026ndash;C, B\u0026ndash;N, and N\u0026ndash;B bonds reduces the number of atoms maintaining CN\u0026thinsp;=\u0026thinsp;3, accompanied by the appearance of atoms with lower coordination number (CN\u0026thinsp;=\u0026thinsp;2, 1, or 0). At a certain temperature, the CN plot exhibits sudden collapse, corresponding to the melted point of the material. By monitoring the evolution of coordination numbers for B, N, and C atoms as functions of temperature, the phase transition can be identified and correlated with other thermodynamic indicators such as the total energy per atom and the heat capacity.\u003c/p\u003e \u003cp\u003eThe transition temperatures and melting points vary systematically with the graphene:hBN ratio. As the number of graphene atoms increases, the critical temperatures shift to higher values, showing that graphene layers act as stabilizing components to the heterostructure.\u003c/p\u003e \u003cp\u003eIn the 32600/25600/32600 configuration, at early temperature points (3750 K to 4750 K), the CN plots of B, N, and C atoms decrease at steady rate (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea). At 3750 K, the CN plots of B and N are both approximately 50% (line with black squares for B, line with red circles for N in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea), much higher than that of C (less than 25% at this temperature point, line with blue triangles for C in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea). At about 4750 K, the CN plots of B and N collapse sharply (line with black squares for B, line with red circles for N in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea), marking the initial disordering of the hBN sublattice. The collapse of C occurs between 6100\u0026ndash;6250 K (line with blue triangles for C in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea), overlapping with the identified melting point at about 6120 K in the heat capacity curve (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea). The fact that the CN plot of C degrades later than that of B and N reflects the delayed disordering of graphene relative to hBN. At higher temperature points (from 6500 K), both B and N atoms reach near-zero coordination (line with black squares for B, line with red circles for N in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea), while C atoms retain a small fraction, but observably greater than 0, of CN\u0026thinsp;=\u0026thinsp;3 (line with blue triangles for C in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ea). After the steep decline of C, the CN plot of C is consistently higher than that of B and N, which is opposite to the observed trend at lower temperatures.\u003c/p\u003e \u003cp\u003eIn the 35860/25600/35860 configuration, the CN plots of B, N, and C have a similar appearance to the CN plots of the 32600/25600/32600 configuration. However, at 3750 K, while the percentage of CN\u0026thinsp;=\u0026thinsp;3 of C is roughly 25% (line with blue triangles for C in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eb), which is about the same that of C in the 32600/25600/32600 configuration, the percentage of CN\u0026thinsp;=\u0026thinsp;3 of B and N is only less than 45% (line with black squares for B, line with red circles for N in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eb). In the 4750 K\u0026thinsp;\u0026minus;\u0026thinsp;5000 K range, the CN plots of B and N drop dramatically from 30% to nearly 0% (line with black squares for B, red circle for N in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eb). The collapse of C\u0026rsquo;s CN plot happens at a higher temperature, about 6100 K (line with blue triangles for C in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eb).\u003c/p\u003e \u003cp\u003eIn the 39120/25600/39120 configuration, at 3750 K, the CN plots of B and N are around 25\u0026ndash;30% (line with black squares for B, red circle for N in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec). Meanwhile, the plot of C is in the vicinity of 20% (line with blue triangles for C in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec). Like the plots of other discussed configurations, at lower temperature points, the plots of B and N are observably higher than that of C. The C atoms in the 39120/25600/39120 configuration also show stronger resilience against heat compared to B and N atoms: the plot of C degrades at higher temperature, about 6250 K (line with blue triangles for C in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec), while the plots of B and N deteriorate at roughly 4750 K (line with black squares for B, line with red circles for N in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec). At temperatures higher than 6500 K, the plot of C continuously decreases as the temperature increases, but still remains well above 0% (line with blue triangles for C in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec); whereas the plots of B and N are close to 0% at these temperature points (line with black squares for B, line with red circles for N in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ec).\u003c/p\u003e \u003cp\u003eIn the 45640/25600/45640 configuration, the CN plots of B, N, and C demonstrate the most delayed transition among all examined configurations. At 3750 K, the CN plot of C does not indicate great differences from the plots of B and N like other configurations (all plots are approximately 8\u0026ndash;12% at this temperature point; line with blue triangles for C, black squares for B, and red circles for N in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ed). The CN plots of B and N begin to collapse above 6300 K, while the CN plot of C plummets at nearly 6700 K (black square for B, red circle for N in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ed). The persistence of C coordination underscores the dominant role of graphene in stabilizing the hybrid heterostructure. Thus, larger graphene layers act as a thermal backbone, resisting complete collapse and ensuring superior stability. At temperatures higher than 6500 K, while the plots of B and N are approximately 0% (line with black squares for B, line with red circles for N in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ed), the plot of C remains at about 3\u0026ndash;4% (line with blue triangles for C in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ed), constantly higher than that of B and N.\u003c/p\u003e \u003cp\u003eA closer comparison among atomic species reveals that B and N atoms lose coordination earlier and more rapidly at lower temperatures than those of carbon. This behavior can be attributed to the relative weakness of B\u0026ndash;N bonds compared to C\u0026ndash;C bonds. As temperature rises, B\u0026ndash;N bonds break first, destabilizing the hBN regions and initiating local disorder. Graphene, with its stronger C\u0026ndash;C bonds, resists this disorder longer, maintaining partial threefold coordination and thereby delaying the global phase transition. The result is a progressive shift of the melting point with increasing graphene content. The coordination number analysis is also consistent with energy and heat capacity calculations.\u003c/p\u003e \u003cp\u003eIn summary, the coordination number analysis demonstrates that melting in graphene/hBN/graphene hybrid heterostructures proceeds via a first-order phase transition. The critical transition temperatures and melting points systematically increase with the graphene:hBN ratio, highlighting the stabilizing role of graphene. The preferential survival of C\u0026ndash;C coordination to higher temperatures further confirms graphene as the primary contributor to thermal resilience. These results underline the importance of atomic composition in dictating the structural integrity of 2D in-plane hybrid graphene/hBN/graphene heterostructures under extreme thermal conditions.\u003c/p\u003e \u003cp\u003eTo analyze structural deformation and clarify the stabilizing role of graphene layers in planar graphene/hBN/graphene heterostructures during heating, the angular distribution of C\u0026ndash;C\u0026ndash;C and B\u0026ndash;N\u0026ndash;B was evaluated for four configurations (32600/25600/32600, 35860/25600/35860, 39120/25600/39120, and 45640/25600/45640). The analysis was performed at two characteristic temperatures: the melting start temperature, determined from the initial abrupt change in total energy per atom (T\u003csub\u003e1\u003c/sub\u003e in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e), and the temperature point just above the melting point (T\u003csub\u003e3\u003c/sub\u003e in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). For ease of reference, the results are presented in Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e to Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e10\u003c/span\u003e.\u003c/p\u003e \u003cp\u003eFirst, for the configuration with the lowest carbon content (the 32600/25600/32600 configuration), at the pre-melting temperature (5850 K), the angular distribution of both C\u0026ndash;C\u0026ndash;C and B\u0026ndash;N\u0026ndash;B remains narrow, exhibiting distinct peaks around 120\u0026deg; (approximately 0.8% for graphene and 0.9% for hBN). This demonstrates that the hexagonal lattice is largely preserved in both regions (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003ea) because at this temperature, thermal vibrations are insufficient to significantly disrupt the sp\u0026sup2; linkage network.\u003c/p\u003e \u003cp\u003eAs the temperature increases to 6320 K, the angular distribution of C\u0026ndash;C\u0026ndash;C expands significantly (spanning approximately 85\u0026deg;\u0026ndash;155\u0026deg;) while retaining a distinct peak at 120\u0026deg; with a decrease in intensity (\u0026asymp;\u0026thinsp;0.6%) (line with red circles for C\u0026ndash;C\u0026ndash;C in Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eb). The persistence of the peak at 120\u0026deg; suggests that the short-range hexagonal order in graphene is not completely destroyed at the melting point. Meanwhile, the B\u0026ndash;N\u0026ndash;B angular distribution becomes significantly less concentrated around 120\u0026deg;. Specifically, the peak intensity decreases to approximately 0.7% and the shape is noticeably distorted (line with black squares for B\u0026ndash;N\u0026ndash;B in Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eb). This marked expansion reflects a more severe angular disorder in the hBN sublattice.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWhen increasing the number of atoms in the graphene layer to 35860 atoms per layer (35860/25600/35860), the melting start temperature increased, implying that the stability of this configuration is significantly improved. Similar to the previous configuration, we will investigate the pre-melting temperature. At 6050 K, the peak at 120\u0026deg; of the C\u0026ndash;C\u0026ndash;C angle distribution is narrow and sharp (peak at approximately 0.9%), whereas the B\u0026ndash;N\u0026ndash;B angle distribution is flatter with more ripples at the top (peaks of these ripples are approximately 0.6%; line with red circles for C\u0026ndash;C\u0026ndash;C and line with black squares for B\u0026ndash;N\u0026ndash;B in Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003ea). The B\u0026ndash;N\u0026ndash;B angle distribution overlaps significantly. Upon reaching 6450 K, the C\u0026ndash;C\u0026ndash;C peak at 120\u0026deg; is broadened, but the peak remains visible (peak reduction to 0.7%, line with red circles for C\u0026ndash;C\u0026ndash;C in Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003eb), while the B\u0026ndash;N\u0026ndash;B begins to significantly deform (peak reduction to 0.5%, line with black squares in Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003eb), indicating a higher degree of structural disorder. Compared to the 32600/25600/32600 configuration, due to the higher amount of C in this configuration, the loss of integrity of the hBN layer occurs more slowly.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eContinuing to increase the carbon content in the graphene layer to 39120 atoms per layer (39120/25600/39120), the melting temperature begins to rise higher than the two previous configurations. At the pre-melting temperature of 6240 K, the angular distribution of C\u0026ndash;C\u0026ndash;C shows a clear and concentrated peak (peak at approximately 0.8%, line with red circles for C\u0026ndash;C\u0026ndash;C in Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003ea). The angular distribution of B\u0026ndash;N\u0026ndash;B is more dispersed, with the distribution at 120\u0026deg; accounts for about 0.5% and the highest peak is about 0.6% (line with black squares for B\u0026ndash;N\u0026ndash;B in Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003ea). Continuing to increase the temperature to 6550 K, the angular distribution of C\u0026ndash;C\u0026ndash;C becomes slightly more dispersed (peak decreases to 0.7%, line with red circles for C\u0026ndash;C\u0026ndash;C in Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003eb). The B\u0026ndash;N\u0026ndash;B angle distribution is flattened (line with black squares for B\u0026ndash;N\u0026ndash;B in Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003eb) but is significantly less affected by temperature than the previous 32600/25600/32600 and 35860/25600/35860 configurations. This implies that the abundance of C molecules in the graphene layer better inhibits hBN structural deformation.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eFinally, when the number of atoms in the graphene layer is 45640 atoms, the best heat resistance is observed. Specifically, the pre-melting temperature is 6430 K. At the initial melting temperature of 6430 K, the angular distribution of C\u0026ndash;C\u0026ndash;C is sharp and narrow (the angular distribution center is approximately 120\u0026deg;, with the peak at over 0.8% (line with red circles for C\u0026ndash;C\u0026ndash;C in Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e10\u003c/span\u003ea). On the other hand, B\u0026ndash;N\u0026ndash;B shows low structural integrity (line with black squares for B\u0026ndash;N\u0026ndash;B in Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e10\u003c/span\u003ea). When the temperature increases to 6790 K (melting point), the C\u0026ndash;C\u0026ndash;C angular distribution is affected but still maintains a peak at 120\u0026deg; (peak at approximately 0.7%, line with red circles for C\u0026ndash;C\u0026ndash;C in Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e10\u003c/span\u003eb). The B\u0026ndash;N\u0026ndash;B angle distribution is wide, but its variation is the least compared to other configurations (line with black squares for B\u0026ndash;N\u0026ndash;B in Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e10\u003c/span\u003eb). This configuration has the highest melting temperature, suggesting that the amount of C in the graphene layer has the strongest impact on hBN stabilization.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eOne can see that graphene contributes significantly to the stabilization of the hBN layer. In all configurations, the change in the B\u0026ndash;N\u0026ndash;B angle distribution is the most sensitive sign reflecting the structural degradation of the hBN layer. As the carbon content increases, this change is delayed, and the intensity is reduced, indicating the protective effect of graphene. At the same time, the C\u0026ndash;C\u0026ndash;C angle distribution still maintains a clear peak at high temperatures, demonstrating that the graphene network has better thermal resistance than hBN. It can be concluded that the angle distribution analysis provides microscopic evidence that graphene not only stabilizes itself but also plays an important role in maintaining the structure of the neighboring hBN layer. Therefore, adjusting the carbon content is an effective strategy to control the thermal stability of 2D heterojunction material systems, aiming at applications in electronic devices or materials operating in high-temperature environments.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eTo characterize the structural evolution of the hybrid graphene/hBN/graphene configuration under increasing temperature, this study employs the Bond-Orientational Order (BOO) parameter, which quantifies the degree of local orientational symmetry in solid-state systems. BOO provides a measure of the bond alignment around each particle, thereby offering insights into phase transitions and structural disorder that may not be fully captured by simple radial distribution functions. In particular, for hexagonal lattices commonly found in two-dimensional materials such as graphene and hBN, the sixth-order BOO parameter (Q₆) is especially relevant, as it reflects the inherent hexagonal symmetry of these systems. The value of Q₆ ranges from 0 to 1, with crystalline phases exhibiting values close to unity, while reductions in magnitude indicate enhanced thermal agitation, amorphization, or melting transitions. In our simulations, monitoring Q₆ as a function of temperature reveals the onset of structural instabilities, providing quantitative evidence of how elevated thermal energy alters the orientational order of bonds within the hybrid graphene/hBN/graphene configuration.\u003c/p\u003e \u003cp\u003eOverall, across all four configurations (Fig.\u0026nbsp;\u003cspan refid=\"Fig11\" class=\"InternalRef\"\u003e11\u003c/span\u003e), the BOO values associated with B, N, and C atoms at the initial temperature of 50 K are close to unity, confirming that the configurations possess nearly perfect orientational symmetry consistent with a hexagonal lattice. However, the BOO parameter subsequently decreases with rising temperature. This trend is consistent with the physical expectation that thermal fluctuations intensify atomic vibrations, thereby displacing atoms from their original positions. Nevertheless, the degree of reduction is not identical across configurations, nor is it uniform among different atomic species within the same configuration. For each configuration, the variations can be described as below.\u003c/p\u003e \u003cp\u003eFor the 32600/25600/32600 configuration (Fig.\u0026nbsp;\u003cspan refid=\"Fig11\" class=\"InternalRef\"\u003e11\u003c/span\u003ea), at 50 K, the BOO values of B and N are nearly identical, both around 0.99, while the corresponding value for C is slightly lower at 0.96. Overall, both the hBN and graphene layers exhibit a high degree of hexagonal orientational symmetry at the initial temperature (50 K). With respect to B and N, their evolution during heating follows almost the same trend: from 50 K to 3000 K, the BOO decreases minimally (reaching 0.96 for both elements at 3258 K); beyond 3500 K, the rate of decrease gradually accelerates until around 6472 K, where the drop becomes most significant with the BOO falling to 0.63, indicating a rapid structural transformation of the hBN layer at this temperature. Thereafter, the decline rate slows down, making the plot reaching a final value of 0.36 at 8000 K.\u003c/p\u003e \u003cp\u003eIn contrast, C shows a distinct three-stage behavior: in the first stage, from 50 K to 5600 K, the BOO decreases linearly from 0.96 to 0.66; in the second stage, between 5600 K and 6300 K, the parameter drops rapidly to 0.34, with the sharpest change observed at 6110 K, reflecting a clear distortion in the original hexagonal structure of the graphene layer, characteristic of the material\u0026rsquo;s phase transition. Finally, in the third stage, the curve continues to decrease linearly, reaching as low as 0.20.\u003c/p\u003e \u003cp\u003eFor the 35860/25600/35860 configuration (Fig.\u0026nbsp;\u003cspan refid=\"Fig11\" class=\"InternalRef\"\u003e11\u003c/span\u003eb), at 50 K, the BOO values for B and N are 0.99, while that for C is 0.96. Similar to the first configuration, the BOO of B and N decreases very slowly from 50 K to 3000 K, after which the rate of decline gradually accelerates until 6343 K, where it reaches 0.55, indicating a rapid structural transformation characteristic of the phase transition. The BOO then decreases more slowly to 0.29 at 8000 K. For C, the reduction occurs in three distinct stages: in the first stage, from 50 K to 5700 K, the BOO decreases linearly from 0.96 to 0.66; in the second stage, from 5700 K to 6300 K, the curve drops sharply to 0.38, confirming the occurrence of a structural phase transition most evident at 6250 K; and in the final stage, it continues to decline linearly, reaching 0.22 at the end of the simulation.\u003c/p\u003e \u003cp\u003eIn the 39120/25600/39120 configuration (Fig.\u0026nbsp;\u003cspan refid=\"Fig11\" class=\"InternalRef\"\u003e11\u003c/span\u003ec), the BOO values for both elements in the hBN layer are 0.99 at the initial temperature, then slightly decrease to 0.96 at 2500 K. From 2500 K to 6443 K, the rate of decrease accelerates, with the BOO reaching 0.45 at the end of this interval, representing the phase transition stage in the hBN layer and subsequently slowing down to a final value of 0.26 at 8000 K. For C, the decline follows a pattern similar to that observed in the first two configurations and can be divided into three distinct stages: stage 1, from 50 K to 6100 K, the BOO decreases linearly from 0.96 to 0.62; stage 2, from 6100 K to 6500 K, the parameter drops rapidly, reaching 0.35 by the end of this stage, indicating the temperature range where the phase transition occurs as the hexagonal structure becomes distorted and loses its initial symmetry, with the most pronounced jump observed at 6342 K; stage 3, from 6500 K to 8000 K, the decline returns to a linear trend, resulting in a final value of 0.22. However, unlike the two previous configurations, where the BOO of B and N remained consistently higher than that of C throughout heating, the third configuration exhibits a reversed trend between 4900 K and 6400 K, during which the BOO of C surpasses that of the other two elements. This behavior arises because the structural degradation in the hBN layer initiates earlier at lower temperatures, whereas the rapid decline in C occurs at higher temperatures compared to the preceding configurations.\u003c/p\u003e \u003cp\u003eIn the 45640/25600/45640 configuration (Fig.\u0026nbsp;\u003cspan refid=\"Fig11\" class=\"InternalRef\"\u003e11\u003c/span\u003ed), the BOO curves of B and N display earlier structural changes compared to the previous configurations. Initially, the BOO is 0.99, then the rate of decrease accelerates rapidly up to around 6000 K, at which point the BOO drops to 0.62, representing a rapid structural deformation characteristic of this phase-transition stage. Thereafter, the decrease slows progressively, reaching 0.24 at 8000 K. The BOO of C continues to follow a three-stage pattern: stage 1, from 50 K to 6200 K, the BOO decreases linearly from 0.96 to 0.65; stage 2, from 6200 K to 6800 K, where a sharp drop is observed particularly at 6648 K to 0.36 at the end of this stage, confirming the phase transition in the graphene layer; stage 3, beyond this point, the curve resumes a linear decrease, reaching a final value of 0.24 at 8000 K. In this configuration, a clear distinction can be observed: for most of the simulated temperature range (from 3800 K to 8000 K), the BOO of C remains higher than that of the other two elements.\u003c/p\u003e \u003cp\u003eIn general comparison, the transformation stages across all four configurations exhibit similar trends. However, as the number of atoms in the graphene layer increases, the BOO variation in the hBN layer occurs at lower temperatures, whereas the BOO variation in the graphene layer shifts to higher temperatures. This leads to noticeable changes in the BOO-temperature distributions among the configurations: in the first configuration, the BOO of the hBN layer consistently remains higher than that of graphene, resulting in two clearly separated curves; in the second configuration, the two curves gradually converge; and in the last two configurations, a distinct temperature interval appears where the BOO of graphene exceeds that of hBN. This trend is most evident in the progressive increase of the temperature at which the sharp BOO drop in graphene occurs, from 6110 K in configuration 32600/25600/32600, to 6250 K in configuration 35860/25600/35860, and 6342 K and 6648 K in configurations 39120/25600/39120 and 45640/25600/45640, respectively. These values are consistent with the preceding peak positions of the Cp curves (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e and Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e), which are predicted to correspond to the phase transition points of the configurations at approximately 6120 K, 6270 K, 6390 K, and 6650 K, respectively. Although the hBN and graphene layers exhibit opposite trends as the number of Carbon atoms increases, the average BOO values of the two layers across all four configurations remain within the range of 0.20\u0026ndash;0.35. This indicates that a significant portion of the initial hexagonal orientational structure is still preserved within the configuration and has not been destroyed even at the final simulation temperature of 8000 K.\u003c/p\u003e"},{"header":"Conclusions","content":"\u003cp\u003eThe melting process of 32600/25600/32600, 35860/25600/35860, 39120/25600/39120, and 45640/25600/45640 configurations of the in-plane hybrid graphene/hBN/graphene heterostructure is investigated using MD simulations. The main results are:\u003c/p\u003e \u003cp\u003e \u003cul\u003e \u003cli\u003e \u003cp\u003eAn analysis of the energy per atom reveals that the four hybrid configurations of graphene/hBN/graphene exhibit a first-order phase transition. Configurations with more graphene atoms exhibit a phase transition at higher temperature ranges, indicating greater thermal tolerance. Additionally, examining the heat capacity of the hybrid system provides insights into how the melting point varies with the number of graphene atoms. The results indicate that as the number of graphene atoms increases, the peak heat capacity also rises, with peaks observed at 6120 K, 6270 K, 6390 K, and 6650 K. This finding highlights the relationship between graphene size and the melting behavior of the hybrid configurations.\u003c/p\u003e \u003c/li\u003e \u003cli\u003e \u003cp\u003eThe coordination number analysis reveals a clear first-order melting transition in the graphene/hBN/graphene hybrid heterostructures under heating. All four configurations exhibit a gradual reduction in threefold atomic coordination at lower temperatures, followed by an abrupt collapse at higher temperatures corresponding to the melting point. The B and N atoms in hBN lose their coordination earlier than the C atoms, confirming weaker B\u0026ndash;N bonding compared to C\u0026ndash;C bonding. Increasing the graphene fraction systematically shifts the transition to higher temperatures, demonstrating that graphene enhances structural stability under extreme thermal conditions. The persistence of CN\u0026thinsp;=\u0026thinsp;3 in carbon atoms to higher temperatures highlights the dominant role of graphene as a thermally resilient backbone. These results are consistent with heat capacity behavior and energy trends, emphasizing that atomic composition strongly governs the thermal resistance of graphene/hBN hybrid systems.\u003c/p\u003e \u003c/li\u003e \u003cli\u003e \u003cp\u003eThrough observations of the angular distributions at the melting onset temperature and the temperature point just above the melting point of four configurations, one can conclude that increasing the carbon content in the graphene layer helps to improve the thermal stability of the graphene/hBN/graphene heterostructure. As the carbon content increases, the change in the B\u0026ndash;N\u0026ndash;B angle distribution is delayed and the peak value of the distribution is less changed, indicating the protective effect of graphene.\u003c/p\u003e \u003c/li\u003e \u003cli\u003e \u003cp\u003eTo comprehensively evaluate the orientational order and the preservation of the hexagonal framework during the heating process of both hBN and graphene layers within the hybrid material, this study employs the BOO parameter. Across all four configurations, the BOO values decrease with increasing temperature due to enhanced atomic vibrations that gradually distort the ideal hexagonal arrangement. As the number of C atoms in the graphene layer increases, an opposite trend arises between the two layers: the hBN layer exhibits a faster BOO reduction, indicating earlier structural transformations, whereas the graphene layer shows a slower variation with higher carbon concentration. This contrasting behavior strongly affects the predicted phase transition temperature, as the abrupt BOO drop in graphene coincides closely with the Cp peaks identified previously.\u003c/p\u003e \u003c/li\u003e \u003c/ul\u003e \u003c/p\u003e"},{"header":"Declarations","content":"\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\u003cp\u003eHang T. T. Nguyen conceived of the presented idea, performed the computations, analyzed the data, supervised the findings of this work, and wrote the manuscript. Dan Lin Lieu analyzed the data, wrote the manuscript, and submitted the manuscript. Nghia Minh Dong analyzed the data and wrote the manuscript. Yen Khanh Nguyen analyzed the data and wrote the manuscript. Tin H. Nguyen analyzed the data and wrote the manuscript.\u003c/p\u003e\u003ch2\u003eAcknowledgement\u003c/h2\u003e\u003cp\u003eWe acknowledge Ho Chi Minh City University of Technology (HCMUT), VNU-HCM for supporting this study.\u003c/p\u003e\u003ch2\u003eData Availability\u003c/h2\u003e\u003cp\u003eData sets generated during the current study are available from the corresponding author on reasonable request.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eN. F. Andrade Neto, A. B. Lima, R. R. Y. O. V. Wilson, T. C. N. Nicacio, M. R. D. Bomio, and F. V. Motta, \u0026ldquo;Heterostructures obtained by ultrasonic methods for photocatalytic application: A review,\u0026rdquo; \u003cem\u003eMaterials Science in Semiconductor Processing\u003c/em\u003e, vol. 139, p. 106311, Mar. 2022, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.mssp.2021.106311\u003c/span\u003e\u003cspan address=\"10.1016/j.mssp.2021.106311\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eH. Wang, F. Liu, W. Fu, Z. Fang, W. Zhou, and Z. Liu, \u0026ldquo;Two-dimensional heterostructures: fabrication, characterization, and application,\u0026rdquo; \u003cem\u003eNanoscale\u003c/em\u003e, vol. 6, no. 21, pp. 12250\u0026ndash;12272, Aug. 2014, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1039/C4NR03435J\u003c/span\u003e\u003cspan address=\"10.1039/C4NR03435J\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eN. Sharma, R. Dev Gupta, R. Chandmal Sharma, S. Dayal, and A. Singh Yadav, \u0026ldquo;Graphene: An overview of its characteristics and applications,\u0026rdquo; \u003cem\u003eMaterials Today: Proceedings\u003c/em\u003e, vol. 47, pp. 2752\u0026ndash;2755, 2021, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/j.matpr.2021.03.086\u003c/span\u003e\u003cspan address=\"10.1016/j.matpr.2021.03.086\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eL. H. Li and Y. Chen, \u0026ldquo;Atomically Thin Boron Nitride: Unique Properties and Applications,\u0026rdquo; \u003cem\u003eAdv Funct Materials\u003c/em\u003e, vol. 26, no. 16, pp. 2594\u0026ndash;2608, Apr. 2016, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1002/adfm.201504606\u003c/span\u003e\u003cspan address=\"10.1002/adfm.201504606\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eI. C. Taskin, O. Sen, M. Emanet, M. Culha, and B. Yilmaz, \u0026ldquo;Hexagonal boron nitrides reduce the oxidative stress on cells,\u0026rdquo; \u003cem\u003eNanotechnology\u003c/em\u003e, vol. 31, no. 21, p. 215101, May 2020, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1088/1361-6528/ab6fdc\u003c/span\u003e\u003cspan address=\"10.1088/1361-6528/ab6fdc\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eF. Kar et al., \u0026ldquo;In Vivo Assessment of the Effect of Hexagonal Boron Nitride Nanoparticles on Biochemical, Histopathological, Oxidant and Antioxidant Status,\u0026rdquo; \u003cem\u003eJ Clust Sci\u003c/em\u003e, vol. 32, no. 2, pp. 517\u0026ndash;529, Mar. 2021, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1007/s10876-020-01811-w\u003c/span\u003e\u003cspan address=\"10.1007/s10876-020-01811-w\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eH. Lin et al., \u0026ldquo;2D Materials and Primary Human Dendritic Cells: A Comparative Cytotoxicity Study,\u0026rdquo; \u003cem\u003eSmall\u003c/em\u003e, vol. 18, no. 20, p. 2107652, May 2022, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1002/smll.202107652\u003c/span\u003e\u003cspan address=\"10.1002/smll.202107652\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003e\u0026Ouml;. Şen, M. Emanet, and M. \u0026Ccedil;ulha, \u0026ldquo;Stimulatory Effect of Hexagonal Boron Nitrides in Wound Healing,\u0026rdquo; \u003cem\u003eACS Appl. Bio Mater.\u003c/em\u003e, vol. 2, no. 12, pp. 5582\u0026ndash;5596, Dec. 2019, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1021/acsabm.9b00669\u003c/span\u003e\u003cspan address=\"10.1021/acsabm.9b00669\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eM. Emanet Ciofani, \u0026Ouml;. Şen, and M. \u0026Ccedil;ulha, \u0026ldquo;Hexagonal Boron Nitride Nanoparticles for Prostate Cancer Treatment,\u0026rdquo; \u003cem\u003eACS Appl. Nano Mater.\u003c/em\u003e, vol. 3, no. 3, pp. 2364\u0026ndash;2372, Mar. 2020, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1021/acsanm.9b02486\u003c/span\u003e\u003cspan address=\"10.1021/acsanm.9b02486\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eS. M. Sharker, \u0026ldquo;Hexagonal Boron Nitrides (White Graphene): A Promising Method for Cancer Drug Delivery,\u0026rdquo; \u003cem\u003eInt J Nanomedicine\u003c/em\u003e, vol. 14, pp. 9983\u0026ndash;9993, 2019, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.2147/IJN.S205095\u003c/span\u003e\u003cspan address=\"10.2147/IJN.S205095\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJ. Wang, F. Ma, and M. Sun, \u0026ldquo;Graphene, hexagonal boron nitride, and their heterostructures: properties and applications,\u0026rdquo; \u003cem\u003eRSC Adv.\u003c/em\u003e, vol. 7, no. 27, pp. 16801\u0026ndash;16822, 2017, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1039/C7RA00260B\u003c/span\u003e\u003cspan address=\"10.1039/C7RA00260B\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eH. Wang, F. Liu, W. Fu, Z. Fang, W. Zhou, and Z. Liu, \u0026ldquo;Two-dimensional heterostructures: fabrication, characterization, and application,\u0026rdquo; \u003cem\u003eNanoscale\u003c/em\u003e, vol. 6, no. 21, pp. 12250\u0026ndash;12272, Aug. 2014, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1039/C4NR03435J\u003c/span\u003e\u003cspan address=\"10.1039/C4NR03435J\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eB. Luan and R. Zhou, \u0026ldquo;Spontaneous ssDNA stretching on graphene and hexagonal boron nitride in plane heterostructures,\u0026rdquo; \u003cem\u003eNat Commun\u003c/em\u003e, vol. 10, no. 1, p. 4610, Oct. 2019, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1038/s41467-019-12584-w\u003c/span\u003e\u003cspan address=\"10.1038/s41467-019-12584-w\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eA. Kiakojouri, I. Frank, and E. Nadimi, \u0026ldquo;In-plane graphene/h-BN/graphene heterostructures with nanopores for electrical detection of DNA nucleotides,\u0026rdquo; \u003cem\u003ePhys. Chem. Chem. Phys.\u003c/em\u003e, vol. 23, no. 44, pp. 25126\u0026ndash;25135, 2021, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1039/D1CP03597E\u003c/span\u003e\u003cspan address=\"10.1039/D1CP03597E\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eF. A. L. De Souza, G. Sivaraman, M. Fyta, R. H. Scheicher, W. L. Scopel, and R. G. Amorim, \u0026ldquo;Electrically sensing Hachimoji DNA nucleotides through a hybrid graphene/ h -BN nanopore,\u0026rdquo; \u003cem\u003eNanoscale\u003c/em\u003e, vol. 12, no. 35, pp. 18289\u0026ndash;18295, 2020, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1039/D0NR04363J\u003c/span\u003e\u003cspan address=\"10.1039/D0NR04363J\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eG. C. Abell, \u0026ldquo;Empirical chemical pseudopotential theory of molecular and metallic bonding,\u0026rdquo; \u003cem\u003ePhys. Rev. B\u003c/em\u003e, vol. 31, no. 10, pp. 6184\u0026ndash;6196, May 1985, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1103/PhysRevB.31.6184\u003c/span\u003e\u003cspan address=\"10.1103/PhysRevB.31.6184\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eD. W. Brenner, \u0026ldquo;Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films,\u0026rdquo; \u003cem\u003ePhys. Rev. B\u003c/em\u003e, vol. 42, no. 15, pp. 9458\u0026ndash;9471, Nov. 1990, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1103/PhysRevB.42.9458\u003c/span\u003e\u003cspan address=\"10.1103/PhysRevB.42.9458\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eM. B\u0026uuml;ttiker, Y. Imry, R. Landauer, and S. Pinhas, \u0026ldquo;Generalized many-channel conductance formula with application to small rings,\u0026rdquo; \u003cem\u003ePhys. Rev. B\u003c/em\u003e, vol. 31, no. 10, pp. 6207\u0026ndash;6215, May 1985, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1103/PhysRevB.31.6207\u003c/span\u003e\u003cspan address=\"10.1103/PhysRevB.31.6207\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJ. Tersoff, \u0026ldquo;New empirical approach for the structure and energy of covalent systems,\u0026rdquo; \u003cem\u003ePhys. Rev. B\u003c/em\u003e, vol. 37, no. 12, pp. 6991\u0026ndash;7000, Apr. 1988, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1103/PhysRevB.37.6991\u003c/span\u003e\u003cspan address=\"10.1103/PhysRevB.37.6991\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eF. H. Stillinger and T. A. Weber, \u0026ldquo;Computer simulation of local order in condensed phases of silicon,\u0026rdquo; \u003cem\u003ePhys. Rev. B\u003c/em\u003e, vol. 31, no. 8, pp. 5262\u0026ndash;5271, Apr. 1985, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1103/PhysRevB.31.5262\u003c/span\u003e\u003cspan address=\"10.1103/PhysRevB.31.5262\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJ. Tersoff, \u0026ldquo;Modeling solid-state chemistry: Interatomic potentials for multicomponent systems,\u0026rdquo; \u003cem\u003ePhys. Rev. B\u003c/em\u003e, vol. 39, no. 8, pp. 5566\u0026ndash;5568, Mar. 1989, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1103/PhysRevB.39.5566\u003c/span\u003e\u003cspan address=\"10.1103/PhysRevB.39.5566\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eA. Kınacı, J. B. Haskins, C. Sevik, and T. \u0026Ccedil;ağın, \u0026ldquo;Thermal conductivity of BN-C nanostructures,\u0026rdquo; \u003cem\u003ePhys. Rev. B\u003c/em\u003e, vol. 86, no. 11, p. 115410, Sept. 2012, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1103/PhysRevB.86.115410\u003c/span\u003e\u003cspan address=\"10.1103/PhysRevB.86.115410\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eD. W. Brenner, \u0026ldquo;Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films,\u0026rdquo; \u003cem\u003ePhys. Rev. B\u003c/em\u003e, vol. 42, no. 15, pp. 9458\u0026ndash;9471, Nov. 1990, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1103/PhysRevB.42.9458\u003c/span\u003e\u003cspan address=\"10.1103/PhysRevB.42.9458\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eS. Plimpton, \u0026ldquo;Fast Parallel Algorithms for Short-Range Molecular Dynamics,\u0026rdquo; \u003cem\u003eJournal of Computational Physics\u003c/em\u003e, vol. 117, no. 1, pp. 1\u0026ndash;19, Mar. 1995, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1006/jcph.1995.1039\u003c/span\u003e\u003cspan address=\"10.1006/jcph.1995.1039\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eH. T. T. Nguyen, \u0026ldquo;Edge effects on the melting process of two-dimensional hexagonal boron nitride,\u0026rdquo; \u003cem\u003eJ Nanopart Res\u003c/em\u003e, vol. 26, no. 8, p. 199, Aug. 2024, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1007/s11051-024-06108-x\u003c/span\u003e\u003cspan address=\"10.1007/s11051-024-06108-x\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eW. Humphrey, A. Dalke, and K. Schulten, \u0026ldquo;VMD: Visual molecular dynamics,\u0026rdquo; \u003cem\u003eJournal of Molecular Graphics\u003c/em\u003e, vol. 14, no. 1, pp. 33\u0026ndash;38, Feb. 1996, doi: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.1016/0263-7855(96)00018-5\u003c/span\u003e\u003cspan address=\"10.1016/0263-7855(96)00018-5\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"molecular dynamics simulations, in-plane Graphene/hBN heterostructure, Thermal stability, Melting behavior, BOO parameter","lastPublishedDoi":"10.21203/rs.3.rs-8822403/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-8822403/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eThe lateral hybridization of graphene is studied using molecular dynamics simulation. The goal is to understand the structural and thermal properties of the in-plane hybrid graphene/hBN/graphene heterostructure, focusing on the dependence of the phase transition on the relative amount of hBN to graphene in the in-plane graphene/hBN/graphene hybrid heterostructure. The Tersoff potential is used to describe covalent bonds between carbon, boron, and nitrogen atoms. Four in-plane graphene/hBN/graphene hybrid configurations were constructed, in which the numbers of atoms in the left graphene, central hBN, and right graphene regions were 32600/25600/32600, 35860/25600/35860, 39120/25600/39120, and 45640/25600/45640, respectively. These configurations were used to examine the effect of the relative graphene-to-hBN ratio on the phase transition behavior. All systems exhibited first-order phase transitions, with broader transition ranges and higher melting points observed as the graphene content increased. Heat capacity peaks shifted from 6120 K to 6650 K, correlating with increased carbon concentration. Coordination number analysis revealed earlier structural breakdown in hBN compared to graphene, reflecting weaker B\u0026ndash;N bonding. Angular distribution studies and bond orientational order parameters confirmed that higher graphene content delays structural distortion and enhances thermal resilience. These results demonstrate the critical role of graphene in preserving structural integrity under extreme thermal conditions, offering valuable insights for the design of high-temperature 2D heterostructures.\u003c/p\u003e","manuscriptTitle":"Thermal Stability and Phase Transition Behavior in In-Plane Graphene/hBN/Graphene Heterostructures: A Molecular Dynamics Study","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2026-03-05 18:34:26","doi":"10.21203/rs.3.rs-8822403/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"94892433-4fc5-4e47-b4e6-3998c46c57fd","owner":[],"postedDate":"March 5th, 2026","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2026-03-09T01:09:23+00:00","versionOfRecord":[],"versionCreatedAt":"2026-03-05 18:34:26","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-8822403","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-8822403","identity":"rs-8822403","version":["v1"]},"buildId":"XKTyCvWXoU3ODBz1xrDgd","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}