Density functional theory-based surrogate kinetic models for heterogeneous reactions of hydrocarbon intermediates on silicon carbide

Abstract

The increasing demand for high-performance materials in advanced technologies highlights the importance of achieving a fundamental understanding and potential control of silicon carbide (SiC) deposition processes. However, existing models often lack sufficient theoretical detail, relying heavily on empirical data and offering limited predictive capability. In particular, the complex surface chemistry governing SiC growth remains poorly understood. This study addresses these challenges by employing density functional theory (DFT) to investigate key heterogeneous reactions involving hydrocarbon intermediates on SiC surfaces, including dehydrogenation, hydrogenation, and carbon deposition. Transition state searches were conducted to identify reaction pathways and energy barriers. While first-principles calculations offer high accuracy, they are computationally intensive. To extend the utility of these first-principles results, vibrational analyses were performed using phonon-based statistical thermochemistry to compute temperature-dependent reaction rates which were used to develop Arrhenius-type surrogate kinetic models. The resulting framework provides a more rigorous, physically grounded basis for integrating atomistic insights into continuum-scale modeling, ultimately enabling improved prediction and optimization of SiC film growth in high-performance material systems.

---

1. Introduction

Ceramic matrix composites (CMCs) are advanced materials composed of ceramic fibers embedded within a ceramic matrix.¹ This results in a material with enhanced strength, toughness, and high-temperature resistance compared to traditional ceramics. Silicon carbide (SiC) is one of the most widely used CMCs in many advanced technologies.² Among the various fabrication techniques, the chemical vapor infiltration (CVI) process stands out as a promising method for producing CMCs. However, challenges related to process optimization and control, and cost remain significant.³ Process modeling plays a crucial role in addressing these issues as it achieves fundamental understanding and thereby enables exploration of optimal processing conditions cost-effectively. By reducing the reliance on extensive experimental trials, modeling facilitates improved process efficiency while ensuring or even enhancing material quality.⁴

Chemical Vapor Infiltration (CVI) was first proposed and developed by Bickerdike in 1962 with the goal of increasing the density of porous carbon materials.⁵ Despite decades of advancement, the reactive gas–surface interactions that drive material deposition in CVI remain poorly understood.⁶ This is mostly due to the complex physico-chemical phenomena involved, which occur across multiple length and time scales.⁷ Existing models of surface chemistry in CVI exhibit significant inaccuracies, primarily due to an incomplete understanding of the coupling between gas and surface reactions.⁸

Methyltrichlorosilane (CH₃SiCl₃), commonly referred to as MTS, is the most widely used precursor for SiC production via CVI.¹⁰ Existing studies on the decomposition of MTS and its subsequent surface deposition during CVI often overlook complex chemical interactions and primarily rely on simplified reaction pathways.¹¹ These models, typically calibrated through experimental data, tend to be case-specific and lack generalizability, frequently leading to inaccurate predictions of deposition rates and material properties. To improve accuracy, it is essential to identify and characterize elementary surface chemical reactions involving a representative set of intermediate species.^11–15 MTS decomposition produces various Si- and C-bearing species that undergo surface reactions, with some pathways involving higher energy barriers and others being more energetically favorable.¹⁶ Through the stepwise reaction mechanism for SiC deposition, it is important to identify and model the rate-controlling steps accurately.¹⁷ In particular, an imbalance in the deposition of silicon and carbon atoms is a major limiting factor in achieving stoichiometric SiC, which requires a 1 : 1 Si-to-C ratio. Several studies have identified the carbon deposition path as the primary limiting reaction set in CVI.¹⁸ Resolving the rate-limiting reaction pathways between hydrocarbon species and SiC is key to gaining critical insight into controlling and enhancing both the deposition rate and quality.

To address these limitations, a carefully validated, bottom-up modeling approach is needed which begins with atomistic, first-principles investigations of surface chemical reaction mechanism pathways, particularly suited for capturing the intricacies and complexities of heterogeneous surface chemistry.⁹ The gas phase decomposition reactions, surface adsorption, surface chemical reactions, and nucleation/growth in CVD and CVI processes can readily be investigated by Density Functional Theory (DFT).^19–21 DFT is a quantum mechanical method used to calculate the electronic structure of many-body systems. DFT calculations enable exploration of complex physico-chemical properties that are difficult or impossible to study by experiments.^22–24 In the current context, accurately accounting for electron correlation remains a significant challenge when describing complex molecules with strong bond potentials or analyzing variations in electron density during bond stretching or breaking. Due to these complexities, DFT-based studies of CVI and CVD are limited, often focusing solely on isolated modeling of gas-phase pyrolysis. Although valuable, only a few studies attempt to explore surface reactions, typically using oversimplified Transition State Theory (TST) models that do not reflect realistic conditions.

First-principles calculations can describe surface chemistry with high accuracy, but they require tremendous computing power. These methods solve complex quantum equations, which become harder and slower as more atoms and reaction steps are involved. Because of this, it's not practical to use them for simulating large systems or exploring many different conditions, like changes in temperature or pressure and mass transport. To overcome this, researchers use reduced-order or surrogate models as simplified versions that still capture the key chemistry but run much faster. These models significantly lower the cost of simulating complex reactions, making them useful for large-scale simulations and real-time optimization. With careful development, they can closely match the accuracy of first-principles methods while being much more efficient. For meaningful progress in advanced manufacturing technologies, it is essential to develop accurate models of surface reactions involving a representative set of key intermediate species, capturing sufficient chemical complexity and structural sensitivity.¹⁶

In this study, we employ DFT to explore the heterogeneous reactions of hydrocarbon intermediates on the SiC surface to uncover key reaction pathways. By developing surrogate kinetic models based on DFT-calculated activation energies and temperature-dependent reaction rates, we seek to bridge atomistic insights with practical deposition models. The following sections present the theoretical background and fundamental principles of DFT. The methodology section details the computational approaches used. We resolved dehydrogenation, hydrogenation, and resulting carbon deposition. Findings are then used to derive reaction rates, which are fitted to an Arrhenius-type surrogate kinetic models for into larger-scale simulations. The paper aims to provide a deeper understanding of the fundamental processes governing SiC deposition.

2. Decomposition of MTS and surface reactions

