Role of generated free radicals in synthesis of amorphous hydrogenated boron carbide from orthocarborane using argon bombardment: a ReaxFF molecular dynamics study

In this study, we modeled and analyzed a critical aspect of the synthesis process for a-BxC:Hy i.e. the argon bombardment from the ortho-carborane precursor. Utilizing the reactive molecular dynamics simulations, we identified and evaluated the formation of free radicals as a result of ion bombardment. We found that increasing kinetic energy of ions releases an increasing amount of hydrogen. Thus, hydrogen content in the product can be tuned by varying the ion energy. Overall, our approach allows for a better understanding of the mechanism of Ar bombardment and the role of radical species toward the formation of ortho-carborane network.


Introduction
High strength ceramics such as boron carbide and silicon carbide can potentially be used for a wide range of structural applications. These ceramics in general have a higher mechanical strength than typical metals [1][2][3][4][5][6]. Boron carbide, for example, is the third hardest material in nature following diamond and cubic boron nitride [1]. Its combined high hardness and very low density offer an attractive materials choice, especially when the strength to weight ratio is taken into consideration. In addition, it is thermally stable with a very high melting temperature, above 2500 K. Combined with its high abrasion resistance [1], the thermal stability and strength of boron carbide amount to an attractive set of physical properties well-suited for body armor components and abrasives. Since boron-10 has a high cross-section for neutron absorption, boron carbide is also an important material for nuclear reactors [7,8] and neutron detection [9,10]. Lastly, within the semiconductor industry, thin-film boron carbide can also be used for insulating low-dielectric-constant (low-k) intra/interlayer dielectrics [11], patterning materials for semiconductor devices, optical films, and other specialized electronic applications [12][13][14][15].
There have been a wide variety of growth techniques developed to produce boron-carbide-based solids [1,2] including hot pressing and pressureless sintering [16][17][18]. These synthesis routes generally require high temperatures and have been used primarily to produce boron carbide for structural/engineering applications such as armor and abrasives. However, these bulk processing techniques can not be easily integrated with semiconductor device processing and oftentimes do not yield desirable electronic properties [19,20]. Alternatively, chemical vapor deposition (CVD) or physical vapor deposition (PVD) techniques have been found to be quite useful in generating B x C thin films with mechanical, thermal, and chemical properties comparable to those of the bulk materials. Furthermore, their electronic properties can be tuned and enhanced [21][22][23][24][25]. In particular, thin films derived from ortho-carborane (o-C 2 B 10 H 12 ) using a plasma-enhanced chemical vapor deposition (PECVD) technique have been demonstrated to possess promising properties for a variety of device applications [26][27][28][29][30][31]. This fabrication approach produces an amorphous hydrogenated boron carbide (a-B x C:H y ) variant, which is quite novel in the field of electronic materials. The electronic, optical, chemical, and mechanical properties of these films have been extensively investigated, and in particular the influence of hydrogen on these properties [32]. However, it is still an open question as to which types of free radicals form during the PECVD process as well as how the ionic bombardment contributes to the overall formation of orthocarborane aggregates leading to the formation of the thin films and the overall hydrogen retained in the product. In this study, we attempt to gain a more fundamental insight into this process through a systematic investigation of the impact of argon bombardment on ortho-carborane in the gaseous phase.

