Pipe and grain boundary diffusion of He in UO2

Molecular dynamics simulations have been conducted to study the effects of dislocations and grain boundaries on He diffusion in UO2. Calculations were carried out for the {1 0 0}, {1 1 0} and {1 1 1} 〈1 1 0〉 edge dislocations, the screw 〈1 1 0〉 dislocation and Σ5, Σ13, Σ19 and Σ25 tilt grain boundaries. He diffusivity as a function of distance from the dislocation core and grain boundaries was investigated for the temperature range 2300–3000 K. An enhancement in diffusivity was predicted within 20 Å of the dislocations or grain boundaries. Further investigation showed that He diffusion in the edge dislocations follows anisotropic behaviour along the dislocation core, suggesting that pipe diffusion occurs. An Arrhenius plot of He diffusivity against the inverse of temperature was also presented and the activation energy calculated for each structure, as a function of distance from the dislocation or grain boundary.


Introduction
An important measure of nuclear fuel performance is the propensity to retain inert gasses produced during irradiation. The gas inventory is largely comprised of Xe and Kr, which are fission products, while He can also be formed via fission. However, He is far more likely to occur through the α-decay of minor actinides which results in it being very important for storage and disposition of spent nuclear fuel, where a large quantity of He will be produced over both short and very long time scales [1]. Furthermore, due to its high thermal conductivity, He is used as a filling gas between the fuel rods and cladding materials, making He the most abundant inert gas species in the fuel rod, if not in the pellet itself. Many studies principally focus on Xe and Kr [2][3][4] as He is produced in smaller quantities within conventional UO 2 pellets during normal operation. However, unlike Xe and Kr, He has an atomic radius that is small compared to the lattice spacing in UO 2 leading to different solubility and diffusivity behaviour than Xe and Kr [5,6]. He tends to form bubbles or diffuse into grain boundaries resulting in fuel swelling [7,8], thereby degrading fuel performance and properties. Moreover, fuel swelling could increase the internal pressure on the cladding that surrounds the fuel rod.
It has been reported that He diffusion can be influenced by defects in UO 2 [9,10]. Grimes et al [10] showed that at vacancy and interstitial sites the solution energy for He is small (∼ −0.1 eV/He atom), while in the perfect lattice on interstitial sites the barrier to diffusion is large (3.8 eV). In radiation environments this diffusion barrier can become lower (∼0.3 eV) through the creation of defects such as cation and anion vacancies. Govers et al [9] took the study of He diffusion in UO 2 further, predicting assisted He diffusion by extrinsic oxygen vacancies for low temperatures. At higher temperatures, or if no extrinsic oxygen vacancies are present, the diffusion mechanism occurs by interstitial hops [9]. Moreover, for low temperatures a decrease in activation energy for He in hypo-stoichiometric fuel is almost proportional to the number of vacancies [11]. These papers provide evidence that He could follow a vacancy pathway, if available in the structure.
In UO 2 , it has been predicted that self diffusion is enhanced by the presence of a dislocation, suggesting pipe diffusion occurs [12]. Also in UO 2 , the activation energy for isolated Xe atoms near dislocations is low relative to the bulk [13]. However, as the isolated Xe atoms diffuse along the core they begin to cluster together, actually blocking the core and inhibiting further pipe diffusion.
Here we investigate the influence of dislocations and special grain boundaries on the transport of He in UO 2 based on a combination of pair and many-body interactions and building upon previous work carried out by Murphy et al [14] and Parfitt et al [15] that investigated the effects of pair potentials on the core structure of dislocations and the methods to perform atomistic simulations on dislocations.

