Models and regressions to describe primary damage in silicon carbide

Silicon carbide (SiC) and SiC/SiC composites are important candidate materials for use in the nuclear industry. Coarse grain models are the only tools capable of modelling defect accumulation under different irradiation conditions at a realistic time and length scale. The core of any such model is the so-called “source term”, which is described by the primary damage. In the present work, classical molecular dynamics (MD), binary collision approximation (BCA) and NRT model are applied to describe collision cascades in 3C-SiC with primary knock-on atom (PKA) energy in the range 1–100 keV. As such, BCA and NRT are benchmarked against MD. Particular care was taken to account for electronic stopping and the use of a threshold displacement energy consistent with density functional theory and experiment. Models and regressions are developed to characterize the primary damage in terms of number of stable Frenkel pairs and their cluster size distribution, anti-sites, and defect type. As such, an accurate cascade database is developed with simple descriptors. One of the main results shows that the defect cluster size distribution follows the geometric distribution rather than a power law.

Silicon carbide (SiC) and SiC/SiC composites are candidate materials for use in fission and fusion applications. SiC is an attractive material due to its favorable neutronic and high-temperature properties. In fusion applications, SiC/SiC composites are considered as structural and/or functional material in the tritium breeding blanket module [1][2][3] . In fission applications, SiC and SiC/SiC composites are considered as cladding material, ranging from coated particle fuels to advanced accident tolerant fuel claddings in light water reactors 1,2 .
Under irradiation, the attractive properties of SiC are known to degrade as the material's microstructure is modified. Depending on the irradiation conditions, point-defects accumulate into defect clusters that may provoke macroscopic swelling and amorphization of the material 2 .
Snead and Katoh 4 developed a comprehensive map of the evolution of microstructural features in SiC under irradiation (both ion and neutron) with dose and irradiation temperature. Initially, damage accumulates in SiC in the form of black spot defects, i.e., defect clusters that are too small for adequate characterization. As either the temperature or dose (dpa) increases, interstitials gather to form faulted Frank loops. When both temperature and dose are increased, voids and unfaulted loops develop.
At nearly all temperatures, SiC swells, despite the lack of cavities at low temperatures and doses 4 . A comprehensive overview of the trend of swelling with temperature is given in 4 . Below a critical temperature (about 150 °C), swelling is primarily caused by amorphization of the SiC material as strain accumulates from irradiationinduced defects from either ion-or neutron irradiation [5][6][7][8] . In the intermediate temperature regime, above the critical amorphization temperature (about 200 °C) and up to the vacancy mobility temperature (about 1000 °C), the swelling of SiC under neutron irradiation appears to saturate with irradiation dose and decreases with irradiation temperature 4 .
The works by Katoh et al. 5 and Kishimoto et al. 9 demonstrated the effects of helium on the swelling in SiC and SiC/SiC composites. In both works, the addition of helium resulted in up to 1.5 times larger amounts of swelling at temperatures up to 1000 °C in SiC. At temperatures higher than 1000 °C, the addition of helium has a negligible effect on the swelling.
Theoretical works employing various density functional theory (DFT) methods have been applied to determine the formation and migration energy of carbon and silicon point-defects [10][11][12][13] . Furthermore, the clustering behavior of H and He, and their interaction with point-defect clusters was also studied by DFT [14][15][16][17][18] . From these works, it was concluded that vacancy (clusters) have a strong impact on the stability and mobility of He and H clusters.
Conversely, the mixed experimental and theoretical works by Daghbouj et al. 19,20 showed that He and H also have an important impact on the formation and stability of vacancy clusters filled with hydrogen/helium. It was found that vacancy clusters grow via Ostwald ripening (vacancy diffusion) and coalescence (because of the SCK CEN, Nuclear Materials Science Institute, Boeretang 200, B-2400, Mol, Belgium. ✉ e-mail: gbonny@sckcen.be open www.nature.com/scientificreports www.nature.com/scientificreports/ overlap of the strain fields around the small voids due to the pressure exerted by He or H 2 gas) where hydrogen and helium help stabilize the vacancy clusters.
Thus, clearly, regardless the irradiation type, the initial damage distribution has an important impact on the long-term damage evolution, especially when H and/or He impurities are present.
To model defect accumulation in SiC and SiC/SiC composites at realistic time and length scales, coarse grain models are necessary. In the literature, cluster dynamics 21 and object kinetic Monte Carlo models 22,23 have been developed to study defect accumulation in SiC. The starting point of any such model is the so-called "source term", i.e., primary damage.
In ceramics, the NRT model 24 and binary collision approximation (BCA) are accepted to be good first approximations to describe primary damage 25 . This is because in ceramics, unlike metals, most defects are formed as isolated Frenkel pairs with limited clustering. However, this approach does not provide information on the type of produced Frenkel pairs, C or Si vacancy/interstitial, number/type of anti-sites, nor size distribution of vacancy/ interstitial clusters. Although little clustering is expected, in the work by Liu et al. 21 , it was shown that a cluster dynamics model employing primary damage described as isolated Frenkel pairs was insufficient to reproduce the experimentally observed defect sizes and densities. It was argued that spatial correlations and clustering are likely necessary to obtain results consistent with experiment.
Therefore, the focus of the present work is to develop models and regressions to provide a detailed characterization of the primary damage in 3C-SiC (β-SiC). The primary damage is generated by simulating collision cascades employing classical molecular dynamics (MD), including electronic stopping, for a primary knock-on atom (PKA) energy in the wide range 1-100 keV. Although studies describing collision cascades in SiC are available in the literature [26][27][28][29][30][31][32] , in none of those electronic stopping is accounted for, which will be shown to play an important role.
Models and regressions are developed to characterize the primary damage in terms of number of generated stable Frenkel pairs and their cluster size distribution, anti-sites, and defect type, i.e., C/Si vacancies V C /V Si , C/Si interstitials I C /I Si and C/Si anti-sites C Si /Si C . Where possible, an in-depth comparison is made with data from the literature as well as with the binary collision approximation (BCA) 25,33 and the NRT model 24 . The latter two methods are computationally efficient and their validity compared to MD is assessed. Thus, the present study allows to fully characterize defect production in 3C-SiC using simple models and regressions.
The paper is organized as follows. Following this introduction, the computational details are summarized in section 2. In section 3, the stopping power is characterized such that electronic stopping is consistent with stopping and range of ions in matter (SRIM) tabulations 25,34 and the threshold displacement energy (TDE) is consistent to density functional theory (DFT) and experimental data. In section 4, the results from our MD simulations are presented and where appropriate compared to the literature, BCA simulations and the NRT model. In section 5, an analytical model based on the geometric distribution is developed to describe the defect cluster size distribution. The paper is finalized with conclusions in section 6.

