Cs diffusion in SiC high-energy grain boundaries

One of the fuel types proposed for the next generation high-temperature nuclear gas reactors is Tristructural-Isotropic (TRISO)-coated fuel particles.¹ Recent TRISO fuels are designed to have a fuel kernel, typically UO₂ but can also be UC or UCO, surrounded by four successive layers consisting of C buffer, inner pyrocarbon, silicon carbide (SiC), and outer pyrocarbon.¹ The SiC layer mainly consists of the cubic 3C-SiC polytype, with a thickness of ∼35 μm. The primary functional of this SiC layer is to provide a structural stability and to be a barrier for fission products (FPs). While this layer in the modern fuels effectively retains most of the FPs under both operating and accident conditions,² some undesirable metallic FPs (e.g., Ag, Cs, and Sr) are historically reported to be released, particularly under accident conditions. Among these FPs, the release of the Cs isotopes has been long studied due to the large quantity of its production and the relatively long radioactive half-life (¹³⁷Cs half-life is 30.2 years). Such release of Cs has significant safety and design implications.

To address issues associated with the Cs release from TRISO particles, it is important to first understand the transport mechanisms through the SiC that is the primary barrier in fuel particles. In the past few decades, studies have been conducted with experimental or computational approaches to elucidate the diffusional rate of Cs and its release mechanism. The previously reported diffusion prefactors and activation energies for Cs are compared in Table I, and these studies can be put into three categories, which will be briefly described in the following.

The first category of diffusion studies is based on integral release (I.R.) measurements. In these measurements, the effective diffusion coefficients are determined by fitting a simple diffusion model to the overall integrated FP release data obtained from a batch of TRISO fuel particles.^3,4 While a few in-reactor data are available,^5,6 most of these experiments measured the release of Cs from TRISO particles that were pre-irradiated and then post-annealed at high temperatures. The effective diffusion coefficient through SiC is estimated by integrating with models for the Fickian diffusion through the other coatings in TRISO and by adjusting the diffusion coefficients in SiC to match the observed FP release.⁷ An interesting observation from combined I.R. measurements in Table I is that the activation energies are in the order of 1–2 eV for low-temperature annealing (i.e., less than about 1400 °C), while they are in the order of 4–5 eV for annealing at higher temperatures (above 1500 °C). The difference in activation energies suggests that there are two mechanisms which become dominate in different temperature regimes.

The second category of diffusion studies is based on ion-implantation (I.I.) measurements, where the diffusion of Cs in single crystalline (sc-SiC) 6H-SiC^8,9 and polycrystalline (pc-SiC) chemical vapor deposition (CVD-) SiC^8,9 were investigated using an ion implantation, followed by heating tests. These surrogate attempts can help understanding the release mechanism of FPs; however, this method is known to create several complications. For instance, in I.I. studies, the diffusion of Ag may be affected by the introduction of trapping sites from the implantation damage¹⁰ and the change in diffusion mechanisms due to the higher concentrations of FPs than the concentrations in TRISO under operation.¹¹ Recently, a novel design for the I.I. experiment was proposed¹² to avoid these complications, by introducing FPs into pyrocarbon adjacent to SiC without causing implantation damage in the SiC and at concentrations that would be relevant to TRISO fuel.

The third category of diffusion studies is based on computational studies, which have been carried out primarily to develop an atomistic level understanding of Cs diffusion. An ab initio study by Shrader et al.¹³ showed that the Cs thermal diffusion in the bulk is slow (∼10⁻³²–10⁻²⁵ m²s⁻¹ at 1200–1600 °C), and these low bulk diffusion rates cannot account for the measurements in I.R. However, Shrader et al.¹³ found that the activation energy for Cs bulk diffusion agrees with the I.R. measurements of temperature above ∼1500 °C, which suggests that Cs release occurs via bulk diffusion when irradiated. Recently, Rabone et al.¹⁴ investigated that Cs diffusion in ∑5 coincides with the site lattice GB (∑5-GB) using an ab initio molecular dynamics at 1227 °C. This study showed that the Cs diffusion in ∑5-GB is 3–4 orders of magnitude faster than the bulk diffusion.

The previous I.R. studies suggest that there are two diffusion mechanisms involved in Cs diffusion.¹⁵ One hypothesizes to explain these trends is that grain boundary (GB) diffusion dominates at low temperature because of having an expected lower activation energy, while bulk diffusion dominates at a higher temperature where the higher temperature makes the migration energy less critical, and the larger accessible region for transportation allows the bulk diffusion to dominate. Another hypothesis, proposed by Malherbe,¹⁵ is that there is de-trapping of Cs atoms from irradiation-induced defects in the bulk at high temperatures, and these de-trapped Cs atoms then diffuse via a volume diffusion. However, the large variations within the I.R. measurements and the difference in their experimental condition (e.g., thermal exposure time) make it hard to determine which transport mechanism dominates and how quickly Cs releases. Furthermore, it remains unclear if there is any radiation effect on the release of Cs, and if yes, a quantitative assessment is necessary to determine how significant the effect is.

