Cite this Article
Sun Yu†, Liu Yi-Jun, Xu Fei. ReaxFF molecular dynamics study on oxidation behavior of 3C-SiC: Polar face effects. Chinese Physics B, 24(9): 096203  
Permissions
ReaxFF molecular dynamics study on oxidation behavior of 3C-SiC: Polar face effects
Sun Yu†a), Liu Yi-Juna),b), Xu Feia)
Institute for Computational Mechanics and Its Applications, Northwestern Polytechnical University, Xi’an 710072, China
Mechanical Engineering, University of Cincinnati, Cincinnati, Ohio 45221-0072, USA

†Corresponding author. E-mail: npuyusun@gmail.com

*Project supported by the 111 Project (Grant No. B07050) and the National Natural Science Foundation of China (Grant No. 11402206).

Abstract

The oxidation of nanoscale 3C-SiC involving four polar faces (, Si(100),, and Si(111)) is studied by means of a reactive force field molecular dynamics (ReaxFF MD) simulation. It is shown that the consistency of 3C-SiC structure is broken over 2000 K and the low-density carbon chains are formed within SiC slab. By analyzing the oxygen concentration and fitting to rate theory, activation barriers for, Si(100),, and Si(111) are found to be 30.1, 35.6, 29.9, and 33.4 kJ·mol−1. These results reflect lower oxidative stability of C-terminated face, especially along [111] direction. Compared with hexagonal polytypes of SiC, cubic phase may be more energy-favorable to be oxidized under high temperature, indicating polytype effect on SiC oxidation behavior.

PACS: 62.25.–g; 81.16.Pr; 61.46.–w
Keyword: molecular dynamics; ReaxFF field; 3C-SiC; oxidation
1. Introduction

To improve the performance of space vehicles for extremely harsh service environment, it is crucial to develop the materials and structures for their thermal protection systems. Due to the advanced features, composite materials based on carbon, silicon carbide (SiC), silica (SiO2), and carbon fibers become the latest choice for the thermal protection system.[1– 3] However, because of the high-temperature and high-pressure working condition, oxidation occurs which can lead to system degradation.[4– 7] Therefore, the understanding of related behaviors involving chemical reaction and microstructure evolution turns out to be a critical issue.

Plenty of experimental efforts have been devoted to explore the oxidation process of these silicon carbide materials, while limitation exists for reaching temperature higher than 2000 K as well as understanding the break and formation of molecular bonding during chemical reaction.[8– 17] Besides, recent studies on compressive properties of SiC composite materials also found temperature-dependent phenomena and attributed them to the oxidation of SiC. Therefore, it is difficult to provide a clear explanation due to the limitation of observation condition.[18] Modeling approaches have their advantages and can play important roles in the study. Previously, calculations based on density functional theory (DFT) were utilized to investigate the oxygen adsorption behavior for different SiC polytypes.[19– 21] Although stable configurations and reaction pathway have been achieved, computational cost was too high to contain sufficient numbers of atoms.[22– 24] Classical molecular dynamics simulation on oxygen diffusion in SiC was also performed without providing information on chemical reaction details.[25] Basically, the oxidation process includes chemical reactions, mass transport, and electrostatic interactions, all of which have its own complexity.[26, 27] To be able to model the profile at a lower calculation cost, a sophisticated bond-order potential, ReaxFF method, was established by the group of Duin.[28] The parameters have been extended for Si/C/H/O system by Newsome et al. They studied the oxidation process of 2H-SiC(0001) by O2 and H2O in a temperature range of 500– 5000 K.[29] Furthermore, by developing the Deal– Grove (D– G) theory, results from MD simulation have been fitted to compute the transport diffusivity at various temperatures and further to obtain the activation barriers.[30] They focused on 2H-SiC(0001) face without providing information on the cubic SiC phase (3C-SiC), which is technologically important related to crystal growth.[9, 13, 19, 31] Besides, the anisotropic oxidation nature of SiC has not been taken into account by MD simulations according to the best knowledge of the authors.

In the present work, we concentrate on the polar face effects, to first investigate the oxidation behavior of , Si(100), , and Si(111) of 3C-SiC by using the ReaxFF MD simulation, and then calculate the Arrhenius parameters for each face following a similar framework of Newsome et al.[29, 30]

2. ReaxFF MD modeling

