Radiation damage effects on helium diffusion in zircon

We report the effects of radiation damage on helium diffusion in zircon using data from molecular dynamics simulations. We observe an increase in activation energy for helium diffusion as a result of radiation damage and increasing structural disorder. The activation energy in a heavily damaged region is smaller than in a completely amorphous system which is correlated with remaining order in the cation sublattices of the damaged structure not present in the fully amorphized system. The reduction of activation energy is related to the disappearance of fast diffusion pathways that present in the crystal. Consistent with the change in activation energy, we observe the accumulation of helium atoms in the damaged structure and discuss the implications of this effect for the formation of helium bubbles and zircon's performance as an encapsulation material for nuclear waste.

We report the effects of radiation damage on helium diffusion in zircon using data from molecular dynamics simulations. We observe an increase in activation energy for helium diffusion as a result of radiation damage and increasing structural disorder. The activation energy in a heavily damaged region is smaller than in a completely amorphous system which is correlated with remaining order in the cation sublattices of the damaged structure not present in the fully amorphized system. The reduction of activation energy is related to the disappearance of fast diffusion pathways that present in the crystal. Consistent with the change in activation energy, we observe the accumulation of helium atoms in the damaged structure and discuss the implications of this effect for the formation of helium bubbles and zircon's performance as an encapsulation material for nuclear waste.

I. INTRODUCTION
The effects of radiation damage on zircon are one of the most well studied and characterized of any material. For many years it has been considered a prime candidate for the use as an encapsulation matrix for radioactive waste 1 . The initial interest was excited when natural zircon samples were found to contain low concentrations of radioactive actinides and as such was considered a 'natural encapsulation matrix'. This has resulted in many years of experimental studies on natural zircon samples examining the direct impact of radiation damage via ion beam irradiation 2-5 . These experiments observed that irradiation results in the amorphization of zircon, which is accompanied by a large volume expansion of approximately 18-20 % 6 . In more recent years the use of Molecular Dynamics (MD) simulations have been carried out to aid the understanding of radiation damage in zircon, predicting the formation of Si-O bridges, the percolation of damage, and density inhomogenities formed as a direct result of simulated collision cascades [7][8][9] .
Crucial for the development of future safe wasteforms is determining the effect of helium buildup produced as a result of the alpha decays of the radioactive actinides present. Helium buildup in wasteforms has been linked to inducing embrittlement of the wasteform, swelling from helium bubble formation and potential pressurization of metal canisters enclosing the waste form 10,11 . Decreased brittleness and microcracking as a result of swelling can lead to increased dissolution rates and influence the movement of actinides 11 . These concerns are important due to helium bubbles being identified in different highly damaged natural minerals [12][13][14] .
Another result of the radioactive elements present in natural zircon is that it is often used as a thermochronometer, due to zircon (U-Th)/He (ZHe) decay 15,16 . Thermochronometery gives us information on the history of the sample, its age and most importantly the rate at which the specific sample has been cooled. However for accurate ages to be calculated, the mechanisms of helium loss need to be well understood, as the age is calculated by the ratio of uranium and thorium to helium and lead. The alpha-decay of actinides in zircon introduces the helium (at high energies) whilst simultaneously changing the local structural environment effecting the rate of helium diffusion as a result of radiation damage. The destruction of crystallographic directions due to even low levels of radiation damage can decrease diffusivity by destroying fast diffusion pathways 17,18 . This is due to the highly anisotropic nature of He diffusion in zircon. The most favourable pathway for He diffusion is along the [001] direction parallel to the c axis 15,16,19,20 . Even small changes in structure as the substitution of rare earth elements instead of Zr and Si produce differences in diffusion rates between natural samples and synthetic analogues 18 . The production of voids induced by radiation damage 21 can result in regions of up to 3 times greater concentration of helium than the surrounding undamaged crystal. This suggests that the rate of diffusion in a damaged, semi-amorphised region of zircon will be lower than that in a crystalline region.
Similar results are found when comparing the effect of radiation dose on diffusion in different samples 19 . Dose is often used to quantify the amount of 'damage' in a sample and is measured as the number of alpha-decays/g. As dose increases up to 1 × 10 16 α/g the activation energy needed for diffusion also increases. However, further increasing the dose to 1 × 10 19 α/g results in a decrease in the activation energy 19 . This indicates that the initial damage to the sample decreases diffusivity until a critical damage point is reached, after which further damage only increases diffusivity. This has been linked to the effect of percolation of damage creating larger regions of overlapped 22 damage, allowing easier He diffusion to take place along these pathways. In minimally damaged natural zircon samples, with doses of 1 × 10 15 α/g, the mean distance between damaged regions can be as little as 5 µm. This means all fast diffusion pathways are interrupted, or completely damaged, thus resulting in a regime where He atoms diffuse alternately between isolated damaged and crystalline regions 23 . Taking this into account it can be difficult to compare diffusion data between different samples, as they all have different ages, doses and thermal histories.
Various computational modelling techniques have been used to study helium diffusion in zircon. Common practise is to use density functional theory (DFT) and in particular nudged elastic band to calculate the activation energy for helium diffusion from one minima to another. 17,22,24,25 . The values for the activation energy vary depending on the parameters used, but in general they agree with each other and experimental work that diffusion parallel to c is more favourable than any other direction. It is also worth noting that calculations of activation energies for helium diffusion are much smaller than experimental studies, due to simulations being conducted on perfect crystals with no elemental impurities, defects or radiation damage which can cut off the c-channel along which He prefers to diffuse. The only computational work we know taking into account the effect of radiation damage on helium diffusion calculates diffusion based on the hopping of helium atoms between interstitial sites 22 , and models the damage by blocking some percentage of the sites depending on the degree of damage.
The aim of this work is to use molecular dynamics (MD) to calculate activation energies for helium diffusion in three different zircon structures; a perfect crystal, a homogeneous amorphous structure and a damaged crystalline structure produced as a result of 6 overlapping collision cascades of a 70 keV U ion. This is the first time such calculations have be performed using large scale molecular dynamics. We will present Arennhius plots of the different structures at a range of temperatures and show that radiation damage results in an increase in activation energy for helium diffusion. We will also show that even a small amount of retained order can lead to a much faster rate of diffusion than a completely amorphous structure. We also show that there is a greater build up of helium in damaged regions over time, and discuss the implications of this effect for helium bubbles and waste form performance.

