Iron phosphate glasses: Structure determination and radiation tolerance Nuclear Instruments and Methods in Physics Research B

Iron phosphate glass (IPG) has gained recent interest for use in encapsulating radioactive waste for long term storage. In this work, we investigate 5 different compositions of iron phosphate glass. We consider amorphous structures of 3 known crystalline phases: Fe 2+ Fe 3 þ 2 (P 2 O 7 ) 2 , Fe 3 þ 4 (P 2 O 7 ) 3 and Fe 3+ (PO 3 ) 3 , and structures of IPG (40 mol% Fe 2 O 3 and 60 mol% P 2 O 5 ), with 4% and 17% Fe 2+ ion concentrations. Using constant volume molecular dynamics (MD), we quench a set of structures for each glass composition, to ﬁnd the optimal density structure. We found that the lowest energy structures of IPG with 4% and 17% concentration of Fe 2+ , have a density of 3.25 and 3.28 g/cm 3 respectively. This is slightly higher than the experimentally measured values of 2.9 and 2.95 g/cm 3 respectively. We also estimate an upper and lower bound on the melting temperatures of each glass, then for each glass, we simulate radiation damage cascades at 4 keV. The cascade structures can be in the form of either a concentrated thermal spike or more diffuse with sub-cascade branching. We found that the glass compositions with a higher Fe/P atomic ratio, contained a greater number of displacements after the cascade. We also found that the IPG with 4% Fe 2+ , contained slightly fewer displacements than the IPG with 17% Fe 2+ . This is consistent with our previous work, which showed that the threshold displacement energies are lower for glasses with a lower Fe 2+ content. In all the simulations, many PO 4 polyhedra are destroyed during the early stages of irradiation, but recover strongly over a time scale of picoseconds, leaving very few over or under co-ordinated P atoms at the end of the ballistic phase. This is in contrast to recent work in apatite. The strong recovery indicates that phosphate glasses with a low Fe 2+ content could be good materials for waste encapsulation. (cid:1) 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
One leading method for long term storage of high-level nuclear wastes (HLW), is to encapsulate it within a glass matrix [1,2]. Borosilicate glasses have previously been studied by us and others [3][4][5], however, there is also interest in iron phosphate glass (IPG) [6][7][8] for use in immobilisation of HLW.
In this work, we consider 5 different compositions of iron phosphate glass; these are summarised in Table 1. The first 3 compositions have the same composition as the known crystalline phases: Fe 2+ Fe 3þ 2 (P 2 O 7 ) 2 , Fe 3þ 4 (P 2 O 7 ) 3 and Fe 3+ (PO 3 ) 3 . We also consider 2 compositions of IPG with 4% and 17% concentration of Fe 2+ . These compositions are comparable to experimentally prepared IPG [9,10].

Methodology
All computer simulations are performed using molecular dynamics (MD). We used our own MD code (LBOMD), which has been developed over many years at Loughborough, see for example [11] for oxides, and [5] for borosilicate glasses.
We use an existing inter-atomic potential from the literature originally developed by Bushra Al-Hasni and Gavin Mountjoy [12] and modified by us [13] for use in high energy radiation damage simulations. A brief overview is presented here, while full details can be found in the original papers [12,13].
The potential contains both two and three-body terms. The two-body interactions are modelled by the Buckingham potential with fixed partial charges (Eq. (1)). This equation also includes a Coulomb term, which models the long-range ionic interactions.
Here, we use the subscripts, i and j, to refer to each ion. r ij is the inter-ionic distance and 0 is the permittivity of free space. The  Table 2. The two-body interactions between the positive ions (Fe 2+ , Fe 3+ and P) are purely Coulombic. q i and q j are the ion charges. Throughout this work we used fixed partial charges for each ion given by: q P ¼ 3:0, q Fe 2þ ¼ 1:2, q Fe 3þ ¼ 1:8 and q O ¼ À1:2. During an MD simulation, the two iron species (Fe 2+ and Fe 3+ ) cannot change their oxidation state. Therefore, the Fe 2+ /Fe total atomic ratio remains constant. The three-body interactions (within the O-P-O and P-O-P triplets), are modelled by the Stillinger-Weber potential (Eq. (2)) [14].
For an ion triplet (j-i-k), H jik is the bond angle formed at the central ion i, by ions j and k. r ij and r ik are the two internal ionic separations. r c is the cut-off radius beyond which the three-body terms are set to zero. The parameters: k, c, H 0 and r c for each triplet, are given in Table 3. Short range interactions, which are important during radiation damage cascades, were modelled with the ZBL potential [15]. The ZBL potential was joined to the Buckingham potential with a spline function given by Eq. (3). The constants a 0 -a 5 are set such that the potential energy and its first and second derivatives are continuous. These parameters for each ionic bond are given in the original paper [13].