The crystal structures of 3C-SiC polar faces including , Si(100), , and Si(111) are shown in Fig. 1. To perform the simulation, a bulk of SiC containing 640 atoms with (100) surface normal to the z axis was set up, the size of which was 17.4 × 17.4 × 20.6 Å . The simulation cell was created by extending the length along the z direction to 52.3 Å . Periodic boundary condition was used for such system and both the C-terminated and Si-terminated face were considered, respectively. Similarly, in order to investigate the effect of different crystal orientations, atomistic models with [111] direction along the z axis with C- and Si-terminated face were set up as well, with the simulation cell having the size of 16.0 × 15.4 × 54.0 Å . Within the rest of the simulation cell, 200 O2 molecules were randomly inserted. Correspondingly, the oxygen concentration was 3.48 × 104 and 4.58 × 104 mol· m− 3 for (100) and (111) models, respectively, which were high-pressure cases following Newsome’ s work.[29] Hydroxyl groups were added bonding with surface atoms to protect the SiC surface which is highly reactive. The cutoff radius was set as 10 Å . The length of SiC slab with [100] and [111] direction was respectively 21 Å and 25 Å , either of which was larger than twice of the cutoff radius. Initial configurations were obtained by minimizing the whole system and then running MD at 100 K for several steps to make O2 molecules well-distributed, as shown in Fig. 2. After that, the temperature was quickly increased to 800, 1000, 1400, 2000, 3000, 4000, 5000 K, controlled by a Berendsen thermostat to simulate the NVT (Canonical) ensemble. The temperature damping constant was set as 100 fs and the time step was adjusted with the value of temperature, 0.25 fs for 2000 K and below, while 0.1 fs for higher temperature. Nose– Hoover thermostat has also been applied for one case to verify the thermostat effect, showing no obvious difference. Furthermore, modeling of SiC slab with thickness of 30 Å was carried on. Lager system has been tested as well which verifies the convergence of the calculations. The initial stage of bonds breaking and forming which are of most interests was similar as in the smaller size model. Since the oxidation process is governed by breaking and forming of covalent bonds which only involve very few angstroms, the structure we used here is large enough to model plenty of oxidizing reactions and the surface oxidation of SiC.[29] All the simulations were conducted using the open-source LAMMPS code.[32, 33]

Fig. 1. Side view of 3C-SiC polar face. In this and all the following figures showing atomistic structures, black, and gray spheres indicate carbon and silicon atoms, respectively. (a) , (b) Si(100), (c) , (d) Si(111).
Fig. 2. Initial configurations for MD. In this and all the following figures showing atomistic structures, red, and white spheres indicate oxygen and hydrogen atoms, respectively. (a) , (b) Si(100), (c) , (d) Si(111).
3. Results and discussion
3.1. Atomistic structure analysis

The configurations with [100] and [111] directions for C- and Si-terminated face oxidized under 1000 K, 2000 K, 3000 K, and 5000 K are presented in Figs.3 and 4, respectively. All of them are observed by ReaxFF MD at 80 ps. It is noticeable that all the four polar faces turn to react very fast with oxygen when the temperature is higher. Above 2000 K, there is a mixture of carbon chain and SiOx in the middle part of simulation cell, and almost all the oxygen molecules have been consumed according to Figs. 3 and 4. Since the melting point of SiO2 is nearly 1800 K given by Castro-Marcano et al., [29, 34] the separated slab is forming a liquid, amorphous silica phase as temperature increases above 2000 K. While comparing with the structures of 2H-SiC oxidation given by Newsome et al., [29] the consistency seems to be broken at earlier stage for 3C-SiC. These results imply that the oxidation degree of 3C-SiC could be greater than that of 2H-SiC. However, by simply comparing the ReaxFF MD configurations of , Si(100), , and Si(111), no visible difference has been shown.

To further understand the structure evolution, the numbers of O2 molecules changing with time at various temperatures have also been tracked during the simulation. The O– O bond consumptions for (100) and (111) surfaces oxidized at 1000 K, 3000 K, and 5000 K are illustrated as Fig. 5. The consumption rates of oxygen increase greatly as temperature increasesfor all cases. C-face consumes O– O bond faster than that of Si-face, with a more evident difference at 1000 K, while less obvious under higher temperature. It may be attributed to the activity enhancement of subsurface as temperature increases, so that the effect of terminated element as C or Si on surface is weakened. By means of oxygen concentration data at 800, 1000, 1400, 2000, 3000, 4000, and 5000 K obtained from MD simulation, analysis on transport rate and energy barriers is conducted in subsection 3.2.

