Accommodation and Diffusion of Nd in Uranium Silicide-U 3 Si 2

Uranium silicide, U3Si2, is considered as an advanced nuclear fuel for commercial light water reactors with improved accident tolerance as well as competitive economics. Nd is employed as a local burnup indicator for conventional oxide fuels due, among other reasons, to its low mobility in the UO2 fuel matrix and its high fission product yield. As part of the studies necessary to determine whether Nd can be considered as a candidate burnup indicator in the U3Si2 concept fuel, we investigate the mobility of Nd in U3Si2. In this work, density functional theory (DFT) calculations are performed to predict the most stable accommodation sites of Nd in U3Si2, found to be within the uranium sublattice. Based on DFT calculations of binding energies and migration activation energies, we investigate Nd diffusion by computing the transport coefficients within the framework of the self-consistent mean-field method. Our calculations predict that the diffusion ratio of Nd to U is smaller in U3Si2 than in UO2. Moreover, at the individual maximum centerline temperature of the fuel, the diffusion of Nd in U3Si2 is much slower than in UO2. From this perspective, Nd represents a good candidate burnup indicator, in similarity to that in UO2.


Introduction
Quantification of fission products is one of the methods used to evaluate the fuel consumption in a nuclear reactor, commonly known as burnup. A well-established method for the determination of burnup is the fission product monitoring method [1]. In this method, a suitable fission product isotope is selected, and its concentration is measured during post-irradiation examination (PIE) and then correlated with residual uranium and plutonium. The isotopes to serve as burnup indicators should satisfy certain requirements such as: being long-lived or stable, having very little mobility in the irradiated fuel matrix, a high and constant fission yield, and negligible transmutation rates by nand γ-reactions [2]. 148 Nd satisfies most of these requirements [3,4] and is well accepted as a burnup indicator for uranium oxide fuels since the issue of ASTM-E321 [5]. It is worth noting that, unlike other requirements related to the intrinsic properties of the isotopes, the diffusion behavior of atomic species in crystalline materials is strongly dependent on the crystal structure. The mobility of Nd in a fuel matrix with a different crystal structure and atomic configurations could thus differ substantially. Therefore, to employ Nd as a burnup indicator for other types of fuel, its mobility must be studied.
Advanced nuclear fuel concepts are being investigated for use in light water reactors in order to improve fuel performance, economics, and safety. U 3 Si 2 is a well-studied example and one of the most promising accident-tolerant fuel candidates for commercial LWRs [6] with both a wide variety of experimental and modeling investigations undertaken. Unlike the conventional fuel system, UO 2 that has a cubic fluorite crystal structure, the U 3 Si 2 crystal structure is tetragonal with the space group P4/mbm. Many crucial properties of U 3 Si 2 , such as thermal conductivity, thermal expansion, mechanical properties, and phase stability have been successfully determined and are impacted by the anisotropic structure of the material [7][8][9][10][11][12][13][14][15]. Recently, self-diffusion and fission gas diffusion in U 3 Si 2 under thermal conditions were investigated by Andersson et. al. [16,17]. U and Si selfdiffusion, as well as Xe diffusion, are predicted to be faster in U 3 Si 2 than U self-diffusion and Xe diffusion in UO 2 . However, to our knowledge, there is no study yet on diffusion of solid fission products (such as Nd and other elements) in U 3 Si 2 .
Diffusion behaviors could be studied using well-established laboratory technologies or state-of-theart first-principle calculations. However, there is no general rule for assessing whether the diffusivity is low or high in absolute terms. Since Nd is widely used as a burnup indicator for UO 2 and its mobility is therefore low enough, we take the mobility of Nd in UO 2 as a reference. The direct comparison of Nd diffusion coefficients in different crystal structures is, however, of little value. Instead, we can compare the diffusion length scale over a given time in U 3 Si 2 and in UO 2 at their respective operating temperatures. Moreover, we can define a ratio of interest, which is the Nd diffusion coefficient over U self-diffusion coefficient to determine the Nd diffusivity with respect to the host atoms and compare this ratio in U 3 Si 2 to that in UO 2 . Previous work [18,19] demonstrated that Nd diffusion in UO 2 is 10 7~1 0 8 times faster than U self-diffusion. This means that the study of diffusion behaviors is not sufficient since, in addition to bulk diffusivity, Nd in the bulk can also form precipitates, interact with dislocations or other traps, be involved in fission spikes, etc.
Nevertheless, computing the diffusion coefficient can provide a first-approximation criterion to evaluate the mobility of Nd and its potential use as a burnup indicator.
In the present work, we study the stability and thermal diffusion of Nd in U 3 Si 2 by combining density functional theory calculations, with a Hubbard correction (DFT+U), and the self-consistent mean field method [20] implemented in the KineCluE code [21]. The Nd diffusion length scale at relevant LWRs operating temperature is calculated. The diffusivity is analyzed and evaluated in relation to self-diffusion in U 3 Si 2 and diffusion of Nd in UO 2 . The aim of this work is to determine whether the thermal diffusion of Nd in U 3 Si 2 is as low as in UO 2 , so that Nd can be considered as a candidate local burnup indicator for U 3 Si 2 , in a general effort of assessing the fuel performance of U 3 Si 2 . The methods and results presented here may also be useful for further studies on fission products behavior in U 3 Si 2 .

