The energies and stresses associated with the decohesion of *β*-SiC in the presence of mobile Pd and Ag impurities are studied from first principles. Density functional theory calculations are parameterized with a generalized cohesive zone model and are analyzed within a thermodynamic framework that accounts for realistic boundary conditions in the presence of mobile impurities. We find that Pd impurities will embrittle SiC when Pd is in equilibrium with metallic Pd precipitates. Our thermodynamic analysis predicts that Pd embrittles SiC by substantially reducing the maximum stress of decohesion as a result of a phase transition between decohering planes involving an influx of Pd atoms. The methods presented in this work can be applied to study the thermodynamics of decohesion of SiC in other aggressive environments containing oxygen and water, for example, and yield environment dependent cohesive zone models for use in continuum approaches to study crack propagation and fracture.

## I. INTRODUCTION

Silicon carbide is a very strong covalent crystal with a high melting point, making it an ideal structural material for applications requiring high strength at high temperatures. SiC also exhibits very low solubilities for impurities so it can serve as a diffusion barrier for metallic elements. These properties have led to its use in ceramic matrix composites and have made it a strong candidate as a cladding material to contain uranium oxide and its fission products in gas cooled nuclear reactor designs. Silicon carbide's electronic properties also make it an attractive compound for electronic devices subjected to extreme environments and it is used for high temperature hydrogen and hydrocarbon sensors based on Shottky diodes consisting of Pd on SiC.^{1} Single crystal silicon carbide has even found application as a substrate for the growth of single graphene layers.^{2}

While considered a promising diffusion barrier, experimental evidence nevertheless suggests that various metals such as Ag are able to permeate through SiC when applied as a coating material in TRISO (TRIstructural-iSOtropic) particles for nuclear reactors.^{3,4} This implies that transition metals pass through SiC either by means of bulk and/or grain boundary diffusion or by a mechanism that emerges upon degradation of the SiC shell.^{5,6} So far, there have been no experiments that have been able to isolate Fickian diffusion of metal impurities in SiC. A consequence of the very high melting temperature of SiC is that impurity solubilities and the concentrations of thermally activated point defects are exceedingly low at typical operating temperatures.^{7} Diffusion couples involving SiC and Ag or Pd, for example, have not led to measurable concentrations of dissolved Ag or Pd in the bulk of SiC.^{8–10} Ion implantation methods are able to yield high nonequilibrium impurity concentrations in SiC along with high concentrations of defects, but these very quickly cluster and precipitate out.^{3,9–12} Recent experiments suggest that the observed migration of Ag through the SiC layer in TRISO particles may occur as a result of corrosive damage to SiC due to the presence of Pd.^{13–16}

In this paper, we investigate impurity embrittlement of *β*-SiC. We focus, in particular, on the role of Pd and Ag in reducing the resistance to the propagation of intragranular cracks in SiC, as Pd and Ag are fission products of uranium and have been observed to pass through the spherical SiC shells of TRISO particles. The thermodynamic framework that we develop here is applicable to any impurity (e.g., oxygen, water, etc.) that may also embrittle SiC.^{17}

## II. COHESIVE ZONE MODEL OF FRACTURE WITH IMPURITIES

SiC is brittle and has little susceptibility to plastic deformation through dislocation glide, although there is evidence that strong electronic excitations can enhance plasticity.^{18,19} The propagation of transgranular cracks should therefore occur with minimal, if any, crack tip blunting due to plasticity. Hence, the geometry of a transgranular crack will be very sharp as schematically illustrated in Fig. 1. While the atomic planes ahead of the crack tip are crystallographically identical to the atomic planes in the crystal adjacent to the crack path, they are nevertheless subjected to a higher normal stress due to the broken symmetry resulting from the presence of the crack. The absence of plasticity during crack propagation makes it possible to analyze the energies and stresses at the crack tip with a cohesive zone model (CZM) describing the response of the crystal along the planes ahead of the crack tip. The total response of the crystal containing a crack can then be modeled by combining the cohesive zone model with a description of bulk elasticity for the regions adjacent to the crack path.^{20} The slow variation in the crack opening behind and ahead of the tip ensures that there is very little inclination between the decohering atomic planes. Decohesion at the atomic scale can therefore be viewed as the separation of essentially parallel crystallographic planes. This makes the parameterization of a cohesive zone model amenable to first-principles electronic structure approaches.

