Stability and Elasticity of Quasi-Hexagonal Fullerene Monolayer from First-Principles Study

: As a newly synthesized two-dimensional carbon material, the stability study of mono-layer fullerene networks or quasi-hexagonal phase fullerenes (qhp-C 60 ) is timely desirable. We have investigated the stabilities of qhp-C 60 , including thermal, structural, mechanical, and thermodynamic stabilities, as well as the bonding characteristics, ductility, and mechanical properties, via ﬁrst-principles calculations. The results show that qhp-C 60 is energetically, mechanically, and thermodynamically stable. The thermodynamic stability of qhp-C 60 at 300 K and 600 K is veriﬁed. The bonding characteristics of qhp-C 60 are analyzed from the bond length, and it has sp 2 and sp 3 hybridization. The Pugh ratio ( B/G ) and Poisson’s ratio ( v ) indicate similar ductility with graphite and graphene. We also found that qhp-C 60 has the lowest hardness and the anisotropy of the material. In addition, the electronic characteristics, including electron localization function (ELF), crystal orbital Hamiltonian population (COHP), and density of states (DOS) at different temperatures, are analyzed to verify the thermal stability of the material. Our results might be helpful in the material design of qhp-C 60 -related applications.


Introduction
Carbon-based two-dimensional (2D) materials have attracted extensive attention because of their unique and rich properties, functions, and applications.As the paradigm of 2D materials, graphene has a perfect sp 2 hybridized carbon structure in a two-dimensional plane.In addition, it has excellent electronic, thermal, optical, and mechanical properties [1][2][3].Graphene has been successfully used to fabricate field-effect transistors (FETs) [4][5][6], transparent conducting electrodes [7][8][9], supercapacitors [10,11], and lithium-ion batteries [12].The successful preparation of graphene [13] has shown that it is possible to exfoliate stable, single-atom, or single-polyhedral-thick 2D materials from van der Waals solids.These materials can exhibit unique and fascinating physical properties.Frindt et al. showed that layered van der Waals materials, such as layered metal dichalcogenides, could be mechanically and chemically exfoliated into few and single layers [14,15].Novoselov et al. isolated individual monolayers from various layered materials, including h-BN and dihydroxides (e.g., MoS 2 and NbSe 2 ) [16].Moreover, fullerenes have interesting electronic and magnetic properties.Superconductivity was also found in fullerene compounds [17][18][19][20].Thus, scientists want to prepare 2D fullerene, further studying its intrinsic properties.Hou et al. prepared a large-size monocrystalline 2D carbon material, a monolayer of quasi-hexagonal crystalline phase fullerene (C 60 ), by interlayer bonding cleavage [21].This 2D material has good thermodynamic stability and high crystallinity.In particular, this monolayer polymer C 60 has a moderate band gap of about 1.6 eV, solving the problem of zero band gap in graphene.
In practical applications, the stability of materials is an essential factor affecting the reliability and lifetime of devices.The structure and properties of 2D materials are susceptible to environmental conditions [22][23][24][25][26][27][28], so the stability of 2D materials is critical.With the development of computer science, first-principles calculation has shown great advantages in the stability of materials.So far, several studies on the stability of carbon materials have been carried out.Sen et al. studied the thermal stability of sp 3 -bonded 2D carbon allotropes using ab initio molecular dynamics [29].This study shows that only low-energy structures can be stable under environmental conditions.Li et al. reported that monolayer h-BN shows strong thermal stability when heated in the air [30].Wen et al. investigated the mechanical stability of simple cubic carbon [31].Simple cubic carbon is mechanically stable above 400 GPa.And an excellent mechanical and thermodynamic stability of carbon can be found above 2700 GPa.Lu et al. discovered that metastable carbon allotrope could withstand temperatures up to 1000 K by first-principles study, and its in-plane stiffness is greater than that of graphene [32].He et al. found that both sp 2 -diamond and cubic-graphite are dynamically stable by phonon spectroscopy.The band gap of sp 2 -diamond is 1.66 eV, and it is a semiconductor [33].Luo et al. found that biphenylene is dynamically stable and has a melting point of up to 4500 K by first-principles calculation [34].Mortazavi et al. showed that the qhp-C 60 has desirable dynamical stability through predicted phonon dispersion [35].However, they explored the unilateral stability of qhp-C 60 only.There is still a lack of systematic research into the stability of qhp-C 60 .In addition, the analysis of the electronic structure is effective for obtaining the interactions between atoms and the properties of the material.Mei et al. found that octagonal-C is metallic and the C-C bond is predominantly covalent by electronic analysis [36].Alborznia et al. showed that in the presence of applied pressure, the band gap of T-Carbon becomes smaller, but it is still a semiconductor material [37].Wang et al. found that sp 3 hybridization facilitates the formation of a symmetrical lattice structure to ensure the high hardness of the carbon material [38].Thus, electronic structure analysis provides insight into the potential applications of the material.
In this paper, we have assessed the stability of the qhp-C 60 structure (structural, mechanical, and thermodynamic stability), bonding characteristics, ductility, and mechanical properties using first-principles calculations.The stability and bonding changes of the materials at different temperatures were analyzed by electron localization function (ELF) and crystal orbital Hamiltonian population (COHP).In addition, the contribution of orbitals at different temperatures is analyzed by using density of states (DOS).The result shows that quasi-hexagonal fullerene monolayer offers an interesting platform for potential application in 2D electronic devices.