computational Details
Classical molecular dynamics (MD) simulations were applied using the large-scale atomic molecular massively parallel simulator (LAMMPS) 35 . To describe the interatomic interactions in 3C-SiC, the potential developed by Erhart and Albe 36 (henceforth ERH) was used. This potential reproduces the DFT computed point defect energetics in SiC reasonably well (see also section 3), and also provides an excellent description of pure Si and graphite. This may be important for consistency with future studies considering SiC/SiC composites.
With the selected potential (ERH), cascade simulations were performed at zero Kelvin to construct a database where only ballistic effects are accounted for. As shown in the work by Samolyuk et al. 30 , the diffusion barriers of point-defects in SiC are not always described accurately by empirical potentials. By performing the simulations at zero Kelvin, inaccuracies on the resulting cascades due to an inaccurate description of the diffusion barriers are minimized.
Consistent with the work by Samolyuk et al. 30 , a variable timestep in the range 0.01-1 fs was used such that the maximum distance a particle can travel within a single step is 0.015 Å. This distance corresponds to the most probable distance travelled by a C atom during 1 fs (typical MD timestep) in a lattice equilibrated at 1600 K following the Maxwell-Boltzmann distribution.
To account for electronic stopping, a friction term consistent with the stopping and range of ions in matter (SRIM) tabulations 25 was added. Thereby communication between atomic and electronic lattice is neglected and only the net result on the atomic lattice is accounted for. The nuclear stopping was guaranteed by the universal screened Coulomb function, so-called "ZBL function", at short distances 25 .
Given the stiffened potential, the threshold displacement energy (TDE) for both Si and C is estimated in the [100], [110], [111] and [111] directions at zero Kelvin. To do so, the PKA energy was varied on a 1 eV grid using ten independent runs in a given direction. For this, a box containing 1728 atoms (6 × 6 × 6 a 0 3 ) was used. The lowest energy in any run that generated a stable Frenkel pair is taken as TDE.
To compute the formation energy of point-defect configurations, the same box size (6 × 6 × 6 a 0 3 ) is used. The configurations were relaxed using the conjugate gradient method at constant volume. The formation energy, E f , of defect configuration D was obtained following 36 , with E D ( ) the total energy of the box containing D, n i the number of atoms of type i, µ SiC 3C the chemical potential of 3C-SiC, µ Si diam the chemical potential of diamond Si and µ C grap the chemical potential of hexagonal graphite.
www.nature.com/scientificreports www.nature.com/scientificreports/ Cascades were simulated at constant equilibrium volume with a Si-PKA energy, E PKA Si , of 1, 5, 10, 50 and 100 keV in a high index direction <135>, consistent with other works 26,27,30 . Such a direction reduces channeling effects that lead to unnecessary complications. During channeling, the PKA and recoils disperse with minimum energy loss, typically along close packed directions. As a result of the channeling, the cascade may self-interact through the periodic boundaries.
Variations of a few percent on this ideal direction were applied to allow for statistics to be collected, i.e., for each E PKA Si , the cascade was run ten times. Simulation boxes of size 80 × 80 × 80 (4.096 × 10 6 atoms) for = E PKA Å the equilibrium lattice parameter following ERH. These simulation boxes proved to be large enough to avoid self-interaction of the cascade through the periodic boundaries (see also Supplementary Material). In all simulations, periodic boundary conditions were applied in three dimensions. The simulations were performed up to 10 ps, within which a stable number of Frenkel pairs were produced.
The simulation results (including TDE calculations) were post-processed to identify point defects using the Wigner-Seitz analysis as implemented in the open visualization tool (OVITO) 37,38 . As such, vacancies, interstitials and anti-sites could be identified. A cluster analysis was performed with a third nearest neighbor distance criterion.
Simulations using the binary collision approximation (BCA) were performed using the transport of ions in matter (TRIM) code 25,34 . The code was executed in "full cascade" mode, such that recoils were also followed in addition to the PKA. Cascades with E PKA Si = 0.5-500 keV were performed in bulk SiC for different TDEs. As TDE, the values obtained from the stiffened potential and DFT computed lower and upper bound 39,40 were used.