A critical step in the CVD and/or CVI of SiC is the breakdown of MTS through a sequence of reactions that reduce it to smaller molecules, ultimately leading to the deposition of SiC on a surface. These processes involve a complex set of coupled homogeneous gas-phase pyrolysis reactions and heterogeneous surface reactions, producing intermediate species that react to form SiC (Fig. 1).²⁵ The decomposition analysis of MTS begins with a thermodynamic evaluation of the gas-phase reactions and the resulting chemical species. The standard stable species are listed in the JANAF Thermochemical Tables and the NBS Tables of Chemical Thermodynamic Properties. Gas chromatography and various forms of spectroscopy experiments are key techniques used to detect many of these species, such as SiCl₄, CH₄, SiCl₃, CH₃, CH₃Cl₂, C₂H₄, SiCl₂, and C₂H₂.

Fig. 1 Illustration of material deposition through homogenous and heterogeneous reactions.

The general breakdown pathway of methyltrichlorosilane (MTS) typically involves intermediates such as methyl (CH₃) and silicon–chlorine (SiCl₃) species. In the most widely accepted descriptions, methyl and silicon–chlorine groups participate in intermediate reactions that lead to Si and C deposition pathways.²⁶ The intermediates formed from homogeneous gas-phase reactions generate the primary gaseous reactant species.²⁷ Studies have proposed a range of intermediate reaction mechanisms of differing complexity for the heterogeneous reactions occurring between Si- and C-bearing species on the surface.^28,29 The key challenge remains in identifying the most dominant and influential intermediate reactions to achieve accurate predictions while maintaining manageable levels of model complexity and computational cost. While a comprehensive review of these mechanisms is provided in literature, most efforts to date have concentrated primarily on detailed gas-phase chemistry. In contrast, studies addressing surface chemistry are relatively scarce, often considering only one or a few reaction steps involving a limited number of intermediate species.

A detailed set of primary deposition precursors was identified by Papasouliotis and Sotirchos et al.¹¹ through comprehensive chemical kinetics modeling. In their study, they proposed 14 heterogeneous reaction mechanisms for SiC deposition involving 13 gaseous precursors. The rate constants for these reactions were expressed using the Arrhenius equation or reactive sticking coefficients for adsorption processes. However, the model parameters were primarily calibrated using data from hot-wall chemical vapor deposition (CVD) reactor experiments. This highlights a broader issue in the current literature: deposition models trained on specific experimental conditions often lack generalizability and fail to capture the complexity of the system. These models, typically calibrated through limited experimental datasets, tend to be case-specific, resulting in poor predictive performance when applied to different operating conditions or reactor configurations. As a consequence, many experimental investigations report low reproducibility in deposition behavior and kinetics.

These issues can be addressed using a bottom-up modeling approach as the intricacies of heterogeneous surface chemistry are captured from the first-principles calculations. While empirical and semi-empirical models offer a starting point, they fall short in representing the true complexity of the surface reactions involved in CVD and CVI processes. In contrast, first-principles methods such as Density Functional Theory (DFT) allow for the simulation of individual reaction steps at the atomic level, including gas-phase decomposition, surface adsorption, bond breaking and formation, and subsequent nucleation and growth.

This approach enables the identification of reaction pathways and transition states, along with the calculation of key parameters (e.g. activation energies) that govern deposition kinetics. Moreover, DFT provides the capability to evaluate a wide range of intermediate species and reaction scenarios under varying local surface environments. By integrating these outputs into meso- or continuum-scale models, researchers can develop more accurate, physically-grounded simulations of SiC growth that overcome the limitations of reactor-specific empirical fits. Ultimately, such comprehensive modeling frameworks are essential to achieving predictive capabilities and enabling the rational design of deposition processes with greater reproducibility and material control.

3. Current methodology

The Vienna Ab initio Simulation Package (VASP) is employed to simulate heterogeneous surface chemistry based on quantum mechanical principles using DFT. It employs a plane-wave basis set, which offers computational efficiency for modeling periodic systems, particularly surfaces. To describe core-valence interactions, VASP implements the Projector Augmented-Wave (PAW) method, which provides higher accuracy and efficiency compared to norm-conserving or ultrasoft pseudopotentials.³⁰ The Perdew–Burke–Ernzerhof (PBE) functional within the Generalized Gradient Approximation (GGA) framework was used to approximate the exchange–correlation energy.³¹ Simulations were performed using the MedeA software interface.^32–35

For the model system, various hydrocarbon molecules (e.g., methane) were adsorbed on the surface of 6H–SiC (0001) using a (3 × 3) supercell, resulting in a 1% lattice mismatch. The SiC surface was constructed from a single-layer crystal of the P6₃mc SiC unit cell with a density of 3.2129 g cm⁻³. The stacking sequence used was 6H–SiC with the (0001) orientation, corresponding to the ABCACB stacking order. The crystal structure belongs to the hexagonal space group P6₃mc. A 10.0 Å vacuum gap was included to eliminate spurious interactions due to periodic boundary conditions. Atom-only structural optimization was subsequently carried out. First, an initial (reactant) and final (product) structure are created using the SiC crystal lattice. The reactant structure created by adding the hydrocarbon molecules of studied heterogeneous reaction. This structure was first geometry-optimized using the GGA–PBE exchange–correlation functional at an initial precision of 273.11 eV plane-wave cutoff energy. All subsequent calculations including final geometry optimizations and transition state searches were performed with a higher plane-wave cutoff of 500 eV to ensure convergence of adsorption and activation energies. Using real space projection operators and the Normal (blocked Davidson) method, the electronic iterations’ convergence is 10⁻⁶ eV. All calculations were spin-polarized. Although the clean 6H–SiC surface is non-magnetic, spin polarization was included to capture any magnetic effects potentially induced by hydrocarbon adsorption. During structural optimizations, the bottommost atomic layer of the slab model was fixed, while the remaining atoms were allowed to relax until the force convergence criterion of 0.02 eV Å⁻¹ was met. A Gaussian smearing with a width of 0.05 eV was applied in all geometry optimization and transition state search calculations. A representative snapshot of simulation domain is given Fig. 2. After geometry optimization and relaxation, the 6H–SiC structure was fully stabilized with lattice parameters of 9.3 Å × 9.3 Å in the periodic directions and 25.17 Å along the surface-normal direction. For all geometry optimizations and transition state searches, PAW–PBE pseudopotentials were employed. The following valence electron configurations were treated explicitly: Si (3s²3p², ZVAL = 4), C (2s²2p², ZVAL = 4), and H (1s¹, ZVAL = 1).

Fig. 2 A representative three-dimensional model system used in periodic DFT simulations.