As mentioned above, the GB diffusion is postulated as a dominant mechanism at low temperatures, but a critical missing part of the present understanding is the diffusion mechanisms and diffusion rates in GBs. Among GBs, the high-energy GBs (HEGBs) are expected to play a significant role in the FP transportation in pc-SiC, for the following reasons: The HEGBs are often highly disordered structures and represent >40% of GBs in TRISO prototype materials with grain sizes of 0.5–1.5 μm.^16,17 The high fraction of HEGBs allows them to provide a percolating path for FP transportation, and HEGBs are one of the two GB types that are present in a high enough concentration to form a percolating path (the other GB type is ∑3-GB).¹⁸ Also, it is often found that disordered (amorphous) materials, which are structurally similar to HEGBs,^19,20 provide a faster transportation pathway for extrinsic defects compared to crystalline materials.^21–23 Furthermore, our recent study has suggested that Ag diffusion via HEGBs is the dominant path for diffusion in unirradiated SiC.²⁴ Therefore, we expect HEGBs to be the most likely GB type to dominate Cs transport and the HEGB Cs diffusion coefficient to be faster than those in other GBs or bulk. In this work, we used an ab-initio based stochastic modeling approach to predict the Cs diffusion in HEGBs. Based on the diffusion coefficient from the model, we will determine under what conditions (e.g., temperature), if any, HEGB D is a dominant pathway in SiC and responsible for the release of Cs.

The approach taken in this work follows that used for a similar study on Ag in SiC from Ko et al.²⁴ We briefly repeat the essential details here, with relevant changes for treating Cs in place of Ag.

Modeling of the HEBGs, often called disordered GBs, has been a great challenge in covalent materials.²⁵ In this work, we follow the approach of Ref. 24 and model the HEGB structure as an amorphous SiC (a-SiC) because (i) the local environments in HEGBs of covalent materials are known to be similar to amorphous phases^19,20 and (ii) a bulk a-SiC is computationally more tractable than modeling a full GB structure. A simulation cell is used that is approximately cubic and contains 128 atoms (64 silicon, and 64 carbon). The cell was trimmed from a bulk a-SiC, which was prepared with the melt-quench method²⁶ using a Tersoff potential.²⁷ This supercell was then fully relaxed (constant zero pressure and zero K temperature) using the density functional theory with the periodic boundary conditions. The relaxed supercell had cell vectors of [11.50, −0.12, 0.04; −0.12, 11.45, −0.10; 0.04, −0.08, 11.51] (non-cubic). The supercell used in this study is shown in Fig. 5 of the  Appendix. The density of SiC in this supercell, after fully relaxed with the density functional theory (DFT), was 2.81 g/cm³, which was the comparable value with the simulated a-SiC densities of 3.057 and 2.896 g/cm³ (Ref. 26) (by Tersoff potential) and was lower than the sc-SiC density of 3.18 g/cm³ (calculated with DFT), as expected.

Ab initio density functional theory (DFT) calculations were performed using the Vienna Ab-Initio Simulation Package (VASP).^28–31 The VASP calculations were performed using the projector-augmented plane-wave (PAW) method^32,33 and the generalized gradient approximation (GGA) exchange-correlation potentials, parameterized by Perdew, Burke, and Ernzerhof (PBE).³⁴ The PAW file electronic configurations used: 3s2 3p2 for Si, 2s2 sp2 for C, and 5s2 5p6 6s1 for Cs. These pseudopotentials have been widely used in SiC studies^35–37 and Cs studies;¹³ therefore, no benchmarks are presented in this study. The energy cut-off was set as 450 eV, and a single Γ-point k-point mesh was used to sample reciprocal space. The plane wave energy cutoff and the k-point mesh density were separately tested for eight cases and were also converged to give a total energy within 180 meV (with a minimum of 32 and a maximum of 175 meV) per Cs defect compared to 3 × 3 × 3 k-point mesh, where this larger mesh is expected to be well-converged.