Materials and Methods
Density functional theory (DFT) calculations are carried out using the Vienna ab initio simulation package (VASP) in the study [39].The electron interaction with ions was replaced by the projector augmented wave (PAW) method, and exchange-correlation potentials between electrons were parameterized by the generalized gradient approximation (GGA) with Perdew-Burke-Ernzerhof (PBE) function [40,41].The qhp-C 60 primitive cell containing 60 carbon atoms was used.The conjugate-gradient algorithm was used to relax all atomic positions and cell size to achieve the minimum total energy of the systems until the average force on each atom is less than 0.01 eV/atom.The plane-wave cutoff energy was specified as 520 eV for all models.The 4 × 4 × 1 k-points mesh obtained by the Monkhorst-Pack method was employed to sample the Brillouin zone.The total energy convergence was set to 10 −6 eV/atom.A 15 Å (1 Å = 0.1 nm) vacuum layer was set along the z-axis direction to avoid the influence of the interaction force between atomic layers.The energy-strain relationship was used to calculate elastic constants.For the quadratic symmetry 2D materials, the orthogonal lattice transformed by selecting the [−110] [−1 −10] [001] direction to redefine the coordinate axis was selected and the energy-strain method was used to calculate elastic constants.The structure's stability was ensured by the energy minimization during the structure relaxation via the conjugate gradient method [42].AIMD simulations were performed to verify the thermodynamic stability of qhp-C 60 at finite temperatures.In doing this, a 2 × 2 × 1 supercell of qhp-C 60 was simulated within the NVT systems by using the Nose-Hoover thermostat for 10 ps with a time step of 0.0001 ps at 300 K and 600 K, and a 1 × 1 × 1 k-point mesh.The crystallographic results of qhp-C 60 are shown in Figure 1a.The structure parameters after relaxation are listed in Table 1, along with the calculated structures of other carbon materials (graphite, graphene, and diamond).The calculated values of our obtained structure parameters are in agreement with those of the references, demonstrating the accuracy of our calculations.
Crystals 2023, 13, x FOR PEER REVIEW 3 of 13 energy-strain relationship was used to calculate elastic constants.For the quadratic symmetry 2D materials, the orthogonal lattice transformed by selecting the [−110] [−1 −10] [001] direction to redefine the coordinate axis was selected and the energy-strain method was used to calculate elastic constants.The structure's stability was ensured by the energy minimization during the structure relaxation via the conjugate gradient method [42].AIMD simulations were performed to verify the thermodynamic stability of qhp-C60 at finite temperatures.In doing this, a 2 × 2 × 1 supercell of qhp-C60 was simulated within the NVT systems by using the Nose-Hoover thermostat for 10 ps with a time step of 0.0001 ps at 300 K and 600 K, and a 1 × 1 × 1 k-point mesh.