Fig. 3. Configurations of (100) face oxidation at 80 ps. oxidation at (a1) 1000 K, (a2) 2000 K, (a3) 3000 K, (a4) 5000 K. Si(100) oxidation at (b1) 1000 K, (b2) 2000 K, (b3) 3000 K, (b4) 5000 K.
Fig. 4. Configurations of (111) face oxidation at 80 ps. oxidation at (a1) 1000 K, (a2) 2000 K, (a3) 3000 K, (a4) 5000 K. Si(111) oxidation at (b1) 1000 K, (b2) 2000 K, (b3) 3000 K, (b4) 5000 K.

Moreover, to investigate the effect of initial oxygen concentration, 300 O2 molecules were inserted and relative simulations at 2000 K and 3000 K have been performed for the model with Si(100) face. The configurations after oxidizing for 80 ps are presented in Fig. 6. It has been shown that there are CO2 molecules formed in the middle part of the slab at 3000 K, where SiC was separating and producing graphitized layer. The analysis of gas concentration changing with time at 2000 K and 3000 K is illustrated in Fig. 7, including O2, CO, and CO2. Obviously, CO2 molecules have been generated when both the initial concentration of O2 and system temperature are high. It seems that excess O2 molecules have passed through the liquid, amorphous silica phase and totally reacted with carbon layers under 3000 K to form CO2. These results could be instructive to realistic situation under high temperature and high pressure to some extent.

Fig. 5. Consumption of O2 with time during the oxidation at 1000 K, 3000 K, and 5000 K. (a) (100) face. (b) (111) face.

Fig. 6. Configurations of SiC oxidation at 2000 K (a) and 3000 K (b) for 80 ps. Systems contain 300 O2 molecules initially.

Fig. 7. Comparison of systems with different initial O2 concentrations at 2000 K and 3000 K. (a) 2000 K, contained 200 O2 molecules initially. (b) 2000 K, contained 300 O2 molecules initially. (c) 3000 K, contained 200 O2 molecules initially. (d) 3000 K, contained 300 O2 molecules initially.

3.2. Estimation of rate constant and activation energy

The transport rate theory to model gas and solid reaction was developed by Newsome et al. It provides a way to quantify the rate parameters, which is difficult to conduct via experiments. The relation between oxygen concentration and reaction time is given by[30]

where and are the concentrations of bulk fluid containing O atoms changing with time and at initial stage, is the initial density of O atoms in the fluid phase of simulation cell, which is calculated by total number of O atoms divided by fluid phase volume, Lf is taken as the length of the simulation cell, Lslab, 0 is the initial thickness of SiO2 slab, which is used as the thickness of hydroxyl layer nearly 2.8 Å , and Dt is the transport diffusivity.

By fitting the MD data shown in Fig. 5 with the concentration-time relationship as Eq. (1), the values of α and β for each temperature are obtained. Then, using Eq. (2) or Eq. (3), two values of Dt could be calculated from α or β . By weighted averaging, these two Dt values to maximize the correlation coefficient of theory and simulation results, final determination of Dt is made. After all the Dt at various temperatures as 800, 1000, 1400, 2000, 3000, 4000, 5000 K for (100) and (111) face are obtained, activation energies and pre-exponential factors are calculated by Arrhenius equation as follows:

Take the logarithm to make a linear form as

where Ea and D0 are the activation energy and pre-exponential factor respectively, and R is the gas constant taken as 8.3145 J· K− 1· mol− 1.