To obtain diffusion pathways and the corresponding energy barriers, the climbing image nudged elastic band (CI-NEB) method was employed.³⁸ For CI-NEB calculations, three images were used and images were linearly interpolated between minima. Following the approach of Ref. 24, in the cases of images where the Cs atom is too close to Si or C atoms, their positions were adjusted to be separated by 75% of their covalent bond length (2.7 Å for Cs-Si and 2.4 Å for Cs-C³⁹). The k-point mesh test on eight barriers showed that the error in the energy barrier with respect to k-points for a Γ-point vs. 2  × 2 × 2 k-point mesh is within an acceptable error range of 220 meV/Cs (with a minimum of 53 and a maximum of 216 meV), which corresponds to about a factor of 5× error at 1500 K.

For the defect formation energy (⁠ $Δ E f$⁠) of defects, we used the following expression:⁴⁰ 

where $E d e f$ and $E undef$ are energies of the defected and the undefected cells, and $Δ n I$ is the change in the number of the atomic species $I$ in the defected cell from the number of the same species in the undefected cell. $μ I$ is the chemical potential of atomic species I relative to its reference states, which are taken as the bulk Si and C VASP energies (which are in turn referenced to the appropriate atomic species energies given as defaults in the pseudopotential files) in their groundstate structures (diamond lattice for Si and graphite for C). These values are $E Si$ = −5.431 and $E C$ = −9.186 eV/atom. Throughout this study, the Si-rich condition chemical potentials are used for silicon and carbon (⁠ $μ Si$ = −5.431 and $μ C$ = −9.619 eV) for consistency. It should be noted that the chemical potential for each element of a binary system is not unique^41–43 and may vary within a range between Si-rich and C-rich conditions. Under C-rich conditions, C based defects will increase their formation energy in the positive direction (toward being less stable) by the calculated heat of formation of SiC (0.44 eV). However, this will not affect our results for Cs diffusion via the interstitial mechanism since interstitial is the most stable form of the Cs defect regardless of the chemical potential choice. The ab-initio formation energy of the bulk solid phase Cs metal (⁠ $μ Cs$ = −0.856 eV) was set to be the chemical potential for Cs. We found that the interstitial diffusion is expected to be the dominant mechanism in the HEGB, and it will be the focus in this study (further discussion in Sec. III A and Table II). The result was also found to be independent of the chemical potentials chosen for Si and C.

It should be noted that no charged supercells or explicitly charged defects were considered in the present study. We believe that this approximation is reasonable as the neutral state for Cs interstitials is the stable charge state for Cs in n-type crystalline SiC, as predicted by Shrader et al.¹³ Also, the effects of the association of Cs defects with intrinsic defects and the interaction between Cs defects had not been considered. The results should be interpreted within these limitations and are an indication of what the behavior would be if only non-interacting defects are present in the material.

To model Cs interstitial (⁠ $Cs I$⁠) diffusion in an a-SiC, the energy landscape of Cs interstitials was first investigated. In an amorphous material, locating extrinsic interstitial defect sites has many complications due to the lack of long-range order. In this paper, we used a simple gridding method, in which we grid the entire supercell with a fine uniform grid and relax a Cs interstitial at every grid point in the cell, to obtain a complete list of possible $Cs I$ sites. This method for finding sites is more computationally expensive but more comprehensive when compared to approaches often taken in previous studies, such as geometry guided guesses^44,45 and free-volume based search.⁴⁶ In this method, the entire supercell is gridded by a 1 Å grid along each axis, and Cs interstitials are placed at every grid points and then relaxed to search all existing sites. The details of the gridding method, including the choice of grid size and an assessment of its ability to identify all interstitial sites, can be found in our previous work on Ag diffusion in HEGB.²⁴ 

The modeling of diffusion in the amorphous systems has been a challenging problem in many applications. One of the most pragmatic and simple methods is molecular dynamics with interatomic potentials; however, potentials for such complex systems (e.g., FPs in SiC) are not developed. There are other plausible methods such as ab-initio molecular dynamics⁴⁷ or potential energy surface sampling approach (e.g., kinetic activation relaxation^48,49 and autonomous basin climbing⁵⁰), which can be used, but the computational demands of such calculations are often prohibitive with full ab-initio methods. Therefore, in this work, we use an approach based on the kinetic Monte Carlo (kMC) model. This section describes the methods used in the KMC model, and the results of using these approaches are given in Sec. III B.

An essential input for the KMC model is the energy landscape for hops. In our previous study of Ag diffusion in HEGB,²⁴ we modeled an effective medium (virtual lattice) that mapped the effective energy landscape of Ag interstitials. However, to avoid additional complexity and uncertainties introduced by the effective medium approach, it was not used in this study. Instead, the diffusion was directly modeled by applying periodic boundary conditions to the simulation cell where the Cs interstitial sites were identified by the DFT calculations. A similar approach in a finite cell with periodic boundary conditions successfully modeled the Li diffusion in amorphous (a-) materials, e.g., a-Si,⁵¹ a-Al₂O₃,⁵² and a-AlF₃,⁵² where they generally showed good agreements with experimental diffusion rates when the comparison was possible.