Results
3.1.Crystal Structure, Thermal, and Structural Stability The crystallographic results of qhp-C60 are shown in Figure 1a.The structure parameters after relaxation are listed in Table 1, along with the calculated structures of other carbon materials (graphite, graphene, and diamond).The calculated values of our obtained structure parameters are in agreement with those of the references, demonstrating the accuracy of our calculations.The results in Figure 1 and Table 1 show dA = 1.64 Å, dB = 1.61 Å, dC = 1.41-1.45Å, dD = 1.44-1.45Å, dE = 1.42 Å, and dF = 1.55 Å. Liu et al. showed that the length of C≡C is the shortest (1.28 Å) and that of C-C is the longest (1.55 Å) under the same system [46].Therefore, by comparing the bond lengths, we find that qhp-C60 has two different carbon  The results in Figure 1 and Table 1 show 42 Å, and d F = 1.55 Å. Liu et al. showed that the length of C≡C is the shortest (1.28 Å) and that of C-C is the longest (1.55 Å) under the same system [46].Therefore, by comparing the bond lengths, we find that qhp-C 60 has two different carbon networks with sp 2 and sp 3 hybridizations.This may lead to qhp-C 60 having different properties from other carbon materials.
In the diamond crystal, the four valence electrons of each carbon atom were sp 3 hybridized to form four identical atomic orbitals, forming covalent bonds (σ bonds) with the four most adjacent carbon atoms [47].In addition, the carbon atoms of graphene were bonded to each other in sp 2 hybridization and formed a sizable π-conjugated system.Kaciulis pointed out that π-type bonds determined the electronic properties, such as electrical conductance and optical gap, whereas the mechanical hardness is defined by the percentage of σ bonds.Therefore, the determination of the sp 2 /sp 3 ratio is decisive for characterizing all carbon allotropes.The study has shown that an appropriate sp 2 /sp 3 (about 1.0-1.25)ratio will strengthen the field emission properties [48].Low-dimensional structure materials with excellent field emission have broad application prospects in the vacuum microelectronic device.It is known that carbon-based materials such as diamond, diamond-like carbon, and carbon nanotubes.According to the results of the bond length, the value of sp 2 /sp 3 in the qhp-C 60 structure is approximately 1.4, which contributes to its application in electronic devices.

Formation and Binding Energy
To verify the structural stability of qhp-C 60 , the formation and binding energy of the qhp-C 60 are calculated using the following equations.
The enthalpy of formation is defined as the energy released (or absorbed) by the system during the formation of the compound from its elementary substance.The calculation formula is Binding energy is the energy released by an atom from forming a compound in a free state, and is calculated as: E tot is the total energy of the qhp-C 60 , E atom is the energy of each carbon atom in a free state, and E soild is the energy of a single carbon atom in the graphite.For stable structures, both the binding energy and the enthalpy of formation are negative.Generally, the smaller the value is, the more stable the system is.As shown in Table 1, the results show that the formation energy of graphene (0.05 eV/atom), diamond (0.27 eV/atom), and qhp-C 60 (0.45 eV/atom) is greater than zero, indicating that graphite is the most stable of the four carbon materials, and the conversion of graphite to other three carbon materials is a heatabsorbing process.The binding energy of qhp-C 60 is −8.75 eV/atom, which is close to that of graphite (−9.20 eV/atom), graphene (−9.15 eV/atom) and diamond (−8.93 eV/atom), respectively, indicating that qhp-C 60 has excellent structural stability.

