Evolution of amorphous structure under irradiation: zircon case study
A. Diver, O. Dicks, A. M. Elena, I. T. Todorov, K. Trachenko
EEvolution of amorphous structure under irradiation: zircon case study
A. Diver, O. Dicks and K. Trachenko
School of Physics and Astronomy, Queen Mary University of LondonSchool of Physics and Astronomy, Mile End Road, London, E1 4NS, UK
A. M. Elena, I. T. Todorov
Daresbury Laboratory STFC UKRI, Scientific Computing Department,Keckwick Lane, Daresbury WA4 4AD, Cheshire, UK
The nature of the amorphous state has been notably difficult to ascertain at the microscopic level.In addition to the fundamental importance of understanding the amorphous state, potential changesto amorphous structures as a result of radiation damage have direct implications for the pressingproblem of nuclear waste encapsulation. Here, we develop new methods to identify and quantifythe damage produced by high-energy collision cascades that are applicable to amorphous structuresand perform large-scale molecular dynamics simulations of high-energy collision cascades in a modelzircon system. We find that, whereas the averaged probes of order such as pair distribution functiondo not indicate structural changes, local coordination analysis shows that the amorphous structuresubstantially evolves due to radiation damage. Our analysis shows a correlation between the localstructural changes and enthalpy. Important implications for the long-term storage of nuclear wastefollow from our detection of significant local density inhomogeneities. Although we do not reachthe point of convergence where the changes of the amorphous structure saturate, our results implythat the nature of this new converged amorphous state will be of substantial interest in futureexperimental and modelling work.
I. INTRODUCTION
Research into materials suitable for the encapsulationof nuclear waste is both a societal and scientific challenge.Part of this challenge is how to deal with existing stock-piles of plutonium. Currently the UK has the largestcivilian plutonium inventory in the world, expected tobe 140 tonnes this year . Two solutions being consideredare reprocessing Pu into mixed oxide fuel or immobiliza-tion by encapsulation in a waste form followed by deepgeological disposal.Many materials have been proposed as potential im-mobilization matrices (waste forms), but no matter howresistant any material is to radiation damage after a fewthousand years any initial crystalline waste form with atypical waste loading will undergo a transformation, be-coming completely amorphized long before it is deemedsafe. In fact, it is typically considered that waste formsneed to be an effective barrier for long-lived nuclear wasteencapsulation for up to a 1,000,000 years. For a modelwaste form, zircon, this transition to a fully amorphousstate has been measured experimentally and occurs atdoses between 10 and 10 alpha decay events/g, whichcorresponds to a waste storage time in the order of 1,700years.Previous modeling work has only simulated high-energy radiation damage effects in the crystalline phasesof materials, which only accounts for a small fraction ofthe waste form’s life cycle. Herein lies the key issue; itcannot be assumed that the amorphous state of zirconis itself unchanging when exposed to radiation (a temp-tation for such as an assumption lies in considering thatintroducing more damage does not increase disorder inan already amorphous system). Therefore, it is neces- sary to study radiation damage effects in the long-termamorphous phase in order to understand the long-termperformance of the waste form.In addition to nuclear waste encapsulation, the ques-tion of the evolution of amorphous structure is of funda-mental importance. Understanding the amorphous stateand ascertaining its main characteristics in the short andmedium range order has a long history (see, e.g., Refs. for review) and is developing as new high-quality datacome in. Understanding the “problem of glass” was in-terestingly regarded as the most profound problem inmodern physics . It is remarkable that despite morethan one hundred years of modern research, there remainfundamental open questions regarding the nature of theamorphous state. For example, although the “continu-ous random network” structure proposed in early daysof glass science is consistent with the wealth of experi-mental data, the presence of quasi-crystalline cybotacticgroups can not be ruled out by the current experimentaldata .Here, we pose a fundamental question of whether anamorphous structure produced by liquid quench (glass)can be made “more amorphous” in a sense of reducing thedegree of order in the short- and medium range? We alsoexplore related and more specific questions: how does theamorphous structure evolve during progressive introduc-tion of radiation damage, does the damaged structureconverge to a distinct amorphous state with new short-and medium range order and what is the nature of thisnew amorphous state?The work done on distinguishing in which range crys-talline and amorphous states differ can be built upon toaid this endeavour. Wright suggests that to answer thesequestion we need to look at the composition of struc- a r X i v : . [ c ond - m a t . d i s - nn ] A p r tural units, their connectivity and fluctuations in the lo-cal density of the glass to determine the differences inoverall network topology between crystalline and amor-phous states. The same observables can also be used todistinguish between evolving amorphous states in a sim-ilar fashion.In this paper, we address these questions using molecu-lar dynamics (MD) simulations. We use zircon (ZrSiO )as a case study, for several reasons. First, it “accepts” alldamage and shows very little recovery to the crystallinestate. Second, it is a well characterized. Third, it wasconsidered as a potential waste form for the encapsula-tion of plutonium because of its high durability .Radiation damage experiments in zircon have alsobeen conducted, including analyzing naturally radiationdamaged zircon samples , Pu doped samples over peri-ods of 6.5 years and heavy ion beam experiments in-volving a range of different ions and energies . Somework has been carried out simulating amorphous zirconstructures . These simulations show Zr and Si becomeunder-coordinated, with a corresponding volume expan-sion, when compared to the crystalline phase. Previousradiation damage simulations in crystalline zircon havefocused on structural effects like polymerization of SiO n units and percolation . However, there have been nohigh energy collision cascades been performed in amor-phous zircon.Here, we perform large-scale MD simulations of high-energy overlapping collision cascades in both crystallineand amorphous zircon in order to ascertain and under-stand the difference in response of the two structuresand develop a new algorithm to detect and quantify theradiation-induced defects. We find that an average prop-erty such as pair distribution function does not detectthe differences between the damaged configurations inthe amorphous structures. On the other hand, our newanalysis method based on local coordinations finds thatthe amorphous structure substantially evolves due to ra-diation damage. Our analysis shows a notable correla-tion between the local structural changes and enthalpy.We also find significant local density inhomogeneities inthe amorphous structure due to radiation damage, whichhave important implications for the long-term encapsu-lation of nuclear waste. Although we do not reach thepoint of convergence where the changes of the amorphousstructure saturate, our results open up the way to studythis process and imply that ascertaining the nature ofthis new converged amorphous state will be of substan-tial interest in future experimental and modelling work. II. METHODSA. Molecular dynamics simulation method
The molecular dynamics package
DL_POLY_4 , whichuses domain decomposition to simulate systems withlarge numbers of atoms in parallel , has been used to simulate the propagation of heavy recoils due to alphadecay of Pu in crystalline and amorphous zircon. Wemodel radiation damage to the heavy recoil nucleus ofan alpha decay, specifically U atoms with a 70 keV recoilenergy. During the collision cascade thousands of atomsare displaced in the system, with the permanently dis-placed atoms causing amorphization of the material. Inorder to analyze the extent of these atomic displacementsand their effect on the material structure this work makesuse of new on-the-fly analytical tools we have developedfor DL_POLY_4.10 (pre-release), including time resolvedcoordination analysis. All collision cascades were car-ried out using Archer high-performance computing facil-ity. Each cascade was simulated in a box with 10125000atoms using 1152 processors. Each simulation took about20 hours and resulted in 180 ps of trajectory time.Interatomic potentials that have been previously usedand well tested for zircon have been used in this study.The Zr-O and O-O interactions are modeled using Buck-ingham potentials, and for the Si-O interactions a Morsepotential is used. These pairwise potentials are thenjoined at a short-range with a repulsive ZBL (ZieglerBiersack Littmark) stopping potential. Simulationboxes containing in the order of 10 million atoms wereused for both the amorphous and crystalline structures.These large cell sizes were necessary to contain multipleoverlapping collision cascades. All collision were carriedout using the NVE ensemble. A variable timestep wasused in these simulations, starting at 1 fs, in order toaccount for the initial high velocities of particles.In order to study the effect of multiple cascades over-lapping, each U recoil ion was placed at the same initialstarting coordinates within the simulation box and givenan impulse along the same initial trajectory, with the onlydifference being the initial structure being taken from theprevious cascade. This was done in order to maximise thedamage in one area of the simulation box and to observethe effects of damage saturation on the local structure. B. Producing amorphous zircon
The amorphous structures of zircon were produced bya slow melt quenching of the same crystal structures usedin this work. First the crystal structure was melted at6500 K for 100 ps using an NPT ensemble and then sub-sequently quenched to 300 K using a quench rate of 10K ps − . This quench rate is much higher than those usedin experiments, but a slower quench rate of 1 K ps − on asmaller system of 100,000 atoms resulted in little differ-ence in overall structure. The volume of the amorphousstructure produced by the melt quench was calculated tobe 16% greater than the crystalline zircon which com-pares well with the experimentally measured volume in-crease due to amorphization in zircon of 18% .We note that as the motivation behind creating theamorphous structures comes from the amorphization ofzircon due to collision cascades one might ask why theamorphous structures were not produced by collision cas-cades. Each individual alpha decay event produces in theorder of a few thousand permanently displaced atoms,as shown from nuclear magnetic resonance (NMR) stud-ies of natural radiation damaged zircons and our ownwork. This means that the amount of collision cascadesnecessary to completely amorphize zircon via moleculardynamics would be around 2,000 MD simulations, if weassume approximately 5,000 atoms are displaced per col-lision cascade. Therefore, the amount of time, and asso-ciated computing cost, for achieving an amorphous statevia simulated cascades is not practical. III. RESULTS AND ANALYSISA. Defining damage in an amorphous structure
Previous works quantify the effects of radiation dam-age cascades by measuring the number of displaced or‘defect’ atoms . An atom is considered displaced if itmoves more than a set cut off radius away from its ini-tial position. Some of these displaced atoms may moveto crystalline positions, resulting in damage recovery.Therefore, a defect is considered to occur if there is eitherno atom within a sphere of a given radius centred on aninitial crystalline atomic position (a vacancy defect) orif there is more than one atom inside that same sphere(an interstitial defect). Both definitions of displacementand defect rely on the initial atomic positions, which isreasonable for crystalline systems but is inadequate fordefining defects in glass or amorphous systems. Indeed,in a crystal the atoms in the surrounding lattice to adisplaced atom relax back to their crystalline positions,with small perturbations, and therefore are not countedas defects. However, in an amorphous system, if a bondis broken in a glass forming unit, such as a tetrahedronor octahedron, and an atom is displaced, then the sur-rounding glass forming units may relax back to positionsdisplaced from their initial coordinates, whilst maintain-ing the same connectivity of the glass network. To ensurethat only the truly “displaced” atoms are counted, andnot all of the atoms in the surrounding polyhedra, a newdefinition for defects and displaced atoms must be con-sidered, using the connectivity of the atoms rather thantheir initial positions, due to the lack of crystalline coor-dinates.Here, we develop a different method which measuresthe change in local coordination of each atom due to thecollision cascades. The nearest neighbour coordination ofan atom is defined by cutoffs between different pairs ofatomic species. In this paper the cutoff radius is definedas the first minimum after the first peak in the individualpair distribution functions, g ij (r), of each pair of species.Our on-the-fly analysis allows two different coordina-tion based quantities to be measured. The first one canbe used to track displaced atoms during a cascade. Wedefine an atom to be displaced if there is (a) a change FIG. 1: The total radial distribution function of crystallinezircon in the ’overlap’ region near the initial recoil site, beforeany cascades, and after the 1 st and 4 th cascade. to the specific atoms it is coordinated to, and (b) it isdisplaced by more than a cutoff distance from its initialposition. A Si atom that remains bonded to the sameO atoms that form its tetrahedral cage is not consid-ered a displaced atom, regardless of any translation, as itdoes not lead to a change in the glass network. However,if the O atoms the Si is coordinated with are replacedby different O atoms, even if the coordination number ismaintained, it is considered displaced if it has also moveda minimum distance from its initial position, as the glassconnectivity has changed.The second quantifiable measure is the overall distri-bution of species-species coordination numbers, and howthey change as a result of radiation damage. The coordi-nation distribution of each species A-B pair is calculatedevery n time steps, and can be plotted as a function oftime. This can be used to measure the change in theoverall number of ’defects’ in the system before, during,and after a cascade. The average change in ’defect’ num-ber as the result of a radiation event is also calculatedin this paper. The initial coordination distribution aver-age is calculated during a simulation run of the structurebefore any collision cascades and then compared to theaverage of each coordination distribution over the last 20ps of the cascade simulation run.Both of these methods can be used to monitor timeresolved changes to the structure of the material duringa cascade, and the overall changes to the structure as aresult of a single, or multiple cascades. FIG. 2: Total radial distribution function for amorphous zir-con in the overlapping region before and after cascade 5 inthe region of the cascades
B. Collision cascades
1. Pair distribution functions in both crystalline andamorphous zircon after collision cascades
We firstly analyze the pair distribution functions(PDF) of the crystalline and amorphous structures dam-aged by collision cascades. The PDFs of the total sys-tems (averaged over 10 million atoms) are insensitive toradiation damage, there is no discernible change to thecrystalline or amorphous PDF before and after the col-lision cascades as the structural changes are small whencompared to the total system size.Some use can be made of the PDF if the regions of max-imum damage are zoomed in on. The PDF was calculatedin the structure containing around 6000 atoms and en-compassing the collision cascades. In Fig. 1 we calculatethe PDF peaks of a region of approximately 6000 atomsaffected by the collision cascade. It can be observed thatthe peaks in the PDF of crystalline zircon decrease andsubsequently disappear as more collision cascades over-lap. This is characteristic of the crystalline-amorphoustransition where only short and medium-range order ispresent in the amorphous state. However, the effect ofradiation damage in the amorphous structure in a sim-ilarly sized region (see Fig. 2) is markedly smaller: theonly observable difference is a decrease in the height ofthe first peak of the PDF.The insensitivity of the PDF to local structuralchanges reinforces earlier discussion , and underscoresthe need for more sensitive and local analysis. This willbecome more apparent in our analysis of coordinationchanges and local density changes discussed later.
FIG. 3: Pair distribution functions, g ij ( r ), of Si-O and Zr-Opairs in crystalline zircon before any radiation damage.
2. Coordination statistics in crystalline zircon due tomultiple cascades
For all collision cascades in the crystal structure thecut offs to build local atomic coordination are 2.0 Å and2.80 Å for Si-O and Zr-O respectively, with these valuestaken from the first minimum after the first maximum ofthe pair distribution functions in Fig. 3. These valuesare then used in our coordination displacement definitionwith the specific atom needing to also be displaced by aminimum 0.75 Å, 1.0 Å and 1.0 Å respectively for Si, Oand Zr atoms. The values for the minimum displacementcome from being slightly less than half the bond lengthbetween atoms. For O we choose the largest bond lengthpair to define this.In Fig. 4 the evolution of coordination displacementsduring the initial cascade simulation (cascade 1) and sub-sequent cascades is shown. We observe that overlappingcascades produce a similar number of displacements tothe previous cascades, and that the displacements perrecoil do not increase when a damaged region under-goes further radiation damage. This result is at oddswith previous work on radiation damage in crystallinezirconolite , where a second overlapping cascade leadsto an increase in displaced atoms in the already damagedregion. However, there are many factors which could ex-plain this discrepancy, the two different materials mayhave different radiation resistance properties especiallyas zirconolite contains more mobile alkali ions. The pre-vious study also used simulations cells approximately a10 th of the size of the simulation cells used here, whichmay have lead to larger temperature increases during thecascade.Although each overlapping cascade produces similarnumbers of total displacements, there is a change in howthe system relaxes from the peak number of displace-ments in the first 5 ps after the recoil (see Fig. 4). Previ- FIG. 4: The time evolution of coordination displacements ofZr, O and Si during collision cascades 1, 3 and 5 in crystallinezircon. ous work looking at naturally damaged zircon samplesestimated around 830 Si atoms are displaced per alphadecay event which is in good agreement with the numberof displaced Si atoms we get after the first cascade be-ing around 1300. As more cascades overlap the amountof recovery of displaced O atoms decreases and becomesinsignificant by the 5 th cascade, meaning once a regionbecomes significantly damaged less recovery occurs dueto amorphization, and atoms are unable to return to theircrystalline lattice positions.To see a measurable difference between each cascade,and to quantify changes to the overall microscopic struc-ture, we follow the change in overall coordination statis-tics as seen in Fig. 5 and Fig. 6. Both the Zr-O and Si-Ocoordination changes reach a steady equilibrium within20 ps of the recoil after an initial spike. It is interesting FIG. 5: The change in Zr-O coordination (where Zr is 5-coordinated Zr-O) during consecutive overlapping U recoilcascades in crystalline zircon, as a function of time. to note that the number of displaced atoms recovers to anew equilibrium in much less time, within 1 ps, but therecovery is much smaller proportionally. This suggeststhat while approximately 10,000 atoms become displacedduring a cascade, in the system fewer than 1,500 atomschange their local coordination, or become topological‘defects’.It is also possible to compare the difference betweenthe number of Si atoms that are displaced from Fig. 4and the overall change in Si-O coordination from Fig.6. After the initial cascade there are roughly 1,200 Siatoms that can be considered to be ‘displaced’, but thereis only an increase of 30 3-coordinated Si atoms fromthe original 4. Even after 3 cascades there is only a totalchange of 120 3-coordinated Si atoms due to the damage,which is orders of magnitude lower than the damage sug-gested by only examining the number of displacements,which is regularly reported. Surprisingly the opposite ef-fect is seen in Fig. 5 where the amount of coordinationchanges is greater than the number of displaced Zr atomsas shown in Fig. 4. The number of Zr displacements isfewer than the number of coordination changes. This islikely to be because Zr will undergo a change in coordina-tion when O ions are displaced, whilst the much heavierand larger Zr ion itself does not move from its originalequilibrium position in the cation sublattice and is there-fore not treated as a displacement by our algorithm. Ocan be considered to be more mobile than Si and Zr.It can also be observed that radiation damage decreasesthe number of 8-coordinated Zr in the crystal with eachsubsequent overlapping cascade, as the damaged regionbecomes amorphized and the local average coordinationdecreases.
FIG. 6: Change in Si-O coordination during consecutive col-lision cascades in crystalline zircon.
3. Coordination statistics in amorphous zircon due tomultiple overlapping cascades
In this section, we analyze the radiation damage in theamorphous zircon structure. The change in local atomiccoordination in the amorphous zircon system is shown inFig. 7, with the cut offs chosen to be 2.0 Å and 2.96 Å forSi-O and Zr-O respectively, where once again these valuesare taken from Fig. 8. To be defined as a displacement,the minimum distances the atom must move, as well asundergoing a coordination change, are 0.75 Å, 1.0 Å and1.0 Å respectively for Si, O and Zr atoms. These werethe same values as used in the crystalline case for directcomparison.Interestingly, there is no sharp recovery peak at 1 psin the amorphous structure Fig.7, in contrast to crys-talline zircon. In fact it is very similar to how the 5 th cascade in the crystal shows no peak in the O displace-ments, because the region of the cascade propagation isalready amorphized. In crystalline materials, the sharprecovery peak was attributed to reversible elastic defor-mation of the structure surrounding the collision cascade:the strongly-anharmonic motion in the cascade core re-sults in the volume expansion (similar to thermal expan-sion due to anharmonicity) which exerts tensile stress onthe surrounding lattice . The associated strain displacesatoms in the surrounding lattice, resulting in the brief in-crease of displaced and defect atoms around 1 ps. Thisdisplacement is elastic and reversible, resulting in a sharppeak. The observed absence of the peak in amorphousstructure interestingly implies that the response of theamorphous structure to radiation damage is more slackin the above sense of reversible elastic deformation.Another notable observation is that the peak numberof coordination displacements produced by cascades inthe amorphous system is approximately 4 times greaterthan that in the crystalline phase and shows an increase FIG. 7: Evolution of coordination displacements of Zr, O, Siduring cascades 1,3,5 in an initial amorphous zircon structure. with cascade number from the initial cascade to the3 rd cascade. Hence the amorphous structure is softer,atoms are more easily displaced and it is more suscepti-ble to damage. However from the 3 rd cascade onwardsthis increase in the peak displacements disappears, whichmay be due to the amorphous system produced by melt-quenching differing from one produced by radiation dam-age. Therefore, saturation in the number of displacedatoms in the heavily damaged structure by multiple over-lapping collision cascades is reached.Although the number of displaced atoms per collisionsaturates, local structural changes still evolve. Indeed,measurable differences of local structural changes can beanalyzed by observing the overall change in coordinationstatistics as shown in Fig. 9 and Fig. 10. Again, as withthe crystalline case, the number of Si atoms identifiedas displaced atoms is much higher, roughly 6,000, but FIG. 8: Pair distribution functions, g ij ( r ), of Si-O and Zr-Opairs in amorphous zircon before any radiation damage.FIG. 9: Change in Si-O coordination during overlapping col-lision cascades in amorphous zircon. the number of newly created 3-coordinated Si-O is onlyaround 200. The majority of displaced Si atoms relaxback into a tetrahedral configuration, merely changingtheir nearest neighbour O atoms. After 5 cascades only600 3-coordinated Si-O are produced, orders of magni-tude fewer than the number displaced.There is a greater degree of fluctuation in the Zr-Ocoordination number, with noise between steps predomi-nantly due to the minimum of the g(r) not going to zeroafter the first peak. Although these fluctuations mightseem rather large, they account for less than 0.03% ofall Zr atoms inside the MD box. The fluctuation of co-ordination numbers is due to vibrational motion causingsome atoms to move in and out of the cutoff radius, thusswitching their coordination. However, it can still beclearly seen that, despite the large variations in coordi-nation number during each run, the cascades lead to a FIG. 10: Change in Zr-O coordination during collision cas-cades in amorphous zircon. substantial change in the average coordination number.The long simulation times (180 ps) show that the averageis stable after each individual cascade.An interesting trend in the data is the decrease of 6coordinated Zr-O accompanied by an increase in 7 and 8coordinated Zr-O after each cascade. The decrease in 6coordinated Zr suggests there may be a saturation limitfor the number of coordination changes produced by acollision cascade. As crystalline zircon undergoes a de-crease in average Zr-O coordination, and the melt-quenchamorphous zircon undergoes an increase, we hypothesizethere is a stable coordination distribution between them,where the structure is healed as quickly as it is dam-aged by radiation. This hypothesis is further supportedby the change in 6 coordinated Zr-O in the amorphoussystem after each cascade, decreasing from 1500 initiallyto roughly 500 from the third cascade onwards. Thissuggests a possible saturation limit is being approachedwhere fewer coordination changes can occur and remainstable after collision cascades.The same effect is seen with a decrease in the num-ber of 4-coordinated Si-O after each cascade, both in theabsolute change in number of 4-coordinated Si from theoriginal structure (Fig. 11) and the change between sub-sequent cascades.
4. Other structural and energetic effects
As well as examining the effects of the recoil nucleuson atomic displacements, changes to other structural andmacroscopic properties, like enthalpy, can be studied.We now look at the density variations in the amor-phous structure due to collision cascades and their over-lap. We investigate this effect by splitting the cell up intoboxes of length 10 Å and calculating the density of eachbox before and after each collision cascade. Examinationof different slices along the z-direction of the simulation
FIG. 11: Total decrease in number of 4 coordinated Si-Ocompared to the initial structure after each collision cascadein amorphous zircon.Change in enthalpy Crystal (eV) Amorphous (eV)Cascade 1 2822 ±
27 1795 ± ±
43 1417 ± ±
37 718 ± ±
47 1141 ± ±
62 824 ± ±
68 1284 ± cells, as shown in Figures 12 and 13, shows that radia-tion damage results in increased density inhomogeneitiesof up to 30% of the initial regions density. The regionsof density inhomogeneity grow in size with the cascadeoverlap as seen in Fig. 13.Local density inhomogeneities seen as the presence ofregions of increased and reduced densities (relative tothe undamaged structure) have important implicationsfor the storage of nuclear waste. Regions with decreaseddensity serve as fast diffusion pathways in the waste form.At large radiation doses, these regions will connect intopercolating clusters where enhanced diffusion is observedon the macroscopic scale .We note that density inhomogeneities at the nanoscalehave been reported before from small-angle x-ray scatter-ing experiments on natural, radiation-damaged zircon ,as well as from previous MD simulations of close overlapsin crystalline zircon . However, the effect is much morepronounced in amorphous zircon, as expected due to theincreased diffusion rates of radiation damaged materialscompared to their undamaged phases .We now address the energetic effects of collision cas-cades. We calculate the enthalpy, a macroscopically mea-surable quantity that can be directly compared with ex-periment. The enthalpy of formation of the amorphousphase compared to the crystal phase of zircon is calcu- FIG. 12: The change in density of amorphous zircon after the1 st collision cascade. As shown by the scale, the more red aregion, the larger the increase in percentage change in densityand the bluer a region the larger the decrease in percentagechange in density. lated from our simulations to be 109 kJ mole − . Thiscan be compared with experimentally measured valuesof metamict zircon, natural samples of zircon with par-tial radiation damage over hundreds of millions of years,where the enthalpy is 59 kJ mole − . This is a similar or-der of magnitude, with potential differences attributableto the variations in structure between a fully amorphizedsample in our simulations and a partially irradiated sam-ple which has also undergone relaxation, as well as theeffects of the specific interatomic potentials used.More specifically, we can quantify the change in en-thalpy between each collision cascade in both the crystaland amorphous systems, as shown in TABLE I. Inter-estingly, although we give the impacting atom 70 keV ofkinetic energy, there is only a total change of 2-3 keV inthe enthalpy of the crystalline system. An even smallerchange in the enthalpy of the amorphous system is ob-served, 0.8-1.7 keV per cascade. The difference in theenthalpy changes between the two systems can be ex- FIG. 13: The change in density of amorphous zircon after the5 th collison cascade. As shown by the scale, the more red aregion, the larger the increase in percentage change in densityand the bluer a region the larger the decrease in percentagechange in density. plained by the difference in the initial entropy of thecrystalline system, which is in the global minimum offree energy, compared to the amorphous system. In crys-talline zircon radiation damage is producing some co-ordination changes, but most of the entropy change isattributed to amorphizing the connectivity of the crys-talline structural units, rather than breaking them. Inthe amorphous structure the majority of the change instructural entropy comes from the breaking of the Si-O,and to a lesser extent the Zr-O, structural units, as thesystem is already fully amorphized.The effect of the difference between breaking struc-tural units and amorphizing the structure is shown inFig. 14. We visualize the position of all the newly cre-ated 3-fold coordinated Si atoms as a result of 6 over-lapping cascades in the amorphous and crystalline zir-con structure. In the crystal, 160 new 3-coordinated Siatoms are created and are predominantly concentratedin the region along the initial trajectory of the U recoil a)b) FIG. 14: Positions of all newly created 3 coordinated Si-Oatoms inside the simulation cell for a) crystalline and b) amor-phous zircon after 6 overlapping collision cascades. nucleus. In the amorphous structure, 640 3-coordinatedSi atoms are produced, suggesting that the lack of crys-tal lattice inhibits the rate of recovery of damaged Si-Otetrahedra. This shows that in the crystalline case thecascades mainly contribute to amorphizing the structure,while producing some breaking of structural units, but inthe amorphous structure cascades instead lead to an in-crease in the breaking of network structural units. Evenin the amorphous case the density of 3-coordinated Si isnot high, demonstrating the difficulty in creating a fullysaturated damaged region in the material.Interestingly, we observe that both the changes in thenumber of 4-coordinated Si-O and the enthalpy of theamorphous structure (see Fig. 15) are highly correlatedbetween subsequent cascades, adding further evidencethat the enthalpy change in amorphous zircon is linkedto the breaking of Si-O bonds. Further, both enthalpyand Si-O coordination are measurable using experimentaltechniques which would allow direct verification of thisprediction. This implies that the macroscopic changein enthalpy could be directly linked to an atomic levelchange in the zircon structure.An important conclusion from our analysis is that thelocal structures in the amorphous system substantially0 a)b)
FIG. 15: Change in both a) the number of 4 coordinated Si-Oand b) the system enthalpy between each collision cascade inamorphous zircon. evolve with each collision cascade but these changes can-not be observed in the pair distribution functions of thematerial where minimal to no variation is observed. Thisresult implies that amorphous structures can significantlyevolve at the local structural level. Apart from the insightinto the nature of the amorphous state itself , this hasimportant consequences for practical applications, suchas nuclear waste storage due to density inhomogeneities.We have not reached a point in our simulations wherethe amorphous structure saturates or converges to anamorphous structure which no longer changes in responseto radiation damage, measured using local coordinationnumbers. This is due to limitations related to compu-tational costs. However, one expects that such a con- vergence eventually takes places due to the need to fulfillchemical and connectivity constraints at the local level .The nature of this converged amorphous state will be ofsubstantial interest in future experimental and modellingeffort, in view of continuing changes of coordination anddensity at the nanoscale. IV. SUMMARY
In summary, we have proposed a new method to detectand quantify damage in amorphous structures. Combin-ing this method and large-scale molecular dynamics sim-ulations in zircon, we have arrived at several insights re-garding the response of crystalline and amorphous struc-tures to radiation damage. Whereas an average, macro-scopic property of the material, such as the pair distri-bution function, does not detect the differences betweenthe damaged configurations in the amorphous structures,our new analysis method based on local atomic coordi-nation finds that the amorphous structure substantially evolves as a result of radiation damage. Our analysis alsoshows a notable correlation between the local structuralchanges and changes in enthalpy. Importantly for thelong-term immobilization of nuclear waste, we find sig-nificant local density inhomogeneities in the amorphousstructure due to radiation damage. Although we do notreach the point of convergence where the changes of theamorphous structure saturate, we have proposed that ourresults open up the way to study this process and im-ply that ascertaining the nature of this new convergedamorphous state will be of substantial interest in futureexperimental and modelling work.
V. ACKNOWLEDGEMENTS N. C. Hyatt, Energy Policy , 303 (2017). W. J. Weber, R. C. Ewing, C. R. A. Catlow, T. D. D. L. Rubia, L. W. Hobbs, C. Kinoshita, H. Matzke, A. T.Motta, M. Nastasi, E. K. H. Salje, and et al., Journal of Materials Research , 1434–1484 (1998). R. Zallen, The physics of amorphous solids (John Wiley &Sons, 2008). A. C. Wright, International Journal of Applied Glass Sci-ence , 31 (2014). P. W. Anderson, Meeting of the American Physical Society,Indianapolis (1992). A. C. Wright, International Journal of Applied Glass Sci-ence https://doi.org/10.1111/ijag.14630 (2019). R. Ewing, W. Lutze, and W. J. Weber, Journal of MaterialsResearch , 243–246 (1995). T. Murakami, B. C. Chakoumakos, R. C. Ewing, G. R.Lumpkin, and W. J. Weber, American Mineralogist ,1510 (1991). W. J. Weber, R. C. Ewing, and L.-M. Wang, Journal ofMaterials Research , 688–698 (1994). J. Yu, R. Devanathan, and W. J. Weber, J. Mater. Chem. , 3923 (2009). K. Trachenko, M. T. Dove, and E. K. H. Salje, Phys. Rev.B , 180102 (2002). K. Trachenko, M. T. Dove, T. Geisler, I. Todorov, andB. Smith, Journal of Physics: Condensed Matter , S2623(2004). I. T. Todorov, W. Smith, K. Trachenko, and M. T. Dove, Journal of Materials Chemistry , 1911 (2006). J. F. Ziegler and J. P. Biersack, Treatise on Heavy-IonScience , 93–129 (1985). I. Farnan and E. K. H. Salje, Journal of Applied Physics , 2084–2090 (2001). K. Trachenko, M. T. Dove, E. Artacho, I. T. Todorov, andW. Smith, Phys. Rev. B , 174207 (2006). H. F. Chappell, M. T. Dove, K. Trachenko, R. E. A. McK-night, M. A. Carpenter, and S. A. T. Redfern, Journal ofPhysics: Condensed Matter , 055401 (2012). I. Farnan, H. Cho, and W. J. Weber, Nature , 190–193(2007). T. Geisler, K. Trachenko, S. Ríos, M. T. Dove, and E. K. H.Salje, Journal of Physics: Condensed Matter , L597(2003). S. Ríos and E. Salje, Applied Physics Letters , 2061(2004). W. J. Weber and R. C. Ewing, MRS Proceedings ,JJ3.1 (2002). S. Ellsworth, A. Navrotsky, and R. C. Ewing, Physics andChemistry of Minerals , 140 (1994). J. C. Phillips, Journal of Non-Crystalline Solids34