Overall, the 55 interstitial sites have 3.4 neighbors within 4.0 Å, which results in a large number of hopping barriers to determine. To reduce the number of calculations, barriers were only calculated until the diffusion coefficient (D) showed convergence with respect to the barrier sampling. Here, we describe the approach used in detail. The $Cs I$ sites were identified first, and then, migration barriers between these sites were calculated. To reduce the number of hops that needed to be considered, it was assumed that the diffusion was likely dominated by the more stable $Cs I$ sites. With this assumption, the most relevant migration barriers were determined as follows: A set of sites was created by taking the n most stable sites. Then, all migration barriers between all the (⁠ $n!/2$⁠) pairs were determined with DFT for any pair of sites within a cutoff distance of 4.0 Å. This was the longest distance over which no intermediate state was found from preliminary test calculations of randomly sampled 25 barriers. This procedure was done first for n = 1 (i.e., the most stable site), then n = 2 (i.e., the 2 most stable sites), and so on. At each n, the migration barriers determined between those n sites were used in the kinetic Monte Carlo simulation to evaluate the Cs diffusion coefficients on the network of those n Cs interstitial sites. This process was done for increasing n until the diffusion coefficient was converged, i.e., stopped changing as n increased.

As a further test of convergence for a given n, migration barriers for all other $Cs I$ sites were estimated and are also included in the model, to see if the estimate for D changed. The approach to estimating the values not calculated with DFT made use of the kinetically resolved activation (KRA) formalism.⁵³ First, we define some notation. We refer to the network of n interstitial sites and their calculated barriers as subsystem S_n. We refer to S_n with additional barriers estimated by the KRA formalism (E_KRA) for all remaining interstitial sites as S_n⁺. We refer to the diffusion coefficient calculated on the subsystem S_n or S_n⁺ by D(S_n) or D(S_n⁺), respectively. In a simple KRA approach, unknown migration energies between sites X and Y are estimated from the energies at sites X and Y by setting the transition state energy between them, E_TS (X,Y), to

where $E X$ and $E Y$ are the site energies of $Xth$ and $Yth$ stable $Cs I$ and E_KRA is approximated as a constant. For E_KRA in the system S_n⁺, we have used the E_KRA value of zero. This approximation sets $E T S (X,Y)=max( E X , E Y )$⁠, which leads to barriers of $E X − E Y$ for hops that increase energy and zero for hops that decrease energy. This approach provides an estimated upper bound on the D one might get using DFT values for all the interstitial sites. When D(S_n) stops changing with n for three consecutive n values and D(S_n) = D(S_n⁺) for all these n values, we conclude that the system size n is sufficiently large to include all essential interstitial sites along the diffusion path, and thus, the D is converged. One problem with this approach is that for a few sites, the very low barriers (which occur when $E X ≈ E Y ,$ as E_KRA = 0, and then leads to a near zero barrier) caused the Cs to become trapped hopping between just two sites during the KMC calculations. For these cases that led to trapping (there were four for n = 23), we calculated and used the DFT value of the barrier, which was always much higher than the lower-bound estimate, and therefore removed the trapping.

The kinetic Monte Carlo (kMC) method was employed with the Bortz-Kalos-Liebowitz⁵⁴ algorithm. The hopping rates (⁠ $Γ$⁠) for Cs atoms were given by the transition state theory as $Γ=ω· e ( − E m / k B T )$⁠, where ω is the attempt frequency and E_m is the migration barrier of the hop. The attempt frequency (ω) can be determined as $ω= ∏ 3 N v equil / ∏ 3 N − 1 v saddle$⁠, where $v equil$ and $v saddle$ are the vibrational frequencies at the equilibrium and saddle point, and N is the number of atoms in an adequately large region around the diffusion path.⁵⁵ We took N to be all atoms within a 4 Å cutoff of the Cs in the transition state of the hop being considered, which was typically 12–16 atoms. Tests for longer cutoff distances and larger N result in a change in ω of less than 5%. From ab-initio calculations for 5 test cases, we found that the attempt frequencies typically range from 10¹² to 10¹³ s⁻¹, with an average of 3.8 × 10¹² s⁻¹. This average value was then used for all hops.