II. METHODS
The molecular dynamics (MD) package DL_POLY_4 (version 4.10 with changes released in version 5.0.0) 26 has been used to model He diffusion in zircon systems, due to its ability to simulate large systems with a large number of atoms, in parallel, by using domain decomposition.
The interatomic potentials used in this work to model zircon have been described in previous papers 7, 9 . The Zr-O and O-O pairwise potentials are in the form of a Buckingham potential, with the Si-O pairwise potential taking the form of a Morse potential. The interaction between He and the constituent atoms in zircon must also be included, and are taken from previous work 24 . We also include a term to model the He-He interaction which has been well used and tested in MD simulations of uranium dioxide 27,28 .
This work uses the radiation damage structures produced in our previous paper 9 . All the systems modelled contain 10 million atoms. There are 3 systems used: a pure crystal structure, a melt-quench amorphous structure produced using a cooling rate of 10 K/ps 9 , and crystalline structure with radiation-damaged regions. The radiation damaged structure is produced as a result of 6 overlapping radiation cascades, each one due to a U ion with 70 keV of kinetic energy with trajectories along the same initial direction from the same initial position. The irregular shape of the damaged region can be seen in Fig. 1, with an estimated volume of 100,000 Å 3 . The damaged region was determined empirically as a parallelepiped to include the full damaged region, such that there is a clear boundary visible between the crystal and the damaged region. The lack of long ranged order in the partial pair distribution functions (see Fig.  2) act as a further check that this region is sufficiently damaged compared to the rest of the structure.
These previously produced structures were then doped with helium atoms amounting to 1% of the total atomic number of the original system. These were randomly placed in the existing simulation boxes. This allowed 100,000 He atoms to be sampled when calculating He diffusion in the amorphous and crystalline zircon structures.
In order to calculate the diffusion separately in both the damaged and undamaged regions, and compare fairly between systems, the mean squared displacements (MSDs) of a subsample of the He atoms within the volume of the damaged region (approximately 250-400 at a time depending on temperature) were calculated. The He atoms only contributed to the MSD of the damaged region whilst within the volume, with its initial position within the region being recorded and these coordinates being where the MSD is calculated from. If the He atom were to leave and re-enter the system it would be considered as a new atom entering the subregion, with a new initial position. This means that any He paths outside the damaged region are not included in the calculation of the diffusion. The diffusion coefficients D plotted are calculated during the diffusive regime of the helium atoms where the square of atomic displacement is linear with time.
The Arrhenius equation (Eq. (1)) relates diffusion to temperature and allows calculation of the activation energy for diffusion.  where D is the diffusion constant, D 0 the pre-exponential factor, E a is the activation energy required for diffusion and T is the temperature. One can take the natural logarithm to yield and calculate E a from the slope of ln(D) vs 1/T . A series of MD simulations (8 in the crystal, 6 in the damaged system and 5 in the amorphous system) were run in each of the three systems at a range of temperatures from 1000 K to a maximum temperature of 2400 K in order to calculate the diffusion coefficients and the activation energy required for helium diffusion. For each temperature the systems were initially equilibrated for 100 ps followed by a 100 ps (80 ps for the damaged structure) run to collect statistics using the NPT hoover ensemble.