Stopping Power and Model Validation
The stopping power of a PKA in a medium consists of the electronic stopping power, S el , due to electronic excitations and nuclear stopping power, S nuc , due to the Coulomb interaction with the nuclei. To illustrate the importance of both contributions, the ratio S S / el nuc as provided by SRIM 25,34 is plotted with E PKA in Fig. 1a. As shown in the figure, ≥ S S el nuc from 10 keV and 70 keV for a C-PKA and Si-PKA, respectively. However, even for E PKA as low as 0.1 keV, S el is 10-20% of S nuc . Thus, for all considered cascade conditions S el is important and should be accounted for.
In the present work, S el was included by introduction of a friction term, such that the force experienced on atom i due to S el is proportional to its velocity, i i with γ a friction coefficient fitted to S el . In Fig. 1b . Within this approximation, communication between atomic and electronic lattice is neglected and only the net result on the atomic lattice is accounted for.
The nuclear stopping power is implemented by merging the screened Coulomb (ZBL) potential by Ziegler et al. 25,34 , V ZBL , to the equilibrium potential by Erhart and Albe 36 , V eq . To ensure a smooth connection between V ZBL and V eq , the scheme implemented in LAMMPS was used such that the resulting effective potential, V eff , is given as, www.nature.com/scientificreports www.nature.com/scientificreports/ r 0 95 C Å. It was verified that this parameter setting did not modify the equilibrium properties of the potential (e.g. defect formation energy).
To validate the simulation set-up and used potential, in the following, point-defect properties and the threshold displacement energy (TDE) are compared to available experimental and DFT data. For later reference, values obtained from the potential by Devanath et al. 26,41,42 (henceforth DEV) are also included in the comparison. In Section 4, cascade results obtained in the present work are compared to cascade results from the literature employing the DEV potential.
In Table 1, the chemical potential, formation energy of vacancies, anti-sites and various interstitial configurations obtained by the different potentials are compared to DFT values. The values for the most stable configurations are indicated in bold. For vacancies and anti-sites, both potentials reproduce the most stable DFT configurations, i.e., V C and C Si . Quantitatively, the values obtained by the potentials are of similar magnitude as the DFT ones.
For the C-interstitial, ERH predicts the C + -C <100> configuration as most stable configuration, consistent with DFT, while DEV predicts the C TSi as most stable. For the Si-interstitial, both ERH and DEV predict the Si + -C < 100> configuration as most stable configuration, which is inconsistent with DFT that predicts Si TC as most stable. Quantitatively, there is a large variation in magnitude between the different DFT data sets. Within this variation, the values obtained by the potentials are of similar magnitude as the DFT ones.
Thus, based on the point-defect stabilities presented in Table 1, ERH is slightly more consistent with DFT data compared to DEV. However, during high energy cascades, defect creation is expected to be dominated by ballistic effects. An important parameter for this process is the threshold displacement energy, which is discussed in the following.
The threshold displacement energy (TDE) for both Si and C is computed in the [ For later reference, the values for the DEV potential 26,41,42 as reported in 47 are also included as some results will be compared to cascade results obtained with that potential. The values of that potential are higher than ERH, but still in a similar range as the DFT data.

Config. DFT/Exp (eV) ERH (eV) DEV (eV)
Chemical potential www.nature.com/scientificreports www.nature.com/scientificreports/ As shown in the figure, the effect of TDE is limited, i.e., less than a factor two in total number of produced Frenkel pairs. The differences between NRT, TRIM and MD are also limited, i.e., less than a factor two. The values from TRIM simulations and NRT formula are essentially equivalent to the ones from the MD simulations. For For comparison, the data from the study by Cowen et al. 47 are added to Fig. 2. In that work, the potential by Devanathan et al. 42 is used and its TDE is summarized in Table 1. As shown in Table 1, the TDE is higher than for ERH and close to the largest values predicted by DFT. Therefore, it is striking that the number of observed Frenkel pairs is systematically higher than our MD, TRIM and NRT results. Clearly, the increased number of observed Frenkel pairs is related to the absence of electronic stopping in the simulation set-up by Cowen et al. 47 . This highlights the importance of accounting for electronic stopping in the simulations set-up.

High energy cascades
In Fig. 3, the ratio of V C /V Si and I C /I Si obtained from MD for all investigated E PKA Si are summarized and where possible, compared to TRIM and the NRT model. It is noted that the large error bar for 1 keV cascades is due to the low statistics, i.e., low number of generated Frenkel pairs. It is observed that for both interstitials and vacancies the C/Si ratios remain constant for all investigated E PKA Si . Also, for both interstitials and vacancies, more C interstitials/vacancies are produced than Si interstitials/vacancies. Ratios of V C /V Si =1.1 and I C /I Si =3.4 are observed.
The C/Si ratio derived from TRIM and NRT model is in excellent agreement with the MD results for vacancies. It is noted that the C/Si ratio for vacancies is equal to Si-TDE /C-TDE. Thus, both TRIM and NRT (or Si-TDE /C-TDE) provide an accurate estimate for V C /V Si , and based on our MD results, I C /I Si is about three times the V C /V Si ratio.
For comparison, the results obtained in the study by Cowen et al. 32 are added to Fig. 3. In that work, the ratios I C /I Si = 15.7 and V C /V Si = 4.7 are at large discrepancy to the values found in the present work. Also, V C /V Si is not  , GP code 40  46  18  45  14  22  38  21  16  38  19   ERH potential  34  19  58  38  31  31  38  58  45  36   DEV potential  36  31  71  38  113  28  39  71  64  40   Table 2. Threshold displacement energy (TDE) as obtained by experiment, density functional theory (DFT) and potentials. www.nature.com/scientificreports www.nature.com/scientificreports/ equal to Si-TDE/C-TDE. However, it is striking that consistent with our results, the I C /I Si ratio is also about three times larger than V C /V Si .
Qualitatively, for both ERH and DEV, I C /I Si > 1 and V C /V Si > 1 is consistent with the difference in formation energy shown in Table 1. The formation energy for the most stable I C and V C configurations is lower than the one for the most stable I Si and V Si configurations, respectively. However, a quantitative correlation could not be found. The quantitative differences are likely the result of non-trivial differences in dynamic behavior between both potentials.
In Fig. 4, the ratio of Frenkel pairs/anti-sites and C Si /Si C anti-sites for all investigated E PKA Si is plotted. These quantities can only be estimated from MD and appear to be constant with E PKA Si . For any E PKA Si , the number of anti-sites is ~2.5 times smaller than the total number of Frenkel pairs, with C Si /Si C = 0.25. For comparison, the same ratios taken from the study of Cowen et al. 47 are added to Fig. 4. While the Frenkel pairs/anti-sites ratio is in good agreement with our results, the C Si /Si C ratio is very different. While in our case the number of Si C anti-sites outnumber the C Si anti-sites, the opposite is true for the results by Cowen et al.
The large discrepancy between the results of both potentials is reflected in their difference in formation energy of C Si and Si C . While for ERH the formation energy of C Si is similar to the one of Si C , for DEV the formation energy of C Si is an order of magnitude lower than the one for Si C . A priori, the values predicted by ERH (used in this work) are more consistent with DFT than the values predicted by DEV (the work by Cowen et al.).
In summary, the number of produced Frenkel pairs can be described by the NRT model and TRIM simulations. Thereby accounting for electronic stopping is essential. The number of anti-sites is about a factor 2.5 lower than the total number of Frenkel pairs and the C/Si ratio for interstitials is about 3 times larger than for vacancies. These results are consistent with the work of Cowen et al. 32 and therefore independent of potential. The following ratios were observed: V C /V Si = 1.1 and C Si /Si C = 0.25, which is at large discrepancy with the work of Cowen et al. 32 .  www.nature.com/scientificreports www.nature.com/scientificreports/ The reason for these discrepancies are at least qualitatively reflected by the differences in formation energy of the respective point-defects by the different potentials used in both works.
The occurrence frequency, f N ( ), of a defect cluster of size N generated in a cascade for given E PKA Si is summarized in Fig. 5, for both vacancy and interstitial clusters. As expected for ceramics, defect cluster sizes are limited and decrease progressively with decreasing E PKA Si . It is found that vacancy clusters are about twice the size of interstitial clusters.
, the largest observed vacancy cluster contains 12 vacancies, while the largest interstitial cluster contains 5 interstitials. In the next section, an analytical model to describe the defect cluster size distribution is developed.

Analytical Model
Based on the cluster occurrence frequency, the occurrence probability, P N ( ), of a point defect cluster of size N is defined as, This is the probability of occurrence of a cluster of size N amongst all observed clusters. The P N ( ) for vacancy and interstitial clusters obtained from our cascade simulations are plotted in Fig. 6a, for all investigated E PKA Si . Interestingly, for both vacancies and interstitials, all data scales to a single curve, independent of E PKA Si . This indicates that the process of cluster creation occurs ballisticly, i.e. with a certain probability, p, without correlations. As such, a higher E PKA Si leads to the occurrence of larger clusters, only because the number of trials is increased, i.e.  www.nature.com/scientificreports www.nature.com/scientificreports/ the Frenkel pair production. The clustering probability, on the other hand, remains constant. Within this scope, P N ( ) can be described by the geometric distribution. The geometric distribution describes the probability that the first failure (the point-defect does not cluster) requires N Bernoulli trials (thus a cluster of size N is obtained) for a given probability p for success (clustering probability) and is given as, The best fit of Eq. 6 to the data for vacancies and interstitials, is obtained for p = 0.45 and = . p 0 14 for vacancy and interstitial clusters, respectively. These curves are added to Fig. 6a, where the data clearly follow a geometric distribution.
To illustrate the generality of this theory, it is also applied to the results obtained by Liu et al. (the data is generated from Fig. 3 of the respective paper) 31 and added to Fig. 6b. As with our data, P N ( ) for different PKA energy scale to a single curve for vacancy and interstitial clusters, respectively. Clearly, these data also satisfy the geometric distribution, even though the authors fitted a power law to their data. It is observed that the independent cluster probability, p, although slightly increased, is in the same range as found for our data. Small differences are to be expected due to differences in simulation set-up, the most important ones being: different interatomic potential and lack of electronic stopping power.
Thus, for 3C-SiC (and likely ceramics) the cluster occurrence probability is described by the geometric distribution, rather than a power law, that is adequate to describe defect clusters in tungsten 49 . The presented distribution function allows to generate cascade libraries beyond the scope of MD simulations by accounting for rare events that are inaccessible on the MD length scale.
It is emphasized that given the generality of the theory, the applicability of the geometric size distribution can probably extend beyond the case of SiC, possibly extending to a wide class of ceramics.

conclusions
Primary damage in 3C-SiC was characterized using classical molecular dynamics (MD), binary collision approximation (BCA) and NRT model for a PKA energy range 1-100 keV. As such, BCA and NRT are benchmarked against MD. Particular care was taken to account for electronic stopping and the use of a threshold displacement energy (TDE) consistent with density functional theory (DFT) and experiment.
It was found that the Frenkel pair production described by BCA and NRT are equivalent with an MD description. The effect of TDE, within the range provided by DFT, affects the total Frenkel pair production by less than a factor two.
The ratio of number of Frenkel pairs over number of anti-sites is about 2.5, C Si /Si C = 0.25, V C /V Si = 1.1 and I C / I Si = 3.4. It is emphasized that these ratios are constant over the range of investigated PKA energy.
The defect cluster size occurrence probability was found to follow the geometric distribution rather than the power law proposed by Sand et al. 49 to describe defect clusters in tungsten. The geometric distribution is characterized by the number of trials N, before first failure, provided each trial is a Bernoulli trial. This view is consistent with ballistic effects that dominate cascade formation in ceramics and 3C-SiC in particular. It is thus important to emphasize that the applicability of the geometric size distribution can probably extend beyond the case of SiC, possibly extending to a wide class of ceramics.
Thus, the present study allows to fully characterize defect production in 3C-SiC using simple models and regressions.
To conclude, a call is launched for experimental validation by in-situ transmission electron microscopy (TEM) under ion irradiation at cryogenic temperatures. Cryogenic temperature may avoid thermally driven recombination and only allow for ballistic effects to occur. If not quantitatively, the form of the distribution function can be verified.