Transition states were determined using the MedeA Transition State Search Module with the Nudged Elastic Band (NEB)^36,37 and Climbing Image NEB (CI-NEB)³⁸ methods. In NEB, the minimum energy path (MEP) between reactants and products is traced using a series of intermediate images connected by virtual springs, which are optimized to locate the approximate transition state. Calculations employed a plane-wave cutoff of 400 eV, a k-point spacing of 0.5 Å⁻¹, and GGA–PBE exchange–correlation functional. The NEB procedure was refined in three sequential steps (NEB0–NEB2), each using five intermediate images to iteratively improve the resolution near the highest energy point, providing a progressively more accurate estimate of the transition state geometry and energy.³⁹ In all calculations, the NEB method was used with a force convergence criterion requiring the maximum force on any image to be less than 1 eV Å⁻¹. Electronic self-consistent field (SCF) convergence was set to 10⁻⁵ eV.

The final transition state was obtained using the CI-NEB method, where the highest-energy image from NEB2 is refined by inverting the force component along the path to drive it toward the true saddle point. This approach ensures accurate determination of the activation energy and MEP. To verify the identified transition state, vibrational frequency analysis was performed using the finite displacement method, confirming the presence of a single imaginary frequency characteristic of a first-order saddle point.

Zone-center phonon (ZCP) calculations were carried out to determine the vibrational frequencies of both reactants and transition states, which are crucial for evaluating zero-point energies, thermal corrections, and reaction rates through transition state theory. The ZCP method involves single-point energy calculations and evaluates vibrational properties at the Γ-point (k = 0), representing the center of the Brillouin zone. To ensure accurate sampling, an odd-numbered k-point mesh centered on the Γ-point was employed. Using the finite displacement method, atoms were displaced by ±0.015 Å to compute forces and energy changes, and a Gaussian smearing width of 0.05 eV was applied to smooth electronic states.

Interatomic interactions were modeled using the Johnson potential (e.g., Johnson–R–A–Wadley–H–N–G–W potential), which provides reliable force predictions for phonon calculations. The resulting vibrational frequencies were used to calculate vibrational entropy, enthalpy, Gibbs free energy, and zero-point energy corrections, enabling accurate thermodynamic and kinetic analyses at various temperatures. Furthermore, gradients and stress tensors derived from the ZCP calculations provided insight into the mechanical stability and equilibrium of the structures. These results were integrated into statistical thermochemistry calculations to obtain zero-point energy and characteristic vibrational temperatures from the computed vibrational frequencies using eqn (1) and (2).

ZPE = hμ/2 (1)

Θ_v,i = hμ/k_B (2)

where ZPE is the zero-point energy, μ is the vibrational frequency of the mode (in Hz) from VASP zone center phonon calculations, h is the Planck's constant (6.626 × 10⁻³⁴ J s), Θ_v,i is the characteristic vibrational temperature of the i-th vibrational mode (in Kelvin), and k_B is the Boltzmann constant (1.381 × 10⁻²³ J/K). Then vibrational entropy is calculated from the summation over all vibrational modes usingwhere S_v is the total vibrational entropy of the molecule (J mol⁻¹ K) representing the entropy contribution arising from all quantized vibrational modes, R is the universal gas constant (8.314 J mol⁻¹ K) and T is the absolute temperature (in Kelvin) setting the thermal energy scale against which the vibrational excitation energy is compared. The first term on the left hand side accounts for the thermal occupation of the vibrational energy levels by quantifying the extent to which vibrational levels are populated at temperature T. The second term comes from the vibrational partition function and reflects the quantum nature of vibrations as it contributes negatively to entropy, capturing the restriction of accessible states at low temperatures. Together, these terms describe how each vibrational mode contributes to the overall entropy of a molecule, accounting for both quantum mechanical energy level spacing and thermal excitation.

Next, the total vibrational internal energy (E_v) is calculated using eqn (4) and it is used to calculate the total vibrational enthalpy for a given structure from eqn (5). Simply, the quantum mechanical ground-state vibrational energy estimated from ZPE is extended by including the thermal vibrational energy.

H = E_v + ZPE (5)

Next, the total vibrational enthalpy and entropy are used to calculate the vibrational contribution to Gibbs free energy (G_v) of the structure using eqn (6). Here, we combine both energetic and entropic contributions from molecular vibrations to determine how favorable a configuration or reaction is at a given temperature.

G_v = H − TS_v (6)

Using the vibrational contribution to Gibbs free energy of activation (G_v-a = G_{v-TransitionState} − G_v-Reactant) and energy barrier of the initial and the transition state structures, we calculate the Gibbs free energy of activation with eqn (7).

ΔG_a = ΔE_a + G_v-a (7)

Next, ΔG_a is used in the Eyring equation (eqn (8)) to calculate the reaction rates based on activation energy with the statistical thermochemistry corrections for temperature dependence.⁴⁰ As an extension of TST, the Eyring equation describes the thermodynamics of the transition state and provides a microscopic view of the reaction mechanism. In contrast, the Arrhenius equation is primarily empirical, derived from experimental observations. The reaction rate constant K can be described by the Eyring equation:

where R is the universal gas constant, T is the absolute temperature, k_B is the Boltzmann constant, and h is Planck's constant. In this way, the Eyring model provides reaction rates using microscopic quantities, offering a molecular-level interpretation of reaction dynamics.

While surface chemistry determined through detailed first-principles calculations offers high accuracy and predictive power, these methods are computationally intensive. Each calculation involves solving complex quantum mechanical equations for electron interactions, which becomes increasingly unmanageable as the number of atoms and possible reaction pathways grow. As a result, applying such high-fidelity calculations to simulate large, realistic systems or to explore a wide parameter space (e.g., varying temperatures, pressures, or flow conditions) becomes practically infeasible. This limitation underscores the need for reduced-order or surrogate models that can retain essential chemical accuracy while enabling efficient simulation of complete, system-scale behavior.

Surrogate modeling of surface chemistry offers a significant advantage by drastically reducing the computational cost of simulating complex reaction mechanisms, making them ideal for large-scale process simulations and real-time optimization. Models based on first-principles calculations can be developed with minimal loss of accuracy. For such a case, the Arrhenius model provides a macroscopic description of temperature dependence, typically expressed in terms of an activation energy, and is often favored as a practical approximation. The Arrhenius equation:

where A is the pre-exponential factor, and ΔE_a is the activation energy. As a last step, an Arrhenius form of surrogate model is developed from DFT calculations. Basically, Arrhenius parameters were calculated from a curve fit from DFT calculated reaction rates. For such a case, we plotted the logarithm of the rate constant versus the inverse temperature to employ the following form in eqn (10).where (E_a)_A is the activation energy for Arrhenius model. The activation energy is found from the slope and A parameter is calculated from the intercept.

