Comparison and Assessment of Different Interatomic Potentials for Simulation of Silicon Carbide

Interatomic potentials play a crucial role in the molecular dynamics (MD) simulation of silicon carbide (SiC). However, the ability of interatomic potentials to accurately describe certain physical properties of SiC has yet to be confirmed, particularly for hexagonal SiC. In this study, the mechanical, thermal, and defect properties of four SiC structures (3C-, 2H-, 4H-, and 6H-SiC) have been calculated with multiple interatomic potentials using the MD method, and then compared with the results obtained from density functional theory and experiments to assess the descriptive capabilities of these interatomic potentials. The results indicate that the T05 potential is suitable for describing the elastic constant and modulus of SiC. Thermal calculations show that the Vashishta, environment-dependent interatomic potential (EDIP), and modified embedded atom method (MEAM) potentials effectively describe the vibrational properties of SiC, and the T90 potential provides a better description of the thermal conductivity of SiC. The EDIP potential has a significant advantage in describing point defect formation energy in hexagonal SiC, and the GW potential is suitable for describing vacancy migration in hexagonal SiC. Furthermore, the T90 and T94 potentials can effectively predict the surface energies of the three low-index surfaces of 3C-SiC, and the Vashishta potential exhibits excellent capabilities in describing stacking fault properties in SiC. This work will be helpful for selecting a potential for SiC simulations.