The first principles analysis of uniform decohesion has been called the *ab initio* tensile test. The method was first introduced to study the effects of defects on the tensile strength of aluminum^{21} and has also been applied to study the mechanical response of grain boundaries^{22–25} and the effects of impurities on decohesion.^{26–31} The competition between elastic and surface energy during uniform decohesion has been analyzed by Nguyen and Ortiz^{32} yielding scaling laws that have found application in other studies of fracture.^{33–35} In this work, we use a thermodynamic formalism that incorporates the effect of impurities at constant chemical potential on the cohesive properties of a crystal.^{30}

Formally, the cohesive zone model should describe the response of the decohering region to externally imposed loads or displacements. This information can be encapsulated within a function that describes the dependence of the internal energy per unit area, *u*, of the decohering zone as a function of a metric of crack opening, or separation, *δ*. The derivative of this energy with respect to separation will then yield a traction-separation curve. In formulating a cohesive zone model, it is necessary to exercise care when defining the cohesive zone energy and the cohesive zone separation *δ*. Here, we restrict ourselves to decohesion between a single pair of atomic planes. Even with this crystallographic simplicity, ambiguity remains. While decohesion occurs between two well-defined crystallographic planes, the atomic relaxations and strain fields in the crystals adjacent to the decohering planes will likely differ from that in the bulk further away from the crack. We can be guided in our choice of definitions for *u* and *δ* by the fact that the cohesive zone model will ultimately be embedded in an elastic medium and that the cohesive zone model together with the elastic description of the adjacent crystal must reproduce the response of the actual crystal containing a crack. This requirement suggests the use of excess quantities, similar to those introduced by Gibbs when formulating a thermodynamic description for grain boundaries and interfaces.^{30,36,37}

Excess quantities account for any deviation from uniform bulk behavior due to local symmetry breaking, either as a result of a grain boundary, an interface, or a propagating crack. An excess energy due to a crack, for example, can be defined as the difference between the total energy of a crystal with a crack and the energy of ideal bulk crystals in the absence of a crack under the same state of local stress. This definition implies that crystalline regions adjacent to the crack path are assigned ideal bulk properties all the way to the atomically sharp crack region, in spite of the fact that non-uniform relaxations may extend into the adjacent crystals. The cohesive zone energy, *u*, set equal to this excess energy, then accounts not only for the energy associated with decohering a pair of crystallographic planes, but also that due to all the deviations from bulk behavior in the adjacent crystalline regions. Similarly, the separation of the cohesive zone *δ* can be set equal to an excess length such that it accounts for the separation of the decohering planes plus any deviation in separation of the adjacent crystals from ideal bulk straining at the same state of externally imposed stress.^{30}

Impurities may be attracted to regions of high stress intensity at the crack tip. They can diffuse to the crack tip from the adjacent bulk crystal. However, in view of the exceedingly low solubility limits of most elements in SiC, a more likely path in this case is through surface diffusion along the newly created surfaces behind the crack tip. The cohesive zone model must account for the effect of impurities in altering cohesive properties in the presence of high impurity concentrations. Restricting the analysis to sub-monolayer impurity concentrations, we can extend the internal energy of the cohesive zone model to also include an impurity concentration dependence, i.e*., u*(*δ*, *θ*), where *θ* is the concentration of impurities per unit area between the decohering planes. While, in general, the impurity concentration in the cohesive zone should also be defined as an excess quantity to account for enhanced impurity concentrations adjacent to the decohering crystallographic planes, the very low impurity solubility in SiC suggests that only appreciable impurity concentrations will be encountered between the decohering planes.