Potential model
Molecular dynamics (MD) simulations, employing a set of interatomic potentials derived previously by Cooper, Rushton and Grimes (CRG) [16][17][18] 3 , are carried out using  the large-scale atomic/molecular massively parallel simulator (LAMMPS) [20]. In this model the potential energy, E i , of an atom i with respect to all other atoms has two components-(a) a pair potential description and (b) a many-body embedded atom method (EAM) contribution using the model of Daw and Baskes [22].
Although He parameters have not been explicitly developed for this UO 2 potential, it has been shown previously [9] that the He-He, He-O and He-U parameters derived by Grimes et al [10] are transferable to other potentials. This was tested by comparing He incorporation energies at interstitial, uranium and oxygen sites in UO 2 . Hence, to validate the combination of the Grimes He potential [10] with the CRG UO 2 [16] potential, incorporation energies were calculated and compared to the literature data using other UO 2 potentials (see section 3.2). For comparison bulk He diffusion was also tested using the Buckingham-Morse description for UO 2 of Basak et al [23].

Dislocation creation
The dislocations considered in this study are straight lines, with a Burgers vector of a 2 1 1 0 / ⟨ ⟩. To avoid interface effects caused by free surfaces (such as surface polarisation), full periodic boundary conditions were used, which impose that the sum of the Burgers vectors for all the dislocations in the supercell must be zero. The supercells used contained dislocation dipoles, with the Burgers vector of the first dislocation opposite to that of the second one. The supercells were about 230 230 130 Å × × (∼500 000 atoms). The dislocations lines were along the z axis with the slip planes for the edge dislocations perpendicular to the y axis. The distance between the two dislocations was about 115 Å.
The creation process for these dipoles is analogous to the Volterra construction [43], and has been used previously for work on dislocations in UO 2 [45]. To create a dipole with the dislocations having a Burgers vector of respectively b and b − , the crystal is cut along a half-plane ending at the dislocation lines. The displacement field caused by the dislocations is calculated in each point of the supercell from anisotropic elasticity theory. This analytic displacement field is then applied to the atoms. Close to the cut plane, the atoms either side are displaced by b 2 / and b 2 / − respectively. For the case of edge dislocations, this causes both sides to overlap, bringing atoms very close to each other. To avoid this, the atoms whose initial position is within b 2 / | | from the cut plane are removed from the configuration.
The displacement field includes a non-periodic component along the direction of the dislocation line, which causes internal stresses. To accommodate those, the shape of the supercell is slightly adjusted, following the procedure published by Cai et al [41,42]. This procedure has been implemented in the Babel code [44], which has been used in this study.
The displacement field from the elasticity theory does not describe appropriately the displacement of the atoms very close to the dislocations cores. Thus, even after adjusting the supercell, some of these atoms can be in unstable positions, and the configurations are not immediately suitable for MD simulations.
Following an initial conjugant gradient energy minimisation during which atomic positions and supercell parameters were relaxed, the system was equilibrated at 1500 K and 0 GPa for 200 ps in the NPT ensemble using Nosé-Hoover thermostat and barostat relaxation times of 0.1 ps and 0.5 ps respectively. The barostat is applied to all supercell parameters independently to allow for asymmetric relaxation. To calculate the dislocation line energy the system temperature is reduced to 300 K over 80 ps and then energy minimised again using the conjugant gradient method. The line energy is the difference in energy between this relaxed cell and the energy of an equivalent number of UO 2 formula units in the perfect cell, per unit supercell length in the direction of the dislocation line.
The resulting core structure are shown in figure 1, and are similar to previously reported structures [14,45]. It must be noted these core structures have not yet been observed experimentally and the influence of the potential on the core structure cannot be disregarded, although Murphy et al [14] showed that there is good agreement between UO 2 potentials on the core structure. The tilt axis is along z, the grain boundary plane indicated by the dashed lines is perpendicular to x and the misorientation angle is Θ 2 .