Computational approach
Since a single molecule of orthocarborane is already comprised of 24 atoms [32,33], modeling an assembly of hundreds of these molecules would be quite difficult if we were to utilize a quantum mechanical (QM) approach [34][35][36][37]. The maximum number of atoms that can be modeled by means of ab-initio calculations is limited to only a few hundred atoms. In order to overcome this limitation while maintaining nearly the same level of accuracy as achieved by QM calculations, we thus opted instead to perform reactive molecular dynamics (MD) simulations by using a carefully parameterized empirical ReaxFF reactive force field method [38,39]. The reactive molecular dynamics simulation can significantly reduce simulation time by several orders of magnitude and is capable of defining continuous bond formation/breaking events during the simulations by monitoring the bond orders and energies of all of the pairs within the simulation [38][39][40][41][42][43]. For this study, we have used the ReaxFF force field developed by An et al, [44] for which parameters have been carefully parametrized to model crystalline and amorphous materials made of boron, carbon, and hydrogen atoms. In An et al, [44] though QM interactions of two C 2 B 10 H 12 molecules (icosahedral ortho-carborane or one of its isomers) were included in the ReaxFF fitting parameters, the application of that force field is shown mainly for boron-carbide-based material. Therefore, before proceeding, we determined the applicability of ReaxFF by comparing some results as calculated by QM and ReaxFF. This includes BH and CH bond dissociation energies in the ortho-carborane molecule as well as the total energies of different common species of boron, carbon, and hydrogen relevant to our study.
To simulate the Ar ion bombardments, Lennard Jones (LJ) potentials have been commonly used, as during bombardment the consideration is mostly the transfer of kinetic energy from high-energy ions to the bombarded species [45,46]. For our studies, we have also used the LJ potential to represent the pair interactions between Ar with the rest of the species. Each of the two LJ interaction parameters (σ, which is correlated to equilibrium distance of paired atoms and ε, which is associated with the minimum energy) for argon-argon and argon-carbon were obtained from Norio Inui et al [46] They modeled the dynamic impacts of argon clusters onto a free-standing graphene monolayer. The optimized LJ parameters from their study are: σ Ar-Ar =3.4 Å, ε Ar-Ar =0.23983 Kcal mol −1 , σ Ar-C =3.385 Å, and ε Ar-C =0.115303 Kcal mol −1 . Boron-boron LJ parameters were obtained from Baowan et al, [47] a study on boron nitride (BN) and carbon boron nitride (CBN) nanocones. Finally, hydrogen LJ parameters were obtained from thermophysical studies of H 2 O and argon plasma [48]. Using the Lorentz-Berthelot combination [49,50], we calculated the LJ parameters, namely σ Ar-B =3.4265 Å, ε Ar-B =0.150885 Kcal mol −1 , σ Ar-H =3.054 Å, and ε Ar-H =0.13283 Kcal mol −1 . According to the Lorentz-Berthelot combination rule, we can define σ 1-2 =(σ 1-1 +σ 2-2 )/2 and ε 12 =(ε 1 ε 2 ) 1/2 . For all MD simulations, we implemented the same and sufficiently large cutoff radius of 10 Å. For the MD simulations, we used the highly efficient Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) MD code [51], which was developed by Sandia National Laboratories and has been successfully used to perform classical MD events on a wide range of materials [52][53][54][55][56][57][58].

Simulation details
3.1. Quantum mechanics calculations As previously mentioned, for determining the applicability of the force field, we calculated the bond dissociation energies of BH and CH in the ortho-carborane molecule and total energies of different relevant species (both chemically active and stable) of B, C, and H such as ortho-carborane, borane, and diborane using both ReaxFF and QM calculations. The QM calculations were performed using the commercial Vienna Ab-initio Simulation Package (VASP) [59]. Plane wave basis set and projector-augmented wave (PAW) based pseudo-potentials are implemented in the VASP software [59][60][61][62]. The ultra-soft pseudopotential with the Perdew-Burke-Ernzerhof (PBE) version of the generalized gradient approximation (GGA) exchange-correlation functional of DFT is adopted in substitution of core nuclei and valence electron exchange-correlation term calculation. We used the Monkhorst-Pack [63] scheme to create the k-space matrix. The plane wave basis set energy cut off was set to 520 eVm, and the tetrahedron method with Blöchl corrections smearing was selected with sigma value 0.2.