Equilibrium point-defect concentration and Nd accommodation
As shown in Fig. 1 (a), the U 3 Si 2 crystal structure has two symmetry-unique uranium sites, the U1 site with Wyckoff position 2a and the U2 site with Wyckoff position 4h, and one silicon site with Wyckoff position 4g. Interstitial positions with Wyckoff positions 2b, 2c, 2d, 4e, 4f, 4g, 4h, 8i, 8j, and 8k are illustrated in Fig. 1.(b). Si at 2b, with fractional coordinates (0.0, 0.0, 0.5) in the reference frame of the unit cell, was found to be the lowest-energy Si interstitial defect structure [13]. U at an 8j site with coordinates (0.90, 0.95, 0.5), pushing away the two neighboring atoms and forming a delocalized split structure involving 3 atoms asymmetrically occupying 2 sites, is the lowest energy structure for U interstitials, see in previous DFT studies [13,16,17] is predicted to be less stable than the U (0.90, 0.95, 0.5) interstitial position considered in this work (the formation energy of the former being 2.5 eV higher than the latter according to our DFT calculations). Diffusion of substitutional Nd atoms is enabled by point defects (either vacancies or self-interstitial atoms). Interstitial Nd atoms can diffuse independently of defects to neighboring interstitial sites.
Consequently, the vacancy and interstitial defect concentrations in U 3 Si 2 are needed in order to compute the U and Si self-diffusion coefficients, which are used in turn to evaluate the mobility of Nd in U 3 Si 2 . The equilibrium point-defect concentrations depend on the Gibbs formation energy G D.f , which is calculated by means of ab initio calculations from the defect formation enthalpy and entropy, . = . − . . The defect formation enthalpy is expressed as: where H tot is the total enthalpy, and superscripts D and 0 represent the supercell with and without defect, respectively. µ i is the chemical potential of species i (U or Si). ∆N i is the number of atoms of the species i that are added (∆N i > 0) or removed (∆N i < 0) from the perfect supercell to create the defect, and the sum runs over all added and removed species. We approximate the chemical potentials µ U and µ Si by assuming equilibrium between two adjacent phases in the phase diagram, with the generic nomenclature U a Si b and U c Si d and solving the following equations [16,17]: The defect formation entropy can be calculated with the approach of Mishin et al. [22] which relates point-defect entropies with lattice vibrations in the harmonic approximation: where ν n are the lattice vibrational frequencies (collected in Table A.1), N±1 and N are the number of atoms of the perfect supercell with and without defect (plus for interstitial and minus for vacancy), respectively. s i is the entropy of either U or Si calculated from: where N is the number of atoms in the cells used to describe the adjacent U a Si b and U c Si d phases.
Three different atomic environments are used to calculate the chemical potentials: Si-rich, U-rich, and perfect stoichiometry. For Si-rich conditions, we assume U 3 Si 2 is in equilibrium with USi. The obtained chemical potentials are applied to the system with either a U vacancy or a Si interstitial.
For U-rich conditions, U 3 Si 2 is in equilibrium with U 3 Si. and The Nd accommodation in U 3 Si 2 can be investigated by computing the incorporation energy E inc and solution energy E sol . The former indicates the possibility to accommodate Nd in a pre-existing vacancy caused by neutron irradiation and can be expressed as: where + represents the total energy of the supercell with a Nd solute and a vacancy, that of a supercell with one vacancy only, and µ Nd the chemical potential of Nd which is derived from the equation: is the total energy of NdSi, i.e., the phase which would form if enough Nd were present [23]. A negative E inc thus indicates that the incorporation of a Nd atom in a pre-existing vacancy is energetically favorable. On the other hand, the solution energy is calculated as the sum of the vacancy formation energy and the incorporation energy of the solute in this vacancy: and reflects the most likely accommodation site for a solute when no defects are already present in the lattice.