Elasticity Constants
Elastic constants characterize the amount of material elasticity and reflect the ability of a crystal to resist elastic deformation.The energy-strain method is used to find the elastic constants [49].Elastic constants of diamond are first calculated to verify the feasibility of the method.The elastic constants, bulk modulus, and shear modulus are shown in Table 2 and compared with the values in the literature.The results show that the elastic constants of the diamond obtained by the energy-strain method are consistent with the results reported in the literature, so the calculation method is available.a is the unit of elastic modulus of graphite and diamond.b is the unit of elastic modulus of graphene and qhp-C 60 .
As shown in Table 2, It is worth noting that our calculation results are consistent with Zhao's [53].For two-dimensional materials, we assess their mechanical stability by the Born-Huang criterion.The two-dimensional materials need to satisfy C 11 > 0, C 11 C 22 -C 2 12 > 0, and C 66 > 0 [54], where C ij is the elastic constants.C 11 and C 22 are the ability of the crystal to resist axial strain along the [100] and [010], respectively.C 66 is the ability to resist shear strain on the (001), and C 12 is the modulus along the [110] that resists shear deformation on the (110).As the results in Table 2 show, the elastic parameters of qhp-C 60 2 also shows the bulk modulus, shear modulus, and Poisson's ratio of the four carbon materials.The results show that qhp-C 60 has the lowest bulk modulus, shear modulus, Young's modulus, and Poisson's ratio.

Ductility
The ductility of a crystalline could be predicted via the Pugh ratio, which is the ratio between the bulk and the shear modulus, B/G [55,56].Pugh ratio is used empirically to predict brittle or ductile behavior of materials.In general, ductility is positively correlated with the B/G value.A larger value implies a more ductile behavior under mechanical loading.The B/G and v of the four carbon materials are compared in Figure 2. Diamond has great hardness and brittleness, so it has smaller B/G and v than other three carbon materials.The B/G and v of qhp-C 60 are very close to those of graphite and graphene, which shows that qhp-C 60 has similar ductility with them.a is the unit of elastic modulus of graphite and diamond.b is the unit of elastic modulus of graphene and qhp-C60.

C11
As shown in Table 2, It is worth noting that our calculation results are consistent with Zhao's [53].For two-dimensional materials, we assess their mechanical stability by the Born-Huang criterion.The two-dimensional materials need to satisfy C11 > 0, C 11 C 22 -C 12 2 > 0, and C66 > 0 [54], where Cij is the elastic constants.C11 and C22 are the ability of the crystal to resist axial strain along the [100] and [010], respectively.C66 is the ability to resist shear strain on the (001), and C12 is the modulus along the [110] that resists shear deformation on the (110).As the results in Table 2 show, the elastic parameters of qhp-C60 are C11 = 140.20 N/m, C22 = 192.62N/m, C12 = 23.91N/m, and C66 = 63.54N/m, respectively, which obviously satisfy the Born-Huang criterion for two-dimensional materials.Thus, qhp-C60 has good mechanical stability.Table 2 also shows the bulk modulus, shear modulus, and Poisson's ratio of the four carbon materials.The results show that qhp-C60 has the lowest bulk modulus, shear modulus, Young's modulus, and Poisson's ratio.

Ductility
The ductility of a crystalline could be predicted via the Pugh ratio, which is the ratio between the bulk and the shear modulus, B/G [55,56].Pugh ratio is used empirically to predict brittle or ductile behavior of materials.In general, ductility is positively correlated with the B/G value.A larger value implies a more ductile behavior under mechanical loading.The B/G and v of the four carbon materials are compared in Figure 2. Diamond has great hardness and brittleness, so it has smaller B/G and v than other three carbon materials.The B/G and v of qhp-C60 are very close to those of graphite and graphene, which shows that qhp-C60 has similar ductility with them.

Hardness
Hardness is a crucial mechanical parameter in the mechanical properties of materials [57].There are many empirical models for calculating the hardness (Hv) of materials.The

Hardness
Hardness is a crucial mechanical parameter in the mechanical properties of materials [57].There are many empirical models for calculating the hardness (H v ) of materials.The selection of a suitable hardness model is determined based on the chemical bonding state of the material.Chen proposed a simple prediction model based on the elastic modulus of the material [58]: where k is the Pratt's modulus ratio, k = G/B.An improved model based on the above model was proposed by Tian, as shown in the following equation [59]: Jiang explained that Young's modulus (E) could be calculated using E = 15.76Hjiang [60].In addition, we use two other formulas to calculate hardness.
As shown in Table 3, the qhp-C 60 has the lowest hardness compared to other carbon materials.a is the unit of elastic modulus of graphite and diamond.b is the unit of elastic modulus of graphene and qhp-C 60 .