The fitting of lnDt with 1000/T as Eq. (5) for (100) and (111) face are shown in Fig. 8. The correlation coefficients of all cases are larger than 0.95. By these fittings, Ea and D0 are calculated for each situation, as displayed in Table 1. The activation barriers of C-face oxidation are lower than that of Si-face for both (100) and (111) crystal planes, which indicates the oxidation facility of C-terminated polar face. Previous results have shown that the oxide growth rate of Si-face was about 1/10th that of C-face, and both rates were much higher than the evaluation given by D– G model at oxidation thickness less than 20 nm. However, almost all the studies only focused on the oxidation properties under temperature lower than 2000 K for hexagonal structures.[14, 35, 36] Here, the small difference between C-face and Si-face presented in Table 1 may result from the influence of the data point at higher temperature. As temperature increases, the polar face effects seem to be reduced. Moreover, at the initial stage involving extremely thin oxide layer, the oxidation rate is high and correspondingly leads to smaller activation barriers. Apart from the difference between C- and Si-face, comparison between (100) and (111) plane reveals slight energy-favorable of (111) plane, while it is not noticeable.

Fig. 8. Fitting of Arrhenius equation for oxidation of 3C-SiC (a) (100) face and (b) (111) face.

Table 1. Activation energy and pre-exponential factor for SiC oxidation with different polar faces.

All the Ea values are less than that of 2H-SiC (see Ref. [30]: 46.3 kJ· mol− 1), revealing higher oxidation rate of 3C-SiC. D0 is similar for each case, nearly 10− 2 cm2· s− 1. It has reported that the oxidation rates of SiC were dependent on crystal polytypes.[9] An experimental oxidation test on β -SiC and the comparison with α -SiC done by the same group revealed that the activation energy for β -SiC were much lower than α -SiC.[37, 38] These seem to be in qualitative agreement with our calculations which show that the value of activation energy of 3C-SiC oxidation is less than that of 2H-SiC. However, since the oxygen concentration is much higher in the present simulation than in experiments, all the values of activation energy obtained here are obviously lower. Additionally, research on the formation of porous SiC ceramics manifested that the cubic phase with oxide was unstable and further converted to hexagonal form.[39] These could also indicate the weaker oxidative stability of 3C-SiC in agreement with our simulation results.

To find the root of polytype effect on oxidation behavior, it is necessary to analyze the differences of their crystal structures. Among SiC polymorphs, 3C-SiC is the only cubic one with stacking sequence as · · · ABCABC· · · along its [111] direction. This layer stacking direction is defined as the c axis along [0001] direction for hexagonal structures, along which 2H-SiC has · · · ABAB· · · sequence. Following the notation rule by Jagodzinski, [40] the letter h is used to describe the stacking position with a 180° rotation, whereas the letter k refers to the stacking without rotation. As shown in Fig. 9, 3C-SiC and 2H-SiC are two typical structures with only k or h stacking sequence. This means that 3C-SiC has only cubic structure component, while 2H-SiC has only hexagonal one. The oxidation behavior of these polytypes may be influenced by the amount of k or h component which causes differences in bonding strength and further induces different adsorption barriers of oxygen. As the results given by density functional theory (DFT) calculation, [20] the adsorption energy of oxygen to Si-terminated surface of 3C-SiC is approximately 0.2 eV less than that of 2H-SiC. The surface configurations are also different depending on the surface sequence. 3C-SiC has larger relaxed bond angle than that of 2H, which means that cubic structure may experience inward relaxation and can be more flexible than 2H while 2H has a closer structure to sp3. This sp3 character gives rise to a stronger bonding between adsorbed oxygen and Si atoms on surface for 2H and makes it more difficult for Si– O bond formation. Such kind of state can facilitate the initial oxygen adsorption to Si-terminated surface of 3C-SiC and further increase the oxidation rate of 3C bulk. As to C-terminated SiC surfaces, all the adsorption energy values are much lower than that of Si-terminated surfaces for both 3C and 2H, which indicates that the polytype effect may diminish whereas the influence of top layer can become the dominating factor.

Apart from the analysis above, the configurations containing only two neighboring bonding tetrahedrons for 3C and 2H stacking are shown in the bottom half of Fig. 9. The adjacent layer of 3C is twisted by 60° with respect to the atom position of 2H. Due to such twisting, the distance between Si atom and C atom belong to the neighboring stacking (shown by the arrow) for 3C is slightly larger than that of 2H. Thus, 3C-SiC may provide a wider space for oxygen, correspondingly holding more oxygen molecules and enhancing the reaction. This can partially lead to the higher oxidation rate of 3C than that of 2H-SiC.

Fig. 9. Comparison between the micro-structures of 3C- and 2H-SiC.

4. Conclusions