Argon bombardment of ortho-carborane molecules
Our primary focus was to study the impact of argon bombardment on ortho-carborane in the gaseous phase. In order to create the gaseous mixture, 1000 ortho-carborane molecules and 25 Ar atoms were packed into a small cubic cell with a unit cell of 102 Å. PACKMOL code [64] was used to randomly populate these species within the box. We also sampled one larger simulation box (8×the size of the initial box) comprising 8000 orthocarborane molecules and 200 Ar atoms to provide better statistics on the gaseous products during the bombardment processes. The kinetic energy of the argon atoms was varied from 100 to 600 eV in 50 eV increments. For our interatomic potentials, 50 eV argon energy is not enough to break ortho-carborane, while 600 eV is set as the upper limit of the theoretical studies of argon energy distribution similar to previous studies [65,66]. The direction of the high-kinetic-energy argon atoms was initially kept random in order to mimic the interaction in the gaseous phase. The bombardment was performed for 5 ps with 0.05 fs time steps using the constant volume and energy simulation method (NVE ensemble). We chose to adopt a relatively small time step [67] due to the high kinetic energy. The whole simulation cell was subsequently equilibrated at 300 K using the constant volume and temperature method (NVT) for 100 ps with 0.25 fs time steps. We imposed periodic boundary conditions in all axes. In the NVT ensembles, the temperature of the box was controlled via a Nose-Hoover thermostat.
In order to analyze the impact of argon bombardment, we started the analysis by examining a pair distribution function (PDF) plots of different atomic pairs, namely BH, HH, and CH. We carried out the time evolution of the species count of those that were created at the initiation of the NVE simulation. To identify and then quatify the number of bond pairs, we used the distance beyond the first peak observed in PDF as the only criteria to mark the cut off radius for bond pairs. Consequently, the distance for a given pair will slightly vary (less than 15%) depending on the simulation ensembles used (e.g. NVT or NVE) and the simulation temperatures.
We found tens of different types of radicals and small clusters during bombardment, including both chemically stable and unstable species. We recorded the statistics of the most common species every 0.5 fs for the first 5 ps of the Ar bombardment. The data was then used to quantify correlations among these. To calculate the correlations we used the machine learning software WEKA [68]. WEKA is an open-source code embedded with several machine learning algorithms. It has been widely used for data mining purposes [69][70][71][72]. It also features a collection of mathematical approaches such as linear regression or multilayer perception to study correlations. We have used the linear regression method for our study.

ReaxFF verification
The reactive force field we opted to use has been previously parameterized through QM calculations to reproduce structure, energy, and reaction barriers [44]. It has been successfully used to analyze the mechanisms of brittle failure of non-hydrogenated crystalline and amorphous boron carbide. Here, we initially begin by testing the applicability of the ReaxFF potential in the analysis of ortho-carborane, a hydrogenated BC-based precursor. Figure 1 shows the bond dissociation energy curves of BH (a) and CH (b) from ReaxFF and DFT calculations. For both bonds the distances were taken up to 5 Å. Figure 1 indicates that the ReaxFF data are in good agreement with the DFT results. Additionally, table 1 compares the formation energy of different species using DFT and ReaxFF. While table 1 is far from complete, it shows that ReaxFF yields results fairly similar to those produced by DFT. However, we do notice that the formation of the B 2 H 6 dimer is much more favorable in the case of DFT. From the DFT calculations, the energy of formation, ΔH, of the reaction 2BH 3 →B 2 H 6 is calculated as -50.54 Kcal mol −1 , whereas for ReaxFF it is calculated as -11.58 Kcal mol −1 . The calculated DFT value of -50.54 Kcal mol −1 is slightly higher than the experimental value of -41.50 Kcal mol −1 [73]. This indicates that ReaxFF would have underestimated the stability of B 2 H 6 , and indeed this was observed in our MD simulation results as discussed in section 4.2.2. Nevertheless, since both BH 3 and B 2 H 6 are gaseous products [74,75] with the same relative hydrogen percentage, this should have a minimal impact on the primary focus of our analysis, which is the hydrogen content in the solid clusters.
In additional testing, we demonstrated that ReaxFF was able to produce a stable structure for the solid crystalline phase of ortho-carborane. As shown in figure 2, we were able to reproduce the FCC crystal structure of ortho-carborane at room temperature with a unit cell dimension of ∼10.9 Å, which is in good agreement with the experimental findings [76].