Diffusion and transport coefficients
The self-consistent mean field (SCMF) theory [20], as implemented in the KineCluE code [21] is an analytical statistical-physics framework that allows for the calculation of transport coefficients at equilibrium, based on the knowledge of atomic jump rates. In this model, the presence of thermodynamic driving forces, e.g., chemical potential gradients (CPG), causes a perturbation of the equilibrium configuration and produces atomic and defect fluxes. The transport coefficients reflect the kinetic response of the system to the driving forces. The Onsager matrix allows for the analysis of flux coupling between different species (vacancies V and solutes S) and its components L ij are defined as: The matrix is perfectly symmetric as long as microscopic detailed-balance conditions are fulfilled.
KineCluE ensures that such conditions are fulfilled at all times for transport-coefficient calculations [21]. In Eq. (14), J i denotes the atomic flux of species i and ∇µ j the CPG of each species j. In our case, the transport coefficients for vacancy-assisted diffusion are , = and . The off-diagonal coefficients and describe the flux of one chemical species induced by forces acting on other species (either vacancies or solute), allowing therefore for an accurate analysis of the flux-coupling tendency and the mutual directions of the two fluxes. For interstitial diffusion, the independent transport coefficient is . Our interstitial transport coefficient calculation is valid as long as possible correlation effects between interstitial solutes and vacancies or other species are negligible. This would not hold in the presence of solute interstitial-vacancy binding interactions, in which case a cross-term would be involved; this is the case for instance for foreign impurities (C, N, O) in Fe alloys [24]. Further studies will be needed to evaluate the effect of such possible correlations.
Under the assumption that cross terms can be neglected, and in the dilute-limit approximation, the tracer diffusion coefficient of a solute that migrates through mechanism δ (δ = vac for the vacancyassisted mechanism and δ = inter for the interstitial mechanism) can be directly derived from the transport coefficient: where C S is the solute concentration, and P δ the probability of the configuration that allows for the migration of solute through the mechanism δ. For the vacancy mechanism, the probability (under certain assumptions) can be approximated as: where [V] is the probability of forming a vacancy (i.e., the vacancy concentration), and Z VS the partition function related to the binding energy of the vacancy-solute pair as = ∑ exp ( ) (index i marks all the possible configurations of the vacancy-solute pair). E b is the vacancy-solute binding energy which can be computed by ab initio methods with the following equation: The terms in Eq. (17) are in turn the total energy of the supercell with a single vacancy, a single solute, a vacancy-solute pair, and the perfect supercell. Positive binding energies stand for attraction between the two defects, and negative binding energies stand for repulsion.
For the simple interstitial mechanism, the solute is the interstitial defect and has a probability equal to 1 of being in an interstitial position, so = . In this work, there are three distinct Nd The flux-coupling character of the vacancy-assisted diffusion can be analyzed by computing the drag ratio = ⁄ and the partial diffusion coefficients (PDC) ratio which is given as [25]: where = − − , and A and B represent the host and foreign atom, respectively. In the absence of attractive interactions between defect and solute, the atomic flux is opposed to the vacancy flux, because the vacancy moves by exchanging with either solute or host atoms. However, when attractive solute-defect interactions are present, the vacancy may remain in proximity of the solute and produce consecutive solute jumps in the same direction: this is usually referred to as vacancy drag. The PDC ratio describes the relative diffusion rate of solute with respect to matrix atoms, which determines the solute segregation or depletion tendency at defect sinks (e.g. dislocation lines, grain boundaries, free surfaces, and so on). Vacancy drag leads to solute enrichment at sinks, in which case the PDC ratio D pd is negative. In the absence of vacancy drag, two regimes are possible: when 0< D pd <1, solutes diffuse slower than host atoms, and enrichment at sinks takes place; on the other hand, if D pd >1 , solutes are faster than host atoms, and depletion thus occurs.

Multi-mechanism diffusion coefficients
Self-diffusion can occur either by vacancy or interstitial mechanism. The total self-diffusion coefficients result from the sum of diffusion coefficients of the two mechanisms: Nd in U 3 Si 2 diffuses through multiple mechanisms: uranium vacancy-assisted migration, silicon vacancy-assisted migration, and interstitial migration. The total Nd diffusion coefficient can be obtained by taking into account the probability of Nd to be located at each trap site (x= U1, U2, Si, inter). It is thus given as: where is the solution energy of Nd at the trap site which can be calculated using Eq. (13).