4. Results and discussion

The literature indicates that the deposition of Si atoms is a favorable reaction pathway, whereas carbon deposition tends to limit surface growth.⁴¹ Therefore, the current study focuses on the heterogeneous surface reactions between hydrocarbon species and the SiC surface to better characterize the rate-limiting reaction pathways. Understanding these interactions is crucial for improving both the deposition rate and the quality of the resulting film. Papasouliotis et al.^6,42 identified CH₂ as key intermediates in hydrocarbon interactions, serving as primary carbon sources and proposed a comprehensive set of reaction pathways for hydrocarbons. However, existing models rely on parameters primarily tuned through experimental data, underscoring the need for a more accurate theoretical treatment. For such a case, we focus on heterogeneous hydrocarbon–SiC surface interactions set provided by Papasouliotis to examine carbon deposition.

Table 1 presents the adsorption reactions, gas-phase species interactions with interfacial species, and deposition reactions involved in the heterogeneous hydrocarbon–SiC surface chemistry. For simplicity, these are collectively referred to as surface reactions (SRs). There are total of five species (CH, CH₂, CH₃, CH₄, H, and H₂). We followed the same notation as the surface species are denoted by subscript s. Assuming that the silicon and carbon sites on the surface are energetically equivalent, [S] given in these equations represents for either of these two types of surface sites.

Table 1 Heterogeneous hydrocarbon–SiC surface interactions for the deposition of SiC

Type of the reaction	#	Reaction
Adsorption reaction	SR1	CH3 + [S] → [CH3]S
Reactions of gas phase species with interface species	SR2	[CH3]S + H → [CH2]S + H2
	SR3	[CH]S + H2 → [CH2]S + H
	SR4	[CH]S + H → [CH2]S
Deposition reactions	SR5	[CH2]S → C[S] + [S] + H2
	SR6	2[CH2]S + [S] → C[S] + [S] + CH4

---

In the following sections, each reaction is introduced and discussed in detail, supported by comprehensive DFT simulations. We present energy profiles along the reaction pathways, including relative energy variations obtained through transition state searches. Simulation snapshots are included to visualize nearly every bond-breaking and bond-forming event along the reaction paths, clearly illustrating the initial, transition, and final states of each reaction. Key thermodynamic properties such as activation energies, enthalpy, and entropy are reported with zero-point phonon corrections. Reaction rates are calculated using the Eyring equation, and generalized surrogate models are developed by fitting the data to Arrhenius-type expressions.

These reactions collectively demonstrate the intricate interactions between gas-phase molecules and surface-bound species, emphasizing the role of the surface in facilitating transformations. The surface provides a site for the adsorption, activation, and transformation of reactants, thereby influencing the reactivity and product distribution in these chemical processes. We start with adsorption reaction SR1 (CH₃ + [S] → [CH₃]_S) represents adsorption of the gas-phase methyl (CH₃) on an available surface site [S] yielding a surface-bound [CH₃]_s species. This is a simple adsorption mechanism without a true transition state, and applying a transition state search can be problematic. Since there is no real energy barrier to overcome, the algorithm may artificially force the system to locate a saddle point, which can introduce nonphysical or artificial energy differences if the calculations are not accurate (e.g., poor convergence, insufficient k-points, or inappropriate initial/final structures). Therefore, using NEB combined with Cl-NEB may be unnecessarily difficult and could lead to artifacts, but it can serve as a test case to evaluate whether the current transition state search calculations are robust.

The initial and final structures are carefully geometry optimized and used for the transition state search. Fig. 3(a) presents close-up snapshots of the initial and final structures. An initial standard NEB TSS calculation is performed between the initial and final structures. The NEB pathway was constructed by linearly interpolating five intermediate images between the initial and final geometries of the CH₃ adsorption process, followed by force relaxation and energy optimization. The resulting energy profile revealed a spontaneous and energetically favorable reaction pathway, characterized by a monotonic decrease in energy along the reaction coordinate (Fig. 3(b)). The reaction was exothermic yielding a net energy release of 395.97 kJ mol⁻¹. A similar adsorption energy reported by literature (∼360 kj mol⁻¹) validating the current calculation procedure.⁴³ Through the TSS steps, capturing a maximum force gradient dropping below 1 eV Å⁻¹, indicating a negligible energy barrier. NEB calculations could accurately swipe the reaction path way and followed relative energy change with no energy peak along the trajectory. During the TSS steps, the maximum force gradient dropped below 1 eV Å⁻¹, indicating a negligible energy barrier. The NEB calculations were able to accurately trace the reaction pathway, showing relative energy changes without any energy peak along the trajectory. Subsequently, the climbing image NEB method and RMM-DIIS optimization were applied to refine the calculations. These converged rapidly, within a single iteration, due to the very low energy difference (−0.34 kJ mol⁻¹), supporting the conclusion that the forward SR1 reaction is effectively barrierless.

Fig. 3 (a) Close-up snapshots illustrating the initial and final structures of SR1 highlighting the major events along the pathway. The initial and final images are labeled accordingly in snapshots. (b) Relative energy profile along the adsorption path obtained from the transition state search.

DFT calculations of charge distribution and stress tensors predict minimal electronic redistribution and no significant structural perturbation of the substrate. Band structure and Fermi energy analysis indicate a narrow width of 0.043 eV, which corresponds to unpaired electrons in the CH₃ radical. Smooth electronic occupancy was ensured using Gaussian smearing with a width of 0.05 eV, which is appropriate for simulations involving radical–surface interactions. Overall, the DFT-VASP results confirm that the SR1 step is spontaneous, exothermic, and barrierless.

Next, we focused on macroscopic modeling of SR1, where two main approaches exist for describing the adsorption mechanism: the sticking coefficient or the Arrhenius rate form. For the sticking coefficient approach, a barrierless adsorption process is typically represented by a sticking coefficient of unity. For the Arrhenius approach, we computed the temperature dependence of the adsorption process by calculating reaction rates and then fitting them with the Arrhenius equation. This yielded a small positive slope, corresponding to an apparent activation energy of ∼6 kJ mol⁻¹ with a pre-exponential factor of ∼1 × 10⁻¹³ s⁻¹. Importantly, this apparent E_a does not represent a true potential energy barrier; rather, it originates from the weak temperature dependence of the prefactor (e.g., partition function scaling, adsorption flux ∼T^−1/2). Such small effective activation energies are commonly observed in barrierless adsorption and diffusion-limited processes and are consistent with transition state theory predictions. The results are summarized in Table 2.