Modelling of glass formation
All the structures are amorphous, there is short range order due to the local bonding, but no long range order. We generate structures for our simulations in the same way as experimental glass is prepared; We create a molten mixture of the required composition, and then slowly quench the system down to room temperature. We quench the system for 1.2 ns using a quench rate of 5 Â 10 12 K=s. To obtain the final 0 K lattices used in this work, we minimise the system using the conjugate gradient method.
This method is limited to lattice sizes of around 3000 atoms, due to the long MD simulation time required (1.2 ns). To quench a system containing 2000 atoms, takes 9 days on the latest 4 GHz Intel Core i7 processor. To generate large system sizes, we first quench lattices of around 2000 atoms, and stack multiple copies to create a super-cell of the desired size.

Density
For each glass composition, we quench a set of structures to obtain 0 K lattices using constant volume MD. We vary the lattice dimensions to produce structures of varying density, from 2.7 to 3.9 g/cm 3 . For the structures with the crystalline composition (Fe 2+ Fe 3þ 2 (P 2 O 7 ) 2 , Fe 3þ 4 (P 2 O 7 ) 3 and Fe 3+ (PO 3 ) 3 ), we used a lattice size of 168, 248 and 260 atoms respectively. However, for the IPG systems with 4% and 17% Fe 2+ , we used much larger system sizes of 2744 and 2380 atoms respectively. This was due to the requirement of a small number of Fe 2+ ions within the structure and the need to maintain charge neutrality of the system. Fig. 1 shows the average total lattice energy per atom at 0 K, plotted against density for each glass. The curves have been offset on the y-axis to allow better viewing within the figure, however, their relative order is preserved. For the crystalline compositions, each data point is an average of 30 quenched structures, however, only one set of quenches were performed for the IPG structures. We fit a 2nd order polynomial to the data to estimate the density of the optimal structure for each glass. The values of the optimal Table 1 The atomic composition of each glass studied in this work.

Glass
Atomic composition % Fe 2+ /Fe total atomic ratio Fe/P atomic ratio  Table 2 The parameters for the two-body interactions used in this work.  Table 3 The parameters for the Stillinger-Weber potential used in this work.  1. For each glass, we plot the average total lattice energy per atom of a set of quenched 0 K lattices for a range of glass densities. The data has been offset on the y-axis to enable viewing of all curves within the same figure, however, the order of the curves has been preserved.
density and the corresponding minimum energy per atom for each glass, are shown in Table 4. We note that these densities are calculated for the 0 K structures of each glass. Since the linear thermal expansion coefficient of glass at room temperature is typically in the range 3 À 10 Â 10 À6 =K, we expect a density change of 1% over the range 0-300 K (assuming a constant thermal expansion coefficient). We also expect little variation in the other glass properties. We see from Table 4, that both the optimal density and the corresponding minimum energy per atom of the 0 K structures, increases with increasing Fe/P atomic ratio. The optimal density found for the IPG systems with 4% and 17% Fe 2+ , was 3.25 and 3.28 g/cm 3 respectively. These values are slightly higher than the experimentally reported values of 2.9 g/cm 3 for IPG with 4% Fe 2+ [16] and 2.95 g/cm 3 for IPG with 17% Fe 2+ [9]. The structures for all glass compositions studied, appear similar to those shown. In Fig. 2a, we draw all phosphorus-oxygen and iron-oxygen bonds. We see the structure is amorphous, and that the iron atoms are distributed evenly throughout the lattice.