Kinetic analysis
The kinetic analysis of this work, i.e., the calculation of transport coefficients, is performed with the KineCluE code [21]. This code implements a cluster-expansion approach to the SCMF where the Onsager matrix is split into cluster contributions. A detailed explanation of the theory and workflow of this code can be found in [21,25]. In the dilute approach applied in this work, two "clusters" are considered for vacancy-assisted diffusion: a mono-vacancy and a vacancy-solute pair, given that an isolated substitutional solute is immobile. Any multi-vacancy and multi-solute effects were neglected, as a first approximation, although their possible contribution to Nd mobility should be further investigated by extending the set of DFT calculations to Nd-Nd, Nd-vacancy, and vacancyvacancy pairs (which is beyond the scope of this work).
In the SCMF theory, one needs to impose a cut-off range for thermodynamic interactions, R th , as well as one for kinetic interactions, R kin . Thermodynamic interactions embodied by the binding energies determine the probability of a certain solute-defect configuration to occur in thermodynamic equilibrium conditions, and the binding energies of solute-defect pairs whose distance is beyond R th are regarded as negligible. Kinetic interactions, on the other hand, are fictitious interactions introduced in the SCMF framework to parameterize the near-equilibrium distribution function. Thus, R kin is the maximum extension of the diffusion trajectories included in the model. Migration paths extending further than R kin are not considered. The transport coefficients converge for longer kinetic radii, but in practice we need to keep R kin as short as possible, taking the computational effort into account. Previous tests [25] confirmed that a kinetic radius of a few Å for vacancy and simple interstitial mechanisms is sufficient to obtain a satisfactory convergence.
Therefore, we set the kinetic radius R kin to 4 a 0 , and the thermodynamic radius R th to √2 2 0 , corresponding to the first nearest-neighbor distance for U1 sites.
In order to construct continuous diffusion trajectories for vacancy-assisted diffusion mechanisms, the migration barriers of three different types of jumps are needed: Nd-vacancy exchanges, vacancy-U and vacancy-Si exchanges with a Nd nearby, as well as vacancy-U and vacancy-Si exchanges in pure U 3 Si 2 . In this work, the migration barriers of Nd-vacancy exchanges ( ) and of vacancy-U (or Si) exchanges in pure U 3 Si 2 ( 0 ) are calculated with DFT (see Table 2), while the migration barriers of vacancy-U (or Si) exchanges with a Nd nearby are calculated automatically by the KineCluE code with the KRA approximation [26]: The latter migration barriers are thus obtained based on the binding energies of the initial ( ) and final ( ) state, and a reference migration barrier corresponding to 0 ( ) or 0 ( ) for vacancy-U and vacancy-Si exchanges, respectively.
Because of the tetragonal crystal structure of U 3 Si 2 , the diffusion behaviors in the a-b plane and along the c-axis are distinct. Therefore, we need to compute the diffusion coefficient in each direction. This is done by setting first the CPG along (110) direction (chosen as a representative direction in the a-b plane), then along the (001) direction. In addition, due to anisotropy, the CPG in the (110) direction may cause a diffusion also in the (001) direction, and vice versa. However, this second-order effect was here found to be several orders of magnitude (10 7 to 10 9 ) lower than the diffusion coefficient in the main directions and was thus neglected. The solute concentration C S is set to 10 -4 in line with the mole fraction of Nd in U 3 Si 2 fuel for a burnup of 52.5 MWd/kg, which is predicted by Serpent calculations [27]. This concentration ensures that Z 0 (pair)*C S < Z 0 (mono) where Z 0 is the count of possible geometric configurations taken by the vacancy-solute pair or mono-vacancy within R kin . Even though the solute concentration is set to this fixed value, both the solute diffusion coefficient and the drag ratios are independent of C S [25].
To limit the computational expense, all the attempt frequencies are set to the Debye frequency of U 3 Si 2 (11.19 THz), which is obtained from our DFT calculations. Due to a current KineCluE limitation that prevents the modeling of anti-sites, we could not consider the migration of a vacancy from the U to the Si sublattice (and vice versa) that involves the creation of an anti-site, although our DFT calculations showed that this jump is possible and has a relatively low migration barrier (see Section 3.3 for further details). In addition, we ignore the possible Nd migration across the two sublattices. In principle, Nd needs only one vacancy to jump from one sublattice to another.
However, once the jump to the other sublattice has occurred, an additional vacancy in that sublattice is needed to ensure a continuous diffusion trajectory. With one vacancy only, Nd would exchange across sublattices back and forth without producing a net migration. In order to take this mechanism into consideration, we should thus consider configurations including two vacancies (one per sublattice), and compute as a consequence multi-defect properties that go beyond the dilute approach targeted in this work. For these reasons, we investigate the diffusion of Nd in the uranium sublattice and in the silicon sublattice separately. On the other hand, the diffusion of Nd through U1 and U2 sublattices can be treated simultaneously. As a final remark, while using the KineCluE code, an energy modification using the most stable configuration as a reference is applied to the nonequivalent configurations in mono-vacancy and vacancy-solute calculations, respectively (more details are provided in Appendix B).