Anisotropy of the Material
The study of mechanical properties can bring a new understanding of the application of the material.On this basis, mechanical anisotropy is also crucial to the practical use of the material, and the generation of microcracks is closely related to anisotropy [61].Young's modulus (E) and Poisson's ratio (v) can be derived from the elastic constant combinations, as is shown in Figure 3. Young's modulus of the qhp-C60 increased from a minimum at 0 • to a maximum at 90 • .Similarly, we note that v exhibits strong anisotropic.
Crystals 2023, 13, x FOR PEER REVIEW 6 of 13 selection of a suitable hardness model is determined based on the chemical bonding state of the material.Chen proposed a simple prediction model based on the elastic modulus of the material [58]: where k is the Pratt's modulus ratio, k = G/B.An improved model based on the above model was proposed by Tian, as shown in the following equation [59]: Jiang explained that Young's modulus (E) could be calculated using E = 15.76Hjiang [60].In addition, we use two other formulas to calculate hardness.
As shown in Table 3, the qhp-C60 has the lowest hardness compared to other carbon materials.a is the unit of elastic modulus of graphite and diamond.b is the unit of elastic modulus of graphene and qhp-C60.

Anisotropy of the Material
The study of mechanical properties can bring a new understanding of the application of the material.On this basis, mechanical anisotropy is also crucial to the practical use of the material, and the generation of microcracks is closely related to anisotropy [61].Young's modulus (E) and Poisson's ratio (v) can be derived from the elastic constant combinations, as is shown in Figure 3. Young's modulus of the qhp-C60 increased from a minimum at 0° to a maximum at 90°.Similarly, we note that v exhibits strong anisotropic.

Thermodynamic Stability
To further validate the thermodynamic stability of qhp-C 60 via ab initio molecular dynamics (AIMD) at finite temperatures [62,63].AIMD simulations were conducted for a 2 × 2 × 1 supercell of qhp-C 60 at the temperature of 300 K and 600 K, as shown in Figure 4a,b, which is drawn by taking a point for each 0.1 ps (100 points in total).To further validate the thermodynamic stability of qhp-C60 via ab initio molecular dynamics (AIMD) at finite temperatures [62,63].AIMD simulations were conducted for a 2 × 2 × 1 supercell of qhp-C60 at the temperature of 300 K and 600 K, as shown in Figure 4a,b, which is drawn by taking a point for each 0.1 ps (100 points in total).The AIMD results illustrated that the total energy and temperature of the system converge with fluctuations after 2 ps relaxation.The small fluctuation of the temperature and energy after convergence at 300 K reflects that the system is well equilibrated and qhp-C60 is thermodynamically stable at ambient conditions.The thermal dynamic equilibrium at 600 K indicates stability at high temperature.
Moreover, as shown in Figure 4c,d, although the single-layer qhp-C60 structure simulated by AIMD has some special bond breaks and changes, it is almost identical to its initial configuration.Our results display that qhp-C60 is thermodynamically stable at both 300 K and 600 K.It suggests the high-temperature applications of the qhp-C60 monolayerbased materials and devices.

Electronic Localization Function Analysis
In quantum chemistry, the electron localization function (ELF) is a function of the possibility of finding electrons in the neighborhood space of a given reference electron with the same spin.It is usually used to analyze the spatial localization of electron charge.In order to further investigate the stability of qhp-C60 at room temperature and high The AIMD results illustrated that the total energy and temperature of the system converge with fluctuations after 2 ps relaxation.The small fluctuation of the temperature and energy after convergence at 300 K reflects that the system is well equilibrated and qhp-C 60 is thermodynamically stable at ambient conditions.The thermal dynamic equilibrium at 600 K indicates stability at high temperature.
Moreover, as shown in Figure 4c,d, although the single-layer qhp-C 60 structure simulated by AIMD has some special bond breaks and changes, it is almost identical to its initial configuration.Our results display that qhp-C 60 is thermodynamically stable at both 300 K and 600 K.It suggests the high-temperature applications of the qhp-C 60 monolayer-based materials and devices.