Impurity induced crack growth is usually a slow process and if impurity mobilities are sufficiently high (e.g., through surface diffusion), it is valid to assume that local equilibrium has been reached along the cohesive zone with respect to stress and impurity concentration. This assumption is particularly relevant to the high-temperature applications of SiC. Mechanical equilibrium between the cohesive zone and the adjacent crystals requires an equality of stresses at the cohesive zone/bulk crystal boundary. If the stress state is predominantly tensile, then the cohesive zone will be subjected to a tensile stress *σ*_{t} imposed by the adjacent crystal. If chemical equilibrium exists, then the chemical potentials of the impurities within the cohesive zone must be equal to the chemical potentials within the impurity reservoir from which impurities are drawn. Equilibrium of a system is determined by the state that minimizes a characteristic potential obtained by applying Legendre transforms to the internal energy for every intensive variable controlled as a boundary condition. For a cohesive zone at constant temperature, stress and impurity chemical potential the characteristic potential to be minimized is the grand-force potential defined as

where *s* is the excess entropy per unit area of the cohesive zone. The intensive quantities *σ* and *μ* are related to *δ* and *θ* by

The effects of temperature can be calculated with statistical mechanical approaches that account for configurational, vibrational, and electronic excitations. Contributions to the grand force potential due to configurational disorder when impurities segregate to the cohesive zone, for example, can be treated with cluster expansion approaches.^{30}

## III. PALLADIUM AND SILVER AS IMPURITIES AND FRACTURE GEOMETRY

The energetics of palladium and silver impurities in *β*-SiC have been extensively investigated using Density Functional Theory (DFT) methods.^{38–41} Pd and Ag are very similar in the way they interact with SiC, both within the bulk crystal and on SiC surfaces. For bulk SiC, the formation energies of defects to accommodate dissolved Pd or Ag are, in general, very high. For instance, the interstitial Pd defect has a formation energy in the range of 7 ∼ 9 eV, while substitutional Pd defects have energies in the range of 4 ∼ 6 eV when using the energies of the elements in their stable bulk form as the reference state.^{38} In the case of Ag, the interstitial defects have energies in the range 10 ∼ 11 eV, while substitutional defects have energies in the range 3 ∼ 7 eV.^{39} Interstitial Pd and Ag are predicted to have high mobilities, having migration barriers of the order of 1 eV. The high formation energies associated with dissolving Pd and Ag in SiC are consistent with the absence of any experimental evidence that these metals can be found in solid solution in SiC in any measurable amount.

Silicon carbide is a covalent crystal with the zincblende structure. C and Si are arranged in the form of two interpenetrated fcc lattices, where each C atom is coordinated by a tetrahedron of Si atoms and each Si atom is coordinated by a tetrahedron of C. There are two available interstitial sites that are tetrahedrally coordinated: TC and TSi. The TC tetrahedral sites are coordinated by carbon, while the TSi tetrahedral sites are coordinated by Si. As interstitials, both Pd and Ag have lower energies in the TC sites. The other interstitial site, TSi, is mechanically unstable and is, in the case of neutral impurities, the saddle point of the interstitial diffusion mechanism between neighboring TC sites.

In this work, we focus on decohesion between a pair of (111) planes in SiC in the presence of Pd and Ag impurities. The (111) planes in SiC are very compact, having a low density of atomic bonds connecting adjacent (111) planes and are thus more likely to undergo cleavage before other families of planes do. There are two types of (111) interplanar spaces: a weak one crossed by one bond per atom and a stronger one crossed by three bonds per atom. Decohesion will occur between the weaker interplanar spaces. We will refer to the region between decohering (111) planes as the cohesive zone. The SiC crystal can be viewed as an alternating sequence of pure-Si (111) and -C (111) planes. Decohesion between (111) planes results in a Si-terminated surface and a C-terminated surface. After decohesion, any impurities within the cohesive zone will have a choice of residing on either the Si-terminated surface or the C-terminated surface. The free surfaces of covalent crystals tend to reconstruct after cleavage, thereby lowering their surface energy.^{42} We neglect the role of surface reconstructions as our analysis indicates that the embrittling effect of Pd and Ag impurities occurs before there is a thermodynamic driving force to reconstruct.

## IV. THE CALCULATION OF EXCESS QUANTITIES