Density functional theory calculations
The DFT calculations in this work were performed with the Vienna Ab Initio Simulation (VASP) package [28][29][30][31]. The projector augmented wave method [32] and the Perdew-Burke-Ernzerhof (PBE) parametrization of the generalized gradient approximation (GGA) [33] were used to describe the wave function and the exchange-correlation function respectively. The potentials used in this work treat 14, 14, and 4 electrons as valence for, respectively, uranium, neodymium, and silicon. To capture the strong correlation effects of the uranium 5f orbitals, the rotationally invariant implementation of the Hubbard-U correction introduced by Dudarev et al. [34] was applied.
Following our previous study [35], an effective Hubbard-U value, U eff =U-J, of 1.5 eV was applied to the ferromagnetic structure of U 3 Si 2 . The Hubbard correction was not applied to the neodymium 4f orbitals in this work since previous research has confirmed that Nd-bearing compounds can be described well without Hubbard correction [36]. Moreover, in this work, the calculation results were obtained as differences with respect to the reference states, thus the relative effect of applying U-correction to the single Nd can be regarded as rather small.
The defect energies, binding energies, and migration barriers were calculated in a 2×2×3 supercell containing 72 uranium atoms and 48 silicon atoms. The Brillouin zone was sampled with 4×4×4 Monkhorst-Pack k-point meshes. The plane-wave cut-off energy was set to 500 eV and the partial occupancies were smeared according to the Methfessel-Paxton method with a smearing width of 0.2 eV. To compute the defect energies, the atomic positions, supercell volume and supercell shape were allowed to fully relax. The convergence criteria for the self-consistent ionic relaxation loop was 10 -4 eV. Migration barriers were calculated using the climbing image nudged elastic band (NEB) method [37][38][39] with three images. The spring constant between the images was set to 5 eV/Å 2 . The NEB calculations were performed with the volume fixed at that of the starting configuration. The force convergence criteria for NEB calculations was 0.05 eV/Å.
The vibrational frequencies at the Γ point were calculated via density functional perturbation theory (DFPT) [40]. To maintain the ferromagnetic structure of the system, the wave function of the ground states was used as the starting point for the DFPT calculations. Considering the high computing cost of DFPT calculations, we used a 1×1×2 supercell of U 3 Si 2 . The k-point grids were chosen to be 8×8×16. The energy convergence criteria were set at 10 -8 eV.
In order to calculate the chemical potentials and the corresponding entropies of U and Si, the energies and vibrational frequencies of the neighboring phases, U 3 Si with space group Fmmm and USi with space group Pbnm, were calculated. The chemical potentials of Nd were derived from NdSi in the iron-boride (FeB) crystal structure. Calculations of these compounds used the same plane-wave cut-off energy as U 3 Si 2 and the density of the k-point grid was chosen to ensure the same accuracy across all systems. All other computational details were the same as for U 3 Si 2 .

Point-defect concentrations in U 3 Si 2
The calculated defect formation enthalpies and entropies are collected in Table 1. The chemical potentials and the corresponding entropies of U and Si were calculated using Eqs. (2)(3) and (5-6) respectively. The defect concentrations calculated using Eqs. Andersson et al. [16,17] are shown for comparison. Here, the "Si-rich" means that either adding one Si interstitial, or removing one U atom (U vacancy), while "U-rich" means either adding one U interstitial or removing one Si (Si vacancy). These results are for most parts in good agreement with those of Anderson et al. [16,17]. The slight discrepancies appearing in Fig. 2 can be explained as follows. For U1 vacancy, Si vacancy, and U interstitial concentrations, they mainly originate from the differences between the defect formation entropies (see Table 1). For U2 vacancy and Si interstitial concentrations, the discrepancy originates mostly from enthalpy differences, although significant entropy differences are also observed.
Compared to Ref. [16,17], where all point defect systems considered as reference states for the chemical potential calculations were approximated as stoichiometric, we treated the systems as nonstoichiometric, i.e., we used U-rich, and Si-rich phases as references, resulting in different chemical potentials of U and Si. Besides, in Ref. [16,17], the finite displacement method in non-cubic 2×2×2 or 1×2×3 supercells were used to compute the vibration frequencies, while in this work the latter were computed with the DFPT method in the cubic-symmetry 1×1×2 supercells. This is the main reason for the substantial differences between our calculated defect entropies and those in Ref. [16,17]. Note that the well-known metastability issue induced by the strong correlation corrections for uranium 5f elections could also be the reason of the discrepancy. In this work, the occupation matrix control scheme (OMC) [41] is applied to find the ground state of the perfect supercell which is crucial for the point defect energy calculations. This was shown in an earlier study to be potentially important since metastable states with total energy differences of nearly 1 eV could appear without this correction [35]. No correction scheme was applied in Ref. [16,17] and their calculations may have converged to metastable states leading to different point defect energies. 3.2 Self-diffusion in U 3 Si 2 The self-diffusion mechanisms by vacancies and interstitials are discussed in Appendix A.
Migration energies 0 calculated with NEBs are collected in Table 2 Fig. 3 (a) shows the total self-diffusion coefficients calculated using Eqs. (19) and (20), Fig. 3 (b) collects the contributions of each diffusion mechanism. As shown in Fig. 3 (a), the diffusion of U is systematically faster than that of Si in both directions. For each species, diffusion in the a-b plane is faster than along the c-axis. This relative relationship is consistent with the calculations of Andersson et al. [16,17], but our calculated self-diffusion coefficients are lower in absolute terms.
The discrepancy lies in the differences between migration barriers and between defect concentrations. Since the discrepancies between the migration barriers listed in Table 3 (left) for Si interstitials and vacancies are much larger than for U interstitials and vacancies, the mismatch is more profound for the Si diffusion coefficient. This discrepancy is also augmented by the difference between the Si interstitial and vacancy concentrations, and can be explained by the same arguments as discussed in Section 3.1 concerning the different DFT calculation parameters. . This combination gives rise to a potentially long-ranged diffusion path with a low total energy barrier, and thus makes the U interstitial diffusion faster than any of the vacancy mechanisms. Si interstitials along the c axis move much slower than in the a-b plane owing to the higher migration barrier of I Si-c , but still dominate the Si diffusion. U-vacancy diffusion is faster than U interstitials along the c-axis, and is thus the prevalent U self-diffusion mechanism along the c axis. This can be attributed to the lower-barrier jump V U1-c and the multi-step mechanism V U2-U1-U2 as illustrated in Fig. A. 1 (b).