Electronic Localization Function Analysis
In quantum chemistry, the electron localization function (ELF) is a function of the possibility of finding electrons in the neighborhood space of a given reference electron with the same spin.It is usually used to analyze the spatial localization of electron charge.In order to further investigate the stability of qhp-C 60 at room temperature and high temperature.The electron localization function (ELF) of qhp-C 60 after the structural optimization at different temperatures is shown in Figure 5. Here, the ELF in the range of 0 to 1 is a space function describing the atomic bonding properties of structures, and the values of 0, 0.5, and 1 correspond to vacuum, homogeneous electron gas, and fully localized states, respectively.temperature.The electron localization function (ELF) of qhp-C60 after the structural optimization at different temperatures is shown in Figure 5. Here, the ELF in the range of 0 to 1 is a space function describing the atomic bonding properties of structures, and the values of 0, 0.5, and 1 correspond to vacuum, homogeneous electron gas, and fully localized states, respectively.As shown in Figure 5a, at a temperature of 0 K, the optimized qhp-C60 was found to have electron localization at the center of the bond, retaining its 2D network structure of covalent bonds, indicating that the covalent bonds between carbon atom pairs dominate.After that, as the temperature rises to room temperature, the ELF changes very little, and the charge localization of C0-C1 and C2-C3 bonds is weakened to a certain extent.At the same time, when the temperature rises to 600 K, the ELF value between C0-C1 and C2-C3 bonds becomes 0. On the other hand, the length of the C0-C1 bond in these three configurations is 1.61 Å, 1.68 Å, and 2.39 Å, while the C2-C3 bond length is 1.61 Å, 1.64 Å, and 2.51 Å, respectively.As the temperature increases, the distance between the C atom and the C atom becomes larger, the bond energy is smaller, and the bond strength also decreases, so the longer the bond, the easier it is to break.It means that the bond is broken at 600 K.However, the structure of qhp-C60 basically remains the original structure.

COHP and ICOHP Analysis
COHP analysis can describe the bonding characteristics between atoms and adjacent atoms, as shown in Figure 6.In general, the positive value of −COHP indicates the bond interaction, and the negative value indicates the anti-bond interaction.As the temperature increases, the interaction of carbon-carbon bonds has been decreasing and is more pronounced at 600 K.As shown in Figure 5a, at a temperature of 0 K, the optimized qhp-C 60 was found to have electron localization at the center of the bond, retaining its 2D network structure of covalent bonds, indicating that the covalent bonds between carbon atom pairs dominate.After that, as the temperature rises to room temperature, the ELF changes very little, and the charge localization of C 0 -C 1 and C 2 -C 3 bonds is weakened to a certain extent.At the same time, when the temperature rises to 600 K, the ELF value between C 0 -C 1 and C 2 -C 3 bonds becomes 0. On the other hand, the length of the C 0 -C 1 bond in these three configurations is 1.61 Å, 1.68 Å, and 2.39 Å, while the C 2 -C 3 bond length is 1.61 Å, 1.64 Å, and 2.51 Å, respectively.As the temperature increases, the distance between the C atom and the C atom becomes larger, the bond energy is smaller, and the bond strength also decreases, so the longer the bond, the easier it is to break.It means that the bond is broken at 600 K.However, the structure of qhp-C 60 basically remains the original structure.

COHP and ICOHP Analysis
COHP analysis can describe the bonding characteristics between atoms and adjacent atoms, as shown in Figure 6.In general, the positive value of −COHP indicates the bond interaction, and the negative value indicates the anti-bond interaction.As the temperature increases, the interaction of carbon-carbon bonds has been decreasing and is more pronounced at 600 K.The ICOHP analysis of orbital bonding strength under the configuration at a time of 10 ps and a temperature of 0 K, 300 K, and 600 K is shown in Figure 7.The integral of crystal orbital Hamilton populations (ICOHP) is obtained by integrating the Hamiltonian population below the Fermi level and can be used to measure the strength of covalent bonds in crystals.As a valuable measure of the covalent bond strength, the lower the value, the stronger the chemical bond, and vice versa [64,65].As shown, we found that the bond lengths of C0-C1 and C2-C3 bonds increased slightly from 0 K to 300 K, and the ICOHP value also increased slightly.Therefore, the bond strength did not change much at room temperature.At this time, no bond breakage and structural changes occurred.As the temperature increases from 300 K to 600 K, the bond length and ICOHP increase dramatically.C0-C1 and C2-C3 are broken.Our first-principles study is self-consistent, which is the same as the ELF analysis.