Table 2 Activation barriers calculated using DFT, along with Arrhenius-type surrogate kinetic model parameters derived from DFT reaction rates estimated by Eyring equation

	DFT	Arrhenius fit
	E a (kJ mol^−1)	A (1 s^−1)	E a (kJ mol^−1)
SR1	0	0.11 × 10^14	6
SR2	366.39	0.96 × 10^14	339.91
SR3	290.33	0.048 × 10^14	273.79
SR4	293.38	0.12 × 10^14	275.33
SR5	357.25	120 × 10^14	376.00
SR6	398.77	4.65 × 10^14	384.08

---

Next, we studied the dehydrogenation reactions SR2 ([CH₃]_S + H → [CH₂]_S + H₂) where a surface-bound methyl species ([CH₃]_s) lose a hydrogen when it reacts with another hydrogen (H), resulting in the formation of a surface-bound methylene species ([CH₂]_s) and the release of hydrogen gas (H₂). The optimized geometries of the initial and final states served as the basis for locating the transition state. In Fig. 4(a), detailed snapshots of these states illustrate the dominant bond dissociation and formation processes. To connect the two structures, a standard NEB transition state search was applied, and the resulting relative energy variation is depicted in Fig. 4(b). The NEB0 method resolves the reaction path using relatively crude convergence, while the transition state is refined through more accurate steps in NEB1 and NEB2.

Fig. 4 (a) Close-up snapshots illustrating the initial and final structures of SR2 highlighting the major bond-breaking and bond-forming events along the reaction pathway. The initial, transition state, and final images are labeled accordingly in snapshots. (b) Relative energy profile along the reaction path obtained from the transition state search. (c) Temperature-dependent reaction rates derived from the DFT-extended Eyring equation and fitted to macroscopic Arrhenius-type surrogate kinetic model.

During the NEB0 calculation, a series of five intermediate images is used to map the reaction pathway between the initial and final structures. For SR2, NEB0 identifies a preliminary transition state, located at the second image point where the energy is highest. Building on this, the search proceeds to the first refinement step (NEB1), where five new intermediate images are inserted between the first and third image points from NEB0. NEB1 offers a higher resolution of the reaction path and again finds a peak energy at the second image. In the third step, NEB2 further refines this region by placing five additional images between the two NEB1 images that bracket the energy maximum. NEB2 yields a more precise estimate of the transition state geometry and its associated energy. Finally, the Climbing Image NEB (Cl-NEB) method uses the NEB2-estimated coordinates to refine the transition state further. In this step, the climbing image is no longer affected by spring forces and instead moves directly along the true energy gradient to converge on the saddle point. This improves the accuracy of both the transition state structure and the calculated activation energy. Based on the maximum energy identified by Cl-NEB, the activation energy for the reaction is determined to be E_a = 366.39 kJ mol⁻¹. Literature suggests activation energies ranging between 300 and 400 kJ mol⁻¹ for similar set of dehydrogenation.⁴⁴ Moreover, Fermi level band crossing shows that the TS maintained its electronic character. Charge analysis and stress tensor results indicate that the TS is further stabilized through local distortions and charge redistribution at the reactive centers. The stability of the reaction mechanism pathway and the corresponding energy profile for SR2 is confirmed by the consistent convergence of the TS energy during climbing-image iterations and the smooth, continuous energy decrease observed throughout the NEB steps.

Next, we focus on reaction rate calculations. Based on the methodology described earlier, we performed ZCP calculations to determine the vibrational frequencies of the reactants and transition state structures (see Methodology, eqn (1)–(4)). From the computed vibrational frequencies, we evaluated the ZPE of each structure using eqn (1) and the characteristic vibrational temperature of each mode using eqn (2). These values provided the foundation for calculating vibrational contributions to the system's entropy and enthalpy. Specifically, the vibrational entropy was calculated using eqn (3), which accounts for the quantized nature of vibrational energy levels and their thermal population at a given temperature. The total vibrational internal energy was evaluated using eqn (4), and the vibrational enthalpy was obtained by summing with the ZPE as shown in eqn (5).

These vibrational contributions were essential for determining the Gibbs free energy of each structure. The vibrational component of the Gibbs energy, computed using eqn (6), captured both the energetic and entropic effects of molecular vibrations, which are crucial for assessing reaction thermodynamics. Using the vibrational Gibbs energies of the transition state and reactant structures, we calculated the Gibbs free energy of activation through eqn (7). This activation barrier, corrected for quantum and thermal vibrational effects, was then used in the Eyring equation (eqn (8)) to estimate the temperature-dependent reaction rates. We observed that the reaction rate increased significantly with temperature, consistent with expected thermally activated behavior. To provide a general kinetic model, we fitted the computed reaction rates to an Arrhenius-type expression as described by eqn (10). The resulting Arrhenius plot (Fig. 4(c)) allowed us to extract the activation energy from the slope and the pre-exponential factor from the intercept. These fitted parameters tabulated in Table 2 provide a surrogate model that simplifies the integration of atomistic reaction kinetics into larger-scale simulations or engineering models.

We continued with SR3 ([CH]_S + H₂ → [CH₂]_S + H) where a surface-bound methylene species ([CH]_s) reacts with a hydrogen molecule (H₂) to form a new surface-bound [CH₂]_s species, with the release of a hydrogen atom (H). This reaction illustrates hydrogen addition to a surface species, as shown in the snapshots in Fig. 5(a). The reactant and product structures were geometry-optimized using consistent interatomic potentials and used to identify the transition state for the SR2 reaction. The relative energy change shown in Fig. 5(b). A preliminary transition state was obtained using NEB0, with the saddle point located between the third and fifth intermediate images. The reaction path was further refined using NEB2 by introducing five additional intermediate images, narrowing the saddle point to between the second and fourth images. NEB2 yielded a more accurate transition state geometry and relative energy, which was then used as input for the Cl-NEB method to precisely determine the transition state coordinate and activation energy. The resulting energy barrier for the reaction was calculated as E_a = 290.33 kJ mol⁻¹. The Fermi energy falls within the band gap, ranging from −0.009 to +0.031 eV, while the transition state preserves its character, as evidenced by band crossings near the Fermi level. Stress tensor calculations indicate local lattice distortions along the pathway, and charge analysis confirms partial charges consistent with the presence of delocalized electrons.

Fig. 5 (a) Close-up snapshots illustrating the initial and final structures of SR3 highlighting the major bond-breaking and bond-forming events along the reaction pathway. The initial, transition state, and final images are labeled accordingly in snapshots. (b) Relative energy profile along the reaction path obtained from the transition state search. (c) Temperature-dependent reaction rates derived from the DFT-extended Eyring equation and fitted to macroscopic Arrhenius-type surrogate kinetic model.