The diffusion coefficients were determined by the Einstein relation, $D=⟨ r 2 t ⟩/2dt$⁠, from the calculated mean square displacement as a function of time,⁵⁶ where d is the dimensionality of the system (here d = 3), t is the time, and  $⟨ r 2 t ⟩$ is the mean square displacement of Cs as a function of time. We evaluated the average using the multiple time origin method.⁵⁷ The details of the kMC can be found elsewhere.²⁴ Simulation at each n was typically performed for 5 × 10⁹ kMC steps to obtain a well-converged diffusion coefficient.

To understand the transportation behavior of Cs in a-SiC, it is important to determine which of the defect types are stable and will contribute to the diffusion. The formation energies (E_f) of the point defects (vacancies, Cs substitutionals on Si and C lattices, and Cs interstitials) were investigated. Each of these defects was sampled over a wide range of local environments of a-SiC. In Table II, the average and standard deviation of E_f values over many defect sites in the a-SiC are shown to represent the spread in values. The formation energies of these defects in the sc-SiC are also presented for a comparison. The sample set includes 15 substitutionals and 55 interstitials, where this latter set is quite comprehensive as they were determined to be the dominant mechanism and studied in detail (see Sec. III A). All the defect types in a-SiC, on average, had lower E_f than those in the sc-SiC. Particularly, Cs interstitials in the a-SiC showed a dramatic difference compared to E_f in sc-SiC. For all identified sites, the mean E_f of Cs interstitials in a-SiC was lower by >13 eV, compared to the most stable Cs interstitial (23.46 eV) in sc-SiC (6–fold coordination for Cs, with three carbon and three silicon atoms). Some of the Cs interstitial formation energies were near zero or negative, i.e., the Cs interstitials are more stable in the a-SiC than in a bulk metallic Cs atom. This result suggests that Cs will segregate strongly to disordered HEGBs, consistent with what was found in the previous GB modeling study.¹⁸ 

Under the dilute Cs condition, the diffusion can take place via either a vacancy-mediated substitutional diffusion or some form of interstitial diffusion. The vacancy-mediated diffusion mechanism, however, requires a migration of a vacancy.⁵⁸ The barriers to migrating a vacancy around Cs in a-SiC were found as 2.6–4.5 eV, and the barriers to migrating a vacancy in a undefected a-SiC were 3.0–4.5 eV (consistent with the previous study²⁴) from tests on five barriers for each case. These values suggested that the vacancy mediated Cs diffusion is expected to be relatively slow compared to interstitial migration, which can avoid moving Si and C atoms. Therefore, in this study, we only considered the interstitial diffusion mechanism for Cs in the a-SiC. Our results below predicted a migration barrier of 2.41 eV through the interstitial mechanism. As this value is smaller than the ranges mentioned above, it strongly supports our hypothesis that interstitial diffusion is dominant for Cs in a-SiC. No interstitialcy (or kickout) mechanisms⁵⁹ were considered.

By applying gridding to the a-SiC sample, we investigated Cs interstitials on 12 × 12 × 12 grid points. The interstitial sites that were both energetically and geometrically close to each other were grouped and represented by the lowest energy site in the group. We have referred this grouping as a hierarchical clustering approach, and the details are available elsewhere.²⁴ After hierarchical clustering, a total of 55 Cs interstitial sites were identified. In Fig. 1, the distribution of formation energies is shown. The formation energies of Cs interstitials in a-SiC had an average value of 0.41 eV, and this value is ∼23.1 eV more stable than the most stable Cs interstitial in sc-SiC, as can be seen in Table II. The values ranged from −2.68 to 8.15 eV, with a standard deviation of 3.40 eV. The Cs, which is an impurity with a relatively larger size than those studied in the similar amorphous systems,^{46,47,51,60–63} turned out to be highly stable as interstitials in the a-SiC. The high stability of Cs defects vs. bulk is possibly due to the large relaxations available to the a-SiC system, analogous to the observation of Ag defects in the HEGB,²⁴ although the detailed mechanism of stabilizing Cs interstitials in a-SiC was not explicitly investigated here.

After identifying interstitial sites, the migration barriers were calculated for the Cs. Initially, we randomly sampled 25 E_m and associated kinetically resolved activation barriers (E_KRA), which can represent migration barriers independent of the initial and final state energies (details are presented in our previous study, Ref. 24). As it can be seen from Fig. 2, more than 40% of E_KRA were found to be less than 1 eV, suggesting that a fast-interstitial diffusion may be possible. To obtain the complete set of barriers and hopping paths, we have sampled migration barriers by increasing the sampling size (n) starting from sites with the lowest energy, as explained in Sec. II B and with the results given in Sec. III B.