Grain boundaries creation
A grain boundary is a surface defect that separates two crystalline grains with different orientations. In this study, symmetrical tilt grain boundaries were generated from perfect UO 2 unit cells, and their properties are summarised in table 1. They represent geometrically simple configurations and have been investigated using the Morelon potential [26,46].
In the case of symmetric tilt boundaries, the crystalline structure from one grain can be obtained from that of the other grain by rotating it with a rotation axis along the grain boundary plane. The procedure to generate them follows this description.
First, the unit cell is rotated by an angle Θ equal to half the misorientation angle (see figure 2). The rotated cell is replicated along the x and y directions, and a supercell is created from this tiling, making sure that the supercell axes are along x, y and z respectively, and are periodic vectors of the rotated structure. A grain is built from this supercell by replicating it along the three directions to provide adequate system size and grain size. It is then combined with its image in a reflection on a (y, z) plane to form an unrelaxed simulation box. Because of the periodic boundary conditions, each simulation box contains two infinite grain boundaries in (y, z) planes, separated by a distance equal to half the box size along the x direction.
Exactly the same heat-quench relaxation method is applied to the grain boundary structures as was applied to the dislocation structures (see section 2.2). This ensures that if the initial grain boundary structure is in a metastable minimum it can reconfigure to the global minimum.

Calculation details
Defect energies were determined by energy minimisation using a supercell of 10 10 10 × × fluorite unit cells, implemented within LAMMPS. This has been shown previously to be of sufficient size to satisfy the criterion for the dilute limit when using the CRG potential [17]. Firstly, the perfect UO 2 structure was minimised using a conjugate gradient method at 0 GPa, ensuring fully relaxed lattice parameters. Subsequently, the defect of interest was introduced to the structure and the atomic positions are relaxed at constant volume using a steepest decent procedure [24]. The defect energy was given by the difference between the defective cell and the perfect cell: Incorporation energies are defined as the defect energy of He substituted at an incorporation site minus the formation energy for that site. Therefore, as the incorporation site for an interstitial is simply the perfect lattice, it is equivalent to the defect energy of a He atom occupying an interstitial site. Molecular dynamics investigations of He diffusivity in bulk UO 2 were carried out using a supercell, as described in section 2.2, with He randomly distributed at 4d Wyckoff interstitial sites. Two concentrations of 0.83 and 4.13 at.% He were To calculate diffusivity as a function of distance from the dislocation core the supercell is divided into two sets of concentric shells centred on each dislocation. For the grain boundary the supercell is split into a set of slabs surrounding the boundary. The colour coding represents the local energy density.  tested using both the Basak and CRG potentials to describe the UO 2 interactions.
Initially the system was energy minimised using the conjugate gradient method followed by 40 ps of equilibration in the isothermal-isobaric (NPT) ensemble at 0 GPa and at the desired temperature. For all MD simulations Nosé-Hoover barostat and thermostat relaxation times of 0.5 and 0.1 ps respectively are used. A timestep of 2 fs is used throughout. Subsequently, the He mean squared displacement (MSD), R He 2 ⟨ ⟩, is calculated and the diffusivity, D, is determined using equation (2), where t and d are the time and the number of dimensions over which diffusion is calculated respectively. The bulk diffusivity values reported in figure 4 were averaged over ten simulations, seeded with random velocities and run for 1 ns each. By plotting the log of diffusivity, D ln , as a function of 1/T and assuming an Arrhenius plot the activation energy for He migration (E a ) can be obtained from the gradient of the graph: where D 0 is the pre-exponential, k B is the Boltzmann constant and T is the temperature.
To calculate diffusivity as a function of distance from the dislocation core, the supercell is divided into two sets of concentric cylindrical shells centered on each dislocation (see figure 3). Each set contains 10 shells that are 5 Å thick, meaning that they cover a cylindrical region between 0 Å and 50 Å from the dislocation centre. The symmetry of the dislocations means that the two sets of shells can be merged into a single set of 10. For grain boundaries, rather than dividing the supercell into sets of concentric shells the supercell was split into a set of slabs around the two grain boundaries. Similarly, each set contains 10 slabs 5 Å thick, meaning that they cover a region between 0 Å and 50 Å from the grain boundary. The supercells were then equilibrated for 35 ps at the target temperature (ranging from 2300 K to 3000 K) before the He MSD was calculated in each shell for 2 ps (the first 200 fs representing the balistic phase are omitted). This short time for calculating the diffusivity is selected intentionally to ensure that He atoms do not diffuse between shells during the simulation. Therefore, to obtain sufficient statistical significance the calculation was repeated 25 times for each temperature with each dislocation structure. Not only was the total MSD in each region calculated but also the individual Cartesian He, 2 ⟨ ⟩. As such, the ratio of diffusivity in line with the dislocation, D z , to the diffusivity in the perpendicular plane, D xy , can be determined as a measure of the anisotropy of the total diffusion. Table 2 shows the predicted line energies calculated over the entire simulation cell. It suggests that the most stable edge dislocation is the {1 0 0} system followed by the {1 1 1} system with an energy difference of ∼0.24 eV Å 1 − . This agrees well with the calculations, using the Morelon potential [47] carried out by Parfitt et al [15], where dislocation {1 0 0} was found to have an energy of 0.25 eV Å 1 − greater than dislocation {1 1 1}. Parfitt's paper [15] also notes that this is consistent with the prediction made by Keller et al [25] in that the dominant slip system in UO 2 is {1 0 0} 1 1 0 ⟨ ⟩ with the {1 1 1} 1 1 0 ⟨ ⟩ system existing as a secondary system. Atomic simulations were also carried out on dislocations in UO 2 by Murphy et al using many different potentials [12,14]. The ordering of the most stable dislocation system in this paper agree with both Parfitt [15] and Murphy [12,14] (The screw 1