The anisotropy property of SiC oxidation is of great importance for both theoretical and technical development, which was seldom studied at a high temperature range of above 2000 K. With the aid of ReaxFF MD, the structure evolution and oxygen consumption at initial oxidation stage for , Si(100), , and Si(111) of 3C-SiC were investigated in this current study. It was found that over 2000 K, all the four cases completely consumed oxygen molecules at 80 ps, producing a mixture of carbon chain and SiOx in the middle part of the configurations. By further data fitting, Arrhenius parameters were obtained. The activation energy of oxidation was the lowest among the four polar faces, which can indicate inferior oxidation resistance. In addition, the values of activation energy for 3C-SiC were computed and they are in the range of 30– 35 kJ· mol− 1, which were lower than that for 2H-SiC. This implies weaker stability of the cubic SiC surrounded by oxygen. It is believed that the polar face effects of 3C-SiC oxidation presented in this paper can find applications in developing future high-temperature materials.

Acknowledgments

We thank Prof. Satoshi Izumi of the University of Tokyo for valuable suggestions.

Reference
1 Park C H, Cheong B H, Lee K H and Chang K J 1994 Phys. Rev. B 49 4485 DOI:10.1103/PhysRevB.49.4485 [Cited within:1]
2 Zhang X D, Cui S X and Shi H F 2014 Chin. Phys. Lett. 31 016401 DOI:10.1088/0256-307X/31/1/016401 [Cited within:1]
3 Dou Y K, Qi X, Jin H B, Cao M S, Usman Z and Hou Z L 2012 Chin. Phys. Lett. 29 077701 DOI:10.1088/0256-307X/29/7/077701 [Cited within:1]
4 Opila E J and Hann R E 1997 J. Am. Ceram. Soc. 80 197 DOI:10.1111/j.1151-2916.1997.tb02810.x [Cited within:1]
5 Opila E J 1999 J. Am. Ceram. Soc. 82 625 DOI:10.1111/j.1151-2916.1999.tb01810.x [Cited within:1]
6 Williams S, Curry D M, Chao D C and Pham V T 1995 J. Thermophys. Heat Transfer 9 478 DOI:10.2514/3.690 [Cited within:1]
7 Li X H, Yan Q Z, Mi Y Y, Han Y J, Wen X and Ge C C 2015 Chin. Phys. B 24 026103 DOI:10.1088/1674-1056/24/2/026103 [Cited within:1]
8 Powell J, Petit J, Edgar J, Jenkins I, Matus L, Choyke W, Clemen L, Yoganathan M, Yang J and Pirouz P 1991 Appl. Phys. Lett. 59 183 DOI:10.1063/1.105960 [Cited within:1]
9 Gupta S K and Akhtar J 2011 Silicon Carbide-Materials, Processing and Applications in Electronic Devices(InTech) pp. 207 230 DOI:10.5772/20465 [Cited within:2]
10 Deng J, Liu W, Du H, Cheng H and Li Y 2001 J. Mater. Sci. Technol. 17 543 [Cited within:1]
11 Amy F and Chabal Y J 2003 J. Chem. Phys. 119 6201 DOI:10.1063/1.1602052 [Cited within:1]
12 Shishkin Y, Oborina E, Maltsev A, Saddow S and Hoff A 2006 J. Phys. D: Appl. Phys. 39 2692 DOI:10.1088/0022-3727/39/13/009 [Cited within:1]
13 Lilov S 2007 Cryst. Res. Technol. 42 1054 DOI:10.1002/(ISSN)1521-4079 [Cited within:1]
14 Yamamoto T, Hijikata Y, Yaguchi H and Yoshida S 2008 Jpn. J. Appl. Phys. 47 7803 DOI:10.1143/JJAP.47.7803 [Cited within:1]
15 Hijikata Y, Yaguchi H and Yoshida S 2009 Appl. Phys. Express 2 021203 DOI:10.1143/APEX.2.021203 [Cited within:1]
16 Zhang Y, Li H, Qiang X and Li K 2010 J. Mater. Sci. Technol. 26 1139 [Cited within:1]
17 Yao X, Li H, Zhang Y and Wang Y 2014 J. Mater. Sci. Technol. 30 123 DOI:10.1016/j.jmst.2013.09.006 [Cited within:1]
18 Suo T, Fan X, Hu G, Li Y, Tang Z and Xue P 2013 Carbon 62 481 DOI:10.1016/j.carbon.2013.06.044 [Cited within:1]
19 Wang J, Zhang L, Zeng Q, Vignoles G L, Cheng L and Guette A 2009 Phys. Rev. B 79 125304 DOI:10.1103/PhysRevB.79.125304 [Cited within:2]
20 Wang J, Zhang L, Zeng Q, Vignoles G L and Cheng L 2010 J. Phys. Condens. Matter 22 265003 DOI:10.1088/0953-8984/22/26/265003 [Cited within:1]
21 Liu Y, Su K H, Zeng Q F, Cheng L F and Zhang L T 2012 Theor. Chem. Acc. 3 1 DOI:10.1021/jp3022503 [Cited within:1]
22 Devynck F, Giustino F, Broqvist P and Pasquarello A 2007 Phys. Rev. B 76 075351 DOI:10.1103/PhysRevB.76.075351 [Cited within:1]
23 Deák P, Gali A, Knaup J, Hajnal Z, Frauenheim T, Ordejón P and Choyke J 2003 Physica B 340 1069 DOI:10.1016/j.physb.2003.09.252 [Cited within:1]
24 He X M, Chen Z M and Li L B 2015 Chin. Phys. Lett. 32 036801 DOI:10.1088/0256-307X/32/3/036801 [Cited within:1]
25 Ye Y, Zhang L, Cheng L and Xu Y 2003 J. Mater. Sci. Technol. 19 29 [Cited within:1]
26 Cheng T, Wen Y and Hawk J A 2014 J. Phys. Chem. C 118 1269 DOI:10.1021/jp409811e [Cited within:1]
27 Sun J B, Tang X Y, Yang Z W, Shi Y and Zhao Y 2014 Chin. Phys. B 23 066103 DOI:10.1088/1674-1056/23/6/066103 [Cited within:1]
28 van Duin A C T, Dasgupta S, Lorant F and Goddard W A 2001 J. Phys. Chem. A 105 9396 DOI:10.1021/jp004368u [Cited within:1]
29 Newsome D A, Sengupta D, Foroutan H, Russo M F and van Duin A C T 2012 J. Phys. Chem. C 116 16111 DOI:10.1021/jp306391p [Cited within:6]
30 Newsome D A, Sengupta D and van Duin A C T 2013 J. Phys. Chem. C 117 5014 DOI:10.1021/jp307680t [Cited within:4]
31 Song Y, Dhar S, Feldman L, Chung G and Williams J 2004 J. Appl. Phys. 95 4953 DOI:10.1063/1.1690097 [Cited within:1]
32 Aktulga H M, Fogarty J C, Pand it S A and Grama A Y 2012 Parallel Comput. 38 245 DOI:10.1016/j.parco.2011.08.005 [Cited within:1]
33 Plimpton S 1995 J. Comput. Phys. 117 1 DOI:10.1006/jcph.1995.1039 [Cited within:1]
34 Castro-Marcano F, Kamat A M, Russo Jr M F, van Duin A C T and Mathews J P 2012 Combust. Flame 159 1272 DOI:10.1016/j.combustflame.2011.10.022 [Cited within:1]
35 Rys A, Singh N and Cameron M 1995 J. Electrochem. Soc. 142 1318 DOI:10.1149/1.2044170 [Cited within:1]
36 Chen X, Ning L, Wang Y, Li J, Xu X, Hu X and Jiang M 2009 J. Mater. Sci. Technol. 25 115 [Cited within:1]
37 Charpentier L, Balat-Pichelin M and Audubert F 2010 J. Eur. Ceram. Soc. 30 2653 DOI:10.1016/j.jeurceramsoc.2010.04.025 [Cited within:1]
38 Charpentier L, Balat-Pichelin M, Glénat H, Bêche E, Laborde E and Audubert F 2010 J. Eur. Ceram. Soc. 30 2661 DOI:10.1016/j.jeurceramsoc.2010.04.031 [Cited within:1]
39 Kim Y, Min K, Shim J and Kim D J 2012 J. Eur. Ceram. Soc. 32 3611 DOI:10.1016/j.jeurceramsoc.2012.04.044 [Cited within:1]
40 Jagodzinski H 1949 Acta Crystallogr. 2 201 DOI:10.1107/S0365110X49000552 [Cited within:1]