Ideally, all the migration between 55 Cs interstitial sites can be sampled to obtain the effective diffusion coefficient in HEGB (D_HEGB), but this would be computationally intractable. As discussed in Sec. II B, we have assumed that Cs visits only the more stable sites during its diffusion. Therefore, the kMC simulations were performed at a set of sites of size n, where n was increased to include less stable sites until D showed convergence. Furthermore, following the approaches discussed in Sec. II B, we have calculated D for the system of sites and barriers associated with the n points (S_n) and a system of all points where additional barriers are calculated by the KRA method (S_n⁺). The resulting D(S_n) and D(S_n⁺) from two methods are shown in Fig. 3. For n ≤ 8, there was no diffusion across the boundary of the simulation cell for both S_n (no diffusion path with small n) and S_n⁺ (Cs is trapped in few n). When a faster diffusion path is added (S_n=9, S_n=12, S_n=18, S_n=9⁺, and S_n=18⁺), a steep increase in D is observed, as shown in Fig. 3. For 9 ≤ n ≤ 20, we have D(S_n⁺) > D(S_n). This is because the sites that were not included in the S_n system, i.e., the (55–n) sites, were playing a role in Cs diffusion as they generally offer hops with low E_ms (they are estimated with E_KRA = 0 eV, which provides very low barriers). For n ≥ 21, both the convergence of D with n, and the equality of D(S_n) = D(S_n⁺) support that n = 23 is large enough to model D for the full set of 55 Cs interstitial sites. The diffusion coefficients reported in the following were calculated with S_n=23.

The diffusion coefficients using n = 23 were simulated for 1200–1600 °C and are summarized in Fig. 4. In this temperature range, the D_HEGB values ranged from 10⁻¹⁹ to 10⁻¹⁷ m²s⁻¹, which is 4–6 orders of magnitude higher than the bulk diffusion coefficient (D_bulk).¹³ The HEGB Cs diffusion exhibited an Arrhenius type diffusion as shown in Eq. (3)

where D₀ is the diffusion prefactor, E_A is the effective activation energy for the interstitial diffusion, and k_b is the Boltzmann constant. This Arrhenius form fitted to the KMC simulations yielded D₀ and E_A values of (3.75  $±$  $3.15)$× 10⁻¹¹ m²s⁻¹ and 2.41  $±$ 0.26 eV, respectively. It is interesting to note that for a simple ideal interstitial diffusion model, we expect the prefactor, $D 0 (ideal)$⁠, to be $D 0 ideal =g l 2 ω$⁠, where $g$ is a geometric factor of order one, $l$ is the hop length, and ω is the attempt frequency. Using the attempt frequency given above for our model and a hop length of 3 Å (similar to our average hop distance) yields an estimate of $D 0 ideal =$(3.4) × 10⁻⁷ m² s⁻¹, about 10⁴ times larger than our fitted value. A similar suppression of the D₀ value vs. a simple model was found for Ag in SiC in Ref. 24. While the D₀ value can be strongly affected by small deviations from Arrhenius behavior when fitting in only a small range at high temperatures, this result also suggests that the Cs diffusion exhibits significant reduction in diffusion relative to a simple model, presumably due to the disordered energy landscape through which it must move.

It should be noted that our Cs HEGB diffusion is modeled by a limited size simulation cell, and this cell is likely not fully capable of representing the a-SiC nature. For this reason, D_HEGB is expected to have some uncertainty originating from the sampling. This uncertainty had been estimated in the previous work on the Ag diffusion in a-SiC.²⁴ By taking the uncertainty estimated in Ref. 24, we expect D_HEGB to have an uncertainty of approximately a factor of 10 (up to 100 at the most) due to a limited structural sampling.

The HEGBs are expected to be the dominant diffusion path for Cs not only because of the diffusion rate but also because the highest fraction of GBs is HEGBs and because Cs clearly has strong segregation tendencies to the HEGB. Here, we discuss each of these properties, their relation to net diffusion in pc-SiC, and the comparison to other Cs release studies.

As mentioned in Introduction, for HEGBs to be a dominant Cs release path, they must enable a connected (percolating) pathway through the polycrystalline material. In pc-SiC, ∑3-GBs and random GBs (HEGBs) are the most abundant type of GBs^16,64 and the HEGBs constitute a majority, more than 40%, of GBs in CVD-SiC.¹⁷ Among the GB types, only HEGBs and ∑3-GBs, with their high fraction of all GBs, are expected to provide percolation paths for FP diffusion in SiC.¹⁸ 