Using the previously described methodology (eqn (1)–(10)), we computed vibrational frequencies and thermodynamic corrections to determine the Gibbs free energy of activation. This was then used in the Eyring equation to calculate temperature-dependent reaction rates, which were further fitted to an Arrhenius model to extract the activation energy and pre-exponential factor. Surrogate model parameters tabulated in Table 2.

Next, we studied SR4 ([CH]_S + H → [CH₂]_S) where a surface-bound methylene species ([CH]_s) reacts with a hydrogen atom (H) to form a new surface-bound [CH₂]_s species. This reaction illustrates hydrogen addition to a surface species, as shown in the snapshots in Fig. 6(a). The reactant and product structures were geometry-optimized using consistent interatomic potentials and used to identify the transition state. The reaction mechanism shown in Fig. 6(b) reveals that the primary transformation involves bond breaking or atom migration from the initial geometry, followed by a secondary rearrangement leading to product formation. As illustrated by the atomic configurations and the corresponding energy profile, the SR4 reaction proceeds via a two-step mechanism with two well-defined transition states along the reaction coordinate. This mechanistic interpretation is corroborated by the energy landscape presented in Fig. 6(b), which displays two distinct energy maxima corresponding to transition states TS1 and TS2. Initial reaction pathways were approximated using NEB0 and NEB1 calculations and subsequently refined via the NEB2 path. The first TS is characterized by a significant energy barrier of approximately 293 kJ mol⁻¹, while the second TS exhibits a relatively lower barrier of approximately 176 kJ mol⁻¹.

Fig. 6 (a) Close-up snapshots illustrating the initial and final structures of SR4 highlighting the major bond-breaking and bond-forming events along the reaction pathway. The initial, transition state, and final images are labeled accordingly in snapshots. (b) Relative energy profile along the reaction path obtained from the transition state search. (c) Temperature-dependent reaction rates derived from the DFT-extended Eyring equation and fitted to macroscopic Arrhenius-type surrogate kinetic model.

The transition state structures obtained through Cl-NEB calculations were further optimized with maximum residual forces slightly above the threshold (approximately 2.3 eV Å⁻¹). Energy plateaus were reached, indicating that the optimized configurations are reasonable approximations of the accurate TS geometries. The overall energy difference between the initial and final states is approximately 107.1 kJ mol⁻¹, confirming that the SR4 reaction is endothermic. Through the two-step mechanism of SR4, the high-energy TS1 governs the reaction kinetics.

We calculated vibrational frequencies and applied thermodynamic corrections to obtain the Gibbs free energy of activation following the methodology outlined in eqn (1)–(10). These values were used in the Eyring equation to compute temperature-dependent reaction rates, which were subsequently fitted to an Arrhenius model as illustrated in Fig. 6(c) to extract the activation energy and pre-exponential factor. The resulting surrogate model parameters are summarized in Table 2. The integration of NEB simulations, TS optimizations, and Arrhenius modeling provides a comprehensive understanding of the SR4 pathway under CVI conditions.

The next reaction is the first deposition reaction of surface-bound methylene species. SR5 ([CH₂]_S → C_[S] + [S] + H₂) represents the decomposition of a surface-bound methylene group ([CH₂]_s) into a surface carbon atom (C_s) releasing a molecular hydrogen (H₂). This is an important interaction depositing C onto surface as shown in the snapshots in Fig. 7(a). A one-step reaction with five intermediate images and a spring constant of 5 eV Å⁻² was found by the NEB calculations using linear interpolation and two refinement steps. Both refined TS by climbing image NEB and final TS confirm that the reaction surmounts a single near-maximum energy barrier of about 357.5 kJ mol⁻¹, as indicated by the energy profile in Fig. 7(b). The activation energy for the dehydrogenation of surface-bound methylene studied on highly active catalytic Fe₂O₃ surface presented as 142 kJ mol⁻¹ that activation on SiC is expected to be between 150 and 350 kJ mol⁻¹.^45,46 The cleavage and C–Si bond formation along the reaction coordinate are highlighted in the reaction sequence, which also clearly illustrates the progressive migration and desorption of surface-bound species. Additionally, the electronic structure analysis of the transition state optimized shows that the system remains its nature, with the Fermi level located within a very close energy window of −0.051 to 0.018 eV, indicating a non-zero density of states at the Fermi level throughout the reaction.

Fig. 7 (a) Close-up snapshots illustrating the initial and final structures of SR5 highlighting the major bond-breaking and bond-forming events along the reaction pathway. The initial, transition state, and final images are labeled accordingly in snapshots. (b) Relative energy profile along the reaction path obtained from the transition state search. (c) Temperature-dependent reaction rates derived from the DFT-extended Eyring equation and fitted to macroscopic Arrhenius-type surrogate kinetic model.

The final reaction is the second deposition reaction of surface-bound methylene species. SR3 (2[CH₂]_S + [S] → C_[S] + [S] + CH₄) reaction involves two surface-bound methylene species ([CH₂]_s) resulting in the formation of a carbon species (C_s), and methane (CH₄) as a product. This reaction involves the coupling of surface species, resulting in carbon deposition and methane formation. Our DFT calculations successfully captured this process, as illustrated by the snapshots in Fig. 8(a).

Fig. 8 (a) Close-up snapshots illustrating the initial and final structures of SR6 highlighting the major bond-breaking and bond-forming events along the reaction pathway. The initial, transition state, and final images are labeled accordingly in snapshots. (b) Relative energy profile along the reaction path obtained from the transition state search. (c) Temperature-dependent reaction rates derived from the DFT-extended Eyring equation and fitted to macroscopic Arrhenius-type surrogate kinetic model.

The reaction pathway, as visualized in Fig. 8(b), shows that the surface mediated coupling begins with the two adjacent [CH₂]_s groups undergo a rearrangement where one carbon donates a hydrogen to the other, facilitating the formation of a CH₄ molecule. The initial step involves hydrogen transfer between the [CH₂] units, forming an intermediate structure resembling CH₃–CH, where one carbon becomes saturated and the other unsaturated. This hydrogen transfer has a moderate energy barrier showing itself as the first transition state. As the CH₄ molecule forms, it detaches from the surface, leaving behind a more reactive carbon species that remains chemisorbed.