Solubility of Nd in U 3 Si 2
The solubility of Nd was modeled considering U1, U2 and Si substitutional sites, the interstitial site with Wyckoff position 2b, and two non-equivalent interstitial sites with Wyckoff position 8j.
Although there are more possible interstitial sites in U 3 Si 2 as illustrated in Fig.1. (b), 2b and 8j sites were observed to be the most stable [42]. The calculated solution energies E sol and incorporation energies E inc are summarized in Table 3. U1 and Si are found to be the most favorable trap sites for Nd, followed by U2 and 8j-2 sites. When vacancies are pre-existing, for instance under irradiation conditions, U2 is the most favorable trap site, followed by U1 and Si. The less preferred sites 2b and 8j-1 sites were not taken into consideration in our mobility study.  Migration of Nd assisted by U1 vacancies in the a-b plane only is improbable due to the high migration energy, 4.92 eV, which essentially originates from the long jump distance. This jump can instead take place with a lower-barrier combination of two steps: first, Nd in a U1 trap site jumps to its nearest U2 vacancy spontaneously, and then jumps to another U1 vacancy with a barrier of 1.20 eV, as indicated by the blue arrows in Fig. A. 2 (a). The diffusion of Nd assisted by U2 vacancies along the c-axis and in the a-b plane could also occur by similar multi-step mechanisms. However, the migration energies of Nd by Si vacancies in the a-b plane are lower than along the c-axis. Nd interstitial diffusion mechanisms have relatively high barriers. The binding energies (E b ) between a vacancy and a Nd atom at each substitutional trap site are plotted in Fig.4. Nd at U2 and U1 trap sites are strongly binding to all kinds of vacancies. Interactions between Si vacancies and Nd atoms located at Si sites are significantly repulsive, which indicates that vacancy-assisted Nd diffusion on the Si sublattice is likely to play a less important role than that on the U sublattice. The total Nd diffusion coefficients in U 3 Si 2 in the a-b plane and along the c-axis were calculated using Eqs. (21)(22)(23)(24). Fig. 5 shows the total Nd diffusion coefficients and the contribution of each mechanism. Nd diffusion across the two sublattices was not taken into account because the contribution is expected to be negligible, based on the fact that Nd at the Si sublattice is more likely to jump back to the U sublattice instead of continuously migrating in the Si sublattice. In the hightemperature range, Nd atoms diffuse at the same rate along the c-axis and in the a-b plane because the dominant diffusion mechanism, uranium vacancy-assisted diffusion, is isotropic. The possible anisotropy depends on the difference between barriers in the a-b plane, and along the c direction.
The larger the difference, the stronger the anisotropy. For Nd migration on U1 or U2 sublattices, the difference is small, as shown in Table 2 (right). In the a-b plane, the direct interstitial diffusion coefficient in the a-b plane has a much lower activation energy compared with the other diffusion coefficients. As a matter of fact, contrary to the other mechanisms, direct interstitial diffusion of Nd atoms does not require the nearby presence of assisting defects, thus its activation energy includes no defect formation energy. So, the total activation energy for interstitial diffusion (0.24 eV, by fitting the points in Fig. 5) is the effective migration energy resulting from the statistical average between Nd inter-ab-1 (2.67 eV) and Nd inter-ab-2 (0.08 eV). With decreasing temperature, uranium vacancy-assisted diffusion becomes slower than the interstitial diffusion. Below 700 K, Nd atoms diffuse mainly through the interstitial diffusion mechanism, and this is reflected in the sudden change of slope appearing in the Nd total diffusion coefficient in Fig. 5. Silicon vacancy-assisted mechanisms contribute less to the total diffusion of Nd since the Si vacancy-Nd interactions are repulsive, while the U vacancy-Nd ones are attractive. We can then conclude that Nd diffusion occurs primarily on the U sublattice.   For silicon vacancy-assisted diffusion along the c-axis, Nd atoms diffuse more slowly than Si selfdiffusion, so that 0<D pd <1 and Nd enrichment at sinks occurs through the inverse Kirkendall mechanism. The Nd enrichment/depletion tendency will be dominated by diffusion on the U sublattice because diffusion on the Si sublattice is much slower.