The energy of a decohering solid, within the approximation that decohering planes remain parallel, can be calculated with density functional theory using the supercell method. A supercell having a long axis perpendicular to the plane of decohesion and two short axes parallel to the decohering plane can then be elongated in the direction of the long axis, with decohesion localized between a particular pair of atomic planes. All DFT calculations in this work were performed within the generalized gradient approximation using the Perdew-Burke-Ernzerhof^{43} (PBE) exchange correlation functional as implemented in the VASP code.^{44,45} A supercell size of 12 Si + 12 C atoms was used. Energy cutoffs were set at 400 eV, while the k-space grid was set with a single point in the decohesion direction and 8 × 8 points in the perpendicular directions. With these parameters, we obtained the equilibrium lattice parameter of SiC *a*_{0} = 4.38 Å.

The geometry of the supercell used in this work is illustrated in Fig. 2. For the particular configuration shown in Fig. 2, the Pd atom is placed in the TC site and is biased to relax towards the Si-terminated side as the degree of decohesion progresses. We analyzed different decohesion geometries with and without impurities. An analysis of uniform decohesion^{46} of a single crystal shows that there are potentially several local minima with respect to internal degrees of freedom for each elongation of the supercell along the decohesion direction. One local minimum corresponds to a uniform separation of all atomic planes. Another local minimum localizes most of the elongation of the supercell between one pair of atomic planes within the supercell. The final relaxed state in the supercell calculations will therefore depend on the choice of the initial coordinates of the relaxing atoms. To ensure that both these local minima were probed, we initialized the coordinates within the supercell such that the elongation of the cell was either entirely localized between the designated decohering planes, or distributed uniformly over all the atomic planes that were allowed to relax. In all the calculations, N = 8 atomic planes were allowed to relax to account for local atomic relaxations surrounding the cleaving planes that deviate from bulk behavior. There are therefore N + 1 = 9 interatomic spaces undergoing elongation. We denote the total elongation of the supercell by *δ _{T}* and the total energy by

*u*.

_{T}A cohesive zone model relates the energy associated with the pair of decohering planes to the separation of those planes. Since relaxations within several atomic layers away from the cleavage planes are allowed in the supercell calculations, the elastic response of the adjacent planes must be subtracted off from the total supercell energy *u _{T}* and total supercell length

*δ*. To avoid ambiguity, we work with thermodynamic excess quantities

_{T}^{30}and subtract off the (non-linear) elastic response of the homogeneous crystal under the same state of stress along the elongation direction (i.e., the [111] direction) according to

The resulting excess energy *u* and excess separation *δ* then define the cohesive zone model (*δ* is here defined as the separation relative to the equilibrium interplanar spacing in the absence of stress). The energy, *u _{ℓ}*(

*δ*(

_{ℓ}*σ*(

_{T}*δ*))), and separation,

_{T}*δ*(

_{ℓ}*σ*(

_{T}*δ*)), are the energy and elongation of the slabs adjacent to the decohering pair of planes when they are part of a homogeneous solid at the same stress

_{T}*σ*(

_{T}*δ*) that is present in the supercell.

_{T}^{30}

There are several ways in which *u _{ℓ}* and

*δ*can be calculated. In our supercell calculations, 8 atomic planes are allowed to relax. The elastic contribution of these 8 interatomic planes can be calculated using the constitutive stress-strain relationship of homogeneously strained SiC. However, we find that it is numerically more consistent to evaluate the elastic response by using the same supercell geometry (without impurities) allowing N − 1 = 7 atomic planes to relax. In this manner, any finite size effects that are present in the supercell calculations due to the discontinuity in the allowed atomic displacements is systematically cancelled when substracting the elastic contribution from

_{ℓ}*u*and

_{T}*δ*. Figure 3 shows this

_{T}*baseline*elastic response of the non-decohering planes. The dependence of this energy,

*u*, on interplanar displacement

_{ℓ}*δ*was fit with a fourth order polynomial to capture the nonlinear elastic behavior

_{ℓ}The stress, *σ _{ℓ}*, between the relaxing planes is then equal to the derivative of this energy with respect to the displacement, which can be inverted to give the displacement

*δ*as a function of

_{ℓ}*σ*. In Fig. 3, these curves are expressed in terms of the strain