These interatomic potentials have different characteristics and applications in simulations.The Tersoff potential accurately describes the covalent bonds of C-Si [19], which are predominant in the bonding character of SiC.Therefore, the Tersoff potentials have been used to simulate mechanical phenomena in SiC, such as tension [20], shear [21], and impact simulations [22].The Vashishta potential is well-suited for accurately modeling the deformation, both bending and stretching, of ionic and covalent bonds in 3C-SiC [15].It is extensively used in simulations involving the impact behavior [23] and nanoindentation [24] of SiC.The EDIP potential can accurately describe the transition from covalent to metallic bonding, making it suitable for simulations of amorphous SiC phases [16].MEAM improves upon the embedded atom method by adding angular-dependent terms to better capture the directional bonding characteristics and has been used in crystal growth [17].The GW potential, similar in expression format to the Tersoff potential and accounting for stable defect configurations [18], is primarily used in cascade collisions [25].
The parameters of the interatomic potential are determined by fitting to certain physical properties.For instance, Vashishta potential fitted the cohesive energy, bulk modulus, and C 11 elastic constant of 3C-SiC [15], and GW potential has performed fittings for the cohesive energy, lattice constant, and bulk modulus [18].However, for physical properties not involved in the parameter fitting, the description of these properties by the interatomic potential remains to be validated.Therefore, it is necessary to assess the descriptive ability of the potentials before applying them.The results calculated using the potentials are usually compared with the experimental or density functional theory (DFT) results.
Researchers have predicted and compared some interatomic potentials of properties related to 3C-SiC [11,16,17,26], such as elastic constants and moduli.As determined, T05 provides good results of the elastic constants of 3C-SiC [11,16,27].Moreover, Vashishta potential shows agreement with experimental results, except for C 44 for 3C-SiC [27].In the field of thermal properties, most studies have focused on the thermal conductivity of 3C-SiC [28][29][30][31][32][33][34][35] and 4H-SiC [35].The research revealed that compared with T89 potential, the thermal conductivity predictions by GW potential closely match the experimental results [32].In addition, studies using interatomic potentials to calculate phonon dispersion in SiC have focused on 3C-SiC [19,35].The results indicate that Tersoff overestimates the optical branch of the dispersion curve [19,36], and compared with GW, MEAM, T05, and T89 potentials, T90 potential is better in describing phonon thermal transport [35].However, there were only few comparisons of the prediction with these interatomic potentials on the elastic mechanics, thermal conductivity, and phonon dispersion curves for hexagonal SiC.
Defect properties have a significant impact on the evolution of materials.For example, the formation energies of defects can help assess their stability, and the migration barriers contribute to understanding their migration mechanisms.Some interatomic potentials have been employed to predict and compare the formation energies of point defects in 3C-SiC [11,[16][17][18].The study demonstrated that the GW interatomic potential closely mirrors the ab initio calculations in describing C dumbbell interstitials and antisite defects [18].The T05 potential tends to underestimate the formation energy of vacancy defects [11], and EDIP potential can better describe C interstitial defects of 3C-SiC [16].As for point defect migration, compared with the T94 potential, the GW potential is more suitable for predicting carbon atom vacancy-interstitial recombination [37].For 4H-SiC, empirical potential has been employed to investigate the threshold displacement energies and displacement cascades using MD simulation [38].However, predictions of properties for defects in hexagonal SiC such as 4H-SiC are rarely mentioned.Regarding planar defects such as generalized stacking fault, Li et al. compared the descriptions of the generalized stacking fault energy (GSFE) using the Vashishta and T05 potentials and concluded that the Vashishta potential is more suitable for describing GSFE in 3C-SiC [27].Compared to T94, T05 exhibits a closer approximation to the DFT calculation in predicting the unstable GSFE of <112>{111} for 3C-SiC [39].Xue et al. calculated the GSFE of 4H-SiC using the Vashishta potential and suggested that the easy slip in the <11 20> direction is caused by its significantly low GSFE [24].The shape of GSFE determines the nature of dislocation cores [40], and thus the descriptions of GSFE in hexagonal SiC by other potentials deserve further explorations.
Motivated by the aforementioned research status, we present predictions of multiple interatomic potentials for the elastic, thermal, and defect properties of four SiC structures with MD methods.These predictions include elastic constants, moduli, Vickers hardness, Debye temperature, vibrational properties, thermal conductivity, point defect formation energies, point defect migration energy barriers, surface energy, and GSFE.In addition, based on DFT and experiments results, a comprehensive comparison of these interatomic potentials in describing the physical properties has been conducted.This evaluation aims to provide insights into the performance of different interatomic potentials in describing properties of SiC.This study provides guidance for the selection of interatomic potentials in SiC simulations.

Mechanical Properties Calculation
The mechanical characteristics are crucial criteria for assessing the quality of interatomic potential.Among these, elasticity is a vital mechanical property, serving as a criterion for mechanical stability, intrinsic ductility/brittleness, and fracture toughness assessment.To address this, we calculate the elastic constants (C ij ) of SiC.For the models of 3C-, 2H-, 4H-, and 6H-SiC, conventional unit cells were used as repeating units, with expansion performed using the following repeat patterns: 8 × 8 × 8, 8 × 8 × 8, 8 × 8 × 8, and 12 × 12 × 2, respectively.The essential parameters, including the bulk modulus (B), shear modulus (G), Young's modulus (Y), and Poisson's ratio (ν), were derived based on the elastic constants.Among these, B and G were obtained using the Hill average method [42].The bulk modulus B, Young's modulus Y, and shear modulus G are obtained using the following formulas [43].
The cubic structure (3C-SiC) is given by and the hexagonal structure (2H-, 4H-, and 6H-SiC) given by The Young's modulus Y and the Poisson's ratio ν are given by The Vickers hardness can be obtained from elastic properties.Because SiC is a ceramic material, Vickers hardness is calculated using the following formula [44]:

Thermal Properties Calculation
In SiC, transverse velocity (v t ), longitudinal velocity (v l ), and average velocity (v m ) are given by where ρ represents the density of the crystal, n is the number of atoms in the molecular formula, N A is Avogadro's number, k B is the Boltzmann constant, and M is the molecular weight.
For the calculation of phonon dispersion curves, the primitive cell, as shown in Figure 1, was employed as the repeating unit, with repeat patterns of 8 × 8 × 8 (3C), 8 × 8 × 8 (2H), 8 × 8 × 8 (4H), and 12 × 12 × 2 (6H), following the method proposed by Kong et al. [45].Periodic boundary conditions were applied in all three directions.The time step was set to 0.0002 ps.The microcanonical ensemble (NVE) was chosen, and the fix phonon [45] command was applied to the entire system during the calculation process.A time step of 0.002 ps was used, and the measurement of Green's function was conducted every 10 MD steps in the following 12,000 ps.After obtaining the dynamic matrix, PHANA [45] and PHONOPY [46] codes were used to calculate the phonon frequency.The phonon dispersion curves were also calculated by DFT with the VASP (version 6.1.0)[47].The generalized gradient approximation (GGA) in the Perdew-Burke-Enzerhof (PBE) [48] exchange-correlation energy functional wa and the cutoff energy was set to 600 eV.The energy convergence standard was se 10 −8 eV.The model was constructed based on protocells, as the repeating unit, a ployed expansion patterns of 3 × 3 × 3 (3C), 3 × 3 × 2 (2H), 3 × 3 × 1 (4H), and 3 × 3 × The 5 × 5 × 5 k-meshes were used in the Brillouin region.
The molar heat capacity C V and entropy S can be calculated by PHONOPY, and the calculation formula is as follows [49]: where ω qj is the frequency calculated by D(q) → e qj = ω 2 qj → e qj and D is the dynamical matrix.Thermal conductivities were calculated using equilibrium molecular dynamics (EMD) based on Green's function approach.The expression for thermal conductivity calculation can be derived from the Green's function as follows: where k B , T, V, and < J(0)J(t) > represent the Boltzmann constant, temperature, volume of the model, and the heat flux autocorrelation function, respectively.For 3C-SiC, thermal conductivities were averaged across the [100], [010], and [001] directions, while for hexagonal structures, thermal conductivities were calculated separately for both in-plane and out-of-plane directions.In the thermal conductivity calculations, the number of atoms in the models for 3C-SiC, 2H-SiC, 4H-SiC, and 6H-SiC were 1728, 1728, 2400, and 2400, respectively.Convergence can typically be achieved to approximately 1500 atoms when calculating thermal conductivity using the EMD method [50].During the calculations, we used a time step of 0.0002 ps.After relaxing for 106 steps (200 ps) in the NPT ensemble, followed by an additional 106 steps at least in the NVE ensemble, data collection was performed in the NVE ensemble.The results are the average of five independent simulations.Owing to the difficulty in achieving convergence for the heat flux autocorrelation function at 300 K, data collection was performed for 12,000 ps when the temperature was at T = 300 K.For temperatures T > 300 K, data collection was performed for 6000 ps.

Defect Properties Calculation
This study primarily focuses on the formation energies of point defects in SiC, including V C (C vacancy), V Si (Si vacancy), substitutional defect C Si (C occupying Si site), substitutional defect Si C (Si occupying C site), SiT C (Si occupying C tetrahedral site), CT Si (C occupying Si tetrahedral site), and eight types of dumbbell interstitials.The conjugate gradient (cg) method was used to minimize the energy.The timestep was set to 0.001 ps.The formation energies of a single vacancy

and substitutional defect E
A B f were given as follows: where A and B represent C or Si, respectively, and E(A) and E(B) are the energy per atom in the bulk phase (diamond structures of C and Si).E d and E perfect represent the energy of the system with defects and without defects, respectively.The migration energies of point defects were estimated by means of the nudged elastic band (NEB) method with 'quickmin' minimization style.The calculation of GSFE refers to the method of Vitek et al. [51].The model of 3C-SiC consists of 10 × 10 × 20 unit cells along X[110], Y[112], and Z[111] directions, respectively.The models of hexagonal SiC consist of 10 × 10 × 20 (2H), 10 × 10 × 10 (4H), and 10 × 10 × 6 (6H) unit cells along X[2110], Y[0110], and Z[0001] directions, respectively.The GSFE of 3C-SiC were calculated for displacements along the [110] and [112] directions in the (111) plane, and hexagonal SiC were calculated for displacements along the [2110] direction in the (0001) plane.The timestep was set to 0.005 ps.
To provide reference results for the FPs calculations of point defect formation energy, defect migration energy, and GSFE of hexagonal SiC, we employed the VASP software with using the GGA in the PBE parameterization.For the calculation of point defect formation energy, the energy cutoff was set to 500 eV [52].When calculating defect migration energy and GSFE, the energy cutoff was set to 600 eV.The force convergence criterion was set to 0.01 eV/Å.The method of climbing image nudge elastic band (CI-NEB) was used to calculate the defect migration energy [53].

Mechanical Properties
Before the examination of mechanical properties, a comprehensive understanding of the fundamental structural attributes of SiC under the influence of diverse interatomic potentials was established.Lattice constants and cohesive energies for the four SiC structures were computed, and the results are presented in Table 1.The cubic phase (3C-SiC) and the hexagonal phases (2H-, 4H-, and 6H-SiC) exhibit similar cohesive energy (E c ) values.Except for T94, other interatomic potentials describe the lattice constants of SiC well.The calculated bulk modulus B values for 3C-SiC, which closely match the DFT and experimental results, are illustrated in Figure 2.This is because most interatomic potentials are constructed using Young's modulus as a fitting parameter.The elastic constants of all interatomic potentials satisfy the following mechanical stability criteria [43]: C 11 > 0, C 11 − C 12 > 0, C 44 > 0, C 12 + 2C 12 > 0. The T05 potential exhibits excellent descriptive capabilities, showing strong consistency with the DFT results.In comparison, T89, T90, and T94 potential slightly overestimate C 44 and the shear modulus G.The GW potential does not perform as well as Tersoff potentials regarding describing elastic constants and moduli.This difference may be attributed to the fact that GW potential was developed primarily to describe the formation energy of point defects [16,25]  to describe the formation energy of point defects [16,25] rather than mechanical properties.The calculated values of C11 and C12 obtained through Vashishta, EDIP, and MEAM potentials show good agreement with experimental results.However, Vashishta, EDIP, and MEAM potentials underestimate C44, which consequently leads to a lower Young's modulus Y.Although 2H-SiC consists of a 100% hexagonal structure, which is fundamentally different from the 100% cubic structure of 3C-SiC, the calculated bulk modulus B for 2H-SiC closely agrees with the experimental and DFT results.This implies that the description of the bulk modulus B by SiC interatomic potentials remains largely unaffected by the Although 2H-SiC consists of a 100% hexagonal structure, which is fundamentally different from the 100% cubic structure of 3C-SiC, the calculated bulk modulus B for 2H-SiC closely agrees with the experimental and DFT results.This implies that the description of the bulk modulus B by SiC interatomic potentials remains largely unaffected by the differences in crystal structures.For 2H-SiC, the mechanical properties calculated using the Tersoff potential show excellent agreement with the experimental and DFT results.
In addition, the MEAM potential tends to overestimate C 12 and C 33 while significantly underestimating C 13 for 2H-SiC.
In the case of 4H-SiC, the moduli are in excellent agreement with the experimental values.Comparing the results for the three hexagonal phases (Figures 3-5), except for the MEAM and EDIP potentials, the elastic properties of 4H-SiC, 6H-SiC, and 2H-SiC are quite similar under the interactions of the other interatomic potentials, aligning well with the DFT results.By comparing the elastic constants of 2H-SiC, 4H-SiC (50% hexagonality), and 6H-SiC (33% hexagonality) obtained using the MEAM and EDIP potentials, it can be observed that, unlike other interatomic potentials, the elastic constants of these three phases are different.This indicates that the mechanical properties are influenced by the different structures when using these two interatomic potentials.differences in crystal structures.For 2H-SiC, the mechanical properties calculated using the Tersoff potential show excellent agreement with the experimental and DFT results.In addition, the MEAM potential tends to overestimate C12 and C33 while significantly underestimating C13 for 2H-SiC.
In the case of 4H-SiC, the moduli are in excellent agreement with the experimental values.Comparing the results for the three hexagonal phases (Figures 3-5), except for the MEAM and EDIP potentials, the elastic properties of 4H-SiC, 6H-SiC, and 2H-SiC are quite similar under the interactions of the other interatomic potentials, aligning well with the DFT results.By comparing the elastic constants of 2H-SiC, 4H-SiC (50% hexagonality), and 6H-SiC (33% hexagonality) obtained using the MEAM and EDIP potentials, it can be observed that, unlike other interatomic potentials, the elastic constants of these three phases are different.This indicates that the mechanical properties are influenced by the different structures when using these two interatomic potentials.To conduct a more in-depth assessment of these interatomic potentials, we calculated Poisson's ratio, B/G ratio, and Vickers hardness based on the elastic constants.Poisson's ratio can be used to predict the shear stability of materials and distinguish whether a material is brittle or ductile, with 0.33 as the boundary point [8].Table 2 shows that the Poisson's ratio obtained from the Tersoff, Vashishta, EDIP, and MEAM potentials indicates that 3C-SiC is a brittle material, and the same result holds for other phases (2H, 4H, and 6H) listed in Tables 3-5.The value of B/G can also characterize the toughness of materials, with 1.75 as the boundary to distinguish between brittle and ductile materials.Materials with a B/G ratio below 1.75 are considered brittle, whereas those above 1.75 are considered ductile [8].The B/G ratio of 3C-SiC is also listed in Table 2. Tersoff, Vashishta, and EDIP consistently suggest that 3C-SiC exhibits brittle behavior, which is consistent with the experimental findings.For 3C-SiC, the experimental values for Vickers hardness fall within the range of 26 ± 2 GPa [59].Tersoff and Vashishta potentials can reasonably predict the Vickers hardness of 3C-SiC, particularly for T05.However, for 2H-, 4H-, and 6H-SiC, Vashishta, GW, and MEAM potentials tend to underestimate the Vickers hardness.To conduct a more in-depth assessment of these interatomic potentials, we calculated Poisson's ratio, B/G ratio, and Vickers hardness based on the elastic constants.Poisson's ratio can be used to predict the shear stability of materials and distinguish whether a material is brittle or ductile, with 0.33 as the boundary point [8].Table 2 shows that the Poisson's ratio obtained from the Tersoff, Vashishta, EDIP, and MEAM potentials indicates that 3C-SiC is a brittle material, and the same result holds for other phases (2H, 4H, and 6H) listed in Tables 3-5.The value of B/G can also characterize the toughness of materials, with 1.75 as the boundary to distinguish between brittle and ductile materials.Materials with a B/G ratio below 1.75 are considered brittle, whereas those above 1.75 are considered ductile [8].The B/G ratio of 3C-SiC is also listed in Table 2. Tersoff, Vashishta, and EDIP consistently suggest that 3C-SiC exhibits brittle behavior, which is consistent with the experimental findings.For 3C-SiC, the experimental values for Vickers hardness fall within the range of 26 ± 2 GPa [59].Tersoff and Vashishta potentials can reasonably predict the Vickers hardness of 3C-SiC, particularly for T05.However, for 2H-, 4H-, and 6H-SiC, Vashishta, GW, and MEAM potentials tend to underestimate the Vickers hardness.Overall, in terms of the mechanical properties of SiC, including elastic characteristics, brittleness, and hardness, Tersoff potentials provide the most accurate descriptions, whether for the cubic or hexagonal phase.The Vashishta potential closely follows as a strong performer.Although the Vashishta potential relies solely on fitted data referencing cohesive energy, bulk modulus, and the C 11 elastic constant from experimental data [15], it is still suitable for describing the mechanical properties of hexagonal SiC.

Debye Temperature
As we know, the Debye temperature is a crucial thermal property of solid materials.It not only provides fundamental characteristics of vibrational spectra but also characterizes properties such as the coefficient of thermal expansion.The calculated sound velocities and Debye temperatures are presented in Table 6.Except for the slightly larger values obtained with T89, all Tersoff potentials, whether in terms of sound velocities or Debye temperatures, are close to the DFT results.Other interatomic potentials generally exhibit slightly smaller values, particularly for the GW potential.

Phonon Dispersion Curve
Phonon dispersion provides crucial insights into the structural and thermal properties of materials.For 3C-SiC, the primary path used is Γ-X-K-Γ-L [36], as show in Figure 6.Overall, none of the interatomic potentials show imaginary frequencies when describing phonon dispersion curves.Some of the outcomes of interatomic potentials exhibit slight variations compared with previous research results [35], which are likely attributed to differences in the selected model sizes.To further investigate the accuracy of the dispersion curves, Table 7 summarizes the phonon frequencies at high-symmetry points for 3C-SiC.Tersoff potentials, namely T05, T94, T90, and T89, provide good descriptions of the longitudinal acoustic (LA) branch.For example, the frequencies of the LA branch at the L point for T05, T94, T90, and T89 are 16.65 THz, 16.65 THz, 17.99 THz, and 17.61 THz, respectively, which agree well with DFT (16.65 THz) and the experimental results (16.95 THz).However, Tersoff potentials do not provide a very accurate description of the transverse acoustic (TA) branch.The frequency of the TA branch at the X point is approximately 27% higher than the experimental results.In addition, the optical branches of their phonon dispersion curves are not well-described.This is because some interatomic potentials, such as Tersoff potentials, only consider the interaction between covalent bonds when proposed.However, SiC also contains 12% ionic bonds; thus, interatomic potentials like Tersoff are unable to accurately describe the optical branches, and the TA branch tends to be overestimated [19].
Previous research found that T90 potential can describe the acoustic branch of phonons well [35], but the acoustic branches obtained using Vashishta, EDIP, and MEAM potentials agree well with the DFT curves.Among them, the frequencies of the optical branches from Vashishta potential are in the same frequency range as those from DFT.This is because both covalent interactions and ionic interactions are taken into account when constructing the Vashishta potentials [15].This contributes to a relatively accurate calculation of phonon dispersion using the Vashishta potential.and MEAM potentials agree well with the DFT curves.Among them, the frequencies of the optical branches from Vashishta potential are in the same frequency range as those from DFT.This is because both covalent interactions and ionic interactions are taken into account when constructing the Vashishta potentials [15].This contributes to a relatively accurate calculation of phonon dispersion using the Vashishta potential.The phonon dispersion curves for 2H-SiC are shown in Figure 7. Similar to the case of 3C-SiC, none of the interatomic potentials can accurately describe the optical branches.The Tersoff potential overestimates the frequencies of the optical branches and provides a relatively good description of the acoustic branches, particularly the LA branches.For the TA branches, particularly along the L-H path, the frequencies are approximately 18% higher than the DFT curves.GW fails to describe the acoustic branch frequencies of 2H-SiC, and the frequencies at various high-symmetry points are generally lower.The results from the Vashishta potential agree well with the DFT results, particularly for the acoustic branches, and the optical branches also fall within the same frequency range as the DFT optical branches.The results from the EDIP and MEAM potentials are similar, with a good description of the acoustic branches but overestimated frequencies for the optical branches.The results of other hexagonal phases (4H and 6H) are similar to those of 2H-SiC, as shown in Figures 8 and 9.The phonon dispersion curves for 2H-SiC are shown in Figure 7. Similar to the case of 3C-SiC, none of the interatomic potentials can accurately describe the optical branches.The Tersoff potential overestimates the frequencies of the optical branches and provides a relatively good description of the acoustic branches, particularly the LA branches.For the TA branches, particularly along the L-H path, the frequencies are approximately 18% higher than the DFT curves.GW fails to describe the acoustic branch frequencies of 2H-SiC, and the frequencies at various high-symmetry points are generally lower.The results from the Vashishta potential agree well with the DFT results, particularly for the acoustic branches, and the optical branches also fall within the same frequency range as the DFT optical branches.The results from the EDIP and MEAM potentials are similar, with a good description of the acoustic branches but overestimated frequencies for the optical branches.
The results of other hexagonal phases (4H and 6H) are similar to those of 2H-SiC, as shown in Figures 8 and 9.

Heat Capacity and Entropy
The behavior of the molar heat capacity C V of 3C-SiC with temperature is depicted in Figure 10a.The trends of the curves for various interatomic potentials are relatively consistent with the DFT results.Each curve exhibits a rapid increase as the temperature rises, and once the temperature surpasses the Debye temperature, they gradually converge to the same value (Dulong-Petit limit).To better illustrate the curve changes, we have zoomed in on the low-temperature region.For 3C-SiC, in the low-temperature range (T < 100 K), Vashishta, EDIP, and MEAM potentials align closely with the heat capacity curve derived from DFT calculations.This phenomenon may be attributed to the fact that at low temperatures, only the low-frequency acoustic branches are excited, which means the heat capacity is primarily influenced by low-frequency vibrations.Combining Figure 6e,g,h, it can be inferred that the acoustic branches of the Vashishta, EDIP, and MEAM potentials match well with the DFT results, resulting in a good agreement in their C V curves with DFT.The acoustic branch frequencies of Tersoff potentials are generally higher than those from DFT.The acoustic branch frequencies of GW potential are lower than the results of DFT, leading to discrepancies in the heat capacity calculated using this potential compared with DFT results.However, at T > 300 K, the GW results exhibit a better agreement with DFT results.

Heat Capacity and Entropy
The behavior of the molar heat capacity CV of 3C-SiC with temperature is depicted i Figure 10a.The trends of the curves for various interatomic potentials are relatively con sistent with the DFT results.Each curve exhibits a rapid increase as the temperature rises and once the temperature surpasses the Debye temperature, they gradually converge t the same value (Dulong-Petit limit).To better illustrate the curve changes, we hav zoomed in on the low-temperature region.For 3C-SiC, in the low-temperature range (T 100 K), Vashishta, EDIP, and MEAM potentials align closely with the heat capacity curv derived from DFT calculations.This phenomenon may be attributed to the fact that at low temperatures, only the low-frequency acoustic branches are excited, which means the hea capacity is primarily influenced by low-frequency vibrations.Combining Figure 6e and Figure 6g,h, it can be inferred that the acoustic branches of the Vashishta, EDIP, and MEAM potentials match well with the DFT results, resulting in a good agreement in thei CV curves with DFT.The acoustic branch frequencies of Tersoff potentials are generall higher than those from DFT.The acoustic branch frequencies of GW potential are lowe than the results of DFT, leading to discrepancies in the heat capacity calculated using thi potential compared with DFT results.However, at T > 300 K, the GW results exhibit better agreement with DFT results.
Similarly, the heat capacity CV of hexagonal SiC increases with increasing tempera tures, as illustrated in Figure 10b-d.Different SiC structures converge to distinct limits in the high-temperature region.When the temperature T < 100 K (low-temperature range) the results from the Vashishta, EDIP, and MEAM potentials are relatively close to the hea capacity curve obtained from DFT calculations.When T > 300 K, the results from the GW potential gradually approach the molar heat capacity curve from DFT.Furthermore, we also calculated the vibrational entropy (S) as a function of temper ature, and the results are shown in Figure 11.For four SiC structures, all potentials exhib consistent trends.In the low-temperature region of the S-T curves, the results of Vash ishta, EDIP, and MEAM potentials closely match the DFT calculation.In particular, th Vashishta potential, when describing S for 2H-SiC, still closely overlaps with the referenc curve even in the high-temperature region.Similarly, the heat capacity C V of hexagonal SiC increases with increasing temperatures, as illustrated in Figure 10b-d.Different SiC structures converge to distinct limits in the high-temperature region.When the temperature T < 100 K (low-temperature range), the results from the Vashishta, EDIP, and MEAM potentials are relatively close to the heat capacity curve obtained from DFT calculations.When T > 300 K, the results from the GW potential gradually approach the molar heat capacity curve from DFT.
Furthermore, we also calculated the vibrational entropy (S) as a function of temperature, and the results are shown in Figure 11.For four SiC structures, all potentials exhibit consistent trends.In the low-temperature region of the S-T curves, the results of Vashishta, EDIP, and MEAM potentials closely match the DFT calculation.In particular, the Vashishta potential, when describing S for 2H-SiC, still closely overlaps with the reference curve even in the high-temperature region.

Thermal Conductivity
The temperature dependence of the thermal conductivity of 3C-SiC is shown in Figure 12.The predictions of the interatomic potentials, except for the T89 and T94 potentials, decrease significantly with increasing temperature, which is in agreement with the results of previous studies [34,64].The decrease in thermal conductivity with increasing temperature is attributed to the reduction in the average phonon mean free path.The thermal conductivity obtained with GW and T89 potentials is in agreement with previous studies [32].In addition, the thermal conductivity calculated using T94 potential closely matches that obtained by previous researchers who employed T94 potential with a Ziegler-Biersack-Littmark short-range correction [34,65].The results from the MEAM and T90 potentials fall between the DFT [66] and experimental values [67].When the temperature exceeds 300 K, the results from all interatomic potentials deviate from the DFT and experimental results.Overall, in comparison, T90 and GW potentials show better agreement with the reference values.T94 and T89 potentials tend to significantly underestimate the thermal conductivity of SiC.

Thermal Conductivity
The temperature dependence of the thermal conductivity of 3C-SiC is shown in Figure 12.The predictions of the interatomic potentials, except for the T89 and T94 potentials, decrease significantly with increasing temperature, which is in agreement with the results of previous studies [34,64].The decrease in thermal conductivity with increasing temperature is attributed to the reduction in the average phonon mean free path.The thermal conductivity obtained with GW and T89 potentials is in agreement with previous studies [32].In addition, the thermal conductivity calculated using T94 potential closely matches that obtained by previous researchers who employed T94 potential with a Ziegler-Biersack-Littmark short-range correction [34,65].The results from the MEAM and T90 potentials fall between the DFT [66] and experimental values [67].When the temperature exceeds 300 K, the results from all interatomic potentials deviate from the DFT and experimental results.Overall, in comparison, T90 and GW potentials show better agreement with the reference values.T94 and T89 potentials tend to significantly underestimate the thermal conductivity of SiC.[68] and DFT [66].
For hexagonal SiC, thermal conductivities in both in-plane and out-of-plane directions were calculated separately, as shown in Figure 13.In comparison, thermal conductivity obtained using T90 potential appears to be closer to the FPs results for the two directions.However, other interatomic potentials tend to significantly underestimate the thermal conductivity of hexagonal SiC, whether in the in-plane or out-of-plane direction.For hexagonal SiC, thermal conductivities in both in-plane and out-of-plane directions were calculated separately, as shown in Figure 13.In comparison, thermal conductivity obtained using T90 potential appears to be closer to the FPs results for the two directions.However, other interatomic potentials tend to significantly underestimate the thermal conductivity of hexagonal SiC, whether in the in-plane or out-of-plane direction.

Point Defect Formation Energy
The results of the formation energy of defects are presented in Table 8.Because of structural instability observed during DFT calculations for Si-C + <110> (the symbol "+" represents the interstitial atom), data on this structure are not listed in Table 8.In addition, when applying the Vashishta potential to compute the cohesive energy of diamond C phase and diamond Si phase, both C and Si exhibited excessively low cohesive energies, approaching zero, making it impossible to calculate the formation energy of point defects in SiC.Thus, we have opted not to include the defect data associated with the Vashishta potential.To describe vacancies, T94, T90, and GW potentials seem to have an advantage.The GW potential is better suited for describing substitutional defects.It should be noted that when the GW potential was proposed, other formulas were employed to calculate the formation energy of point defects.To ensure a meaningful comparison, we used Equations ( 20)-( 22) to recalculate the formation energy.The GW potential does not obtain a good interstitial defect formation energy because the GW potential cannot describe the energies of pure silicon and pure carbon.Regarding interstitial defects located in tetrahedral sites, only T05 can achieve the best results for CT Si and CT C .In the case of dumbbell interstitial defects, T89, T90, and MEAM potentials can describe C-C + <100> well.T05 and T89 potentials can predict the formation energy of Si-C + <100>.Furthermore, T90 and GW potentials provide accurate formation energies for Si-Si + <100> and C-Si + <100>.Overall, it appears that none of these interatomic potentials can effectively predict the formation energies of all dumbbell interstitials.However, based on the description of more easily generated point defects, T94 should be the optimal choice.Compared with other interaction potentials, the formation energies of point defects obtained using GW potential are generally lower.This leads to a higher number of point defects in cascade collisions obtained using the GW potential than those obtained with the T94 potential [37].The results of formation energies for point defects in hexagonal SiC are presented in Figure 14 and also listed in Tables S1-S3.The structural reference for interstitial defects is taken from a previous study [72].The DFT results show that the difference in structure has no significant effect on the point defect formation energy in 2H and 6H-SiC.Each has its advantages in describing these defects.The T90, T94, and GW potentials perform well in describing V C , V Si , and Si C .The MEAM and T89 potentials provide formation energies for Si C that are in excellent agreement with DFT results.However, from the perspective of the energy trends, only the results obtained from the T05 and EDIP potentials align with those from DFT.Among these, EDIP potential appears to provide a better description of the defect formation energies in hexagonal SiC. for SiC that are in excellent agreement with DFT results.However, from the perspective of the energy trends, only the results obtained from the T05 and EDIP potentials align with those from DFT.Among these, EDIP potential appears to provide a better description of the defect formation energies in hexagonal SiC.
(a) (b) (c) for 2H-, 4H-, and 6H-SiC, respectively.The DFT results for 4H-SiC are referenced from a previous study [72].The h represents the hexagonal environment, and k represents the cubic environment.

Point Defect Migration Barrier
The nudged elastic band (NEB) method helps in exploring the diffusion mechanisms of vacancies and interstitials [73].Through the evaluation of migration pathways and energy barriers associated with point defects, including vacancies and interstitials, the precision of interatomic potentials in characterizing atomic mobility and energy variations can be corroborated.
Here, we calculate the migration of vacancies in 3C-SiC along three paths, namely VC → VC, VSi → VSi, and VC-CSi → VSi, respectively.The DFT results are also calculated, as shown in Figure 15.The migration barriers for these three paths obtained from DFT calculations are 3.64 eV, 3.51 eV, and 2.41 eV, respectively, which is consistent with the values reported by Bockstedte et al. [74].For interatomic potentials, in the VC → VC pathway, only the GW potential can obtain relatively good results, followed by the T89 potential.In the VSi → VSi pathway, GW and EDIP potentials yield relatively good results for migration barriers.Interstitial migration outcomes are illustrated in Figure 16, with (a) and (b) denoting the migration for carbon interstitial (IC) and silicon interstitial (ISi), respectively.For IC, the most favorable migration pathway is Csp <100> → Csp <100> [75].For ISi, the two most

Point Defect Migration Barrier
The nudged elastic band (NEB) method helps in exploring the diffusion mechanisms of vacancies and interstitials [73].Through the evaluation of migration pathways and energy barriers associated with point defects, including vacancies and interstitials, the precision of interatomic potentials in characterizing atomic mobility and energy variations can be corroborated.
Here, we calculate the migration of vacancies in 3C-SiC along three paths, namely V C → V C , V Si → V Si , and V C -C Si → V Si , respectively.The DFT results are also calculated, as shown in Figure 15.The migration barriers for these three paths obtained from DFT calculations are 3.64 eV, 3.51 eV, and 2.41 eV, respectively, which is consistent with the values reported by Bockstedte et al. [74].For interatomic potentials, in the V C → V C pathway, only the GW potential can obtain relatively good results, followed by the T89 potential.In the V Si → V Si pathway, GW and EDIP potentials yield relatively good results for migration barriers.
for SiC that are in excellent agreement with DFT results.However, from the perspective of the energy trends, only the results obtained from the T05 and EDIP potentials align with those from DFT.Among these, EDIP potential appears to provide a better description of the defect formation energies in hexagonal SiC.

Point Defect Migration Barrier
The nudged elastic band (NEB) method helps in exploring the diffusion mechanisms of vacancies and interstitials [73].Through the evaluation of migration pathways and energy barriers associated with point defects, including vacancies and interstitials, the precision of interatomic potentials in characterizing atomic mobility and energy variations can be corroborated.
Here, we calculate the migration of vacancies in 3C-SiC along three paths, namely VC → VC, VSi → VSi, and VC-CSi → VSi, respectively.The DFT results are also calculated, as shown in Figure 15.The migration barriers for these three paths obtained from DFT calculations are 3.64 eV, 3.51 eV, and 2.41 eV, respectively, which is consistent with the values reported by Bockstedte et al. [74].For interatomic potentials, in the VC → VC pathway, only the GW potential can obtain relatively good results, followed by the T89 potential.In the VSi → VSi pathway, GW and EDIP potentials yield relatively good results for migration barriers.Interstitial migration outcomes are illustrated in Figure 16, with (a) and (b) denoting the migration for carbon interstitial (IC) and silicon interstitial (ISi), respectively.For IC, the most favorable migration pathway is Csp <100> → Csp <100> [75].For ISi, the two most Interstitial migration outcomes are illustrated in Figure 16, with (a) and (b) denoting the migration for carbon interstitial (I C ) and silicon interstitial (I Si ), respectively.For I C , the most favorable migration pathway is C sp <100> → C sp <100> [75].For I Si , the two most favorable interstitial sites are the carbon tetrahedral site and the <110> dumbbell interstitial [74].In addition, Si sp <110> → Si TC is an essential constituent unit for the long-range migration of neutral silicon interstitial [76].We calculated the energy barriers for both C sp <100> → C sp <100> and Si sp <110> → Si TC migration pathways.In the C sp <100> → C sp <100> pathway, the energy barrier calculated through DFT is found to be 0.89 eV, which is very close to the value of 0.88 eV in a previous study [75].The results show that the GW potential describes the migration of I C well, particularly in terms of the barrier height.The results for I Si migration from Si sp <110> to Si TC also show that the empirical potentials differ significantly from the DFT calculations, particularly in the description of the saddle point barrier.During migration, some of these interatomic potentials yield negative values, indicating that the structures of the initial and final states of migration are in metastable states.
favorable interstitial sites are the carbon tetrahedral site and the <110> dumbbell interstitial [74].In addition, Sisp <110> → SiTC is an essential constituent unit for the long-range migration of neutral silicon interstitial [76].We calculated the energy barriers for both Csp <100> → Csp <100> and Sisp <110> → SiTC migration pathways.In the Csp <100> → Csp <100> pathway, the energy barrier calculated through DFT is found to be 0.89 eV, which is very close to the value of 0.88 eV in a previous study [75].The results show that the GW potential describes the migration of IC well, particularly in terms of the barrier height.The results for ISi migration from Sisp <110> to SiTC also show that the empirical potentials differ significantly from the DFT calculations, particularly in the description of the saddle point barrier.During migration, some of these interatomic potentials yield negative values, indicating that the structures of the initial and final states of migration are in metastable states.In 2H-, 4H-, and 6H-SiC, we primarily investigate the VC and VSi in a hexagonal structural environment.The DFT reference results for 2H-, 4H-, and 6H-SiC are calculated by the current study.For the migration barriers of VC and VSi in 4H-SiC, the DFT results are 3.30 eV and 3.90 eV, respectively, consistent with the values obtained in previous research using the GGA and HSE hybrid functional, which were 3.31 eV for VC and 4.32 eV for VSi [77].The results in Figure 17 indicate that the empirical potential exhibits good symmetry for VC migration, but only GW potential provides accurate results for the energy barrier calculation in hexagonal SiC.When describing the migration behavior of VSi in hexagonal SiC, both the GW and EDIP potentials demonstrate excellent descriptive capabilities.The main features are that the migration path exhibits good symmetry, and the energy barrier of the saddle point matches the DFT value.In contrast, other interatomic potentials tend to overestimate or underestimate the migration energy barrier.Overall, GW potential is well-suited for describing vacancy migration in hexagonal SiC.In 2H-, 4H-, and 6H-SiC, we primarily investigate the V C and V Si in a hexagonal structural environment.The DFT reference results for 2H-, 4H-, and 6H-SiC are calculated by the current study.For the migration barriers of V C and V Si in 4H-SiC, the DFT results are 3.30 eV and 3.90 eV, respectively, consistent with the values obtained in previous research using the GGA and HSE hybrid functional, which were 3.31 eV for V C and 4.32 eV for V Si [77].The results in Figure 17 indicate that the empirical potential exhibits good symmetry for V C migration, but only GW potential provides accurate results for the energy barrier calculation in hexagonal SiC.When describing the migration behavior of V Si in hexagonal SiC, both the GW and EDIP potentials demonstrate excellent descriptive capabilities.The main features are that the migration path exhibits good symmetry, and the energy barrier of the saddle point matches the DFT value.In contrast, other interatomic potentials tend to overestimate or underestimate the migration energy barrier.Overall, GW potential is well-suited for describing vacancy migration in hexagonal SiC.

Surface Energy and GSFE
Surface is also a very important defect in materials [78].Furthermore, in nanoscale indentation and similar processes, surface energy can also serve as a measure of deformation on low-index crystallographic planes [27].The calculated surface energy results for 3C-SiC are shown in Figure 18.It can be seen that the surface energies of (100) and (110) surfaces calculated by T94 and T90 potentials are in excellent agreement with DFT results.In the case of the (111) surface, the result obtained through the MEAM potential is the closest to DFT, followed by those from T90 and T94 potentials.Overall, T90 and T94 potentials exhibit good capabilities in describing the formation energies of the three low-index surfaces in 3C-SiC.

Surface Energy and GSFE
Surface is also a very important defect in materials [78].Furthermore, in nanoscale indentation and similar processes, surface energy can also serve as a measure of deformation on low-index crystallographic planes [27].The calculated surface energy results for 3C-SiC are shown in Figure 18.It can be seen that the surface energies of (100) and (110) surfaces calculated by T94 and T90 potentials are in excellent agreement with DFT results.In the case of the (111) surface, the result obtained through the MEAM potential is the closest to DFT, followed by those from T90 and T94 potentials.Overall, T90 and T94 potentials exhibit good capabilities in describing the formation energies of the three lowindex surfaces in 3C-SiC.The GSFE is one of planar defects that control fracture behavior and determine dislocation properties [51,78].Furthermore, it represents the energy required for lattice movement along a specific plane and can serve as an indicator for evaluating dislocation slip [24].This study evaluates these interatomic potentials for stacking fault energy by calculating the GSFE of SiC.Considering the slip systems of 3C-SiC, calculations were performed on the (111) plane along both the [110] and [112] directions.The results in Figure 19 show that the displacements are normalized.It can be observed that the empirical potentials tend to overestimate the GSFE in the [110] direction.In the [112]   The GSFE is one of planar defects that control fracture behavior and determine location properties [51,78].Furthermore, it represents the energy required for latt movement along a specific plane and can serve as an indicator for evaluating disloca slip [24].This study evaluates these interatomic potentials for stacking fault energy calculating the GSFE of SiC.Considering the slip systems of 3C-SiC, calculations w performed on the (111) plane along both the [110 ] and [112] directions.The result Figure 19 show that the displacements are normalized.It can be observed that the em ical potentials tend to overestimate the GSFE in the [110 ] direction.In the [112] direct the T05 and Vashishta potentials exhibit energy values consistent with those obtai from DFT, with Vashishta potential being notably consistent.In hexagonal SiC, when nanoindentation occurs on the (0001) basal plane, slip is served along the [1 12 0 ] direction [79].Thus, for hexagonal SiC, we primarily calcula the variation of GSFE with displacement for the slip system along the a-axis (   The GSFE is one of planar defects that control fracture behavior and determine dislocation properties [51,78].Furthermore, it represents the energy required for lattice movement along a specific plane and can serve as an indicator for evaluating dislocation slip [24].This study evaluates these interatomic potentials for stacking fault energy by calculating the GSFE of SiC.Considering the slip systems of 3C-SiC, calculations were performed on the (111) plane along both the [110 ] and [112] directions.The results in Figure 19 show that the displacements are normalized.It can be observed that the empirical potentials tend to overestimate the GSFE in the [110 ] direction.In the [112]  In hexagonal SiC, when nanoindentation occurs on the (0001) basal plane, slip is observed along the [1 12 0 ] direction [79].Thus, for hexagonal SiC, we primarily calculated the variation of GSFE with displacement for the slip system along the a-axis (  In hexagonal SiC, when nanoindentation occurs on the (0001) basal plane, slip is observed along the [1120] direction [79].Thus, for hexagonal SiC, we primarily calculated the variation of GSFE with displacement for the slip system along the a-axis ([1120] direction) during rigid body movement.The DFT results indicate that the GSFE for 2H-, 4H-, and 6H-SiC gradually increases, as shown in Figure 20.The stacking fault energy of 4H-SiC calculated using the Vashishta potential is in close agreement with the calculations reported by previous researchers [24].The results obtained for the three phases using empirical potentials exhibit remarkable similarity.Overall, the Vashishta potential is consistent in the shape of GSFE curves.The structure of the dislocation core is closely associated with the shape of GSFE [27,40].Therefore, the Vashishta potential is often used for the nanoindentation simulation of hexagonal SiC.The GSFE is closely related to the mechanism of dislocation motion in materials.For instance, when calculated using the T89 potential, the GSFE results for 4H-SiC and 6H-SiC are highly similar.Consequently, under the interactions of this interatomic potential, the stress-strain curves obtained from nanoindentation simulations for 4H-SiC and 6H-SiC would exhibit a close resemblance [79].
for the nanoindentation simulation of hexagonal SiC.The GSFE is closely related to the mechanism of dislocation motion in materials.For instance, when calculated using the T89 potential, the GSFE results for 4H-SiC and 6H-SiC are highly similar.Consequently, under the interactions of this interatomic potential, the stress-strain curves obtained from nanoindentation simulations for 4H-SiC and 6H-SiC would exhibit a close resemblance [79].

Conclusions
To assess the ability of multiple interatomic potentials to describe the properties of four common SiC structures, elastic constants, moduli, vibrational properties, thermal conductivities, point defect formation energies, point defect migration energies, surface energies, and GSFE were calculated using MD methods.The results show that the T05 potential agrees with both the DFT and experimental results in terms of elastic constant, modulus, brittle/toughness, and Vickers hardness, which implies that this interatomic potential is more suitable for mechanical simulations of SiC.The better prediction of the elastic constants allows the T05 potential to obtain more accurate Debye temperatures.The Vashishta, EDIP, and MEAM potentials effectively describe the lattice vibrational properties of SiC.At the macroscopic level, the thermal conductivity of SiC obtained from the T90 potential agrees well with the FPs and the experimental results.In terms of point defects of 3C-SiC, all interatomic potentials can only describe the formation energies of some point defects.The EDIP potential provides more accurate point defect formation energies in hexagonal SiC.For point defect migration, the GW potential is more advantageous.Thus, the GW potential appears to be more suitable for the cascaded collision simulation of hexagonal SiC.The T90 and T94 potentials describe the surface energy of the three lowindex surfaces of 3C-SiC well.The results of the GSFE with the Vashishta potential are consistent with the DFT, which indicates that this potential is also suitable for the mechanical simulation of the impact and polishing of hexagonal SiC.This study is helpful for selecting interatomic potentials in SiC simulation.

Figure 2 .
Figure 2. Elastic constants and moduli for 3C-SiC calculated using different interatomic potentials.

Figure 2 .
Figure 2. Elastic constants and moduli for 3C-SiC calculated using different interatomic potentials.

Figure 3 .
Figure 3. Elastic constants and moduli for 2H-SiC calculated using different interatomic potentials.

Figure 3 .
Figure 3. Elastic constants and moduli for 2H-SiC calculated using different interatomic potentials.

Figure 3 .
Figure 3. Elastic constants and moduli for 2H-SiC calculated using different interatomic potentials.

Figure 4 .
Figure 4. Elastic constants and moduli for 4H-SiC calculated using different interatomic potentials.

Figure 5 .
Figure 5. Elastic constants and moduli for 6H-SiC calculated using different interatomic potentials.

Figure 12 .
Figure 12.Thermal conductivity for 3C-SiC.The pink and purple lines represent the results of the experiment[68] and DFT[66].

Figure 12 .
Figure 12.Thermal conductivity for 3C-SiC.The pink and purple lines represent the results of the experiment[68] and DFT[66].

Figure 13 .
Figure 13.Temperature-dependent thermal conductivity for 2H-, 4H-, and 6H-SiC.The in-plane and out-of-plane thermal conductivities of 2H-SiC are shown in (a) and (b), respectively.The in-plane and out-of-plane thermal conductivities of 4H-SiC are shown in (c) and (d), respectively.The inplane and out-of-plane thermal conductivity of 6H-SiC are shown in (e) and (f), respectively.The experimental and FPs results are from a previous study [69].

Figure 13 .
Figure 13.Temperature-dependent thermal conductivity for 2H-, 4H-, and 6H-SiC.The in-plane and out-of-plane thermal conductivities of 2H-SiC are shown in (a) and (b), respectively.The in-plane and out-of-plane thermal conductivities of 4H-SiC are shown in (c) and (d), respectively.The in-plane and out-of-plane thermal conductivity of 6H-SiC are shown in (e) and (f), respectively.The experimental and FPs results are from a previous study [69].

Figure 14 .
Figure 14.Defect formation energies for different point defects.(a), (b), and (c) represent the results for 2H-, 4H-, and 6H-SiC, respectively.The DFT results for 4H-SiC are referenced from a previous study [72].The h represents the hexagonal environment, and k represents the cubic environment.

Figure 14 .
Figure 14.Defect formation energies for different point defects.(a-c) represent the results for 2H-, 4H-, and 6H-SiC, respectively.The DFT results for 4H-SiC are referenced from a previous study [72].The h represents the hexagonal environment, and k represents the cubic environment.

Figure 14 .
Figure 14.Defect formation energies for different point defects.(a), (b), and (c) represent the results for 2H-, 4H-, and 6H-SiC, respectively.The DFT results for 4H-SiC are referenced from a previous study [72].The h represents the hexagonal environment, and k represents the cubic environment.

Figure 15 .
Figure 15.Vacancy migration barrier for 3C-SiC.(a-c) depict the migration energy barriers for V C → V C , V Si → V Si , and V C -C Si → V Si in 3C-SiC, respectively.

Figure 16 .
Figure 16.Interstitial migration energy barrier for 3C-SiC.(a,b) depict the migration of I C along the C sp <100> → C sp <100> path and I Si along the Si sp <110> → Si TC path, respectively, with the DFT results in (b) referenced from a previous study[76].The "sp" represents split interstitial.

Figure 17 .
Figure 17.Vacancy migration barrier for hexagonal SiC.(a,b), (c,d), and (e,f) illustrate the migration energy barriers for V C → V C and V Si → V Si in 2H-, 4H-, and 6H-SiC, respectively.
direction, the T05 and Vashishta potentials exhibit energy values consistent with those obtained from DFT, with Vashishta potential being notably consistent.

Figure 18 .
Figure 18.Surface energies of three low-index surfaces for 3C-SiC.The DFT, T05, and Vashishta are taken from a previous study [27].

Figure 19 .
Figure 19.Generalized stacking fault energy for 3C-SiC.(a,b) show the GSF energies for movem along the [110] and [112] directions on the (111) plane, with data from DFT, T05, and Vashi taken from a previous study [27].
[1 12 0 ]  rection) during rigid body movement.The DFT results indicate that the GSFE for 2H-, , and 6H-SiC gradually increases, as shown in Figure20.The stacking fault energy of

Figure 18 .
Figure 18.Surface energies of three low-index surfaces for 3C-SiC.The DFT, T05, and Vashishta data are taken from a previous study [27].

Figure 18 .
Figure 18.Surface energies of three low-index surfaces for 3C-SiC.The DFT, T05, and Vashishta data are taken from a previous study [27].

Figure 19 .
Figure 19.Generalized stacking fault energy for 3C-SiC.(a,b) show the GSF energies for movement along the [110] and [112] directions on the (111) plane, with data from DFT, T05, and Vashishta taken from a previous study [27].
[1 12 0 ] direction) during rigid body movement.The DFT results indicate that the GSFE for 2H-, 4H-, and 6H-SiC gradually increases, as shown in Figure20.The stacking fault energy of 4H-SiC calculated using the Vashishta potential is in close agreement with the calculations reported by previous researchers[24].The results obtained for the three phases using empirical potentials exhibit remarkable similarity.Overall, the Vashishta potential is

Figure 19 .
Figure 19.Generalized stacking fault energy for 3C-SiC.(a,b) show the GSF energies for movement along the [110] and [112] directions on the (111) plane, with data from DFT, T05, and Vashishta taken from a previous study [27].

Table 1 .
Lattice parameters a and c (Å) as well as cohesive energies E c (eV) calculated using interatomic potentials.The experimental and DFT results are also listed.
rather than mechanical properties.The calculated values of C 11 and C 12 obtained through Vashishta, EDIP, and MEAM potentials show good agreement with experimental results.However, Vashishta, EDIP, and MEAM potentials underestimate C 44 , which consequently leads to a lower Young's modulus Y.

Table 2 .
Calculated elastic constants C ij (GPa), bulk modulus B (GPa), shear modulus G (GPa), Young's modulus Y (GPa), Poisson's ratio ν, the rate of B/G, and Vickers hardness H V (GPa) for 3C-SiC.The results of the experiment and DFT calculation are listed as well.

Table 3 .
Calculated elastic constants C ij (GPa), bulk modulus B (GPa), shear modulus G (GPa), Young's modulus Y (GPa), Poisson's ratio ν, the rate of B/G, and Vickers hardness H V (GPa) for 2H-SiC.The results of the experiment and theoretical calculation are listed as well.

Table 4 .
Calculated elastic constants C ij (GPa), bulk modulus B (GPa), shear modulus G (GPa), Young's modulus Y (GPa), Poisson's ratio ν, the rate of B/G, and Vickers hardness H V (GPa) for 4H-SiC.The results of the experiment and theoretical calculation are listed as well.

Table 5 .
Calculated elastic constants C ij (GPa), bulk modulus B (GPa), shear modulus G (GPa), Young's modulus Y (GPa), Poisson's ratio ν, the rate of B/G, and Vickers hardness H V (GPa) for 6H-SiC.The results of experiment and theoretical calculation are listed as well.

Table 7 .
Phonon frequencies ν (THz) in high-symmetry points of the Brillouin zone calculated using different potentials and results of DFT and experiments for 3C-SiC.

Table 7 .
Phonon frequencies ν (THz) in high-symmetry points of the Brillouin zone calculated using different potentials and results of DFT and experiments for 3C-SiC.

Table 8 .
Formation energy (eV) of different point defects for 3C-SiC, calculated by DFT and empirical potentials.The "-" represents an unstable structure.