Discussion
Once the diffusion coefficients are known, we need a criterion for assessing whether Nd mobility is low enough to be considered as practically immobile. Since Nd has been licensed as a burnup indicator in UO 2 , its mobility can be taken as a reference. Here, we decide to evaluate Nd mobility both in absolute terms and relatively to self-diffusion because the fuel temperature would not be identical between oxide and silicide fuel. According to the technical report of Westinghouse Electric Company, for a peak linear power of 49.9 kW/m, the maximum fuel centerline temperature is 1066 °C for U 3 Si 2 /SiC fuel and 908 °C for U 3 Si 2 /Zr-alloy fuel, against 2058 °C for UO 2 /Zr-alloy fuel [43].
The Nd diffusion and U self-diffusion in U 3 Si 2 and UO 2 are shown in Fig. 7. The U self-diffusion coefficients in UO 2 were experimentally observed by Matzke et. al. [18] and the Nd diffusion coefficient under thermal condition was measured by Han et. al. [44]. The Nd diffusion coefficient in UO 2 at 1070 K is similar to the Nd diffusion coefficient in U 3 Si 2 at 800 K. Assuming that Nd diffusion has similar activation energies in UO 2 and U 3 Si 2 , D Nd in UO 2 at 2058 °C is close to D Nd in U 3 Si 2 at 2328 °C, which is clearly significantly larger than D Nd in U 3 Si 2 at 1066 °C. This makes Nd a promising burnup indicator in U 3 Si 2 concept fuels because it has a slower migration in U 3 Si 2 than in UO 2 at the relevant fuel temperature, and should thus remain in the grain where it first appears.
In addition, if Nd has a similar or slower diffusivity as U, it is reasonable to assume that it evolves with the matrix and with limited probability to make it to the grain boundaries. We compare Nd diffusion with U self-diffusion rather than Si self-diffusion because Nd diffusion on the Si sublattice is significantly slow and can therefore be neglected. Similarly, in UO 2 , there is no interaction between Nd and O sublattices, and Nd diffusion occurs mainly on U sublattice [45]. One can note that the diffusivity of Nd is 7~8 orders of magnitude higher than U in UO 2 at 1070 K, where there is experimental data, but this ratio is expected to be smaller at higher temperature. In U 3 Si 2 , Nd diffusion is significantly slower than U in the a-b plane in the whole temperature range. Along the c-axis, Nd diffusion is slightly faster than U and the ratio is about 7 at 1300 K, which is significantly lower than the ratio in UO 2 . So, if Nd mobility is low enough in UO 2 , it is certainly low enough in U 3 Si 2 . However, the large gap between the diffusivity of Nd and U in UO 2 implies that low diffusivity is not the only possible condition to ensure that Nd is immobilized in the fuel matrix. Analysis of other microstructure phenomena involving Nd (for instance, clustering and formation of precipitates, interaction with dislocations) and the calculation of Nd diffusivity in UO 2 must be done in order to draw a stronger conclusion in this respect. The Xe bulk diffusion coefficients in U 3 Si 2 calculated by Andersson et al. [16,17] are illustrated in Fig. 7 for comparison. One can see that Xe diffusion coefficients are close to our calculated Nd diffusion coefficients. In Ref. [16,17], the Xe diffusion coefficient was calculated using a ratelimiting model that does not take into account kinetic correlations between solutes and defects. With this model, the total diffusion coefficient along a specific direction can be only affected by the migration mechanisms along this direction, while in the framework applied in this work, migration mechanisms in all directions contribute to the total diffusion coefficient. In order to compare the two models, we computed the diffusion coefficient of Xe in Si vacancies using the KineCluE code, selecting the required DFT results from [16,17], and compared the results to those by Andersson et al. in Fig. 8. The rate-limiting model underestimates the diffusion coefficient by 2 to 3 orders of magnitude. Taking this discrepancy into account, Xe diffuses slightly faster than Nd in U 3 Si 2 .
However, this similarity between the diffusion coefficients of Xe and Nd in U 3 Si 2 does not necessarily entail that the two species have similar mobilities overall because, as previously mentioned, bulk diffusion is not the only phenomenon at play. It is known, for instance, that diffusivity in microstructural defects such as dislocations, grain boundaries or cracks provides the main contribution to fission gas release: in UO 2 , the pipe diffusion coefficients of fission gases, such as Xe, along dislocations are 10 13 ~ 10 15 times higher than the bulk diffusion coefficients [46].
An analogous mechanism is not known for Nd in UO 2 . It is thus reasonable to assume that a similar argument would apply to U 3 Si 2 , in which case the similar bulk diffusivities of Xe and Nd would not be relevant to judge whether Nd can be regarded as a slow diffuser in U 3 Si 2 . Nd could also precipitate in the matrix and therefore immobilize, or be trapped by other impurities, but these two effects are not considered in this work. Fig. 8 The diffusion coefficients of Xe in Si vacancies computed with the rate-limiting model (as reported in Ref. [16,17]) and the SCMF theory (computed with the KineCluE code), both based on the DFT results of Ref [16,17].