If HEGBs form a percolating network, the diffusion rates in HEGBs must be compared to other possible paths in pc-SiC to determine which is the fastest. The present work showed that HEGBs can provide much faster diffusion paths compared to the bulk,¹³ suggesting that the GB diffusion dominates over bulk diffusion in most cases. The study on Cs diffusion in a low angle ∑5-GB by Rabone et al.¹⁴ predicted a range of D_∑5-GB to be 2–4 orders lower than D_HEGB as shown in Fig. 4. Unfortunately, studies on the GB diffusion are very limited in a number, and it is not possible to compare diffusion rates in a range of different GB types. In particular, the diffusion data in ∑3-GBs, which is also likely to form a percolating path, need to be determined to conclude that the HEGB is indeed the fastest and dominant pathway.

Even though there are no D_∑3-GB data available, we hypothesize that HEGBs are likely to provide the fastest path among the GBs. This hypothesis is reasonable because the diffusion barriers of Cs in ∑3-GB are expected to be greater than those in HEGBs. Assuming that the diffusion prefactors in the Arrhenius equation are in a similar order of magnitude for both the ∑3-GB and HEGB, the effective activation barriers contribute dominantly to their diffusion coefficient. The effective activation barrier in HEGBs is expected to be lowest, as the diffusion barriers are, in general, inversely proportional to the free volume of the system. This idea is supported by the Cs migration barrier of 4.65 eV in the ∑5-GB,¹⁴ whereas the effective barrier in the HEGBs was 2.41 eV. Also, in the similar study of Ag diffusion in unirradiated SiC, the D_∑3-GB¹⁸ was found to similar or lower than D_HEGB.²⁴ Analogous trends are expected for Cs, perhaps favoring HEGBs even more, as Cs is a much larger atom than Ag, and therefore, the movement of Cs will be sluggish in the more confined geometry of the ∑3-GB.

Assuming that HEGBs are the dominant GB diffusion pathways, it is important to consider the effect of microstructural features on the effective diffusion we expect in pc-SiC. The net diffusion coefficient in pc-SiC (⁠ $D n e t$⁠) allows the comparison to the diffusion coefficients measured from both I.I. and I.R. in pc-SiC. Diffusion coefficients obtained from I.R. measurements are effective values that combine contributions from all transport mechanisms available to Cs. The D_net value is given by the modified Hart equation:⁶⁵ 

where s is the segregation factor, d is the grain size, and δ is the GB width. For D_GB, we approximate it to be D_HEGB. The segregation factor (s) can be calculated as

where C_GB and C_Bulk are the equilibrium Cs concentrations at the GB and bulk. We assume that C_GB is approximately equal to C_HEGB. The equilibrium concentrations can be calculated from formation energies of Cs defects in the GB and bulk, as shown in Eq. (6)

where i is the type of defect, $ρ i$ is the site density of defects of type i per unit volume, and $E f i$ is the formation energy of the defect of type i. With this simple non-interacting model for the Cs solubility, we predict the solubility limit in the HEGB (⁠ $C HEGB Cs$⁠) to be 2.25 × 10²⁶ m⁻³ at 1200 °C. The solubility limits in the sc-SiC can be calculated using values from Ref. 13, as $C bulk Cs$ = 1.9 × 10⁻⁹m⁻³. From Eq. (5), the segregation factor is found as s ∼ 10³⁵ at 1200 °C and s ∼ 10²⁷ at 1600 °C. Although there is no experimental evidence showing a Cs segregation to GB (while Ag atoms are detected in GBs and triple junctions⁶⁴), these large values of s clearly suggest a strong segregation of Cs to GBs (HEGBs).