Line and grain boundary energies
The energy of each grain boundary is calculated over the simulation cell and shown in table 3. The most stable grain boundary is the Σ5. The Σ5 symmetrical is one of the most studied grain boundaries both experimentally and theoretically  which has several studies on fluorite-like structures [26][27][28][29][30][31][32][33][34]. Grain boundary energies reported in this paper for Σ13 and Σ25 are lower than those by Van Brutzel et al [26] and the Σ5 higher (this is most likely due to the fact their grain boundaries were calculated in a different size simulation cell), however, the trend is the same. Experimental and theoretical studies carried out by Nerikar et al [30] investigated grain boundary structures in UO 2 . Using the Basak potential [23], the Σ5 tilt grain boundary was found to have a grain boundary energy of 1.58 J m 2 − , a value slightly higher than the energy predicted in this paper (see table 3).

Incorporation energy
The energy for He incorporation into a uranium vacancy, oxygen vacancy and two possible interstitial sites as predicted using different UO 2 potentials is shown in table 4. For all calcul ations the He-He, He-O and He-U interactions of Grimes et al [10] were used. The extent to which the incorporations energies agree with the Grimes UO 2 potential, for which the He interactions were originally developed, gives an indication of the transferability of the He potential. The Basak potential and CRG potential both exhibit excellent agreement with the values of Grimes et al [10] for all incorporation sites considered. However, the Morelon potential yields markedly different results for the 24d interstitial site. As such, the Basak and CRG potential have been selected for further comparison in MD for bulk He diffusivity at two concentrations: 0.83 at.% and 4.13 at.%.