Argon bombardment
As illustrated in figure 3, which represents a trajectory of the ortho-carborane molecule almost immediately after its collision with energetic argon, the kinetic energy from Ar is transferred into the molecule and as a consequence, we observe the formation of a variety of species. In figure 3, for example, we observe the formation of a cluster consisting of seven B atoms, one C atom, and six H atoms, a BH 3 molecule, BH radicals, and isolated carbon and hydrogen radicals.

Pair distribution function analysis
In order to analyze the formation of species during argon bombardment in greater detail, we took the structure data of the simulation cell at different times and performed a pair distribution function (PDF) analysis. In figure 4, we present the PDF analysis of BH pairs at different times of the simulation with 25 argon atoms of 600 eV kinetic energy. The dominant peaks at 0 ps are from the different orders of the nearest neighbors from the initial ortho-carborane molecule. Due to the presence of carbon atoms in the molecule, the structure is slightly distorted from a regularly shaped B-only icosahedron because of the assymetry in bond length generated by the C atoms. In addition, there are some slight variations in the BH and CH bond lengths at the initial setting prior to bombardment noted in figure 4. Because of Ar bombardments, some of the icosahedral structures would fracture causing a wider range of interatomic distances for each pair as evident from the gradual attenuation/ broadening of peaks, also seen in figure 4. PDF analysis of HH pairs further supports this observation as we can again see the gradual diminishing of initially sharp peaks, as shown in figure 5. These initial peaks originated from adjacent H atoms within ideal icosahedral structures, initially bound to either B or C atoms within the ortho-carborane. Within the initial crystalline structure, no H-H bonds were present, as shown in figure 2. After argon bombardment, we observe two newly formed peaks at ∼0.75 Å and ∼2.1 Å, which can be attributed to the formation of hydrogen molecules and BH 3 molecules, respectively [77]. At the end of the NVE simulations, the  temperature remained high and thus the PDF remain fairly diffused. However, after quenching the structure to room temperature by means of NPT canonical ensembles, we would then observe one dominant and fairly sharp peak below 2 Å indicating the formation of stable hydrogen molecules.