_{ℓ}*ϵ*relative to the length of the interatomic planes at equilibrium. Figure 3 also shows the elastic constitutive curves for homogeneously strained SiC. The discrepancy between the two curves arises from finite-size effects in the supercell calculations due to a transition between planes that are allowed to relax and planes that are held fixed.

In Fig. 4(a), we illustrate the construction of a CZM for (111) decohesion of pure SiC. The green circles relate the total energy *u _{T}* to the total elongation

*δ*of the supercell when the atomic planes that are free to relax are initialized to deform uniformly. The yellow circles relate

_{T}*u*to

_{T}*δ*when the elongation of the supercell is concentrated within the cohesive zone, with local relaxations allowed in the 8 atomic planes surrounding the cohesive zone. Note that due to the existence of two local minima, the supercell energy

_{T}*u*is not necessarily a single valued function of

_{T}*δ*. Substracting the elastic contribution from the adjacent relaxing planes as quantified by Eq. (4) from

_{T}*u*to

_{T}*δ*according to Eq. (3) yields the

_{T}*u*versus

*δ*for the cohesive zone.

The resulting cohesive zone model *u*(*δ*) is shown in Fig. 4 as red diamonds. We parameterized this relationship with an xUBER expression, an extension of the universal binding curve for metals,^{47,48} recently developed for covalent crystals.^{46} The xUBER expression is constructed in a way that allows for a more complex nonlinearity in the elastic response, and has the general form:

where *u*_{0} is the energy at the equilibrium separation *δ*_{0} (relative to the equilibrium spacing in SiC in the absence of impurities), 2*γ* is the energy at infinite separation, and *λ* is a parameter that has dimensions of distance. The parameters *α*_{0} and *α*_{1} are equal to 1 to ensure consistency with the universal binding curve for metals, while the remaining dimensionless parameters *α _{i}* depend on the chemistry and crystal structure. All parameters of the xUBER are fit to the DFT calculated

*u*versus

*δ*relationship.

As an illustration of the effects of local atomic relaxations on the CZM model, we compare the excess energy *u* to the cohesive zone energy calculated in the absence of any relaxations in the adjacent planes (corresponding to the rigid decohesion of SiC along its weak (111) plane) in Fig. 4(b). While for small separations, the curves with and without local relaxations are similar, they start to deviate at larger separations, indicating that relaxations within the planes adjacent to the cleavage planes do affect the response of the cohesive zone. At large separations, the binding energy curve tends to 2*γ*, the work required to create the free surfaces. We obtain values of 2*γ* = 8.37 J/m^{2} for the rigid separation of the crystal and 2*γ* = 7.59 J/m^{2} when atomic relaxations are allowed. Experimental values of the surface energy of SiC are not available. Density functional theory calculations performed by Pastewka *et al*.^{49} under similar conditions (by creating both a Si- and C-terminated surface) yield the value 2*γ* = 8.34 J/m^{2} for the unrelaxed surfaces.

Figure 5 shows a similar construction of a CZM for the decohesion of a pair of (111) SiC planes but now saturated with Pd in the TC-interstitial sites. The energy in Fig. 5 is relative to an equivalent amount of bulk SiC and fcc Pd. These reference states set the zero of the Pd chemical potential to that of fcc Pd. The separations *δ* and *δ _{T}* are defined relative to the equilibrium atomic separation of SiC. The energy minimum at

*δ*

_{0}≈ 1.9 Å represents the equilibrium separation between a pair of (111) SiC interatomic planes saturated with interstitial Pd. The Pd impurities were biased to reside on the Si-terminated side upon decohesion. The same supercell was used as for pure SiC and, as before, eight interatomic planes surrounding the cohesive zone were allowed to relax. As with the CZM for pure SiC, the elastic contribution of the non-decohering planes was subtracted off from the energy and separation of the supercell calculations. The resulting cohezive zone energy versus separation,

*u*(

*δ*), (red squares) was parameterized with an xUBER curve and is shown in Fig. 5 (solid curve). The same approach was repeated for Pd biased towards the C rich surface upon decohesion and for Pd occupying the TSi site and residing on either the Si or C rich surfaces. The difference between the surface energies at large separations with and without impurities equals the adsorption energy of one Pd per surface unit cell. We obtain values of 0.80 eV for TC and 0.91 eV for TSi, which are in agreement with the values previously reported in the literature by Shuck and Stoller