DOS Analysis
In order to better understand the differences in orbital contributions at different temperatures.We study the projected density of states (DOS) at different temperatures.The results are shown in Figure 8.The lower valence band below the Fermi level and the above part of the conduction band are mainly contributed by the s-orbital.The outermost electron valence layer of the p-orbital is more concentrated at the Fermi level and is also the main contributor to the bonding.The contributions of the three p-orbit components (px, py, pz) are almost the same.The ICOHP analysis of orbital bonding strength under the configuration at a time of 10 ps and a temperature of 0 K, 300 K, and 600 K is shown in Figure 7.The integral of crystal orbital Hamilton populations (ICOHP) is obtained by integrating the Hamiltonian population below the Fermi level and can be used to measure the strength of covalent bonds in crystals.As a valuable measure of the covalent bond strength, the lower the value, the stronger the chemical bond, and vice versa [64,65].As shown, we found that the bond lengths of C 0 -C 1 and C 2 -C 3 bonds increased slightly from 0 K to 300 K, and the ICOHP value also increased slightly.Therefore, the bond strength did not change much at room temperature.At this time, no bond breakage and structural changes occurred.As the temperature increases from 300 K to 600 K, the bond length and ICOHP increase dramatically.C 0 -C 1 and C 2 -C 3 are broken.Our first-principles study is self-consistent, which is the same as the ELF analysis.The ICOHP analysis of orbital bonding strength under the configuration at a time of 10 ps and a temperature of 0 K, 300 K, and 600 K is shown in Figure 7.The integral of crystal orbital Hamilton populations (ICOHP) is obtained by integrating the Hamiltonian population below the Fermi level and can be used to measure the strength of covalent bonds in crystals.As a valuable measure of the covalent bond strength, the lower the value, the stronger the chemical bond, and vice versa [64,65].As shown, we found that the bond lengths of C0-C1 and C2-C3 bonds increased slightly from 0 K to 300 K, and the ICOHP value also increased slightly.Therefore, the bond strength did not change much at room temperature.At this time, no bond breakage and structural changes occurred.As the temperature increases from 300 K to 600 K, the bond length and ICOHP increase dramatically.C0-C1 and C2-C3 are broken.Our first-principles study is self-consistent, which is the same as the ELF analysis.

DOS Analysis
In order to better understand the differences in orbital contributions at different temperatures.We study the projected density of states (DOS) at different temperatures.The results are shown in Figure 8.The lower valence band below the Fermi level and the above part of the conduction band are mainly contributed by the s-orbital.The outermost electron valence layer of the p-orbital is more concentrated at the Fermi level and is also the main contributor to the bonding.The contributions of the three p-orbit components (px, py, pz) are almost the same.Furthermore, we can see that qhp-C60 has a small bandgap, and its value is approximately 0.7 eV.It is well known that the standard DFT exchange-correlation function greatly underestimates the band gap.For insulators, this is due to the discontinuous existence of energy to electron derivatives.The difference between the calculated and the experimental band gap can reach 0.5 eV [66].The experimental and calculated values of the band gap of qhp-C60 are 1.6 eV and 1.2 eV.Based on these facts, the calculated band gap of qhp-C60 (0.7 eV) is in good agreement with their results.It is worth noting that the band gap values at 300 K and 600 K are 0.4 eV and 0.5 eV, respectively, with a slight decrease, indicating that the temperature rise does not seriously affect the band gap of the Furthermore, we can see that qhp-C 60 has a small bandgap, and its value is approximately 0.7 eV.It is well known that the standard DFT exchange-correlation function greatly underestimates the band gap.For insulators, this is due to the discontinuous existence of energy to electron derivatives.The difference between the calculated and the experimental band gap can reach 0.5 eV [66].The experimental and calculated values of the band gap of qhp-C 60 are 1.6 eV and 1.2 eV.Based on these facts, the calculated band gap of qhp-C 60 (0.7 eV) is in good agreement with their results.It is worth noting that the band gap values at 300 K and 600 K are 0.4 eV and 0.5 eV, respectively, with a slight decrease, indicating that the temperature rise does not seriously affect the band gap of the structure.The phenomenon provides guidance for the preparation of devices serving at hightemperature conditions.