The D_net value can be estimated from Eq. (4) by using s values calculated from above and using experimentally measured values of d and δ, which are taken as 1 μm and 0.5 nm, respectively. Due to high segregation factors at all relevant temperatures, D_net becomes approximately equal to D_HEGB (the difference is less than 0.001% of D_HEGB, at most). In other words, the GB diffusion governs the total diffusion. While it is hard to predict the real diffusion process in the radiated material, a simple diffusion kinetic model involving both the bulk and GB (types A and B, Harrison's model⁶⁶) is unlikely based on the analysis presented here.

The contribution of D_HEGB to the net diffusion in pc-SiC can be determined by comparing this value to previously reported I.I. and I.R. measurements. For the ion implantation study, Audren et al.⁶⁷ and Friedland, van der Berg, Hlatshwayo, Kuhudzai, Malherbe, Wendler, and Wesch⁸ have investigated Cs diffusion in the ion implanted SiC. In both studies, no Cs diffusion was observed below the certain temperature (around 1050 and 1200 °C, respectively), and Cs started to become mobile above those temperatures. Friedland et al.,⁸ postulated that Cs impurities are trapped by the implantation damage, and thus, no movement of Cs is observed. The implantation damage is a disadvantage in the ion implantation studies, while it is not a concern in the TRISO case. To resolve the issue, Dwaraknath and Was¹² have developed a novel design for I.I., where the FPs are implanted in a thin pyrocarbon layer which is in between two SiC layers. In this way, no direct ion implantation damage is done on the SiC layer. Dwaraknath and Was⁶⁸ showed that Cs diffusion indeed takes place in a temperature range of 900–1300 °C, if there is no implantation damage. The calculated D_GB from Ref. 68 showed a good agreement with the D_HEGB calculated in this work, giving measured values about an order of magnitude lower over a wide temperature range, as it can be seen in Fig. 4. It is worth noting that a factor of 10 is within the uncertainty we expected due to the sampling from a finite cell. Furthermore, assuming an Arrhenius behavior, the activation energy from Dwaraknath and Was⁶⁸ is 2.6  $±$ 0.6 eV, which compares very well to the values predicted in this work of 2.41  $±$ 0.26. The agreement strongly supports that the HEGB is the dominant diffusion path for Cs in non-irradiated SiC. However, it is worth noting that the bulk diffusion coefficient determined by Dwaraknath and Was is orders of magnitude higher than that predicted from ab-initio methods by Shrader et al.,¹³ as well as having a much lower activation energy (1  $±$ 0.1 eV from Dwaraknath and Was vs. 5.14 eV from Shrader et al.). The origin of this discrepancy is still unclear and suggests a need for further study of possible bulk mechanisms for Cs transport.

The diffusion coefficients extracted from I.R. measurements (D_IR) are greater than the D_HEGB values in all temperature ranges, as it can be seen in Fig. 4. Particularly, in the relevant temperature ranges of 1200–1600 °C, there are approximately 1–3 orders of magnitude difference between the predicted D_HEGB and the D_IR. Compared to I.R. measurements, the ion implantation measurements by Dwaraknath and Was⁶⁸ also exhibit a low value, quite consistent with our calculations. The I.I. experiments and our calculations are under conditions quite different from I.R., but perhaps the most obvious difference is that I.R. measurements all involved irradiation. Therefore, we speculate that the predicted D_HEGB is representative of the effective diffusivity in unirradiated pc-SiC but that irradiation accelerates the Cs diffusion. In general, radiation enhanced diffusion (RED) has been widely observed in metallic nuclear materials and semiconductors.^69–71 More specifically, a recent study on Cs diffusion in SiC by Dwaraknath and Was⁷² observed RED of Cs in Si⁺⁺ irradiated pc-SiC, and the enhancement was about 2 orders of magnitude in Cs diffusion, even when the diffusion was measured after the irradiation had stopped. At this point, it is not clear how the irradiation might enhance the Cs transport and whether GB diffusion is the active mechanism, perhaps with RED, with irradiation. There is no detailed determination of the evolution of the HEGB structure in SiC under irradiation of which we are aware. However, given that the grains are stable after irradiation,⁷³ we do not expect the HEGB to completely transform or be removed from the system. Nonetheless, as we conclude that irradiation will change the Cs diffusion, there may be structural changes in HEGBs, which lead to these transport changes. Furthermore, there are other amorphous systems that show structural changes under irradiations.⁷⁴ Overall, the present study supports that HEGB is the dominant pathway for Cs release in unirradiated SiC, and enhancement is expected due to the radiation.

In summary, Cs defects in a-SiC (used to represent a HEGB) were studied with DFT calculations and Cs diffusion in the HEGB was modeled using the kMC method. The interstitial was found as not only the most stable defect type but also the fastest diffusion mechanism for Cs in a-SiC. The low formation energies of Cs defects indicated a strong segregation to GB. The Cs interstitial diffusion coefficients were calculated for 1200–1600 °C, over which range they exhibited Arrhenius type diffusion. The diffusion prefactor and effective barrier values were calculated to be (3.75 $±$ 3.15) × 10⁻¹¹ m² s⁻¹ and 2.41  $±$ 0.26 eV, respectively. We predicted that D_HEGB is 4–6 orders higher than D_bulk in the operating temperature range. A good agreement with D_GB from a recent ion implantation study suggested that the HEGB is the dominant pathway for Cs diffusion in unirradiated 3C-SiC. However, the discrepancy between D_HEGB and D from ion implantation studies with irradiation and integral release measurements suggested that radiation enhanced diffusion may occur for Cs in SiC.

This material is based upon the work supported by the U.S. Department of Energy office of Nuclear Energy's Nuclear Energy University Programs under Grant No. 12–2988. Computations in this work benefitted from the use of the Extreme Science and Engineering Discovery Environment (XSEDE), which was supported by the National Science Foundation Grant No. OCI–1053575.