^{40}(0.93 eV for TC and 1.04 eV for TSi). We also performed similar calculations for Ag TC interstitials, as this type of defect has the lowest free energy up to the maximum theoretical stress.

## V. EQUILIBRIUM DURING DECOHESION

Equilibrium at a constant impurity chemical potential *μ* and applied stress *σ*, is determined by minimizing the grand-force potential $\varphi $, Eq. (1). In this section, we compare the relative stability of a cohesive zone without Pd impurities (*θ* = 0) and a cohesive zone saturated with interstitial Pd impurities (*θ* = 1/*A*_{0}, where *A*_{0} is the area per two-dimensional unit cell of the decohering plane). A more complete description would also consider intermediate levels of impurity segregation using, for example, a cluster expansion approach,^{30} but our analysis shows that the fully saturated extreme is thermodynamically viable at realistic Pd chemical potentials.

Figure 6 shows the grand-force potential as a function of applied stress for pure SiC without any impurities and for a cohesive zone with all TC tetrahedral sites filled with Pd. The chemical potential of Pd is set equal to that of its fcc metallic form (*μ*_{Pd} = 0 for the reference states used here). For comparison, we also show a grand-force potential for a (111) cohesive zone saturated with Ag impurities (*μ*_{Ag} = 0 as well). It is convenient to interpret the variation of the grand-force potential with stress by considering the separation *δ* as an internal parameter. As an example, consider the grand-force potential of pure SiC without any impurities. Starting from the equilibrium separation, *δ*_{0} at *σ* = 0, the grand force potential of the cohesive zone first decreases, having a negative curvature until the theoretical maximum stress is reached at the cusp. Upon a further increase of *δ*, the corresponding stress *σ* decreases from its maximum value, while the grand force potential increases with positive curvature. The curvature of the grand force potential is related to minus the stiffness of the cohesive zone. A negative curvature indicates a regime where the cohesive zone is mechanically stable, while a positive curvature indicates a negative stiffness and thus a cohesive zone that is mechanically unstable. In the limit of infinite separation, we return to the *σ* = 0 axis and the difference between the two branches of the grand-force potentials at *σ* = 0 is the work performed in creating the two free surfaces.

The equilibrium state of the cohesive zone is determined by a global minimum among all the grand-force potentials. Whether the cohesive zone is saturated with impurities or not depends on which grand-force curve is minimal. Figure 6 shows that at zero stress, the state of the cohesive zone with the lowest grand force potential corresponds to pure SiC, or *θ* = 0. However, the slope of the grand force potential for the cohesive zone saturated with interstitial Pd, or *θ* = 1/*A*_{0}, is more negative. In fact, it can be shown that the slope is given by the negative of the equilibrium separation,−*δ*_{0}. This causes the grand force potential for the Pd saturated cohesive zone to cross that of pure SiC at a stress substantially below the maximum stress of the pure SiC cohesive zone. This occurs both for TC and TSi- type interstitials, although TC interstitials have a lower grand-force potential for stresses below the maximum. The crossing of free energies for distinct states (or phases) signifies a first-order thermodynamic phase transformation. In the present example, the crossing of the grand force potentials at *σ*_{tr} indicates a thermodynamic driving force for a stress induced phase transformation within the cohesive zone.

The equilibrium traction curve of the cohesive zone at constant Pd chemical potential *μ*_{Pd} = 0 (the Pd chemical potential in the presence of fcc Pd precipitates), as determined by the minimum envelope of the grand force potential, is shown in Fig. 7. Once the transformation stress *σ*_{tr} of 16.3 GPa is reached, corresponding to the crossing of the pure SiC and Pd saturated grand-force potentials, the separation increases discontinuously at a constant stress *σ*_{tr} to the traction curve for the Pd saturated cohesive zone. This discontinuous increase in separation is accompanied by an abrupt change in equilibrium Pd concentration within the cohesive zone from *θ* = 0 to *θ* = 1/*A*_{0}. As is clear from Fig. 7, the maximum stress for decohesion of the Pd saturated cohesive zone is substantially lower than that of pure SiC. Hence the stress induced transformation within the cohesive zone, requiring an influx of a monolayer of Pd between the decohering planes, leads to an embrittlement of SiC. In contrast, Fig. 6 indicates that such embrittlement mechanisms will not occur in the presence of Ag impurities at the chemical potential of its fcc bulk phase, *μ*_{Ag} = 0 (when using fcc Ag as our reference state).

## VI. DISCUSSION

Impurities can alter the fracture toughness of a material by directly altering the bonding between the atoms of the crystal. Here, we have developed a first-principles thermodynamic description of the cohesive zone ahead of a sharp crack tip subjected to boundary conditions imposed by the adjacent elastic medium and an impurity reservoir. We have applied it to analyze decohesion of SiC along a pair of (111) planes in the presence of Pd and Ag at constant chemical potential, the thermodynamic boundary condition in the limiting case of a very mobile supply of Pd and Ag. Using DFT input, our analysis has shown that for Pd chemical potentials equal to that of the crystalline metallic state, there is a thermodynamic driving force for Pd to saturate the interstitial sites between separating (111) planes of SiC and thereby lower the maximum stress of decohesion. This suggests that SiC in contact with precipitates of Pd is susceptible to embrittlement provided Pd can diffuse to the crack tip. A similar driving force does not exist in the presence of metallic Ag.

The impurity embrittlement mechanism predicted for Pd is the result of a first order phase transition within the cohesive zone. It emerges due to a crossing of two distinct free energies and involves a discontinuous change in the extensive quantities *δ* and *θ*. For such a process to occur there must be a kinetic pathway that is likely accompanied by energy barriers. As with other first-order phase transitions, this embrittlement mechanism will likely occur through a nucleation and growth mechanism.

The mechanism with which Pd impurities embrittles SiC is by reducing the maximum stress for decohesion, which occurs before free surfaces are formed. As a result, surface reconstructions, while lowering the energy of decohesion, do not participate in this particular embrittlement mechanism. In fact, we find that the grand force potentials for cohesive zones having low energy surface reconstructions only become thermodynamically favored at separations that are much larger than the separations at which the Pd reduces the maximum stress for decohesion. We are therefore justified in neglecting the energetic constributions of surface reconstructions in our thermodynamic analysis.

The tendency of an impurity to embrittle a crystal depends on the value of its chemical potential. In the current treatment, embrittlement occurs when the grand force potential of the impurity saturated cohesive zone crosses the grand force potential of pure SiC. The grand force potential of the impurity saturated cohesive zone is very sensitive to the impurity chemical potential, rigidly shifting down as the chemical potential increases. While the grand force potential for a Pd saturated cohesive zone crosses that of pure SiC when in equilibrium with metallic Pd, that of a Ag saturated cohesive zone in contact with metallic Ag does not (Fig. 6). Under normal thermal equilibrium conditions, it is unlikely that a Ag reservoir having a chemical potential substantially above that of metallic Ag can be devised such that Ag would have a thermodynamic driving force to embrittle SiC. In nuclear applications, however, high concentrations of Ag may be implanted in SiC at concentrations far above its solubility limit, thereby generating a source of Ag with an effective chemical potential far above that of metallic Ag precipitates. Under these conditions, the possibility exists that even Ag may embrittle SiC.

Both Pd and Ag are fission products of the uranium oxide constituting the core of TRISO particles. These fission products are implanted in the inner graphitic shell covered by a SiC shell. Due to its low solubility, supersaturated Pd will form metallic precipitates, which will act as a reservoir of Pd at a chemical potential coinciding with that of metallic Pd. As our first-principles analysis predicts, Pd at this chemical potential will embrittle stressed (111) planes ahead of pre-existing cracks by lowering the critical stress for decohesion. Given that Pd is known to form silicides,^{9} it is also possible that a stressed SiC crystal will be susceptible to stress induced phase transformations involving any of the Pd-Si compounds. This, however, will require the concerted diffusion of not just Pd but also Si and C, and should therefore occur more slowly. Nevertheless, in TRISO particles, irradiation could enhance the mobilities of Pd, Si, and C such that silicide formation may compete with Pd embrittlement. The creation of fresh new SiC surfaces as cracks advance due to Pd embrittlement, though, are likely the most favorable sites for Pd-silicide nucleation and growth.

Our predictions suggest a mechanism by which fission products such as Ag, Cs, Sr, and I among others can escape the SiC shell of TRISO particles. While solubilities and mobilities of these fission products in SiC are predicted to be exceedingly low, transport along micro and nano-cracks will occur much more rapidly. Pd is always present as a result of fission, and our thermodynamic predictions here indicate that Pd can embrittle SiC by substantially reducing the critical stress for decohesion. A pair of SiC layers can locally be under high levels of tensile stress, either from residual growth stresses or due to the presence of thermal gradients. While these stresses might be below the critical stress for decohesion of pure SiC, they may nevertheless exceed the critical stress when Pd segregates ahead of crack tips. Depending on the distribution of internal stresses within the SiC shell, Pd embrittlement could open up a short-circuit diffusion path for other fission products.

We have restricted our analysis to decohesion between (111) planes in SiC. Similar embrittlement mechanisms could occur between other crystallographic planes and along grain boundaries. Our prediction that Pd embrittles SiC for (111) decohesion at realistic Pd chemical potentials suggests that it is a likely mechanism of SiC degradation and failure in TRISO particles. It is possible that other crystallographic planes and grain boundaries decohere at even lower stresses in the presence of Pd. If this is the case, it would make this mechanism of embrittlement even more likely than is suggested by the results reported on in this study. Our thermodynamic analysis was also restricted to charge neutral cohesive zones. Since SiC is an insulator, internal defects, including grain boundaries and impurity clusters segregated to stress concentrations, may attract excess charge from the bulk crystal depending on the Fermi level. To account for the effect of excess charge within the cohesive zone, it would be necessary to apply an additional Legendre transform to the grand force potential with respect to the electron chemical potential (Fermi level) and number of excess electrons within the cohesive zone. Supercell calculations would also need to be performed having different states of excess charge. An analysis of equilibrium within the cohesive zone would then occur by minimizing this new grand force potential not just as a function of stress and impurity chemical potential, but also as a function of Fermi level. Since charged impurities add a new degree of freedom to the cohesive zone, its inclusion in a thermodynamic analysis of equilibrium can only result in predictions of embrittlement occurring at lower stresses and impurity chemical potentials than already predicted in this study. Our study also neglects contributions of entropy to the grand force potentials and therefore corresponds to the low temperature limit. The effect of impurity configurational disorder within the cohesive zone, for example, can be accounted for with cluster expansion techniques^{30} and is likely to favor impurity segregation to regions of high stress concentration, thereby lowering the impurity chemical potentials at which embrittlement can occur.

The analysis presented here can also be applied to study the tendency of other impurities such as oxygen, hydrogen, or water molecules to embrittle SiC. This type of environmental damage has particular relevance to the design of SiC fiber composites.^{50}

## VII. CONCLUSION

We have investigated the decohesion of SiC in the presence of mobile Pd and Ag impurities by combining first-principles density functional theory calculations with a thermodynamic analysis of decohesion at constant impurity chemical potential. We find that Pd impurities can substantially reduce the critical stress of decohesion of SiC when in equilibrium with metallic Pd. The embrittlement mechanism relies on a stress induced phase transition requiring an influx of Pd between the decohering planes.

Our thermodynamic analysis is based on a cohesive zone model that is parameterized with density functional theory calculations of supercell geometries. The cohesive zone model is defined in terms of thermodynamic excess quantities, whereby the elastic response of the atoms adjacent to the decohering planes is subtracted from the total response of the decohering solid. We have restricted our focus to decohesion between (111) planes of SiC, however, our prediction that Pd can reduce the maximum stress for decohesion at realistic chemical potentials indicates that other crystallographic planes and grain boundaries may embrittle at even lower stresses and impurity chemical potentials than predicted in this study.

## ACKNOWLEDGMENTS

This work was funded by DOE under NEUP of the Idaho National Laboratory, Award No. 103195, Project No. 10-924.