Bulk diffusivity
The bulk He diffusivity is reported in figure 4 using both the Basak and the CRG descriptions of the UO 2 interactions. For both potential sets the trends are similar, in that the gradient becomes flatter at high temperatures. This indicates a change in the activation energy and is not dissimilar to oxygen diffusivity, so may be linked to the change undertaken by the lattice during the superionic transition for reasons similar to those given in previous work [17,18,35,36]. Nonetheless, this result is an important reminder to focus consideration on temperatures below the superionic transition (<2700 K) to identify behaviour that may be extrapolated to lower temperatures that are more relevant to fuel either in reactor or in storage. Additionally, both potentials show very little difference between the two concentrations. The higher concentration of He has a slightly lower diffusivity and is outside the bounds of the error bars at lower temperatures, which can be attributed to He clustering. Furthermore, as the CRG potential describes a wider range of thermophysical properties than the Basak potential [16] and so that the results are applicable to lower He concentrations, the remainder of this work will focus on enhanced pipe and GB diffusion using the CRG potential with 0.83 at.% He between 2300 K and 3000 K.   Figure 5 indicates that all dislocations studied exhibit enhanced diffusion in regions closer to the dislocation core. The enhancement is clearest for lower temperatures-for 2300 K the enhancement is just under 2 orders of magnitude for the three edge dislocations and about 1 order of magnitude for the screw dislocations. However, at 3000 K there is no clear enhanced diffusivity for any dislocations studied. If the trend of greater enhancement at lower temperatures is continued there should be an even greater effect at temperatures relevant to spent nuclear fuel (SNF) storage. Moreover, radiation damage will enhance diffusion by creating dislocations. However, further studies are required to rule out He clustering at lower temperatures which may act to limit diffusion. This is particularly important for dislocations pinned to fission product precipitates as this may enhance the rate of growth for these bubbles. Experimental evidence for precipitate-bubble association is given by Baker [37] and Turnbull [38]. Figure 6 depicts diffusivity using a log scale as a function of 1/T for each region from 2.5 Å to 47.5 Å from the dislocation core. Again enhanced diffusivity in regions closest to the dislocation core is predicted. The spread in the diffusivity at lower temperatures also demonstrates that this behaviour may extend and be further enhanced at lower temperatures. Helium activation energy as a function   of distance from the dislocation core (see figure 7) is calculated by taking the gradient of figure 6 between 2300 K and 2600 K (below the superionic transition temperature). Close to the dislocation core the activation energy is quite low and gradually increases further from the core for all dislocations studied here. When the distance is far enough away from the influence of dislocation core it tends towards the bulk value (however there is scatter around the bulk value due to the other dislocation in the system causing a small perturbation at large distances due to elastic strains, see section 2.4). For comparison to the bulk activation energy (see figures 7 and 11), using the nudged elastic band (NEB) technique [39], the energy barrier for a He hop between two interstitial sites in bulk UO 2 was calculated. The energy barrier for this process was calculated at 3.96 eV (which is consistent with previous calculations conducted by Grimes et al [10]). The transport of He via interstitial sites would occur through a non-defective lattice. However, it should be noted that at higher temperatures (where the lattice can contain defects) He migration may be vacancy assisted. According to Grimes et al [10] the energy barriers associated with He migration via oxygen and uranium vacancies is much lower than in the non-defective lattice. There is a slight discrepancy between the energy barrier calculated via NEB and the activation energy for He diffusion in bulk UO 2 (figures 7 and 11). While the NEB technique calculates the energy needed for a He atom to hop between two interstitial sites, the activation energy is calculated from the Arrhenius plot ( figure 3). This entails not only the hopping energy but the energy needed to create defects. Although only one separation between the dislocations has been chosen for this study, the effect of separation is not thought to be important unless the dislocations are much closer together. There is little effect due to the dislocations after 2 nm, considering the activation energy tends to bulk values at this distance. The low activation energy provides evidence for enhanced diffusion in the core. Figure 8 shows the ratio of diffusivity in line with the dislocation core (D z ) against the diffusivity in the perpendicular plane (D xy ) as a function of distance from the dislocation core. This indicates that diffusivity is greatest along the dislocation core for the {1 0 0} edge dislocation, which is also the most stable dislocation system (see section 3.1). There is a slight enhancement in diffusion along the core for the {1 1 0} and {1 1 1} edge dislocations and none at all for the screw dislocation. As there is enhanced diffusivity for all the dislocations (see figure 5) the discrepancy between the extent of anisotropic behaviour for each dislocation could be due to the difference of He diffusing in the directions of the tensile region caused by the dislocation compared to the proportion of He diffusing in the core direction. If this is the case, it gives insight into the effect of the It is also interesting to note, there is virtually no void-like core structure in the screw and so there is also no anisotropy. Conversely, the 3 edge dislocations exhibit a more open core structure (to a greater or lesser extent). The {1 0 0} edge dislocation which exhibits the largest anisotropy has the largest core region (see figures 1(a) and 8), whereas the {1 1 0} edge dislocation which has the smallest core region exhibits limited anisotropy (see figures 1(b) and 8). Although we have predicted a correlation between core structure and aniso tropy, future work should focus on examining the exact mech anism of diffusion in the core and the effect on anisotropic diffusion. The anisotropic behaviour properly channels the He diffusion and will lead to long distance diffusion in UO 2 , whereas isotropic enhancement will only create regions of enhanced diffusion.