Structural properties
For all the glass compositions we studied, we found that all the phosphorus atoms had 4 neighbouring oxygen atoms, and hence, formed PO 4 tetrahedra. This is in contrast to a 1.9% concentration of 5-fold coordinated phosphorus found previously [12]. Here, the authors used a quench rate of 10 13 K=s, which is twice as fast our procedure. Fig. 2b shows only the phosphorus-oxygen bonds, to enable better viewing of the PO 4 tetrahedra. Here, we see that there are bridging oxygen atoms that are shared by two PO 4 tetrahedra. We found that all the glass structures contained mostly long chains of PO 4 tetrahedra, however, all quenched structures also contained a few isolated PO 4 tetrahedra (an example can be seen in Fig. 2c).
In Table 4, we report some mechanical properties of each optimal glass structure. The values of the Young's moduli and bulk moduli were computed using GULP [17]. For the bulk modulus, we report the average of the three methods: Reuss, Voigt and Hill,  From Table 4, we see that both the bulk modulus and Young's modulus slightly decrease with increasing Fe/P atomic ratio. The values of the Young's moduli, E, and bulk moduli, K B , reported, were found to be higher than the experimental values of 70 and 47 GPa respectively [18,19]. However, if the structures with the lower experimentally observed density are examined, the calculated elastic properties are in good agreement with the experimental values [13].

Estimation of the melting point
We now perform some simulations to estimate the melting temperatures of the structures with the optimal density. For each glass, we thermalise a small system at a constant temperature for 10 ps, then run an NVE MD simulation for around 100 ps. We then compute the time average of the temperature and total lattice energy of the whole system. We do this for a range of temperatures from 0 to 6000 K. In Fig. 3a, we plot the total lattice energy against temperature. We can clearly see that above 3000 K, the plot is linear, indicating that the system is fully melted above 3000 K.
The point at which the system departs linear behaviour, indicates that the system is no longer exploring the entire phase space available, and that it is becoming more viscous. We take the point, where the difference between the linear fits, and the curves in Fig. 3a exceeds 0.01 eV, as the upper bound on the melting temperature. This can also be seen in the plot of the specific heat capacity against temperature shown in Fig. 3b, where the curves plateau. We see from the table in Fig. 3c, that this upper bound melting temperature decreases with increasing Fe/P atomic ratio.
We also estimate a lower bound of the melting temperature as the point that the specific heat begins to increase rapidly. These values are also shown in the table in Fig. 3c. The lower bound melt- ing temperatures of IPG are found to be around 1800 K, this is greater than the experimental value of 1198 K [20].