The transition state for this reaction is determined using a stepwise NEB approach, starting with NEB0 to identify a rough energy maximum along the reaction path. Successive refinements using NEB1 and NEB2 added additional intermediate images to increase resolution. Two transition states with a reaction intermediate in between determined. Climbing Image NEB (Cl-NEB) was then applied to both transition state points that second one representing CH₄ desorption gave the highest energy barrier as E_a = 398.779 kJ mol⁻¹. Electronic structure analysis based on the Fermi energy reveals that the system remains with the Fermi level positioned within a narrow band that intersects the conduction states. Both deposition barrier values are consistent with the experimental findings of Kojima et al.,⁴⁷ who reported activation energies of ∼515 kJ mol⁻¹ for SiC decomposition.

From our DFT calculations, this transformation proceeds through a relatively well-defined transition state, characterized by the simultaneous breaking of a C–H surface bond and the formation of a C–H bond in CH₄. The snapshots from the NEB calculation reveal that the transition state corresponds to the point where the CH₄ molecule is nearly formed but still weakly interacting with the surface. This reaction illustrates two key surface phenomena: hydrocarbon coupling and carbon deposition. Temperature dependent reaction rates determined by including phonon calculations. The model parameters for the Arrhenius-type surrogate kinetic model are determined from Fig. 8(c) and tabulated in Table 2.

5. Conclusion

This work presents a comprehensive investigation of hydrocarbon–SiC surface reactions critical to carbon deposition during SiC film growth. Building on the framework proposed by Papasouliotis et al., we analyzed three key surface reactions through DFT calculations, identifying transition states and calculating activation energies with zero-point phonon corrections. The reaction energy landscapes revealed that SR6, representing methylene coupling and carbon deposition, has the highest energy barrier, making it the most kinetically limiting step. Surface dehydrogenation (SR2) and hydrogenation (SR3 and SR4) were found to play essential roles in forming the reactive [CH₂]_s intermediate that facilitates carbon incorporation.

Temperature-dependent reaction rates were computed using transition state theory and vibrational thermodynamics, then fitted to Arrhenius-type surrogate models to enable integration into mesoscale simulations. These models captured the expected thermally activated behavior and provide useful parameters for reactor-scale process modeling. Collectively, this study highlights the mechanistic complexity of carbon incorporation on SiC surfaces and emphasizes the importance of surface hydrogen dynamics and methylene availability in controlling film growth rates and quality. The approach demonstrated here serves as a template for future kinetic modeling of other heterogeneous systems relevant to materials synthesis.

We observe key differences from the mechanisms and rates reported by Papasouliotis et al., where some reactions were treated as barrierless or assigned high-rate constants. While our DFT-derived kinetics suggest mechanistic pathways and rate parameters that differ from the prevailing semi-empirical understanding in the literature, a direct and quantitative comparison requires system-level validation. To assess their impact, furnace-level validation is required. By coupling these first-principles kinetics with reactor-scale simulations, a more transferable and physically grounded framework can clarify the limitations of earlier semi-empirical approaches.

The findings underscore the utility of atomistic modeling in isolating rate-limiting surface events that are difficult to determine experimentally. By explicitly resolving transition state geometries and vibrational contributions, we provide an accurate, first-principles thermodynamic and kinetic description that overcomes limitations of typically employed, purely empirical models. Importantly, the identified energy barriers and reaction pathways help explain why carbon buildup may inhibit surface growth, especially in hydrogen-deficient environments where CH₂ recombination dominates. This study offers both mechanistic insight and practical tools for guiding process optimization in chemical vapor deposition of SiC.

Conflicts of interest

There are no conflicts to declare.

Notes