Species calculation
We showed earlier that, when a high-energy argon atom collides with an ortho-carborane molecule, different kinds of species and clusters form. Here we quantitatively and qualitatively investigate these species. We will refer to each simulation by the kinetic energy of the argon atoms set at the beginning of the bombardment. For instance, the particular simulation for which the bombardment was started with assigning 400 eV kinetic energy to all of the 25 argon atoms will be referred to as a 400 eV simulation. We have calculated bond order values to determine the presence of chemical bonds. The bond order calculations were performed using the built-in method of ReaxFF in LAMMPS (reax/c/species). In the beginning, we put 1000 ortho-carborane molecules in all of our simulations. Figure 6 shows the most common species found including small clusters, free radicals, ortho-carborane, and bigger clusters at different times for the 600 eV simulation. With the the use of small clusters, the dynamics involving the statistics of the main species can be assessed. It does not however necessarily capture the statistical dynamics of the smaller species population. Using 600 eV as the highest kinetic energy we have used for argon, we found a considerable amount of hydrogen generated in the simulation cell after 1 ps of the simulation time. Figure 6 indicates that most of the hydrogen was produced within 1 ps and it grew slowly but consistently over time. More than 300 molecules of hydrogen (more than 5% of the total hydrogen) were found after 5 ps of bombardment and the quantity was unchanged after equilibration at room temperature. Initially, lots of H radicals were formed, which are chemically very active and unstable. The number of H atoms decreased drastically and almost disappeared by the end of the simulation. This can be explained by the fact that two hydrogen atoms easily react to produce H 2 molecules. Other active BH and BH 2 free radicals were found to form after 1 ps of simulation due to the fragmentation of ortho-carborane molecules by the energetic Ar impact. In the experimental Ar plasma environment, free radicals such as H, BH, CH, and CH 2 have been reported as the main products during the plasma compositional analysis of PECVD with trimethylboron [78] and triethylboron [79] precursors. In our modelling process using ortho-carborane, the dominant BH, BH 2 , CH, and atomic H species are similarly found to be formed right after argon ion impact. The BH and BH 2 species, however, are found to subsequently disappear after equilibration, which can be explained by the fact that some of them readily form bonds with H radicals and transform into BH 3 gas. This is consistent with a previous Ar bombardment study [80], which concluded that the formation of BH 3 gas was the favorable pathway for PECVD chemistry in both Ar and H 2 atmospheres with triethylboron as precursor. The formation of these gaseous species could be the possible reason for the lack of any extraicosahedral B environments detected in the PECVD experimentally derived thin films [81]. This may imply that while favorable subsequent reactions involving BH, BH 2 and then BH 3 gaseous species should indeed readily take place in the presence of the H radicals, these remaining volatile species will likely be removed from the thin films. Thus, the resulting thin film would be slighly depleted from B which was quite consistent woth the experimental finding. However, we should however note that there have been reports of both B and C types of extraicosahderal environment found in films derived from orthocarborane precursors subject to different processing environments [82,83].
The trend for the amount of BH 3 molecules was found to be very similar to that of the H 2 molecules, whereas the variation in the concentration of BH and BH 2 species tends to follow the trend of the variation in the H concentration. Metastable species like BCH, B 2 CH 2 , B 3 CH 3 , B 2 CH 4 , B 3 CH 4 , and others are found in small numbers. Most likely, clusters like BCH and BCH 2 were formed from direct collisions between Ar and orthocarborane. Others may have been created from isolated B, C, and H atoms/radicals in the vicinity of BCH and BCH 2 . In analyzing the trend in behavior of the bigger clusters, we observed that the number of orthocarborane molecules dropped to around 33% of the initial count within 1 ps. Deprotonation of the orthocarborane molecules was observed with the evolution of icosahedral clusters including B 10 C 2 H 11 , B 10 C 2 H 10 , and B 10 C 2 H 9 . The formation of these partially dehydrogenated species of ortho-carborane is in good agreement with the PECVD radicals formed during the Ar plasma environment [84]. In addition, some of the free radicals generated by bombardments such as B, BH, C, and CH are relatively unstable and thus readily react with H radicals in the vicinity of the surrounding ortho-carborane molecules. Figure 7 shows that for the 100 eV simulation, around 10% of the ortho-carborane molecules were fragmented after 5 ps, whereas more than 800, i.e. 80%, were fragmented for the 600 eV simulation. Figure 8 shows the time evolution of different species created for different simulation energies. For the 600 eV simulation, the highest number of BH and CH free radicals was observed within 200-300 fs. BH was found to be the most common type of free radical at the beginning, aside from H. After peaking at around ∼200-300 fs, the quantity decreased dramatically, presumably due to its reaction with the surrounding species. The amount of BH 2 radicals, which also initially peaked, shows a similar decrease with time, albeit slightly more gradual. Presumably, this is due to the difficulty, over time, to find remaining H radicals to transition to BH 3 . In general, the amount of CH n is low. This may indicate the difficulty in stabilizing CH n species relative to their BH n  counterparts. In addition, we observed smaller amounts of CH 3 and CH 4 species relative to CH and CH 2 . Presumably, this is due to the majority amount of C-related radicals generated by bombardment have been interconnected with the surrounding ortho-carborane molecules. This may also be due to the relatively smaller amount of CH (relative to BH), limiting the possibily for further reactions with H radicals. Our finding of C atoms in extra-icosahedral environment in the form of CH 2 and CH 3 has been reported in a number of PECVD film depositions [81,[85][86][87][88]. BH 3 gas molecules, on the other hand, were highest among all these species at the end of the bombardment simulations and consistently increased over time. This can be understood by the fact that BH 3 is a relatively stable molecule, forming a 2D planar structure [77]. Without a sufficient external kinetic disturbance, these molecules should remain stable once formed. In figure 9, the amount of hydrogen produced over time is shown. The trend is similar with BH 3 since H 2 is also a stable molecule and, once formed, will remain chemically stable. Our observation of the formation of stable BH 3 and H 2 is in agreement with the plasma reaction using the single solid precursor triethylboron [80]. Another previous study by Imam et al [79] showed the dominant formation of BH and H 2 in the plasma. This was proposed to result from the dehydrogenation of BH 3 in the presence of excess H radicals. In our MD simulation study, since the H radicals have been primarily consumed by H 2 gas formation at the early stage of the simulated Ar bombardment, we did not observe the proposed chemical reaction.
Overall, the correlation between the initial kinetic energy of argon and the number of species generated can be quantitatively evaluated through our MD simulations. In addition, the hydrogen content in the formed clusters is in good agreement with previous studies [32,81,89].

Species correlation: regression analysis
Based off MD simulation results using 1000 ortho-carborane and 25 Ar atoms, we were able to identify the main species formed and evolved during the processes. To improve the statistics, we made a new simulation cell with 8000 ortho-carborane molecules and 200 argon atoms. The density of the gaseous mixture (i.e., the ratio of orthocarborane to Ar) was kept the same as before. For the sake of representation, we categorized the important species in four different groups based on the combination of types of atoms. Figure 10 shows the time evolution of these groups of species for 600 eV argon energy. The outcomes from this new structure are consistent with those produced by smaller scale simulations. The number of active species like B, BH, H, C, CH, CH 2 , and BCH were grown at the initiation of the ion bombardment. However, their number declined quickly as they started creating bonds with other species and clusters in their vicinity. Less active and metastable species like BH 2 , CH 3 , BCH 2 declined comparatively slowly. Gradual protonation showed the following pattern: B>BH>BH 2 >BH 3 , C>CH>CH 2 >CH 3 >CH 4 , BCH>BCH 2 >BCH 3 >BCH 4 and so on. From group D the trend of originated clusters from ortho-carborane is observed. It should be noted that the Y axis for group D is in log scale. Mostly, orthocarborane deprotonated and generated clusters such as B 10 C 2 H 11 , B 10 C 2 H 10 , and B 10 C 2 H 9 .
In figures 11(a) and (b) the actual and predicted data of BH 3 and isolated H atoms are shown respectively. We chose BH 3 as stable species and H as active species. The correlations between stable species with BH 3 , and active species with H were analyzed extensively. From the visualization, we see the predicted data, using the linear correlation, is in good agreement with the actual simulated data. The linear regression models to predict the amount of BH 3 and H radicals are as follow:   Table 2 shows a summary of the statistical outcomes of the regression models. If we look carefully to the regression model and the summary in the table, stable species BH 3 is highly correlated with other stable species. More importantly, this correlation is linear. Similarly, active isolated hydrogen atoms are linearly correlated with other active species like BH, CH, and BCH. Other correlations can be found from the available data. It should be noted that there can be other species not specified here which might not be significant individually but altogether could have a nontrivial impact. Further detailed study will be needed using larger simulation scales. Nevertheless, these types of correlations might be useful to deduce some of the unmeasurable (low count) species based of the quantity of the measurable species.
Overall, our systematic modeling study confirms some of the observations obtained experimentally [32]. For example, the previous experimental results did indicate that when the percentage of hydrogen content is relatively high, the density is normally low and vice versa [79,87]. Additionally, we have also been able to show the evolution of the radical or stable species as a function of the Ar bombardment. These modeling results are  useful for further strategizing the synthesis process so as to maximize the density and quality of the boron carbide thin films.

Conclusion
We have successfully modeled the argon bombardment process of ortho-carborane using a reactive force field (ReaxFF) molecular dynamics simulation. By varying the energy of argon atoms, we systematically determined the amount of hydrogen gas and various types of gas species produced. The time evolution of PDF analysis allowed us to monitor the decomposition of ortho-carborane molecules as a result of bombardment. In addition, we have successfully identified the dynamics in the free radical formation.

Acknowledgments
We are grateful for the support of NSF (DMREF) Grant no 1729176. The computational work was performed at the National Energy Research Scientific Computing Center (NERSC) under ERCAP0007777 computational project.