A. Diffusion in zircon
The Arrhenius plots produced for pure crystal, amorphous and radiation damaged zircon can be seen in Fig. 3. From the gradient of the relevant slopes in the Arrhenius plots a value for the activation energy of helium diffusion for the pure crystal, amorphous and damaged zircon structures are calculated (see Table I ). The R 2 -values quantifying the fit of the Arrhenius plots to the straight line are in the range 0.96-0.99. The slight scatter seen in the damaged region is due to the smaller number of helium atoms sampled (250-400) in the damaged structure as compared to the 100,000 helium atoms sampled for the crystalline and amorphous simulations, because the radiation-damaged region is small compared to the total system size.
As expected the crystal has the smallest activation energy due to the crystallographic pathways being both intact and not blocked at any points, allowing easy diffusion of the helium atoms. The activation energy of He diffusion in the damaged region of the crystalline structure is greater than that in the pure crystal, but is around a fourth of the activation energy required for He diffusion in the completely amorphous structure.
In our previous work 9 we showed that the pair distribution function (PDF) of the most highly overlapped damaged region defined by a regular shape of a box containing 6,000 atoms was very close to the PDF of the fully amorphous structure. In order to study He diffusion in the damaged structure, we have used a larger region damaged by overlapping collision cascades. This region has an irregular shape as shown in Fig. 1a. To make a consistent comparison of the structural correlations in different regions measured by a distribution function, we have selected atoms in the melt-quench and crystalline structure bound by the same surface as in the radiation-damaged region (see Fig. 1b and 1c). We calculated the correlation function as n(r)/r 2 , where n(r) is the number of atoms in a bin, and refer to this function as "pseudo-PDF" because it is calculated for a finite non-periodic collection of atoms where the usual normalisation does not apply.
The pseudo PDFs for all three of those regions are shown in Fig. 2. We observe the decrease and widening of PDF peaks in the disordered structures. We also observe a difference between the radiation-damaged structure and the fully amorphous one. For example, the radiation-structure shows fairly small peaks for O-O, Si-Si, Zr-Si and Zr-Zr sublattices in the medium-range order (e.g. second and third peaks) where no discernible peaks are seen in the fully amorphous structure. This is in agreement with our previous work 9 in which it was shown that even after 6 collision cascades further structural damage and amorphization was still possible. As well as continual change to the local PDFs a continued increase in the number of coordination defects was also observed in the overlapped damaged region. The persistence of PDF peaks at the cations' crystalline peak positions have been also observed in amorphous zirconolite 29 . We will return to this point below when we discuss the variation of activation energies with the degree of introduced disorder.
Our values for the activation energies in crystalline and amorphous zircon (see Table I) are the first we know of from large scale molecular dynamics simulations calculating bulk diffusion from a range of temperatures. As seen in Table  FIG. 1: The atomic structures of a) the amorphised radiation damaged region in the pure crystalline zircon, b) a region in the melt-quench amorphous structure that is the same shape as the damaged region, c) a region in the crystalline zircon also of the same shape. These regions are used to calculate and fairly compare He diffusion between the various systems.
I the calculated activation energy in the perfect crystal is 8.44 ± 0.26 kJ/mol. This is lower than experimentally measured values which are 10-20 times greater 15,16,19,20 . However, the activation energy barriers for He diffusion calculated using a variety of computational methods, including DFT and molecular dynamics, are significantly lower than those measured in experiments. The modelling results also show a large variation of activation energy depending on the direction. For example, DFT nudged elastic band calculations of the activation energy barrier along the c channel in crystalline zircon have reported values of 21.3 kJ/mol 25 and 42 kJ/mol 17 . Similar calculations using DFT parameterized empirical potentials calculate the energy activation barrier along [001] to be 13.4 kJ/mol 24 . Saadoune et al. 25 demonstrate that a single oxygen vacancy can lower the energy barrier for a jump from 21.3 kJ/mol to as little as 5.9 kJ/mol. The caveat is this subsequently acts as a trap for He, increasing the energy required for the next jump to 30 kJ/mol. Although the exact values for the He activation energy differ, most likely due to use of different  Fig. 1 functionals, step sizes and relaxation parameters for the surrounding lattice, they all calculate the activation energy barrier perpendicular to the c-direction (i.e. along [100]) to be much higher, 255 17 -259 25 kJ/mol. In an atomistic simulation study using pair wise potentials instead of DFT Saadoune et al. 30 showed the fastest diffusion pathway [100] to be along tetrahedral interstial sites in the perfect zircon lattice. This corresponds to an energy barrier for helium diffusion as low as 5 kJ/mol. This agrees well with our results, that in crystalline zircon the activation energy is low as diffusion is dominated by He pathways along the c channel, and therefore increases as this channel is destroyed by radiation damage and subsequent complete amorphization.
Other experiments similarly report larger activation energies parallel to the c-channel as compared to simulations 15,16 . In one study 19 where 7 samples are measured, activation energies parallel to the c-direction (in samples that are still crystalline) are reported to be between 138.22 kJ/mol and 166.19 kJ/mol, and those perpendicular to c range from 106.52-169.75 kJ/mol. Two of the samples studied have no clear c direction and are considered amorphized due to long term radiation damage and have activation energies of 145.96 and 70.76 kJ/mol. In another experimental study 20 where He is implanted into polished zircon crystals via ion beam, the measured activation energies along each direction are almost equal, 148 kJ/mol parallel to the c-direction and 146 kJ/mol perpendicular. This discrepancy between experiment and the large number of different simulations carried out using various methodologies and techniques can be attributed to the large effects that even a small number of impurities, defects and radiation damage has on the diffusion of He along the c-channel. In the simulation studies zircon is a perfect crystal, which is impossible in the natural samples used experimentally where the zircon will contain Pb, Th and U impurities and will have experienced some level of radiation damage 19 . The use of ion beams to implant He into zircon also adds the possibility of atomic displacements and defects in the structure 20 . Saadoune et al. 25 show that the introduction of even a single oxygen vacancy along the channel can form a defect trap for He requiring over 30 keV to escape. Recent experimental work using laser depth profiling showed that the presence of said trace elements can have a significant effect on helium diffusion 31 . An interesting experimental observation is that as the dose rate of natural samples goes from 1 × 10 16 α/g to 1 × 10 19 α/g, the resulting activation energy of helium diffusion falls from 166.19 kJ/mol to 70.7 kJ/mol due to the effect of damage percolation 19 . This is due to zircon being rendered totally amorphous at between 1 × 10 18 and 1 × 10 19 α/g 11 . The value calculated from our MD simulations for the fully amorphous structure (114 kJ/mol) is found to be between the activation energies of the two amorphous samples, 19 145.96 kJ/mol and 70.7 kJ/mol respectively, which is in the range expected once zircon has been rendered completely amorphous by radiation damage.
The increase in activation energy due to small amounts of radiation damage or defects are corroborated by the only (to the authors' knowledge) computational work that models the effect of damage on diffusion in zircon 22 . Gautheron et al. 22 use a kinetic Monte-Carlo (KMC) model based on DFT calculations of activation energy and attempt rates where the helium is modelled jumping between interstitial sites. The effect of "damage" is modelled by randomly blocking some percentage of the interstitial sites. Even after blocking only 1 percent of interstitial sites they calculate a large increase in the activation energy from 23 kJ/mol in the perfect crystal to around 60 kJ/mol, around the half value we calculated for the completely amorphous structure and within 4 kJ/mol of the activation energy along the perpendicular directions to c. Additional blocking of sites, up to 20%, does not increase the activation energy calculated using KMC. Whilst we calculated the activation energy in the damaged region after 6 radiation cascades to be 30.87 ±3.93 kJ/mol, half that of the value Gatheron et al. 22 calculated for 'damaged' zircon, consistently we also see a three-fold increased in the value of activation energy in the damaged region compared to the crystal. The lower estimate of the activation energy calculated using MD compared to DFT may include different factors. For example, DFT does not account for temperature and collective modes facilitating diffusion. Another difference between this and previously mentioned modelling studies on one hand and our current simulations on the other hand is that we do not specify a certain direction or pathway of diffusion: our activation energies should be viewed as effective values with contributions from all possible pathways. Nevertheless, the same trends are observed in both studies, that radiation damage significantly increases the activation energy in zircon, even before the complete loss of the cation sublattice.
In the above discussion, decreasing crystallinity results in the disappearance of fast diffusion pathways. This is consistent with (a) our finding that the activation energy is larger in the fully amorphous structure as compared to  radiation-damaged one and (b) our earlier observation that the radiation-damaged structure has more order than the fully amorphous systems: recall that the radiation-structure shows fairly small peaks for O-O, Si-Si, Zr-Si and Zr-Zr sublattices in the medium-range order (e.g. second and third peaks) where no discernable peaks are seen in the fully amorphous structure.
To summarise, the main result of this section is the increase of the activation energy with the amount of introduced disorder.

B. Helium accumulation in damaged zircon
Another effect of radiation damage to consider is the accumulation of helium inside the damaged region. The diffusion rate in the crystalline region is much higher, and as such He will arrive at the damaged region much faster than it will exit, leading to a build up. Fig. 4 shows the number of helium atoms inside the damaged region during the temperature runs, as well as the number density of He in the damaged region, against time. The He number density in the damaged region at each temperature can be compared to the average number density of the whole system at each temperature shown in Table II, where the number density ranges from 0.825-0.848×10 −3 Å −3 in the crystal and 0.694-0.743×10 −3 Å −3 in amorphous zircon. When comparing the number densities of helium atoms in the damaged region to the whole system it is evident there is significantly more helium in that region in Fig. 4 than in both the amorphous and crystalline structures in Table II, and that the number density continues to increase with time. We note that we do not observe saturation of the number of He atoms in the damaged region due to limitations related to simulating large system sizes for long times.
Our finding indicates that during the process of damage and amorphization helium atoms will be mainly trapped in the damaged region and likely remain trapped there for longer due to increased activation energies. This raises the possibility that He will become inhomogeneously stored in zircon, becoming concentrated in regions that are damaged first.
Although from irradiation experiments on zircon no helium bubble formation was reported from electron microscope analysis 32 , it is known that the formation of helium bubbles occurs in other materials studied as potential waste forms. This is particularly apparent in euxenite, were helium bubble formation occurs in high enough concentrations that it can be related to swelling and eventually cracking 14 . Recent work on zirconolite showed that helium bubble formation occurs in the amorphous damaged phase and leaching would be expected to occur in a helium bubble containing zirconolite waste form 33 . Our simulations suggest that He build up occurs in radiation-damaged parts of a waste form, a prediction that can be studied experimentally.

IV. CONCLUSIONS
Using molecular dynamic simulations, we have modelled the diffusion of helium atoms in crystalline and amorphous zircon and radiation damaged crystalline zircon. Examination of the PDFs of the radiation damaged zircon structure shows that there is some retained order in the cation sublattice. This retained 'ordering' plays a large part in the diffusion of helium, leading to a lower activation energy in a radiation damaged zircon region compared to that of a fully amorphous structure generated using the melt-quench method.
We have also for the first time reported activation energys for helium diffusion using molecular dynamics taking into account a 1% doping of helium atoms in a crystal, damaged and amorphous zircon structures, which is consistent with results from different computational methods used previously. We also show good agreement with experimentally measured activation energies in amorphous zircon 22 .
The large difference in diffusion rates and activation energies calculated in crystalline and damaged zircon is further evidenced by the build up of He that we predict in the damaged region via our MD simulations. Understanding how this build up may effect bubble formation or cracking in zircon is of importance for predicting waste form performance in long term.