Grain boundary diffusion
3.5.1. Diffusivity versus distance. Similar to section 3.4, He diffusivity as a function of distance from the grain boundary is calculated for the Σ5, Σ13, Σ19 and Σ25 tilt grain boundaries (depicted in figure 9). Again, enhanced diffusivity is observed for each grain boundary structure particularly at lower temperatures. Unlike some of the dislocations, isotropic diffusion is predicted for all grain boundaries and anisotropy will not be further discussed. Figures 10 and 11 report diffusivity using a log scale as a function of 1/T for each grain boundary for every region from 2.5 Å to 47.5 Å and a plot of He activation energy as a function of distance from the grain boundary. This result is similar to that attained for the dislocations (see figure 6) in that there is an enhancement of diffusivity of two orders of magnitude. As with the dislocations, the activation energy (figure 11) is calculated by taking the gradient of figure 10 between 2300 K and 2600 K. Close to the grain boundary the activation energy is low and gradually grows as the distance from the grain boundary is increased until it reaches the bulk value. This again is similar to the dislocation activation energy where there is a deviation in activation energy of ∼4 eV at the grain boundary compared to the bulk and that the deviation begins at 20 Å. As such, there is an increase in diffusivity (visualised in figure 9) close to the grain boundary (similar in terms of magnitude to that of the dislocation, seen in figure 5). Ronchi et al [40] carried out experiments on He diffusion in UO 2 and concluded that 'atomic diffusion is controlling short range migration to still unidentified traps from which the gas subsequently migrates and escapes with a very low activation enthalpy' and that the sink strength of the traps depends on the level of lattice damage. It is possible that grain boundary dislocation interactions can be a release mechanism for He, so that He could diffuse through the bulk quickly via dislocations (created by radiation damage) to grain boundaries and then from grain boundaries  [40] also suggested that migration on dislocations and grain boundaries is the possible cause of the low activation enthalpy associated with one of the He release stages.

Conclusions
In this work, MD calculations are performed to predict the influence on He diffusion in UO 2 of dislocations ({1 0 0}, {1 1 0}, {1 1 1} edge dislocations and screw) and grain boundaries (Σ5, Σ13, Σ19 and Σ25 tilt) over the temperature range 2300-3000 K. It is clear that dislocations and grain boundaries enhance the diffusivity of the He atoms in UO 2 by up to two orders of magnitude at 2300 K. There is also evidence for pipe diffusion of He in dislocations. However, not all dislocations are equal in their influence on He transport in the z direction. In particular, the {1 0 0} edge dislocation exhibits enhanced diffusion by a factor of 2 in the z direction compared to the xy plane (see figure 8). This is not the case for the {1 1 0} and {1 1 1} edge dislocations and the screw dislocation. The activation energy as a function of distance is also investigated for each dislocation and grain boundary. It is seen that for each dislocation and grain boundary the activation energy decreases by ∼4 eV from the bulk value as the distance to the dislocation or grain boundary decreases. As mentioned in section 3.5.2, grain boundary dislocation interaction can be a release mechanism for He from UO 2 . This can be problematic as the release of He can influence fuel performance and cause pressurisation problems between the fuel and cladding.  Sigma 19 Sigma 25 Figure 11. The activation energy for helium diffusion as a function of distance from the Σ5, Σ13, Σ19 and Σ25 tilt grain boundaries. The activation energy was calculated from 2300 K to 2600 K using figure 10. Activation energy (eV) Distance from grain boundary (Å) Bulk Moreover, with respect to waste management, a change in transport of this magnitude should be taken into account when considering time scales over which containers may be subject to increased He pressurisation.