Conclusion
Density functional theory (DFT) calculations have been used to investigate Nd stability and diffusion in U 3 Si 2 , which is one of the main accident tolerant fuel (ATF) candidates. Nd is predicted to preferentially occupy the uranium sublattice. The thermal stability of vacancy-Nd pairs was estimated by calculating the binding energies. Our results show attraction for uranium vacancies with Nd on the U sublattice, and repulsion for silicon vacancies with Nd on the Si sublattice.
Migration energies were then computed between the various configurations. Migration through alternated jumps on U1 and U2 sites was proven to have the highest probability. Diffusion coefficients of Nd in U 3 Si 2 were then obtained within the SCMF framework. The diffusivity of Nd showed only a small anisotropic character. The fastest diffusion mechanism for Nd is the uranium vacancy-assisted diffusion, which is isotropic. On the other hand, direct interstitial diffusion in U 3 Si 2 shows strong anisotropy.
For vacancy-assisted diffusion, the flux-coupling characters of the vacancy-Nd clusters were analyzed based on the transport coefficients obtained with the KineCluE code. In the temperature range of interest, no vacancy drag occurs in any cases, which means that Nd and vacancies migrate in opposite directions. The computed PDC ratios imply that, for silicon vacancy-assisted diffusion in the a-b plane and uranium vacancy-assisted diffusion in both directions, Nd depletion at sinks is expected. This depletion tendency increases the probability for Nd to stay within the matrix. In any case, any radiation-induced segregation/depletion tendency should occur too slowly to be observed, given the low mobility of Nd.
Nd diffuses slower than U in U 3 Si 2 in the a-b plane and slightly faster than U along the c-axis.
Comparing to the situation in UO 2 that Nd diffusion is significantly faster than U self-diffusion, Nd can be regarded as a slow diffuser. At relevant maximum centerline temperatures of the fuel, the Nd diffusion rate in U 3 Si 2 is significantly smaller than that in UO 2 . Since the mobility of Nd is low enough to satisfy the requirement of burnup indicator for UO 2 fuel, Nd with even lower mobility in U 3 Si 2 can be regarded as a strong candidate burnup indicator for U 3 Si 2 concept fuel. Further studies on its behavior in the irradiated microstructure are needed to confirm its suitability as local burn-up indicator.  Fig. A. 1 (a-c). Fig. A. 1 (d) shows that U interstitial pushes its neighboring U2 atom away from the equilibrium lattice site. It migrates by taking the U2 position and spontaneously kicking the original U2 atom to an interstitial site, which is named I U-ab-1 and is marked with the red arrows. The green arrow is the migration of interstitial U to its neighboring equivalent position, named I U-ab-2 . The U interstitial migration along the c-axis (I U-c ) also involves multiple steps. As illustrated in Fig Fig. A. 2 (a). Nd jumps to a U1 vacancy along the c-axis and in the a-b plane, respectively named Nd U1U1-c and Nd U1U1-ab , are marked with the red arrows. Nd exchanging positions with a U 2 or Si vacancy is illustrated with the blue or green arrows and named as Nd U1U2 or Nd U1Si . Note that Nd U1U2 (or Nd U1Si ) and its reverse mechanism Nd U2U1 (or Nd SiU1 ) are non-equivalent, see Table 2 for their migration barriers. The long distance migration Nd U1U1-ab requiring a very high migration energy can be replaced by a multi-step migration mechanism, by which Nd pass through a U2 site (blue arrows) and then jump to the destination U1 site (blue dashed arrows). Fig. A. 2 (b) shows the migration mechanisms of Nd at U2 sublattice. It can migrate to a U2 vacancy along the c-axis (red dashed arrows) or in the a-b plane (red arrows).

Appendix B. Reference binding energy for configurations in different sublattices
For the mono-vacancy calculation in the U sublattice, there are two non-equivalent configurations, depending on whether the vacancy is in the U1 or U2 sublattice. In order to take into account the difference in formation energy, we need to provide a suitable energy for each configuration that correctly describes this difference. This is done by modifying the binding energies in the KineCluE code, by using the most stable configuration as a reference. In our case, a U1 vacancy having lower formation energy is more stable than a U2 vacancy. Therefore, we add a binding energy of -|E vacU1.f -E vacU2.f | to the U2 vacancy and leave the binding energy of the U1 vacancy at zero (taken therefore as the reference).
For the vacancy-solute pair calculation, in addition to the difference in formation energy, we need to take the difference of solution energy into account. As it is reported in Table 3, the most stable substitutional site of Nd is on the U1 sublattice. Therefore, when computing the binding energies of vacancy-solute pairs with Eq. (17), terms and are fixed to the total energies of the supercells with a single U1 vacancy and with a single Nd on U1 sublattice, respectively.
No modifications are needed for Si self-diffusion and Nd diffusion in Si sublattice because there is only one type of Si sublattice.

Appendix C. Defect formation entropy and vibration frequencies
The defect formation entropy can be obtained by S D and S 0 are the vibrational entropy of the supercell with and without defect respectively. The vibrational entropies were calculated according to the approach of Mishin et al. [22], which approximates the entropy of crystalline solids at temperatures higher than the Debye temperature as = − ∑ ln( ℎ ) 3 −3 =1

+ (3 − 3)
A. (2) where ν n are the normal vibrational frequencies of the crystal, N the number of atoms in the investigated supercell. Eq. (4) is the results of combining the equations above.
Vibrational frequencies of the perfect U 3 Si 2 supercell 0 , and of the supercell with point defects (D = VacU1, VacU2, VacSi, InterU, and InterSi) calculated using the DFPT method are provided in

Declaration of interests
☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.