This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a nonexclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (https://energy.gov/downloads/doe-public-access-plan).

Data availability

The code for VASP can be found at https://www.vasp.at/ with https://doi.org/10.1007/978-1-4615-5943-6_10. The version of the code employed for this study is version VASP.6.5.0.

Acknowledgements

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Award Number DE-SC0024510. The authors acknowledge the support of the University of Tennessee at Chattanooga; research reported in this publication was partially supported by the FY2025-28 Center of Excellence for Applied Computational Science competition.

References

1. I. Donald and P. McMillan, Ceramic-matrix composites, J. Mater. Sci., 1976, 11, 949–972.
2. J. Binner, et al., Selection, processing, properties and applications of ultra-high temperature ceramic matrix composites, UHTCMCs–a review, Int. Mater. Rev., 2020, 65(7), 389–444.
3. P. Delhaes, Chemical vapor deposition and infiltration processes of carbon materials, Carbon, 2002, 40(5), 641–657.
4. H. Li, et al., Numerical simulation of chemical vapor infiltration of propylene into C/C composites with reduced multi-step kinetic models, Carbon, 2005, 43(14), 2937–2950.
5. R. Bickerdike, et al., The deposition of pyrolytic carbon in the pores of bonded and unbonded carbon powders, in Proceedings of the Fifth Conference on Carbon, Pergamon Press, London, 1962.
6. G. D. Papasouliotis and S. V. Sotirchos, Heterogeneous kinetics of the chemical vapor deposition of silicon carbide from methyltrichlorosilane, MRS Online Proc. Libr., 1993, 334, 111.
7. B. W. Sheldon and T. M. Besmann, Reaction and diffusion kinetics during the initial stages of isothermal chemical vapor infiltration, J. Am. Ceram. Soc., 1991, 74(12), 3046–3053.
8. J. K. Nørskov, et al., Density functional theory in surface chemistry and catalysis, Proc. Natl. Acad. Sci. U. S. A., 2011, 108(3), 937–943.
9. A. Bruix, et al., First-principles-based multiscale modelling of heterogeneous catalysis, Nat. Catal., 2019, 2(8), 659–670.
10. K. Guan, et al., Numerical modeling of SiC by low-pressure chemical vapor deposition from methyltrichlorosilane, Chin. J. Chem. Eng., 2020, 28(6), 1733–1743.
11. G. D. Papasouliotis and S. V. Sotirchos, Experimental study of atmospheric pressure chemical vapor deposition of silicon carbide from methyltrichlorosilane, J. Mater. Res., 1999, 14, 3397–3409.
12. J. Xue, et al., Thermodynamic Analysis on the Codeposition of SiC–Si 3 N 4 Composite Ceramics by Chemical Vapor Deposition using SiCl 4–NH 3–CH 4–H 2–Ar Mixture Gases, J. Am. Ceram. Soc., 2013, 96(3), 979–986.
13. S. H. Mousavipour, V. Saheb and S. Ramezani, Kinetics and mechanism of pyrolysis of methyltrichlorosilane, J. Phys. Chem. A, 2004, 108(11), 1946–1952.
14. Y. Ge, et al., Theoretical study of the pyrolysis of methyltrichlorosilane in the gas phase. 2. Reaction paths and transition states, J. Phys. Chem. A, 2007, 111(8), 1475–1486.
15. T. Ogawa, et al., CFD simulation of CVD reactors in the CH₃SiCl3 (MTS)/H2 system using a two-step MTS decomposition and one-step SiC growth models, Heliyon, 2023, 9(4), e15061.
16. X. Wang, et al., Initial decomposition of methyltrichlorosilane in the chemical vapor deposition of silicon-carbide, Comput. Theor. Chem., 2011, 967(2–3), 265–272.
17. J. Wang, et al., The rate-limiting step in the thermal oxidation of silicon carbide, Scr. Mater., 2010, 62(9), 654–657.
18. Y. M. Manawi, et al., A review of carbon nanomaterials’ synthesis via the chemical vapor deposition (CVD) method, Materials, 2018, 11(5), 822.
19. S. Mahamulkar, et al., Formation and oxidation/gasification of carbonaceous deposits: a review, Ind. Eng. Chem. Res., 2016, 55(37), 9760–9818.
20. S. Chabi and K. Kadel, Two-dimensional silicon carbide: emerging direct band gap semiconductor, Nanomaterials, 2020, 10(11), 2226.
21. P. Sukkaew, O. Danielsson, O. Kordina, E. Janzén and L. Ojamae, Ab initio study of growth mechanism of 4H–SiC: adsorption and surface reaction of C2H2, C2H4, CH4, and CH₃, J. Phys. Chem. C, 2017, 121(2), 1249–1256.
22. P. Geerlings and F. De Proft, Conceptual DFT: the chemical relevance of higher response functions, Phys. Chem. Chem. Phys., 2008, 10(21), 3028–3042.
23. T. Ziegler, Approximate density functional theory as a practical tool in molecular energetics and dynamics, Chem. Rev., 1991, 91(5), 651–667.
24. H. B. Schlegel, Exploring potential energy surfaces for chemical reactions: an overview of some practical methods, J. Comput. Chem., 2003, 24(12), 1514–1527.
25. H.-C. Liu, et al., Stepwise Reaction for Chemical Vapor Deposition of Stoichiometric SiC Films Using Methyltrichlorosilane and Hydrogen as Reactants, Ind. Eng. Chem. Res., 2024, 63(36), 15800–15808.
26. B. Heidenreich, C/SiC and C/C–SiC Composites, in Ceramic matrix composites: materials, modeling and technology, 2014, pp. 147–216.
27. C. Lu, et al., Kinetics of chemical vapor deposition of SiC from methyltrichlorosilane and hydrogen, Appl. Surf. Sci., 2009, 255(17), 7495–7499.
28. G. S. Fischman and W. T. Petuskey, Thermodynamic Analysis and Kinetic Implications of Chemical Vapor Deposition of Sic from Si–C–C1–H Gas Systems, J. Am. Ceram. Soc., 1985, 68(4), 185–190.
29. T. Sen and M. Barisik, Internal surface electric charge characterization of mesoporous silica, Sci. Rep., 2019, 9(1), 137.
30. J. Hafner, Ab–initio simulations of materials using VASP: Density–functional theory and beyond, J. Comput. Chem., 2008, 29(13), 2044–2078.
31. G. Kresse and J. Furthmüller, Efficiency of ab initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci., 1996, 6(1), 15–50.
32. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54(16), 11169.
33. M. Iskin and C. Sá de Melo, BCS-BEC crossover of collective excitations in two-band superfluids, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 72(2), 024512.
34. MedeA®: Comprehensive materials design and simulation software, Materials Design Inc, Santa Fe, New Mexico, USA, 2013.
35. H. G. Ozcelik, Y. Sozen, H. Sahin and M. T. Barisik, Parametrizing nonbonded interactions between silica and water from first principles, Appl. Surf. Sci., 2020, 504, 144359.
36. H. Jónsson, G. Mills and K. W. Jacobsen, Nudged elastic band method for finding minimum energy paths of transitions, in Classical and quantum dynamics in condensed phase simulations, World Scientific, 1998, pp. 385–404.
37. G. Henkelman and H. Jónsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, J. Chem. Phys., 2000, 113(22), 9978–9985.
38. G. Henkelman, B. P. Uberuaga and H. Jónsson, A climbing image nudged elastic band method for finding saddle points and minimum energy paths, J. Chem. Phys., 2000, 113(22), 9901–9904.
39. H. G. Ozcelik, E. Satiroglu and M. Barisik, Size dependent influence of contact line pinning on wetting of nano-textured/patterned silica surfaces, Nanoscale, 2020, 12(41), 21376-91.
40. T. Sen and M. Barisik, Internal surface electric charge characterization of mesoporous silica, Sci. Rep., 2019, 9(1), 137.
41. M. D. Allendorf and C. F. Melius, Theoretical study of thermochemistry of molecules in the silicon-carbon-chlorine-hydrogn system, J. Phys. Chem., 1993, 97(3), 720–728.
42. G. D. Papasouliotis and S. V. Sotirchos, On the homogeneous chemistry of the thermal decomposition of methyltrichlorosilane: Thermodynamic analysis and kinetic modeling, J. Electrochem. Soc., 1994, 141(6), 1599.
43. P. Sukkaew, O. Danielsson, O. Kordina, E. Janzen and L. Ojamae, Ab initio study of growth mechanism of 4H–SiC: adsorption and surface reaction of C2H2, C2H4, CH4, and CH₃, J. Phys. Chem. C, 2017, 121(2), 1249–1256.
44. P. Sukkaew, O. Danielsson and L. Ojamae, Growth mechanism of SiC CVD: surface etching by H2, H atoms, and HCl, J. Phys. Chem. A, 2018, 122(9), 2503–2512.
45. Z. Cheng, L. Qin, M. Guo, J. A. Fan, D. Xu and L. S. Fan, Methane adsorption and dissociation on iron oxide oxygen carriers: the role of oxygen vacancies, Phys. Chem. Chem. Phys., 2016, 18(24), 16423–16435.
46. N. Chen and F. Sun, Atomic–Level Investigation of CHx and C2Hx Adsorption on β–SiC (111) Surface for CVD Diamond Growth from DFT Calculations, J. Nanomater., 2011, 2011(1), 910372.
47. K. Kojima, S. Kuroda, H. Okumura and K. Arai, Effect of Additional Silane on In Situ H2 Etching Prior to 4H-SiC Homoepitaxial Growth, Mater. Sci. Forum, 2007, 556–557, 85–88.

---