Radiation damage cascades at 4 keV
For each glass, we generate a large system of 450,000 atoms using the methods above. All MD simulations begin with a 0 K lattice. We perform radiation damage cascade simulations, by choosing a primary knock-on atom (PKA), and injecting 4 keV of kinetic energy in a random direction. We run each simulation for 4 ps using constant energy and volume MD (NVE ensemble). This was found to be sufficient time for the systems to reach a steady state in all cases. Since we are using the NVE ensemble, during the simulation the temperature rises to 32 K. Fig. 4 shows plots of the number of displacements of each atomic species against simulation time for the 3 glass structures with the crystalline composition. Here, we consider an atom as displaced if it has moved at least 1.2 Å compared to its original position. We clearly see that in each case, there were fewer displacements when choosing a phosphorus atom as the PKA compared to a Fe 3+ atom. For the Fe 2+ Fe 3þ 2 (P 2 O 7 ) 2 glass, which contains both Fe 2+ and Fe 3+ atoms, the displacement cascades were very similar for the Fe 2+ and Fe 3+ PKAs, as can be seen in Fig. 4a.
The Fe 2+ Fe 3þ 2 (P 2 O 7 ) 2 glass contains twice as many Fe 3+ atoms as Fe 2+ , yet there were the same number of displacements of each atom after a cascade. This means that, a greater fraction of Fe 2+ atoms were displaced compared to Fe 3+ . This was observed for all glass compositions that contain both iron species. Therefore, it appears that the Fe 2+ atoms are less resistant to radiation damage than Fe 3+ atoms. This agrees with the findings of our previous work [13], where we found that the average threshold energy for displacement of Fe 2+ atoms (21 eV) was significantly lower than it was for the Fe 3+ atoms (36 eV). Fig. 5 shows the total number of displacements for each glass during a 4 keV cascade using a phosphorus atom as the PKA. We see that the number of displacements after the cascade, generally decreases with decreasing Fe/P atomic ratio. The Fe/P atomic ratio of the IPG compositions and the Fe 3þ 4 (P 2 O 7 ) 3 blend are similar, while the fraction of Fe 2+ atoms varies. Therefore, we consider these cascades to investigate the effect of the Fe 2+ concentration. While the cascades in IPG with 4% Fe 2+ had fewer total displacements than the IPG with 17% Fe 2+ , a greater number of displacements were observed in the Fe 3þ 4 (P 2 O 7 ) 3 glass, which contains no Fe 2+ atoms.
Since these glass structures are amorphous, we found that the best method for identifying defects, is to consider the coordination number of the atoms. Any phosphorus atom that does not have 4 neighbouring oxygen atoms, is considered to be a defect. Fig. 6,  shows a plot of the total number of phosphorus atoms in each coordination state against simulation time, after a 4 keV cascade for the Fe 2+ Fe 3þ 2 (P 2 O 7 ) 2 glass. We found that all the glass compositions we studied, exhibited the same behaviour. There is a large spike of up to 45 atoms with 3-fold coordination, which decays to zero within 1 ps. There are also smaller spikes of under coordinated atoms, which decay much faster. A small spike of up to 15 atoms with 5-fold coordination occurs in all glass compositions. This typically decays to less than 5 within 300 fs. At the end of the simulation there were no under coordinated phosphorus atoms remaining, however, for some cascades, a maximum of four 5-fold coordinated P atoms would sometimes be observed. Fig. 7 shows two images of the final state of a radiation damage simulation. Fig. 7a shows the typical behaviour, where the PKA is displaced by up to 150 Å and causes isolated regions of damage along its trajectory. This image also shows that a secondary collision occurred, leading to two distinct branches of displaced atoms. We typically see less than 10 atoms with large displacement of up to 150 Å while the majority of atoms move less than 10 Å.
In some cases, the PKA loses its energy in collisions over a much shorter range (within 50 Å), resulting in one or two large regions of damage. An example is shown in Fig. 7b. Generally, for all glass compositions studied, the lighter oxygen and phosphorus PKAs, result in mostly long range damage trajectories (>100 Å), while iron PKAs cause shorter damage trajectories (<100 Å).

Conclusion
For each glass composition, we performed a series of quenching simulations to generate amorphous structures for a range of glass densities. For the IPG systems, the structures with the lowest lattice energy, had densities that exceed the experimental values by 12%.
For the IPG structures, reducing the fraction of Fe 2+ caused fewer displacements after the 4 keV cascades, however, the Fe 3þ 4 (P 2 O 7 ) 3 glass, which contains no Fe 2+ atoms, had a greater number of displacements. However, we found that the fraction of Fe 2+ had no effect on the amount of defects remaining in the system after a radiation damage cascade. In general, the cascades leave very few defects in all the glass compositions, as measured by the number of non 4-fold coordinated P atoms. This is in contrast to cascades performed in fluorapatite [21], where the authors find around 60 PO 5 polyhedra remaining at the end of their 1 keV cascade simulations. We believe this occurs due to the formation of chains of PO 4 tetrahedra, in which two phosphorus atoms share a bridging oxygen atom. Since all phosphorus atoms begin as isolated PO 4 tetrahedra in this system, a liberated oxygen atom must be accommodated by the lattice. In fluorapatite, the authors found that the calcium meta-prisms remained stable, and so the liberated oxygen atom has nowhere else within lattice to go.
Whereas in our case the initial configuration contains mostly chains of PO 4 tetrahedra with very few isolated groups. Furthermore, the iron atoms can be displaced and change their coordination to accommodate varying local oxygen concentration. This allows IPG to recover far more defects than fluorapatite.