Conclusions
We have investigated the stability, bonding characteristics, and mechanical properties of qhp-C 60 by means of first-principles calculation.The formation and binding energy are calculated to verify the structural stability of qhp-C 60 .The Born-Huang criterion is satisfied and qhp-C 60 has good mechanical stability.The thermodynamics stability is further ensured via ab initio molecular dynamics (AIMD) at 300 K and 600 K.The analysis result of the bond length demonstrates the coexistence of two different carbon networks with sp 2 and sp 3 hybridization in qhp-C 60 .The value of sp 2 /sp 3 is approximately 1.4, implying that it can be helpful in electronic devices.Both Pugh ratio and Poisson's ratio show that qhp-C 60 has approximate mechanical properties with graphite and graphene.In addition, qhp-C 60 has the lowest hardness, and the values of Young's modulus and Poisson's ratio in different directions show that the material has pronounced anisotropy.The results of ELF analysis show that the bond breaks when the temperature rises to 600 K, which is also verified by COHP analysis, but the structure at this temperature still maintains the basic initial structure.The analysis of the electronic density of states shows that the contribution of the orbital does not change with the increase in temperature, but the band gap decreases.

Figure 1 .
Figure 1.Crystal structure (a-c) qhp-C 60 , (d) graphite, and (e) diamond.Here, d A -d D , d E , and d F are the bond length of qhp-C60, graphite, graphene, and diamond.

Figure 2 .
Figure 2. B/G (a) and Poisson's ratio (b) for the four carbon materials.

Figure 2 .
Figure 2. B/G (a) and Poisson's ratio (b) for the four carbon materials.

Figure 4 .
Figure 4. Total energy (left) and temperature (right) of qhp-C60 at (a) 300 K and (b) 600 K.The snapshot of the configuration at time of 10 ps and temperature of 300 K (c) and 600 K (d).

Figure 4 .
Figure 4. Total energy (left) and temperature (right) of qhp-C 60 at (a) 300 K and (b) 600 K.The snapshot of the configuration at time of 10 ps and temperature of 300 K (c) and 600 K (d).

Figure 7 .
Figure 7.The ICOHP value of C 0 -C 1 and C 2 -C 3 bonds of at 0 K, 300 K, and 600 K.

3. 6 .
DOS AnalysisIn order to better understand the differences in orbital contributions at different temperatures.We study the projected density of states (DOS) at different temperatures.The results are shown in Figure8.The lower valence band below the Fermi level and the above part of the conduction band are mainly contributed by the s-orbital.The outermost electron valence layer of the p-orbital is more concentrated at the Fermi level and is also the main contributor to the bonding.The contributions of the three p-orbit components (p x , p y , p z ) are almost the same.

Funding:
S.T. gratefully acknowledges the financial support of the National Natural Science Foundation of China (Grant Nos.52175293 and 51774083).Q.P. and X.C. would like to acknowledge the Shenzhen Science and Technology Program (Grant No. KQTD20200820113045081). Q.P. would like to acknowledge the support provided by the National Natural Science Foundation of China (Grant No. 12272378) and the LiYing Program of the Institute of Mechanics, Chinese Academy of Sciences (Grant No. E1Z1011001).J.J. acknowledges the Fundamental Research Funds for the Central Universities (Grant Nos.N2007002 and N2007011), 111 Project (Grant No. B20029).

Table 2 .
Elastic constants, bulk modulus, shear modulus, Young's modulus, and Poisson's ratio of the four carbon materials.

Table 3 .
The hardness of the four carbon materials.

Table 3 .
The hardness of the four carbon materials.